优先流对土壤水文及滑坡触发的影响

陈思婕，李美骏，罗雅婕，邵 伟*，杨宗佶

(1.南京信息工程大学水文与水资源工程学院，南京 210044; 2.中国科学院水利部成都山地灾害与环境研究所，成都 610041)

滑坡地质灾害一般由地震、极端降雨等因素触发[1-4]。地震在短时间内可以导致边坡荷载的剧烈变化，改变山区岩土特性，造成滑坡或促使大量的边坡转变为滑坡隐患点。在震后期，降雨触发的次生滑坡灾害，具有突发性强、破坏力强等特点，对山区人民财产和生命安全造成极大的威胁，是世界上最主要的地质灾害之一[5-6]。

降雨渗流可造成孔隙内部侵蚀，细颗粒发生迁移[7-8]。这种微观潜蚀机制不仅影响孔隙结构和渗透性，也会减小岩土黏聚力和抗剪强度。经过长期微观潜蚀，裂隙、大孔隙、软弱结构面可能相互贯通，并形成潜在滑动面[9-11]。高强度降雨造成的入渗水流可以绕过致密基岩和土壤基质，在裂隙和大孔隙中快速流动，急剧增大潜在滑动面上的孔隙水压力，并降低抗剪强度[12-13]。因此，在震后极端降雨触发的次生滑坡灾害中，非饱和条件下优先流的快速入渗成为地震滑坡堆积体最显著的影响因素。中外许多学者基于优先流影响下的滑坡预警监测开展了大量的工作，如杨宗佶等[14]采用人工降雨模型试验模拟了典型砾石土滑坡的破坏过程，揭示了降雨条件下滑坡堆积体形成优先流作用下的变形破坏的过程。李鑫等[15]总结了优先流触发滑坡致灾的4个发展阶段，并以典型滑坡作为验证，分析了优先流触发滑坡作用研究机理。潘网生等[16]研究发现黄土孔隙、裂隙中优先流对于滑坡等地质灾害的研究有重要的理论意义和实用价值。王力等[17]采用各向异性渗透材料在降雨作用下对滑坡进行渗流场及位移场的模拟，并与常规模拟及监测数据进行对比，验证了优先流在滑坡事件中的存在。Ronchetti等[18]通过采用染色剂示踪剂法描述斜坡中地下水流运动路径，结果表明，存在两种类型的地下水流，即优先流和基质流，可为滑坡预警提供新的贡献。Chen等[19]通过进行一系列的人工降雨试验，对水文响应进行表征，结果表明大部分降雨渗入滑坡体，受优先流的影响，孔隙水压力和含水量对模拟降雨时间响应迅速。

但目前为止，中国仍缺乏对于优先流影响下的水-力耦合机理的定量研究[20-21]。现行的滑坡预警系统常常忽略优先流对边坡稳定性的影响，造成强降雨条件下滑坡灾害风险的低估[3,6,22-24]。因此，现构建并完善考虑优先流的渗流-边坡稳定模型，将野外实验与数值模拟技术结合构建滑坡预警系统,基于单渗透模型和双重渗透模型，探寻降雨条件下的滑坡堆积体水文响应过程以及失稳过程，并将其应用到地震扰动区的潜在滑坡体的野外实验以及灾害监测对滑坡触发的影响，为提高滑坡预警系统准确度提供新的思路。

1 研究区概况

研究区位于四川省都江堰市白沙河流域银洞子沟(经度103°40′22.57″E；纬度31°9′36.15″N)，距市区约30 km。银洞子沟所在的都江堰市四季分明，属于四川盆地中亚热带湿润季风气候区。降雨充沛但多集中在6—9月的雨季，最低月平均降雨量出现在1月底。同时降雨强度差异极大，降雨空间分布不均匀，东南部平原区大于西北部地区，且降雨量表现为随地势逐渐升高而增加的趋势。该地区降雨为流域泥石流爆发的主要驱动因素。研究区的银洞子沟面积约2.2 km2，主沟整体长约 2.5 km，最大沟壑宽约1.0 km，主沟平均纵坡降 310‰。该流域最高海拔高程2 050 m，最低海拔高程 1 070 m，相对高差1 020 m。山体走向为北东南西向，山脊狭窄，地形为上陡下缓，沟谷深切且狭窄，沟床陡直,多为“V”形谷，其全景图如图1所示。

2008年汶川地震导致银洞子沟两侧山体多处产生崩塌滑坡，形成大量的崩塌堆积体。银洞子沟滑坡为典型的震后滑坡，滑坡堆积体在降雨激发条件下发生再次浅层滑动、拉槽，并转化为泥石流。滑坡堆积体结构松散，主要由砾石土组成。震后降雨作用下，堆积体表面受地表径流冲刷并形成细沟侵蚀，其沟道受到坡面径流的冲刷和渐进浅层滑塌破坏而不断扩大，滑坡体前缘发生局部破坏并不断向后延伸，同时失稳垮塌的滑坡物质，成为后续极端降雨条件下发生次生地质灾害的泥石流物源。

2 野外实验验证

2.1 野外实验系统

人工降雨实验在都江堰市白沙河流域银洞子沟地震滑坡堆积体现场开展，实验系统由人工降雨、摄像记录、地表及地下位移监测组成。坡体的探测装置主要由雨量计、孔隙水压力计、土壤体积含水量探头、数据接收器和倾斜式位移计组成，如图2所示。

图2 实时监测预警系统仪器组成

现场实验布置如图3所示，土壤含水量和孔隙水压力传感器埋藏在滑坡体的关键位置，以便测量含水量和基质吸力的动态变化。翻斗式雨量计可实时收集瞬态降雨量。地表位移和倾斜数据通过GPS(global positioning system)位移计及地表位移倾斜计观测，从而建立滑坡触发与降雨的相关关系。

图3 现场试验方案

2.2 数学模型

2.2.1 瞬态孔隙水压力响应

针对边坡土体瞬时孔隙水压力的变化，单渗透模型采用考虑滑坡坡度的Darcy-Richards方程[24]，即

(1)

双重渗透模型利用两个Darcy-Richards方程分别模拟基质流和优先流[25-26]，即

(2)

(3)

式中：下标m和f分别表示基质流和优先流；wm和wf为基质流和优先流的体积含水量；Γw是两域水分交换量[13]，其计算公式为

(4)

式(4)中：αw为水交换项系数。

采用van Genuchten-Mualem模型描述基质和大孔隙域的土壤水动力特性[27]，即

(5)

(6)

K(Θ)=KsΘ0.5[1-(1-Θ1/m)m]2

(7)

式中：Θ为有效饱和度；θ为体积含水量；θs为饱和含水量；θr为残余含水量；α、n和m(m=1-1/n)为拟合参数；Ss为储水系数；Ks为饱和导水系数。

双重渗透模型基于连续模型概念，土体总含水量θ为基质含水量θm和大孔隙域含水量θf加权平均，即

θ=wfθf+wmθm

(8)

而饱和导水率Ks为

Ks=wfKsf+wmKsm

(9)

式(9)中：Ksf和Ksm分别是大孔隙域和基质域的饱和导水率。

2.2.2 无限边坡稳定性分析

无限边坡稳定性分析方法将浅层滑坡体的滑动面视为与软弱面平行的规则滑动面，稳定性可用安全系数Fs表示为

(10)

式(10)中：Fs为安全系数；c′为有效黏聚力；φ′为有效内摩擦角；G为土体重度；σs为土体基质吸力。

σs=χpw=χγwh

(11)

式(11)中：pw为孔隙水压力；γw为水的容重；χ为基质吸力系数，可用有效饱和度近似计算[28]。考虑到研究区的岩土体渗透系数极高，双重渗透模型中，式(11)中的h可用hm替代进行计算。

2.3 模型设定及参数化

2.3.1 土壤水分特征曲线

首先将Darcy-Richards方程的表面边界设为通量型(flux boundary)，以计算非饱和土表面边界的降雨入渗。模型计算区域深度设为1.0 m，土壤水力学参数分为三层进行设定(0～0.55 m，0.55～0.8 m，0.8～1.0 m)。将26 h的数值模拟周期分为率定和验证两部分。在率定周期内(即0～6 h)对导水率进行率定。根据实测结果，确定了基质与大孔隙域体积比相关参数。大孔隙域体积比设定为25%，表示土壤中大孔隙、裂缝等优先流路径的体积比。测定的饱和土壤含水量范围为0.4～0.5，其下限与基质域孔隙度相关，而上限则与大孔隙所占体积有关。因此，在单渗透模型和双重渗透模型中，残余土壤含水量θr和饱和土壤含水量θs分别为0.1和0.425。

土壤水力学特征曲线中，参数α与进气吸力有关，而n为拟合参数。其中，大孔隙域的n，各土层深度设定为4.0；而基质域的n利用MATLAB软件中的“Isqcurvefit”函数，基于非线性最小二乘算法拟合参数，以描述土壤水分特征曲线。单渗透模型拟合后的土壤水分特征曲线如图4(a)、图4(c)、图4(e)所示，双重渗透模型确定的复合土壤水分特征曲线如图4(b)、图4(d)、图4(f)所示。

图4 单渗模型和双重渗透模型中土壤水分特征曲线

双重渗透模型的复合土壤水分特征曲线可以用来描述土壤整体的持水特性，与单渗透模型和实测结果基本吻合。双重渗透模型的土壤水分特征曲线采用两组参数描述基质流和优先流的水力特性。依据脱湿过程的土壤含水量和孔隙水压力监测数据，将其作为单渗透模型和双重渗透水分特征曲线参数，其参数如表1所示。

表1 单渗透模型、双重渗透模型的基质和大孔隙域对应的水力学参数

2.3.2 土力学相关参数

土力学参数(即黏聚力、摩擦角、比重)和边坡几何关系(厚度和长度)如表2所示。基于非饱和土直剪试验测试了土体的有效黏聚力和有效内摩擦角。土体力学参数将用于非饱和条件下的边坡稳定性分析，土力学参数由文献[21]进行室内实验测试而来。

表2 滑坡体土体参数

3 模拟结果

3.1 数值模拟和实测结果的比较

三种不同深度(0.50、0.75、1.0 m)土壤含水量模拟结果如图5所示。开始时，各深度土壤含水量与数值分析中指定的残余土壤含水量θr非常接近(如表1所示)。率定期间降雨强度较为均匀，平均102 mm/h，在5 h过程中，土壤含水量由于重力排水而逐渐降低，达到稳定值，在0.50、0.75、1.0 m的土层深度分别降低0.12、0.18、0.25。在边坡失稳破坏期间，4次降雨强度均引起土壤含水量的显著变化，降雨强度变化条件下，相对复杂的降雨输入导致入渗和排水过程显著变化。

图5 不同土层深度在率定期、排水期、边坡破坏过程中实测和模拟土壤含水量

由图5可知，在0.50 m深度处，单渗透模型和双重渗透模型基质域土壤含水量模拟值均与实测值吻合较好。各深度实测土壤含水量对降雨的响应非常迅速，数值模拟土壤含水量在0.75 m和1.0 m深度存在明显的滞后。在1.0 m时，单渗透模型的结果为显著延迟0.33～0.5 h，双重渗透模型基质域土壤水文响应仅滞后0.17～0.33 h。在模拟和实测的土壤水分对比中，模拟的土壤水动力过程表明，入渗主要发生在基质域，而大孔隙域的水文响应在达到较高的饱和条件后才会显著，故在1.0 m处双重渗透模型能更好地模拟真实土壤含水量。双重渗透模型能体现基质域和大孔隙域的水力特性，比单渗透模型能更好地揭示滑坡体中土壤含水量的快速响应。

非饱和土孔隙水压力的负值可以定义为土体吸力，其在三个不同深度的值如图6所示。0.50、0.75、1.0 m深度处土体初始吸力P分别为14.7、4.5、1.9 kPa。降雨期间各深度土体吸力显著降低至1.0 kPa。间歇排水期间，土壤含水量降低，土体吸力逐渐增大。在边坡失稳破坏过程中，由于强降雨作用，土体吸力下降至最低值。

图6 不同土层深度率定期、排水期、边坡破坏期中实测和模拟土壤孔隙水吸力

以上对比表明，单渗透模型和双重渗透模型都可获得较为合理的土壤含水量和基质吸力。但双重渗透模型，可更准确地反映快速水文响应。

3.2 边坡稳定分析

现场实测边坡破坏滑移面深度为0.5 m，因此将边坡稳定性分析的深度定为该深度以模拟滑移面的稳定系数。深度为0.5 m处的安全系数计算如图7所示。初始条件为稳定(Fs>1.0)，当安全系数Fs<1.0时，则发生边坡失稳破坏。单渗透模型模拟的初始安全系数为1.115，然而双重渗透模型模拟的安全系数更高为1.120。试验期间模拟的安全系数随降雨作用逐渐降低。单渗透模型模拟的安全系数下降至1.005，双重渗透模型模拟的安全系数下降至1.000，双重渗透模型计算的安全系数的波动大于单渗透模型的安全系数。在排水过程中，由于孔隙水压力的增大，两种渗透模型模拟的安全系数均增大，分别达到1.025(单渗透模型)和1.045(双重渗透模型)。

图7 单渗透和双重渗透模型模拟的边坡稳定性

在高强降雨作用下，边坡共发生了4次滑坡破坏事件。在此期间，土体始终处于非饱和状态，因此滑坡全部发生在非饱和条件下。第一次滑坡(t=22 h)的触发可通过两种模型来进行预警:单渗透模型的安全系数下降到0.995，双重渗透模型模拟的安全系数下降到0.99。随后，第二次和第三次(t=23.5 h和t=24 h)失稳破坏发生于低雨强(20 mm/h)降雨过程，单渗透模型和双重渗透模型结果均高估了边坡稳定性。第四次(t=24.7 h)降雨作用持续时间较长，大约为40 min，单渗透模型模拟的安全系数仍然接近1.0，而双重渗透模型模拟的安全系数为0.99。因此，双重渗透模型可以更好地揭示边坡失稳破坏过程，而单渗透模型则因无法量化优先流对滑坡触发的影响而高估边坡稳定性。

4 结论

研究以汶川地震扰动区银洞子沟滑坡易发区的野外人工降雨滑坡实验为例，采用单渗透模型和双重渗透模型模拟滑坡体的土壤水动力过程，并分析降雨入渗和滑坡失稳机理，得到如下结论。

(1)单渗透模型无法准确表征优先流影响下的水力特征，只能模拟基质流，导致模拟土壤含水量与实测之间存在明显的滞后效应。相反，双重渗透模型既可以模拟优先流，也可以模拟基质流，可以更准确地模拟土壤含水量和孔隙水压力的变化。

(2)将双重渗透模型与边坡稳定性计算相耦合，可以量化土体非匀质性对边坡稳定性的影响，为滑坡泥石流灾害预警提供更可靠的结果。