基于DAN3D的五花海泥石流运动过程模拟研究

盛 豪， 胡 卸 文, 2， 熊 冲 冲

(1.西南交通大学 地球科学与环境工程学院，四川 成都 610031; 2.西南交通大学高速铁路运营安全空间信息技术国家地方联合工程实验室，四川 成都 610031)

1 概 述

2017年8月8日，九寨沟(33.20°N，103.82°E)发生M7.0级强烈地震，地震诱发境内多处崩滑、滑坡、泥石流等次生地质灾害，严重破坏了九寨沟世界自然遗产的景观面貌[1-2]。五花海泥石流位于九寨沟核心景区，受地震影响，沟内发育多处崩滑、震裂物源，松散堆积体大大增加，为泥石流发育提供了充足物源条件。在暴雨条件下，泥石流一旦冲出，会对沟口五花海景观造成严重破坏。

国内外对于震后泥石流的研究很多，取得了较为丰硕的成果。Alexander J. Horton[3]利用Massflow对汶川震后红椿沟进行了运动过程反演后认为，利用voellmy铲刮模型能够较好适用于震后泥石流的模拟；胡旭东[4]选取了九寨沟震后五个集水区，综合采用泥沙增加量计算、粒子模型及数值计算的方法，对九寨沟震后泥石流危害性进行了评价；胡卸文[5-6]针对汶川震后泥石流冲出规模远大于规范计算结果，提出了震后泥石流普遍存在的一类新物源—震裂物源，并给出其启动模式；张向营[7]对地震扰动区泥石流启动阈值进行了探讨，认为对比震后泥石流启动激发雨量远远低于震前，仅为震前的三分之一左右；顾文韬[8]在调查统计了汶川震后多条泥石流物源基础上，总结了极震区震后泥石流形成过程及启动方式，得出了震后泥石流动储量与流域面积、高差的计算关系，为震后泥石流动储量估算及工程治理提供思路。综上所述，震后泥石流区别于常规泥石流，其具备临界雨量小、物源储量多、一次性冲出规模大等特点。因而，对其深入研究能够为治理工程设置及灾后恢复提供新的思路及方法。

笔者在对五花海泥石流野外调查的基础上，获取了地形地貌、地层岩性、物源条件及水源条件等相关形成条件及发育特征，计算了泥石流流速和一次冲出规模等动力参数。同时，运用DAN3D软件对震后泥石流运动过程进行模拟研究，对比分析两种结果，得出公式计算结果普遍偏小的结论，为防治措施提供了参考意见。

2 五花海泥石流发育特征

2.1 地形条件

五花海泥石流沟域形态呈栎叶形，无支沟发育，流域沟道纵向长度3.5 km，平均宽度0.85 km，流域面积2.6 km2，平均纵坡降434‰。流域内水系沿主沟道分布，最高点位于流域西侧山脊部位，高程3 981 m，最低点位于五花海，高程2 460 m，相对高差1 521 m，具有较显著的清水区、形成流通区和堆积区划分,见图1。

清水区平均纵坡475‰，沟谷呈深“V”字形，两侧纵坡陡峻，植被发育良好，少见崩塌滑坡等不良地质现象；形成流通区平均纵坡412‰，两侧岸坡不良地质现象发育，沟道内堆积震后垮塌滚落块石，受地震影响，树木多被折断；堆积区平均纵坡170‰，沟形不明显，植被覆盖良好，未见新近泥石流堆积。

2.2 物源条件

经野外调查，五花海泥石流沟域在地震前两侧斜坡植被发育较好，水土流失轻微，基本未见斜坡变形破坏现象。九寨沟“8·8”地震之后，五花海泥石流新增大量松散固体物源，沟内物源类型主要为崩滑物源、沟道堆积物源、坡顶震裂物源及坡面物源。震裂物源主要分布于沟道上游两侧山体顶部，主要为“8·8”地震产生。崩滑物源、沟道堆积物源主要分布于沟道中游区域沟道内和沟道右岸斜坡中下部。坡面物源大部分分布于谷坡中上部，储量一般。物源分布情况见图2。此外，在沟道中上游两侧坡面存在大量震倒树木，可能成为潜在漂木物源，一旦冲出，会对下游防护措施产生巨大冲击。据现场调查统计，崩滑物源总量约为31.2×104m3，沟道堆积物源总量约为2.5×104m3，坡面物源总量约为2.1×104m3，震裂物源总量约为25.4×104m3，共计松散固体物源量61.2×104m3，可能参与泥石流活动的物源量7.34×104m3(图3)。

图1 五花海泥石流主沟纵剖面

2.3 水源条件

五花海泥石流沟的水源主要来源于大气降水。鉴于既往研究，九寨沟区域泥石流均发生于雨季，春季冰雪融水一般不会成为引发泥石流的水源。此外，沟域内中上游地下水不丰富，不构成引发泥石流的主要水源，且沟域内没有水库、湖泊等集中的地表水体。因而，暴雨形成的地表径流是引发泥石流的主要水源，暴雨是泥石流的主要激发因素。泥石流沟域整体形态呈栎叶形，清水区沟道纵坡大，岸坡陡峻，有利于地表降水汇集，为五花海泥石流的发育提供有利的水源条件。

3 泥石流运动特征参数计算

泥石流运动特征参数是数值模拟求解设置和分析对比的依据，所需确定参数主要为泥石流重度、流速及暴雨洪峰流量、泥石流峰值流量，可采用雨洪修正法计算获得。

(1)泥石流重度。现场配浆法与查表法所得泥

图2 五花海泥石流平面图

图3 五花海泥石流物源统计直方图

石流重度综合取值15.1 kN/m3，属于稀性泥石流。

(2)泥石流流速。五花海泥石流属稀性泥石流，流速计算公式采用西南地区(铁二院陈光曦)公式。

(1)

式中Vc为泥石流流速；H为泥石流平均泥深；I为泥位纵坡降； 1/n为沟床糙率系数。

沟口附近泥石流流速计算结果见表1。

表1 五花海泥石流沟口流速计算结果

(3)暴雨洪峰流量和泥石流峰值流量。泥石流峰值流量采用雨洪修正法进行计算。

QC=Dc(1+φ)Qp

(2)

Qp=0.278ψAS/tn

(3)

式中QC为泥石流断面峰值流量(m3/s)；φ为泥沙修正系数；Qp为暴雨洪峰流量；Dc为堵塞系数；ψ为洪峰径流系数；S为单位历时的暴雨平均强度(mm/h)；A为流域面积(km2)；t为汇流时间(h)；n为暴雨递减指数。

沟口泥石流峰值流量计算结果见表2。

表2 五花海泥石流沟口流量计算结果

4 DAN3D数值模拟

4.1 DAN3D软件

DAN(Dynamic Analysis)方法由Hungr(1995)[9]提出，最早为解决高速远程滑坡运动过程复杂、组成物质多样、传统理论模型无法很好适应等问题。该方法的核心思想是等效流体假设与连续性假设。DAN模型将运动中的复杂非均值的滑体或泥石流浆体视为某种理想的、具有简单的流变关系的流体，并将流体视为连续性介质，以连续的场代替小范围内离散物质相互作用。因而，物质的变形及运动仅取决于内部及底部的流变准则。目前国内利用DAN3D对高速远程滑坡、碎屑流等灾害进行模拟反演已经取得了比较好的成果，如深圳“12.20”渣场滑坡、汤家沟滑坡～碎屑流等案例[10-11]。国外采用DAN3D对泥石流进行模拟的成功案例比较多，如皮茨河谷地区泥石流、上阿诺谷泥石流[12-14]等案例。

4.2 Voellmy模型的基本原理

DAN3D内置Frictional、Plastic、Newtonian、Bingham及Voellmy五种流变模型，使用者可根据实际问题选择恰当的模型进行数值模拟分析。多项研究结果表明，对泥石流进行动力特性分析的时候，Bingham及Voellmy模型能够较好反映泥石流运动历史。由于五花海泥石流为稀性泥石流，Voellmy模型更能考虑到流体运动过程中的湍流效应与速度效应。因而，本文采用Voellmy模型进行模拟。

Voellmy模型认为，流体在运动过程中所受到的抗剪应力应为流体所受摩擦力和湍流流动共同产生的阻力之和，其模型公式如下所示：

式中τ为滑体所受抗剪应力；f为摩擦项，等同于tanφ；ξ为紊流项，与流体速度和密度相关。

5 五花海泥石流数值模拟

5.1 参数设置

利用五花海泥石流1∶1000地形图进行沟域DEM三维地形建模，并采用Kriging插值法进行网格化处理，生成沟域网格模型。同时，根据现场调查情况，将物源厚度建立与沟域网格模型相对应的网格模型，生成的DEM模型见图4、5。

图4 五花海泥石流沟域DEM模型

本文采用Voellmy摩擦模型作为五花海泥石流流体运动准则，其中最重要的参数为摩擦系数μ、湍流系数ξ及泥石流流体容重。摩擦系数及湍流系数的取值范围，在前期的许多文献中有过研究。一般认为，Voellmy模型中摩擦系数为0.05～0.8，湍流系数30～800 m/s2。因此，本文在前人研究基础上选取不同参数进行模拟，最终确定摩擦系数选用0.3，湍流系数选用50 m/s2。

图5 五花海泥石流物源DEM模型

5.2 模拟结果分析

根据上述设定参数，利用DAN3D对五花海泥石流在P=2%及P=1%暴雨时的运动情况进行数值模拟，得到泥石流堆积厚度模拟结果，见图6。

图6的模拟结果表明：50 a一遇泥石流发育过程中，崩滑物源首先启动，向坡脚沟道堆积；位于分水岭的震裂物源渐次运动，进入沟道或沿坡面向下运动，形成沟道物源及坡面物源，参与到泥石流活动中。从t=30～150 s这段时间的云图可以看出，泥石流运动到沟道宽缓(海拔2 675～2 710 m段)位置时，流体逐渐停止运动，在沟道内堆积，仅部分能够冲出继续向下游流动。t=300 s的云图可以看到，沟口有泥石流堆积，但是方量很小。t=300～600 s云图可以看出，这段时间泥石流冲出量极少，小于公式计算结果，前缘漫流至五花海栈道，对栈道及游客安全构成一定威胁。同时，从图8可以看出：泥石流在沟口最大速度也小于公式计算结果。因此，可以推测，海拔2 675～2 710 m处的宽缓平台为泥石流起到了天然停淤作用，很好地减缓了泥石流冲出速度，滞留了大量固体物质。

图7的模拟结果表明：100年一遇泥石流运动及物源启动方式与50年一遇的情况相一致。但是由于泥石流流经宽缓沟道时流速较大，因而，大量泥石流进一步向下游冲出，逐渐堆积形成扇形堆积扇。根据图9可以看出，泥石流在沟口最大速度约为6.5 m/s，与公式计算结果符合较好。泥石流一次冲出总量为5.7×104 m3，远大于公式计算结果，与震后泥石流冲出规模远远大于常规泥石流的认知一致。

图6 五花海泥石流P=2%堆积厚度云图

图7 五花海泥石流P=1%堆积厚度云图

5.3 影响范围预测及防治建议

五花海泥石流沟口正对五花海及栈道，因此，利用模拟结果对泥石流冲出危害进行一定的预测，在后期对五花海泥石流提出针对性的建议很有必要。50年一遇工况下，虽然物源启动量很大，但大部分在流经宽缓沟道处便开始减速，并逐步停止运动，仅少量冲出沟口，漫流至栈道处。此时，泥石流流速不大，最大约3 m/s，对栈道会构成一定威胁，但对五花海景观影响不大；100年一遇工况下，尽管有部分泥石流在宽缓沟道处停淤，但仍有大量泥石流冲出沟口。堆积扇前缘延伸至五花海，沟口最大速度可达6.5 m/s，可能冲毁栈道。同时，大量泥沙可随之裹挟进入五花海，对五花海景观造成毁灭性打击。

图8 P=2%泥石流最大速度云图

图9 P=1%泥石流最大速度云图

根据模拟结果可知：五花海泥石流爆发主要物源来源于崩滑物源，可采取支挡措施避免崩滑物源垮塌进入沟道参与泥石流活动；同时，泥石流在宽缓沟道处会大量堆积，可在此设置停淤围堰，提高其停淤能力；为避免泥沙进入五花海，对五花海水质景观造成破坏，可在临近沟口处设置拦水堤，利用天然土体对含砂水流进行过滤。

6 结 语

五花海泥石流位于九寨沟核心景区境内，沟域形态呈栎叶形，流域面积为2.6 km2，主沟长0.85 km，主沟平均纵坡降434‰。“8·8”地震后，沟域内地质灾害发育松散物源丰富，为泥石流的形成提供了地形、物源和水源条件。

结合泥石流基本发育特征和动力学参数计算，可得出九寨沟泥石流容重15.1 kN/m3，50 a、100 a一遇降雨频率下沟口流速分别5.68 m/s、6.21 m/s，一次过流总量为0.92 m3、1.11 m3。

通过DAN3D软件对泥石流进行数值模拟，预测了泥石流50 a与100 a一遇的运动过程，可以得知50 a一遇工况下泥石流不会大量冲出，少量冲出泥石流会对栈道构成一定威胁。100 a一遇工况下，泥石流浆体大量涌出，泥沙大量冲入五花海，会对五花海景观水质构成严重威胁。

根据模拟结果，提出了“拦挡物源、围堰停淤、拦水过滤”的防治方案。