侯圣山,曹 鵬,陳 亮,馮 振,王立朝,李 昂,劉軍友,李陽(yáng)光,鄭 浩
(1.中國(guó)地質(zhì)環(huán)境監(jiān)測(cè)院(自然資源部地質(zhì)災(zāi)害技術(shù)指導(dǎo)中心),北京 100081;2.中國(guó)地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院,北京 100083)
泥石流是指由強(qiáng)降雨誘發(fā),在地形陡峭、巖性松散的溝谷或者斜坡上形成的特殊洪流,表現(xiàn)為固體、液體、氣體三相混合流體,一般具有歷時(shí)短、強(qiáng)度高、破壞性大等特點(diǎn)[1]。山區(qū)泥石流災(zāi)害防控越來(lái)越受到重視,而泥石流危險(xiǎn)性評(píng)價(jià)是泥石流災(zāi)害防控的關(guān)鍵環(huán)節(jié)[2]。泥石流危險(xiǎn)性評(píng)價(jià)方法諸多,傳統(tǒng)泥石流危險(xiǎn)性評(píng)價(jià)方法主要有層次分析法、模糊綜合評(píng)價(jià)、邏輯回歸模型、物源可拓模型等。層次分析法確定指標(biāo)權(quán)重具有很大的主觀性;其它評(píng)價(jià)模型在泥石流危險(xiǎn)性區(qū)劃中以流域?yàn)榛締卧?,無(wú)法確定流域中精確地點(diǎn)的危險(xiǎn)程度,也無(wú)法對(duì)泥石流的運(yùn)動(dòng)進(jìn)行精確刻畫;而數(shù)值方法評(píng)價(jià)的重點(diǎn)在于對(duì)泥石流流動(dòng)過(guò)程的刻畫,可識(shí)別流動(dòng)速度、堆積厚度的空間分布及其隨時(shí)間的變化,非常適合小流域中人口聚集區(qū)泥石流的危險(xiǎn)性評(píng)價(jià)。宋兵等[3]運(yùn)用RAMMS 模型模擬了白沙溝泥石流運(yùn)動(dòng)特征;Hungr[4]使用DAN3D 模型對(duì)泥石流的流動(dòng)過(guò)程進(jìn)行了模擬;Magirl 等[5]使用Laharz模型對(duì)圣卡塔利納山脈7 條泥石流的運(yùn)動(dòng)距離和堆積范圍進(jìn)行了模擬,并與實(shí)際情況相驗(yàn)證;Horton 等[6]將二維動(dòng)態(tài)泥石流模型Massflow 運(yùn)用于汶川地震震中地區(qū)流域尺度的災(zāi)害圖編制;Lin 等[7]采用二維洪水與土石流數(shù)值模擬軟件FLO-2D 對(duì)臺(tái)風(fēng)誘發(fā)的松河地區(qū)泥石流進(jìn)行了分析;常鳴等[8]運(yùn)用FLO-2D 對(duì)都江堰八一溝泥石流進(jìn)行了多種降雨頻率條件下的危險(xiǎn)性評(píng)價(jià)。上述數(shù)值分析方法在泥石流過(guò)程模擬方面均取得了較好的效果,尤其是FLO-2D 模型在泥石流危險(xiǎn)性模擬方面潛力巨大。
泥石流災(zāi)害的致災(zāi)因子主要為強(qiáng)降雨,在相對(duì)脆弱的自然環(huán)境條件下,短時(shí)強(qiáng)降雨極易引發(fā)群發(fā)性泥石流[9]。2012年5月10日,甘肅岷縣近一個(gè)小時(shí)大范圍冰雹及強(qiáng)降雨誘發(fā)耳陽(yáng)河流域突發(fā)泥石流,沖擊溝口及河道兩側(cè)村落、公路、輸電設(shè)備,造成18 人死亡,直接損失4 800 余萬(wàn)元[10]。經(jīng)過(guò)與岷縣多個(gè)氣象站點(diǎn)歷史降水記錄相對(duì)比,確定誘發(fā)本次泥石流的降水強(qiáng)度約為百年一遇。
耳陽(yáng)河流域泥石流溝眾多,主溝平均縱比降為235‰,主溝和支溝塌岸嚴(yán)重。區(qū)內(nèi)大量種植黃芪、黨參、玉米、胡豆等作物。每年5—6月為農(nóng)耕時(shí)節(jié),農(nóng)業(yè)生產(chǎn)后土壤松散,抗侵蝕能力差,一旦出現(xiàn)強(qiáng)降雨,容易發(fā)生群發(fā)性泥石流災(zāi)害。泥石流的威脅巨大,識(shí)別危險(xiǎn)地段、預(yù)測(cè)致災(zāi)過(guò)程至關(guān)重要。因此,本文運(yùn)用FLO-2D 模型,選取耳陽(yáng)河流域6 條具有代表性的泥石流,模擬了2012年5月10日實(shí)際降雨(1%降雨概率)情況下的泥石流過(guò)程,并與地質(zhì)災(zāi)害現(xiàn)場(chǎng)調(diào)查成果對(duì)比,以驗(yàn)證模擬的可靠性。進(jìn)而,運(yùn)用相同的評(píng)價(jià)模型,模擬上述泥石流在2.0%和0.2%降雨頻率條件下的堆積范圍、深度、流速,并進(jìn)行了泥石流危險(xiǎn)性分析,為泥石流災(zāi)害防治提供依據(jù)。
岷縣位于甘肅南部,屬于西秦嶺北支中段,緊鄰長(zhǎng)江流域和黃河流域的分界線,洮河穿流而過(guò),為典型的中高山侵蝕地貌,是甘肅省泥石流活躍區(qū)之一[11]。
研究區(qū)位于岷縣洮河右岸支流耳陽(yáng)河流域(圖1),地處西秦嶺山地和黃土高原的交匯區(qū),新構(gòu)造運(yùn)動(dòng)活躍,地貌上呈現(xiàn)為年輕的構(gòu)造侵蝕中山地貌類型。耳陽(yáng)河流域面積約63 km2,主溝長(zhǎng)度約15.8 km,流域內(nèi)最高海拔為東山頂?shù)? 130 m,最低海拔2 310 m,位于耳陽(yáng)河匯入納納河的溝口。耳陽(yáng)河的主溝和支溝多呈“V”型谷,溝壁陡峭,臨空面發(fā)育,滑坡塌岸多見;溝內(nèi)耕地眾多,土地松散,物源豐富。區(qū)內(nèi)植被覆蓋總體上差,耳陽(yáng)河兩岸巖土裸露,僅有農(nóng)作物和零星灌木,上游局部可見人工林。
圖1 耳陽(yáng)河流域地形圖Fig.1 Topographic map of the Eryang river watershed
流域內(nèi)出露的地層主要是泥盆系、侏羅系的板巖、頁(yè)巖和砂巖等,經(jīng)過(guò)構(gòu)造作用和風(fēng)化作用,巖體破碎,穩(wěn)定性較差;區(qū)內(nèi)第四紀(jì)沉積物主要為黃土、沖洪積物、泥石流堆積物、殘坡積物,耳陽(yáng)河下游東岸(右岸)剖面圖見圖2。
圖2 耳陽(yáng)河右岸地質(zhì)剖面圖Fig.2 Geological section of the right bank of the Eryang River
研究區(qū)為高寒濕潤(rùn)氣候帶,多年平均降水量為560.8 mm,6—8月降水占總降水量的 60%以上。受地形及海拔高度影響,降水量空間分布差異較大,呈現(xiàn)出由河谷到山頂降水量逐漸增大的趨勢(shì)。
耳陽(yáng)河屬于黃河水系,為黃河的三級(jí)支流:耳陽(yáng)河在茶埠鎮(zhèn)溝門村匯入納納河,之后向西流約500 m,在茶埠鎮(zhèn)陽(yáng)坡匯入洮河。經(jīng)過(guò)野外調(diào)查,耳陽(yáng)河流域內(nèi)共發(fā)育23 條泥石流溝,其中左岸9 條,右岸14 條,泥石流溝線密度為1.45 條/km。
本研究選取的6 條泥石流溝對(duì)稱分布于耳陽(yáng)河兩岸(表1),均處于活躍期,呈“V 形谷”,主溝縱比降為189.7‰~316.9‰,平均坡度35°~50°。松散物質(zhì)類型主要為坡面的殘坡積物以及耕作土,其次為溝道兩岸的小型崩塌、滑坡及溝岸坍塌,其余為溝底再搬運(yùn)形成的松散物質(zhì)。6 條溝谷中都有大面積耕地,占溝谷總面積的51%,義子溝、接哈溝、扎龍溝上游農(nóng)田已經(jīng)退耕,牌嘴堡北溝(牌嘴溝)和接哈溝上游有人工林。6 條溝溝底情況相似,可見泥石流相砂礫,碎石以及泥球,碎石粒徑多在2~10 cm,最大可見1.3 m,母巖主要為砂巖、砂質(zhì)板巖、灰?guī)r、灰質(zhì)板巖。6 條泥石流溝的補(bǔ)給段長(zhǎng)度占各自溝道總長(zhǎng)度的61.8%~69.2%。接哈溝堆積扇擠壓耳陽(yáng)河河道,使耳陽(yáng)河出現(xiàn)一處牛軛狀河曲(圖3)。
表1 泥石流溝特征參數(shù)Table 1 Property list of the key debris flow gullies
圖3 本研究選取的六條泥石流平面分布圖Fig.3 Distribution of the six gullies in this study
2012年“5·10”特大暴洪泥石流發(fā)生之后,在牌嘴溝中修建了10 道攔擋壩。2019年9月現(xiàn)場(chǎng)調(diào)查發(fā)現(xiàn),前6 道攔擋壩已被泥石流堆積物、兩側(cè)崩滑堆積物淤滿、掩埋,后4 道攔擋壩仍具有攔擋作用。在扎龍溝溝口修建了約400 m 長(zhǎng)的排導(dǎo)渠,寬10~15 m,深3 m。
FLO-2D 模型由O’Brien[12-13]提出,采用有限差分計(jì)算垂向深度和流動(dòng)速度,預(yù)測(cè)泥石流流動(dòng)和堆積范圍,可用于泥石流、洪水、潰壩、城市淹沒(méi)等過(guò)程的模擬,并可用于災(zāi)害危險(xiǎn)性評(píng)價(jià)。
FLO-2D 軟件在模擬泥石流運(yùn)動(dòng)過(guò)程中,將地形分為若干等大網(wǎng)格,每一個(gè)網(wǎng)格中,其曼寧系數(shù)和高程值都是唯一的,通過(guò)運(yùn)動(dòng)方程以及連續(xù)方程,可以計(jì)算出每個(gè)網(wǎng)格中流體深度、流量,進(jìn)而得知流體的運(yùn)動(dòng)范圍,通過(guò)動(dòng)量方程計(jì)算出相鄰網(wǎng)格間流體的速度變化。FLO-2D 模型方程如下:
式中:t—泥石流運(yùn)動(dòng)時(shí)間/s;
x—X軸方向距離/m;
y—Y軸方向距離/m;
h—流體深度/m;
I—有效降雨強(qiáng)度/(mm·h-1);
u—X軸方向流速/(m·s-1);
v—Y方向流速/(m·s-1)。
連續(xù)方程控制了泥石流運(yùn)動(dòng)時(shí)每個(gè)網(wǎng)格內(nèi)質(zhì)量守恒。
式中:g—重力加速度/(m·s-2);
Sfy、Sfx—X、Y方向摩擦坡降;
Sox、Soy—X、Y方向床底坡降。
X、Y方向的運(yùn)動(dòng)方程保持了運(yùn)動(dòng)平衡。
式中:Sf—摩擦坡降;
Sy—屈服坡降;
Sv—黏性坡降;
Std—分散坡降;
η—黏滯系數(shù);
k—層流阻滯系數(shù);
γm—泥石流容重/(kN·m-3);
ntd—等效曼寧系數(shù)。
流變方程考慮了泥石流運(yùn)動(dòng)時(shí)顆粒之間的碰撞對(duì)泥石流流動(dòng)阻力的影響。
2.1.1 地形數(shù)據(jù)處理
地形數(shù)據(jù)來(lái)自1∶5 萬(wàn)數(shù)字高程模型(DEM)。將高程數(shù)據(jù)轉(zhuǎn)化為ASCII 格式,導(dǎo)入FLO-2D,將之剖分為10 m×10 m 的評(píng)價(jià)單元,通過(guò)插值計(jì)算,確定每個(gè)評(píng)價(jià)單元的高程,完成地形數(shù)據(jù)處理。
用ArcGIS 將研究區(qū)流域邊界、出水點(diǎn)、攔擋壩等主要地形要素矢量化,導(dǎo)入FLO-2D。
2.1.2 確定泥石流參數(shù)
郭富赟等[14]對(duì)耳陽(yáng)河流域泥石流流體重度進(jìn)行了測(cè)算,為17.0~17.6 kN/m3。本文采用文獻(xiàn)[14]的泥石流流體重度測(cè)算結(jié)果,并結(jié)合堆積物平均重度,計(jì)算出泥石流泥沙體積濃度為46.67%~50.67%,取中間值48.67%。
高濃度泥石流在屈服應(yīng)力的作用下出現(xiàn)層流特征,層流之間的摩擦作用即可用層流阻滯系數(shù)(k)表示,本文參考前人研究實(shí)例,用工程地質(zhì)類比法,k取值定為2 250。
黏滯系數(shù)(η)和屈服應(yīng)力(τy)的大小主要取決于泥沙體積濃度,關(guān)系式如下:
式中:α1、α2、β1、β2—試驗(yàn)系數(shù);
Cv—體積濃度。
根據(jù)泥石流堆積物、流體的容重和物質(zhì)組成,判斷研究區(qū)泥石流為黏性,且為中阻性,泥沙比(Rns)為0.75。王裕宜等[15]統(tǒng)一的泥沙比-體積濃度-流變參數(shù)關(guān)系式如下:
式中:nc—曼寧系數(shù);
h—泥深/m(本文以平均堆積厚度表示)。
將各泥石流溝泥深帶入式(7)即可得到拉龍溝、扎龍溝、拉路溝、接哈溝、義子溝、牌嘴溝的曼寧系數(shù),分別為0.019 3,0.022 4,0.004 5,0,0,0.022 4。將式(5)(6)(8)(9)聯(lián)立,求解得到α1、α2、β1、β2分別為0.000 16,0.002 30,17.41,18.28。
根據(jù)野外調(diào)查和訪問(wèn)得知,耳陽(yáng)河流域泥石流于2012年5月10日17 時(shí)50 分許開始爆發(fā),18 時(shí)左右達(dá)到洪峰,21 時(shí)逐步回落,泥石流過(guò)程歷時(shí)約3 h。降雨量數(shù)據(jù)選用岷縣麻子川自動(dòng)站降雨觀測(cè)數(shù)據(jù),如圖4所示[16]。流量過(guò)程曲線計(jì)算選用《四川省中小流域暴雨洪水計(jì)算手冊(cè)》中川西地區(qū)水文模型,本次研究區(qū)位于甘肅省南部,和川西地區(qū)距離較近,地質(zhì)、地貌、氣候均有相似性,此計(jì)算手冊(cè)中的參數(shù)在本區(qū)降雨-徑流分析中具有良好的適用性。運(yùn)用推理公式法概化各泥石流溝清水流流量過(guò)程曲線,并乘以放大系數(shù)BF(BF=1/(1-Cv)=2.10),得到泥石流流量過(guò)程曲線,見圖5。
圖4 岷縣麻子川自動(dòng)站降水曲線圖(2012年5月10日)Fig.4 Precipitation curve of the Mazichuan monitoring station(May 10,2012)
從圖5可以看出,各泥石流溝于爆發(fā)15~30 min后達(dá)到洪峰,約3 h 后流量逐步回落;這與現(xiàn)場(chǎng)調(diào)查和訪問(wèn)了解到的“5·10”當(dāng)時(shí)的情況基本相符。
圖5 清水流量過(guò)程曲線(a)和泥石流流量過(guò)程曲線(b)Fig.5 Flow process curve of(a)clear water and(b)debris flow
模擬結(jié)果顯示:泥石流流體在流通區(qū)運(yùn)動(dòng)速度較快,最大速度可達(dá)7 m/s;在溝口附近速度迅速下降,大部分泥石流匯入耳陽(yáng)河,部分區(qū)域堆積深度較大,堵塞河道(圖6)。
圖6 “5·10”泥石流深度(a)和流速(b)Fig.6 “5·10” debris flow(a)depth and(b)velocity
通過(guò)FLO-2D 模擬,得到了6 條泥石流溝的堆積扇面積和泥石流堆積平均厚度,如表2所示。實(shí)地調(diào)查測(cè)量得到的面積和厚度同時(shí)列于表2,并將模擬結(jié)果和實(shí)際面積進(jìn)行對(duì)比。可以看出:除拉龍溝和拉路溝外,其他各溝堆積扇面積誤差率在-6.2%~28.5%之間、堆積扇均厚誤差率在-7.5%~23.3%之間。模擬堆積厚度大部分低于實(shí)際厚度,主要原因是自“5·10”泥石流爆發(fā)后,每年有少量松散物質(zhì)沖出溝口,堆積扇厚度逐年加大,實(shí)測(cè)的厚度比“5·10”泥石流實(shí)際堆積厚度大。拉龍溝與拉路溝堆積面積誤差率分別達(dá)到189.8%和2 827.5%,模擬堆積扇面積遠(yuǎn)遠(yuǎn)大于實(shí)際面積,主要原因?yàn)榇迕裨跍峡诟浇ǚ?、開墾土地,對(duì)泥石流堆積扇進(jìn)行了人為改造。
表2 模擬結(jié)果與調(diào)查結(jié)果對(duì)比表Table 2 Comparison of the simulation and survey results
當(dāng)前泥石流危險(xiǎn)性分級(jí)指標(biāo)主要有流速和泥深[17]、沖擊力和堆積厚度[18]、頻率和強(qiáng)度[19-20]等,耳陽(yáng)河流域泥石流的危害性主要在于泥石流攜帶的固體物質(zhì)淤埋農(nóng)田和房屋,以及中高速黏性洪流對(duì)沿岸的沖擊和侵蝕。因此,本次選取流速、泥深等指標(biāo)進(jìn)行危險(xiǎn)性分級(jí),采用如表3所示的邏輯關(guān)系,確定泥石流危險(xiǎn)性的空間分布。
表3 泥石流危險(xiǎn)性分區(qū)指標(biāo)Table 3 Risk classification of debris flow
危險(xiǎn)性評(píng)價(jià)結(jié)果見表4、圖7?!?·10”泥石流總危險(xiǎn)區(qū)面積為57 900 m2,占總流通區(qū)面積的60.2%,其中高危險(xiǎn)區(qū)面積占總危險(xiǎn)區(qū)面積比例較高的是接哈溝和扎龍溝,分別為34.23%和38.60%。
表4 泥石流危險(xiǎn)性分區(qū)統(tǒng)計(jì)表Table 4 Statistics of debris flow hazard zoning
由于“5·10”泥石流過(guò)后重新規(guī)劃了房屋建筑,目前耳陽(yáng)河沿岸的大部分建筑在1%降雨頻率條件下較為安全,但扎龍溝、接哈溝溝口和義子溝沿岸的部分房屋建筑易受到泥石流威脅,共涉及約21 處房屋(見圖7中的局部放大)。
圖7 “5·10”泥石流危險(xiǎn)性分區(qū)圖Fig.7 Hazard zoning map of “5·10”debris flow
通過(guò)對(duì)“5·10”過(guò)程進(jìn)行模擬,并與實(shí)際情況進(jìn)行對(duì)比,驗(yàn)證了模型和參數(shù)的可靠性,說(shuō)明FLO-2D 技術(shù)能夠較好地重現(xiàn)泥石流過(guò)程、模擬出泥石流流速、堆積深度的時(shí)空分布。使用同樣的參數(shù)及方法,對(duì)P=2.0%及P=0.2%降雨頻率條件下泥石流情況進(jìn)行了模擬,并開展了上述兩個(gè)工況的泥石流危險(xiǎn)性評(píng)價(jià)。P=2.0%及P=0.2%降雨參數(shù)取值見表5,模擬得到的流量過(guò)程曲線見圖8。
表5 不同頻率的降雨參數(shù)Table 5 Rainfall amounts of different frequencies
圖8 不同降雨頻率條件下泥石流流量過(guò)程曲線Fig.8 Flow process curve of debris flow in different rainfall frequencies
利用表3中的邏輯運(yùn)算方法,得到了P=2.0%及P=0.2%降雨概率的泥石流危險(xiǎn)性分區(qū)圖(圖9、表6)。
圖9 2.0%(a)和0.2%(b)降雨頻率條件下泥石流危險(xiǎn)性分區(qū)圖Fig.9 Hazard zoning map of debris flow in(a)2.0% and(b)0.2% rainfall frequencies
表6 2.0%和0.2%降雨頻率條件下泥石流堆積區(qū)危險(xiǎn)性分區(qū)統(tǒng)計(jì)表Table 6 Statistics of hazard zoning of debris flow accumulation area in 2.0% and 0.2% rainfall frequencies
模擬結(jié)果顯示:在2.0%降雨頻率條件下,泥石流總危險(xiǎn)面積為37 900 m2,牌嘴溝、拉路溝、義子溝無(wú)高危險(xiǎn)區(qū),加之牌嘴溝已修建道攔擋壩,危險(xiǎn)性不高。低危險(xiǎn)區(qū)占總危險(xiǎn)區(qū)的72.2%。扎龍溝高危險(xiǎn)區(qū)占比最高,為13.95%,主要分布于溝口,距離左岸居民建筑最近約30 m,威脅較小。但義子溝沿岸仍有4 戶居民處于中危險(xiǎn)區(qū)內(nèi),接哈溝溝口仍有5 戶居民處于中危險(xiǎn)區(qū)內(nèi)。
在0.2%降雨頻率條件下,泥石流危險(xiǎn)區(qū)總面積為60 500 m2,與百年一遇條件下相比,危險(xiǎn)區(qū)總面積僅增加了4.5%(2 600 m2);但中-高危險(xiǎn)區(qū)面積占比增大,其中接哈溝堆積區(qū)高危險(xiǎn)面積約4 000 m2,占總面積的34.78%;其次為義子溝和扎龍溝,高危險(xiǎn)區(qū)面積分別為3 100 m2和2 500 m2,占比分別為12.82%和35.82%??紤]到“5·10”后,在牌嘴溝已增設(shè)攔擋壩,牌嘴溝堆積區(qū)無(wú)高危險(xiǎn)區(qū)。除拉路溝外,其他各溝溝口和溝道附近居民均受到威脅,共約38 戶,其中中-高危險(xiǎn)區(qū)居民總戶數(shù)約22 戶。
(1)利用自動(dòng)站降雨觀測(cè)數(shù)據(jù),運(yùn)用FLO-2D 模擬重現(xiàn)了“5·10”泥石流運(yùn)動(dòng)和堆積特征,得到了泥石流堆積范圍、堆積深度、發(fā)展過(guò)程。與野外調(diào)查所得的現(xiàn)場(chǎng)情況進(jìn)行了對(duì)比,排除人為大規(guī)模改造作用,堆積扇均厚模擬結(jié)果誤差在-7.5%~23.3%之間,堆積扇面積模擬結(jié)果誤差在-6.2%~28.5%之間,數(shù)值模擬結(jié)果和實(shí)際情況基本相符,模擬效果較好。
(2)通過(guò)模擬,以泥石流堆積深度和沖擊速度為標(biāo)準(zhǔn),劃分了不同降雨頻率場(chǎng)景下的泥石流危險(xiǎn)區(qū);在此基礎(chǔ)上圈定了易于受泥石流威脅的房屋、財(cái)產(chǎn)等潛在受災(zāi)對(duì)象,為泥石流風(fēng)險(xiǎn)防控及治理提供了參考。
(3)FLO-2D 模型能夠較好地進(jìn)行泥石流過(guò)程重現(xiàn)和危險(xiǎn)性區(qū)劃,通過(guò)不同降雨頻率下的數(shù)值模擬,可以得到泥石流流體深度和速度的空間分布和變化過(guò)程,能夠?qū)δ嗍鬟^(guò)程及沖擊場(chǎng)景進(jìn)行直觀表達(dá)。