強震區(qū)泥石流堵江-洪水災害鏈式過程模擬研究
——以汶川縣下莊溝為例

羅玉婷,王立娟,馬 松,唐 堯

(1.重大危險源測控四川省重點實驗室,四川 成都 610045;2.四川安信科創(chuàng)科技有限公司,四川 成都 610045;3.四川省安全科學技術研究院,四川 成都 610045)

汶川地震是新中國成立以來最具破壞性的地震,產(chǎn)生了大量的崩塌和滑坡[1-2]。由于大量同震滑坡物質和極端天氣條件,泥石流成為汶川震區(qū)常見的地質災害之一[3-5]。在震后十余年,汶川強震區(qū)泥石流發(fā)生的頻率和規(guī)模都急劇增加,對山區(qū)人民的生命財產(chǎn)、基礎設施和重建工作構成了嚴重威脅。強震區(qū)的山地災害一般以災害鏈的形式發(fā)育,泥石流則是鏈式災害效應中的重要環(huán)節(jié)[6],與單一災害相比,鏈式災害具有時間尺度長、范圍廣、破壞力大的特點,造成的人員傷亡與經(jīng)濟損失遠大于各單一災害的破壞效應之和[7]。汶川地震大幅降低了泥石流的降雨閾值,在地震后,汶川、北川、清平等地區(qū)多次因強降雨引發(fā)了大規(guī)模的泥石流災害。泥石流堵江是其致災過程中較為關鍵的環(huán)節(jié),直接影響到其致災的嚴重程度。汶川7.9級大地震后,次生地質災害進入一個長達25年的活躍期,受地震影響區(qū)域將頻繁發(fā)生泥石流災害鏈[8]。

泥石流災害鏈的形成過程極為復雜且具有突發(fā)性,人們往往難以對災害發(fā)生及時作出反應。近年來,隨著遙感觀測和計算機應用技術的發(fā)展,研究人員開發(fā)了許多監(jiān)測和模擬泥石流災害鏈的技術手段,其中,數(shù)值模擬軟件成為了重現(xiàn)并預測泥石流災害鏈的有效工具[9-10]。Cao等[11]通過水流、泥沙輸移和河床演變的耦合,提出了一個淺水流體動力學模型來描述潰壩水流在河床上的運動。江中舟[12]采用SPH的物理方法,提出一種基于粒子的固液兩相泥石流侵蝕模型來模擬泥石流災害過程。楊春陽[13]利用FLO-2D和RAMMS 2種軟件,分別模擬不同暴雨條件下泥石流的活動特征和啟動臨界降雨量。然而,上述研究所涉及模型都只適用于單一災害,難以模擬由幾種具有不同物理機制的次生災害組成的災害鏈。目前,泥石流數(shù)值模擬方面的研究仍在不斷進步、不斷完善。雖然許多軟件平臺可以重現(xiàn)泥石流的運動過程,但大部分都忽視了泥石流所帶來的鏈式效應,無法對災害鏈所帶來的危害進行有效的評估。

在目前的模型方法中,與泥石流有關的過程,如水文過程、滑坡和泥石流運動,大多是單獨模擬的,忽略了不同過程之間所產(chǎn)生的相互作用或反饋,模型的預測能力受到極大的限制。為了更有效地評估災害鏈產(chǎn)生的風險問題,應考慮采取一種綜合集成的方法,同時耦合滑坡、泥石流運動、洪水運動各個過程。鏈生的災害過程以及泥石流與洪水之間的相互作用很可能會造成新的致災效應,會超過它們各自獨立相加的影響,將2種過程分開考慮通常會嚴重低估其產(chǎn)生的危害。因此,在災害鏈調(diào)查和風險評估中需要實現(xiàn)多災害模擬,而數(shù)值模擬軟件OpenLISEM在邊坡破壞、泥石流沖出預測以及洪水運動模擬等方面都表現(xiàn)出優(yōu)越的模擬效果。本文以汶川縣2019年下莊溝泥石流為對象,通過利用OpenLISEM這一多災害耦合模型來實現(xiàn)泥石流堵江災害鏈的數(shù)值模擬研究。

1 研究區(qū)概況

下莊溝溝口地形溝道狹窄呈“V”型,岸坡陡峻,兩岸坡度40°以上的面積占研究區(qū)總面積的49%,大量崩滑體松散固體物質堆積于坡腳及溝道內(nèi),在降雨作用下極易啟動轉化為泥石流。下莊溝最大高程為4 120 m,溝口地勢最低處高程為1 450 m,相對高差達2 670 m,溝床平均縱坡降為262.4‰,高陡的地形有利于快速匯水,為堆積在溝谷兩岸的崩滑體物源的侵蝕提供了充足的動力條件。

2019年8月20日,汶川縣境內(nèi)突降暴雨,導致岷江、雜谷腦河沿岸數(shù)十條溝暴發(fā)泥石流災害。其中,下莊溝于8月20日凌晨暴發(fā)泥石流,泥石流發(fā)生前期累積降雨量達27.9 mm,誘發(fā)此次泥石流的小時雨強為12.7 mm/h,由于溝口處的雜谷腦主河寬度約90 m,此次暴發(fā)泥石流沖出的固體物質達10.05萬m3,大量碎屑物質搬運出溝口后,在溝口形成大型泥石流堆積扇,堵塞雜谷腦河并形成堰塞湖,致使上游水位抬升近5 m,上游大量房屋與基礎設施被淹沒。此外,泥石流致使溝口G317國道被於埋近180 m,同時損壞下莊村房屋51間,掩埋蓉昌高速公路的多個橋墩,對上下游居民的生命安全構成巨大威脅。

2 模型理論

OpenLISEM模型中集成了滑坡、泥石流和山洪等多種災害的理論過程,主要通過輸入柵格圖層的屬性求解單元網(wǎng)格特定過程和控制流的微分方程,該模型可將水文過程和地表運動過程結合起來,在流域尺度上實現(xiàn)基于降雨事件的徑流、洪水、滑坡和泥石流運動模擬。

a)地表流。在OpenLISEM模型中,地表流是地表徑流和洪水的組合,其運輸以及動力學特征可以用描述流體的微分方程來模擬,所有的流體運動計算都基于淺表流的圣維南二維非恒定流運動方程。

(1)

(2)

(3)

式中h——流深,m;u——流速,m/s;R——降雨量,m;I——降雨入滲量,m;g——重力加速度,m/s2;S——摩擦阻力項,m/s2;Sf——摩擦動力項,m/s2;n——曼寧摩擦系數(shù)。

b)邊坡穩(wěn)定性。邊坡的穩(wěn)定性評估是基于無限邊坡模型,該方法計算了下滑力和抗滑力,其中下滑力的計算是基于破壞面平行于表面的假設。在OpenLISEM模型中可以提供動態(tài)的降雨輸入選項,邊坡破壞的觸發(fā)機制為降雨事件中濕潤鋒下滲,濕潤鋒可通過增加土體總重量和提高地下水位來影響邊坡穩(wěn)定性。若濕潤鋒下滲未達到地下水位,則只增加土體總重量。濕潤鋒對邊坡穩(wěn)定性的影響可以見式(4):

(4)

式中 SF——邊坡安全系數(shù);c′——土壤有效黏聚力,kPa;c——根系的有效黏聚力,kPa;γ——土壤單位重度,kN/m3;γs——飽和土單位重度,kN/m3;γw——水的重度,kN/m3;β——邊坡坡度,(°);φ′——有效內(nèi)摩擦角,(°);Z——土壤厚度,m;Zw——初始地下水水位,m。

c)泥石流運動方程。為了模擬洪水和泥石流的動力學特征和相互作用,OpenLISEM模型中運用了一組使用廣泛的兩相泥石流方程,該方程基于物理的二相動量守恒定律,除包含壓力和重力外,還包括黏滯力、非牛頓流體黏滯力、固相阻力和摩爾-庫侖摩擦力[14]。這種方法可以在非黏性流、高濃度流和泥石流之間實現(xiàn)平穩(wěn)過渡,解決了不同類型的流體之間的相互作用。

(5)

(6)

(7)

(8)

式中αs——固體體積分數(shù);αf——液體體積分數(shù);Pb——基面壓強,kPa;NR——雷諾數(shù);NRA——準雷諾數(shù);CDG——阻力系數(shù);ρs——固體密度,kg/m3;ρf——液體密度,kg/m3;γ——液體固體之間的密度比;x——流體垂直剪切的速度,m/s;ε——模型的高寬比;ξ——固體體積分數(shù)的垂直分布。

3 數(shù)據(jù)獲取與參數(shù)厘定

3.1 數(shù)據(jù)獲取

為探究下莊溝泥石流形成、運動以及堵江形成堰塞湖的整個災害鏈過程,收集了下莊溝高精度的地形數(shù)據(jù)、無人機正射影像、衛(wèi)星遙感影像以及下莊溝鄰近雨量站的實時降雨監(jiān)測數(shù)據(jù),以此重現(xiàn)降雨事件下泥石流堵江災害鏈的動態(tài)情境。在OpenLISEM數(shù)值模擬軟件中,主要輸入數(shù)據(jù)分為4類:高程數(shù)據(jù)(DEM)、土地植被數(shù)據(jù)、物源數(shù)據(jù)和降雨數(shù)據(jù)。

a)地形數(shù)據(jù)。本文的數(shù)值模擬基于1∶1 000等高線地形數(shù)據(jù),經(jīng)過與無人機航拍地形數(shù)據(jù)合成校正后,利用Arcgis平臺生成研究區(qū)的數(shù)字高程模型(DEM),最終得到5 m×5 m的柵格數(shù)據(jù)(圖1a)。

a)數(shù)字高程模型

b)物源數(shù)據(jù)?；诟叻直媛蔬b感影像及無人機航拍影像,對研究區(qū)進行精細的遙感解譯以此獲取泥石流物源分布,并根據(jù)后期野外現(xiàn)場測量數(shù)據(jù)進行校準,物源厚度是泥石流數(shù)值模擬中極為重要的參數(shù)之一。地震滑坡為震后泥石流主要物源,其厚度計算依據(jù)前人研究提出的滑坡平均厚度與滑坡面積之間的函數(shù)關系,這里運用式(9)計算滑坡物源厚度[14],溝道物源厚度則根據(jù)野外侵蝕斷面測量數(shù)據(jù)進行估算。

SD=1.105lnA-4.795

(9)

式中 SD——物源厚度,m;A——滑坡面積,m2。

計算后通過Arcgis平臺對研究區(qū)物源的厚度屬性進行賦值,并將物源矢量文件轉換為柵格文件,導入QGIS軟件后轉換為模擬所需的Map文件(圖1c)。

c)土地植被數(shù)據(jù)。基于遙感觀測和野外調(diào)查,下莊溝流域內(nèi)植被發(fā)育,主要以灌木為主,次為喬木以及草本植物,覆蓋率達70%以上。為反映植被覆蓋情況并準備基礎數(shù)據(jù),綜合影像質量、成像時間和影像可使用率等方面選擇了空間分辨率30 m的2020年7月27日Landsat7影像作為數(shù)據(jù)源,利用ENVI提取并估算了研究區(qū)的植被覆蓋度(圖1b)。不同的土地利用類型在災害形成過程中,往往表現(xiàn)出不同的力學特征?；趯Σ煌恋乩妙愋偷倪b感解譯分類,根據(jù)OpenLISEM使用手冊(2018版)對相應土地利用類型的多個屬性進行賦值[16](表1)。

表1 土地利用類型相關參數(shù)取值(OpenLISEM使用手冊,2018)

d)降雨數(shù)據(jù)。水文過程包括降雨、地表蓄水、截留、滲透、蒸發(fā)蒸騰和徑流,模型中利用時間和空間的降雨強度值來模擬降雨。為了還原下莊溝“8·20”泥石流的災害過程,需要輸入觸發(fā)泥石流的降雨事件的實時監(jiān)測降雨過程,包括降雨歷時與泥石流激發(fā)雨強。2019年8月20日凌晨3:00左右下莊溝暴發(fā)泥石流,本文收集了下莊溝鄰近雨量站的降雨數(shù)據(jù),并將泥石流暴發(fā)過程的雨量數(shù)據(jù)作為降雨事件輸入軟件進行模擬,見圖2。

圖2 降雨過程曲線

3.2 參數(shù)厘定

泥石流與洪水模型預測中有一個共同的難題,即輸入數(shù)據(jù)和參數(shù)的準確性。只有在輸入大量足夠精確參數(shù)的情況下,才能使用綜合方法來模擬邊坡破壞、泥石流和洪水這一系列災害過程。然而,由于近年來科學技術的發(fā)展,可獲得的精密數(shù)據(jù)越來越多,數(shù)據(jù)獲取問題已變得不那么困難。而本文為了進行更精準合理的模擬預測,針對研究區(qū)開展了多方面的現(xiàn)場災害調(diào)查以及室內(nèi)外試驗(圖3),并結合相關研究厘定模型中所需要的特征參數(shù)[15],見表2。

a)現(xiàn)場調(diào)查

表2 下莊溝OpenLISEM模擬主要參數(shù)

4 模擬結果分析與驗證

4.1 結果分析

利用對應模型和相應參數(shù)對下莊溝“8·20”泥石流運動進行模擬,分析其影響范圍及運動特征。此次模擬總時長為120 min,通過OpenLISEM模型計算,重現(xiàn)了下莊溝“8·20”泥石流啟動—搬運—侵蝕—堆積整個過程的運動特征,圖4分別是不同步長(1 step=45 s)的泥石流分布。從模擬結果可以看出,降雨事件的前期階段t=step30時,隨著降雨強度逐漸增大,松散土體達到飽和并受到徑流的強烈沖刷,數(shù)個邊坡逐漸開始失穩(wěn)運動,此刻還未啟動形成泥石流。t=step40時,大量靠近溝道的岸坡失穩(wěn)破壞并快速向地勢低洼處匯集,主要集中在下莊溝溝段下游部分,成為補給泥石流的主要來源之一。t=step50時,持續(xù)的降雨形成地表徑流,松散固體物質不斷在溝道匯集,在洪流強大的水動力作用下形成陣流,不斷向下游推進。t=step100時,由于1號滑坡失穩(wěn),形成平均厚度約7 m、寬約35 m、長約80 m的滑坡壩,后由于上游匯聚的攜沙水流的不斷沖擊,堰塞體發(fā)生潰決,流量及流速產(chǎn)生放大效應,高速向下游溝口推進。

a)step30

圖5、6表明了監(jiān)測點分布和泥石流沿程最大流深及流速的變化,距離溝口5 500 m的上游溝道處流深及流速較低,最大流速約為0.79 m/s,流深約2.62 m,這一階段處于徑流匯聚及能量積聚過程。由于重力勢能不斷轉化為動能,流速隨時間逐漸變大,但在局部彎道處,流速有所減緩,但整體上表現(xiàn)出增加的趨勢,最大達3.39 m/s。隨著動力作用的增強,泥石流流體的侵蝕作用及攜沙能力也得到大幅提高,刮鏟溝床及溝岸,使溝床普遍加寬,溝床下切深度加大,由于物源的不斷匯入,泥石流規(guī)模變大。當泥石流流體達到距溝口1 300 m處的1號滑坡堵潰點,由于堵塞作用流速急劇下降,約為0.58 m/s,泥石流深度約為8.32 m,在泥石流流體持續(xù)的沖擊作用下,滑坡壩失穩(wěn),產(chǎn)生潰決放大效應,導致流速急劇增加至4.73 m/s,沖出溝口,在溝口的最大堆積厚度為7.39 m。

圖5 泥石流動力過程監(jiān)測點位置分布

圖6 泥石流沿程最大流深及流速變化曲線

t=step120時,泥石流繼續(xù)向前推進至G317國道距離雜谷腦河道約10 m處(圖7),到達溝口堆積扇頂部,此時泥石流流深普遍在2.7～3.4 m范圍內(nèi),流速在4.2～4.7 m/s。t=step130時,泥石流進入雜谷腦河,由于雜谷腦河河床縱坡降較緩,泥石流整體推進速度降低至1.5～2.7 m/s,沖出的固體物質基本堵塞河道,形成長約146 m,寬約51 m,平均厚5.5 m的堰塞體。t=step140及t=step150時,泥石流完全堵塞主河,分別形成了長、寬、厚為236.0、74.0、6.2 m和312.0、79.0、7.1 m的堰塞體(圖7),因雜谷腦河對堆積扇的擾動以及沖刷作用,泥石流堆積扇平面形態(tài)呈偏向下游的長條形。主河完全被堵塞后,造成上游水位抬升,居民區(qū)房屋及公路被淹沒。壩后及監(jiān)測點A、B、C的洪水水深時間曲線見圖8,上游居民區(qū)淹沒深度為3～5 m,與實際調(diào)查情況較為吻合。

圖8 監(jiān)測點洪水過程曲線

a)step120

4.2 模擬結果驗證

前期野外實地測量得知,下莊溝“8·20”泥石流堆積范圍為2.17×104m2,堆積扇平均厚度5.2 m左右,沖出方量約10.5×104m3。為驗證模擬結果的準確性,以野外實際調(diào)查情況與模擬計算堆積結果對比驗證(表3),經(jīng)過對比計算,模擬堆積范圍與實測范圍重疊區(qū)域面積2.04×104m2,形態(tài)大小基本一致。同時,由式(10)計算得到模擬精度,經(jīng)過計算,模擬精度達到84.9%,符合數(shù)值模擬精度要求。

表3 模擬結果與野外實測驗證對比

(10)

式中A——泥石流數(shù)值模擬精度,%;S0——模擬與實測堆積范圍重疊區(qū)域,m2;SM——野外實測堆積范圍,m2;SN——模擬堆積范圍,m2。

5 結論

震后泥石流的形成是地震與降雨共同作用的結果,研究其運動過程和災害鏈效應,有助于實現(xiàn)泥石流災害鏈風險評價的精準性和有效性。通過對2019年8月20日下莊溝泥石流堵江災害鏈的調(diào)查,首次利用OpenLISEM重現(xiàn)了下莊溝“8·20”泥石流從啟動、搬運、侵蝕、沖出堵江到上游洪水淹沒整條災害鏈的動態(tài)過程,揭示了泥石流堵江災害鏈的形成過程與運動特征,得出以下認識。

a)作為震后泥石流堵江-洪水災害鏈的典型案例,下莊溝泥石流堵江災害鏈的形成主要是流域內(nèi)滑坡體失穩(wěn)造成溝道堵塞,滑坡壩失穩(wěn),產(chǎn)生潰決放大效應,導致流速急劇增加至4.73 m/s,沖出溝口物質達10.05萬m3。隨著泥石流形成的堆積扇漸進向河道推進,主河過流斷面逐漸縮減,造成上游的洪水淹沒危險。

b)OpenLISEM模型應用于下莊溝“8·20”泥石流災害,模擬重現(xiàn)了泥石流的形成—運動—沖出—堵江—洪水動態(tài)過程,不僅確定了泥石流的危險范圍,同時獲得了堵江后洪水淹沒深度動態(tài)圖,可以有效評估后期該區(qū)域可能存在的災害風險,為當?shù)仫L險管控提供依據(jù)。通過將模擬結果與野外實測數(shù)據(jù)對比,模擬的準確度達84.9%,表明該模型對滑坡-泥石流-洪水的災害鏈過程模擬具有較好的適用性,尤其在同震滑坡物源充分發(fā)育的汶川震區(qū)具有極大的意義。

c)災難性的地震對環(huán)境的影響是一個長期的過程,汶川震后的十余年來,群發(fā)性泥石流事件一次又一次對山區(qū)人民的生存環(huán)境造成毀滅性的打擊。在對泥石流災害事件的模擬中,如果只單獨考慮單一災害過程的模擬,忽略其他次級過程相互作用機制,這會對模擬結果的預測精度和應用造成較大影響。另外,泥石流的發(fā)生往往會伴隨洪水,將這2個過程分開考慮通常會嚴重低估其產(chǎn)生的風險及危害。