胡海軍， 李博鵬， 田堪良， 巴亞東， 崔玉軍

(1.西北農(nóng)林科技大學 水利與建筑工程學院， 陜西 楊凌 712100; 2.中國科學院水利部水土保持研究所， 陜西 楊凌 712100; 3.法國國立路橋大學 納維爾研究所巖土力學實驗室,巴黎 77455)

非飽和黃土具有顯著的水敏性力學特性，即其強度和變形模量與含水率密切相關(guān).黃土滑坡與降雨、灌溉和渠道滲漏等具有很大的關(guān)聯(lián)性.建于黃土地層之上的構(gòu)筑物在浸水如降雨入滲情況下會產(chǎn)生沉降變形，而計算浸水環(huán)境作用下地層水分遷移過程是計算邊坡穩(wěn)定性和地基沉降的前提.降雨入滲是發(fā)生最為頻繁的地層浸水情況.積水入滲是降雨入滲的極端情況而通常被室內(nèi)試驗所采納.有關(guān)實際降雨情況下的水分入滲，已有大量的試驗研究.如,李毅等[1]通過室內(nèi)試驗研究了不同降雨強度下黃土地層的水分遷移規(guī)律.基于近似物理模型如Green-Ampt模型的解析方法和求解Richard方程的數(shù)值方法是計算積水和降雨入滲條件下土體水分遷移過程的兩類基本方法.

基于近似物理模型的解析方法計算簡便，特別是考慮浸潤鋒處吸力合理計算后的Green-Ampt修正模型能較好地模擬積水情況下累積入滲量和時間的關(guān)系[2]，卻不能較好地計算水分剖面以及實際浸潤鋒遷移過程.為了模擬積水情況下水分剖面變化過程，王文焰等[3]在Green-Ampt模型基礎(chǔ)上提出了考慮積水入滲過程中原狀黃土水分剖面形狀的入滲模型.本文應用該模型預測重塑黃土積水入滲情況水分遷移過程時，發(fā)現(xiàn)由于其不合理地考慮浸潤區(qū)吸力作用而過快地預測了浸潤鋒遷移速率；結(jié)合Green-Ampt修正模型，本文提出了能合理考慮浸潤鋒處吸力作用和水分剖面形狀的改進模型，并基于Mein等[4]降雨分析理論，將該模型引入到降雨情況下的水分剖面計算.

求解Richard方程的數(shù)值方法能嚴謹求解水分入滲過程，然而非飽和滲透參數(shù)的準確選取是影響該方法預測準確度的重要因素[2]，有時不能采用間接法如VG-Mualem模型獲得的非飽和滲透系數(shù)參數(shù)而需要反演[5].基于近似物理模型的解析方法同樣存在該問題，鑒于浸潤鋒前進法[6]或改進的浸潤鋒前進法[7]能較精確測試非飽和土的滲透系數(shù)函數(shù)，故本文應用改進的浸潤鋒前進法來獲得非飽和滲透系數(shù)函數(shù).

1 模擬方法

1.1 改進模型

1.1.1積水入滲情況

對于積水入滲情況，Green-Ampt模型采用的水分剖面如圖1所示，其將實際入滲深度為Zf的浸潤體簡化為長度為Zs的飽和體，水力梯度采用式(1)計算，根據(jù)達西定律得到入滲量增量，該值等于飽和體長度增量下土體水量的增量如式(2)所示，由式(1)和式(2)得到飽和體深度Zs和時間的關(guān)系如式(3)所示.

(1)

Ksidt=(θs-θi)dZs

(2)

(3)

式中：i為水力梯度;H為積水水頭；Zs為Green-Ampt模型的浸潤鋒深度(也即飽和體深度)；Sm為Green-Ampt模型浸潤鋒Zs處的吸力水頭(取為正值)；Ks為飽和滲透系數(shù)；θi為土體初始體積含水率；θs為土體飽和體積含水率.

圖1 積水入滲下Green-Ampt模型、文獻[3]模型以及實際的水分剖面示意圖

Fig.1 Water content profiles adopted by the Green-Ampt model, the model proposed by Wang Wenyan et al. and the actual water content profile in ponding infiltration

該模型浸潤鋒處吸力水頭Sm并非土體初始含水率θi對應的吸力水頭.Philip[8]認為Sm是一個數(shù)學量，沒有實際的物理意義，該值遠小于浸潤鋒下土體初始含水率對應的吸力水頭Si.Bourwer[9]最早給出了Sm的計算公式如式(4)所示，并認為Sm可近似取為增濕過程中進水值,即增濕過程中氣不連續(xù)時對應的吸力值，然而該值不容易測定，其建議取為飽和樣(對應于抽真空飽和樣)脫濕過程中進氣值的50%，該確定方法比較簡單可行[10].在浸潤鋒以上為飽和的假設(shè)下，Mein等[4]應用式(5)獲得Sm，該值在初始含水率變化較大范圍內(nèi)相差很小，其對不同初始含水率和不同降雨強度下采用相同的值.Mein等[11]經(jīng)過推導認為應用式(4)獲得Sm是合理可靠的.Neuman[12]通過推導得到了浸潤鋒處吸力值表達式與式(4)相同.Brakensiek[13]通過比較，表明上述各種方法可得相近的結(jié)果.對于采用上述等效吸力的模型稱為Green-Ampt修正模型.已有大量研究表明， Green-Ampt修正模型在模擬入滲量和時間的關(guān)系上具有較高的準確度[2,14]，說明上述確定Sm的方法有足夠精度.從公式(4)可見，Sm小于初始含水率對應的吸力水頭Si.

(4)

(5)

式中：Si和Ki分別為初始含水率對應的吸力水頭和滲透系數(shù);S和K分別為吸力水頭和滲透系數(shù).

Green-Ampt模型及修正模型僅能計算如圖1所示的等效飽和體深度.王文焰等[3]根據(jù)大量原狀黃土現(xiàn)場積水入滲試驗中實測水分剖面，假定入滲深度為Zf的區(qū)域近似由上部長度為Zf/2的飽和區(qū)(嚴格來講是傳導區(qū)，因為飽和區(qū)和過渡區(qū)的區(qū)間都很小[15]，其中傳導區(qū)和過渡區(qū)為接近飽和的區(qū)域，所以這里仍用飽和區(qū)稱之)和下部長度為Zf/2、水分剖面為橢圓的非飽和浸潤區(qū)組成(見圖1)，提出了能計算不同時刻水分剖面的模型，然而在計算浸潤鋒處等效吸力Sm上卻沒有采用上述確定方法，其采用下部Zf/2非飽和浸潤區(qū)平均初始含水率對應的吸力Si施加于上部飽和區(qū)底部，得到上部Zf/2飽和區(qū)的水力梯度如式(6)所示，根據(jù)達西定律和水量平衡原理，得到式(7)，將式(6)代入式(7)整理可得入滲深度Zf和時間t的關(guān)系如式(8)所示.如上文所述，作用于長度Zs飽和體的等效吸力Sm小于Si，而將Si作用于長度為Zf/2的飽和區(qū)，可見過大估計了下部吸力的作用.該模型[3]計算所得現(xiàn)場積水入滲試驗中的水分剖面運移過程較為符合實測，可能在于原狀黃土常存在節(jié)理和大孔隙，節(jié)理[16]或大孔隙[17]存在導滲的作用，加快了浸潤鋒遷移速率，對此問題應采用相應模型[16-17]進行分析，而不應單獨在模型中增加吸力作用.

(6)

Ksidt=0.5(θs-θi)dZf+0.125π(θs-θi)dZf

(7)

(8)

(9)

(10)

式中：ε為浸潤區(qū)占入滲區(qū)的比例，其隨Zf的增加而線性增加，可表示為ε=aZf+b.

(11)

(12)

式中：Zs為如圖1所示的飽和體長度，按浸潤區(qū)為橢圓過渡，經(jīng)過換算Zs=(1-ε+0.25πε)Zf.

Zf=

(13)

(14)

1.1.2降雨入滲情況

對于不同降雨強度下的水分入滲模擬，Mein等應用Green-Ampt修正模型進行了降雨入滲分析[4]，得到與求解Richard方程的數(shù)值計算方法相近的入滲率和時間關(guān)系，但其并沒有引入水分剖面形狀來預測實際水分剖面.本文應用上述改進模型，獲得實際水分剖面.當降雨強度小于飽和滲透系數(shù)時，不會發(fā)生徑流，入滲率f(t)一直等于降雨強度P.當降雨強度大于飽和滲透系數(shù)時，在地表發(fā)生徑流前，入滲率等于降雨強度；發(fā)生徑流時，降雨強度P等于入滲強度Ksi，據(jù)此可得徑流發(fā)生時等效飽和體長度Zs，進而根據(jù)式(14)得到此時的實際入滲深度Zf；根據(jù)此時入滲量Fp等于等效飽和體范圍內(nèi)土體水量變化，可得Fp如式(15)所示，進一步可得徑流發(fā)生的時刻tp如式(16)所示.徑流發(fā)生后，假定水及時排走沒有積水水頭，即發(fā)生積水水頭為0的積水入滲，入滲量F和時間的關(guān)系可按式(17)計算，入滲深度Zf由式(18)和式(14)計算得到.

(15)

(16)

(17)

(18)

1.2 基于求解Richard方程的數(shù)值方法

Hydrus-1d軟件應用迦遼金有限元法求解給定初始條件和邊界條件下的一維Richard方程，獲得非飽和土水分遷移過程.一維Richard方程如式(19)所示.對于本文積水和降雨入滲中初始含水率沿深度均相等情況，初始條件如式(20)所示，下邊界條件為自由排水邊界如式(21)所示.對于積水入滲情況，上邊界條件如式(22)所示；對于降雨入滲情況，徑流發(fā)生前，上邊界條件如式(23)所示，當計算到上邊界水頭為0的時刻，即為徑流時刻，從該時刻開始，將上邊界條件定義為水頭為0的邊界，即應用式(22)所示邊界條件并將H設(shè)置為0；另外還分析了停雨后1 h的水分再分布，該過程將上邊界設(shè)置為流量為0的邊界,即將式(23)所示P設(shè)置為0.求解時需要h和θ的關(guān)系，以及K與h或θ的關(guān)系；然而有時不能采用間接法(如由VG模型持水曲線函數(shù)獲得的非飽和滲透系數(shù)函數(shù))而需要反演計算[5]，本文采用改進的浸潤鋒前進法[7]進行求取，具體參數(shù)確定過程見下文.

(19)

θ(z,0)=θi

(20)

(21)

h(0,t)=H

(22)

(23)

式中：z坐標取向上為正，地表處z=0；此處h為吸力水頭，當孔隙水壓力為負時取為負值;P為降雨強度.

2 兩類模擬方法的驗證與分析

2.1 積水入滲情況水分運移計算的驗證與結(jié)果分析

2.1.1試驗情況及參數(shù)確定

為檢驗上述兩類模擬方法在模擬重塑黃土積水情況下水分遷移的可靠性，以延安治溝造地工程建設(shè)中開挖邊坡土料制取的兩種干密度重塑黃土土柱積水入滲試驗[19]為對象進行模擬.兩種土樣干密度分別為1.35 g·cm-3和1.53 g·cm-3，土柱初始質(zhì)量含水率均約為12.5%，初始體積含水率分別為0.174和0.194，積水水頭為2 cm，土柱高220 cm.在不同高度處布置水分儀獲得了積水過程中的含水率變化，土柱底部為透氣板，且入滲速率較慢，浸潤鋒下的氣壓阻滲作用可以不考慮.

為了能用上述兩類方法模擬試驗中水分運移過程，本文測試了飽和滲透系數(shù)和持水曲線，各分析方法所需參數(shù)如表1所示.其中,間接法1采用室內(nèi)制取相同初始狀態(tài)的土樣進行增濕或脫濕測試的持水曲線(圖2)，應用VG模型如式(24)擬合；對制成相同初始狀態(tài)的試樣進行浸水飽和，應用變水頭法測得的飽和滲透系數(shù)Ks，應用VG-Mualem模型[20]即式(25)間接獲得非飽和滲透系數(shù)函數(shù)，進而獲得Sm.間接法2除飽和滲透系數(shù)應用改進的浸潤鋒前進法根據(jù)積水入滲土柱試驗獲得的飽和滲透系數(shù)外，其他參數(shù)與間接法1相同.直接法則采用改進的浸潤鋒前進法[7]根據(jù)積水入滲土柱試驗獲得的非飽和滲透系數(shù)函數(shù).

(24)

K=KsΘ0.5[1-(1-Θ1/m)m]2

(25)

圖2 實測的持水曲線及VG模型擬合函數(shù)

Fig.2 Measured water-retention curves along with VG water-retention functions fitted to independent data

表1 各模擬分析中所需參數(shù)

注：① 兩種干密度抽真空飽和試樣滲透系數(shù)分別為6.50×10-5cm·s-1和3.32×10-5cm·s-1，由于浸水飽和試樣與一維土柱試驗入滲后飽和度接近，滲透系數(shù)均采用浸水飽和試樣的滲透系數(shù);② 根據(jù)浸水飽和試樣及一維土柱試驗入滲后土柱含水率確定θs，該含水率對應土柱入滲試驗中傳導區(qū)含水率[20],并不是完全飽和的含水率;③ 根據(jù)式(5)計算得到兩種干密度試樣Sm分別為65.6 cm和80.1 cm，鑒于式(4)更為嚴格，這里僅采用式(4)計算值.

其中應用改進的浸潤鋒前進法獲得非飽和滲透參數(shù)的過程如下：根據(jù)土柱上各測點含水率開始變化的時間即浸潤鋒達到該處的時間，按冪函數(shù)關(guān)系擬合得到浸潤鋒遷移距離和時間的關(guān)系，對于兩種干密度試樣，所得結(jié)果分別如式(26)和(27)所示.應用該關(guān)系很容易得到不同時刻浸潤鋒遷移速率vZf，以其中一個測點A為例，按式(28)計算得到t1～t2時間段通過A斷面的水流速度v，浸潤鋒前進法[12]采用式(29)計算相應水力梯度i，依據(jù)改進的浸潤鋒前進法[7]相應采用式(30)計算t1～t2時間段A斷面的水力坡降i，根據(jù)達西定律便可得到該時間段平均吸力下的滲透系數(shù)，選取不同的時間段便可得到不同吸力下的滲透系數(shù).圖3給出了應用改進的浸潤鋒前進法所得不同吸力下的滲透系數(shù)結(jié)果.應用VG模型即式(24)，對所得非飽和滲透系數(shù)和基質(zhì)吸力的關(guān)系進行擬合，擬合得到的飽和滲透系數(shù)在抽真空飽和試樣和浸水飽和試樣所測滲透系數(shù)之間.可見,雖然浸水飽和試樣和土柱入滲試驗后的飽和度接近，但在滲透系數(shù)上仍存在試樣不對等性.文獻[21]則建議采用入滲試驗后的土柱進行滲透試驗以獲得飽和滲透系數(shù).由于應用式(23)擬合非飽和滲透系數(shù)時，所得a和n與原持水曲線擬合參數(shù)不同，因此需要在新的a和n下，重新對持水曲線擬合，得到殘余含水率θr分別為0.114和0.138，雖然該參數(shù)有較大改變，但前后兩次所得持水曲線在大于初始含水率時都很接近測試點，而入滲過程是增濕過程，所得持水曲線和非飽和滲透系數(shù)函數(shù)在含水率大于初始含水率時具有足夠的精度.

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

Zf=1.053t0.622

(26)

Zf=0.692t0.592

(27)

[vZf(t1)+vZf(t2)]

(28)

(29)

i(zA,t3)]

(30)

式中：θ(zA,t2)為時刻t2、深度zA處測點A的體積含水率；Ψ為基質(zhì)吸力，由前面所測持水曲線擬合函數(shù)根據(jù)含水率反算得到；vZf為t1～t2時間段浸潤鋒遷移速度.

2.1.2模擬結(jié)果分析

圖4給出了各方法所得浸潤鋒入滲深度和時間的關(guān)系以及與實測結(jié)果的對比.總體上，王文焰等提出的模型[3]過快地預測了入滲過程，Green-Ampt修正模型較慢地預測了入滲過程，而采用相同參數(shù)情況下改進模型較Green-Ampt修正模型提高了浸潤鋒遷移速率，相比前兩種方法均更加接近實測值.針對采用不同方法確定的參數(shù)，采用間接法1所得參數(shù)誤差較大，而采用間接法2或直接法所得參數(shù)預測結(jié)果均較為接近實測值，說明飽和滲透系數(shù)的準確性在很大程度上決定了預測的準確度.另外整體上，對干密度1.53 g·cm-3土柱預測較為準確，這可能與1.35 g·cm-3土柱具有導滲的大孔隙有關(guān).

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

圖5給出了土柱不同深度測點體積含水率隨時間變化的實測與預測結(jié)果.本文改進模型和Hydrus軟件數(shù)值方法都能模擬出淺部測點體積含水率變化快，深部測點體積含水率變化稍慢的特點.總體上，對于干密度1.53 g·cm-3土柱預測較好，而對干密度1.35 g·cm-3土柱預測稍差.

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

圖6給出了浸潤鋒到達各測點時的水分分布預測值和實測結(jié)果的對比.本文改進模型和Hydrus軟件數(shù)值方法所得結(jié)果很接近，只是在浸潤鋒位置處稍有差異.有限的實測點結(jié)果接近于兩種方法所得的水分分布線，說明兩種方法的可靠性.另外,根據(jù)Hydrus軟件所得結(jié)果，分析水分剖面從飽和到非飽和的過渡點，得到飽和區(qū)(嚴格來講為傳導區(qū))所占比例隨著入滲深度的增加由0.3變化到0.6.

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

圖6 積水入滲過程浸潤鋒到達各測點時的水分分布預測值和實測結(jié)果對比

Fig.6Comparison of predicted and measured water content profile when wetting front arrived at measurement points in ponding infiltration

2.2 降雨入滲情況的驗證與分析

為了獲得積水入滲和降雨入滲的差別以及本文改進模型相對于Green-Ampt修正模型在降雨入滲過程分析的適宜性，對比了Green-Ampt修正模型、本文改進模型和Hydrus軟件數(shù)值方法所得降雨情況下浸潤鋒入滲規(guī)律和水分剖面，并與積水情況進行了對比.

結(jié)合取土地區(qū)的降雨情況資料[22]，本文選取特大暴雨1.170 cm·h-1、大雨0.208 cm·h-1和中雨0.104 cm·h-13種雨強下降雨24 h和停雨24 h作為分析情況.對于Green-Ampt修正模型和本文改進模型，按式(16)計算得到干密度1.35 g·cm-3和1.53 g·cm-3土柱僅特大暴雨情況下分別在3.4 h和0.73 h發(fā)生徑流，其他兩種雨強在24 h內(nèi)均未發(fā)生徑流；Hydrus分析結(jié)果表明24 h內(nèi)也僅特大暴雨下發(fā)生了徑流，徑流時間分別為4.5 h和0.8 h，這與Green-Ampt修正模型和本文改進模型所得結(jié)果相近.

圖7給出了積水和3種雨強條件下浸潤鋒入滲深度和時間的關(guān)系.對于Green-Ampt修正模型和本文改進模型，徑流之前，由于土柱表面處于非飽和狀態(tài)不能應用改進模型預測，圖中僅給出了徑流發(fā)生后的入滲深度和時間的關(guān)系.從結(jié)果上可見，本文改進模型比Green-Ampt修正模型更加接近Hydrus所得結(jié)果.對比3種雨強和積水入滲結(jié)果，特大暴雨情況比較接近于積水入滲情況，特別是干密度為1.53 g·cm-3土柱.這也說明了積水入滲可以作為降雨入滲的極端條件予以研究.

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

圖7 降雨24 h過程中入滲深度和時間的關(guān)系以及與積水入滲的對比

Fig.7 Infiltration depth versus time during 24 h rainfall compared with that under ponding condition

圖8給出了特大暴雨情況下，在降雨24 h和停雨24 h時間段，兩種干密度土柱水分入滲及水分再分布過程水分剖面圖.Hydrus分析結(jié)果表明，初始入滲時土柱表面由非飽和狀態(tài)逐漸過渡到飽和狀態(tài)，Green-Ampt修正模型和本文改進模型不能預測徑流發(fā)生前的水分分布變化；發(fā)生徑流時及降雨24 h時，本文改進模型所得水分分布與Hydrus所得水分分布接近，特別是干密度1.53 g·cm-3土柱，而Green-Ampt修正模型預測水分剖面精度較差.相對于降雨24 h，發(fā)生徑流時Green-Ampt修正模型和本文改進模型所得結(jié)果與Hydrus分析結(jié)果差異較大.這主要是由于徑流發(fā)生的時間不同導致兩者入滲量不同，另外如前面分析表明入滲深度淺時，飽和區(qū)所占比例小，改進模型在此階段采用了相對較大的飽和區(qū)比例.降雨結(jié)束后，改進模型不能預測水分重分布，Hydrus分析得到土柱表面含水率減少并且更深處水分增加，這與已有實測規(guī)律[1]相符.需要說明的是，該分析過程宜采用增濕后脫濕的持水曲線函數(shù)，另外停雨后一般涉及到蒸發(fā)，需要測試或計算蒸發(fā)強度E，此時可將式(23)中P替換為-E的上邊界條件進行計算.

a 干密度1.35 g·cm-3土柱

b 干密度1.53 g·cm-3土柱

3 結(jié)論

在Green-Ampt修正模型和王文焰等提出的模型[3]基礎(chǔ)上，提出了可較為合理預測入滲過程中水分剖面及遷移過程的改進模型，應用各模型和求解Richard方程的數(shù)值方法模擬了室內(nèi)重塑黃土積水入滲試驗，并對比了特定雨強下本文改進模型和數(shù)值方法所得的水分剖面.主要結(jié)論如下：

(1) 本文改進模型合理考慮了浸潤鋒處吸力作用和水分剖面形狀，相比Green-Ampt修正模型和王文焰等提出的模型，能更好地模擬室內(nèi)非飽和重塑黃土積水入滲試驗中水分遷移過程.

(2) 根據(jù)土柱入滲試驗各測點水分變化，采用改進的浸潤鋒前進法[7]獲得的飽和滲透系數(shù)和非飽和滲透系數(shù)函數(shù)時，能較好地預測水分遷移過程，而采用室內(nèi)制樣獲得的飽和滲透系數(shù)和間接法獲得的非飽和滲透系數(shù)函數(shù)時預測誤差較大.

(3) 根據(jù)文獻[4]提出的降雨入滲理論，將改進模型推廣到降雨情況下徑流發(fā)生后的水分剖面預測，所得結(jié)果與求解Richard方程的數(shù)值方法分析結(jié)果接近，特別是在試樣干密度比較大的情況，改進模型相對數(shù)值方法計算簡便、易于應用.另外,應用Hydrus軟件可以分析得到徑流發(fā)生前和停雨后的水分遷移變化過程，結(jié)果與已有試驗規(guī)律相符，顯示出較強的模擬能力.