李明忠,林 妮

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

0 引 言

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

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

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

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

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

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

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

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

2 三維滲流模型

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

本次模擬選用穩(wěn)定滲流,堰塞湖順河長(zhǎng)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 邊界條件給定

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

3 滲流模擬分析

對(duì)上述4種水位條件下堰塞湖滲流場(chǎng)進(jìn)行模擬分析,見圖2～5。

圖2 水位1 705 m滲流場(chǎng)Fig.2 Seepage field diagram at 1 705 m

圖3 水位1 710 m滲流場(chǎng)Fig.3 Seepage field diagram at 1 710 m

圖4 水位1 715 m滲流場(chǎng)Fig.4 Seepage field diagram at 1 715 m

圖5 水位1 720 m滲流場(chǎng)Fig.5 Seepage field diagram at 1 720 m

3.1 滲流特性分析

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

3.2 滲透破壞臨界水位確定

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

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

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

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

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

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

4 結(jié) 論

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

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