基于Visual Modflow的關(guān)門山溝堰塞壩滲流穩(wěn)定模擬研究

李明忠,林 妮

(長江勘測規(guī)劃設(shè)計研究有限責(zé)任公司,湖北 武漢 430010)

0 引 言

滑坡、崩塌、泥石流堵塞天然河道,形成壅塞體,上游來水在庫內(nèi)蓄積,形成堰塞湖[1]。堰塞湖水位不斷上漲,將首先造成淹沒災(zāi)害,水位繼續(xù)上漲,滲流作用加強,可能出現(xiàn)管涌,進一步發(fā)展可能導(dǎo)致堰塞湖潰決。潰決后,庫水宣泄,下游洪水可相當(dāng)于千年一遇至萬年一遇[2],將毀壞沿途城鎮(zhèn)、耕地,造成大量生命財產(chǎn)損失,同時,潰壩洪水將可能引起新的滑坡、崩塌等次生災(zāi)害,進一步加劇損失[3]。

“5·12”汶川地震以來,國內(nèi)學(xué)者開始關(guān)注堰塞湖的滲流問題,徐文杰等研究了肖家橋堰塞湖滲透穩(wěn)定性,通過出逸流量及高程差,反算流量系數(shù),推求了肖家橋堰塞湖水位在743 m和761 m時的滲透坡降,二者均小于材料允許比降,在缺少地質(zhì)勘察的情況下,初步判斷肖家橋堰塞壩不會出現(xiàn)滲透破壞[4]。嚴祖文等[5]計算了唐家山堰塞湖非穩(wěn)定滲流和穩(wěn)定滲流兩種工況的滲流穩(wěn)定性,得出了唐家山堰塞壩滲透破壞可能性較小的結(jié)論。石振明等[6]從滲流的基本規(guī)律、數(shù)值方法、模型試驗等方面總結(jié)了堰塞壩組成材料滲流的機理和研究方法,指出土體滲透變形和滲透穩(wěn)定性由幾何條件和水力條件共同決定。胡卸文等[7]對唐家山堰塞湖滲流問題進行了數(shù)值模擬,得出結(jié)論:下游壩腳處某些部位會發(fā)生零散滲透破壞,但對于堰塞壩的整體穩(wěn)定性影響不大。徐軼等[8]對堰塞湖應(yīng)急處置工程措施及典型案例進行分析,較全面地總結(jié)梳理了堰塞湖應(yīng)急處置工程措施的經(jīng)驗及技術(shù)發(fā)展趨勢。

本文基于關(guān)門山溝堰塞壩實測土層物理參數(shù),采用Visual Modflow可視化三維地下水模型軟件模擬滲流場,分析滲透坡降變化規(guī)律,判斷在各模擬水位條件下,堰塞壩是否發(fā)生滲透破壞及滲透破壞形式。根據(jù)滲透坡降變化規(guī)律,提出發(fā)生滲透破壞時堰塞湖臨界水位,對堰塞湖滲流穩(wěn)定性研究具有理論和實際參考價值。

1 關(guān)門山溝堰塞湖概況

關(guān)門山溝堰塞湖位于都江堰市白沙河上游關(guān)門山溝,“5·12”汶川特大地震時右岸百余米高的山體在強震作用下形成滑坡,滑坡體快速沖擊左岸基巖,淤積河道形成堰塞湖。堰塞壩下游河底高程1 648.0 m,壩頂高程1 728.0 m,上游水位1 705.0 m,下游壩腳水位1 649.0 m。堰塞壩高80 m,寬230 m,順河長約500 m,不計算右岸滑坡體部分,估算土石方量2.7×107m3,蓄水量約370萬m3。堰塞壩主要由碎石土組成,結(jié)構(gòu)松散,滑坡陡壁高150 m,坡度50°～60°,右岸邊坡尚存大量物源,有再次下滑的可能。堰塞壩以上主河道長約11 km,河道坡降13.5%,堰塞湖以上集雨面積約56 km2。根據(jù)SL 450-2009《堰塞湖風(fēng)險等級劃分標(biāo)準》,關(guān)門山溝堰塞湖風(fēng)險等級為Ⅱ級。

關(guān)門山溝堰塞湖所在區(qū)域海拔高程在1 200～3 750 m之間,地貌特征主要為構(gòu)造剝蝕型,其次為侵蝕堆積地貌。堰塞湖處于揚子地臺與松潘-甘孜地槽區(qū)之間的構(gòu)造過渡帶上,位于龍門山構(gòu)造帶北段,龍門山主中央斷裂帶附近。根據(jù)四川防震減災(zāi)信息網(wǎng)公布的地震烈度分布圖,“5·12”汶川地震對工程場地的影響烈度為Ⅹ度。堰塞壩主要由碎石土構(gòu)成,通過鉆孔及物理指標(biāo)試驗,獲取碎石土物理參數(shù)見表1。

表1 堰塞壩物理參數(shù)

對于原河道覆蓋層,由于沒有進行物理指標(biāo)試驗,參考白沙河上游地區(qū)相關(guān)文獻確定其參數(shù)。根據(jù)太沙基公式,綜合考慮計算成果及實驗測量,確定關(guān)門山溝堰塞壩允許坡降取值0.40～0.55,河道覆蓋層取值0.47～0.60。

2 三維滲流模型

Visual Modflow由加拿大Waterloo水文地質(zhì)公司開發(fā)研制,基于有限差分法理論進行計算,作為三維地下水水流和溶質(zhì)數(shù)值模擬評價的標(biāo)準可視化專業(yè)軟件系統(tǒng),運用廣泛[9]?；诔Ｃ芏鹊牡叵滤S流動基本方程:

(1)

式中:Kxx,Kyy,Kzz分別為滲透系數(shù)在x,y,z方向的分量;H為作用水頭;w為單位時間內(nèi)單位體積流進或流出的水的體積通量;Ss為儲水系數(shù),表示含水層地下水頭變化時,由于含水層垂向壓縮和地下水彈性膨脹,單位體積含水層變化的水的體積;t為時間。

(1) 第一類邊界條件。在邊界上,各時刻每一點水頭是給定的:

(2)

式中:H(x,y,z,t)為邊界S1上點(x,y,z)在t時刻的水頭,f1(x,y,z,t)為S1上的已知函數(shù)。

(2) 第二類邊界條件。已知某一部分邊界單位面積上流入(流出)的流量,表示為

(3)

式中:n為邊界S2的外法相方向,q1為邊界z單位面積的補給量,為已知函數(shù)。

Modflow采用迭代法求解滲流差分方程,在計算之前,需要給定一個初始水位值,并求解下一個時間段的水位值;所求得的水位值作為下一步初始水位,再次迭代,直到兩次迭代結(jié)果相差足夠小為止[10]。

采用Visul Modflow,建立關(guān)門山溝堰塞湖三維滲流模型。取兩岸山體及河床下伏基巖為相對不透水層,上下游作為自由流動面,根據(jù)滲透系數(shù)對堰塞壩進行分區(qū),分區(qū)內(nèi)考慮土體性質(zhì)均一,在滲流模擬過程中,土體物理力學(xué)性質(zhì)不發(fā)生變化,假定堰塞湖上游水位和下游出水口處水位不變,并作為定水頭邊界處理。

2.1 網(wǎng)格劃分及參數(shù)給定

本次模擬選用穩(wěn)定滲流,堰塞湖順河長500 m,寬230 m,堰塞壩三維模型見圖1。模型平面上劃分為200×100的網(wǎng)格,將堰塞壩概化為均質(zhì),縱向劃分為堰塞壩及原河道覆蓋層兩層[11]。第一層滲透系數(shù)賦值為6.5×10-5m/s,第二層滲透系數(shù)賦值為1.3×10-5m/s。第一層貯水率為1×10-4,重力給水度、有效孔隙度、總孔隙度分別為0.2,0.18和0.36;第二層貯水率取為0.5×10-5,重力給水度、有效孔隙度、總孔隙度分別取為0.03,0.085和0.12。

圖1 堰塞壩三維模型Fig.1 Three-dimensional model of barrier dam

2.2 邊界條件給定

堰塞壩地下水運動形式主要是上游湖水補給,下游出口排泄。假設(shè)在模擬時間段內(nèi),水位變化微小,上下游均采用定水頭邊界[12],堰塞壩與兩岸、下伏基巖接觸部位不設(shè)定邊界條件,系統(tǒng)默認其為隔水邊界,順河方向上,上下游壩坡作為自由臨空面[13]。上游水位分別取1 705,1 710,1 715 m和1 720 m,下游水位取1 649 m。

3 滲流模擬分析

對上述4種水位條件下堰塞湖滲流場進行模擬分析,見圖2～5。

圖2 水位1 705 m滲流場Fig.2 Seepage field diagram at 1 705 m

圖3 水位1 710 m滲流場Fig.3 Seepage field diagram at 1 710 m

圖4 水位1 715 m滲流場Fig.4 Seepage field diagram at 1 715 m

圖5 水位1 720 m滲流場Fig.5 Seepage field diagram at 1 720 m

3.1 滲流特性分析

對4種水位條件下堰塞壩滲流場進行分析,當(dāng)水位為1 705 m時,堰塞壩浸潤線較低,滲水主要由堰塞壩與覆蓋層分界處溢出,與現(xiàn)場踏勘情況相符。壩軸線以上,滲透坡降較小且變化也較小,壩軸線以下,滲透坡降變化加大,在壩腳附近達到最大。當(dāng)水位上升到1 710 m時,堰塞壩浸潤線提升,溢出點向上游移動,最大滲透坡降增大。隨著水位繼續(xù)上升、浸潤線將進一步升高,溢出點進一步向上游移動,最大滲透坡降持續(xù)增大。

3.2 滲透破壞臨界水位確定

堰塞壩發(fā)生滲透破壞的條件是滲透坡降達到允許坡降。堰塞壩在不同水位條件下出口處滲透坡降見表2。從表2可知,滲透坡降與水位呈正相關(guān)性,最大滲透坡降出現(xiàn)在下游出水口附近,堰塞壩平均坡降及最大坡降均大于覆蓋層,且堰塞壩層的滲透坡降變化速度大于覆蓋層。

表2 堰塞壩出口處滲透坡降

當(dāng)堰塞湖水位在1 705 m和1 710 m時,兩土層的最大滲透坡降均小于允許坡降,堰塞壩處于穩(wěn)定狀態(tài)。當(dāng)堰塞湖水位達到1 715 m時,原河道覆蓋層的滲透坡降仍小于允許坡降,但堰塞壩層最大滲透坡降大于允許坡降下限,一些較細的顆粒將被帶走,堰塞壩局部出現(xiàn)滲透破壞[14]。當(dāng)水位達到1 720 m時,堰塞壩層最大滲透坡降超過允許坡降上限,覆蓋層滲透坡降達到允許滲透坡降上限,在壩腳處將出現(xiàn)較大范圍的滲透破壞。如果水位進一步上升,將引起更大范圍的滲透破壞,甚至形成貫穿性通道,危及堰塞壩整體安全。

最大滲透坡降與允許坡降對比見圖6,根據(jù)滲透趨勢,隨著堰塞湖水位的上升,兩土層的滲透坡降均逐漸增大。在一定范圍內(nèi),水位上升主要影響堰塞壩層,當(dāng)堰塞湖水位達1 712 m時,堰塞壩層最大滲透坡降達到其允許滲透坡降下限,堰塞壩局部細顆粒將被水流帶走;當(dāng)水位達到1 714 m時,覆蓋層滲透坡降達到其允許滲透坡降下限,覆蓋層局部細顆粒也將移動,水位繼續(xù)上升;當(dāng)水位達到1 719 m時,堰塞壩層最大滲透坡降達到其允許滲透坡降上限,堰塞壩層將出現(xiàn)更大范圍的滲透破壞;當(dāng)水位達到1 722 m時,覆蓋層滲透坡降達到其允許滲透坡降上限,覆蓋層也將出現(xiàn)更大范圍的滲透破壞,嚴重影響堰塞壩安全。因此,將1 719 m確定為關(guān)門山溝堰塞湖滲透破壞臨界水位。

圖6 最大滲透坡降與允許坡降對比Fig.6 Comparative analysis of the maximum seepage gradient and the allowed gradient

隨著堰塞湖水位的抬升,關(guān)門山溝堰塞壩表層碎石土發(fā)生滲透變形,穩(wěn)定性降低,但壩體整體穩(wěn)定。預(yù)測堰塞壩潰決模式為:下游側(cè)碎石土層因滲透發(fā)生破壞,出水口甚至發(fā)生局部坍塌,但由于大塊石骨架的存在,由滲流引起整體潰決可能性不大。隨著堰塞湖蓄滿,漫頂導(dǎo)致碎石土被侵蝕、淘刷,水流速度加大,帶動進一步?jīng)_刷下切,形成以漫頂沖刷為主的大范圍潰決。

4 結(jié) 論

(1) 本文采用Visual Modflow三維地下水可視化軟件,建立關(guān)門山溝堰塞壩三維模型。將滲流區(qū)分成堰塞壩層和原河道覆蓋層,取兩岸山體及下伏基巖作為相對不透水層,假設(shè)堰塞湖水位和下游出水口水位不變并作為定水頭邊界處理。模擬了堰塞湖在1 705,1 710,1 715 m和1 720 m四種水位條件下的滲流場。隨著水位的上升,堰塞壩滲透坡降逐漸增大,且最大滲透坡降出現(xiàn)在下游出水口附近。

(2) 基于滲流模擬成果,得出4種水位條件下出水口平均滲透坡降和最大滲透坡降,并與筑壩材料允許坡降進行比對。通過滲透坡降變化趨勢,確定1 719 m為堰塞壩滲透破壞臨界水位。