基于Massflow模型的西藏仁布杰仲溝泥石流運動特征分析

段學(xué)良，馬鳳山，郭 捷，孫琪皓

(1.中國科學(xué)院地質(zhì)與地球物理研究所中國科學(xué)院頁巖氣與地質(zhì)工程重點實驗室，北京 100029；2.中國科學(xué)院地球科學(xué)研究院，北京 100029；3.中國科學(xué)院大學(xué)，北京 100049)

0 前言

泥石流是山區(qū)發(fā)育的一種自然地質(zhì)現(xiàn)象，泥石流具有突然性以及流速快，流量大，物質(zhì)容量大和破壞力強等特點[1-3]。發(fā)生泥石流常常會沖毀公路鐵路等交通設(shè)施甚至村鎮(zhèn)等，造成巨大損失。2018年7月15日，在距離研究區(qū)20 km的G318國道西藏日喀則市仁布段4 739 km+200 m處發(fā)生泥石流災(zāi)害，導(dǎo)致交通中斷，并造成經(jīng)濟(jì)損失。所以，對泥石流災(zāi)害進(jìn)行模擬分析，并預(yù)測泥石流的流動距離和泛濫范圍，從而提出防治措施，對于研究區(qū)內(nèi)的國道通行以及人員安全有著重要的意義。

一般研究泥石流的方法是基于水源、物源、流通條件區(qū)域調(diào)查數(shù)據(jù)的統(tǒng)計評價和動力分析方法，這種方式結(jié)果可信度高，但耗費大、操作繁瑣，經(jīng)過參數(shù)校正的數(shù)值模擬也可以提供合理的分析，且成本低廉，是一種重要的研究方式[4]。因此，通過數(shù)值分析方法建立泥石流的數(shù)值模型，對泥石流的啟動過程、流動過程和堆積過程進(jìn)行模擬，得到的結(jié)果借助數(shù)據(jù)、圖形等方式將泥石流的運動過程呈現(xiàn)出來，能夠為泥石流災(zāi)害防治規(guī)劃和設(shè)計作為參考[5-7]。胡明鑒等[8]利用二維顆粒流程序(PFC2D)分析了降雨作用下松散碎屑物質(zhì)啟動形成泥石流的過程。王純祥等[9]采用基于GIS的深度積分的二維數(shù)值模型，模擬了泥石流的運動規(guī)律。HUNGR等[10]使用基于半經(jīng)驗方法和等效流體的DAN3D模型模擬了泥石流的運動過程。Peng等、Chang等、楊濤等使用基于有限差分原理的FLO-2D 模型對泥石流進(jìn)行模擬，并進(jìn)行了危險性評價[11-13]。

SAVAGE等[14]提出了基于深度積分的連續(xù)介質(zhì)力學(xué)模型，并利用該模型成功模擬了干性沙土的流動過程。該模型的優(yōu)勢是簡化了Navier-Stokes模型，大大降低了方程的求解難度，有效地提高了計算效率[15]。OUYANG等提出了基于深度積分的Massflow 數(shù)值分析方法，并利用Massflow 軟件對滑坡、泥石流等地質(zhì)災(zāi)害進(jìn)行了數(shù)值模擬，不僅很好吻合理論解和實驗結(jié)果，而且還原了真實災(zāi)害的發(fā)生過程[16-18]。

以往研究，雖然在整個流域范圍內(nèi)對泥石流進(jìn)行了分析或評價，但選取的DEM精度都較低，因此不能精細(xì)的刻畫泥石流的運動規(guī)律。為探明研究區(qū)流域上游的潛在物源在極端條件下發(fā)生泥石流災(zāi)害，對溝口建筑及公路的影響，本研究利用無人機技術(shù)獲取研究區(qū)的高精度DEM數(shù)據(jù)。并采用基于泥石流運動深度積分的Massflow數(shù)值分析方法模擬了杰仲溝泥石流的運動過程。

1 研究區(qū)地質(zhì)特征

研究區(qū)杰仲溝位于西藏自治區(qū)日喀則市仁布縣西北方向16 km處(圖1)，G318國道從溝口穿過，地理位置為東經(jīng)89°41′10″～89°43′30″，北緯29°16′00″～29°19′27″，屬于深切構(gòu)造高中山地形，流域總體呈樹葉狀，呈南北向展布，表現(xiàn)為南高北低形態(tài)。杰仲溝流域長約6 km,平均寬約0.3 km，匯水面積約27 km2。杰仲溝流域最高點海拔5 394 m，溝口海拔3 794 m，相對高差1 600 m，流域平均縱坡306‰，流域地貌深切，呈“U”型寬谷，溝谷兩側(cè)斜坡坡度多在30°左右，物源區(qū)坡度約30°。雅魯藏布江由西向東從堆積區(qū)上切過。

圖1 研究區(qū)及斷裂位置分布Fig.1 The location of the study area and the distribution of the fault

泥石流流域內(nèi)主要分布紫紅、灰色礫巖、砂巖、頁巖夾凝灰?guī)r。區(qū)域內(nèi)活動構(gòu)造發(fā)育，主要受第四系活動斷裂達(dá)機翁-彭錯林-朗縣斷裂影響較大。仁布縣屬于地震強烈活動區(qū)，杰仲溝的地震峰值加速度為0.2g，對應(yīng)的地震基本烈度為Ⅷ度。仁布縣屬于溫帶半干旱高原季風(fēng)氣候，年平均降水量450 mm，雨季集中在7、8月份，降水量占全年的95%。在雨季尤其是暴雨時期，研究區(qū)內(nèi)會增加大量的崩塌、滑坡現(xiàn)象，為泥石流提供了松散固體物源。

2 泥石流形成條件

2.1 地形地貌條件

杰仲溝泥石流流域?qū)儆谏钋袠?gòu)造高中山地貌，相對高差1 600 m。上游形成區(qū)地形三面環(huán)山，出口為漏斗狀，地形比較開闊，周圍山高坡陡、植被生長不良，這樣的地形有利于降雨的匯集。流通區(qū)狹長，溝谷縱坡降大，為松散固體物質(zhì)參與泥石流活動提供了有利條件。下游堆積區(qū)的地形為開闊平坦的河谷。較大的高差和匯水面積使得雨水快速而大量匯集，并侵蝕地表和帶動松散固體物質(zhì)的運動，從而極易引發(fā)泥石流的形成。

2.2 物源條件

研究區(qū)地質(zhì)構(gòu)造復(fù)雜，區(qū)內(nèi)地震烈度較高，且發(fā)育有新構(gòu)造斷裂，如圖1所示的達(dá)機翁-彭錯林-朗縣斷裂從研究區(qū)穿過(在圖1中簡稱為D-P-L斷裂)。流通區(qū)溝道兩側(cè)的坡體上未發(fā)現(xiàn)有大量的破碎松散體，流域內(nèi)主要松散體物源位于上游的構(gòu)造斷裂附近(圖2)，在強烈的構(gòu)造作用下, 巖石裂隙發(fā)育, 巖體破碎，同時研究區(qū)地處高寒山區(qū)，冬季長，氣溫差異較大，尤其是在春融季節(jié)早晚溫差更大，加之植被稀少，使基巖的物理以及化學(xué)風(fēng)化作用強烈。長期積累的松散固體物質(zhì)儲量大, 分布集中, 有利于泥石流的補給。

圖2 杰仲溝流域全貌Fig.2 The panoramic of Jiezhonggou watershed

2.3 水源條件

研究區(qū)年平均降雨量450 mm，7、8月份降雨最為集中，降水量占年降水量的95%, 大雨、暴雨出現(xiàn)頻率較高, 加之氣溫升高, 上游的冰雪消融, 泥石流活動頻率, 規(guī)模隨之增大, 9月降水逐漸減少, 泥石流活動進(jìn)入低潮期?？偟膩碚f，冬春季節(jié)降水量少，夏季降水多，豐富的降雨為泥石流的形成提供了水源和動力條件。

3 Massflow 模型

3.1 模型控制方程

對流體動力學(xué)三維Navier-Stokes方程進(jìn)行深度積分，推導(dǎo)出質(zhì)量和動量控制方程。深度積分是指將控制方程中的各種變量沿深度方向進(jìn)行積分，這里將對質(zhì)量守恒方程和動量守恒方程進(jìn)行深度積分，再運用萊布尼茲法則和動力學(xué)邊界條件對方程化簡，最終得到控制方程如下，公式(1)為質(zhì)量守恒方程，公式(2)、(3)為動量守恒方程。

(1)

(2)

(3)

式中：ρ——流體密度/(kg·m-3)；

t——時間；

u，v——x方向和y方向上的流體速度(m/s)；

gx，gy，gz——各坐標(biāo)軸上的重力分量；

Zb——河床標(biāo)高(m)；

h——泥石流體的泥深(m)；

τb——底部的剪應(yīng)力(Pa)；

kap——土壓力系數(shù)。

Massflow 在控制方程的求解方面采用MacCormack-TVD 有限差分法。有限差分法在離散點上直接用差分代替微分，相比于有限元法以及有限體積法的優(yōu)點是更加簡單成熟及可以構(gòu)造高精度格式，其數(shù)值方法具體介紹見文獻(xiàn)[16]。

3.2 模型參數(shù)選取

3.2.1地形條件

以實地?zé)o人機拍攝的照片數(shù)據(jù)為基礎(chǔ)(拍攝區(qū)域如圖2中的紅色陰影所示)，建立數(shù)字高程模型(DEM)，進(jìn)而生成研究區(qū)的地形信息。首先，利用攝影三維建模軟件PhotoScan，把研究區(qū)的無人機航拍照片進(jìn)行分析處理，并生成研究區(qū)的DEM數(shù)據(jù)(分辨率為1 m)。然后，運用GlobalMapper、Surfer軟件，對DEM數(shù)據(jù)文件進(jìn)行裁剪、生成矢量格式等操作，得到流域和物源坐標(biāo)及高程信息的柵格數(shù)據(jù)。最后，利用ArcGIS將得到的流域和物源的柵格數(shù)據(jù)投影，并轉(zhuǎn)換為Massflow軟件能識別的ASCⅡ文件格式。

3.2.2物源條件

在流域的上游存在松散破碎區(qū)，如圖2中黃色虛線所圈的區(qū)域。由前文分析可知，該破碎區(qū)是一個潛在的泥石流物源，在極端條件下如暴雨、地震等，發(fā)生泥石流的概率較高。通過在地圖上測量，得到潛在物源區(qū)的面積約為6×105m2。

3.2.3泥石流的運動條件

根據(jù)現(xiàn)場實地調(diào)查，并參考相關(guān)研究[19]，本模擬選取泥石流容重為1 800 kg/m3。考慮到本模擬中，流體的初始運動位置不是實際的起點位置。因此，通過在低精度DEM(30 m分辨率)的條件下進(jìn)行模擬，估算出流體的初始速度為10 m/s，所以在進(jìn)行高精度DEM模擬時，給流體賦一個沿主溝方向的大小為10 m/s的初速度。

分別選用Manning模型和Voellmy摩擦模型，進(jìn)行初步模擬計算，以選取適用于該研究區(qū)的摩擦模型和相應(yīng)的運動參數(shù)。

(1)Voellmy模型

對于Voellmy模型在庫倫摩擦模型的基礎(chǔ)上增加了湍流項，認(rèn)為運動的流體受到的抗剪應(yīng)力為流體所受摩擦力和湍流流動產(chǎn)生的額外阻力之和，因此該模型由摩擦項和湍流項兩部分組成，需確定摩擦系數(shù)和湍流系數(shù)，其表達(dá)式如下：

(4)

式中：τb——底部的剪應(yīng)力/Pa;

σ——正應(yīng)力/Pa;

μ——摩擦系數(shù)；

ρ——泥石流容重/(kg·m-3)；

ξ湍流系數(shù)/(m·s-2)。

根據(jù)相關(guān)研究[20]并結(jié)合現(xiàn)場調(diào)查情況，將摩擦系數(shù)μ設(shè)定為0.1～0.3，湍流系數(shù)ξ設(shè)定為150～250 m/ s2。

(2)Manning模型

對于Manning模型，需要確定的計算參數(shù)是基底糙率即曼寧系數(shù)。本研究根據(jù)王裕宜等[21]提出的不同阻力介質(zhì)狀態(tài)下統(tǒng)一的阻力糙率系數(shù)表征公式計算，公式如下：

nc=0.033Rns-0.51exp(0.34Rns0.17)lnh

(5)

式中：nc——曼寧系數(shù)；

Rns——泥砂比；

h——泥深/m。

根據(jù)現(xiàn)場調(diào)查情況，取Rns=0.75，h為2～10 m，代入公式(4)得到杰仲溝的曼寧系數(shù)取值范圍0.04～0.12。

為了選取適于杰仲溝泥石流數(shù)值模擬的摩擦模型和相應(yīng)運動參數(shù)，對兩種摩擦模型Voellmy和Manning模型分別設(shè)定不同的運動參數(shù)，得到了不同的泥石流泥深、流動強度以及堆積范圍，對比已發(fā)生的泥石流，選擇與實際情況最匹配的模型及運動參數(shù)。對于Voellmy模型將摩擦系數(shù)μ設(shè)定在0.1～0.3，湍流系數(shù)ξ設(shè)定為150～250 m/ s2。對于Manning模型，將曼寧系數(shù)設(shè)定在0.04～0.12。泥石流容重取1 800 kg/m3，流體體積取6×104m3。對比不同運動參數(shù)的模擬結(jié)果，最終選取μ=0.1，ξ=200 m/s2，nc=0.08。如圖3和圖4分別是Voellmy模型和Manning模型的模擬結(jié)果。

圖3 泥石流最大泥深和流動強度分布, Voellmy模型Fig.3 Maximum depth and flow intensity distribution of the debris flow, Voellmy model

圖3(a)和圖4(a)分別是兩種模型整個運動過程中最大泥深的分布，圖3(b)和圖4(b)是最大泥深與泥石流運動速度的乘積，代表了泥石流的流動強度。從泥深分布圖中可以看出，兩種模型的泥石流的堆積范圍模擬結(jié)果與實際已發(fā)生的泥石流幾乎相同，但是Manning模型的模擬結(jié)果更接近實際情況，而Voellmy模型模擬的泥石流運動范圍誤差相對較大，尤其在流通區(qū)內(nèi)。同樣在流動強度分布圖中也可看出，Manning模型比Voellmy模型的模擬結(jié)果更精確。雖然兩種模型在精度上有一定差距，但是總的規(guī)律是相同的。圖3和圖4的(a-1)和(a-2)是實際被泥石流沖毀的兩處公路石墩護(hù)欄。圖3、圖4的(a-1)中的兩條黑色虛線之間的公路西側(cè)的石墩護(hù)欄被沖毀，(a-2)中被黑色虛線圈起的公路西側(cè)的石墩護(hù)欄也被沖毀。兩種模型的模擬結(jié)果都與實際情況相匹配，尤其是Manning模型，從這兩個圖中可以看出，泥深在這兩處都達(dá)到了4 m左右，這說明模擬結(jié)果與實際匹配的較好。此外，在圖4(b)泥石流流動強度分布圖中可以看出，泥石流初始運動階段流動強度迅速增大，隨后逐漸減小，但是在上述兩處石墩護(hù)欄沖毀處，流動強度表現(xiàn)為增大。因此，也與實際情況匹配良好，這說明所選取的運動參數(shù)能夠較真實的描述研究區(qū)內(nèi)泥石流的運動狀態(tài)。綜合比較上述兩種摩擦模型，最終選取Manning模型對該泥石流溝做進(jìn)一步的模擬研究。

圖4 泥石流最大泥深和流動強度分布，Manning模型Fig.4 Maximum depth and flow intensity distribution of the debris flow, Manning model

4 模擬結(jié)果

根據(jù)估算的流域上游存在1.0×105～6.0×105m3的固體物源，考慮到運動過程中流體的損失量，本次模擬選取流體體積為3.0×105m3。利用已選模型和相應(yīng)參數(shù)對此體積的泥石流運動進(jìn)行模擬，并分析其影響范圍和危險性。

此次模擬共進(jìn)行了600 s。如圖5和圖6分別是不同時刻的泥深分布圖和整個泥石流運動過程中最大泥深及流動強度的分布圖?？梢钥闯瞿嗍鲃傞_始運動30 s最大泥深就達(dá)到了6 m，90 s時達(dá)到了8 m，且流動速度快，已達(dá)到甘丹桑阿曲林寺前，此后流動速度變慢，并出現(xiàn)堆積。這是因為甘丹桑阿曲林寺南側(cè)的溝道較為平直，而在寺前以及其北側(cè)，溝道出現(xiàn)了兩次較大的彎曲，從而使流體速度減慢并發(fā)生堆積。模擬計算進(jìn)行到180 s時，泥石流體穿過了溝口處的G318國道，并開始擴散流入堆積區(qū)。進(jìn)行到270 s時，泥石流體開始流入雅江，并在寺北側(cè)第二個彎道處發(fā)生堆積回流，流入西側(cè)的支溝中。450 s時，大部分流體流進(jìn)雅江，最大泥深減小到4 m，600 s時泥石流運動基本結(jié)束。此次模擬，整個泥石流運動過程中出現(xiàn)的最大泥深高達(dá)9.6 m。

圖5 不同時刻的泥深分布Fig.5 Mud depth distribution of debris flow at different times

圖6 泥石流最大泥深和流動強度分布Fig.6 Maximum depth and flow intensity distribution of the debris flow

綜合最大泥深和流動強度分布圖，定義了研究區(qū)內(nèi)四處危險區(qū)，分別在圖6中標(biāo)為A、B、C、D，圖7是其放大圖。最大泥深出現(xiàn)在圖6(a)中的A和B處，A處在溝道第一個彎道前，流體進(jìn)入彎道后流速減小，發(fā)生堆積；B處是石墩護(hù)欄被沖毀的地點，此處位于兩個彎道之間，所以泥石流也易在此堆積；另一個護(hù)欄被沖毀的地點C處位于第二個彎道處，最大泥深達(dá)到了8.7 m。G318國道南側(cè)溝道D處的最大泥深為7.8 m，國道上的最大泥深達(dá)到了5 m。此次模擬中，流動強度的最大值為104 m2/s，流動強度大的區(qū)域集中在第一個彎道南側(cè)，即A處以南，這部分區(qū)域中溝道較為平直，泥石流流速快，沖擊力強，但A處之后流動強度有所減小，但在B處，流動強度再次出現(xiàn)一個較大值，達(dá)到了75 m2/s，在B處危險區(qū)邊緣有一通向上游村莊的公路，從泥石流溝上架橋而過。G318國道處，達(dá)到36 m2/s，泥石流穿過國道后流入堆積區(qū)，流動強度逐漸減小。

圖7 危險區(qū)放大圖Fig.7 Enlarged map of risk zones

5 結(jié)論

對實際已發(fā)生的泥石流進(jìn)行了數(shù)值模擬，模擬結(jié)果與實際情況匹配良好，從而確定了適合研究區(qū)杰仲溝泥石流模擬的摩擦模型和運動參數(shù)。

在選取了摩擦模型和相應(yīng)運動參數(shù)的基礎(chǔ)上，又對可能發(fā)生的更大體積的泥石流災(zāi)害進(jìn)行了模擬，對泥石流從發(fā)生到休止，整個動態(tài)過程進(jìn)行了分析。主要描述了不同時刻的泥深分布情況，以及對整個運動過程中的最大泥深分布和流動強度分布進(jìn)行了分析，并劃定了研究區(qū)內(nèi)的四處危險區(qū)，從而得出以下結(jié)論。

(1)泥深和泥石流流動強度較大的區(qū)域，主要分布危險區(qū)A及以南的區(qū)域，即溝道第一個彎道南部，但處于這一范圍的大部分區(qū)域人類活動較少，僅在彎道處有一條通向寺廟的公路。

(2)危險區(qū)B在泥石流溝的第一個彎道與第二個彎道之間，危險區(qū)C在第二個彎道處，并且這兩處的模擬結(jié)果與實際情況較為符合，公路旁的石墩護(hù)欄被泥石流沖毀，并且危險區(qū)B內(nèi)有一通向上游村莊的公路，從泥石流溝上架橋而過。因此，在地震(地震烈度為Ⅷ度)和極端氣候條件如暴雨的情況下，這兩個區(qū)域的危險性相對較高，需對其加強防護(hù)。

(3)G318國道在溝口穿過，危險區(qū)D位于國道南側(cè)附近，雖然流動強度相對其他危險區(qū)不大，但是流體在這個區(qū)域容易發(fā)生堆積，最大泥深達(dá)到了7.8 m。因此，大量堆積的泥砂、石塊等容易造成對國道的沖埋，嚴(yán)重影響交通通行。所以，應(yīng)提高此處公路下涵洞的排導(dǎo)能力，預(yù)防災(zāi)害的發(fā)生。