陳 祖 煜，張 強(qiáng)，侯 精 明，王 琳，馬 利 平

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710048; 2.中國(guó)水利水電科學(xué)研究院 巖土工程研究所, 北京 100048)

2018年10月10日22：06，西藏自治區(qū)昌都市江達(dá)縣波羅鄉(xiāng)白格村和四川省甘孜藏族自治州白玉縣交界處金沙江河道右岸發(fā)生大規(guī)模山體滑坡，堵塞金沙江干流河道，形成白格堰塞湖。12日17：15堰塞湖開(kāi)始漫溢自然過(guò)流，13日14：30基本退至基流，22：00險(xiǎn)情解除。此次泄流，西藏共轉(zhuǎn)移安置群眾13 637人，云南共轉(zhuǎn)移群眾11 550人，無(wú)人員傷亡，下游在建的葉巴灘、蘇洼龍水電站導(dǎo)流洞過(guò)流，造成一定損失，但影響可控。

堰塞湖潰決洪水過(guò)程線通常是通過(guò)量測(cè)水位降落的過(guò)程獲得。但是，由于“10·10”堰塞湖在兩天后即發(fā)生了漫頂潰決，相應(yīng)的水位量測(cè)設(shè)備未能準(zhǔn)確捕捉到本次堰塞湖潰決的水位下降過(guò)程線，因而無(wú)法通過(guò)庫(kù)容水位關(guān)系曲線計(jì)算“10·10”白格堰塞湖潰決洪水過(guò)程線。幸運(yùn)的是，在堰塞湖下游54 km葉巴灘水文站和137 km的拉哇水文站均獲得了實(shí)測(cè)洪水過(guò)程線。因此，應(yīng)用洪水波演進(jìn)的原理和相關(guān)的程序可以對(duì)壩址(0 km)處的潰壩洪水過(guò)程線進(jìn)行反演。本文應(yīng)用中國(guó)水利水電科學(xué)研究院潰壩洪水分析軟件DB-IWHR和西安理工大學(xué)侯精明研發(fā)的GST洪水演進(jìn)模型，通過(guò)參數(shù)敏感性分析反演獲得堰塞湖的潰決洪水。

1 基礎(chǔ)資料

“10·10”白格滑坡體積約3 500萬(wàn)m3，滑入河道體積約2 500萬(wàn)m3，堰塞體順河長(zhǎng)約2 km，寬約450～700 m，堰塞體整體呈左岸高右岸低之勢(shì)，左岸最高點(diǎn)高程為3 005.0 m，右岸埡口高程2 931.4 m，堰塞體高度約61～100 m。堰塞湖從12日17：15開(kāi)始漫溢自然過(guò)流，13日00:45堰塞湖壩上水位達(dá)到最高值2 932.69 m，對(duì)應(yīng)蓄水量2.9億m3。13日01:00過(guò)流明顯增加，06：00過(guò)流流量達(dá)到最大值約10 000 m3/s，此后開(kāi)始消退，13日14：30基本退至基流[1-2]。

在堰塞湖發(fā)生后第一時(shí)間進(jìn)行了無(wú)人機(jī)拍攝，三維點(diǎn)云圖如圖1所示。主堰塞體長(zhǎng)度為961.89 m，在潰壩發(fā)生后經(jīng)過(guò)水流沖刷，形成長(zhǎng)1 622 m、底寬80～120 m的泄流槽，泄流槽口高程約2 888 m、出口底板高程2 872 m，對(duì)應(yīng)堰塞湖蓄水量約0.5億m3。堰塞湖庫(kù)容水位關(guān)系曲線如圖2所示。

圖1 “10·10”白格堰塞體Fig.1 “10·10” Baige Landslide Dam

圖4 潰口的側(cè)向崩塌過(guò)程等效簡(jiǎn)化Fig.4 Equivalent simplifcation of lateral enlargement process

圖2 白格堰塞湖庫(kù)容-水位關(guān)系曲線Fig.2 Relation of reservoir storage capacity and water level of Baige barrier lake

2 潰壩洪水水力學(xué)基礎(chǔ)

2.1 潰口洪水分析

DB-IWHR潰壩洪水分析模型主要包括以下3個(gè)部分[3]。

2.1.1 潰口水力學(xué)條件

潰口斷面處的流量計(jì)算采用寬頂堰公式，如式(1)～(2)所示。

(1)

(2)

式中，C為綜合流量系數(shù)，理論值為1.7 m1/2/s[4]，建議取值1.3～1.7[5]；mb,mq分別為寬頂堰的流量系數(shù)和側(cè)收縮系數(shù)；H為庫(kù)水位；z為渠底高程。

水流經(jīng)過(guò)堰口跌落后水深為h，如圖3所示。流量系數(shù)m計(jì)算方法如式(3)所示，建議取值為0.8～0.6之間。

(3)

式中，h為潰口水深。

圖3 寬頂堰泄流示意Fig.3 Discharge of broad crest weir

2.1.2 沖刷模型

DB-IWHR采用雙曲線模型，其形式如式(4)。

(4)

式中，γ為扣除臨界剪應(yīng)力后的剪應(yīng)力。

γ=k(τ-τc)

(5)

2.1.3 潰口側(cè)向擴(kuò)展模型

潰口橫向擴(kuò)展采用了巖土工程中的圓弧滑動(dòng)面分析方法，并且比較合理地模擬了巖坡崩塌過(guò)程。巖坡崩塌過(guò)程如圖4(a)所示。但是在實(shí)際操作時(shí)，需輸入每一步的圓心坐標(biāo)、半徑和渠底坐標(biāo)，尋找臨界滑裂面，過(guò)程太過(guò)繁瑣。在實(shí)際應(yīng)用時(shí)，簡(jiǎn)化為一系列梯形[3]，如圖4(b)所示。

根據(jù)近期研究分析，我們提出了用雙曲線模型來(lái)計(jì)算梯形側(cè)面傾角增量Δβ。

(6)

(7)

式中，1/m1和1/m2分別表示雙曲線的初始切線和漸近線。文獻(xiàn)[6]詳細(xì)介紹了此雙曲線模型。

2.2 洪水波演進(jìn)模型

2.2.1 控制方程

洪水演進(jìn)采用GAST模型，它運(yùn)用平面二維淺水方程(SWES)模擬洪水演進(jìn)過(guò)程，控制方程如下[7]：

(9)

式中，q為變量矢量，包括水深h，x、y方向的單寬流量qx、qy；g為重力加速度，u和v分別為x、y方向的流速；F、G為x、y方向的通量矢量；S為源項(xiàng)矢量；zb為河床底面高程，Cf=gn2/h1/3為謝才系數(shù)，n為曼寧系數(shù)。

2.2.2 數(shù)值方法

GAST模型相較于經(jīng)典的水動(dòng)力學(xué)數(shù)值模型，其數(shù)值求解格式為Godunov類(lèi)型的有限體積法，它可以很穩(wěn)健地解決不連續(xù)問(wèn)題，并嚴(yán)格保持物質(zhì)守恒[7]。網(wǎng)格中的水、泥沙和動(dòng)量的數(shù)值通量應(yīng)用HLLC近似黎曼求解器計(jì)算[8]。干濕邊界處負(fù)水深問(wèn)題通過(guò)靜水重構(gòu)法解決。底坡源項(xiàng)和摩阻源項(xiàng)分別使用底坡通量法和半隱式法處理[9]。時(shí)間積分采用二階顯示Runge Kutta方法進(jìn)行，構(gòu)造具有二階時(shí)空精度的MUSCL型格式，能有效解決復(fù)雜地形引起的計(jì)算發(fā)散和物質(zhì)動(dòng)量的不守恒問(wèn)題，提高了計(jì)算的穩(wěn)定性[10]。并且，GST模型采用了GPU加速計(jì)算技術(shù)，大大提高了計(jì)算效率[11]。

3 潰壩洪水反演分析

金沙江電站分布如圖5所示。根據(jù)實(shí)測(cè)和團(tuán)隊(duì)前期研究相關(guān)數(shù)據(jù)，“10·10”白格堰塞體潰壩計(jì)算參數(shù)如表1所示。以往對(duì)唐家山、易貢、小崗劍等堰塞湖反演工作表明[12-14]，沖刷侵蝕系數(shù)b是其中敏感性較高的參數(shù)。其它參數(shù)對(duì)洪峰計(jì)算影響僅在5%左右。因此，我們?nèi)為0.000 4,0.000 5和0.000 6進(jìn)行潰壩洪水計(jì)算，結(jié)果獲得的潰口處的流量過(guò)程線如圖6所示。

圖5 金沙江部分電站分布Fig.5 Distribution of some power stations on Jinsha river

名稱符號(hào)數(shù)值地理學(xué)參數(shù)庫(kù)容系數(shù)p1,p2,p3,Hr0.11, 4.77, 90.20, 2910m初始庫(kù)水位H02932.50m初始底高程z02929m入流流量q800m3/s庫(kù)容水位關(guān)系h2910.00,2915.00,2935.00 mW90.20,116.80,278.40億m3水力學(xué)參數(shù)寬頂堰系數(shù)C1.43跌落系數(shù)m0.8啟動(dòng)流速Vc4.0侵蝕參數(shù)a, b1.1000,0.0006巖土力學(xué)參數(shù)壩體材料參數(shù)c,φ,γ50kPa, 45°,18.5kN/m3潰口雙曲線模型參數(shù)m1, m20.95,0.04初始潰口側(cè)向傾角 β0129°

注：其中庫(kù)容與水位關(guān)系用公式表示W(wǎng)=p1(H-Hr)2+p2(H-Hr)+p3

由圖6可以看出，當(dāng)b為0.000 6時(shí)，潰口發(fā)展到6.2 h達(dá)到洪峰流量10 882.78 m3/s，最終潰口水面寬度為99.6 m；b為0.000 5時(shí)，5.65 h達(dá)到洪峰12 261.51 m3/s；b為0.000 4時(shí)，洪峰流量最大，5.16 h即達(dá)到洪峰流量14 219.86 m3/s。不同沖刷侵蝕參數(shù)b對(duì)潰口洪峰流量計(jì)算值的影響敏感度約為30%。

GST模型主要輸入?yún)?shù)資料為地形資料、入流資料、模擬參數(shù)，地形資料采用白格至拉哇下游處的數(shù)字地形高程數(shù)據(jù)(DEM)，全長(zhǎng)約137 km，共計(jì)187萬(wàn)網(wǎng)格。上游邊界為入流邊界，采用DB-IWHR模型計(jì)算結(jié)果，下游邊界為開(kāi)邊界的自然出流，其余邊界為封閉邊界。曼寧系數(shù)選用0.02，庫(kù)郎數(shù)為0.5[15]，模擬時(shí)間為27 h，計(jì)算步長(zhǎng)0.01 s。

下游模擬流量與實(shí)測(cè)結(jié)果如圖7所示。

圖7 模擬流量與實(shí)測(cè)過(guò)程對(duì)比Fig.7 Simulation and measured discharges of the dam break flood propagation

可以發(fā)現(xiàn)，圖7(a)的計(jì)算成果與下游葉巴灘、拉哇的實(shí)測(cè)值最為接近。相應(yīng)此工況，b為0.000 6。13日14:00洪水演進(jìn)至拉哇，洪峰流量為6 786.24 m3/s，比實(shí)測(cè)推遲1 h左右，洪峰流量誤差約為8 %。為此，我們將該工況相應(yīng)潰口洪峰流量10 882.78 m3/s以及圖6中b=0.000 6的過(guò)程線作為“10·10”白格堰塞湖潰決的反演成果。

4 結(jié) 論

本文以下游在建水電站處實(shí)測(cè)的洪水過(guò)程為依據(jù)，分別采用沖刷侵蝕參數(shù)b為0.000 6,0.000 5和0.000 4，通過(guò)DB-IWHR潰壩洪水分析程序和GST洪水演進(jìn)模型，對(duì)“10·10”白格堰塞湖漫頂自然泄流過(guò)程進(jìn)行了反演分析。

(1) 當(dāng)沖刷侵蝕參數(shù)為a=1.100 0、b=0.000 6時(shí)，潰決洪水模擬至葉巴灘和拉哇的模擬結(jié)果和實(shí)測(cè)洪水流量過(guò)程最為接近。據(jù)此判斷“10·10”白格堰塞湖歷時(shí)6.2 h潰決到達(dá)洪峰流量10 882.78 m3/s，最終潰口水面寬度為99.66 m。

(2) 從不同沖刷侵蝕參數(shù)的潰壩洪水演進(jìn)模擬結(jié)果可以發(fā)現(xiàn)：在潰決洪峰流量相差較大的情況下，洪水演進(jìn)至拉哇處的洪峰流量相差較小。由此說(shuō)明，洪峰流量越大，洪水演進(jìn)過(guò)程中坦化作用越明顯，洪峰流量到達(dá)時(shí)間提前。

(3) DB-IWHR是基于EXCEL程序VBA開(kāi)發(fā)的潰壩洪水分析軟件，在輸入合適的參數(shù)后，即可完成計(jì)算，得到潰口流量過(guò)程。GST模型采用了GPU加速技術(shù)，在白格至拉哇137 km的大尺度空間，共計(jì)187萬(wàn)網(wǎng)格，計(jì)算時(shí)間僅半小時(shí)，計(jì)算效率大大提高。DB-IWHR潰壩程序結(jié)合GST模型在堰塞湖應(yīng)急搶險(xiǎn)的過(guò)程中，可以實(shí)現(xiàn)快速、精準(zhǔn)的預(yù)測(cè)。