薛聯(lián)青，廖淑敏，劉遠洪，魏 卿，符芳兵

(1.河海大學水文水資源學院，江蘇 南京 210098；2.石河子大學水利建筑工程學院，新疆 石河子 832003;3.濟寧市水務綜合執(zhí)法支隊，山東 濟寧 272000)

自20世紀70年代大西海子水庫建成后，塔里木河下游河道斷流，不合理的水資源開發(fā)和極端干旱氣候加劇了生態(tài)環(huán)境惡化[1-3]。為改善塔里木河生態(tài)環(huán)境，自2000年起進行間歇性生態(tài)輸水，至2018年1月輸水總量達70億m3[4]。塔里木河下游植被得到恢復，地下水埋深抬升至6 m以上[5-6]。國內(nèi)外有關地下水對生態(tài)輸水的響應研究多集中在模型建立、地下水運動方程的開發(fā)上[7-9]，對輸水帶來的累積效應研究有待進一步加強。認識輸水對地下水的累積效應可為輸水策略的中長期規(guī)劃提供科學指導，對節(jié)約干旱區(qū)水資源和塔河下游生態(tài)環(huán)境的恢復重建具有重要意義。

本文基于塔里木河下游2011—2017年進行的第12～18次生態(tài)輸水，建立阿拉干斷面局部三維地下水數(shù)值模型模擬7次間歇性輸水條件下地下水的時空分布格局，明晰間歇性輸水對地下水的有效影響范圍，評估現(xiàn)有生態(tài)輸水對地下水產(chǎn)生的累積效應，確定輸水對地下水恢復的有效影響范圍。

1 數(shù)據(jù)來源及模型建立

塔里木河下游大西海子水庫以下共有英蘇、老英蘇、博孜庫勒、喀爾達依、阿拉干、依干不及麻、庫爾干和臺特瑪湖8個監(jiān)測斷面，其中，英蘇和喀爾達依監(jiān)測斷面位于其文闊爾河河道，老英蘇、博孜庫勒監(jiān)測斷面位于老塔里木河河道，阿拉干監(jiān)測斷面位于兩河道交匯處。由于各監(jiān)測斷面所處的地形地勢有較大差別，只有分區(qū)構建模型才能較為準確地模擬塔里木河下游的地下水動態(tài)變化。以平原區(qū)的阿拉干監(jiān)測斷面為例，根據(jù)阿拉干監(jiān)測斷面的水文地質條件，利用實測地下水位數(shù)據(jù)和水文氣象資料，以ArcGIS 10.4和GMS 10.4.5為工具，對2011—2017年第12～18次輸水期間地下水在阿拉干斷面的非穩(wěn)定運動的狀態(tài)進行模擬。

1.1 數(shù)據(jù)來源

1.1.1 氣象水文數(shù)據(jù)

2001—2019年鐵干里克氣象站的逐日數(shù)據(jù)集合從中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)下載，包括降水量、平均風速和蒸發(fā)皿蒸發(fā)等。生態(tài)輸水量來自大西海子水庫下泄水量，研究區(qū)有2條支流——北側的其文闊爾河(約185 km)和南側的老塔里木河(約145 km)，2條支流在阿拉干監(jiān)測斷面處匯合，最終注入臺特瑪湖。生態(tài)輸水數(shù)據(jù)由新疆塔里木河流域管理局提供。地下水數(shù)據(jù)使用的是塔里木河管理局提供的2011年1月至2017年12月各監(jiān)測斷面生態(tài)監(jiān)測井每隔4 h記錄一次的地下水埋深數(shù)據(jù)。

1.1.2 地質高程數(shù)據(jù)

輸水過程中對地表水、地下水、植被響應等進行了野外監(jiān)測，并對模擬區(qū)域開展了相關水文地質調(diào)查工作。模型建立采用的空間分辨率為30 m×30 m的ASTER GDEM數(shù)字高程數(shù)據(jù)來源于地理空間數(shù)據(jù)云，是目前唯一覆蓋全球陸地表面的高分辨率高程影像數(shù)據(jù)。

1.2 地下水數(shù)值模擬模型

MODFLOW是模擬和預測地下水狀況以及地下水/地表水相互作用的國際標準模型，被廣泛應用于建立三維地質模型，模擬地下水運動。GMS(groundwater modeling system)是模擬地下水水環(huán)境的軟件系統(tǒng)，支持USG、LGR、NWT、2000、2005等MODFLOW版本。

1.2.1 間歇性輸水下局部地下水模擬

1.2.1.1 概念模型導入有限差分模型

由于河流雙側地下水的對稱性及模擬區(qū)域的典型性，模擬區(qū)域的邊界以典型平原區(qū)斷面阿拉干監(jiān)測斷面處的河道為中心，向河道兩側各擴展5.0 km，由ArcGIS裁剪緩沖區(qū)。該區(qū)域包含老塔里木河(河道長度19.29 km)、其文闊爾河(河道長度15.97 km)、塔里木河(河道長度14.38 km)，為不規(guī)則區(qū)域，如圖1所示。該區(qū)域有6個生態(tài)監(jiān)測井(H1、H2、H3、H4、H5和H6)，距離河道右岸垂直距離分別為50 m、150 m、300 m、500 m、750 m和1 050 m。

圖1 塔里木河下游模擬區(qū)域示意圖Fig.1 Sketch map of simulated area in the lower Tarim River

a.邊界條件設置：模擬區(qū)域取自距河道兩側5 km、距阿拉干監(jiān)測斷面兩側5～7 km的流域，輸水對河道兩側5 km外的區(qū)域影響甚微[10]，將模擬區(qū)域的四周邊界線概化為通用水頭邊界(GHB)，模擬區(qū)域內(nèi)的河道在過水時均被概化為第三類邊界[11]。

b.高程數(shù)據(jù)插值：選擇自然鄰域法進行插值。

c.初始水位計算[12]：2010年12月底的水位應作為模型模擬的初始水位，由于2010年11月12日至2011年1月6日沒有進行輸水，2011年1月初的水位與2010年12月底水位并無明顯差異，取2011年1月1日監(jiān)測井平均地下水埋深6.5 m計算的地下水位作為初始水位。

1.2.1.2 模擬區(qū)域水文地質條件概化

由于整個塔里木河下游的地下水縱向運動滯緩且與河道走向一致[13]。概化的大厚度潛水含水層垂向厚度取為40 m，分為4層。塔里木河下游灘槽高差1～3 m，河床寬約20～30 m，坡度約1∶4 500至1∶7 500[14]。根據(jù)各含水層的巖性、導水性和相關試驗，對模擬區(qū)域參數(shù)賦初始值。結合已有研究與試驗[11,13-14]，模擬區(qū)域水平滲透系數(shù)取2.5 m/d，垂直滲透系數(shù)取2 m/d；水平各向異性系數(shù)取2.0；單位儲水量取0.000 1/m；給水度取0.24。模擬區(qū)域331.5 km2，由GMS自動生成250 m×250 m正方形網(wǎng)格，有效單元格5 304個。模擬時間為2011年1月1日至2018年1月4日，共2 561 d，應力期根據(jù)輸水期進行劃分，分為22個應力期。

1.2.1.3 源匯項計算

研究區(qū)域多年平均潛在蒸散發(fā)量約為多年平均降水量的70倍[15-16]，降水在入滲補給地下水前就已消耗殆盡，垂向的補給為輸水期河道的滲漏。模擬區(qū)域遠離居民聚集地，無農(nóng)業(yè)生產(chǎn)生活用水，僅考慮河道補給。主要排泄方式為潛水蒸發(fā)、植被蒸騰和地表水蒸發(fā)[17]，由于模擬期間阿拉干監(jiān)測斷面的地下水埋深小于4.5 m，此時潛水蒸發(fā)趨近于零，實際蒸散發(fā)采用基于互補相關理論的平流-干旱模型進行計算。河流數(shù)字化后提取河道高程，輸入如表1所示的斷面流量作為河道流量。

1.2.2 不同應力期最大蒸散發(fā)速率

根據(jù)鐵干里克水文站的實測日氣象數(shù)據(jù)，采用基于互補相關理論平流-干旱模型[18]計算鐵干里克站2011—2019年實際蒸散發(fā)量[19]。不同應力期的最大蒸散發(fā)速率按照距離河道遠近進行賦值[20]。

(1)

(2)

(3)

式中：Ep——潛在蒸散發(fā)量；Ew——濕潤環(huán)境蒸散發(fā)量；Ea——實際蒸散發(fā)量；Δ——溫度-飽和水汽壓曲線斜率；Rn——地表凈輻射；G——土壤熱通量；γ——干濕表常數(shù)；α——經(jīng)驗系數(shù)；E——充分供水時的實際蒸散發(fā)量；e——水汽壓；Eab——干燥力。Ep通過EToClaculator計算，G≈0，Δ、Rn、γ、e、Eab參數(shù)的詳細計算見文獻[18，21]。

1.2.3 參數(shù)調(diào)整及靈敏度分析

監(jiān)測井水位觀測值、模擬值見圖2。監(jiān)測井H1、H2、H3、H4、H5和H6的觀測水位和模擬水位的平均誤差σ分別為0.36 m、-0.27 m、-0.42 m、0.56 m、0.65 m和0.05 m，均方根誤差分別為1.097 m、0.418 m、0.389 m、0.620 m、0.787 m和0.328 m。其中，H2、H3模擬水位比實測水位低，H4、H5模擬效果較差，H6模擬效果最好，整體模擬結果可行，模型識別比較成功。

圖2 阿拉干監(jiān)測斷面監(jiān)測井觀測水位與模擬水位對比Fig.2 Simulation level and monitoring level comparison of monitoring wells in Alagan section

為對已識別過模型的不確定性進行量化，需要對模型參數(shù)進行靈敏度分析。分別對河道水位Hr、水平滲透系數(shù)KH、垂直滲透系數(shù)KV、給水度SY、單位儲水量SS和水平各向異性系數(shù)A這6個參數(shù)進行分析，取參數(shù)增加、減少20%和10%分別與參數(shù)無變化時的地下水位做比較，對地下水位變化量取絕對值繪制監(jiān)測井的水位變幅，如圖3所示。由圖3可知，大幅度的參數(shù)變化對地下水變化影響較小，變幅20%時水位變化量小于0.4 m，采用的水文地質參數(shù)可用于地下水數(shù)值模擬。

圖3 參數(shù)變化下各監(jiān)測井地下水位變幅Fig.3 Variation amplitudes of groundwater level of each monitoring well under different parameters

2 生態(tài)輸水對地下水的累積效應

輸水對阿拉干斷面地下水的累積效應通過地下水輸水前后的時空響應規(guī)律揭示，通過含水層等水位線時空分布格局和地下水變化量時空分布格局分析地下水的響應范圍，包括地下水距河道的響應距離和受水面積。為表達方便，將老塔里木河河道稱作南側河道，其文闊爾河河道稱作北側河道，塔里木河河道稱為塔河河道。

2.1 地下水距河道的響應距離

運行所得模擬區(qū)域的含水層等水位線分布如圖4所示，圖中的日期選擇輸水前一天和輸水最后一天。等水位線越密集表示水位變化較快，輸水補給作用強烈；越稀疏表示水位變化緩慢，輸水補給作用弱化。對模擬的7次輸水來說，輸水前后地下水位基本以河道為中心，向河道兩側遞減，形成以河道為脊的分水嶺。

圖4 第18次輸水前后含水層等水位線分布Fig.4 Water level distribution of water aquifer before and after 18th ecological water conveyance

輸水后南北河道兩側1 000 m的響應距離內(nèi)地下水位抬升1～2 m(表2)，地下水位為812 m、814 m、816 m距河道的最大響應距離分別為1 100 m、800 m、600 m(表3)。隨著輸水的持續(xù)進行，地下水的抬升幅度和抬升范圍均有增大，輸水前的高水位覆蓋范圍也逐漸擴大，體現(xiàn)在高水位等水位線距河道寬度變大，范圍擴大，直觀地展示了輸水帶來的累積效應，如圖4所示。

表2 輸水前后距河道不同距離的地下水位及其抬升量

表3 輸水后地下水的有效響應距離 單位：m

值得注意的是，由于河道滲漏和蒸散發(fā)影響，兩河交匯處的塔里木河在模擬區(qū)域下游盡頭由河道損耗導致的水位降低十分明顯，輸水后812 m以上的地下水距河道的響應距離分布在距河道兩側400～1 000 m內(nèi)，816 m以上的地下水距河道的響應距離分布在距河道兩側150～500 m內(nèi)，越靠近尾閭湖范圍越小。由多階段(第12、14、18次)輸水結果得知，由于Ⅰ階段河道過水補充了含水層且輸水間隔小于2個月，減緩了后續(xù)階段輸水前地下水位的下降，說明小間隔多階段的持續(xù)輸水可以保持地下水位有效抬升，也能繼續(xù)維持地下水的有效響應距離。

2.2 地下水受水面積

由圖5可知，第18次輸水后，低水位(低于807 m的水位)的受水面積由250.1 km2縮減至127.0 km2，減少了49.22%；中水位(介于807～813 m的水位)受水面積由81.4 km2增加至151.8 km2，增長了86.49%；高水位(高于813 m的水位)受水面積由0增加至52.7 km2。地下水低水位受水面積減小，中高水位受水面積增大，中水位受水面積逐漸超過低水位受水面積，占模擬區(qū)域的45.78%。第12次Ⅰ階段輸水后中高水位受水面積為低水位受水面積的0.32倍，第18次輸水后中高水位受水面積達到低水位受水面積的1.61倍。

圖5 輸水后不同區(qū)間地下水的 受水面積變化Fig.5 Variation of affected areas of groundwater in different zones after water conveyance

2.3 地下水位變化量

為量化輸水后地下水位變化幅度、變化范圍，由圖4得到輸水后相對輸水前地下水位變化量，變化量構成的等值線沿河道對稱分布，河道附近的變化量等值線與輸水后含水層等水位線圖所展示的等水頭線分布高度重合。模擬期的地下水位變化量在-2～6 m以內(nèi)，其絕對值呈距河道由近及遠逐漸減小的規(guī)律。

表4為各階段輸水后與輸水前地下水位變化量距河道的響應距離，由表4可知，多階段輸水的第12、14、18次輸水的第一階段地下水位變化量為1～3 m，距河道最大響應距離由第12次輸水的300 m增加至第18次輸水的600 m；第二階段地下水位變化量增加至1～4 m，距河道最大響應距離可達1 100 m，響應距離比第一階段擴大。地下水位變化量距河道的最大響應距離呈波動性減小，累積效應對地下水的最大有效影響范圍有所收窄，是由于隨著輸水量的累積地下水位逐漸抬升至趨于穩(wěn)定狀態(tài)，累積效應變化速率由大到小逐漸趨于穩(wěn)定。

表4 輸水后地下水位變化量距河道的響應距離

3 結 論

a.模型對阿拉干斷面的地下水位模擬效果較好，平均誤差為 0.162 m，平均均方根誤差 0.607 m，水文參數(shù)合理，能夠用于間歇性輸水條件下的地下水模擬。

b.阿拉干斷面含水層地下水等水位線時空分布格局呈現(xiàn)明顯的以河道為中心，向河道兩側遞減，形成以河道為脊的分水嶺特點。

c.累積效應對地下水的有效影響范圍包括地下水的響應距離、受水面積2個指標：輸水后 812 m以上的地下水距河道的響應距離分布在距河道兩側 400～1 100 m內(nèi)；第18次輸水后中高水位受水面積是低水位受水面積的1.61倍；輸水的有效響應距離主要分布于距河道 1 000～2000 m范圍內(nèi)。