1. 引言

随着全球气候变暖的加剧，极端天气事件频发，降雨对边坡稳定性的影响更为复杂和显著，冻土区由于季节性的冻融循环，土壤的物理力学性质在冻融过程中发生显著变化，这进一步增加了滑坡发生的可能性，泥石流会迅速冲向下方的建筑物，对山坡、河岸或不稳定地带的建筑物带来强大的冲击力，导致建筑物遭到严重的损坏，甚至倒塌。因此，研究季节性冻土区在降雨作用下的滑坡稳定性，对加强滑坡风险区域的监测和预警、预防地质灾害、保障人民生命财产安全具有重要意义。

郑等人[1]揭示了降雨强度动态变化特征对滑坡发生运移和危险性动态演进的影响机制，并为区域滑坡预警体系的建立提供科学依据。吴等人[2]分析了滑坡前的地表形变过程，提出了基于动力学基础的多年冻土区小型滑坡破坏机理概念模型。唐等人[3]分析了不同降雨形式对边坡稳定性的影响，指出前期降雨对边坡稳定性影响显著。Rahardjo等人[4]模拟了不同降雨模式下孔隙水压力的变化，并指出前期降雨强度对土坡稳定性有重要影响。张等人[5]通过对多处滑坡的分析，认为前期降雨对边坡稳定性有一定的影响，并且研究了前期降雨对边坡稳定性影响的机理。

当前，针对滑坡稳定性的研究虽取得了显著进展，但针对季节性冻土区降雨作用下滑坡稳定性的研究仍存在诸多不足。一方面，传统的极限平衡法在分析过程中往往将土体假设为刚体，忽略了土体的实际变形情况，导致计算结果不够准确；另一方面，虽然数值模拟方法如FLAC 3D等能够更真实地反映土体的应力、应变和位移等特性，但在冻土区降雨作用下的具体应用研究仍较为匮乏，特别是在复杂冻融循环和降雨耦合作用下的滑坡稳定性方面，仍然需要进一步研究。

因此本文旨在基于FLAC 3D有限差分软件对季节性冻土区降雨作用下的滑坡进行数值分析，通过构建合理的数值模型，考虑冻土在自然条件和降雨条件下的物理力学特性变化，模拟滑坡的形成和发展过程。同时，利用FLAC 3D内置的强度折减法，计算自然条件和降雨条件下滑坡的安全系数，分析降雨对滑坡稳定性的影响机制。

2. 研究区域

2.1. 滑坡地形地貌

果洛藏族自治州玛沁县位于青海省东南部，地处青藏高原腹地的巴颜喀拉山和阿尼玛卿山之间，地理坐标介于东经97˚54'~101˚50'，北纬32˚31'~35˚40'之间。果洛藏族自治州平均海拔4200米以上，海拔4000~5000米的地区约占全区面积的80%左右，其中位于果洛藏族自治州玛沁县的某滑坡如图1所示。滑坡位于古沟槽内，两侧与近南北向延展的基岩脊相连，边界清晰可辨。在路基的外边坡前沿，明显呈现出两个主要的滑动出口，一浅一深，主要由紫红色的砂粘土与砾石土交织而成，其中角砾岩的含量占10%~40%，紫红色的泥岩受中厚层状构造的严格控制，形成了特定的地质结构。受断裂的影响，坡面发育构造裂缝，泥岩破碎，容易被地下水渗入整个坡体内部，对滑坡的稳定性构成了潜在威胁。

Figure 1. Engineering geological map of landslide area

图1. 滑坡区工程地质图

边坡结构如图2所示，根据地质勘探结果，滑坡体被划分为自上而下堆积层、强风化泥岩和中等风化泥岩三个主要部分。堆积层滑坡明确界定为沿崩积层潜水线发生的滑动现象，滑坡体沿接触界面滑动下伏基岩，后拉强烈风化泥岩随之滑动。滑坡总长设计为250 m，高度为53 m，宽度为50 m。

2.2. 滑坡地区气候分析

果洛藏族自治州玛沁县属大陆性寒润性气候，东西部差异较大，西北部寒冷湿润，东南部由寒温潮温逐渐到冷温湿润。如图3所示，从2014~2023年每月平均最低气温为−12˚C~7˚C，年平均气温为−3.8˚C~3.5˚C，气温低，日温差较大。如图4、图5所示，每年降水量多集中在6~9月份，年平均降水量在10~460 mm之间，其中2016年6~9月份的月平均气温相对较高，达到25.7˚C，2016年年平均降水量尤为显著，达到452.28 mm，全年日照时间为2313~2607小时，相对日照45%~63%，一年之间无明显四季之分，冬季寒冷而漫长，时间长达8~9个月，春季干旱多风，夏秋季短而多雨，并常伴有暴雨和冰雹。

Figure 2. Geological structure profile of landslide

图2. 滑坡地质结构剖面图

Figure 3. Monthly temperature from 2014 to 2023

图3. 2014~2023年每月气温图

Figure 4. Monthly precipitation from 2014 to 2023

图4. 2014~2023年每月降水量图

Figure 5. Annual average precipitation from 2014 to 2023

图5. 2014~2023年年平均降水量图

3. 数值模拟

3.1. 理论分析

3.1.1. 降雨入渗强度理论

自然状态下，有饱和及非饱和状态常同时存在于坡内，降雨入渗过程就是非饱和–饱和的转化过程[6]。

饱和土的抗剪强度通过Mohr-Coulomb破坏准则[7]和Terzaghi有效应力原理[8]，可得：

$\tau =c^{\prime }+\left(\sigma -{\mu }_w\right)\tan {\phi }^{\prime }$ (1)

其中： $\tau$ ——破坏时破坏面上的剪应力；

$c^{\prime }$ ——有效黏聚力；

${\phi }^{\prime }$ ——有效内摩擦角；

${\mu }_w$ ——孔隙水压力。

而非饱和土的抗剪强度：

$\tau =c^{\prime }+\left(\sigma -{\mu }_{\alpha }\right)\tan {\phi }^{\prime }+\left({\mu }_{\alpha }-{\mu }_w\right)\tan {\phi }^b$ (2)

其中： ${\mu }_a$ ——孔隙中的气压；

${\mu }_a-{\mu }_w$ ——破坏时破坏面上的基质吸力，当气压为0时，吸力为负的孔隙水压力；

${\phi }^b$ ——抗剪强度随基质吸力而增加的速率。

二维饱和–非饱和的渗流控制方程可根据达西定律得：

$\frac{\partial \left(k_x\frac{\partial h}{\partial x}\right)}{\partial x}+\frac{\partial \left(k_z\frac{\partial h}{\partial z}\right)}{\partial z}+S=\frac{\partial \theta }{\partial t}$ (3)

式中： $k_x,k_z$ ——x和z方向的渗透系数；

h——某点测压管水头；

S——源汇项；

$\theta$ ——含水率；

t——时间。

孔隙水压力，也称静水压力，是降雨诱发滑坡的主要机制[9]。降雨入渗和地下水的流动增加孔隙水压力，降低有效应力，最终会削弱岩土体的物理力学指标，改变岩土体的力学性质，降低边坡稳定性。

3.1.2. Morh-Coulomb本构模型

该模型的破坏包络线对应于一个莫尔–库仑准则(剪切屈服函数)，张拉截止值(张拉屈服函数)，该包络线上应力点的位置由非关联法则控制，为剪切破坏，与此相关的规则为张拉破坏。FLAC 3D中的Mohr-Coulomb准则[10]-[12]是以主应力表示的 ${\sigma }_1$ 、 ${\sigma }_2$ 和 ${\sigma }_3$ ，为该模型广义应力向量的3个分量(n = 3)，相应的广义应变向量的分量为主应变 ${\epsilon }_1$ 、 ${\epsilon }_2$ 和 ${\epsilon }_3$ 。

1) 增量弹性定律

胡克定律关于广义应力和应力增量的增量表达式具有形式：

$\begin{array}{l}\Delta {\sigma }_1={\alpha }_1\Delta {\epsilon }_1^e+{\alpha }_2\left(\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e\right)\\ \Delta {\sigma }_2={\alpha }_1\Delta {\epsilon }_2^e+{\alpha }_2\left(\Delta {\epsilon }_1^e+\Delta {\epsilon }_3^e\right)\\ \Delta {\sigma }_3={\alpha }_1\Delta {\epsilon }_3^e+{\alpha }_2\left(\Delta {\epsilon }_1^e+\Delta {\epsilon }_2^e\right)\end{array}$ (4)

式中， ${\alpha }_1=K+\left(4/3\right)G$ ； ${\alpha }_2=K-\left(2/3\right)G$ 。

$S_i$ 在增量公式中：

$\begin{array}{l}{S}_1=\left(\Delta {\epsilon }_1^e,\Delta {\epsilon }_2^e,\Delta {\epsilon }_3^e\right)={\alpha }_1\Delta {\epsilon }_1^e+{\alpha }_2\left(\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e\right)\\ S_2=\left(\Delta {\epsilon }_1^e,\Delta {\epsilon }_2^e,\Delta {\epsilon }_3^e\right)={\alpha }_1\Delta {\epsilon }_2^e+{\alpha }_2\left(\Delta {\epsilon }_1^e+\Delta {\epsilon }_3^e\right)\\ S_3=\left(\Delta {\epsilon }_1^e,\Delta {\epsilon }_2^e,\Delta {\epsilon }_3^e\right)={\alpha }_1\Delta {\epsilon }_3^e+{\alpha }_2\left(\Delta {\epsilon }_1^e+\Delta {\epsilon }_2^e\right)\end{array}$ (5)

2) 复合失效准则及流动法则

模型中采用的破坏准则为张拉截止的复合Mohr-Coulomb准则，对3个主应力进行标记时，令 ${\sigma }_1\le {\sigma }_2\le {\sigma }_3$ (以拉应力为正，压应力为负)。

在应力空间和 $\left({\sigma }_1,{\sigma }_3\right)$ 平面的破坏准则可表示为图6的形式。由 $f\left({\sigma }_1,{\sigma }_3\right)=0$ 的定义，自A点到B点根据摩尔–库仑失效准则 $f^s=0$ 的破坏包络线为：

$f^s=-{\sigma }_1+{\sigma }_3{N}_{\varphi }-2c\sqrt{N_{\varphi }}$ (6)

Figure 6. Mohr-Coulomb model and failure criterion of geomaterials

图6. 岩土材料 Mohr-Coulomb模型及破坏准则

由B点到C点拉应力屈服函数的定义为：

$f^t={\sigma }_3-{\sigma }^t$ (7)

$N_{\varphi }=\frac{1+\sin \left(\varphi \right)}{1-\sin \left(\varphi \right)}$ (8)

${\sigma }_{\max }^t=\frac{c}{\tan \varphi }$ (9)

式中， $\varphi$ ——摩擦角；

C——凝聚力；

${\sigma }^t$ ——抗拉强度，材料强度不超过 ${\sigma }_{\max }^t$ 。

剪切势函数用两个函数来描述，分别用来定义剪切塑性流动和拉伸塑性流动，该函数对应于一个非关联的流动法则，具体形式为：

$g^s=-{\sigma }_1+{\sigma }_3{N}_{\psi }$ (10)

$N_{\psi }=\frac{1+\sin \left(\psi \right)}{1-\sin \left(\psi \right)}$ (11)

式中， $\psi$ ——剪胀角。

势函数 $g^t$ 对应于拉应力破坏的相关联流动法则：

$g^t=f^t={\sigma }_3-{\sigma }_t$ (12)

对于剪切应力和拉应力处于边界的情况，可由摩尔–库仑流动法则，并通过定义三维应力空间中边界附近的混合屈服函数进行计算。定义函数 $f\left({\sigma }_1,{\sigma }_3\right)=0$ 用以表示 $\left({\sigma }_1,{\sigma }_3\right)$ 平面中 $f^s=0$ 和 $f^t=0$ 所代表曲线的对角线，该函数的表达式为：

$h={\sigma }_3-{\sigma }_t+{\alpha }^p\left({\sigma }_3-{\sigma }^p\right)$ (13)

式中， ${\alpha }^p$ 和 ${\sigma }^p$ 为两个常量， ${\alpha }^p=\sqrt{1+N_{\varphi }^2}+N_{\varphi }$ , ${\sigma }^p={\sigma }^t{N}_{\varphi }-2c\sqrt{N_{\varphi }}$ 。

弹性假设和破坏准 ${\alpha }^p$ 则不一样，分别在 $\left({\sigma }_1,{\sigma }_3\right)$ 平面中位于1区域和2区域(对应于 $h=0$ 区域内或 + 区域)，如图7所示。如果位于1区，则属于剪切破坏，应用由势函数 $g^s$ 确定的流动准则，应力回归到 $f^s=0$ 的曲线上；如果位于2区，则属于拉应力破坏，应用由势函数 $g^t$ 确定的流动准则，应力回归到 $f^t=0$ 的曲线上。

Figure 7. Region definition of flow criterion

图7. 流动准则的区域定义

3) 剪切塑性校正

首先，考虑剪切破坏，方程(10)的偏微分得到：

$\begin{array}{l}\frac{\partial g^s}{\partial {\sigma }_1}=-1\\ \frac{\partial g^s}{\partial {\sigma }_2}=0\\ \frac{\partial g^s}{\partial {\sigma }_3}=N_{\psi }\end{array}$ (14)

将 $\frac{\partial g^s}{\partial {\sigma }_1}$ 、 $\frac{\partial g^s}{\partial {\sigma }_2}$ 和 $\frac{\partial g^s}{\partial {\sigma }_3}$ 用 $\Delta {\epsilon }_1^e$ 、 $\Delta {\epsilon }_2^e$ 和 $\Delta {\epsilon }_3^e$ 替代，分别在等式(5)中得到：

$\begin{array}{l}{S}_1\left(\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}\right)=-{\alpha }_1+{\alpha }_2{N}_{\psi }\\ S_2\left(\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}\right)={\alpha }_2\left(-1+N_{\psi }\right)\\ S_3\left(\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}\right)={\alpha }_1{N}_{\psi }-{\alpha }_2\end{array}$ (15)

用 $f=f^s$ (见公式(6))得：

$\begin{array}{l}{\sigma }_1^N={\sigma }_1^I+{\lambda }^s\left({\alpha }_1-{\alpha }_2{N}_{\psi }\right)\\ {\sigma }_2^N={\sigma }_2^I+{\lambda }^s{\alpha }_2\left(1-N_{\psi }\right)\\ {\sigma }_3^N={\sigma }_3^I+{\lambda }^s\left(-{\alpha }_1{N}_{\psi }+{\alpha }_2\right)\end{array}$ (16)

和

${\lambda }^s=\frac{f^s\left({\sigma }_1^I,{\sigma }_3^I\right)}{\left({\alpha }_1-{\alpha }_2{N}_{\psi }\right)-\left(-{\alpha }_1{N}_{\psi }+{\alpha }_2\right)N_{\varphi }}$ (17)

3.2. 边界条件和材料特性

模拟岩土材料的非线性力学性能时，采用了各项同性弹塑性本构模型，结合了拉伸屈服准则与经典的Mohr-Coulomb剪切屈服准则进行计算，明确定义了两种极端条件——自然状态与饱和状态。自然状态是指岩土材料在经历了夏季和秋季强烈的蒸发过程之后，其内部的水分含量减少到最低水平，岩土材料的强度和硬度会有所提高。而饱和状态是指冻土在热流和水分的共同作用下完全融化的状态，材料内部的孔隙被水完全占据，这会导致材料的强度下降，同时其变形的能力会得到增强。

模型的南北两端和东西两端被施加了固定约束条件，底部边界被设定为完全固定约束，模型的顶部设定为自由状态。根据相关领域文献数据集和实验室数据结果，得到其土体物理力学参数值见表1，暴雨条件下，(1)、(2)、(3)土工参数根据《岩土工程勘察规范GB500212009》的标准，根据实际情况进行调整，相应作减量处理。采用弹性解法计算初始应力场，模型一共划分为9369个网格单元和9505个节点数量，在5 m宽的三维网格模型中(如图8所示)，计算自然状态和饱和状态边坡变形，设置初始位移和初始速度为零，降雨量取最高日平均降雨量44.92 mm计算，共设置4个监测点。

Figure 8. Grid model

图8. 网格模型

Table 1. Physical and mechanical parameters of soil

表1. 土体物理力学参数

4. 分析结果

4.1. 自然条件

自然条件下FLAC3D数值模拟结果如图9所示，玛沁县滑坡总体位移方向沿坡面向下。根据图9(a)可知，x方向最大位移发生在滑坡中下部强风化泥岩和中等风化泥岩交界处，发生轻微变形范围约为0.48 m。由图9(b)可知，z方向的最大位移为0，滑坡处于比较稳定状态。最大剪切应变–应力的不断增加可以反映滑坡内部最不稳定的区域，此区域容易发生破坏。由图9(c)可知，研究区剪切应变最大增幅约为8.627 × 10⁻³ Pa，主要集中在滑坡后缘中下位置，存在剪切破坏的可能。由图9(d)可知，最大主拉应力发生在滑坡后缘中间位置与坡体表面堆积层，最大拉应力增幅约为4.948 × 10⁴ Pa，最容易发生破坏。从图9(e)可以看出，滑坡后缘中间部位和滑坡中上堆积层表面监测点1与监测点2附近同时存在拉升破坏和剪切破坏，且沿坡体向下垂直方向分布，增加了形成连续滑动面的可能性；图中显示自然条件下滑坡体的安全系数为2.34，坡体处于比较稳定状态。从图9(f)可以看出，自然条件下坡体整体处于未饱和状态。

Figure 9. Landslide analysis results under natural conditions. (a) x direction displacement nephogram; (b) Displacement cloud diagram in z direction; (c) The maximum shear strain cloud diagram; (d) Maximum shear stress nephogram; (e) Plastic zone cloud diagram and safety factor; (f) Natural condition saturation cloud map

图9. 自然条件下滑坡分析结果。(a) x方向位移云图；(b) z方向位移云图；(c) 最大剪切应变云图；(d) 最大剪切应力云图；(e) 塑性区云图及安全系数；(f) 自然条件饱和度云图

4.2. 降雨条件

降雨条件下FLAC 3D数值模拟结果如图10所示，玛沁县滑坡中上段堆积层界面处x方向位移最大。根据图10(a)可知，在降雨条件下，滑坡沿坡面总体位移方向保持向下，x方向最大位移达到25.666 m，与自然条件相比位移增加了25.186 m，表明在降雨条件下玛沁县滑坡有更大的位移和变形。从图10(b)可以看出，降雨过程中，坡面z方向最大沉降位移为0.52 m，位于坡底，与自然条件相比位移有所增加，表明在降雨条件下z方向滑动的可能性比自然条件下更大。根据图10(c)可知，与自然条件下相比剪切应变增量增加，

Figure 10. Landslide analysis results under rainfall conditions. (a) x direction displacement nephogram; (b) Displacement cloud diagram in z direction; (c) The maximum shear strain cloud diagram; (d) Maximum shear stress nephogram; (e) Plastic zone cloud diagram and safety factor; (f) Natural condition saturation cloud map

图10. 降雨条件下滑坡分析结果。(a) x方向位移云图；(b) z方向位移云图；(c) 最大剪切应变云图；(d) 最大剪切应力云图；(e) 塑性区云图及安全系数；(f) 自然条件饱和度云图

最大值为1.1046 Pa，表明玛沁县滑坡中上部的堆积层界面监测点3附近是最容易发生裂缝和变形破坏的区域。图10(d)显示，拉伸屈服主要分布在坡体后缘中间位置以及表面堆积层，最大拉应力增量值约为6.3741 × 10³ Pa，而剪切屈服则更广泛，表明存在崩塌变形的潜力。从图10(e)可以看出，滑坡中上堆积层界面同时存在拉升破坏和剪切破坏，且沿坡体向下形成连续滑动面，表明降雨条件下坡体容易发生失稳；图中显示降雨条件下滑坡体的安全系数为1.27，与自然条件相比坡体安全系数显著降低，导致坡体稳定性降低容易发生滑坡。从图10(f)可以看出，降雨条件下坡体堆积层表面处于饱和状态，饱和度为1。

根据上述设置的边坡监测点对自然条件和降雨条件边坡变形演化过程进行跟踪监测，绘制如图11所示各监测点x方向随时间变化曲线。根据图11(a)可知，坡体土体密实度处于比较稳定状态。根据图11(b)可以看出，四个监测点x方向初始阶段坡体处于平稳状态，随降雨时间变化监测点x方向位移缓慢增加，其中监测点2和监测点3的x方向位移增加尤为显著，位移最大达到25.666 m，堆积层表面有效应力总体出现减小趋势，可以推测坡体堆积层表面土体中存在许多微裂缝，渗透性强，会加速雨水渗透导致坡体失稳。监测点2处形成裂隙与坡体堆积层表面贯通，雨水径流汇入裂隙并转为有压入渗，导致边坡土体含水量增大，孔隙压力加速累积，抗剪强度减小，土体相对松散，表层不断破碎滑落，随着降雨时间持续增加导致边坡失稳。

Figure 11. Displacement evolution of monitoring points with time. (a) Natural conditions; (b) Rainfall conditions

图11. 监测点位移随时间发展演化情况。(a) 自然条件；(b) 降雨条件

5. 结论

对发生在玛沁县的滑坡事件进行了滑坡研究，得出以下结论：

1) 在地下水和降雨渗透的共同作用下，斜坡经历了多个阶段的滑坡过程，形成了一个最大变形区域，其总位移达到27.68米。这种滑坡的特点是在滑坡的后部有推移的动作，而在前部则表现为牵引的力量，整体上可以被归类为牵引型滑坡。

2) 在斜坡变形和破坏的过程中，斜坡内部的软弱层首先开始发生缓慢的变形，形成了剪切破坏面。同时，斜坡顶部的堆积层地面出现了微小的裂缝，这些裂缝为地表水和雨水的渗透提供了路径，从而加速了斜坡的变形。

3) 降雨导致沉积物的含水量迅速上升，这使得斜坡的稳定性系数显著下降，并且沉积物发生了软化、增重和水力作用。随后，整个沉积体开始整体滑动，带动了背后的泥岩也一起滑动。在滑动过程中，孔隙水压力对斜坡土体的失稳起到了关键作用。随着含水率的不断累积，斜坡内部土体的有效应力减小，斜坡堆积层表面的中部位移和剪切应变比上部更大，这与实际的滑动特征相吻合。

参考文献