-
开放科学（资源服务）标识码（OSID）：
-
极端降雨引发的边坡失稳与路基塌陷是岩土工程领域持续关注的重大问题[1-2]，四川汉源山洪灾害和广东梅大高速路基塌陷事件进一步凸显了降雨诱发滑坡机理研究的紧迫性。降雨型滑坡的本质是非饱和渗流场动态演化诱发有效应力场劣化，进而触发渐进式失稳，因此合理量化渗流场的变化是揭示灾变机理的基础性工作。传统非饱和渗流场的计算一般采用非饱和达西方程[3-4]，边界条件根据土表孔隙水压力的正负来选用Dirichlet边界或Neumann边界[5]，但这种简化的处理方法未考虑不同边界之间的转化过程[6]。通过施加“空气单元”或添加“半透水层”可以改进传统方法的不足，实现不同边界之间的转化[7-8]。
值得注意的是，岩土材料显著的地域性和变异特征决定了不同土质的渗流计算理论无法完全借鉴。以我国东南沿海省份广泛分布的花岗岩残积土为例，受物理化学风化作用，其粒度两头多、中间少，呈现出典型的大孔隙特征[9-10]，该介质中以优先流主导的渗流过程普遍采用双重渗透模型来表征[11-12]。窦智等[13]将降雨补给量和实际入渗量作为控制变量，提出了双重渗透模型降雨边界设置方法。Nian等[14]以孔隙水压力作为控制变量，提出双重连续介质降雨边界模型，为方程求解提供了新思路，但有压入渗阶段的积水深度主要根据经验来取值[15]。
在基于渗流-应力耦合计算边坡稳定性的研究方面，尽管极限平衡法与强度折减法被广泛用于获取边坡的整体最小点稳定系数，但无法有效刻画稳定性的空间分布与浅层失稳的时序演化过程。实际上，降雨触发型滑坡多以浅层渐进破坏模式发展为主，Lu等[16]突破传统滑面预设约束提出了局部稳定系数场(LFS)法，特别适用于浅层滑坡动态演化分析；Zhu等[17]通过场分布理论验证了LFS法可定量表征潜在破坏域的空间梯度特征，对揭示边坡灾变的空间分异规律具有重要理论价值。
综上所述，目前坡面积水的动态变化过程在大孔隙流控制方程中未显式表征，并且对考虑坡面积水的大孔隙边坡浅层失稳机制认识不清晰，基于此，本文提出考虑径流-入渗耦合作用的大孔隙边坡点稳定系数计算方法，揭示不同典型参数对边坡浅层失稳机制的影响，以期进一步丰富残积土边坡灾害治理的理论体系。
HTML
-
传统的均质土(无大孔隙)非饱和流控制方程一般采用Richards模型描述，在此基础上，大孔隙流控制方程采用双重渗透模型描述[18]，其核心思想是将土壤中复杂的孔隙网络概化为2个相对独立的系统，分别为基质域和大孔隙域，且两域之间通过水分交换项耦合，如式(1)和式(2)所示：
其中：Cm、Cf分别为基质域、大孔隙域比水容量(m-1)；hm、hf分别为基质域、大孔隙域压力水头(m)；Km、Kf分别为基质域、大孔隙域非饱和导水系数(m/s)；ωm、ωf分别为基质域、大孔隙域体积占比，满足ωm+ωf=1；Gw为水分交换项(s-1)；t为时间(s)。
上述参数Cm、Cf、Km、Kf是与hm、hf分别相关的变量，具体关系参考VG(Van Genuchten)模型[12]。
水分交换项以两域水头差作为驱动力，如式(3)所示：
式中：αw为一阶水分交换系数(m-2)。
-
采用Kinematic wave方程描述坡面径流过程[19]，如式(4)所示：
式中：H为坡面积水深度(m)；q为坡面单宽流量(m2/s)；x为沿着坡面向下的坐标轴；R为降雨强度(m/s)；α为坡角(°)；i为入渗率(m/s)；n为曼宁粗糙系数(m1/3/s)，取0.035 m1/3/s。
-
点稳定系数最早由Lu等[16]提出，与传统边坡稳定系数计算方法相比，点稳定系数法突破传统滑面几何预设的约束，可以刻画降雨入渗过程中应力场的空间梯度，其数学定义表征如式(5)所示[20]：
式中：LFS为点稳定系数；c′为有效黏聚力(kPa)；φ′为有效内摩擦角(°)；σ′1、σ′3分别为第一、第三有效主应力(kPa)。
1.1. 大孔隙流控制方程
1.2. 坡面径流控制方程
1.3. 点稳定系数方程
-
对于双重渗透模型，总降雨强度R与两域等效降雨强度Rm、Rf之间的关系如式(6)所示：
式中：Rm、Rf分别为基质域、大孔隙域等效降雨强度(m/s)。
总入渗率i与两域实际入渗率im、if之间的关系如式(7)所示：
式中：im、if分别为基质域、大孔隙域入渗率(m/s)。
土壤入渗能力随着降雨不同历程而发生变化，根据降雨强度、土壤入渗能力、饱和渗透系数等参数将降雨过程划分为以下3个阶段。
1) 降雨初始阶段，土壤表层含水率较低，入渗能力大于降雨强度；雨水完全入渗，两域入渗率等于降雨强度，如式(8)所示：
此时，两域边界控制方程如式(9)所示：
式中：ρw为水的密度(kg/m3)。
2) 当降雨历时达到基质域超渗产流临界状态时，两域流体动态重分布机制导致坡面径流向大孔隙域补给。此时基质域入渗率受控于土壤入渗能力，基于Green-Ampt入渗理论框架[17]，基质域边界控制方程如式(10)所示：
式中：icm为基质域入渗能力(m/s)；Fm为单位面积的入渗量(m)；sm为基质域湿润锋处平均吸力水头(m)，取0.2 m；Δθm为基质域饱和含水率与初始含水率之差(%)；Ksm为基质域饱和导水系数(m/s)。
大孔隙域等效降雨强度如式(11)所示[13-14]，边界控制方程如式(12)所示：
3) 当降雨历时进一步增大，大孔隙域无法接收全部降雨和径流补给量，坡面产生积水，此时方程边界条件为有压边界，如式(13)和式(14)所示：
其中：icm-pond为积水状态下基质域入渗能力(m/s)，相比式(10)中的icm，计算icm-pond时需考虑坡面积水深度；if-pond为积水状态下大孔隙域入渗率(m/s)；Ksf为大孔隙域饱和导水系数(m/s)。
上述边界条件通过图 1所示流程表示。
-
由于边界条件处于动态变化过程中，需根据降雨强度、两域饱和渗透系数等参数来确定，压力边界通过径流-入渗耦合作用来实现。基于COMSOL Multiphysics多物理耦合求解器，并通过二次开发自适应边界算法，实现径流-入渗耦合作用的压力边界自迭代修正。数值模型中由上至下设置3个域点探针(1#、2#、3#)用于监测体积含水率，基于室内模型试验验证数值模型的有效性，如图 2所示。
监测点体积含水率时变响应特征如图 3所示。由图 3a可知，基质域体积含水率模拟值和试验值变化趋势具有显著一致性，值得注意的是，由于试验过程大孔隙长度只有60 cm，3#测点大约在3.75 h后试验值略低于模拟值。由图 3b可知，在大孔隙域内，两者变化曲线也保持较高的吻合度，但1#测点存在微小偏差，可能原因是大孔隙通道上部石英砂填充密实度不足导致孔隙增大，在非饱和渗流过程中持水能力下降[21]，导水性增强，水分流经该区域快速补给下部土体使体积含水率减小[22]，但总体偏差小于8%。综上所述，室内试验结果与数值计算结果可以互相验证。
-
理论上，降雨作用下渗流场与应力场是一个双向耦合的关系。需要说明的是，本文点稳定系数求解方法源自线弹性本构框架[16]，不适用于塑性及以上阶段，鉴于线弹性阶段土壤介质变形量级处于可忽略范围，因此仅考虑水分变化对应力的影响。根据式(5)可知，求解点稳定系数需要先求得土体任一点的应力，具体求解流程如图 4所示。
文献[23]对比分析了点稳定系数法与传统极限平衡法，结果发现，高边坡工况下点稳定系数法与传统极限平衡法的计算结果具有良好的一致性；而在缓坡地形，前者在分析边坡局部失稳状态上优势明显，而后者忽略应力场梯度特征导致安全系数预测趋于保守，且无法捕捉潜在的局部破坏区域。这本质源于点稳定系数法突破了滑面几何预设约束，建立了应力场驱动的稳定性连续场表征机制，故在浅层破坏分析中更具理论优势。
2.1. 考虑坡面积水动态变化的双重渗透模型边界条件
2.2. 大孔隙流控制方程求解
2.3. 点稳定系数方程求解
-
以宁武高速公路福建政和段某路基典型残积土边坡为地质原型，坡高为20 m，坡比为1∶1.5~1∶2。该边坡于2021年6月汛期受持续强降雨影响触发多点群发性浅层滑坡，如图 5所示。根据气象统计，数值模拟采用重现期匹配的降雨强度极值，分别取1.5 cm/h和6 cm/h。边坡几何模型概化后如图 6所示，为了提高模拟结果的精度，将土层由上往下设置加密区网格和细化区网格。上部采用降雨入口边界，两侧和底部为无流量边界。数值计算中，大孔隙占比ωf取0.05，一阶水分交换系数αw取4×10-5 cm-2，吸力水头平均值sm取20 cm，其余计算参数如表 1所示。
-
均质土边坡(不考虑大孔隙)与大孔隙边坡点稳定系数对比如图 7、图 8所示，将点稳定系数取值1.2作为失稳区划分界限。根据图 7a、图 7b可知，低雨强工况(1.5 cm/h)下，2种边坡均发生浅层失稳，其中均质土边坡最小点稳定系数约为0.94，大孔隙边坡约为1.07。根据图 7c、图 7d可知，高雨强工况(6.0 cm/h)下，均质土边坡最小点稳定系数减小至0.92，大孔隙边坡最小点稳定系数减小至0.86，表明大孔隙结构在极端降雨下显著加剧灾变风险。需特别指出，云图显示坡体深部区域稳定性的异常下降源于本构模型适用的边界问题，现行理论基于线弹性假设，而深部土体在较大应力作用下进入塑性屈服阶段，其稳定状态不再适用点稳定系数法来判断，后文中云图深部土体的稳定性异常亦源于此。
由图 8a可知，低雨强工况(1.5 cm/h)下，大孔隙边坡失稳区面积始终小于均质土边坡，两者最大相差3倍。由图 8b可知，高雨强工况(6.0 cm/h)下，大孔隙边坡的失稳区面积在40 h后开始快速增大，并于44 h后超过均质土边坡。此后，由于降雨快速入渗导致土体抗剪强度大幅下降，两者失稳区面积差距逐渐扩大，最终大孔隙边坡失稳区面积达到均质土边坡的2.4倍。
图 7和图 8综合表明大孔隙对边坡灾变演化具有显著的雨强-历时双阈值效应，均质土边坡稳定性对雨强的敏感性显著低于大孔隙边坡。降雨历时越大，大孔隙对边坡稳定性的不利作用相比均质土边坡更明显，其原因是降雨强度或降雨历时增大到一定程度，来水无法完全流入基质域，大孔隙流发挥作用，触发更深土层的孔隙水压力，孔压快速响应引起有效应力减小，这也解释了图 5案例中持续性降雨诱发大孔隙边坡发生多点群发性浅层失稳的原因。由此表明，对于大孔隙边坡，本文提出的计算方法适合降雨强度较大或降雨历时较长的情形，若采取传统均质土边坡稳定性的计算方法，计算结果则偏保守。
3.1. 数值模型
3.2. 对比分析
-
大孔隙占比ωf对点稳定系数的影响如图 9-图 11所示。由图 9可知，ωf取值为0.05、0.10、0.20时，最小点稳定系数分别为0.87、0.76、0.61，点稳定系数显著减小；最大失稳层深度分别为50、80、150 cm。同时还可以发现，ωf由0.05增加到0.20时，失稳区域进一步扩展，尤其在坡顶和坡趾附近向外延伸约6 m，此时边坡土表出现较大范围的浅层失稳破坏。
图 10是在图 9基础上对点稳定系数空间变化关系更详细的描述。由图 10a可知，坡脚剖面点稳定系数随着深度增大而线性增大，ωf取值为0.05、0.10、0.20时，失稳层深度分别为3、19、47 cm，最小点稳定系数减小了28.1%，ωf对边坡浅层稳定性具有显著的负面作用。由图 10b可知，ωf取值为0.05、0.10、0.20时，失稳层深度分别为21、49、87 cm，最小点稳定系数为0.88、0.78、0.63，说明坡中剖面灾变敏感性相比坡脚更为显著。
由图 11可知，失稳区面积曲线呈“阶梯型”变化，不同ωf取值下增长速率具有一定差异，ωf取值为0.05和0.10时，失稳区面积增长速率先小后大，当ωf由0.05增至0.20时，最终失稳区面积增长了128.0%。上述分析均表明，增大ωf对边坡稳定具有显著的负面影响，其原因是增大ωf导致大孔隙域深部土壤孔隙水压力快速响应，基质吸力迅速下降使得边坡失稳时间提前和失稳区面积变大。同时，增大ωf使两域水头差变大，水分交换加强增加了基质域入渗深度，边坡土体抗剪强度减小导致浅层稳定性降低。
-
水分交换系数αw对点稳定系数的影响如图 12-图 14所示。由图 12可知，不同αw取值下，滑裂面与坡面大致平行，其深度随αw增大而略微增大。αw由4×10-5 cm-2增至2.5×10-2 cm-2时，失稳深度由80 cm增至100 cm，最小点稳定系数由0.76减小至0.73，表明增大αw在一定程度上不利于大孔隙边坡稳定。
由图 13可知，αw分别取值为4×10-5、1×10-3、2.5×10-2 cm-2时，坡脚表层点稳定系数分别为0.85、0.84、0.82，尽管坡脚与坡中剖面的失稳层深度随αw增大略有增加，但增幅有限，这表明αw对点稳定系数的空间分布影响相对较弱。对比图 13a、图 13b可知，相同αw取值下，坡中剖面失稳层深度明显大于坡脚剖面，表明大孔隙边坡坡中处的稳定性更差。
由图 14可知，不同αw取值下失稳区面积曲线在40~60 h时间范围内区分度较明显，且随着αw的增大，失稳区面积呈现增大的趋势。αw分别为4×10-5、1×10-3、2.5×10-2 cm-2时，失稳区面积最终依次为32.09、32.64、33.01 m2，最大增幅为2.9%，表明失稳区面积与αw呈现一定的正相关性。这是因为降雨至一定时间后，大孔隙域湿润锋深度超过基质域，水分交换总体上从大孔隙域向基质域迁移。因此，水分交换系数越大，意味着基质域湿润锋深度也越大，体积含水率也随之升高。由于导入应力场的水分场由两域体积含水率加权求和得出，基质域权重显著大于大孔隙域，因此，基质域体积含水率的升高使整体含水率增大，导致土体有效应力减小，从而削弱边坡的稳定性。
-
导水系数比μ对点稳定系数的影响如图 15-图 17所示。由图 15可知，潜在失稳深度随μ增大呈现正相关增长，而最小安全系数则表现出显著负相关特性。当μ从20增至100时，失稳深度由70 cm延伸至90 cm，失稳区域在纵深和水平方向同步扩展，其中坡脚失稳区外延长度从1.5 m扩展至4.5 m。这表明当两域导水性能差异显著时，需重视坡脚处的失稳破坏现象。
由图 16可知，当μ从20增至100时，坡体中部剖面表层点稳定系数从0.83降至0.74，降幅达10.8%；而在50 cm深度处安全系数由1.06微降至1.02，仅减少3.8%。进一步分析图 16b中破坏面发育特征，μ取值为20、50、100时对应的临界滑动面深度分别为38、44、48 cm，呈现显著的正相关增长趋势，这表明μ增大对边坡整体稳定性构成不利影响。
由图 17可知，当μ取20和50时，降雨初期边坡失稳区域扩展规律呈现高度相似性，而μ为100时失稳区域面积显著大于前2种工况。整体观察发现，大孔隙边坡失稳面积呈现典型的“阶梯型”增长特征，且随μ提升而逐级扩大。当μ从20增至100时，破坏区域面积增幅达3.7%。这种演变机制源于：μ增大意味着大孔隙域渗透性能增强而基质域渗透能力减弱，导致大孔隙域入渗通量增速及增量显著超过基质域，不仅极大地加速了优先流现象的产生，还加快了两域间的水分交换频率，使得土体整体含水率升高，进而对边坡稳定状态产生不利影响。
4.1. 大孔隙占比对点稳定系数的影响
4.2. 水分交换系数对点稳定系数的影响
4.3. 导水系数比对点稳定系数的影响
-
1) 降雨强度为1.5 cm/h时，大孔隙边坡失稳区面积小于均质土边坡，总体稳定性相比均质土边坡更好；降雨强度为6.0 cm/h时，大孔隙边坡最大失稳层深度、失稳区面积大约是均质土边坡的6倍、2.4倍，整体稳定性更差。大孔隙边坡在强降雨或持续性降雨条件下更易发生浅层失稳破坏。
2) 不同大孔隙占比ωf下，边坡均表现为浅层失稳破坏。ωf越大，土壤深部孔隙水压力响应越显著，致使失稳层深度增加，失稳破坏区域由坡顶和坡脚向外扩展，失稳区面积曲线呈“阶梯型”增长态势，ωf由0.05增加到0.20时，最小点稳定系数减小了28.1%，失稳区面积增加了128.0%。
3) 水分交换系数αw对大孔隙边坡点稳定系数的影响相比大孔隙占比更小。随着αw增大，最大失稳层深度增加了25.0%，不同αw取值下失稳区面积在降雨至40~60 h范围内差异最为显著，最终失稳区面积增加了2.9%。
4) 随着两域导水系数比μ增大，大孔隙边坡失稳层深度增大，最小点稳定系数减小，失稳区范围由坡脚向外延伸，面积增加了3.7%。当大孔隙域与基质域导水系数之比较大时，需重视坡脚处的失稳破坏现象。
DownLoad: