“引哈濟黨”工程對敦煌盆地地下水位影響的數值模擬研究

何劍波，李玉山，胡立堂，尹 政，胡彥斌

（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）地下水側向徑流量是維系西湖自然保護區核心區植被生態生長的重要水源，人類對黨河的利用以及在黨河灌區的人類活動造成黨河洪積扇地下水水位下降，從而減少了地下水向西湖保護區的側向徑流量，這種變化的時間響應是緩慢的，而且預測分析的“引哈濟黨”工程對西湖自然保護區地下水水位影響也是緩慢的，但這種調水工程最終會恢復黨河沖洪扇區虧空的地下水儲量，進而維持地下水向西湖保護區側向徑流量的穩定回升，保證植被生態增長的地下水水位要求。

本次重點研究的是 “引哈濟黨”工程實施前后的地下水水位變化規律。由于研究區西部區域偏遠，只有間斷的和統測的地下水水位資料，缺乏長期地下水水位動態監測井，難以建立地下水與植被生態的動態關系，仍需進一步開展相關工作，以支撐敦煌規劃中適宜調水量的論證。

致謝：感謝甘肅酒泉市黨河流域水資源管理局在調查和研究過程中給予的幫助和支持！