熊朝正,吉 鋒,石豫川

(1.四川藏區(qū)高速公路有限責(zé)任公司,四川 成都 610041;2.成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,四川 成都 610059)

近年來(lái)隨著全球氣候變化,臺(tái)風(fēng)對(duì)中國(guó)東南沿海地區(qū)造成了嚴(yán)重影響[1-2]。伴隨臺(tái)風(fēng)而產(chǎn)生的暴雨是一種復(fù)雜多變的降雨類型,具有雨量大、降雨集中,極端雨強(qiáng)大、雨型單峰變化、持續(xù)時(shí)間短等特點(diǎn)[3-7]。由臺(tái)風(fēng)暴雨直接誘發(fā)的臺(tái)風(fēng)暴雨型泥石流是一類特殊的災(zāi)害類型,此類泥石流憑借其顯著的突發(fā)性、群發(fā)性和破壞性,愈發(fā)引起人們的重視。開(kāi)展臺(tái)風(fēng)暴雨型泥石流運(yùn)動(dòng)特征研究,具有十分重要的理論和現(xiàn)實(shí)意義[8]。

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

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

本文運(yùn)用Massflow軟件進(jìn)行臺(tái)風(fēng)暴雨型泥石流運(yùn)動(dòng)特征數(shù)值模擬,對(duì)比分析泥石流的運(yùn)動(dòng)堆積過(guò)程、泛濫范圍及堆積厚度等動(dòng)力學(xué)特征。研究成果不僅能給臺(tái)風(fēng)影響地區(qū)的地質(zhì)災(zāi)害防治提供參考,同時(shí)也能為今后臺(tái)風(fēng)誘發(fā)泥石流的深入研究提供計(jì)算方法及手段參考。

1 研究區(qū)概況

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

圖1 青龍溝及支溝分布

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

圖2 青龍溝溝道改變

2 Massflow數(shù)值模擬工作

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

2.1 地形數(shù)據(jù)處理

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

圖3 青龍溝程序計(jì)算模型

2.2 泥石流特征值設(shè)定

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

圖4 模擬降雨工況設(shè)計(jì)

根據(jù)研究區(qū)泥石流溝谷實(shí)際匯流情況,運(yùn)用實(shí)際降雨資料及當(dāng)?shù)睾樗?jì)算手冊(cè),計(jì)算概化泥石流流量曲線,結(jié)合已有學(xué)者的研究,本文采取“概化五邊形模型”來(lái)進(jìn)行相關(guān)計(jì)算[26]。根據(jù)青龍溝泥石流匯流特征,分別設(shè)定泥石流啟動(dòng)點(diǎn)1(Q1)和啟動(dòng)點(diǎn)2(Q2),流量過(guò)程曲線的持續(xù)時(shí)長(zhǎng)為900 s。各啟動(dòng)點(diǎn)在不同工況下的泥石流“概化五邊形”流量曲線計(jì)算結(jié)果見(jiàn)圖5、6。

圖5 啟動(dòng)點(diǎn)1泥石流流量曲線

圖6 啟動(dòng)點(diǎn)2泥石流流量曲線

2.3 基底摩擦模型及參數(shù)

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

3 泥石流運(yùn)動(dòng)特征對(duì)比分析

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

圖7 堆積范圍對(duì)比

圖8 調(diào)查處泥位線高度

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

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

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

在非臺(tái)風(fēng)暴雨工況下,泥石流分別在700、200 s到達(dá)監(jiān)測(cè)點(diǎn)j1、j2,約在1 700 s時(shí)兩處泥石流匯流,匯流后的泥石流在1 800 s運(yùn)動(dòng)至j3,2 000 s時(shí)運(yùn)動(dòng)至j4,2 200 s時(shí)到達(dá)j5并在溝口停留堆積。在臺(tái)風(fēng)暴雨型降雨工況下,泥石流分別在600、100 s到達(dá)監(jiān)測(cè)點(diǎn)j1、j2,約在1 400 s時(shí)兩處泥石流匯流,匯流后的泥石流在1 500 s運(yùn)動(dòng)至j3,1 600 s時(shí)運(yùn)動(dòng)至j4,1 900 s時(shí)到達(dá)j5并在溝口停留堆積。泥石流到達(dá)各監(jiān)測(cè)點(diǎn)的時(shí)間分布見(jiàn)圖11,由此可見(jiàn)臺(tái)風(fēng)暴雨型泥石流暴發(fā)后,在更少的時(shí)間內(nèi)到達(dá)沿途監(jiān)測(cè)點(diǎn)位置,能更快對(duì)沿途溝道及建筑物造成破壞,居民的逃生時(shí)機(jī)更短,這也是臺(tái)風(fēng)暴雨型泥石流通常會(huì)造成更大生命財(cái)產(chǎn)損失的重要原因之一。

圖11 到達(dá)監(jiān)測(cè)點(diǎn)時(shí)間對(duì)比

取2種工況下流深區(qū)別較大的j1、j2、j5進(jìn)行統(tǒng)計(jì),見(jiàn)圖12,各監(jiān)測(cè)點(diǎn)在相同時(shí)刻泥石流流深以及最終流深,臺(tái)風(fēng)暴雨型泥石流均大于非臺(tái)風(fēng)暴雨型泥石流。在監(jiān)測(cè)點(diǎn)j5表現(xiàn)最為明顯,2 500 s時(shí)工況1最終流深5.86 m,工況2最終流深7.47 m。

圖12 監(jiān)測(cè)點(diǎn)流深對(duì)比

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

圖13 堆積區(qū)最大流深

圖14 堆積區(qū)面積

綜上分析,臺(tái)風(fēng)暴雨型泥石流在整個(gè)運(yùn)動(dòng)過(guò)程中具有更快的前進(jìn)速度,龍頭能在更短的時(shí)間內(nèi)對(duì)沿途溝道及建筑物造成影響和破壞,此外,其泛濫過(guò)程中的垂直影響高度和平面影響范圍也都大于非臺(tái)風(fēng)暴雨型泥石流,說(shuō)明臺(tái)風(fēng)暴雨型泥石流具備更大的泥石流物質(zhì)量,其攜帶能力及破壞能力都遠(yuǎn)強(qiáng)于非臺(tái)風(fēng)暴雨型泥石流,在沒(méi)有防治措施的情況下對(duì)工程建筑安全及居民的生命財(cái)產(chǎn)安全具有很大的威脅能力。

4 結(jié)論

a)臺(tái)風(fēng)暴雨型泥石流及非臺(tái)風(fēng)暴雨型泥石流運(yùn)動(dòng)方向及整體運(yùn)動(dòng)過(guò)程基本一致。啟動(dòng)點(diǎn)1與2處泥石流匯合,表現(xiàn)出疊加增強(qiáng)效果,運(yùn)動(dòng)至溝口后快速停留堆積,形成上庫(kù)庫(kù)區(qū)、青龍溝下游溝道、青龍溝溝口3處主要堆積區(qū)。

b)按50年一遇降雨強(qiáng)度設(shè)計(jì)的非臺(tái)風(fēng)暴雨條件下,泥石流分別在700、200、1 800、2 000、2 200 s運(yùn)動(dòng)至監(jiān)測(cè)點(diǎn)j1—j5,堆積區(qū)最大堆積厚度分別為2.65、3.77 、6.06 m,堆積面積分別為33 378.4、3 349.8、18 824.9 m2。

c)在臺(tái)風(fēng)暴雨型降雨條件下,泥石流分別在600、100、1 500、1 600、1 900 s運(yùn)動(dòng)至監(jiān)測(cè)點(diǎn)j1—j5,堆積區(qū)最大堆積厚度分別為3.08、3.80、7.75 m,堆積面積分別為36 324.4、3 700.4、32 279.5 m2。

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