李 宁,刘冠麟,许建聪,张利伟,唐诗晴,吴 数

(1.上海理工大学环境与建筑学院,上海 200093；2.同济大学地下建筑与工程系,上海 200092；3. 四川省城乡规划设计研究院,四川 成都 610081)

滑坡已成为仅次于地震的全球第二大自然灾害,其每年在我国造成的经济损失高达100亿元以上[1～2]。统计资料表明,降雨是导致滑坡的最主要和最普遍的环境因素[3]。因此,对降雨条件下的边坡稳定性进行研究具有重要的现实意义。

刘艳辉等[4]、赵国通等[5]采用统计方法对降雨滑坡进行了研究；张泰丽等[6]采用物理模型试验对台风暴雨条件下边坡变形演化规律进行了研究。数值模拟作为一种有效的研究方法也被广泛应用于边坡稳定性的研究中,尤其是有限元强度折减法,因其为边坡稳定系数及其滑面的确定提供了一条有效途径,而被越来越多的学者应用于降雨条件下边坡稳定性的研究中。Cai等[7]采用有限元渗流分析与强度折减法相结合,对降雨作用下水平排水对边坡稳定性的影响进行了研究。随后,Cai等[8]又采用相同的方法对降雨对边坡稳定性的影响进行了较系统的研究。于玉贞等[9]将有限元渗流分析与基于Bishop非饱和土抗剪强度公式的强度折减有限元法相结合,对Okura等[10]试验边坡的稳定性进行了评价,以此来探讨将边坡稳定分析应用于滑坡发生的时机预测可能面临的问题。荣冠等[11]根据饱和非饱和降雨入渗有限元模拟,并结合有限元强度折减法对元磨公路K261段典型滑坡体工程边坡,进行了降雨过程边坡变形情况及稳定性变化分析。林鸿州等[12]利用降雨条件下的土坡渗流有限元分析与强度折减有限元法相结合,对滑坡模型试验进行了数值分析。王协群等[13]同时采用强度折减法和极限平衡法中的Bishop法,分析了不同降雨入渗深度对非饱和路堤稳定性和滑弧变化特征的影响规律。

但目前采用有限元强度折减法确定边坡稳定系数时,往往要根据不同的折减系数,不断进行试算,过程较繁琐[14～15]。尤其是对于降雨滑坡的计算,对于每一时间步,都需要进行反复试算,才能得到这一时间步对应的稳定系数。而要得到边坡稳定性随降雨历时的变化过程,就必须得到每一时间步所对应的稳定系数,其计算量是非常大的,同时,对降雨边界的处理更加大了降雨滑坡在有限元分析计算方面的复杂性。因此,目前从计算角度来看,降雨滑坡的强度折减法有限元分析还比较繁琐,其计算工作量巨大[16],这也在一定程度上阻碍了该方法在降雨诱发滑坡中的推广应用。

为了能够将有限元强度折减法更加准确、高效地应用于降雨滑坡计算中,本文以ABAQUS软件为基础开发平台,首先基于GeoStudio软件采用迭代算法来对降雨边界进行处理,从而克服ABAQUS软件不能合理考虑降雨边界为不定边界的情况[17],然后根据总应力平衡方程并结合Lu[18]和Likos[19]提出的饱和-非饱和土的普遍有效应力公式,采用Fortran编制程序来对ABAQUS中的有限元节点力进行修正,最后结合改进的有限元强度折减法[20],基于ABAQUS提供的场变量对其进行二次开发,从而建立降雨条件下边坡有限元强度折减法计算平台。然后结合典型边坡案例,对其在降雨条件下的稳定性进行计算,并将其与传统极限平衡法得到的结果进行比较,从而对建立的计算平台的可行性进行验证。最后采用该平台对降雨滑坡机理进行研究,研究结果可为边坡防减灾提供参考。

1 基于ABAQUS软件降雨诱发滑坡有限元强度折减法计算平台开发

1.1 饱和-非饱和渗流控制方程

根据水体质量守恒,模拟降雨入渗的饱和-非饱和渗流控制方程可表示为如下形式[21]：

(1)

式中：c(h)——容水度；

β——参数,非饱和区β=0,饱和区β=1；

Ss——单位储水系数,对非饱和体,Ss=0,对饱和体Ss为常数；

t——时间变量；

kr(h)——相对渗透系数；

xi，xj——坐标；

h——压力水头。

1.2 降雨入渗边界的处理

对于降雨入渗边界,以往的数值模拟中多采用定边界进行处理[22],而朱伟等[23]与陈学东[24]通过土柱降雨入渗试验,指出降雨入渗边界是随降雨时间而不断变化的动边界,其具体形式要根据地表径流情况来确定。

(1)未出现径流时

此时降雨将全部入渗,因此其边界流量应满足：

(2)

式中：ni——坡面法线向量；

p——降雨强度。

其边界上的水头应满足：

h≤0

(3)

(2)出现径流时

此时降雨不能全部入渗,未入渗部分将在坡面形成径流,由于坡面具有一定坡度,假定径流可迅速流走而未形成积水,则此时边界上水头应满足：

h=0

(4)

其边界上的流量应满足：

(5)

采用迭代算法对以上降雨入渗边界进行处理：

(1)首先将上一时步降雨边界所处的状态(有无径流)作为下一时步降雨边界迭代的初始条件。

(2)假如处于无径流状态,则首先以式(2)为定解条件进行计算；然后将得到的边界上节点的水头计算结果代入式(3)进行校核,如全部满足则转入下一时步计算；若不满足则转入步骤(3)。

(3)假如处于有径流状态,则首先以式(4)作为定解条件进行计算；然后将得到的边界上节点的流量计算结果代入式(5)进行校核,如全部满足则转入下一时步计算；若不满足则转入步骤(2)。

(4)重复步骤(2)、(3),直至降雨边界上所有节点均满足相应的定解及校核条件。

1.3 水压形成的等效节点力

由饱和-非饱和渗流计算得到的水压,根据Cai等[8]提出的方法,可求得相应的等效节点荷载：

(6)

式中：[B]——应变-位移矩阵,{I1}T={1,1,1,0,0,0}；

uw——水压。

χ可以按照Lu[18]和likos[19]提出的饱和、非饱和相统一的土体有效应力公式选取：当土饱和时取为1,当土体非饱和时取为相对饱和度Se。

1.4 基于ABAQUS软件降雨诱发滑坡有限元强度折减法计算平台开发

由于目前ABAQUS软件将降雨边界作为定边界处理(流量边界)[17],并不能实现1.2节中所述的迭代型降雨入渗边界,因此首先采用Geostudio进行迭代边界的处理,随后利用ABAQUS提供的二次开发功能实现降雨过程中强度折减法的自动完成,从而简化降雨滑坡强度折减法的计算过程。其具体开发思路如下：

(1)采用Geostudio/Seep,通过迭代方法来确定降雨边界,进行降雨条件下的非稳定渗流分析,得到降雨过程中每一时刻的边坡水压分布。

(2)根据(1)中得到的边坡水压分布,根据式(6)利用Fortran编制程序求得由水压形成的等效结点力fu。

(3)采用Fortran编制接口程序,将(1)中Geostudio/Seep渗流计算的有限元信息及(2)中得到的水压形成的等效结点力fu导入ABAQUS软件中。

(4)根据李宁等[20]提出的基于场变量有限元强度折减法,采用USDFLD.f对ABAQUS进行二次开发,对降雨边坡进行强度折减有限元分析,得到降雨过程中每一时刻的边坡稳定系数。其计算流程见图1。

图1 计算流程Fig.1 Calculation process

2 计算平台的可行性验证

为了检验所开发的计算平台的可行性,利用该计算平台对Cai等[8]文中的典型边坡在降雨条件下的稳定性进行计算,并将计算结果与采用传统的刚体极限平衡法的计算结果进行比较,从而对开发的计算平台可行性进行验证。

2.1 计算条件及材料参数

采用文献[8]中的匀质土坡模型进行分析,边坡坡度为1∶1.5,坡高为10 m,几何模型尺寸及有限元模型网格划分见图2。其中,AF、EF边界为不透水边界,AB、BC为水头边界,CD、DE为降雨入渗边界,施加的降雨强度为10 mm/h,持续72 h。此外,AB、EF为水平约束边界,AF为固定约束边界,其他为自由边界。初始时刻,地下水位与BC面等高,地表相对饱和度Se=0.72,并随高程线性变化至地下水位处。

图2 几何参数及有限元网格Fig.2 Geometric parameters and finite element mesh

边坡土体参数亦按文献[8]进行选取,其中非饱和土水力特性采用Van Genuchten(VG)模型[25]来描述,具体见表1所示。

表1 土体材料参数Table 1 Soil material parameters

2.2 与刚体极限平衡法的对比分析

将计算平台得到的结果与采用Geostudio进行刚体极限平衡法计算得到的结果进行对比分析,这里只列举了降雨结束时段(降雨72 h)时的对比结果,其他各降雨时段的对比结果与此类似,限于文章篇幅,在此不再列出。

采用本文开发的计算平台进行计算,可得持续降雨72 h后的孔压水头分布图(图3),由图3可知,随着降雨的进行,地下水位将不断抬高,土坡上层的非饱和区将不断减少。根据图3中得到的水压,利用该计算平台可进一步确定此时的稳定系数为0.864；同时还可得到降雨72 h后边坡破坏时的等效塑性应变云图,据此可以确定滑面的位置(图4)。

图3 持续降雨72 h后孔压水头分布图(m)Fig.3 Distribution of pore pressure head after 72 hours of continuous rainfall(m)

图4 降雨72 h后边坡破坏时的等效塑性应变云图Fig.4 Equivalent plastic strain cloud image of slope failure after 72h rainfall

为了验证本文开发的计算平台的可靠性,采用传统刚体极限平衡法中的简化Bishop法和Morgenstern-Price法,采用Geostudio分别对降雨72 h后的边坡稳定性进行计算,可以得到边坡稳定系数为0.854和0.857；同时还可得到相应的滑面(图5和图6)。

图5 Bishop法得到的稳定系数与相应的滑面Fig.5 The safety factor obtained by Bishop method and the corresponding slip surface

图6 Morgernstern-Price法得到的稳定系数与相应的滑面Fig.6 The safety factor obtained by Morgernstern-Price method and the corresponding slip surface

将两者得到的结果进行比较可以看出,本文开发的计算平台得到的降雨72 h后边坡的稳定系数为0.864；采用传统的简化Bishop法和Morgerstern-Price法得到的边坡稳定系数分别为0.854与0.857,分别比计算平台得到的结果小1.2%和0.8%。但总体来看,它们之间的误差并不大,这说明本文开发的计算平台与传统极限平衡法的计算结果具有较好的一致性。从图4、图5、图6可以看出,计算平台所确定的滑面位置与刚体极限平衡法所得到的滑面位置也是基本一致的。

通过对比分析可以看出,采用开发的计算平台得到的滑面位置与传统方法确定的滑面位置是一致的,稳定系数的相对误差最大为1.2%,说明所开发的计算平台是稳定、可靠的,可以用于降雨诱发滑坡的计算。

3 降雨条件下土坡破坏机理及破坏过程研究

采用本文开发的计算平台,仍以Cai等[8]文中的典型边坡为基础,进一步对其在降雨条件下的破坏过程及破坏机理进行研究。

采用开发的计算平台,可以得到稳定系数随降雨时间的变化图(图7),从图7中可以看出,随着降雨时间的增加,边坡的稳定性将不断降低；当降雨时间持续24 h后,边坡的稳定系数为0.997,此时边坡发生破坏。以下将进一步对降雨过程中该边坡的孔隙水压力、有效应力、等效节点力及塑性应变的变化情况进行研究,从而阐明该边坡在降雨条件下的破坏机理及破坏过程。

图7 稳定系数随降雨时间变化Fig.7 Variation of safety factor with rainfall

3.1 孔隙水压力变化情况

图8为降雨各时段边坡内孔压水头的分布,由此可以看出,随着降雨时间的增加,边坡内的浸润线逐渐升高,坡体内的饱和区不断扩大,非饱和区不断减小,饱和区的水压不断增加；由于初始时刻地表的相对饱和度Se=0.720,并随高程线性变化到地下水位处；所以边坡非饱和区的水压较小,其随降雨时间的变化也没有饱和区明显,在图中显示的非饱和区的水压变化并不明显,所以这里选取非饱和区中的A1点，A2点(图8c)为代表,对非饱和区内的水压随降雨时间的变化进行研究。

图8 各降雨时段孔压水头分布图(m)Fig.8 Distribution of pore pressure head at various rainfall time (m)

图9显示了边坡非饱和区中的A1，A2两点的孔隙水压力随降雨时间的变化过程,从图9中可以看出,随着降雨时间的增长,雨水不断入渗进入边坡,非饱和区的基质吸力不断减小。

图9 非饱和区水压随降雨时间的变化Fig.9 Variation of water pressure in unsaturated zone with rainfall time

3.2 有效应力变化情况

图10为降雨各时段边坡内竖向有效正应力的分布,由此可以看出,在整个降雨过程中,边坡土体均处于受压状态,且其分布规律基本上是相同的；但是随着降雨时间的增长,有效应力将逐渐减小,这一点可以通过前面对降雨过程中水压变化规律的分析来加以解释：伴随着降雨过程,边坡内的浸润线逐渐升高,饱和区的水压不断升高,而非饱和区的基质吸力不断减小,从而造成边坡内有效应力逐渐减小。而有效应力的减小,又将造成土体抗剪强度的减小,因此在整个降雨过程中,土体的抗剪强度将逐渐降低。

图10 各降雨时段竖向有效正应力云图(kPa)Fig.10 Vertical Effective Normal Stress at various rainfall time (kPa)

3.3 等效节点力变化情况

图11为坡脚处水平向等效节点力随降雨时间的变化,其中B1点的具体位置如图8c所示,水平等效结点力以朝向坡内为正。由图11可知,在整个降雨过程中,坡脚处受到朝向边坡临空面方向的水平渗流力作用,且其值随着降雨时间的增长而不断增大。这主要是因为随着降雨的进行,边坡地下水位不断抬升,边坡后缘与坡脚处的水头差越来越大,造成坡脚处的水力梯度也越来越大,所以坡脚处的水平向渗流力也越来越大；由于边坡后缘的水位高于坡脚处水位,因此雨水朝向坡脚处流动,从而造成坡脚处受到朝向边坡临空面方向的水平渗流力作用,这对于边坡稳定性是不利的。

图11 坡脚处水平向等效节点力随降雨时间的变化(kN)Fig.11 Changes in equivalent nodal force at the toe of the slope with rainfall time

3.4 塑性应变的变化情况

图12显示了降雨作用下边坡塑性区的发展过程,由此可以看出,在降雨初期,坡脚部位首先出现塑性区,随着降雨的进行,坡脚处的塑性区进一步向坡顶发展,在降雨24 h后,整个塑性区从坡脚到坡顶贯通,边坡发生破坏,此时的稳定系数为0.997。

图12 各降雨时段边坡塑性应变云图Fig.12 Plastic strain of slope at various rainfall time

图13 各降雨时段边坡内流速矢量图Fig.13 Velocity vector in the slope at various rainfall time

图13分别为降雨2 h和24 h后的边坡内的流速矢量图,从图13可以看出,坡体内部的雨水通过向坡脚处汇集进而排出坡体,因而坡脚处的流速较大,也即是说在坡脚处具有较高的水力梯度,水压的变化率较大,因此产生的渗流作用力也较强,所以在坡脚处首先出现塑性区；随着降雨时间的增长,坡体内部的雨水不断增加,地下水位不断抬升,坡脚处的渗流作用力也进一步增加,而坡脚处土体已处于塑性状态,其抗剪强度已得到充分发挥,无法继续承担不断增加的渗流作用力,因而会进一步牵动周围土体来发挥抗剪强度,所以塑性区会进一步向上发展,使更多的土体发挥抗剪强度，直到最后塑性区从坡脚发展到坡顶,塑性区贯通,造成边坡的整体破坏。

3.5 降雨作用下边坡破坏机理及破坏过程分析

在降雨作用下,边坡内的浸润线逐渐抬高,坡体内的饱和区逐渐扩大,非饱和区逐渐缩小,饱和区的水压不断升高,非饱和区的基质吸力不断减小,从而使坡体内的有效应力减小,导致土体的抗剪强度降低。由于边坡内地下水位不断抬升,使边坡后缘的地下水位高于坡脚处水位,雨水朝向坡脚处流动,对坡脚形成朝向边坡临空面方向的水平渗流力作用,且其值随降雨时间不断增加。土体的抗剪强度不断减小,而朝向边坡临空面方向的水平渗流力作用不断加大,从而导致了边坡的最终破坏。

在降雨作用下,坡体内部的雨水通过向坡脚处汇集进而排出坡体,因而坡脚处具有较高的水力梯度,水压的变化率较大,因此产生的渗流作用力也较强,所以在坡脚处首先出现塑性区。随着降雨时间的增长,地下水位不断抬升,使边坡后缘与坡脚处的水头差越来越大,从而使坡脚处的渗流作用力进一步增加,而坡脚处土体已处于塑性状态,其抗剪强度已得到充分发挥,无法继续承担不断增加的渗流作用力,因而会进一步牵动周围土体来发挥抗剪强度,所以塑性区会进一步向上发展,使更多的土体发挥抗剪强度。直到最后塑性区从坡脚发展到坡顶,塑性区贯通,形成从坡顶到坡脚的一个贯通的滑动面,从而造成边坡的整体破坏。

香港是一个降雨滑坡频发的地区,GEO (Geotechnical Engineering Office) 通过对香港地区大量降雨滑坡实例的调查发现,大部分降雨滑坡首先是在坡脚处发生破坏,然后形成整体破坏[26]。采用开发的计算平台所得到的边坡破坏过程与这一实际观察到的破坏过程是一致的,这也从侧面说明该计算平台可以用于降雨滑坡破坏机理及破坏过程的研究。

4 结论

(1)以ABAQUS软件为基础开发平台,通过二次开发,建立了降雨条件下边坡有限元强度折减法计算平台。该计算平台较合理地考虑了降雨边界为不定边界的情况,并简化了有限元强度折减法的计算过程。通过与传统极限平衡法进行比较,证明所开发的计算平台是稳定、可靠的,可以用于降雨诱发滑坡的计算。

(2)通过该计算平台对案例边坡在降雨作用下的破坏机理进行研究,发现该边坡在降雨作用下的破坏主要由两方面因素造成,其一：降雨导致地下水位抬升,从而使土体抗剪强度减小；其二：由于边坡内地下水位不断抬升,从而对坡脚形成朝向边坡临空面方向的水平渗流力作用,且其值随降雨时间不断增加。这两个因素的共同作用,最终导致边坡的破坏。

(3)通过该计算平台对案例边坡在降雨作用下的破坏过程进行研究,发现该边坡在降雨条件下的破坏具有典型的牵引式滑坡特点：在降雨作用下,坡脚处首先发生破坏,并不断牵动邻近的上部土体发生破坏,最终形成从坡顶到坡脚的一个贯通的滑动面,从而造成边坡的整体破坏。该破坏过程与香港地区实际观察到的降雨条件下边坡的破坏过程是一致的。

参考文献：

[1] Dai F C, Lee C F, Ngai Y Y. Landslide risk assessment and management: an overview[J]. Engineering Geology, 2002, 64(1): 65-87.

[2] 刘新喜,晏鄂川,王鹏飞,等.滑坡稳定性评价的非线性方法[J]. 中国安全科学学报,2003,13(1):34-36.[LIU X X, YAN E C, WANG P F,etal. Nonlinear method for evaluating the landslide stability[J]. China Safety Science Journal, 2003, 13(1): 34-36. (in Chinese)

[3] Ng C W, Shi Q. A numerical investigation of the stability of unstaturated soil slopes subjected to transient seepage[J].Computer and Geotechanics, 1998, 22(1):1-28.

[4] 刘艳辉,方志伟,温铭生,等.川东北地区强降雨诱发崩滑流灾害分析[J]. 水文地质工程地质,2014,41(2): 111-115.[LIU Y H, FANG Z W, WEN M S,etal. Analysis of geo-hazards and rainfall in northeast of Sichuan Province, China[J]. Hydrogeology & Engineering Geology, 2014, 41(2):111-115. (in Chinese)

[5] 赵国通,朱小龙,马佰衡.河北省崩塌滑坡滞后于降雨的地质背景分析[J]. 水文地质工程地质,2015,42(2):151-155.[ZHAO G T, ZHU X L, MA B H. Geological background analysis of the collapse and landslide lagging behind rainfall in Hebei Province[J]. Hydrogeology & Engineering Geology, 2015, 42(2): 151-155. (in Chinese)

[6] 张泰丽,周爱国,施斌,等.台风暴雨条件下滑坡变形特征物理试验研究[J]. 水文地质工程地质,2016,43(6): 127-132.[ZHANG T L, ZHOU A G, SHI B,etal. Physical experiment research on landslide deformation characteristics under the condition of the typhoon heavy rain[J]. Hydrogeology & Engineering Geology, 2016, 43(6):127-132. (in Chinese)

[7] Cai F, Ugai K, Wakai A,etal. Effect of horizontal drains on slope stability under rainfall by three-dimensional finite element analysis[J]. Computers and Geotechnics, 1998, 23:255-275.

[8] Cai F, Ugai K. Numerical analysis of rainfall effects on slope stability[J]. International Journal of Geomechanics, 2004, 4(2):69-78.

[9] 于玉贞,林鸿州,李广信.边坡滑动预测的有限元分析[J]. 岩土工程学报,2007,29(8):1264-1267.[YU Y Z, LIN H Z, LI G X. Analysis of finite element method for early warning of landslide[J]. Chinese Journal of Geotechnical Engineering, 2007, 29 (8):1264-1267. (in Chinese)

[10] Okura Y, Kitahapa H, Ochiai H,etal. Landslide fluidization process by flume experiments[J]. Engineering Geology, 2002, 66(1): 65-78.

[11] 荣冠,王思敬,王恩志,等.强降雨下元磨公路典型工程边坡稳定性研究[J].岩石力学与工程学报,2008, 27(4): 704-711. [RONG G, WANG S J, WANG E Z,etal. Stability analysis of Yuanmo highway slope under intensive rainfall[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(4), 704-711.

[12] 林鸿州,于玉贞,李广信,等.强度折减有限元法在滑坡特性预测的应用探讨[J]. 岩土工程学报,2009,31(2):229-233.[LIN H Z, YU Y Z, LI G X,etal. ON Application of soli-water characteristic curves to landslide forecast[J]. Chinese Journal of Geotechnical Engineering,2009, 28(12): 2569-2576. (in Chinese)]

[13] 王协群,张有祥,邹维列,等.降雨入渗条件下非饱和路堤变形与边坡的稳定数值模拟[J].岩土力学,2010,31(11): 3640-3655.[WANG X Q, ZHANG Y X, ZOU W L,etal. Numerical simulation for unsaturated road-embankment deformation and slope stability under rainfall infiltration[J].Rock and Soil Mechanics, 2010,31(11):3640-3655. (in Chinese)]

[14] Dawson E M, Roth W H, Drescher A. Slope stability analysis by strength reduction[J]. Geotechnique, 1999, 49(6): 835-840.

[15] Xu Q J, Yin H L, Cao X F,etal. A temperature-driven strength reduction method for slope stability analysis[J]. Mechanics Research Communications, 2009, 36:224-231.

[16] 贾苍琴.考虑非饱和非稳定渗流和剪胀性的土坡稳定分析[D].上海：同济大学,2006.[JIA C Q. Soil slope stability analysis with consideration of the saturated unsaturated seepage and dialatancy characteristic of the soil [D]. Shanghai: Tongji University, 2006. (in Chinese)]

[17] 章正,张本卓. 应用ABAQUS模拟边坡降雨入渗的评价性分析[J]. 中国科技信息,2012, 89:893-900. [ZHANG Z, ZHANG B Z. Evaluation of ABAQUS’s performance in simulating rainfall infiltration into slope [J]. China Science and Technology Information, 2012, 89: 893-900.

[18] Lu N, Likos W J. Unsaturated Soil Mechanics [M]. New York：John Wiley and Sons,2004: 132-133.

[19] Likos W J, Lu N. Suction stress characteristic curve for unsaturated soil [J]. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(2), 131-142.

[20] 李宁,许建聪.基于场变量的边坡稳定分析有限元强度折减法[J].岩土力学,2012, 33(1): 314-318.[LI N, XU J C. Strength reduction FEM for slope stability analysis based on field variable[J]. Rock and Soil Mechanics, 2012, 33(1): 314-318. (in Chinese)]

[21] Neuman S P. Saturated-unsaturated seepage by finite elements[C]//Proceedings of ASCE Hydrology, 1973, 99(12):2233-2250.

[22] Li J H, Du L, Chen R,etal. Numerical investigation of the performance of covers with capillary barrier effects in south china [J]. Computers and Geotechnics, 2013, 48: 304-315.

[23] 朱伟,程南军,陈学东,等.浅谈非饱和渗流的几个基本问题[J].岩土工程学报,2006, 28(2):235-240.[ZHU W, CHENG N J, CHEN X D,etal. Some fundamental problems of unsaturated seepage[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(2): 235-240. (in Chinese)]

[24] 陈学东.浅层非饱和带降雨入渗规律的试验与数值研究[D]. 南京：河海大学,2005.[CHEN X D. Experiments and numerical study on rainfall infiltration into shallow unsaturated soils [D].Nanjing: Hohai University, 2005. (in Chinese)]

[25] Van GENUCHTEN M T H. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils [J]. Soil Science Society of American Journal, 1980, 44(5):892-898.

[26] Geotechnical Engineering Office. Geotechnical Manual for Slopes[M]. Hongkong: The government of the Hong Kong special Administrative Region, 2011.