基于改进积分型Richards方程的地下水开采方案环境影响评估研究

孔令尧

(辽宁省阜新水文局，辽宁 阜新 123000)

对于干旱半干旱区域而言，地下水占据水资源的比重均在50%以上，是区域水资源总量的重要组成。随着当前水资源矛盾的日益紧张，对于区域地下水的开采量也逐步较大，而由于不当的地下水开采，对区域地下水环境产生的影响也越来越显著。为此对不同开采方案下的地下水环境进行有效评估，确定适宜的地下水开采方案对于区域地下水环境保护十分重要。近些年来，对于地下水开采下环境影响评估的研究也逐步成为国内学者的研究热点[1- 5]。采用的方法主要为两种，一种是观测试验方式，这种方式较为直接，但这种试验方式需要大量的人力和物力；另一种方式为采用数值模型，对区域地下水进行数值模拟，分析不同地下水开采方案下的地下水模拟结果，来对环境进行评估[6- 10]。这种方式较为快捷，这种方式在地下水模拟试验中应用较多。辽西为辽宁省地下水开采量最大区域，年地下水开采量在900万m3左右，而地下水开采量的加剧，势必造成区域地下水环境的影响，为对区域地下水环境进行保护，本文结合改进的积分型Richards方程对辽宁西部某区域不同地下水开采方案下的环境进行评估，提出适宜的地下水开采方案。

1 改进的积分型Richards方程计算原理

改进的积分型Richards方程计算每个计算节点的控制水量，各节点控制水量的计算方程为：

(1)

式中，ΔWi—第i个节点控制水量，m3；Bβ—β单元的控制面积，m；Aβ—β单元的控制深度，m；ΔHi—第i个节点水头，m；Ss—节点给水率。在进行数值求解时，需要对能量方程的源汇项进行确定，源汇项方程为：

(2)

式中，Qsi—源汇等效流量，m3/s；Si—控制单元i的源汇项。在进行源汇项计算后，对控制单元的垂向流量进行计算，方程为：

(3)

(4)

(5)

(6)

式中，Δβ—控制单元β的平均水力传导度；Hs—节点s的水头，m；bi、bs、ci、cs—水力传导系数。此外，在控制单元不同方向流量计算的基础上，对各节点进行均衡分析，分析方程为：

(7)

式中，各变量同上述方程中的变量含义，Qsi和Qvi采用隐式方程进行求解。

2 实例应用

2.1 区域概况

本文以辽宁西部某区域为研究区域，区域地下基岩含水层的富集程度分布十分不均匀，地下水主要分布在风化岩裂隙以及线性充水带中，上、中、下更新层的砂砾石分布不均，水文地质较为复杂，区域单井涌量变化较为稳定，稳定在50～150m3/d，整个地下系统主要为冲击岩组成。受地下水开采量逐年增加的影响，区域地下水环境造成一定程度的影响。为对区域不同地下水开采方案进行环境影响评估，对区域内布设6组地下水观测井。各观测点的布设位置如图1所示。

图1 观测井分布位置

2.2 模型参数动态识别及验证

为对模型进行参数动态识别并对模型进行验证，结合6个区域观测点的组水文地质参数进行设定，并对比分析改进前后模型的计算质量及时间，结果见表1、2及图1。

表1 模型参数动态识别结果

表2 模型计算质量及计算时间对比

*Em表示计算质量；tc表示网格控制单元计算时间。

表1中给出了改进的积分型Richards方程下区域不同水文地质下的主要敏感参数。结合参数设置结果对比分析改进前后模型的计算质量和计算时间，从6组样本序列计算质量和计算时间的对比结果可以看出，相比于传统模型，改进模型下各组样本序列计算时间缩短约11.5倍，格网点计算质量整体提高22倍。

2.3 不同开采方案下的潜水位动态变幅分析

结合不同地下水开采方案(见表3)，应用改进的积分型Richards方程模拟分析了不同开采方案下的地下水潜水位的变幅，如图2、3所示。

表3 不同地下水开采方案

图2中方案下浅层地下水前期变幅波动性较大，后期区域稳定，各观测点整体下降率平均为1.2～1.5m，后期区域地下水变化较为平稳，在在深层区域，由于深层地下水压采量的减少，深层地

图2 方案1下地下水变幅结果

图3 方案2下地下水变幅结果

图4 不同开采方案下的区域地下水盐碱化风险评估结果

下水水位整体抬升。而在方案2下，由于同时开出深层和浅层地下水，浅层地下水水位变幅明显高于方案1下的水位变幅，各观测点整体下降水位为2.5～4.8m，而受到深层地下水开采影响，方案2下深层地下水变幅要低于方案1，不同方案下的地下水潜水位降低和抬升幅度主要受到各方案的压采率密切相关。

2.4 不同开采方案下的地下水盐碱化风险评估

结合改进积分型Richards方程实现了区域地下水的三维数值模拟，并结合地下水盐碱化风险评估指标对区域不同开采方案下的地下水盐碱化风险进行评估，各开采方案下的评估结果如图4所示。

从区域地下水盐碱化风险评估结果可看出，发生地下水盐碱化风险值较高区域地下水水位埋深为3.5m，而在潜水埋深2.5～3.5m之间的区域较易发生次生的地下水盐碱化风险。从方案对比结果可以看出，方案1下区域发生地下水盐碱化风险程度高于方案2，建议区域主要采用方案2进行地下水的开采。从区域地下水盐碱化风险分布可看出，区域的中南部较易发生地下水盐碱化风险，应尽量减少该区域的地下水开采，维持区域地下水环境健康。

3 结论

(1)改进的积分型Richards方程较传统方法的计算时效性改善程度较高，可用于非饱和水流的数值模拟。

(2)只开采深层地下水，其水位较同时开采浅层和深层地下水的方案下的水位变幅波动性更大，建议应采取浅层和深层地下水同时开采的方案，适当调整方案的压采率，保持区域地下水的平稳变化。

(3)方案1下区域发生地下水盐碱化风险程度高于方案2，区域中南部较易发生地下水盐碱化风险，应尽量减少该区域的地下水开采。

[1] 段艳华. 浅层地下水系统中砷富集的季节性变化与机理研究[D]. 中国地质大学, 2016．

[2] A. H. 阿辛顿, 钱卓洲, 付湘宁. 地下水生态系统的分析与评估[J]. 水利水电快报, 2015, 36(12): 13- 16+26．

[3] Hussein Ibrahim Ahmed. 中国北方大型灌区地下水数值模拟和预测[D]. 西北农林科技大学, 2014．

[4] 魏钰洁. 海水入侵区的地下水开采控制方法及应用研究[D]. 郑州大学, 2013．

[5] 项国圣. 黑河中游张掖盆地地下水开发风险评价及调控[D]. 兰州大学, 2011．

[6] 苗添升, 卢文喜, 欧阳琦, 等. 地下水数值模拟的不确定性分析在水质预测中的应用[J]. 水电能源科学, 2016, 34(08): 20- 23+44．

[7] 陈南祥, 杨杰, 屈吉鸿. 中牟县地下水生态水位研究[J]. 华北水利水电大学学报(自然科学版), 2016, 37(01): 84- 88.

[8] 张丽, 温斯钧, 田晓龙. 巩义市城区地下水水质评价分析[J]. 华北水利水电大学学报(自然科学版), 2016, 37(05): 68- 71.

[9] 张丽, 田晓龙. 巩义市城区地下水动态特征分析[J]. 华北水利水电大学学报(自然科学版), 2015, 36(03): 24- 27.

[10] 严锋, 郭玉法, 刘波, 等. GIS技术在地下水资源研究中的应用现状及前景[J]. 水利技术监督, 2005(05): 63- 65．