考慮土壤膨脹性的流域水文模型應用研究

劉海瀅，甘永德，賈仰文，仇亞琴，牛存穩

(1.中國水利水電科學研究院流域水循環模擬與調控國家重點實驗室，北京 100038；2.青海大學省部共建三江源生態與高原農牧業國家重點實驗室，西寧 810016)

目前，水資源短缺問題越來越突出，為了適應變化環境和人類活動影響下的水文水資源研究，考慮水文參數和水文過程的空間異質性的分布式水文模型成為研究流域水文循環的一種有效工具，被認為是探索和認識復雜水文循環過程和機理的有效手段。隨著遙感技術(RS)、地理信息系統(GIS)和全球定位系統(GPS)的緊密結合，分布式和半分布式流域水文模型不斷開發出來，例如USGS-MMS模型[1]、SWMM模型[2]、SHE/MIKESHE模型[3]、TOPMODEL模型[4]、OHyMoS模型[5]、SWAT模型[6]、WEP[7]等。

我國是膨脹性土壤分布最廣泛的國家之一，膨脹性土壤問題又不斷影響水文過程，現成為水文學研究的熱點。土壤膨脹是由于土壤膠體，尤其是黏粒部分的水化所引起的[8]。土壤水分運動受土壤吸水膨脹、失水收縮，濕脹干縮過程影響，土壤膨脹變形主要與初始含水量和上覆荷載等有關，膨脹力和膨脹變形隨土壤增濕程度增加而增加。在吸水膨脹過程中，土壤膨脹變形多反映為垂直向上。國內外眾多學者對膨脹性土壤入滲進行了相應研究。McGarry等研究了土壤含水量與土壤變形量間關系，并提出了用于描述土壤濕脹干縮變化特征的三直線模型[9]。蘇寧虎采用分數階偏微分方程建立了膨脹性土壤入滲方程[10]。甘永德采用土壤膨脹特征曲線和土壤應力-應變關系曲線建立了膨脹性土壤飽和水分運動參數計算模型[11]。然而，在進行流域水文過程模擬時對已有分布式水文模型還未考慮土壤膨脹性的影響。

由此，基于前期研究成果，對傳統降雨入滲產流模型TGAM(Traditional Green-Ampt mode)模型進行修正，引入膨脹性土壤導水系數和膨脹性土壤飽和含水量，提出考慮土壤膨脹性的降水入滲產流模型GJGAM(Gan and Jia propose a modified Green-Ampt model)。進一步地，將改進后的入滲產流模型應用到流域分布式水文模型WEP-L模型中，并選擇南小河溝流域對改進前后的WEP-L模型模擬效果進行對比分析，驗證模型改進后的改進效果及適用性。

1 分布式水文模型WEP-L模型改進

1.1 分布式水文模型WEP-L模型介紹

WEP-L(Water and Energy transfer Processes in Large river basins)模型是在WEP模型基礎上，舍棄柵格模式，采用子流域套等高帶作為最小計算單元，形成WEP-L模型[12]。2006-2010年，又同多目標決策分析模型(DAMOS)和水資源配置模型(ROWAS)耦合形成流域二元水循環模型，用于強人類活動影響下的流域水循環過程模擬。

1.2 考慮土壤膨脹性影響的WEP-L模型改進

(1) TGAM模型修正假設。目前的WEP-L模型中的降雨入滲計算模塊采用的是傳統沒考慮土壤膨脹性Green-Ampt模型，其未考慮土壤膨脹性變形引起的土壤導水系數和飽和含水量的變化，這種認識無法真實地模擬降雨入滲過程，進而影響WEP-L模型的模擬精度。為此，對Green-Ampt模型做如下修正。

取地面為參照面，向下為正(圖1)。膨脹性土壤吸水膨脹變形，土壤變形受土壤膨脹力和自重應力影響，導致土壤濕潤區剖面飽和導水系數和飽和含水量均隨深度變化而變化，為了便于分析，做如下假設：①變形前土壤為均質土壤；②土壤變形為彈性，即變形無滯后性；③土壤變形只引起土壤孔隙度變化；④入滲過程中存在明確的濕潤鋒面，濕潤鋒面將濕潤區和未濕潤區截然分開，濕潤區土壤達到飽和，未濕潤區土壤含水量為初始含水量。為了便于描述導水系數隨深度的變化特性，引入膨脹性土壤導水系數描述濕潤區導水系數；引入考慮膨脹性土壤飽和含水量描述濕潤鋒以上飽和含水量。與傳統沒考慮土壤膨脹性的Green-Ampt模型相比，模型考慮了土壤膨脹變形引起的土壤導水系數和飽和含水量的變化。

圖1 土壤入滲過程示意圖Fig.1 Schematic diagram of the infiltration in soil

土壤吸水膨脹變形是土壤膨脹力和自重應力共同作用的結果，假設土壤膨脹變形是由土壤孔隙度的變化引起的，則當土壤飽和時，土壤膨脹力引起的孔隙度變化量可以表示為：

(1)

式中：ρd為土粒密度，g/cm3；e0為土壤初始孔隙度，cm3/cm3；ew為由土壤吸水膨脹導致的孔隙度變化量，cm3/cm3；ρsw為由土壤吸水膨脹導致的容重變化量，可以采用土壤膨脹特征曲線計算(三直線模型結構段計算)：

(2)

式中：c和α3均為三直線模型參數；U為土壤質量含水量，g/g。

同理，土壤自重應力引起的孔隙度變化量可以表示為：

(3)

式中：ep為由土壤自重應力導致的孔隙度變化量，cm3/cm3；ρsp為由土壤自重應力導致的容重變化量，可以采用土壤應力----應變關系曲線計算：

ρsp=A+Bln (γZ)

(4)

式中：γ為土壤濕比重，N/cm3；Z為土壤深度；A和B均為參數。

土壤膨脹變形是膨脹力和自重應力共同作用的結果，則合力導致的土壤空隙變化量可以表示為：

de=dew+dep=Δew+Δep

(5)

當土壤飽和時，土壤孔隙被水分充滿，即土壤飽和含水量等于孔隙度，則土壤剖面飽和含水量總量θT可以表示為：

式中：θT為土壤深度Z以上區域的飽和含水量，cm3/cm3；Z為土壤深度，cm。

受土壤膨脹變形影響，土壤飽和導水系數隨深度的變化而變化。針對土壤飽和導水系數，采用改進的Lambe模型計算：

Ks(e)=K010m(ez-e0)

(7)

(8)

式中：Ks(e)為孔隙度為ez時土壤飽和導水系數，cm/min；ez為土壤深度為Z時土壤孔隙度；K0為孔隙度為e0時土壤飽和導水系數，cm/min；m為與土壤孔隙度性質有關的參數。

由于降雨過程中，雨強非恒定，因此，將降雨過程劃分為x個時段，每個時段內降雨強度恒定，x(x∈ 0,1,…,n)，n為時段數。同一降雨時段(tx-1～tx)內，降雨入滲特性由本時段的降雨強度I、時段初的積水深度h0和潛在入滲強度fpt決定。由于積水深度h0與雨強單位不同，模型計算時，將h0除以相應時段，轉換為雨強P′。根據時段內降雨強度、時段初積水深和潛在入滲強度，時段內入滲過程可以分為以下情景(圖2)。 ①h0=0,I>fpt≥Ks(e)，其中為膨脹性土壤濕潤區導水系數；I為雨強；h0為時段初積水深度；fpt為積水入滲率。這種情況下，隨著降雨的持續，地表開始積水，土壤入滲過程可以分為非積水入滲過程和積水入滲過程；②h0>0,P′+I0,P′+I≥fpt≥Ks(e)。 這種情況下，土壤繼續進行積水入滲過程；④h0=0,I

圖2 積水過程和非積水過程轉換情景示意圖Fig.2 Scenario diagram of the Infiltration process for ponding condition and non-ponding condition

根據達西定理有：

積水前：

fnpt=I

F=Fx-1+(t-tx-1)I

(9)

積水后：

(10)

忽略地表積水：

(11)

式中：fnpt為積水前入滲強度，cm/min；I為雨強，cm/min；fp為積水后土壤入滲率，cm/min；SW為濕潤鋒土壤水吸力，cm；H0為積水深度，cm；Z為濕潤鋒距離，cm。

由水量平衡原理，可以得出某一時刻t的累計入滲量F可以表示為：

F=θT-θ0Z=

令：

(13)

則：

(14)

式中：θT為濕潤鋒以上土壤飽和含水量，cm3/cm3；F為土壤累計入滲量，cm。

(15)

積分得：

(16)

A=SWΔθ

(17)

積水時刻確定：

(18)

(19)

式中：fnpt為非積水時段入滲強度,cm/min；Fp積水發生時刻土壤累計入滲量，cm；F為土壤累計入滲量，cm；tp為土壤表層積水發生時間，min；Ip為土壤表層積水發生時的時段降雨強度，cm/min；I為x時段降雨強度，cm/min；fpt為積水時段土壤入滲率，cm/min；t為時間，min；θ0為土壤初始含水量，cm3/cm3；θT為濕潤鋒以上土壤飽和含水量，cm3/cm3；SW為濕潤鋒平均吸力，cm；A為參數；Ks(e)為土壤飽和導水系數，cm/min。

(2) 模型參數。甘永德等[13]基于室內試驗分析，給出了膨脹性土壤飽和含水量和飽和導水系數，以及這些參數計算方法。通過模型應用，以及與室內試驗對比表明，該方法可以給出較好的模擬結果。

在實際應用中，這些參數值因土壤類型不同而不同，通過測定土壤膨脹特征曲線、水分特征曲線、土壤應力----應變關系曲線、土密度、土壤孔隙度與飽和導水系數間關系曲線、K0及其對應孔隙度等獲得相應參數。

2 模型應用

應用改進前后的WEP-L模型模擬南小河溝流域的降雨-徑流過程，檢驗模型改進效果。

2.1 研究流域概況

南小河溝(圖3)是涇河支流蒲河左岸的一條支溝，屬于黃土高原溝壑典型區，位于東經107°30′～107°37′，北緯35°41′～35°44′，流域面積36.3 km2。流域長約13.6 km，平均寬度為3.4 km，高程1 081～1 430 m。流域主要地貌單元為塬、坡和溝。塬面平緩，坡度一般5°以下；梁峁坡為破面與塬邊間的緩沖帶，坡度一般在10～20°之間；梁峁坡以下為溝谷，呈“V”字形，坡度一般25°以上。流域年平均降水量520.0 mm，降雨年內不均勻，多集中在7-9月。流域平均氣溫8.7 ℃，年積溫2 700～3 300 ℃，年日照2 454.1 h。流域土壤主要為黑壚土和黃綿土，其中黑壚土廣泛分布于塬面，黃綿土分布于塬面以下溝坡部位。植被類型以暖溫帶森林草原為主。

圖3 南小河溝流域概況圖Fig.3 Research of Nanxiaohegou river basin map

2.2 計算單元劃分

以數字高程模型和天然河網水系為依據，提取南小河溝流域的數字河網水系，并生成流域產匯流分區。采用全球SRTM90m分辨率的DEM數據，劃分河網閾值南小河溝流域為10 km2。依據Pfafstetter流域編碼規則將南小河溝流域劃分為29個天然子流域。進一步地，在各子流域單元內部依照高程劃分等高帶，以此得到流域水文模擬的基本計算單元。南小河溝流域基本計算單元為228個，等高帶平均面積為0.229 km2(見圖4)。

圖4 南小河溝流域子流域與等高帶劃分Fig.4 Sub-basin and contour belt division in Nanxiaohegou river basin

2.3 參數率定與模型驗證

WEP-L模型的參數包括3類：地表面及河道系統參數、植被參數以及土壤與含水層參數。所有參數均有物理意義，可根據實驗觀測數據或遙感數據進行測算。然而，由于參數存在空間變異性，模擬計算時一般采用其在基本計算單元內的概化均值。因此，一些關鍵參數仍需要結合模型檢驗，根據流量實測數據進行調整。模型中需要率定的關鍵高敏感參數有：不同土地利用類型的洼地儲留深、3層土壤厚度、滲透系數K0、地下水透水系數和產水系數、河床透水系數以及河道與坡面的Manning糙率。

改進入滲模型應用中，將南小河溝流域簡單劃分為兩個土壤類型區，即塬面區(黑壚土)和溝道坡面區(黃綿土)，并分別給定兩種土壤膨脹特征曲線、土密度、應力----應變關系曲線和容重與飽和導水系數間關系曲線等。針對滲透系數K0進行了率定。傳統模型(不考慮土壤膨脹性)應用中，只對包括導水系數進行了率定。濕潤鋒吸力采用進氣值一半代替，兩種土壤均為20 cm(見表1)。

表1 改進入滲模型兩種土壤所需水分特征參數Tab.1 Improving the moisture characteristics of two soils in the infiltration model

針對南小河溝流域，模型率定期為2007-2008年，驗證期為2009-2010年。模擬采用變時間步長，即對降雨強度超過10 mm以上的入滲過程采用1 h，坡地與河道匯流采用6 h，其余采用1 d。選擇Nash-Sutcliffe效率系數和相對誤差作為模擬精度表征指標，運用試錯法對模型進行校準，使得模擬期年均徑流量的相對誤差RE盡可能接近0，而其Nash-Sutcliffe效率系數NSE盡可能接近1。

(1)Nash-Sutcliffe效率系數：

(20)

(2)相對誤差：

(21)

式中：RE為模擬月徑流總量的相對誤差，%。

2.4 結果分析

由于南小河溝流域十八畝臺測站具有較完整的序列資料，因此選取十八畝臺斷面利用改進前、后的WEP-L模型進行率定和驗證，模擬計算結果見表2。十八畝臺斷面逐月流量過程模擬如圖5所示，2010年兩種模擬的逐月流量過程差別如圖6所示。

表2 改進前后的WEP-L模型模擬結果Tab.2 Simulation results of the traditional and modified WEP-L model

圖5 考慮與不考慮土壤膨脹性的逐月流量過程模擬Fig.5 Simulation of monthly flow process considering or without consideration of soil swelling

圖6 2010年兩種模擬的逐月流量過程差別Fig.6 Differences in monthly flow processes between the two simulations in 2010

驗證結果顯示，通過考慮土壤膨脹性因素，率定期內十八畝臺斷面月徑流量模擬值與實測值的相對誤差變化不大，Nash-Sutcliffe效率系數由0.33提高到0.57；驗證期內相對誤差由-27.5%減小到22.37%，Nash-Sutcliffe效率系數由0.37提高到0.49。由此可知，改進后的WEP-L模型對南小河溝流域水文過程的模擬效果具有明顯改善。

3 結 語

在利用分布式水文模型開展流域水文過程模擬時，需要考慮土壤膨脹性對降雨入滲產流過程的影響。本文對傳統降雨入滲產流模型進行修正，通過引入膨脹性土壤導水系數和膨脹性土壤飽和含水量，進而對分布式水文模型WEP-L模型進行改進。將考慮土壤膨脹性的改進后的WEP-L模型應用于南小河溝流域。根據相關算法，將流域劃分為不同計算單元，并將模型所需基礎資料包括水文氣象數據、土壤類型數據、水文地質參數、土地利用類型數據、社會用水數據等展布到了相關計算單元上，并生成了模型輸入文件。然后，采用相對誤差和Nash-Sutcliffe效率系數對南小河溝流域十八畝臺測站月徑流量進行了率定和驗證，結果表明通過考慮土壤膨脹性對入滲產流過程影響，不僅改善了模型模擬失真問題，而且可以有效提高模型計算精度。

土壤膨脹性變形在垂直方向上是自重應力和膨脹力共同影響的結果，土壤在兩合力的影響下，導致土壤剖面基質區土壤水分特征參數隨深度變化而變化。土壤膨脹性變形在水平方向上是土壤膨脹力影響的結果，導致水平方向產生大量干縮裂隙。本次模型只考慮了垂直方向變形，而忽略了水平方向變形，進而可能導致在模擬徑流時，考慮與不考慮膨脹土影響的徑流模擬結果差別不大，因此未來還需要考慮水平方向變形影響。

□