降雨入渗下非饱和土边坡临界稳定性分析

王叶娇，曹 玲，徐永福

(上海交通大学船舶海洋与建筑工程学院，上海 200240)

天然边坡大部分为非饱和土边坡，降雨入渗是诱发非饱和土边坡失稳的重要原因。边坡失稳给社会带来较大的财产损失以及人身安全威胁。目前，国内外学者对非饱和土边坡的降雨入渗机理及其稳定性变化作了很多研究［1－5］。Rahardjo 等［5］对影响降雨过程中非饱和土边坡稳定性的因素及其变化规律作出了研究。Collins等［6］通过计算降雨入渗下边坡内部孔隙水压力的分布，并结合边坡极限平衡分析条分法，提出了边坡稳定性临界深度的计算公式。徐永福［7］研究了分形理论在非饱和土力学中的应用，提出了非饱和土的土水特征曲线与导水系数曲线的分形模型以及非饱和土抗剪强度理论。

本文运用分形模型预测非饱和土的土水特征曲线与导水系数曲线，结合孔隙气流运动，求解孔隙吸力的分布。并且分析了降雨入渗作用下非饱和土边坡体内吸力随入渗深度的变化，根据非饱和土抗剪强度的分形模型，推导非饱和土边坡在降雨入渗下的临界稳定深度公式，结合降雨渗流计算，分析了非饱和土边坡临界稳定性的影响因素。

1 计算模型与参数

1.1 一维非饱和渗流分析与气体流动分析

降雨引起的非饱和土边坡失稳多发生在边坡浅层部位，滑面大致与坡面平行，假设边坡长度与垂直于坡面的深度之比大于10∶1。在此假设下，三维Richard渗流方程可简化成一维渗流模型［8］。降雨入渗在边坡坡面是均匀分布的，取坡体内任一竖直土条进行一维渗流分析，一维渗流模型如图1所示。

在土体顶部建立如图1所示的x－z坐标系，降雨从顶部入渗，设入渗深度为din，在入渗土体部分产生一渗流力S，忽略各土条之间的作用力，设土条在渗流力S，重力W，底部作用力N与TN作用下平衡。当入渗深度din达到一临界入渗深度dcr时，此时土坡的安全系数为1，即非饱和土边坡在降雨入渗下达到临界稳定状态。

一维渗流偏微分控制方程为［9］

式中:Hw为水头;Pa为孔隙气压;kw为y方向上的导水系数;Qw为边界流量;mw为土水特征曲线的斜率;γw为水的重度;t为时间。

考虑降雨入渗下非饱和土体孔隙气体的流动，假设环境温度恒定，气体流动控制方程［10］为

图1 一维渗流计算模型示意图Fig.1 Sketch of the onedimensional infiltration model

式中:ka为气体渗透系数;θa为气体体积含量;ρa为气体密度;γoa与ρoa分别为初始气体重度与初始密度;R为理想气体常数;T为温度;Hw为水头;Pa为孔隙气压。

运用SEEP/W(2007)软件耦合Air/W(2007)软件，求解上述控制方程(1)和方程(2)。

1.2 非饱和土分形模型

非饱和土土体质量分布具有分形特征，根据土体孔隙分形特征与Young-Laplace方程可以得到非饱和土的土水特征曲线与导水系数曲线。

根据分形模型，非饱和土的土水特征曲线方程为［7］

式中:Se是有效饱和度，Se=(θ－θr)/(θs－θr)，θ是体积含水量，θs是饱和体积含水量，θr是残余体积含水量;δ=D－3，D为孔隙分布的分维，对于一般土质，D 介于2.0到3.0之间;ψ 是基质吸力;ψe是进气值。因此，运用孔隙分布的分维和进气值可以预测非饱和土的土水特征曲线。

非饱和土的相对导水系数可表示为

式中:η=3D－11;kr是相对导水系数。若饱和导水系数设为ks，则导水系数k可以表示为

非饱和土的抗剪强度为

式中:τ为非饱和土的抗剪强度;c'为有效黏聚力;(σ －ua)为总正应力;ø'为有效内摩擦角;ξ=D －2。

根据Collins等［6］提出的降雨入渗下非饱和土边坡一维渗流临界深度方程，结合非饱和土的分形模型，临界深度dcr的确定可以运用基质吸力ψ来表示

式中:hp是孔隙正水压;γ与γw分别为土的重度与水的重度;β为边坡坡角。文中不考虑土体抗剪强度的变化，统一取 c'=5 kPa，φ'=15°，β =40°，γ =19.6 kN/m3，γw=10 kN/m3。临界深度曲线与土体吸力分布曲线的交点，即为对应的非饱和土边坡的稳定性临界值，对应的临界深度为边坡发生失稳时的深度值。但临界深度曲线与初始吸力分布值的交点意义仅在于初始时刻非饱和土边坡是稳定的，则该交点在之后的降雨入渗过程中不体现任何意义。

此外，计算孔隙气体渗流运动采用的孔隙气体导水系数与土体的体积含水量的关系如图2所示。在以后的计算中，暂不考虑该类土体性质的变化，计算参数均采用图2所示的数据。

图2 孔隙气体导水系数与体积含水量关系［10］Fig.2 Pore air hydraulic conductivity vs.volumetric water content［10］

2 计算结果

2.1 分形模型参数的影响

2.1.1 分维

取进气值ψe=10 kPa不变，土的分维D分别取为2.1，2.35，2.78。饱和导水系数均取为 1 × 10－5m/s。根据分形理论绘制相应的非饱和土导水系数曲线与土水特征曲线，如图3、图4所示。图中土质类型标为S10，2.1，－5等，其中S是指土质类型，数字分别指进气值 ψe=10 kPa，分维D=2.1，饱和导水系数 ks=1 ×10－5m/s。

图3 不同分维土质的导水系数曲线Fig.3 Hydraulic conductivity functions for soils of different fractal dimensions

计算模型为1 m×4 m的均质非饱和土土条，底部为初始地下水位线，降雨从土条顶部入渗，降雨强度统一设为36 mm/h，孔隙气体初始边界条件设为0 kPa，即与大气压一致，随着降雨入渗，气压逐渐增大至进气值。

SEEP/W(2007)软件与AIR/W(2007)软件耦合计算可以求出各个时刻的孔隙水头Hw与孔隙气压Pa，进而准确求出土体内的吸力值ψ。3种不同分维土质的计算结果如图5所示。

图4 不同分维土质的土水特征曲线Fig.4 Soil-water characteristic curves for soils of different fractal dimensions

图5 降雨入渗下不同分维的土体内部吸力变化Fig.5 Variation of soil suction for soils of different fractal dimensions in the presence of rainfall infiltration

图5 中，标题如 S10，2.1，－5，R36mm/h是指进气值 ψe=10 kPa，分维 D=2.1，饱和导水系数ks=1×10－5m/s的土质类型，R36mm/h是指降雨强度Ir为 36 mm/h。分维 D=2.78的土体(图 5(c))，由于初始时刻土体顶部的导水系数大于D=2.1，2.35的土体的初始导水系数，所以随着降雨入渗，分维较大的土体雨水入渗量较大，土体含水量较快增大，从而导致土中吸力快速下降。另一方面，相同进气值、不同分维土体相应的临界深度曲线形状也不同，从图5中可以看出，分维较高的土体其临界深度曲线斜率较大，且在低分维值的情况下，临界曲线有较明显的拐点。在饱和导水系数与进气值相同的情况下，分维越大，在降雨入渗下非饱和土边坡稳定性下降较快，稳定性临界深度较大。

2.1.2 进气值

取分维D=2.35不变，进气值ψe分别取为1，10，20 kPa 3个值。饱和导水系数均取为1×10－5m/s。图6、图7为3种土质的导水系数曲线与土水特征曲线。

图6 不同进气值土质的导水系数曲线Fig.6 Hydraulic conductivity functions for soils of different air-entry values

图7 不同进气值土质的土水特征曲线Fig.7 Soil-water characteristic curves for soilsof different air-entry values

图8 降雨入渗下不同进气值的土体内部吸力变化Fig.8 Variation of soil suction for soils of different air-entry values in the presence of rainfall infiltration

计算模型保持不变，降雨强度均设置为36 mm/h，且各边界条件与初始条件均保持不变。计算结果如图8所示。不同的进气值的土质在降雨作用下各时刻的吸力分布曲线形态也不同。进气值为1的土条，在降雨达32 h时土体上部产生孔隙正水压(吸力为负值)。临界深度曲线与吸力分布曲线的交点横坐标为负值，即该非饱和土边坡由于孔隙产生正压力而发生失稳。进气值为20的土条在降雨4 h后底部产生正水压力，顶部吸力达到进气值，非饱和区域仅出现在土条中间一小部分。随着进气值的增加，边坡稳定性破坏的降雨入渗时间逐渐减小。主要由于进气值较小的土体顶部导水系数较小，降雨入渗量小且缓慢;而进气值较大的土体，初始时刻土体顶部的导水系数较大，雨水快速入渗，顶部与底部土体均达到饱和，非饱和区域快速减小，非饱和土边坡稳定性也较快下降。

2.2 饱和导水系数与降雨强度的影响

降雨入渗下非饱和土边坡的稳定性与土的饱和导水系数及降雨强度有较大的影响。饱和导水系数决定了土体的入渗能力，降雨强度影响了雨水的入渗量。分别选取 ks=1 ×10－5，1 ×10－6，1 ×10－7m/s 3种不同的饱和导水系数，降雨强度分别取为100，36，15，9，3，0.3 mm/h。

由图9(a)，当降雨强度远大于土体的饱和导水系数时，顶部土体在较短时间内快速饱和，在1.5 h时吸力已接近0，随着降雨继续入渗，雨水入渗深度的增加，在土体底部出现了正的吸力值，这是由于顶部土体快速饱和，封闭了气体，随着孔隙气压的增加，底部土体出现去饱和的现象。而当降雨强度小于饱和导水系数时，如图9(c)所示，顶部土体吸力逐渐减小，但由于降雨量较小，顶部未达到饱和状态，降雨入渗深度也较小。

图9 不同降雨强度下土体内的吸力变化Fig.9 Variation of soil suction in the presence of different rainfall intensities

将饱和导水系数 ks=1 ×10－5，1 ×10－6，1 ×10－7m/s的3类非饱和土边坡出现临界稳定状态的时刻ts与降雨强度的关系曲线绘制于对数坐标系中，如图10所示。曲线呈指数型减小趋势，且不同量级的饱和导水系数对边坡出现失稳所需的时间ts有较大的变化，饱和导水系数ks=1×10－5的边坡一般需数10 h，而ks=1×10－6的边坡在各降雨强度下需10～100 h，ks=1×10－7的边坡需要更长的降雨持时。

图10 不同饱和导水系数的非饱和土坡在降雨入渗下失稳时间与降雨强度的关系Fig.10 The time of soil slope failure with different saturate hydraulic conductivity vs.rainfall intensity

2.3 考虑气体与不考虑气体流动情况下的对比

在降雨入渗非饱和土边坡的过程中，存在孔隙气体的流动，尤其在降雨强度远远大于土体的饱和导水系数时，由于表层土体快速达到饱和，土体内孔隙气体受到封闭，产生气体压力，从而影响孔隙吸力值。

对比图11与图9(a)，未考虑气体流动的土体顶部逐渐由负水压变化至零，而雨水入渗深度以下部分的土体均保持初始状态;若考虑土中孔隙气体的流动，如图9(a)，随着雨水入渗，入渗深度以下部分土体的吸力值逐渐增大，主要由于雨水快速入渗，封闭了一部分气体，孔隙气压逐渐增大，底部饱和土体逐渐开始产生正的吸力值，即出现非饱和区域。

图11 不考虑气体流动时的土体孔隙水压变化Fig.11 Variation of pore water pressure in the absence of air flow

3 结论

(1)运用分形模型预测不同的分维与进气值的非饱和土土水特征曲线及导水系数曲线，并将非饱和土强度分形模型结合非饱和土边坡临界深度的计算，预测降雨入渗作用下非饱和土边坡的稳定性变化。

(2)非饱和土的饱和导水系数与进气值不变时，不同的分维影响降雨入渗速度与入渗深度，分维较大的情况下，雨水入渗速度较快，相应的临界深度曲线斜率较大。当非饱和土的进气值变化时，降雨入渗下土体内不同时刻的吸力变化曲线形态不同。

(3)饱和导水系数与降雨强度对非饱和土边坡在降雨作用下的稳定性影响较大，饱和导水系数直接影响了非饱和土边坡达到临界稳定状态的时间;降雨强度大于饱和导水系数时，表层土体快速到达饱和状态，入渗边界转化成定水头边界，降雨强度较小时，雨水入渗量也较小，入渗深度较浅。

(4)考虑了非饱和土孔隙气体在降雨过程中的流动，分别计算降雨过程中非饱和土边坡体内孔隙气压与孔隙水压，可以准确求出各吸力值。当降雨强度远远大于土体饱和导水系数时，气体受到封闭，产生气压，改变土体的吸力分布。

［1］刘 洋，周 健，王国强.膨胀土基坑边坡降雨入渗的一维数值模拟［J］.同济大学学报，2005，33(2):158－161.(LIU Yang，ZHOU Jian，WANG Guo-qiang.Numerical Method of One-dimension for Rainfall Infiltration in Expansive Soil Slope of Excavation［J］.Journal of Tongji U-niversity，2005，33(2):158 －161.(in Chinese))

［2］RAHARDJO H.Rainfall-Induced Slope Failures［R］.Singapore:Nanyang Technological University，2002.

［3］RAHARDJO H，NIO A S，LEONG E C，et al.Effects of Groundwater Table Position and Soil Properties on Stability of Slope During Rainfall［J］.Journal of Geotechnical and Geoenvironmental Engineering，2010，136(11):1555 －1564.

［4］刘兴旺，戴源夏.基于非饱和土理论的降雨入渗路基稳定性分析［J］.公路交通技术，2009，(3):1 －3.(LIU Xing-wang，DAI Yuan-xia.Analysis of Roadbed Stability after Rain Seepage Based on Unsaturated Soil Theory［J］.Technology of Highway and Transport，2009，(3):1 －3.(in Chinese))

［5］RAHARDJO H，ONG T H，REZAUR R B，et al.Factors Controlling Instability of Homogeneous Soil Slopes under Rainfall［J］.Journal of Geotechnical and Geoenvironmental Engineering，2007，133(12):1532 －1543.

［6］COLLINS B D，ZNIDARCIC D.Stability Analyses of Rainfall Induced Landslides［J］.Journal of Geotechnical and Geoenvironmental Engineering，2004，130(4):362－372.

［7］XU Y F.Fractal Approach to Unsaturated Shear Strength［J］.Journal of Geotechnical and Geoenvironmental Engineering，2004，130(3):264－273.

［8］RICHARD M IVERSON.Landslide Triggering by Rain Infiltration［J］.Water Resources Research，2000，36(7):1897－1910.

［9］Geo-slope International Ltd.Seepage Modeling with SEEP/W2007(Third Edition)［M］.Canada:Geo-slope International Ltd.，2008.

［10］Geo-slope International Ltd.Air Flow Modeling with AIR/W2007(Third Edition)［M］.Canada:Geo-slope International Ltd.，2008.