1. 引言
水库建设及运营期间水位的涨落变化，实则是库岸边坡边界条件的变化并间接导致坡脚区域组成物质力学性质的改变，这些因素的变化打破了库岸边坡岩土体原有的应力平衡，应力将发生重新分布，使坡体某些局域的应力集中或安全可靠度的下降，最终造成边坡区域的破坏和失稳[1]，如意大利的瓦伊昂水库在蓄水3年后发生滑坡，在大坝附近2亿多m3的岩土体迅速滑入水库，库水涌浪翻过大坝流向下游致使2600人丧生[2]以及国内的三峡水库蓄水期间发生的著名的千将坪滑坡，致使24人遇难，1100多人无家可归，334省道被毁，交通被迫中断[3]。因此，研究库水位升降条件下岸坡稳定性变化是十分必要的，目前国内外针对库水位升降作用对边坡稳定性的影响作了大量的研究工作，张夏冉以重庆万州下坪滑坡为例，将暂态孔隙水压力分布应用到滑坡极限平衡分析中，得到了滑坡在不同库水位降速、不同降雨量以及不同降雨时长下的稳定性系数和破坏概率[4]，Manzoor通过详细现场调查和反演分析研究了蓄水作用下滑坡复活的成因机制[5]，Guo采用实证模型研究了库水位升降作用下的滑坡稳定性变化规律[6]，Kmac和Benoit采用现场监测方法研究了库水位升降条件下库岸滑坡水力触发因素[7] [8]。堆积体滑坡在库水浸泡下，岩土体的物理力学性质不断恶化[9] [10]，受库水浮托力和坡体内渗透力的影响，坡体稳定性不断改变。此外，滑坡是复杂的多相系统，内部渗流场和应力场相互作用、互相影响[11] [12]，分析库水位波动条件下的滑坡渗流场特征和稳定性变化规律，是研究水库蓄水运行对库岸堆积体滑坡稳定性影响的重要内容[13]。
龙羊峡水库区属高原干旱–半干旱大陆性气候，多年平均降水量243.0~404.0 mm，多年平均蒸发量1327.2~1378.5 mm，最大季节性冻土深度为2.4 m。降雨多集中在夏季，夏季降雨约占全年降雨量的75%以上。研究区地质环境条件复杂，崩塌、滑坡、塌岸等边坡失稳灾害发育，尤以滑坡最为严重。水电站自1986年建成投产以来，库水位的长期波动变化对库岸边坡稳定性造成了重大影响，库区地质灾害不断凸显，如2002年6月汪什科滑坡东侧发生了600 × 104 m3巨型滑坡[14]；2005年11月在高库水位持续阶段，沙沟加达口–文昌庙南侧西岸产生了约450 × 104 m3的高速滑坡[15]，以及白刺滩滑坡、农场滑坡、水磨沟北侧滑坡、龙西滑坡、查西滑坡等等[16]。近年来黄河上游汛期来水较丰，龙羊峡水库水位持续上升，水库库岸存在潜在的不稳定高边坡，在库水位的反复浸润、浪蚀作用下局部地段已有变形，这些边坡一旦发生大规模的滑坡，对居住于沿岸的人民群众生命财产的安全将构成重大危害。因此，研究龙羊峡库区边坡稳定性随水位波动的变化规律有助于迅速、科学地评估龙羊峡水电站库区边坡的潜在重大危险性。
2. 龙羊峡库岸边坡工程地质概况
2.1. 区域地质背景
龙羊峡水电站库区位于黄河上游干流青海省海南州境内，库区在区域构造上属于青藏高原内部南祁连–西秦岭海西褶皱带共和盆地东部的菱形拗陷块体中。此外，龙羊峡水电站所处的共和盆地沉积着大量半成岩土，系上新世、更新世早中期内陆湖相沉积的产物，受青藏高原晚第三纪以来大幅抬升和差异性升降运动的控制，在黄河凿穿龙羊峡口后，由于区域性的强烈抬升运动和黄河的深切作用，在近坝库区南岸形成了一系列广范围的高陡边坡，这些高陡边坡的孕育为库区大型滑坡的发生奠定了必要的坡形基础。
2.2. 地层岩性
库区出露的地层主要为第四系河湖相松散堆积物，库位零星分布有二叠系下统、中统碳酸盐岩、三叠系中下统碎屑岩以及印支期侵入岩，第四系地层在库区广泛分布，成因类型复杂，岩性多样。库区高陡边坡内的半成岩土性质与内陆湖相沉积和当地干旱少雨的气候密不可分，有学者经过土化学试验后表明，半成岩土中可溶性钙盐和钠盐含量较高，一般高于20%，其中难溶性钙盐所占比重极大，岩土中的微碎屑和粘胶微粒，经微晶碳酸钙胶结成集粒，与老黄土的物质组成颇相同(吴其伟，王成华[17]，1986)，此外这些半成岩土流塑性上的中、低液限性质和遇水易于崩解的特征即是其反映，所以这类岩土容易受水的影响。
2.3. 水文地质条件
龙羊峡水库自1986年建成蓄水以来，水库水位从7月份开始上升，到11月份达到峰值，随后开始逐渐下降，直到次年春季降到最低点。库区地下水埋藏深，且流量小，主要埋藏于湖相地层下部I大层以及III大层下部砂层及砂壤土中。向岸内300 m处地下水位高程多在2510 m以下，低于坡顶400~500 m。地下水属孔隙–裂隙水，向黄河排泄，水力坡度为3%~8%。由于本区干旱的气候环境，地下水动态变化受大气降水的影响很小，以接受远处地下水的稳定补给为主[18]。经学者多年来的监测调查结果表明，同期地下水位实际上升是库水位上升幅度的7%~36%，水位年变幅一般小于5 m，地下水位的上升幅度还受库/地下水位差的影响，当库水位持续上升，库/地下水位差大于30~40 m时，地下水位明显上升，反之，地下水位无明显变化，反映了河湖相地层低渗透性的特征[19]。
2.4. 岸坡滑动历程
水库运行以来的主要滑坡隐患来源于库区的古滑坡以及老滑坡，在近四十年的浸泡以及库水位波动作用下，库岸边坡经历三个过程(如图1)。
A. 建坝初期水库蓄水导致库岸边坡坡脚浸润，而近坝库岸的高陡边坡上部覆盖层为遇水易软化的半成岩土，在水的作用下，经历漫长的坡脚冲刷、卸荷、应力集中，直至潜在滑动面贯通，发生滑坡，上部覆盖层滑落，滑坡方量大的巨型滑坡易产生堵江、涌浪等危害。
Figure 1. Slope change process
图1. 岸坡变化过程
B. 在上部覆盖层滑落后，下部较硬的岩石层出露，但很快被滚落下来的碎石土覆盖，边坡处于一个暂时稳定的阶段，随后受库水位升降影响，坡脚继续受到冲刷，致使覆盖碎石被冲走，岩层浸水不断软化。
C. 受库水冲刷，于坡脚岩体处形成掏蚀穴，在岩体临空面增加，切岩体浸水软化的影响下，导致坍岸发生，为后续边坡破坏或二次破坏孕育条件。
3. 有限元模型建立
3.1. 计算模型及边界条件
龙羊峡库区潜在滑坡主要以古、老滑坡为主，近年来除了部分滑坡发生如后缘拉张裂缝的延展、滑体四周发生小型坍塌、滑体前缘坡脚处出现土体凸起现象以及滑坡位移监测数据增大等迹象以外，并没有大型滑坡的发生。因此，本文选择距坝6.5 km的查纳以及查东滑坡为典型案例分析(如图2)，其中查纳滑坡曾于1943年发生大规模滑动，下滑量达到1.6亿m3，滑体前缘向前推进了约3 km，埋没了黄河岸边的树林带以及查纳村，将其作为老滑坡典型案例；查东滑坡位于查纳滑坡后缘东侧，坡向40˚，体量相对查纳滑坡略小，全新世以来并无滑动历史。地下水位于模型原点坐标以下，初始应力的产生与坡体自重及地下水有关。
Figure 2. Study the relative position of slope
图2. 研究边坡相对位置
根据研究成果及地勘数据资料，本文借助二维有限元软件进行数值模拟分析，所建模型如下图(如图3)：查纳滑坡模型高365 m，长1050 m，对Xmin位置以及Ymin位置处设置固定约束，划分单元网格1682个，节点13,845个；查东滑坡模型高320 m，长900 m，对Xmin位置以及Ymin位置处设置固定约束，划分单元网格2018个，节点16,437个。
Figure 3. Two-dimensional finite element model ((a) Chana landslide and (b) Chadong landslide)
图3. 二维有限元模型((a) 查纳滑坡与(b) 查东滑坡)
3.2. 库水位变动特征
龙羊峡水库具有多年调节性能，早年水库长期在低水位2575 m以下运行，在2005年11月19日库水位达到历史最高值2597.62 m之后，水库进入高水位运行阶段。图4为利用水利部黄河水利委员会公布龙羊峡水库水位数据绘制的近五年来水位高程变化曲线。从图4可以看出，库水位在一个水文周期年内，1~6月份呈下降趋势，其中4、5、6三个月份下降趋势较快，并在6月份达到全年库水位最低值，随后开始呈上升趋势，7、8两个月上升趋势较快，到10、11月份水位上升至峰值，并保持稳定直到下一个水位年。
3.3. 岩土体强度参数选取
研究区地层岩性多为第四纪中、早更新世湖相地层，可划分为七大层(如图3)，坡体上部主要为砂性土，中、下部主要由密实度很高的、超固结呈半成岩状的粘性土夹薄层砂土类地层组成，岸坡下部主要为粘性土，厚度约250 m，内部有多次沉积韵律的变化，夹有多层0.30~3.0 m的薄砂层，是库水浸泡的主要地层，上部主要为砂性土，总厚150~300 m，根据地勘资料及土工试验数据，将各层岩土体强度参数取值见表1。
Figure 4. Changes of reservoir water level and elevation of Longyangxia Hydropower Station in recent five years
图4. 龙羊峡水电站近五年库水位高程变化情况
Table 1. Strength parameters of rock and soil mass
表1. 岩土体强度参数表
地层
划分 |
天然容重(kN/m3) |
饱和容重(kN/m3) |
弹性模量(MPa) |
孔隙比 |
天然状态内摩擦角Φ’ |
饱和状态内摩擦角Φ’ |
粘聚力C’ (kPa) |
渗透系数 |
Ⅰ |
21.1 |
20.9 |
479 |
0.452 |
35˚ |
33˚ |
323 |
8.64 × 10−5 |
Ⅱ |
18.7 |
20.6 |
468 |
0.554 |
37˚ |
35˚ |
317 |
8.58 × 10−5 |
Ⅲ |
21.1 |
21.8 |
347 |
0.459 |
35˚ |
33˚ |
88 |
8.36 × 10−2 |
Ⅳ |
18.7 |
20.59 |
325 |
0.554 |
37˚ |
35˚ |
245 |
6.96 × 10−5 |
Ⅴ |
22.0 |
22.8 |
339 |
0.355 |
32˚ |
30˚ |
79 |
6.73 × 10−5 |
Ⅵ |
17.3 |
20.3 |
222 |
0.573 |
38˚ |
- |
39 |
0.865 |
Ⅶ |
17.5 |
20.5 |
243 |
0.594 |
38˚ |
- |
36 |
4.6 |
3.4. 计算工况选取
根据2021年一整年水位变化趋势(如图5)，可细分为平稳运行阶段(0 m/d)、水位慢速上升阶段(上升速率+0.03~+0.04 m/d)、水位中速上升阶段(上升速率+0.05~+0.08 m/d)、水位快速上升阶段(上升速率+0.10~+0.15 m/d)、水位慢速下降阶段(下降速率−0.02~−0.05 m/d)、水位中速下降阶段(下降速率−0.06~−0.09 m/d)，水位快速下降阶段(下降速率−0.12 m/d)共7种类型24个阶段。根据每个阶段水位和升降速率，分17种工况进行计算(如图5)。
Figure 5. Water level fluctuation and working condition fitting of Longyang Gorge reservoir in 2021
图5. 2021年龙羊峡库水位波动及工况拟合
4. 数值模拟结果及分析
4.1. 增量位移计算结果
位移增量云图是最后一个计算步相对上一个计算步的增量，它的分布能够揭示边坡正在破坏的机制以及滑弧(潜在滑动面)的位置[20]，本文通过分析滑坡增量位移云图来揭示极限状态下的滑弧，两种滑坡的增量位移如下图(如图6)：
Figure 6. Landslide incremental displacement cloud map ((a) Chana landslide and (b) Chadong landslide)
图6. 滑坡增量位移云图((a) 查纳滑坡与(b) 查东滑坡)
查纳滑坡不同工况下的潜在滑动面位置差别不大，岸坡内部位移由内到外逐渐变小，呈高–中–低的分布梯次，大位移产生位置集中在滑坡前缘中下部，相比查纳滑坡，查东滑坡属于古滑坡，并无滑动历史，上部覆盖较厚，变形区域较大，滑弧较长，且岸坡内部变形分布不均，呈低–高–中的分布梯次，最大增量位移区位于滑弧中下部，接近坡脚。本文以24个阶段计算得到增量位移最大值为依据，绘制了两种滑坡的增量位移变化曲线(如图7)：
Figure 7. Incremental displacement curves at each stage
图7. 各阶段增量位移变化曲线
查纳滑坡的增量位移变化曲线显示了其内部变形的不均匀性，阶段1~3属于水位下降阶段，下降速度为慢速~中速，在图中并没有显示出较大的位移量差别，但是从阶段4开始，位移增量发生大幅度变化，直到阶段6结束才恢复变化前水准，这说明岸坡内部的变形存在滞后和累计性，这两种特性，在随后的阶段8~阶段17同样有所体现，库水位在阶段7发生了最快速度的下降，并且持续时间长，通过观察水位波动图(如图5)可知库水位在阶段8~阶段17处于一个相对平稳的变化阶段，既有水位上升也有水位下降，甚至还存在相当长一段时间的稳态水位，因此之前累计的变形在这一时期得到了释放，由于水位发生了大幅度的变化，造成的影响就是岸坡内部变形随之增大，图中表现位大幅度的增量位移变化，并一度达到了整个变化过程中的峰值，随后在阶段17~22，库水位开始了整体上升阶段，上升速度由快到慢，水位同样发生了较大幅度的变化，但增量位移虽然也发生了不同程度的变化并没有表现出水位下降时期的高度，这说明水位下降阶段对岸坡位移的影响大于水位上升阶段。
而查东滑坡同样表现出了坡体内部变形的不均匀性，不同于查纳滑坡在阶段4就发生明显的上升，查东滑坡则是在阶段8才发生增量位移的大幅上升，阶段1~7呈现出上下波动的特征，说明在前7个阶段，水位的变化并未对岸坡造成明显的影响，亦或是查东的滑坡变形滞后性更甚于查纳滑坡，其内部坡体结构较查纳滑坡更加稳定，因此直到阶段8才发生明显的变化，这一轮变化的持续时间也比较短，紧随而来的是第二轮大幅度变化，持续时间为阶段11~14，这一轮变化的幅度较小，但持续时间较长，在第二轮大幅变化后，同查纳滑坡相似，在阶段15~18开始了变形孕育期，这一阶段的增量位移变化不大，此后查东滑坡又经历了两轮大幅变化过程。相较于查纳滑坡，查东滑坡增量位移波动周期多，多达四轮，结合查东滑坡的滑弧范围，可见虽然查东滑坡虽然没有滑动历史，但其内部变化并非想象中的那么稳定，这也与调查人员在查东滑坡后缘发现拉张裂缝的现实情况相吻合。
4.2. 稳定性计算结果
两种滑坡在不同工况下的稳定性计算结果见图8，查纳滑坡总体上水位上升阶段的稳定性系数(安全系数)要高于水位下降阶段的稳定性系数，且查纳滑坡在所有工况下都处于比较稳定的状态，在阶段21达到了峰值1.833，最低点位于阶段8和阶段10，两种阶段由于计算精度的原因稳定性系数都达到了1.789，，关于阶段8和阶段10达到最低稳定性系数的情况，增量位移变化曲线图(如图7)已经给出了解释，阶段8处于上一累计变形释放的结束以及新一轮累计变形的开始，而阶段10则达到了增量位移的峰值。纵观查纳滑坡在水位整体性上升阶段的安全系数，则没有先前阶段的“纠结”，呈现出比较简单的变化规律，即随着水位上升和下降，稳定性系数先增大后减小，且上升速率越快，稳定性系数增长得越快。
查东滑坡稳定性变化曲线与查纳滑坡的变化曲线呈现出了镜像的变化规律，查纳滑坡稳定性系数大幅度变化存在于阶段19以后，而查东滑坡在库水位变化前中期就凸显出了大幅变动的特征，稳定性系数最高位于阶段9，达到了1.816，最低位于阶段7，低至1.781，两种边坡截然不同的稳定性系数变化曲线，反映了即便是临近滑坡，地层基本一样，也会出现不同的内在变化规律，查东滑坡的稳定性系数变化主要集中在库水位快速下降阶段以及中期的库水位反复升降阶段，不同工况下的稳定性系数没有明显变化规律，也说明查东滑坡目前还处于孕灾阶段，内部并没有产生贯通的滑动面，受内外在因素所影响的滑动力与抗滑力仍处于拉锯状态，岸坡稳定性系数并没有随库水位波动表现出明显规律。
Figure 8. Variation curve of landslide stability coefficient
图8. 滑坡稳定性系数变化曲线
4.3. 计算结果分析
通过观察上文两种边坡的位移增量云图可以发现，对于库水位不同下降速度的工况，最终获得潜在滑动面的形状、范围、深度差别不大，考虑到原因是由于计算中设置的阶段时间较短，仅仅显示出了变形趋势，对滑弧位置的影响还未来得及体现，此外，由于查东滑坡坡形向外凸起，查纳滑坡存在滑动历史，表层发生过滑移，而显得略微凹陷，从而查东滑坡的滑弧不论是深度、大小都远大于查纳滑坡，总体趋势是库水位下降速率越大，滑弧深度越大，此时比较容易在深层发生滑动；而在水位上升阶段，随着水位上升速率的增加，坡体变形集中区由滑弧面逐渐向外部推移，且从坡脚向坡肩变形逐渐变小，缓慢上升阶段的下降梯度较为明显，推测由于坡体内地下水位随库水位不断抬升，内部出现孔隙水压力且方向由坡外指向坡内，且快速上升阶段的孔隙水压力消散速度慢于缓慢上升阶段，导致坡内变形传递较快[21]。
将库水位变化散点图与两种边坡不同阶段的稳定性系数变化汇总到一起(见图8)，进行综合分析，可见查纳的滑坡稳定系数变化规律与库水位的变化规律相当符合，总体上的水位降低阶段的稳定性系数低于水位上升阶段，在水文周期前期，库水位以不同的速率逐渐降低，查纳滑坡的稳定系数也相应的呈下降趋势，水位快速下降阶段的稳定性系数小于中速下降阶段小于慢速下降阶段，其原因是在库水位下降过程中，坡内的渗流方向是由坡内指向坡外，导致动水渗透压力指向坡外，从而增大了坡体下滑力，以及岸坡浸水范围增大，岩土体遇水软化导致物理力学参数降低的结果[22]。此外，库水位下降速率越快，坡体稳定性达到最差的时间也就越短，稳定值也最小，当库水位下降速率较大时，坡体稳定性以静水压力为主导作用，随着库水位的不断降低，动水压力不断增大，边坡的稳定性也会不断降低，如果库水位下降速率较小，相对的动水压力的增大速率也会略低，导致坡体稳定性整体高于快速下降阶段；当水位进入中期的平稳运行阶段，库水位在这一时期的变化幅度不大，既包含上升阶段，也包含下降阶段，还存在水位保持阶段，相应的查纳滑坡的稳定性系数在这一时期也变幅不大，其原因在于库水位在短期内变幅小、频率高，坡体内部渗流场变化存在滞后作用，当水位下降结束并开始保持不变或者开始上升时，坡内浸润线还未来得及变化，坡外水位就已经发生变化，导致坡内浸润线高度下降幅度无法与之前相比，从而主导这一时期稳定性系数变化的是前期积累的位移；当水位进入整体上升阶段时，查纳滑坡稳定性系数也随之提高，当库水位以不同速率上升过程中，边坡的稳定性先迅速增大，再缓慢减小，最后达到逐渐稳定的状态，分析其原因，可能是上升过程中的渗流方向是由坡外指向坡内，导致动水渗透压力指向坡内，从而增加坡体的抗滑力，但是如果库水位以较小的速率上升，坡体则会受到静水压力、动水压力、浮托力以及岩土体的物理力学性质的变化的共同影响，加之孔隙水压力有相对充足的消散时间，导致坡体抗滑力的增大不明显，当水位逐渐稳定时，库水会不断渗入滑坡体内部，并造成滑带部位岩土体物理力学参数降低，这也是稳定性先增大后减小的原因[23]。反观查东滑坡的稳定性系数变化情况，与库水位变化关联性不大，查东滑坡的稳定性系数与库水位阶段性的升降、平稳趋势不甚相同，除了前中期存在一个骤降和骤升的结果，整体处于一个相对稳定的阶段，其骤降存在于库水位最快速下降阶段，其原因与查纳滑坡一致，这里不再赘述，而随后的骤升过程，则可能是查东滑坡独特的坡体结构在坡内历经大幅度的变形之后所进行自我修复，这种机制的存在也是查东滑坡迄今为止并没有发生滑动的原因，但两种临近岸坡、地层相仿却展现出了截然不同的变化规律，由专家曾指出：边坡形态对坡体整体稳定性的影响也至关重要，边坡角和边坡高度是边坡形态最主要的两大要素[24]，本文以查纳和查东两种坡形不同的边坡为例，参照综合对比图(图8)对稳定性计算结果进行对比分析，查纳滑坡相对于查东滑坡坡高较高，但坡角略缓，查东滑坡上部覆盖层未曾滑落，导致上部凸起，根据计算结果显示，当库水位较高时，查纳滑坡稳定性高于查东滑坡，推测是由于库水位较高时坡面浸润位置临近上部覆盖层底部，致使覆盖层产生变形，而查纳滑坡由于覆盖层已经滑落，影响相对较小，反观库水位较低时，两种边坡稳定性差距不大互有高低，究其原因可能是两种边坡中下部结构类似，坡形差距主要体现在边坡上部，故而产生了这种情况[25]。
5. 结论
本文以龙羊峡水电站的两个边坡为研究对象，利用有限元分析软件，结合地质勘察数据，对库水位波动过程中滑坡体内部的渗流场、变形场变化规律进行了分析，研究结论显示：
1) 当库水位呈上升状态时边坡稳定性较高，水位上升速率越大，边坡稳定性系数越高；当库水位呈下降状态时边坡稳定性较低，稳定性系数随下降速率的提高而降低。
2) 当水库以低水位状态运行时边坡稳定性较低，而以高水位状态运行时，由于动水渗透压力起到一个增大抗滑力的效果，加之水对岩石软化的作用还未产生质变，因而边坡稳定性较高，边坡稳定性也会随着水位下降速率的提高而降低。
3) 查纳滑坡的稳定性随水位波动均匀变化，而查东滑坡稳定性保持相对稳定，其稳定性系数与库水位波动关联性不大。
两种类型的边坡稳定性计算结果不论在什么工况下都是 > 1的，与龙羊峡水库蓄水以来的实际情况相符合，近年来并没有滑动历史，但是随着库水位的变化，库水位升降会导致边坡稳定性系数有不同的变化特征，根据学者研究龙羊峡其他滑坡的孕育、诱发情况，滑坡的发生需要相当长的周期，产生因素和诱发因素众多，不单单是库水位波动这一个因素所能概括的，库水位波动使坡内渗流场发生变化，坡内变形增大，存在累积效应，不断为边坡耐久度降低积累因素，最终在多种因素的耦合作用下，边坡累计位移达到阈值，坡内滑动面贯穿，滑坡发生，对周围人民的生命财产安全造成威胁，后续应该在本文研究内容的基础上，多角度，多层次，多深度的分析致使边坡稳定性发生变化的因素，进而通过综合分析，结合边坡监测数据，及时掌握滑坡体的变形动态，以便对其变形趋势进行判断和处置。
基金项目
第二次青藏高原综合科学考察研究(2019QZKK0904)资助。