基于Green-Ampt模型的多層結(jié)構(gòu)邊坡降雨入滲改進計算方法及穩(wěn)定性影響研究

宋宜祥，尹子航，黃 達

（河北工業(yè)大學土木與交通工程學院, 天津 300401）

土質(zhì)邊坡的穩(wěn)定性分析是巖土工程中十分重要的研究課題。非飽和土質(zhì)邊坡在天然狀態(tài)下常具有較高的穩(wěn)定性，但在降雨的作用下，邊坡表層土體的飽和度增加，邊坡內(nèi)部的孔隙水壓力增大，基質(zhì)吸力有所減少，邊坡穩(wěn)定性將大大降低[1]。降雨強度、歷時、土體力學性質(zhì)等多種因素都會對邊坡的穩(wěn)定性產(chǎn)生影響[2-3]。根據(jù)調(diào)查，在全國290 個縣發(fā)生的地質(zhì)災(zāi)害中，90%的滑坡由降雨誘發(fā)[4]。因此對降雨作用下非飽和土質(zhì)邊坡穩(wěn)定性影響的研究具有重要意義。

Green 等[5]在考慮薄層積水的情況下，基于毛細管理論提出一種新的入滲物理模型—Green-Ampt（GA）模型。G-A 模型很好的預(yù)測了濕潤鋒的發(fā)展，但實際上濕潤鋒以上的土體并非完全飽和?；诖?，彭振陽等[6]根據(jù)Richards 方程和入滲實驗對入滲過程的精確模擬，將入滲面以下的均質(zhì)土體分為飽和層、過渡層和初始層；并發(fā)現(xiàn)過渡層飽和含水率變化至初始含水率符合橢圓形曲線分布。蘇永華等[7]、李誠誠[8]基于上述含水率的分層假定改進G-A 模型，提出均質(zhì)邊坡降雨入滲深度及其穩(wěn)定性的計算方法，并給出過渡層厚度與濕潤層厚度之間的關(guān)系。胡慶國等[9]基于降雨入滲機理及均質(zhì)土體雨水入滲深度理論計算公式，提出多層結(jié)構(gòu)土質(zhì)邊坡降雨入滲模型。

雖然上述學者已經(jīng)將G-A 模型推廣至多層結(jié)構(gòu)土質(zhì)邊坡，但模型精度不足。為此很多學者對滲流過程進行進一步研究，并將其細分為不同階段。如王繼華[10]分析均質(zhì)土體中雨水的入滲過程，將非飽和土邊坡降雨入滲過程分為2 個階段，分別為降雨強度控制階段和土體滲透系數(shù)控制階段，2 個階段的交點為積水點。積水點形成時刻的研究對于多層結(jié)構(gòu)土質(zhì)邊坡穩(wěn)定性的評價格外重要[11-12]。但已有的計算方法（多層結(jié)構(gòu)土）卻忽略上述事實，對計算模型的精度產(chǎn)生較大影響。

上述學者在對雨水入滲深度的研究中，取得了一系列的成果。在此基礎(chǔ)上，學者對邊坡穩(wěn)定性進行分析。 Fredlund 等[13]提出非飽和土體抗剪強度理論，將雙應(yīng)力狀態(tài)變量與多相連續(xù)介質(zhì)學進行結(jié)合。方薇[14]提出非飽和土的非線性抗剪強度包絡(luò)殼模型，根據(jù)土水特征曲線參數(shù)計算該模型抗剪強度。黃夢宏等[15]在利用極限平衡法求解穩(wěn)定性系數(shù)時，通過引入載荷系數(shù)將其轉(zhuǎn)化為求解一系列線性規(guī)劃問題。

學者已認識到過渡層存在的事實，并給出過渡層厚度與濕潤層厚度間的關(guān)系，將入滲理論推廣至多層結(jié)構(gòu)土。但這些研究很少考慮降雨作用下多層邊坡土體中基質(zhì)吸力隨土體含水率變化的事實，同時對邊坡潛在滑動面位置隨降雨歷時的變化考慮不足。

基于此，本文將在考慮基質(zhì)吸力變化與層間積水形成的情況下，結(jié)合飽和層內(nèi)沿邊坡內(nèi)部滲流管網(wǎng)向坡底流失部分雨水的事實，將入滲過程分解為若干個子過程，并基于G-A 入滲模型對胡慶國等[9]提出的傳統(tǒng)多層結(jié)構(gòu)土質(zhì)邊坡降雨入滲深度計算方法提出改進并對子過程進行求解，最后將其合并為整體入滲過程的解。在此基礎(chǔ)上對層間積水點的形成時刻進行預(yù)測，進而分析雨水入滲深度與時間的關(guān)系，并研究隨降雨歷時增長，降雨強度與雨水入滲深度對邊坡不同位置處（濕潤鋒、飽和層）穩(wěn)定系數(shù)的影響。以達到對潛在滑動面位置隨降雨歷時動態(tài)變化分析的目的，進一步提高傳統(tǒng)多層結(jié)構(gòu)邊坡入滲計算方法的精確度與適用范圍。

1 非飽和-飽和滲流分析

1.1 降雨作用下土體含水率變化規(guī)律

G-A 模型（圖1）對濕潤鋒的發(fā)展進行了很好的預(yù)測，其入滲模式可以表示為:

圖1 G-A 入滲模型Fig.1 G-A infiltration model

式中：i—雨水入滲強度/（m·s-1）

ks—土壤飽和滲透系數(shù)/（m·s-1）；

sf—濕潤鋒處基質(zhì)吸力/m；

hd—濕潤鋒深度/m。

彭振陽等[6]在G-A 模型基礎(chǔ)上提出含水率的分層假定模型（圖2），蘇永華等[7]、李誠誠[8]給出飽和層厚度（hs）與濕潤層厚度（hd）之間的關(guān)系，如式（2）所示。

圖2 降雨入滲過程剖面圖Fig.2 Profile of rain fall infiltration process

式中：a、b—各層土的相關(guān)系數(shù)。

胡慶國等[9]提出多層結(jié)構(gòu)土質(zhì)邊坡降雨入滲模型（圖3）；圖3 中，h1、h2分別為第1、2 層土體厚度，θi1、θi2、θi3分別為不同土層初始飽和度，θs1、θs2、θs3分別為不同土層飽和度。第1 層土的雨水滲流強度為i1且濕潤鋒位于該土層時，濕潤鋒到達第1、2 層分界面時間（t1）為：

圖3 多層結(jié)構(gòu)土質(zhì)邊坡含水率分布示意圖Fig.3 Schematic diagram showing water content distribution of soil slope with the multi-layer structure

式中：Δθ1—第1 層土體飽和度差值；

n1—第1 層土的孔隙率。

降雨入滲的深度（hd1）為：

1.2 基于G-A 模型改進多層結(jié)構(gòu)邊坡入滲理論

Mein 等[16]將G-A 模型應(yīng)用到降雨入滲分析，并考慮到降雨入滲過程分為雨水入滲強度控制階段和土體入滲能力控制階段的情況（圖4）。

圖4 入滲率-時間關(guān)系曲線Fig.4 Relation between infiltration rate and time

基于此，本文提出一種新型多層結(jié)構(gòu)土質(zhì)邊坡降雨入滲模型（圖5）。如圖5 所示，計算模型分為3 層，橫坐標為含水率的變化，縱坐標表示雨水入滲深度。hd1、hd2、hd3分別為第1、2、3 層土體中濕潤鋒的深度，hd1（1）、hd1（2）分別為第1 層土體中第1 階段與第2 階段濕潤鋒深度，hd2（1）、hd2（2）分別為第2 層土體中第1 階段與第2 階段濕潤鋒深度。假設(shè)邊坡傾角為α，降雨強度為q，降雨初期，表層土體處于非飽和狀態(tài)，雨水入滲強度小于土體入滲能力，雨水全部滲入土體。此時，降雨在垂直于坡面方向的入滲強度（i1）為：

圖5 多層結(jié)構(gòu)土質(zhì)降雨入滲過程剖面圖Fig.5 Section of soil rainfall infiltration process of the multi-storey structure

隨著降雨歷時的增長，土體入滲能力隨著地表土體的飽和而減小，直至小于降雨強度。此時雨水入滲強度由土體入滲能力控制，考慮不斷變化的基質(zhì)吸力，垂直于坡面方向的入滲強度（i12）為：

式中：sf1—第1 層土濕潤鋒處平均吸力/m；

k1—第1 層土體飽和滲流系數(shù)/（m·s-1）。

由降雨入滲強度連續(xù)性可知存在臨界時刻tp1（積水點）使得i1=i12，則tp1時雨水入滲深度（hdp1）為:

其中：

將式（8）代入式（3）可得臨界時刻tp1為：

濕潤鋒到達第2 層土體的時刻為t1，若tp1＞t1，取tp1=t1。故當降雨歷時t＜tp1時，將式（6）代入式（4）可得此時濕潤鋒深度（hd1）為：

當降雨歷時tp1＜t＜t1時，伴隨著入滲過程的發(fā)展，土質(zhì)邊坡開始出現(xiàn)積水點（tp1），過渡層作用于飽和層的吸力效能不斷減少，加之受邊坡幾何條件與重力的影響，假設(shè)該部分雨水（飽和層）沿邊坡內(nèi)部滲流管網(wǎng)平行于邊坡表面向坡底流動（圖5）。

由Darcy 定律分析各向同性邊坡體可得，第1 層土體飽和層排出雨水的強度（ip1）為：

式中：hs1—第1 層土體中飽和層厚度/m。

將式（2）代入式（11）可得：

式（7）與式（12）的差值即為第1 層土體第2 階段的有效降雨入滲強度（i1(2)）：

第1 層土體的滲流速率（v1）為：

降雨入滲深度（hd1）為：

式中：Δt—計算時間步長/h；

j—所取計算時間步長數(shù)量。

計算時，首先根據(jù)式（9）計算積水點形成的時刻tp1，并根據(jù)tp1將入滲過程分為2 個階段，分別用式（6）、式（13）計算雨水入滲強度控制階段的入滲強度、土體入滲能力控制階段的入滲強度（考慮變化的基質(zhì)吸力）。

當降雨歷時t1＜t＜tp2時,其中tp2為第2 層土體中積水點的形成時刻，滲流初期，第2 層土體入滲能力大于雨水的入滲強度，層間雨水全部滲入土中，此時:

式中：i2—第2 層土體雨水入滲強度/（m·s-1）。

當降雨歷時tp2＜t＜t2時，其中t2為雨水到達第2 層土與第3 層土交界面處的時刻，伴隨著入滲過程的發(fā)展，土體入滲能力越來越小，直至小于雨水滲流強度，此時層間開始出現(xiàn)積水點。將第2 層土體參數(shù)代入式（13），得到：

式中：sf2—第2 層土濕潤鋒處的平均吸力/m；

i2(2)—第2 層土體的有效雨水入滲強度/（m·s-1）；

k2—第2 層土體飽和滲透系數(shù)/（m·s-1）；

hd2—第2 層土體中濕潤鋒的深度/m；

ip2—第2 層土體飽和層排出雨水的強度/（m·s-1）。

將式（17）與式（18）依次代入式（14）、式（15）與式（16），得到hd2?；诖耍瑢⒋朔椒ㄍ茝V至多層結(jié)構(gòu)土質(zhì)邊坡。

（1）根據(jù)降雨時間選擇合適的時間步長，并大致判斷出雨水入滲深度，并根據(jù)土層的入滲率，飽和度以及其他擬合參數(shù)計算濕潤鋒深度。

假設(shè)降雨強度為0 層土的滲透率，記為k0。首先將第1,2,3,···,n層土的滲流過程分為雨水入滲強度控制階段和土體入滲能力控制階段，并計算第r層土的臨界時刻tpr。若tpr＞tr，取tpr=tr，表示第r層土體只存在雨水入滲強度控制階段。然后分別計算每層土體中各個階段的雨水滲流強度，并根據(jù)式（3）分別計算雨水在不同土層中的入滲時間t1,t2,t3,···,tn-1，進而得出濕潤鋒到達第n與n+1 土層交界面的時間為：

根據(jù)式（19）得到濕潤鋒到達第n與n+1 土層交界面時間，據(jù)此初步估計濕潤鋒所在土層。

（2）根據(jù)上述計算結(jié)果推測濕潤鋒所在土層。結(jié)合降雨歷時，選擇合適的時間步長Δt，計算此時的濕潤鋒深度：

在推導(dǎo)多層結(jié)構(gòu)土質(zhì)邊坡降雨入滲深度理論公式時，并沒有考慮邊坡表面產(chǎn)生厚積水情況，這是因為當降雨強度大于表層飽和滲透系數(shù)k1時，雨水的入滲強度最大為k1，邊坡與地面存在一定的角度，邊坡表面無法聚集雨水，僅存在一層薄積水[17]。同時，僅表層土體在降雨入滲過程中可達到完全飽和。降雨后濕潤區(qū)內(nèi)土體的體積含水率大約為0.8～0.9 倍的飽和土體含水率[18]。

2 非飽和土邊坡的穩(wěn)定性分析

2.1 邊坡穩(wěn)定性的極限平衡法

降雨入滲條件下，非飽和多層結(jié)構(gòu)土質(zhì)邊坡多發(fā)生平行于邊坡表面的淺層滑坡（濕潤鋒處破壞）[19-20]。此外，當飽和層到達土層交界面時，易在土層交界面處形成積水，甚至產(chǎn)生平行于坡面的滲流，故飽和層與過渡層交界面亦可能成為潛在的滑動面[21]。以單位寬度、長度為L的邊坡體作為研究對象，分別將飽和層面與濕潤鋒面作為潛在滑動面，采用有限元軟件分析降雨入滲過程中滑動面隨降雨歷時增長的動態(tài)變化過程[22]。

采用飽和層處穩(wěn)定系數(shù)（Fsm）與濕潤鋒處穩(wěn)定系數(shù)（Fsf）評價邊坡的穩(wěn)定性。當潛在滑動面位于濕潤鋒處時，對胡慶國等[7]提出的關(guān)于土質(zhì)路塹邊坡的穩(wěn)定系數(shù)計算方法進行改進并計算。胡慶國等[7]給出土質(zhì)路塹邊坡穩(wěn)定性系數(shù)（Fs）的計算公式：

式中：c—土體有效黏聚力/kPa；

φ—土體有效內(nèi)摩擦角/（°）；

ρ—水的密度/（kg·m-3）；

H—路塹邊坡的高度/m；

y—未被雨水浸濕的土體厚度/m。

為將式（21）推廣至普通邊坡，本文將水的重度（ρg）替換為對應(yīng)土體的飽和重度（γs），并引入濕潤鋒深度（hd）這一變量。此時濕潤鋒影響下的邊坡穩(wěn)定性計算公式為：

式中：r—第r層土體；

γsr—第r層土體飽和重度/（kN·m-3）；

hdr—第r層土體中浸潤鋒深度/m。

潛在滑動面位于飽和層處時，對李誠誠[6]提出的以交界面為潛在滑動面時邊坡穩(wěn)定系數(shù)計算方法進行改進與求解，此時飽和層影響下的邊坡穩(wěn)定性計算公式為：

式中：c0—飽和土有效黏聚力/kPa；

σn0—飽和層面作潛在滑動面時土體正應(yīng)力/kN；

φ0—飽和土有效內(nèi)摩擦角/（°）；

τ0—飽和層土體重度沿坡面的分力/kN；

j—作用于飽和層的滲透力/kN。

降雨歷時t＜tp1時，飽和層土體重度沿坡面方向的分量與飽和層滲透力的計算公式[5-6]為：

式中：γw—水的重度/（kN·m-3）；

hs1—第1 層土飽和層厚度/m；

γs1—第1 層土的飽和重度/（kN·m-3）。

故當t＜tp1時，飽和層處穩(wěn)定系數(shù)（Fsm）為：

降雨歷時tp1＜t＜tp2時

降雨歷時t＞tp2時，土層交界面處形成積水，此時積水點的發(fā)展強度為：

式中：kr—第r層土體的飽和滲透系數(shù)/（m·s-1）。

將式（29）代入式（4）可將積水強度換算為等價的濕潤層厚度：

將式（30）代入式（2）可將積水強度換算為等價的飽和層厚度：

此時飽和層內(nèi)平行于坡面的滲流力為：

飽和層土體重度沿坡面方向的分量為：

某一降雨歷時下，邊坡整體穩(wěn)定系數(shù)Fs為Fsf和Fsm較小值即：

利用式（22）、（34）、（35）分析潛在滑動面隨降雨歷時增長的動態(tài)變化過程，進而對邊坡穩(wěn)定性進行評估。該方法既考慮了降雨過程中土層交界面處存在積水點的情況，又考慮了不同的潛在滑動面，使得邊坡穩(wěn)定性評價更加完善、更加符合實際。

3 模型驗證與分析

3.1 模型建立與計算參數(shù)的確定

選擇文獻[7]中的邊坡計算模型與實驗數(shù)據(jù)進行驗證，該邊坡分別由坡積土、風化土、基巖構(gòu)成（圖6）。如圖6 所示，模型坡比為1∶1.5，邊坡高度為5 m，坡積土與風化土體的厚度均為1 m。滲流計算時，將模型底部以及兩側(cè)設(shè)為不透水邊界，并將降雨強度邊界轉(zhuǎn)化為單位流量邊界施加于計算模型表面。邊坡的力學參數(shù)如表1 所示。

圖6 計算模型Fig.6 Calculation model

表1 土體力學參數(shù)Table 1 Mechanical parameters of soil

本次計算降雨強度為8×10-6m/s。由于邊坡土體的壓實度[23]與粒度的分布[24]均會影響土層滲透系數(shù)的變化，為此坡積土、風化土、基巖的飽和滲透系數(shù)分別取為5.0×10-5，1.39×10-7，1.2×10-8m/s。

3.2 邊坡降雨入滲深度與穩(wěn)定性分析

3.2.1 降雨入滲深度

基于有限元軟件GEO-Studio 進行數(shù)值模擬，并利用式（20）計算不同時刻的雨水入滲深度，計算結(jié)果如圖7 所示。由圖7 可知，降雨歷時5.5 h 時，雨水入滲深度為1 m，雨水到達坡積土與風化土交界面，本文所提計算模型與數(shù)值模擬結(jié)果基本一致，證實了本文模型的正確性。

圖7 降雨入滲深度理論的數(shù)值模擬結(jié)果、本文結(jié)果、文獻[7]結(jié)果曲線Fig.7 Numerical simulation results of the theory of rainfall infiltration depth, results in this paper, and results in literature [7]

降雨歷時7 h 時，雨水入滲深度為1.12 m，此時本文改進理論模型計算值與數(shù)值模擬值結(jié)果開始出現(xiàn)偏差，但總體趨勢基本吻合。

將本文改進理論模型的計算結(jié)果與文獻[7]的結(jié)果進行對比，二者在降雨5.5 h 之前較吻合，這是由于降雨強度小于坡積土飽和滲透系數(shù)時，tp1＞t1，此時取tp1=t1，改進理論模型僅存在第1 階段，退化為文獻[7]的模型。但在降雨5.5 h 之后，文獻[7]的結(jié)果與模型的數(shù)值模擬結(jié)果開始出現(xiàn)偏差，而本文改進理論模型出現(xiàn)偏差的時刻較晚，約在7 h 左右。造成這種現(xiàn)象的原因是，降雨5.5 h 時，濕潤鋒到達交界面處，此時第2 層土體入滲能力大于雨水的入滲強度，交界面雨水全部滲入土中，積水點并未形成。改進的理論值更好的預(yù)測了積水點形成于坡積土與風化土滲透能力相等時，而非濕潤鋒到達土層交界面的時刻。

綜合對比分析2 條理論值曲線總體趨勢，可以得出由于改進后的理論值考慮隨土體含水率不斷變化的基質(zhì)吸力，并較好的預(yù)測了積水點形成時刻，因而其與數(shù)值模擬計算值的相對誤差更小。并且，降水15 h 后，本文提出的理論模型與數(shù)值模擬的結(jié)果非常接近，這表明改進的理論模型計算精度更高。

3.2.2 穩(wěn)定系數(shù)

圖8 為濕潤鋒處的穩(wěn)定系數(shù)（Fsf）與飽和層處穩(wěn)定系數(shù)（Fsm）隨降雨時間（t）的變化曲線。由圖8 可知，穩(wěn)定系數(shù)的變化可分為2 個階段。

圖8 不同位置處邊坡穩(wěn)定系數(shù)變化曲線Fig.8 Slope stability coefficient variation curve at different positions

第1 階段：t＜5 h，此時Fsf＜Fsm，濕潤鋒處穩(wěn)定系數(shù)更小，坡體在此處偏于不安全。若土體的抗剪強度較小，濕潤鋒處極有可能轉(zhuǎn)化為滑動面。這是由于降雨初期，分層假設(shè)下飽和帶土體厚度較小，此時處于降雨入滲第1 階段（降雨強度控制階段），雨水全部滲入土體，且坡體飽和層內(nèi)滲透力較小。

第2 階段：t＞5 h，此時Fsf＞Fsm，潛在滑動面為飽和層面，主要是由于降雨后期，分層假設(shè)下飽和帶土體抗剪強度比濕潤鋒處低；隨著飽和帶土體厚度不斷發(fā)展，坡體飽和層內(nèi)產(chǎn)生的平行于坡表的滲透力也在逐漸發(fā)展，導(dǎo)致非飽和土邊坡發(fā)生平行于邊坡表面的較深層滑坡時，往往在飽和層面處發(fā)生滑動。

降雨5.5 h 時，濕潤鋒到達土層交界面，此時Fsf出現(xiàn)明顯波動。結(jié)合文獻[7]可知，當雨水入滲到土層交界面并滲入風化土后，交界面土體體積含水率增長較快，土體抗剪強度迅速減小，且此時需考慮風化土（第2 層土）的土體性質(zhì)，導(dǎo)致Fsf出現(xiàn)明顯的波動。

降雨6.5 h 時，F(xiàn)sm出現(xiàn)陡降，隨后繼續(xù)相對平緩的下降。這主要是由于隨著降雨歷時的增加（7 h 左右），雨水在分界面處聚集，形成積水點，導(dǎo)致此時該處土體體積含水率增長速度較快，土體抗剪強度迅速減小且基質(zhì)吸力消失，甚至產(chǎn)生于平行于坡表的滲透力，進而導(dǎo)致Fsm陡降。

對比Fsf與Fsm2 條曲線發(fā)展規(guī)律可知：當t＜5 h時，濕潤鋒面處穩(wěn)定系數(shù)更小，為潛在滑動面；而當t＞5 h 時，模型則以飽和層面為潛在滑動面；即計算邊坡穩(wěn)定性系數(shù)時邊坡不安全面的相對位置會隨降雨歷時改變。相較于文獻[7]中分析邊坡穩(wěn)定性時僅以濕潤鋒處為潛在滑動面位置，本文所提邊坡穩(wěn)定性評價方法更為全面。

4 結(jié)論

（1）土體中基質(zhì)吸力對雨水入滲深度的發(fā)展具有重要影響。在考慮基質(zhì)吸力隨降雨入滲過程不斷變化情況下提出的多層非飽和土邊坡雨水入滲深度的分析方法更好的反映了多層結(jié)構(gòu)土層邊坡降雨入滲的過程，得到的理論計算結(jié)果與數(shù)值模擬結(jié)果更加吻合。

（2）在降雨入滲過程中，隨著降雨入滲深度的增加，濕潤鋒處的穩(wěn)定系數(shù)與飽和層處穩(wěn)定系數(shù)不斷降低，濕潤鋒到達土層交界面時，濕潤鋒處的穩(wěn)定系數(shù)出現(xiàn)明顯的波動，積水點在土層交界面形成時，飽和層處穩(wěn)定系數(shù)出現(xiàn)陡降。表明土體的體積含水率隨著雨水到達土層交界面時迅速增大，甚至產(chǎn)生平行于坡面的滲流，黏聚力與內(nèi)摩擦角降低，土體抗剪強度迅速下降。

（3）潛在滑動面相對位置隨降雨歷時的增長而改變，降雨初期，濕潤鋒處穩(wěn)定系數(shù)相對更小，若土體抗剪強度較小或降雨強度較大則此時濕潤鋒處極有可能轉(zhuǎn)化為潛在滑動面，隨著降雨歷時的增長，潛在滑動面轉(zhuǎn)移到飽和層面處。