非饱和土边坡饱和区降雨强度-时间临界曲线

王述红， 何 坚， 刘 欢， 韩文帅

(东北大学 资源与土木工程学院， 辽宁 沈阳 110819)

山体滑坡造成全球巨大的生命财产损失，降雨是引起滑坡最重要的自然因素.降雨型滑坡主要以浅层滑坡为主，降雨引起坡面浅层含水率迅速升高，形成的饱和区将软化和削弱土体强度，诱发浅层滑坡[1].坡面浅层产生的饱和区是导致斜坡浅层破坏的主要原因[2].随着数值计算技术的发展，边坡稳定性数值演化分析得以推动[3].

降雨形成饱和区过程复杂，特别在考虑边坡入渗过程中的边界效应[4]，即在底部边界排水较差时，较小的降雨强度，持续时间足够长，底部水位抬升，斜坡浅层仍能形成饱和区.饱和区的成形涉及到降雨强度、土体渗透系数的动态变化，目前缺少获取饱和区临界降雨条件的有效方法.

Mein等研究发现降雨过程中土体表层渗透系数与降雨强度存在定量关系[5]，基于此，提出了以土体饱和度表征降雨强度的计算方法.并以土体表层饱和度首次达到饱和度0.9后，表层饱和度持续稳定在0.9±0.005内，为土体表面产生饱和区的临界判据，数值计算出了饱和区产生的临界降雨强度与时间.通过创新性构建的饱和区降雨强度-时间临界曲线模型，实现了对饱和区产生的临界降雨阈值的判定.

1 非饱和土渗流基本理论

非饱和土中渗流运动可使用Richards方程进行描述，通过质量守恒定律获得适用于非饱和土瞬态流的基本控制方程[6]:

(1)

式中：H为总水头；kx为x方向的渗透系数；ky为y方向的渗透系数；Q为施加在边界的流量；mw为土水特征曲线的线性段斜率；rw为水容重；t为时间.

非饱和土的渗透系数与土壤含水率有关，通过SWCC(土水特征曲线)可间接求解.经过大量研究和实践发现，SWCC和土体渗透系数能被Mualem-van Genuchten(MVG)模型准确表征，SWCC的表达式为[7]

(2)

式中:θ为体积含水率；h为压力水头；θs为饱和含水率；θr为残余含水率；α和n为土水特征参数；m=1-1/n.

渗透系数的表达式为[7]

(3)

式中，Ks为土体饱和时的渗透系数.MVG模型为Richards基本控制方程的求解提供必要条件.

2 饱和区临界条件分析方法

2.1 饱和区临界条件分析的基本原理

含不透水基岩的土体纵向含水率剖面如图1所示.根据Coleman等提出的均质土层入渗的含水率剖面图，在入渗过程中仅表面浅层土体才能达到饱和，其余土体介于饱和与非饱和之间[8]，土体饱和度0.9以上区域近似为饱和区[9].

典型非饱和土入渗过程如图2所示.当降雨强度q≤Ks，Mein等发现土体表面入渗率K0与降雨强度数值相等，即土体表面的渗透系数等于降雨强度[5].因此，以土体饱和(饱和度0.9)对应渗透系数作为降雨强度(渗透系数与降雨强度同量纲，均为m/s)，一定降雨时间后，土体表面饱和度将趋近饱和值并保持稳定，此时的降雨强度为临界降雨强度.当降雨强度低于临界降雨强度时，土表面渗透系数将低于临界降雨强度值，即表面含水率低于对应饱和含水率，土表面将无饱和区产生(不考虑渗流边界效应).

在下部排水条件较差的情况下(图1)，由于湿润峰极易下渗到不透水边界，导致雨水在底部蓄积并抬升地下水位，产生渗流边界效应.在较小的降雨强度下，持续足够长时间，坡面仍能产生饱和区.因此，饱和区的临界降雨强度仅在湿润峰未抵达底部不透水边界时才存在，且有意义.

2.2 土体饱和度表征降雨强度的计算方法

首先假设:①降雨强度小于土体饱和渗透系数；②湿润峰未下渗到底部不透水边界或地下水位处，即无渗流边界效应；③土体是均质的.将降雨强度以土体饱和度所对应的渗透系数的方式进行表示，使土体表层饱和度到达预设数值.考虑到残余含水率，这里使用有效饱和度Se表征土体含水程度:

(4)

式中，S为饱和度；Sr为残余饱和度.在初始入渗阶段，土体均为非饱和土，h<0.根据MVG数学模型，联立式(2)～式(4)获得饱和度与渗透系数关系:

(5)

根据降雨入渗规律，一定入渗时长后，表层饱和度所对应的渗透系数与降雨强度相等，降雨强度在数值上等于浅层土体渗透系数，即

(6)

进一步考虑土体坡度β，由式(6)获得边坡降雨强度与表层渗透系数的关系:

(7)

2.3 饱和区降雨强度-时间临界曲线模型

对饱和区临界曲线模型提出以下假设条件:①非饱和土为均质土；②降雨强度小于土体饱和渗透系数(同量纲m/s)；③降雨过程中，湿润峰未到达底部不透水边界，无渗流边界效应.饱和区降雨强度-时间临界曲线如图3所示，降雨强度为纵坐标，对应的降雨时间为横坐标，利用临界曲线将区域划分为不产生饱和区与产生饱和区.

其中，饱和区临界曲线中存在一临界点，临界点对应临界降雨强度，即在降雨强度低于临界降雨强度时，若不考虑渗流边界效应影响，低于该降雨强度的条件下，土体浅层无饱和区产生.因此，采用分段函数的数学模型表达饱和区临界曲线:

(8)

式中：φ与ε分别为拟合参数；q的单位为m/s，在曲线段使用指数函数并引入一个修正系数γ对曲线进行修正.

3 土柱饱和区临界曲线

3.1 土柱模型建立

土柱模型如图4所示，底部为地下水位线，在深度0，0.3，0.8，1.5，2.0 m处设置监控点，模型高度7 m，顶部为入渗边界，其余各边界为不透水边界.采用易产生饱和区的粉质黏土，土水参数见表1.使用ABAQUS软件模拟，采用CPE4P单元，单元长度为0.1 m.首先进行稳态土应力平衡分析获得初始含水分布，后续使用瞬态分析步做四组降雨分析，根据式(6)以预设的饱和度设置对应降雨强度，见表2.

表1 粉质黏土土水参数

3.2 运算结果分析

降雨36 h后，湿润峰最大下移深度为1.4 m，如图5所示.方案中湿润峰均未达到底部不透水边界，未产生渗流边界效应.边坡浅层最先形成饱和区，方案中饱和区最大深度为1.105 m，饱和区内土层整体饱和度与土体表面饱和度相同.

表2 土柱模型计算方案

降雨强度小于土体饱和度0.93对应的降雨强度(9.474×10-7m/s)如图6所示，计算的土体表面饱和度与预设饱和度的平均误差较大，为-0.013，计算的土体表面渗透系数平均误差为-1.10×10-7m/s；在降雨强度大于9.474×10-7m/s时，计算饱和度与预设饱和度平均误差为-0.003，计算的渗透系数与预设的降雨强度平均误差为-2.08×10-7m/s.

表层土体达到饱和后，雨水以恒定速度渗入土内，饱和区内的饱和度与土体表面饱和度相同，而该饱和度所对应的渗透系数与降雨强度一致.采用预设饱和度对应的渗透系数值作为降雨强度进行数值计算，计算的土体表面饱和度和渗透系数对比方案预设值的误差最大为-1.74%.

通过模型表面监控点，获得土体表面饱和度随时间的变化曲线如图7所示，饱和度在降雨开始时迅速增大，并达到稳定.在饱和度0.915对应渗透系数的降雨强度下，土体表层饱和度达到0.9后，稳定在0.9015，在0.9±0.002以内，认为该降雨强度为产生饱和区的临界降雨强度，此时对应的降雨时间为临界降雨时间.

获取各降雨强度及对应的降雨时间，结合临界降雨强度与降雨时间，进行数据拟合获得一维入渗饱和区降雨强度-时间临界曲线，见图8.在降雨强度大于临界降雨强度时，粉质黏土产生饱和区耗时急剧减少.

4 边坡饱和区临界曲线

4.1 边坡模型及计算方案

为考虑坡度对饱和区临界曲线的影响，建立边坡模型，如图9所示.以坡面中部监控点的饱和度判断浅层坡面饱和区的产生情况.

具体为:ef与gh平行，eg与fh平行；ef为入渗边界，其余各边界为不透水边界，土柱入渗中湿润峰最大入渗深度为1.105m，边坡模型3m土厚有效排除渗流边界效应；采用ABAQUS软件进行数值模拟，网格单元为CPE4P，长度为0.2 m；首先稳态平衡分析获得初始含水分布，稳定后改变边界gh的初始孔隙水压力，使边坡表面饱和度在降雨前与土柱模型顶部保持相同表面饱和度，饱和度误差为±0.005；采取垂直降雨的模式，垂直坡面入渗的降雨作为有效降雨强度，平行坡面的降雨则沿坡面流出.由式(7)获得预设饱和度所对应的实际降雨强度，由于临界降雨强度产生饱和区时间差异大，对临界降雨强度计算时间会加长，见表3.

表3 边坡模型计算方案

4.2 坡度对饱和区临界曲线的影响

设置不同坡度，并保持坡面相同初始饱和度，探究坡度和垂直于坡面入渗的有效降雨强度，对饱和区有效降雨强度-时间临界曲线的影响见表4.

50°和60°边坡的有效临界曲线在A点相交，如图10所示.当qe>1×10-6m/s，60°边坡产生饱和区速度比50°边坡快，而当qe<1×10-6m/s，则反之.40°和60°边坡的有效临界曲线在B点相交，当qe<7.78×10-7m/s，40°边坡产生饱和区耗时比60°边坡要少，依此地，30°和60°边坡的有效临界曲线在C点相交(即qe<7.25×10-7m/s)，同时在C点后坡度越小饱和区产生时间越短.

表4 饱和区临界曲线拟合方程

50°边坡有效临界曲线依此在D点和E点，分别与40°，30°边坡有效临界曲线相交，当降雨强度小于交点处的强度时，坡度较小边坡产生饱和区时间越短.40°与30°的有效临界曲线始终近似平行，qe一定时，坡度越大，饱和区产生时间越短.

表4中所有坡面初始饱和度相差较小，坡面基质吸力可认为相同，因此主要分析重力对饱和区产生时间的影响.坡度较大时，雨水受重力沿坡面向下分力越大，但qe较大时，后续入渗雨水迅速补给，重力沿坡面分力使雨水沿坡面下渗的影响相对迅速，垂直坡面入渗的雨水可忽略.同时，50°与60°边坡内雨水所受重力垂直坡面分力较30°与40°边坡小，相同qe，坡度越大，雨水垂直下渗速度越慢，坡面易聚水产生饱和区.综上，qe较大时，坡度越大，坡面浅层产生饱和区越快.

随qe减小，50°与60°边坡内渗入雨水受重力沿坡面向下分力的作用突显.当qe<1×10-6m/s时，60°边坡雨水受重力沿坡面向下分力影响较大，雨水沿坡面下渗较大，60°边坡产生饱和区速度较50°边坡慢.然而30°与40°边坡入渗雨水受重力沿坡面向下分力十分接近，因此30°与40°边坡中入渗雨水的重力效应无法在曲线中体现.

4.3 边坡饱和区临界曲线运用

4.3.1 诱发浅层滑坡的坡度预测

通过国家气象局制定的降雨强度等级，对临界曲线上的实际降雨强度进行划分，如图11所示.相同降雨强度，坡度越大产生饱和区耗时越长.降雨强度越大，50°与60°边坡产生饱和区时间差越小，饱和区形成时间与30°和40°边坡逐渐趋近.

不同坡度边坡的饱和区临界曲线反映出降雨过程中边坡产生饱和区存在“危险坡度”.30°与40°边坡在不同降雨强度下产生饱和区耗时均较少，能迅速产生饱和区.根据大量数据统计，香港地区的天然场地发生滑坡的坡角大都在35°～40°之间，重庆地区滑坡中30°～40°的斜滑坡最为发育[11].因此，结合不同坡度边坡的饱和区临界曲线，推测潜在浅层滑坡的危险坡度为30°～40°，并与实际统计数据有较好一致性.

4.3.2 诱发浅层滑坡的降雨强度预测

在图11的暴雨等级中，各坡度产生饱和区时间长短不一致，很难形成大量的边坡饱和区，大范围边坡滑坡的潜在危险较小.而降雨强度大于147 mm/d(大暴雨级)时，30°～50°边坡迅速在6.73 h内形成饱和区，具有大范围发生的潜在风险.因此，多种坡度的饱和区临界曲线结合当地气象统计资料，可推测区域边坡发生大范围饱和区的危险降雨强度范围，应用潜力巨大.

5 结 论

1) 使用土体饱和度表征降雨强度的方法计算土体表面饱和度，对比预设值的平均误差为-0.008.

2) 40°与30°边坡在相同有效降雨强度作用下，坡度越大，坡面产生饱和区速度越快；坡度50°与60°的边坡，有效降雨强度较大时，60°的边坡产生饱和区速度比坡度50°，40°，30°的边坡要快；随着有效降雨强度的降低，60°与50°的边坡产生饱和区速度会依此低于40°与30°边坡.

3) 相同实际降雨强度，30°～40°的边坡形成饱和区速度最快，形成浅层滑坡风险较大.

4) 通过多种坡度的饱和区降雨强度-时间临界曲线，可获得区域内边坡发生潜在大量饱和区的危险降雨强度范围.