雙溝、灘山溝泥石流災害數值模擬及危險性評價

柳 軍，年 峰，鄒轉妮

(1.甘肅林業職業技術學院，甘肅 天水 741020； 2.甘肅省小隴山林業科學研究所，甘肅 天水 741020)

危險性指受災體接觸危險因子時遭受損失的可能性大小[1]。泥石流危險性評價是對泥石流發育的地質環境及災害可能造成的不良效應進行綜合判定。

19世紀后期，B.H.斯塔科夫斯基在土建工程中開展了泥石流危險性研究[2]。歐莫里[3]等對泥石流的研究為后來以規模和頻率為危險性評價因子奠定了基礎。戈麥斯[4]等利用GIS技術、結合人工智能，對委內瑞拉某地泥石流進行了危險性研究。安德魯[5]將GIS和RS結合，對泥石流堆積扇進行了危險性評價。譚炳炎[6]于1986年將影響泥石流溝的因素分為15個項目，根據每個項目的重要性分別賦值，計算泥石流的嚴重程度。劉希林[7]等提出了泥石流危險區劃的定量指標，結合灰色關聯度法評判了泥石流的危險性并通過實例進行了驗證。韋方強[8]等根據四川省泥石流數量眾多且地質環境條件復雜、無法用直接指標對其危險性進行精確區劃的特點，利用模糊信息模型，對泥石流危險性開展了區劃研究。寧娜[9]等采用將信息量模型和敏感性方法相結合的方式對武都地區的泥石流進行了危險性分析。徐惠梁[10]等將云模型和層次分析法相結合，對某一小流域泥石流的危險性進行了評價。牛全福[11]等利用灰色關聯度法、粗糙依賴度法和將兩者相融合的模糊綜合方法，對蘭州地區泥石流進行了危險性評價，結果表明可信度較高。辛聰聰[12]等利用單溝泥石流法，對芍藥溝泥石流進行了研究分析，確定了泥石流的危險程度等級。

基于有限差分數值方法結合GIS技術、利用高精度遙感影像資料，對雙溝、灘山溝泥石流進行了危險性評價，該流域面積小，模擬精度高，危險性評價效果顯著。

1 研究區概況

研究區為雙溝、灘山溝流域，地處麥積區琥珀鎮，屬渭河一級支流。自1998年以來，流域內發生過5次規模較大的泥石流災害及多次小規模泥石流，平均每兩年發生一次。2018年7月10日，流域內又一次暴發災害性泥石流，造成溝口排洪渠嚴重堵塞，隨后溢出并掩埋多處民房，毀壞大量通村公路、果園和耕地，造成了嚴重破壞。

流域全長1.12 km，面積1.11 km2，高差363 m。主溝平均縱坡降比70‰，溝道兩側山坡坡度30°～60°，呈“V”型溝谷，植被覆蓋率約40%，主要以喬木及草地為主。流域地勢陡峻，植被覆蓋率低，降雨、地表徑流、地下水豐富，固體物源松散且規模大，人類工程活動劇烈，具備暴發泥石流的有利條件。

2 研究目的與方法

2.1 研究目的

根據流域內常見的洪水重現期，利用有限差分數值模型，模擬雙溝、灘山溝流域不同頻率(P=1%、P=2%和P=5%)運動狀態，得到最終流速和流深，根據結果開展危險性評價，為當地泥石流災害防治提供參考。

2.2 研究方法

以有限差分數值模型對研究區進行數值模擬。該模型是奧布萊恩[13]利用非牛頓體物理原理和有限差分法數值計算建立的二維洪水災害數值模擬軟件，主要利用流體運動的質量守恒定律、動量守恒定律和連續性原理，建立固液二相流運動模型[14]，可以求解泥石流的流動速度、流動深度和堆積范圍，結合GIS技術，可進行泥石流危險性分析，開展危險性分區研究，從而為泥石流監測預警、預測預防和治理措施的制定提供有力支持。模型運動控制方程為：

連續性方程：

(1)

式中：h為泥深，單位m；I為水力坡降；u為水平方向平均流速，單位m/s；v為垂直方向平均流速，單位m/s。

運動方程：

(2)

(3)

式中：Sox和Soy分別為x和y方向的溝床坡降，Sfx和Sfy分別是x和y方向的摩擦坡降。

流變方程(據O′Brien-Julien模式[12]，以應力形式表示的流變方程)：

τ=τc+τmc+τv+τt+τd

(4)

式中：τ為泥石流流體總應力，單位Mpa；τc為黏聚性屈服應力，單位Mpa；τmc為庫倫剪應力，單位Mpa；τv為黏滯剪應力，單位Mpa；τt為紊流剪應力，單位Mpa；τd為擴散剪應力，單位Mpa。

在式(4)的基礎上，O′Brien和Julien提出了用坡降形式表達的流變方程：

(5)

式中：Sf為摩擦坡降；Std為紊流分散坡降；Sv為黏性坡降；Sy為屈服坡降；τy為屈服應力，單位Mpa；γm為流體比重，單位g/cm3；k為層流阻滯系數；h為泥深，單位m；η為流體黏滯系數；n為曼寧系數；v為流速，單位m/s。

3 數值模擬結果及分析

3.1 模擬參數選取

基于現場調查，并參考祝恩珍[15]在天水大砂溝泥石流研究中采用的方法，確定了本研究區體積濃度CV；利用軟件自帶的使用手冊確定了曼寧系數n；參考Woolhiser的研究成果[16]，確定了層流阻滯系數k；現場配漿確定了泥石流土石比重。選用的模擬參數如表1所示。

根據麥積區雨量站多年觀測數據，結合研究區地質災害實際，分別對研究區進行不同頻率(P=5%、P=2%、P=1%)數值模擬，得到自然條件下兩條泥石流溝的流速和流深(圖1)。自然條件下不同降雨頻率模擬結果見表2。

表1 研究區模擬參數取值表Tab.1 Simulation parameter value of the study area

由模擬結果可知，隨著降雨頻率降低，最終流速呈增大趨勢。不同頻率下，兩條溝的最終流速分布特征具有高度相似性，總體表現為從上游到下游逐漸降低，雙溝流速最大處位于形成區主溝道內，灘山溝則位于形成區主溝道及兩條支溝交匯處。兩條溝流通區溝道兩側及堆積區最大流速都小于1.5 m/s，泥石流做低速運動，此區域縱坡普遍小于50‰，溝谷開闊，泥石流勢能小，但動能還是較大的，具有較強的破壞力。兩條溝流通區主溝道最大流速為1.54～3.81 m/s，泥石流做中高速運動，此區域沒有住戶，但有一些耕地，因此會對當地居民的生產活動造成危害。3種不同頻率下都出現了溝口右上角堆積扇深入到居民區，對當地住戶的生產生活有一定的威脅，這與現場調查情況一致。2018年發生50年一遇的泥石流時，對溝口右上角村民房屋、300 m通村硬化道路及溝口105畝果園與耕地造成掩埋破壞，模擬結果基本符合現實情況。

圖1 自然條件下兩條泥石流溝的流速和流深Fig.1 Flow velocity and depth of two debris flow gullies under natural conditions

表2 自然條件下不同降雨頻率模擬結果Tab.2 Simulation results of different rainfall frequencies under natural conditions

隨著降雨頻率降低，最終流深呈增大趨勢。不同頻率下，溝道內最終流深表現為從上游向下游逐漸減小，堆積區后緣流深大而前緣逐漸減小的特征。雙溝最大流深出現在溝道內兩條支溝交匯處上游，灘山溝最大流深出現在主溝道形成區與流通區交匯處，主要原因是從上游急陡下沖的泥石流具有很高的動能，而此處溝道突然變緩開闊，大量泥石流動能轉化為勢能，流深增加，滑坡等松散物源加入，導致出現流深峰值。大量泥石流沿溝道傾瀉而下，到達溝口時縱坡接近水平面，泥石流便在溝口堆積，導致流深較大，后緣由于補給物源豐富而出現第二個流深峰值。在3種不同頻率下，溝口右上角沖入居民區的泥石流最大泥深分別為0.47 m、0.91 m、1.26 m，最大流速分別為0.5 m、1.3 m、1.5 m，可見隨著頻率降低，泥石流危害性增強。

3.2 自然狀態下泥石流危險性評價

根據數值模擬結果，考慮泥石流暴發頻率和強度之間的聯系[17]，分別將3種不同降雨頻率(P=1%，P=2%和P=5%)下的流速圖和流深圖加載到有限差分軟件中，得到2種不同工況下流域危險性分布圖，將雙溝、灘山溝劃分為高危險區、中危險區和低危險區。由于兩溝溝口緊臨，危險性分區存在重疊，故進行統一劃分(圖2)。

圖2 自然條件下泥石流危險性分區Fig.2 Zoning of debris flow risk under natural condition

自然條件下流域內高危險區、中危險區和低危險區分別為4 922 m2、7 358 m2和6 484 m2，各占流域總面積的0.44%、0.66%和0.58%。由于溝口地勢低矮、地形平坦，3種頻率下泥石流沖出范圍都比較大，直接沖入居民區，對當地居民正常的生產生活造成巨大危害，因此有必要對兩條泥石流溝進行防治。

4 結束語

利用有限差分數值模型對雙溝、灘山溝流域進行數值模擬，得到研究區3種不同降水頻率的最終流速和最終流深。結果顯示，隨著頻率降低，泥石流一次最大沖出量、堆積面積、流速和堆積厚度均呈增大趨勢。

基于數值模擬結果對研究區進行危險性評價，結果顯示，高危險區主要位于主溝道及堆積區左下角，面積4 890 m2；中危險區主要位于溝道兩側及下游堆積區，面積7 318 m2；低危險區主要位于下游溝道兩側邊緣及堆積區上緣部位，面積3 428 m2。總體而言，主溝道及堆積區危險性較高，有必要進行防治。