引發某近壩庫岸滑坡失穩臨界蓄水位研究

余政興，孫 寧

(中國電建集團中南勘測設計研究院有限公司，湖南 長沙 410014)

自1963年意大利瓦依昂滑坡發生后，庫岸滑坡，尤其是近壩庫岸滑坡一直是國內外學者研究的重點[1]。水庫蓄水后，近壩庫區水位變化幅度最大，滑坡體前緣被水淹沒后穩定性降低，若發生失穩，對水庫大壩的安全運行和下游城鎮均會造成巨大威脅。

長期以來，大量學者在庫岸滑坡的形成條件和預測評價等方面進行了大量研究。Zou Zongxing等[2]基于對三峽庫區滑坡監測數據分析，提出了水庫滑坡變形機制的解釋模型；李卓等[3]對近壩庫岸滑坡進行了模型試驗研究，揭示了水庫運營情況下滑坡土壓力、孔隙水壓力、基質吸力等變化，對此類滑坡破壞特點進行了總結；張永剛[4]基于流固耦合的三維數值模擬方法，對三峽庫區不同降水速率滑坡破壞過程進行了研究；吳泉澳[5]、王灝[6]等基于FLAC3D的強度折減法，對不同降雨量、不同蓄水位下滑坡的穩定性進行了研究。

以上研究均建立在水庫建設期或運營期的庫水位調度對滑坡影響的基礎上，而在蓄水位尚未確定的規劃階段，針對不同設計蓄水位對庫岸滑坡的影響少有研究。若在水庫或流域梯級規劃階段就考慮對近壩庫岸不良地質體問題，則可通過調整設計方案等方式提前預防，能極大減輕后期治理工程的工作量。本文在前人研究的基礎上，基于某水電站近壩庫岸大型滑坡，在假定蓄水位尚未確定的情況下，對引發滑坡失穩的臨界蓄水位進行研究，以期為水庫規劃及后期防治工作提供參考。

1 研究區概況

1.1 工程概況

某水電站地處藏區，為I等大(1)型工程，正開展前期勘察工作。該工程在流域規劃階段曾采用310m的高壩方案，蓄水位2605.00m，但考慮到高壩建設難度大等原因，目前設計采用低壩方案，正常蓄水位2535.00m，死水位2530.00m。由于區域內工程地質條件復雜，河谷狹窄，岸坡陡峻，在近壩庫岸形成了多處不良地質體。其中Z1滑坡體位于擬選壩址上游左岸，距離壩址約1.7Km，為壩址區內規模最大、危害性最大的滑坡體。

1.2 基本地質條件

1.3 氣象水文

研究區地處高山峽谷區，為干熱河谷氣候，氣溫高，降雨少，由于河谷內發育局部小氣候，與附近氣象站差異較大，故對區內氣象進行實測。2014—2020年實測多年平均氣溫16.0℃，最高氣溫36.7℃，多年平均降水338.3mm，最大日降雨量33.5mm。研究區內地下水主要為基巖裂隙水，鉆孔揭露地下水位較低，基本與河水位平齊(河面常水位2367m)，監測顯示，暴雨期間地下水位上升幅度最大可達15m，但仍位于滑坡體以下。

1.4 滑坡體發育特征

Z1滑坡體前沿高程2500m，后緣最高約3110m，最大高差達610m，平面分布面積約35.6萬m2，滑體一般厚度28m左右，總方量約996.8萬m3，為大—巨型滑坡。滑坡中后部滑體較厚，揭露最大厚度可達57m，為中—深層推移式滑坡。滑面地形整齊、坡度35°左右，以中間沖溝為界將滑坡體分為Ⅰ、Ⅱ兩區，其中Ⅱ區后緣見少量裂縫。滑體、滑帶組成物質以碎石土為主，滑床為堅硬的花崗巖，滑坡體邊界清晰，自然條件下基本穩定。

1.5 成因機制分析

滑坡體目前處于蠕滑變形階段，后緣見少量拉裂縫，其產生變形的原因主要為滑體規模較大，結構松散，在重力作用下，滑體與基巖面產生相對位移。除此之外，滑坡的形成條件還包括了地形地貌、地質構造、巖土性質、水文地質等內在因素，以及地震、降雨、地下水位變化、河水位變化、人類活動等外在誘因。

針對Z1滑坡，在天然條件下，自重應力對滑坡體的形成起主導作用；而降雨、地震、庫水位作用等是促使滑坡體失穩的誘發因素。由于庫區降雨稀少，附近不具備發育6級以上地震條件，周圍亦多見有垂直自穩高度超10m的堆積體，故地震與降雨對滑坡體影響有限，僅在極端情況下考慮。水電站運營蓄水后，庫水位變化與滑坡體穩定性關系最密切，其主要體現在：庫水位的上升使坡體巖土體浸水軟化，物理力學強度降低；庫水位上升對坡體產生浮托力；同時庫水位變化使坡體內地下水位發生改變，進而使得坡體內產生滲流作用，影響庫岸滑坡的穩定性。

故本文主要從庫水位變化對庫岸滑坡形成的影響角度出發進行研究。

2 不同蓄水位條件下滑坡體穩定性評價

2.1 滲流求解

滲流求解的主要目的是確定某一時刻的浸潤線，浸潤線確定后可以通過水頭求解坡體內孔隙水壓力的變化，為坡體后期的強度、變形及穩定性分析奠定基礎。浸潤線的求解方法主要有Boussinesq方程法[7]、Forchheimer方程法[8]、實際監測法、以及基于SEEP/W模塊的滲流—應力非完全耦合法[9- 10]。其中Boussinesq方程法、Forchheimer方程法求解較為困難，故本文首先基于監測數據推測穩定蓄水后的坡體的浸潤線，再通過Geostudio的SEEP/W模塊，對庫水位變化后的滲流過程進行二維數值計算，進而得到庫水位變化后的浸潤線。

由于地區降雨量稀少，實測發現地下水位基本與河水位平齊，故假設正常蓄水穩定后坡體地下水位仍與庫水位平齊，而在水庫運營時地下水隨庫水位的變化而改變。據統計，約60%的滑坡發生在庫水位驟降時期，其余40%發生在水位上升或穩定蓄水時期[11]，通常又以水位驟降時滑坡穩定系數最低[12]，故本文僅對不同設計蓄水位下水位驟降時滑坡的穩定性進行對比。由于該水電站是以發電為主的日調節水庫，可調節水位5m左右，故對穩定蓄水后水位驟降5m(驟降速率為5m/d)的浸潤線進行求解，考慮到滑坡Ⅰ區位置相對較低，更易受到庫水位影響，故取A—A＇剖面作為計算剖面，蓄水位取2535m時初始孔隙水壓力分布如圖1所示。

圖1 初始孔隙水壓力分布

進行瞬態分析前需要定義土水特征曲線，即水力傳導率函數曲線和體積含水量函數曲線，其中土體參數設置如圖2所示。根據滲透室內試驗和壓水試驗結果設置滲透系數，土體飽和滲透系數Ks=4.8m/d，基巖Ks=0.01m/d。

圖2 土體滲流參數設置

定義臨水坡面為水頭邊界，底部為不透水邊界，左右垂直端為零流量邊界，對模型進行瞬態滲流求解(降水5m，速率為5m/d，時間24h)，滲流進行24h后浸潤線變化情況如圖3所示。可見坡體內部水位有所降低，而基巖、滑帶處水位降低不明顯。

圖3 降水24h浸潤線

2.2 穩定系數計算

在滲流計算完成的基礎上，采用傳遞系數法對滑坡整體的穩定性進行計算，基于室內試驗、參數反演、經驗等取值，計算參數見表1，其中浸潤線以上取天然參數，以下取飽和參數。由于該水電站曾分別考慮過高壩方案(蓄水位2605m)和低壩方案(蓄水位2535m)，現以10m梯度分別求解蓄水位2535～2605m之間庫水位驟降5m時滑坡整體的穩定性。

表1 滑坡穩定性計算參數

計算結果如圖4所示，在蓄水位為2535～2605m時，滑坡穩定系數隨水庫蓄水位的增高而降低，平均蓄水位每升高10m，滑坡穩定系數降低約0.017。在目前采用的低壩方案(蓄水位2535m)中，水庫蓄水后滑坡穩定系數為1.071，處于基本穩定狀態；而若采用高壩方案(蓄水位2605m)，蓄水后滑坡穩定系數僅為0.948，極可能發生失穩。故僅從對近壩庫岸穩定性影響考慮，采取低壩方案更有利。當蓄水位位于2465～2575m時，滑坡穩定系數在0.997～1.012之間，處于臨界穩定狀態，內插可得滑坡失穩的臨界蓄水位約2571m。

圖4 不同蓄水位穩定性計算結果

3 基于FLAC3D的滑坡三維數值模擬

由于滑坡規模較大，實際地形地貌與地質條件較復雜，滑坡變形破壞模式還不明確，而二維計算僅能代表滑坡剖面附近的穩定性情況，有較大的的局限性，故需采用三維數值模擬軟件對二維計算結果進行驗證和推廣。本文采用FLAC3D有限差分數值軟件，對不同蓄水條件下滑坡穩定性進行評價，并對滑坡破壞過程、破壞模式等進行演化。

3.1 模型建立

采用外接軟件建模后導入FLAC3D進行計算，模型如圖5所示，分為滑體、滑床、基巖3個部分，劃分網格數554244個，對滑帶、滑體內的網格進行了加密處理。

圖5 三維計算模型

滑坡的荷載組合從兩個方面考慮，即自重力和水作用力，其中水作用力考慮水庫水位變化對滑體穩定性產生的影響，浸潤線參照2.1節設置，滑坡迎水面施加等效于庫水壓力的應力。分別對未蓄水及蓄水位2535、2571、2605m四種情況下水位驟降5m工況進行計算，除表1所給的巖土參數外，其余計算參數見表2。

表2 數值模擬計算參數

本次計算模型選用莫爾-庫侖模型，計算滿足莫爾—庫侖破壞準則。

3.2 計算結果分析

FLAC3D計算完成的標志是最大不平衡力趨于0，如圖6所示，當最大不平衡力減小至1x10-6Pa以下，說明在天然重力場條件下，模型通過自身的調整可以達到自平衡，天然應力場形成。

圖6 最大不平衡力曲線

首先對蓄水前滑坡應力、應變場情況進行分析。未蓄水時滑坡的應變云圖如圖7所示，可見當前滑坡變形主要集中在坡體中后部，前緣為阻滑段，與推移式滑坡的特點相吻合。中部變形最大為0.2m，考慮到中部土體較厚，主要為土體沉降產生的垂直位移，而沿滑動方向變形較小，說明此時滑坡尚處于蠕滑階段。對B—B′位置進行切片后，調出最大剪應力云圖，如圖8所示，可見滑體、滑帶內最大剪應力基本小于50kPa，尚未達到坡體的抗剪強度，故滑坡目前基本穩定。

圖7 蓄水前滑坡位移云圖

圖8 最大剪應力分布云圖

FLAC3D中通常以塑性區的貫通性作為滑坡是否失穩判據。各工況滑坡體的塑性區分布如圖9所示，由圖9可知，未蓄水時，滑坡僅Ⅱ區中后部側邊界有少量塑性變形，變形體長度約40m，寬度約2m，與實際中滑坡出現裂縫位置基本對應，此時滑坡基本穩定；蓄水位為2535m時，Ⅱ區側邊界塑性區范圍擴大，有逐漸連通趨勢，Ⅰ區邊界開始出現少量塑性變形；蓄水位為2571m時，滑坡后緣塑性區基本連通，I區坡腳等位置出現少量塑性變形，實際Ⅰ區內已經出現連通的塑性變形面，此時滑坡處于欠穩定—臨界穩定狀態；蓄水位為2605m時，滑坡邊界位置塑性區完全貫通，滑帶內出現剪應力塑性變形。此時對B—B′剖面位置進行切片處理，滑坡體內最大剪應變增量分布如圖10所示，可見滑坡剪應變增量主要集中在滑帶位置，前、后緣已完全貫通，滑坡沿滑帶整體剪切破壞，坡體已發生失穩。對蓄水前后坡體中部的位移進行監測，如圖11所示，蓄水至2605m時坡體最大位移已增加至12m，說明滑動后部分坡體已產生堆積。

圖9 塑性區分布云圖

表3 穩定系數計算結果

圖10 最大剪應變增量云圖

圖11 坡體中部位移監測曲線

不同蓄水位下滑坡的破壞狀態也反應了滑坡失穩的演化過程，具體為：滑體前緣為阻滑段，下滑力主要來自于滑坡中后部土體，隨著庫水位升高，前緣被水淹沒，地下水滲透入滑帶中，前緣滑帶力學性質降低，使阻滑作用降低。此時在后緣坡體巨大推力作用下，前緣逐漸被推動導致坡體整體破壞。其破壞過程可總結為：蠕滑、拉裂、剪切、滑動4個階段。

Flac3D內置了穩定系數計算方法，即強度折減法，計算中通過不斷降低抗剪強度折減系數Fs，將折減后的力學參數重新帶入模型計算，直到模型發生破壞，此時的抗剪強度折減系數Fs就是滑坡的安全系數F。該功能可通過調用FOS命令進行計算，計算結果見表3。與基于傳遞系數法的二維計算結果相比，三維數值模擬結果偏大。除了與內置算法不同外，這主要由于二維計算僅考慮了單條剖面位置的局部失穩，而三維計算更注重滑坡的整體性穩定性，故三維計算的結果與實際情況更相符合。

4 結論

(1)該近壩庫岸滑坡為重力主導型滑坡，水庫蓄水是其失穩主要誘發條件。在設計蓄水位為2535(低壩方案)～2605m(高壩方案)時，蓄水位高度與滑坡穩定性呈負相關，具體為：在未蓄水時，天然條件下滑坡處于基本穩定；在蓄水位為2535m時，滑坡穩定系數有所降低，但仍處于基本穩定；當蓄水位達到2571m時，滑坡處于臨界穩定；蓄水位達到高壩方案設計蓄水位2605m時，滑坡發生失穩破壞。

(2)三維數值模擬結果證明了該滑坡為上部土體擠壓下部產生變形的推移式滑坡，滑動發育過程主要包括了蠕滑、拉裂、剪切、滑動4個階段，蓄水后加快了滑坡變形發育的速度。

(3)三維計算滑坡穩定系數稍高于二維計算結果，這是由于三維數值模擬相比于二維計算更著重于滑坡的整體穩定性，其結果與實際情況也更相符。結合兩種分析方法，可對滑坡的成因模式和發育趨勢進行相對準確的預測。

(4)僅從對近壩庫岸的穩定性方面考慮，該工程目前所采用的低壩方案(蓄水位2535m)更有利，但蓄水后滑坡穩定性仍有降低，需做好長期監測工作。本文使用的分析方法也可推廣到類似工程中，可為水庫、流域梯級等規劃及后期防治工作提供參考。