“引哈济党”工程对敦煌盆地地下水位影响的数值模拟研究

何剑波，李玉山，胡立堂，尹 政，胡彦斌

（1.甘肃省地矿局水文地质工程地质勘察院，甘肃 张掖 734000；2.甘肃省地下水工程及地热资源重点实验室，甘肃 兰州 730050；3.北京师范大学水科学研究院，北京 100875）

敦煌盆地是我国典型的内陆河盆地，是古陆上丝绸之路的重要途经之地，气候极其干旱，生态脆弱。近半个世纪以来，敦煌绿洲区已经出现湿地萎缩和泉水衰减等生态环境问题。针对敦煌水资源开发利用和生态环境存在的突出问题，2011年6月国务院批准了《敦煌水资源合理利用与生态保护综合规划（2011—2020）》（简称《敦煌规划》）[1]。其中规划的“引哈济党”工程是保证《敦煌规划》成功的重要骨干工程之一，即引与敦煌盆地相邻的苏干湖盆地大哈尔腾河补给党河[2]，并在党河洪积扇进行生态补水。随着《敦煌规划》的实施，党河灌区地下水水位出现上升，而且2016—2018年西湖湿地可见水面的面积增大，“引哈济党”工程实施的必要性受到质疑。

敦煌盆地地下水的研究一直受到政府和国内学者关注。李世明等[3]在国家“九五”攻关项目中对水资源合理利用和生态环境保护问题进行了总结。2003年国土资源部开展的第二轮地下水资源评价[4]和地质大调查[5]，以及中国工程院重大咨询研究[6]，均对地下水资源进行了评价。其中，月牙泉区作为著名景点，泉水湖面面积的变化是地下水系统衰减与否的直接体现[7]。敦煌西湖湿地是国家级自然保护区，是典型的干旱区湿地生态系统和荒漠生态系统[8]，其中水资源是制约西湖湿地的重要和关键因素。学者对地下水补给及湿地变化规律及驱动因素[9−10]关注较多，包括同位素技术评价地下水补给特征[11−12]、土地利用与地下水和生态关系[13−14]、河水和地下水转化关系等[15]，但对规划实施以来地下水对生态的支撑作用分析得不足，尚没有专门针对“引哈济党”工程对地下水水位影响分析的论文。本研究组近几年开展了水文地质补充调查和地下水动态监测，为深入分析敦煌盆地地下水储量变化和“引哈济党”工程对地下水水位的影响提供了条件。

敦煌盆地除敦煌绿洲外，大部分地区是戈壁，地下水动态监测信息缺乏，建立可靠的地下水模型是重要的研究基础[16]。本研究旨在综合利用已有的和补充调查的地下水水位动态、统测地下水水位以及水文地质勘探信息，建立地下水数值模型，模拟区域地下水流场，并分析2010年以来敦煌盆地地下水储量的时空变化规律，重点估算盆地范围、党河灌区和西湖湿地的地下水储量变化，预测“引哈济党”工程不同调水量情景下敦煌盆地地下水水位的动态，该结果将为进一步论证“引哈济党”工程提供技术支撑。

1 研究区概况

敦煌盆地位于甘肃河西走廊西端，地势南部高，北部低，海拔为840～1 700 m。东部山前洪积平原坡度较大；西部洪积平原和中部冲洪积平原的地势比较低缓；西部边界地势最低，约为840 m，见图1（a）。其为典型的大陆性气候，降水稀少，蒸发强烈。盆地多年平均降水量为40.1 mm，多年平均潜在蒸发量为2 486 mm。

图1 研究区平面（a）和剖面图（b）Fig.1 (a)Map of the study area and(b)a profile

区内主要河流为疏勒河和党河。疏勒河自昌马水库流出后经玉门市和瓜州、敦煌等市县，再流经哈拉湖，向西注入尾闾哈拉齐，地质时期曾流入罗布泊。根据水文站监测资料，疏勒河双塔站平均流量为24.64 m3/s，在双塔水库的下游基本断流。根据1955—2017年沙枣园水文站资料，党河进入敦煌盆地多年平均径流量为3.17×108m3。除党河外，敦煌盆地还有南湖泉水、崔木土沟、多坝沟、东水沟、西水沟等小沟小河从盆地南部进入，疏勒河在盆地东部瓜州与敦煌交界处有很小的潜流进入，在盆地西部有很小潜流流出。敦煌西湖国家级自然保护区位于疏勒河下游故河道的湾窑—马迷兔—玉门关一带，是一个以湿地生态系统为主要保护对象的区域。根据疏勒河河道疏浚之前的实地调查，疏勒河干流在后坑以东为间歇性河道，在后坑以西为干涸河道，双塔水库建成运行后，玉门关一带已很少有河水流入。杨根生等[17]遥感影像研究显示，1973—2007年，西湖自然保护区沼泽面积从184.2 km2减小到107.5 km2。同时根据监测，月牙泉湖水位从1960年的1 139.55 m 降落到2002年7月份的1 133.19 m，平均每年下降0.15 m，1999年出现干涸。

2 研究方法和数据

2.1 水文地质概念模型

根据程旭学等[5]研究和本项目组补充的地质钻孔与水文地质试验孔，进而将含水层概化为1 个潜水含水层和2 个承压含水层，见图1（b）。平面上，模型范围为敦煌盆地。垂向上，以潜水面为顶面，以半胶结的下更新统的底面作为模型的底边界，厚度一般约50～300 m。在党河洪积扇顶部，出山口河水渗漏补给，地下水经水平运动，再以垂向为主地向上流动。区内的地下水在洪积扇群带为单一的潜水，地下水接受垂直向下的河流、渠系和灌溉回归水等的入渗补给，然后逐渐在细土平原带转为多层的含水系统，地下水在洪积扇群带前缘形成浅埋区。在党河洪积扇接受地表水入渗补给后，地下水向洪积扇外呈放射状的径流运动。研究区地下水概化为含6 个模拟层（种植土层、潜水含水层、弱透水层、第一承压含水层、弱透水层和第二承压含水层）的地下水三维不稳定流模型。取火焰山和卡拉塔什塔格山山前基岩与第四系松散沉积物的分界线为模型的南边界，将其概化为弱透水边界。取北山山脚为模型的北边界，该处地下水径流很弱，处理为零流量边界。取瓜州与敦煌市行政交界处（即瓜敦公路）为模型的东边界，是流量边界类型。关于西边界地下水的流向问题学术界一直存在争议，历史上的疏勒河是出甘肃省西界流入新疆的罗布泊，尽管疏勒河河水早已断流，而且缺乏实际的勘探资料，研究过程中曾打浅层井进行地下水流向判断和流量估算，发现地下水仍向西流出边界，因此取敦煌盆地西侧库穆塔格沙漠东界为模型的西边界，将其作为第二类边界。地下水的主要补给项为降水入渗、党河入渗、渠系和灌溉水入渗补给；主要地下水排泄项为潜水蒸发、地下水开采和泉水溢出。

水文地质模型的数学模型为：

式中：H—含水层或弱透水层的水头函数/m；

H0，H1—研究区初始和第一类水头函数/m；

OUTsp—泉口处位置；

Zsp—泉口标高/m；

Kze—垂向等效渗透系数/（m·d−1）；

Kh、Kz、Kn—含水层或弱透水层的水平和垂直、二类边界法线方向的渗透系数/（m·d−1）；

μs—含水层或弱透水层的单位储水系数/m−1；

μd—无压含水层的重力给水度；

Qw、Vw—开采井的开采量和井孔工作段的体积/m3；

W—大气降水、渠系和水库等滞后入渗补给强度的代数和/（m·d−1）；

Q—研究区第二类边界已知单位面积流量函数/（m·d−1）；

B1—研究区第一类边界；

B2—研究区第二类边界；

D—研究区的分布范围。

2.2 地表水系的概化与地下水的交互

由于研究区降水量小，地表水入渗是地下水的重要补给来源，因此合理描述地表水和地下水的关系非常重要。地表水系主要包括党河、疏勒河、多坝沟、崔木土沟、南湖泉水（山水沟和西头沟）和东水沟。模型将地表水系概化为河段或引水渠道的连线对节点连接组成拓扑网络[18]，按水量平衡计算，地表水向地下水的渗漏量。河道和渠系概化为线状源汇项，田间灌溉回归、地下水开采等概化为面状源汇项。党河灌区是主要用水区，党河经党河水库调蓄后，沿渠道供给农业、工业和生活用水，党河冲淤水沿党河主河道入渗补给地下水。南湖泉水在供给南湖灌区后，其余部分消耗于蒸发和补给地下水。

2.3 潜水蒸发的模拟

在盆地北部的大片地区，地下水水位埋深在1 m以下，潜水蒸发量大，对其合理模拟是模型仿真的关键。根据清华大学玉门均衡试验场3 种土质（砂土、壤土和黏土）的潜水蒸发系数（潜水蒸发与水面蒸发的比值）与地下水水位埋深的关系（图2），在模型中根据幂函数型的拟合公式计算不同土质的潜水蒸发量。

图2 玉门试验站潜水蒸发系数与地下水位埋深的关系Fig.2 Relationship between the evaporation coefficient and depth to groundwater level in Yumen

2.4 数据

收集和利用1∶20 万地质普查报告（《敦煌幅》（1972）和《玉门关幅》（1975））、1∶20 万区域地质和水文地质普查报告（《敦煌幅》（1982）和《玉门关幅》（1982））、《河西走廊西端水文地质调查报告》（1996）、《河西走廊疏勒河流域地下水资源合理开发利用调查评价》（2008）等研究资料，新增地质钻孔15 个和水文地质钻孔36 个，参考《甘肃省河西走廊西段安敦盆地水文电测深普查成果报告》（1989），完成了三维地质结构的描述，利用克里格插值法生成各个格点的地面高程和模拟层厚度。南、东和西侧弱透水边界的流量，根据收集的资料作为估计，并在模拟识别过程中加以调整。降水入渗补给系数与地下水埋深和包气带岩性相关，一般在0.06～0.50 之间[5]。研究区共有22 眼长观孔和16 眼2016—2019年勘探的观测孔，均为潜水观测孔，见图1（a），其水位数据用于模型参数率定和模型验证。2018—2019年，研究组对敦煌盆地地下水水位进行了统测。收集到2000—2018年月尺度的党河下泄水量数据，下泄水量总体呈增加趋势。收集了2016—2018年双塔水库下泄水量数据，近3年，年均下泄水量达到2.1×108m3。

3 模型识别与验证

模型计算采用了中国地质大学（武汉）陈崇希等[19]开发的多边形有限差分数值模型（PGMS），并对地表水的水量平衡模块进行修改。PGMS 模型曾用于疏勒河流域地下水开发利用分析[5]，相关的模型数据可以参考。模拟识别时段为2010—2017年，将2018年作为验证时段，并取时间步长为1 个月，共计108 个时间步长。每层剖分为2 089 个结点，4 002 个辅助三角形，共6 个模拟层。模拟区域总面积约为1.08×104km2。

3.1 水文地质参数分区

对于党河洪积扇中上部潜水含水层，顶部上更新统为单一的砂砾石层；党河洪积扇至北山山前逐渐变薄，党河洪积扇一带大于200 m；中更新统属单一砂砾石层，颗粒变得稍细，到扇的中部逐渐变成砂砾石、含砾砂、亚砂土、亚黏土互层的多层结构。冲积扇的中上部由于水位埋深大，上更新统多为透水不含水的砂砾石层，地下水多赋存于中更新统砂砾石层中。含水层岩性为砂砾石夹薄层泥质砂砾石及含砾中粗砂，由南向北颗粒变细，砂层增多。在北部，含水层厚度一般小于30 m。水位埋深在洪积扇的上部大于100 m，中部为50～100 m；含水层渗透性强，渗透系数在50 m/d以上。研究区北部的地下水埋深较浅，导水性较差，西端疏勒河排泄带渗透系数仅有1.45 m/d。根据地质和水文地质条件以及前人研究成果，拟定了水文地质参数分区，然后从参数初始估计值出发，拟合96 个时段的各观测孔水位。在平面上，从党河洪积扇向北部和东西部共划定了7 个水文地质参数分区，水平渗透系数范围为0.056～60 m/d，垂向渗透系数取为水平渗透系数的1/5，给水度范围为0.002～0.25，单位储水系数为0.48×10−6～3.00×10−6m−1。

3.2 模型识别验证

2010年1月—2017年12月和2018年1—12月分别为模型识别和验证期。典型观测孔（B8、B17、B18、B14 孔）观测与模拟水位的对比曲线如图3，两者水头的动态趋势一致。识别验证期共1 501 个观测值的均方根误差为2.45 m。分别计算22 个长观孔水位拟合的纳什效率系数[20]、观测井观测和模拟数据的变异系数见图4。10 眼长观孔水位拟合的纳什效率系数在0.85 以上，同时22 眼孔的平均纳什效率系数为0.61。观测井（B5、B6、B7、B13、B14、B20、B21）观测数据的变异系数（CV值/%）比模拟结果的大2～16倍（图4），地下水水位变化幅度偏小的主要原因是开采井空间分布和数量的信息不准确，模型将党河灌区的开采井概化为面井，反映的是区域开采条件下地下水水位的动态，因此模拟的水位变化幅度偏小（如B14孔）。研究组根据经验和实测结果绘制了2018年12月潜水位等值线，与模拟值对比，见图5，两者在趋势上有较好的一致性，反映了地下水的流向。

图3 4 个典型观测孔观测和模拟水位对比曲线Fig.3 Comparison of the observed and simulated water levels for the four typical wells

图4 模型识别和验证期22 眼长观孔的纳什效率系数和数据变异系数比较Fig.4 Compassion of the Nash–Sutcliffe Efficiency Index and coefficient of data variation for 22 wells

图5 2018年12月模拟和观测的潜水位等值线图Fig.5 Contour map of the observed and simulated water table at the end of 2018

4 结果分析和讨论

4.1 地下水均衡和储量变化

根据2000—2018年地下水均衡模拟计算结果可知：降水入渗补给量约为0.33×108m3/a，河流入渗量达到1.87×108m3/a，渠系和田间灌溉回归水入渗为1.10×108m3/a，地下水侧向边界补给量约0.30×108m3/a，主要的入渗补给项为河流、渠系和田间灌溉水回归，约占总补给量的83%。地下水侧向流出量约为0.47×108m3/a，泉水溢出量约为0.03×108m3/a。潜水蒸发量（2.35×108m3/a）和地下水开采量（0.84×108m3/a）是主要排泄项（占总排泄量的87%以上），潜水蒸发量约占总排泄量的64%。2010—2018年平均地下水储量空间变化见图6（正值表示储量增加，负值表示储量减小）。地下水储量亏空最大的区域在党河洪积扇，党河灌区南区和北区的地下水储量变化有显著差异，灌区南部的地下水储量增加，但灌区北部的地下水储量减少，其主要原因是党河灌区开采量的影响。值得注意的是，在西湖自然保护区及研究区南部洪积扇的地下水储量仍呈小幅度减少趋势。盆地南部洪积扇地下水储量减少，主要是因为党河向地下水的补给量减少，而且灌区地下水开采也引起洪积扇地下水水位的下降。西湖自然保护区地下水储量减小是因为南部洪积扇地下水水位下降，导致西湖保护区的侧向径流量下降，而且地下水潜水蒸发量也呈现缓慢的下降。

图6 2010—2018年平均地下水储量变化空间分布Fig.6 Contour map of the yearly average groundwater storage changes in the period from 2010 to 2018

敦煌盆地2010—2018年地下水的主要均衡项变化见图7。地下水开采量由2010年的0.93×108m3下降到2018年的0.73×108m3，减少了约22%。党河入渗补给量自2010—2014年间呈现减少趋势，减小了约0.18×108m3/a；在2014—2018年呈增加趋势，增加了约0.23×108m3/a。2010—2018年潜水蒸发量的变化不大，敦煌盆地地下水储量变化为−0.76×108～0.16×108m3，年平均地下水的储量减少约0.40×108m3。自2014年以来，党河径流量的增加使得其入渗量逐渐增加，区

图7 敦煌盆地主要地下水均衡项2010—2018年变化趋势Fig.7 Changes in the main groundwater components in the Dunhuang Basin from 2010 to 2018

域地下水储存量（负值表示亏损，正值表示增加）呈增加趋势。在党河灌区，2010—2018年年均灌溉水回归和渠系的入渗水量约1.10×108m3，地下水侧向径流补给量约0.35×108m3，降水入渗补给量约为36×104m3，地下水开采量为0.84×108m3，潜水蒸发量约为0.64×108m3，年均地下水储量减少约2.62×106m3。2000—2018年，西湖自然保护区核心区地下水侧向径流补给量为0.85×108m3/a，降水入渗补给量为0.56×106m3/a，地下水侧向流出量约1.73×106m3/a，潜水蒸发量约0.93×108m3/a，地下水储量减少约9.99×106m3/a。在敦煌盆地其他地区，地下水储存量减少约0.27×108m3/a，其主要在党河洪积扇，见图6。

4.2 引哈济党不同调水情景的地下水水位动态预测

以平水年党河来水条件和2018年降水、地下水开采及灌溉水利用方式为基础，预测“引哈济党”不同调水量的情景。模拟未来50年地下水水位动态不同调水量的方案分别为0.80×108m3/a（S1 情景）、0.90×108m3/a（S2 情景）、1.00×108m3/a（S3 情景）和1.20×108m3/a（S4情景）。模拟未来50年地下水水位动态。考虑“引哈济党”工程的调配方案，调水量在损失16.5%后进入党河水库，之后通过多级拦水坝、渗水池和渗水坑槽等工程措施在七里镇以上16 km 的河段拦蓄补给地下水。

“引哈济党”调水情景下，50年后地下水水位上升预测值等值线见图8。“引哈济党”工程调水的补水量越大，地下水水位回升0.5 m 的范围有明显增加；在河道补水地段地下水水位回升值最大，之后地下水水位回升值由党河洪积扇向西北侧的外围逐渐减小，党河洪积扇水位回升在5.0～20.0 m。4 种调水情景均可发现党河灌区的地下水水位普遍回升，回升值从河流上游的数十米到灌区北部的0.5 m，在月牙泉区地下水水位回升7.0～15.0 m。4 种调水情景下，南湖灌区附近地下水水位有明显恢复，为1.0～6.0 m，地下水恢复0.5 m 的范围还没有到达北部的疏勒河河道附近，西湖自然保护区地下水水位回升在0.5 m 以内。

图8 4 种调水情景预测的50年后敦煌盆地上升地下水位等值线（a，b，c 和d 分别为S1、S2、S3 和S4 情景）Fig.8 Contour map of the predicted increased groundwater levels after 50 years under four water transferring scenarios.a,b,c and d represent scenarios S1,S2,S3 and S4,respectively.

模拟的典型4 眼观测井的地下水水位变化曲线，见图9。“引哈济党”调水补水量越大，地下水水位回升值越大。4 种调水情景下，50 年后，位于洪积扇的B18 孔地下水水位上升了2 .67～3.93 m，位于研究区北侧的B17 孔地下水上升约0.42 m，党河灌区南部的B8孔和北部的B14 孔分别上升1.80～3.24 m 和约0.02 m。

图9 4 种调水情景预测的4 个典型观测孔地下水位变化曲线Fig.9 Changes in the predicted water levels for four typical wells under four water transferring scenarios

S1、S2、S3 和S4 情景下年均地下水储量分别增加约0.58×108，0.66×108，0.73×108，0.88×108m3。储量恢复最大的区域为党河洪积扇，与图6 反映的情况一致。在月牙泉附近地下水水位有着明显回升。西湖自然保护区是研究区最重要的生态保护目标之一。虽然研究区地下水储量均出现增加，但50年的情景预测西湖自然保护区核心区边界内地下水水位恢复缓慢。西湖自然保护区核心区（图1）的补给水源主要包括降水量、疏勒河径流量和地下水侧向径流量。2016年疏勒河河道疏浚后，河道沿程水头损失减少，疏勒河进入西湖自然保护区的水量主要取决于双塔水库泄水量。西湖自然保护区的地下水侧向径流量是主要的补给来源，也受人类活动影响大。根据2016—2018年模型计算，西湖自然保护区核心区的年均降水量、疏勒河径流量和地下水侧向径流量分别约为0.60×108，1.00×108，0.84×108m3。在“引哈济党”调水进行生态补水之前，党河洪积扇地下水水位下降，地下水向西湖自然保护区流动的水力梯度减小，因此侧向径流量缓慢减小，保护区的地下水水位呈现约0.01～0.05 m/a的下降趋势。模拟“引哈济党”调水，党河洪积扇地下水水位逐渐恢复，地下水向保护区的侧向径流量也逐渐增大，预测结果地下水水位恢复约为0.02 m/a，说明该保护区地下水水位变化是缓慢的，时间尺度至少在50年以上，而地下水水位的变化将对植被生态有显著影响。“引哈济党”调水工程的实施将使地下水向西湖自然保护区的侧向径流量保持稳定，从而维系西湖植被生态生长的地下水水位。

5 结论

（1）根据区域水文地质条件建立了地下水三维流数值模型，22 眼长观孔模拟水位的平均纳什效率系数为0.61，10 眼长观孔水位的纳什效率系数高于0.85，观测井水位流场的拟合情况说明模型较好反映了区域地下水水位动态规律。

（2）模拟结果显示2010—2018年区域地下水储量亏空0.40×108m3/a，党河灌区地下水储量损失约2.62×106m3/a，西湖自然保护区的核心区地下水储量损失约9.99×106m3/a，而且党河洪积扇的地下水储量衰减最大，党河灌区南部的地下水储量增加但在其北部的地下水储量减少。

（3）模拟了调水情景分别为0.80×108，0.90×108，1.00×108，1.20×108m3/a 的地下水水位动态，发现了50年后地下水水位在区域上有差异性回升，党河洪积扇地下水水位恢复在5.0～20.0 m，月牙泉区地下水水位回升7.0～15.0 m，而西湖自然保护区地下水水位恢复在0.5 m 以内。

（4）地下水侧向径流量是维系西湖自然保护区核心区植被生态生长的重要水源，人类对党河的利用以及在党河灌区的人类活动造成党河洪积扇地下水水位下降，从而减少了地下水向西湖保护区的侧向径流量，这种变化的时间响应是缓慢的，而且预测分析的“引哈济党”工程对西湖自然保护区地下水水位影响也是缓慢的，但这种调水工程最终会恢复党河冲洪扇区亏空的地下水储量，进而维持地下水向西湖保护区侧向径流量的稳定回升，保证植被生态增长的地下水水位要求。

本次重点研究的是 “引哈济党”工程实施前后的地下水水位变化规律。由于研究区西部区域偏远，只有间断的和统测的地下水水位资料，缺乏长期地下水水位动态监测井，难以建立地下水与植被生态的动态关系，仍需进一步开展相关工作，以支撑敦煌规划中适宜调水量的论证。

致谢：感谢甘肃酒泉市党河流域水资源管理局在调查和研究过程中给予的帮助和支持！