基于Massflow的达摩沟泥石流模拟与危险性分析

杨鑫 孟华君 路璐 王宇 吴季寰

摘要：

达摩沟小流域历史上曾发生泥石流灾害事件，造成较大人员伤亡和财产损失。目前，该小流域主沟堆积大量煤矸石，主沟两侧坡面堆积有较多崩、坡积物。为分析泥石流再次发生的可能性与危险性，在无人机摄影测量和野外现场调查的基础上，详细分析达摩沟小流域泥石流发育条件，基于改进MacCormack-TVD有限差分法，采用Massflow软件对10，20，50，100 a一遇4种降雨概率条件下达摩沟暴发泥石流的工况进行数值模拟，得到达摩沟小流域在不同降雨概率下的泥石流危险性。结果表明：①  在此4种降雨概率下，大部分区域的泥石流流速维持在0～3 m/s，泥深保持在0～4 m的水平，而峰值流量最大可达19 m/s，泥深最大可达3.85 m；②  在100 a一遇的降雨强度下，达摩沟高危险性地区占16.54%，中风险地区占49.27%，低风险地区占34.19%；③  模拟结果较好地符合泥石流的时空发育特征，展现了降雨强度对泥石流运动的影响规律，可为泥石流的预警监测提供参考。

关 键 词：

泥石流； 危险性； Massflow； 数值模拟； 达摩沟

中图法分类号： P642.23

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2023.08.020

0 引 言

1867～1999年间，北京地区共发生超过50次泥石流灾害事件，造成2 400多人死亡，经济损失达数亿元［1］。近年来，随着极端天气频发和人类工程活动增强，北京地区先后发生2012年“7·21”和2018年“7·16”泥石流灾害群发事件［2］，虽处置得当，避免了较大人员伤亡，但仍造成了较大财产损失和社会影响。

长期以来，大量专家学者基于野外调查、监测预警等取得的一手资料，对北京市山区泥石流进行了大量研究。在20世纪80年代，洪惜英、李容全、张仲德等就开始研究北京地区泥石流的发育特征、分布规律等基础性工作［3-5］，旨在通过对泥石流灾害的科学认识，有效降低地质灾害造成的经济、财产损失，提出防治对策、保障社会安全。20世纪90年代，相关研究工作已经从泥石流事件的现场调查和资料统计分析逐渐转向成因机理分析及防治对策设计［6-7］和泥石流危险性区划［8-10］，并尝试开展北京地区泥石流预警预报工作［11-15］。21世纪后，随着仿真模拟技术、信息传输技术的快速发展，基于泥石流启动物理过程的泥石流灾害精细化风险管理成为可能，不断出现精度更高的预警模型。例如，白利平等［16-17］根据力学平衡条件，尝试推出松散固体物质起动判别式和雨强表达式。赵忠海［1］综合考虑地形地貌、地质条件、土壤类型以及降雨情况等因素，建立了泥石流是否发生的判断公式。王海芝［18］通过整理1949年以来北京地区发生的泥石流与其发生时期的降水数据，建立了临界雨量的定量估算I-D模型等。这些工作为北京市泥石流灾害管控提供了理论指导和技术支撑。然而，随着京津冀协同一体化国家发展战略的实施，北京市对泥石流灾害管控提出了更高的要求，逐渐由灾害点管控向“隐患点+风险区”双管控转变。新的管控策略客观上要求对隐患点泥石流形成、运动及承灾机理进行深入研究［19］，而数值模拟方法是解决该问题的较好手段。通过数值模拟，可以再现泥石流运动演化过程，得到用于定量评价泥石流有关动力学、运动学要素和防治工程的参数［20］。

目前，国内外开展泥石流数值模拟主要采用FLO-2D、Fluent、OpenLISEM及PFC等成熟商业软件，研究结果显示数值模拟方法在模拟泥石流运动过程、冲出量、堆积等方面具有较好的效果［21-24］。鉴于此，本文采用数值模拟方法，选择北京市达摩沟作为研究区，利用更专业的泥石流模拟软件Massflow模拟达摩沟泥石流运动过程，通过研究泥石流隐患点在不同降雨工况下的泥石流运动及沉积过程，厘清达摩沟泥石流危险性，为该小流域泥石流隐患点风险管理提供依据。

1 研究区概况

达摩沟位于北京市门头沟区清水镇东南方向 3 km 处，是清水河支沟，位于东经115°37′～115°39′，北纬39°52′～39°56′（见图1）。该泥石流沟沟域面积为3.96 km2，相对高差647 m，主沟沟向325°，断面大多呈“V”形，长2 680 m，沟宽30～100 m，平均15 m，沟床平均坡降230%。陡坎、卡口不多，河段宽窄均匀，主支沟交角大，沟道较顺直平整。沟床两岸坡面较陡，坡度30°～60°，平均40°。植被覆盖率为70%。

达摩沟内出露地层有侏罗系下统窑坡组（J1y）、龙门组（J1l），侏罗系中统九龙山组（J2j）、髫髻山组（J2t）及第四系全新统（Q4）。岩性主要为砂岩、粉砂岩、安山岩、玄武岩、凝灰质砂岩、粉砂岩、页岩及煤系地层等，第四系全新统主要由砾石、岩屑、黏土质等残、坡积物组成［25］。

研究区地处北京西部山区，属中纬度大陆性季风气候，年平均气温10.2 ℃，极端最高气温37.6℃，极端最低气温西部-22.9 ℃。降水量自东向西逐渐减少，受中纬度大气环流的不稳定性和季风影响，降水量年际变化大。

达摩沟流域历史上曾多次暴发泥石流，最近一次发生于1950年8月4日，当日最大降雨量達190.1 mm，最大雨强为56.5 mm/h。持续的暴雨造成达摩庄后港沟和王家港沟泥石流爆发，造成84人死亡、24人重伤，基础设施损毁严重，给当地人民造成了大量的财产损失［26］。

2 泥石流形成条件

2.1 物源条件

（1） 受达摩沟历史煤矿开采影响，主沟及支沟沟道堆积大量煤矸石（见图2（a）～（b）），其中仅少部分煤矸石堆有浆砌石护坡（图2（a）），绝大部分煤矸石堆没有任何防护措施，稳定性较差，是达摩沟未来暴发泥石流的主要固体物源。

（2） 沟道两侧岸坡浅表层堆积大量坡积物。这些松散堆积物主要分布于坡面中下部浅表层，厚度不等，约30～80 cm（见图2（c）～（d）），植被覆盖较好。然而，根系延伸同时也使堆积更加松散，空隙发育，加之下伏基岩具有倾向主沟的优势结构面（图2（d）），当降雨量较大时，雨水在缓慢入渗的同时极易沿大空隙灌入坡积物，抵达土岩接触面，在该面形成积水，从而引起崩坡积物失稳补给泥石流。图2（c）～（d）可见，坡面下部已经缺失了一部分浅表堆积物，下伏基岩面直接出露。

（3） 在坡脚分布有零星崩滑堆积体（图2（e））和历史泥石流堆积体（图2（f））。其中，泥石流堆积体可见明显的定向排列，大颗粒石块长轴方向基本一致。这些堆积体规模大小不等，固结程度较高，不易失稳，但在极端降雨工况下，受上游准泥石流体冲击侵蚀，会受到强烈侧蚀作用，可能作为泥石流补给。

2.2 水文条件

达摩沟泥石流一般由强降雨激发。该地区降雨分布具有年际间不均衡、年内变化大、降雨集中、降雨强度大等特点。该沟一般年降水量为550～650 mm，平均为575 mm，汛期平均雨量为400～500 mm；历史年最大降水量932 mm，年最小降水量329 mm。最大日降雨量为224.8 mm，最大小时降雨量达56.5 mm。每年雨季降水量占全年降水量的70%左右，且多以大暴雨或特大暴雨形式出现。

3 研究方法

数值模拟是解决单沟泥石流隐患点危险性评价的有效方法。本文采用专业泥石流模拟Massflow软件对研究区泥石流进行运动过程模拟。

Massflow是中国科学院成都山地灾害与环境研究所欧阳朝军研究员研发的一款数值模拟软件。该软件基于深度积分连续介质力学理论，利用改进MacCormack-TVD有限差分方法，兼顾复杂地形地貌，是具有二阶精度和自适应求解域特征的山地灾害动力学高效计算模拟软件。与国内外成熟软件相比，该软件不仅提供了丰富的物理模型以提高模拟准确性，还能够极大提升计算效率，缩短模拟时间。在甘肃岷县二马沟泥石流模拟［27］、金沙江白格滑坡模拟［28］以及泥石流防治工程防治效果模拟［29］中，Massflow均取得了良好的效果。

3.1 Massflow模型控制方程

4 模拟结果与分析

4.1 计算模拟结果

通过对10%，5%，2%，1%降雨频率下达摩沟爆发泥石流的过程进行模拟，可以得到不同降雨条件下达摩沟泥石流的运动过程和不同时刻的泥深以及流速分布，并在此基础上开展达摩沟泥石流隐患点在不同降雨频率下的危险性分析评价。部分模拟结果如图3，4所示。

4.2 模拟结果分析

4.2.1 流速变化

根据计算结果（见表2），不同降雨工况下达摩沟泥石流流速变化具有以下3个特点。

（1） 总体流速较为稳定。在10，20，50，100 a一遇降雨概率下，整个主沟道流速总体普遍维持在0～4 m/s的水平；只有下游的阶地以及陡峭区域流速急剧增大，局部平均流速达到8 m/s以上（见图3）。

（2） 峰值流速随降雨强度增加，呈先增大后逐渐趋于稳定的趋势，见图5。其中，10 a一遇和20 a一遇降雨频率下的峰值流速差异较大，而20 a，50 a，100 a一遇降雨的峰值流速差异较小，总体维持在18～20 m/s的水平。

（3） 随着降雨强度增大，泥石流高速流动区的距离明显变长。对比50 a一遇和100 a一遇降雨工况下高速流动区的沟道长度，发现100 a一遇降雨条件下的峰值流速略小于50 a一遇降雨条件下的峰值流速，但其高速流动区更长。

以上3点表明：①  达摩沟泥石流峰值流速受限于降雨强度和沟道地质条件两个方面。当降雨强度小于20 a一遇时，峰值流速会随降雨强度增加而增大；当降雨强度达到20 a一遇时，降雨强度继续增大，并不会继续提高峰值流速，而会延长高速流动的距离。②  降雨强度增大导致泥石流流速和高速流动距离增加，本质上是提高了泥石流体在流动过程中的冲击速度和冲击距离，即提高了侵蚀效率，大大提高了泥石流体的破坏能力。这既有利于沿程物源进一步补给泥石流体，又进一步增加了泥石流体的重度，提高了破坏力。

4.2.2 运动过程泥深变化

泥石流堆积深度是评估泥石流危险性的重要指标。图4，6是不同降雨工况下泥石流堆积情况统计趋势。总体来看，达摩沟泥石流堆积特征主要可以概括如下：①  泥石流最大泥深随降雨强度增加呈先增大后减小的趋势；②  达摩沟在不同降雨条件下均不会形成明显的堆积扇；③  泥石流在沟道狭窄处泥深将会增大，偏向流通区上游的部分堆积深度较小，而流通区中下游和堆积区等泥石流前缘部分堆积深度较大。

造成以上泥石流堆积特征的原因主要有：①  在降雨强度较低的情况下，水流对堆积物的侵蚀能力和搬运能力均较弱，泥石流流量较小；随着雨量的增加，产汇流增加，水流的侵蚀搬运能力逐渐增强，泥石流流量变大，导致泥石流的堆积厚度增大；当降雨强度增加到一定程度后，泥石流流量增加、流速变大，其裹挟能力更强，更容易将固体物质搬运至更远的地方，因此高强度降雨条件下沟道内的最大堆积厚度反而有所减小。②  达摩沟狭窄的沟道不具备形成堆积扇的地形条件。主沟下游地势相对平坦，有利于泥石流运动减速开始堆积，但历史泥石流在流通区下游和近沟口处有堆积，导致目前流通区末端和近沟口处地势仍然相对较高，不利于泥石流淤积蠕动形成扇体，而有利于未来泥石流活動进一步向清水河流动；此外，在流通区下游以及堆积区南侧有部分居民区并修建有道路，导致达摩沟主要为带状堆积。③  达摩沟与清水河交接处沟道比降仍较大，不利于泥石流减速沉积，更可能出现的情况是泥石流体高速冲入清水河，因此，达摩沟沟口没有明显的堆积扇。④  由于沟道较窄处泥石流的运动受到阻碍，泥石流流速变缓，导致其更易在此处堆积下来。⑤  在泥石流堆积过程中，随着泥石流不断向周围扩散堆积，其流速、动能不断降低，泥石流边运动边堆积，最终整体停止运动，所以堆积区中心厚度最大，而发散区域堆积厚度逐渐降低。

4.2.3 危险性评价

泥石流的强度可以依据泥石流的流速与泥深的强度进行划分。本文参照常鸣［32］提出的泥石流强度影响以及危险性分区标准对模拟结果进行危险性评价。泥石流强度影响划分标准如表3所列。

根据表3，利用Arcgis对泥石流流速与泥深进行栅格叠加计算，将计算结果分类为高危险、中危险、低危险，最终得到100 a一遇降雨工况下的达摩沟泥石流危险性分区如图7所示。

在图7中，高危险性地区占16.54%，中危险性地区占49.27%，低危险性地区占34.19%。低中危险地区主要分布在主沟边缘地带，主沟绝大部分区域处于中高危险区。考虑到达摩村位于主沟沟口处，涉及到大量人、财、物的潜在损失，说明达摩沟是一条风险性较大的泥石流沟，因此需要加强对达摩沟的预警监测。

5 结论与建议

（1） 达摩沟小流域在强降雨条件下有发生泥石流的可能，模拟结果显示：泥石流流速主要分布在0～4 m/s，随着降雨强度的增大，泥石流平均流速逐渐增大，但泥石流峰值流量随着降雨强度的增大呈现先增大后减小的趋势。达摩沟泥石流堆积泥深主要分布在0～3 m之间，最大泥深可达3.85 m，平均泥深随着降雨强度的增大而增大，峰值泥深则呈现先增大后减小的趋势。

（2） 依据泥石流流速与泥深，利用Arcgis对达摩沟泥石流危险性进行分区，结果显示：高危险区占16.54%，中危险区占49.27%，低危险地区占34.19%。然而，中高危险区仍有居民居住，因此风险较高，建议加强监测。

（3） 泥石流流速及泥深与降雨强度的统计关系表明：泥石流运动不是简单受降雨强度影响，还受到地形条件、植被覆盖等多因素约束，但是不同工况下各因素贡献有差异。因此，在对泥石流进行预警监测时，需针对不同小流域进行具体分析，建立不同的监测方案。

参考文献：

［1］ 赵忠海.北京泥石流灾害预警研究现状与技术方法探讨［J］.城市地质，2015，10（增1）：78-84.

［2］ 北京市地质研究所，中国地质科学院地质力学所，北京林业大学.北京市沟道侵蚀泥石流多元预警阈值模型研究及示范工程［R］.北京：北京市地质研究所，中国地质科学院地质力学所，北京林业大学，2021.

［3］ 洪惜英，王礼先，朱靖才，等.北京市密云县汗峪沟泥石流［J］.北京林学院学报，1983（4）：38-45.

［4］ 李容全，徐振溥.北京山区泥石流分布特点［J］.山地研究，1983（3）：42-48.

［5］ 张仲德.北京山区泥石流分布和成因分析［J］.北京师院学报（自然科学版），1984（2）：77-89.

［6］ 魏永明.1991年北京北部山区泥石流成因及防灾对策［J］.山地研究，1993（3）：187-192.

［7］ 张世清，毕小刚，刘振国.北京山区泥石流防治对策探讨［J］.水土保持通报，1992（3）：46-51.

［8］ 张春山.北京北山地区泥石流灾害危险性评价［J］.北京地质，1996（2）：11-20.

［9］ 張春山.北京地区泥石流灾害危险性评价［J］.地质灾害与环境保护，1995（3）：33-40.

［10］ 赵波，董桂芝.北京西山清水河流域泥石流及危害预测［J］.北京地质，1996（3）：1-14.

［11］ 谢又予，崔之久.北京山区泥石流的防治与预报［J］.贵州科学，1992（3）：132.

［12］ 吴正华，储锁龙.北京泥石流暴雨基本特征［J］.大气科学，1992（4）：476-481.

［13］ 韦京莲，董桂芝，赵波.北京山区近代泥石流活动规律及暴发周期分析［J］.中国地质灾害与防治学报，1994（4）：48-53.

［14］ 韦京莲，赵波，董桂芝.北京山区泥石流降雨特征分析及降雨预报初探［J］.中国地质灾害与防治学报，1994（增1）：45-52

［15］ 吴正华.北京泥石流灾害及其降水触发条件［J］.水土保持研究，2001（1）：67-72.

［16］ 白利平.北京山区泥石流灾害临界雨量研究［J］.地质灾害与环境保护，2006（4）：101-104.

［17］ 白利平，孙佳丽，南赟.北京地区泥石流灾害临界雨量阈值分析［J］.地质通报，2008（5）：674-680.

［18］ 王海芝.北京地区暴雨泥石流预警阈值研究［J］.第四纪研究，2020，40（5）：1371-1380.

［19］ 黄晓虎.泥石流降雨启动机制及早期预警模型的研究：以达摩沟为例［D］.长春：吉林大学，2016.

［20］ 涂剑，马超，杨海龙.北京山区暴雨泥石流激发雨量条件［J］.中国水土保持科学，2017，15（5）：103-10.

［21］ 方群生，陈志和，唐川，等.基于FLO-2D数值模拟震区急陡型泥石流冲出量：以瓦窑沟为例［J］.中山大学学报（自然科学版），2020，59（3）：23-31.

［22］ 赵峥，刘铁骥，陈亮，等.基于Fluent强震区王家沟泥石流运动特征研究［J］.灾害学，2021，36（2）：208-12.

［23］ 陈明，唐川，王飞龙，等.基于OpenLISEM模型的泥石流启动冲出预测［J］.泥沙研究，2020，45（4）：59-65.

［24］ 胡明鉴，汪稔，陈中学，等.泥石流启动过程PFC数值模拟［J］.岩土力学，2010，31（增1）：394-397，434.

［25］ 慎乃齐，韦京莲，赵波，等.京西达摩沟泥石流形成条件与防治［J］.探矿工程（岩土钻掘工程），2003（增1）：96-98.

［26］ 齊干，张长敏.达摩沟泥石流形成的物质条件分析及防治对策［J］.水文地质工程地质，2011，38（5）：102-109.

［27］ 乔渊，刘铁骥，陈亮，等.基于Massflow模型的甘肃省岷县二马沟泥石流危险性评价［J］.水利水电技术，2020，51（4）：184-192.

［28］ 曹水和，吴新明，钟东.基于Massflow的金沙江滑坡数值模拟及堵江风险预测［J］.地质灾害与环境保护，2021，32（4）：3-7.

［29］ OUYANG C，HE S，XU Q，et al.A MacCormack-TVD finite difference method to simulate the mass flow in mountainous terrain with variable computational domain ［J］.Computers & Geosciences，2013，52：1-10.

［30］ HUNGR O，MCDOUGALL S.Two numerical models for landslide dynamic analysis ［J］.Computers & Geosciences，2009，35（5）：978-92.

［31］ 周凯琦.基于单层单相流和多层多相流的灾害动力学数值模拟［D］.武汉：华中科技大学，2017.

［32］ 常鸣.基于遥感及数值模拟的强震区泥石流定量风险评价研究［D］.成都：成都理工大学，2014.

（编辑：高小雲）

Abstract：

Mudslide disaster events had occurred in the Dharma Gully sub-basin area historically and caused a large casualties and property losses.At present，a large amount of coal gangue was piled up in the main ditch of the sub-basin，and more landslide deposits were accumulated on both slopes sides of the main ditch.In order to analyze the possibility and risk of mudslides reoccurring，this paper analyzed the development of mudslides in Dharma Gully sub-basin area in detail through UAV photogrammetry and field investigation.Based on the improved MacCormack-TVD finite difference method，the Massflow software was used to simulate the mudslide conditions and obtained the mudslide risk in Dharma Gully under 4 rainfall probabilities of 10years，20years，50years，and 100years.The results showed that： ① Under the above 4 rainfall probabilities，the flow velocity of debris flow in most areas was 0～3 m/s，the mud depth was 0～4 m.The maximum flow rate could reach 19 m/s，and the maximum mud depth could reach 385 m；② Under the rainfall of 100-year，the high-risk areas in Damogou accounted for 1654%，the medium-risk areas accounted for 4927%，and the low-risk areas accounted for 3419%；③ The simulation results reproduced the temporal and spatial development characteristics of debris flow well，showed the influence of rainfall intensity on the movement of debris flow，which could provide some references for the early warning and monitoring of debris flow.

Key words：

debris flow；risk；Massflow；numerical simulation；Dharma Gully