强降水条件下豫北平原地下水动态响应研究

朱洪生，王继华，豆敬峰，朱瑞霞

（1.河南省地质环境监测院，河南郑州450016；2.河南省地质环境保护重点实验室，河南郑州450016；3.平舆县杨埠镇农业服务中心，河南平舆463400）

近年来，随着全球气候变暖，水循环速度加快，汛期极端强降水事件增多，引起洪涝灾害频发。我国北方地区受地下水资源开采影响，地下水位多呈下降趋势。降水作为地下水的主要补给来源，影响着地下水资源的动态变化特征，极端强降水的输入导致地下水产生变化。目前，研究地下水动态与其影响因素的理论方法包括多元回归分析法[1-2]、时间序列分析预测模型[3-4]、灰色系统模型[5-6]和地质统计学方法[7-8]等。平建华等[9]总结回顾了地下水动态预测模型，将各种模型分为确定性模型和随机模型，并指出了各种模型的适用条件；陶虹等[10]研究了关中城市群50 a地下水动态变化及其影响因素，指出该地区地下水位整体呈下降趋势，地下水开采是其主要影响因素；危润初等[11]采用标准化降水指数研究了黑龙江建三江地区地下水位变化的突变点，指出地下水位下降受水稻种植面积以及降水影响；迟宝明等[12]利用基于遗传算法的BP神经网络模型预测了元宝山露天矿区地下水动态特征；张文化等[13]研究了气候变化与人类活动对石羊河流域地下水动态变化的影响，指出人类活动在地下水动态变化中占主导因素。上述主要是对降水的周期性、频率，以及降水的变化趋势、空间分配等的研究，而对于典型区域强降水条件下地下水动态响应的研究较少，单次强降水引起的地下水动态响应和预测研究更少。

豫北平原位于河南省北部，地下水过量开采严重，造成地下水埋深很深，包气带厚度增大。近些年，该地区汛期强降水频发，特别是2016年7月，华北地区出现5次明显的降水过程，河北北部、北京、河南等地出现了50 mm以上的暴雨天气。笔者利用豫北地区水文地质、地下水位统测资料和水质监测资料等，建立研究区地下水数值模型，对强降水条件下地下水动态响应进行研究。

1 研究区概况

研究区位于黄河以北的河南省北部。地形西高东低，高程42～1 725 m。研究区受地质构造控制明显，地貌类型分为山高谷深、山势陡峻雄伟的断块山地、堆积地貌，局部发育侵蚀、剥蚀型地貌。

研究地层主要有太古界、中太古界、古生界、中生界与新生界的古近系、新近系和第四系地层。研究区地下水类型可分为松散岩类孔隙水、碎屑岩类裂隙孔隙水和碳酸盐岩类裂隙岩溶水，其中：松散岩类孔隙水在垂向上划分为浅层地下水系统和深层地下水系统，在平面上分为漳卫河地下水系统（Ⅰ）、黄河地下水系统（Ⅱ）；碎屑岩类裂隙孔隙水可划分为安阳丘陵区裂隙孔隙水系统（Ⅰ）和汤阴火龙岗裂隙孔隙水系统（Ⅱ）；碳酸盐岩类裂隙岩溶水系统可分为焦作九里山泉岩溶水（Ⅰ）、辉县百泉岩溶水系统（Ⅱ）、鹤壁许家沟泉岩溶洞水系统（Ⅲ）和安阳小南海—珍珠泉岩溶水系统（Ⅳ）。

2 数据来源和研究方法

2.1 数据来源

地下水位监测包括2016年1月—2017年12月豫北地区浅层地下水动态监测点34个，深层地下水监测点4个；水质监测点41个，2016年4月共采集41组水样，2016年7月、9月10个加密点各采集1次共20组水样。

2.2 地下水数值模拟

为了研究2016年7—9月强降水对地下水位的影响，建立地下水数值模拟模型，对地下水位变化趋势进行研究。

2.2.1 水文地质概念模型

模拟区为整个豫北平原，总面积约21 460 km2。模拟计算的目的层为第四系全新统-下更新统孔隙含水层，含水层概化为2层，包括潜水含水层和承压含水层。模型上部边界为潜水面，为水量交换边界，包括降水入渗、潜水蒸发、河渠入渗和灌溉水回渗等，下部多为第三系碎屑岩类，处理为隔水边界。模型西部以山区与平原自然分界线为边界，为地下水流入边界。模型北侧以行政区边界为流量边界，侧向流入流出量根据达西定律计算得出。模型其他边界为河流边界，主要由黄河、金堤河和徒骇河组成。地下水接受大气降水和地下水侧向补给，整体流向为由西南流向东北，以蒸发、人工开采和地下水侧向流出为主要排泄方式。

2.2.2 地下水流数学模型

模拟区第四系松散沉积物岩性在水平、垂直方向上都有较大变化，在同一点，水平方向渗透系数在各个方向大小一致，而垂直方向渗透系数小于水平方向渗透系数，因此将计算区域含水层概化为非均质各向异性介质。模拟区地下水分为潜水和深层承压水，因此将研究区地下水流系统概化为非均质各向异性准三维非稳定地下水流系统，可用下面偏微分方程的定解问题来描述：

式中：K为渗透系数；S为自由面以下含水层的贮水率；W为源汇项；μ为潜水含水层在潜水面上的重力给水度；P 为潜水面上的降水入渗和蒸发量等；h（x，y，z，t）为水头；h0（x，y，z）为初始压力水头；h1（x，y，z，t）为第一类边界上的压力水头；f1（x，y，z，t）为第二类边界上的水分通量；Kn为边界法线方向的渗透系数；t为时间；Ω 为研究范围；→n 为有效孔隙度；Γ0为潜水面；Γ1为第一类边界；Γ2为第二类边界。

2.2.3 模型的识别和验证

模拟软件采用Modflow软件，初始（2016年6月）流场见图1，第一含水层渗透系数、给水度分区和第二含水层渗透系数、贮水率分区见图2。根据模拟区地下水位资料，将2016年6—12月作为模型识别验证期。验证期末刻（2016年12月）地下水流场拟合见图3。验证期模拟水位与观测水位拟合较好，说明含水层结构、边界条件的概化、水文地质参数的选取是合理的，所建立的数学模型能较真实地刻画研究区地下水系统特征，识别后含水层参数见表1。

图1 地下水初始流场（单位：m）

图2 模拟区分区

图3 2016年12月地下水流场拟合（单位：m）

表1 含水层参数

2.2.4 预测情景

2016年保持降水量、蒸发量、开采量等源汇项不变；其他年份保持降水量总量不变，7—9月增加降水量30%，其他月份降水量减少30%，其他源汇项保持不变。预测未来10、20 a地下水位变化情况，以及典型漏斗区地下水位变化情况。

3 强降水条件下地下水动态响应

3.1 地下水位的动态响应

利用6—8月、6—12月地下水位等值线得到地下水位变幅，将研究区浅层地下水划分为基本平衡区（-0.5～0.5 m）、缓慢下降区（-0.5～-2.0 m）、急剧下降区（＜-2.0 m）、缓慢上升区（+0.5～2.0 m）和急剧上升区（＞2.0 m）5个水位变幅分区。

2016年6月浅层地下水平均水位为61.11 m，平均变幅0.75 m，最大降幅（-12.50 m）位于濮阳市胡村乡北豆门村，最大升幅（13.81 m）位于鹤壁市浚县淇门水文站。2016年6—8月豫北地区浅层地下水位变幅分区分布见图4，其中：急剧上升区面积2 121.89 km2，占10.65%；缓慢上升区面积6 891.28 km2，占34.58%；基本平衡区面积9 505.54 km2，占47.69%，缓慢下降区面积682.53 km2，占2.90%；急剧下降区面积577.63 km2，占4.18%。与2016年6—12月豫北地区地下水位变幅分区（见图5）相比较，急剧上升区面积减少200.6 km2，缓慢上升区面积增大1 373.97 km2，基本平衡区面积减少905.65 km2。2016年6—8月浅层地下水位急剧上升区面积与缓慢上升区面积之和比2016年6—12月浅层地下水位急剧上升区面积与缓慢上升区面积之和多1 173.37 km2，说明7—8月强降水过后水位埋深明显减小。

3.2 地下水水质的动态响应

图4 豫北地区浅层地下水位变幅分区（2016年6—8月）

图5 豫北地区浅层地下水位变幅分区（2016年6—12月）

浅层地下水各单项组分变化情况见表2。根据2016年4—7月浅层地下水水质数据分析水化学类型和单个化学组分的变化情况。2016年4月HCO3-Ca·Mg·Na型水占80%；2016年7月HCO3-Ca·Mg·Na型水占70%，且与2016年4月相比井位有所改变。2016年4月、7月HCO3·SO4-Ca·Mg·Na型水都占20%，但2016年7月井位有所改变。SO4·Cl-Ca·Mg·Na型水仅在2016年7月测试中有1眼井为此类型。

单个化学组分变化情况：pH值为 7.3～8.3，与2016年4月相比，7月pH值增高井数为5眼，减少井数为2眼，不变井数为3眼。溶解性总固体含量7月低于4月的有9眼。与地下水质量标准Ⅲ类水相比，2016年4月和7月，氯化物含量均不超标，最高含量为231.49 mg/L；硫酸盐超标率均为10%；硝酸盐超标率均为50%，并且是相同井；氨氮均不超标；铁含量超标率均为60%，最高含量为5.09 mg/L；耗氧量4月超标率为10%，7月不超标。与2016年4月相比，2016年7月指标含量减少的占比相对较多，其原因是，4月为相对枯水期，到7月受降水量增加影响，浅层地下水接受降水补给，地下水受到稀释，指标含量减小。指标含量增加或不变的浅层地下水，可能是浅层地下水埋深大，降水入渗量小所致。

表2 浅层地下水各单项组分变化情况

3.3 地下水位动态响应预测

通过预测得出10 a和20 a后不同含水层地下水流场，见图6。研究区整个地下水流场发生了一定程度的变化，整体地下水位有所下降，不同地区水位下降程度不同，但整体地下水流向未发生明显变化，尤其是深层地下水。对于浅层地下水，在安阳市西北侧，10 a后出现一定面积的疏干，但深层水未疏干；在安阳市东侧和濮阳市北部，水位下降比较明显，形成较大的地下水漏斗区；在模拟区南部，水位下降不明显。研究区整体地下水位的下降说明地下水的补采不平衡，开采量大于补给量。对于当前存在的降落漏斗，如果按照当前开采量继续开采，则20 a后漏斗范围将持续扩大，地下水埋深将越来越大。

图6 含水层地下水水位等值线（单位：m）

3.4 典型漏斗区地下水位的动态响应预测

为了详细分析漏斗区地下水位变化情况，将安阳市东侧和濮阳市大部分区域的漏斗区作为典型漏斗区进行分析。将典型漏斗区作为一个均衡区，同时在漏斗区的上、中、下3个地方安置水位观察点，分别位于内黄县、濮阳市和留固镇，通过分析地下水位随时间的变化了解漏斗区水位的变化。由图7可知，不同观察点浅层地下水位和深层地下水位均出现不同程度的下降，且下降速率不同。濮阳市浅层地下水位20 a下降了1.6 m，下降速率比较平稳；深层地下水位下降了3.3 m，前2 a下降速率较快，后期下降速率比较稳定。内黄县浅层地下水位20 a变化不大，总体下降幅度不超过1.5 m；深层地下水位下降3.6 m，前期水位下降迅速，后期较为稳定。留固镇浅层地下水位20 a下降了1.8 m，深层地下水位下降了3.7 m。按照当前开采速率，典型漏斗区浅层和深层地下水位将持续下降。

图7 不同观察点地下水位随时间变化曲线

4 结 论

受强降水影响，研究区地下水位变化可划分为基本平衡区、缓慢下降区、急剧下降区、缓慢上升区和急剧上升区5个水位变幅分区。受强降水影响，地下水化学类型和水化学组分均发生了变化，水质超标率减小。保持当前地下水开采量不变，降水量总量保持不变，7—9月增加降水量30%，其他月份减少降水量30%情景下，研究区地下水位将下降，不同地区水位下降程度不同，但地下水流向未发生明显变化，尤其是深层地下水；在渗透性较好的地区，10 a后将出现一定面积的疏干，但深层地下水未疏干，在安阳市东侧和濮阳市北部形成地下水漏斗区。典型漏斗区保持现状开采，20 a后浅层地下水位和深层地下水位均将出现不同程度的下降，濮阳市浅层地下水位下降1.6 m，深层地下水位下降3.3 m；内黄县浅层地下水位下降幅度不超过1.5 m，深层地下水位下降3.6 m；留固镇浅层地下水位下降1.8 m，深层地下水位下降3.7 m。