基于DAN3D的五花海泥石流運動過程模擬研究

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

(1.西南交通大學 地球科學與環境工程學院，四川 成都 610031; 2.西南交通大學高速鐵路運營安全空間信息技術國家地方聯合工程實驗室，四川 成都 610031)

1 概 述

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

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

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

2 五花海泥石流發育特征

2.1 地形條件

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

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

2.2 物源條件

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

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

2.3 水源條件

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

3 泥石流運動特征參數計算

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

(1)泥石流重度。現場配漿法與查表法所得泥

圖2 五花海泥石流平面圖

圖3 五花海泥石流物源統計直方圖

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

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

(1)

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

溝口附近泥石流流速計算結果見表1。

表1 五花海泥石流溝口流速計算結果

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

QC=Dc(1+φ)Qp

(2)

Qp=0.278ψAS/tn

(3)

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

溝口泥石流峰值流量計算結果見表2。

表2 五花海泥石流溝口流量計算結果

4 DAN3D數值模擬

4.1 DAN3D軟件

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

4.2 Voellmy模型的基本原理

DAN3D內置Frictional、Plastic、Newtonian、Bingham及Voellmy五種流變模型，使用者可根據實際問題選擇恰當的模型進行數值模擬分析。多項研究結果表明，對泥石流進行動力特性分析的時候，Bingham及Voellmy模型能夠較好反映泥石流運動歷史。由于五花海泥石流為稀性泥石流，Voellmy模型更能考慮到流體運動過程中的湍流效應與速度效應。因而，本文采用Voellmy模型進行模擬。

Voellmy模型認為，流體在運動過程中所受到的抗剪應力應為流體所受摩擦力和湍流流動共同產生的阻力之和，其模型公式如下所示：

式中τ為滑體所受抗剪應力；f為摩擦項，等同于tanφ；ξ為紊流項，與流體速度和密度相關。

5 五花海泥石流數值模擬

5.1 參數設置

利用五花海泥石流1∶1000地形圖進行溝域DEM三維地形建模，并采用Kriging插值法進行網格化處理，生成溝域網格模型。同時，根據現場調查情況，將物源厚度建立與溝域網格模型相對應的網格模型，生成的DEM模型見圖4、5。

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

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

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

5.2 模擬結果分析

根據上述設定參數，利用DAN3D對五花海泥石流在P=2%及P=1%暴雨時的運動情況進行數值模擬，得到泥石流堆積厚度模擬結果，見圖6。

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

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

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

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

5.3 影響范圍預測及防治建議

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

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

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

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

6 結 語

五花海泥石流位于九寨溝核心景區境內，溝域形態呈櫟葉形，流域面積為2.6 km2，主溝長0.85 km，主溝平均縱坡降434‰。“8·8”地震后，溝域內地質災害發育松散物源豐富，為泥石流的形成提供了地形、物源和水源條件。

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

通過DAN3D軟件對泥石流進行數值模擬，預測了泥石流50 a與100 a一遇的運動過程，可以得知50 a一遇工況下泥石流不會大量沖出，少量沖出泥石流會對棧道構成一定威脅。100 a一遇工況下，泥石流漿體大量涌出，泥沙大量沖入五花海，會對五花海景觀水質構成嚴重威脅。

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