基于数值模型的肥城市地下水源地保护区划分

王晓玮，王延岭，赵志伟，陈伟清

(1.中国地质大学(北京)水资源与环境学院，北京 100083；2.山东省第五地质矿产勘查院，山东 泰安 271000)



基于数值模型的肥城市地下水源地保护区划分

王晓玮1，王延岭2，赵志伟2，陈伟清2

(1.中国地质大学(北京)水资源与环境学院，北京 100083；2.山东省第五地质矿产勘查院，山东 泰安 271000)

在北方隐伏岩溶分布地区，许多城市开发利用岩溶地下水资源作为饮用水水源，但利用数值模型方法进行岩溶水饮用水源保护区划分的工作较少。以山东省肥城市城区水源地为例，采用数值模型方法，构建地下水流数值模型，并采用质点追踪技术对岩溶水饮用水源保护区进行了划分。一级保护区面积1.53km2，二级保护区面积28.4km2，准保护区面积约272.9km2，全部位于肥城市辖区内。此次划定的保护区范围比之前划定的保护区范围更合理、准确，有利于加强地下水资源的管理和保护。

岩溶水；饮用水源保护区；数值模型；质点追踪；肥城市

0 引言

北方地区广泛分布的隐伏岩溶地下水水资源是很多城市重要甚至唯一的饮用水水源。由于岩溶含水介质的隐蔽性和分布的复杂性，污染的途径和运移时间长，一旦污染不易被及时发现[1]。因此该类水源地地下水资源的管理应以预防为主，划定合理的保护区范围，建立完善的监督和管理措施。根据国家环境保护总局发布的《饮用水水源保护区划分技术规范》(HJ/T338-2007)[2]，岩溶水饮用水水源保护区的划分办法包括经验公式法和数值模型法。与经验公式法相比，数值模型方法可以刻画复杂的水文地质条件，且能够较为真实的还原地下水系统的特征和运动状态[3]，是较为科学的水源地保护区划分方法[4-6]，通过数值模拟刻画岩溶地下水流是可行的[7-9]，但在对北方岩溶地下水源地保护区划分的研究中仅有少数采用了数值模型方法[10-12]。

肥城盆地是北方岩溶地区典型的向斜-盆地型岩溶水系统[13]，盆地内肥城市城区以地下水作为唯一的供水水源地，之前未针对饮用地下水水源地开展过数值模拟研究[14-15]，只有有关部门曾采用经验法对水源地保护区范围进行了划分[16]。采用数值模型方法，尝试建立肥城市城区水源地三维非稳定流地下水数值模型，通过质点追踪技术划定水源地开采井保护区范围，对有效防止岩溶水源污染，保障城市供水安全具有借鉴意义，可以促进水源地保护和城市建设共同发展。

1 水源地概况

1.1 地质及水文地质条件

肥城盆地位于山东省中部偏西，泰山西麓，肥城市北部，是以肥城市北部平原及其周边山体为主体形成的独立地质构造单元，在区域构造上隶属华北型沉降。区域内在原有的变质岩基底上沉积了浅海相、陆相和海陆交互相的地层。肥城盆地内断裂构造发育，肥城断裂、马山断裂、桃园断裂及其伴生的网状断层，使得大地构造更加复杂。

肥城市城区水源地位于肥城断陷盆地东部(图1)。属暖温带大陆性半湿润季风气候区，四季分明，多年平均降水量633mm，汛期降水量占全年降水量的70%以上，多年平均气温为12.9℃。

1—行政界；2—肥城盆地边界；3—断层；4—模拟第一含水层范围：5—模拟第二含水层范围；6—城区位置；7—水源地开采井图1 模拟区位置图

肥城市城区水源地主要取水层为碳酸盐岩类裂隙岩溶水含水岩组，包括碳酸盐岩裂隙岩溶水含水层和碳酸盐岩夹碎屑岩岩溶裂隙水含水层。碳酸盐岩类裂隙岩溶水分布在水源地大部分地区，在城区一带隐伏于第四系以下，在南部地区则出露于地表，岩性为奥陶系马家沟群、寒武系三山子组灰岩、白云质灰岩，单井涌水量一般1000～3000m3/d；碳酸盐岩类夹碎屑岩岩溶裂隙水分布在水源地南部丘陵区以南，岩性为寒武系炒米店组、崮山组页岩夹薄板状灰岩、竹叶状灰岩，单井涌水量一般500～1000m3/d。含水层主要接受大气降水入渗补给和侧向渗流补给；水源地地下水由南向北径流至北部受石炭二叠系弱透水层的阻挡后折向西南；人工开采已成为水源地范围内主要的岩溶水排泄方式。水源地范围内发育有马山断裂、桃园断裂等倾向E—SE的高角度正断层，多具导水性，对地下水含水系统起一定的控制作用。

1.2 地下水开发利用情况

肥城市城区内居民密集，工业发达。肥城市自来水公司在城区水源地范围内现有生产井和备用井32眼，企业自备井数以百计。根据肥城市环境保护局统计结果，目前肥城市城区水源地工业和生活用自来水集中供水开采量约为3.0万m3/d，而2012年工业自备井取水量则达到了1.67万m3/d。区内地下水动态类型为降水入渗-开采型，近年来年际水位变化总体呈下降趋势，在城区范围内产生了大范围的降落漏斗。2015年肥城市城区水源地地下水总体质量均可达到《地下水质量标准》(GB/T14848-1993)的Ⅲ类水以上标准，但与往年相比总体水质呈下降趋势，特别是硝酸盐、硫酸盐等离子含量一直呈上升趋势。水源地主要污染源包括工业、农业、生活、固体废弃物、加油站及油库等，漏斗区的形成加重了地下水污染[17]。肥城市曾于2000年开展了城区水源地保护区划，并于2001年发布实施了《肥城市城区水源地保护区管理办法》，但城市化进程的加快，经济社会的持续发展，必然导致水资源需求量的增长，进一步加剧了水资源供需矛盾。

2 数值模型方法的应用

2.1 水文地质概念模型与数学模型

该次研究划定的模拟区范围包括肥城市城区水源地补给和径流区，东北和东部至肥城-王庄弧形断裂，南至地表分水岭以南寒武系灰岩与泰山岩群地层界线，北至肥城煤田边界，西南至孙牛公路，西部以行政区划为界。

垂向上将模拟区内含水层概化为两个含水层，第一含水层为灰岩裸露区及第四系，概化为潜水，第二含水层为隐伏灰岩，概化为承压水。一层北部山前冲洪积平原概化为定流量流入边界；西部沿河道方向为定流量流出边界；东部和南部为天然隔水边界；西南部根据流场，判断为隔水边界。二层东部边界定为肥城盆地的东部零流量边界；南部边界处主要为寒武系灰岩，赋水性较差，岩溶连通性一般，概化为零流量边界；中北部边界为老城镇的煤田南部边界，概化为定流量流出边界，模型中分为3段；西部边界本身为行政区边界及其延伸，且为径流带的流出边界，概化为定流量边界。含水层之间存在一定的水量交换。

整个地下水系统接受大气降水入渗、农业回灌、侧向流入的补给，主要排泄方式是人工开采，以及少量的侧向流出和排泄到地表河流。将模拟区含水系统概化为非均质、垂向各项异性、三维非稳定地下水流系统。其控制方程为：

潜水：

式中：Kx，Ky，Kz分别为潜水含水层水平x，y方向和垂向渗透系数(m/d)，μ为给水度，h为含水层的水位标高(m)，ε为含水层的源汇项(L/d)。

承压水：

式中：Kx，Ky，Kz—分别为承压含水层水平x，y方向和垂向渗透系数(m/d)，S为贮水系数(L/m)，H为含水层的水位标高(m)，其他参数意义同前。

2.2 模型的识别和验证

模拟区面积357.8km2，剖分为100m×100m的单元格，选取2011年6月1日至2015年3月31日为该次模型研究的识别阶段，以2015年4月1日至5月31日为验证阶段。以模拟区2011年6月的地下水位实测值为二层承压水模拟水位的初始数据，以2012年6月的地下水位实测值近似替代一层潜水模拟水位的初始数据。

采用2014年4月19日的Landsat8遥感数据，使用ENVI5.1软件进行解译，采用基于CART的自动决策树分类方法，确定了工作区范围内的地物类型，结合岩性特征，综合确定了一层的参数分区，分为12个区。二层参数分区主要考虑不同岩性区域的赋水性差异、区内马山断裂构造的影响等，分为5个区。

该次模拟的源汇项主要包括降水入渗、农业灌溉开采、企业自备井开采和水源地开采。将月降雨量转化为平均降雨入渗强度，根据前述遥感解译结果和地形资料，将降水入渗补给的范围，根据补给强度的差异分为4个区。区内农灌井的数量众多，无详细井位和开采量数据，因此确定4个主要的井灌区，按照200m×200m间距布置开采井，根据《山东省主要农作物灌溉定额》(DB37/T1640-2010)中的规定，核算不同降水保证率下亩均灌溉定额，根据灌溉制度分解为每月的地下水开采量，核算出模拟期内每个月的农灌开采量用于计算。将收集到的工业和水源地开采量进行处理后，分配到逐月用于模型计算。

根据模型最终调试结果分析，识别后的模型概化和参数符合实际的水文地质条件。模拟的地下水位变化过程与趋势同实际的地下水动态过程基本一致，模拟的地下水流场基本可以反映区域内地下水开采形成的降落漏斗(图2、图3)。

图2 观测孔BH25水位识别期过程线拟合情况

图3 模型模拟的识别期末二层地下水流场拟合情况

2.3 模型的预测

以数值模型识别阶段取得的末流场为初始流场，进行一个水文年的地下水流场的预测。预测过程中，边界条件不发生变化；在源汇项中，选择保证率为50%的年份的降雨量作为降雨入渗和农业开采的计算依据，工业开采量保持不变，水源地开采量除备用井外，选择2014年的开采量数据作为水井开采量，目前尚未完全使用的华美热电1-3号井、西水厂1号井、象山井、桃园街小学井等，类比附近开采井的开采量设定。取得的模拟结果作为模拟水源井迹线的依据。

3 水源地保护区划分

3.1 划分技术方法

对三维稳定流，modpath的质量平衡方程可表示为：

(3)

式中：Vx，Vy，Vz—位线形流动流速矢量在各坐标轴方向的分量，N为含水层有效孔隙率，W—由含水层单位体积源、汇产生的水量，由GMS依据公式(1)(2)计算得出。对于三维非稳定流,可以将其视为由一系列的稳定流时间步长组成，在每一个时间步长内按稳定流方法计算质点运移矢量，通过矢量累加得到非稳定流一定时间段内的质点运移轨迹[12]。

以modflow中各个水井所在网格的中心点作为水质点的起始位置。利用modpath的反向示踪功能按100d，1000d确定水质点位置。连接质点位置构成地下水流截获区。截获区主要沿地下水流向展布，质点在水源地上游运移迹线较长，呈条带状展布；在水源地中游运移迹线较短，呈扇形展布。

由于肥城城区水源地属于群井水源地，根据技术规范对集中式供水水源地保护区范围的规定，井群内井间距小于或等于保护区半径的2倍时，以外围井的外界多边形为边界，向外径向距离为保护区半径的多边形区域。对modpath确定的地下水流截获区进行修正。并结合水源地保护区附近的地标、地物特点，充分利用具有永久性的明显标志，最终确定各级保护区的界线。

3.2 保护区划分结果

保护区划定结果为：一级保护区面积合计1.53km2，二级保护区面积合计28.4km2，准保护区面积约272.9km2。以二级保护区为例，计算得到的1000d地下水流截获区面积为3.36km2。考虑到井群间距，向外径向出保护区半径，并结合附近地标地物，最终确定以北部康汇河河道和南部的道路为保护区边界，圈定保护区面积28.4km2。

图4 1000d截获区范围及划定的二级保护区示意图

3.3 分区可靠性评价

与前人成果相比[16]，该次研究得到的保护区类型从6个减少为3个，面积变化明显，一级保护区面积减少为1.53km2，二级保护区面积从99.2km2减少为28.43km2，准保护区面积从44.6km2扩大为272.9km2。一、二级保护区的范围主要位于建城区，通过合理计算得到较小的保护区面积，有利于地方政府和管理部门安排开展有关建设和保护工作，减少城市发展和水源地保护的冲突。

4 结论

(1)该次研究重点针对肥城市城区水源地开展数值模拟工作，为了弥补资料欠缺，参考了遥感影像、灌溉定额等数据，对缘汇项及水文地质参数进行厘定，对水文地质概念模型进行了较为详细准确的刻画，具有一定的先进性。

(2)数值模型建立的关键在于岩溶含水介质的合理概化以及资料的详细程度。该次研究较为准确的概化了城区水源地的水文地质情况，但在水文地质参数、开采量、水位观测数据等资料获取方面存在欠缺，影响了模型评价的精度。另外现有的等效多孔介质模型，在岩溶地区使用时忽略了岩体的非均质性，本身极易造成较大的计算误差。以上因素降低了分区的可靠性。未来岩溶地区数值模型应用中引入不确定性分析，可在一定程度上提高评价结果的可靠性。

(3)建议在保护区范围内建立完善的地下水监测网络，结合当前开展的国家地下水监测工程和省级地下水监测工程等项目，对水源地水位及水质实现动态观测。按照《饮用水水源保护区污染防治管理规定》相关要求，结合肥城市城区水源地环境现状，针对不同级别的保护区，采取不同的防范措施和管理措施。大力控制地下水开采量，多方寻找替代水源，提高工业的单位用水效率，提倡居民生活节约用水，避免因过度开采引起地下水动力条件发生进一步的变化，对水源地水质造成污染。

[1] 郝永艳.三给地垒岩溶地下水系统及水源地保护区划分研究[D].太原：太原理工大学,2011.

[2] HJ/T338-2007.饮用水水源地保护区划分技术规范[S].

[3] 刘志涛,周群道,杨建华.地下水溶质运移数值法和解析法预测结果对比分析——以沾化电厂为例[J].山东国土资源,2016,32(7):78-82.

[4] 赵红梅,肖杰.公式法与数值模拟法在地下水饮用水源保护区划分中的应用——成都平原某水源地为例[J].四川环境,2013(S1):60-64.

[5] 邓媛媛,胡立堂,高童,等.吴忠市金积地下水饮用水源地保护区划分[J].南水北调与水利科技,2013(01):127-131.

[6] 汤新梅.北京市水源八厂水源地保护区划分研究[D].长春：吉林大学,2012.

[7] 付延玲.邯邢水文地质南单元岩溶地下水系统数值预报[J].中国岩溶,2002,21(4):269-275.

[8] 钱家忠,吴剑锋,董洪信.徐州市张集水源地裂隙岩溶水三维等参有限元数值模拟[J].水力学报,2003(3):37-41.

[9] 刘猛,束龙仓,刘波.地下水数值模拟中的参数随机模拟[J].水利水电科技进展,2005,25(6):25-27.

[10] 卜华,陈占成,张良鹏.山东羊庄岩溶水系统饮用水水源地保护区划分探讨[J].科技创新导报,2008(20):110-112.

[11] 翟立娟.岩溶水饮用水水源保护区划分技术方法——以邯郸市羊角铺水源地为例[J].中国岩溶,2011(01):47-52.

[12] 李星宇,南天,王新娟,等.数值模拟方法在隐伏岩溶水源地保护区划分及污染治理中的应用[J].中国岩溶,2014(03):280-287.

[13] 梁永平,王维泰.中国北方岩溶水系统划分与系统特征[J].地球学报,2010,31(06):860-868.

[14] 魏晓燕.肥城盆地地下水动态数值模拟研究[D].济南：济南大学,2014.

[15] 王玮,李云峰,侯东辉,等.肥城盆地石横电厂水源地岩溶水系统地下水管理模型[J].地球科学与环境学报,2004,26(3):32-39.

[16] 王万喜,王启田,刘福臣,等.水源地保护区划与防护措施——以肥城盆地为例[J].地下水,2008,30(04):42-44.

[17] 崔素芳,张保祥,范明元,等.肥城盆地地下水水化学演变规律研究[J].人民黄河,2015,37(3):75-79.

Partition of Groundwater Source Protection Zones of Feicheng City Based on Numerical Model

WANG Xiaowei1，WANG Yanling2，ZHAO Zhiwei2，CHEN Weiqing2

(1. School of Water Resource and Environment, China University of Geosciences (Beijing), Beijing 100083，China；2. No.5 Exploration Institute of Geology and Mineral Resources, Shandong Tai'an 271000, China)

Karst groundwater had been developed as the drinking water source in many cities in the areas of covered karst areasof northern China. There remains lots of work to do in demarcating karst water drinking source protection zones with the technique of numerical model. This research conducted the partition of karst groundwater drinking source protection zones of Feicheng City in Shandong Province by building a groundwater numerical model and executing particle tracing. The square of the first level protection zone for the drinking water source is 1.53km2, the second level is 28.4km2, and the square beyond protection zone is about 272.9km2. It is showed that protection zones divided by numerical model are reliable and rational compared with the former division. It is suitable for management and protection of groundwater resources by local government.

Karst water; protection of drinking water source; numerical model; particle tracing; Feicheng city

2016-11-21；

2017-01-20；编辑：曹丽丽 作者简介：王晓玮(1982—)，男，山东泰安人，博士研究生，主要研究方向为地下水资源评价和管理；E-mail：wangxw@cugb.edu.cn

P641.8

B

王晓玮，王延岭，赵志伟，等.基于数值模型的肥城市地下水源地保护区划分[J].山东国土资源，2017,33(4)：29-33.WANG Xiaowei，WANG Yanling，ZHAO Zhiwei，etc. Partition of Groundwater Source Protection Zones of Feicheng City Based on Numerical Model[J].Shandong Land and Resources, 2017,33(4)：29-33.