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

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

目前，水資源短缺問題越來越突出，為了適應(yīng)變化環(huán)境和人類活動影響下的水文水資源研究，考慮水文參數(shù)和水文過程的空間異質(zhì)性的分布式水文模型成為研究流域水文循環(huán)的一種有效工具，被認(rèn)為是探索和認(rèn)識復(fù)雜水文循環(huán)過程和機(jī)理的有效手段。隨著遙感技術(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é)研究的熱點(diǎn)。土壤膨脹是由于土壤膠體，尤其是黏粒部分的水化所引起的[8]。土壤水分運(yùn)動受土壤吸水膨脹、失水收縮，濕脹干縮過程影響，土壤膨脹變形主要與初始含水量和上覆荷載等有關(guān)，膨脹力和膨脹變形隨土壤增濕程度增加而增加。在吸水膨脹過程中，土壤膨脹變形多反映為垂直向上。國內(nèi)外眾多學(xué)者對膨脹性土壤入滲進(jìn)行了相應(yīng)研究。McGarry等研究了土壤含水量與土壤變形量間關(guān)系，并提出了用于描述土壤濕脹干縮變化特征的三直線模型[9]。蘇寧虎采用分?jǐn)?shù)階偏微分方程建立了膨脹性土壤入滲方程[10]。甘永德采用土壤膨脹特征曲線和土壤應(yīng)力-應(yīng)變關(guān)系曲線建立了膨脹性土壤飽和水分運(yùn)動參數(shù)計(jì)算模型[11]。然而，在進(jìn)行流域水文過程模擬時(shí)對已有分布式水文模型還未考慮土壤膨脹性的影響。

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

1 分布式水文模型WEP-L模型改進(jìn)

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

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

1.2 考慮土壤膨脹性影響的WEP-L模型改進(jìn)

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

取地面為參照面，向下為正(圖1)。膨脹性土壤吸水膨脹變形，土壤變形受土壤膨脹力和自重應(yīng)力影響，導(dǎo)致土壤濕潤區(qū)剖面飽和導(dǎo)水系數(shù)和飽和含水量均隨深度變化而變化，為了便于分析，做如下假設(shè)：①變形前土壤為均質(zhì)土壤；②土壤變形為彈性，即變形無滯后性；③土壤變形只引起土壤孔隙度變化；④入滲過程中存在明確的濕潤鋒面，濕潤鋒面將濕潤區(qū)和未濕潤區(qū)截然分開，濕潤區(qū)土壤達(dá)到飽和，未濕潤區(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)土壤飽和時(shí)，土壤膨脹力引起的孔隙度變化量可以表示為：

(1)

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

(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)系曲線計(jì)算：

ρsp=A+Bln (γZ)

(4)

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

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

de=dew+dep=Δew+Δep

(5)

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

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

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

Ks(e)=K010m(ez-e0)

(7)

(8)

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

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

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

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

積水前：

fnpt=I

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

(9)

積水后：

(10)

忽略地表積水：

(11)

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

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

F=θT-θ0Z=

令：

(13)

則：

(14)

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

(15)

積分得：

(16)

A=SWΔθ

(17)

積水時(shí)刻確定：

(18)

(19)

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

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

在實(shí)際應(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)用改進(jìn)前后的WEP-L模型模擬南小河溝流域的降雨-徑流過程，檢驗(yàn)?zāi)Ｐ透倪M(jìn)效果。

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 計(jì)算單元劃分

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

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

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

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

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

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

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

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

(20)

(2)相對誤差：

(21)

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

2.4 結(jié)果分析

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

表2 改進(jìn)前后的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

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

3 結(jié) 語

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

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

□