1. 引言
为了贯彻落实习近平总书记关于生态文明建设以及保障水安全重要指示精神，全面解决华北平原因地下水超采引发的生态环境危机，水利部联合河北省人民政府组织开展地下水超采综合治理专项活动，联合制订了《华北地下水超采综合治理河湖地下水回补试点方案(2018~2019年)》，依据“重点突破、示范先行”的实施方案原则，在超采严重程度最突出、水源补给条件最好的典型区域选择具备条件的河段，开展目标河段地下水人工回补试验工程，可以为后期华北平原地下水超采综合治理活动全面开展积累一定的实践经验，提供借鉴。目前，河南、河北两省境内的部分河道通过南水北调中线工程沿线退水闸开展生态补水工作。
南水北调中线工程是缓解北方水资源短缺的重大战略性基础设施，滏阳河流域邯郸段平原区作为其受益区域，生态环境和水资源状况备受关注。邯郸市平原区地表水资源匮乏，农业和生活用水主要依赖地下水开采。然而，随着人口增长和城市化进程加快，地下水过度开采引发了一系列环境地质问题[1]。滏阳河作为邯郸市的母亲河，源自太行山东麓，属海河南系子牙河水系，是当地唯一的常年性河流。流域地处北温带季风气候区，降水年内分配极不均衡，全年降水量的70%集中于夏季，由此形成春旱秋涝的典型气候特征。同时，滏阳河来水受东武士水库调控，非汛期流量较小，下游断流现象频发[2]。
2015年12月起南水北调中线启动向滏阳河生态补水；2018年9月13日滏阳河纳入南水北调中线生态补水示范区，随后围绕滏阳河开展一系列综合整治，不断促进流域可持续利用。研究表明南水北调中线工程(邯郸段)引江水替代受水区生产用水，补充生态用水，对区域地下水位的回升、促进地下水环境恶化趋势的扭转起到了重要作用[3]。
在地下水的研究与应用中，数值模拟受到广泛关注[4] [5]，段旭东等结合自身团队优势，应用数值模拟技术对临汾盆地这一较为复杂的地下水位动态变化形式进行了全方位系统的模拟和预测研究，量化了影响地下水位变化的各项因素[6]。
在水资源短缺与生态环境需水矛盾日益突出的当下，生态补给地下水成为改善区域水文生态的关键举措，而数值模拟则是深入探究其过程的有力工具。胡立堂等以永定河流域(北京段)为研究区，利用生态补水期间日尺度监测数据，揭示地下水位恢复和动态响应的机制，为制定科学的生态补水方案提供技术参考[7]。王树芳等通过分析南水进京前后地下水水位的变化规律，研究了南水北调中线一期工程通水后对北京地下水的影响[8]。边雨佳以永定河流域平原区的北京段为例，主要应用地下水数值模拟方法，结合永定河生态补水方案，系统探究了生态补水条件下地下水的动态响应特征[9]。阿德·阿布拉依据2000~2020年塔里木河下游段地下水监测数据，拟合生态输水前后地下水水面线方程。基于水均衡原理，分析21次输水后地下水埋深时空变化，估算地下水补给量，为后续生态输水方案优化提供支撑[10]。马尧等人开展2019年永定河北京段生态补水对区域地下水补给分析时，采用点线面结合的方式，选取典型地下水位监测井、监测断面及研究区，运用水均衡法计算生态补水期间山峡段、平原段不同区段的地下水补给量[11]。孙浩收集了生态补水期间水文站和地下水监测井的观测资料，利用水均衡法估算出各河段的渗漏量，并参考62个地下水监测井的水位变动数据进行了对比分析[12]。黎光和等人在永定河(北京段)进行生态补水研究时，探讨了生态补水环境下，区域内地表水与地下水之间的转换规律，以及对玉泉山泉水补给过程的影响[13]。
因此，开展南水北调中线生态补水地下水数值模拟研究对其作用机制的进一步认识和合理的管理方法的提出都具有重要的意义。
2. 区域概况
滏阳河作为邯郸市的母亲河，源自太行山东麓，属海河南系子牙河水系，是当地唯一的常年性河流。流域地处北温带季风气候区，降水年内分配极不均衡，全年降水量的70%集中于夏季，由此形成春旱秋涝的典型气候特征。同时，滏阳河来水受东武士水库调控，非汛期流量较小，下游断流现象频发。
研究区域为滏阳河流域邯郸段平原区，沿滏阳河径流路径展布，研究区面积2691 km2，该地区地势总趋势为西、南高，东、北低。
3. 方法
3.1. Visual MODFLOW
Visual MODFLOW是Waterloo水文地质公司基于MODFLOW框架研发的三维地下水模拟平台，采用有限差分法求解非均质含水层中的非稳定流运动。该软件系统具有以下技术优势：1) 集成化的前处理模块，支持通过AutoCAD、GIS等空间分析工具进行数据预处理；2) 高效的数值求解引擎；3) 多维可视化分析功能。其建模流程包括网格剖分、参数赋值、数值求解及成果输出四个核心环节，各模块通过标准化接口实现无缝衔接，构成完整的数值模拟流程。
3.2. 有限差分法
有限差分法：常用有限差分离散偏微分方程，主要适用于含有时间和空间变量的方程。该方法常用在水资源研究中描述水流及水质问题。它将水流空间划分成网格，并求解网格上的差分方程式，计算出各个网格点的水流参数，进而得到水文过程的数值解答。有限差分法为水资源管理决策提供模拟预测工具。
4. 模型识别与验证
4.1. 含水层结构及参数划分
研究区含水介质以第四系和新近系松散沉积层为主，构成典型的孔隙型地下水系统。本次研究聚焦浅层地下水动态，考虑到其水位对气象因素(降水–蒸发)的敏感响应，将含水层结构简化为单层三维非均质非稳定流模型。
4.2. 网格剖分和边界条件
依据地下水模拟的基本流程，模型构建首先需要对研究区进行单元剖分。在地理空间数据云下载DEM高程数据，选择分辨率为90 M高程数据产品，进行下载。本研究采用ArcGIS软件绘制滏阳河流域邯郸段平原区区划图作为背景底图，并将其导入Visual MODFLOW软件中。采用规则矩形网格对研究区进行1 km × 1 km的剖分，共划分为78行、85列，垂向设置为单层结构，总计生成6630个单元格，其中有效单元格2691个，对应研究区面积2691 km²。通过ArcGIS软件提取研究区高程数据，并按软件要求格式输入模型，采用克里金插值法确定地面高程分布。研究区剖分结果如图2所示，其中灰色区域为活动单元，绿色区域为非活动无效单元。
研究区西侧受太行山、山前冲积平原天然边界的影响，作为流入边界，接受山前侧向地表径流补给，概化为流量边界；东边作为流量边界，依据达西公式计算该边界流量；南边作为流量边界，流量依据达西公式计算；北边设置为弱透水边界，概化为隔水边界；西北与东南边界为河流补给的定流量边界，研究区整体构成一个独立且完整的水文地质单元。
含水层上边界为潜水自由水面，实现与外界垂向水交换(包括降水入渗、河渠渗漏及灌溉回归等)。下边界以深层含水层顶板为界，考虑到第四系浅层与深层含水系统间存在稳定隔水层，故简化为隔水边界。
4.3. 补给和排泄
研究区含水层地下水的补给项主要有降水入渗、灌溉回归、河渠渗漏与山前侧向流入，排泄项包括人工开采、侧向流出与蒸发。
4.4. 模型识别与验证
根据建模规范，系统完成边界条件与源汇项的预处理及模型输入。数值求解器配置如下图1：采用非稳定流模拟方案，将2023年设为率定期，按月份划分为12个应力期，每期设置3个等比例时间步长(步长系数 = 1)。模型校验采用PEST自动优化与人工干预相结合的方法：1) PEST程序执行22轮参数反演(模型调用1076次)，参数取值范围限定为[1 × 10−1, 1 × 10−10]；2) 基于抽水试验数据对反演结果进行人工校核，确保参数组合符合实际水文地质特征。模拟结果通过实测值与计算值的残差分析进行可靠性评估。
Figure 1. Subdivision map of the study area
图1. 研究区剖分图
依据滏阳河流域邯郸段平原区的实际水文地质状况，综合本次研究过程中所采集获取的数据资料，同时借鉴前人在这一区域所进行的水文地质勘探相关成果，本研究将该区域划分为7个水文地质参数区，在这7个分区中，给水度分区与渗透系数分区保持一致。为确定各参数区渗透系数(K)和给水度(μ)的初始值，研究过程中系统梳理了区域内既有的水文地质资料，结合各种岩性的经验水文地质参数值，以研究区抽水试验实测资料相结合，多类资料综合对比后确定区内各参数取值区的初始参数值为宜。校核完成后的水文地质参数取值见表1。
Table 1. Hydrogeological parameter zoning table
表1. 水文地质参数分区表
参数/分区 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
Kx (m/d) |
5.37 |
0.55 |
6.88 |
2.88 |
8.86 |
0.80 |
2.30 |
Kx (m/d) |
5.37 |
0.55 |
6.88 |
2.88 |
8.86 |
0.80 |
2.30 |
Kx (m/d) |
0.54 |
0.06 |
0.07 |
0.29 |
0.89 |
0.08 |
0.23 |
μ |
0.065 |
0.041 |
0.045 |
0.043 |
0.047 |
0.05 |
0.043 |
为验证模型准确性，采用2023年1月至2024年12月8日观测井的地下水位监测数据作为测试数据。经过多次参数调整，模拟值呈散点分布(图2)，研究结果表明，根据模拟水头计算结果绘制的地下水水位变化曲线，与实测水位变化曲线呈现一致的变化趋势，模型表现出良好的拟合性能。通过对比分析模拟水位与实测水位的动态变化特征，有效验证了该模型的稳定性和可靠性，且模拟误差较小。
(a) 观测井1-YN-RZ (b) 观测井2-YN-LMG
(c) 观测井3-YN-MJY (d) 观测井12-QZ-DST
(e) 观测井19-HDX-HC (f) 观测井21-HDX-ZZQ
(g) 观测井47-CX-TP (h) (h) 观测井48-CX-CG
Figure 2. Observation hole fitting during the verification period
图2. 验证期观测孔拟合
Figure 3. 95% Confidence interval
图3. 95%置信区间图
经过多次参数调整及运行模型后，从图中模拟值散点分布可以看出，模型模拟值的点位集中分布在计算值–观测值对角线上，散点越靠近斜线，表明模拟值和观测值之间的误差越小，拟合程度越高，反之则表明拟合效果不理想。验证期末绝对平均残差为0.61 m，最大残差2.46 m，最小残差0.031 m，均方根误差为0.37 m，标准化均方根为2%，相关系数为0.99。个别偏差在95%置信区间内(图3)。模型模拟值与实际观测值之间的误差处于允许范围之内，且变化趋势基本一致。
4.5. 区域地下水均衡
基于研究区地下水数值模型，运用Visual MODFLOW的Mass Balance和Zone Budget模块，对2023年1月~2024年12月模拟期的地下水均衡量进行了系统计算，验证期的水均衡结果详见表2。
Table 2. Calculation results of regional water balance. Unit: 108 m3
表2. 区域水均衡计算结果。单位：108 m3
项目 |
均衡项 |
2023 |
2024 |
总计 |
|
|
补给量 |
比例(%) |
补给量 |
比例(%) |
补给量 |
比例(%) |
补给项 |
降水入渗 |
2.20 |
50.09 |
1.44 |
40.65 |
3.64 |
45.88 |
灌溉回归 |
1.15 |
26.18 |
1.14 |
32.18 |
2.29 |
28.86 |
河渠渗漏 |
0.81 |
18.44 |
0.73 |
20.61 |
1.54 |
19.41 |
侧向流入 |
0.23 |
5.29 |
0.23 |
6.56 |
0.46 |
5.85 |
合计 |
4.39 |
100 |
3.54 |
100 |
7.93 |
100 |
排泄项 |
开采量 |
3.19 |
97.35 |
2.58 |
97.14 |
5.77 |
97.25 |
侧向流出 |
0.03 |
0.92 |
0.03 |
1.13 |
0.06 |
1.01 |
蒸发 |
0.06 |
1.74 |
0.05 |
1.73 |
0.10 |
1.74 |
|
合计 |
3.28 |
100 |
2.66 |
100 |
5.93 |
100 |
均衡差 |
1.12 |
|
0.89 |
|
2.00 |
|
根据计算结果，研究区地下水系统总补给量为7.93 × 108 m3，总排泄量为5.93 × 108 m3，均衡差为2.00 × 108 m3。由于地下水系统补给量超过消耗量，表明研究区地下水系统处于正均衡状态。从补给来源分析，降水入渗补给在平原区占主导地位。在总排泄量中，人工开采是主要排泄方式。
模型输出与研究区水文地质条件高度吻合，模拟水头与实测值变化趋势一致。这表明该模型能够准确反映研究区地下水运动状况，为后续地下水运动趋势预测等相关研究提供了可靠的工具。
5. 生态补给对地下水的影响
5.1. 生态补水方案设计
本文利用水资源管理模型，以确保研究区基本用水需求为前提，在此基础上提出不同方案在满足其目标前提下的地下水位抬升目标，经过模型模拟，在每个方案既定目标前提下，分配水资源量，然后从整体上评估各生态补水方案对研究区地下水的影响，为后期科学制定滏阳河流域生态补水方案提供重要的支撑。
为深入探究生态补给对滏阳河流域邯郸段平原区地下水的影响，拟定以下六种方案，并依据相关规划生态补水量。所有方案中，河流地表水量、降水量、蒸发量取流域多年平均值。
方案一：汛期(6~9月)集中补水8000万m³
汛期(6~9月)集中补水8000万m³，其中汛期占60% (4800万m³，每月1200万m³)，非汛期占40% (3200万m³，每月400万m³)。
方案二：非汛期(10~5月)集中补水8000万m³
非汛期(10~5月)集中补水8000万m³，其中非汛期占70% (5600万m³，每月700万m³)，汛期占30% (2400万m³，每月600万m³)。
方案三：平均分配8000万m³
全年平均分配8000万m³，每月666.6万m³，日均22.2万m³。
方案四：汛期(6~9月)集中补水9000万m³
汛期集中补水9000万m³，其中汛期占60% (5400万m³，每月1350万m³)，非汛期占40% (3600万m³，每月450万m³)。
方案五：非汛期(10~5月)集中补水9000万m³
非汛期集中补水9000万m³，其中非汛期占70% (6300万m³，每月787.5万m³)，汛期占30% (2700万m³，每月675万m³)。
方案六：平均分配10,000万m³
年平均分配10,000万m³，每月833万m³，日均27.78万m³。
5.2. 生态补给对地下水影响的模拟结果
各方案地下水水位抬升表见表3。
Table 3. Groundwater level rise for each scheme. Unit: Meters
表3. 各方案地下水位抬升表。单位：米
监测井 |
补水前 水位 |
方案一 水位抬升 |
方案二 水位抬升 |
方案三 水位抬升 |
方案四 水位抬升 |
方案五 水位抬升 |
方案六 水位抬升 |
1-YN-RZ |
13.52 |
0.52 |
0.85 |
0.88 |
0.86 |
0.89 |
0.89 |
2-YN-LMG |
34.52 |
0.42 |
0.46 |
0.27 |
0.67 |
0.88 |
0.52 |
3-YN-MJY |
24.86 |
0.57 |
0.50 |
0.53 |
0.59 |
0.60 |
0.46 |
12-QZ-DST |
26.17 |
0.55 |
0.37 |
0.24 |
0.61 |
0.60 |
0.71 |
19-HDX-HC |
45.23 |
0.64 |
0.55 |
0.29 |
0.68 |
0.59 |
0.33 |
21-HDX-ZZQ |
50.03 |
0.47 |
0.48 |
0.21 |
0.62 |
0.55 |
0.58 |
47-CX-TP |
50.92 |
0.65 |
0.39 |
0.11 |
0.52 |
0.82 |
0.98 |
48-CX-CG |
62.13 |
0.57 |
0.43 |
0.10 |
0.50 |
0.51 |
0.64 |
平均抬升 |
0.55 |
0.51 |
0.31 |
0.63 |
0.68 |
0.64 |
生态补给措施实施之后，地下水位整体上呈现出上升的趋势。其中方案五水位抬升最优，抬升了0.68米。这一结果进一步证实了生态补给措施的积极影响。
5.3. 河道渗漏量
河道渗漏量是表征地表水与地下水相互作用的关键指标，其数值直接反映了河水对地下水的补给强度。该数值可作为评估人工补水方案效果的重要依据，检验生态补给地下水的效果。各方案河道渗漏量见表4。
方案一汛期(6~9月)集中补水8000万m³，河道渗漏量为4738万m³/年；方案二非汛期(10~5月)集中补水8000万m³，河道渗漏量为4686万m³/年；方案三平均分配8000万m³，河道渗漏量为4488万m³/年；方案四汛期(6~9月)集中补水9000万m³，河道渗漏量为4808万m³/年；方案五非汛期(10~5月)集中补水9000万m³，河道渗漏量为4869万m³/年；方案六平均分配10,000万m³，河道渗漏量为4878万m³/年。
Table 4. Ecological water replenishment capacity of each scheme
表4. 各方案河道渗漏量
方案 |
河道渗漏量(万m³/年) |
方案一 |
4738 |
方案二 |
4686 |
方案三 |
4488 |
方案四 |
4808 |
方案五 |
4869 |
方案六 |
4878 |
6. 结论与建议
本文以滏阳河流域邯郸段平原区为研究区域，在分析研究区水文地质概况和含水层结构特征的基础上，运用Visual MODFLOW软件建立了研究区的地下水数值模型，采用有限差分法的求解方法，对模型参数进行了率定和验证。根据假定的研究区生态补水方案，对生态补给地下水动态变化进行预测。结果表明，根据地下水位上升和河道渗漏量综合考虑，方案五非汛期(10~5月)集中补水9000万m³使得非汛期地下水水位得到有效提升，地下水位抬升0.68米。若要探究“南水北调中线”来水对邯郸市平原区乃至更广阔区域地下水的影响，建议扩大模型覆盖范围，搭建大区域地下水流模型，以此系统分析生态补水工程对地下水的影响。
基金项目
河北省高等学校科学技术研究项目资助，CXY2024037。
NOTES
*通讯作者。