东北大学学报:自然科学版  2020, Vol. 41 Issue (10): 1452-1458  
0

引用本文 [复制中英文]

王述红, 何坚, 刘欢, 韩文帅. 非饱和土边坡饱和区降雨强度-时间临界曲线[J]. 东北大学学报:自然科学版, 2020, 41(10): 1452-1458.
[复制中文]
WANG Shu-hong, HE Jian, LIU Huan, HAN Wen-shuai. Critical Curve of Rainfall Intensity-Time for Saturated-Zone in Unsaturated Soil Slope[J]. Journal of Northeastern University Nature Science, 2020, 41(10): 1452-1458. DOI: 10.12068/j.issn.1005-3026.2020.10.014.
[复制英文]

基金项目

国家自然科学基金资助项目(51474050，U1602232)；中央高校基本科研业务费专项资金资助项目(N170108029)；辽宁省自然科学基金资助项目(20170540304，20170520341)；中建股份科技研发课题(CSCEC-2016-Z-20-8)；辽宁省科学技术计划重大项目(2019JH2/10100035)

作者简介

王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。

文章历史

收稿日期：2020-01-02
Contents              Abstract              Full text              Figures/Tables              PDF
非饱和土边坡饱和区降雨强度-时间临界曲线
王述红 , 何坚 , 刘欢 , 韩文帅     
东北大学 资源与土木工程学院，辽宁 沈阳 110819
收稿日期：2020-01-02
基金项目：国家自然科学基金资助项目(51474050，U1602232)；中央高校基本科研业务费专项资金资助项目(N170108029)；辽宁省自然科学基金资助项目(20170540304，20170520341)；中建股份科技研发课题(CSCEC-2016-Z-20-8)；辽宁省科学技术计划重大项目(2019JH2/10100035)。
作者简介：王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。
摘要：为研究降雨入渗边坡产生饱和区的临界条件，使用基于MVG模型的土体饱和度表征降雨强度的计算方程.当土表面饱和度首次达到0.9以上时，饱和度持续稳定在0.9±0.005为饱和区出现的临界判据，数值计算得出土柱和边坡的饱和区临界降雨强度与降雨时间.通过数据拟合构建了饱和区降雨强度-时间临界曲线模型，分析了不同坡度、有效降雨强度对临界曲线的影响.结果表明:饱和度表征降雨强度的计算方法可使土体表面饱和度达到预定值，平均精度为-0.008；有效降雨强度较大时，坡度为50°~60°的边坡，坡度越大产生饱和区的耗时越短，反之耗时越长；坡度30°~40°为易发浅层滑坡的“危险坡度”.
关键词：降雨入渗    饱和区    坡度    临界降雨强度    Mualem-van Genuchten模型    
Critical Curve of Rainfall Intensity-Time for Saturated-Zone in Unsaturated Soil Slope
WANG Shu-hong , HE Jian , LIU Huan , HAN Wen-shuai     
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: HE Jian, E-mail: hejian95@foxmail.com.
Abstract: In order to study the critical condition of rainfall-induced saturated-zone forming, the soil saturation based on the MVG model was used to characterize the calculation equations of rainfall intensity. When the surface saturation of the soil reaches above 0.9 for the first time, the saturation is continuously stable at 0.9 ± 0.005 as the critical criterion of the saturated-zone. The numerical calculation method was used to obtain the critical rainfall intensity and rainfall time in the saturated-zone of the soil column and slope. Through the data fitting, a rainfall intensity-time critical curve model of the saturated zone was constructed, and the influence of different slope and effective rainfall intensity on the critical curve was analyzed. The results showed that the calculation method of saturation characterizing rainfall intensity can make the soil surface saturation reach a predetermined value, with an average accuracy of -0.008. When the effective rainfall intensity is large, the greater the slope, the shorter the time consuming to produce the saturated zone for the slope of 50°~60°, on the contrary, the longer it takes. A slope of 30°~40° is a "dangerous slope" prone to shallow landslides.
Key words: rainfall infiltration    saturated-zone    slope angle    critical rainfall intensity    Mualem-van Genuchten model    

山体滑坡造成全球巨大的生命财产损失，降雨是引起滑坡最重要的自然因素.降雨型滑坡主要以浅层滑坡为主，降雨引起坡面浅层含水率迅速升高，形成的饱和区将软化和削弱土体强度，诱发浅层滑坡[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].

图 1 含不透水基岩的土体纵向含水率剖面 Fig.1 Longitudinal water content profile of soil with impervious bedrock

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

图 2 典型非饱和土入渗过程[5] Fig.2 Typical infiltration process in unsaturated soil[5]

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

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

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

(4)

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

(5)

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

(6)

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

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

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

图 3 饱和区降雨强度-时间临界曲线 Fig.3 The critical curve of rainfall intensity and time for the saturated-zone

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

(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.

图 4 土柱模型 Fig.4 Soil column model
表 1 粉质黏土土水参数 Table 1 Soil-water characteristic parameters of silty clay
表 2 土柱模型计算方案 Table 2 Soil column calculation scheme
3.2 运算结果分析

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

图 5 降雨36 h粉质黏土饱和度随深度变化 Fig.5 Saturation changes of silty clay with depth after 36 h rainfall

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

图 6 数值计算与方案预设值对比曲线 Fig.6 Contrast curves of numerical result and scheduled values (a)—饱和度；(b)—渗透系数.

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

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

图 7 降雨36 h粉质黏土表层饱和度随时间变化 Fig.7 Surface saturation of silty clay changes with time after 36 h rainfall

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

图 8 一维入渗饱和区降雨强度-时间临界曲线 Fig.8 One-dimensional infiltration saturated zone rainfall intensity-time critical curve
4 边坡饱和区临界曲线 4.1 边坡模型及计算方案

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

图 9 边坡模型 Fig.9 Slope model

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

表 3 边坡模型计算方案 Table 3 Calculation scheme for slope model
4.2 坡度对饱和区临界曲线的影响

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

表 4 饱和区临界曲线拟合方程 Table 4 Fitted equation of the critical curves for the saturated-zone

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

图 10 不同坡度边坡的饱和区有效降雨强度-时间临界曲线 Fig.10 The critical curves of effective rainfall intensity and time for the saturated-zone in slope with different slope angles

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

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

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

4.3 边坡饱和区临界曲线运用 4.3.1 诱发浅层滑坡的坡度预测

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

图 11 不同坡度边坡的饱和区降雨强度-时间临界曲线 Fig.11 Critical curves of rainfall intensity-time for the saturated-zone in slope with different slope angles

不同坡度边坡的饱和区临界曲线反映出降雨过程中边坡产生饱和区存在“危险坡度”.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) 通过多种坡度的饱和区降雨强度-时间临界曲线，可获得区域内边坡发生潜在大量饱和区的危险降雨强度范围.

参考文献
[1]
Matsushi Y, Hattanji T, Matsukura Y. Mechanism of shallow landslides on soil-mantled hillslopes with permeable and impermeable bed-rocks in the Boso Peninsula, Japan[J]. Geomorphology, 2006, 76(1/2): 92-108.
[2]
王述红, 何坚, 王帅, 等. 上覆强透水层的双层土质边坡降雨入渗特征[J]. 东北大学学报(自然科学版), 2020, 41(3): 438-444.
(Wang Shu-hong, He Jian, Wang Shuai, et al. Rainfall infiltration characteristics of double-layered soil slope covered by a highly permeable stratum[J]. Journal of Northeastern University (Natural Science), 2020, 41(3): 438-444.)
[3]
Wang S H, Ni P P. Application of block theory modeling on spatial block topological identification to rock slope stability analysis[J]. International Journal of Computational Methods, 2014, 11(1): 903-914.
[4]
Ali A, Huang J, Lyamin A V, et al. Boundary effects of rainfall-induced landslides[J]. Computers and Geotechnics, 2014, 61: 341-354.
[5]
Mein R G, Larson C L. Modeling infiltration during a steady rain[J]. Water Resources Research, 1973, 9(2): 384-394. DOI:10.1029/WR009i002p00384
[6]
Fredlund D G, Rahardjo H. Soil mechanics for unsaturated soils[M]. New York: John Wiley & Sons, 1933.
[7]
Schaap M G, Van Genuchten M T. A modified mualem-van genuchten formulation for improved description of the hydraulic conductivity near saturation[J]. Vadose Zone Journal, 2006, 5(1): 27-34. DOI:10.2136/vzj2005.0005
[8]
Bodman G B, Coleman E A. Moisture and energy condition during downward entry of water into soil[J]. Soil Science Society of America Journal, 1944, 8: 116-122. DOI:10.2136/sssaj1944.036159950008000C0021x
[9]
Shakoor A, Smithmyer A J. An analysis of storm-induced land-slides in colluvial soils overlying mud rock sequences, Southeastern Ohio, USA[J]. Engineering Geology, 2005, 78(3/4): 257-274.
[10]
李宁, 许建聪, 钦亚洲. 降雨诱发浅层滑坡稳定性的计算模型研究[J]. 岩土力学, 2012, 33(5): 1485-1490.
(Li Ning, Xu Jian-cong, Qin Ya-zhou. Research on calculation model for stability evaluation of rainfall-induced shallow landslides[J]. Rock and Soil Mechanics, 2012, 33(5): 1485-1490.)
[11]
李焯芬, 汪敏. 港渝两地滑坡灾害的对比研究[J]. 岩石力学与工程学报, 2000, 19(4): 493-497.
(Li Chao-fen, Wang Min. Comparison of landslide hazards between Hong Kong and Chongqing[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 19(4): 493-497.)