基于DAN3D的五花海泥石流運(yùn)動(dòng)過程模擬研究

盛 豪， 胡 卸 文, 2， 熊 沖 沖

(1.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 610031; 2.西南交通大學(xué)高速鐵路運(yùn)營(yíng)安全空間信息技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，四川 成都 610031)

1 概 述

2017年8月8日，九寨溝(33.20°N，103.82°E)發(fā)生M7.0級(jí)強(qiáng)烈地震，地震誘發(fā)境內(nèi)多處崩滑、滑坡、泥石流等次生地質(zhì)災(zāi)害，嚴(yán)重破壞了九寨溝世界自然遺產(chǎn)的景觀面貌[1-2]。五花海泥石流位于九寨溝核心景區(qū)，受地震影響，溝內(nèi)發(fā)育多處崩滑、震裂物源，松散堆積體大大增加，為泥石流發(fā)育提供了充足物源條件。在暴雨條件下，泥石流一旦沖出，會(huì)對(duì)溝口五花海景觀造成嚴(yán)重破壞。

國(guó)內(nèi)外對(duì)于震后泥石流的研究很多，取得了較為豐碩的成果。Alexander J. Horton[3]利用Massflow對(duì)汶川震后紅椿溝進(jìn)行了運(yùn)動(dòng)過程反演后認(rèn)為，利用voellmy鏟刮模型能夠較好適用于震后泥石流的模擬；胡旭東[4]選取了九寨溝震后五個(gè)集水區(qū)，綜合采用泥沙增加量計(jì)算、粒子模型及數(shù)值計(jì)算的方法，對(duì)九寨溝震后泥石流危害性進(jìn)行了評(píng)價(jià)；胡卸文[5-6]針對(duì)汶川震后泥石流沖出規(guī)模遠(yuǎn)大于規(guī)范計(jì)算結(jié)果，提出了震后泥石流普遍存在的一類新物源—震裂物源，并給出其啟動(dòng)模式；張向營(yíng)[7]對(duì)地震擾動(dòng)區(qū)泥石流啟動(dòng)閾值進(jìn)行了探討，認(rèn)為對(duì)比震后泥石流啟動(dòng)激發(fā)雨量遠(yuǎn)遠(yuǎn)低于震前，僅為震前的三分之一左右；顧文韜[8]在調(diào)查統(tǒng)計(jì)了汶川震后多條泥石流物源基礎(chǔ)上，總結(jié)了極震區(qū)震后泥石流形成過程及啟動(dòng)方式，得出了震后泥石流動(dòng)儲(chǔ)量與流域面積、高差的計(jì)算關(guān)系，為震后泥石流動(dòng)儲(chǔ)量估算及工程治理提供思路。綜上所述，震后泥石流區(qū)別于常規(guī)泥石流，其具備臨界雨量小、物源儲(chǔ)量多、一次性沖出規(guī)模大等特點(diǎn)。因而，對(duì)其深入研究能夠?yàn)橹卫砉こ淘O(shè)置及災(zāi)后恢復(fù)提供新的思路及方法。

筆者在對(duì)五花海泥石流野外調(diào)查的基礎(chǔ)上，獲取了地形地貌、地層巖性、物源條件及水源條件等相關(guān)形成條件及發(fā)育特征，計(jì)算了泥石流流速和一次沖出規(guī)模等動(dòng)力參數(shù)。同時(shí)，運(yùn)用DAN3D軟件對(duì)震后泥石流運(yùn)動(dòng)過程進(jìn)行模擬研究，對(duì)比分析兩種結(jié)果，得出公式計(jì)算結(jié)果普遍偏小的結(jié)論，為防治措施提供了參考意見。

2 五花海泥石流發(fā)育特征

2.1 地形條件

五花海泥石流溝域形態(tài)呈櫟葉形，無(wú)支溝發(fā)育，流域溝道縱向長(zhǎng)度3.5 km，平均寬度0.85 km，流域面積2.6 km2，平均縱坡降434‰。流域內(nèi)水系沿主溝道分布，最高點(diǎn)位于流域西側(cè)山脊部位，高程3 981 m，最低點(diǎn)位于五花海，高程2 460 m，相對(duì)高差1 521 m，具有較顯著的清水區(qū)、形成流通區(qū)和堆積區(qū)劃分,見圖1。

清水區(qū)平均縱坡475‰，溝谷呈深“V”字形，兩側(cè)縱坡陡峻，植被發(fā)育良好，少見崩塌滑坡等不良地質(zhì)現(xiàn)象；形成流通區(qū)平均縱坡412‰，兩側(cè)岸坡不良地質(zhì)現(xiàn)象發(fā)育，溝道內(nèi)堆積震后垮塌滾落塊石，受地震影響，樹木多被折斷；堆積區(qū)平均縱坡170‰，溝形不明顯，植被覆蓋良好，未見新近泥石流堆積。

2.2 物源條件

經(jīng)野外調(diào)查，五花海泥石流溝域在地震前兩側(cè)斜坡植被發(fā)育較好，水土流失輕微，基本未見斜坡變形破壞現(xiàn)象。九寨溝“8·8”地震之后，五花海泥石流新增大量松散固體物源，溝內(nèi)物源類型主要為崩滑物源、溝道堆積物源、坡頂震裂物源及坡面物源。震裂物源主要分布于溝道上游兩側(cè)山體頂部，主要為“8·8”地震產(chǎn)生。崩滑物源、溝道堆積物源主要分布于溝道中游區(qū)域溝道內(nèi)和溝道右岸斜坡中下部。坡面物源大部分分布于谷坡中上部，儲(chǔ)量一般。物源分布情況見圖2。此外，在溝道中上游兩側(cè)坡面存在大量震倒樹木，可能成為潛在漂木物源，一旦沖出，會(huì)對(duì)下游防護(hù)措施產(chǎn)生巨大沖擊。據(jù)現(xiàn)場(chǎng)調(diào)查統(tǒng)計(jì)，崩滑物源總量約為31.2×104m3，溝道堆積物源總量約為2.5×104m3，坡面物源總量約為2.1×104m3，震裂物源總量約為25.4×104m3，共計(jì)松散固體物源量61.2×104m3，可能參與泥石流活動(dòng)的物源量7.34×104m3(圖3)。

圖1 五花海泥石流主溝縱剖面

2.3 水源條件

五花海泥石流溝的水源主要來(lái)源于大氣降水。鑒于既往研究，九寨溝區(qū)域泥石流均發(fā)生于雨季，春季冰雪融水一般不會(huì)成為引發(fā)泥石流的水源。此外，溝域內(nèi)中上游地下水不豐富，不構(gòu)成引發(fā)泥石流的主要水源，且溝域內(nèi)沒有水庫(kù)、湖泊等集中的地表水體。因而，暴雨形成的地表徑流是引發(fā)泥石流的主要水源，暴雨是泥石流的主要激發(fā)因素。泥石流溝域整體形態(tài)呈櫟葉形，清水區(qū)溝道縱坡大，岸坡陡峻，有利于地表降水匯集，為五花海泥石流的發(fā)育提供有利的水源條件。

3 泥石流運(yùn)動(dòng)特征參數(shù)計(jì)算

泥石流運(yùn)動(dòng)特征參數(shù)是數(shù)值模擬求解設(shè)置和分析對(duì)比的依據(jù)，所需確定參數(shù)主要為泥石流重度、流速及暴雨洪峰流量、泥石流峰值流量，可采用雨洪修正法計(jì)算獲得。

(1)泥石流重度?，F(xiàn)場(chǎng)配漿法與查表法所得泥

圖2 五花海泥石流平面圖

圖3 五花海泥石流物源統(tǒng)計(jì)直方圖

石流重度綜合取值15.1 kN/m3，屬于稀性泥石流。

(2)泥石流流速。五花海泥石流屬稀性泥石流，流速計(jì)算公式采用西南地區(qū)(鐵二院陳光曦)公式。

(1)

式中Vc為泥石流流速；H為泥石流平均泥深；I為泥位縱坡降； 1/n為溝床糙率系數(shù)。

溝口附近泥石流流速計(jì)算結(jié)果見表1。

表1 五花海泥石流溝口流速計(jì)算結(jié)果

(3)暴雨洪峰流量和泥石流峰值流量。泥石流峰值流量采用雨洪修正法進(jìn)行計(jì)算。

QC=Dc(1+φ)Qp

(2)

Qp=0.278ψAS/tn

(3)

式中QC為泥石流斷面峰值流量(m3/s)；φ為泥沙修正系數(shù)；Qp為暴雨洪峰流量；Dc為堵塞系數(shù)；ψ為洪峰徑流系數(shù)；S為單位歷時(shí)的暴雨平均強(qiáng)度(mm/h)；A為流域面積(km2)；t為匯流時(shí)間(h)；n為暴雨遞減指數(shù)。

溝口泥石流峰值流量計(jì)算結(jié)果見表2。

表2 五花海泥石流溝口流量計(jì)算結(jié)果

4 DAN3D數(shù)值模擬

4.1 DAN3D軟件

DAN(Dynamic Analysis)方法由Hungr(1995)[9]提出，最早為解決高速遠(yuǎn)程滑坡運(yùn)動(dòng)過程復(fù)雜、組成物質(zhì)多樣、傳統(tǒng)理論模型無(wú)法很好適應(yīng)等問題。該方法的核心思想是等效流體假設(shè)與連續(xù)性假設(shè)。DAN模型將運(yùn)動(dòng)中的復(fù)雜非均值的滑體或泥石流漿體視為某種理想的、具有簡(jiǎn)單的流變關(guān)系的流體，并將流體視為連續(xù)性介質(zhì)，以連續(xù)的場(chǎng)代替小范圍內(nèi)離散物質(zhì)相互作用。因而，物質(zhì)的變形及運(yùn)動(dòng)僅取決于內(nèi)部及底部的流變準(zhǔn)則。目前國(guó)內(nèi)利用DAN3D對(duì)高速遠(yuǎn)程滑坡、碎屑流等災(zāi)害進(jìn)行模擬反演已經(jīng)取得了比較好的成果，如深圳“12.20”渣場(chǎng)滑坡、湯家溝滑坡～碎屑流等案例[10-11]。國(guó)外采用DAN3D對(duì)泥石流進(jìn)行模擬的成功案例比較多，如皮茨河谷地區(qū)泥石流、上阿諾谷泥石流[12-14]等案例。

4.2 Voellmy模型的基本原理

DAN3D內(nèi)置Frictional、Plastic、Newtonian、Bingham及Voellmy五種流變模型，使用者可根據(jù)實(shí)際問題選擇恰當(dāng)?shù)哪Ｐ瓦M(jìn)行數(shù)值模擬分析。多項(xiàng)研究結(jié)果表明，對(duì)泥石流進(jìn)行動(dòng)力特性分析的時(shí)候，Bingham及Voellmy模型能夠較好反映泥石流運(yùn)動(dòng)歷史。由于五花海泥石流為稀性泥石流，Voellmy模型更能考慮到流體運(yùn)動(dòng)過程中的湍流效應(yīng)與速度效應(yīng)。因而，本文采用Voellmy模型進(jìn)行模擬。

Voellmy模型認(rèn)為，流體在運(yùn)動(dòng)過程中所受到的抗剪應(yīng)力應(yīng)為流體所受摩擦力和湍流流動(dòng)共同產(chǎn)生的阻力之和，其模型公式如下所示：

式中τ為滑體所受抗剪應(yīng)力；f為摩擦項(xiàng)，等同于tanφ；ξ為紊流項(xiàng)，與流體速度和密度相關(guān)。

5 五花海泥石流數(shù)值模擬

5.1 參數(shù)設(shè)置

利用五花海泥石流1∶1000地形圖進(jìn)行溝域DEM三維地形建模，并采用Kriging插值法進(jìn)行網(wǎng)格化處理，生成溝域網(wǎng)格模型。同時(shí)，根據(jù)現(xiàn)場(chǎng)調(diào)查情況，將物源厚度建立與溝域網(wǎng)格模型相對(duì)應(yīng)的網(wǎng)格模型，生成的DEM模型見圖4、5。

圖4 五花海泥石流溝域DEM模型

本文采用Voellmy摩擦模型作為五花海泥石流流體運(yùn)動(dòng)準(zhǔn)則，其中最重要的參數(shù)為摩擦系數(shù)μ、湍流系數(shù)ξ及泥石流流體容重。摩擦系數(shù)及湍流系數(shù)的取值范圍，在前期的許多文獻(xiàn)中有過研究。一般認(rèn)為，Voellmy模型中摩擦系數(shù)為0.05～0.8，湍流系數(shù)30～800 m/s2。因此，本文在前人研究基礎(chǔ)上選取不同參數(shù)進(jìn)行模擬，最終確定摩擦系數(shù)選用0.3，湍流系數(shù)選用50 m/s2。

圖5 五花海泥石流物源DEM模型

5.2 模擬結(jié)果分析

根據(jù)上述設(shè)定參數(shù)，利用DAN3D對(duì)五花海泥石流在P=2%及P=1%暴雨時(shí)的運(yùn)動(dòng)情況進(jìn)行數(shù)值模擬，得到泥石流堆積厚度模擬結(jié)果，見圖6。

圖6的模擬結(jié)果表明：50 a一遇泥石流發(fā)育過程中，崩滑物源首先啟動(dòng)，向坡腳溝道堆積；位于分水嶺的震裂物源漸次運(yùn)動(dòng)，進(jìn)入溝道或沿坡面向下運(yùn)動(dòng)，形成溝道物源及坡面物源，參與到泥石流活動(dòng)中。從t=30～150 s這段時(shí)間的云圖可以看出，泥石流運(yùn)動(dòng)到溝道寬緩(海拔2 675～2 710 m段)位置時(shí)，流體逐漸停止運(yùn)動(dòng)，在溝道內(nèi)堆積，僅部分能夠沖出繼續(xù)向下游流動(dòng)。t=300 s的云圖可以看到，溝口有泥石流堆積，但是方量很小。t=300～600 s云圖可以看出，這段時(shí)間泥石流沖出量極少，小于公式計(jì)算結(jié)果，前緣漫流至五花海棧道，對(duì)棧道及游客安全構(gòu)成一定威脅。同時(shí)，從圖8可以看出：泥石流在溝口最大速度也小于公式計(jì)算結(jié)果。因此，可以推測(cè)，海拔2 675～2 710 m處的寬緩平臺(tái)為泥石流起到了天然停淤作用，很好地減緩了泥石流沖出速度，滯留了大量固體物質(zhì)。

圖7的模擬結(jié)果表明：100年一遇泥石流運(yùn)動(dòng)及物源啟動(dòng)方式與50年一遇的情況相一致。但是由于泥石流流經(jīng)寬緩溝道時(shí)流速較大，因而，大量泥石流進(jìn)一步向下游沖出，逐漸堆積形成扇形堆積扇。根據(jù)圖9可以看出，泥石流在溝口最大速度約為6.5 m/s，與公式計(jì)算結(jié)果符合較好。泥石流一次沖出總量為5.7×104 m3，遠(yuǎn)大于公式計(jì)算結(jié)果，與震后泥石流沖出規(guī)模遠(yuǎn)遠(yuǎn)大于常規(guī)泥石流的認(rèn)知一致。

圖6 五花海泥石流P=2%堆積厚度云圖

圖7 五花海泥石流P=1%堆積厚度云圖

5.3 影響范圍預(yù)測(cè)及防治建議

五花海泥石流溝口正對(duì)五花海及棧道，因此，利用模擬結(jié)果對(duì)泥石流沖出危害進(jìn)行一定的預(yù)測(cè)，在后期對(duì)五花海泥石流提出針對(duì)性的建議很有必要。50年一遇工況下，雖然物源啟動(dòng)量很大，但大部分在流經(jīng)寬緩溝道處便開始減速，并逐步停止運(yùn)動(dòng)，僅少量沖出溝口，漫流至棧道處。此時(shí)，泥石流流速不大，最大約3 m/s，對(duì)棧道會(huì)構(gòu)成一定威脅，但對(duì)五花海景觀影響不大；100年一遇工況下，盡管有部分泥石流在寬緩溝道處停淤，但仍有大量泥石流沖出溝口。堆積扇前緣延伸至五花海，溝口最大速度可達(dá)6.5 m/s，可能沖毀棧道。同時(shí)，大量泥沙可隨之裹挾進(jìn)入五花海，對(duì)五花海景觀造成毀滅性打擊。

圖8 P=2%泥石流最大速度云圖

圖9 P=1%泥石流最大速度云圖

根據(jù)模擬結(jié)果可知：五花海泥石流爆發(fā)主要物源來(lái)源于崩滑物源，可采取支擋措施避免崩滑物源垮塌進(jìn)入溝道參與泥石流活動(dòng)；同時(shí)，泥石流在寬緩溝道處會(huì)大量堆積，可在此設(shè)置停淤圍堰，提高其停淤能力；為避免泥沙進(jìn)入五花海，對(duì)五花海水質(zhì)景觀造成破壞，可在臨近溝口處設(shè)置攔水堤，利用天然土體對(duì)含砂水流進(jìn)行過濾。

6 結(jié) 語(yǔ)

五花海泥石流位于九寨溝核心景區(qū)境內(nèi)，溝域形態(tài)呈櫟葉形，流域面積為2.6 km2，主溝長(zhǎng)0.85 km，主溝平均縱坡降434‰?！?·8”地震后，溝域內(nèi)地質(zhì)災(zāi)害發(fā)育松散物源豐富，為泥石流的形成提供了地形、物源和水源條件。

結(jié)合泥石流基本發(fā)育特征和動(dòng)力學(xué)參數(shù)計(jì)算，可得出九寨溝泥石流容重15.1 kN/m3，50 a、100 a一遇降雨頻率下溝口流速分別5.68 m/s、6.21 m/s，一次過流總量為0.92 m3、1.11 m3。

通過DAN3D軟件對(duì)泥石流進(jìn)行數(shù)值模擬，預(yù)測(cè)了泥石流50 a與100 a一遇的運(yùn)動(dòng)過程，可以得知50 a一遇工況下泥石流不會(huì)大量沖出，少量沖出泥石流會(huì)對(duì)棧道構(gòu)成一定威脅。100 a一遇工況下，泥石流漿體大量涌出，泥沙大量沖入五花海，會(huì)對(duì)五花海景觀水質(zhì)構(gòu)成嚴(yán)重威脅。

根據(jù)模擬結(jié)果，提出了“攔擋物源、圍堰停淤、攔水過濾”的防治方案。