初始地下水对浅层边坡降雨入渗及稳定性影响

赵 宇,马朋博,梁 腾,刘 准,吕 庆

(1.浙江大学 建筑工程学院,杭州 310058;2.软弱土与环境土工教育部重点实验室(浙江大学),杭州 310058;3.浙江大学 超重力研究中心,杭州 310058)

近年来,中国滑坡灾害数量一直呈上升趋势,其中50%以上由极端降雨直接诱发[1],浙江省作为中国降雨型滑坡最频发的地区之一,地下水资源丰富且分布复杂[2];由降雨引发的滑坡深度一般较浅,滑面通常为与坡面平行的湿润锋面[3]。一般认为,降雨入渗会使边坡湿润层土体基质吸力消失、重度增加,进而导致土体抗滑力减小、下滑力增大;当土体抗滑力难以维持边坡稳定时即会发生失稳[4-5]。

目前,对于非饱和土入渗过程的分析主要有以Kostiakov模型和Horton模型等为代表的经验方程模型[6]、以Iverson等研究为代表的Richard渗流理论模型[7-8]和Green等[9]基于毛细理论提出的薄层积水入渗模型(Green-Ampt模型)。其中,Green-Ampt(GA)模型假设土体初始含水率均匀分布,且入渗时存在明确的湿润锋面,锋面上下两侧土体含水率分别为饱和含水率和初始含水率,忽略土体从饱和含水率到初始含水率的过渡区域;由于GA模型物理意义明确且形式简单,近年来,大量学者针对其适用范围及主要参数确定等方面开展研究。Chen等[10]在GA模型的基础上推导了考虑坡度的边坡降雨入渗模型,该模型假设土体初始含水率均匀分布且湿润层土体为饱和状态;Gavin等[11]的改进GA模型认为降雨后边坡湿润层土体含水率并未饱和并假定土体基质吸力沿深度线性变化;李宁等[12]认为在产生径流前湿润层土体未饱和,提出采用非饱和土Van Genuchten(VG)模型确定湿润层土体含水率;马世国等[13]研究了强降雨和地下水共同作用下的入渗模型,但未考虑坡角和弱降雨工况对于湿润锋推移过程的影响;张洁等[14]推导了斜坡初始吸力为三角形和梯形分布的改进GA模型,其基于Runge-Kutta原理的数值求解过程复杂;简文星等[15]研究了坡体含水率沿深度近似为直线——指数曲线分布的降雨入渗模型。可见,现有研究一般采用纯数学方程描述模型边坡土体初始含水率/基质吸力的分布,难以准确反映地下水富集地区坡体初始含水率真实分布情况;同时,稳定性分析时一般认为湿润锋处为潜在滑动面,忽略了由于初始地下水影响提前在基岩面处失稳的可能性,难以准确评价浙江省山地丘陵地区基岩型浅层滑坡灾害致灾机制。

综上,首先基于浅层边坡地下水分布特征及土体非饱和特性,提出考虑初始地下水影响的坡体含水率指数分布模型;基于GA入渗模型及水量平衡原理推导了考虑初始地下水影响的浅层边坡在弱/强两种降雨工况下入渗方程解析解;通过与传统模型和数值模拟结果对比,验证了本文坡体含水率分布模型和降雨入渗方程解析解的有效性;最后,结合浅层无限边坡模型,对比分析了边坡在湿润锋和基岩面两个潜在危险滑动面处的稳定性演化规律。

1 考虑初始地下水影响的浅层边坡降雨入渗过程分析

1.1 考虑初始地下水影响的坡体含水率分布

浅层边坡下覆弱透水基岩层,基岩上部土体中大多存在稳定的初始地下水,非饱和土毛细作用使土体含水率分布呈现从坡面到基岩面逐渐增大的趋势(图1)。Baum等[16]研究表明,当存在地下水时,坡面到地下水位之间的土体含水率分布可以参照Gardner等[17]提出的非饱和土土-水特征曲线(SWCC)模型进行改进,其模型为

图1 考虑初始地下水影响的坡体含水率分布

θ=θr+(θs-θr)e-αψ*

(1)

式中:θ为体积含水率;θr为残余体积含水率;θs为饱和体积含水率;α为对SWCC曲线拟合得到的参数;ψ为基质吸力水头,且ψ*=ψ-ψ0,令ψ0=1/α以使该SWCC模型具有良好拟合效果[18],其中,1/α表示地下水位线以上的毛细饱和带高度。

假定边坡土体基质吸力水头ψ与垂直于坡面方向的深度z存在线性关系[19],即

ψ=Hw-z

(2)

式中Hw为地下水深度。

因此,坡面到地下水毛细饱和带之间的土体体积含水率分布可表示为指数分布形式,即

θz,i=θA+(θs-θA)e(1-αHw)eαz

(3)

式中θA为根据现场实测数据或数值模拟确定的拟合参数。

1.2 考虑初始地下水影响的降雨入渗解析解

边坡降雨入渗是一个二维问题,但在浅层边坡的降雨入渗及稳定性分析时,为求解方便,常将该过程简化为一维问题[8,11,14]。假定浅层边坡坡度为β,基岩面以上土层存在稳定的初始地下水,且地下水位与基岩面平行,如图2所示。下面根据降雨强度q与土体饱和渗透系数ks的大小关系分为弱降雨和强降雨工况分别进行分析。

图2 考虑初始地下水影响的浅层边坡降雨入渗模型

1.2.1 弱降雨工况(q≤ks)

当降雨强度q小于或等于土体饱和渗透系数ks时,降雨全部入渗到坡体内部,此时降雨入渗率i(t)为

i(t)=qcosβ

(4)

由水量平衡原理[20],降雨时间t内坡面法向累计入渗量I等于坡体在该方向体积含水率的增量,即

(5)

式中:zw为降雨时间t对应的湿润锋推移深度,θw为湿润层土体体积含水率。根据李宁等[12]提出的方法:由非饱和土渗透系数方程与降雨入渗率的对应关系求得湿润层土体基质吸力水头ψ,将ψ代入土体含水率-基质吸力关系曲线模型可得降雨强度q对应的湿润层土体体积含水率θw,该方法中描述非饱和土渗透系数和计算θw的方程[21]如下:

(6)

(7)

式中:αvg、n和m为与SWCC曲线有关的拟合参数,m=1-1/n。

将式(3)代入式(5),可得

(8)

对式(8)积分即得当存在初始地下水时弱降雨工况下降雨时间t与湿润锋推移深度zw的关系,为

(9)

1.2.2 强降雨工况(q>ks)

当降雨强度q大于土体饱和渗透系数ks时,降雨过程应根据降雨入渗率i(t)的变化分为两个阶段:降雨强度控制阶段和土体入渗能力控制阶段,如图3所示。

图3 强降雨工况下入渗率与降雨时间关系曲线

降雨初期(0≤t≤tp),降雨强度q小于土体入渗能力,即降雨强度控制阶段,此阶段降雨入渗率i(t)同式(4),降雨时间t与湿润锋推移深度zw的关系同式(9)。

随着降雨持续进行,土体入渗能力不断减小至小于降雨强度后(t≥tp),坡面产生径流,即土体入渗能力控制阶段,此阶段降雨入渗率i(t)为

(10)

式中sf为湿润锋处基质吸力水头,可按下式计算[22]:

(11)

式中:ψi为土体体积含水率θi对应的基质吸力水头,k(ψ)为渗透系数方程即式(6)。

坡面开始产生径流时,降雨强度q沿垂直于坡面z方向的分量刚好等于该方向土体入渗能力大小,令开始产生径流时的湿润锋推移深度为zp,则有

(12)

求解式(12),可得

(13)

令开始产生径流时的累计入渗量为Ip,则有

(14)

将式(3)代入式(14)并积分,可得

(15)

则坡面产生径流时间tp为

(16)

将式(15)代入式(16),可得

(17)

当t>tp时,由水量平衡原理可得

(18)

对于浅层边坡,参照文献[10]的假定,可认为

(19)

将式(19)代入式(18),并对式(18)积分后整理可得当存在初始地下水时,强降雨工况下土体入渗能力控制阶段降雨时间t与湿润锋推移深度zw的关系为

(20)

综上,当存在初始地下水时,强降雨工况下降雨时间t与湿润锋推移深度zw的关系为

(21)

2 考虑初始地下水影响的浅层边坡稳定性分析

当浅层边坡下覆基岩处土层存在初始地下水,且降雨为诱发滑坡主要因素时,根据文献[13]的结论,应同时考虑基岩面处与湿润锋处两个潜在危险滑动面(图4)。下面将上文改进降雨入渗模型与浅层无限边坡模型结合,分别分析两个潜在危险滑动面的安全系数随降雨时间的变化关系。

图4 浅层边坡稳定性分析

2.1 湿润锋处稳定性分析

湿润锋位置zw随降雨时间t的增加而不断向下推移,考虑湿润锋处作为一个潜在危险滑动面,如图5(a)所示,采用安全系数Fs评价边坡稳定性,其定义为土体抗滑力与下滑力之比,其中,抗滑力采用Fredlund 等[23]提出的土体抗剪强度理论求解,下滑力为湿润层土体重度沿坡面方向的分量,则湿润锋处安全系数可表示为

图5 单宽土条分析

(22)

式中:τf为土体抗滑力;τm为土体下滑力;c′和φ′分别为土的有效黏聚力和有效内摩擦角;σn为单宽土条底部正应力,σn=γtzwcosβ;γt为湿润层土体重度;ua为孔隙气压力,假定为0;uw为孔隙水压力;(ua-uw)为土体基质吸力,当土体饱和时,该项为0;tanφb为抗剪强度随基质吸力(ua-uw)增加而增加的速率,根据Vanapalli等[24]提出的φb与土体体积含水率间的经验模型可表示为

(23)

2.2 基岩面处稳定性分析

当浅层边坡基岩面以上土层存在稳定的初始地下水时,由非饱和土毛细作用导致的毛细饱和带至基岩面之间的土体呈近饱和状态,基质吸力基本消失;对基岩面以上单宽土条沿深度分为3部分进行分析,如图5(b)所示。

第1部分:坡面到湿润锋面之间(0≤z≤zw)

(24)

式中γsat为土体饱和重度。

第2部分:湿润锋面到毛细饱和带之间(zw≤z≤Hw-1/α),在其任意深度z处取dz段分析,则dz段自重应力可表示为

(25)

式中γz,i为深度z处土体重度。

由土体基本物理特性可推得

(26)

γz,i=ρz,ig

(27)

式中:ρz,i为深度z处土体密度;ρw为水的密度;ρs为土粒密度;C/V为土体孔隙体积与总体积之比,假设入渗过程中土体含水率的变化不引起土体体积的变化及土颗粒骨架的变形[23]。

将式(26)、(27)代入式(25),得

(28)

将式(3)代入式(28),并积分可得

(29)

第3部分:毛细饱和带到基岩面之间(Hw-1/α≤z≤H)

(30)

式中γ′为土体有效重度。

基岩面以上单宽土条底部正应力与切应力可分别表示为

σn=∑Wicos2β

(31)

τm=∑Wicosβsinβ

(32)

式中∑Wi=W1+W2+W3。

因此,基岩面处安全系数可表示为

(33)

3 模型验证与分析

首先将改进模型计算结果分别与传统模型和数值模拟结果对比,验证改进模型的有效性;然后利用改进模型分别计算不同降雨强度下湿润锋推移深度与降雨时间的关系;最后对存在初始地下水影响的浅层边坡在降雨条件下的稳定性进行了讨论。

浙江省山地丘陵地区降雨型滑坡基岩上覆土层深度大多在4～8 m,且坡度集中在25°～35°[25]。因此,假定一浅层边坡(图6),坡长为65 m,坡度为30°,垂直坡面土层厚度H为6.5 m,地下水深度Hw为5.8 m。选取浙江省山地丘陵地区边坡常见残积土为上覆土体[26],算例中土体参数见表1,土-水特征曲线和渗透系数曲线分别如图7、8所示。通过计算,降雨导致的湿润层土体体积含水率θw为0.37,此时湿润层土体重度γt为19.5 kN/m3。

表1 土体参数

图6 均质浅层边坡分析示意

图7 土-水特征曲线

图8 渗透系数曲线

3.1 改进模型有效性验证

3.1.1 改进模型与传统模型计算结果对比

当坡体含水率分布与Chen等[10]假设(即初始含水率沿深度均匀分布)相同,改进模型可简化为

(34)

利用式(34)与Chen-Young模型[10]分别计算当降雨强度为13 mm/h时(假定坡体初始体积含水率为0.18)湿润锋推移深度zw与降雨时间t的关系,如图9所示。可以看出,两者的湿润锋推移曲线比较接近,降雨持续60 h,两条曲线几乎重合,当湿润锋推移深度达到396 cm时,其相对最大误差小于3%,且改进模型计算结果偏于安全。可见,当假定初始含水率均匀分布时,改进模型与Chen-Young模型是等效的,改进模型是对传统模型在应用范围上的扩展。

图9 含水率均匀分布时改进模型与传统模型计算结果对比

3.1.2 改进模型与数值模拟计算结果对比

采用SEEP/W数值模拟软件对浅层边坡存在初始地下水时的含水率分布及降雨入渗过程进行分析,以进一步验证改进模型的有效性。分析模型如图10所示,其中,虚线为边坡初始地下水位线,其深度与图6一致,并预先进行稳态渗流分析。

图10 数值模拟分析模型

首先利用式(3)对SEEP/W软件中稳态渗流分析获得的坡体体积含水率沿深度分布数据点进行拟合(图11),其中,拟合得到的参数θA为0.10,拟合优度达到0.96,可见采用式(3)描述存在初始地下水的浅层边坡土体体积含水率沿深度分布是可行的。

图11 坡体初始含水率沿深度分布曲线

然后将雨强定义为流量边界进行降雨入渗过程分析,通过定义孔隙水压力值获得降雨时间t对应的湿润锋推移深度zw[8],并与本文改进模型结果进行比较。如图12所示,两种方法计算得到的降雨强度为11 mm/h时湿润锋推移深度zw与降雨时间t对应关系比较接近,但数值模拟得到的湿润锋推移深度要比改进模型结果稍大,在降雨持续95 h内最大偏差为12 cm,相对误差为2.71%,这是因为SEEP/W在降雨入渗分析时土体湿润层与初始层之间存在一个较窄的过渡段,该段即为两个模型计算结果的偏差。

图12 改进模型与数值模拟分析结果对比

3.2 降雨入渗过程分析

3.2.1 弱降雨工况(q≤ks)

采用改进模型(即考虑地下水影响时初始含水率指数分布)与Chen-Young模型(取改进模型坡面体积含水率0.10为初始均匀含水率值)分别计算降雨强度为5、7、9、11 mm/h时湿润锋推移深度与降雨时间的关系,结果如图13所示。初始含水率均匀分布时,湿润锋推移速度随降雨强度增加而加快,且同一降雨强度在入渗过程中其推移速度保持不变;考虑地下水影响即初始含水率呈指数分布时,由于坡体初始含水率沿深度逐渐增大,单位深度土体达到近饱和状态时所需雨量逐渐减小,因此同一降雨强度在入渗过程中湿润锋推移速度逐渐加快。

图13 弱降雨工况下降雨时间与湿润锋推移深度关系曲线

3.2.2 强降雨工况(q>ks)

采用改进模型与Chen-Young模型分别计算降雨强度为13、15、18、21 mm/h时,考虑地下水和不考虑地下水影响时(含水率分布与弱降雨工况相同)湿润锋推移深度与降雨时间的关系,结果如图14所示。初始含水率均匀分布时,坡面产生径流的时间受降雨强度影响明显,坡面产流后不同降雨强度下,降雨时间达到130 h时对应的湿润锋推移深度最大差距仅为2.7 cm,此时,降雨入渗率受土体自身入渗能力控制;考虑地下水影响即初始含水率呈指数分布时,强降雨工况下增加降雨强度导致湿润锋推移速度虽仍有增加,但相比弱降雨工况增加速度变缓。同一降雨强度下湿润锋推移速度呈现先基本保持不变然后逐渐加快的趋势,此时,降雨入渗率受土体自身入渗能力和坡体初始含水率共同控制,在靠近毛细饱和带附近土体时,湿润锋推移速度达到最大。

图14 强降雨工况下降雨时间与湿润锋推移深度关系曲线

利用改进模型计算降雨强度为13、15、18、21、24、27 mm/h时对应的坡面产流临界时间tp和对应的湿润锋推移深度zp。如图15所示,强降雨工况下提高降雨强度会明显缩短坡面产生径流的临界时间,此后一部分降雨将沿坡面流走。

图15 降雨强度与坡面产流临界时间/深度关系

3.3 考虑初始地下水影响的浅层边坡降雨作用下稳定性分析

3.3.1 坡体初始含水率分布对失稳时间和失稳位置的影响

分别分析坡体在初始含水率均匀分布(假定初始体积含水率为0.18)和考虑地下水影响时初始含水率指数分布两种情况下降雨型浅层滑坡失稳时间和失稳位置,当降雨强度为7 mm/h时,浅层边坡湿润锋处和基岩面处的安全系数与降雨时间的关系分别如图16和17所示。根据唐扬等[27]的结论,当安全系数小于1.05时,边坡即将失稳,应进行滑坡灾害风险预警预报。

图16 初始含水率均匀分布时安全系数与降雨时间关系

由图16可知,当降雨时间达到115 h时,湿润锋处安全系数已减小至1.05,此时基岩面处的安全系数为1.14,说明该工况下将在湿润锋处发生失稳;当浅层边坡基岩处土层存在初始地下水时,如图17所示,降雨时间达到77 h,基岩面处的安全系数降低至1.05,此时湿润锋处安全系数为1.15,说明该工况下失稳位置将发生在基岩面处,且基岩面处初始安全系数仅为1.08,即该类边坡在降雨前即处于欠稳定状态。

图17 考虑初始地下水影响的安全系数与降雨时间关系

3.3.2 降雨强度对失稳时间和失稳位置的影响

利用改进模型计算降雨强度为5、7、9、11、13、15、18、21、24、27 mm/h时存在初始地下水的浅层边坡失稳时间和失稳位置,结果如图18所示。

图18 不同降雨强度下失稳时间和失稳位置关系

由图18可以看出,当降雨强度q小于ks时,降雨强度从5 mm/h提高到11 mm/h,边坡触发失稳时间从108 h呈近指数关系快速缩短至45 h;当降雨强度q大于ks时,降雨强度从15 mm/h增大到27 mm/h,失稳时间仅从42 h缩短到30 h,失稳时间缩短速度趋缓,此时大多降雨以径流的形式沿坡面流走,当不考虑坡面径流对入渗的促进作用时,影响边坡稳定性的降雨入渗量趋于稳定;当基岩面处土层存在稳定的初始地下水时,失稳位置往往发生在基岩面处而非传统预警预报模型中的湿润锋处,因此,重点关注该类边坡可以提高滑坡降雨阈值预测的准确性。

4 结 论

1)浙江省山地丘陵地区浅层边坡土层通常存在稳定的初始地下水,本文基于土体非饱和特性,提出了边坡土体体积含水率指数分布模型,基于GA入渗模型及水量平衡原理推导了存在初始地下水的浅层边坡在弱降雨和强降雨两类工况下的入渗过程解析解。

2)将改进模型与传统模型和数值模拟计算结果对比,验证了改进模型的有效性,将其应用于基岩处土层存在初始地下水的浅层边坡降雨入渗及稳定性分析,拓展了传统模型的应用范围,且解析解模型更有利于工程上的应用。

3)当浅层边坡存在初始地下水时,边坡在降雨前即处于欠稳定状态,当降雨发生后,该类边坡失稳位置将发生在基岩面处,此时湿润锋处还偏于安全,且失稳时间提前于传统模型的预测结果,须重点关注。

4)弱降雨工况下,降雨强度提高,边坡失稳启动时间呈近指数关系快速缩短;强降雨工况下,降雨强度提高,边坡失稳启动时间缩短速度趋缓,此时大多数降雨以径流的形式沿坡面流走,当不考虑坡面径流对入渗的促进作用时,影响边坡稳定的降雨入渗量趋于稳定。