金沙江“10·10”白格堰塞湖潰壩洪水反演分析

陳 祖 煜，張 強，侯 精 明，王 琳，馬 利 平

(1.西安理工大學 省部共建西北旱區生態水利國家重點實驗室, 陜西 西安 710048; 2.中國水利水電科學研究院 巖土工程研究所, 北京 100048)

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

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

1 基礎資料

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

在堰塞湖發生后第一時間進行了無人機拍攝，三維點云圖如圖1所示。主堰塞體長度為961.89 m，在潰壩發生后經過水流沖刷，形成長1 622 m、底寬80～120 m的泄流槽，泄流槽口高程約2 888 m、出口底板高程2 872 m，對應堰塞湖蓄水量約0.5億m3。堰塞湖庫容水位關系曲線如圖2所示。

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

圖4 潰口的側向崩塌過程等效簡化Fig.4 Equivalent simplifcation of lateral enlargement process

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

2 潰壩洪水水力學基礎

2.1 潰口洪水分析

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

2.1.1 潰口水力學條件

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

(1)

(2)

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

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

(3)

式中，h為潰口水深。

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

2.1.2 沖刷模型

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

(4)

式中，γ為扣除臨界剪應力后的剪應力。

γ=k(τ-τc)

(5)

2.1.3 潰口側向擴展模型

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

根據近期研究分析，我們提出了用雙曲線模型來計算梯形側面傾角增量Δβ。

(6)

(7)

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

2.2 洪水波演進模型

2.2.1 控制方程

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

(9)

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

2.2.2 數值方法

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

3 潰壩洪水反演分析

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

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

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

注：其中庫容與水位關系用公式表示W=p1(H-Hr)2+p2(H-Hr)+p3

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

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

下游模擬流量與實測結果如圖7所示。

圖7 模擬流量與實測過程對比Fig.7 Simulation and measured discharges of the dam break flood propagation

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

4 結 論

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

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

(2) 從不同沖刷侵蝕參數的潰壩洪水演進模擬結果可以發現：在潰決洪峰流量相差較大的情況下，洪水演進至拉哇處的洪峰流量相差較小。由此說明，洪峰流量越大，洪水演進過程中坦化作用越明顯，洪峰流量到達時間提前。

(3) DB-IWHR是基于EXCEL程序VBA開發的潰壩洪水分析軟件，在輸入合適的參數后，即可完成計算，得到潰口流量過程。GST模型采用了GPU加速技術，在白格至拉哇137 km的大尺度空間，共計187萬網格，計算時間僅半小時，計算效率大大提高。DB-IWHR潰壩程序結合GST模型在堰塞湖應急搶險的過程中，可以實現快速、精準的預測。