(吉林大学 新能源与环境学院，吉林 长春 130021)

饮用水水质直接关系到人们的生命健康，习近平总书记“4·13”重要讲话中提到：全面开展集中式饮用水源地专项整治，加强近岸海域污染治理，尽快实现城镇污水管网全覆盖、全收集、全处理[1]。地下水水源地保护区划分的主要目的是为了避免水源地附近的人类活动及其所排放污染物对水源造成污染而划定的保护区域，因此在理想条件下，应在人类活动区域与地下水水源地之间预留足够的空间，在人为和自然条件下使地下水得到适当处理，从而达到人们生活和饮用的标准。此外，在水源地出现意外污染时，也可作为处理和修复的区域。建立水源保护区，是为了更好地保护水源地水质安全和应对突发情况做准备[2]。

水源地保护区划分的常用方法有公式法、简化图形法、分析法、解析模型法和数值法等[3-6]，其中数值法通过构建地下水模型，并计入溶质运移模型，以达到划分保护区的目的。已有的地下水保护区划分数值法方面的研究[7-10]大多只建立了水流模型，然后利用粒子示踪来划分保护区范围，比如王涵等[10]利用MODFLOW建立地下水渗流场，利用MODPATH进行粒子反向示踪模拟，对呼和浩特市城市的水源地进行了保护区划分；韩京龙等[11]使用同样的方法对磐石市地下水水源地进行了保护区划分。考虑到污染物在随地下水进行迁移的过程中，由于自身分子扩散的作用，污染物的运移速度要快于地下水流的运移速度，如果仅进行地下水流模拟，并不能真正模拟出污染物的运移范围。因此，以山东省莒县城区相近的三个水源地为研究对象，根据研究区的水文地质条件建立地下水的水流模型，在水流模型的基础上进行溶质运移模拟[11-18]，得到三个水源地的保护区范围，为莒县城区的地下水水源地保护和开发利用提供了可靠的依据。

1 研究区概况

研究区位于山东省东南部日照市莒县，南北长16 km，东西宽10 km，面积约160 km2，如图1所示。地形较为简单，东、西部为丘陵，中部为冲积平原区，地形总趋势为北部高、南部低。平原地形标高95～110 m，地形坡度2～14°。地下水主要以松散岩类孔隙裂隙水为主，由第四系砂及砂砾石组成，厚度20～30 m，含水量大，单井出水量最低可达1 000 m3/d以上，浅层地下水的埋深平水期为3～5 m，丰水期更浅，为2～4 m。包气带岩性以粉质粘土为主，并有粉土、粉砂夹层。地下水主要为降水入渗补给和河流渗漏补给，农业开采和工业开采为主要排泄项，以地下水作为饮用水源，地下水水源地主要有3个，由北向南分别为净水厂、第一水厂和文心水厂。过净水厂由南向北的剖面图如图2所示。三个水厂水井井径约8 m，井深60 m，动态水位范围为80～100 m。

2 水文地质概念模型

根据研究区基础水文地质条件、水动力场、水化学场分析确定概念模型。整个莒县沭河盆地的上层除河谷地区外均为弱透水性的黏土层，中间层为第四系孔隙潜水砂层，底部为隔水的厚层砂岩。

2.1 确定模型模拟区域

确定模拟区域要把可能影响到的区域尽量包括进来，故模拟区东到第四系与基岩的分界处，西到西部沭河和柳清河之间的分水岭处，北到具有相对稳定地下水位的定水头边界处，南到沭河及其一条支流的定水头边界处。

2.2 地层概化

根据研究区的水文地质概况，整个莒县沭河盆地的上层除河谷地区外为弱透水性的黏土层，中间层为第四系孔隙潜水砂层，底部为隔水的厚层砂岩。因此，把地层概化为3层：第一层为透水性弱的第四系孔隙潜水层，第二层为透水性强的第四系孔隙潜水层，第三层为隔水层。

3 地下水水流数值模拟

根据莒县的水文地质条件对研究区进行模型概化，利用GMS软件中的MODFLOW模块建立地下水水流的数值模型。

图1 研究区地理位置图Fig. 1 Geographical location of the study area

图2 地层剖面图Fig. 2 Stratigraphic profile

3.1 地下水水流数学模型

根据地质条件和水文地质条件分析，莒县城区地下水流可视作平面二维流，对应的地下水流数学模型为：

(1)

其中：h—地下水水头，m；Kx、Ky—x、y方向渗透系数，m/d；h1—含水层第一类边界上的水头，m；q—含水层第二类边界上的流量，m3/d；h0—含水层初始水头，m；W—源汇项强度(包括开采强度等)，m3/d；Σ1—含水层第一类边界；Σ2—含水层第二类边界；D—研究区；μ—给水系数。

3.2 边界条件处理

研究区存在一类水头边界和二类流量边界。南部为一条河流，常年有水，与地下水存在水力联系，故处理为水头边界。北部有多个地下水位观测孔，根据长期水位观测资料显示，在沭河流域进行地下水开采对该处地下水位无明显影响，因此可处理为水头边界。在研究区的西部有一地表分水岭，位于沭河和柳清河之间，两侧的地下水水位高于河流两侧水位，所以该处可作为研究区的地下隔水边界。研究区东部的定流量边界是第四系和基岩的交界线，由于基岩区地势高且透水性较差，大气降水沿基岩表面经地表径流汇入研究区内的含水层内，故定义为给定流量边界。

3.3 模型区域剖分

根据研究区的水文地质条件，研究区在垂向上分为3层，水平方向分为128×187的网格，每个网格代表100×100的区域。

3.4 源汇项的处理

研究区内潜水含水层的补给方式主要包括降水入渗补给、农业水灌溉补给和河流侧向径流补给，排泄方式主要为蒸发排泄和人工开采。人工开采按不同开采分区的开采总量平均分布到每个开采分区上，开采分区图如图3所示，每个分区的开采量如表1所示；根据研究区降雨情况，降水入渗系数取0.125～0.3，得到不同时间的降雨补给强度，见图4；灌溉回渗强度为0.000 022 m/d；多年平均蒸发为1 410 mm，蒸发系数为0.04，蒸发强度为0.000 016 m/d，蒸发极限埋深为2.5 m。

图3 研究区开采分区图Fig. 3 Zoning map of mining in the study area

图4 不同时间段下的降雨补给强度Fig. 4 Rainfall recharge intensity at different time periods

表1 不同开采分区的开采强度

Tab. 1 Mining intensity of different mining zones m/d

分区一区二区三区四区开采量0.000 024 70.000 140 00.000 411 00.000 685 0

3.5 水文地质参数

主要为渗透系数K、潜水的给水度、承压水的贮水系数和降雨入渗系数α等，各水文地质参数按相应试验、地层岩性和参考经验值等途径综合给定初始值，并通过模型调参，进而得到研究区模型的水文地质参数。

3.6 模型的识别验证

在上述条件下运行模型，可得到研究区的模拟地下水流场图。通过与相同时期的真实地下水流场图进行对比拟合，对相关水文地质参数进行合理化调整，使模型与实际情况更加契合。该模型进行拟合的主要原则为：在水文地质参数满足研究区水文地质条件的基础上使模型的等水位线和水均衡中各源汇项与实际相符。取两个长期观测点对2016年的观测水位进行模型识别，取2017年丰水期的地下水等水位线进行验证。经过识别和验证，得到研究区水文地质参数值如表2所示。其中第一层的渗透系数(岩性)分区如图5所示，第二层和第三层渗透系数未分区，第二层、第三层均为砂岩，每个分区的渗透系数见表2。

取2017年丰水期的观测水位进行模型验证。图6为模型计算的2017年丰水期流场与实测流场的对比，模型模拟的水位和模拟的流场与实际水位和实际流场也比较接近。经过模型验证，表明本次模拟具有较高的精度和可信性。

表2 研究区地层水文地质参数

单位：k—m/d；μ—m3/d

图5 第一层含水层渗透系数分区图Fig. 5 Zoning map of permeability coefficient of the first aquifer

图6 模型验证期实测水位与 计算流场的对比图Fig. 6 Comparison of measured water level and calculated flow field during model validation period

图7 保护区范围图Fig. 7 Scope map of the protected areas

4 溶质运移模拟

在上述地下水水流模型的基础上，应用GMS软件中的MT3DMS模块进行地下水的溶质运移模拟，用以计算污染物的捕获区范围为保护区范围。由于NaCl性质稳定，一般不参与物理化学反应，本身对环境影响小，在保护区划分过程中更常用。依据《饮用水水源保护区划分技术规范》(HJ/T338—2018)[4],模型以地下水开采井为中心，NaCl做为溶质，以NaCl迁移100 d的范围做为水源地一级保护区范围；在一级保护区外，以NaCl迁移1 000 d的范围做为二级保护区范围。

4.1 地下水溶质运移数学模型

地下水溶质运移的对流-弥散方程：

(2)

其中：Rd—阻滞因子；c—地下水因子浓度，mg/L；t—时间，d；xi—坐标轴方向上的距离，m；Dij—地下水水动力弥散系数；vi—地下水渗流速度，m/d；qs—源和汇的单位流量，m3/d；cs—源和汇的浓度，mg/L；θ—含水层孔隙率，%。

4.2 模拟过程

保护区范围应以开采井为中心，地下水在开采过程中，周围污染物向开采井中心移动，以100 d和1 000 d污染物能运移到开采井中心的位置点组成的范围为保护区边界，但在实际应用过程中不可能对每个点都进行模拟预测来寻找保护区边界，因此在应用过程中转换思路，将水源地的开采井变换为注水井，注水量与开采量相同，给开采井设置一个初始污染物浓度，模拟污染物从开采井向四周运移的过程，以污染物从开采井向四周运移100 d和1 000 d形成的污染晕为保护区边界，本次模拟过程中注水量设置为10 000 m3/d，污染源浓度设置为10 000 mg/L,根据场地的水文地质情况，参照吉林省双阳区弥散试验场地所得出的结果，对比地层岩性、水流速度等参数，取此研究区弥散系数为50 m，模拟总时长为1 440 d(4年)。本次研究中的三个水源地均为孔隙潜水型水源地，主要取水层位均为第二层的砂层，因此主要关注污染物在模型第二层中的运移范围。最终模拟的溶质运移100 d和1 000 d的污染晕范围如图7和表3所示。

表3 水源地溶质运移各方向位移表

由模拟结果可以看出：随着时间的增加，污染物的污染晕逐渐扩大，而且沿水流方向污染晕的迁移速度更快，说明污染物的运移扩散主要受地下水流动的影响，同时在其他方向污染物也发生了迁移，说明污染物的运移扩散除了受水流作用的影响外还受污染物自身分子扩散作用的影响。

在保护区实际管理操作中，考虑到三个水源地相距较近，参照计算结果，可以将独立的三个水源地二级保护区适当外扩边界，整体划入一个沿河谷的条带区域内(可作为准保护区)进行保护管理。

5 结论

以莒县城区的三个水源地为研究对象，利用GMS中MODFLOW和MT3DMS模块，对水源地进行保护区划分，得出结论如下：

1) 首先应用GMS软件的MODFLOW模块建立了莒县城区的地下水水流数值模型，通过等水位线的拟合情况验证了模型的合理性，说明建立的模型能够准确刻画保护区的范围。

2) 在水流模型的基础上再通过MT3DMS模块建立了研究区的溶质运移模型，得到溶质运移100 d和1 000 d的污染晕范围：依据《饮用水水源保护区划分技术规范》(HJ/T338—2018)，将溶质运移100 d和1 000 d的范围划分为一级保护区和二级保护区，净水厂、第一水厂和文心水厂的一级保护区和二级保护区的污染晕范围面积分别为34 617和150 800 m2、54 467和203 000 m2、32 200和164 000 m2。

3) 在实施保护区划分时，考虑到三个水源地相距较近，参照计算结果，实际操作时，可以将三个独立的二级保护区适当外扩将三个水源地整体划入一个沿河谷的条带区域内进行保护管理。