库岸堆积体斜坡降雨过程稳定性研究
吴海真
河海大学水利水电工程学院
江西省水利科学研究院
在对降雨入渗及其引发边坡失稳机制的研究基础上，基于饱和-非饱和渗流数学模型，运用延伸的摩尔-库伦强度理论，对某库岸堆积体斜坡进行了降雨过程的饱和-非饱和渗流场计算和稳定分析。结果表明降雨是导致库岸堆积体滑坡的主要外动力因素之一，并得出了在采用实际可能降雨特征时，库岸堆积体的最小安全系数发生在雨后一定时间内的结论。这与实际工程情况相吻合，可更合理的解释因降雨引发的库岸堆积体滑坡问题。
1 概述
库岸堆积体是一类由滑坡堆积、崩塌堆积、残积层、冰溃堆积、坡积物等组成的松散体，即主要受重力作用堆积而形成的第四系覆盖层斜坡，在我国尤其是三峡库区和西南地区分布广泛。其物质组成为典型的土石混合体，是介于岩体和土之间的一类特殊地质体，使得人们较难把握其稳定性及其变形发展趋势，故在国内水电建设工程中发生滑坡并酿成重大灾害的情况时有发生。研究表明，降雨是影响库岸堆积体变形发展的主要外动力因素［1-2］。因此，研究因降雨入渗引发库岸堆积体滑坡的机制和规律具有重要意义［3］。在对降雨入渗及其引发边坡失稳机制的研究基础上，对影响边坡降雨入渗过程的条件及入渗参数进行分析研究，并基于饱和-非饱和渗流数学模型和实测土水特征曲线，在延伸摩尔-库仑强度理论基础上对某堆积体进行了降雨过程的饱和-非饱和渗流计算和稳定分析，得到了一些具有实用价值的结论。
2 降雨在边坡中的入渗过程分析
降雨入渗过程实质上是入渗水分在非饱和区运动的过程，即水在下渗过程中驱替空气的过程［4］。研究表明，雨水渗入边坡土体的强度主要取决于降雨的方式和强度以及土体的渗水性能。如果土体渗水性能较强，大于外界供水强度，则入渗强度主要决定于外界供水强度，在入渗过程中边坡表面含水率随入渗而逐渐提高，直至达到一稳定值；如果降雨强度较大，超过了岩土体的渗水能力，入渗强度就决定于岩土体的入渗性能，这样就会形成地表径流积水。这两种情况可能发生在入渗过程的不同阶段。土体的入渗速率随时间而变化。一般来说，开始入渗阶段，边坡表层的水势梯度较陡，所以入渗速率较大，但随着雨水渗入土体中，土体的基质吸力下降，湿润层下移使基质吸力梯度减小，入渗速率逐渐减小最后接近于一常量而达到稳定入渗阶段。实际模拟计算中，降雨入渗边界虽然只是一个流量边界，但这个流量并不是不变的，在计算过程中需要根据含水率的变化不断的调整入渗流量，从而实现对降雨入渗问题的数值模拟。
据以上分析，可以确定降雨强度R（t）与实际入渗率ir（t）的关系。设入渗面的外法线方向为n=（n1，n2，n3），则降雨在入渗面法向上的入渗分量为qn（t）=R（t）·n3，由入渗面法向入渗能力
为整体坐标系中x，y，z 3个方向的入渗能力］，可确定入渗强度ir（t）与qr（t）及qn（t）的关系：
入渗过程中，累积入渗量、入渗率和岩土体含水率随时间的变化与地表处水的施加方式有关。为了求出入渗过程中土壤含水率分布，以及入渗率等随时间变化的定量结果，可以在一定的初始含水率分布条件下，根据入渗边界条件，求解水流运动基本方程。
本文发表于2008年。
3 雨水入渗引发边坡失稳的物理机制
如前所述，滑坡的发生与降雨关系密切早为人知，降雨入渗及地下水状态对边坡稳定的影响也是显而易见的，其作用主要有［2］：
（1）降低岩土体特别是有泥质充填物的断层、层间错动带以及以节理裂隙等软弱结构面为滑裂面的岩土体强度。
（2）地下水的静水压力一方面降低了滑裂面上的有效法向应力，从而降低了滑裂面上的抗滑力，另一方面岩土体非饱和区暂态附加水荷载又增加了边坡岩土体的下滑力，非饱和区饱和度增加，基质吸力（即毛细压力）降低，从而使边坡的稳定条件恶化，这也是雨季边坡失稳的重要原因［4］。
为定量研究因降雨入渗而导致的岩土体抗剪强度降低问题，引入非饱和土的抗剪强度理论［4-5］，该理论为定量计算因水分入渗而引起的岩土体软化的强度提供了一种计算方法，即Fredlund于1993年提出的岩土体吸水软化抗剪强度公式：
式中：c′，φ′分别为土的有效黏聚力和有效内摩擦角；（ua-uw）为孔隙气压与孔隙水压力之差；（ua-uw）tanφb为与基质吸力直接相关的抗剪强度，称为基质吸力的附加强度；φb为因基质吸力上升而引起抗剪强度增加的曲线倾角。
从式（2）可以看出，非饱和土的抗剪强度除了与c′、φ′及正应力有关外，还与基质吸力（ua-uw）有关。当土体饱和时，可以认为φb=φ′，退化为传统的摩尔-库仑公式，故式（2）又称延伸的摩尔-库仑公式。降雨时，岩土体饱和度和孔隙水压力上升，基质吸力（ua-uw）减小，抗剪强度明显减小，从理论上解释了降雨条件下滑坡发生的机理。
4 降雨入渗条件下的饱和-非饱和渗流场研究
目前国内外进行岩土体渗流场以及渗流应力耦合计算分析时，常根据实际工程水文地质条件复杂程度而采用不同的计算模型，主要有［6］：①等效连续介质模型；②裂隙网络非连续介质模型；③岩土体裂隙双重介质模型；④离散介质-连续介质耦合模型。大量工程研究和应用表明，尽管依据裂隙实测资料的统计数字由计算机生成等效网络已有可能，但三维裂隙网络水力学分析仍很难进行，故用离散网络模型来研究裂隙岩体非饱和渗流并不合适。采用等效连续介质模型。
4.1 堆积体非饱和水力参数的确定
对于饱和岩土类多孔介质，目前已有成熟的试验方法来确定其饱和渗透张量。对于非饱和岩土类多孔介质，关键是确定非饱和相对渗透系数与基质吸力、体积含水量与基质吸力的对应关系曲线，这两个对应关系曲线一般通过试验获取［7］。基于对江西某大（2）型水库库岸堆积体斜坡（崩坡积体）基质吸力的现场观测和观测竖井中的碎（块）石黏性土土原状样（取样的深度为1.0～5.0m）试验成果，得到该堆积体的试验参数如下：饱和含水量ωs=19.24%；天然含水量ωs=18.44%；天然状态饱和度S=95.84%；饱和容重γ饱和=22.5kN/m3；干容重为γd=17.5kN/m3；天然容重γ天然=21.5kN/m3。堆积体的土水特征曲线根据测试结果依据V-G模型拟合得到［7］，非饱和渗透系数依据土水特征曲线按下式确定［5］：
式中：k（θ）i为用第i个间断体积含水量θ确定的渗透系数；Ad为调整常数，一般取1；ks为饱和渗透系数；ksc为计算饱和渗透系数，计算公式为：
（ua
，i=1，2，…，m；ua，uw分别为孔隙气压和孔隙水压力。
由此得到非饱和渗透系数与基质吸力的关系，见图1（对数坐标）和图2。
图1 堆积体基质吸力—渗透系数关系图
图2 堆积体基质吸力—体积含水量关系图
对于基岩裂隙介质非饱和等效水力参数，可采用均化方法求取［7］，关于应用裂隙网络方法来确定裂隙介质的非饱和水力参数的方法可参见文献［8-9］。采用均化方法确定部分岩层的非饱和等效水力参数。其中假设：①裂隙充分发育，裂隙介质存在表征单元体且体积不是太大；②流动随时间变化缓慢，也就是指岩块饱和度变化不大，岩块和裂隙间的水量交换瞬时完成。
取一体积为V的表征单元体，其中裂隙体积V1，岩块体积V2，V=V1+V2。假定垂直于水流方向的平面上裂隙内水头和岩块内水头相等，则按照流量等效和水头近似等效的原则，可依据裂隙和岩块各自的饱和及非饱和水力参数、各自所占的体积，通过体积加权平均来确定裂隙岩土体的等效非饱和水力参数，则裂隙介质的等效相对渗透系数Kr、等效比容水度C及等效单元贮存量Ss分别为
其中
式中：
，
，
，
，
，
分别为裂隙网络和岩块饱和渗透张量的3个主值；
，
分别为裂隙和岩块的非饱和相对渗透系数；C1，C2分别为裂隙网络和岩块的比容水度；
，
分别为裂隙网络和岩块的单元贮存量。
4.2 渗流场数学模型的建立及有限元求解
4.2.1 饱和非饱和渗流数学模型
基于地下水运动的连续性方程和Darcy定律，略去推导过程，在不考虑骨架及流体压缩作用下建立的非稳定饱和非饱和渗流基本微分方程为
式中：h为总水头；hc为压力水头；
为饱和渗透张量；ki3为饱和渗透张量中仅和第3坐标轴有关的渗透系数值；kr为相对透水率，在非饱和区0<Kr<1，在饱和区kr=1；C为比容水度，在正压区C=0；β为饱和-非饱和选择常数，在非饱和区等于0，在饱和区等于1；Ss为弹性贮水率，在非饱和土体中Ss=0，当忽略土体骨架及水的压缩性时对于饱和区也有Ss=0；Q为汇源项。
4.2.2 定解条件
（1）初始条件
（2）边界条件
式中：ni为边界面外法线方向余弦；t0为初始时刻；Γ1为已知结点水头边界；Γ2为流量边界；Γ3为饱和逸出面边界。
4.2.3 有限元求解
对于上述数学模型，微分方程和定解条件往往都是非线性的，通常采用数值法求解。应用Galerkin加权余量法求解上述数学模型。
将计算空间域Ω离散为有限个单元，对于每个单元，选取适当的形函数Nm（xi），满足式中：hcm（t）为结点压力水头值。
将hc（xi，t）分别代入微分方程及边界条件，一般不能精确满足，分别会产生一定的误差，记残差值（余量）分别为R和
。按照加权余量法原理，选择权函数W（xi）及
，取W（xi）=Nn（xi），在边界上，取
，略去中间推导过程，对于离散后空间域Ω有：
令式（11）中
则式（11）变为
采用向后有限时间差分格式对时间域的积分进行求解，其基本方程为
根据求解的稳定性，采用tk时刻的压力水头值
近似作为Δtk时段的平均压力水头值，可得式（13）的隐式差分格式
5 工程应用
5.1 基本情况
江西省某大（2）型水库库岸堆积体位于坝址左岸上游约70～300m地段，边坡底部横向宽度约220.0m，纵向最大高度120m，厚一般12～15m，最厚近20m，总体积约30万m3。经勘探查明，堆积体物质组成以碎（块）石为骨架，内充填黏性土，碎（块）石、砾含量达40%～50%。堆积体下伏基岩顶部为石炭系岩层，一般为弱～新鲜岩体，岩层反倾；中下部为震旦系变余细砂岩夹板岩，其中顺坡向层间挤压带是边坡稳定的控制性结构面，见图3。
为定量分析边坡在强降雨条件下的稳定性，采用前述理论和方法，对其进行了强降雨过程的渗流场计算和稳定分析，其中降雨强度和雨型有：雨型1为降雨强度30mm/h，降雨持时24h；雨型2的选取系参照文献［10］中的降雨区划而定，降雨持时24h，其中前12h的降雨强度见图4，后12h的降雨强度和雨型同前12h，之后为无雨状态，相应库水位均取正常蓄水位180.0m，初始稳定渗流场系根据钻孔稳定地下水位确定，见图5。
图3 堆积体斜坡典型剖面图
图4 堆积体强降雨条件下的雨强和雨型分布图
基岩的非饱和水力参数参照4.1节所述方法，同时参照地质勘察中压水试验成果和相关地质描述确定，用于渗流计算分析的物理力学参数见表1。抗滑稳定计算采用可考虑基质吸力、适用于任意滑裂面的Morgenstern-Price法。
表1 堆积体斜坡抗滑稳定计算物理力学参数选取表
5.2 计算结果分析
限于篇幅，仅列出雨型2情况下降雨初始时刻的稳定渗流场和强降雨36h（雨后12h）时刻的暂态渗流场，见图5和图6。图7为堆积体斜坡安全系数随降雨过程的变化情况。由图可看出：
图5 堆积体斜坡初始稳定渗流场
图6 雨后12h时刻堆积体斜坡的暂态渗流场
（1）随着强降雨的进行，雨水不断入渗到堆积体中，非饱和带饱和度上升，暂态水荷载增加，基质吸力降低，同时地下水位略有上升，使滑裂面上的孔隙水压力增加，斜坡的安全系数随之逐渐减小。降雨停止后，因堆积体内的一部分孔隙水和基岩裂隙水缓慢排泄于库区和低洼地带，其安全系数有所回升。
（2）计算表明，在雨型1条件下，堆积体斜坡在强降雨24h时刻即降雨停止时刻安全系数达到最小值1.005，较降雨前降低6.34%；而雨型2条件下则是在强降雨36h时刻（即雨后12h）安全系数达到最小值0.989，较降雨前降低7.83%，这一现象与多数自然边坡在雨后一段时间产生滑动失稳的情况相吻合。可见降雨特征（雨强、雨型）对堆积体的降雨过程稳定性有一定影响，故采用雨型2时，计算结果更符合工程实际。
图7 堆积体斜坡降雨过程稳定性分布图
6 结语
（1）库岸堆积体斜坡在我国许多水电站地区分布广泛，其稳定性问题关系到近坝库岸和大坝的安全。降雨是导致库岸堆积体变形发展和滑动的主要外动力因素之一，研究因降雨入渗引发库岸堆积体滑坡的机制和规律具有重要意义。
（2）基于饱和-非饱和渗流数学模型和实测土水特征曲线，运用延伸的摩尔-库伦强度理论，对某库岸堆积体进行了降雨过程的饱和-非饱和渗流场计算和稳定分析，得出其最小安全系数发生在雨后一定时间内的结论，这与多数自然边坡在雨后一段时间产生滑动失稳的实际情况相吻合，较为合理地阐述了因降雨引发的边坡失稳问题。同时研究表明，在模拟因降雨入渗引发边坡失稳过程时，应优先采用实际可能降雨强度和雨型，降雨特征（雨强、雨型）对边坡稳定的影响不容忽视。
参考文献
［1］田陵君，王兰生，刘世凯.长江三峡工程库岸稳定性［M］.北京：中国科学技术出版社，1992.
［2］孙广忠.中国典型滑坡［M］.北京：科学出版社，1988.
［3］吴宏伟，陈守义，庞宇威.雨水入渗对非饱和土坡稳定性影响的参数研究［J］.岩土力学，1999，20（1）：1-14.
［4］包承纲.非饱和土的性状及膨胀土边坡稳定问题［J］.岩土工程学报，2004，26（1）：1-15.
［5］Fredlunnd D G.非饱和土土力学［M］.陈仲颐，译.北京：中国建筑工业出版社，1997.
［6］戚国庆，黄润秋，速宝玉，等.岩质边坡降雨入渗过程的数值模拟［J］.岩石力学与工程学报，2003，22（4）：625-629.
［8］胡云进，速宝玉，詹美礼.裂隙岩土体非饱和渗流研究综述［J］.河海大学学报，2000，28（1）：45-51.
［9］王恩志，王洪涛，孙役.三维裂隙网络渗流数值模型研究［J］.工程力学，1997（增2）：520-525.
［10］江西省山洪灾害防治规划领导小组办公室.江西省山洪灾害防治规划报告［R］.2004，12.
本文发表于2008年。