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

第54卷第8 期2023年8 月文章编号：10 0 1-417 9(2 0 2 3)0 8-0 146-0 7引用本文：杨鑫，孟华君，路璐，等.基于Massflow的达摩沟泥石流模拟与危险性分析J.人民长江,2 0 2 3,54（8）：146-152.人民长江YangtzeRiverVol.54,No.8Aug.，2 0 2 3基于 Massflow 的达摩沟泥石流模拟与危险性分析杨鑫1-2,孟华君1.3.4,路璐,王字,吴季裳（1.中国地质科学院地质力学研究所，北京10 0 0 8 1；2.北京科技大学土木与资源工程学院，北京10 0 0 8 3；3.自然资源部新构造运动与地质安全重点实验室，北京10 0 0 8 1；4.自然资源部陕西宝鸡地质灾害野外科学观测研究站，北京10 0 0 8 1；5.北京市地质研究所，北京10 0 0 8 1）摘要：达摩沟小流域历史上曾发生泥石流灾害事件，造成较大人员伤亡和财产损失。目前，该小流域主沟堆积大量煤研石，主沟两侧坡面堆积有较多崩、坡积物。为分析泥石流再次发生的可能性与危险性，在无人机摄影测量和野外现场调查的基础上，详细分析达摩沟小流域泥石流发育条件，基于改进MacCormack-TVD有限差分法，采用Massflow软件对10,2 0,50,10 0 a一遇4种降雨概率条件下达摩沟暴发泥石流的工况进行数值模拟，得到达摩沟小流域在不同降雨概率下的泥石流危险性。结果表明：在此4种降雨概率下，大部分区域的泥石流流速维持在0 3m/s，泥深保持在0 4m的水平，而峰值流量最大可达19m/s，泥深最大可达3.8 5m；在10 0 a一遇的降雨强度下，达摩沟高危险性地区占16.54%，中风险地区占49.2 7%，低风险地区占34.19%；模拟结果较好地符合泥石流的时空发育特征，展现了降雨强度对泥石流运动的影响规律，可为泥石流的预警监测提供参考。关键词：泥石流；危险性；Massflow；数值模拟；达摩沟中图法分类号：P642.230引言18671999年间，北京地区共发生超过50 次泥石流灾害事件，造成2 40 0 多人死亡，经济损失达数亿元。近年来，随着极端天气频发和人类工程活动增强，北京地区先后发生2 0 12 年“7 21”和2 0 18 年“716”泥石流灾害群发事件2 ,虽处置得当，避免了较大人员伤亡，但仍造成了较大财产损失和社会影响。长期以来，大量专家学者基于野外调查、监测预警等取得的一手资料，对北京市山区泥石流进行了大量研究。在2 0 世纪8 0 年代，洪惜英、李容全、张仲德等就开始研究北京地区泥石流的发育特征、分布规律等收稿日期：2 0 2 2-0 6-2 2基金项目：国家自然科学重点基金项目（42 130 7 2 0）；全国地质灾害风险区划技术方法研究（DD20221738）；北京市科技计划课题（Z191100001419015）；第二次青藏高原综合科学考察项目（2 0 19QZKK0902）作者简介：杨鑫，男，硕士研究生，主要从事地质灾害数值模拟及危险性研究工作。E-mail:通信作者：孟华君，男，副研究员，博士，主要从事地质灾害和地质环境等方面的研究工作。E-mail:文献标志码：AD01:10.16232/ki.1001-4179.2023.08.020基础性工作3-5】,旨在通过对泥石流灾害的科学认识，有效降低地质灾害造成的经济、财产损失，提出防治对策、保障社会安全。2 0 世纪90 年代，相关研究工作已经从泥石流事件的现场调查和资料统计分析逐渐转向成因机理分析及防治对策设计6-7 和泥石流危险性区划8-10 ,并尝试开展北京地区泥石流预警预报工作-15。2 1世纪后，随着仿真模拟技术、信息传输技术的快速发展，基于泥石流启动物理过程的泥石流灾害精细化风险管理成为可能，不断出现精度更高的预警模型。例如，白利平等16-17 根据力学平衡条件，尝试推出松散固体物质起动判别式和雨强表达式。赵忠海综合考虑地形地貌、地质条件、土壤类型以及降第8 期雨情况等因素，建立了泥石流是否发生的判断公式。王海芝18 通过整理1949年以来北京地区发生的泥石流与其发生时期的降水数据，建立了临界雨量的定量估算ID 模型等。这些工作为北京市泥石流灾害管控提供了理论指导和技术支撑。然而，随着京津冀协同一体化国家发展战略的实施，北京市对泥石流灾害管控提出了更高的要求，逐渐由灾害点管控向“隐患点+风险区”双管控转变。新的管控策略客观上要求对隐患点泥石流形成、运动及承灾机理进行深入研究19,而数值模拟方法是解决该问题的较好手段。通过数值模拟，可以再现泥石流运动演化过程，得到用于定量评价泥石流有关动力学、运动学要素和防治工程的参数2 0 目前，国内外开展泥石流数值模拟主要采用FLO-2D、Flu e n t、O p e n LI SEM 及PFC等成熟商业软件，研究结果显示数值模拟方法在模拟泥石流运动过程、冲出量、堆积等方面具有较好的效果2 1-2 4。鉴于此，本文采用数值模拟方法，选择北京市达摩沟作为研究区，利用更专业的泥石流模拟软件Massflow模拟达摩沟泥石流运动过程，通过研究泥石流隐患点在不同降雨工况下的泥石流运动及沉积过程，厘清达摩沟泥石流危险性，为该小流域泥石流隐患点风险管理提供依据。1研究区概况达摩沟位于北京市门头沟区清水镇东南方向3km处，是清水河支沟，位于东经11537 11539，北纬3952 3956（见图1）。该泥石流沟沟域面积为3.9 6 km,相对高差6 47 m，主沟沟向32 5,断面大多呈V形，长2 6 8 0 m，沟宽30 10 0 m，平均15m，沟床平均坡降2 30%。陡坎、卡口不多，河段宽窄均匀，主支沟交角大，沟道较顺直平整。沟床两岸坡面较陡，坡度30 6 0，平均40。植被覆盖率为7 0%。达摩沟内出露地层有侏罗系下统窑坡组（Jy）、龙门组（J,1）,侏罗系中统九龙山组（J2j）、馨髻山组（J,t)及第四系全新统（Q4）。岩性主要为砂岩、粉砂岩、安山岩、玄武岩、凝灰质砂岩、粉砂岩、页岩及煤系地层等,第四系全新统主要由砾石、岩屑、黏土质等残、坡积物组成(2 3。研究区地处北京西部山区，属中纬度大陆性季风气候，年平均气温10.2，极端最高气温37.6，极端最低气温西部-2 2.9。降水量自东向西逐渐减少,受中纬度大气环流的不稳定性和季风影响,降水量年际变化大。达摩沟流域历史上曾多次暴发泥石流，最近一次发生于19 50 年8 月4日，当日最大降雨量达19 0.1杨鑫，等：基于Massflow的达摩沟泥石流模拟与危险性分析2泥石流形成条件2.1物源条件（1）受达摩沟历史煤矿开采影响，主沟及支沟沟道堆积大量煤研石（见图2（a）（b）,其中仅少部分煤研石堆有浆砌石护坡（图2（a），绝大部分煤研石堆没有任何防护措施，稳定性较差，是达摩沟未来暴发泥石流的主要固体物源。（2）沟道两侧岸坡浅表层堆积大量坡积物。这些松散堆积物主要分布于坡面中下部浅表层，厚度不等，约30 8 0 cm（见图2（c）（d），植被覆盖较好。然而，根系延伸同时也使堆积更加松散，空隙发育，加之下伏基岩具有倾向主沟的优势结构面（图2（d），当降雨量较大时，雨水在缓慢入渗的同时极易沿大空隙灌入坡积物，抵达土岩接触面，在该面形成积水，从而引起崩坡积物失稳补给泥石流。图2（c）（d）可见，坡面下部已经缺失了一部分浅表堆积物，下伏基岩面直接出露。（3）在坡脚分布有零星崩滑堆积体（图2（e）和历史泥石流堆积体（图2（f)）。其中，泥石流堆积体可见明显的定向排列，大颗粒石块长轴方向基本一致。这些堆积体规模大小不等，固结程度较高，不易失稳，但在极端降雨工况下，受上游准泥石流体冲击侵蚀，会受到强烈侧蚀作用，可能作为泥石流补给。2.2水文条件达摩沟泥石流一般由强降雨激发。该地区降雨分布具有年际间不均衡、年内变化大、降雨集中、降雨强度大等特点。该沟一般年降水量为550 6 50 mm，平均为57 5mm，汛期平均雨量为40 0 50 0 mm；历史年147mm,最大雨强为56.5mm/h。持续的暴雨造成达摩庄后港沟和王家港沟泥石流爆发，造成8 4人死亡、2 4人重伤，基础设施损毁严重，给当地人民造成了大量的财产损失2 6 北京市Fig.1 Location of Damogou gully图1达摩沟地理位置148坡面石堆积人民长江3.1NMassflow模型控制方程沟道研石堆积质量运动的特征可以用NavierSt o k e s 方程来描述，该方程以微分守恒形式写成：+V p u=02023年(1)+V.p=V.T+pgat(2)式中：p为质量密度；g为重力加速度；u为质量速度矢(a)坡面研石堆积(a)Slopegangueaccumulation浅表层坡积物(c)浅表层坡积物A(c)Superficial slopedeposits A历史崩滑堆积物(e)历史崩塌堆积物(e)Historicalcollapseaccumulation图2 达摩沟泥石流物源Fig.2Material source of Damogou debris flow最大降水量9 32 mm，年最小降水量32 9mm。最大日降雨量为2 2 4.8 mm，最大小时降雨量达56.5mm。每年雨季降水量占全年降水量的7 0%左右，且多以大暴雨或特大暴雨形式出现。3研究方法数值模拟是解决单沟泥石流隐患点危险性评价的有效方法。本文采用专业泥石流模拟Massflow软件对研究区泥石流进行运动过程模拟。Massflow是中国科学院成都山地灾害与环境研究所欧阳朝军研究员研发的一款数值模拟软件。该软件基于深度积分连续介质力学理论，利用改进Mac-Cormack-TVD有限差分方法,兼顾复杂地形地貌,是具有二阶精度和自适应求解域特征的山地灾害动力学高效计算模拟软件。与国内外成熟软件相比，该软件不仅提供了丰富的物理模型以提高模拟准确性，还能够极大提升计算效率，缩短模拟时间。在甘肃岷县二马沟泥石流模拟2 7 、金沙江白格滑坡模拟2 8 以及泥石流防治工程防治效果模拟2 9 中，Massflow均取得了良好的效果。(b）沟道煤研石堆积(b)Trenchgangueaccumulation浅表层坡积物光滑下垫面(d)浅表层坡积物B(d)Superficial slopedeposits B历史泥石流堆积物(f)历史泥石流堆积物(f)Historical debrisflowaccumulations量,其分量位于坐标轴x,y,z三个轴上；T与重力和垂直方向对齐。可以通过在垂直深度上积分，将它们简化为二维浅水形式。如果质量在水平方向的运动远远超过在垂直方向的运动,这种简化是有效的。对质量守恒方程和动量守恒方程进行深度积分，再运用莱布尼兹法则和动力学边界条件化简方程，最终得到控制方程如下：ah+(hu)+a(hv)(3)+ata(hu?+a(hu)+at=gxh-kiapg.hapa(hv)a(h uv)+atx=g,h-hapg.hay式中：h,u和v分别为从基础底面z,到自由表面z的总质量高度、x方向上流体速度和y方向上流体速度；t表示时间；T为底部的剪应力;kap为侧向土压力系数；g，g，g.为各坐标轴上的重力分量30 。3.2参数选取3.2.1摩擦模型山地灾害数值模拟主要使用的摩擦模型有4种：Coulomb-viscous模型、Manning模型、Voellmy模型和Bingham模型,分别适用于不同的灾害类型。Coulomb-viscous模型主要适用于滑坡、岩崩等碎屑流灾害；Manning模型则适用于溃坝、洪水等水力学灾害；Voellmy模型适用于泥浆、泥石流灾害；Bingham模型广泛运用在尾矿、黏性泥石流等灾害模拟中31。本文根据达摩沟历史泥石流特征和现有物源特性，选用Voellmy模型对达摩沟泥石流进行模拟计算。Massflow软件中Voellmy模型表达式如下：T=ou+Pgu$式中T为底部的剪应力,Pa;为正应力,Pa;u为摩擦=0 x+2at(T）xPa(hv+kapg,h/2)ay(T.）Pay+a(huv)at(4)(5)(6)第8 期系数;p为泥石流密度,kg/m;为流系数。3.2.2参数选取使用Massflow软件模拟泥石流运动过程主要关注3个参数：摩擦系数、流系数和泥石流流量。摩擦系数和端流系数通常采用2 个步骤来确定：参考历史泥石流事件确定参数取值范围。以往研究成果中,摩擦系数的取值范围为0.1 0.3,端流系数的取值范围为10 0 30 0。基于第一步确定的摩擦系数和端流系数来反演确定具体参数值。本文根据现场古泥石流堆积情况，进行了一系列数值模拟反演计算后，通过整体比较，发现当摩擦系数=0.3,瑞流系数=32 0 时的模拟结果与现场古泥石流堆积情况最为吻合。另外，Massflow软件不能直接输入降雨量，因此，需要将降雨量转换为流量进行计算。本文根据现场实地调查资料，并参考历史泥石流资料，计算得到达摩沟泥石流密度为1540 kg/m。在此基础上，根据预设的4种不同降雨工况（10,2 0,50 10 0 a一遇），采用雨洪法计算不同降雨工况下泥石流流量。根据DZ/T02392 0 0 4泥石流灾害防治工程设计规范中的雨洪计算法，可以计算得到不同降雨频率下的泥石流流量：Q。=(1+)Q,D式中：Q，为清水峰值流量；为泥沙修正系数;D。为堵塞系数。雨洪法计算中，泥石流清水流量的计算公式为Q,=0.278 FT式中：0.2 7 8 为单位换算系数；F为流域面积,h,p为净雨雨量；T为流域汇流时间。达摩沟流域清水流量及泥石流流量计算结果如表1 所列。表1不同降雨条件下流量计算值Tab.1 Calculated flow under different rainfall conditions降雨频率P/泥石流峰值流量Q。/%(m:s-)1017.88522.39227.68134.374模拟结果与分析4.1计算模拟结果通过对10%，5%，2%，1%降雨频率下达摩沟爆发泥石流的过程进行模拟，可以得到不同降雨条件下达摩沟泥石流的运动过程和不同时刻的泥深以及流速杨鑫，等：基于Massflow的达摩沟泥石流模拟与危险性分析h清水流量Q/149分布，并在此基础上开展达摩沟泥石流隐患点在不同降雨频率下的危险性分析评价。部分模拟结果如图3,4所示。流速/(ms)0.00-1.691.69-3.183.18-7.497.49-9.699.69-13.32流速/（m：s)0.00-2.492.49-4.764.76-8.0618.06-12.96112.96-18.96图3不同降雨频率条件下泥石流流速分布Fig.3 Distribution of debris flow velocity under different(7)(8)(m3:sl)41.3051.7263.9379.40流速/（ms)0.00-2.412.41-4.594.59-7.927.92-13.4613.96-18.32(a)P=10%流速/(ms)0.00-2.552.55-4.774.77-8.108.10-13.1013.10-18.32(c)P=2%rainfall frequencies泥深/m10.00-0.530.53-0.9010.90-1.301.30-1.8511.85-3.10泥深/m10.00-0.740.74-1.321.32-1.961.96-2.632.63-3.74(c)P=2%图4不同降雨频率条件下泥石流泥深分布Fig.4Depth distribution of debris flow under differentrainfall frequencies(b)P=5%(d)P-1%泥深/m10.00-0.770.77-1.321.32-1.981.98-2.6112.61-3.86(a)P=10%泥深/m10.00-0.7310.73-1.291.29-1.9211.92-2.562.56-3.65(d)P=1%(b)P=5%1504.2模拟结果分析4.2.1流速变化根据计算结果（见表2），不同降雨工况下达摩沟泥石流流速变化具有以下3个特点。表2 不同降雨条件下流速占比分布Tab.2Flow velocity distribution under differentrainfall conditions频率04m/s48m/s812m/s1216m/s1620m/s1092.5581.1278.5178.4（1）总体流速较为稳定。在10,2 0,50,10 0 a一遇降雨概率下，整个主沟道流速总体普遍维持在0 4m/s的水平；只有下游的阶地以及陡峭区域流速急剧增大，局部平均流速达到8 m/s以上（见图3）。（2）峰值流速随降雨强度增加，呈先增大后逐渐趋于稳定的趋势,见图5。其中,10 a一遇和2 0 a一遇降雨频率下的峰值流速差异较大，而2 0 a,50a,100a一遇降雨的峰值流速差异较小，总体维持在18 2 0m/s 的水平。2018(_5.u)/率16141210图5不同降雨概率工况下峰值流速趋势Fig.5 Peak flow under different rainfall probability conditions（3）随着降雨强度增大，泥石流高速流动区的距离明显变长。对比50 a一遇和10 0 a一遇降雨工况下高速流动区的沟道长度，发现10 0 a一遇降雨条件下的峰值流速略小于50 a一遇降雨条件下的峰值流速，但其高速流动区更长。以上3点表明：达摩沟泥石流峰值流速受限于降雨强度和沟道地质条件两个方面。当降雨强度小于20a一遇时,峰值流速会随降雨强度增加而增大；当降雨强度达到2 0 a一遇时，降雨强度继续增大，并不会继续提高峰值流速，而会延长高速流动的距离。降雨强度增大导致泥石流流速和高速流动距离增加，本质上是提高了泥石流体在流动过程中的冲击速度和冲击距离，即提高了侵蚀效率，大大提高了泥石流体的破坏能力。这既有利于沿程物源进一步补给泥石流体，又进一步增加了泥石流体的重度，提高了破坏力。人民长江4.2.2运动过程泥深变化泥石流堆积深度是评估泥石流危险性的重要指标。图4，6 是不同降雨工况下泥石流堆积情况统计趋势。总体来看，达摩沟泥石流堆积特征主要可以概括如下：泥石流最大泥深随降雨强度增加呈先增大后减小的趋势；达摩沟在不同降雨条件下均不会形成%明显的堆积扇；泥石流在沟道狭窄处泥深将会增大，偏向流通区上游的部分堆积深度较小，而流通区中5.51.914.82.916.33.616.33.720重现期/a2023年0.100.000.700.500.121.480.121.4850100下游和堆积区等泥石流前缘部分堆积深度较大。4.03.8F/3.6F3.4F3.23.010图6不同降雨概率工况下峰值泥深趋势Fig.6Peak mud depth under different rainfallprobability conditions造成以上泥石流堆积特征的原因主要有：在降雨强度较低的情况下，水流对堆积物的侵蚀能力和搬运能力均较弱，泥石流流量较小；随着雨量的增加，产汇流增加，水流的侵蚀搬运能力逐渐增强，泥石流流量变大，导致泥石流的堆积厚度增大；当降雨强度增加到一定程度后,泥石流流量增加、流速变大，其裹挟能力更强，更容易将固体物质搬运至更远的地方，因此高强度降雨条件下沟道内的最大堆积厚度反而有所减小。达摩沟狭窄的沟道不具备形成堆积扇的地形条件。主沟下游地势相对平坦，有利于泥石流运动减速开始堆积，但历史泥石流在流通区下游和近沟口处有堆积，导致目前流通区末端和近沟口处地势仍然相对较高，不利于泥石流淤积蠕动形成扇体，而有利于未来泥石流活动进一步向清水河流动；此外，在流通区下游以及堆积区南侧有部分居民区并修建有道路，导致达摩沟主要为带状堆积。达摩沟与清水河交接处沟道比降仍较大，不利于泥石流减速沉积，更可能出现的情况是泥石流体高速冲人清水河，因此，达摩沟沟口没有明显的堆积扇。由于沟道较窄处泥石流的运动受到阻碍，泥石流流速变缓，导致其更易在此处堆积下来。在泥石流堆积过程中，随着泥石流不断向周围扩散堆积，其流速、动能不断降低，泥石流边运动边堆积，最终整体停止运动，所以堆积区中心厚度最大，而发散区20重现期/a50100第8 期域堆积厚度逐渐降低。4.2.3危险性评价泥石流的强度可以依据泥石流的流速与泥深的强度进行划分。本文参照常鸣32 提出的泥石流强度影响以及危险性分区标准对模拟结果进行危险性评价。泥石流强度影响划分标准如表3所列。表3泥石流强度影响划分Tab.3Strength division of the debris flow泥石流强度泥深H/m高H2.5中0.5H2.5低0H0.5注：OR表示或条件,2 个条件（H与VH）满足其中一个即可;AND表示与条件，必须同时满足2 个条件。根据表3,利用Arcgis对泥石流流速与泥深进行栅格叠加计算,将计算结果分类为高危险、中危险、低危险，最终得到10 0 a一遇降雨工况下的达摩沟泥石流危险性分区如图7 所示。危险性低中高图7 达摩沟泥石流危险性分区Fig.7Hazard zones of Damogou debris flow在图7 中，高危险性地区占16.54%，中危险性地区占49.2 7%，低危险性地区占34.19%。低中危险地区主要分布在主沟边缘地带，主沟绝大部分区域处于中高危险区。考虑到达摩村位于主沟沟口处，涉及到大量人、财、物的潜在损失，说明达摩沟是一条风险性较大的泥石流沟，因此需要加强对达摩沟的预警监测。5结论与建议（1）达摩沟小流域在强降雨条件下有发生泥石流的可能，模拟结果显示：泥石流流速主要分布在0 4m/s,随着降雨强度的增大，泥石流平均流速逐渐增大，但泥石流峰值流量随着降雨强度的增大呈现先增大后减小的趋势。达摩沟泥石流堆积泥深主要分布在03m之间，最大泥深可达3.8 5m，平均泥深随着降杨鑫，等：基于Massflow的达摩沟泥石流模拟与危险性分析计算关系泥深与最大流速（V)乘积/(ms-)ORANDAND151雨强度的增大而增大，峰值泥深则呈现先增大后减小的趋势。（2）依据泥石流流速与泥深，利用Arcgis对达摩沟泥石流危险性进行分区，结果显示：高危险区占16.54%，中危险区占49.2 7%，低危险地区占34.19%。然而，中高危险区仍有居民居住,因此风险较高,建议加强监测。（3）泥石流流速及泥深与降雨强度的统计关系表明：泥石流运动不是简单受降雨强度影响，还受到地形VH2.5条件、植被覆盖等多因素约束，但是不同工况下各因素0.5VH2.50VH0.50530106015902120im贡献有差异。因此，在对泥石流进行预警监测时,需针对不同小流域进行具体分析，建立不同的监测方案。参考文献：1赵忠海.北京泥石流灾害预警研究现状与技术方法探讨J.城市地质,2 0 15,10（增1)：7 8-8 4.2 北京市地质研究所，中国地质科学院地质力学所，北京林业大学。北京市沟道侵蚀泥石流多元预警阅值模型研究及示范工程R.北京：北京市地质研究所，中国地质科学院地质力学所，北京林业大学,2 0 2 1.3洪惜英，王礼先，朱靖才，等.北京市密云县汗峪沟泥石流J.北京林学院学报，198 3（4）：38-45.4李容全，徐振涛.北京山区泥石流分布特点J.山地研究，198 3(3):42-48.5张仲德.北京山区泥石流分布和成因分析J.北京师院学报（自然科学版）,198 4（2）：7 7-8 9.6魏永明.19 9 1年北京北部山区泥石流成因及防灾对策J.山地研究,19 9 3(3):18 7-19 2.7张世清，毕小刚，刘振国.北京山区泥石流防治对策探讨J.水土保持通报,1992(3）：46-51.83张春山.北京北山地区泥石流灾害危险性评价J北京地质，1996(2):11-20.9张春山.北京地区泥石流灾害危险性评价J.地质灾害与环境保护,19 9 5(3):3340.10 赵波，董桂芝.北京西山清水河流域泥石流及危害预测J.北京地质,19 9 6(3):1-14.11 谢又予，之久.北京山区泥石流的防治与预报.贵州科学，1992(3):132.12吴正华，储锁龙.北京泥石流暴雨基本特征J.大气科学，1992(4):476 481.13韦京莲，董桂芝，赵波.北京山区近代泥石流活动规律及暴发周期分析J.中国地质灾害与防治学报，19 9 4（4）：48-53.14韦京莲，赵波，董桂芝.北京山区泥石流降雨特征分析及降雨预报初探J.中国地质灾害与防治学报，19 9 4（增1）：45-5215吴正华.北京泥石流灾害及其降水触发条件J.水土保持研究，2001(1):67-72.16白利平.北京山区泥石流灾害临界雨量研究J.地质灾害与环境保护,2 0 0 6(4)：10 1-10 4.17白利平，孙佳丽，南赞.北京地区泥石流灾害临界雨量闵值分析J.地质通报,2 0 0 8(5）：6 7 4-6 8 0.18王海芝.北京地区暴雨泥石流预警阅值研究J.第四纪研究，2020,40(5):1371-1380.15219 黄晓虎.泥石流降雨启动机制及早期预警模型的研究：以达摩沟为例D.长春：吉林大学,2 0 16.20 涂剑，马超，杨海龙.北京山区暴雨泥石流激发雨量条件J.中国水土保持科学,2 0 17,15（5）：10 3-10.21方群生，陈志和，唐川，等.基于FLO-2D数值模拟震区急陡型泥石流冲出量：以瓦窑沟为例J.中山大学学报（自然科学版）,2 0 2 0,59(3）：2 3-31.22 赵峰，刘铁骥，陈亮，等.基于Fluent强震区王家沟泥石流运动特征研究J灾害学,2 0 2 1,36(2)2 0 8-12.23 陈明，唐川，王飞龙，等.基于OpenLISEM模型的泥石流启动冲出预测J.泥沙研究,2 0 2 0,45（4）：59-6 5.24 胡明鉴，汪稔，陈中学，等.泥石流启动过程PFC数值模拟J.岩力学,2 0 10,31（增1）：394-397,434.25 慎乃齐，韦京莲，赵波，等.京西达摩沟泥石流形成条件与防治J.探矿工程（岩土钻掘工程）,2 0 0 3（增1)：96-98.26 齐干，张长敏.达摩沟泥石流形成的物质条件分析及防治对策J.水文地质工程地质,2 0 11,38（5）：10 2-10 9.人民长江27乔渊，刘铁骥，陈亮，等.基于Massflow模型的甘肃省岷县二马沟泥石流危险性评价J.水利水电技术,2 0 2 0,51（4)：18 4-192.28曹水和，吴新明，钟东.基于Massflow的金沙江滑坡数值模拟及堵江风险预测J.地质灾害与环境保护,2 0 2 1,32（4)：3-7.29OUYANG C,HE S,XU Q,et al.A MacCormack-TVD finite differ-ence method to simulate the mass flow in mountainous terrain withvariable computational domain J.Computers&Geosciences,2013,52:1-10.30HUNGR O,MCDOUGALL S.Two numerical models for landslide dy-namic analysis J.Computers&Geosciences,2009,35(5):978-92.31周凯琦.基于单层单相流和多层多相流的灾害动力学数值模拟D.武汉：华中科技大学,2 0 17.32 常鸣.基于遥感及数值模拟的强震区泥石流定量风险评价研究D.成都：成都理工大学,2 0 14.（编辑：高小雲）2023年Numerical simulation and risk analysis of mudslide in Dharma Gullysub-basin based on MassflowYANG Xin,MENG Huajun-,LU Lu,WANG Yu,WU Jihuan(1.Institute of Geomechanics,Chinese Academy of Geological Sciences,Bejing 100081,China;2.School of Civil and ResourcesEngineering,University of Science and Technology Bejing,Bejing 100083,China;3.Key Laboratory of Neotectonic Movementand Geohazard,Beijing 100081,China;4.Observation and Research Station of Geological Disaster in Baoji,Shaanxi Province,Ministry of Natural Resources,Beijing 100081,China;5.Beijing Geological Research Institute,Beijing 100083,China)Abstract:Mudslide disaster events had occurred in the Dharma Gully sub-basin area historically and caused a large casualtiesand property losses.At present,a large amount of coal gangue was piled up in the main ditch of the sub-basin,and more land-slide deposits were accumulated on both slopes sides of the main ditch.In order to analyze the possibility and risk of mudslides re-occurring,this paper analyzed the development of mudslides in Dharma Gully sub-basin area in detail through UAV photogram-metry and field investigation.Based on the improved MacCormack-TVD finite difference method,the Massflow software was usedto 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 debrisflow 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 muddepth could reach 3.85 m;U n d e r t h e r a i n f a ll o f 10 0 -y e a r,t h e h i g h -r i s k a r e a s i n D a m o g o u a c c o u n t e d f o r 16.54%,t h e m e d i-um-risk areas accounted for 49.27%,and the low-risk areas accounted for 34.19%;The simulation results reproduced thetemporal and spatial development characteristics of debris flow well,showed the influence of rainfall intensity on the movement ofdebris 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