降雨对多年冻土区浅层滑坡失稳的影响研究

周 保， 魏 刚， 魏赛拉加， 沈凌铠， 张明哲， 常文斌， 臧佳园

（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）。现场调查表明，滑坡高程范围为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 才逐渐增大，且流速增量很小；降雨停止后12 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）极限状态下斜坡位移演化规律印证了冻土浅层滑坡是活动层沿多年冻土层滑脱的过程，土体位移变化与雨水入渗形成的渗流场关系密切，低强度的持续降雨导致坡体中下部变形最显著。