刘晶, 左双英
(1.贵州大学资源与环境工程学院, 贵阳 550025; 2.贵州大学喀斯特地质资源与环境教育部重点实验室, 贵阳 550025)
中国国土辽阔,地理条件复杂,易发生滑坡灾害,滑坡规模巨大、损坏性强、形成机制和诱发因素众多,且灾后影响区域广,治理难度高。影响涉水岸坡稳定性的因素分为内因和外因,外界因素中以水位升降[1]和降雨[2]为主,两者的共同作用改变了坡体的渗流场、应力场和位移场,易使坡体发生明显变形,从而诱发滑坡灾害。
中外诸多学者对水位升降引起的库岸边坡的稳定性进行了深入研究,肖诗荣等[3]、Tang等[4]针对三峡库区数千余处涉水滑坡,将库区复活型滑坡分为浮托减重型、动水压力型和浸泡软化型。浮托型滑体大多透水性弱,一般发生在库水骤升时期,库水对阻滑段的浮托减抗作用是坡体变形破坏的主要因素。汪发武等[5]以千将坪滑坡为例,分析了顺层结构边坡在水位上升过程失稳并产生高速滑动的机理。肖捷夫等[6]以具有典型多期次滑动特征的藕塘滑坡为对象,制作了大型物理模型试验,模拟库水涨落和降雨工况,通过试验获得了坡体位移、孔隙水压力和土压力全时程曲线,得出该滑坡发生多期次变形的原因。周永健等[7]从定性与定量的角度分析了木鱼包滑坡位移与库水位高程之间的相关性,发现变形对库水位变动的响应较为明显,并且在时间上存在约40 d的滞后期。动水压力型滑坡坡体一般为多孔隙介质,透水性微弱,往往发生在库水骤降时,孔隙水压力的消散滞后于库水消落,形成局部动水压力产生滑坡,如新滩滑坡、树坪滑坡等。殷跃平等[8]对巫山塔坪滑坡进行研究,认为库水抬升导致的孔隙水压力以及库水位下降时滑坡体内的渗透力是导致水库滑坡失稳的主要因素,并分析了多年平均库水位升降工况下堆积层滑坡的浸润线和安全系数。谭淋耘等[9]、朱权威等[10]研究表明,库水位快速下降使得滑坡地下水位调整滞后,指向坡外的渗透压力增大导致滑坡稳定性变差。江强强等[11]通过不同组合条件下的物理模拟,研究库岸滑坡影响因素、变形演化规律及失稳条件。库水浸泡软化型滑坡,一般滑带含有亲水性黏土矿物,浸水易软化,力学强度大幅降低,导致滑坡失稳,最典型的是意大利瓦伊昂水库滑坡[12]。
库岸边坡除了受水位升降影响外,降雨也是导致滑坡复活的重要因子[13],因此探讨坡体在两者共同作用下的动态响应尤为重要。Gu等[14]、Yao等[15]、Shang等[16]通过对库水位、降雨和滑坡变形过程三者之间的相关性分析,研究了其主要成因机制,建立了滑坡变形与外界影响因素之间的对应关系。Jia等[17]以树坪滑坡为例,通过分解滑坡位移,选取粒子群算法优化的最小二乘支持向量机模型对基于降雨量和库水位的滑坡位移预测进行预测。近年来,数值模拟[18-23]也逐渐运用到对滑坡的研究中,以分析库水位与降雨联合作用下坡体的运动特征及变形机理。
综上所述,目前对于涉水滑坡的研究主要集中在基于多年调节水库(如三峡库区)大型滑坡的变形破坏特征和稳定性变化方面,降雨、水位涨落、流量几乎呈周期性变化,坡体变形规律性强。但针对山区高弯陡坡河岸涉水滑坡的研究却不多见,此类河流的特点是降雨引起的水位涨落快、周期短、水位变幅大,与一般人工调节的大库区别较大。而且以上研究虽然运用现场监测、理论分析及数值模拟等各手段分析降雨和库水对滑坡稳定性的影响,但是几乎都是在固定水位变速和雨量条件下进行研究分析,考虑的组合条件过于简单,难以重现实际情况。鉴于此,以贵州平塘六硐南岸滑坡为例,通过工程地质勘查获取了滑坡体的工程地质、环境地质条件,并通过监测手段定量确定坡体变形基本特征,应用有限元Geo-studio软件模拟6—10 月水位和雨量持续波动下的坡体变形,与现场监测进行对比分析,同时探究多种组合(不同水位升降速率、不同降雨强度及同一滞后时间)条件下的滑坡稳定情况,并进行渗流场-应力场耦合分析,探究坡体在不同影响因素下的变形迹象和演化规律,为河岸斜坡在水位变化和降雨影响下的稳定性分析提供一定依据。
滑坡位于贵州省平塘县,属亚热带湿润季风气候区,溶蚀性峰丛槽谷地貌,地形较复杂,地表水系发育。据气象资料统计,年平均气温17 ℃,年内降雨分布不均匀,集中在5—10月,全县降雨量最多年为1 564.4 mm,最少年为998.4 mm。滑坡前缘的六硐河汛期水位为685 m,浸泡滑坡前缘,当恢复到常水位673 m时,暴露整个滑坡体前缘。
滑坡区地处六硐背斜的东翼,构造形迹主要为喜马拉雅山运动发生的褶皱、断裂,呈单斜构造。岩体节理裂隙发育,裂隙发育部位易风化,抗剪强度低。边坡地形坡度约30°,坡体较陡,滑坡体位于河谷南侧山体。滑坡平面形态呈半圆形,后缘、左侧及右侧均以陡崖为界,前缘及中部较宽、后缘窄,坡脚为陡坎。滑坡全貌如图1所示。
HT01~HT04为滑塌01~滑塌04;ZG1~ZG4为GPS监测点
滑体最大纵长约160 m,中前缘宽约500 m,中后部宽约200 m,平均宽约350 m,滑体体积约为5.48×105m3,主滑方向为31°。滑体最小厚度6.2 m,最大厚度13.1 m,主要由第四系冲洪积粉质黏土组成,该层结构复杂、松散,岩土体孔隙度大和可压缩性大,粉质黏土呈褐黄、褐黑色,磨圆好,厚度0.3~13.1 m,分布于滑坡前缘至六硐河河床底部。滑床岩性为二叠系中统茅口组P2m强风化~中风化石灰岩,岩体结构致密,细晶结构,地层产状140°∠13°,滑坡剖面图如图2所示。
图2 工程地质剖面图Fig.2 Engineering geological profile
2012年7月16日在特大暴雨影响下发生滑塌,导致六硐河此段河道变得狭窄,致使县城1/4被淹,经勘察设计论证后对滑坡采取“抗滑桩+截排水沟+河道清理”的方式进行了治理。但在2020年雨季,经特大暴雨入渗冲刷、坡脚浸泡侵蚀作用促使滑坡复活,滑坡前缘出现垮塌及大量裂缝,两侧则沿着岩土界线向下滑动,先后多次发生浆砌石护坡垮塌,抗滑桩呈现不同程度位移,明显垮塌的有4处(图3)。滑塌01(HT01):长16.9 m,宽9.0 m,坡脚掏蚀严重,坡面裂缝长17 m,宽10 cm,可见深度20 cm;滑塌02(HT02):长45.5 m,宽32.6 m,坡脚掏蚀严重,坡面已全部垮塌,部分堆积于河道内,后侧桩体出露;滑塌03(HT03):长90.2 m,宽12.3 m,坡脚掏蚀严重;滑塌 04(HT04):长109.2 m,宽31.5 m,坡脚掏蚀严重,坡面下部已全部垮塌。右侧坡体后缘存在多条裂缝走向110°,长度100 m,后侧裂缝较宽,为10~20 cm,最大可见深度1 m。
图3 滑坡变形破坏迹象Fig.3 Signs of landside deformation and damage
岸坡上共布设4个GPS监测点ZG1、ZG2、ZG3、ZG4(图1),监测时间2020年6月1日—10月31日,定期收集数据。绘制滑坡变形监测数据、6—10月日降雨分布及河流水位变化曲线,如图4所示。
图4 滑坡累积位移时间序列监测曲线Fig.4 Time series monitoring curve of landslide cumulative displacement
由图4可知,坡体上监测点累积位移呈现“阶跃状”特征。2020年6—7月中旬,坡体无明显变形,在经历多次降雨后,监测点ZG1和ZG3变形剧烈,水平和垂直累积位移量均出现两次较明显的上升,于7月14—24日出现第一次小幅度增长,后增长较为缓慢,第二次阶跃上升出现在9月8—20日,该时间段内,河流水位降落加上暴雨作用,使坡体位移量增长变化明显,两监测点水平累积位移分别为123.9 cm和67.3 cm,垂直累积位移达到151.8 cm和52.8 cm。其余两个监测点位移变化幅度平缓,监测点ZG2和监测点ZG4水平累积位移13.1 cm和4.1 cm,垂直方向基本无明显变形。
结合日降雨量分布与河流水位波动数据看,山区河流水位一般在强降雨之后1~2 d上升,保持一定时间,如果没有降雨或其他补给来源水位开始缓慢下降。雨强越大或持续天数越多,水位升得相对更快、更高;连续几天不降雨或降雨强度不大,水位消落迅速。纵观6—10月,此段河流水位在673~685 m波动,水位升降与坡体变形的关系较明显,即水位有变动的时候,变形也随之增加,尤其是水位变动大随之强降雨,会产生较大的变形增量。
由此可见,该类滑坡主要是岩土体内部结构及外在不利因素共同作用引起。在雨季由于上部岩土堆积体松散,降水易于入渗到坡体内部,而坡底又遭遇水位升降动态影响,地下水位线持续变动,饱和区、非饱和区变换频繁,孔隙水压力不稳定,内部岩土体形成渗流场,促使坡体稳定性朝着不利的方向发展。因此,只有深入分析该滑坡的复活机理,才能采取更加科学、有效的措施进行治理。
3.1.1 饱和-非饱和渗流理论
采用seep/w模块进行滑坡渗流场的模拟研究,水位升降和降雨的耦合作用使土体在饱和与非饱和之间转化,地下水位也随之发生规律性的变化。将达西定律和连续方程相结合,得到控制方程为
(1)
式(1)中:H为总水头;kx、ky分别为x、y方向的渗透系数;Q为施加的边界流量;mw为水容比;ρw为水的密度;t为时间。
3.1.2 稳定性分析理论
以极限平衡法为基本原理,采用Morgenstern-Price法,将seep/w模块的渗流计算结果耦合到稳定性计算slope/w模块,在考虑孔隙水压力的情况下进行稳定性模拟分析。
3.2.1 几何模型
利用CAD建立剖面计算模型,高100 m,长156 m,模型左上顶部以滑坡后缘为界,右下底部以河床右缘为界,将建立完成的剖面模型输入有限元软件中。计算区域包括滑体和滑床,上层为粉质黏土,下层为灰岩。利用软件自带的三角形和四边形的方式,对模型进行网格划分,共划分为2 481个网格,2 591个节点。同时,为监测滑坡土体渗流场、 位移场、应力场的变化,在滑面上、中、下部共设3个监测点,如图5所示。
图5 滑坡有限元计算模型Fig.5 Finite element calculation model of landslide
3.2.2 计算参数
(1)计算模型和参数。滑带土抗剪强度的取值对滑坡稳定性分析结果影响较大,釆用传递系数法建立折线形滑面力学结构模型,计算时把整个断面上的滑坡体适当划分成6 个条块(①~⑥),如图6所示。
图6 滑坡断面计算模型简图Fig.6 Calculation model diagram of landslide section
(2)敏感性分析。基于室内试验,分析获得的抗剪强度参数分布区间,滑带土黏聚力c为31.8~37.3 kPa,统计标准值为35.85 kPa,内摩擦角φ为9.7°~12.6°,统计标准值为11.83°,各取差值1作为计算步距进行计算[24],同时绘制黏聚力c和内摩擦角φ与稳定性系数的关系曲线,如图7所示。
图7 抗剪强度参数与稳定性系数关系曲线Fig.7 Relationship curve between shear strength parameters and stability coefficient
由图7可知,稳定性系数Fs与c、φ呈线性相关,从曲线斜率(图7)即可判断,在相同系数区间内,φ对Fs的影响程度明显大于c,表明φ对Fs的敏感性较强,滑坡稳定性与滑动面的内摩擦角关系更为密切。
(3)反演结果。c的敏感度较低,故对其可通过试验来确定其数值,取c=35.85 kPa,而对于敏感性较高的参数φ,在滑坡的极限平衡状态及边界条件下,进行滑带抗剪强度参数反演[25]。计算得到φ=10.84°,所得强度参数介于试验值范围附近,说明本次反演结果较可靠。通过对滑坡土工试验数据统计、工程类比法以及反演计算综合分析,采用的岩土物理力学参数如表1所示。
表1 岩土体物理力学参数
3.2.3 计算工况
为探究不同水位升降速率及降雨强度下,滑坡的渗流场和稳定性变化规律,结合滑坡的工程地质特征和场区河流水位、降雨量变化情况,设计以下计算工况和荷载组合,如表2所示。
表2 有限元模拟计算工况
3.3.1 不同降雨强度对坡体稳定性的影响
参照平塘县多年降雨及6—10月日降雨资料,设置40、75、150 mm/d 3种降雨强度来模拟入渗条件,持续降雨历时3 d、滞后5 d,即工况1、2、3,分析同一降雨历时和滞后时间下、不同降雨强度对滑坡渗流场及稳定性的影响,计算结果如图8所示。
图8 降雨作用下坡体稳定系数变化Fig.8 Variation of slope stability coefficient under rainfall
由图8可知,坡体稳定性随降雨持续入渗逐渐降低,降低速率随降雨强度的增大而加快。岸坡初始状态时,稳定系数为1.050,属于基本稳定状态。降雨强度为150 mm/d时,坡体稳定系数降落最快,达到最低值1.011,随后岸坡稳定性逐渐恢复,该雨强下坡体稳定性回升速率最慢,最终仅升至1.015;降雨强度为75 mm/d时,稳定系数降至最低值1.030,最终在滞后5 d时上升至1.036;降雨强度为40 mm/d时,稳定系数变化平缓,降低至1.043,后缓慢上升至1.048。
根据瞬态渗流场的计算结果,得到不同降雨强度下,滑坡监测点1、2、3的孔隙水压力变化规律及坡体内浸润线的分布,如图9和图10所示。
图9 降雨作用下最大孔隙水压力变化Fig.9 Variation of maximum pore water pressure under rainfall
图10 降雨作用下坡体浸润线变化Fig.10 Variation of slope infiltration line under rainfall
由图9可知,3个监测点的最大孔隙水压力变化趋势遵循一定规律,降雨持续作用时变化明显,滞后时间段内变化平缓。该过程中,土体基质吸力逐渐减小,土体饱和区随降雨的持续入渗向上移动,土体饱和度增大,饱和区扩大,非饱和区逐渐转变为饱和区。监测点1位于水位之上,呈现负孔隙水压力,土体孔压变化速率随雨强的增大而增大。监测点2和监测点3位于水位之下,监测点2距离地下水位线较近,孔隙水压较小,监测点3布设于基岩内,深度较大,故孔隙水压值较大。
由图10可知,各雨强下的地下水位线渗入点一致,整体均呈现出略上凹的趋势,当降雨强度较大时,随着降雨持续入渗,地下水位不断上升,雨强越大,地下水位抬升越高。
3.3.2 河流水位升降速率对坡体稳定性的影响
分别以0.5、1.0、1.5 m/d的河流水位升降速率模拟计算滑坡渗流场的变化情况,即水位上升工况4~6和水位下降工况7~9,水位在673~685 m变化,得到滑坡的稳定系数和浸润线变化图,分别如图11和图12所示。
图11 河水位升降下坡体稳定系数变化Fig.11 Variation of slope stability coefficient under water level fluctuation
图12 河水位升降下坡体浸润线变化Fig.12 Variation of slope infiltration line under water level fluctuation
由图11可知,当水位以不同的速率升降时,坡体稳定系数的变化趋势基本一致,变化速率随水位变化速率的增大而增大。岸坡在水位763 m时稳定系数为1.048,工况4、工况5和工况6下水位每上升1 m,稳定系数平均上升0.021 4、0.021 9、0.025 1,最终分别升至1.305、1.312和1.349。水位上升时,坡体的弱透水性使河水入渗坡体较慢,坡体地下水位的上升速率小于河水位,河水补给地下水,形成指向坡内的动水压力,反作用于坡体,对坡体的稳定性有利,滑坡稳定性增强。当河水位为685 m时,岸坡稳定性系数为1.192,工况7、工况8和工况9下稳定系数分别降至1.119、1.112和1.097。河水位下降时,地下水位降落较慢,动水压力指向坡外,作用于滑坡前缘,对滑坡稳定性不利。
由图12可知,坡体内地下水位随水位上升而上升,呈右弯下凹的趋势。因坡体渗透系数较小,河水位和地下水位之间形成指向坡体内的动水压力,河水位上升速率越大,地下水位线所处位置越低、变化趋势越明显,与河水位形成的水头差越高,产生的动水压力越大,岸坡稳定性增加越明显;河水位下降,地下水位线的变化趋势与水位上升时相反,其所处位置与河水位下降速率有关,降落速率越大,地下水位线所处位置越高。水位下降时,岸坡内地下水缓慢渗出,滞后于水位的降落,动水压力指向临空面,坡体下滑力增大,岸坡的稳定性降低。
3.3.3 河水位下降联合降雨对坡体稳定性的影响
(1)稳定系数。分析可知,水位上升会使坡体稳定性上升,降雨、水位下降分别作用会使坡体稳定性降低,故最不利组合因素为降雨与水位下降联合作用,以0.5、1.0、1.5 m/d的水位下降速率及75 mm/d雨强入渗3 d的条件进行模拟,即工况10~工况12,与仅在水位下降作用下(工况7~工况9)的结果进行对比分析,如图13所示。
图13 坡体稳定系数变化规律比较Fig.13 Comparison of slope stability coefficient
工况10~工况12下,降雨叠加河水位降落作用时,滑坡稳定系数依次从1.192降低到1.062、1.056、1.034;工况7~工况9下,仅在降雨作用下,稳定系数分别从1.192降落到1.119、1.112和1.096,由此可知,滑坡在两者的共同作用下,稳定系数降落速率和幅度更大,且水位下降速率越大,滑坡稳定性下降得越快。降雨入渗会使滑坡土体由非饱和变为饱和,坡体下滑力增大,并降低岩土体抗剪强度,使滑坡体的抗滑力减小,同时产生指向坡外的动水压力,两者的叠加作用使坡体更易失稳。
(2)位移场和应力场。河水位下降及降雨的作用不仅影响滑坡体的渗流场,而且使得滑坡体的应力场、位移场发生相应的变化。
工况10和工况7计算得到的坡体水平、竖直方向的最大拉应力变化规律比较如图14(a)所示。两种工况下,水平和竖直方向最大拉应力随着水位的下降而增大,且工况10增大得更快,水平方向由0升至424.49 kPa,竖直方向由0升至427.62 kPa。因土体承受拉力的能力较弱,故坡体表面易发生拉裂破坏,产生裂隙。
图14 坡体应力场和位移场变化规律比较Fig.14 Comparison of slope stress field and displacement field variation rules
从图14(b)可以看出,随着水位的下降,工况7和工况10坡体最大水平位移分别达到8.9 cm和22.6 cm,最大竖直位移分别达到为12.8 cm和29.5 cm,降雨和水位下降的耦合作用(工况10)使坡体变形增大较快,水平、垂直位移增长速率分别达到0.79 cm/d和0.99 cm/d。
3.3.4 模拟结果与实际情况分析对比
通过对降雨和河流水位变化过程对滑坡变形的影响分析,表明滑坡变形与水位升降速率和降雨强度等因素有密切关系。因此,很有必要讨论6—10月实际气候变化对滑坡变形的具体影响,并将模拟结果与实际情况进行对比分析。以收集到的降雨量和河流水位数据为水力条件,进行为期153 d的模拟监测,期间最大日降雨量为106.8 mm/d,河流水位经历两次较大幅度的“下降-上升”,坡体总位移变化分布如图15所示。
图15 坡体位移分布云图Fig.15 Distribution cloud diagram of slope displacement
据现场监测数据,监测点ZG1累积位移变化最大,水平累积位移为123.9 cm,垂直累积位移为151.8 cm。有限元模拟结果中监测点ZG1总位移为115.7 cm,略小于现场位移变化,其原因可归结于:岸坡在河水的长期浸泡和冲刷作用下,岩土体的强度降低,坡体整体稳定性变差;而数值模拟难以还原坡体数年以来所受的水-岩作用,仅模拟153天岸坡在降雨和河水位波动耦合作用下的变形趋势。
由图15可知,滑坡发生时,整体向前移动,坡体挤压河道,使河道鼓起,最大变形量发生在坡体前缘,从前缘往坡内位移逐渐变小,基岩几乎不产生变形量。综上,当滑坡变形产生时,其前缘部分在降雨入渗、河流冲刷和水位升降的综合作用下,破坏最先发生,变形量最大,整体表现为前缘发生拉剪破坏,呈典型的牵引式滑坡变形特征。
综合以上分析,确定该滑坡形成条件主要为内因和外因,内在因素有地形地貌、岩土体结构、地质构造等,外在因素有连续降雨、河水位降落、人为工程活动等,此处主要讨论外部影响机制。
据资料查阅,可知全县多年平均降水量为1 259 mm,雨量充沛;且在坡体产生两次大变形之前的1~2 d有持续强降雨。在强降雨作用下,岩土体重度变大,下滑力增加,抗剪强度降低,抗滑力减小,同时雨水沿岩土分界面渗流,增加了动水压力。
山区河流水位升降受降雨和地表径流影响大。坡体发生两次大变形分别在7月和9月,而8月降雨很少,河流波动小,给水位下降与强降雨叠加作用提供了条件。
(1)动水压力作用。水位上升时,河水缓慢渗入坡体内部,河水位与坡内地下水位形成水头差,动水压力指向坡内,并产生垂直于滑面表面的水压力,反压滑坡体阻滑段,有利于滑坡的稳定;水位下降时,滑体渗透性较差,坡内地下水无法及时排出,具有一定滞后性,坡内地下水和河水位之间形成指向坡外的动水压力,降低滑坡的稳定性。
(2)冲刷作用。水位升降产生的纵向冲刷,带走前缘土体中的细小颗粒,发生局部渗透破坏,产生淤积体;产生的横向冲刷,使岸坡土颗粒被搬运,使河岸线发生后移,且冲刷作用强弱与水流速成正比。在横纵冲刷的共同作用下,滑坡前缘第四系冲洪积粉质黏土长期受到冲刷淘蚀,岸坡发生局部塌方,使得坡脚临空面高度逐渐增加,加剧了滑坡的形成。
坡体上属于耕地,为村集体所有,植被覆盖率较差,使得坡体表面土层结构更为松散,利于雨水渗入,增加土坡重量,减小抗剪强度,不利于滑坡的稳定性。
结合滑坡工程地质条件,分析滑坡的变形特征和形成诱因,通过数值模拟分析滑坡在不同工况下的稳定性变化规律,得出如下主要结论。
(1)该滑坡累计位移曲线呈现“阶跃状”特征,二次明显变形与河水位升降、强降雨在时间上存在对应关系。
(2)该滑坡灾害发生的内因为不良的地形地貌及岩土体组成,外因为降雨和河水位升降的共同作用。以饱和-非饱和渗流理论及极限平衡原理为依据,利用Geo-studio软件模拟不同工况下坡体的响应及其稳定性变化趋势。
(3)水位下降,坡体内地下水的排出滞后于河水位的降落,动水压力指向坡外,岸坡稳定性降低,稳定系数变化速率随水位下降速率的增大而增大;河水位上升,滑坡岩土体所受压力增大,且动水压力指向坡内,滑坡稳定性增加,水位上升速率越大,岸坡稳定性增加越明显。叠加降雨作用后,雨水渗入坡体,产生不利作用,坡体稳定性系数进一步降低。
(4)通过数值模拟分析岸坡在不同工况下,地下水位线、稳定系数、应力场和位移场的变化规律,得到滑坡对外界不同影响因素的响应规律。可知坡体受降雨和水位下降的共同作用时,渗流场和应力场发生极大变化,稳定性最差。
(5)据153 d的二维变形破坏分析可知,坡体ZG1监测点总位移变化为115.7 cm,与现场监测数据相近,滑坡后缘及中部未发现变形,坡体变形破坏主要发生在坡体前缘,该滑坡属于典型的牵引式滑坡。2020年7—9月,正是因为坡体历经2次“河水位下降+强降雨”工况,导致滑坡复活并且变形加剧,结论可为山区河岸滑坡防灾减灾提供一定参考。