崩滑堰塞湖的形成-孕災-致災機理與模擬方法

鐘 啟 明，錢 亞 俊，單 熠 博

(1.南京水利科學研究院，江蘇 南京 210029； 2.水利部土石壩破壞機理與防控技術重點實驗室，江蘇 南京 210029)

1 研究背景

堰塞湖是一定量固體物質阻塞山區(qū)河谷或河道導致上游壅水形成的具有一定庫容的水體，而堵塞河谷或河道的固體物質被稱為堰塞體[1-2]。堰塞湖在世界范圍內廣泛分布，依據成因，堰塞湖可分為熔巖、滑坡、崩塌、泥石流和冰磧堰塞湖等5類[3]，全球350個大型堰塞湖分布如圖1所示[4-5]。世界范圍內1 393個堰塞湖案例的統計數據表明[6]，形成堰塞湖的誘因依次是地震(50.5%)、降雨(39.3%)、融雪(2.4%)、人為原因(2.2%)、火山噴發(fā)(0.9%)，其他未知原因的占4.7%。由此可以看出，地震和降雨是形成堰塞湖的主導因素，兩種成因的堰塞湖約占總數的90%。

圖1 全球350個大型堰塞湖分布圖[4-5]Fig.1 The distribution map of 350 large dammed lakes around the world[4-5]

美國地質調查局Costa和Schuster對全球73個堰塞湖的壽命統計發(fā)現[1]，85%的堰塞湖壽命小于1 a。我國學者Peng和Zhang[7]、石振明等[8]、Shen等[6]也分別通過對全球204，276個和352個堰塞湖的壽命統計得出了類似結論。堰塞湖的破壞模式包括水流漫頂沖刷、滲透破壞或邊坡失穩(wěn)，其中89%為水流漫頂沖刷潰決，10%為滲透破壞潰決[1]。

近年來，受地形地貌、地質構造及氣象水文等條件綜合作用，我國堰塞湖事件呈多發(fā)、頻發(fā)態(tài)勢。據作者統計，21世紀以來(2000年1月至2019年12月)，我國發(fā)生有文獻記載的堰塞湖案例362個。2000年4月9日，西藏波密縣易貢鄉(xiāng)發(fā)生巨型山體滑坡，形成易貢堰塞湖[2]，并于同年6月10日潰決，21億m3洪水下泄，潰口峰值流量達到12.4萬m3/s，導致我國墨脫、波密、林芝3縣90余鄉(xiāng)近萬人受災，印度布拉馬普特拉河沿岸7個邦94人死亡，250萬人無家可歸[5，9]；2008年，“5·12”汶川地震形成了257處崩滑堰塞湖[10]，其中唐家山是集雨面積最廣、蓄水量最大、威脅最嚴重的堰塞湖，在人工干預下于2008年6月7日應急泄流，由于泄流前緊急轉移下游風險人口約25萬人，所幸未造成人員傷亡[11]；2018年10～11月，我國金沙江和雅魯藏布江各接連發(fā)生兩次滑坡事件，形成了白格和加拉堰塞湖，并在短期內發(fā)生潰決，對西藏、四川和云南的人民群眾生命財產安全構成巨大威脅[12-13]。

與人工填筑的土石壩不同，堰塞體一般由崩滑土石料快速堆積而成，沒有經過充分壓實，結構較為松垮、組成物質雜亂，局部存在由大顆粒骨架架空形成的高滲透區(qū)域，滲流和力學穩(wěn)定性較差[14]，且缺乏必要的洪水溢流設施，容易發(fā)生潰決造成嚴重的洪水災害，對下游公眾生命財產和基礎設施安全構成巨大威脅。

因此，快速地評估堰塞體的穩(wěn)定性、準確地預測潰口流量過程對堰塞湖的科學應急處置至關重要。本文針對崩塌型和滑坡型這兩類發(fā)生頻率最高的堰塞湖，重點圍繞崩滑堰塞湖形成-孕災-致災機理與模擬方法開展研究，總結其形成機理，剖析影響其穩(wěn)定性的關鍵特征參數，并介紹作者提出的堰塞體穩(wěn)定性快速評價方法和堰塞湖潰決過程數值模擬方法，選擇我國21世紀形成的3個典型崩滑堰塞湖案例，利用作者的研究成果對其災害鏈進行模擬分析。

2 崩滑堰塞湖形成機理與堰塞體顆粒組成

2.1 崩滑堰塞湖形成機理

研究表明，崩滑堰塞體的幾何形態(tài)、粒徑分布和材料力學性質等重要特征取決于其形成過程，并極大地影響了堰塞湖的孕災和致災過程[15-16]。根據國內外常用的分類方法[17-19]，并考慮崩滑體的運動特征，可將崩滑堰塞湖的形成原因細分為3類：崩塌、滑坡和碎屑流。

崩塌形成堰塞湖是指岸坡上部巖土體被裂縫切割、拉裂后，在外荷載作用下失去穩(wěn)定，脫離母巖急劇向下翻滾、跳躍，從而堵塞河道[19-20]?；滦纬裳呷侵赴镀略谕夂奢d作用下，當坡內軟弱結構面上的剪應力超過了抗剪強度，滑坡體沿滑裂面發(fā)生整體滑動。一般可將滑坡形成堰塞湖的過程概括為岸坡上部拉裂、下部滑移隆起、中段快速剪斷及整體的高速滑動[21]。碎屑流形成堰塞湖是指在滑坡或崩塌過程中，坡體物質在高速差異滑移速度和相互碰撞作用下產生平移剪切運動、跳躍及滾動等綜合流動形式，將高速滑坡碎屑流沖到河對岸并阻擋停積于河床上，從而堵塞河道[18-19]。

2.2 堰塞體顆粒分布特征

目前，國內外對崩滑堰塞湖基本特性的研究主要集中在堰塞湖事件的統計與類型識別、堰塞體的形態(tài)特征、堰塞湖形成的影響因素等方面[16]。調查研究發(fā)現[2，5]，堰塞體的顆粒組成對其穩(wěn)定性和潰決過程有重要影響。一般而言[16]，巖質整體滑坡形成的堰塞體穩(wěn)定性好，潰決時潰口發(fā)展速度慢、潰決程度低；少量崩塌塊石堆積形成的堰塞體可達成入流和滲流平衡而避免潰決；崩滑碎屑流形成的堰塞體易于快速潰決，危害性更大。

我國學者Fan等[15，22]基于崩滑發(fā)生區(qū)域的地質條件、堰塞體的地貌特征和材料物理力學特性，根據崩滑堰塞體的顆粒分布特征將2008年“5·12”汶川地震形成的堰塞體分為3類：第1類，堰塞體由于深部巖體的滑塌導致，內部結構可分為2層或3層，底部由較為完整的巖層構成，頂部為巖石碎屑和細顆粒土體，如圖2(a)所示；第2類，堰塞體由邊坡崩塌滑落的石塊形成，內部結構可分為兩層，底部為高度破碎的巖石碎屑和細顆粒土體，頂部由顆粒較大的碎石組成，如圖2(b)所示；第3類，堰塞體由遠程高速滑坡的碎屑流構成，由于滑坡體經過較長的移動，堰塞體顆粒較為松散且含有較多細顆粒，如圖2(c)所示。

圖2 基于顆粒分布特征的堰塞體分類Fig.2 Classification of landslide dams based on particle size distribution characteristics

3 崩滑堰塞體穩(wěn)定性快速評價方法

崩滑堰塞體形成后，其在外荷載作用下的穩(wěn)定性是揭示堰塞湖孕災機理的核心。1999年，Casagli和Ermini[23]最早提出了堆積指標法(BI)來判斷堰塞體的穩(wěn)定性，并改進提出了無量綱堆積體指標法(DBI)[4]。隨后，國內外學者基于已潰和未潰的堰塞湖資料，采用統計學的方法提出了一系列崩滑堰塞體穩(wěn)定性的快速評價方法[24]，如Korup[25]、Dong等[26]和Stefanelli等[27]提出的評價方法，這些評價方法大多基于堰塞體的地貌學指標和堰塞湖的水動力條件，未考慮堰塞體的顆粒組成。目前國內外常用的堰塞體穩(wěn)定性快速評價方法如表1所列。

表1 國內外常用堰塞體穩(wěn)定性快速評價方法對比Tab.1 Comparison of rapid evaluation methods for stability of landslide dams at home and abroad

由于堰塞湖的成因各異，堰塞體的顆粒分布特征存在巨大差異，例如主要由大塊石組成的堰塞體較土質堰塞體具有更好的穩(wěn)定性，但目前尚缺乏能合理考慮堰塞體顆粒組成的評價方法。作者通過考慮堰塞體的形態(tài)特征、顆粒組成及被阻塞河道的水動力條件，采用邏輯回歸的數值計算方法，建立了一套新的堰塞體穩(wěn)定性快速評價方法[28]。根據可獲取的堰塞體顆粒組成信息的多寡，該快速評價方法又可分為精細化和簡化方法。

3.1 考慮顆粒組成的堰塞體穩(wěn)定性精細化快速評價方法

若堰塞體材料的顆分情況已知，則選擇高度、寬度(即順河向距離)和體積來表征堰塞體的形態(tài)特征，選擇堰塞湖流域面積來表征被阻塞河道的水動力特征，選擇d90，d60，d30和d5等代表粒徑表征顆粒組成特征，Shan等建立了堰塞體穩(wěn)定性精細化快速評價方法，數學表達式如下[28]：

LsIVAS=-0.264lgI+1.166lgVd

-1.551lgAb-0.168lgSd-4.847

(1)

式中：Ls(IVAS)為堰塞體穩(wěn)定性精細化評價指標，當Ls(IVAS)>0時，堰塞體是穩(wěn)定的，當Ls(IVAS)<0時，堰塞體是不穩(wěn)定的；I為堰塞體高度與寬度的比值，即I=Hd/W；Sd為精細化顆粒組成指標，Sd=(d90-d60)/(d30-d5)。

3.2 考慮顆粒組成的堰塞體穩(wěn)定性簡化快速評價方法

若堰塞體材料的顆粒特征僅為定性描述，則同樣選擇高度、寬度和體積來表征堰塞體的形態(tài)特征，選擇堰塞湖流域面積來表征被阻塞河道的水動力特征，選擇顆粒特征參數來表征顆粒組成特征，作者建立了堰塞體穩(wěn)定性簡化快速評價方法，數學表達式如下[28]：

LsIVAM=-0.198lgI+1.387lgVd

-1.432lgAb+4.169Mi-8.674

(2)

式中：Ls(IVAM)為堰塞體穩(wěn)定性簡化評價指標，當Ls(IVAM)>0時，堰塞體是穩(wěn)定的；當Ls(IVAM)<0時，堰塞體是不穩(wěn)定的。Mi為簡化顆粒組成指標，當堰塞體材料以塊石為主時(粒徑>200 mm的顆粒重量超過50%，且粒徑>2 mm的顆粒重量超過80%)，Mi取0.75～1.00；當堰塞體材料為塊石加粗粒土時(粒徑>200 mm的顆粒重量不超過50%，粒徑>2 mm的顆粒重量超過80%)，Mi取0.50～0.75；當堰塞體材料以粗粒土為主時(粒徑>2 mm的顆粒重量占20%～80%)，Mi取0.25～0.50；當堰塞體材料以細粒土為主時(粒徑≤2 mm的顆粒重量超過80%)，Mi取 0～0.25；當顆粒含量位于某一取值范圍之內，Mi采用線性插值法進行計算。

4 崩滑堰塞湖潰決機理

近年來，由于崩滑堰塞湖案例頻發(fā)，國內外學者們圍繞堰塞體的材料沖蝕特性、潰口形態(tài)和流量的演化規(guī)律開展了一系列的研究[29-30]。主要的研究手段包括：現場測試、小尺度模型試驗和離心模型試驗。

2008年“5·12”汶川地震形成的唐家山堰塞湖和2018年“11·03”白格堰塞湖都是在人工開挖泄流槽的情況下發(fā)生了潰決，科技人員通過各種手段獲取了現場實測資料[11，31]，為崩滑堰塞湖潰決機理的研究提供了第一手的寶貴資料。以“11·03”白格堰塞湖泄流過程為例，可將潰決過程分為3個階段[32]：

(1) 均勻沖蝕階段。在沖蝕的初期階段，由于潰口(初始泄流槽)水位較低、流速較小，僅堰塞體表層細顆粒被帶走，潰口未發(fā)生明顯下切現象，漫溢水流表現出均勻沖蝕的特征，此階段潰口出流量小于來流量，堰塞湖水位繼續(xù)抬升，如圖3(a)所示。

(2) 溯源沖蝕階段。隨著上游水位的不斷抬升，流速持續(xù)增大，由于堰塞體下游坡腳處的流速最大，因而率先發(fā)生沖蝕而導致結構性破壞，隨后堰塞體發(fā)生溯源沖蝕直至達到上游坡頂部，如圖3(b)所示。

(3) 沿程沖蝕階段，當溯源沖蝕發(fā)展到潰口頂部后，由于潰口底高程的下降，水頭突然增大，導致流量猛增直至峰值流量，潰口不斷下切和展寬，并伴隨潰口邊坡的失穩(wěn)，由于潰口流量大于來流量，堰塞湖水位不斷下降直至漫溢水流無法沖蝕堰塞體材料為止，如圖3(c)所示。

圖3 2018年“11·03”白格堰塞湖潰決過程Fig.3 Breach process of the Baige landslide dammed lake at November 3，2018

另外，國內外學者開展了一系列針對堰塞湖潰決機理的小尺度模型試驗和離心模型試驗[33-37]，通過對潰口在橫斷面和縱斷面的演化過程，可以總結出與現場觀測資料類似的堰塞體潰決機理如下：堰塞體頂部潰口的向下沖蝕和下游坡面的溯源沖蝕→堰塞體下游坡角減小→溯源沖蝕發(fā)展至上游坡頂部→潰口快速下切和展寬→潰口邊坡失穩(wěn)(破壞面為平面)→潰口沖蝕結束。

值得一提的是，與均質土壩潰決機理不同，堰塞體潰決過程中下游壩坡逐漸變緩，且潰口的最終深度受堰塞體顆粒分布特征的影響。

5 崩滑堰塞湖潰決過程數值模擬方法

5.1 模型基本假設

基于崩滑堰塞湖的潰決機理，考慮堰塞體顆粒分布特征，在堰塞湖潰決過程數值模擬中采用以下假設：① 為簡便起見，將堰塞體的橫斷面和縱斷面分別視為倒梯形和正梯形；② 基于崩滑類型，結合地質調查獲得的堰塞體顆粒分布特征，對崩滑堰塞體進行分層，且各層呈水平狀分布；③ 各層頂部和底部的沖蝕系數由原位試驗或經驗公式確定，且各層的沖蝕系數由上至下呈線性變化；④ 在堰塞體潰決過程中，潰口邊坡坡角在失穩(wěn)前一直保持不變；⑤ 由于各層土體力學性質不同，潰口邊坡失穩(wěn)后的坡角由潰口所在位置處的土體特性確定；⑥ 潰口邊坡失穩(wěn)時沿平面滑動破壞，且滑落塊體瞬間被沖走；⑦ 下游坡的溯源沖蝕速率等于潰口底部的向下沖蝕速率除以該時刻下游坡的坡比(垂直/水平)；⑧ 下游坡坡腳處的沖蝕深度為0，沖蝕深度沿下游坡腳向上至下游坡頂部呈線性增大，根據潰口底部和溯源沖蝕的深度確定下游坡坡角的減小量。

5.2 數值模擬方法

基于數值模擬的假設，建立了一個可描述崩滑堰塞湖潰決過程的數值模擬方法。該模擬方法主要包括3個組成部分[32，38]：水動力模塊、材料沖蝕模塊和潰口發(fā)展模塊。采用按時間步長迭代的計算方法模擬崩滑堰塞湖潰決過程中的水土耦合作用，輸入初始參數，設置計算時長tc和時間步長Δt，計算流程如圖4所示。圖4中，zb為潰口底部高程，zs為堰塞湖水位，t為時間，kd為沖蝕系數，τb為頂部潰口底床水流剪應力，τc為土顆粒臨界剪應力，m為潰口邊坡系數(水平/垂直)，H為堰塞體頂部潰口處水深(H=zs-zb)，Qb為潰口流量，Qin為入流量，As為堰塞湖湖面面積，B為潰口頂寬，b為潰口底寬，nloc為潰口位置參數(nloc=1表示潰口只能向一側發(fā)展，nloc=2表示潰口可向兩側發(fā)展)，Fd為驅動力，Fr為抵抗力，α為失穩(wěn)后潰口邊坡的坡角，β為潰口邊坡坡角，φ為土體的內摩擦角，C為土體的黏聚力，γs土體的容重，Hs為潰口邊坡高度。

圖4 崩滑堰塞湖潰決過程數值模擬計算流程Fig.4 Flow chart of numerical simulation of landslide dammed lake breach process

6 我國21世紀典型崩滑堰塞湖案例分析

為了對本文提出的崩滑堰塞湖形成-孕災-致災模擬方法的可行性進行驗證，針對不同的崩滑方式，分別選擇2008年“5·12”汶川地震形成的唐家山堰塞湖[5，11]、小崗劍堰塞湖[5，39-40]和2018年“11·03”白格堰塞湖[12，41]等3個典型崩滑堰塞湖，開展堰塞體穩(wěn)定性快速評價和堰塞湖潰決過程模擬。

按照崩滑堰塞湖的形成機理和地質調查資料獲取的堰塞體顆粒分布特征，根據前文2.2中的分類標準，唐家山屬于第1類堰塞體，小崗劍屬于第2類堰塞體，“11·03”白格屬于第3類堰塞體。

6.1 典型堰塞體穩(wěn)定性快速評價

分別采用作者提出的考慮顆粒組成的堰塞體穩(wěn)定性精細化(式(1) )和簡化快速評價方法(式(2) )對3個崩滑堰塞體的穩(wěn)定性進行評價，并選擇國內外常用的堰塞體穩(wěn)定性評價方法，如Casagli和Ermini[23]提出的BI法、Ermini和Casagli[4]提出的DBI法、Korup[25]提出的Is法和Ia法、Dong等[26]提出的Ls(AHWL)和Ls(AHV)法、Stefanelli等[27]提出的HDSI法，與作者提出的方法進行比較。堰塞體穩(wěn)定性快速評價計算參數如表2所列，各種評價方法的對比情況如表3所列。

表2 各評價方法計算參數Tab.2 The calculation parameters for each evaluation method

表3 各評價方法計算結果對比Tab.3 Comparison of calculation results of each evaluation method

調研結果表明，唐家山、小崗劍和“11·03”白格堰塞湖均為高危堰塞湖，若不開挖泄流槽也會自然潰決，因此其堰塞體均為不穩(wěn)定堰塞體。計算結果表明：本文提出的Ls(IVAS)法和Ls(IVAM)對3個堰塞體的穩(wěn)定性給出了準確的評價結論，同樣，DBI法、Ia法、Ls(AHWL)法和Ls(AHV)法也給出了準確的評價結論；但是，BI法對唐家山和“11·03”白格堰塞體給出了錯誤的評價結論，Is法無法對3個堰塞體的穩(wěn)定性給出明確的評價結論，HDSI無法對唐家山堰塞體的穩(wěn)定性給出明確的評價結論。

6.2 典型堰塞湖潰決過程分析

唐家山堰塞湖通過開挖泄流槽的方式，于2008年6月10日泄流；小崗劍堰塞湖通過人工爆破的方式，于2008年6月12日泄流?！?1·03”白格堰塞湖亦采用人工開挖泄流槽的方式，于2018年11月12日泄流。地質調查資料顯示，唐家山堰塞體自上而下大致可分為3層，頂部以碎石土為主，第2層為塊碎礫石土，第3層為似層狀碎裂巖體；小崗劍堰塞體主要以孤石和塊石為主，含少量碎石，細粒土填充于塊石骨架之間，且頂部含有較多大粒徑石塊，自上而下大體可分為兩層；白格堰塞體是由高位滑坡形成，滑坡土體經過長距離的運移顆粒發(fā)生了崩解，以砂礫石夾碎石土為主，整體顆粒較小，自上而下大體可分為兩層。唐家山、小崗劍和“11·03”白格堰塞湖的潰決過程計算參數如表4和表5所列。

表4 我國3個典型堰塞湖物理力學參數Tab.4 Physical and mechanical parameters of three typical landslide dammed lakes in China

表5 我國3個典型堰塞體各層計算參數Tab.5 Calculation parameters of each layer of three typical landslide dams in China

為了驗證模型的合理性，提取計算結果中的潰口峰值流量(Qp)、潰口最終頂寬(Bf)、潰口最終底寬(bf)、潰口最終深度(Df)以及潰口峰值流量出現時刻(Tp)等堰塞湖潰決特征參數與實測值進行比對(見表6)。

由表6可以看出：唐家山堰塞湖潰決時的潰口峰值流量、潰口最終頂寬、底寬和深度的計算誤差控制在±5%以內，峰值流量出現時刻的計算誤差控制在±10%以內；小崗劍堰塞湖潰決時的潰口峰值流量的計算誤差控制在±5%以內，潰口最終底寬和深度、峰值流量出現時刻的計算誤差控制在±10%以內；白格“11·03”堰塞湖潰決時的潰口峰值流量、潰口最終頂寬和峰值流量出現時刻的計算誤差控制在±5%以內，潰口最終深度的計算誤差控制在±15%以內，僅潰口最終底寬的誤差較大，超過±30%。綜上可知，崩滑堰塞湖潰決過程數值模擬方法可較好地反饋分析3個堰塞湖的潰決過程，比對結果驗證了數值模擬方法的合理性。

表6 我國3個典型堰塞湖潰決特征參數計算值與實測值比較Tab.6 Comparison of the calculated and measured breach characteristic parameters of three typical landslide dammed lakes in China

7 結論與建議

本文圍繞崩滑堰塞湖形成機理、堰塞體穩(wěn)定性快速評價方法、堰塞湖潰決機理和潰決過程數值模擬方法進行了系統的闡述，介紹了作者提出的崩滑堰塞湖形成-孕災-致災過程數值模擬方法并驗證了其合理性。主要結論和建議如下。

(1) 崩滑堰塞湖風險評估對于應急搶險具有重要的指導意義，可圍繞其形成-孕災-致災的災害鏈開展，其中的重要環(huán)節(jié)就是對堰塞體穩(wěn)定性的快速評價和對堰塞湖潰決過程的模擬。研究表明，崩滑堰塞體的形態(tài)特征、顆粒組成及被阻塞河道的水動力條件是其穩(wěn)定性的決定因素；另外，顆粒分布特征決定了堰塞體在潰決洪水作用下的材料沖蝕特性和潰口演化規(guī)律。

(2) 本文介紹的崩滑堰塞湖形成-孕災-致災模擬方法以崩滑堰塞湖形成機理為切入點，重點分析堰塞體材料顆粒組成和分布特征，在此基礎上通過考慮堰塞體的形態(tài)特征、顆粒組成及被阻塞河道的水動力條件，建立了其穩(wěn)定性快速評價方法；通過顆粒分布特征將堰塞體分層，建立了基于水動力模塊、材料沖蝕模塊和潰口發(fā)展模塊的堰塞湖潰決過程模型，并通過典型案例驗證了其合理性。

(3) 崩滑堰塞湖的災害鏈牽涉到材料和荷載不確定性問題，以及復雜的水土耦合問題，應在材料測試設備與方法、沖蝕試驗系統與技術和潰決模擬理論與算法等方面開展深入研究，更加科學地評估其形成-孕災-致災過程。