干旱灌区长时间尺度地下水数值模拟及预测

2022-12-24王志勇蒋梦瑶

水利技术监督 2022年12期

王志勇，蒋梦瑶

(1.黄河水利委员会河南水文水资源局，河南郑州 450000;2.兰州大学资源环境学院，甘肃兰州 730000)

E-mail：2072802457@qq.com

20世纪以来随着世界经济的快速发展，许多干旱地区地下水开采量越来越大，以满足日益增长的用水需求。2012年12月12日，国务院发布并实施了最严格水资源管理制度[1]。数值模拟作为地下水资源研究领域中，定量分析水量、水质不可或缺的手段，被广泛应用于最严格水资源管理制度的落实和预测[2-4]。本文运用目前国际上应用最广泛的三维地下水模拟系统软件GMS(Groundwater Modeling System)，建立了焉耆盆地地下水数值模型，并利用15个观测井的地下水位监测数据及遥感反演数据[5]进行验证和校准，分析了现状地下水均衡状态。最后，采用减少地下水开采量、提高地表水引水量的处理方案，设计了5种不同的情景模式，并依据数值模拟结果提出焉耆盆地灌溉农业的可持续开采量。

1 研究区概况

焉耆盆地位于中国西北部的新疆，如图1所示，地理坐标为41°23′—43°31′N，86°39′—88°20′E，隶属焉耆县、和静县、和硕县、博湖县以及8个兵团[9]，平均高程约为1050～1200m，占地总面积为1.32×106ha2，盆地最洼处为中国最大的内陆淡水湖—博斯腾湖，面积约为1.0×105ha2[6]。

作为典型的中温带大陆性气候区，年蒸发量约为2000～2450mm，而年降雨量只有50～80mm。开都河是盆地内最大的河流，也是最主要的地表水补给来源，多年平均径流量为36.4×108m3[7]。焉耆盆地地下水开采量的85.08%用于农业灌溉[8]，且用水效益低，属于典型的干旱绿洲灌区。随着高效节水技术的深入推广，焉耆盆地高效节水面积已从2000年的160ha2上升为2016年的1.47×105ha2，其中高效节水灌区每增加10000ha2，地下水位降低0.25m[9]，盆地内地下水年开采量达6.76×108m3，处于过度开发状态[10]。

焉耆盆地是一个完整的地下水系统，四周环山，西北、北部山地断块隆起强烈，形成高山区。盆地内可大致划分为三大地貌单元：洪积一冲积倾斜平原、开都河三角洲以及博斯腾湖湖积平原[11]。

2 建立MODFLOW模型

本文基于GMS软件建立干旱绿洲灌区地下水流数值模拟模型。模拟面积为10880km2，包括了盆地内的5个灌区和博斯腾湖泊。

2.1 水文地质概念模型

2.1.1含水层概化

根据地质钻孔资料、第四系岩性结构特征以及含水层埋藏条件，将研究区内模型模拟深度定为150～200m，模拟深度以上的含水介质分为四层，区内地下水主要开采层在第四层。地下水流系统通过模型上边界，即潜水含水层自由水面与外部发生垂向水量交换，模型的底部为第四系松散堆积物的下边界。含水层岩性呈现非均匀的变化特征，地下水水力坡度较大，渗流速度随着方向不同而发生变化，且具有垂直分量；地下水流的各运动要素随时间而发生变化。综上所述，研究区概化为四层非均质各向异性的三维非稳定流数值模型。

图1 研究区地理位置图

2.1.2边界条件和源汇项处理

依据模拟区内的水文地质条件资料，模型底边界为隔水边界；上边界参与入渗和蒸发过程，概化为开放边界。补给方式主要有山前侧向补给、河道渗漏、渠系渗漏、灌溉入渗和少量降水入渗。主要排泄方式有农业灌溉开采、蒸发蒸腾、排水渠排水、河流排泄和向博斯腾湖排泄。区内的北部和西部为接受山区侧向径流补给的山前地带，其边界概化为变流量边界(GHB)；东部和南部以被看作等势体的博斯腾湖为界，将其边界概化为已知水头边界，其水头值由博斯特湖实测水位指定，如图2所示。

图2 模型侧边界条件示意图

2.2 数学模型

根据上述地下水流场、地下水运动定律及水文地质条件等的系统分析，本文所建立的四层结构的非均质三维非稳定流数值模型的数学表达式见式(1)。

(1)

式中，Kxx、Kyy、Kzz—X、Y、Z方向的渗透系数，m/d；Kxx=Kyy；H—水头值，m；H0—初始水位，m；ε—源汇项，1/d；S—给水度，1/m；潜水含水层指重力给水度μd，第四层含水层在发生地下水开采活动时，会出现由承压转无压的状况，当呈现承压含水层时，模型计算使用弹性释水系数，而含水层在无压状态下时，模型计算使用重力给水度，其它含水层则使用弹性释水系数μe；Ω—模型模拟区域的范围；n——边界面的外法线方向；Г—侧边界；B—底部边界。

2.3 地下水数值模型

根据研究区现有资料将模拟时间定为2000—2016年，划分为17个应力期，时间步长为1天。在模型区域空间离散化的基础上，采用GMS中1000m×1000m的三维网格将模拟区划分为80个东西向和136个南北向的网格，如图2所示，上下分共4层，共计有效计算单元数43520个。求解数值模型的过程中，根据研究区相关勘察报告、水文地质条件和已有的抽水试验数据赋值模型中各分区参数的初值。在瞬态模型中，初始水头等于2000年1月稳定流模型生成的解。

2.4 模型识别与验证

以2001年1月1日到2007年12月31日作为模型参数的率定期，通过参数区划和试点法相结合的方法，依据不同时段的水位实测值对模型含水层参数反复调整，见表1，得到含水层的渗透系数范围为0～90m/d，给水度范围为0.1～0.25。以2007年1月1日到2016年12月31日作为模型参数的验证期。本文选取图2所示的15个观测井进行参数的率定及验证。选用3眼典型的长观井，对比分析其模拟水位与实测水位，如图3所示，各个观测点地下水位模拟值与实测值拟合程度较高，相关系数为0.89，水位动态变化趋势基本一致。且依据相关勘察报告可知，自2000年以后，地下水位不断下降，降幅度在2～3m之间，部分地区水位下降超过5m，与模拟结果相符。可见所建模型能较为真实地反映焉耆盆地地下含水层的实际情况，可以用来进行预测分析。

表1 含水层特性及水文地质参数取值

图3 部分观测井水位动态拟合图

图4 实测水位与模拟将结果拟合图

3 地下水开采方案及预测

3.1 地下水量均衡分析

灌区地下水的过度开采将会对农业生产和生态产生不利影响，破坏水资源平衡。在保持其他源汇项不变条件下，以现状年(2016年)6.75×108m3/a的开采量运行MODFLOW模型至2050年，得到地下水变化均衡表，见表2。

3.1.1现状水均衡

区内模拟所得现状年地下水总补给量为143658×104m3/a，入渗补给作为主要补给来源的补给量是77862.6×104m3/a，其他补给来源为河流渗漏补给、山前侧向补给、和湖泊补给，补给量分别为16664.3、46114.2和3016.8×104m3/a；区内地下水排泄量为144466×104m3/a，其中地下水开采量为67566.0×104m3/a，占总排泄量的46.2%，排水渠排泄量为7926.0×104m3/a，河流排泄量为12535.0×104m3/a，蒸发蒸腾量为56439.0×104m3/a。该数值模拟结果与已有的焉耆盆相关工作报告中均衡法计算的水均衡结果和补排关系相一致。

3.1.2水均衡量变化

运行到2050年，发现盆地内的博斯腾湖北部和开都河南部部分区域出现疏干现象，地下水位不断下降，地下水储量也减少。这也是为什么2025年的地下水开采量还没有发生变化，但到2050年已经减少到了64258.0×104m3；蒸散发量增加至2025年后也开始减小。

由此表明，焉耆盆地地下水处于严重超采的状态。为防止研究区内地下水位进一步下降，减少地下水开采量是最有效的手段。

3.2 开采情景下的地下水模拟

考虑高效节水灌溉工程的实施、控制地下水超采和用水红线指标[12]3个方面，通过减少地下水开采，结合地表水置换地下水的原理设置了5个开采方案，各开采方案的年开采量分别为6.38×108、5.35×108、4.52×108、4.29×108、4.12×108和3.83×108m3，并以2016年为现状年，将该方案及现状年的开采规模输入识别后的模型，在保持其他源汇项不变的情况下，根据各方案开采量仅对Well和Recharge参数作相应修改，运行数值模拟至规划年(2050年)，得到各开采方案下的地下水流场图，如图5所示。

表2 焉耆盆地地下水均衡表单位：104m3

图5 2050年各开采方案下地下水流场图

本文将焉耆盆地合理地下生态水位埋深下限设为4.5m，上限埋深定为3.2m[13]。按照方案1的开采规模运行到2050年，有超过65%的区域地下水埋深大于生态水位下限；方案2和方案3在规划年的地下水埋深有较为明显的上升，但只有超过45%左右的区域地下水埋深介于合理的生态水位范围内。方案4和方案5在规划年的地下水埋深60%以上都在合理生态水位之间，且与初始地下水流向及流态基本保持一致，地下水系统的流入流出量基本均衡。但考虑到国务院下达的最严格用水制度对地下水和地表水利用指标的要求，认为在有限的水资源情况下应最大限度的利用地表和地下水资源，方案4开采规模为4.12×108m3/a时，地表水引水量达12.63×108m3，下水位不会再持续降低，且逐渐恢复稳定，满足焉耆盆地用水红线指标。