地下水水源地保护区划分方法研究
——以大汶口盆地为例

王 昕 翌，王 敏，冯 建 国，高 宗 军

(山东科技大学 地球科学与工程学院，山东 青岛 266510)

地下水是我国重要的供水资源，地下水水质的优劣直接关系到人类的饮水安全和健康状况。目前，我国居民饮水安全状况不容乐观，饮用地下水水源污染已成为一个较突出的社会问题，地下水保护的一个直接有效的途径就是建立地下水水源地保护区[1-4]。国内外已有许多关于地下水水源地保护区划分方面的研究与讨论，如田文英[5]通过解析模型法对辽宁省柳河平原区水源地进行了保护区划分；王丽娟等[6]选用稳定流和非稳定流模型法进行了保护区划分。但是，地下水水源地的复杂性导致运用单一的方法进行区划划分会存在划分结果与实际不符、划分精度不高等问题。因此，学者们对区划方法做出了许多改进：如赵红梅等[7]选用公式法和数值模拟法2种方法对成都平原某水源地进行保护区划分，并比较2种方法的区划结果，认为在实际应用中2种方法相结合所得结果更为合理；江广长等[8]对多种保护区划分方法进行对比分析，认为各类划分方法的不足之处是未能综合考虑水源地内各项影响因素，区划方法与地下水水源地间的联系不甚紧密，并提出在数值模拟法中引入不确定分析的观点。

空间叠加技术在水文地质研究中的广泛应用，也为改进区划方法提供了思路。如Dawoud[9]通过空间叠加手段建立了覆盖埃及全国的地下水质量监测网，将多组包含了不同含水层类型和水化学特征的空间数据图层叠加，实现了如人口密度、水文地质条件、盐渍化风险等众多影响因素的信息整合。将空间叠加技术应用到区划工作中，有助于提高划分精度，使区划结果能够全面地反映出水源地的特征。

结合前人的研究，本文以大汶口盆地为例，首先选用数值模拟法对盆地内水源地进行保护区划分，在此基础上，应用空间叠加分析，将数值模拟法与地下水污染风险评价、防污性能评价结合，修正数值模拟法的区划结果，并将2种方法所得结果进行对比讨论，以期完善保护区划分方法，为地下水安全保护提供依据。

1 大汶口盆地概况

大汶口盆地位于山东省西南部，地跨泰安、新泰、肥城等县市，城镇众多，交通便利，农业发达，矿产丰富。盆地内地表水系发育，区域性河流为汶河水系，位于盆地中南部，主要支流有大汶河和漕河。整体地势东高西低，南高北低，总体形态为单断凹陷，是一个较完整的水文地质单元。

盆地内地下水类型有岩浆岩风化裂隙水、碎屑岩类孔隙裂隙水、松散岩类孔隙水和碳酸盐岩裂隙岩溶水4种类型。地下水主要依靠大气降水补给，汶河的侧向补给和渗漏对大汶口至西界一带的孔隙水及岩溶水起着重要作用。由于盆地受构造控制，各时代沉积物差异很大，地下水的开发利用受地层及富水性控制明显，地下水开采以岩溶水为主，其次为孔隙水和裂隙水。

盆地内有东武、楼德2处地下水水源地，是当地居民的主要供水水源地。2处水源地均以开采地下岩溶水为主，水位动态主要受大气降水、汶河径流量以及地下水开采量的影响，排水途径均为补给上覆第四系孔隙水，侧向径流排泄及人工开采，二者间的水力联系十分密切。其中东武地下水水源地分布在东武驾庄-大汶口镇一带，开采井深度为100～200 m，单井涌水量一般在5 000 m3/d以上，地下水开采量2.2万m3/d；楼德地下水水源地分布在楼德-甘露村-西杜村一带，开采井深度为130～300 m，单井涌水量一般在3 000 m3/d 以上，地下水开采量为2.0万m3/d。自1989年以来，东武、楼德2处水源地长期处于开采状态，2006年12月至2018年12月，东武地下水水源地平均水位下降了41.58 m，楼德地下水水源地平均水位下降了38.97 m，岩溶水水位呈逐渐下降趋势。同时工业废水、农药化肥等污染源不断增多，渗入地下导致地下水污染，地下水水质急剧恶化。根据1996年水质监测资料显示，东武、楼德水源地水化学类型均以HCO3—Ca型水为主，经2018年调查，东武水源地水化学类型以HCO3·SO4·Cl—Ca型水为主，可溶性固体平均为1 291 mg/L，总硬度平均为889 mg/L，SO42-浓度平均为314 mg/L；楼德地下水水源地水化学类型为HCO3·SO4—Ca型水，总硬度平均为330 mg/L，SO42-浓度平均为30 mg/L，HCO3-浓度为230 mg/L。与开采之初相比，2处水源地地下水各项指标均有所增加，SO42-浓度明显升高，水质状况以Ⅲ类较差水和Ⅳ类极差水为主，地下水环境保护刻不容缓。

2 地下水水源地保护区划分方法

2.1 划分方法与思路

地下水水源地并非独立存在的，开展地下水水源地保护区划分工作，应在全面了解整个水文地质单元特征的前提下，以地质单元内现有水源地为核心进行区划，将保护区的范围涉及整个水文地质单元。传统的数值模拟法虽然能够较客观地刻画含水层结构及水文地质条件，但依然存在参数的不确定性、模型概化过程中易出现误差等问题[8]。因此，本文首先应用传统的数值模拟法，在建立研究区地下水水流模型和水质模型的基础上，依据各级保护区对污染溶质运移时间的不同要求，来初步划定不同级别保护区的边界范围。在此基础上，将地下水防污性能评价和污染风险评价引入区划方法中，把研究区的地质条件、人类活动、污染源分布情况等因素纳入区划标准中，综合评估整个水文地质单元的地下水防污能力和污染风险，然后利用空间叠加技术，整合数值模拟法和2项评价所得的空间数据图层，同时结合研究区自然地理和水文地质条件对保护区边界进行调整，得出最终的地下水水源地保护区范围(见图1)。

图1 保护区划分方法流程示意Fig.1 Structure diagram of protection zone division method

2.2 数值模拟法—保护区初步划分

应用数值模拟法划分地下水水源地保护区范围的原理，是根据研究区的含水层结构、地质构造等基础数据，建立研究区的水文地质概念模型，在此基础上，利用模拟软件如GMS(Groundwater Modeling System)、Visual Modflow等来模拟研究区的地下水流场。根据研究区内污染源的空间分布情况和水质资料，进行地下水溶质运移模拟，对地下水水源地周边的污染物入渗点位置进行预测，在水源地周边不同位置设置入渗点，其中那些污染溶质运移100 d和1 000 d后刚好到达水源地的入渗点，分别为水源地一级和二级保护区的边界点[10-11]。

2.3 地下水防污性能评价

目前国内外一般采用DRASTIC法进行地下水防污性能评价。DRASTIC法通常选取地下水埋深D、地下水的净补给量R、含水层介质A、土壤类别S、地形坡度T、包气带介质I和含水层水力传导系数C这7项影响因子，通过对指标进行量化评分来代表各项因子的数值范围或类别(见表1)，并给予各项因子不同的权重参数[12]，DRASTIC地下水防污性能DI值由公式

表1 评价因子评分体系Tab.1 Evaluation factor scoring system

(1)确定：

DRASTIC=5D+4R+3A+2S+1T+5I+3C

(1)

7项因子的权重赋值可根据研究区实际情况进行调整，使评价结果更为合理。在一般条件下，地下水防污性能指数为23～226，根据计算后所得DI值，将地下水防污性能指数DI共分为5个等级(见表2)，DI值越高，则地下水防污能力越差。

表2 地下水防污性能分级标准Tab.2 Classification standards of groundwater anti-fouling performance

2.4 地下水污染风险评价

地下水污染风险评价是在污染载荷风险(P)评价、含水层固有脆弱性(U)评价、污染危害性(V)评价的基础上进行的。污染源荷载风险(P)的评价指标为：研究区内污染源种类K、污染物产生量Q、污染物释放可能性L和调查点与污染源间的距离D，并依据研究区实际情况对4项评价指标分别打分[13-14]，按公式(2)计算单个污染载荷P′的值。

P′=KQLD

(2)

将单个污染载荷P′的计算结果按照表3的分级标准重新分类，即为污染载荷风险(P)评价的结果。

表3 污染载荷风险评价结果重新分级标准Tab.3 Pollution load risk assessment results reclassification criteria

含水层脆弱性的高低与地下水防污能力的强弱密不可分，地下水防污能力强，则表示脆弱性低；地下水防污能力弱，即脆弱性高[15]。故含水层固有脆弱性评价(U)将地下水防污性能评价的结果进行重新分级并赋值，赋值标准见表4。

污染危害性(V)一般根据研究区内地下水的用途来进行分级，因此本文据此将污染危害性分为3级(见表5)。最终按照表6的评价标准得到污染风险评价结果R，其中“0”代表低污染风险；“1”表示中等风险；“2”表示高污染风险[14]。

表4 含水层固有脆弱性重新分级标准Tab.4 Reclassification criteria for the inherent vulnerability of aquifers

表5 污染危害性分级标准Tab.5 Pollution hazard classification standards

表6 地下水污染风险评价标准Tab.6 Risk assessment criteria for groundwater pollution

2.5 空间叠加分析

通过对不同图层的叠加分析，可以整合多层数据间的联系与变化特征[9，16]。应用GIS软件如MapGIS、ArcGIS平台上的空间分析功能，首先进行空间叠加求交运算，分离出地下水防污能力低且污染风险高的区域和地下水防污能力低且污染风险中等、地下水防污能力中等且污染风险高的区域；再通过空间叠加求并运算，将上述区域分别纳入数值模拟法拟定的一、二级保护区范围内。空间叠加分析分级标准列于表7。

3 大汶口盆地地下水水源地保护区划分

本文将大汶口盆地整个水文地质单元作为研究区，以盆地内东武、楼德2处地下水水源地为核心，根据2018年调查获得的45个承压水水质监测点及盆地基础数据，利用数值模拟法和基于数值模拟法的空间叠加方法对大汶口盆地进行保护区划分，并对2种方法所得结果进行对比讨论。

表7 空间叠加分析分级标准Tab.7 Spatial superposition analysis division criteria

3.1 大汶口盆地数值模拟结果及保护区初步划分

根据大汶口盆地水文地质条件，盆地内地下水开采以岩溶水为主，且岩溶含水层与上覆孔隙含水层间有稳定分布的巨厚第三系泥页岩阻隔，故本次模拟以岩溶含水层为目的含水层，将其概化为非均质各向同性的承压含水层，水流运动为平面二维流。模拟区包括一类水头边界和二类流量边界，其中盆地北侧为多组断裂视为隔水边界，大汶河流域处理为水头边界。模拟时仅考虑对流、弥散作用对于溶质浓度变化的影响，依据盆地经验数据，孔隙度取值为0.07，弥散系数取值为15 m2/d。

数值模拟模型包括水流模型和水质模型2个部分，模拟选用GMS6.0软件中的MT3D模块，在运用ModFlow模块模拟出地下水流场的基础上运用MT3D模块进行地下水溶质运移模拟。模拟区内地下水流特征用如下数学模型表示：

(3)

式中：K为岩溶含水层渗透系数，m/d；W为含水层的源汇项，m2/d；Ss为承压含水层的贮水率；τ为流量边界；h0为初始水位。

模拟区地下水溶质运移用如下数学模型表示：

(4)

结合大汶口盆地内污染源的空间分布情况，分别在东武和楼德水源地周围不同位置设置污染物入渗点进行试算，其中那些污染质运移100 d和1 000 d后刚好分别到达东武、楼德地下水水源地的入渗点，分别为其一级和二级保护区的边界点(见图2)。将这些边界点圈闭，即为水源地一级、二级保护区的理论范围(见图3)。其中，东武水源地一级保护区面积为17.2 km2，楼德水源地一级保护区面积为15.3 km2，模拟结果显示2处水源地的二级保护区范围互相交叠，故将二者的二级保护区范围划定为一个整体，二级保护区面积为456.3 km2。

图2 溶质运移模拟Fig.2 Solute transport simulation diagram

图3 数值模拟法初步划定保护区范围Fig.3 The scope of protected areas is preliminarily determined by numerical simulation method

3.2 大汶口盆地地下水防污性能评价

由于大汶口盆地内农田大量喷洒农药，所以为贴合盆地实际情况，对地下水防污性能评价的权重稍作调整，适当增加土壤类别S、地形坡度T和包气带介质I这3项影响因子的权重，调整后的地下水防污性能DI值由公式(5)确定：

DRASTIC=5D+4R+3A+5S+3T+4I+3C

(5)

根据盆地内各岩溶水点的调查数据，结合表1的评分体系和表2的分级标准，将大汶口盆地地下水防污性能分为四级(见表8)。防污性能高和非常高的地区呈带状分布在盆地北部、东北部和南部山区，富水性小，第三系盖层较厚，地下水防污能力强；防污性能低区主要分布在大汶口镇、楼德镇周边地区以及大汶河富水性大的地段，河流对地下水的补给量大，局部地区第三系盖层缺失，地下水易受地表水的影响；防污性能中等的地区主要分布在防污性能低区的外围。大汶口盆地地下水防污性能评价空间分布如图4所示。

表8 地下水防污性能评价综合成果Tab.8 Comprehensive evaluation results of groundwater anti-fouling performance

图4 大汶口盆地地下水防污性能评价空间分布Fig.4 Spatial distribution of groundwater anti-fouling performance evaluation

3.3 大汶口盆地地下水污染风险评价

综合大汶口盆地内各调查点污染载荷风险评价(P)、含水层固有脆弱性评价(U)、污染危害性评价(V)的最终评分，对应表6的评价标准，将大汶口盆地地下水污染风险分为三级(见表9)。污染风险高的地区主要分布在汶河北岸、大汶口镇和楼德镇周边地区，污染渠道较多且人口众多，地下水开采量大且多用于饮用，另外汶河南岸第三系缺失地区，地下水受污染河流直接补给，地下水受污染可能性极大；污染风险中等的地区主要分布在高污染风险地区的外围；污染风险低的地区主要分布在盆地北部、东部和南部的基岩山区，人类活动稀少，对地下水的开采利用较少，并且远离污染源。大汶口盆地污染风险评价空间分布如图5所示。

表9 地下水污染风险评价综合成果Tab.9 Comprehensive results of groundwater pollution risk assessment

图5 大汶口盆地地下水污染风险评价空间分布Fig.3 Spatial distribution of groundwater pollution risk assessment

3.4 空间叠加及保护区范围确定

在MapGIS软件中的空间分析平台上，将地下水防污性能评价结果(见图4)、地下水污染风险评价结果(见图5)和数值模拟法划定的保护区范围(见图3)进行空间叠加(见图6)。依据表7的叠加分级标准，把大汶口盆地内地下水防污能力低且污染风险高的区域加入一级保护区范围内；把地下水防污能力中等且污染风险高、地下水防污能力低且污染风险中等的区域加入二级保护区范围内。

图6 空间叠加过程示意Fig.6 Schematic diagram of spatial superposition process

根据叠加整合结果，虽然盆地西部汶河下游区域，地下水防污能力弱且污染风险高，但距离东武、楼德水源地极远，是地下水的排泄区，故不纳入地下水水源地一、二级保护区范围内。最后参照盆地的自然环境和水文地质条件对保护区边界进行适当调整：尽量利用盆地内的交通道路、隔水断层等界线作为保护区边界，结合盆地的地形地貌对保护区边界进行适当增加或缩减，准保护区范围扩展至整个盆地，即以水文地质单元的地表分水岭为边界。叠加后东武水源地一级保护区面积为33.6 km2，位于盆地中部；楼德水源地一级保护区面积为30.4 km2，位于盆地东南部；二级保护区面积为569.5 km2(见图7)。

图7 大汶口盆地地下水水源地保护区Fig.7 Groundwater source protection area of Dawenkou basin

与数值模拟法的划分结果相比，东武水源地一级保护区边界向东北方向有延伸，增加面积为13.2 km2；楼德水源地一级保护区边界向东南方向有延伸，增加面积为15.1 km2；二级保护区向西侧、南侧以及东南侧延伸，面积增加113.2 km2，一级、二级保护区边界均向地下水补给源方向扩展。根据结果显示，经过空间叠加后得到的大汶口盆地地下水水源地保护区，包括了盆地内地下水补给区中污染风险高和防污能力弱的区域，得到的结果更为贴近大汶口盆地实际情况。

2种方法的划分结果存在差异，主要原因是传统的数值模拟法虽然考虑了地下水径流方向等因素的影响，但由于数值模拟过程中存在的种种误差，使得划分出的保护区与实际情况存在出入。而在对盆地进行地下水污染风险和防污性能评价时，能够充分利用渗透系数、包气带厚度等地质点量化信息，同时亦将盆地内污染源的空间分布及地下水用途等因素纳入了区划标准，客观反映了盆地的水文地质条件和地下水环境状况，最终通过空间叠加手段，将综合了多项指标的地下水污染风险和防污性能评价结果整合入数值模拟法中，优化数值模拟法的区划结果，弥补了数值模拟法的不足。

基于数值模拟的空间叠加法的优点是思路清晰，适用性强，可以根据研究区的实际情况来选择评价方法和需要叠加整合的要素，结果贴合实际；本次针对实际水源地的区划过程也反映出该方法工作量较大且需要足够多的基础资料支撑，需要以研究区准确的水文地质信息及充足的调查资料为依托。本次区划由于盆地资料限制，仅考虑了地下水受污染的可能性和防污性能等要素，而如地下水水质、水量的动态变化等因素未能综合处理，在后期工作中，可以将更多可能对区划结果造成影响的因素纳入考量，优化叠加分级标准，进一步完善数字化的综合评价体系。该方法适用于水文地质背景复杂的研究区，尤其对于以开采岩溶水为主的大型水源地效果更佳，在实际应用中，可针对不同研究区的具体情况来调整评价手段和打分体系，并与其他区划方法和评价手段相结合，使保护区划分结果更加准确合理。

4 结 论

(1) 在数值模拟法所得保护区范围基础上，通过空间叠加技术整合地下水防污性能和污染风险评价的结果，得到大汶口盆地东武地下水水源地一级保护区面积为33.6 km2，楼德地下水水源地一级保护区面积为30.4 km2，二级保护区面积为569.5 km2，准保护区面积为1 490.0 km2。区划结果更为贴合大汶口盆地实际情况，证明基于数值模拟的空间叠加法具有可行性。

(2) 地下水防污性能和污染风险评价可以较全面地综合研究区的水文地质特征、人类活动及污染源的空间分布等环境要素，而空间叠加技术有效地将2项评价结果与数值模拟法的区划结果结合，修正传统数值模拟法存在的误差，使区划结果能更加全面地反映出水源地的特征，为保护区内污染防治、土地利用、设施建设等做出相关规定提供依据，可更好地实现地下水资源的保护，提高地下水资源的利用效率。