降雨對多年凍土區淺層滑坡失穩的影響研究

周 保， 魏 剛， 魏賽拉加， 沈凌鎧， 張明哲， 常文斌， 臧佳園

（1. 青海省地質環境監測總站，青海 西寧 810007； 2. 青海省環境地質勘查局，青海 西寧 810007； 3. 上海交通大學，上海 200240）

0 引言

活動層沿多年凍土層滑脫引起的凍土淺層滑坡廣泛分布于青藏高原多年凍土地區，是該區域熱融災害的主要類型之一［1-2］。近年來青藏地區氣候變暖、降雨量增大趨勢顯著［3-4］，直接影響多年凍土穩定性，導致凍土淺層滑坡災害頻發［5-6］，造成高寒生態系統破壞，對社會經濟發展產生威脅［7］。鑒于降雨是誘發凍土淺層滑坡的主要外部因素之一［8］，研究降雨對凍土淺層滑坡失穩的影響對青藏地區熱融災害的預防工作具有重要意義。

隨著高寒地區重點工程的修建，降雨引發的多年凍土地區斜坡穩定性問題受到了廣泛關注，已有學者通過模型試驗、數值分析和現場監測等手段開展了降雨對凍土斜坡穩定性影響的研究。Vidie等［9］通過模型試驗研究不同降雨強度對凍土斜坡失穩的影響，提出中等降雨會引起凍脹蠕變和小規模滑坡，強降雨會誘發規模大、速度快的滑坡和泥流；段東明等［10］采用數值模擬對多年凍土斜坡路基進行降雨入滲分析，得出降雨導致斜坡穩定性的減小存在一定滯后；Niu 等［11］基于無限斜坡模型評估多年凍土斜坡穩定性，計算夏季降雨觸發凍土淺層滑坡的臨界地下水位；趙文等［12］開展模型試驗研究凍土斜坡在反復降雨-日曬-凍融循環條件下的演化特征，認為凍融循環引起坡面土體密實度降低，松散土體在降雨過程中流失導致斜坡逐年后退。

凍土淺層滑坡失穩是一個水熱力動態耦合過程，雨水入滲和氣溫變化均會影響凍土斜坡內部的水熱遷移［13-14］，但現有的研究主要關注降雨條件下凍土斜坡失穩過程及穩定性變化規律，氣溫變化作用下降雨對多年凍土斜坡水熱力影響機制尚不明確。本文通過對青海省祁連地區凍土淺層滑坡開展現場調查，總結了當地凍土淺層滑坡發育的地質環境特征與氣候條件，以典型凍土淺層滑坡災害為研究對象，基于現場勘察和原位鉆孔資料，根據當地夏季降雨特征［15］，采用凍土水熱力耦合數值方法，系統研究了氣溫變化條件下低強度連續降雨對多年凍土斜坡水熱力演化的影響，研究結果為青藏地區凍土淺層滑坡災害防治提供了理論依據。

1 青海省祁連縣凍土淺層滑坡

1.1 地質環境特征

基于多源遙感數據調查青海省多年凍土區熱融災害分布，共解譯凍土淺層滑坡災害290處，在青海省祁連縣、治多縣和曲麻萊縣廣泛分布，結合現場調查驗證，祁連縣共發育有凍土淺層滑坡災害54處，占青海省該類型災害總數的18.6%。采用無人機航拍、動力觸探和鉆探等方式對祁連縣凍土淺層滑坡開展了現場調查。

青海省祁連縣位于祁連山中段腹地，地形呈北高南低、西高東低的特點，地貌類型較豐富。經統計，祁連縣凍土淺層滑坡發育的海拔范圍為3 550～3 850 m，主要集中在3 750 m；發育位置多為匯水溝中上部，斜坡坡度集中在10°～20°，整體上陡下緩，發育坡向多為陰坡；斜坡表層植被以異針茅和高山蒿草為主，上部植被較下部稀疏，植被總體覆蓋率約50%～80%，滑體表層植被破壞嚴重；滑體表面由第四紀上更新統坡（洪）積物（Q3dl+pl）組成，下伏基巖主要為三疊系（T）砂巖，活動層土體主要為細粒黏土和泥炭，受凍融循環影響土體力學性質較差；后緣往往發育小型不規則陡坎，底部常有地下水匯集或地下冰出露；滑體向溝口流動堆積，前緣多發育有魚鱗狀草皮。

1.2 氣候條件

青海省祁連縣屬于典型的高原大陸性氣候，地形條件復雜和海拔梯度懸殊導致水熱條件呈顯著的區域性差異和垂直差異。東部降雨多、氣溫較高，西部降雨少、氣溫較低，低海拔地區的氣溫較高海拔地區高；年內氣溫變化符合正弦函數分布，1 月平均氣溫達到最低值-11 ℃，7月平均氣溫達到最高值15 ℃；研究表明當地多年平均氣溫呈波動上升趨勢，氣溫傾向率達0.33 ℃·（10a）-1［16］。

根據當地氣象站資料可得：從祁連縣降雨量月際分配來看，5—9 月多年平均降雨量達到324 mm，占全年總量的84.6%，7 月往往達到降雨量最大值；從祁連縣降雨量的季節分配來看，夏季多年平均降雨量為234 mm，占全年降雨量的57.8%，秋季和春季平均降雨量分別占20.9%和17.8%，冬季僅占3.5%。研究表明近年來青海省祁連山地區降雨量呈波動上升趨勢，尤其進入21 世紀后，降雨量呈快速增大趨勢，降雨量傾向率為46.3 mm·（10a）-1［17］。

1.3 滑坡概況

2018 年8 月末青海省祁連縣默勒鎮海浪村（100°49′58″ E， 37°43′47″ N）發生凍土淺層滑坡（圖1）?，F場調查表明，滑坡高程范圍為3 763～3 800 m，平面呈不規則舌狀，位于陰坡中下部，斜坡整體坡度約15°。根據運動及堆積特征將滑坡分為滑源區和流動堆積區，滑源區長32 m，流動堆積區長58 m。根據2020 年10 月10 日鉆孔資料獲取滑坡地層信息（圖2），地表以下0～1.6 m 為活動層，由粉質黏土、砂質黏土和泥炭組成，其中0～1.4 m 土層較松散，1.4～1.6 m土層相對致密；1.6～12.7 m為多年凍土層，由黏土和泥巖組成，含有大量冰晶；12.7～15.0 m 為基巖層，主要為砂礫巖。結合現場調查和鉆孔資料獲取主滑方向地質剖面如圖3所示。

圖1 祁連縣凍土淺層滑坡Fig. 1 Scene photograph of the active layer detachment in Qilian

圖2 鉆孔地層信息Fig. 2 Stratum information based on borehole

圖3 滑坡體主滑剖面Fig. 3 Section of the main sliding mass

圖4 為根據氣象站觀測資料所得祁連縣降雨概況。根據祁連縣2008—2020 年的年降雨量統計圖［圖4（a）］，2018 年降雨量達574 mm，為2008—2020年間降雨量峰值，2018年6—8月降雨量占全年總降雨量67.6%。根據2018 年6—8 月的日降雨量統計圖［圖4（b）］，當地6 月22 日至7 月9 日（凍土淺層滑坡發生前）連續18 d 總降雨量為162 mm，日均降雨量為9 mm·d-1，最大日降雨量為33 mm·d-1。

圖4 青海省祁連縣降雨概況Fig. 4 Rainfall situation of Qilian， Qinghai Province： annual rainfall in Qilian from 2008 to 2020 （a）；daily rainfall in Qilian from June to August in 2018 （b）

2 水熱力耦合數值模擬

2.1 計算原理

基于有限元軟件COMSOL Multiphysics 求解多物理場偏微分方程，實現多年凍土斜坡氣溫變化和降雨過程中的水熱力演化數值仿真，并用強度折減法計算斜坡穩定性［18］。鑒于氣溫變化和降雨對凍土淺層滑坡的影響是復雜的水熱力動態耦合過程，為簡化這一過程，假設入滲雨水溫度與地表溫度一致，溫度場變化受熱傳遞、水分運移和冰水相變影響；水分場變化由基質吸力驅動，并由孔隙冰阻隔作用控制；水熱過程單向影響土體應力和位移。相應物理場控制方程如下：

（1）溫度場方程

凍土二維傳熱中，忽略對流傳熱、考慮相變的溫度場方程可寫為［19］

式中：ρ為土體的密度（kg·m-3）；C(θ)為土體的熱容量（J·kg-1·℃-1）；λ(θ)為土體導熱系數（W·m-1·℃-1）；T為土體溫度（℃）；L為水凍結成冰的相變潛熱，取335 kJ·kg-1，ρI為冰的密度（kg·m-3）；θI為體積含冰量。

（2）水分場方程

凍土中水分遷移規律與非飽和土中水的運移相似，用Richards方程描述［19］

式中：D為未凍土的擴散率，未凍土的擴散率可寫為，凍土中孔隙冰的存在削弱了水分的遷移效應，因此，引入阻抗系數考慮孔隙冰的隔水作用［20］，凍土擴散率寫為，I=1010θI；θu為體積含水率；ρw為水的密度（kg·m-3）；k為導水率（m·s-1）；c為比水容量（1·m-1）；kg為重力方向導水率（m·s-1）。

（3）應力場方程

土體任意微元的平衡微分方程和土體位移-應變關系表達式可寫為［21］

式中：[?]T為微分算子矩陣；{σ}為應力張量，{σ}={σxσyσxy}T；{F}為單位土體體積力，{F}={FxFy}T；{ε} 為應變，{ε}={εxεyεz}T；{u} 為位移，{u}={uxuy}T。

2.2 網格模型和計算參數

根據地層信息（圖1）和主滑剖面（圖2）建立二維有限元網格計算模型如圖5 所示，斜坡模型長90 m，坡度為15°，地表以下0～1.6 m為活動層，1.6～12.7 m 為多年凍土層，12.7 m 以下為基巖層；采用自由三角形網格，由于活動層水熱力演化劇烈，故將活動層網格細化。

圖5 有限元計算模型Fig. 5 Finite element model

凍土比熱容和導熱系數與未凍水含量和含冰量的關系［22-23］、抗剪強度參數與土溫的關系［24］參考相關研究設置，土水特征曲線和滲透系數采用Van Genuchten 提出的模型［25-26］。土骨架的比熱容和導熱系數分別取1.4×106J·m-3·℃-1和1.3 W·m-1·℃-1；活動層深度0～1.4 m 土層的飽和滲透系數取2×10-6m·s-1，1.4～1.6 m 土層飽和滲透系數取4×10-9m·s-1［27］；其余數值模擬所需參數根據鉆孔取樣進行室內試驗的結果設置如表1所示。

表1 地層物理力學參數Table 1 Physical and mechanical parameters of the formation

2.3 計算方案及約束條件

建立不考慮降雨的模型一和考慮降雨的模型二，對照分析降雨對凍土淺層滑坡的影響：模型一僅施加氣溫變化工況模擬斜坡10 月10 日至次年12月31 日的水熱力演化；模型二以模型一7 月1 日的水熱力計算結果為初始狀態，在氣溫變化的基礎上施加強度9 mm·d-1、歷時18 d 的降雨工況，模擬斜坡7 月1 日至7 月18 日降雨過程中和7 月19 日至7月30日降雨停止后12 d內的水熱力演化。

根據當地氣溫變化情況和氣候變暖特征［28］，結合附面層理論［29］將地表溫度視為平均氣溫并簡化成三角函數形式作為模型上表面溫度邊界條件，表達式如下：

式中：t為時間（s），模型下表面溫度邊界條件為熱流邊界條件［30］，根據深層地熱梯度與土體導熱系數的乘積進行換算，取0.06 W·m-2；模型兩側為絕熱邊界；將鉆孔實測地溫設為初始值，下表面采取熱流邊界條件，計算時間設為100年，使地溫分布基本穩定，結果作為模型一2020 年10 月10 日溫度場初始狀態。模型兩側和下表面為不透水邊界條件；不考慮降雨時上表面不透水邊界條件，考慮降雨時上表面為流量邊界；如表1所示，土體初始體積含水率作為水分場初始狀態。模型上表面為自由位移邊界條件，兩側水平位移為0，下表面為固定位移邊界條件。

3 模擬結果與分析

3.1 降雨對溫度場的影響

2018年7月1日至18日模型一和模型二地表以下0.4 m 和1.6 m 處地溫變化分別如圖6（a）和圖6（b）所示：僅施加氣溫變化時，7 月18 日0.4 m 深度地溫升高0.90 ℃，1.6 m 深度地溫升高0.087 ℃；在氣溫變化基礎上施加降雨工況后，7 月18 日0.4 m深度地溫升高1.84 ℃，1.6m深度地溫升高0.18 ℃。降雨導致0.4 m 深度地溫增長速率隨時間減小，1.6 m 深度地溫增長速率隨時間增大。活動層基底存在0.2 m厚的相對致密層，阻礙雨水入滲，降雨0～13 d對1.6 m深度地溫不產生影響。

圖6 降雨導致不同深度地溫隨時間變化Fig. 6 Ground temperature at different depths varies with time caused by rainfall： ground temperature at a depth of 0.4 m varies with time （a）； ground temperature at a depth of 1.6 m varies with time （b）

3.2 降雨對水分場的影響

圖7 為7月18日和7月30日模型一和模型二的含冰量云圖，結果表明：7 月18 日模型一融深位于1.35 m，模型二融深位于1.49 m，降雨18 d 導致模型二融深下降14 cm，多年凍土上限位置含冰量較模型一減小2.94%；7月30日模型一融深位于1.43 m，模型二活動層已完全融化，融深位于1.60 m，可見雨水滲流對土體含冰量產生持續影響，降雨停止后12 d 模型二融深下降17 cm，模型二多年凍土上限位置含冰量較模型一減小40.74%。

圖7 7月18日和7月30日模型一和模型二含冰量Fig. 7 Ice content of model 1 and model 2 on July 18th and July 30th

圖8 為7月18日和7月30日模型一和模型二的飽和度云圖。7 月18 日模型一0～1.4 m 為飽和度約0.58 的融土層，活動層底部20 cm 為飽和度約0.45的凍土層；7月30日模型一融深增大，凍土層厚度減小為10 cm，飽和度增至0.48。雨水入滲導致活動層飽和度大大增加，7 月18 日模型二融土層的飽和度約0.88，活動層底部13 cm 為飽和度約0.45 的凍土層；由于水分下滲和凍土進一步融化，7月30日活動層完全融化，飽和度約0.85，活動層基底以下形成10 cm厚飽和度約0.90的富水層。

圖8 7月18日和7月30日模型一和模型二飽和度Fig. 8 Saturation of model 1 and model 2 on July 18th and July 30th

根據不同深度滲流速度隨時間變化曲線［圖9（a）］和降雨過程中流速分布［圖9（b）］，分析降雨條件下滲流速度變化規律如下：7 月1 日1.2 m 深度內土體已完全融化，滲流速度隨深度減小，此時融水沿坡面法向遷移，降雨0～3 d 雨水入滲導致流速迅速增加，此時雨水滲透深度有限，流速增大速率隨深度增大而減小；隨著雨水繼續下滲，自降雨3 d起，流速增大速率隨深度增大而增大；1.4 m 深度土體的孔隙冰在降雨9 d融化，流速開始增大；降雨9 d雨水沿重力方向豎直下滲，在孔隙冰的阻隔作用下，降雨18 d 隨土層深度增加，滲流方向愈發接近順坡；降雨停止后，1.4 m 深度范圍內土體流速下降，此時流速隨深度增大而增大；由于1.4～1.6 m 深度為滲透性較差的相對致密層，1.6 m 深度處流速在降雨停止后2 d 才逐漸增大，且流速增量很?。唤涤晖Ｖ购?2 d，滲流方向完全轉變為順坡。

圖9 降雨條件下滲流速度變化規律Fig. 9 Variation of seepage velocity under rainfall condition： seepage velocity at different depths varies with time （a）；distribution of seepage velocity at different stages （b）

3.3 降雨對斜坡穩定性的影響

圖10（a）為模型一年內的安全系數變化曲線，可以得出斜坡穩定性與氣溫均呈正弦函數變化，且兩者呈負相關，穩定性相對氣溫具有一定滯后，最高氣溫和最低氣溫出現在7 月和1 月，而8 月和9 月斜坡穩定性最差，2 月穩定性最好，年內氣溫升高導致安全系數下降約12.1%。

圖10（b）為模型二7 月1 日至7 月30 日安全系數變化曲線，結果表明：降雨導致安全系數下降約49.1%。降雨0～6 d安全系數緩慢下降，6～15 d安全系數下降趨勢顯著增大，安全系數從4.30 下降至2.61，降雨15 d 安全系數下降趨勢逐漸減小，于降雨停止后10 d 達到最小值2.22，此時斜坡仍處于穩定狀態。

圖11 （a）直觀反映了不同時刻斜坡極限狀態下的位移。降雨條件下斜坡位移主要集中于活動層，降雨0、9、18 d 活動層最大位移分別為5.2 mm、8.4 mm 和35.9 mm，降雨停止后第12 d 達到19.1 cm，最大位移所在位置從坡頂逐漸向坡腳移動，降雨停止后12 d 最大位移位于斜坡中部和坡腳間。根據斜坡極限狀態下不同位置的位移隨時間變化曲線［圖11（b）］，降雨期間位移緩慢增大，坡頂的位移較大，坡腳位移較??；降雨入滲至第15 d，坡表變形迅速發展，其中坡中和坡腳位移增大趨勢尤為顯著，坡頂位移增幅相對較小。

圖11 極限狀態下模型二位移變化規律Fig. 11 Variation of displacement of model 2 under limit state： distribution of slope displacement under limit state （a）；variation of slope displacement at different positions under limit state （b）

4 討論

降雨18 d 導致0.4 m 深度地溫升高0.94 ℃、1.6 m深度地溫升高0.093 ℃，可見雨水入滲導致斜坡淺層地溫升高，降雨對地溫的影響隨深度增加而減小。地溫升高促進凍土融化，同一時刻降雨條件下活動層融深較大，降雨停止后12 d 多年凍土上限附近含冰量減小40.74%。雨水入滲和凍土融化導致活動層飽和度提高77%～95%，降雨停止后12 d活動層以下出現了10 cm 厚富水層。分析認為，降雨導致活動層融土含水量增加、容重增大，水分聚集對土顆粒產生動態浮托力［31］，土體力學性質下降；7月1 日至18 日地表溫度達到最大值，地溫隨深度增加而降低，雨水下滲過程中內能向溫度較低的土體轉移，且活動層飽和度上升提高了土的導熱性，有利于大氣與土體之間的熱傳導；降雨條件下融深可能進一步增大，活動層以下細粒富冰凍土層融化產生富水層，孔隙水壓難以消散，土體有液化的可能［32］，且水分的潤滑作用削弱了凍融交界面的抗滑力［33］，發生凍土淺層滑坡的風險大大增加。

極限狀態下斜坡的位移集中在活動層，潛在滑面位于活動層和多年凍土交界面，符合凍土淺層滑坡變形特征。降雨初期水分沿重力方向下滲，此時位移緩慢增大；當雨水入滲量較大時，地下冰起隔水作用，導致滲流方向逐漸轉變為順坡，流速隨深度增大而增大，同一深度斜坡中部流速最大，此時位移顯著增大，最大位移所在位置由坡頂附近轉移至坡中和坡腳之間。由此可見，降雨條件下融土達到一定飽和度后，活動層形成順坡方向的滲流場，產生順坡方向的滲透力［34］，增大了上覆融土的下滑力，此時土體飽和度較高，力學性質降低，導致斜表變形加劇，增大了活動層沿多年凍土層滑脫的可能性。

斜坡穩定性與氣溫均以正弦函數形式變化，兩者呈負相關，相關研究表明，氣溫對多年凍土斜坡穩定性的影響需要進一步考慮短時間尺度下的極端高溫［35-36］和長時間尺度下的凍融循環作用［8］。施加降雨工況后，斜坡安全系數從4.30 下降至2.22，下降幅度達48.4%，降雨入滲前3 d 安全系數下降速率緩慢，3～15 d安全系數迅速下降，第15 d起安全系數下降速率逐漸降低，于降雨停止后10 d 達到最小值，此時斜坡仍處于穩定狀態，推測后續的降雨將進一步降低斜坡安全系數導致斜坡在8 月末失穩。綜上所述，降雨對多年凍土斜坡淺層穩定的影響顯著，雨水入滲對多年凍土斜坡淺層水熱力產生持續性的影響，在低強度、長時間的降雨工況下，斜坡最小安全系數出現時間滯后于降雨數天，斜坡可能在降雨停止后的數天失穩，推測此時斜坡的破壞形式為坡腳附近活動層最先沿凍融交界面薄弱層破壞，從而牽引上部土體下滑引起凍土淺層滑坡。

5 結論

通過COMSOL Multiphysics 建立了多年凍土斜坡模型，在氣溫變化的基礎上模擬了9 mm·d-1、持續18 d 降雨條件下多年凍土斜坡水熱力演化過程，基于模擬結果探討了氣溫變化條件下青藏地區低強度、長持時降雨對凍土淺層滑坡失穩的影響，得出以下結論：

（1）夏季降雨擾動斜坡淺層溫度場，影響了活動層的凍融過程，可能導致活動層以下含冰量較高的細粒土融化形成富水層，增大了斜坡因活動層基底孔隙水壓力積累而失穩的風險。

（2）雨水入滲導致活動層飽和度大幅提高，增大了土體容重、削弱了土體力學性質，在多年凍土層隔水作用下水分沿順坡方向遷移，產生順坡方向的滲透力，加劇活動層變形，對多年凍土斜坡淺層穩定性產生威脅。

（3）降雨對多年凍土斜坡淺層穩定性產生的不利影響較明顯，雨水入滲一段時間后斜坡穩定性迅速下降，且雨水入滲對活動層水熱力的演化產生持續影響，安全系數最小值滯后于降雨數天。

（4）極限狀態下斜坡位移演化規律印證了凍土淺層滑坡是活動層沿多年凍土層滑脫的過程，土體位移變化與雨水入滲形成的滲流場關系密切，低強度的持續降雨導致坡體中下部變形最顯著。