毛烏素沙地南緣降水入滲滯后補給與模型參數敏感性分析

王 鍇, 楊澤元, 袁 悅, 陳志軍

(1.長安大學 環境科學與工程學院, 陜西 西安 710054; 2.旱區地下水文與生態效應教育部重點實驗室,陜西 西安 710054; 3.陜西省地下水與生態環境工程研究中心, 陜西 西安 710054)

毛烏素沙地位于半干旱氣候區，境內包含陜北和內蒙能源化工基地，隨著該地區工農業發展與人口增長，用水量短缺形勢越發嚴峻。該地區主要用水來源為地下水，降水入滲是地下水補給的主要來源，研究該地區降水入滲補給是其地下水資源評價的重要基礎工作[1]。前人[2-3]常采用年降水量與區域入滲補給系數計算區域性降水入滲補給量，很少考慮入滲過程的滯后性。地下水補給總是滯后于降水事件，某一時段內補給量實際由該時段之前的時段內降水所決定。同時，由于補給滯后的存在，含水層對于人類活動與氣候變化的響應往往不會立即發生[4]，這給地下水可持續開發帶來挑戰。因此研究地下水的補給滯后現象對于提升地下水資源評價精度與地下水資源合理利用具有重要意義。

1 材料與方法

1.1 研究區概況

研究區位于陜西省榆林市補浪河鄉，地處毛烏素沙地南緣。年均降水量365 mm，年均蒸發量1 246 mm。主要潛水含水層為第四系風積孔隙含水層，巖性為中砂，不含黏土(表1)。地下水常年埋深1.3～1.5 m，主要補給來源為大氣降水，占天然補給總量85%以上[16]，主要排泄去路為蒸發排泄。

表1 研究區不同層土壤質地

1.2 研究方法

在原國土資源部地下水與生態陜西榆林野外觀測基地建立原位試驗場。監測地下5，10，30，60，90 cm深處含水率，同時獲取場地氣象數據。包氣帶含水率和溫度由5 TM傳感器監測(美國Decagon公司)，采集頻率為10 min。氣象數據由試驗場內的波紋比系統(05130-5RMYong，美國Campbell公司)自動監測。地下水位由MiniDiver自動監測并記錄數據，氣象數據和地下水位的采集頻率均為1 h。

以野外原位監測數據為基礎研究場地內包氣帶含水率與地下水位對降水入滲的滯后響應。根據場地條件建立水文地質概念模型并利用Hydrus-1D軟件建立降水入滲模型。采用野外監測數據進行模型參數的識別與優化。

采用相對平均誤差(AVRE)、相關系數(R)作為指標進行識別期模型誤差分析。

AVRE=∑(VC1/VO1)/n

(1)

(2)

式中：VCi——模擬值； VOi——原位觀測值；n——數據總數量。

通過上述過程識別獲取基本參數，基于原位監測含水率，負壓，溫度數據建立反解模型(共14 d，4 695個數據點)。反解計算采用Hydrus-1D的Inverse solution模塊進行，對van-Genuchten方程中的各參數θr(殘余含水率)，α,n(水土特征曲線方程的參數)和Ks(飽和滲透系數)進行反解，反解方法為最小目標函數法[17]。選取SSR，r2最優對應的參數作為參數優化結果。其中SSR為反解過程中目標函數值。r2為擬合值與觀測值的相關系數，r2為1表示擬合最佳。

(3)

注意到,熱帶西北太平洋位于東亞季風區,而熱帶東南印度洋位于南亞季風區以南(Wang and Linho,2002),環流結構有所不同。同時,在氣候學上夏季位于對流層低層起自索馬里至西太平洋熱帶地區的西風氣流和對流層上層的東風氣流將西太平洋熱帶地區與熱帶印度洋地區聯系起來。對印度洋與太平洋上環流變化聯系的研究僅限于赤道或近赤道的緯向區域,其聯系的重要紐帶是Walker環流,而對于赤道外地區兩大洋上空環流變化的聯系卻未能有較好的研究。那么,位于海洋性大陸東南邊緣和西北邊緣的兩個地區上的環流變化是否存在聯系呢?這種聯系是否存在年代際變化呢?其原因是什么呢?本文將回答這些問題。

參數獲取完畢后，以累計底部通量為指標，采用單參數擾動方法進行敏感性分析，重點關注不同參數改變對于地下水補給滯后的影響。采用輸入輸出變化率(ROV)來分析模型計算對于參數的敏感性[15]。

(4)

式中：C(t)——模型實際計算結果(改變某一參數值的結果)；Cref(t)——模型參考計算結果(所有參數保持不變的結果)，此處計算結果均為累計底部通量；P——實際輸入參數；Pref為參考輸入參數；t——時間。ROV(t)的絕對值越大，則模型對該參數越敏感。由于介質被劃為多層，因此輸入輸出變化率均采用各層平均值。

2 結果與分析

2.1 包氣帶含水率響應滯后

研究區降水年際變化大，降水主要集中在每年的6—9月份，為豐雨期。由于每年的11月至次年3月是該區的冰凍期，只發生少量降雪，研究時段確定為2014年3—10月并將時段內的降水依據國標劃分為小雨、中雨、大雨與暴雨[18]。包氣帶深度越淺，降水后含水率響應越劇烈，部分降水事件中地表最大含水率甚至接近飽和；同一地下水埋深下，降水量越大，包氣帶響應深度越大。以8月26日至9月5日為例(圖1)，8月27日降水0.82 cm(小雨)，僅5 cm深度有明顯響應，其他深度含水率均無明顯變化；8月30日降水7.1 cm(暴雨)，包氣帶所有深度含水率都明顯上升。選取研究時段內主要降水事件，通過相關性分析研究降水量與入滲響應深度的關系。如圖1所示，入滲響應深度與降水量呈顯著線性相關，相關系數R為0.94>R(16，0.01)=0.59。隨著包氣帶深度增加，含水率響應存在明顯滯后；以7月8日降水為例，10 cm處滯后時間為13 h，30 cm為19 h，90 cm為29 h。

圖1 入滲響應深度與降水量相關關系

注意到7月3日降水(降水量2.12 cm)歷時4 h，最大雨強14 mm/h，5，10，30 cm深度含水率均明顯響應，但90 cm處含水率只微弱變化；7月8日(降水量4.1 cm)歷時27 h，最大雨強僅3 mm/h，90 cm處含水率明顯上升。前者降水強度更大，后者降水量更大，這說明降水量對于響應深度的影響大于降水強度(圖2)。

圖2 降水量與含水率變化時間序列

2.2 地下水滯后補給

無降水的時段，地下水位保持相對穩定(圖3)，小雨和中雨后地下水位均未響應，如6月25日降水5.2 mm，水位并無變化。大雨暴雨均會對地下水形成有效補給(如6月29日、7月3日、7月9日的降水)，使得地下水位上升。試驗時段內對地下水形成有效補給的降水統計結果詳見表2。從表2可知，降水量越大，地下水位上升幅度越大。相近降水事件如7月3日降水(降水量2.12 cm)與6月29日降水(共2.24 cm)下，水位埋深與地下水位升幅反相關。這是由于在蒸發條件相似時，水位埋深越大，包氣帶水分虧缺越大，降水后補給量也越小。降水之后，地下水位并不立即上升，一般把從降水開始到地下水位上升的時間間隔稱為補給滯后時間；由表2可知，6月29日降水(共2.24 cm)補給滯后時間為4 h，7月3日降水(共2.12 cm)補給滯后時間為5 h；注意到這一滯后時間甚至小于30 cm處含水率響應時間，表明降水后存在優勢流通道，降水入滲通過優勢流通道提前進入地下水，引起地下水位響應。歷次降水后地下水補給滯后時間均隨著降水強度的增大而減小，但減小幅度逐漸放緩。這是由于隨著降水強度的增大，入滲率也隨之增大；但當降水強度逐漸接近或超過土壤入滲能力后，入滲率不再變化，相應的入滲補給滯后時間也趨于一定值。需要注意的是，9月10日降水3.06 cm，較9月22日降水5.46 cm差別并不大，但后者地下水位上升劇烈，這是由于后者的前期降水量更大，前期降水有效補給了包氣帶水分，較高的包氣帶含水率促進了后期降水入滲。

圖3 地下水位與降水量變化時間序列

雨 型大 雨暴 雨降水量/mm21.241.430.671.654.622.4降水強度/(mm·h-1)4.241.531.717.92.14.48地下水位升幅/cm1.31.62.99.19.91.8補給滯后時間/h5118364 地下水埋深/cm雨前129.4129.5130.6145.3133.1121.6雨后128.1127.9127.7136.2123.2119.8

2.3 降水入滲數值模型建立與參數敏感性分析

2.3.1 模型建立 以地表為坐標原點，向上為正。建立水文地質概念模型，地表作為模型的上邊界。地下埋深150 cm處作為模型的下邊界，忽略水平徑流，由此可以概化為變飽和帶一維非穩定流模型。含水介質依據顆粒分析和容重概化為3層，其中0—20 cm為第1層，20—80 cm為第2層，80—150 cm為第3層。在Hydrus中選用大氣邊界作為水分運移的上邊界條件，可變壓力水頭作為水分運移的下邊界條件。熱運移的上下邊界條件分別為上下邊界溫度。選用Van-Genuchten方程作為水力特征曲線，考慮持水曲線的滯后效應。含水率初始條件采用與初始時間地下水位相平衡的含水率分布。溫度初始條件設置為初始時刻監測值。模型的初始土壤水力學參數詳見表3。采用Hydrus中基于顆粒分析和容重數據的神經網絡預測結果。

表3 經神經元預測分析得到的土壤水力學參數

注：θs，α，n,，Ks，θr均為van-Genuchten方程中參數。θr為殘余含水率；α和n為水土特征曲線方程的參數；Ks為飽和滲透系數；θs為飽和含水率。

2.3.2 模型參數獲取 選取2014年7月2—15日作為模型參數的識別期。識別期包氣帶深度5，30，60 cm處含水率、溫度模擬值與觀測值對比如圖4所示。

圖4 模型識別期包氣帶不同深處含水率與溫度

模型在識別期的計算誤差結果詳見表4—5，除埋深60 cm處誤差較大外，總體計算結果較好。經過參數反演優化，選取SSR=1.854，r2=0.926時對應的參數作為優化結果，最終確定的水、熱力學參數詳見表6—7。

表4 識別期模型計算含水率誤差結果

注： AVRE為模型計算值與監測值的相對平均誤差，R為相關系數。

表5 識別期模型計算溫度誤差結果

2.3.3 van-Genuchten方程參數敏感性分析 由于調參過程中涉及的參數較多，對van-Genuchten主要的土壤水力學參數θr，α，n，Ks進行了敏感性分析。通過對某一參數進行增減，維持其余參數不變，確定參數敏感性[15]。敏感性分析中參數取用情況詳見表8。

表6 模型識別及優化后所得土壤水力學參數

表7 模型識別及優化后所得土壤熱力學參數

注：b1,b2,b3為熱傳導系數λ(θ)的參數；Cn為固態物質的體積熱容；Cw為水的體積熱容；Di為熱擴散系數；So為固態物質體積百分比。

表8 參數敏感性分析中所用土壤水力學參數

注：考慮到n最小取1，故n的變化為在原始值與1作差的基礎上再增減25%。

采用累積底部通量誤差進行參數敏感性分析。時間起點為7月2日0：00，圖中下標為0的表示該參數的初始值。-25表示減小25%的參數值，+25表示增加25%的參數值(圖5)。各參數的輸入輸出變化率詳見表9。從表9可知，底部通量對θs最為敏感，對Ks最不敏感。α與n的敏感性介于θs與Ks之間，其中n比α更敏感。若以底部通量代替地下水補給量，則地下水補給量的大小主要由θs控制，θs越大，入滲過程中包氣帶持水能力越強，地下水補給量越小。從圖5中可以看出，θs與Ks對于降水補給的滯后影響較大，當Ks減小或θs增大時，累計底部通量拐點時間(指示著補給開始時間)明顯滯后；其中θs增大25%，滯后時間增長10～16 h；Ks減小25%，滯后時間增長6～8 h。值得注意的是，若包氣帶剖面處于蒸發狀態，Ks的改變幾乎不引起底部通量的變化，此時Ks的敏感性很弱。

圖5 不同參數變化對模擬底部通量的影響

項目α-25n?-25Ks-25θs-25α+25n?+25Ks+25θs+25ROV(tend)1.222.930.442.21.783.760.254.44ROV(tmax)9341550241517136018401322630ROV(tmin)1.212.880.310.00751.782.110.00134.42

注：tend表示結束時刻，tmax為ROV(t)最大時刻，tmin為ROV(t)最小時刻； 下標為-25表示在原參數基礎上減小25%的參數值，+25表示增加25%的參數值。

3 討 論

3.1 降水入滲滯后補給的影響因素探討

由于土地利用變化、氣候變化和工農業發展，世界不同地區的土壤水分和地下水正越來越多地受到影響但呈現出不同規律。本文確定的地下水滯后補給時間為4～11 h，小于相近水位埋深下三江平原地區滯后時間(28～50 d)[19]，這是由于后者包氣帶巖性更細，入滲速率較慢。在大埋深地區，由于入滲補給路徑很長，地下水補給滯后時間甚至能達數十至數百年[4，20]。土壤水分分布同樣是影響入滲補給的重要因素，按照活塞流入滲理論，包氣帶含水率響應滯后時間與土壤水分分布的關系可以表示為[21]：

I(t)=(θ0-θi)Zf

(5)

式中：I(t)——累計入滲量；Zf——入滲深度；θ0——入滲面上部含水率；θi——剖面初始含水率。

由公式(5)得知，包氣帶含水率響應滯后時間隨包氣帶初始含水率增大而減小，而初始含水率反映了前期降水的影響,因此入滲補給除受到本次降水的影響，還受到前期降水的控制。

3.2 多參數敏感性分析

本研究采用的是單參數敏感性計算，能夠較為精確得出各個參數的敏感性。但實際的土壤參數之間存在相關性且總是處于一個合理范圍。比如θs增大導致降水入滲滯后時間增大，但是θs增大同樣會導致Ks增大，后者能促進降水入滲；因此，對于Ks與θs等相關性較強的參數采用多因素分析能更為準確地探討其敏感性。

3.3 未來研究

由于研究區的地下水埋深淺，降水入滲補給速率相對較快，包氣帶淺部補給滯后不明顯，通過監測數據確定補給滯后時間的精度并不高。此外，研究時段內降水頻次高。因此后續研究應當考慮累計降水對于降水入滲的影響并提高淺部包氣帶監測精度。研究區處在固定沙丘上，地表平緩，側向徑流微弱，故采用一維滲流模型；但場地內稀疏分布著沙柳，干旱時期沙柳吸收水分較強，對土壤水分分布存在影響[22]。同時地下水補給滯后時間小于包氣帶含水率響應滯后時間，這表明場地中存在不可忽視的優勢流通道，可能的原因是場地內包氣帶巖性不均一或植被根系吸水[23]。因此后續研究也應當針對植被吸水做進一步分析并考慮使用雙滲透模型解決優勢流問題。

4 結 論

結合原位監測數據，分析了毛烏素沙地地下水淺埋區降水入滲補給規律，利用Hydrus-1D軟件建立了降水入滲數值模型，獲取了模型水力學參數和熱力學參數；分析了van-Genuchten方程主要水力學參數敏感性；為研究包氣帶和地下水位對降水入滲的響應機制、分析降水入滲補給機理奠定了基礎。

(1) 包氣帶含水率響應深度隨著降水量增大而增大，其中小雨型響應深度為3～10 cm，中雨型為30～60 cm，大雨型為60～90 cm，暴雨型均大于90 cm。入滲響應深度與降水量的關系可表示為：D=35.5P-17.1，其中D為入滲響應深度(cm)，P為降水量大小。隨著前期累計降水量增大，降水入滲深度也增大。

(2) 研究區內地下水補給滯后時間約為4～11 h。包氣帶水分和地下水滯后補給受地下水埋深、包氣帶巖性和含水率分布共同影響，其中包氣帶水分滯后補給時間隨著初始含水率增大而減小。

(3) 建立了該地變飽和帶水汽熱耦合運移模型。采用模型識別與驗證獲取了包氣帶水力學參數與熱力學參數，以此作為參數反演優化的初始輸入值；參數優化最佳結果θr為0.03，θs為0.34，Ks為23 cm/h,各層n為2.0～2.2，α為0.060～0.074。

(4) 與數值模型有關的各項水力學參數中，敏感性自強到弱的排序是θs，n，α，KS，底部通量對于θs最為敏感；降水補給滯后時長與θs呈正相關關系。