曹 磊, 吳 清, 趙 廣

(1.成都市城市安全與應(yīng)急管理研究院,四川 成都 610023;2.四川省公路規(guī)劃勘察設(shè)計(jì)研究院有限公司,四川 成都 610041;3.中國科學(xué)院地理科學(xué)與資源研究所,北京 100101)

0 引言

降雨誘發(fā)土質(zhì)滑坡是我國西南山區(qū)最為普遍的災(zāi)害之一,也是災(zāi)害防治的重點(diǎn)[1-2]。降雨誘發(fā)土質(zhì)滑坡的過程可劃分為邊坡失穩(wěn)階段和滑坡運(yùn)動(dòng)階段,這2個(gè)階段之間存在緊密關(guān)聯(lián)。針對上述2個(gè)階段,有研究利用室內(nèi)實(shí)驗(yàn)[3-4]、物理模型[5-6]、數(shù)值模擬[7-8]等方法揭示其形成機(jī)理。對于土質(zhì)滑坡而言,降雨入滲是坡體從非飽和狀態(tài)到飽和狀態(tài)的過程,該過程中坡體內(nèi)部物理力學(xué)性質(zhì)發(fā)生顯著改變,從而導(dǎo)致坡體失穩(wěn)破壞。以抗剪強(qiáng)度為例[2,5],降雨入滲導(dǎo)致坡體內(nèi)部含水率增加與基質(zhì)吸力下降,使得土體抗剪強(qiáng)度下降,進(jìn)而造成邊坡失穩(wěn)。

為定量描述降雨入滲導(dǎo)致土體抗剪強(qiáng)度的改變,Lu等[9]將廣義有效應(yīng)力的概念引入淺層邊坡穩(wěn)定性計(jì)算中,并廣泛應(yīng)用于局部飽和條件下的邊坡穩(wěn)定性分析。國內(nèi)外將非飽和土入滲理論引入淺層邊坡穩(wěn)定性研究并取得長足進(jìn)展[10-11]。當(dāng)邊坡失穩(wěn)后,如何定量描述滑坡運(yùn)動(dòng)過程是評價(jià)降雨誘發(fā)滑坡危害的重要前提。一方面,通過物理實(shí)驗(yàn)可測量直接反映滑坡運(yùn)動(dòng)特征的滑坡厚度和滑坡速度等[12-13]參數(shù),隨著測量技術(shù)的發(fā)展,坡體內(nèi)部孔隙水壓力、有效壓應(yīng)力以及地震動(dòng)信號等參數(shù)也間接被用于分析滑坡運(yùn)動(dòng)特征[14-15];另一方面,數(shù)值計(jì)算方法成為滑坡運(yùn)動(dòng)過程定量評價(jià)常用方法之一[16-17],如Savage等[18]首次提出基于深度平均理論的滑坡物理模型,并成功應(yīng)用于顆粒流、滑坡、泥石流等物質(zhì)的運(yùn)動(dòng)模擬中[17-18];此外,侵蝕效應(yīng)[19]、剪脹效應(yīng)[20]等顯著影響滑坡運(yùn)動(dòng)特征的現(xiàn)象也被考慮。盡管針對降雨誘發(fā)土質(zhì)滑坡過程單一階段的研究已較為深入,但有關(guān)降雨誘發(fā)土質(zhì)滑坡過程多階段的耦合研究較為缺乏,可同時(shí)模擬邊坡失穩(wěn)和滑坡運(yùn)動(dòng)過程的數(shù)值算法還有待建立。

因此,本文以降雨誘發(fā)土質(zhì)滑坡過程為研究對象,建立可反映邊坡失穩(wěn)與滑坡運(yùn)動(dòng)2個(gè)階段演化機(jī)理的物理模型,明確2個(gè)階段轉(zhuǎn)化銜接的關(guān)鍵因子,提出適用于2個(gè)階段物理模型同步計(jì)算的數(shù)值方法。搭建考慮水力耦合作用下從邊坡失穩(wěn)到滑坡運(yùn)動(dòng)全過程的數(shù)值模擬流程,實(shí)現(xiàn)降雨誘發(fā)土質(zhì)滑坡過程的定量描述,在此基礎(chǔ)上開展數(shù)學(xué)模型的參數(shù)敏感性分析,旨在為降雨誘發(fā)土質(zhì)滑坡提供一種考慮水力耦合作用的可定量化評價(jià)方法,為此類滑坡防災(zāi)減災(zāi)提供理論與技術(shù)支撐。

1 考慮水力耦合作用的降雨誘發(fā)滑坡耦合數(shù)學(xué)模型

1.1 邊坡失穩(wěn)數(shù)學(xué)模型

水力耦合作用對土質(zhì)邊坡失穩(wěn)階段的影響主要體現(xiàn)在降雨入滲下邊坡內(nèi)部土體含水率與吸應(yīng)力的改變。為描述這一過程,本文采用基質(zhì)勢形式的Richards模型[8],模型中涉及的土—水特征曲線和非飽和滲透系數(shù)等采用Van-Genuchten公式[21]擬合。模型方程表達(dá)如下:

(1)

(2)

(3)

式中:

t—土壤吸應(yīng)力變化的時(shí)間,s;

z—土壤垂向厚度,m;

ψ—土體基質(zhì)吸力,Pa;

θ—土體含水率;

θr—土體殘余含水率;

θs—土體飽和含水率;

C—比水容重,C=dθ/dh,其中,h表示邊坡厚度,m;

Se—土體飽和度,Se=(θ-θr)/(θs-θr);

α,n,m—與邊坡物理特征相關(guān)的擬合參數(shù);

K—非飽和滲透系數(shù);

Ks—飽和滲透系數(shù)。

考慮非飽和條件下的降雨誘發(fā)土質(zhì)邊坡穩(wěn)定安全系數(shù)通常使用基于無限邊坡假設(shè)的一維極限平衡模型進(jìn)行評價(jià)[9],方程表達(dá)如下:

(4)

式中:

FoS—邊坡穩(wěn)定安全系數(shù);

φ—邊坡內(nèi)摩擦角,°;

β—邊坡斜面傾角,°;

cp—邊坡內(nèi)部有效黏聚力,Pa;

ρs—邊坡固相密度,kg/m3;

ρf—邊坡孔隙水密度,kg/m3。

1.2 滑坡運(yùn)動(dòng)數(shù)學(xué)模型

降雨誘發(fā)土質(zhì)邊坡失穩(wěn)后通常處于松散及臨近飽和狀態(tài),變形破壞時(shí)由于剪縮作用會(huì)產(chǎn)生超孔隙水壓力,可降低滑坡底部有效摩擦阻力,提升滑坡流動(dòng)性;相反,對于密實(shí)飽和土體,變形破壞時(shí)由于剪脹作用導(dǎo)致孔隙水壓力消散,可增加滑坡底部有效摩擦阻力,降低滑坡流動(dòng)性。因此,為考慮水力耦合作用對滑坡運(yùn)動(dòng)階段的影響,本文采用可考慮剪脹效應(yīng)的Iverson-George模型[20],該模型由質(zhì)量守恒方程、動(dòng)量守恒方程、體積含水率方程和孔隙水壓力方程組成,方程表達(dá)如下:

(5)

(6)

(7)

(8)

(9)

式中:

u,v—滑坡沿x和y方向的運(yùn)動(dòng)速度,m/s;

gx,gy,gz—重力加速度沿x,y和z方向的分量,m/s2;

pb—滑坡內(nèi)部孔隙水壓力,Pa;

kap—滑坡內(nèi)部側(cè)向土壓力系數(shù),通?？稍O(shè)為1[17];

D—滑坡剪脹率,D=-2K(pb-ρfgzh)/μh;

μ—滑坡內(nèi)部孔隙水黏度,Pa·s;

c—滑坡固相體積分?jǐn)?shù),c=1-θ;

ρ—滑坡密度,kg/m3,ρ=ρfθ+ρsc;

ω—滑坡固相壓縮率;

κ—滑坡剪脹角,°,κ=atan(c-ceq),其中,ceq表示滑坡固相臨界體積分?jǐn)?shù);

χ,ξ—無量綱系數(shù),χ=(ρs+3ρf)/4ρ,ξ=3/(2ωh)+gzρf(ρ-ρf)/4ρ。

上述模型將剪脹理論耦合至滑坡運(yùn)動(dòng)模型,基于滑坡所處狀態(tài)(固相體積分?jǐn)?shù)大小)來判斷其內(nèi)部發(fā)生剪脹或剪縮效應(yīng),當(dāng)滑坡固相體積分?jǐn)?shù)大于臨界值時(shí)表明滑坡處于剪脹狀態(tài),反之則處于剪縮狀態(tài),從而判定孔隙水壓力消散或增加。

1.3 數(shù)學(xué)模型耦合

通過上述數(shù)學(xué)模型分析可見,含水率是同時(shí)影響邊坡失穩(wěn)和滑坡運(yùn)動(dòng)2個(gè)階段的關(guān)鍵參數(shù),可以作為實(shí)現(xiàn)階段耦合的銜接因子,如圖1。通過邊坡失穩(wěn)模型可獲得邊坡內(nèi)部含水率、基質(zhì)吸力、穩(wěn)定安全系數(shù)、失穩(wěn)體積等變量狀態(tài)。當(dāng)滑坡開始運(yùn)動(dòng)后,將邊坡失穩(wěn)體積和含水率作為初始條件代入滑坡運(yùn)動(dòng)模型,可獲得滑坡運(yùn)動(dòng)過程的厚度、速度等變量狀態(tài)。值得注意的是,由于Iverson-George模型是由Navier-Stokes方程基于深度平均理論簡化而來,無法考慮滑坡垂直方向上的孔壓分布,故滑坡物質(zhì)含水率可取垂直方向上的平均值代入。

圖1 考慮水力耦合作用的降雨誘發(fā)土質(zhì)滑坡全過程物理模型

2 數(shù)學(xué)模型求解算法

由于上述數(shù)學(xué)模型方程的結(jié)構(gòu)特征存在差異,且包含變量的演化方向也有所不同,故分別采取相應(yīng)的算法進(jìn)行數(shù)值求解。

針對Richards模型方程,為保持計(jì)算穩(wěn)定性,采用隱式格式離散,具體離散形式如下:

(10)

式中:

k—時(shí)間迭代次數(shù);

i—垂向計(jì)算網(wǎng)格位置;

ψi—非飽和區(qū)基質(zhì)勢在垂向i網(wǎng)格處取值,Pa;

Ci—比水容重在垂向i網(wǎng)格處取值;

Δt—單次迭代時(shí)間步長;

Δz—垂直方向網(wǎng)格長度。

可見,方程(10)屬于典型的三對角矩陣形式,利用Tridiagonal Matrix Algorithm(TDMA)方法[22]求解。

針對Iverson-George模型方程,同樣采用隱式格式離散,具體離散格式如下:

(11)

式中:

j—平面計(jì)算網(wǎng)格位置;

Δx—水平x方向網(wǎng)格長度;

Δy—水平y(tǒng)方向網(wǎng)格長度;

U—變量矩陣;

F,G—通量矩陣;

S—源項(xiàng)矩陣。

可見,方程(11)屬于非線性雙曲格式,利用Finite Volume Method(FVM)方法求解[17]。值得注意的是,為保證2個(gè)階段模型同步計(jì)算時(shí)間的一致性,時(shí)間步長取決于2個(gè)階段中偏小的一方。綜上所述,搭建考慮水力耦合作用的降雨誘發(fā)土質(zhì)滑坡過程數(shù)值流程,如圖2。

圖2 考慮水力耦合作用的降雨誘發(fā)土質(zhì)滑坡過程數(shù)值計(jì)算流程

3 數(shù)學(xué)模型及數(shù)值算法驗(yàn)證

為驗(yàn)證所用各階段物理模型及數(shù)值算法的適用性,首先,對楊詩秀等[23]開展的室內(nèi)一維降雨入滲非飽和土壤試驗(yàn)進(jìn)行模擬,該試驗(yàn)選擇60cm高的砂壤土柱,分2種工況開展,一種將土柱水平放置,一端表面以飽和含水率作為進(jìn)水條件(第一類邊界條件),經(jīng)150min后測定土柱內(nèi)部含水率分布;一種將土柱垂直放置,頂端以0.02427cm/min的供水強(qiáng)度噴灑在土柱表面(第二類邊界條件),經(jīng)180 min后測定土柱內(nèi)部含水率分布。土柱初始含水率θ為0.025,飽和含水率θs為0.4,殘余含水率θr為0.01,滲透系數(shù)K=1.42θ10.24為土壤飽和度函數(shù),擬合參數(shù)取值α=0.025,n=1.391。網(wǎng)格大小為1cm條件下的數(shù)值結(jié)果如圖3中黑色虛線,表明模型數(shù)值數(shù)據(jù)與試驗(yàn)實(shí)測數(shù)據(jù)是一致的。

圖3 2種條件下入滲監(jiān)測數(shù)據(jù)與計(jì)算數(shù)據(jù)對比

其次,對2015年深圳光明新區(qū)紅坳村土質(zhì)滑坡運(yùn)動(dòng)過程開展模擬?，F(xiàn)場調(diào)查顯示,該滑坡發(fā)生時(shí)處于飽和松散狀態(tài),受剪脹效應(yīng)影響具有滑動(dòng)迅速、滑距長等特點(diǎn)[24],計(jì)算模擬參數(shù)通過相關(guān)文獻(xiàn)資料[24-25]獲取,見表1。網(wǎng)格大小為5m條件下的模擬結(jié)果,如圖4。圖4給出了滑坡運(yùn)動(dòng)4個(gè)時(shí)刻的厚度分布,時(shí)間為0.041667h時(shí)刻線條代表實(shí)際滑坡波及范圍,滑坡堆積邊緣處存在幾棟被波及的建筑。結(jié)果表明計(jì)算滑坡堆積范圍與實(shí)際滑坡堆積范圍吻合度較好。

表1 模型計(jì)算所需參數(shù)取值

圖4 深圳滑坡運(yùn)動(dòng)過程不同時(shí)刻滑坡形態(tài)分布

4 數(shù)學(xué)模型參數(shù)敏感性分析

為深入研究水力耦合作用在降雨誘發(fā)土質(zhì)滑坡失穩(wěn)及運(yùn)動(dòng)過程中的影響,利用上述耦合物理模型開展全過程模擬及關(guān)鍵參數(shù)敏感性分析。由于缺乏降雨誘發(fā)土質(zhì)滑坡從邊坡失穩(wěn)到滑坡運(yùn)動(dòng)全過程的實(shí)測資料,本文采取一維數(shù)值模擬方案進(jìn)行研究。具體方案設(shè)置,如圖5。初始邊坡位于坡度39°的斜坡上,邊坡水平長5m,斜坡水平長10m,邊坡初始含水率0.1,飽和含水率0.4,殘余含水率0.05,飽和滲透系數(shù)7×10-9m/s,雨強(qiáng)50mm/h,滑坡固相臨界體積分?jǐn)?shù)0.64,滑面摩擦系數(shù)0.65,其余模型參數(shù)與上述模擬案例一致。

圖5 降雨誘發(fā)土質(zhì)滑坡初始條件示意圖

邊坡失穩(wěn)過程計(jì)算結(jié)果,如圖6。隨著降雨入滲,坡體內(nèi)部含水率由上到下逐漸提升,相應(yīng)的坡體穩(wěn)定性系數(shù)隨之下降,當(dāng)降雨至5h,坡面向下1m處發(fā)生失穩(wěn)破壞,坡體失穩(wěn)部分趨于飽和,平均含水率為0.38。將上述計(jì)算結(jié)果代入滑坡運(yùn)動(dòng)模型模擬滑坡運(yùn)動(dòng)過程(如圖7),由于滑坡初始固相體積分?jǐn)?shù)小于臨界固相體積分?jǐn)?shù),滑坡初始滑動(dòng)階段處于剪縮狀態(tài),孔隙水壓力增大,從而減少滑坡基底摩擦阻力,增加滑坡流動(dòng)性,滑坡最大速度可達(dá)10m/s。當(dāng)滑坡進(jìn)入平面后,摩擦阻力增大,導(dǎo)致滑坡減速堆積,最大水平滑距為32m。

圖6 降雨入滲下邊坡含水率及穩(wěn)定安全系數(shù)變化

圖7 滑坡運(yùn)動(dòng)過程中不同時(shí)刻滑坡形態(tài)

在上述案例的基礎(chǔ)上,選擇滲透系數(shù)及土體含水率作為研究水力耦合作用的關(guān)鍵因子開展敏感性分析。該參數(shù)對降雨入滲過程的影響分析已有較多研究成果[5-6,26],因此本文僅針對滑坡運(yùn)動(dòng)過程進(jìn)行研究。在其余模型參數(shù)相同的條件下選取K=7×10-8,7×10-9,7×10-10進(jìn)行模擬,在相同時(shí)間內(nèi)的滑坡形態(tài),如圖8(a)。從圖8(a)可知,滲透系數(shù)較小,水力耦合作用對滑坡運(yùn)動(dòng)能力影響明顯,這是由于當(dāng)滲透系數(shù)較大時(shí)滑坡內(nèi)部孔隙水壓力耗散迅速,導(dǎo)致滑坡底部摩擦阻力減小幅度變小,從而減少了滑坡運(yùn)動(dòng)距離。同樣,在其余模型參數(shù)相同的條件下選取不同滑坡固相體積分?jǐn)?shù)c=0.61,0.62,0.64進(jìn)行模擬,在相同時(shí)間內(nèi)的滑坡形態(tài),如圖8(b)。從圖8(b)可知,滑坡固相體積分?jǐn)?shù)小于臨界體積分?jǐn)?shù)時(shí),滑坡運(yùn)動(dòng)能力越強(qiáng),滑坡固相體積分?jǐn)?shù)的細(xì)微差異可對滑坡運(yùn)動(dòng)能力帶來顯著變化。綜上可知,滑坡土體物質(zhì)特征(如滲透系數(shù)、含水率等)均對滑坡運(yùn)動(dòng)能力存在著關(guān)鍵影響作用。

圖8 不同條件下的滑坡運(yùn)動(dòng)形態(tài)

5 結(jié)論

針對降雨誘發(fā)土質(zhì)滑坡形成過程,將其劃分為2個(gè)階段:邊坡失穩(wěn)階段和滑坡運(yùn)動(dòng)階段。通過分析水力耦合作用對降雨誘發(fā)土質(zhì)滑坡的影響機(jī)理,采用Ricahrds模型與Iverson-George模型相耦合,以及針對方程結(jié)構(gòu)特征采用不同數(shù)值方法求解,搭建降雨誘發(fā)土質(zhì)滑坡過程數(shù)值模擬流程,實(shí)現(xiàn)該過程的定量計(jì)算及參數(shù)敏感性分析,得到以下結(jié)論。

(1)水力耦合作用是影響降雨誘發(fā)土質(zhì)滑坡過程的重要因素,且對不同階段影響的方式存在差異。邊坡失穩(wěn)階段,主要通過入滲影響邊坡含水率變化;滑坡運(yùn)動(dòng)階段,主要通過剪脹影響內(nèi)部孔隙水壓力進(jìn)而改變滑坡基底摩擦阻力。

(2)通過模擬物理實(shí)驗(yàn)及實(shí)際案例,驗(yàn)證了構(gòu)建的物理模型及數(shù)值算法的適用性。數(shù)值模擬分析表明,土體滲透系數(shù)及含水率均對滑坡運(yùn)動(dòng)能力存在顯著影響。同一初始條件下,土體滲透系數(shù)越小,滑坡運(yùn)動(dòng)能力越強(qiáng),土體趨于飽和時(shí)含水率越大,滑坡運(yùn)動(dòng)能力越強(qiáng),相比于滲透系數(shù)而言土體含水率對滑坡運(yùn)動(dòng)影響程度較大。