間歇性生態輸水塔里木河下游斷面地下水位變化模擬

劉遷遷,古力米熱·哈那提,王光焰,蘇里坦,張 音

1 中國科學院新疆生態與地理研究所,荒漠與綠洲生態國家重點實驗室,烏魯木齊 830011 2 中國科學院大學資源與環境學院,北京 100049 3 新疆水利水電科學研究院水資源研究所,烏魯木齊 830049 4 新疆塔里木河流域干流管理局,庫爾勒 841000

作為干旱區生態系統的重要組成部分,植被主要依賴于地下水和地表水,而在塔里木河下游幾乎無降水,地下水則成為維系下游生態系統的唯一水源[1- 2]。20世紀70年代初期,塔里木河下游修建人工水庫,致使下游321 km河道徹底斷流,地下水位大幅度下降。為解決下游地下水位過低問題,有關部門利用博斯騰湖持續高水位的有利時機,對塔里木河下游實施應急生態輸水工程[3]。然而,由于實際條件的限制,無法實現輸水的連續性,因此充分了解在間歇輸水情況下河道地下水動態變化規律是輸水生態效益評價的關鍵[4]。

地下水水位的估算能夠很好地描述地下水流的動態變化,對水文地質研究、地下水管理等具有重要作用,常常有必要使用確定性模型重建地下水變化過程[5- 6]。塔里木河下游河道作為典型的間歇性生態輸水河道,在輸水-斷流交替出現的過程中,其地下水滲流場運動過程中各個運動要素(水位,流速,流向等)隨時間產生一定變化,具有非穩定流特性,可根據非穩定流理論進行求解。自1935年C. V. Theis[7]提出了地下水非穩定流解析解以來,國內外學者基于非穩定流原理對地下水變化模型構建的研究逐漸增多,王玉林等[8]運用非穩定流原理對抽-灌同軸非完整井承壓層地下水變化過程進行了模型構建;陳志強等[9]基于非穩定流原理對巖體滲透系數進行了求解;Liggett等[10]基于邊界積分方程對承壓含水層非穩定流理論進行了研究;Koch等[11]基于非穩定流理論構建了耦合地表水、地下水運動關系的冰川融水徑流模型;另外,Jha等[12]、Merkhali等[13]將非穩定流理論運用到河道泥沙運動過程中,得到了較好的擬合結果。

對于塔里木河下游河道,基于非穩定流理論,前人也進行了地下水變化模型構建研究[14- 16],并取得了一定進展。然而,以往的地下水模型主要適用于河岸附近,忽略了地下水運動過程中河岸距離的增加引起的地下水變化的時間滯后效應,另外忽略了潛水蒸發的影響,使得地下水變化擬合精度隨河岸距離增加而逐漸降低。本研究基于前人研究成果,以塔里木河下游英蘇斷面地下水動態變化為研究內容,基于非穩定流理論,以水位邊界條件作為初始求解條件,以余誤差函數erfc(x)為求解函數,綜合考慮地下水變化滯后效應及潛水蒸發作用,構建了間歇性輸水河道地下水模型,并運用英蘇斷面C3、C4、C5、C6、C7監測井2011—2015年河岸1050 m范圍內逐月地下水埋深監測數據進行了定量分析,以期構建適用于間歇性輸水河道附近地下水變化過程模型,合理估算地下水變化趨勢,一方面為指導塔里木河下游輸水方案及生態效益評價提供科學依據,另一方面為間歇性輸水河道地下水位變化研究提供模型參考。

1 地下水非穩定流計算方法

1.1 流量邊界條件已知的情形

河道充水后,若河道滲漏量隨時間變化值已知,則輸水引起的兩岸地下水位的上升值,可根據地下水動力學原理[17]進行求解。

(1)

式中,tj為時間軸刻度,其中j(j= 1,2,3,…,k)為間歇輸水次數,t0= 0;S(x,tj)為離河渠x距離的點在tj時刻地下水位變化值;qj為tj時刻河道滲漏量;ierfc(λ)為補余誤差函數;a為含水層壓力傳導系數;k為含水層滲透系數;?為含水層平均厚度;S(x,tj)為離河渠x距離的點在tj時刻地下水位變化值(自初始水位算起),S(x,t0) = 0。

1.2 水位邊界條件已知的情形

河道充水后,如果河岸x處水位變化值已知,通常把水位變化曲線按時間軸剖分成n(n≥ 0)個小段,相鄰之間的變化可看作是水位定值時的狀況,從而把整個過程理解為各個微小段的疊加,利用特定水位邊界時的水位計算方法[17- 18],就可以確定充水時所引起的兩岸地下水位的上升值S(x,tj),即是

(2)

式中,Hj為tj時刻河岸x(x=0)處地下水位;erfc(λ)為余誤差函數。

2 間斷輸水河道地下水交互式模型構建

影響地下水數值模擬結果不確定性的因素很多,如何從眾多的因素中篩選出影響模型的主控參數是模型分析的關鍵[19]。基于生態組成部分相互影響,地下水變化需要綜合考慮水的補給(生態輸水、自補給)和消耗(地下水流動)之間的平衡[20]。其中,河道水位和流量的變化是影響地下水位動態的重要因素,一般情況下,河流水位高于附近地下水位時,將補給地下水;河流水位低于地下水位時,河渠就成為地下水的排泄出路。塔里木河下游作為典型的間歇輸水河道,其地下水位變化包括兩個階段：一、生態輸水階段(圖1)。此時河道滲漏主要經歷以下兩個過程,當地下水位低于河道底部高程時,河道產生自由滲漏;當地下水峰上升至河道底部后,地下水與河水中的地表水連成一體,也即是地下水水位邊界條件已知,在這種情況下,河道將產生頂托滲漏,由于其輸水特性,河道滲漏過程為自由滲漏與頂托滲漏方式交互進行。二、河道斷流階段(圖1)。此時地下水缺乏持續有效的地表徑流補給,地下水變化主要經歷以下兩個過程,地下水位側向自補給過程,以及潛水蒸發作用引起的地下水位降低的過程。

圖1 英蘇斷面輸水河道生態輸水-斷流交替階段地下水變化示意圖Fig.1 Ground water level changes diagram in ecological water conveyance-cutoff alternating period of Yingsu-sectionH: 地下水位;D: 地下水埋深;C3—C7: 監測井;x: 河岸距離;s: 水位變化值;S: 坐標軸; Eg: 潛水蒸發;h0: 初時水位

然而,在塔里木河下游已有的間歇性輸水過程中,河道生態輸水流量年際變化較大,導致不同輸水階段流量邊界分布不穩定,使得實際測量過程中,生態輸水流量邊界的確定比較困難,工作量較大,不易求解。然而,河道附近設置有監測井,可較好地反映河道附近水位變化,易于界定水位邊界。整個輸水及斷流過程中,在橫向上僅存在側向自補給,即高水位區向低水位區形成的側向流,可根據公式(2)對英蘇斷面地下水位進行模擬,然而,隨著河道距離x的增加,余誤差函數erfc(x)計算結果將會逐漸減小,并且小于實際減小幅度,導致計算誤差較大。因而,為提高擬合精度,需綜合考慮河道距離所產生的水位變化時間誤差,則有

(3)

式中,tl為距離引起的滯后時間,x為河岸與河道間的距離,u為地下水的實際流動速度。

根據達西定律,實際流動速度u=v/n=ki/n,其中,v為地下水滲流流速,n為含水介質的孔隙度,k為含水層滲透系數,i為地下水運動的水力梯度,i=Ha/xa,xa為非穩定流相鄰兩個滲流階段的河岸距離差,Ha為xa距離內的水位差。然而,達西實驗是基于恒定均勻滲流發生的,對于英蘇監測斷面非穩定地下流,水力坡度i隨河岸距離呈非穩定變化,此時u已不再是斷面平均流速,而是滲流斷面中任一點的流速[21]。根據研究斷面地下水非穩定流運動規律,將河岸帶地下水滲流過程分為m個部分,此時地下水實際流速為

(4)

那么,式(3)中不同河岸距離所引起的時間滯后差則可表示為

(5)

在此情況下,地下水位變化計算公式為

(6)

式中,(Hk-H1)為模擬期間水位變化總值。

為了便于測算,將地下水水位(H)變化轉換為地下水埋深(D)變化,則計算公式可表示為

(7)

式中,D為地下水埋深值,(DK-D1)為模擬期間埋深變化總值;其中S計算結果為負值表征地下水埋深減小,水位升高,正值表征地下水埋深增加,水位下降。

對于豎直方向上,潛水蒸發作用在一定程度上影響地下水埋深變化[22],此時,由于蒸發作用所引起的地下水位的變化量T為：

(8)

式中,Eg為實際蒸發量,b、c為蒸發相關系數,由實驗資料確定,E0為水面蒸發強度,D0為初始地下水埋深,其中j(j≥ 1)為天數,T(0)=0。

綜合公式(7)、(8),最終地下水埋深模擬值為

D(x,tj)=D0+S(x,tj)+T(x,tj)

(9)

式中,D0為初始埋深值,S與T分別為非穩定流、潛水蒸發作用引起的地下水埋深波動值。

3 研究區概況

研究區英蘇監測斷面(C),距離大西海子水庫約61 km,位于若羌縣鐵干里克鄉境內的其文闊爾河上,地處北緯40°25′52.3″,東經87°56′18.7″,是下游輸水監測的核心位置(圖2)。區域內干旱少雨而蒸發強烈,大氣降水對地下水幾乎不產生補給,地下水的主要補給來源是區域內的河道。已有研究表明[23],生態輸水期間,塔里木河下游生態輸水過程對單側河岸影響寬度超過1000 m,并且隨著輸水量的增加,其影響寬度呈現增加的趨勢,而河道斷流期間,由于河道地下水自補給,地下水埋深波動范圍也超過1000 m。因此,為全面掌握英蘇監測斷面處地下水埋深的動態變化情況,分別利用距離生態輸水河道59、300、500、750 m及1050 m范圍內的C3、C4、C5、C6、C7監測井,對河岸約1050 m范圍的地下水埋深進行動態變化監測。

圖2 英蘇斷面地理位置圖Fig.2 Yingsu section location map

作為間歇性輸水河道,塔里木河下游自2011年1月初至2015年12月底,共進行了五個周期共計八次生態輸水工作(圖3),其他時間階段均為斷流階段。基于河道內地下水埋深在生態輸水期間易達到頂托滲漏,不適宜作為水位邊界求解條件,因而,本研究選取河岸59 m處C3監測井作為水位邊界逐月埋深變化值初始求解條件,如圖3,為C3監測井2011—2015年地下水埋深值變化圖。其中,C3、C4、C5、C6、C7監測井地下水埋深初始值D0分別為3.47,4.80,5.04,5.91,7.48 m;英蘇監測斷面地處塔克拉瑪干沙漠及庫魯克沙漠的交界地帶,其巖性主要為粉砂、細砂、粉土及粉質粘土的互層,河岸附近含水介質孔隙度n約為35 %,地下水含水層滲透系數k在1—5 m/d之間;地下水壓力傳導系數a為575.62 m2/d;i為階段性水力梯度,大小由不同河岸距離滲流范圍及該滲流范圍內非穩定流滲流變化規律決定;另外,水面蒸發強度E0為0.0075 m/d,蒸發相關系數b為0.0195,c為0.0011。

圖3 2011—2015年間歇性輸水過程圖及河道附近地下水埋深值變化圖Fig.3 Maps of intermittent water delivery process and groundwater depth change near the river during 2011—2015

4 結果與分析

4.1 地下水埋深值變化滯后期分析

由于地下水補給量的差異,隨著間歇性生態輸水過程的持續,塔里木河下游河道地下水位的總趨勢是上升的,但離河道不同距離處的地下水位上升幅度有所不同,且波動規律具有一定的差異。研究表明,在河道輸水期間,由于含水層滲透系數的限制,在短期內,臨近河道地下水埋深波動較大,然而遠離河道地下水埋深波動并未受到較大的影響,生態輸水作用影響河岸附近地下水埋深變化會產生一定時段的滯后期,且隨著距離河道的增加,滯后期增長。本研究所涉及的C3、C4、C5、C6、C7監測井分布在距河道1050 m范圍內,處于河道輸水影響緩沖區內,在生態輸水期間,由于滯后效應,地下水埋深變化在短期內并不能對生態輸水產生響應,因此在運用非穩定流理論對地下水模擬過程中,除了考慮河道附近水位邊界變化之外,還需綜合考慮滯后效應所引起的地下水埋深滯后變化。

通過河道距離x及含水層滲透系數k,計算出不同監測井地下水埋深變化滯后期,以2011—2015年生態輸水期為時間起點,綜合滯后期,得到地下水埋深變化滯后統計結果(表1)。由表1可知,在不同年份生態輸水期間,C3、C4、C5、C6、C7監測井由于河岸距離的不同,地下水埋深變化滯后期達到1—13月之間,其中,2015年輸水期間,C5、C6、C7監測井滯后期過長,超出計算年限,因此忽略了此段時期滯后效應計算。

表1 2011—2015年不同監測井輸水前后地下水埋深變化滯后期統計結果

4.2 地下水埋深值擬合結果分析

經過初始地下水埋深和水位邊界條件設定、水文地質參數賦值、不同河岸距離地下水埋深變化時間滯后效應界定、潛水蒸發條件設定、模型校驗、參數微調、公式計算,最終完成了英蘇監測斷面地下水埋深變化模型構建。以C4、C5、C6、C7監測井初始埋深值、2011—2015年河道附近C3監測井地下水埋深變化值為初始求解條件,運用公式(5)、(6)非穩定流計算模型,潛水蒸發計算模型進行模擬計算,得到C4、C5、C6、C7監測井2011—2015年逐月埋深變化及潛水蒸發變化模擬結果(圖4)。

圖4 2011—2015年非穩定流及蒸發作用引起的地下水埋深值變化模擬值Fig.4 Simulation of groungwater depth changes caused by unstedy flow and evaporation during 2011—2015

由圖4可知,生態輸水作用對地下水埋深月變化及年變化影響較大,引起的地下水埋深變化值較大,并且波動明顯;基于綜合考慮生態輸水期間地下水埋深變化滯后效應,不同河岸距離不同監測井地下水埋深變化波動趨勢具有一定的差別,但總體變化趨勢一致。潛水蒸發作用在短期內對地下水埋深變化影響作用較小,然而在多年時間序列上對地下水埋深變化總影響值不可忽視,在2011—2015年期間,潛水蒸發作用引起河道附近地下水埋深變化累加平均值超過0.1 m,在一定程度上對地下水埋深產生了影響。其中距離河道1050 m處C7監測井潛水蒸發作用變化累加值為0.066 m,距離河道750 m處C6監測井潛水蒸發作用變化累加值為0.111 m,距離河道500 m處C5監測井潛水蒸發作用變化累加值為0.118 m,距離河道300 m處C4監測井潛水蒸發作用變化累加值為0.122 m,表明隨著河岸距離的減小,潛水蒸發作用對地下水埋深的影響逐漸增強,且隨著時間的累積,以及生態輸水引起的河道附近生態恢復效應的突顯,該影響作用將持續增強。

運用公式(7),將C4、C5、C6、C7監測井地下水埋深初始值與地下水埋深變化模擬值(圖4)累加計算,得到2011—2015年地下水埋深最終模擬結果。綜合英蘇斷面不同監測井2011—2015年逐月地下水埋深實測數據,得到地下水埋深變化擬合-實測對照圖(圖5)。由圖5可知,對于各監測井埋深變化實測值,在時間序列上波動變化明顯,在間歇性輸水過程中地下水埋深值呈現階段性增減,總體表現為生態輸水階段地下水埋深減小,斷流階段地下水埋深增加。對于各監測井埋深變化模擬值,總體變化趨勢與實測值一致,在河道斷流期間,由于地下水自補給及潛水蒸發作用,模擬得出的地下水埋深值逐漸增加;在生態輸水期間,由于河道內地表徑流對地下水補給作用,河道附近地下水埋深值總體呈現減小的趨勢。綜合考慮表1滯后效應所對應的滯后期進行埋深變化求解,有效的獲取了生態輸水期間不同監測井地下水埋深滯后變化規律,得到較好的擬合變化結果,地下水埋深模擬值在輸水月份及其對應的滯后月份均符合實測變化趨勢。

圖5 2011—2015年地下水埋深值擬合結果Fig.5 Simulation of groungwater depth during 2011—2015

4.3 地下水埋深擬合精度分析

為研究地下水擬合效果,對模擬結果進行精度分析。圖6為地下水埋深觀測值與模擬值檢驗結果圖,可以反映模擬值偏離實測值的程度。由圖可知,數據主要集中在擬合曲線周圍,沒有明顯偏離的情況,表明塔里木河下游斷面尺度地下水埋深觀測值與模擬值表現出較好的相關性。

圖6 2011—2015年地下水埋深觀測值與模擬值有效性檢驗Fig.6 Validity test of observation value and simulation value of groundwater depth during 2011—2015

表2為地下水埋深擬合精度統計結果,其中確定系數R、相關系數R2反映擬合精度,平均絕對誤差(mean absolute error,MAE)表征模擬值偏離實測值的大小,平均相對誤差(mean relative error,MRE)為平均絕對誤差與平均模擬值的比值,用以表征擬合結果的準確性水平[24- 25]。對于C4、C5、C6、C7監測井地下水埋深模擬值與實測值相關關系,在P<0.01置信區間內均達到顯著性水平;對于不同監測井其確定系數R值均大于0.791、相關系數R2均大于0.625;另外,對于平均相對誤差MRE,均不超過13%。結果表明,基于地下水非穩定流理論及潛水蒸發模型擬合結果的相關性較好,表明以河道為基點的地下水非穩定流運動及潛水蒸發作用是影響沿岸地下水變化的諸多環境因子中最敏感的因子。

表2 C4—C7監測井地下水埋深值擬合精度統計結果

*表示在P<0.01水上平顯著相關

另外,隨著距離河道的增加,確定系數R值、相關系數R2逐漸減小,地下水模擬精度逐漸降低,其中,C5、C7監測井擬合精度較C4、C6監測井低,主要原因是隨著距離河道的增加,影響地下水埋深變化的隨機性因素逐漸增多。其中C5監測井分布于距河道約500 m處,為植物分布比較集中的位置,植物生長及蒸散發作用對地下水埋深變化產生一定的影響;C7監測井分布于距河道約1050 m處,距離河道較遠,受到多元方向的地下水補給與排泄作用,影響模擬精度。但是總體上,對于河道附近1050 m范圍內地下水埋深模擬,已達到了所需精度要求,對于間歇性輸水河道地下水埋深變化狀況模擬具有重要意義。

5 結論與討論

本研究采用非穩定流理論,以水位邊界條件為自變量,并綜合考慮潛水蒸發及生態輸水期間地下水埋深變化在時間及距離尺度上的滯后效應,對塔里木河下游斷面尺度單側河岸1050 m范圍內地下水變化進行模擬,得出以下結論：

(1) 生態輸水期間,在短期內,臨近河道地下水埋深波動較大,而遠離河道地下水埋深波動并未受到較大的影響,生態輸水作用影響河岸附近地下水埋深變化會產生一定時段的滯后期,且隨著距離河道的增加,滯后期增長,其中C4監測井滯后期最短,C7監測井滯后期最長。

(2) 非穩定流作用對地下水埋深月變化及年變化過程中影響較大,引起的地下水埋深變化值較大,并且波動明顯;潛水蒸發作用在短期內對地下水埋深變化影響作用較小,然而在多年時間序列上對地下水埋深變化總影響值不可忽視。綜合分析表明,以河道為基點的地下水非穩定流運動及潛水蒸發作用是影響沿岸地下水變化的諸多環境因子中最敏感的因子。

(3) 通過本研究所應用的非穩定流理論及潛水蒸發模型,在河岸300、500、750、1050 m處擬合精度R2分別達到0.701、0.654、0.701、0.625,達到了地下水模擬精度的要求,對塔里木河下游斷面尺度上地下水恢復狀況研究具有重要的理論意義與參考價值。

輸水效益的顯現是一個漫長的過程,地下水的響應和下游植被的生態響應在一個大的空間和時間尺度上將逐步顯現[26]。英蘇斷面作為塔里木河下游核心監測帶,自河道處開始至1050 m范圍內,依次分布有胡楊-檉柳-草甸混合帶、草甸帶、草甸-荒漠混合帶、半荒漠帶,且不同植被類型所需地下水埋深值有一定的差異。其中,當地下水埋深3—3.5 m時,胡楊-檉柳林分布最大[27];當地下水埋深在3—5 m時,最適合草本植物生長;當地下水位低于6 m時,物種的多樣性均勻度豐富度迅速降低[28]。運用非穩定流理論對地下水位變化模擬,在已知河道附近地下水埋深變化值、間歇性輸水時間以及初始地下水埋深情況下,可對河道附近地下水影響范圍內任意位置地下水埋深變化進行模擬,結合不同植被類型生態響應期,進而對生態效益進行評估,在一定程度上可提高生態評價的效率。

然而,在實際計算中,由于觀測不準或極端氣候因素、多元地下水補給及排泄等外因的緣故,地下水模擬結果與實際測量結果存在一定的偏差[29]。本研究側重測算間斷歇性輸水過程及潛水蒸發作用對地下水位變化影響,忽略了其他因素對地下水位的影響,同時本研究是基于斷面尺度研究,在一定程度上限定了區域內地下水補給及排泄的邊界。因此在今后的研究中需要進一步細化模型,增加更多的影響變量以提高模型的精度,同時伴隨今后監測精度的提高,獲得相關參數的精確性會更高。另一方面,下一步研究過程中,將在以地下水模擬為主的生態水文研究的基礎上,增加地下水-植物耦合關系,以更好的建立生態輸水-地下水-植物生態效益快速評價體系。