1、第47 卷第2 期2023年6 月doi:10.3969/j.issn.1674-3636.2023.02.012雅鲁藏布江色东普沟冰崩一堵江一溃决灾害数值模拟田杨杨，姜亮，郭江（1.西藏地勘局区域地质调查大队，西藏拉萨8 51414；2.四川省地质环境调查研究中心，四川成都6 10 0 36)摘要：为了揭示雅鲁藏布江色东普沟2 0 18 年10 月17 日冰崩一堵江一溃决灾害链的动力演化过程，基于Massflow数值模拟仿真平台，使用Fortran编程语言，根据研究区域地质条件特征对程序进行二次开发以优化Voellmy模型，模拟冰崩一泥石流动力过程；将模拟泥石流得到的堰塞坝体嵌入地形中，运用
2、ArcGIS计算堰塞湖范围及体积，通过Manning模型模拟堰塞湖溃决洪水动力过程。采用分段模拟法再现冰崩一泥石流一堵江一堰塞湖一溃坝的完整动力过程，对泥石流运动过程中的流速、流深，坝体高度，溃决洪水的流深、流速等参数进行定量化研究，为色东普流域的防灾减灾工作提供有效支撑。关键词：冰崩；泥石流；堰塞湖；灾害链;Massflow数值模拟；色东普沟;雅鲁藏布江中图分类号：P642.230 引 言堵江事件是指高海拔斜坡上的冰川在重力作用下断裂解体后迅速滑塌、铲刮沟道并裹挟碎屑物堵塞河道形成堰塞坝，在坝体溃决后又引发洪水的灾害,属于典型的冰崩一泥石流一堵江一溃决链式灾害。为揭示该链式灾害的动力过程，对
3、灾害链进行模拟反演。基于动力过程的数值模型主要有连续介质力学模型、离散介质模型及混合介质模型（王钟文，2 0 2 0）。Savage等（19 8 9）最早提出了深度积分连续介质力学的方法概念，并在干性沙土的流动试验中得到了验证。连续介质运动方向特征尺度远大于垂直方向特征尺度，使得复杂的三维地表过程可通过二维计算求解。该方法提出以来,崩滑碎屑流、洪水、泥石流等地表过程动力学数值模拟研究均取得了突破性进展。深度积分连续介质力学方法在泥石流等山地灾害动力学过程分析中具有独特优势：可输出多种泥石流特征参数，不仅能为定量化评价提供更多的参考依据,而且可以与已有的野外监测数据、大规模试验数据进行对比分析,
4、提高模型的准确性(Iverson收稿日期：2 0 2 2-0 9-30；修回日期：2 0 2 2-11-0 9；编辑：陈露，朱明君基金项目：西藏自治区自然资源项目“监测示范工程（SCXX-XZCG-20098）作者简介：田杨杨（19 8 1一），男，工程师，主要从事地质矿产、地质灾害动力演化过程研究，E-mail：35155447 q q.c o m通信作者：姜亮（19 9 4），男，助理工程师，硕士，主要从事地质灾害动力演化特征研究，E-mail：1130 18 32 55 q q.c o m地质学刊Journal of Geology文献标识码：A文章编号：16 7 4-36 36(2 0
5、 2 3)0 2-0 19 6-0 7etal.，2 0 10,2 0 14);将上下基底位置演变的影响直接计人控制方程，减少了单独处理边界的麻烦(Biscarini et al.,2010;George et al.,2014);极大减少了控制方程中未知量的个数,结合适当的数值求解方法和网格划分方法,能够在大区域高地形分辨率的前提下快速求解，满足地质灾害应急抢险救灾的要求(Berger et al.,2011;George,2011)。目前，对泥石流灾害强度及范围的获取主要有经验公式法、统计类比法和基于动力过程的数值分析法。经验公式法及统计类比法虽然能计算泥石流的流速及运动距离（余斌，2 0
6、 0 8），但误差通常较大。针对单次泥石流的运动距离，经验公式法的误差可能 50%(Guo et al.,2014)。O u y a n g 等(2 0 15)从广义的深度积分方法出发，运用MacCormack-TVD有限差分格式求解开发的Massflow软件,计算效率高，支持二次开发，被应用于实际地质灾害的数值分析中。以雅鲁藏布江大峡谷段的色东普流域为研究区,采用连续介质法，基于Massflow软件平台模拟2018年10 月17 日冰崩导致泥石流一堵江一溃决洪水灾害的动力演化过程，为揭示灾害动力演化特征、提升青藏高原地区防灾减灾能力提供理论依据。Vol.47 No.2Jun.2023第47
7、卷第2 期1区域环境背景1.1流域特征色东普流域位于西藏林芝东侧约55km，加拉白垒峰西坡,雅鲁藏布江左岸。该流域面积约6 6.8 9km,主沟长约7.4km,冰川发育，岩土体风化严重，侵蚀剥蚀作用强烈。流域内高程最高点为加拉白垒峰主峰7 2 9 4m，最低点为色东普沟口2 7 46 m,高差4548m，平均高程4540 m，平均坡度约34.8 9。流域内总体坡度较陡,坡度 30 的区域占整个流域的50%以上,在流域中部有一平均坡度 15的平坦区域,为冰川和冰物主要覆盖区域，由于冰川侵蚀和山体岩崩，该区域覆盖大量冰物等碎屑物质,为泥石流提供了丰富的物质条件。1.2地层岩性研究区出露地层主要有前
8、寒武、泥盆、石炭、二叠、三叠、侏罗、白垩系和第四系。除侏罗、第四系外,其他地层的岩石均有不同程度的变质。岩性主要为混合片麻岩、斜长片麻岩、石英质糜棱岩等，岩石总体较坚硬,抗风化能力较强。但岩体中片麻理、片理发育，同时岩体受北东、北西、北北西、近东西向等多组断裂的切割破坏,岩体中节理发育,具有发育崩塌、滑坡的构造条件，也为冰崩碎屑流的产生提供了丰富的松散固体物质。此外，第四系发育于各级河流、谷地、山麓及其边缘,成因有崩积与滑坡堆积，这些成因复杂的第四系松散物质在受到冰崩碎屑流的铲刮后更容易导致大规模堵江事件（黄文星等，2013)。1.3地质构造色东普沟分水岭地带常年冰雪覆盖，多条分支冰川活动造成
9、沟谷深切，谷坡陡峻，冰物丰富。雅鲁藏布江大拐弯地处印度板块与欧亚板块碰撞形成的东喜马拉雅构造结,是挤压、旋转走滑和隆升的强烈地壳活动区。区内主要发育北东向墨脱断裂、米林断裂、里龙断裂,以及北西向嘉黎断裂、察隅断裂、阿帕隆断裂、雅鲁藏布江断裂等,断裂相互切割贯通导致山体表面碎裂化。崩滑-泥石流多发生在月均温度高、日均温差大的时间段。色东普沟域位于阳坡，日夜温差大,冰雪融水较多,冰川破裂退缩，物理风化加剧,坡面水动力作用强烈，冻融作用为冰物长期饱水液化、蠕田杨杨等：雅鲁藏布江色东普沟冰崩一堵江一溃决灾害数值模拟更易发生渗流滑动作用。1.4气象水文色东普流域属高原温带半湿润性季风气候，雨量充沛，日照
10、长、辐射高，昼夜温差大，无霜期短（每年约17 0 d），降水量垂直变化大，干湿季节分明。区内气候较湿润，无明显夏季，多年平均气温为8.2，全年平均日照时数为140 0 2 0 0 0 h；降水相对较多,多年平均降水量为7 13.1mm,6一9 月降水量占全年总降水量的8 0%以上。研究区内河流均属雅鲁藏布江水系，地表水系发达。雅鲁藏布江两侧次级水系发育，由南向北或自北向南汇人雅鲁藏布江,平面上呈树枝状展布。研究区内的雅鲁藏布江段水面落差和平均坡降大，多年平均径流深度为1 5003 000 mm。2冰崩泥石流事件2.1典型事件色东普沟在19 6 9 2 0 18 年至少发生了8 起冰崩碎屑流事件
11、（杨超，2 0 19），其中2 0 18 年发生了2起典型冰崩事件：2 0 18 年10 月17 日，发生大规模冰崩碎屑流,碎屑物质在沟口堆积形成堰塞坝，上游回水15km，水位上涨40 m，库容约1.5亿m；2 0 18 年10 月2 9 日，坝体表面堆积大量新的碎屑物质,堵塞了色东普沟口的狭窄流水通道,形成堰塞湖，堵江坝体宽约310 6 2 0 m，顺河长度约2.3km，高9 0 m，估算堰塞体体积约6 6 0 0 万m。2.2原因分析2.2.1地形条件色东普流域的地形呈典型的阶梯状,利于泥石流的形成。上部陡峭区域易引发岩体和冰雪的崩塌,使得碎屑物不断崩落;中部平缓区有利于累积崩塌物质，为泥
12、石流提供了大量的物源；下部较陡区域为泥石流的高速下滑提供了地形条件。该流域内典型的陡一缓一较陡的阶梯状地形特征,是泥石流频繁发生的基础（王磊等,2 0 14）。2.2.2沟口碎屑物堵塞河道状态2 0 17 年12 月2 1日发生的泥石流规模较大，在色东普沟口形成了扇形堆积体，一直未被冲毁,堵塞了雅鲁藏布江河道宽197动、流动提供了条件，形成高含水的滑坡-泥石流或冻融泥石流。近年来,色东普沟降水量增加,冰川裂隙渗流增大了冰川底部静水压力，降低了底部摩擦力，导致冰川易发生垮塌；冰崩冲击使得冰物底部198度的大部分区域,且之后的泥石流事件十分频繁,因此色东普沟口雅鲁藏布江河道大部分一直处于堵塞状态,
13、即使发生较小规模的泥石流也可造成雅鲁藏布江堵塞回水。2.2.3冰川类型色东普沟冰川为海洋型冰川,与大陆型冰川相比更易受到全球气温升高的影响而加速消融，运动速度更快，运动距离更长。2.2.4冰裂隙色东普流域内的冰川表面裂隙较多,冰川较为破碎，在冰川表面覆盖有较厚的表物,冰舌表面的色调与土体极为相似,流域内丰富的冰物为冰崩泥石流的发生提供了物质基础。同时,大型裂隙的存在使降水更易下渗,降低了冰川底部的摩擦力，也更易导致冰川在重力作用下沿裂隙发生崩塌。2.2.5林芝地震2 0 17 年11月18 日林芝发生地震，震中距色东普流域约2.7 km，震级Mw6.4,震源深度8 km。地震对该流域造成严重扰
14、动,破坏了冰川及沟谷内堆积物和冰物的完整性,使末端有临空面及后缘有宽大裂隙的冰川变得更不稳定，沟谷内堆积物和冰物更松散，为后来连续发生冰崩灾害提供了物质基础（普布次仁等,2 0 2 2）。2.2.6降水色东普流域2 0 17 年降水量偏高,导致色东普沟冰川累积区增大,加之冰川表物大且厚，降水后重力进一步增加,致使冰川下滑力增大，容易发生垮塌（童立强等,2 0 18）。3模型构建色东普沟冰崩一堵江一溃决灾害链本构模型基于Voellmy和Manning模型,其中Voellmy模型模拟冰崩形成泥石流导致堰塞坝的动力过程，Manning模型模拟堰塞湖溃决洪水动力演化过程。灾害链各阶段的衔接过程：将模拟
15、泥石流得到的堰塞坝体模型嵌人地形模型中,采用ArcGIS计算堰塞湖范围及体积,将其作为起始条件模拟堰塞湖的溃决动力过程。采用12.5m12.5m的DEM数据,通过ArcGIS将DEM栅格数据导出为ASCII文件，再将ASCII文件和模型参数导人Massflow中构建三维地形模型。Massflow以深度积分连续介质力学法为理论基础，采用MacCormack-TVD有限差分法求解,实现地表灾害动力过程的数值分析。地质学刊3.1泥石流动力过程模型构建将泥石流视为不可压缩流体，密度恒定，根据质量守恒定律及动量守恒定律，在笛卡尔直角坐标系中建立示意图,使之与平均坡度角平行。分别围绕X、Y轴旋转角度、0,
16、；流动层上、下部临界表面分别用z,z,表示,流动层厚度h=z,-zb;u、w 分别为流体速度在X、Y、Z轴上的分量（图1）。ZZ二Y0图1泥石流流动示意图Fig.1 Schematic diagram of debris flow流动层需满足质量守恒方程（1)以及X、Y、Z方向上的动量守恒方程（2）、（3）、（4）：p,(pu),(pv)(pw),atdxdy(pu),(pu)(pw),(puw)atxaT+za(pv)a(puw),(pu),(pvw)atxaydTxOTTXy十y(pw),a(puw)(pvw)(pw)otxT2ayaz式中,p为流体的密度，t为时间，T,为不同方向的应力分
17、量（如T为X方向的切应力，T,为XY方向的切应力），gvgyvg.分别为重力加速度g在X、Y、Z轴上的分量。矩阵形式的表达为：2023年Whdydy流动层侵蚀沉积层X=0(1)=pg+zz=pg,+=pg,+z(2)(3)(4)第47 卷第2 期cos 0,0-sin,y8y0(8.）sin 0,0cos 010cos 0.0sin 0LO-sin(cos x(g根据动力学边界条件（Gray，2 0 0 1），流体上部和下部临界面需分别满足方程（6）、（7）：+Votxazb+b+Vbtxay3.2堰塞湖溃坝模型构建针对最经典的流体动力学三维Navier-Stokes方程进行深度积分，推导出溃
18、决洪水必须满足的质量和动量方程：=0(pu)(pu)a(puxU)a(puv,),ypgx+x(pu,)a(puu,)a(pu,aT.T,yypgx+xa(pu,)(pv.v.)a(pv,.)OTaTaT十Pgx+xayz基底阻力模型采用Manning阻力模型：pgnT3.3参数选取参数选取见表1。泥石流较大的内摩擦角会产生更大、更快的能量耗散，使运动过程整体滞后；较高的流系数增加了流体层之间的动量交换强度，增强了与周围介质混合的能力，对周围介质的卷吸作用增大；沟道摩擦系数与灾害体流速、移动距离，以及堆积区体积、面积关系很密切,在很大程度上影响泥石流和溃决洪水的致灾范围和规模。田杨杨等：雅鲁藏
19、布江色东普沟冰崩一堵江一溃决灾害数值模拟100W,=0w,=0ay十ay十z199表1模型和计算参数Table 1IModel and calculation parameters容重内摩擦角角摩擦湍流灾害类型模型/(kN/m)/（）系数n系数m泥石流Voellmy(5)堰塞湖溃决Manning注：数据据李昆仲等，2 0 2 1；空白表示无此项4模拟结果(6)根据数值分析结果，泥石流从启动点至雅鲁藏布江以沟道流动为主,沟道狭窄,泥石流并无明显的(7)扩散流动堆积现象。泥石流的流深介于7 2 10 6 m之间（图2）,运动速度介于2 0 2 5m/s之间,运动速度快,对沟道的基底侵蚀作用较大，主
20、要原因是此段地势陡立,大比例的势能转化为动能（图3），泥石流在演进过程中对基底快速铲刮、侵蚀，并夹带沟床松散沉积物和松散坡积物,规模逐渐增加,从而增大了(8)下游的致灾规模。模拟结果与现场调查结果基本一致。泥石流在流通区的流动深度一般为8 0 9 0 m，宽度一般为50 10 0 m,流速为15 18 m/s。泥石流在沟道下游流速逐渐降低，堆积厚度逐(9)渐增加,靠近沟口的运动速度逐渐降低,因地势相对(pu,v,)平坦且沟道开阔,在沟口呈扩散堆积,最后泥石流堆积扇堵塞雅鲁藏布江,形成泥石流堰塞坝，坝体宽度600m,高度一般介于8 0 9 0 m之间,最高处达9 2(10)m,与现场调查结果基本
21、吻合。a(pu图4为泥石流堵塞河道导致形成堰塞湖的过一.2110程,泥石流堰塞坝形成2 7.5h后，湖水在泥石流堰塞坝右岸自然漫顶过流,水流冲蚀泥石流堰塞坝,坝(11)体迅速溃决（钟华邦，2 0 10）。其中，溃口最大流速为2 1 m/s,溃决洪水向下游快速演进,水位逐渐恢复到灾前水位，图5、图6 为堰塞湖溃决过程的流(12)深、流速演化过程。5 结 论通过资料收集、遥感解译、野外调查等方法对色东普冰崩灾害进行研究,并结合数值模拟方法对灾害链进行反演，得出以下结论。(1)泥石流在流通区的流深介于7 2 10 6 m之间,但在该区域的运动速度及对沟道的基底侵蚀作240.0512000.04200
22、地质学刊2023年N流深50s72m6km0(a)50 s流深(a)Flow depth at the 50th second;(b)Flow depth at the 250h second;(c)Flow depth at the 500t second流深250s106m6km0(b)250 s流深图2 泥石流流深Fig.2Flow depth of the debris flow流深500s89m6km0(c)500s流深流速50523m/sokm0(a)50 s流速(a)Flow velocity atthe50t second;(b)Flowvelocity atthe250t s
23、econd;（c)Fl o w v e l o c i t y a t t h e 50 0 s e c o n d用较大；泥石流的运动速度介于2 0 2 5m/s之间，主要原因是此段地势陡立，大比例的势能转化为动能。模拟结果与现场调查结果基本吻合。(2)色东普沟冰崩形成泥石流,堵江形成堰塞坝，模拟结果显示坝体宽度约为6 0 0 m，顺河长度约为2 km,高度最高为9 2 m,与实地状况基本一致。（3）通过对色东普沟2 0 18 年10 月17 日冰崩灾害链进行模拟反演,验证了模型的适用性,为未来冰崩灾害链的防灾减灾工作提供了依据。流速250s24m/sokm0(b)250 s流速图3泥石流流
24、速Fig.3 Flow velocity of the debris flow参考文献黄文星，王国灿，王岸，等，2 0 13.雅鲁藏布江大拐弯地区河流形态特征及其意义J.地质通报，32(1：130-140.李昆仲，张明哲，邢爱国,2 0 2 1.雅鲁藏布江色东普沟崩滑-碎屑流过程模拟及运动特征分析J.中国地质灾害与防治学报,32(1):18-2 7.普布次仁,巴桑多吉,刘松超,2 0 2 2.西藏仁布县巴曲沟泥石流发育特征及防治建议J.四川地质学报，42（1)：100-103.流速500s一6 m/s6km0(c)500 s流速第47 卷第2 期田杨杨等：雅鲁藏布江色东普沟冰崩一堵江一溃决灾害
25、数值模拟201流深33h32m6km0(a)33h后堰塞湖深度及范围Fig.4,Formation process of the barrier lake(a)Depth and extent of the barrier lake after 33 hours;(b)Depth and extent of the barrier lake after 57 hours;c)Depth and extent of the barrier lake after 85 hours流深57h62m6km0(b)57h后堰塞湖深度及范围图4堰塞湖形成过程N流深85h92mokm0(c)85h后堰塞湖深
26、度及范围V流深150s76mkm(a)150s溃决洪水流深(a)Flow depth of the outburst flod at the i5oh second;(b)Flow depth of the outburst flood at the 450h second;(c)Flow depth of the outburst flood at the 90oth second流深450s62m0(b)450s溃决洪水流深图5堰塞湖溃决过程流深Fig.5Flow depth in the barrier lake outburst processN流深900 s55mKm0(c)900s
27、溃决洪水流深流速150s21m/s6km0(a)150s溃决洪水流速Fig,6Flow velocity in the barrier lake outburst process(a)Flow velocity of the outburst flood at the 150i second;(b)Flow velocity of the outburst flood at the 450 second;(c)Flow velocity of the outburst flood at the 90o second流速450s18m/s0(b)450s溃决洪水流速图6 堰塞湖溃决过程流速流速9
28、00s17m/sKm0(c)900s溃决洪水流速202童立强，涂杰楠，裴丽鑫，等，2 0 18.雅鲁藏布江加拉白垒峰色东普流域频繁发生碎屑流事件初步探讨J.工程地质学报,2 6(6):1552-156 1.王磊，冯毅，2 0 14.四川茂县宗渠沟低危险堰塞湖溃决泥石流发育特征及工程对策J.四川地质学报,34（增刊2)：85-87.王钟文，2 0 2 0.基于动力过程的山区小流域泥石流灾害风险评估研究D.成都：成都理工大学.余斌,2 0 0 8.粘性泥石流的平均运动速度研究J.地球科学进展，2 3(5):52 4-532.杨超,2 0 19.雅鲁藏布江色东普沟冰崩碎屑流演化特征及成因机理研究D.
29、成都：成都理工大学.钟华邦,2 0 10.地质素描：西藏巴松错堰塞湖J.地质学刊，34(1):103.詹威威，黄润秋，裴向军,等,2 0 17.沟道型滑坡-碎屑流运动距离经验预测模型研究J.工程地质学报，2 5（1)：154-163.Biscarini C,di Francesco S,Manciola P,2010.CFD modellingapproach for dam break flow studies J.Hydrology andEarth System Sciences,14(4):705-718.Berger M J,George D L,LeVeque R J,et al.
30、,2011.The GeoClawsoftware for depth-averaged flows with adaptive refinementJ.Advances in Water Resources,34(9):1195-1206.Gray J M N T,2001.Granular flow in partially filled slowly rota-ting drums J/OL.Journal of Fluid Mechanics,441:1-29.https:/sci-hub.se/10.1017/s0022112001004736.(1.Regional Geologi
31、cal Survey Team of Tibet Geological Exploration Bureau,Lhasa 851414,Tibet,China;2.Geological Environment Survey and Research Center of Sichuan Province,Chengdu 610036,Sichuan,China)Abstract:To reveal the dynamic evolution process of the glacier avalanche-river blocking-bursting disaster chain in Sed
32、ongpu gully,Yarlung Zangbo River on October 17,2018,Massflow numerical simulation platform and Fortran language were used to simulate dy-namic process of glacier avalanche and debris flow by conducting secondary development and optimization of the Voellmy model accord-ing to the characteristics of r
33、egional geological conditions.The barrier lake dam obtained from simulated debris flow model was embed-ded into the terrain model,its scope and volume were calculated by ArcGIS,and the dynamic process of its bursting flood was simula-ted by Manning model.The complete dynamic process of glacier avala
34、nche-debris flow-river blocking-barrier lake-dam bursting wasreproduced by piecewise simulation method,and the parameters of flow velocity,flow depth and dam height in the process of debrisflow,and flow depth and flow velocity of bursting flood were quantitatively studied,to provide effective suppor
35、t for disaster preventionand reduction in Sedongpu Basin.Key words:glacier avalanche;debris flow;barrier lake;disaster chain;Massflow numerical simulation;Sedongpu gully;YarlungZangbo River地质学刊George D L,2011.Adaptive finite volume methods with well-balanced Riemann solvers for modeling floods in ru
36、gged ter-rain:Application to the Malpasset dam-break flood(France,1959)J.International Journal for NumericalMethods in Fluids,66(8):1000-1018.George D L,Iverson R M,2014.A depth-averaged debris-flowmodel that includes the effects of evolving dilatancy.II.Nu-merical predictions and experimental tests
37、 J/OL.Pro-ceedings of the Royal Society A:Mathematical,Physicaland Engineering Sciences,470(2170):1-31.https:/royalsocietypublishing.org/doi/10.1098/rspa.2013.0820.Guo Deping,Hamada M,He Chuan,et al.,2014.An empiricalmodel for landslide travel distance prediction in Wenchuanearthquake area J.Landsli
38、des,11(2):281-291.Iverson R M,Logan M,LaHusen R G,et al.,2010.Theperfect debris flow?Aggregated results from 28 large-scaleexperiments J/OL.Journal of Geophysical Research:Earth Surface,115:1-29.https:/ Chaojun,He Siming,Tang Chuan,2015.Numericalanalysis of dynamics of debris flow over erodible beds
39、 inWenchuan earthquake-induced area J.EngineeringGeology,194:62-72.Savage S B,Hutter K,1989.The motion of a finite mass ofgranular material down a rough incline J.Journal ofFluid Mechanics,199:177-215.Numerical simulation of glacier avalanche-river blocking-outburstin Sedongpu gully of the Yarlung Zangbo RiverTian Yangyang,Jiang Liang,Guo Jiang2023年
