基于MODFLOW 模型滹沱河傍河地下水源地保护区划分

冷文鹏，陶 亚，孙若涵，夏建新

(1. 中央民族大学 生命与环境科学学院，北京 100081；2. 生态环境部环境规划院，北京 100012)

我国饮用水源地面临水质污染超标、突发水污染事件频发等安全隐患，特别是具有水量丰富、水质比较稳定等特点[1]的地下水水源地被污染后不易修复，因此，要加强对地下水饮用水源地的保护。为保护地下水水质，建立地下水源地保护区是必不可少的措施[2]。但是，地下水流动比较复杂，地下水源地保护区划分存在一定的技术困难。利用数值模拟法确定地下水水源中的污染物运移方向和距离，可为保护区划分提供依据。邓媛媛等[3]利用FEFLOW 软件建立了吴忠市金积傍河型地下水水源地的三维地下水数值模型，将水源地的12 眼水井概化成“一眼、四眼、十二眼”，模拟了3 种地下水开采情景下地下水的运动，简化复杂的边界条件，并利用反向示踪粒子跟踪技术，建立了三级保护区。王涵等[4]在内蒙古呼和浩特市城市水源地运用MODFLOW 建立了该水源地的地下水流场，根据148 口抽水井单独设置了148 个反向示踪粒子点位，计算出反向示踪粒子的运移迹线，确定了水源地各级保护区的界线。肖杰等[5-9]运用Visual MODFLOW 中的modflow 和modpath 模块分别对崇州市地下水源地、渭河傍河水源地、关中盆地水源地进行了数值模拟，按照抽水井位置设置单个示踪粒子，对其进行保护区的划分。目前边界条件的相关研究方向集中在单井概化，缺少对井群概化的研究和应用。边界条件处理方法的差异可能导致模拟结果不同，按照单井概化处理的研究结果可能使保护区范围过小，带来污染风险。因此，通过研究两种不同的边界条件，分析其计算结果，对于提高保护区划分结果的科学性和准确性具有重要意义。

本文主要以石家庄市新调整后的滹沱河应急后备地下水源地为例，设定了单井和井群两种边界条件，基于MODFLOW 软件计算分析地下水流场，运用数值模拟法计算不同工况下反向示踪粒子100 和1 000 d的运移距离，并分析入渗场使用情况下的示踪粒子运移规律。

1 研究区概况

滹沱河应急后备地下水源地位于石家庄市新华区，西望太行山，北依滹沱河。按照供水计划该水源地要保持供水规模10×104m3/d，布设开采井17 眼，单井开采量6 000 m3/d，井眼内直径800 mm，利用段为第四系第Ⅱ+Ⅲ含水组，井距800 m，水源地面积约为10 km2。

石家庄市西部地处太行山中段，东部为滹沱河冲积平原。区内面积为1 521.08 km2，地势为西北高东南低，地面坡降(1.6～2.5)‰。滹沱河流域的地下水主要赋存在第四系松散岩类孔隙中，含水层多由亚砂土、砂、卵砾石组成，粒度粗、厚度大，水动力特征为潜水、微承压水。水源地拟取水的含水层为第Ⅱ含水组和第Ⅲ含水组。第Ⅰ含水组埋藏深度为15～20 m，该层沉积较薄，颗粒较细，岩性为粉土、细、中、粗砂及含砾石，由于当地长期取水，该层基本疏干。第Ⅱ含水组底板埋藏深度约为100 m，含水层厚度30～50 m，该层沉积厚度大，含水层颗粒较粗，且磨圆度较好；主要岩性为砂砾、卵砾石，透水性及富水性好。该层分为上、下两段，尤以下段含水层最为丰富，单位涌水量30～40 m3/(h·m)，渗透系数一般为37～145 m/d；地下水质良好，矿化度小于0.5 g /L。第Ⅲ含水组底板埋藏深度约为220 m，自西北向东南倾斜，含水层厚度大于50 m，岩性含砾卵石、砂砾夹砂质黏土，其中砂卵石、砂砾石分选较差，该层在经济技术开发区以西遭受了不同程度的风化，透水性和富水性均较差；开发区以东富水性较好，受地方开采井连通影响，本区水动力特征属潜水-微承压水，矿化度小于0.5 g/L。

2 数值模拟法划分保护区

2.1 软件模拟过程概要

MODFLOW 软件由20 世纪80 年代美国地质调查局（USGS）的Mc Donald 和Harbaugh 研发，基于有限差分法原理用于孔隙介质中地下水流动的数值模拟。本研究采用MODFLOW 6.0.3 计算地下水流场，并通过MODPATH（粒子跟踪后处理软件程序）计算示踪粒子的运移迹线。首先按照特定的格式整理空间展布资料（地表高程、分层数据、边界位置）、含水层参数、初始水位、源汇项（大气降水入渗系数、蒸发排泄系数、人工开采量）及水位动态观测资料。随后，根据水文地质剖面图确定含水层数，根据以上资料输入新建模型中的各含水层底板高程数据、初始水流场数据、含水层参数。其中，人工开采量和侧向径流补给量分别以抽水井、注水井的方式确定[10]，分别按照所需格式输入数据。输入11 个水文地质参数分区数值，使用PEST 模块调整水文地质参数，使得模拟值和实测值的差值在误差允许范围之内，运行MODFLOW，得出符合实际条件的地下水流场数值。最后通过MODPATH 模块设置反向示踪粒子，计算100 和1 000 d 粒子运移的距离。

2.2 地下水流动模型

以滹沱河傍河地下水井群为中心，将研究区边界条件、径排补给、各区渗透系数、给水度等要素进行概化，确定研究区域范围为30 km×45 km。石家庄市区17 口应急后备水井位于滹沱河河心岛处，取水第Ⅱ、Ⅲ含水层即为潜水和微承压水，因此本次模拟的主要含水层为微承压水层。

MODFLOW 共有自然邻点插值法、反距离插值法和Kriging 空间插值法等3 种插值方法。本次模拟采用Kriging 空间插值。刘光孟等[11]在几种插值方法的比较研究中，选取疏密差异较大的均匀和不均匀离散点数据进行试验，得出在少量的不均匀分布数据中，Kriging 空间插值法具有较好的精度。因此，基于本研究区域有限的地表数字高程资料，选择Kriging 空间插值法计算得到模拟区域地表高程模型的效果最优。在垂向上根据含水组特点划分为3 层，通过Kriging 空间插值对数据进行初步处理后，获得区域内第Ⅰ含水组底板高程（38.6～118.5 m）、第Ⅱ含水组底板高程（−36.2～83.7 m）及第Ⅲ含水组底板高程（−96.3～71.6 m）。其中，第Ⅰ含水组层厚为10～49 m、平均层厚19.3 m；第Ⅱ含水组层厚为10～76 m、平均层厚43.4 m；第Ⅲ含水组层厚为10～90 m，平均层厚53.6 m。

在黄壁庄水库西侧（山区与平原衔接处）和东北侧（滹沱河冲洪积扇与磁河冲洪积扇交接处）分别存在一条水文地质分区界线，将这两条分区界线作为西侧和东北侧研究区边界。西侧边界处在山区与平原衔接处，作为弱透水的定流量边界。西北部边界设在黄壁庄水库大坝位置，作为定水位边界。由历年地下水流场图得知：研究区的南部边界和北部边界等水位线形状变化不大，且边界与等水位线基本垂直，因此模型的南部边界和北部边界可作为零流量边界处理。研究区东部边界设在地下分水岭处，以藁城区岗上镇至栾城区城关镇为界，该处水位变化稳定，作为定水位边界。

2.3 数学方程

结合石家庄市区的水文地质条件，将研究区模型概化为非均质非稳定三维流地下水系统。

2.4 计算参数确定

滹沱河地下水源地地层多为砂、砂卵砾石层，结构松散，孔隙度大，一般厚40～50 m，随着地下水位逐年下降，目前水源地含水层厚20～30 m，西厚东薄；区域内富水性及导水性良好，渗透系数K 为300～400 m/d，导水系数T 一般在15 000 m2/d 左右，为强富水含水砂层。由于各含水层的水文地质参数在平面和空间上存在区域差异性，因此本文结合前人水文地质研究成果[12]进行参数率定，经过模型调试和识别，13 个水文地质参数分区参数值如表1 所示。

表1 研究区水文地质参数分区Tab. 1 Division of hydrogeological parameters in the study area

2.5 模型识别及结果

根据研究区的地下水位观测资料，去除观测序列不完整或者波动异常的水位观测数据，选用永安村、省二院、省委党校、西兆通村、凌透村、宋营村6 个水位监测井2016 年1 月至2017 年12 月的观测数据进行模型率定识别，率定期内以30 d 为1 个时间步长，选取2017 年1—12 月作为模型的验证期。本次模拟参数的率定采用试估-校正法，结合实际的水文地质条件和前人的研究成果，不断对模型进行约束性的调试。采用置信区间、均方根、相关系数等进行误差分析，计算值和观测值的拟合效果见图1。

图1 监测井计算值与实测值的拟合Fig. 1 Fitting of calculated value and measured value of monitoring wells

从图1 可见，监测井计算值和观测值的最大拟合误差是凌透村的监测点位0.934 m，最小拟合误差是省委党校的监测点位0.098 m，监测的误差范围为0.098～0.934 m，其中拟合误差绝对值小于0.5 m 的监测井数占总监测井的83%，所有监测点位在满足置信度为95%时的相关系数为0.999。由以上分析可得，模型的整体拟合效果颇为理想，满足校正标准。

3 模拟结果分析

从计算流场的结果看，总体流向从西北向东南，最大流动速度小于1.0 mm/s。在滹沱河傍河地下水井群处，地下水位高程介于45～55 m，如图2所示。

图2 水源地区域点位及地下水流场（单位：m）Fig. 2 Point location and groundwater flow field in water source area (unit: m)

根据《饮用水水源地保护区划分技术规范》[13]（HJ 338—2018）可知，地下饮用水源地一级保护区是以取水井为中心，溶质质点运移100 d 的距离所圈定的范围。二级保护区是在一级保护区外，溶质质点运移1 000 d 的距离所圈定的范围。准保护区是包括水源地的补给区和径流区。在本研究中，溶质质点运移是利用MODPATH 粒子示踪程序来进行计算可视化反向示踪粒子的运移流线。

工况1：水源地的17 口水井单独设置示踪粒子，计算其运移距离。经过计算，示踪粒子在第100 d 时运移距离为0.10～0.83 km，平均长度为0.54 km，其运移流线呈现长方喇叭状；示踪粒子在第1 000 d 时运移距离为3.0～7.4 km，平均长度为6.1 km，其运移流线同样呈现长方喇叭状，如图3（a）所示。

图3 不同工况下反向示踪粒子的运移流线和保护区示意图Fig. 3 Flow line of reverse tracer particles and schematic diagram of protected area under different working conditions

工况2：将水源地的17 口抽水井概化成1 个水井群，将反向示踪粒子设置在井群的边界处（井群边界即应急后备水源地的边界位置），计算其运移距离。经过计算，示踪粒子在第100 d 时运移距离为0.12～0.83 km，平均长度为0.49 km；示踪粒子在第1 000 d 时运移距离为3.1～7.4 km，平均长度为5.6 km，如图3（b）所示。

由反向示踪粒子流线图可知：抽水井群的西北边界100 和1 000 d 运移流线的长度相对于东南边界较长。经分析可知：粒子轨迹方向与地下水位落差方向基本平行，抽水井群东南边界及该区域的水井受到井群内部取水造成的井群内地下水位坡度平缓的影响，其粒子运移距离相对地下水井群西北边界及该区域水井的粒子运移距离较短。

此外，通过入渗场补给水源时，计算过程同工况1，各单井设置反向示踪粒子，计算其运移距离。滹沱河在南水北调干渠与滹沱河相交断面上游建设长约5.5 km、宽约0.6 km 的入渗场，以补给地下水，入渗场的补给水量为2×108m3/a。滹沱河傍河拟建水源井分布在石家庄地下水库范围内，假设区域自黄壁庄水库方向没有地下水量补充，而是由滹沱河内入渗场作为主要的地下水补给来源。在此情景下滹沱河傍河地下水井的粒子运移流线将发生改变，经过反向示踪粒子法计算，各拟建抽水井1 000 d 运移距离为2.4～3.5 km，平均长度为3.0 km。在滹沱河入渗场补给水量为2×108m3/a 的影响下，地下水井群局部地下水流场及水力坡度将会发生剧烈变化，此种情境下17 口拟建抽水井的补水来源主要是滹沱河入渗场，相应的反向示踪粒子的流线方向也将变更为入渗场水源补给水源地的方向。此时，抽水井群1 000 d 的粒子流线范围在地下水井群南部边界与滹沱河入渗场北部边界之间。同理，反向示踪粒子在100 d 的运移流线也是在入渗场范围内，如图3（c）所示。

4 地下水保护区范围的确定

在将地下水井按照设定的不同工况条件概化时，通过反向示踪粒子可得到污染物的迁移距离，见表2。

表2 不同工况下的模拟计算结果Tab. 2 Simulation results under three working conditions

在工况1 中，当抽水井作为基准点单独设置示踪粒子进行计算时，100 d 的粒子运移流线长度范围为0.10～0.83 km。又因各抽水井彼此相距0.80 km，部分0.10 km 的迹线范围可能导致相邻水井的一级保护区之间会有闲置区域。如果闲置区域存在人为活动，那么会增加水源地受到污染的风险，没有起到相应的保护作用。在工况2 中，由于示踪粒子的运移是随着地下水流场的流动而进行的，运移流线的不规则分布和长度的差异可以直观地表现出抽水时水源地对附近地下水流场的影响；同时，示踪粒子反向运移的轨迹几何形状较工况1 更加规则，且迁移距离相对较短，能够有效保证土地资源的利用。在入渗场补给水源地的条件下，滹沱河傍河地下水井的粒子运移流线将发生改变，100 和1 000 d 的运移流线都在入渗场范围内。综上所述，选取工况2 的模拟情景，并将特殊补给情况下的入渗场区作为一级保护区的补充范围，可以做到科学合理划分保护区，满足水源保护的要求。

因此，根据工况2 模拟计算结果，将水源地一级保护区的边界范围设为：西侧上游从水源地边界外推0.49 km，东侧下游外推0.22 km。北侧外推到河南岸大堤外坡脚，南侧从边界外推0.12 km；同时，考虑特殊补水情况，故入渗场整个区域也要作为一级保护区的范围。污染物在1 000 d 内运移到水源地的最长距离为7.4 km，而滹沱河地下水源地上游边界距离黄壁庄水库大坝只有15.9 km，因此，二级保护区从一级保护区上游边界直接延伸到黄壁庄水库坝下。两侧及下游二级保护区范围主要参考数值模拟计算结果确定，上游外推到黄壁水库坝下，下游以一级保护区边界外推1.0 km；北部边界考虑河堤的物理边界作用，从河流北岸外推50 m；南部边界从井群边外推3.1 km。水源地保护区示意图如图3（d）所示。

5 结 语

以石家庄市滹沱河应急后备水源地为例，采用数值模拟法设置了2 种工况，对水源地污染物运移轨迹进行模拟计算，进而划分一级、二级保护区。主要结论如下：

（1）利用MODFLOW 和MODPATH 等程序，建立了研究区的水文地质模型，模拟出地下水流场，每口抽水井单独设置反向示踪粒子和17 口水井概化成1 个水井群两种不同的工况来计算反向示踪粒子的流线，通过比较2 种工况下100 和1 000 d 粒子反向示踪运移流线可知，将水源地概化成井群时反向示踪粒子运移轨迹更加精准，运移距离相对较短，因此选择工况2 并考虑入渗条件，更利于水源地的环境保护。

（2）考虑入渗条件，反向示踪粒子在1 000 d 的运移距离完全在入渗场的场区范围内，建议当地相关主管部门加强对入渗场的保护，采取相应措施制定水质安全突发事件的预案，以保证应急后备水源地安全，真正起到应急和后备的作用。