考虑水文效应的不同植被根系护坡效果数值分析

廖田婷，蒋关鲁，张崇磊，2，单源椿，叶雄威

(1.西南交通大学 土木工程学院，成都 610031； 2.中国科学院水利部成都山地灾害与环境研究所，成都 610041)

1 研究背景

降雨诱发边坡失稳是工程中常见的重要地质灾害之一[1-2]。生态护坡是目前国际前沿的研究焦点，植被对边坡的加固作用已被工程界广泛认识[3]。植被根系固土的力学效应体现在根土界面的摩擦力使根类似于抗滑桩或锚杆的作用，形成根-土复合体，提高土体的抗剪强度；植被冠层具有良好的保水功能，可以有效拦截雨水，减小坡体浅层冲刷；植被根系除加筋作用外，还具有一定的蒸腾作用，能增大土体的基质吸力，增强坡体强度[4]。

以往学者多集中于针对不同土体通过室内模型试验来探讨和总结降雨后坡体内部水分运移规律与物理力学性质的变化[5-8]，但针对植被根系固土效果研究的室内模型试验还较少[9-11]，原因是真实生态根系难以模拟。不少学者采用简化根系，Büsgen等[12]根据根的几何形状和分布特征将三维形态树根分为最具代表性的heart型、plate型和tap型3类。Ng等[13]采用与真实根力学性质相近的高聚物空心管来模拟不同类型的根系，还通过抽气系统来模拟蒸腾作用。除此之外，还有学者选择菩提木棒和3D打印ABS塑料来模拟真实根系[10-11，14]。室内模型试验难以充分考虑每种工况以及不同植被根系，数值分析软件的渗流模块和应力分析模块则解决了这一问题。王一冰等[15]通过VADOSE/W软件分析了植筋带(延长根)加固的边坡和普通植被边坡在降雨条件下的水力特性，但未能从应力场来分析其加固效果。

本文基于FLAC3D软件对饱和-非饱和渗流模块和应力分析模块进行二次开发，在现有学者研究的基础上提出了一种新颖的建立初始渗流场的方法，并对雨水入渗引起有效应力和土体抗剪强度的变化值进行了修正[16]，使用了简化方法来模拟植被根系的加筋和蒸腾作用。为验证本程序的合理性，将素坡(未加固边坡)与tap根型、plate根型和heart根型3种植被类型根系加固边坡在长达46.872 d雨季的大暴雨工况下的数值模拟结果与降雨离心模型室内试验进行对比，探讨不同植被类型根系加固边坡的坡体内部水分运移规律和浅层土体加固效果。

2 基于FLAC3D二次开发非饱和降雨入渗及应力修正模块

2.1 FLAC3D渗流计算基本方程

FLAC3D采用达西渗透定律来表征流体在固体颗粒间的流动准则，即

(1)

对于小变形，液体质量平衡方程为

(2)

式中：qi，i为通量强度(L/s)；qv为流体体积源强度；ζ为含水率的变化；t为时间(s)。

FLAC3D的数值计算方法依赖于流体连续性方程的有限差分节点公式，孔隙压力的变化通过区域分配给节点。FLAC3D在进行渗流计算时，默认为饱和渗流，即饱和度等于1，可通过设置流体抗拉强度使得负孔隙水压存在，也就是说，饱和渗流是通过有限差分迭代节点孔隙水压力来进行流固耦合计算的，即

(3)

式中：M为比奥模量；α为与土体压缩性相关的系数；ε为体积应变。

在非饱和渗流情况下，若想设置饱和度不为1，FLAC3D会强制节点孔隙水压为0，但实际上可通过土-水特征曲线反算出节点孔隙水压，此时FLAC3D是通过有限差分迭代节点饱和度来进行非饱和渗流计算的，反算出的节点孔隙水压并未影响这一时刻的内部渗流计算，即

(4)

式中n为孔隙率。

饱和渗流可以看作饱和度为1的非饱和渗流，上述两种算法均使得运用FLAC3D进行非饱和渗流计算成为可能。

2.2 饱和-非饱和渗流应力及抗剪强度修正模块

Bishop非饱和土有效应力原理为：

σ′=σ-ua+Sr(ua-uw) ；

(5)

σ′=σ-ua+1(ua-uw)-(1-Sr)(ua-uw)。

(6)

式中：σ为总应力(kPa)；ua为孔隙气压力(kPa)；Sr为非饱和土饱和度；uw为孔隙水压力(kPa)；(ua-uw)为基质吸力，一般取正值；σ′为有效应力。

当Sr=1时，即为饱和土的有效应力原理(Terzaghi有效应力原理)[17]，FLAC3D使用Terzaghi有效应力原理进行计算，仅考虑饱和度为1，在进行抗剪强度计算前，采用式(6)进行内置有效应力模块的修正，只需要在总应力上减去(1-Sr)(ua-uw)项。

根据弗雷德隆德等[18]的非饱和土抗剪强度理论有

τ=c′+(σ-ua)ftanφ′+(ua-uw)ftanφb。

(7)

式中：τ为抗剪强度(kPa)；c′为有效黏聚力(kPa)；(σ-ua)f为破坏时在破裂面上的净法向应力(kPa)；φ′有效内摩擦角(°)；(ua-uw)f为破坏时在破裂面上的基质吸力(kPa)；φb为抗剪强度随基质吸力变化的内摩擦角(°)(本文由试验获得)。

FLAC3D在进行应力计算时，未考虑非饱和土体的黏聚力随含水率变化值，需采用式(7)进行黏聚力的修正。在每个力学计算时步，都使用c′+(ua-uw)ftanφb替换初始值，实现实时修正。

2.3 饱和非饱和渗流计算模块的二次开发

边坡降雨入渗是一个非饱和到饱和的过程，降雨强度与土壤入渗能力关系可分为“降水模型”和“积水模型”。“降水模型”中，降雨强度小于边坡表层土体单元的饱和渗透系数，入渗量取决于降雨强度；“积水模型”中，降雨强度大于边坡表层土体单元的饱和渗透系数，入渗量取决于饱和渗透系数，且坡表产生径流。为实现上述两种模型，设置动态边界条件，地表产流前，施加稳定流量边界，浅层土体饱和后，地表雨水产生径流，施加压力边界，即取表面压力为0。

本文为实现降雨入渗饱和-非饱和渗流，总结出两个关键技术。FLAC3D在进行渗流计算时均采用饱和渗透系数，非饱和渗流计算关键在于实时更新非饱和土的渗透系数，这是本文二次开发的关键技术之一。生成正确的初始渗流是水力特性的重要前提，大多数学者均通过设置流体抗拉强度和地下水位线自动生成负孔隙水压力，与真实基质吸力差距较大。即使经过稳态渗流，浅层水分在地下水和降雨的补给条件下，基质吸力难以与软件自动生成的基质吸力相近。部分学者利用FISH函数直接赋予正确的节点孔隙水压力，但不开启渗流模式，降雨初期内部渗流和降雨入渗同时进行，基质吸力变化幅度大。本文的另一个关键技术在于开发了一种简便的基于土-水特征曲线利用有限差分迭代节点饱和度来生成正确的初始渗流场的方法。

对于非饱和土，根据Van Genuchten模型[19]，饱和度和基质吸力的关系式为

(8)

式中：sr为残余饱和度；a、n、m为拟合参数。

非饱和土渗透系数表达式为

(9)

式中：ks为饱和渗透系数；k为非饱和渗透系数。

利用FISH语言在FLAC3D中实现基于正确初始渗流场的降雨入渗非饱和到饱和渗流步骤(图1)如下：

(1)初始渗流平衡阶段，打开渗流模式，关闭力学计算进程。使用initialize命令初始化设置土体每部分的节点饱和度，自定义FISH语言，获取每时刻节点饱和度，通过土-水特征曲线反算出节点孔隙水压力，赋给节点内置变量gp_pp(pnt)。利用8个

图1 初始渗流平衡-降雨-固结阶段步骤Fig.1 Steps of initial seepage balance stage，rainfall stage，and consolidation stage

节点饱和度体积平均计算单元饱和度，进一步求得相对介质渗透率和单元渗透系数，再赋给该单元，利用自带FISHCALL命令实现实时更新节点孔隙水压力和单元渗透系数。

(2)降雨阶段。利用FLAC3D的节点渗流边界条件apply pwell施加降雨边界，一旦边坡表层节点达到饱和状态，则更改边界条件为压力边界，即固定孔隙水压力为0。在上述获取正确的初始渗流场后，彻底删除上述FISH，在本步骤中采用有限差分迭代节点孔隙水压力来实时更新单元渗透系数。首先设置节点饱和度为1，使用gp_pp(pnt)提取通过饱和度反算求得的正确的节点孔隙水压力，利用土-水特征曲线求得节点饱和度(注意这里的节点饱和度不可覆盖原有节点饱和度)，可以使用gp_extra(pnt)存储节点饱和度，利用8个节点饱和度体积平均求得单元饱和度，进一步求得单元渗透系数，再实现实时更新单元渗透系数。

(3)固结阶段。水分在重力的作用下继续向坡体内部渗流，而在边坡表面，坡体内的水分既可以向边界外溢出，也可以在重力作用下运移[20]。因此设置表层节点最大孔隙水压力为0，并且阻止坡体内部水分不断往外溢出。

2.4 FLAC3D模拟根系固坡原理

FLAC3D内置线弹性本构模型和设置平面的渗流边界条件为实现有效模拟根系固土的力学效应和蒸腾作用奠定了良好的基础。考虑到植被根系的复杂性，为简化计算，在构建模型时设置土体为摩尔-库伦(Mohr-Coulomb)模型，对根系另行分组并赋予线弹性模型，模拟根系的锚固作用，设置相应的参数并进行力学计算。为模拟植被根系的蒸腾作用，在渗流阶段，采用apply discharge命令设置根单元表面为渗漏边界，在整个降雨及雨后固结过程中，坡体内水分会不断从根系与土体的接触面部分往外溢出。需要注意的是，实际情况中植被根系为不透水或微透水，但在若设置根系为不透水模型，根系所在位置会出现异常的堵截水行为，这与实际不符，因此设置根系仍为透水模型，忽略植被根系的微弱拦水行为。

3 数值计算模型及离心模型试验对比

3.1 降雨入渗离心模型试验简介

为验证本数值模拟计算程序的正确性和可靠性，将结合降雨离心模型室内试验的结果进行对比分析。

离心模型设计尺寸为0.7 m×0.6 m×0.41 m，以相似比尺15换算原型尺寸为10.5 m×9 m×6.15 m的边坡，采用自主研发且适用于本模型箱的降雨系统(图2)。为模拟长达46.872 d的间隔降雨大暴雨(24 h降雨量145.5 mm)工况，实施间隔降雨，设计参数见表1。

图2 模型箱及降雨装置Fig.2 Model box and rainfall device

表1 降雨工况

将ISO标准砂和伊利石粉按质量比8∶2混合作为试验用土，通过一定的缩尺对比分析，最终采用PA管(尼龙)制作模型根，并在模型根周围放置空心高进气值陶土头。外部抽气系统利用内外压强差将只吸水不进气的陶土头中的水分用导管抽出，从而模拟植被根系的蒸腾作用。具体模型的设计参数如表2所示(根-土接触面参数由直剪试验测得)。

表2 模型设计参数

3.2 计算模型

数值模拟计算与室内降雨离心模型试验按相似比尺15换算的原型相同。边坡所采用的土-水特征曲线根据试验所测数据使用Van Genuchten模型(简称VG模型)进行拟合(如式(8)、式(9)所示)，结合非饱和土三轴试验的结果进行抗剪强度以及有效应力的修正(如式(6)、式(7)所示)，对应拟合土-水特征曲线参数Sr、a、n分别为0.065 6、6.405 02 kPa、1.622 67，非饱和土抗剪强度参数c′、φ′、φb分别为6.5 kPa、28.11°、5°。为与室内降雨离心模型试验进行对比，在坡体内部布置张力计的位置布置测点T1、T0和T6，而测点H9、H3和H4布置位置与含水率计相同。具体边坡尺寸和监测点布置如图3所示。

图3 边坡尺寸及监测点布置Fig.3 Slope dimensions and monitoring point layout

3.3 素坡与heart根型加固边坡水分运移规律

图4为素坡与heart根型加固边坡两阶段降雨后和固结后的饱和度等值线云图。以0.4作为初始饱和度，则可明显观察到湿润锋的变化路径。降雨阶段，素坡(图4(a))呈现出从浅层逐渐向坡体内部入渗的形态。由非饱和土渗透定律可知，土体渗透系数大小与饱和度呈正相关。上部暂态饱水区水分在重力作用下逐渐下移的过程中表现为近坡表侧下渗速度快于远坡表侧。固结阶段，随着水分逐渐下移，湿润锋不断推进，坡顶浅层再次出现干燥区。坡脚处在近坡表侧和右侧边坡底部水分下渗的双重作用下，水位上升最快。

heart根型加固边坡(图4(b))的水分运移路径与素坡截然不同。降雨阶段，根系蒸腾作用的影响使得坡顶远坡表侧雨水入渗量高于近坡表侧(近根侧)，根系附近出现饱和度<0.4的干燥区。坡体内部远坡表侧渗透系数大于近坡表侧，水分下移速率高于近坡表侧。固结阶段，远坡表侧逐渐湿润，且坡后水分增多。整个过程中，根系附近出现干燥区且不消散，说明植被根系能通过自身体积阻挡部分雨水入渗和蒸腾作用(吸收部分水分)这两方面来影响内部水分运移规律。

图4 素坡与heart根型加固边坡降雨后和固结后的饱和度等值线云图Fig.4 Saturation contours of plain slope and heart-root slope after rainfall and consolidation

3.4 4种工况初始和降雨结束渗流场

基于修正初始渗流场的饱和-非饱和渗流分析模块，经过10 d稳态渗流后，4种工况的初始渗流场孔隙水压力分布如图5(a)所示。由于不同根型的比表面积的差异，在设置等值流量溢出边界的前提下，显示出不同的蒸腾作用效果。坡面顶部的水分在蒸腾作用和重力作用下向下迁移，边坡底部水分不断增加，逐渐分层。以素坡最大基质吸力42 kPa划分根系加固边坡的蒸腾作用影响区域。heart根型不仅最大基质吸力最大，且蒸腾作用影响区域最广，plate根型的蒸腾作用稍优于tap根型。

经过长达46.872 d的间隔降雨后，结束渗流场孔隙水压力分布如图5(b)所示(限于篇幅，未给出3种根型每个阶段降雨前后的孔隙水压力分布)。3种根型通过其蒸腾作用和浅层阻水作用都能非常有效地减小和阻止坡体内部水分入渗与运移。与初始状态相比，边坡底部出现一定水位高度，素坡底部坡后水位线约1.73 m，tap根型加固边坡约0.66 m，plate根型加固边坡约0.09 m，heart根型加固边坡未出现水位线。heart根型加固边坡坡体内部水分最少，plate根型加固边坡次之。

图5 4种工况的初始渗流场及结束渗流场孔隙水压力分布Fig.5 Pore water pressures in initial seepage field and end flow field under four working conditions

3.5 4种工况下监测点H9的基质吸力与饱和度

图6为4种工况下监测点H9(近根处)在整个过程中的基质吸力及饱和度变化趋势。采用植被根系加固的3种根型的边坡显著减少了坡体内部的水分，说明植被根系加固边坡的有效性。四阶段降雨结束后，监测点H9处plate根型和heart根型加固边坡水分相等。第四次降雨过程中，plate根型与heart根型两者存在汇合点，其原因是尽管heart根型蒸腾作用效果最好，但是plate根型宽大的两翼更能有效阻挡浅层雨水冲刷与渗透。

图6 监测点H9基质吸力和饱和度随时间的变化趋势Fig.6 Variation trend of matric suction and saturation at monitoring point H9 with time

3.6 4种工况的水平向位移云图对比

在进行应力场及抗剪强度修正后， 可以得到每阶段降雨后及固结后的位移云图。 图7为4种工况结束后的累计水平向位移等值线云图。 据图7分析可得， tap、 plate和heart三种根型加固边坡最大水平向位移百分比分别减少52.22%、 73.68%和70.24%。 结合离心模型试验的滑坡位移现象分析， 水平向位移>1.2 cm范围的区域为潜在滑移区。 3种根型加固的边坡都能有效减小潜在滑移区范围， plate和heart根型表现为浅层冲刷导致的局部破坏。 前者冲刷范围稍大于后者， 但后者由于根型复杂产生应力集中现象， 使得最大水平向位移大于前者。 从整体分析， heart根型固土效果最好， tap根型最差。

图7 4种工况的累计水平向位移云图Fig.7 Cumulative horizontal displacements in four working conditions

3.7 tap根型加固边坡基质吸力和含水率的数据验证

图8为tap根型加固边坡的室内降雨离心模型试验的基质吸力和体积含水率测量值与数值模拟监测值对比。基质吸力与体积含水率数值模拟与试验测试结果趋势吻合度良好，但存在部分差异。考虑由初始状态引起，进行离心模型试验时期为夏季，配好的边坡土体在静置达到初始平衡状态期间存在水分蒸发，致使初始饱和度偏低，基质吸力偏大。本数值模拟不仅对分析边坡降雨入渗土体从非饱和到饱和的变化趋势有一定参考价值，还能有效模拟生态植被根系的蒸腾作用和固土效果，为分析更为复杂的根系奠定了基础。

图8 tap根型加固边坡基质吸力和体积含水率随时间的变化趋势Fig.8 Variation trend of matric suction and volumetric moisture content of tap root slope with time

4 结 论

(1)heart根型加固边坡的水分运移路径与素坡截然不同。素坡表现为近坡表侧水分渗透速率快于远坡表侧，heart根型加固边坡与之相反。heart根型加固边坡近根处出现干燥区且不消散。

(2)初始渗流场显示heart根型蒸腾作用最佳。经过长达46.872 d的间隔降雨后，heart根型加固边坡坡体内部水分最少，plate根型加固边坡次之。由监测点H9的基质吸力和饱和度可知，heart根型蒸腾作用最好，但是plate根型加固浅层土体最能阻止雨水冲刷，有效阻挡坡面水分渗透作用。

(3)从最大水平向位移等值线云图分析，3种根型加固边坡的浅层固土作用都是非常显著的。plate根型最能减小水平向最大位移，heart根型减小量比plate根型少3.44%，heart根型减小量比tap根型大18.02%。heart根型加固边坡存在应力集中现象。

(4)降雨离心模型室内试验和数值模拟的基质吸力与体积含水率变化趋势拟合程度较好，说明本程序能够很好地解决降雨入渗及固结过程中土体从非饱和到饱和的动态变化趋势的问题。但结果存在部分差异，原因可能是夏季水分强蒸发作用使得坡体内部水分较少，初始状态存在一些差异，下一步可考虑结合大气蒸发作用进行进一步的完善与分析。