任思宇， 姜 亮， 翟勝?gòu)?qiáng)， 葉曉華， 郭 江

（1. 四川省地質(zhì)環(huán)境調(diào)查研究中心，四川 成都 610036； 2. 四川省地質(zhì)災(zāi)害防治工程技術(shù)研究中心，四川 成都 610036）

0 引言

受全球氣候變暖的影響，喜馬拉雅山地區(qū)冰川急劇退縮，冰川消融與后退不僅會(huì)影響冰川徑流變化［1］，而且會(huì)導(dǎo)致由冰磧物攔擋的冰川湖的形成［2］，由此引發(fā)的冰湖潰決災(zāi)害也隨之增加［3］。冰湖潰決洪水具有突發(fā)性強(qiáng)、規(guī)模大、破壞力強(qiáng)和危害范圍廣等特點(diǎn)［4］。潰決洪水在演進(jìn)的過(guò)程中，由于海拔落差大，潰決洪水對(duì)溝床沉積物進(jìn)行沖刷侵蝕，進(jìn)而形成泥石流。在泥石流形成、運(yùn)動(dòng)和堆積的過(guò)程中又會(huì)激發(fā)其他次生災(zāi)害，從而形成災(zāi)害鏈，因此致災(zāi)能力非常強(qiáng)［5］。波曲河流域內(nèi)“極高”危險(xiǎn)的冰湖有10 個(gè)，“高”危險(xiǎn)的有22 個(gè)［6］，其下游為聶拉木縣及連接中國(guó)和尼泊爾的重要通道，該區(qū)域受冰湖影響危險(xiǎn)性很高，波曲河流域中下游左側(cè)的一級(jí)支流章藏布于1964 年、1981 年和1983 年先后三次暴發(fā)大規(guī)模冰湖潰決災(zāi)害［7］。其中1981 年的次仁瑪錯(cuò)冰湖潰決最為嚴(yán)重［8］，并引發(fā)泥石流、滑坡等次生災(zāi)害。潰決洪水沖毀下游50 km 范圍內(nèi)的中尼公路和交通設(shè)施，誘發(fā)了多處古滑坡復(fù)活，并沖毀曲鄉(xiāng)附近的建筑及尼泊爾境內(nèi)的遜科西水電站，導(dǎo)致尼泊爾200 人死亡［9］。近年來(lái)流域氣候由干冷向濕熱轉(zhuǎn)化，章藏布冰川融水不斷補(bǔ)給次仁瑪錯(cuò)冰湖，導(dǎo)致冰湖面積不斷擴(kuò)張，除此之外章藏布冰川存在多組裂隙，融水匯入裂隙，增加了冰崩冰滑坡的風(fēng)險(xiǎn)，因此對(duì)冰湖潰決洪水防災(zāi)減災(zāi)志在必行。

因此，要定量評(píng)估冰湖潰決洪水的危險(xiǎn)性，就需要對(duì)其演化過(guò)程及動(dòng)力特征進(jìn)行量化分析，從而制定合理有效的防災(zāi)減災(zāi)措施［10］。近年來(lái)，有關(guān)冰湖潰決的研究正逐步從以分析冰湖潰決災(zāi)害案例中的泥石流成因［11］、冰湖潰決模式［12］、機(jī)理分析［13］等為主的定量分析階段向以物理實(shí)驗(yàn)［14］、數(shù)值模擬［8］為主的冰湖潰決洪水災(zāi)害定量預(yù)測(cè)階段發(fā)展。

目前國(guó)內(nèi)外學(xué)者對(duì)冰湖潰決洪水的模擬演進(jìn)進(jìn)行了大量的研究，如楊瑞敏等［15］利用BREACH 模型和SMPDBK 模型估算和模擬來(lái)古冰湖潰決洪水并對(duì)其進(jìn)行災(zāi)害預(yù)警分析，BREACH 模型可以模擬潰決時(shí)壩體潰口特征及潰口處洪峰流量，SMPDBK 模型可以模擬潰決洪水的洪峰演進(jìn)及洪水深度，但不能反映洪水演進(jìn)過(guò)程中洪水水淹范圍、洪水速度，及洪水對(duì)溝床侵蝕等特征，在潰決洪水災(zāi)害評(píng)估中受到了極大限制；舒有鋒［16］運(yùn)用FLDWAV模型模擬了桑旺錯(cuò)潰決洪水，模型可以考慮多種潰決洪水的啟動(dòng)方式和壩體潰決破壞過(guò)程，但模型不能揭示潰決洪水在演進(jìn)過(guò)程中流深的變化及洪水溝床的侵蝕特征；Sattar 等［17］利用HEC-RAS 一維和二維模型對(duì)南洛納克湖進(jìn)行潰決模擬，模型可以考慮不同潰壩模式，可以采用一維模型模擬潰口處流量隨時(shí)間的變化趨勢(shì)，并將一維模型計(jì)算的結(jié)果導(dǎo)入二維模型，模擬潰決洪水的演進(jìn)過(guò)程，模型能反映潰決洪水演進(jìn)過(guò)程流深和流速變化，能較好地對(duì)冰湖潰決洪水災(zāi)害進(jìn)行定量分析，但模型難以反映潰決洪水侵蝕演化動(dòng)力過(guò)程。歐陽(yáng)朝軍采用連續(xù)介質(zhì)開(kāi)發(fā)的Massflow 數(shù)值平臺(tái)可以根據(jù)用戶需求對(duì)本構(gòu)模型進(jìn)行二次開(kāi)發(fā)，用戶可根據(jù)自身需求獲得相應(yīng)模擬參數(shù)，在滑坡、泥石流、洪水領(lǐng)域得到了廣泛應(yīng)用。

本文基于地表過(guò)程數(shù)值仿真軟件Massflow，根據(jù)實(shí)地調(diào)查和冰磧土力學(xué)參數(shù)對(duì)潰決洪水模型進(jìn)行二次開(kāi)發(fā)，考慮了冰湖潰決洪水在演進(jìn)的過(guò)程中對(duì)溝床及溝岸松散沉積物及坡積物的侵蝕。模型將泥沙輸移模型與水動(dòng)力模型進(jìn)行耦合，使模型更符合冰湖潰決洪水動(dòng)力演化特征研究。通過(guò)建立冰湖潰決洪水動(dòng)力模型，模擬反演次仁瑪錯(cuò)冰湖潰決洪水的演進(jìn)過(guò)程，著重分析潰決洪水在下游基礎(chǔ)設(shè)施處的洪峰流量、流速、流深、侵蝕、沉積及潰決洪水的傳播等特征，將模擬結(jié)果與實(shí)地調(diào)查結(jié)果進(jìn)行對(duì)比分析，驗(yàn)證了模型的可行性與適用性。對(duì)冰湖潰決洪水危險(xiǎn)性進(jìn)行定量評(píng)估，為冰湖潰決洪水防災(zāi)減災(zāi)提供了理論依據(jù)。

1 研究區(qū)域概況

1.1 地質(zhì)背景

次仁瑪錯(cuò)冰湖（86° 03′35″～86° 4′23″ E，28°03′28″～28°4′10″ N）位于波曲河流域章藏布支溝（85°44′06″～85°59′12″ E，28°06′30″～28°19′18″ N）內(nèi)。冰湖從西向東延伸，位于支溝的一個(gè)小冰川侵蝕的支流溝中，北部?jī)H3.97 km2。在這個(gè)小支流中，由于冰川周期性退縮和發(fā)展形成了次仁瑪錯(cuò)冰湖。冰湖距離中尼友誼橋約24 km，溝道平均縱比降為77.3‰。支溝干流長(zhǎng)8.5 km，海拔從6 109 m 到3 168 m，溝道平均縱比降為176.0‰，匯水面積約50.5 km2。支溝內(nèi)巖性主要為云母石英片巖，黑云母花崗巖片麻巖和第四紀(jì)沖積土。支溝內(nèi)有六條主要支流，河床下切深度4～10 m。在雨季（6—9月），章藏布支溝干流流量為14 m3·s-1。圖1 顯示了次仁瑪錯(cuò)冰湖和章藏布支溝的位置及物源條件，圖2為潛在冰崩冰滑坡航拍圖片。

圖1 次仁瑪錯(cuò)冰湖位置及章藏布支溝地層巖性及物源條件Fig. 1 The location of Cirenmaco glacial lake （a） and the formation lithology （b），and materials conditions of Zhangzangbu branch ditch （c）

圖2 潛在冰崩冰滑坡Fig. 2 Potential ice avalanche and landslide

研究區(qū)域分布兩條斷層如圖1 所示，北部斷層穿過(guò)章藏布。在這樣的地質(zhì)背景下，第四紀(jì)層廣泛分布于溝谷，因此很容易發(fā)生滑坡和塌陷。此外，冰磧臺(tái)地分布在溝谷兩岸，冰磧土分布于冰舌前端，冰磧土顆粒粗，孔隙率高，顆粒級(jí)配級(jí)不均勻（圖3）。次仁瑪錯(cuò)冰湖是由冰川作用形成的冰磧湖，下游由疏松的第四紀(jì)沉積土所覆蓋。

圖3 章藏布冰磧土顆粒級(jí)配Fig. 3 Grain gradation of moraine in Zhangzangbu branch ditch

章藏布支溝具有豐富的物源，為泥石流的形成和發(fā)育提供了有利條件，除此之外溝道兩岸分布眾多潛在滑坡，如707滑坡、樟木隧道滑坡和樟木口岸滑坡群，自尼泊爾4·25（2015 年4 月25 日）地震以來(lái)，古滑坡變形速率加快，次仁瑪錯(cuò)冰湖一旦潰決，其導(dǎo)致的直接災(zāi)害及次生災(zāi)害將嚴(yán)重威脅著樟木口岸的正常運(yùn)營(yíng)（圖4～7）。

圖4 707滑坡Fig. 4 Landslide of 707

圖5 樟木隧道滑坡Fig. 5 Landslide of Zhangmu tunnel

圖6 章藏布支溝物源Fig. 6 Materials of Zhangzangbu branch ditch

圖7 樟木口岸滑坡群Fig. 7 Landslide group of Zhangmu port

1.2 氣象條件

聶拉木氣象站位于章藏布支溝溝口北部10 km2處，本文統(tǒng)計(jì)了氣象站1981—2017 年的年均降雨、年均氣溫和年均濕度數(shù)據(jù)如圖8 和圖9 所示，流域內(nèi)呈年均氣溫逐漸升高，降雨逐漸減小的趨勢(shì)，氣象條件由濕冷轉(zhuǎn)為干熱，冰川融水加劇，冰湖面積加速擴(kuò)張，增加了冰湖潰決洪水的風(fēng)險(xiǎn)。

圖8 聶拉木氣象站1981—2017年的年均降雨和年均氣溫?cái)?shù)據(jù)Fig. 8 Average annual rainfall and average annual temperature data of Nyalam Meteorological Station from 1981 to 2017

圖9 聶拉木氣象站1981—2017年年均濕度數(shù)據(jù)Fig. 9 Annual average humidity data of Nyalam Meteorological Station from 1981 to 2017

1.3 水文條件

波曲流域氣候跨越了亞熱帶、暖溫帶、寒溫帶、高山苔原帶和雪山冰漠帶。受喜馬拉雅山脈地殼隆升的影響，河谷深切、河道彎曲，呈現(xiàn)典型的山區(qū)河流特點(diǎn)：河床坡降大，平均縱坡降達(dá)81‰。中、上游形成相對(duì)切深達(dá)1 000～1 500 m 的高山深谷，多呈“V”侵蝕型谷，一般寬在60～80 m；而下游峽谷深切，相對(duì)切深達(dá)1 500～3 000 m，其橫斷面多呈“V”型，兩岸懸崖峭壁，谷坡坡度大于35°，極有利于雨洪的集流和冰雪運(yùn)移。波曲河流域年平均降水量582.9 mm，最大日降水量107.6 mm，雨季主要集中在6—8 月，降雪時(shí)間長(zhǎng)達(dá)6個(gè)月以上。洪水期為5月中旬至9月中旬，7—8月為洪峰期，主河年均徑流流量31.7 m3·s-1，平均徑流深度約500 mm，平均徑流模數(shù)15.9 L·s-1·km-2，出境年徑流量1×108m3（友誼橋），主河和各主要支溝的徑流變幅較小。

2 冰湖潰決洪水動(dòng)力演化過(guò)程的模型構(gòu)建

2.1 模型構(gòu)建

模型采用全局笛卡爾坐標(biāo)。模型分為上層運(yùn)動(dòng)流流體和底部可侵蝕沉積物，上層自由表面定義為Z1top點(diǎn)，上層和下層之間的邊界為Z1bot，h1=Z1top-Z1bot為流動(dòng)深度。（u1，v1，w1）表示三個(gè)獨(dú)立坐標(biāo)軸上深度平均的速度分量［18］。流動(dòng)層質(zhì)量和動(dòng)量守恒方程可以表示如下［19］：

式中：τ1zxbot和τ1zybot為上部流動(dòng)層在底部靜態(tài)層上施加的基礎(chǔ)剪應(yīng)力；E1bot為流體的侵蝕率；u1（Z1bot）和v1（Z1bot）為進(jìn)入流動(dòng)層材料的邊界速度。u1（Z1bot）和v1（Z1bot）用平均速度分量u1和v1表示。質(zhì)點(diǎn)跳躍條件ρ1E1bot=ρ2E1top，其中E1top下層頂部的侵蝕率，ρ2為下層的密度［20］。

對(duì)方程（1）～（3）進(jìn)行推導(dǎo)，其表達(dá)式為：

式中：E和D為侵蝕和沉積公式；P為溝床沉積物孔隙度。侵蝕和沉積方向垂直于基底表面，當(dāng)采用第1 層的基礎(chǔ)剪應(yīng)力τ1zxbot和τ1zybot時(shí)，在動(dòng)量守恒方程右邊存在動(dòng)量交換項(xiàng)ρ2u1(E-D)和ρ2v1(E-D)和Ouyang［18］。

方程式的左邊去除ρ1，則質(zhì)量和動(dòng)量方程式為：

侵蝕沉積方程通過(guò)刪減冗余項(xiàng)，得到如下方程：

考慮到實(shí)際情況的復(fù)雜性，可將溝床沉積物的漲縮特征以及輕微影響項(xiàng)進(jìn)行忽略，將溝床演化方程簡(jiǎn)化為：

式中：E和D為侵蝕和沉積速率；p為河床沉積物孔隙度。

式中：α為經(jīng)驗(yàn)參數(shù)，為近床濃度和深度平均濃度之間的差異，深度平均濃度由min［m，（1-p）/c］；c1是流動(dòng)中的固體濃度；m是經(jīng)驗(yàn)指數(shù)；ws0是單個(gè)顆粒在平靜水中的沉降速度。它可以表示為［21］：

式中：υ為水的運(yùn)動(dòng)黏度；d為沉積物粒徑；s=ρs/ρw-1；g是重力加速度，取9.8 m·s-2；ce為推移質(zhì)泥沙輸移能力，可表示為：

式中：θc為決定沉積物運(yùn)動(dòng)開(kāi)始的希爾茲參數(shù)的臨界值［22］；θ=(sgd)為屏蔽參數(shù)；u*為摩擦速度。剪應(yīng)力τ1zxbot和τ1zybot表示為［18］：

其中n是曼寧粗糙度。

2.2 求解方法

將控制方程（8）至（17）進(jìn)一步量化為向量模型：

采用MacCormack-TVD 有限差分法來(lái)求解以上方程。與解傳統(tǒng)的淺水方程一樣，首先利用算子分裂法將方程（18）分解為兩個(gè)一維方程，表示如下：

然后，(n+ 1)Δt時(shí)刻由以下方程求解：

執(zhí)行預(yù)測(cè)步驟p：

執(zhí)行修正步驟c：

執(zhí)行平均步驟：

式中：

函數(shù)G的表達(dá)式為：

函數(shù)φ（x）代表最小流體通量限制函數(shù)表達(dá)式為：

式中：x代表自變量，變量C的表達(dá)式示為：

式中：Cr為當(dāng)?shù)氐膸?kù)朗數(shù)（Courant number）表達(dá)式如下：

2.3 模型參數(shù)選取

DEM 數(shù)據(jù)源自ALOS 衛(wèi)星數(shù)據(jù)。采用原始地形精度12.5 m×12.5 m的網(wǎng)格進(jìn)行計(jì)算。通過(guò)實(shí)地調(diào)查和遙感影像（高分二號(hào)）對(duì)河道進(jìn)行分析，模擬范圍河床糙率系數(shù)為0.07。冰磧土作為特殊的巖土材料，其性質(zhì)與溝床沉積物區(qū)別大［23］，由于潰決洪水侵蝕夾帶物源主要為冰磧土，因此我們采集了次仁瑪錯(cuò)冰湖下游冰磧土體，進(jìn)行土工實(shí)驗(yàn)（圖3），由于流域范圍較大，不同河段溝床堆積碎屑顆粒級(jí)配相差巨大，難以在模型中實(shí)現(xiàn)區(qū)別對(duì)待，因此本文采用章藏布流域分布最廣泛的冰磧土作為模型的特征參數(shù)，冰磧物容重為1.6～2.0 g·cm-3，潰決洪水容重1.0 g·cm-3。

表1 次仁瑪錯(cuò)冰湖潰決洪水?dāng)?shù)值模型參數(shù)Table 1 Numerical model parameters of the outburst flood of the Cirenmaco glacial lake

3 模型驗(yàn)證

3.1 次仁瑪錯(cuò)冰湖潰決誘因及危害

次仁瑪錯(cuò)冰湖曾在1964 年、1981 年和1983 發(fā)生過(guò)多次潰決，其中規(guī)模最大的一次為1981 年，冰川融水加速湖水位上升，導(dǎo)致冰磧壩體發(fā)生管涌潰決，潰決歷時(shí)60 min，最大潰決流量為16 000 m3·s-1，潰決洪水在壩體上沖刷了一條長(zhǎng)50 m，底寬40～60 m 的沖溝，潰決洪水在章藏布支溝演進(jìn)的過(guò)程中沖刷侵蝕約2×106m3固體物質(zhì)，其中1.6×106m3的固體物質(zhì)形成泥石流沖入波曲河，成為這次泥石流的主要物源。

潰決洪水在演進(jìn)過(guò)程中，對(duì)河床及溝岸殘坡積物進(jìn)行沖刷側(cè)蝕，導(dǎo)致大規(guī)模的崩塌和滑坡，潰決洪水演進(jìn)至友誼橋附近時(shí)，由于其運(yùn)動(dòng)受阻而雍高達(dá)20 m，比洪峰水位高約5 m，沖毀橋梁及東岸附近的全部建筑物，并引起東岸坡積、殘積層大規(guī)模的崩塌和滑坡［24］。泥石流過(guò)后，友誼橋上游段河谷中留下較大范圍的泥石流沉積。自友誼橋至水電站約8.3 km 的河段，除沿途繼續(xù)強(qiáng)烈的側(cè)向侵蝕外，大量泥石流物質(zhì)在這里沉積，河床顯著升高，次仁瑪錯(cuò)的潰決給下游造成了巨大的經(jīng)濟(jì)損失。

3.2 次仁瑪錯(cuò)冰湖潰決洪水災(zāi)害反演結(jié)果

徐道明［24］對(duì)1981 年次仁瑪錯(cuò)冰湖進(jìn)行了實(shí)地調(diào)查，推測(cè)了潰口處的流量隨時(shí)間的演化特征和洪水在演進(jìn)過(guò)程中洪峰流量隨距離的變化特征，因此在反演過(guò)程中以徐道明實(shí)地調(diào)查數(shù)據(jù)為依據(jù)，采用在潰口處添加流量曲線作為模擬的啟動(dòng)條件。

潰決洪水在下游的演進(jìn)過(guò)程中呈衰減趨勢(shì)，如圖10～11所示為實(shí)測(cè)洪峰流量、流深、流速與模擬洪峰流量、流深、流速隨距離的演化過(guò)程對(duì)比圖，模擬結(jié)果與實(shí)測(cè)結(jié)果誤差范圍小于20%，模擬結(jié)果與實(shí)測(cè)結(jié)果擬合度較好。

圖10 模擬洪峰流量與實(shí)測(cè)流量Fig. 10 Simulated peak flow and measured flow

圖11 模擬流深、流速與實(shí)測(cè)結(jié)果對(duì)比Fig. 11 Comparison of simulated flow depth and velocity with measured results

4 預(yù)測(cè)分析

次仁瑪錯(cuò)冰湖從（0.10±0.08） km2（1988年）擴(kuò)大至（0.35±0.04） km2（2021年），最大水深為（115±2） m，庫(kù)容為1.8×107m3，通過(guò)歷史遙感影像觀測(cè)，冰湖有明顯的擴(kuò)張趨勢(shì)，流域氣象條件變化顯著，由濕冷轉(zhuǎn)向干熱（圖8～9），冰川融水流入冰川冰裂隙中，一方面導(dǎo)致冰川強(qiáng)度降低，另一方面起到潤(rùn)滑作用，冰崩冰滑坡發(fā)生風(fēng)險(xiǎn)逐漸增加。

由于冰湖在潰決的過(guò)程中，洪水流態(tài)及洪水對(duì)壩體的沖刷機(jī)理極為復(fù)雜。模型難以考慮潰口的演化過(guò)程，然而潰口假設(shè)法與實(shí)際情況偏差較大，得出的結(jié)果誤差率高。通過(guò)大量的文獻(xiàn)查閱，基于前人對(duì)大量冰湖潰決案例調(diào)查數(shù)據(jù)的基礎(chǔ)上，本文推導(dǎo)得到經(jīng)驗(yàn)公式法，擬合潰口處的流量曲線，以此作為冰湖潰決洪水啟動(dòng)的參考方法。經(jīng)驗(yàn)公式法需要的參數(shù)較少，便于獲取，通過(guò)經(jīng)驗(yàn)公式（表2）計(jì)算潰決洪水的洪峰流量，潰口寬度及潰決歷時(shí)，擬合最危險(xiǎn)工況（全潰）時(shí)的流量曲線，如圖12 所示，其中V為冰湖庫(kù)容，Hw為冰湖水深，K0為經(jīng)驗(yàn)參數(shù)取值為1。

表2 冰湖潰決洪峰流量、潰口寬度和潰決歷時(shí)經(jīng)驗(yàn)公式Table 2 Empirical formula of glacial lake outburst peak discharge， outburst width and outburst duration

圖12 潰口流量曲線Fig. 12 Flow curve at the breach

圖13 為全潰條件下不同剖面點(diǎn)處的模擬洪峰流量［圖13（a）］、流深和流速［圖13（b）］。潰決洪水對(duì)潰口（No.1）進(jìn)行下切侵蝕，侵蝕深度21～29 m（勘查最大侵蝕深度28 m），洪水模擬最大深度達(dá)16 m，洪峰流量時(shí)刻流速11.4 m·s-1，冰磧壩址（No.2）處地形開(kāi)闊，洪水在該處發(fā)散，流速15 m·s-1、流深13 m，洪水裹挾的大顆粒土體在壩址處堆積，形成沖積扇，洪峰流量從1.80×104m3·s-1降至1.64×104m3·s-1。No.7～No.12 段，河床受冰川融水沖刷而下切侵蝕，河床及溝岸巖體破碎，河床縱比降較大，潰決洪水在該段對(duì)河床基底及坡岸進(jìn)行強(qiáng)烈的沖刷侵蝕，潰決洪水逐漸演化為稀性泥石流，洪峰流量不斷增大。并在章藏布支溝溝口（No.12）和下游形成11 m高泥石流堰塞壩短暫堵塞主河。No.13～No.21 段，稀性泥石流逐漸進(jìn)入峽谷，地形相對(duì)較緩，流深有所增加，流速有所放緩，洪峰流量過(guò)水?dāng)嗝婕眲〗档停榉辶髁吭谘葸M(jìn)過(guò)程中不斷衰減，在No.21～No.22段，泥石流對(duì)滑坡群坡腳進(jìn)行侵蝕夾帶，洪峰流量有微弱的增長(zhǎng)，下游溝床相對(duì)開(kāi)闊，泥石流逐漸沉積，洪峰流量逐漸衰減，泥石流演進(jìn)至水電站時(shí)洪峰流量為9 520 m3·s-1，對(duì)下游仍具有較強(qiáng)的破壞性。

圖13 全潰條件下不同剖面點(diǎn)處的模擬洪峰流量（a）、流深和流速（b）Fig. 13 Simulated peak flow （a）， flow depth and velocity （b） at different profile points under the condition of total collapse

圖14 潰決洪水演進(jìn)過(guò)程沿線侵蝕深度（a）和沉積深度（b）Fig. 14 Erosion depth （a）and siltation depth （b）along the evolution of the outburst flood

707 老橋（No. 11）距離潰口6 440 m，潰決洪水到達(dá)該位置用時(shí)3 300 s，洪峰到達(dá)老橋用時(shí)4 500 s，洪峰流量為1.41×104m3·s-1，洪峰到達(dá)老橋時(shí)流深16.3 m，流速13 m·s-1，稀性泥石流對(duì)溝道左岸707滑坡坡腳進(jìn)行沖刷側(cè)蝕［圖15（c）］，最大侵蝕深度約為9 m，威脅707 滑坡穩(wěn)定性。稀性泥石流沖出章藏布支溝，沖擊波曲河右岸，對(duì)右岸土體進(jìn)行沖刷侵蝕，侵蝕深度約9 m［圖15（c）］，隨后稀性泥石流在章藏布支溝溝口下游沉積（No.12），形成泥石流堰塞壩［圖15（d）］，短暫堵塞波曲河。對(duì)比No.11～No.12段的模擬結(jié)果和實(shí)地調(diào)查結(jié)果，模擬結(jié)果能較為準(zhǔn)確地揭示潰決洪水流深、流速及溝床及溝岸松散沉積物的沖淤特征。因此能較為全面地對(duì)冰湖潰決洪水造成的危害進(jìn)行定量評(píng)估。

圖15 章藏布支溝溝口處洪峰流量時(shí)刻流深（a）、流速（b）和最終侵蝕深度（c）和沉積深度（d）Fig. 15 Flow depth （a）， flow velocity （b）， final erosion depth （c） and silting depth （d）at peak discharge time at the mouth of Zhangzangbu branch ditch

滑坡群（No. 21～No. 23）位于樟木口岸附近，受尼泊爾2015 年4 月25 日8.1 級(jí)地震影響，古滑坡群復(fù)活，通過(guò)模擬結(jié)果可知洪峰到達(dá)該位置用時(shí)3.4 h，由圖16（a）和圖16（b）可以得出洪峰到達(dá)滑坡群時(shí)的最大流深25 m，最大流速為16.8 m·s-1，潰決洪水在該段對(duì)滑坡群坡腳進(jìn)行沖刷掏蝕，最大侵蝕深度達(dá)12.0 m。流體中固體顆粒沉積于溝床，沉積深度高達(dá)11.0 m。潰決洪水對(duì)滑坡群坡腳沖刷側(cè)蝕，引起潛在滑坡群坡腳大規(guī)模崩塌，可能增加滑坡群變形速率或失穩(wěn)概率進(jìn)而誘發(fā)更大規(guī)模的次生災(zāi)害。

圖16 滑坡群洪峰流量時(shí)刻流深（a）、流速（b）和最終侵蝕深度（c）和沉積深度（d）Fig. 16 Flow depth （a）， flow velocity（b）， final erosion depth （c） and siltation depth （d）of landslide group at flood peak discharge time

水電站（No. 24）位于尼泊爾境內(nèi)，距離潰口25.0 km，洪峰到達(dá)該位置用時(shí)4.8 h，洪峰流量為9 520 m3·s-1，流深17 m，流速8 m·s-1，圖17為水電站位置洪峰流量時(shí)刻的流體深度和流體速度及最終的侵蝕深度和淤埋深度。水電站高度約為15 m，稀性泥石流淤滿水電站，沉積深度最大約17 m，導(dǎo)致水電站失效。稀性泥石流漫過(guò)水電站繼續(xù)向下游演進(jìn)，對(duì)下游基礎(chǔ)設(shè)施仍具有較高的破壞性。

圖17 水電站洪峰流量時(shí)刻流深（a）、流速（b）和最終侵蝕深度（c）和沉積深度（d）Fig. 17 Flow depth （a）， flow velocity （b）， final erosion depth （c） and siltation depth （d）at peak discharge time of hydropower station

5 討論

冰湖潰決洪水在演化過(guò)程中由清水流轉(zhuǎn)換為高含沙洪水，再到泥石流［28］。前人采用BASEMENT［29］，HEC-RAS［8］，F(xiàn)LO-2D［30］等模型模擬高原冰湖潰決事件。然而，這些模型無(wú)法進(jìn)行冰湖潰決洪水對(duì)沿途基底侵蝕夾帶的模擬。研究發(fā)現(xiàn)冰湖潰決誘發(fā)的泥石流的峰值流量比僅計(jì)算洪水的結(jié)果要高得多［28］。MASSFLOW 模型則可以模擬潰決洪水對(duì)基底的侵蝕夾帶效應(yīng)，進(jìn)而模擬冰湖潰決泥石流形成演化的動(dòng)力過(guò)程，模型基于深度積分的二維數(shù)值模型，可根據(jù)實(shí)際需求對(duì)模型進(jìn)行改進(jìn)，進(jìn)而考慮流體對(duì)基底的侵蝕效應(yīng)［18］。除此之外，高精度的DEM 對(duì)于準(zhǔn)確模擬冰湖潰決泥石流至關(guān)重要，它可以更加真實(shí)地反映潰決洪水在溝道的侵蝕和沉積狀況，由于模擬范圍較大，普通計(jì)算機(jī)難以實(shí)現(xiàn)高精度模擬，通過(guò)模擬結(jié)果與實(shí)地調(diào)查對(duì)比分析12.5 m 精度地形用于流域風(fēng)險(xiǎn)評(píng)估和危險(xiǎn)性分析可以被接受。

冰湖潰決事件模擬中的輸入?yún)?shù)的選擇具有極大的挑戰(zhàn)性。相對(duì)于經(jīng)驗(yàn)?zāi)Ｐ偷膮?shù)校準(zhǔn)，物理模型則利用質(zhì)量和動(dòng)量守恒等物理原理表示過(guò)程。然而，即使是復(fù)雜的物理模型在應(yīng)用于真實(shí)事件時(shí)也有很大的不確定性，因此模型校準(zhǔn)是必不可少的［31］。在現(xiàn)實(shí)應(yīng)用中，受限冰湖潰決事件高質(zhì)量的現(xiàn)場(chǎng)數(shù)據(jù)獲取，有時(shí)必須做出相應(yīng)的妥協(xié)。模型校準(zhǔn)的關(guān)鍵問(wèn)題在于選擇合適的調(diào)整參數(shù)和定義一個(gè)能被接受的參數(shù)調(diào)整范圍［28］，而敏感性分析可以識(shí)別模型校準(zhǔn)中使用的敏感輸入?yún)?shù)，這需要建模者證明標(biāo)準(zhǔn)參數(shù)值的偏差。

在模型敏感性分析中，通過(guò)調(diào)整曼寧系數(shù)和潰決水量來(lái)觀察模擬結(jié)果的差異，其中曼寧系數(shù)控制著流體的紊流運(yùn)動(dòng)。敏感性分析表明，曼寧系數(shù)變化，泥石流淹沒(méi)面積、流深、流速等參數(shù)變化較小，在宏觀上這些誤差可以被接受。但流體在演化過(guò)程中的時(shí)間具有明顯差異，且泥石流侵蝕距離也顯著延長(zhǎng)，但在泥石流平均流速上差異很?。?0］。同時(shí)，潰決水量的增加將加快了泥石流演進(jìn)的速度。較大的曼寧系數(shù)將導(dǎo)致泥石流到達(dá)時(shí)間延遲，并增加泥石流深度和侵蝕距離，這與前人研究成果一致［32］。分析結(jié)果對(duì)區(qū)域風(fēng)險(xiǎn)評(píng)估、國(guó)土空間規(guī)劃及防災(zāi)減災(zāi)工作的開(kāi)展具有一定的參考意義，但由于冰湖潰決洪水演進(jìn)過(guò)程中的復(fù)雜性質(zhì)和流動(dòng)行為可能發(fā)生的變化，在解釋模型結(jié)果時(shí)仍然需要謹(jǐn)慎處理。

6 結(jié)論

通過(guò)資料收集、遙感解譯、野外調(diào)查等方法對(duì)西藏波曲河流域次仁瑪錯(cuò)冰湖潰決洪水動(dòng)力過(guò)程進(jìn)行研究，對(duì)1981 年洪災(zāi)進(jìn)行模擬反演，并對(duì)次仁瑪錯(cuò)冰湖再次潰決后的洪水演進(jìn)動(dòng)力過(guò)程進(jìn)行預(yù)測(cè)，得到以下結(jié)論：

（1）對(duì)次仁瑪錯(cuò)1981年冰湖潰決洪水進(jìn)行模擬反演，對(duì)比洪峰流量、流深、流速等特征參數(shù)，得出模擬結(jié)果與實(shí)測(cè)結(jié)果具有良好的一致性，驗(yàn)證了模型的適用性和可行性。

（2）次仁瑪錯(cuò)冰湖面積從自1988年以來(lái)急速擴(kuò)張，最大水深為（115±2） m，庫(kù)容為（1.8×107） m3，流域氣候由濕冷轉(zhuǎn)向干熱，冰川融水流入冰川冰裂隙中，增加冰崩冰滑坡發(fā)生風(fēng)險(xiǎn)，進(jìn)而導(dǎo)致冰湖潰決風(fēng)險(xiǎn)增加。

（3）次仁瑪錯(cuò)冰湖一旦再次潰決，洪水對(duì)下游冰磧物及松散坡積物和沉積物進(jìn)行侵蝕夾帶，逐漸演化為稀性泥石流，對(duì)下游潛在滑坡坡腳進(jìn)行沖刷側(cè)蝕，加速潛在滑坡變形。對(duì)下游友誼橋構(gòu)成巨大威脅，稀性泥石流在水電處進(jìn)行沉積，導(dǎo)致水電站失效。

（4）潰決洪水在演進(jìn)的過(guò)程中（No. 7～No. 12）對(duì)冰磧土及溝床沉積物進(jìn)行沖刷侵蝕，洪水逐漸演化為泥石流，洪峰流量逐漸增強(qiáng)。稀性泥石流在（No.12～No.20）段，一方面對(duì)溝床及岸坡松散沉積物進(jìn)行沖刷侵蝕，另一方面稀性泥石流在溝床開(kāi)闊處形成沉積，基本上為沖淤平衡。稀性泥石流到達(dá)潛在滑坡群后（No.20～No.21），對(duì)潛在滑坡群坡腳進(jìn)行沖刷側(cè)蝕，洪峰流量有微弱增長(zhǎng)。泥石流到達(dá)水電站處，洪峰流量仍然有9 520 m3·s-1，對(duì)下游仍具有較強(qiáng)的危害性。

（5）通過(guò)對(duì)次仁瑪錯(cuò)冰湖潰決洪水的動(dòng)力演化過(guò)程分析，潰決洪水對(duì)章藏布支溝松散物質(zhì)進(jìn)行侵蝕夾帶，逐漸演化為稀性泥石流，對(duì)下游707滑坡和樟木口岸滑坡群坡腳進(jìn)行沖刷侵蝕，易誘發(fā)滑坡變形，導(dǎo)致次生災(zāi)害的發(fā)生，除此之外，潰決洪水引發(fā)的泥石流將對(duì)中尼友誼橋和下游水電站的安全運(yùn)營(yíng)構(gòu)成威脅。