考慮土壤膨脹性的流域水文模型應(yīng)用研究

劉海瀅，甘永德，賈仰文，仇亞琴，牛存穩(wěn)

(1.中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；2.青海大學(xué)省部共建三江源生態(tài)與高原農(nóng)牧業(yè)國家重點實驗室，西寧 810016)

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

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

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

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

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

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

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

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

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

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

土壤吸水膨脹變形是土壤膨脹力和自重應(yīng)力共同作用的結(jié)果，假設(shè)土壤膨脹變形是由土壤孔隙度的變化引起的，則當(dāng)土壤飽和時，土壤膨脹力引起的孔隙度變化量可以表示為：

(1)

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

(2)

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

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

(3)

式中：ep為由土壤自重應(yīng)力導(dǎo)致的孔隙度變化量，cm3/cm3；ρsp為由土壤自重應(yīng)力導(dǎo)致的容重變化量，可以采用土壤應(yīng)力----應(yīng)變關(guān)系曲線計算：

ρsp=A+Bln (γZ)

(4)

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

土壤膨脹變形是膨脹力和自重應(yīng)力共同作用的結(jié)果，則合力導(dǎo)致的土壤空隙變化量可以表示為：

de=dew+dep=Δew+Δep

(5)

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

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

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

Ks(e)=K010m(ez-e0)

(7)

(8)

式中：Ks(e)為孔隙度為ez時土壤飽和導(dǎo)水系數(shù)，cm/min；ez為土壤深度為Z時土壤孔隙度；K0為孔隙度為e0時土壤飽和導(dǎo)水系數(shù)，cm/min；m為與土壤孔隙度性質(zhì)有關(guān)的參數(shù)。

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

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

根據(jù)達西定理有：

積水前：

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積水發(fā)生時刻土壤累計入滲量，cm；F為土壤累計入滲量，cm；tp為土壤表層積水發(fā)生時間，min；Ip為土壤表層積水發(fā)生時的時段降雨強度，cm/min；I為x時段降雨強度，cm/min；fpt為積水時段土壤入滲率，cm/min；t為時間，min；θ0為土壤初始含水量，cm3/cm3；θT為濕潤鋒以上土壤飽和含水量，cm3/cm3；SW為濕潤鋒平均吸力，cm；A為參數(shù)；Ks(e)為土壤飽和導(dǎo)水系數(shù)，cm/min。

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

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

2 模型應(yīng)用

應(yīng)用改進前后的WEP-L模型模擬南小河溝流域的降雨-徑流過程，檢驗?zāi)Ｐ透倪M效果。

2.1 研究流域概況

南小河溝(圖3)是涇河支流蒲河左岸的一條支溝，屬于黃土高原溝壑典型區(qū)，位于東經(jīng)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，降雨年內(nèi)不均勻，多集中在7-9月。流域平均氣溫8.7 ℃，年積溫2 700～3 300 ℃，年日照2 454.1 h。流域土壤主要為黑壚土和黃綿土，其中黑壚土廣泛分布于塬面，黃綿土分布于塬面以下溝坡部位。植被類型以暖溫帶森林草原為主。

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

2.2 計算單元劃分

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

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

2.3 參數(shù)率定與模型驗證

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

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

表1 改進入滲模型兩種土壤所需水分特征參數(shù)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效率系數(shù)和相對誤差作為模擬精度表征指標(biāo)，運用試錯法對模型進行校準(zhǔn)，使得模擬期年均徑流量的相對誤差RE盡可能接近0，而其Nash-Sutcliffe效率系數(shù)NSE盡可能接近1。

(1)Nash-Sutcliffe效率系數(shù)：

(20)

(2)相對誤差：

(21)

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

2.4 結(jié)果分析

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

表2 改進前后的WEP-L模型模擬結(jié)果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

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

3 結(jié) 語

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

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

□