水热变化对季节冻土区牵引式滑坡变形的作用机制——以2020年8月11日西和县大型滑坡为例

梁小鹏， 陈虹举， 燕强珍， 杨佳亮， 马敏

(1. 甘肃工程地质研究院， 兰州 730030； 2. 中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室， 兰州 730000; 3. 东北林业大学寒区科学与工程研究院， 哈尔滨 150040)

冻融引发的滑坡是季节冻土区滑坡的一种主要形式，该类型滑坡的主要特点是滑动速率较小且变形具有季节特性[1-3]。其主要变形机制是季节性冻融引发的水分迁移导致边坡下部产生了水分集聚，从而在集聚区附近形成软弱层[4-5]。在外部荷载或者土体自重作用下，软弱层可能会演化成滑动面，并最终产生滑坡[6-7]，这是目前季节冻土区滑坡发生机制的主要理论成果。然而，该理论仅仅来自定性分析[8]，并无定量的解释和说明，因此难以指导工程实践。

要对该类型滑坡进行定量分析，必须要建立相应的水热力耦合数值模型[9-11]。目前，饱和状态下冻融土的水热耦合、水热力耦合模型已经取得了一定的成果，并且主要用于各种工程建设的工后评价，如青藏公路和青藏铁路多年冻土融化后的路基沉降问题等[12-16]。其中，Zhang等[8]运用饱和状态下冻土的水热模型计算了滑坡发生过程中的水热变化，然后运用单独的本构模型对该滑坡变形进行了模拟，进而对该类型滑坡的变形机制做了初步解释。然而，该模型假定土体饱和，并且运用两个单独的模型分别模拟了水热和变形过程，因此无法反映水热变化对变形的作用机制。即便如此，也极少有研究运用该饱和土模型对冻融引发的滑坡进行数值分析[4,6,17-18]。主要原因在于饱和状态的假设太过于绝对，不符合冻土中水分迁移的实际情况[19-22]。另外，目前对各类型的滑坡来说，数值研究的主要作用都是对滑坡灾害的评价，极少研究其诱发机制，这也源于理论模型的欠缺[9,23-24]。

由于缺少完整的适用于滑坡变形的非饱和冻融土的水热力耦合数值模型，结果导致该类型滑坡的变形机制只能停留在定性分析层面。比如，一种解释认为季节冻融作用会导致边坡产生季节性冻结滞水促滑效应，这种作用会从3个方面影响边坡稳定性：一是冻结会导致边坡下部水分集聚；二是水分增加会降低土体强度；三是水分集聚还会增大边坡的静、动水压力[25]。上述分析从概念模型的角度解释了季节冻土区滑坡产生的原因，具有一定的合理性。但是，水分如何作用于滑体、水分的来源和消散过程并不清楚[26-27]。此外，水分富集的位置及其原因，水分迁移如何影响变形等原因尚不清楚，这也导致现有的边坡处置措施一直未能起到良好的效果[26,28-31]。如果能够建立一种符合季节冻土区边坡水热和变形相互作用过程的理论模型，并结合上述定性分析结果，对冻融作用引发滑坡的过程、变形特点和机制进行定量研究，那么可以对定性分析的结果进行验证和补充，从而为该类型滑坡处置措施的探索提供更加充分的理论指导。

为此，现从微观角度首次基于热力学观点，建立季节冻土区非饱和土水分、温度和变形相互耦合的数值模型，阐述水分对边坡变形的作用机制。该模型的优点是能够反映冻土边坡中冻融作用引发的水分迁移过程和水分不均匀分布对土体强度的影响，以此来分析边坡的变形特性。并以2020年8月11日甘肃陇南市白冯村发生的大型滑坡H3为例，分析水热变化对季节冻土区牵引式滑坡的作用机制。对滑坡处置前、后边坡的温度场、水分场和变形场进行计算，并对滑坡产生的原因和处置前后边坡的变形特点进行分析。最后讨论目前滑坡治理存在的问题，并结合研究结果给出该类型滑坡可能的治理方式。

1 模型与方法

水分作为滑坡产生的主要原因，水分的来源和去向至关重要。季节冻结会导致水分向冻结锋面迁移从而形成水分聚集。当冰透镜体融化后水分液态水增加必然导致活动层土体强度大幅降低，从而产生蠕动。由此，基于质能平衡的非等温水热气耦合模型，建立了季节冻土区边坡的水分迁移方程如下。

(1)

式(1)中：θw、θv和θi分别为液态水、气态水和冰体积含量，如无特殊标注，下标w、i、v和s分别指的是液态水、冰、气态水和土的性质；ρw、ρv和ρi分别为液态水、气态水和冰的密度；h为水头高度；Kwh和KwT分别为基质势和温度梯度引起的等温和非等温水力传导率；Kvh和KvT分别为等温和非等温气态水传导率。

根据Van Genuchten模型，非饱和土的土水特征曲线表达式为

Θ=[1+(-αh)n]-m

(2)

式(2)中：Θ为有效饱和度；α、n和m为拟合参数。

等温条件下，非饱和土渗透系数[32]为

Kwh=10-ΩθiKsΘl[1-(1-Θ1/m)m]2

(3)

式(3)中：Ω为量纲为1的经验参数；l为Mualem模型调整参数；Ks为饱和渗透系数。

非等温液态水渗透系数、温度梯度引起的气态水渗透系数和基质势梯度引起的气态水渗透系数分别为KwT、KvT和Kvh，表达式分别为

(4)

式(4)中：GwT为评价土水特征曲线受温度影响的经验参数；γ为表面张力，γ0为25℃时的表面张力；M和g分别为水的摩尔质量和重力加速度；R为气体常数；Hr为相对湿度；η为水气扩散增强因子；D为水气在土中的扩散度。

温度梯度作为水分迁移的另一个主要驱动力，冻融产生的温度梯度会导致水分向活动层迁移，从而产生水分聚集。根据能量守恒，传热方程[32-33]为

(5)

式(5)中：C为土体的等效比热；λ为土体等效导热系数；Lv为水的蒸发潜热。其余模型参数表达式见表1。

表1 模型相关参数表达式Table 1 Expression of relevant parameters of the model

牵引式滑坡的变形主要为蠕变变形，本构方程必须考虑水分变化的影响，这也是本模型的一个主要创新点。基于黏塑性理论，建立季节冻土区边坡的本构方程为

Δσ=DT(Δε-Δεvp-Δεv)

(6)

(7)

式中：DT为与温度和含水量相关的弹性矩阵；Δε为总应变增量矢量；Δεvp为黏塑性应变增量矢量；Δεv为水分相变引起的应变增量矢量；ET为弹性模量；vT为泊松比。其他主要力学参数表达式如下。

(8)

(9)

(10)

(11)

式中：Tp为土体冻结温度，273.15 K；其余参数均为试验拟合结果。

式(1)～式(11)组成了非饱和冻土水-热-力耦合数值模型。

2 案例分析——以2020年8月11日西和县大型滑坡为例

2.1 研究区地质概况

研究区位于陇南市西和县，该地区处于两组大断裂带的交汇处，其中一组断裂带经过成县-西和-礼县，横穿西礼盆地南部，呈北西走向，另一组断裂带走向北东，即为西礼断陷的东缘主断裂，该断裂带向南西方向延伸，与武都化马一带的北东向断裂相接。从这两组大断裂对新生代地层的控制和对现代地形地貌的影响可以推断，它们均属于活动性断裂构造。主要地貌有构造侵蚀低中山地貌和侵蚀堆积河谷地貌。该区域受温带大陆性季风气候影响，年平均气温为8.4 ℃，标准冻深0.42 m，最大冻深约为1.0 m。

2.2 滑坡介绍

滑坡发生后，甘肃工程地质研究院立即对该滑坡进行全面勘测，查明了新老滑坡周界(图1)，钻探确定了滑动面深度，对岩土体结构特性进行了实验分析。经研究发现，H3为老滑坡，该滑坡自2019年8月2日以来一直处于极慢速滑动状态。本次因连续强降雨影响，于2020年8月11日在H3滑坡体上又发育1处次级滑坡H3-1，该滑坡截止到2020年8月26日一直处于蠕变滑动状态(表2)。

表2 模型相关参数表达式Table 2 Expression of relevant parameters of the model

图1 白冯村滑坡全貌图Fig.1 Overall view of Baifeng Village landslide

H3老滑坡(东经105°16′50″～105°17′44″，北纬33°59′54″～34°00′21″)长约390 m，宽约320 m，滑坡前缘高程为1 604 m，滑坡后缘高程为1 729 m，相对高差为125 m，滑向为140°，钻孔显示滑体平均厚度为11.5 m，滑坡体积143.52×104m3，属于大型滑坡。

H3-1滑坡长约132 m，宽约125 m，滑坡前缘高程为1 603 m，滑坡后缘高程为1 646 m，相对高差为43 m，滑向为139°，根据钻孔显示滑体平均厚度为9.2 m，滑坡体积15.18×104m3，属于小型滑坡。

2019年8月2日—2020年8月11日期间该边坡一直处于蠕变滑动状态，H3-1滑坡中部发育多条横向张拉裂缝，截止到2020年8月26日，H3滑坡导致白冯村三社9户居民75间房屋倒塌，136间房屋出现不同程度竖向裂缝，通村硬化道路中断损坏约240 m，经济损失约460万元，灾情等级为中型。H3滑坡威胁西和县白村民三社69户村民361人的生命财产安全，同时威胁十天高速西和隧道。

2020年10月28日，甘肃省自然资源厅立项对该滑坡进行治理，具体措施包括：在H3滑坡前部布设29根2 m×3 m的C30砼抗滑桩，中心间距5 m，抗滑桩走向垂直于滑坡滑动方向，局部根据地形调整，同时辅已截排水和局部挡墙加固措施。

2.3 几何模型

为了分析滑坡治理前后边坡的稳定性，选择1-1断面进行分析(图2)，该断面的几何尺寸如图3所示。钻孔数据显示，边坡浅层为厚度变化较大的马兰黄土层，土质疏松、颗粒较均匀，以粉砂为主，呈块状，大孔隙显著，垂直节理发育，一般具有湿陷性。马兰黄土下部为弱风化泥岩层，岩芯呈短柱状，采取率高，干钻困难，遇水易软化，风干易崩解。根据现场勘查和钻孔取样分析结果，勘察单位给出了滑移线的大体位置(图3)。

图2 滑坡平面布置图Fig.2 Landslide layout plan

AB和BC边界为绝热、不可渗透边界；AED为水分可渗透边界图3 1-1断面几何模型及计算域划分Fig.3 Geometric model of the section 1-1 and division of calculation domain

2.4 土层参数

根据原位检测和室内试验结果，土层的物理参数见表3[34-36]，力学参数见表4[37-39]。

表3 土层物理参数Table 3 Physical parameters of the soil layers

表4 土层力学参数Table 4 Mechanical parameters of the soil layers

2.5 边界条件和初始条件

边坡温度边界条件表达式为

T=T0+Asin(2πt+α0)

(12)

式(12)中：α0为初始相位角，由滑坡开始滑动时间决定。年平均地温T0为10.5 ℃，年平均气温变化幅度为14.7 ℃。AED边界的水平位移和竖向位移均不受限制，AB和CD边界的竖向位移不受限制，其他边界及其位移均受限制。

假设2019年8月2日之前边坡的水、热和变形状况稳定，根据边坡初始边界条件对边坡进行稳定性计算，当水热和变形不再发生变化，此时边坡的水热状况将作为模型起算的初始条件。

3 结果和讨论

由于新滑坡是在老滑坡的基础上继续滑动的结果，因此必须对两者同时进行稳定性分析。根据分析结果选择在滑坡前缘布设抗滑桩的方式来控制边坡的变形。

3.1 滑坡处理前边坡稳定性分析

2019年8月2日—2020年8月11日期间该边坡一直处于极慢速滑动状态。图4(a)～图4(d)为2020年8月11日该边坡的温度、水分和变形状况。由图4(a)可见，推测滑移线以上温度梯度较大，推测滑移线以下温度梯度快速减小，温度梯度突变位置与推测滑移线位置接近。这是因为温度梯度作为水分迁移的主要驱动力，其突变位置的水分含量必然会产生明显变化。图4(b)中推测滑移线附近的含水率最大值约为0.37，大于周围土体的含水率，这说明水分在此处发生了集聚(滞水层)。由于含水率越大，泥岩的强度越低，越容易发生滑移破坏。因此，含水率最大的位置与推测滑移线的位置相吻合。图4(c)、图4(d)为边坡垂直变形和水平变形状况，其中最大垂直变形为72.41 cm，最大水平变形为68.57 cm。其中，蓝色线代表0-位移线与地表的距离，垂直变形和水平变形的0-位移线距离地表的距离分别为11.11 m和11.02 m，其中推测滑移线距离地表的距离为10.67 m。另外，根据现场勘查资料显示，滑坡前缘垂直位移大约为0.4 m，图4(c)中滑坡前缘最大垂直位移为0.42 m，与实测结果接近。

图4 2020年8月11日(滑坡处理前)边坡的温度、水分和变形状况Fig.4 Temperature, moisture and deformation of the slope on August 11, 2020

上述结果表明水分和变形计算结果与实际情况相吻合，这也验证了该数值模型的准确性。

需要强调的是，该边坡一直处于蠕变状态，主要是由于水分不均匀分布引起的。这与水分在不同季节的迁移方式不同有关，主要过程包括两个方面：一方面，在冷季，地表温度降低导致下层滞水在温度梯度作用下向活动层内迁移并冻结成冰，造成地表附近体积含冰量和含水量增大[40]；另一方面，在暖季，地表温度上升，活动层内冰晶融化导致体积液态水含量迅速增大。其中小部分液态水通过地气交换排放到大气中，大部分则在温度梯度和基质势作用下向下迁移并在水分迁移驱动力最小的地方发生集聚[35]，导致该处土体的弹性模量、黏聚力和内摩擦角等力学指标大幅降低，加之在上层土体的重力势作用下，势必导致斜坡沿力学强度薄弱位置产生滑移破坏。

3.2 设置抗滑桩后边坡稳定性分析

针对该滑坡的特点，选择在滑坡前缘设置抗滑桩的方式来减小边坡变形。为此，对设置抗滑桩后边坡的变形进行稳定性计算，当边坡达到稳定状态时，其水热和变形状况如图5所示。

图5(a)、图5(b)为设置抗滑桩后边坡的水热和变形状况，与未设置抗滑桩时相比，设置抗滑桩后边坡的温度场和水分场略有不同。首先，边坡整体温度略微升高，最明显的是2 ℃等温线范围有所增加。这是因为黄土的孔隙度较大，导热系数较低；当抗滑桩的阻力有效阻止滑体向下移动时，水分迁移的结果会使得黄土发生湿陷下沉而非滑移。此时黄土的密实度增加导致其导热系数增大，从而使得进入边坡内的热量增加。因此，边坡的温度出现了一定程度的升高。此外，30%含水率等值线高度略有上升，这是因为黄土层密实度增加致使表层土体的温度梯度增大，进而使温度梯度突变点所在高度增加，最终导致暖季水分集聚点位置升高。

图5(c)、图5(d)为设置抗滑桩后边坡的垂直和水平变形，其中，最大垂直变形为2.60 cm，最大水平变形为2.72 cm。与未设置抗滑桩时相比，变形明显减小。由此可见，设置抗滑桩可以有效提升边坡稳定性。此外，垂直变形和水平变形的0-位移线距离地表的距离分别为9.06 m和8.99 m，和图4(c)、图4(d)相比，0-位移线分别提升了2.05 m和2.03 m。由此可见，设置抗滑桩后边坡的0-位移线更高，这说明抗滑桩还可以提升滑移线高度，从而减小了潜在滑体体积，这对提升边坡稳定性有很大帮助。

图5 滑坡治理后边坡的温度、水分和变形状况Fig.5 Temperature, moisture and deformation of the slope after landslide treatment

总的来说，抗滑桩主要是通过增大滑体阻力来提升边坡的力学稳定性[41-42]。除此之外，由于水热和变形之间的相互耦合机制，设置抗滑桩会使得边坡的地温升高、最大含水量所在位置提升，最终导致边坡滑移线升高。这在一定程度上减小了潜在滑体的体积，从而减小了边坡的垂直和水平变形，进而提高了边坡稳定性。

然而，事实证明通过设置抗滑桩的方式对滑坡的治理效果非常有限。毕竟滑移线始终存在，这也是该类型滑坡复活的主要原因。在外界条件(包括暴雨、地震和人为因素)影响下，这种方式难以杜绝滑坡复活。但是，通过详细论述水热变化，尤其是水分不均匀分布对边坡变形的影响，并对牵引式滑坡产生的机理进行了解释，发现可以从减小潜在滑体体积(或者说提升潜在滑移线)的角度提升边坡稳定性，这也为滑坡治理提供了一种新思路，也是下一步要从事的研究。

4 结论

建立了季节冻土区非饱和土的水热力耦合模型，以2020年8月11日甘肃陇南市白冯村发生的大型滑坡为研究对象，对滑坡处理前后边坡中的水热变化和变形进行了分析，讨论了边坡中水分迁移和水分聚集的原因和来源及其对滑坡变形的影响。得到以下主要结论。

(1)详细论述了水热变化，尤其是水分不均匀分布对边坡变形的影响，并对牵引式滑坡产生的机理进行了解释。结果表明，冷季水分向活动层迁移，在冰透镜体作用下水分在活动层下部发生集聚，暖季活动层融化导致液态水大量增加，导致冰透镜体所在位置变成潜在滑移面。在外部荷载(暴雨等)作用下滑体会沿着该滑移面产生滑动，这便是季节冻土区边坡中的水分迁移过程及其对滑体滑动的作用机制。

(2)黄土地区牵引式滑坡复发的主要原因在于滞水层的存在，而季节冻融作用加剧了水分不均匀分布和滞水效应。因此，在外界条件(暴雨和地震等)影响下，原本相对稳定或者已经处于不稳定状态的老滑坡会出现明显滑动。

(3)从力学层面来说，抗滑桩能明显改善滑坡的稳定性。滑坡处理前最大垂直变形和最大水平变形分别为72.41 cm和68.57 cm。在滑坡前缘设置抗滑桩后，边坡达到稳定状态时的最大垂直变形和最大水平变形分别为2.60 cm 和2.72 cm，稳定性显著提升。

(4)抗滑桩只是从力学角度提升了边坡稳定性，但是潜在滑移线仍然存在，在其他外部外荷载作用下，滑坡可能会复活。发现可以从减小潜在滑体体积(或者说提升潜在滑移线)的角度提升边坡稳定性，这也为滑坡治理提供了另一种思路。

此外，目前对滑坡的变形机理研究尚未有比较具体的量化研究结果，尤其是地基中的水分变化非常难确定，从而导致理论研究停滞不前。首次从微观角度基于热力学原理分析了水分变化对季节冻土区滑坡变形的作用机制，这对滑坡变形的防治和量化研究来说具有理论借鉴意义和实践指导价值。