地下水滲流場對庫區(qū)滑坡穩(wěn)定性影響的數(shù)值模擬
——以白馬庫區(qū)羊角灘滑坡為例

田 凱， 姚品品， 鐵永波， 徐 偉

(1.中國地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，四川 成都 610081; 2.中國地質(zhì)大學(武漢)工程學院工程學院，湖北 武漢 430070)

0 引言

降雨入滲與水庫蓄水引起庫水位的抬升及水庫調(diào)度，導致庫水位大幅度升降，使水庫沿岸岸坡地下水滲流場發(fā)生重大改變，是水庫滑坡等地質(zhì)災害發(fā)生的主要誘因[1]。降雨入滲導致水庫岸坡地下水滲流場變化的數(shù)值模擬問題本身就是水文地質(zhì)工程地質(zhì)的難題，而水庫蓄水引起庫水位的抬升及水庫調(diào)度導致庫水位大幅度升降使這一難題更加復雜[2-3]。因此，模擬降雨入滲與水庫蓄水引起庫水位的抬升及水庫調(diào)度導致庫水位大幅度升降聯(lián)合作用下的滑坡地下水滲流場，成為當前水文地質(zhì)工程地質(zhì)界研究的一個熱點內(nèi)容[4]。

地下水滲流場數(shù)值模擬的主要任務是求解水頭函數(shù)，確定滲流自由面和滲流量等滲流狀態(tài)，其方法主要有有限單元法、邊界元法和有限差分法[5-6]。用有限元方法分析地下水滲流問題時，若計算涉及到滲流自由面以及溢出面情況，由于滲流自由面(浸潤線)所處的位置事先不確定，滲流流域范圍便具有不確定性，這將導致邊值問題中具有未定的邊界，使得求解地下水滲流場問題由原本一個簡單的線性問題變成事先部分邊界條件不確定的復雜邊界線性問題[7-9]。本文主要方法是通過迭代運算不斷修改滲流自由面的位置，直至計算結果滿足滲流自由面上的邊界條件，獲得足夠接近實際情況的解答，最終確定滲流自由面位置，計算工作量很大是這種通過多次迭代運算最終確定自由面位置方法的明顯缺陷[10]。

本文以重慶烏江白馬庫區(qū)上游羊角灘滑坡為例，結合區(qū)域水文地質(zhì)工程地質(zhì)資料[8]，通過調(diào)查試驗得到模擬計算所需參數(shù)，從飽和/非飽和滲流的基本理論出發(fā)，討論了滲流模擬有限元法的基本原理，給出了降雨入滲與水庫水位升降聯(lián)合作用下的滲流控制方程及有限元求解方法，利用Geostudio軟件SEEP模塊建立水文地質(zhì)結構概念模型進行模擬，并對模擬結果進行分析，為今后類似工程實際應用提供一定的參考依據(jù)。

1 研究區(qū)概況

羊角灘滑坡位于重慶市武隆區(qū)白馬鎮(zhèn)烏江白馬庫區(qū)上游，庫區(qū)白馬航電樞紐工程已于2019年正式開工建設，滑坡位于該樞紐工程上游6.7 km，其穩(wěn)定性直接影響著白馬航電樞紐工程的修建(圖1)。因此，本文選取羊角灘滑坡進行降雨和庫水位變化條件下的地下水滲流場模擬分析，對滑坡穩(wěn)定性進行定量評價。

圖1 羊角灘滑坡區(qū)域位置

羊角灘滑坡平面呈長舌形，主滑方向NE14°，后緣高程367 m，前緣高程159 m，縱向長1 066 m，平均寬372 m(后緣最小寬度51 m，前緣最大寬度541 m)，面積33.9×104m2，平均厚度19.9 m，體積672×104m3(圖2)。據(jù)該滑坡的長期水文監(jiān)測資料，庫水位變化對滑坡后緣影響較小，因此選取滑坡距烏江0～725 m的范圍進行模擬。圖3為滑坡地質(zhì)剖面。

圖2 羊角灘滑坡地質(zhì)簡圖

通過已有勘查資料進行斜坡地質(zhì)分層： ①分布于滑坡上部的崩坡積層(Qhdl+col)，主要成分為灰?guī)r碎塊石土，最大厚度約41 m； ②分布于滑坡體內(nèi)部的滑坡堆積層(Qhdel)，主要為黏土夾少量碎石，最大厚度29.5 m； ③滑坡前緣分布于滑坡剪出口外側河漫灘，為河流沖積堆積(Qhal+dl)，成分以粉細砂為主，細—中粒結構，較松散； ④發(fā)育志留系下統(tǒng)龍馬溪組(S1l)淺黃、灰綠色砂質(zhì)頁巖。YJ01為鉆孔，YJ03 、YJ06、YJ09為水文長期觀測孔。

圖3 滑坡地質(zhì)剖面

2 滑坡計算模型的建立

本滑坡的剖分是根據(jù)具體邊界條件，在滿足一定精度要求的情況上進行的細致剖分，其中滑體的剖分較詳細，而滑床的剖分較稀疏?；陆＜坝邢拊W(wǎng)格剖分見圖4,斜坡地質(zhì)剖面分為4層。

(1)滑體表部中等滲透帶。位于滑坡中部的坡體內(nèi)部，為滑體表部強滲透帶內(nèi)的透鏡體，下接弱滲透的黏土層，巖土體滲透性一般。

圖4 滑坡模型的材料分區(qū)

(2)滑體下部弱滲透帶。為后期形成的羊角灘滑坡的次級滑動帶，位于滑坡的前緣，分布范圍為海拔175～255 m，以類似透鏡體的形式位于滑體下部。弱滲透帶北端為滑坡前緣的剪出口，南端與滑體表部中等滲透帶相連接，下伏基巖為滑床。巖土體的滲透性較差。

(3)滑坡前緣的強滲透帶。分布于滑坡剪出口外側河漫灘，為河流沖積堆積層，成分以粉細砂為主，層薄，下伏基巖，斜坡地下水主要通過該層向烏江進行排泄，地下水埋深相對較淺。

(4)基巖(不透水層)。表1為各個地層滲流模擬巖土體的飽和滲透系數(shù)。

表1 滲流模擬巖土體飽和滲透系數(shù)取值

2.1 初始條件和邊界條件

在計算降雨前滑坡體的地下水位線時，首先需要根據(jù)勘察資料確定出初始水位。當烏江水位為161 m時，滑坡后部地下水位線傾角為9°，由此可以反算出滑坡左邊界的水頭邊界為249 m。因此，在降雨前可以設定初始條件： 滑坡體左邊界為249 m水頭邊界，滑坡右側為烏江水位。然后，再模擬計算滑坡體降雨入滲期間引起的瞬態(tài)滲流及地下水位變化，通過在滑坡表面SEEP中定義邊界函數(shù)來模擬降雨過程。

(1)入滲邊界。假設降雨區(qū)域覆蓋滑坡體表面，通過函數(shù)定義為單元流量邊界。通過設計單位節(jié)點流量—時間的函數(shù)來確定降雨的大小和時間，通過定義某一時刻的單元流量來實現(xiàn)降雨強度的變化，動態(tài)模擬降雨入滲的過程。

(2)模型兩側。在滑坡右側按定水頭邊界條件確定水位高程，左側也按定水頭邊界條件水位高程。

(3)模型底面。假設滑床基巖的滲透性很小，滑床底部的邊界條件可以設為不透水的邊界條件。

2.2 模型參數(shù)的選取

根據(jù)Geostudio軟件SEEP/W自帶的函數(shù)，用土水特征曲線來研究非飽和土的含水率和土體基質(zhì)吸力的關系，依據(jù)軟件提供的土的滲透性函數(shù)來定義非飽和土的滲透性。在各種工況下，所研究滑坡共有3種巖土類型，所以選擇碎塊石土、黏土、粉細砂3種土的特性函數(shù)。本文選用Fredlund&Xing方法函數(shù)。根據(jù)現(xiàn)場地表注水試驗和鉆孔注水實驗結果，可以得到模擬區(qū)域3種材料的土-水特征曲線及滲透函數(shù)曲線(圖5)。

(a) 滑體碎塊石土

(b) 滑體下部黏土

(c) 滑體前緣粉細砂

3 數(shù)值計算結果分析

3.1 模擬工況

選取YJ1、YJ3 、YJ6、YJ9為觀測孔，與正常水位工況模擬計算結果作對比，然后逐步調(diào)整各滲透帶的水文地質(zhì)參數(shù)，直至滿足模擬精度為止。據(jù)重慶武隆1951～2018年實測降水量資料統(tǒng)計，歷年最大日降水量為157.54 mm(2003年6月24日)，以此作為特殊暴雨工況。根據(jù)水庫實際營運條件，滑坡滲流計算采取下列6種計算工況。

(1)工況1。當前天然161 m河水位，該工況下前緣庫水位為161 m，水位線下邊界為固定水頭邊界，水位線以下為零水頭邊界，后緣為249 m水頭邊界，基巖面為零流量邊界。

(2)工況2。暴雨、久雨+161 m河水工況(以161 m河水位為初始，降雨4 d)，該工況下前緣庫水位為161 m，水位線下邊界為固定水頭邊界，水位線以下為零水頭邊界，后緣為249 m水頭邊界，基巖面為零流量邊界，上表面庫水位線以上為降雨邊界。

(3)工況3。183 m正常庫水位工況，該工況下前緣庫水位為183 m，水位線下邊界為固定水頭邊界，水位線以下為零水頭邊界，后緣為249 m水頭邊界，基巖面為零流量邊界。

(4)工況4。暴雨、久雨+193 m洪水位工況(以193 m河水位為初始，降雨4 d)，該工況下前緣庫水位為193 m，水位線下邊界為固定水頭邊界，水位線以下為零水頭邊界，后緣為249 m水頭邊界，基巖面為零流量邊界，上表面水位線以上為降雨邊界。

(5)工況5。182 m水位升至193 m洪水位工況(以182 m河水位為初始，4 d后升至193 m)，該工況為暴雨情況下，庫水位在4 d時間內(nèi)從183 m洪水位降至182 m水庫正常蓄水位。

(6)工況6。193 m洪水位驟降至183 m水位工況(以193 m河水位為初始，4 d后降至183 m)，該工況為庫水位在4 d時間內(nèi)從193 m洪水位降至183 m水庫正常蓄水位。

本文主要研究的是降雨和庫水位這2個因素對滑坡地下水滲流的影響。通過對模擬結果的分析可知，庫水位變化主要影響滑坡前緣，庫水位的高低控制著地下水位的高低，庫水位的升降控制著地下水位的升降，且地下水位是隨庫水位變化而迅速發(fā)生變化的； 降雨對整個滑坡的水位變化都有一定的影響，而且滑坡中部地下水位變化影響趨于最大，前后部影響趨于較小。

滑坡地下水滲流場的變化會引起其穩(wěn)定性系數(shù)的變化，通過利用地下水滲流場模擬的結果得出最危險滑動面，用畢肖普法對羊角灘滑坡穩(wěn)定性進行了簡單計算，并分析滲流場變化對滑坡穩(wěn)定性的影響。物理力學參數(shù)選取見表2。

表2 滑坡巖土物理力學參數(shù)綜合建議值

3.2 模擬結果分析

采用畢曉普法進行計算，只考慮地下水位的不同對滑坡穩(wěn)定性的影響。經(jīng)計算，穩(wěn)定性系數(shù)計算結果見表3和圖6。

表3 各工況穩(wěn)定性計算結果

(a) 工況1(穩(wěn)定性系數(shù)=1.410) (b) 工況2(穩(wěn)定性系數(shù)=1.222)

(c) 工況3(穩(wěn)定性系數(shù)=1.274) (d) 工況4(穩(wěn)定性系數(shù)=1.225)

(e) 工況5(穩(wěn)定性系數(shù)=1.275) (f) 工況6(穩(wěn)定性系數(shù)=1.158)

圖6 不同工況條件下最危險滑面及穩(wěn)定性系數(shù)

通過對表3和圖6的分析，可以得出以下結論。

(1)不同庫水位對滑坡穩(wěn)定性的影響不同，隨著庫水位的上漲，滑坡穩(wěn)定性系數(shù)呈現(xiàn)先減小后增大的趨勢。分析其原因，是因為滑坡穩(wěn)定性的變化主要受到浮托作用、壓坡作用和滲流力作用的影響。該實例中，庫水位的變化主要對滑坡前緣影響較大，前緣為阻滑段，當上部土層重力越大時，抗滑力與下滑力的差值越大，阻滑作用越強，滑坡越穩(wěn)定。庫水位變化對滑坡穩(wěn)定性的影響正是通過這種方式進行的。初始階段，隨水庫庫水位的上升，水面以下土層比例將增大，水的浮托作用同比增大，在計算公式中體現(xiàn)為滑面反力減小，阻滑作用力降低，滑坡穩(wěn)定性趨于減??； 當達到一個特定值后，如果庫水位繼續(xù)增加，庫水便會起到一定的壓坡作用，使阻滑作用得到增強，滑坡的穩(wěn)定性將會增加。同時，水庫庫水位的上升還會減小滑坡中的水力的梯度值，使地下水滲流力趨于減小，有利于滑坡的穩(wěn)定性。

(2)在降雨條件下，滑坡穩(wěn)定性隨著降雨歷時的增加而趨于變小，到達一定階段后，滑坡穩(wěn)定性系數(shù)不再變化。本文對161 m水位(工況2)和193 m水位(工況4)進行了歷時4 d的穩(wěn)定性計算。由表3(工況2和工況4)結果可以看出，初始階段穩(wěn)定性系數(shù)減小幅度不大，隨著降雨時間的增加，穩(wěn)定系數(shù)將迅速減小，之后慢慢地趨于穩(wěn)定。分析其原因是： 降雨入滲補給地下水需要一定的時間，當降雨入滲補給地下水時，水面上升，地下水滲流速度加大，滑坡的穩(wěn)定性將迅速降低，之后地下水滲流情況趨于穩(wěn)定，穩(wěn)定性系數(shù)也就不再發(fā)生變化。

(3)庫水位從183 m上升到193 m的過程中，隨著庫水位的上漲，滑坡的穩(wěn)定性有所增加。分析其原因是： 183 m水位時，滑坡前緣阻滑段已經(jīng)被庫水淹沒，而隨著水位上漲，浮托作用將不再增大，而壓坡作用卻不斷增大，同比水位的抬升，滑坡內(nèi)水力梯度趨于減小，滲流力趨于減小，接近庫岸處的部位甚至出現(xiàn)水體的倒流； 這些因素都會使得滑坡穩(wěn)定性有所增大。

(4)庫水位從193 m下降至183 m的過程中，滑坡的穩(wěn)定性明顯趨于減小。分析其原因是： 水庫水位下降時，水庫水壓坡作用趨于減小，同時滑體內(nèi)水力梯度增大，滲流力逐漸增大； 這些因素都將導致滑坡的穩(wěn)定性變差，因此滑坡穩(wěn)定性迅速趨于減小。

4 結論

本文對6種工況進行了數(shù)值模擬，并在模擬基礎上對滑坡的穩(wěn)定性進行了分析，得到的主要結論如下：

(1)該滑坡最危險潛在滑動面位置受地形地貌的影響較大，各種工況下，范圍和位置變化相對較小。

(2)不同的水位情況下，隨著降雨時間的增加，滑坡穩(wěn)定性逐漸趨于減小，最后趨于一個定值; 降雨入滲需要一定量的時間，因此最開始滑坡穩(wěn)定性減小相對較小。

(3)隨著庫區(qū)水位的升高，滑坡穩(wěn)定性趨于減小，當水位到達一定值的時候，滑坡穩(wěn)定性將隨著水庫的水位升高而趨于增大。

(4)地下水位的驟升對滑坡穩(wěn)定性的影響要視具體情況而定。地下水位的驟升對浮托作用并沒有太大影響，但對壓坡作用影響比較大，同時降低了地下水滲流力的影響，因此，地下水位驟升不僅沒有使滑坡失穩(wěn)的可能性進一步增大，反而對滑坡穩(wěn)定起到了積極的作用。

(5)當?shù)叵滤E降時，滑坡穩(wěn)定性系數(shù)趨于最小。地下水位的驟降增大了其滲流力，降低了前緣庫水位的壓坡作用，但并未減弱地下水本身的浮托作用，因此在此工況下，滑坡失穩(wěn)的可能性趨于最大。