不同含水率条件下非饱和边坡的稳定性

鲜 丹, 杨仲康

(1.西北大学 城市与环境学院, 陕西 西安 710100; 2.兰州大学 资源环境学院, 甘肃 兰州 730000)

在黄土高原地区，降雨是滑坡灾害最为重要和常见的诱发因素之一[1]，降雨入渗导致边坡失稳主要表现在：一方面是雨水入渗产生渗流作用，土体容重增高，边坡下滑力增加；另外一方面是雨水入渗，软化土体，使岩土体抗剪强度降低。众多研究表明[2]，入渗致滑的因素是多方面的，Rahardjo等[3]认为在饱和—非饱和渗流稳定分析中，对初始条件的正确考量是解决问题的关键之一。

初始条件主要指初始地下水位及初始孔隙水压力或含水量的分布。饱和状态时孔隙水压力可以较方便准确的测定，而实际情况中获取更为普遍的非饱和状态负孔隙水压力则较为困难。为了正确考虑初始含水率条件，朱伟等[4]提出了一种以平均渗透强度来确定初始含水率的计算方法，袁俊平等[5]通过现场实测和数值模拟方法总结了初始孔隙水压力确定的一般方法并分析各自的适用范围。而在考虑初始条件进行非饱和渗流稳定分析的研究中，付建新等[6]基于有限差分方法，直接设置不同含水率作为初始条件进行二维入渗分析，认为土体初始含水率条件对入渗和稳定性的影响是随降雨强度的增加而降低的。Rahimi等[7]则在降雨边界上施加一个长期的小强度降雨来确定初始条件，詹良通等[8]在数值解析中取平均降雨量。此外，更为普遍的方法是通过考虑前期降雨实现对初始条件的考量，唐栋等[9]通过设置初始最大孔隙水压力为-75，-50，-25 kPa作为初始条件考虑对稳定性影响，并认为当考虑前期降雨天数30 d以上时，初始条件对安全系数的影响十分微弱。罗红明[10]通过考虑不同前期降雨过程作用下形成的含水量分布作为后期降雨的初始条件，将降雨过程分为前期降雨和主降雨，分析前期降雨持时、强度和与主降雨的时间间隔与稳定性关系。其中直接设置坡体含水率简单有效，但数值计算中设置土体含水率处处相等并设定入渗深度，忽略了入渗的物理过程，解析结果仍需进一步验证。而通过施加前期降雨获取稳态渗流结果作为初始条件，事实上，在不同地区，不同地层条件下前期降雨的影响十分复杂，Segoni等[11]研究也认为利用前期降雨指数来考虑初始水率条件，获取结果与实际观测值存在一定出入。

本文以甘肃省天水市南部廖集村滑坡为例，综合现有评价方法，基于现场实测土体含水率将其分为天然和湿润状态，通过建立渗流与应力的流—固耦合物理力学模型，力图改进现有初始条件对非饱和渗流与稳定性影响的评估方法，着重探讨不同初始含水率条件下降雨入渗对边坡稳定性的作用机制，并初步尝试提出滑坡预警的建议，进而应用于天水南部地区降雨诱发型滑坡的预警和防灾减灾工作。

1 研究区概况

2013年6月19日至7月25日，甘肃省天水市出现多次罕见暴雨天气，由李子气象站记录共有4次强降雨天气，6.20，7.8，7.22，7.25日降雨量分别为198.2，128.6，22.7，29.8 mm，共降雨379.3 mm，此次连续性暴雨天气共诱发滑坡415处，造成24人死亡，1人失踪，简称“7.25群发性滑坡灾害”[12]。经初步勘察发现，前两次突发性暴雨天气诱发了大量浅表层土质滑坡。后两次降雨量虽然较小，但在前期丰富降雨作用下土体含水率大幅度提高，呈湿润状态，滑坡事件依然发生，这些滑坡往往为黄土滑坡，滑带埋藏较深，虽然数量小但滑坡规模大，潜在危害更大。

2 流—固耦合数值模型

2.1 计算原理

降雨入渗导致斜坡内渗流场的变化，土体容重增加，下滑力增大，抗剪强度降低，其中所涉及的物理力学机制利用数值模型模拟由降雨引发的孔隙水压力重分布和滑坡稳定性评估已被广泛应用[13]。具体的数值分析过程为：基于饱和—非饱和渗流理论，获取入渗条件下斜坡土体的瞬态含水率分布，并根据非饱和抗剪强度理论，利用技术比较成熟的极限平衡法进行边坡稳定性分析[14]。

将达西定律和质量守恒定律相结合就可推导出二维情况下以水头H为因变量的非饱和渗流控制方程，即Richards方程，微分方程为：

(1)

式中：H——全水头，H=h+z(压力水头+位置水头)；Kx，Ky——土体沿X，Y方向的渗透系数；θ——体积含水率(%)；t——时间(h)。

对控制方程采用Galerkin有限元格式，将孔隙水压力作为节点自由度进行空间离散，获取区域渗流分析结果。为了考虑负孔隙水压力对抗剪强度的影响，Fredlund & Rahardjo引入材料变量φb来定量分析非饱和土抗剪强度随基质吸力增加而增加的非线性关系，公式即为：

τf=c′+(σ-uw)tanφ′+(ua-uw)tanφb

(2)

式中：c′——有效黏聚力(N)；φ′——有效内摩擦角(°)；φb——对应基质吸力的内摩擦角(°)，反映抗剪强度增长率与基质吸力的关系；ua，uw——孔隙气压力和水压力(kPa)；σ——条块地面法向正应力(N)。

2.2 建构模型

2.2.1 设置初始条件 非饱和土是由固体颗粒、水及空气构成的三相体，土体中水的流动是由土水势所决定的，基质势的空间分布显著影响土水势，基质势即为基质吸力(即非饱和土的负孔隙水压力)，基质吸力与土体体积含水率密切相关，Ning Lu等[15]认为非饱和土中孔隙压力水头的大小很大程度上取决于土体饱和度(含水量)与土体类型等。故尝试引入表达体积含水率与孔隙水压力关系的Boltzman模型[16]，由此孔隙水压力的设置就有了含水率的意义。Boltzman模型为：

θ(uw)=θseαuw

(3)

式中：θ(uw)——对应孔隙水压力时的体积含水率(%)；θs——饱和体积含水率(%)；α——饱和系数，本文取值0.01。

吴玮江等[17]研究表明，天水地区考虑灾害发生前7 d的有效降雨量，对分析降雨触发滑坡过程最有意义。故在7 d内未降雨条件下测得土体平均含水率为17.63%，在7 d内发生降雨且降雨总量达到80 mm测得土体平均含水率为19.06%。依照土力学对饱和度分类标准将其分别定义为天然和湿润状态，基于实测结果，对孔隙水压力分布进行概化处理，借鉴包承纲等[18]在研究膨胀土边坡时的做法，假定地下水位以上孔隙水压力随深度呈线性分布，在岩土模拟Geostudio软件SEEP/W模块中中设置不同含水率条件下孔隙水压力的分布，首先，模拟多日未降雨情形，设置干燥状态，表层孔隙水压力控制在102 kPa附近；然后，模拟存在前期降雨情形，设置为湿润状态，孔隙水压力设置在-63 kPa左右，具体如表1所示。

表1 不同含水率下材料参数

此外，土体含水率不同，土体抗剪强度也会产生深刻影响，随着含水量的增高，抗剪强度会降低，尤其是在黄土地区[19]。所以为了更加准确进行稳定性评估，采集野外滑带样品，室内利用VJT环剪仪，重塑环状土样测定各对应含水率下的剪切强度，因为廖集村滑坡已显著位移，现阶段属于峰值后的变形与破坏[20]，故取残余剪切强度作为材料初始强度参数(表1)，其中φb取值为2/3φ。

2.2.2 模型参数和降雨工况 土—水特征曲线(SWCC)是用于描述土吸力与含水量之间本构关系的函数，是模型进行非饱和渗流分析所需的重要参数。通过颗粒级配曲线和密度等参数对比，选取李萍等[21]在类似黄土地区实测基质吸力和含水率数据(张力计法)，并用Fredlund & Xing模型进行拟合，将Boltzman模型计算的孔压与其实测结果验证，具有较高吻合度。Fredlund & Xing拟合模型公式为：

(4)

式中：ψ——基质吸力(N)；a，b，c——拟合曲线形状参数，拟合参数取值分别为11.4，2.0，0.9。

非饱和渗透系数函数用于描述渗透系数对土体系统内孔隙水含量的依赖程度，但其变化范围很大且不易测量，通常做法从土水特征曲线获取。利用SEEP/W模块内置程序，采用Fredlund & Xing方法，基于实测土水特征曲线，进行估算。其中饱和渗透系数1.97×10-3cm/s[22]，饱和含水率47.93%，残余饱和度为10%。

研究区多年年降雨总量在600～800 mm，降雨主要集中夏季7—8月，多以暴雨形式出现。故优先考虑暴雨工况，降雨总量取50年一遇为480 mm，降雨强度为10 mm/h，持时48 h。此外，另设置连阴雨天气工况，着重考虑小强度降雨入渗，考察土体吸湿过程和降雨历时的影响，降雨强度为3 mm/h，持时160 h。

进行应力分析时，主要利用Geostudio软件SIGMA模块，通过应力重分布和应力—孔隙水压力耦合来计算，边界条件处理为两侧水平约束和底部完全固定边界[23]。为提高模型计算精度，二维模型网格划分为四边形和三角形单元，为减少高度非线性系统中经常出现的数值震荡和发散问题，非饱和区域单元进行细化处理，共7 445个节点7 430个单元。

3 结果分析与讨论

3.1 非饱和渗流过程

只有计算过程收敛所得到的瞬态渗流分析结果才能合理有效，差值迭代图是十分直观有效的收敛判别依据[24]。如图1差值虽然震荡，但数值递减，最终计算收敛。

图1 瞬态渗流过程差值迭代图

基于SEEP/W模块数值解析结果，得到瞬态孔隙水压力场。通过追踪孔隙水压力变化获取降雨入渗信息。首先，获取横向入渗信息(图2)，随着入渗的开始，孔隙水压力迅速升高，基质吸力锐减，直至稳定在较低值，反映了由浸润渗透转化为稳定渗透的非饱和—饱和的入渗过程。

其中，天然状态时，土体间含水率梯度值大，基质势能高，土体能更快速持水，在12 h内孔隙水压力就由-100 kPa激增至-6 kPa，基质吸力基本消失，达到稳定渗流。通过提取剖面AB处孔隙水压力随深度变化结果，获取垂向入渗信息(图3—4)，天然状态时，垂向渗透过程明显滞后于湿润状态。暴雨天气时，降雨6 h，天然状态入渗深度不到1 m，湿润状态已达到2 m，降雨24 h后，天然状态下入渗深度基本停留在3.2 m左右，而湿润状态继续入渗，深度达到4 m(图3)。连阴雨天气时，由于降雨历时的延长，入渗深度差异更加显著，湿润状态时-5 m以下土体呈过饱和状态，下渗到了更深部位可能已经引起了下部浸润面的上升(图4)。

图2 坡面上一点孔隙水压力与降雨持时的关系

图3 暴雨天气下剖面AB各深度孔隙水压力与降雨历时的关系

图4 连阴雨天气剖面AB各深度孔隙水压力与降雨历时的关系

此外，在不同含水率条件下各深度孔隙水压力变化也呈现不同的特征，在天然状态下，边坡表层1.5 m深度内孔隙水压力被显著提升到-6 kPa，接近正值，然而1.5 m深度以下孔隙水压力改变有限，土体仍能保持较强吸力。湿润状态时，相比较坡体表层，深度在2 m以下孔隙水压力被明显提升，曲线接近垂直，说明该部位渗流作用十分显著。

3.2 滑坡稳定性结果

SEEP/W模块非饱和渗流模拟得到瞬态孔隙水压力场分布，将其导入到SLOPE/W模块进行耦合稳定性计算。结果如表2所示，首先，天然状态下初始稳定系数为1.34，而湿润状态下仅为1.16，初始稳定系数存在较大差异。说明考虑不同含水率条件下力学参数的变化而进行不同含水率状态下的剪切试验是十分有必要的。在不同含水率条件下对各雨型的响应是不同的，天然状态下对暴雨天气响应更加敏感，最小安全系数下降幅度为5.22%，湿润状态下对连阴雨天气响应更加积极，下降幅度高达6.9%。细观安全系数随降雨历时演化曲线(图5)，发现不同初始含水率条件下，安全系数随降雨历时具有不同的阶段性变化特征。以暴雨天气为例，天然状态时在降雨前12 h内Fs迅速下降，12 h以后Fs下降速率明显放缓，衰减特征为先快后慢，湿润状态时在降雨前18 h内Fs缓慢变化，然而在18 h以后开始迅速下降，并且呈现降雨之后继续下降的趋势。

表2 稳定性计算结果

图5 不同雨型条件下各含水率状态安全系数与降雨历时的关系

3.3 降雨入渗致滑作用机制分析

天然状态时，初始含水率仅为17.63%，基质吸力强，持水迅速，仅降雨12 h土体孔隙水压力激增95%，饱和度升高，基质吸力锐减，所以此过程中稳定系数迅速下降。经历短暂的浸润渗透进入到稳定渗透，此阶段之后表层1—2 m内土体孔隙水压力接近正值，不再变化，而下渗基本停滞，其余水分沿着坡表面排出，对深部土体孔隙水压力影响微弱，使其仍保持较强的吸力，所以稳定系数开始下降缓慢。从物理过程而言，干燥坡体遇突发性暴雨，迅速持水，在坡表形成暂态饱和区，入渗使下滑力增大，坡体受力而产生超静孔隙水压力[25]，由于持续入渗而超静孔隙水压力不能被快速消散，土体孔隙被水充填超过其临界孔隙率[26]，最终土体产生类似于“液化”的塑性滑动。在研究区普遍发育，为浅表层土质滑坡，规模小，发生时间短多在降雨中或降雨后迅速发生，滑坡过程多呈“溜滑”形式[27]。

湿润状态时，初始含水率高，坡体内含水率梯度小，土体吸湿持水慢，所以在降雨前12 h内稳定系数变化缓慢，但随着降雨历时的延长，雨水能够持续向深部渗透并使其基质吸力锐减，此时稳定系数开始迅速下降。具体演化过程通过SIGMA模块对滑动面进行应力分析(图6)，在降雨12 h以前，入渗还停留在浅表层1—2 m处，滑动面附近孔隙水压力尚未受到影响，抗剪强度远大于剪应力，稳定系数变化不显著。降雨12 h以后，入渗深度达3—4 m，接近滑动面附近，此时孔隙水压力被显著抬升，土体被软化，如图中24 h时曲线所示，抗剪强度下降，剪应力上升，此时稳定系数开始迅速下降，并且这种应力消涨变化在0—40 m(坡脚)变化最为剧烈。

48 h时，剪应力大范围超过抗剪强度，稳定系数持续下降濒临失稳。暴雨天气时，观测时间短且稳定性变化随降雨具有滞后性，湿润状态时稳定系数下降幅度比天然状态时小，当连阴雨天气时，降雨历时长，降雨入渗对深部滑带软化作用则更加显著，稳定系数下降幅度为6.9%，Fs为1.08，由此可判定此时已达到滑坡失稳临界值。这类滑坡往往为深层滑坡，滑带深度在6～12 m左右，滑坡显著滞后于降雨，存在多级滑动，率先从坡脚开始，滑坡过程为“后退式渐进破坏”。

图6 滑面上各节点应力状态与降雨历时的关系

4 结 论

(1) 在非饱和渗流分析中，基质吸力起着关键作用，通过引入Boltzman模型使不同孔隙水压力的设置具有了含水率的意义，从而实现在土体天然和湿润2种初始条件下分析问题，计算结果也证明在对应含水率下进行剪切试验作为模型参数是很有必要的。

(2) 在天然状态下，土体初始含水率低，基质吸力强，持水快，但入渗仅局限于浅表层1—2 m，所以入渗一旦开始，稳定系数迅速下降，随着入渗停止下移，稳定系数下降放缓。入渗造成边坡表层高饱和状态，使其滑坡特征为类似于液化的“溜滑”运动而形成浅表层滑坡。

(3) 在湿润条件下，含水率梯度小，入渗缓慢，但能入渗到边坡更深部位，并使其孔隙水压力发生显著变化。所以入渗之初，稳定系数无明显变化，只有入渗下移到3～4 m时，靠近滑带部位，稳定系数才开始迅速下降，往往滞后于降雨而形成深部滑坡。

天水市南部地区滑坡预警工作我们建议应该注意以下几点：在多日未降雨土体较干燥时，尤其注重对突发性暴雨天气的防范，这种条件下，滑坡发生迅速，滑动面浅，数量大，单体监测困难，雨量预警则是较高效的方法，可以借鉴彭建兵等[26]在该地区提出的降雨阈值前7 d内降雨总量超过239 mm，同时最大最强达到20 mm/h就可判断会有滑坡发生；当处于持续的连阴雨天气，土体呈湿润状态时，入渗能够到较深部位，此情形下应当防范滑坡的滞后效应，往往发生在降雨数十小时以后，重点监测厚层黄土分布区深度在4 m以下孔隙水压力的变化，同时关注斜坡前缘位移状态，共同作为预防此类深部滑坡发生的措施。