基于改进Green-Ampt降雨入渗模型的花岗岩风化带边坡稳定性分析

蒋天娇，简文星，姚远，田朋飞，张树坡

(中国地质大学(武汉) 工程学院，湖北 武汉 430074)

花岗岩是我国分布最为广泛的岩体之一，约占国内山地总面积的15%，在广东、福建、江西、桂东南和湘南地区一带分布较为集中，其中闽粤地区的花岗岩出露面积可占其总面积的30%～40%[1]。多数花岗岩节理发育，存在良好的物理化学风化入侵通道，受南方地区多雨潮湿气候的影响，风化剥蚀强烈，往往在山谷地段形成高达数十米甚至百余米的棕黄色-红褐色花岗岩风化带[2]。天然的花岗岩风化带边坡大多分带特征明显、结构松散，且含有较多的原生、次生裂隙以及残余结构面[3]，获取此类边坡的分带模式、力学强度、渗流特征和失稳破坏模式，对于花岗岩分布区边坡灾害防治和工程安全保障具有重要意义。

相关研究表明，花岗岩风化带岩土体边坡在切坡卸荷、持续降雨等因素影响下，极易发生失稳破坏，主要包括坡脚渗流冲刷导致的牵引式滑移、坡顶孔隙水压力升高引起的近垂直滑塌、坡面充水强度衰减形成的浅层滑动等，有时甚至发生暴雨溃决而成为大型泥石流的物质来源，给所在地区的人民生命和财产以及工程安全造成严重的危害[4-7]。降雨入渗对花岗岩风化带岩土体边坡的影响涉及边坡各风化带非饱和-饱和的浸润过程，主要表现为岩土体含水率升高、容重增大以及强度衰减，继而导致边坡破坏[8]。结合降雨入渗模型，获取不同降雨特征下花岗岩风化带岩土体边坡的降雨浸润规律并进行边坡稳定性评价，是分析降雨诱发型边坡失稳破坏规律的有效方法[9-10]。

基于毛细理论的Green-Ampt降雨入渗模型(以下简称GA模型)，模型参数少、物理意义明确且适用类型广，其以湿润锋的推进来直观地求解降雨诱发型边坡的失稳破坏模式[11-12]，多年来一直被国内外学者不断发展、改进，并广泛用于边坡降雨入渗问题的研究。如Cho[13]通过改进GA模型中的浅层不透水边界条件，探讨了降雨强度变动下残积土边坡潜在滑移面的变化规律及其稳定性；Dorofki等[14]、Liu等[15]考虑坡度和初始含水率因素，提出了适用于不同坡度坡面、非均质土层的GA模型;马娟娟等[16]基于GA模型，引入湿润区的平均含水率概念，建立了3种变水头积水降雨入渗模型；彭振阳等[17]通过引入湿润层的等效导水率并给出锋面吸力的计算方法，对传统GA模型进行了改进，改进后模型的计算精度有显著提高；唐扬等[18]在改进GA模型的基础上，以非饱和土VG模型研究了初始含水率对坡体浸润稳定性的影响。上述研究结果表明，结合实际边坡特征适当改进GA模型，对于研究花岗岩风化带岩土体边坡降雨诱发失稳的演变过程和致灾机理等问题至关重要。

赣南山区位于江西省南部地区，与湖南、广东、福建等地相接壤，属亚热带季风气候区，年平均降雨量高达1 605 mm，降雨时间集中在每年的3～8月份。研究区内岩浆活动频繁，其中花岗岩出露区面积约占全区面积的1/3，且多形成巨厚层状的花岗岩风化带。区内公路、铁路、市政等工程建设中，常常形成大量挖方边坡，坡体在连续降雨入渗的作用下极易发生失稳破坏。但与沿海地区花岗岩风化带岩土体边坡已有的丰富的科学研究和大量的工程防治实践相比，针对该区域花岗岩风化带岩土体边坡在连续降雨条件下稳定性的研究尚较少。为了研究连续降雨对赣南花岗岩风化带岩土体边坡稳定性的影响，本文依托江西安远至定南段高速公路项目，通过建立花岗岩风化带路堑边坡地质模型，采用资料收集、室内试验等手段获取边坡岩土体物理力学参数，并基于改进的Green-Ampt降雨入渗模型，采用极限平衡法和强度折减法数值模拟方法对不同降雨工况下花岗岩风化带岩土体边坡的降雨浸润规律、稳定性和变形破坏模式进行分析。该研究结果可为赣南及我国南方其他花岗岩分布区风化带岩土体边坡降雨入渗稳定性分析和边坡安全保障提供重要参考。

1 花岗岩风化带路堑边坡地质模型建立

本研究所选取的边坡为赣南安远至定南新建高速公路定南段的K196+380左侧花岗岩风化带路堑边坡，在分析时将该边坡的坡面进行简化处理，未考虑后续边坡支护情况，其工程地质剖面图见图1。

图1 赣南地区花岗岩风化带路堑边坡工程地质剖面图Fig.1 Engineering geological profile of the cutting slope of granite weathering zone in Southern Jiangxi Province of China

该路堑边坡在研究时已开挖完成，边坡上部全风化覆盖层较厚，岩体结构面出露不明显。据勘察报告钻孔资料以及野外地质调查情况，该边坡岩土体的物质组成从上往下可分为3层：①全风化花岗岩层，黄褐色夹麻灰色，原岩结构基本被破坏，可定为黏土质砂,土样中夹少量石英颗粒和强风化碎块，厚度大;②强风化花岗岩层：灰褐色，岩体极破碎，发育有节理裂隙，呈碎块状，厚度因所处的环境不同而有较大的差异;③中风化花岗岩层，灰白色，岩体节理裂隙不发育，岩质较硬，厚度大。

2 边坡岩土体参数获取

本文通过室内基本土工试验、三轴剪切试验、岩体压缩试验等方法，结合线路工程勘察报告，获取了所选取的花岗岩残积土风化带边坡不同风化带岩土体的物理力学参数，具体如下：

(1) 花岗岩全风化带：该层厚度大，具有天然固结的围压条件，开挖后排水措施滞后，取该边坡花岗岩全风化带岩土体样品在室内利用TKA三轴剪切渗流仪进行不同围压下的固结不排水试验，以获取花岗岩全风化带岩土体的抗剪强度参数，并通过施加渗流水头的方法获取花岗岩全风化带岩土体的饱和渗透系数，其水力特征参数见表1。

表1 赣南地区花岗岩全风化带岩土体的水力特性参数

(2) 花岗岩强风化带：由工程地质勘察资料统计分析获取此路段花岗岩强风化带岩土体强度参数的标准值，包括重度、黏聚力、内摩擦角、弹性模量和泊松比等。

(3) 花岗岩中风化带：取该边坡下部花岗岩中风化带岩土体样品在室内RMT-150C岩石力学实验系统上进行单轴压缩试验(见图2)，获取花岗岩中风化带岩土体的单轴抗压强度、弹性模量和泊松比，然后再通过三轴压缩试验获取岩土体黏聚力和内摩擦角。

图2 中风化花岗岩岩土体单轴压缩试验Fig.2 Uniaxial compression test of moderately weathered granite

通过对以上花岗岩残积土风化带边坡各层岩土体的强度参数进行整理与整合，见表2。

表2 赣南地区花岗岩残积土风化带岩土体的强度参数

3 考虑排水的改进Green-Ampt降雨入渗模型

在持续降雨作用下，花岗岩风化带边坡表层土体将处于饱和状态，并且饱和带伴随降雨的连续性逐步往下延伸。为了研究赣南花岗岩风化带岩土体边坡在强降雨作用下的稳定性，本文引入Green-Ampt降雨入渗模型。该模型可以用来描述地表水在重力和基质吸力作用下的入渗规律，见图3和图4。

图3 Green-Ampt降雨入渗模型下地表水的入渗规律Fig.3 Surface water infiltration under Green-Ampt model注：t为选定时刻；Z为地表往下的距离；Zw为选定时刻湿润锋的深度；θs为全风化花岗岩的饱和含水量；θi为全风化花岗岩的初始含水量；ψw为湿润锋处基质吸力水头。

图4 Green-Ampt边坡模型Fig.4 Green-Ampt slope model

Green-Ampt降雨入渗模型的基本理论公式如下：

(1)

式中:iw为降雨入渗率；Ks为全风化花岗岩的饱和渗透系数。

(2)

其中，Zw在边坡中取湿润锋垂直于坡表的深度。此公式适用于湿润锋形成过后的情况，所以应把边坡降雨入渗过程分为湿润锋形成前、后两部分。随着降雨的持续进行，会出现湿润锋形成的临界时刻tp，降雨强度为q，边坡入渗速率iw为qcosα，tp时刻的浸润锋深度为Zp，将iw=qcosα和Zw=Zp代入公式(2)，可得：

(3)

则该时刻的累计入渗量Ip和临界时刻tp为

(4)

(5)

式中:θs为全风化花岗岩饱和状态的体积含水率;θi为全风化花岗岩天然状态的体积含水率。

为了考虑降雨与排水的共同作用，可结合达西定律对GA模型进行改进。边坡湿润锋的下降速率等于降雨入渗率减去饱和带平行于坡表的水分排出率，即：

(6)

式中:L为斜坡沿坡表方向的长度。

初始条件为t=tp、Zw=Zp，湿润锋平均基质吸力水头为ψw，则花岗岩全风化带湿润锋基质吸力水头为8.17 cm[19]。

4 降雨入渗下花岗岩风化带岩土体边坡的稳定性分析

4. 1 降雨工况选取

影响赣南花岗岩风化带岩土体边坡稳定性的主要因素是该地区的持续强降雨天气，定南县近30年的月平均降雨量数据统计结果，见图5。

图5 定南县近30年月平均降雨量数据统计Fig.5 Monthly average rainfall in Dingnan County in recent 30 years

根据中国气象数据网、定南县人民政府网数据，定南县月平均降雨量最大值发生在2007年6月，其值为617.4 mm，单日最大降雨量发生在2010年5月7日，其值为335 mm。通过对定南县降雨资料进行分析，设置了如下两种强降雨工况：

(1) 工况一：选取月平均降雨量为617.4 mm/月，为2007年6月整月降雨工况；

(2) 工况二：选取单日最大降雨量为335 mm/d，为2010年5月7日整日降雨工况。

此外，根据公式(3)和公式(5)，分别计算出两种工况下赣南地区花岗岩风化带岩土体边坡湿润锋形成的临界深度和临界时刻如下：

(1) 工况一：Zp=0.43 m,tp=7.56×104s；

(2) 工况二：Zp=0.006 m,tp=24.5 s。

4. 2 理论计算与数值模拟

本文采用FLAC3D软件自带的强度折减模块进行了该花岗岩风化带岩土体边坡的稳定性计算，计算模型选取摩尔-库仑模型。在不同时间节点下，边坡湿润锋以上选取饱和状态下边坡岩土体的重度和抗剪强度参数，边坡湿润锋以下选取天然状态下边坡岩土体的重度和抗剪强度参数进行数值模拟计算，得到在不同降雨条件下赣南地区花岗岩风化带岩土体边坡的折减系数、破坏方式和潜在滑动面，以强度折减系数作为该边坡的稳定性系数。

4.2.1 工况一条件下边坡的稳定性系数计算

(1) 根据赣南地区花岗岩全风化带岩土体的水力特性参数(见表1)和公式(6)，可得到工况一条件下该赣南地区花岗岩风化带岩土体边坡湿润锋深度随降雨时间的变化曲线，见图6。

图6 工况一条件下赣南地区花岗岩风化带岩土体边坡 湿润锋深度随降雨时间的变化曲线Fig.6 Variation curve of slope wetting front depth of the rock and soil in grantic weathing zone with rainfall time under working condition 1 in Southern Jiangxi Province of China

(2) 根据赣南地区花岗岩风化带岩土体的强度参数(见表2)和工况一条件下该边坡湿润锋深度随降雨时间的变化曲线(见图6)，选取降雨时间分别为1 d、3 d、5 d、10 d、15 d和20 d时的边坡湿润锋深度，将湿润锋处作为潜在的滑动面，滑动面以上边坡岩土体的重度和强度参数选取饱和状态参数，利用刚体极限平衡法计算工况一条件下该边坡的稳定性系数，并绘制边坡稳定性系数随降雨时间的变化曲线，见图7。

图7 工况一条件下极限平衡法计算该边坡稳定性 系数随降雨时间的变化曲线Fig.7 Variation curve of stability coefficient calculated by limit equilibrium method with rainfall time under working condition 1

由图7可见，降雨第10 d时该边坡的稳定性系数为1.08，而降雨第15 d时边坡的稳定性系数为0.97，以边坡稳定性系数小于1作为边坡失稳的评判标准，该赣南花岗岩风化带岩土体边坡在此时工况下会发生失稳破坏。

本文建立的FLAC3D边坡模型包括11 670个单元和13 970个节点，左右边界采取水平方向0位移边界，底部边界采取水平与垂直方向均为0位移边界，拉伸方向即Y轴方向亦采取0位移边界。为了更好地展示边坡的破坏模式，本文所有剪应变增量云图均在水平方向即X轴方向旋转20°进行展示，这样可以方便观察到边坡坡脚以及台阶处的变形破坏情况。

本文选取降雨时间分别为0 d、1 d、3 d、5 d、10 d、15 d和20 d时的边坡湿润锋深度，湿润锋以上分别选取饱和状态下边坡岩土体的重度和抗剪强度参数，湿润锋以下分别选取天然状态下边坡岩土体的重度和抗剪强度参数，利用FLAC3D软件采用强度折减法进行数值模拟计算，得到工况一条件下不同降雨时长时该边坡的剪应变增量云图和边坡稳定性系数随降雨时间的变化曲线，见图8和图9。

通过分析图8和图9可知，赣南地区花岗岩风化

图8 工况一条件下不同降雨时长时该边坡的剪应变增量云图Fig.8 Cloud map of shear strain increment of the slope with different rain fall time under working condition 1

图9 工况一条件下该边坡稳定性系数随降雨时间的 变化曲线Fig.9 Variation curves of slope stability coefficient with rainfall time under working condition 2

带岩土体边坡在工况一的强降雨入渗作用下，在降雨20 d时边坡的稳定性系数小于1，会发生失稳破坏；极限平衡法与强度折减法计算得到的该边坡稳定性系数基本一致，但在降雨初期，由于极限平衡法选取湿润锋处为边坡的潜在滑动面会导致边坡稳定性系数偏高，使其与强度折减法计算所得的边坡稳定性系数有一定的差异，体现为极限平衡法在初期的计算结果相对偏高。

此外，根据图8和图9的模拟结果可见，未降雨时，该花岗岩风化带岩土体边坡的初始稳定性系数为2.72，该边坡应变增量主要在坡脚出现集中；随着降雨的持续进行，该边坡的潜在滑动面在坡体表面逐渐贯通，并在降雨3 d左右时完全贯通，边坡的稳定性系数明显下降；降雨期间该边坡的剪应变增量最大处始终出现在坡脚位置，破坏区域呈现逐步向上的发展趋势；在降雨5 d后该边坡台阶处开始呈现高剪应变增量区，边坡各台阶处很可能会出现剪切破坏；在降雨20 d时，该边坡的强度折减系数低于1，认定为该边坡出现破坏，破坏面在浅层坡表贯通，此时很可能发生坡面冲刷破坏，而高剪应变增量出现在边坡浅层台阶和坡脚处，可认为边坡台阶处、坡脚处以及整个坡表都是容易发生滑移破坏的地方。

综上所述，赣南地区花岗岩风化带岩土体边坡在工况一的降雨作用下，坡表很可能发生冲刷破坏，而坡脚处和台阶处均呈现一定程度的剪切滑移，且坡脚处的破坏程度更大；该边坡的稳定性系数和强度折减系数也随着降雨的持续而下降，直到边坡发生失稳破坏(边坡稳定性系数和强度折减系数均小于1)。

4.2.2 工况二条件下边坡的稳定性系数计算

(1) 根据赣南地区花岗岩全风化带岩土体的水力特性参数(见表1)和公式(6)，可得到工况二条件下赣南地区花岗岩风化带岩土体边坡湿锋峰深度随降雨时间的变化曲线，见图10。

图10 工况二条件下赣南地区花岗岩风化带岩土体边坡 湿润锋深度随降雨时间的变化曲线Fig.10 Variation curve of slope wetting front depth of the rock and soil in grantic weathing zone with rainfall time under working condition 2 in Southern Jiangxi Province of China

(2) 根据赣南地区花岗岩风化带岩土体的强度参数(见表2)和工况二条件下该边坡湿润锋深度随降雨时间的变化曲线(见图10)，选取降雨时间分别为1 h、3 h、6 h、12 h、18 h和24 h时的边坡湿润锋深度，将湿润锋处作为潜在的滑动面，滑动面以上边坡岩土体的重度和强度参数选取饱和状态参数，利用极限平衡法计算工况二条件下该边坡的稳定性系数，并绘制边坡稳定性系数随降雨时间的变化曲线，见图11。

图11 工况二条件下极限平衡法计算该边坡稳定性系 数随降雨时间的变化曲线Fig.11 Variation curve of slope stability coefficient calculated by limit equilibrium method with rainfall time under working condition 2

由图11可见，在降雨12 h时该边坡的稳定性系数达到1.03，在降雨18 h时该边坡的稳定性系数为0.98，即小于1，说明赣南地区花岗岩风化带岩土体边坡在此工况下会发生失稳破坏。

本文选取降雨时间分别为1 h、3 h、6 h、12 h、18 h和24 h时的边坡湿润锋深度，湿润锋以上分别选取饱和状态下边坡岩土体的重度和抗剪强度参数，湿润锋以下分别选取天然状态下边坡岩土体的重度和抗剪强度参数，利用FLAC3D采用强度折减法进行数值模拟计算，得到工况二条件下不同降雨时长时该边坡的剪应变增量云图和边坡稳定性系数随降雨时间的变化曲线，见图12和图13。

图12 工况二条件下不同降雨时长时该边坡的剪应变增量云图Fig.12 Cloud curves of slope shear strain increment with rainfall time under working condition 2

图13 工况二条件下该边坡稳定性系数随降雨时间的 变化曲线Fig.13 Variation curves of slope stability coefficient with rainfall time under wording condition 2

通过分析图12和图13可知，赣南地区花岗岩风化带岩土体边坡在工况二的强降雨入渗作用下，边坡坡脚处剪应变增大较快，边坡的稳定性系数逐渐降低并在降雨18 h之后边坡稳定性系数小于1，可认为边坡即将发生失稳破坏；与工况一相似的是，极限平衡法与强度折减法计算得到的该边坡稳定性系数基本一致，但在降雨初期，由于极限平衡法选取湿润锋处为边坡潜在的滑动面会导致边坡稳定性系数偏高，使其与强度折减法所得的边坡稳定性系数有一定的差异，体现为极限平衡法在初期的计算结果偏高。

除参数外，FLAC3D计算前后的处理方式如工况一，由图12和图13的模拟结果可见，未降雨时，该边坡的初始强度折减系数也为2.72，随着降雨的进行，该边坡的潜在滑动面从坡脚应变集中逐渐由坡脚向上扩展；在降雨达到3 h左右时边坡的潜在滑动面完全贯通，边坡强度折减系数也快速下降，且边坡剪应变增量最大处出现在坡脚位置，破坏区域也随着降雨时长的增加而增大，边坡各台阶处也均出现了潜在的剪切破坏区；在降雨达到24 h时，边坡的稳定性系数仅为0.92，可认为此时边坡已发生失稳破坏。

相对于工况一，工况二的降雨速率更大，因此该边坡的湿润锋下移速率加快，边坡相应的浅层破坏区域快速成型贯通，且剪应变增量较大，边坡有很大的整体失稳风险。整体而言，在工况二条件下，赣南地区花岗岩风化带岩土体边坡的破坏模式为从坡脚往上延伸的浅层剪切滑移破坏以及各个台阶处的局部剪切破坏，各破坏区逐渐贯通后，最终引起该边坡的整体失稳破坏。

5 结论与建议

本文以安远至定南高速公路定南段的K196+380左侧花岗岩风化带路堑边坡为例，采用现场工程地质调查、室内试验等方法，对研究区的花岗岩风化带岩土体进行了风化分带，整理得到赣南地区花岗岩风化带岩土体强度特征，再通过边坡地质模型的建立与降雨资料的搜集整理，采用考虑排水的改进Green-Ampt降雨入渗模型，分别采用极限平衡法和强度折减法数值模拟方法，计算得到在强降雨作用下该边坡的稳定性系数和变形破坏模式随降雨时间的变化情况，得到如下主要结论：

(1) 已经开挖的K196+380花岗岩风化带岩土体边坡可分为三层，从地表往下依次为全风化花岗岩层、强风化花岗岩层和中风化花岗岩层。通过室内TRA三轴剪切渗流试验和岩体压缩试验，结合工程勘察资料，得到了赣南地区花岗岩各风化带岩土体的力学强度和水力特征参数。

(2) 研究区定南县暴雨频繁，近30年的月平均降雨量最大值为617.4 mm(工况一)，单日最大降雨量为335 mm(工况二)。引入Green-Ampt降雨入渗模型，并考虑降雨与排水的共同作用，参考达西定律对该模型进行了改进，使边坡湿润锋的下降速率等于降雨入渗率减去饱和带平行于坡面的水分排出率。

(3) 理论计算与数值模拟结果表明：随着降雨的持续，该边坡的稳定性系数逐步降低，体现在工况一条件下降雨20 d后边坡的稳定性系数小于1，工况二条件下降雨18 h后边坡的稳定性系数小于1；该边坡的坡体破坏从坡脚逐步向上延伸，并在坡体表面形成潜在的滑动面，而后台阶处也逐渐破坏；整体而言，工况一条件下该边坡的破坏模式以坡面冲刷和浅层剪切滑移破坏为主；工况二条件下，该边坡的高应力区集中在坡脚处和台阶处，边坡的破坏模式由浅层剪切滑移逐步转为整体失稳。

基于上述研究结果，建议雨季在该区域进行花岗岩风化带岩土体边坡开挖施工时，须考虑极端降雨条件的影响，并做好坡面冲刷和坡体滑移的防护工作。