基于FLO-2D数值模拟震区急陡型泥石流冲出量——以瓦窑沟为例*

方群生，陈志和，唐川，徐惠梁

(1.中山大学地理科学与规划学院，广东 广州 510275； 2.中山大学土木工程学院，广东 广州 510275； 3.华南地区水循环与水安全广东省普通高校重点实验室，广东 广州 510275； 4.地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059； 5.中国建筑西南勘察设计研究院有限公司，四川 成都 610031)

汶川8.0级地震后，震区山体破碎严重，震后灾区地质环境变得更为脆弱，泥石流灾害处于新的活跃期[1-4]。汶川地震后，震区泥石流活动具有易发、频发、数量多、规模大、临界雨量阈值降低等特点，且多为黏性泥石流[5-7]。2010年“8.13”、2013年“7.10”，震区暴发了致灾更严重的群发性泥石流[8-9]。其中有一类运动、成灾规模较为典型的急陡沟道型泥石流，此类泥石流具有较小流域面积、较为陡峭的地形、较好物源条件、较好沟谷汇水条件、水动力大、冲蚀能力强、致灾严重等特征[10-12]。

目前，相关研究对急陡沟道型泥石流涉及较少，以美国科罗拉多州的泥石流形成的研究为背景，美国联邦地质调查局(United States Geological Survey, USGS) 发现在急陡沟道条件下，泥石流的物源起动中的“消防管效应”较为明显。Cannon等[13]和Godt & Coe[14]认为松散堆积物在山谷汇流作用下起动形成泥石流过程中会形成“消防管效应”，唐川等[15]提出在陡峭流域上游暴雨产生沟道径流导致水流快速集中，这个过程像“消防水管”效应。Yu等[16]通过野外实地调查选取台湾陈有兰溪流域的117 条泥石流，构建消防管效应统计模型。之后，狭陡型泥石流、窄陡沟道型泥石流、急陡型泥石流等对急陡沟道型泥石流的定义被相继提出[10,17-18]。陈明等[10]定义急陡沟道型泥石流为流域面积<5 km2、沟道平均纵比降>30%、流域完整系数<0.4 的区域，沟谷处于发育壮年期，汇水集中，侵蚀能力大，物源丰富，冲出量相对较大，地形条件上相似于沟谷型泥石流，形成、运动上又相似于坡面型泥石流的这一类泥石流。杨东旭等[17]定义狭陡沟道型泥石流为流域面积<3 km2、沟道平均纵比降>30%、流域完整系数<0.4 的区域，流域邻脊(谷) 狭窄、沟床纵坡陡峻、地表物源丰富的发展期山坡泥石流。韩玫等[18]提出窄陡沟道型泥石流为流域面积相对较小、纵坡陡、沟道狭窄，容易瞬间汇流形成突发性泥石流，具有流速快、沟口冲出规模大、迅速冲淤、迅速堵江而形成一系列重大危害的特点，时常导致沟口建筑物遭受毁灭性影响的泥石流。屈永平等[19]通过野外调查和遥感解译，定义了汶川强震区急陡沟道泥石流的平均纵比降为400‰的沟道泥石流。龚凌枫等[11]通过实验研究表明，急陡沟道泥石流堆积体的起动是急陡沟道泥石流冲出量巨大的主要原因。韩玫等[18]以磨子沟泥石流为例，采用CFX软件模拟再现50年一遇暴雨频率下窄陡型泥石流的动力学过程，分析其危险范围和破坏性影响。胡卸文等[20]研究震后泥石流的首期治理工程发现，对于沟道纵坡降较大、沟道狭窄、水动力作用较强的这类泥石流，以往的治理工程难以起作用，往后应有效地组合防治工程方案。陈明等[10]采用FLO-2D数值模拟瓦窑沟急陡型泥石流基本特征及沿途动力学特征，结果表明FLO-2D数值模型对急陡型泥石流仍具有较强的有效性和实用性。唐得胜[21]选取龙溪河的12条泥石流，采用FLO-2D进行自然和工程条件下泥石流危险性分区评价，结果显示，总冲出量越小防治效果越好。

目前，急陡沟道型泥石流的预警和工程控制难度较大，关于这类泥石流的研究虽已有所涉及，但数值模拟未考虑工程措施影响。本文以汶川震区映秀镇急陡型泥石流瓦窑沟为实例，结合野外调研数据，采用FLO-2D模拟4种不同降雨频率下的瓦窑沟泥石流在非工程和工程措施下的冲出规模，为深入了解汶川震区急陡沟道型泥石流规模、频率和危险性等提供参考。

1 典型案例——瓦窑沟概况

瓦窑沟流域汶川县绵虒乡高店村，沟口地理坐标31°20′5.39″N，103°28′58.97″E，沟口前缘有绵虒服务区配电房及国道G213经过。流域面积约1.21 km2，沟道纵向长度约2.78 km，平均纵比降约为48.9%，坡度一般为35°～65°。沟道地形陡峻，临空条件发育，且具有中上游宽，下游窄的特点。沟道的斜坡上部堆积体以松散的粉质粘土夹碎块石土为主，沟道两侧坡脚地堆积大量不稳定堆积体。2013年7月10日，在特大暴雨作用下，诱发了瓦窑沟泥石流，大量松散物质转化为泥石流快速冲出沟道。从7月10日上午9时至7月11日凌晨2时不断有泥石流冲出沟道，累计造成高店村9户民房完全被掩埋，18户民房部分被淤埋，绵虒服务区配电房部分被淤埋，停车场被淤积约3 000 m2，国道213线及老213线各被掩埋330 m。导致绵虒服务区停止运营，国道213线及老213线中断，因当地政府及时组织人员撤离所幸无人员伤亡，经济损失约1 000万元。本次泥石流灾害冲出量淤积在沟道出山口与河道之间的宽阔平坝上，泥石流堆积体前缘未到达岷江河道，未造成堵河。形成的堆积扇最大冲出长度约143 m，最大冲出宽度约307 m，堆积面积约3.9×104m2，累积3次冲出量约22×104m3，如图1。根据绵虒镇区域气象站资料，结合分析本次发生”7.10”泥石流的降雨量为20年一遇降雨强度(参考最近的羊店站点)，但考虑到前期连续降雨的影响及其叠加作用(达67.3 mm)，确定本次“7.10”泥石流属于50年一遇泥石流。

图1 瓦窑沟“7.10”泥石流冲出量淤积Fig.1 Accumulation fan of Wayao gully debris flow runout volume in the 7.10 event

2 数值模拟参数率定

2.1 数据来源

DTM是FLO-2D模拟的基础地形模型。依据瓦窑沟1∶10 000等高线的地形数据。模拟前，本文采用10 m×10 m的网格模拟，计算流域边界，并插值计算网格的高程点。工程情景模拟时，将实际的工程位置导入到地形数据，先创建坝体，并设置属性，赋值给坝体的高度、方向，导入防治工程参数(图2)。

防治工程模拟具体步骤：① 确定泥石流入流点的位置；② 通过数值高程模型、FLO-2D 用户手册、前人研究和野外调查准备数值模拟中所需的各种参数数据；③ 在入流点处设置流量过程线；④ 若高程数据为未发生泥石流或修建工程前的数据，需要对其进行校正，以便达到较好的模拟效果；⑤ 在 FLO-2D 软件中形成的地形模型中，设置防治工程相关参数；⑥ 进行模拟，获取结果。

通过野外的调查与测量，得到瓦窑沟流通区格栅坝的高、宽及设计库容量等坝体参数，坝体的高、宽及设计库容量分别为9 m，2 m，9.78×104m3。瓦窑沟沟口设置有排导槽工程，相对于流通区的格栅坝，其对泥石流冲出量的排导、消能、停淤的作用较小，因此本次数值模拟只考虑起主要拦挡作用的格栅坝，这也是泥石流淤积范围及淤积量与实际有一定差距的原因之一。

2.2 参数率定

结合FLO-2D 手册，对瓦窑沟数值模拟的参数进行率定。假设泥石流的降雨频率同于暴发频率，然后获取一次泥石流清水流量Qp，结合洪峰流量Qs，采用雨洪法获取瓦窑沟不同降雨频率下起始点的流量见表1。

1)重度及体积浓度

Fang等[7]在华溪沟泥石流沟口堆积扇上随机取样，采用现场配浆计算泥石流的容重，结合FLO-2D参考手册，综合选取体积浓度，模拟精确度约达89%。本文参考其方法，通过野外实地调查及现场试验，选取平均重度为2.3 g/m3，密度为2.65 t/m3，综合选取体积浓度为CV=60%。

图 2 FLO-2D 数值模拟软件中的坝体设计窗口Fig.2 Dam design window in FLO-2D numerical simulation software

2)屈服应力及粘滞系数

王纳纳[22]利用FLO-2D对汶川震区羊岭沟和码头沟冲出量进行模拟，根据震区泥石流特性，结合FLO-2D使用手册和詹钱登等[23]的实验，选取屈服应力(τy)及粘滞系数(η)的经验系数α1=0.811、α2=0.004 62、β1=13.72、β2=11.24，码头沟的模拟结果精确度达74%，羊岭沟的模拟结果精确度达 88%。由于本文选取的瓦窑沟位于汶川县绵虒乡，其地质地貌条件和泥石流特性与羊岭沟和码头沟相似，因此参考其方法取值。屈服应力(τy)及粘滞系数(η)的大小主要取决于流体的体积浓度，相应的求解如下：

η=α1exp(β1cν)

(1)

τy=α2exp(β2cν)

(2)

式中，α1、α2、β1、β2的取值为经验系数，取瓦窑沟的屈服应力及粘滞系数经验参数α1=0.811、α2=0.004 62、β1=13.72、β2=11.24。综合野外调查，因此参考其方法取瓦窑沟粘滞系数η=3 449.75。

3)层流阻滞系数(K)

根据FLO-2D手册，定义泥石流流域内稀疏植被所采用的层流阻滞系数为1 000～4 000，唐得胜[21]运用FLO-2D，选取层流阻滞系数K=2 280，对汶川震区茶马古道沟等12条泥石流沟进行模拟，模拟结果平均精确度达70.42%。由于瓦窑沟位于汶川震区，其地质地貌条件和泥石流特性相似，综合野外调查，因此参考其方法取瓦窑沟层流阻滞系数K=2 285。

4)曼宁粗糙系数

曼宁粗糙系数(c2)取决于王裕宜等[24]的公式(3)和FLO-2D手册，常鸣[25]根据王裕宜[24]公式和FLO-2D手册，结合现场野外调查确定泥石流堆积扇的堆积厚度分别计算出汶川震区黄央沟等12条泥石流沟的曼宁粗糙系数，模拟结果平均精确度达70.42%。参考其方法，本文根据现场调查，确定瓦窑沟流域内只存在无碎石分布的矮草原植被，且堆积扇的平均厚度约为3 m，故取值曼宁粗糙系数为0.1。

(3)

式中，Rns为体积浓度，h为泥深。

对比分析陈明等[10]设置降雨频率为2%，瓦窑沟泥石流在非工程措施条件下的数值模拟的参数，本文在曼宁粗糙系数略有差异，其他参数均保持一致，因此降雨频率为2%，非工程措施下的模拟结果和精度均非常接近。但本文同时还模拟了降雨频率为10%、20%、100%时，非工程措施条件下，瓦窑沟泥石流冲出量情况，以及考虑治理工程条件下，分别模拟瓦窑沟在不同降雨频率下的冲出量情况，对其进行工程措施影响效果分析。

沟道的侵蚀/堆积的分界点为本文模拟的入流点[26]，模拟时间按照泥石流实际历时确定，瓦窑沟泥石流按照“7.10”泥石流总历时约0.6 h。为对比非工程措施和工程措施下发生的泥石流，将两种情景下的泥石流入流点选取同一处。

3 急陡沟道型泥石流冲出量模型预测

3.1 非工程措施下模拟结果及分析

不同降雨频率会影响泥石流的冲出范围，数值模拟瓦窑沟在非工程措施下4种不同降雨重现周期(P=1%，2%，5%和10%)的泥石流冲出量。如图3所示，随着重现期的增加，冲出影响范围、平均深度、总冲出量均出现减少，重现期为10、20、50和100年一遇时，均不会堵塞岷江河道，但重现期超过50年一遇时，会对国道G213线及绵虒服务区造成不同程度的损坏，具有较大的危险性。这与瓦窑沟“7.10”泥石流事件的结果保持一致，重现期为50年一遇时，冲出量和冲出规模均不形成堵江。P=1%的冲出量分别是P=2%、P=5%的1.69倍、2.81倍，如表2所示。

3.2 非工程措施下模拟验证

根据野外调查和相关文献资料，瓦窑沟“7.10”泥石流事件是降雨频率为P=2%下暴发的泥石流，运用基于模拟堆积区域与实测冲出量范围叠加分区的精度评估模型[27]，对比验证实际发生的冲出规模与模拟结果，模拟的最大冲出长度、最大冲出宽度、冲出范围分别为156 m、259 m和2.77×104m2，实际发生的最大冲出长度、最大冲出宽度、冲出范围分别为143 m、307 m和2.9×104m2，冲出范围出现重叠的部分为2.45×104m2，模拟精度达到86.4%，模拟结果与实际发生泥石流事件的结果较为吻合。

3.3 工程措施下泥石流冲出量数值模拟

急陡型泥石流沟道修建了工程措施，比无工程措施下的泥石流冲出规模大很多。本文采用FLO-2D准确设置工程位置，以及将防治工程的拦挡效果反应出来。通过分析急陡沟道型泥石流在工程措施下的冲出规模，检验工程措施对泥石流的防治效应。

表1 瓦窑沟泥石流不同降雨频率起始点的流量表Table 1 Discharges of starting points in the catchments of Wayao gully at different rainfall frequency

表2 非工程措施下瓦窑沟泥石流的模拟结果Table 2 Simulation results of the Wayao gully debris flow under nonengineering measures scenarios

通过野外现场调查，瓦窑沟只在流通区修建了一道格栅坝的防治工程，坝后没有松散堆积物质，库容量非常大。拦挡坝主要作用是阻止泥石流进一步侵蚀流通区沟道，消耗泥石流能量，降低流速，减少冲出量冲出沟口。数值模拟时主要考虑这一道拦挡坝(图4)，其坝体有效高度为9 m，坝体淤埋深度为1.2 m，坝体轴线长度为20 m，坝体宽度为2 m，库容量为9.78×104m3。

模拟模型在非工程措施的数值模拟基础上加入一道拦挡坝，模拟结果显示，随着降雨重现期(P=1%，P=2%，P=5%，P=10%)的不断增加，冲出量不断减少，在重现期不超过P=5%时，防治工程防治效果更为明显。P=10%和P=5%时，拦挡坝挡住大部分泥石流，使其停淤在坝后。当重现期为P=2%和P=1%时，部分泥石流越过防治工程冲出沟口，模拟具体结果如图5和表3所示。

排导槽一般都修建于泥石流沟口位置，虽对泥石流掩埋的范围有较大的限制，会影响泥石流冲出量的长度和宽度，但排导槽对泥石流冲出量的排导、消能、停淤的作用较小。同时本文数值模拟工程措施时，瓦窑沟的冲出量在10年、20年、50年一遇降雨频率下并未冲出沟口，即排导槽工程没有影响到冲出量规模，而100年一遇降雨频率下冲出量大部分亦未冲出沟口，排导槽工程对冲出量的影响较小，因此本文没有考虑排导槽工程。

3.4 非工程与工程措施下模拟结果分析

相同降雨重现期，非工程措施下的冲出量远大于工程措施下的冲出量。当重现期为P=10%和P=5%时，由于防治工程的作用，泥石流未冲出沟口。非工程措施下，大量泥石流冲出沟口，冲入岷江，P=5%时，冲入岷江河道的泥石流物质变大，冲出物只淤积在沟道出山口与河道之间的宽阔平坝上，未形成堵塞河道，这与“7.10”瓦窑沟泥石流实际发生的冲出规模较一致。P=1%时，冲出量损毁部分国道G213道路，很大部分冲出量进入岷江河道。

P=2%时，在工程治理效果下，泥石流少量冲出沟口，比非工程措施下的冲出量减少6.47×104m3，减少83.16%，冲出影响范围减少2.1×104m2，减少75.81%，平均堆积深度减少0.85 m，减少30.25%。P=1%时，泥石流冲出量显著降低，比非工程措施下的冲出量减少8.39×104m3，减少63.70%，冲出影响范围减少2.01×104m2，减少50%，平均深度减少1.03 m，减少31.60%，未对国道G213道路形成威胁，且较少冲出量进入岷江河道(图3和图5)。工程措施的拦挡坝受泥石流物质堆积后，堆积深度加大，不断降低沟道纵坡降，致使泥石流冲出动能不断变小，有效控制了泥石流的冲出规模。

图4 瓦窑沟治理工程现场照片(摄于2018年1月8日)Fig.4 The site photos of treatment engineering in Wayao gully (Photographed on Jan 8, 2018 )

表3 工程措施下瓦窑沟泥石流的模拟结果
Table 3 Simulation results of the Wayao gully debris flow under engineering measures scenarios

P/%冲出量/(104m3)平均堆积深度/m堆积扇面积/(104m2)最大冲出长度/m最大冲出宽度/m14.722.232.0113221221.311.960.6789148

图5 工程措施下不同降雨频率的瓦窑沟泥石流冲出量模拟结果Fig.5 Simulation results of the runout volume in Wayao gully at different rainfall frequency under engineering measures scenarios

4 结 论

急陡型泥石流的形成、运动及成灾过程具有特有的规律，研究此类泥石流冲出规模，可以更深入掌握汶川震区急陡沟道型泥石流的风险、危害程度。通过对汶川震区的典型急陡沟道型泥石流瓦窑沟的现场调查分析，以及采用FLO-2D数值模拟非工程和工程措施两种情景下的泥石流冲出量，得到以下结论：

1)以汶川震区的典型急陡沟道型泥石流瓦窑沟为例，采用FLO-2D模拟重现期50年一遇时，瓦窑沟泥石流冲出量情况，模拟精度达到86.4%。模拟结果与瓦窑沟“7.10”泥石流事件的冲出规模基本保持一致，这表明模拟方法具有很高的准确性和有效性。

2)模拟重现期50年一遇时，工程治理措施下比非工程措施下的冲出量减少83.16%，冲出影响范围减少75.81%，平均堆积深度减少30.25%。模拟重现期100年一遇时，冲出规模显著降低，工程治理措施下比非工程措施下的冲出量减少63.70%，冲出影响范围减少50%，平均深度减少31.60%。