非飽和滲流Richards方程數(shù)值求解的欠松弛方法

陳 曦，于玉貞，程勇剛

（1. 北京交通大學(xué) 土木建筑工程學(xué)院，北京 100044；2. 清華大學(xué) 水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084;3. 武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072）

1 引 言

非飽和土滲流理論是非飽和土力學(xué)理論體系的一個(gè)重要組成部分，Richards方程是非飽和土滲流理論的基本方程，在土石壩滲流、污染物傳輸、凍土滲流相變和邊坡穩(wěn)定性分析等領(lǐng)域有著廣泛的應(yīng)用，例如，污染物在飽和與非飽和土中遷移的基本方程可以通過 Richards方程和對(duì)流擴(kuò)散方程導(dǎo)出[1]。對(duì)于滲流引起的邊坡失穩(wěn)問題，目前主要采用滲流與變形非耦合的分析方法[2]，研究結(jié)果表明，非飽和滲流場(chǎng)和侵潤(rùn)面的計(jì)算對(duì)邊坡的穩(wěn)定性和安全系數(shù)等結(jié)果具有顯著的影響[3－4]，開展Richards方程數(shù)值求解方法的研究具有十分重要的意義。

針對(duì) Richards方程數(shù)值求解過程中的各種問題，研究人員已經(jīng)開展了相關(guān)的研究，如吳夢(mèng)喜等[5]研究了非飽和滲流求解過程中的數(shù)值彌散現(xiàn)象（即土體飽和度比較低處孔隙負(fù)壓計(jì)算的結(jié)果分布規(guī)律較亂，出現(xiàn)與實(shí)際物理表現(xiàn)不符的現(xiàn)象），為消除這一現(xiàn)象提出了變坐標(biāo)的特征有限元法。Richards方程空間離散和時(shí)間差分后，一般可采用非線性迭代方法來求解每一時(shí)間步所對(duì)應(yīng)的方程，常用的方法有 Picard迭代法和 Newton迭代法，Mehl[6]通過比較Picard迭代法和Newton迭代法，得出與Newton迭代法相比，Picard的結(jié)論通常更加簡(jiǎn)單有效。非飽和土的水力傳導(dǎo)系數(shù)（即滲透系數(shù)）具有較強(qiáng)的非線性特征，是導(dǎo)致常規(guī)Richards有限元方程的收斂性較差的主要原因之一，一些研究人員嘗試采用變量變換方法來改善Richards有限元方程的求解性能。Williams等[7]認(rèn)為，非飽和滲流變量（如壓力水頭）在較小空間和較短時(shí)間范圍內(nèi)的快速變化，是引起Richards有限元方程非線性較強(qiáng)和求解困難的主要原因，基于常規(guī)的方程格式和求解方法未必有效，建議針對(duì)變量進(jìn)行變換來改善Richards有限元方程的非線性和收斂性，并研究不同的變換方法。Miller等[8]提出基于誤差控制的自適應(yīng)策略，將其同時(shí)應(yīng)用于Richards方程的空間離散和時(shí)間差分，提出空間和時(shí)間自適應(yīng)求解方案。

在Richards方程的數(shù)值求解過程中，某些參數(shù)尤其是水力傳導(dǎo)系數(shù)的計(jì)算，需要采用欠松弛（under-relaxation）法,不同的欠松弛法對(duì)非飽和滲流的數(shù)值求解精度和計(jì)算效率具有顯著的影響[9－10]。變量變換技術(shù)的性能也是依賴于欠松弛方法[11]，從基本變量進(jìn)行欠松弛方法的研究尤為重要。本文通過數(shù)值實(shí)例指出了現(xiàn)有欠松弛法的局限性，提出一種新的混合欠松弛方法，并對(duì)其實(shí)用性和可靠性進(jìn)行了驗(yàn)證。

2 非飽和滲流Richards有限元方程

經(jīng)過多年的發(fā)展，非飽和土滲流Richards方程已衍生出3種基本格式，即壓力水頭格式（h-form）、混合格式（mixed form）和體積含水率格式（θ-form）。盡管混合格式的Richards方程被證實(shí)具有嚴(yán)格質(zhì)量守恒的特性，壓力水頭格式的Richards方程在實(shí)際工程中的應(yīng)用仍十分廣泛。

混合格式的Richards方程表達(dá)式為

式中：θ 為體積含水率，θ = nSf，n為多孔介質(zhì)的孔隙率，Sf為流體的飽和度；s為源匯項(xiàng)；K =KsKr(Θ)，Kr(Θ)為相對(duì)滲透張量，為有效飽和度Θ的函數(shù)；Ks為飽和滲透系數(shù)張量；(x,t)∈Ω ×(0,T ]。非飽和滲流數(shù)值分析通常采用 Mualem[12]定義的水力傳導(dǎo)系數(shù)，即

式中：ks為飽和水力傳導(dǎo)系數(shù)；有效飽和度或標(biāo)準(zhǔn)化含水率Θ可由van Genuchten[13]定義的土-水特征曲線計(jì)算，即

式中：h為壓力水頭；θ、θs、θr分別為體積含水率、飽和含水率和殘余含水率；av、nv、mv分別為經(jīng)驗(yàn)系數(shù)，mv通常簡(jiǎn)化為mv= 1-1/nv。

根據(jù)比存儲(chǔ)量概念（即土-水特征曲線的導(dǎo)數(shù)）：

式中：mw為土-水特征曲線的斜率；γw為水的單位重度。

利用式（4），將式（1）表示為h-form的Richards方程，即

對(duì)式（5）應(yīng)用Galerkin加權(quán)殘量法，經(jīng)過空間離散和時(shí)間差分后可得h-form的Richards有限元方程：

式中：下標(biāo)n表示時(shí)間步指標(biāo)；M、C分別為等價(jià)質(zhì)量矩陣和水力傳導(dǎo)特征矩陣，由相應(yīng)單元矩陣Me、Ce組裝而成：

式中： N、B分別為形函數(shù)矩陣及其導(dǎo)數(shù)矩陣。右端矢量為b=s+g+f，其中：

式中：z為節(jié)點(diǎn)高程矢量；n為邊界單位法線矢量。

式（6）的增量迭代形式又稱為修正Picard（這里記為mPicard）迭代，表示如下：

當(dāng)式（9）中第一式中 JP替換為準(zhǔn)確雅克比矩陣（JN= J時(shí)），mPicard迭代法成為Newton迭代法。

3 非線性迭代收斂震蕩與欠松弛方法

在迭代過程中，式（2）中第一式水力傳導(dǎo)系數(shù)的計(jì)算需要輸入適當(dāng)?shù)膲毫λ^：

根據(jù)Tan等[9]的研究，當(dāng)=hn+1,m時(shí)，水力傳導(dǎo)系數(shù)采用當(dāng)前時(shí)間步、當(dāng)前迭代步的壓力水頭進(jìn)行計(jì)算，這種未使用欠松弛的方法記為UR0。然而，直接使用當(dāng)前迭代步的壓力水頭來計(jì)算水力傳導(dǎo)系數(shù)可能會(huì)引起如圖 1所示迭代收斂震蕩。Phoon等[10]對(duì)這種迭代收斂震蕩現(xiàn)象進(jìn)行了分析，認(rèn)為上述震蕩現(xiàn)象是由水力傳導(dǎo)系數(shù)相對(duì)滲流量計(jì)算不協(xié)調(diào)造成。為了降低迭代收斂震蕩對(duì)求解精度和計(jì)算效率的影響，可采用適當(dāng)?shù)那匪沙诜椒?。?dāng)采用前一時(shí)間步結(jié)束時(shí)的壓力水頭與當(dāng)前時(shí)間步當(dāng)前迭代步的壓力水頭的均值來計(jì)算水力傳導(dǎo)系數(shù)時(shí)，即

這種方法稱為UR1方法。由于UR1方案簡(jiǎn)單有效，商業(yè)軟件GeoStudio中的SEEP/W模塊采用了這一方法。Tan等[9]和 Phoon等[10]分別研究了 3種不同的欠松弛方法，即上述的UR0、UR1和下面的UR2方法：

這里，UR2方法采用當(dāng)前時(shí)間步最近兩次迭代步的壓力水頭的均值來計(jì)算水力傳導(dǎo)系數(shù)。他們通過一維算例得出“UR2方案在網(wǎng)格較粗的情況下能夠給出更加精確的壓力水頭場(chǎng)，而UR1盡管收斂較快，結(jié)果卻明顯偏離正確解”的結(jié)論[9-10]。采用同一算例對(duì)1 m厚、初始條件為干燥的土層進(jìn)行有表面雨水入滲的模擬，具體參數(shù)見章節(jié) 5.1算例。圖 1給出 10單元網(wǎng)格一維非飽和滲流的計(jì)算結(jié)果，顯示了 UR1只需少量迭代步數(shù)便可達(dá)到收斂，但其收斂值（約為-8.0 m）明顯偏離了解析解-0.0216 m。UR2和UR0都能收斂到較為精確的解，但UR2相對(duì) UR0需要較少的迭代步數(shù)。圖 1的局部放大圖顯示UR2和UR0欠松弛方法作用下迭代法仍存在圍繞正確解的迭代收斂震蕩現(xiàn)象，正是由于這種迭代收斂震蕩導(dǎo)致了Picard或Newton等非線性迭代方法收斂緩慢和數(shù)值求解精度降低等問題。

圖1 3種欠弛松方法作用下Picard方法的迭代收斂行為Fig.1 Convergence behaviors of Picard iteration

4 一種混合欠松弛方法

理論上，上述的欠松弛方法可以統(tǒng)一在下面的廣義欠松弛方法的框架內(nèi)，其表達(dá)形式為

式（14）表示水力傳導(dǎo)系數(shù)計(jì)算所采用的壓力水頭為前一時(shí)間步結(jié)束時(shí)的壓力水頭和當(dāng)前時(shí)間步所有迭代步已知壓力水頭的線性組合。很明顯，采用廣義欠松弛方法來構(gòu)造近似壓力水頭主要存在以下兩個(gè)問題：①隨著迭代步數(shù)的增加，廣義欠松弛方法需要存儲(chǔ)前面所有迭代步的壓力水頭；②公式右端各項(xiàng)壓力水頭的系數(shù)很難確定。針對(duì)上述問題，可只保留式（14）右端的2～3項(xiàng)。根據(jù)對(duì)圖1中欠松弛方法的觀察，做出如下推測(cè)：①前一時(shí)間步結(jié)束時(shí)的壓力水頭hn在欠松弛迭代法中具有阻尼收斂震蕩的作用，如UR1；②為了獲得較為精確的收斂結(jié)果，應(yīng)盡可能保留最近的迭代結(jié)果，即hn+1,m，如UR0和UR2；③增加前幾步迭代的結(jié)果也有阻尼收斂震蕩的作用，如UR2中的hn+1,m-1，但阻尼效果較弱。

基于上面的理論分析、數(shù)值觀察和推測(cè)，提出一種混合欠松弛方法，即

式中：α∈[0,1]和β∈[0,1]為欠松弛經(jīng)驗(yàn)參數(shù)，由于hn的阻尼作用，α也稱為阻尼系數(shù)。

此短項(xiàng)欠松弛法可以看作 UR0、UR1和 UR2的混合方法，記為 mUR(α,β)，而 UR0、UR1和UR2法也可以看作是這種混合欠松弛方法的特例，即當(dāng)α= 0，β= 0，mUR退化為UR0；當(dāng)α= 0.5，β= 0，mUR退化為UR1；當(dāng)α= 0，β= 0.5，mUR成為UR2。根據(jù)前面的Richards有限元方程和欠松弛方法編制了非飽和滲流計(jì)算程序和可視化界面，通過具體實(shí)例指出現(xiàn)有欠松弛方法的局限性，并對(duì)所提出的混合欠松弛方法的優(yōu)越性進(jìn)行驗(yàn)證。

5 數(shù)值評(píng)價(jià)

5.1 一維均質(zhì)土的非飽和滲流

算例1. 對(duì)一維均質(zhì)土的非飽和非穩(wěn)態(tài)滲流問題進(jìn)行模擬[10]。模型飽和體積含水率θs和殘余體積含水率θr分別為0.363 和0.186，van Genuchten模型中的參數(shù)av、nv分別為1 m-1和1.53，飽和滲透系數(shù)為ks= 1×106m/s。模型中，1 m厚土層劃分為10單元有限元網(wǎng)格（即Δz = 0.1 m）。在均勻干燥的初始條件下（設(shè)為h(z,0) = -8.0 m），上表面有水逐漸滲入時(shí)的邊界條件設(shè)為h(0,t) = -8.0 m，h(1.0,t)=0.0。計(jì)算t = 46800 s的壓力水頭分布時(shí)，隨網(wǎng)格和時(shí)間步長(zhǎng)逐漸加密，壓力水頭逐漸收斂，如圖2所示。采用UR2欠松弛方法，當(dāng)網(wǎng)格密度達(dá)到100單元時(shí)，壓力水頭分布已經(jīng)足夠精確，與1000個(gè)單元網(wǎng)格對(duì)應(yīng)的曲線非常接近。采用1000個(gè)單元時(shí)，在0.8 m的位置處，壓力水頭的計(jì)算結(jié)果為-0.0208 m，接近解析解為-0.0216 m[10]。

仍采用10單元有限元網(wǎng)格，令混合欠松弛方法中參數(shù)。β取固定值 0.5，觀察α變化對(duì)求解精度和效率的影響。圖3為mUR(α,0.5)作用下mPicard迭代法在t = 46800 s時(shí)的壓力水頭分布。這里，mUR(0,0.5)變?yōu)閁R2欠松弛方法，其結(jié)果用十字標(biāo)記表示。由圖可見，隨著α的增大，壓力水頭分布逐漸偏離紅色虛線表示的密網(wǎng)格解。當(dāng)α= 0.5時(shí)，壓力水頭分布與三角形標(biāo)記表示的 UR1方法給出的結(jié)果幾乎一致，意味著β參數(shù)影響甚微。

圖2 壓力水頭分布隨網(wǎng)格和時(shí)間步長(zhǎng)加密的收斂情況Fig.2 Convergence of pressure head distribution with the refinement of finite element mesh and time step

圖3 mUR中參數(shù)α 對(duì)壓力水頭分布計(jì)算精度的影響Fig.3 Effects of α in mUR on solution accuracy

通過上述分析可知，當(dāng)網(wǎng)格劃分較粗時(shí)，α應(yīng)足夠小（例如α∈(0,0.01]）才可確保收斂到較為精確的解。相應(yīng)地，將 mUR(α,0.5)作用下 mPicard迭代法在z = 0.8 m高程處解的迭代收斂過程繪制在圖4中，這里選取部分具有代表性的α值。有趣的是，隨著α值的增加，迭代求解的精度逐漸下降，mPicard迭代求解的速度卻在增加，且當(dāng)α值足夠小時(shí)，計(jì)算精度下降不是很明顯。但α= 0.01與α=0（即UR2方法）相比，mPicard所需的迭代數(shù)分別為68步和112步，迭代步減少約為40%；當(dāng)α較大時(shí)，迭代收斂盡管很快，結(jié)果卻偏離正確解較大，這一點(diǎn)與UR1方法相似，印證了前面的推測(cè)①的合理性，證實(shí)了混合欠松弛方法 mUR(α,β)中 hn對(duì)迭代震蕩具有較好的阻尼作用，且此算例中參數(shù)α益取一個(gè)較小的值。

圖4 mUR中參數(shù)α 對(duì)mPicard迭代性能的影響Fig.4 Effects of α in mUR on the performance of mPicard

確定 mUR(α,β)中的參數(shù)α后（設(shè)定為α=0.01），還需要進(jìn)一步研究和評(píng)價(jià)β變化對(duì)求解精度和計(jì)算效率的影響。根據(jù)數(shù)值結(jié)果，mUR(α,β)中β變化時(shí)所得到的壓力水頭分布與圖3中α= 0.01對(duì)應(yīng)的曲線基本上一致（因此這里沒有繪出），表明此算例中求解精度主要由α控制。圖5顯示β參數(shù)對(duì)mPicard的計(jì)算性能的影響。由圖可見，與UR2的112步（圖中紅色虛線標(biāo)定）相比，β取0.5以下一定范圍的值都可以導(dǎo)致約 50%的迭代步數(shù)節(jié)省，表明混合欠松弛方法在保證精度的同時(shí)，可有效提高迭代法的計(jì)算效率。

圖5 mUR中參數(shù) β 對(duì)mPicard迭代性能的影響Fig.5 Effect of β in mUR on the performance of mPicard

5.2 二維非均質(zhì)壩體的非飽和滲流

算例2. 參照于玉貞等的文獻(xiàn)[4]，壩體和2層堤基的非飽和土參數(shù)分別為：堤身土，(θs,θr,av,nv,ks) = (0.441,0.106,1.9 m-1，1.31,6.7×10-7m/s)；堤基黏土層：(θs,θr,av,nv,ks) = (0.490,0.123,0.8 m-1,1.09,1.8×10-7m/s)；堤基砂土層：(θs,θr,av,nv,ks) =(0.411,0.049,7.5 m-1,1.89,2×10-5m/s)。采用 3 組有限元網(wǎng)格（記為A，B和C），分別含有160、320和2008個(gè)4節(jié)點(diǎn)四邊形流體單元。對(duì)于最密的網(wǎng)格劃分C，時(shí)間步長(zhǎng)設(shè)為Δt = 900 s（288個(gè)時(shí)間步），為了提高計(jì)算效率，這里采用稀疏格式的對(duì)稱連續(xù)超松弛（SSOR）預(yù)處理的共軛梯度（CG）迭代法并結(jié)合 Eisenstat技巧[14]來求解每一非線性迭代步對(duì)稱正定線性方程組。壩體一側(cè)初始水位為 17 m標(biāo)高，另一側(cè)保持恒定水位為10 m標(biāo)高，在此水位差條件下達(dá)穩(wěn)態(tài)滲流。以此穩(wěn)態(tài)滲流狀態(tài)為初始條件，壩體一側(cè)水位從17 m標(biāo)高處以1 m d-1的速度開始下降，經(jīng)過3 d，水位下降3 m后壩體內(nèi)非穩(wěn)態(tài)滲流壓力水頭分布如圖6所示。

圖6 水位下降3 m后的壓力水頭分布（單位：m）Fig.6 Distribution of pressure head after 3 m water level drawdown (unit: m)

圖6根據(jù)網(wǎng)格C和UR1方法求出的非穩(wěn)態(tài)滲流壓力水頭場(chǎng)，在網(wǎng)格C條件下，采用UR1、UR2和 mUR幾種欠松弛方法所得到的壓力水頭分布基本一致，因此，可將其用于評(píng)價(jià)粗網(wǎng)格的計(jì)算精度。圖7為根據(jù)網(wǎng)格A和網(wǎng)格B所獲得的0壓力水頭線，時(shí)間步長(zhǎng)分別設(shè)為Δt = 10800 s（24個(gè)時(shí)間步）和Δt = 5400 s（48個(gè)時(shí)間步）。 由圖可見，對(duì)于2個(gè)粗網(wǎng)格，藍(lán)色實(shí)線對(duì)應(yīng) UR1方法所獲得的計(jì)算精度（其累計(jì)迭代數(shù)分別為3612和8207）要明顯好于紅色實(shí)線表示的UR2方法（其累計(jì)迭代數(shù)分別為8579和18179），而且隨著網(wǎng)格加密，UR1的計(jì)算結(jié)果較快收斂于正確解（這里指虛線對(duì)應(yīng)的加密網(wǎng)格獲得的0壓力水頭線）。需要強(qiáng)調(diào)的是，比較上述兩個(gè)算例可知，欠松弛方法的性能會(huì)因問題不同而呈現(xiàn)顯著差別，Tan等[9]和 Phoon等[10]的數(shù)值研究并不全面。

參照UR1方法（即β= 0），固定β= 0，α從0.6變化到1.0，計(jì)算得到的0壓力水頭線與圖7中UR1給出的0壓力水頭線基本一致，表明上述參數(shù)情況下計(jì)算精度變化甚微。圖8為相應(yīng)的累計(jì)迭代數(shù)變化曲線，括號(hào)中的數(shù)字為mPicard未收斂的累計(jì)時(shí)間步數(shù)，這里 mPicard采用的收斂允許值為1×10-5，最大允許迭代數(shù)為500。如圖所示，保持β=0，α從0.5變化到1.0時(shí)，累計(jì)迭代數(shù)基本呈下降趨勢(shì)，表明在這種情況下，采用較大的α可以顯著減少計(jì)算量。對(duì)于網(wǎng)格 A和網(wǎng)格 B，mUR(1.0,0)作用下的累計(jì)迭代數(shù)分別是2772和5770，與UR1作用下的累計(jì)迭代數(shù)（分別為3612和8207）相比，計(jì)算量分別節(jié)省了23%和30%。

圖7 水位下降3 m后的0壓力水頭線Fig.7 The 0 pressure head curves after 3 m water level

圖8 mUR中參數(shù) α 對(duì)計(jì)算量的影響Fig.8 Effects of α in mUR on computing cost

6 結(jié) 論

（1）非飽和滲流數(shù)值求解時(shí)，水力傳導(dǎo)系數(shù)相對(duì)滲流量計(jì)算不協(xié)調(diào)會(huì)造成迭代法收斂震蕩，欠松弛方法的主要作用是對(duì)這種迭代震蕩施加一種阻尼。通過數(shù)值試驗(yàn)可知，對(duì)于迭代法的收斂震蕩，前一時(shí)間步結(jié)束時(shí)的壓力水頭hn的阻尼抑制效果明顯。

（2）算例演示了現(xiàn)有兩種常用欠松弛方法計(jì)算性能的顯著差異，表明單一欠松弛方法所存在的局限性。在一維粗網(wǎng)格入滲算例中，UR1方法收斂較快，但收斂結(jié)果卻偏離正確解較大，UR2方法作用下非線性迭代法可收斂到較為正確的解，但收斂速度相對(duì)較慢。二維堤壩滲流算例與一維入滲算例結(jié)論不同，UR1方法收斂較快且計(jì)算精度也明顯好于UR2。

（3）基于廣義欠松弛方法，發(fā)展出的混合欠松弛方法mUR可以看作是UR1和UR2的泛化方法。通過 mUR(α,β)中的系數(shù)可調(diào)整阻尼程度和欠松弛程度。數(shù)值試驗(yàn)表明文中所提出的混合欠松弛方法克服了單一欠松弛方法的局限，具有更好的適用性和計(jì)算性能。

[1]李錫夔. 飽和－非飽和土壤中污染物運(yùn)移過程的數(shù)值模擬[J]. 力學(xué)學(xué)報(bào),1998,30(3): 321－331.LI Xi-kui. Numerical modelling of pollutant transport in unsaturated/saturated soils[J]. Acta Mechanica Sinica,1998,30(3): 321－332.

[2]陳曦,劉春杰. 逐步完善的非飽和土邊坡穩(wěn)定性有限元分析方法[J]. 巖土工程學(xué)報(bào),2011,33(增刊1): 380－384.CHEN Xi,LIU Chun-jie. Staged development of finite element methods for stability of unsaturated soil slopes[J].Chinese Journal of Geotechnical Engineering,2011,33(Supp. 1): 380－384.

[3]黃潤(rùn)秋,戚國(guó)慶. 非飽和滲流基質(zhì)吸力對(duì)邊坡穩(wěn)定性的影響[J]. 工程地質(zhì)學(xué)報(bào),2002,10(4): 343－348.HUANG Run-qiu,QI Guo-qing. The effect of unsaturated soil suction on slope stability[J]. Journal of Engineering Geology,2002,10(4): 343－348.

[4]于玉貞,林鴻州，李榮建,等. 非穩(wěn)定滲流條件下非飽和土邊坡穩(wěn)定分析[J]. 巖土力學(xué),2008,29(11): 2892－2898.YU Yu-zhen,LIN Hung-chou,LI Rong-jian,et al.Stability analysis of unsaturated soil slope under transient seepage flow state[J]. Rock and Soil Mechanics,2008,29(11): 2893－2897.

[5]吳夢(mèng)喜,高蓮士. 飽和-非飽和土體非穩(wěn)定滲流數(shù)值分析[J]. 水利學(xué)報(bào),1999,30(12): 38－42.WU Meng-xi,GAO Lian-shi. Saturated unsaturated unsteady seepage numerical analysis[J]. Journal of Hydraulic Engineering,1999,30(12): 38－42.

[6]MEHL S. Use of Picard and Newton iteration for solving nonlinear ground water flow equations[J]. Ground Water,2006,44(4): 583－594.

[7]WILLIAMS G A,MILLER C T,KELLEY C T.Transformation approaches for simulating flow in variably saturated porous media[J]. Water Resources Research,2000,36(4): 923－934.

[8]MILLER C T,ABHISHEK C,FARTHING M W. A spatially and temporally adaptive solution of Richards’equation[J]. Advances in Water Resources,2006,29(4):525－545.

[9]TAN T S,PHOON K K,CHONG P C. Numerical study of finite element method based solutions for propagation of wetting fronts[J]. Journal of Geotechnical andGeoenvironmental Engineering,ASCE,2004,130(3):254－263.

[10]PHOON K K,TAN T S,CHONG P C. Numerical simulation of richards equation in partially saturated porous media: Under-relaxation and mass balance[J].Geotechnical and Geological Engineering,2007,25(5):525－541.

[11]CHENG Y G,PHOON K K,TAN T S. Unsaturated soil seepage analysis using a rational transformation method with under-relaxation[J]. International Journal of Geomechanics,ASCE,2008,8(3): 207－212.

[12]MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research,1976,12(3): 513－522.

[13]VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America,1980,44(5):892－898.

[14]CHEN X,TOH K C,PHOON K K. A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equations[J]. International Journal for Numerical Methods in Engineering,2006,65(6): 785－807.