论文部分内容阅读
摘要： 堰塞坝冲刷溃决及溃决洪水演进过程十分复杂，其溃决洪水对下游人民生命财产构成巨大威胁。利用数值分析方法对大型滑坡堰塞坝的溃决演进过程进行模拟和重演，对堰塞湖下游的避险与防灾减灾具有重要指导意义。以2018年金沙江“10·11”白格滑坡堰塞湖为例，基于无人机获取的地形数据，建立白格滑坡堰塞坝的三维数值模型，采用Flow-3D软件对堰塞坝的自然泄流冲刷溃决过程进行模拟，分析泄流槽内的流速、冲淤变化特征以及下游溃口处的洪峰流量演变过程。模拟结果表明：堰塞坝漫顶冲刷可以划分为溃决冲刷前、溃口快速拓展阶段、洪峰时刻、溃口稳定发展阶段4个时间段;溃决泄流过程中，泄流槽斜坡道上的水流流速较大，冲刷深度最大，堰塞坝下游出现明显淤积;白格堰塞湖溃决过程中出现了明显的溯源侵蚀现象，在泄流槽不断下切的过程中，泄流槽跌坎不断向上游移动。模拟结果有助于进一步深化对金沙江“10·11”白格滑坡堰塞坝冲刷溃决过程和机理的认识，对于堰塞湖应急处置措施和科学避险方案的制定具有一定的参考价值。
关 键 词： 堰塞坝; 冲刷溃决; 洪峰流量; 溃口演变; 数值模拟; 白格堰塞湖
中图法分类号： P642
文献标志码： A
DOI： 10.16232/j.cnki.1001-4179.2021.08.004
0 引 言
堰塞坝是山体滑坡、崩塌、泥石流等堵塞河道而形成的一种天然坝体[1]。堰塞坝由于形成时间短，坝体结构较为松散，在水流的淘蚀和冲刷作用下容易出现突然溃决，这将导致堰塞湖内大量蓄水迅速下泄而引发洪水灾害，严重威胁着下游人民的生命和财产安全[2-3]。堰塞坝溃决事件在世界各地时有发生，往往造成严重的经济损失和人员伤亡。据记载，1933年10月9日，四川叠溪地震形成的“海子”堰塞坝溃决，洪水倾泻而下，水头超过60 m，冲毁下游农田超过3 300 hm2，冲没房屋、牲畜不计其数，总死亡人数高达2 500余人[4]。1986年9月，巴布亚新几内亚贝拉曼河上的堰塞坝发生漫顶溃决，约1.2亿m3的泥石流下泄39 km汇入所罗门海，导致入海口水位较正常水位高出8 m，貝拉曼村被淹没[5]。2008年汶川地震共计诱发了256处滑坡堰塞湖，其中唐家山堰塞湖严重威胁下游绵阳市和遂宁市两座城市城区的安全，绵竹市绵远河上的小岗剑堰塞湖对下游数万人的生命财产安全构成极大威胁[6]。
开展堰塞坝溃决过程及机理研究对于制定科学合理的应急抢险方案具有重大意义，可最大程度降低下游的灾害损失。目前，已有众多研究者针对这一热点问题开展了大量研究。张婧等[7]通过堰塞坝泄流冲刷水槽试验研究了溃坝过程，探讨了不同内坡坡度和不同坝体级配对洪峰流量的影响。王光谦等[8]建立了耦合边界跟踪技术的平面二维水沙数学模型，并考虑了河床变形对泄流的影响。朱勇辉等[9]将土石坝溃坝模型BRESZHU应用于唐家山堰塞湖的泄流冲刷模拟，所计算的数十种工况结果为制定抢险方案提供了技术支持。目前研究溃坝过程的数值模型大多为二维模型，忽略了河道、复杂地形的平面弯曲和横截面的变化，且对溃口形状做了一定的假设，无法较为真实地描述溃坝水流和溃口发展的复杂性。
2018年10月11日和11月3日，西藏江达县波罗乡白格村与四川甘孜州白玉县交界处先后发生两次山体滑坡，滑体堵塞金沙江并形成堰塞湖。白格堰塞湖泄流时出现较大洪峰，堰塞坝下游在建的叶巴滩、苏洼龙等大型水电站均遭受较大影响。灾害发生后，许多研究者对白格堰塞湖进行了现场调查，并开展了广泛的研究。蔡耀军等[10]基于现场调查数据介绍分析了白格堰塞坝的成因、形态特征及溃流过程，研究表明白格堰塞湖库容、来水量、堰塞体溃流段长度及抗冲性是影响溃决特征值的关键因素。陈祖煜、侯精明等人[1112]运用DB-IWHR溃坝洪水分析程序和GST洪水演进模型比较准确地反演了白格堰塞坝溃口流量过程曲线及洪水演进过程。Zhang等[13]通过野外调查确定了白格堰塞坝坝体材料的可侵蚀性参数，并采用经验方程和漫顶溃坝数值模型DABA对白格堰塞湖两次溃决过程的库水位变化、洪水流量以及溃口尺寸做出了快速合理的预测。Fan等[14]利用数值分析方法对白格堰塞湖的滑坡堵江过程及溃决洪水演进过程进行了模拟。
白格堰塞湖现有研究在模拟溃坝过程时均没有采用实际地形，而是对坝体形态做了一定程度的概化，其采用的溃坝模型限定了溃口发展变化，不能较为真实地描述溃口变化的复杂性。本文基于无人机获取的高精度地形数据，构建了“10·11”白格堰塞坝的三维数值模型，采用Flow-3D软件模拟了白格滑坡堰塞坝自然泄流的冲刷溃决过程，深入分析了溃口处的流速特征及溃口冲刷情况，以期进一步揭示溃口发展规律，为今后堰塞坝的应急处置及下游避险提供技术参考。
1 金沙江“10·11”白格滑坡堰塞湖概况
2018年10月11日凌晨，西藏江达县波罗乡白格村与四川甘孜州白玉县交界处金沙江河道右岸（西藏岸）发生大规模山体滑坡，滑体堵截金沙江干流，形成堰塞湖。白格滑坡属高位高速基岩滑坡，滑源区顶部高程高达3 700 m左右，坡脚金沙江河床高程约2 900 m。破碎的岩体以巨大的势能俯冲而下，沿途铲刮斜坡表面的覆盖层和岩体，并在滑入金沙江后继续向左岸（四川侧）冲击，受左岸山体抵挡最终在河道里堆积，形成堰塞坝。滑坡体瞬间滑入金沙江时激起巨大涌浪，涌浪挟带河床砂卵石冲刷席卷左岸，大片植被和表土层被冲刷殆尽，形成涌浪铲刮区[15]。白格堰塞坝顺河道方向长约1 100 m，宽约400～500 m，总体积约2 795万m3，主要物质组成为碎石土夹块石，土石比约为7 ∶3。堰塞体地形起伏较大，总体呈左侧高，右侧低，左岸最高高程2 985.0 m，右岸垭口最低高程2 931.4 m，最小坝高61 m[16]。
堰塞体形成后，金沙江上游水位急剧壅高，迅速逼近垭口高程，至2018年10月12日17：15（北京时间），湖水开始漫顶溢流。10月13日00：45，堰塞体上游来水量和下泄流量达到平衡，水位达到最高值2 932.69 m，库容达到最大值2.9亿m3，此后水位逐渐消退。10月13日06：00，溃口流量达到峰值约10 000 m3/s，10月13日14：00，溃口流量降至与上游来水量持平[10]。堰塞体溃决后，在靠近右岸处形成宽约80～120 m的天然泄流槽，堰塞坝下游坝体附近出现明显的泥沙淤积，在紧接左岸残留坝体下方形成冲刷淤积区，如图1所示。 2 数值模拟方案
Flow-3D是一款功能强大的计算流体力学软件，它以连续性方程和Navier-Stokes方程作为控制方程，采用基于结构化矩形网格的有限差分法和有限体积法进行求解。Flow-3D运用FAVORTM技术解决复杂几何问题，采用VOF方法追踪自由表面，具有模拟复杂地形、精确处理干湿边界和急缓流交替水流运动的能力，可以较为真实地模拟流体运动现象和计算流场的各种性质。
当上游河道水位不断增高，堰塞体逐渐开始过流，水流会不断冲蚀松散的堰塞堆积体，导致堰塞体发生溃决，堰塞体的溃决过程实质上是溃坝水流对堆积物颗粒的冲刷过程。Flow-3D自带的泥沙冲刷模型通过计算泥沙颗粒的起动、沉降以及推移质和悬移质的运输来描述泥沙运动，这一模型与堰塞坝溃决机理相适应，即每层堰塞体材料的冲蚀具体表现为该层泥沙颗粒运动情况。此算法对溃口发展过程不作任何限制，堰塞体是否冲刷完全取决于泥沙输移规律，以计算时间中止为最终冲刷结果，使模拟更加合理。此外，Flow-3D自带的RNG k-ε模型能精细地模拟低强度紊流和具有强剪切应力区域的水流流动，在计算复杂地形的泥沙冲刷问题时具有良好的适用性，因此能较好地模拟溃坝水流的复杂流态。
本文根据现场无人机勘测地形数据，按实际地形1 ∶1比例建立三维实体模型，并导入Flow-3D软件中。为了在不影响研究区范围的前提下提高计算效率，本文计算域选取包含堰塞体在内的长约1 170 m，宽约830 m，高约240 m的计算范围，如图2所示。本文采用结构化正交网格对计算区域进行划分，网格尺寸为5 m×5 m×5 m，在溃口局部区域进行网格加密，加密网格尺寸为2.5 m×2.5 m×2.5 m，计算网格总数约为360万个。
上游来水量、堰塞坝的库容和溃口发展情况是影响溃坝洪水过程线的重要因素，本文通过严格控制上游水位来消除前两个因素的影响。模型中将入口边界设定为压力边界，水位下降过程的设置参考现场数据如下：10月13日0：45，上游水位达到最大值2 932.69 m;02：45，上游水位下降0.88 m;06：30，水位累计下降8.51 m;08：00，上游水位累計下降12.26 m;09：00，水位累计下降20 m[10]。出口边界设定为自由出流边界，壁面和底面边界均采用固壁边界，自由液面上方设为一个标准大气压。堰塞体上游初始水位设为2 932.69 m。堰塞体物质组成直接影响了坝体的抗冲性和溃口发展过程，目前尚无“10·11”堆积物颗粒的详细资料，但由于“10·11”与“11·3”两次滑坡堰塞体堆积材料物质组成接近，本文参考“11·3”堰塞体颗粒级配曲线（见图3），据此确定数值计算模型的相关参数[10]。本文设定堰塞堆积体为泥沙区域，设置泥沙颗粒密度为1 980 kg/m3，临界Shields数为0.35，挟带系数为0.013，其他参数如表1所列。
3 计算结果分析
本文在堰塞坝下游设置流量监测断面，如图2所示。计算结果显示10月13日凌晨5：10达到了洪峰流量9 719 m3/s，野外实测数据显示溃口在10月13日6：00达到了洪峰流量10 000 m3/s，模拟结果与现场数据相比洪峰提前了50 min，峰值流量相差较小，相对误差小于3%，表明本文数值计算的可靠性和合理性。数值计算得到的溃决洪峰流量与实测结果对比如表2所列。
白格堰塞坝自然泄流漫顶冲刷是个非常复杂的过程，如图4所示。根据地面高程变化率及过流量大小可以将其划分为溃决冲刷前、溃口快速拓展阶段、洪峰时刻、溃口稳定发展阶段4个时间段。本文地面高程变化率指单位时间内河床高程的绝对变化，负值表示高程降低，反之表示高程增加，是一个定量表征洪水对河床冲刷效应的指标。从开始泄流到洪峰时刻，这一阶段的泥沙冲刷率较高，地形变化快，溃口快速拓展。洪峰过后，随着时间的推移，水位和流量逐渐降低，水流冲刷能力降低，溃口发展速度缓慢，地形逐渐趋于稳定。本文的计算由于模拟时长限制虽未能达到冲淤平衡，但仍为揭示堰塞坝的漫顶冲刷规律提供了参考。
3.1 流速分布
图5为不同时期漫坝水流的流速分布图。如图5所示，在形成新的流道后流速一直保持在较高的水平，在流道末段流速最大，最大值约为16 m/s。因流速较快，水流冲刷能力较强，导致溃口不断下切。溃口快速拓展阶段流道内溃决洪水流速最大值与洪峰时刻流道内最大值接近。洪峰过后，流道被显著下切变窄，且流道坡降变小，与此同时，随着泄流过程的持续，堰塞坝上游水位大幅下降，流速呈逐渐减小趋势，在这一阶段流速最大值约为13 m/s，流速最大值出现在流道中段和末段。
从图5可以看出，不同时期的流速分布规律具有一定的相似性：溃决水流在进入流道前，流速为0，整个库区处于静压平衡状态;水流进入流道后，泄流槽进口水流在自重作用下流速呈现明显递增。在溃口快速拓展阶段和洪峰时刻，泄流槽中段流速分布不均匀，流速值相对较小，水流挟沙能力较弱。在泄流槽末段，泄流槽斜坡道的水流在自重作用下速度急剧增大，水流冲刷能力强，挟带着该区域的泥沙形成高含沙水流下泄，泄流槽不断被侵蚀下切，至溃口稳定发展阶段，泄流槽宽度显著降低。在泄流槽不断下切的过程中，泄流槽跌坎不断向上游移动，泄流槽的中段和末段泄流洪水流速均较快。
3.2 溃口拓展过程
3.2.1 冲刷深度变化
水流沿泄流槽流动时对河床上的泥沙颗粒产生拖曳和上举作用，当水流强度超过泥沙临界起动条件时，泥沙开始运动，产生冲刷现象，堰塞体堆积高程随之减小，如图6所示。图6中，冲淤深度为正表示淤积，冲淤深度为负表示冲刷下切。随着冲刷过程的持续，整个堆积体的冲刷深度不断增加，冲刷深度大于10 m的范围也逐步增大。在溃口快速拓展阶段，流道连通并不断扩展，水流挟沙力强且含沙量大，该过程堰塞体的冲刷深度增幅大，冲刷速度快，流道底部下切和侧向侵蚀明显。溃口稳定发展阶段，泄洪流量较小，泄流槽内水面宽度和流深均大幅减小，堰塞体的冲刷深度增长较慢，流道底部和侧壁处于均匀冲刷态势，该现象与流速变化趋势相吻合。这一阶段泄流槽两侧已形成具有一定高度和坡度的侧壁，新流道河谷呈“U”型。泄流槽不同位置处的冲刷情况也有差异，流经泄流槽入口和中部的水流流速小，该区域河床出现了缓慢的冲刷。泄流槽斜坡道上的水流流速大，水流带动泥沙形成高含沙水流冲向下游，导致该区域的冲刷深度最大，堆积体高程显著降低。 在堰塞坝下游，地形坡度减缓，地势平坦开阔，被冲起的部分泥沙由于自重和床面摩擦力作用而沉降，出现泥沙淤积。在泄流初期，水流携着槽内泥沙颗粒输移到下游，且由于流速较小泥沙颗粒堆积在近堰塞体河段。如图6所示，发展至溃口快速拓展阶段，泄流槽一侧流出的水流扩散覆盖了整个河道，河道内有大范围的泥沙淤积，淤积深度达11 m。随着溃口和水流条件的不断变化，近堰塞体河段也会遭受冲刷，初期堆积的泥沙颗粒被冲向下游，河道的底坡更加平直，水流更加通畅。如图6所示，洪峰时刻较前一阶段的淤积范围更小，淤积深度变浅，冲刷程度更强。泄流后期，形成了较为稳定的新流道，流道内泥沙冲淤变化不明显，回水区泥沙不断淤积。白格堰塞体的泄流槽位于堰塞坝右侧，在堰塞坝下游河道左侧，溃决洪水形成涡旋，流速降低，洪水携带的大量泥沙逐渐沉降淤积，导致这一区域的河床明显抬高，最大淤积深度达到27.52 m。
3.2.2 溃口断面历时变化
为了探究白格堰塞体溃决过程，本文截取了泄流槽的纵、横剖面，以进一步揭示漫顶水流的冲刷规律。选取的截面位置如图7所示，其中A-A′截面与泄流槽中心截面大致重合，纵坐标为y=269 m;B-B′截面位于泄流槽斜坡道上，横坐标为x=709 m。
图8展示了堰塞体纵剖面（A-A′截面）不同时期的形态变化过程。如前文所述，堰塞坝整体范围内均有不同程度的冲刷，其中泄流槽斜坡道的冲刷现象最为显著。通过对比发现，白格堰塞湖第一次溃决过程中出现了明显的溯源侵蚀现象。漫顶水流对斜坡道冲蚀形成跌坎，水流流经跌坎处，一部分势能转化为动能，流速突然增大，水流挟沙力增大，导致跌坎处的溃决速度加快，泄流槽断面突增，泄流槽深度和宽度都较上游断面更大，在三维空间上形成一个朝向下游侧的喇叭状形态。跌坎上游，流态稳定，流速较小，冲刷强度低;跌坎下游，水流呈复杂的湍流形式，流速较大，水流对跌坎冲刷强烈，跌坎处砂石输移效率大，物质不断流失，在视觉上跌坎不断向上游区域推进。
泄流槽斜坡道的天然地形是产生溯源冲刷的有利条件，斜坡道的水面比降大，水位突然下降导致流速和紊动强度增大，水流侵蚀作用加强，对河床侵蚀并不断向上游发展。溯源侵蚀在短时间内产生强烈冲刷，一旦发生则发展迅速。相关研究表明，溯源侵蚀是堰塞体溃决的重要冲刷模式，是流道下切和拓展的主要形式[17]。
图9展示了溃口横断面（B-B′截面）不同时期的形态变化过程。在水流作用下，溃口处的冲刷深度不断增加。溃口快速拓展阶段的水位较高，流速较大，水流冲刷能力强，在溃口快速拓展阶段的末段，冲刷深度达19.9 m，溃口形态基本稳定。之后，溃口侧壁和底面呈均匀冲刷态势，流道发展缓慢。
白格堰塞湖现场调查数据显示，10月13日堰塞体溃决后形成新的流道，新流道位于堰塞体右侧凹槽处，紧贴右岸坡脚，宽约80～120 m[10]。如图9（d）所示，本文模拟的流道位置、宽度与实际情况相符。本文模拟中未考虑流道下切过程中两侧侧壁的垮塌，溃口发展以表面侵蚀、冲刷掏深为主，泥沙床面较为平整，横向变化不明显，而实际溃坝过程中，水流在卷起、携带砂石向下游运动的同时不断掏蚀泄流槽坡脚，泄流槽侧壁出现坍塌，坍塌体又被水流冲走，泄流槽宽度因此不断扩展。侧向侵蚀和侧壁坍塌是引起泄流槽横向拓宽的主要原因（图10，11）。
4 结 论
本文基于Flow-3D軟件对金沙江“10·11”白格堰塞湖自然泄流溃决过程进行了重演和分析，复演了白格堰塞坝的第一次溃决过程，模拟的洪峰流量过程和溃口横断面形态与实测值较为符合。本文研究取得如下结论。
（1） 金沙江“10·11”白格堰塞坝漫顶冲刷是一个复杂的时空变化过程，可以划分为溃决冲刷前、溃口快速拓展阶段、洪峰时刻、溃口稳定发展阶段4个时间段。
（2） 溃决泄流过程中，泄流槽斜坡道上的水流流速大，冲刷深度最大，堰塞坝下游出现明显淤积。
（3） 金沙江“10·11”白格堰塞湖溃决过程中出现了明显的溯源侵蚀现象，在泄流槽不断下切的过程中，泄流槽跌坎不断向上游移动。
参考文献：
[1] 任强，陈生水，钟启明，等.堰塞坝的形成机理与溃决风险[J].水利水电科技进展，2011，31（5）：30-34.
[2] 王杨科，孙欢，代娇娇.我国堰塞湖溃坝特点及其防治[J].科技信息，2013（19）：485-486，514.
[3] 匡尚富，汪小刚，黄金池，等.堰塞湖溃坝风险及其影响分析评估[J].中国水利，2008（16）：17-21.
[4] 高继宗.历史回顾：1933年四川叠溪地震[J].中国减灾，2008（6）：58-59.
[5] 徐方军，韩丽宇.部分国家和地区地震诱发堰塞湖及其抢险概况[J].水利发展研究，2008（6）：7-9.
[6] 崔鹏，陈晓清，张建强，等.“4.20”芦山7.0 级地震次生山地灾害活动特征与趋势[J].山地学报，2013，31（3）：257-265.
[7] 张婧，曹叔尤，杨奉广，等.堰塞坝泄流冲刷试验研究[J].四川大学学报（工程科学版），2010，42（5）：191-196.
[8] 王光谦，钟德钰，张红武，等.汶川地震唐家山堰塞湖泄流过程的数值模拟[J].科学通报，2008（24）：3130-3136.
[9] 朱勇辉，范北林，卢金友，等.唐家山堰塞湖溃坝洪水分析及泄流冲刷模拟[J].人民长江，2008，39（22）：79-82.
[10] 蔡耀军，栾约生，杨启贵，等.金沙江白格堰塞体结构形态与溃决特征研究[J].人民长江，2019，50（3）：15-22.
[11] 陈祖煜，张强，侯精明，等.金沙江“10·10”白格堰塞湖溃坝洪水反演分析[J].人民长江，2019，50（5）：1-4，19. [12] 侯精明，馬利平，陈祖煜，等.金沙江白格堰塞湖溃坝洪水演进高性能数值模拟[J].人民长江，2019，50（4）：8-11，70.
[13] ZHANG L M，XIAO T，HE J，et al.Erosion-based analysis of breaching of Baige landslide dams on the Jinsha River，China，in 2018[J].Landslides，2019，16：1965-1979.
[14] FAN X M，YANG F，SUBRAMANIAN S S，et al.Prediction of a multi-hazard chain by an integrated numerical simulation approach：the Baige landslide，Jinsha River，China[J].Landslides，2019，19：1313-1317.
[15] 许强，郑光，李为乐，等.2018年10月和11月金沙江白格两次滑坡-堰塞堵江事件分析研究[J].工程地质学报，2018，26（6）：1534-1551.
[16] LI H B，QI S C，CHEN H，et al.Mass movement and formation process analysis of the two sequential landslide dam events in Jinsha River，southwest China[J].Landslides，2019，16：2247-2258.
[17] 牛志攀，许唯临，张建民，等.堰塞湖冲刷及溃决试验研究[J].四川大学学报（工程科学版），2009，41（3）：90-95.
（编辑：郑 毅）
引用本文：
谢忱，陈骎，侯奇东，等.
金沙江白格堰塞坝自然泄流冲刷溃决过程数值模拟
[J].人民长江，2021，52（8）：22-29.
Numerical simulation of natural erosion and breaching process
of “10·11” Baige landslide dam on Jinsha River
XIE Chen1，CHEN Qin1，HOU Qidong1，HU Yuxiang1，ZHOU Jiawen1，2，FAN Gang1，2
（ 1.College of Water Resource and Hydropower，Sichuan University，Chengdu 610065，China; 2.State Key Laboratory of Hydraulics and Mountain River Engineering，Sichuan University，Chengdu 610065，China ）
Abstract：
Outburst flood of the landslide dam poses a great threat to the life and property of the people downstream.Also，the processes of dam erosion and flood evolution are very complicated.Thus，using numerical methods simulate and repeat the typical landslide dam breaching has a great significance for disaster prevention and mitigation in the downstream.This study focuses on the “10·11” Baige landslide dam on the Jinsha River in 2018.A three-dimensional numerical model of the Baige landslide dam is established with the terrain data acquired by UAV，and the processes of dam discharge，scour and breaching are simulated using the Flow-3D software.The velocity and sediment scour in the natural spillway are analyzed，and the peak discharge process at the downstream is measured.The results show that overtopping and scouring of the Baige landslide dam can be divided into four phases，i.e.before the outburst，rapid expansion phase，peak discharge phase and breach steadily-developing phase.The velocity of flow on the ramp of spillway was large during the releasing period，which led to a maximum scour depth，while deposition occurred at the downstream of the dam body.There was an obvious headward erosion in the breaching process，the spillway bed was incised continuously which led channel bed drop-off moving to the upstream.This study is significant for understanding the process and mechanism of the “10·11” Baige landslide dammed breaching，which could provide references to the emergency treatment and rescue plan formulation.
Key words：
landslide dam;erosion and breaching;peak discharge;breach evolution;numerical simulation;Baige landslide dammed lake
关 键 词： 堰塞坝; 冲刷溃决; 洪峰流量; 溃口演变; 数值模拟; 白格堰塞湖
中图法分类号： P642
文献标志码： A
DOI： 10.16232/j.cnki.1001-4179.2021.08.004
0 引 言
堰塞坝是山体滑坡、崩塌、泥石流等堵塞河道而形成的一种天然坝体[1]。堰塞坝由于形成时间短，坝体结构较为松散，在水流的淘蚀和冲刷作用下容易出现突然溃决，这将导致堰塞湖内大量蓄水迅速下泄而引发洪水灾害，严重威胁着下游人民的生命和财产安全[2-3]。堰塞坝溃决事件在世界各地时有发生，往往造成严重的经济损失和人员伤亡。据记载，1933年10月9日，四川叠溪地震形成的“海子”堰塞坝溃决，洪水倾泻而下，水头超过60 m，冲毁下游农田超过3 300 hm2，冲没房屋、牲畜不计其数，总死亡人数高达2 500余人[4]。1986年9月，巴布亚新几内亚贝拉曼河上的堰塞坝发生漫顶溃决，约1.2亿m3的泥石流下泄39 km汇入所罗门海，导致入海口水位较正常水位高出8 m，貝拉曼村被淹没[5]。2008年汶川地震共计诱发了256处滑坡堰塞湖，其中唐家山堰塞湖严重威胁下游绵阳市和遂宁市两座城市城区的安全，绵竹市绵远河上的小岗剑堰塞湖对下游数万人的生命财产安全构成极大威胁[6]。
开展堰塞坝溃决过程及机理研究对于制定科学合理的应急抢险方案具有重大意义，可最大程度降低下游的灾害损失。目前，已有众多研究者针对这一热点问题开展了大量研究。张婧等[7]通过堰塞坝泄流冲刷水槽试验研究了溃坝过程，探讨了不同内坡坡度和不同坝体级配对洪峰流量的影响。王光谦等[8]建立了耦合边界跟踪技术的平面二维水沙数学模型，并考虑了河床变形对泄流的影响。朱勇辉等[9]将土石坝溃坝模型BRESZHU应用于唐家山堰塞湖的泄流冲刷模拟，所计算的数十种工况结果为制定抢险方案提供了技术支持。目前研究溃坝过程的数值模型大多为二维模型，忽略了河道、复杂地形的平面弯曲和横截面的变化，且对溃口形状做了一定的假设，无法较为真实地描述溃坝水流和溃口发展的复杂性。
2018年10月11日和11月3日，西藏江达县波罗乡白格村与四川甘孜州白玉县交界处先后发生两次山体滑坡，滑体堵塞金沙江并形成堰塞湖。白格堰塞湖泄流时出现较大洪峰，堰塞坝下游在建的叶巴滩、苏洼龙等大型水电站均遭受较大影响。灾害发生后，许多研究者对白格堰塞湖进行了现场调查，并开展了广泛的研究。蔡耀军等[10]基于现场调查数据介绍分析了白格堰塞坝的成因、形态特征及溃流过程，研究表明白格堰塞湖库容、来水量、堰塞体溃流段长度及抗冲性是影响溃决特征值的关键因素。陈祖煜、侯精明等人[1112]运用DB-IWHR溃坝洪水分析程序和GST洪水演进模型比较准确地反演了白格堰塞坝溃口流量过程曲线及洪水演进过程。Zhang等[13]通过野外调查确定了白格堰塞坝坝体材料的可侵蚀性参数，并采用经验方程和漫顶溃坝数值模型DABA对白格堰塞湖两次溃决过程的库水位变化、洪水流量以及溃口尺寸做出了快速合理的预测。Fan等[14]利用数值分析方法对白格堰塞湖的滑坡堵江过程及溃决洪水演进过程进行了模拟。
白格堰塞湖现有研究在模拟溃坝过程时均没有采用实际地形，而是对坝体形态做了一定程度的概化，其采用的溃坝模型限定了溃口发展变化，不能较为真实地描述溃口变化的复杂性。本文基于无人机获取的高精度地形数据，构建了“10·11”白格堰塞坝的三维数值模型，采用Flow-3D软件模拟了白格滑坡堰塞坝自然泄流的冲刷溃决过程，深入分析了溃口处的流速特征及溃口冲刷情况，以期进一步揭示溃口发展规律，为今后堰塞坝的应急处置及下游避险提供技术参考。
1 金沙江“10·11”白格滑坡堰塞湖概况
2018年10月11日凌晨，西藏江达县波罗乡白格村与四川甘孜州白玉县交界处金沙江河道右岸（西藏岸）发生大规模山体滑坡，滑体堵截金沙江干流，形成堰塞湖。白格滑坡属高位高速基岩滑坡，滑源区顶部高程高达3 700 m左右，坡脚金沙江河床高程约2 900 m。破碎的岩体以巨大的势能俯冲而下，沿途铲刮斜坡表面的覆盖层和岩体，并在滑入金沙江后继续向左岸（四川侧）冲击，受左岸山体抵挡最终在河道里堆积，形成堰塞坝。滑坡体瞬间滑入金沙江时激起巨大涌浪，涌浪挟带河床砂卵石冲刷席卷左岸，大片植被和表土层被冲刷殆尽，形成涌浪铲刮区[15]。白格堰塞坝顺河道方向长约1 100 m，宽约400～500 m，总体积约2 795万m3，主要物质组成为碎石土夹块石，土石比约为7 ∶3。堰塞体地形起伏较大，总体呈左侧高，右侧低，左岸最高高程2 985.0 m，右岸垭口最低高程2 931.4 m，最小坝高61 m[16]。
堰塞体形成后，金沙江上游水位急剧壅高，迅速逼近垭口高程，至2018年10月12日17：15（北京时间），湖水开始漫顶溢流。10月13日00：45，堰塞体上游来水量和下泄流量达到平衡，水位达到最高值2 932.69 m，库容达到最大值2.9亿m3，此后水位逐渐消退。10月13日06：00，溃口流量达到峰值约10 000 m3/s，10月13日14：00，溃口流量降至与上游来水量持平[10]。堰塞体溃决后，在靠近右岸处形成宽约80～120 m的天然泄流槽，堰塞坝下游坝体附近出现明显的泥沙淤积，在紧接左岸残留坝体下方形成冲刷淤积区，如图1所示。 2 数值模拟方案
Flow-3D是一款功能强大的计算流体力学软件，它以连续性方程和Navier-Stokes方程作为控制方程，采用基于结构化矩形网格的有限差分法和有限体积法进行求解。Flow-3D运用FAVORTM技术解决复杂几何问题，采用VOF方法追踪自由表面，具有模拟复杂地形、精确处理干湿边界和急缓流交替水流运动的能力，可以较为真实地模拟流体运动现象和计算流场的各种性质。
当上游河道水位不断增高，堰塞体逐渐开始过流，水流会不断冲蚀松散的堰塞堆积体，导致堰塞体发生溃决，堰塞体的溃决过程实质上是溃坝水流对堆积物颗粒的冲刷过程。Flow-3D自带的泥沙冲刷模型通过计算泥沙颗粒的起动、沉降以及推移质和悬移质的运输来描述泥沙运动，这一模型与堰塞坝溃决机理相适应，即每层堰塞体材料的冲蚀具体表现为该层泥沙颗粒运动情况。此算法对溃口发展过程不作任何限制，堰塞体是否冲刷完全取决于泥沙输移规律，以计算时间中止为最终冲刷结果，使模拟更加合理。此外，Flow-3D自带的RNG k-ε模型能精细地模拟低强度紊流和具有强剪切应力区域的水流流动，在计算复杂地形的泥沙冲刷问题时具有良好的适用性，因此能较好地模拟溃坝水流的复杂流态。
本文根据现场无人机勘测地形数据，按实际地形1 ∶1比例建立三维实体模型，并导入Flow-3D软件中。为了在不影响研究区范围的前提下提高计算效率，本文计算域选取包含堰塞体在内的长约1 170 m，宽约830 m，高约240 m的计算范围，如图2所示。本文采用结构化正交网格对计算区域进行划分，网格尺寸为5 m×5 m×5 m，在溃口局部区域进行网格加密，加密网格尺寸为2.5 m×2.5 m×2.5 m，计算网格总数约为360万个。
上游来水量、堰塞坝的库容和溃口发展情况是影响溃坝洪水过程线的重要因素，本文通过严格控制上游水位来消除前两个因素的影响。模型中将入口边界设定为压力边界，水位下降过程的设置参考现场数据如下：10月13日0：45，上游水位达到最大值2 932.69 m;02：45，上游水位下降0.88 m;06：30，水位累计下降8.51 m;08：00，上游水位累計下降12.26 m;09：00，水位累计下降20 m[10]。出口边界设定为自由出流边界，壁面和底面边界均采用固壁边界，自由液面上方设为一个标准大气压。堰塞体上游初始水位设为2 932.69 m。堰塞体物质组成直接影响了坝体的抗冲性和溃口发展过程，目前尚无“10·11”堆积物颗粒的详细资料，但由于“10·11”与“11·3”两次滑坡堰塞体堆积材料物质组成接近，本文参考“11·3”堰塞体颗粒级配曲线（见图3），据此确定数值计算模型的相关参数[10]。本文设定堰塞堆积体为泥沙区域，设置泥沙颗粒密度为1 980 kg/m3，临界Shields数为0.35，挟带系数为0.013，其他参数如表1所列。
3 计算结果分析
本文在堰塞坝下游设置流量监测断面，如图2所示。计算结果显示10月13日凌晨5：10达到了洪峰流量9 719 m3/s，野外实测数据显示溃口在10月13日6：00达到了洪峰流量10 000 m3/s，模拟结果与现场数据相比洪峰提前了50 min，峰值流量相差较小，相对误差小于3%，表明本文数值计算的可靠性和合理性。数值计算得到的溃决洪峰流量与实测结果对比如表2所列。
白格堰塞坝自然泄流漫顶冲刷是个非常复杂的过程，如图4所示。根据地面高程变化率及过流量大小可以将其划分为溃决冲刷前、溃口快速拓展阶段、洪峰时刻、溃口稳定发展阶段4个时间段。本文地面高程变化率指单位时间内河床高程的绝对变化，负值表示高程降低，反之表示高程增加，是一个定量表征洪水对河床冲刷效应的指标。从开始泄流到洪峰时刻，这一阶段的泥沙冲刷率较高，地形变化快，溃口快速拓展。洪峰过后，随着时间的推移，水位和流量逐渐降低，水流冲刷能力降低，溃口发展速度缓慢，地形逐渐趋于稳定。本文的计算由于模拟时长限制虽未能达到冲淤平衡，但仍为揭示堰塞坝的漫顶冲刷规律提供了参考。
3.1 流速分布
图5为不同时期漫坝水流的流速分布图。如图5所示，在形成新的流道后流速一直保持在较高的水平，在流道末段流速最大，最大值约为16 m/s。因流速较快，水流冲刷能力较强，导致溃口不断下切。溃口快速拓展阶段流道内溃决洪水流速最大值与洪峰时刻流道内最大值接近。洪峰过后，流道被显著下切变窄，且流道坡降变小，与此同时，随着泄流过程的持续，堰塞坝上游水位大幅下降，流速呈逐渐减小趋势，在这一阶段流速最大值约为13 m/s，流速最大值出现在流道中段和末段。
从图5可以看出，不同时期的流速分布规律具有一定的相似性：溃决水流在进入流道前，流速为0，整个库区处于静压平衡状态;水流进入流道后，泄流槽进口水流在自重作用下流速呈现明显递增。在溃口快速拓展阶段和洪峰时刻，泄流槽中段流速分布不均匀，流速值相对较小，水流挟沙能力较弱。在泄流槽末段，泄流槽斜坡道的水流在自重作用下速度急剧增大，水流冲刷能力强，挟带着该区域的泥沙形成高含沙水流下泄，泄流槽不断被侵蚀下切，至溃口稳定发展阶段，泄流槽宽度显著降低。在泄流槽不断下切的过程中，泄流槽跌坎不断向上游移动，泄流槽的中段和末段泄流洪水流速均较快。
3.2 溃口拓展过程
3.2.1 冲刷深度变化
水流沿泄流槽流动时对河床上的泥沙颗粒产生拖曳和上举作用，当水流强度超过泥沙临界起动条件时，泥沙开始运动，产生冲刷现象，堰塞体堆积高程随之减小，如图6所示。图6中，冲淤深度为正表示淤积，冲淤深度为负表示冲刷下切。随着冲刷过程的持续，整个堆积体的冲刷深度不断增加，冲刷深度大于10 m的范围也逐步增大。在溃口快速拓展阶段，流道连通并不断扩展，水流挟沙力强且含沙量大，该过程堰塞体的冲刷深度增幅大，冲刷速度快，流道底部下切和侧向侵蚀明显。溃口稳定发展阶段，泄洪流量较小，泄流槽内水面宽度和流深均大幅减小，堰塞体的冲刷深度增长较慢，流道底部和侧壁处于均匀冲刷态势，该现象与流速变化趋势相吻合。这一阶段泄流槽两侧已形成具有一定高度和坡度的侧壁，新流道河谷呈“U”型。泄流槽不同位置处的冲刷情况也有差异，流经泄流槽入口和中部的水流流速小，该区域河床出现了缓慢的冲刷。泄流槽斜坡道上的水流流速大，水流带动泥沙形成高含沙水流冲向下游，导致该区域的冲刷深度最大，堆积体高程显著降低。 在堰塞坝下游，地形坡度减缓，地势平坦开阔，被冲起的部分泥沙由于自重和床面摩擦力作用而沉降，出现泥沙淤积。在泄流初期，水流携着槽内泥沙颗粒输移到下游，且由于流速较小泥沙颗粒堆积在近堰塞体河段。如图6所示，发展至溃口快速拓展阶段，泄流槽一侧流出的水流扩散覆盖了整个河道，河道内有大范围的泥沙淤积，淤积深度达11 m。随着溃口和水流条件的不断变化，近堰塞体河段也会遭受冲刷，初期堆积的泥沙颗粒被冲向下游，河道的底坡更加平直，水流更加通畅。如图6所示，洪峰时刻较前一阶段的淤积范围更小，淤积深度变浅，冲刷程度更强。泄流后期，形成了较为稳定的新流道，流道内泥沙冲淤变化不明显，回水区泥沙不断淤积。白格堰塞体的泄流槽位于堰塞坝右侧，在堰塞坝下游河道左侧，溃决洪水形成涡旋，流速降低，洪水携带的大量泥沙逐渐沉降淤积，导致这一区域的河床明显抬高，最大淤积深度达到27.52 m。
3.2.2 溃口断面历时变化
为了探究白格堰塞体溃决过程，本文截取了泄流槽的纵、横剖面，以进一步揭示漫顶水流的冲刷规律。选取的截面位置如图7所示，其中A-A′截面与泄流槽中心截面大致重合，纵坐标为y=269 m;B-B′截面位于泄流槽斜坡道上，横坐标为x=709 m。
图8展示了堰塞体纵剖面（A-A′截面）不同时期的形态变化过程。如前文所述，堰塞坝整体范围内均有不同程度的冲刷，其中泄流槽斜坡道的冲刷现象最为显著。通过对比发现，白格堰塞湖第一次溃决过程中出现了明显的溯源侵蚀现象。漫顶水流对斜坡道冲蚀形成跌坎，水流流经跌坎处，一部分势能转化为动能，流速突然增大，水流挟沙力增大，导致跌坎处的溃决速度加快，泄流槽断面突增，泄流槽深度和宽度都较上游断面更大，在三维空间上形成一个朝向下游侧的喇叭状形态。跌坎上游，流态稳定，流速较小，冲刷强度低;跌坎下游，水流呈复杂的湍流形式，流速较大，水流对跌坎冲刷强烈，跌坎处砂石输移效率大，物质不断流失，在视觉上跌坎不断向上游区域推进。
泄流槽斜坡道的天然地形是产生溯源冲刷的有利条件，斜坡道的水面比降大，水位突然下降导致流速和紊动强度增大，水流侵蚀作用加强，对河床侵蚀并不断向上游发展。溯源侵蚀在短时间内产生强烈冲刷，一旦发生则发展迅速。相关研究表明，溯源侵蚀是堰塞体溃决的重要冲刷模式，是流道下切和拓展的主要形式[17]。
图9展示了溃口横断面（B-B′截面）不同时期的形态变化过程。在水流作用下，溃口处的冲刷深度不断增加。溃口快速拓展阶段的水位较高，流速较大，水流冲刷能力强，在溃口快速拓展阶段的末段，冲刷深度达19.9 m，溃口形态基本稳定。之后，溃口侧壁和底面呈均匀冲刷态势，流道发展缓慢。
白格堰塞湖现场调查数据显示，10月13日堰塞体溃决后形成新的流道，新流道位于堰塞体右侧凹槽处，紧贴右岸坡脚，宽约80～120 m[10]。如图9（d）所示，本文模拟的流道位置、宽度与实际情况相符。本文模拟中未考虑流道下切过程中两侧侧壁的垮塌，溃口发展以表面侵蚀、冲刷掏深为主，泥沙床面较为平整，横向变化不明显，而实际溃坝过程中，水流在卷起、携带砂石向下游运动的同时不断掏蚀泄流槽坡脚，泄流槽侧壁出现坍塌，坍塌体又被水流冲走，泄流槽宽度因此不断扩展。侧向侵蚀和侧壁坍塌是引起泄流槽横向拓宽的主要原因（图10，11）。
4 结 论
本文基于Flow-3D軟件对金沙江“10·11”白格堰塞湖自然泄流溃决过程进行了重演和分析，复演了白格堰塞坝的第一次溃决过程，模拟的洪峰流量过程和溃口横断面形态与实测值较为符合。本文研究取得如下结论。
（1） 金沙江“10·11”白格堰塞坝漫顶冲刷是一个复杂的时空变化过程，可以划分为溃决冲刷前、溃口快速拓展阶段、洪峰时刻、溃口稳定发展阶段4个时间段。
（2） 溃决泄流过程中，泄流槽斜坡道上的水流流速大，冲刷深度最大，堰塞坝下游出现明显淤积。
（3） 金沙江“10·11”白格堰塞湖溃决过程中出现了明显的溯源侵蚀现象，在泄流槽不断下切的过程中，泄流槽跌坎不断向上游移动。
参考文献：
[1] 任强，陈生水，钟启明，等.堰塞坝的形成机理与溃决风险[J].水利水电科技进展，2011，31（5）：30-34.
[2] 王杨科，孙欢，代娇娇.我国堰塞湖溃坝特点及其防治[J].科技信息，2013（19）：485-486，514.
[3] 匡尚富，汪小刚，黄金池，等.堰塞湖溃坝风险及其影响分析评估[J].中国水利，2008（16）：17-21.
[4] 高继宗.历史回顾：1933年四川叠溪地震[J].中国减灾，2008（6）：58-59.
[5] 徐方军，韩丽宇.部分国家和地区地震诱发堰塞湖及其抢险概况[J].水利发展研究，2008（6）：7-9.
[6] 崔鹏，陈晓清，张建强，等.“4.20”芦山7.0 级地震次生山地灾害活动特征与趋势[J].山地学报，2013，31（3）：257-265.
[7] 张婧，曹叔尤，杨奉广，等.堰塞坝泄流冲刷试验研究[J].四川大学学报（工程科学版），2010，42（5）：191-196.
[8] 王光谦，钟德钰，张红武，等.汶川地震唐家山堰塞湖泄流过程的数值模拟[J].科学通报，2008（24）：3130-3136.
[9] 朱勇辉，范北林，卢金友，等.唐家山堰塞湖溃坝洪水分析及泄流冲刷模拟[J].人民长江，2008，39（22）：79-82.
[10] 蔡耀军，栾约生，杨启贵，等.金沙江白格堰塞体结构形态与溃决特征研究[J].人民长江，2019，50（3）：15-22.
[11] 陈祖煜，张强，侯精明，等.金沙江“10·10”白格堰塞湖溃坝洪水反演分析[J].人民长江，2019，50（5）：1-4，19. [12] 侯精明，馬利平，陈祖煜，等.金沙江白格堰塞湖溃坝洪水演进高性能数值模拟[J].人民长江，2019，50（4）：8-11，70.
[13] ZHANG L M，XIAO T，HE J，et al.Erosion-based analysis of breaching of Baige landslide dams on the Jinsha River，China，in 2018[J].Landslides，2019，16：1965-1979.
[14] FAN X M，YANG F，SUBRAMANIAN S S，et al.Prediction of a multi-hazard chain by an integrated numerical simulation approach：the Baige landslide，Jinsha River，China[J].Landslides，2019，19：1313-1317.
[15] 许强，郑光，李为乐，等.2018年10月和11月金沙江白格两次滑坡-堰塞堵江事件分析研究[J].工程地质学报，2018，26（6）：1534-1551.
[16] LI H B，QI S C，CHEN H，et al.Mass movement and formation process analysis of the two sequential landslide dam events in Jinsha River，southwest China[J].Landslides，2019，16：2247-2258.
[17] 牛志攀，许唯临，张建民，等.堰塞湖冲刷及溃决试验研究[J].四川大学学报（工程科学版），2009，41（3）：90-95.
（编辑：郑 毅）
引用本文：
谢忱，陈骎，侯奇东，等.
金沙江白格堰塞坝自然泄流冲刷溃决过程数值模拟
[J].人民长江，2021，52（8）：22-29.
Numerical simulation of natural erosion and breaching process
of “10·11” Baige landslide dam on Jinsha River
XIE Chen1，CHEN Qin1，HOU Qidong1，HU Yuxiang1，ZHOU Jiawen1，2，FAN Gang1，2
（ 1.College of Water Resource and Hydropower，Sichuan University，Chengdu 610065，China; 2.State Key Laboratory of Hydraulics and Mountain River Engineering，Sichuan University，Chengdu 610065，China ）
Abstract：
Outburst flood of the landslide dam poses a great threat to the life and property of the people downstream.Also，the processes of dam erosion and flood evolution are very complicated.Thus，using numerical methods simulate and repeat the typical landslide dam breaching has a great significance for disaster prevention and mitigation in the downstream.This study focuses on the “10·11” Baige landslide dam on the Jinsha River in 2018.A three-dimensional numerical model of the Baige landslide dam is established with the terrain data acquired by UAV，and the processes of dam discharge，scour and breaching are simulated using the Flow-3D software.The velocity and sediment scour in the natural spillway are analyzed，and the peak discharge process at the downstream is measured.The results show that overtopping and scouring of the Baige landslide dam can be divided into four phases，i.e.before the outburst，rapid expansion phase，peak discharge phase and breach steadily-developing phase.The velocity of flow on the ramp of spillway was large during the releasing period，which led to a maximum scour depth，while deposition occurred at the downstream of the dam body.There was an obvious headward erosion in the breaching process，the spillway bed was incised continuously which led channel bed drop-off moving to the upstream.This study is significant for understanding the process and mechanism of the “10·11” Baige landslide dammed breaching，which could provide references to the emergency treatment and rescue plan formulation.
Key words：
landslide dam;erosion and breaching;peak discharge;breach evolution;numerical simulation;Baige landslide dammed lake