河北平原地下水流系统演变规律研究

高永华

(河北省衡水水文水资源勘测局，河北 衡水 053000)

地下水作为人们生产生活中最重要的供水源，为工业、农业以及生活用水提供了便利[1-2]。然而大规模的开发利用地下水，造成了地下水存储量减小、水质恶化等一系列环境问题，严重威胁了人类的生存[3]。研究地下水流系统的演变规律成为亟待解决的科学问题。张翼龙[4]等对地下水流速、补给速率以及水流路径等重要信息做了详细的阐述，并给出了地下水年龄的测定方法；郭纯青[5]等认为地下水年龄数据可用来帮助认识地下水流畅，校正地下水模型，评估地下水的可更新能力；张光辉[6]认为水层底界至潜水面的非饱和带内水流过程是可忽略的，并假设地下水补给量即为根植层底部的渗出量。由于大陆冰盖融化水强烈垂向补给地下水，河北平原地下水在晚更新世和全新世发生强烈交替。地下水排泄基准面下降后，导致地下径流加快，新渗入水取代了部分深部含水层的老水。采用同位素技术研究发现，河北平原部分地区的深层淡水既不是现代降水也不是同生水，而是晚更新世到全新世更替产生的淡水。目前，由于人类的强烈活动，河北平原地下水开采天然均衡状态已经打破，地下水流模式已经发生改变。本文结合河北平原地下水的流动模式，利用地下水年龄数值模拟技术，对地下水年龄的结构进行了系统分析。通过对河北平原非饱和带中地下水滞留时间的研究，揭示了影响地下水流系统各要素之间的紧密联系，查明了不同自然条件对地下水补给的影响，为今后地下水的质量演变提供参考依据。

1 河北平原天然状态地下水流动模式

1.1 地下水年龄计算

地下水年龄分布模拟控制方程为：

(1)

式中：A为地下水年龄；t为时间；D为弥散系数；i,j为取1,2,3…，n；x为任意观测点；q为流量。

方程(1)的求解反应了扩散、对流以及弥散等地下水的混合作用，但由于古水文信息无法获得，所以一般只关心方程的稳态解。令方程(1)解为0，即：

(2)

假设公式(2)中的弥散系数D为0，即得对流年龄分布：

(3)

式中：x为任意观测点，x0为补给区起始位置。

为了尽可能祛除井管内与年轻水间的混合造成的影响，收集采用80多个同位素数据，利用14C浓度计算地下水年龄，公式为：

Age=(5730/ln(2))ln(A0/A)

(4)

式中：5730为14C半衰期；A为14C浓度；A0为14C初始浓度。

1.2 地下水年龄模拟计算

河北平原的开采初期使用MODFLOW-2000[7]建立稳定流模型，地表高程数据采用SRTM数据，分辨率为90 m，各含水层地板标高经数字化后输入LPF程序包。模型网格剖分为350行，310列，网格大小为3×3 km。网格由下到上共分为10层，最底层代表浅层含水层(第I含水层)，2～5层代表第II含水层，6～8层代表第III含水层，9～10层代表第IV含水层。采用水头边界代表东部沿海边界，假设侧向补给只发生在浅层含水层，采用流量边界定义其他边界。由于缺乏历史时期的流量数据，通常假设地下水系统处于动态平衡状态，利用未开采前的状态代表动态平衡条件，如表 1 观测井信息及13C、14C浓度所示，采用MT3DMS计算得到不同井深下的13C、14C含量值，分析方法中，LSC为液体闪烁计数器，AMS为加速质谱。

表1 观测井信息及13C、14C浓度

有效孔隙度和弥散参数是对年龄分布进行模拟时两个重要的参数，弥散度越大，说明地下水年龄差异越大，如表 2 不同孔隙度和弥散参数下计算地下水年龄统计值可知，垂向弥散度αv设在0.1～10量级上，纵向弥散度αL设在10～100 m之间， 除去极端情况其余所有弥散度下计算的地下水平均年龄值差别不大。

表2 不同孔隙度和弥散参数下计算地下水年龄统计值

2 地下水开采对区域水均衡要素的影响

2.1 河北平原地下水储存量变化

利用Zonebudget中的MODFLOW计算浅层地下水和深层地下水存储量，由图1 1975年至2015年地下水存储量变化和总补给量所示，1975年至2015年间浅层含水层地下水储存量减小速率约为30mm/yr，减少量约占总储量的15%，并且浅层含水层的变化量是深层地下水存储量的10倍，说明河北平原地下水存储量的减少主要是由于浅层地下水存储量的减少。

图1 1970年至2010年地下水存储量变化和总补给量

地下水补给量最少的年份发生在1980年和1985年，对应典型的干旱年，最大地下水补给发生在2000年。模型经校正后得到河北平原年地下水补给量为120 mm，灌溉以及降水补给为100 mm，渗流以及山前侧向补给为20 mm，降水量的变化趋势与地下水补给趋势保持一致，相关系数为0.92，地下水补给有25%来自降水和灌溉。

2.2 模型参数及校正

经原有模型参数计算得蒸散发量如表 3 原有模型参数计算得ETα以及水均衡计算蒸散发对比所示，三模型计算所得年实际蒸散发量与水均衡计算实际蒸散发量变化规律保持一致，使用模型计算的蒸散发量值与实际测量值相差不大，产生的误差值分别为24.62%，15.98%以及20.12%。三个模型在2007年出现最大相对误差，分别为39.54%，34.40%和35.43%，且降水量在小于400 mm的2001、2003、2007以及2007年，相对误差要大于平均值。

表3 原有模型参数计算得ETα以及水均衡计算蒸散发对比

图2 不同土壤地质类型地下水分布直方图

利用参数自动校正软件PEST对模型参数进行校正，固定一个参数保持不变，通过计算每年的水均衡得到的蒸散发数据来校正参数的动态变化，然后利用该数据对模型参数进行校正。如表 4 利用校正后参数计算ETα，以及与水均衡计算蒸散发对比所示，经模型参数校正之后的实际蒸散发量与水均衡计算的实际蒸散发量之间的相对误差分别为1.35%，1.87%和0.39%，相对误差大于平均值的年份，发生于1999、2005、2007以及2011年。当使用PEST软件进行自动调参时，要保证参数数目不多于观测数目，且当模型计算蒸散发与水均衡计算蒸散发之间的误差最小时，模型参数最佳。

2.3 不同土壤地质类型地下水分布

经计算后可知，年平均地下水补给量为166 mm，占总的灌溉量以及降水量的20%。土壤地质类型对地下水补给量造成较大的影响，由图2不同土壤地质类型地下水分布直方图所示，粒径大的土壤地区地下水补给量要高于其他地区：壤土分布地区地下水补给量约为5～300 mm/yr；粘性土分布地区地下水补给量约为0.3～250 mm/yr；壤砂土分布地区地下水补给量约为10～350 mm/yr。计算地下水补给量时，由于不同地区降水量及灌溉量的不同，使得不同土壤分布区域之间都存在着明显的重叠，说明土壤地质类型是影响地下水补给量的重要因素。

地下水对地表条件改变响应的快慢程度很好地反映了包气带中地下水的运动时间，包气带的厚度、地下水渗漏量以及有效孔隙度等可直接影响运动时间的大小。

表4 利用校正后参数计算ETα，以及与水均衡计算蒸散发对比

3 结语

通过建立数值分析模型，对影响地下水补给以及存储量变化的因素进行了研究，并利用地下水的年龄分布特点，在水流数值模拟基础上，对河北平原地下水的年龄分布及其影响因素进行分析，取得了一定成果，通过将土壤水均衡模型、非饱和带水流模型以及饱和带水流模型相耦合的方式，统计了不同土壤地质分布区域的降水量，并通过数据说明土壤地质类型是影响河北平原地下水补给量的重要因素，为今后地下水补给量的进一步研究提供依据。

[1]潘玥,刘勇,曾献奎,等.徐州东部废弃矿井地下水流场演化模拟研究[J]. 水文地质工程地质.2017.44(2):52-56.

[2]许乃政,刘红樱,魏峰,等. 江苏洋口港地区地下水的环境同位素组成及其形成演化研究. 环境科学学报.2015.35(12):3862-3871.

[3]王伟,黄光球. 基于格子Boltzmann的矿井巷道突水计算机仿真[J]. 煤矿安全.2011.42(5):40-44.

[4]张翼龙,王丽娟,王文中. 呼和浩特市浅层地下水水化学特征演变规律[J]. 南水北调与水利科技.2010.8(6):14-17.

[5]郭纯青,方荣杰.中国南方岩溶区岩溶地下河系统复杂水流运动特征[J]. 桂林理工大学学报.2010.30(4):507-512.

[6]张光辉,费宇红,张行南,等. 滹沱河流域平原区地下水流场异常变化与原因[J]. 水利学报.2008.(6):747-752.

[7]梁世川, 徐明, 王磊,等. GMS在地下水数值模拟及断层处理中的应用——以盖孜河水源地为例[J]. 地下水.2013.35(6):53-54.