基于Massflow模型的青龍溝臺風暴雨型泥石流運動特征研究

熊朝正,吉 鋒,石豫川

(1.四川藏區高速公路有限責任公司,四川 成都 610041;2.成都理工大學 地質災害防治與地質環境保護國家重點實驗室,四川 成都 610059)

近年來隨著全球氣候變化,臺風對中國東南沿海地區造成了嚴重影響[1-2]。伴隨臺風而產生的暴雨是一種復雜多變的降雨類型,具有雨量大、降雨集中,極端雨強大、雨型單峰變化、持續時間短等特點[3-7]。由臺風暴雨直接誘發的臺風暴雨型泥石流是一類特殊的災害類型,此類泥石流憑借其顯著的突發性、群發性和破壞性,愈發引起人們的重視。開展臺風暴雨型泥石流運動特征研究,具有十分重要的理論和現實意義[8]。

泥石流作為一種非均勻混合介質,具有復雜的物理過程和動力學特征[9],隨著計算機技術、數值算法、本構模型等方面的發展成熟,數值模擬已成為現在泥石流運動特征研究的主要方法,其中應用較多的軟件主要有FLO-2D、DAN、MassMov2D、RAMMS、EDDA、Massflow等。唐得勝[10]將FLO-2D應用于都江堰龍溪河流域泥石流溝數值模擬,結果顯示重現期為20年一遇時的模擬結果與“8.13特大泥石流”對比精確度在60%～87%。Cheon等[11]將DAN-3D軟件與機器學習方法相結合,對韓國Umyeon山泥石流進行泥石流寬度、體積及沖擊力的計算,在生成危險性評估圖的基礎上選擇最佳攔擋設施放置位置。Beguería[12]應用MassMov2D模型分別對奧地利蒂羅爾和法國阿爾卑斯山2個泥石流進行反分析,模擬泥石流堆積物的最終分布和厚度。宋兵等[13]探討了RAMMS在泥石流模擬中的運用方法,通過模擬北川縣鄉白沙溝泥石流,得出泥石流沖出量與雨洪法計算量(20年一遇)誤差不超過5%,可信度較高。Chen等[14]提出一種三維綜合數值模型EDDA,該模型考慮了泥石流過程中密度、屈服應力和動力黏度的變化,驗證了該模型在流域尺度模擬中的性能。

Massflow基于深度積分的連續介質力學方法[15],能揭示滑坡、泥石流、碎屑流、山洪等山地災害的時空演化全過程,把三維計算問題簡化為二維,有效地提高了計算效率。丁邦政[16]運用Massflow對尾礦庫潰壩泥石流進行了數值模擬,計算泥石流最大沖起高度、沖擊力、流深、流速等數據,在此基礎上對橋梁進行泥石流沖擊作用響應分析和風險評估。段學良等[17]利用無人機航拍三維建模生成DEM,結合Massflow軟件,比選后采用更精準的Manning摩擦模型對西藏仁布杰仲溝泥石流進行災害模擬評估,將研究區劃分為4個危險區,為泥石流防治提供了重要依據。

本文運用Massflow軟件進行臺風暴雨型泥石流運動特征數值模擬,對比分析泥石流的運動堆積過程、泛濫范圍及堆積厚度等動力學特征。研究成果不僅能給臺風影響地區的地質災害防治提供參考,同時也能為今后臺風誘發泥石流的深入研究提供計算方法及手段參考。

1 研究區概況

研究區擬建抽水蓄能電站由下水庫(壩)、上水庫(壩)、地下廠房、輸水系統等主要建筑物組成。青龍溝位于上水庫區,其中部為上水庫蓄水區,擬建水工大壩橫跨青龍溝將其截斷,為更好地評價青龍溝對庫區的影響,以壩址為界,將青龍溝分為上下兩段。青龍溝流域面積約1.84 km2,長度約2.86 km,高差為580 m,整體縱坡降為203‰,流域內植被發育,植被覆蓋率大于60 %,其中以喬木、竹為主。青龍溝地形總體上屬于構造剝蝕低山地貌,兩側支溝和小型的沖溝很發育。據野外地質調查,調查段共發育有13條主要的一級支溝(Qs1—13)及若干的小型沖溝(圖1),兩岸地形總體較陡,坡度多為30～50°。通過調查青龍溝物源點有13處,共有松散物源總量約185.34萬m3,不穩定物源量約17.56萬m3。因此青龍溝具有對泥石流暴發有利的地形和物源條件。

圖1 青龍溝及支溝分布

2019年8月10日在“利奇馬”引發的臺風暴雨作用下,當日下午16時左右(即在臺風暴雨短歷時超強降雨發生時刻),啟動青龍溝暴發了大規模泥石流。通過訪問居住在青龍溝附近的居民,得知青龍溝泥石流暴發持續時長約為15 min。由于泥石流流量過大,無法沿原始溝道正常運動排泄,最終沖出原始溝道,致使溝道改變,在左側形成新的泥石流溝道,并大量擴散堆積(圖2),青龍溝泥石流在上水庫區造成了嚴重的破壞,摧毀部分居民房屋,中斷工程進度。通過收集研究區泥石流暴發過程實際臺風暴雨資料,結合文獻分析[18-20],臺風暴雨具有突出的短歷時強降雨特征,能在短強降雨時刻引發泥石流,青龍溝泥石流的暴發也符合這一特征。臺風對泥石流的影響表現于通過樹木松動表層物源,由于臺風作用難以定量,該耦合過程較為復雜,且主要以臺風暴雨作用為主,本文主要考慮臺風暴雨作用下泥石流運動特征。

圖2 青龍溝溝道改變

2 Massflow數值模擬工作

針對青龍溝泥石流的運動堆積特征,采用中科院山地所研究開發的一款地表動力過程數值模擬軟件Massflow[21-22],來進行臺風暴雨型泥石流運動特征研究,能取得良好的效果。該軟件以深度積分連續介質力學方法為理論基礎,運用MacCormack-TVD有限差分法來求解以上方程[23-24],采用Fortran和C#編程語言,可跨越Windows和Linux平臺,并結合MPICH分布式和OpenMP共享內存并行計算技術,有效地提高了計算效率,減少了計算時間。Massflow通過對大區域切割和重點區域網格細化的方式提高了計算規模;此外還提供了修改物理模型及微分方程的源代碼,支持用戶通過二次開發編入自定義模型[25]。

2.1 地形數據處理

本次研究收集到青龍溝實測1∶2 000地形等高線和高精度遙感影像圖。通過ArcGIS軟件將原始高程數據進行轉化,獲取青龍溝DEM,本文根據實際情況建立4 m×4 m的計算網格數據。進一步運用ArcGIS軟件轉換工具中柵格轉ASCII功能,將上一步獲得的DEM數據轉換為ASCII格式(txt文件)的高程文件,即可導入Massflow軟件進行模擬計算。計算模型見圖3。

圖3 青龍溝程序計算模型

2.2 泥石流特征值設定

本文采用基于泥石流溝易發程度數量化評分的泥石流容重確定方法,確定青龍溝泥石流重度為16.34 kN/m3,相應1+φ取值1.637(φ為泥石流泥沙修正系數)。擬設計非臺風暴雨(工況1)及臺風暴雨(工況2)兩種工況進行對比模擬分析。如前文所述,研究區實際臺風降雨過程具有突出的臺風暴雨特征(即短歷時強降雨特征),是青龍溝泥石流的暴發的關鍵因素,因此,工況2采用實際降雨資料進行計算模擬。同時收集研究區往年雨季非臺風暴雨資料與臺風暴雨相比較,峰值特征不突出,過程雨量分布較為均勻,因此工況1將降雨過程概化為雨強均勻分布,根據青龍溝泥石流防治工程體系50年一遇的設計標準,計算工況1的24 h暴雨量并24 h平均,得到設計雨強為10.5 mm/h。2種工況的降雨過程見圖4。

圖4 模擬降雨工況設計

根據研究區泥石流溝谷實際匯流情況,運用實際降雨資料及當地洪水計算手冊,計算概化泥石流流量曲線,結合已有學者的研究,本文采取“概化五邊形模型”來進行相關計算[26]。根據青龍溝泥石流匯流特征,分別設定泥石流啟動點1(Q1)和啟動點2(Q2),流量過程曲線的持續時長為900 s。各啟動點在不同工況下的泥石流“概化五邊形”流量曲線計算結果見圖5、6。

圖5 啟動點1泥石流流量曲線

圖6 啟動點2泥石流流量曲線

2.3 基底摩擦模型及參數

結合本文泥石流災害特點,參考已有學者的研究成果,選取Voellmy模型作為模擬的基底摩擦模型。結合青龍溝8月10日泥石流暴發后的實際泥位線及堆積范圍,進行大量反演工作,最終確定基底摩擦系數取值0.11,湍流系數取值260 m/s2。

3 泥石流運動特征對比分析

工況1條件下模擬結果顯示,庫區內的模擬堆積范圍與實際堆積范圍基本吻合,僅在北西側邊緣略大于實際堆積范圍(圖7),分析認為是由于該處居民修建住房時對局部地形進行了一定的堆積抬升,導致實際泥石流的泛濫受到一定限制。在青龍溝下段溝道調查處模擬最大泥深為173 cm,與實際泥位線高度170 cm僅相差3 cm(圖8),整體上模擬結果準確性較高。

圖7 堆積范圍對比

圖8 調查處泥位線高度

根據青龍溝泥石流在非臺風暴雨型及臺風暴雨型降雨工況下運動特征模擬結果,泥石流運動方向及整體運動過程基本一致。啟動點1處泥石流首先向庫區中心運動,隨后與啟動點2處泥石流匯合。交匯后的泥石流都表現出疊加增強效果,加快了泥石流在青龍溝下段溝道中的前進速度。當運動至溝口后開闊的地形使得泥石流開始快速停留堆積,均形成了典型的扇狀堆積區。2種工況條件下,泥石流的最終堆積區也基本一致,主要在上庫庫區、青龍溝下游溝道、青龍溝溝口3處,分別為堆積區1、堆積區2和堆積區3(圖9、10)。

圖9 工況1條件下泥石流流深分布

圖10 工況2條件下泥石流流深分布

在非臺風暴雨工況下,泥石流分別在700、200 s到達監測點j1、j2,約在1 700 s時兩處泥石流匯流,匯流后的泥石流在1 800 s運動至j3,2 000 s時運動至j4,2 200 s時到達j5并在溝口停留堆積。在臺風暴雨型降雨工況下,泥石流分別在600、100 s到達監測點j1、j2,約在1 400 s時兩處泥石流匯流,匯流后的泥石流在1 500 s運動至j3,1 600 s時運動至j4,1 900 s時到達j5并在溝口停留堆積。泥石流到達各監測點的時間分布見圖11,由此可見臺風暴雨型泥石流暴發后,在更少的時間內到達沿途監測點位置,能更快對沿途溝道及建筑物造成破壞,居民的逃生時機更短,這也是臺風暴雨型泥石流通常會造成更大生命財產損失的重要原因之一。

圖11 到達監測點時間對比

取2種工況下流深區別較大的j1、j2、j5進行統計,見圖12,各監測點在相同時刻泥石流流深以及最終流深,臺風暴雨型泥石流均大于非臺風暴雨型泥石流。在監測點j5表現最為明顯,2 500 s時工況1最終流深5.86 m,工況2最終流深7.47 m。

圖12 監測點流深對比

圖13為兩工況條件下各堆積區最大流深統計,可見工況2條件下,3個主要堆積區的泥石流流深均大于工況1,在堆積區3最為明顯,工況1泥石流流深為6.06 m,工況2為7.75 m,約為前者1.3倍,表明臺風暴雨型泥石流具有更大垂直影響高度的特征。圖14為兩工況條件下各堆積區泛濫面積統計,在堆積區1,工況2略大于工況1,堆積區2泛濫范圍基本一致,最主要差別在溝口堆積區3,工況1泛濫范圍為18 824.9 m2,工況2為32 279.5 m2,約為前者1.7倍。結合兩工況條件下最終泛濫范圍對比(圖15)不難看出,臺風暴雨型泥石流在溝口堆積區的影響范圍明顯大于非臺風暴雨型泥石流,具有更大平面影響范圍的特征。

圖13 堆積區最大流深

圖14 堆積區面積

綜上分析,臺風暴雨型泥石流在整個運動過程中具有更快的前進速度,龍頭能在更短的時間內對沿途溝道及建筑物造成影響和破壞,此外,其泛濫過程中的垂直影響高度和平面影響范圍也都大于非臺風暴雨型泥石流,說明臺風暴雨型泥石流具備更大的泥石流物質量,其攜帶能力及破壞能力都遠強于非臺風暴雨型泥石流,在沒有防治措施的情況下對工程建筑安全及居民的生命財產安全具有很大的威脅能力。

4 結論

a)臺風暴雨型泥石流及非臺風暴雨型泥石流運動方向及整體運動過程基本一致。啟動點1與2處泥石流匯合,表現出疊加增強效果,運動至溝口后快速停留堆積,形成上庫庫區、青龍溝下游溝道、青龍溝溝口3處主要堆積區。

b)按50年一遇降雨強度設計的非臺風暴雨條件下,泥石流分別在700、200、1 800、2 000、2 200 s運動至監測點j1—j5,堆積區最大堆積厚度分別為2.65、3.77 、6.06 m,堆積面積分別為33 378.4、3 349.8、18 824.9 m2。

c)在臺風暴雨型降雨條件下,泥石流分別在600、100、1 500、1 600、1 900 s運動至監測點j1—j5,堆積區最大堆積厚度分別為3.08、3.80、7.75 m,堆積面積分別為36 324.4、3 700.4、32 279.5 m2。

d)臺風暴雨型泥石流具有很快的前進速度,能在短時間破壞沿途溝道及建筑物,居民的逃生時機較短,同時還具有高度影響范圍和平面影響范圍較大的特征,破壞范圍廣,因此對于臺風暴雨型泥石流更需注重防治措施的修建。