堰塞体溃决流量与溃口形态演化数值模拟

Numerical Simulation of Breach Hydrograph and Morphology Evolution during Landslide Dam Breaching

摘要

堰塞体是山区常见的地质灾害，一旦发生溃决，将对下游人民的生命财产安全构成严重威胁.在应急响应时需要对堰塞体溃口流量过程和溃口形态演化进行快速准确的预测，但目前的溃坝数学模型大多未充分考虑堰塞体的地貌学特征，无法合理反映复杂地形下堰塞体的溃决过程.采用雷诺平均Navier-Stokes方程和湍流重正化群k-ε模型相结合的数值方法，对复杂地形下的溃决水流进行模拟，并利用可考虑推移质和悬移质输移的冲蚀公式模拟溃口冲蚀过程，选择拥有详细勘测资料和水文数据的“11•03”白格堰塞体溃决案例进行反演分析.对比计算和实测的溃口流量过程线、溃决过程水动力学特征及最终溃口断面形态发现，模拟结果可较好地反映堰塞体的实际溃决过程，验证了模型的合理性.

关键词

堰塞体 / 溃决过程 / 数值模拟 / 溃口流量 / 溃口形态

Key words

landslide dam / breach process / numerical simulation / breach hydrograph / breach morphology

0 引言

近年来，受气候变化和人类活动影响，世界各地极端天气和地质灾害的发生频率显著增加.堰塞体作为一种山区常见的地质灾害，一般是由于降水或地震导致山体产生滑坡，堆积物阻塞河谷后形成.由于缺乏泄洪排水设施，大部分堰塞体在形成后不久即发生溃决.美国地质调查局对73座堰塞体的寿命进行统计发现（Costa and Schuster， 1988），85%的堰塞体在形成1 a内会发生溃决，27%的堰塞体1 d内就会溃决.此外，有记录的堰塞体溃决事故中，90%以上都是由于水流漫顶诱发的（Zhong et al.， 2021）.一旦发生溃决，短时间内释放的巨大洪水将对下游居民的生命财产造成严重威胁，而溃口流量过程直接决定了溃决洪水的灾害程度（Fan et al.， 2019； Luo et al.， 2019）.2008年，“5.12”汶川地震形成了257座堰塞体，其中规模最大的是唐家山堰塞体，满蓄容积可达3.16亿m³，对下游170万人民造成了巨大的威胁（Fan et al.， 2018）（图1a）.通过人工干预，堰塞体于2008年6月10日溃决，溃口洪峰流量为6 500 m³/s，在堰塞体溃决之前，下游约25万居民被疏散.2018年10~11月，四川省与西藏自治区交界处的白格村同一位置连续发生两次滑坡，均直接截断金沙江，形成的堰塞体均在短期内发生溃决，两次灾害共造成西藏、四川、云南等地10.2万人受灾，8.6万人紧急转移安置（Zhong et al.， 2020；曹鹏等， 2021）（图1b）.鉴于堰塞体溃决产生的巨大威胁，如何科学预测堰塞体溃决过程是应急处置时面临的关键难题.

为了模拟堰塞体漫顶溃决过程，近几十年来学者们开展了一系列的研究，主要围绕堰塞体的溃决机理与过程试验和溃决过程数值模拟（Peng and Zhang， 2012； Chen et al.， 2020）.对于溃决机理与过程试验，一般采用的方法可分为小型水槽模型试验（Zhu et al.， 2021）、大型现场模型试验（Zhang et al.， 2021）、离心模型试验（赵天龙等， 2016）和实体原位试验（Cai et al.， 2020）.此外通过量化各要素明确其对堰塞体溃决过程的影响，影响要素主要包括：堰塞体形态（Walder et al.， 2015）、坝料的物理力学性质（Zhu et al.， 2021）以及堰塞湖水动力条件（Jiang and Wei， 2020）.

对于溃决过程数值模拟，一般采用的模型可分为参数模型、基于物理过程的简化数学模型和基于物理过程的精细化数学模型（ASCE/EWRI Task Committee on Dam/Levee Breaching， 2011）.参数模型通常是结合历史溃决案例，运用代表性地貌学参数，采用回归方法建立基于经验参数的溃决模型，得到堰塞体溃决特征参数（如溃口峰值流量、溃口尺寸、溃决历时（Costa and Schuster， 1988； Peng and Zhang， 2012； 石振明等， 2014）），但无法提供溃口洪水流量过程.之后学者们又从溃决机理出发，提出了基于物理过程的堰塞体溃决过程简化数学模型.目前此类简化模型应用最为广泛，如香港科技大学的DABA模型（Chang and Zhang， 2010； Zhang et al.， 2019）、中国水利水电科学研究院的DB-IWHR模型（Chen et al.， 2015； Cai et al.， 2020）和南京水利科学研究院的DB-NHRI模型（Zhong et al.， 2020； Mei et al.， 2021）等.这类模型一般假定堰塞体为规则形状（如纵断面为正梯形、横断面为倒梯形），假定初始溃口也为规则形状（如三角形或倒梯形），在溃决过程溃口仅发生尺寸变化，但形状始终保持规则；采用宽顶堰流量公式模拟溃口水流，采用土体冲蚀率公式以及边坡极限平衡分析法分别模拟坝料的连续性侵蚀和溃口边坡的间歇性失稳坍塌.上述简化模型虽然计算效率较高，并在一定程度上反映了堰塞体的溃决机理，但在溃口几何形态、坝料物理力学特性、溃口流量表达和溃口边坡稳定性分析等方面引入了大量假设，无法真正反映复杂地形下堰塞体的溃决特征.随着泥沙科学和计算流体力学的发展，学者们基于水动力和泥沙输移方程，在浅水和静水压力分布假设下，建立了一系列一维、二维和三维数值模型，对溃坝过程进行了精细化模拟（陈生水等， 2019；李坤等， 2022）.这类模型一般包括3个模块，如水动力模块（清水或浑水的质量守恒和动量守恒方程）、坝料冲蚀模块（平衡或非平衡冲蚀方程）和溃口形态演化模块（挟带和沉积质量守恒方程、边坡失稳方程）.根据溃坝模型中冲蚀方程的不同，精细化模型大体可分为平衡模型、非平衡模型、两相流模型和两层流模型（Guan et al.， 2014）.

国内外学者虽然针对堰塞体的溃决机理和溃决过程模拟开展了大量的研究工作，但在溃口演化规律上尚存在争议，而可模拟堰塞体溃决过程的精细化模型较少，且大多为一维或二维模型（Marsooli and Wu， 2015）.本文基于堰塞体溃决机理与溃决过程的研究现状，梳理总结堰塞体溃决过程中的溃口演化规律，基于此选择合适的冲蚀公式和溃坝精细化模型，并选择具备详细勘察资料的金沙江“11•03”白格堰塞体溃决案例，考虑实际地形特征，进行三维精细化数值模拟.通过实测和计算的溃口流量过程与最终溃口形态的对比，验证了数值模型的准确性，研究成果可为堰塞体溃决应急抢险方案的制定提供科学依据.

1 堰塞体溃决特征与溃决过程

在开展堰塞体溃决过程数值模拟之前，深入认识其溃决过程特征及溃决机理对数学模型的构建尤为重要.由于不同的形成条件，堰塞体堆积物通常表面起伏不平，材料颗粒级配宽泛，包含大块石以及碎石土混合物等（Mei et al.， 2021）.漫顶溃决通常发生在坝顶高程较低部位或下游坝坡表面的薄弱部位，经水流侵蚀形成初始溃口.随着流量增加，溃口持续发生下切和扩展（图2）.地质调查和冲蚀试验研究表明堰塞体材料的冲蚀特性和临界剪应力沿深度方向均会发生变化（Zhang et al.， 2019），溃坝时溃口存在材料逐渐粗化的现象，在实际溃决后可能存在残留坝高.

基于模型试验和现场观测，堰塞体的漫顶溃决过程大致可分为初始、加速和稳定3个阶段（Jiang et al.， 2020； Zhong et al.， 2021； Zhu et al.， 2021）.考虑到溃口流量和溃口形态在纵横向的变化，溃决过程可描述如下：（1）表层冲蚀和溯源冲蚀阶段.当持续上涨的水位漫过坝顶，下游坝坡会形成大量细小冲沟.此时由于流速较小，表层冲蚀占据主导地位，仅表层细颗粒被水流带走.随着堰塞湖水位不断抬升，漫顶水流在重力势能的作用下，使得下游坡脚的冲沟迅速下切，形成初始溃口.随后，冲蚀作用逐渐由下游坡脚延伸至上游坝顶溃口，沿下游坡脚向上溃口的深度和宽度逐渐增大.（2）沿程冲蚀阶段.当溯源冲蚀发展到坝顶后，溃口底部高程陡增，漫顶水头突然增加，导致溃口流量迅速增大，漫顶水流顺流而下，溃口不断下切并导致溃口边坡发生失稳坍塌，进一步使溃口尺寸和流量迅速增加.这一阶段溃口发展速度最快，并出现溃口峰值流量.（3）溃口稳定阶段.随着溃口不断加深扩宽、溃口坝料的粗化和上游堰塞湖水位的下降，溃口纵向下切逐渐停止，以横向缓慢扩展为主.当漫顶水流的水动力衰弱到一定程度后，溃口处坝料将无法输移，悬浮物质逐渐沉降到堰塞体下游，使得残余堰塞体坡度逐渐平缓.最后当溃口流量与上游来流量相等时，整个溃决过程结束.

图3和图4分别给出了唐家山堰塞体和“11•03”白格堰塞体溃决前后的地形对比（Fan et al.， 2018； Cai et al.， 2020）.通过两座实际堰塞体溃决前后纵断面形态可知，由于坝料受到水流冲刷和重力沉降的联合作用，溃决后两座堰塞体的下游坡比均相对平缓，这也间接验证了上述的溃决机理.

2 三维堰塞体溃决过程数值模型

堰塞体具有复杂的外部形态，且溃口水流呈现出湍流的特点.因此，采用雷诺平均Navier-Stokes（RANS）方程和湍流重正化群k-ε（RNG k-ε）模型相结合的方法，对复杂地形下溃口水动力特征进行分析.此外，利用可考虑推移质和悬移质输移的坝料冲蚀方程表示溃口形态演化过程，通过流体体积（VOF）法追踪自由液面流动，应用有限体积法离散控制方程进行求解.

2.1　水动力模型

基于笛卡尔坐标系，采用RANS方程描述三维不可压缩流体运动，连续性方程如下：

\frac{\partial {\rho }_{\mathrm{w}}}{\partial t} V_{\mathrm{F}} + \frac{\partial }{\partial x} {\rho }_{\mathrm{w}} u A_x + \frac{\partial }{\partial y} {\rho }_{\mathrm{w}} v A_y + \frac{\partial }{\partial z} {\rho }_{\mathrm{w}} w A_z = \mathrm{0},(1)

式中：ρ _w为清水密度；t表示时间；V _F为流体通过单元的体积与总体积的比值；u、v、w分别为x、y、z方向上的流速；A_x 、A_y 、A_z 分别为x、y、z方向上水流流过的区域占比.

根据连续介质力学原理，动量方程可以表示为（Movahedi et al.， 2018）：

\left\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}(u{A}_x\frac{\partial u}{\partial x}+v{A}_y\frac{\partial u}{\partial y}+w{A}_z\frac{\partial u}{\partial z})=-\frac{\mathrm{1}}{{\rho }_w}\frac{\partial P}{\partial x}+G_x+f_x\\ \frac{\partial v}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}(u{A}_x\frac{\partial v}{\partial x}+v{A}_y\frac{\partial v}{\partial y}+w{A}_z\frac{\partial v}{\partial z})=-\frac{\mathrm{1}}{{\rho }_w}\frac{\partial P}{\partial y}+G_y+f_y\\ \frac{\partial w}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}(u{A}_x\frac{\partial w}{\partial x}+v{A}_y\frac{\partial w}{\partial y}+w{A}_z\frac{\partial w}{\partial z})=-\frac{\mathrm{1}}{{\rho }_w}\frac{\partial P}{\partial z}+G_z+f_z\end{array}\right.，(2)

式中：P为压强；G_x 、G_y 、G_z 分别为x、y、z方向上的质量加速度；f_x 、f_y 、f_z 分别为x、y、z方向上的黏性加速度，其中黏性加速度f_i 的张量形式可以表示为f_i =［τ _{b ，i} /ρ _w-∂（A_jS_ij ）/∂x_j ］/V _F，其中τ _b为壁面剪应力；应变率张量S_ij 可表示为S_ij =-（μ+μ_t ）（∂u_i /∂x_j +∂u_j /∂x_i ），μ为湍流黏性系数；μ _t为动态湍流黏性系数.

当溃坝洪水在堰塞体下游斜坡蔓延冲蚀时，水流变化剧烈，导致流态复杂.而标准的k-ε模型不能解决具有高剪切力和曲折流线的流体运动的计算，因此本文选用RNG k-ε湍流模型来模拟湍流对时均动量输运的影响，使用湍流黏性系数和动态湍流黏性系数将非线性雷诺应力与平均流速梯度联系起来，其基本控制方程为（Yakhot et al.， 1992）：

\frac{\partial ({\rho }_{\mathrm{w}}k)}{\partial t} + \frac{\partial ({\rho }_{\mathrm{w}}k{u}_i)}{\partial x_i} = \frac{\partial }{\partial x_j} \left[\left(\mu +\frac{{\mu }_{\mathrm{t}}}{{\sigma }_{\mathrm{k}}}\right)\frac{\partial k}{\partial x_j}\right] + G_k - {\rho }_w \epsilon，(3)

\frac{\partial ({\rho }_{\mathrm{w}}\epsilon )}{\partial t} + \frac{\partial ({\rho }_{\mathrm{w}}\epsilon u_i)}{\partial x_i} = \frac{\partial }{\partial x_j} \left[\left(\mu +\frac{{\mu }_{\mathrm{t}}}{{\sigma }_{\mathrm{\epsilon }}}\right)\frac{\partial \epsilon }{\partial x_j}\right] + C_{\mathrm{1}} G_k \frac{\epsilon }{k} - C_{\mathrm{2}} {\rho }_{\mathrm{w}} \frac{{\epsilon }^{\mathrm{2}}}{k}，(4)

式中：k为湍流动能；ε为湍流耗散率；此时μ _t=C_μk ²/ε，其中C_μ =0.085；G_k 为由速度梯度产生的湍流动能；C ₁=1.42，C ₂=1.68，σ_k =σ_ε =0.719 4.

2.2　坝料冲蚀模型

在堰塞体溃决过程中，坝料的运动状态可以在推移质和悬移质之间进行转换.图5给出了坝料冲蚀模型的示意图.其中水流中含沙量的变化会导致其密度和黏滞系数的变化，从而影响冲蚀过程.悬移质通常浓度较低，可以通过流体流动进行输送.由于相邻颗粒的限制，推移质往往不易发生移位，堆积的颗粒可以通过滚动、跃移、滑动或拉拽以推移质的形式移动（Samma et al.， 2020）.当沉积速度超过冲蚀速度时，悬移质可以转化为推移质.溃口底床的形态演化可通过挟带和沉积质量守恒方程来表示：

\varphi \frac{\partial \mathrm{z}}{\partial t} = ( \frac{\partial q_{bx}}{\partial x} + \frac{\partial q_{by}}{\partial y} + D - E )，(5)

式中：z为底床高程；q_bx 和q_by 分别为x和y方向单宽体积推移质输移率；ϕ为最大堆积分数（一般取0.64）；D为坝料向下沉积的速度；E为坝料向上挟带的速度.

E = {\alpha }_i n_s d_{\mathrm{*}}^{\mathrm{0.3}} ( {\theta }_i - {\theta }_{cr,i} {)}^{\mathrm{1.5}} {\left[\mathrm{g}(\frac{{\rho }_{\mathrm{s}}-{\rho }_{\mathrm{w}}}{{\rho }_{\mathrm{w}}})d_i\right]}^{\mathrm{0.5}} (6)

由于单个颗粒周围的水动力和界面处的边界层很难计算，因此采用经验模型来进行坝料冲蚀模拟.其中坝料的挟带和沉积可以看作是两个相互对立的微观过程，综合计算可得到推移质和悬移质之间的净交换率.坝料的向上挟带速度E可用来计算推移质转化为悬移质的量（Mastbergen and van den Berg， 2003）：

式中：α_i 为不同坝料的挟带率系数；n _s为充填溃口底床的法向矢量，d _*为无量纲不同坝料的粒径参数，d _*=d_i ［ρ _w（ρ _s-ρ _w）g/μ ²］^1/3；d_i 为不同坝料的平均粒径；ρ _s为坝料密度；θ_i 为无量纲希尔兹数，θ_i =τ_c /［gd _i （ρ _s-ρ _w）］，由剪切应力τ _c计算得到，τ _c=ρgu ²/C _D，其中C _D是与平均粒径相关的拖曳系数，C _D=5.75g^0.5log（2h/d ₅₀），h为局部水深.采用Soulsby-Whitehouse方程计算得到无量纲临界希尔兹数θ_cr，i，θ_cr，i =0.3/（1+1.2d _*）+0.055［1-e^{-0.02 d *}］（van Rijn， 2020）.

坝料冲蚀造成的溃口形态在三维空间的演化会导致溃口边坡的失稳，本文模型采用坝料的休止角来分析边坡的稳定性.在进行溃口形态演化计算时，当相邻网格中心点的连线形成的倾角大于坝料休止角时，溃口边坡即发生失稳，失稳后的倾角恢复到坝料休止角.另外，模型基于坝料的休止角对无量纲临界希尔兹数进行了修正（Samma et al.， 2020）：

{{\theta }_{cr,i}}^{\mathrm{\text{'}}} = {\theta }_{cr,i} \frac{\mathrm{c}\mathrm{o}\mathrm{s}\psi \mathrm{s}\mathrm{i}\mathrm{n}\beta +\sqrt[]{\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\beta \mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}{\phi }_i-\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\psi \mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta }}{\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_i}，(7)

式中：θ_{cr ， i} ʹ为修正后的无量纲临界希尔兹数；β为溃口底床坡角；φ_i 为不同坝料的休止角，ψ为水流与上游坡方向的夹角.

坝料向下沉积速度D是颗粒因重量而从悬移质沉降在推移质中的输移过程，可定义为有效沉积速率ω与悬移质浓度c的乘积：

D = {\omega }_i c_i，(8)

式中：ω_i 为不同坝料的有效沉积速度，ω_i =ν_f ［（10.36²+1.049d _* ³）^0.5-10.36］/d_i，其中ν _f为流体的运动黏度；c_i 为不同坝料的悬移质浓度.

推移质输移公式可分为经验和半经验公式，其中经验模型大多是在试验的基础上采用孤立因子法建立（Meyer-Peter and Muller， 1948），而半经验模型通常根据某种物理假设来确定公式的基本结构，公式中的一些参数可能需要利用实测数据来标定（Bagnold， 1966）.根据钱宁（1980）的研究，可采用下式表示无量纲推移质输沙率：

\mathrm{\Phi } = \left\{\begin{array}{cc}K{\mathrm{\Theta }}^a(\mathrm{\Theta }-{\mathrm{\Theta }}_c{)}^b{(\sqrt[]{\mathrm{\Theta }}-\lambda \sqrt[]{{\mathrm{\Theta }}_c})}^c& \mathrm{\Theta }>{\mathrm{\Theta }}_c\\ \mathrm{0}& \mathrm{\Theta }\le {\mathrm{\Theta }}_c\end{array}\right.，(9)

式中：Φ为无量纲推移质冲蚀率，Φ=q_b （ρ _w/（ρ _s-ρ _w）gd_i ）^0.5/ρ _sg；Θ为希尔兹系数的另一种表达；Θ _c 为颗粒起动系数，Θ _c =0.047；K为推移质系数；a，b，c和λ为相关参数，a+b+c=1.5.

本文选取了得到广泛引用的推移质实测数据，包括Meyer-Peter and Muller（1948）、钱宁（1980）等经典试验数据，以及Roseberry et al.（2012）通过高速摄影技术获得的推移质输沙率数据.这些实测数据已经过诸多学者分析验证，具有较高的可信度.颗粒粒径范围为0.785~28.650 mm，比重范围为1.25~4.22.为了选择合适的推移质公式用于计算，选择了这些实测数据对具有代表性的4个推移质公式进行了对比分析，如图6所示.由图可知，Meyer-Peter和Muller公式在弱输沙阶段（即水流参数ψ=1/Θ>2）预测效果较为准确，而在强输沙阶段（即水流参数ψ=1/Θ<2），也表现出较好的适应性.考虑到堰塞体的宽级配特征，选取Meyer-Peter和Muller的经验公式表示推移质单宽输沙率（Meyer-Peter and Muller， 1948）：

q_b = K {\left({\theta }_i-{{\theta }_{cr,i}}^{\mathrm{\text{'}}}\right)}^{\mathrm{1.5}} {\left[\mathrm{g}\left(\frac{{\rho }_s-{\rho }_{\mathrm{w}}}{{\rho }_{\mathrm{w}}}\right)d_i^{\mathrm{3}}\right]}^{\mathrm{0.5}} .(10)

要计算推移质输移情况，还需要估算推移质层厚度（van Rijn， 1984），将单宽输沙率q _b转换为推移质速度u _{bedload， i} ：

\frac{{\delta }_b}{d_i} = \mathrm{0.3} d_{\mathrm{*}}^{\mathrm{0.7}} {\left(\frac{{\theta }_i}{{{\theta }_{cr,i}}^{\mathrm{\text{'}}}}-\mathrm{1}\right)}^{\mathrm{0.5}},(11)

u_{\mathrm{b}\mathrm{e}\mathrm{d}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d},i} = \frac{q_b}{{\delta }_b{c}_{b,i}{f}_b},(12)

式中：δ _b为推移质层厚度；c_b，i 为不同坝料的体积分数；f _b为坝料的临界堆积分数，这里假定推移质速度与溃口底床界面附近流体流动方向相同.悬移质浓度可以通过对流-扩散方程来计算（Samma et al.， 2020）：

\frac{\partial C_{\mathrm{s},i}}{\partial t} + \nabla \cdot ( u_{s,i} C_{s,i} ) = \nabla \cdot \nabla ( \xi C_{s,i} )，(13)

式中：C_s，i 为不同坝料的悬沙质量浓度，可表示为流体-沉积物混合物的单位体积沉积物质量；ξ为方向扩散系数；u_s，i 为悬移质流速，u_s，i =u _m+ω_ic_s，i，其中u _m为流体-沉积物混合物的速度，c_s，i 为悬移质体积浓度，c_s，i =C_s，i /ρ_i .

为了避免数值计算过程中出现不稳定的情况，坝料冲蚀进行一个时间步长的限制，不允许在一个时间步长内跨多个计算单元输运.而基于欧拉模型的VOF法可以精确描述水和沉积物之间界面（Marsooli and Wu， 2015），通过流体体积与单位体积的比值来追踪自由表面流动，结合FAVOR（分数面积/体积比）方法可准确地解析复杂几何固体的数据（Liang et al.， 2019）.为求解控制方程，采用有限体积法将数值模拟离散成三维矩形网格，并将计算变量定义到网格边界的中心点.此外，由于具有收敛速度快、计算精度高等优点，广义最小残差GMRES算法被选择用于求解方程（Samma et al.， 2020）.

3 “11•03”白格堰塞湖溃决过程概述

2018年10月10日22∶06，中国西藏自治区昌都市江达县和四川省甘孜藏族自治州白玉县交界的波罗乡白格村发生山体第一次滑坡（位于东经98°42′17.98"，北纬31°4′56.41"，如图7所示），堵塞金沙江干流河道形成堰塞体.根据现场实测（Cai et al.， 2020），“10.10”白格堰塞体在未经人工干预的情况下于10月12日自然漫顶溃决.2018年11月3日17∶40，原山体滑坡点发生二次滑坡，滑坡形成的高速碎屑流再次堵塞河道形成堆积体（图8a，图8b）.本文主要针对拥有完整地质勘察资料的“11•03”白格堰塞体，采用地理信息系统（GIS）和数字高程模型（DEM）相结合的快速空间信息处理技术，获取了堆积体的形态特征（图8c）.“11•03”白格堰塞体滑坡区域高程范围从山顶的3 800 m到金沙江河床的2 870 m，堆积体地形起伏较大，垭口高程为2 966 m，最小坝高约为96 m（图9）.地质勘察资料表明（Zhang et al.， 2019），堰塞体的颗粒总体较细，主要物质为碎石土.结合堰塞体的高度、上游来水量以及水位上涨速度等监测信息，现场采取开挖泄流槽导流的方法进行除险，从而有效减少堰塞湖蓄水量.泄流槽开挖长度约220 m，最大顶宽42 m，底宽3 m，平均深度11.5 m，入口处底部高程为2 952.50 m（Cai et al.， 2020）.堰塞体形成时的上游水位为2 892.84 m，随后水位迅速抬升，于11月12日04∶45上升到泄流槽底坎高程（此时水位2 952.52 m），湖水开始进入泄流槽.11月13日13∶45，水位达到最高值2 956.40 m，堰塞湖容积达到最大值5.78×10⁸ m³，18∶00溃口达到峰值流量约31 000 m³/s；当天20∶00，溃口流量迅速下降至7 700 m³/s.11月14日8∶00，溃口流量与来水量持平，水位下降至2 905.75 m，堰塞体溃决过程结束.

4 反演分析

4.1　数值模型建立

基于“11•03”白格堰塞体溃决前的三维扫描地形数据（Cai et al.， 2020），通过Civil3D软件结合实际地貌等高线数据等比例建立三维实体STL（Stereolithography）模型.为了便于后续网格处理，截取求解域成规则长方形区域，尺寸为1 150 m×1 050 m×210 m（图10a）.采用4.5 m×4.5 m×4.5 m的网格对模型进行离散，共计278万个网格单元.考虑到实际溃决时间和计算效率，模拟时间设定为2018年11月12日04∶45至11月13日22∶00.

为了精确模拟堰塞体的实际溃决过程，模型将入口边界设置为水位压力边界，具体水位升降情况参照实际堰塞湖水位变化，出口边界设置为自由出流边界，模型两侧和底部设置为刚体无滑移边界，顶部设置为标准大气压边界.由于白格堰塞体整体颗粒较均匀，综合考虑坝料的颗粒级配曲线和地勘资料（Zhang et al.， 2019； Cai et al.， 2020），采用平均粒径d ₅₀作为堰塞体材料的计算粒径，模型的初始水位设置为2 952.52 m，与泄流槽底坎高程基本齐平，其余输入参数见表1.根据野外调查和现场勘察资料，可以估算出坝料平均粒径、密度以及堆积体的休止角，同时通过敏感性分析可以获取挟带率系数和推移质系数的最佳值（Kaurav and Mohapatra， 2019）.在数值模型中设置了4个监测点，对不同位置的溃决水流特征进行分析，其中点1和点2分别设置于泄流槽入口和转折处，点3设置在堰塞体下游斜坡处，点4设置于堰塞体下游坡脚处（图10b）.此外，在堰塞体泄流槽入口处设置了流量监测断面，以获取堰塞体溃决过程中溃口流量的变化（图10b）.

4.2　溃口流量过程

图11为“11•03”白格堰塞体溃口流量过程计算值与现场监测数据的对比，对比结果表明：模型计算得到的溃口峰值流量为29 518 m³/s，与实测值（31 000 m³/s）的相对误差为-4.78%；计算得到的溃口峰值流量出现时间为35.06 h，与实测值（溃决后37.25 h）相比，两者的相对误差为-5.88%.由此可得，精细化数值模型较好地重现了堰塞体的溃决过程，溃口流量过程与实测值基本吻合.

4.3　溃决过程特征

根据影像观测资料（Cai et al.， 2020），图12为现场观测与数值模拟溃决不同阶段的溃口形态对比图，可以看出三维数值模拟得出的“11•03”白格堰塞体溃决过程与实测过程基本一致.

为了进一步分析堰塞体溃决过程中各断面的水动力学特征，提取4个监测点的数值模拟结果，各监测点自由液面高程、水深和平均流速随时间的变化情况分别如图13所示.由于堰塞体溃决初期水动力条件较弱，各监测点数值差异较小，图13中的起始时间设定为2018年11月12日18∶00.从图中可以看出，在溃决第1阶段的表层冲蚀期，沿水流方向各监测点随时间增加开始有水流，但水深和平均流速维持在较低水平；随后，入口处（点1）的自由液面和水深上升最明显，点3和点4的平均流速在此阶段增长最快，验证了溯源冲蚀的特征.特别是计算得到的最大自由液面高程为2 958.60 m，略高于实测值2 956.40 m（蔡耀军等， 2019）.随着溃决进入第2阶段，受泄流槽持续纵向下切的影响，各监测点的自由液面高程迅速下降；在这一阶段产生了各个监测点的最大水深和平均流速，且沿水流方向各监测点的最大水深依次减小，并出现了图11对应的溃口峰值流量；点1的平均流速增加最为明显，这是由于溃口持续下切和侧向冲刷造成的，同时，其他监测点的最大平均流速与实测值10 m/s也非常接近（蔡耀军等， 2019）.随着自由液面高程进入稳定期，可以推测出溃口的冲蚀已经从纵向下切转变为横向展宽，此时溃决进入第3阶段，各监测点的平均流速基本保持一致，自由液面和流体深度也趋于统一.由此可以看出，本文模型获得的“11•03”白格堰塞体溃决过程水动力特征符合溃决过程各阶段的实际情况，并与测量值基本吻合.

4.4　溃口形态演化过程

基于无人机勘测技术，Cai et al. （2020）对“11•03”白格堰塞体进行了实地测绘，并对水下地形进行了测量，得到溃决后残余堰塞体的轮廓和详细的河道剖面信息，为溃口形态演化过程的研究提供了重要的实测数据.为了验证数值模拟获得的溃口最终断面的合理性，从已知断面（图4a）中选取泄流槽入口处（Ⅱ-Ⅱ'）、中间拐点处（Ⅴ-Ⅴ'）、下游扩展区（Ⅶ-Ⅶ'）和出流处（Ⅺ-Ⅺ'）4个典型断面与模型对应位置的断面（1-1'、2-2'、3-3'和4-4'）进行对比，断面在建立的数学模型中的具体位置如图14所示.

图15给出了4个典型断面计算和实测的最终溃口形态，通过对比可以看出，模型预测结果与实测数据基本吻合.其中，泄流槽转折处（断面2-2'）最终溃口形态与实测数据之间差距相对较大，这是由于此处水流情况较为复杂（图15b）.由于溃口河道两侧陡峻河谷的限制，下游坝坡断面3-3'比其他断面侧向冲蚀更弱（图15c）.此外，堰塞体下游坡脚处（断面4-4'）最终溃口高程因溃决期间水流的持续冲蚀而相对较低（图15d）.另外，实测和计算的顺河向堰塞体最终形态对比如图16所示，可以看出堰塞体溃决后残余坝体的下游坡比减小，这与模型试验和现场实测的溃决过程描述一致.

5 结论与展望

漫顶诱发堰塞体溃决过程往往受到许多因素影响，而在应急抢险中需要进行科学的预测才能有效地指导实践.为了深入探究溃口形态演化及关键溃决参数，基于各类试验揭示的堰塞体溃决机理和溃决过程，建立了三维精细化溃决过程数学模型.选取拥有实际地形数据和翔实地勘资料的“11•03”白格堰塞体，利用快速空间信息处理技术重建了实际三维地形模型，基于有限体积法的三维数值模拟结果与实测数据有良好的一致性.模拟结果与实测数据的对比表明，模型可较好地重现“11•03”白格堰塞体的溃口流量过程，与实测值相比，溃决流量峰值和峰值流量出现时间的相对误差均小于±10%.对比模型结果和现场拍摄实测视频，并提取溃口纵断面上4个监测点获取的水动力学特征值（自由液面高程、水深、平均流速），验证了数值模拟得出的溃口演化规律与实测堰塞体溃决过程表现出的“泄流槽表层和溯源冲蚀—沿程冲蚀—溃口稳定”等3个发展阶段较为吻合.此外，选取堰塞体溃决后残留坝体在泄流槽入口处、中间拐点处、下游坝坡和出流处等4个典型横断面，对比发现溃口各断面最终形态的计算值与实测值基本一致.上述研究成果验证了数学模型的合理性，表明该数值模拟方法可为科学制定堰塞体应急响应方案提供技术支撑.

由于堰塞湖溃决过程牵涉复杂的水土耦合效应，目前大多采用河道泥沙的输移公式描述溃口冲蚀过程，尚缺乏反映高速溃决水流作用下宽级配堰塞体材料输移特性的冲蚀公式；此外，精细化模型虽然可以有效提高数值模拟的精度，但是在堰塞湖溃决过程大尺度模拟时计算效率仍相对较低.未来应基于大比尺冲蚀试验，提出合理的宽级配坝料冲蚀公式，并在数值模拟计算速度提升、计算规模扩大等方面开展深入研究.

参考文献

[1] ASCE/EWRI Task Committee on Dam/Levee Breaching, 2011. Earthen Embankment Breaching. Journal of Hydraulic Engineering, 137(12): 1549-1564. https://doi.org/10.1061/(asce)hy.1943-7900.0000498

[2] Bagnold, R.A., 1966. An Approach to the Sediment Transport Problem from General Physics. U.S. Geological Survey Professional Paper, 422(1): 231-291.

[3] Cai, Y. J., Cheng, H. Y., Wu, S. F., et al., 2020. Breaches of the Baige Barrier Lake: Emergency Response and Dam Breach Flood.Science China Technological Sciences, 63(7): 1164-1176. https://doi.org/10.1007/s11431-019-1475-y

[4] Cai, Y.J., Luan, Y.S., Yang, Q.G., et al., 2019. Study on Structural Morphology and Dam-Break Characteristics of Baige Barrier Dam on Jinsha River. Yangtze River, 50(3): 15-22 (in Chinese with English abstract).

[5] Cao, P., Li, Y.S., Li, Z.L., et al., 2021. Geological Structure Characteristics and Genetic Mechanism of Baige Landslide Slope in Changdu, Tibet. Earth Science, 46(9): 3397-3409 (in Chinese with English abstract).

[6] Chang, D. S., Zhang, L. M., 2010. Simulation of the Erosion Process of Landslide Dams Due to Overtopping Considering Variations in Soil Erodibility along Depth. Natural Hazards and Earth System Sciences, 10(4): 933-946. https://doi.org/10.5194/nhess-10-933-2010

[7] Chen, C., Zhang, L. M., Xiao, T., et al., 2020. Barrier Lake Bursting and Flood Routing in the Yarlung Tsangpo Grand Canyon in October 2018. Journal of Hydrology, 583: 124603. https://doi.org/10.1016/j.jhydrol.2020.124603

[8] Chen, S.S., Chen, Z.Y., Zhong, Q.M., 2019. Progresses of Studies on Failure Mechanism and Numerical Dam Failure Model of Earth-Rockfill Dam and Landslide Dam. Water Resources and Hydropower Engineering, 50(8): 27-36 (in Chinese with English abstract).

[9] Chen, Z. Y., Ma, L. Q., Yu, S., et al., 2015. Back Analysis of the Draining Process of the Tangjiashan Barrier Lake. Journal of Hydraulic Engineering, 141(4): 05014011. https://doi.org/10.1061/(asce)hy.1943-7900.0000965

[10] Costa, J. E., Schuster, R. L., 1988. The Formation and Failure of Natural Dams. Geological Society of America Bulletin, 100(7): 1054-1068. https://doi.org/10.1130/0016-7606(1988)1001054: tfafon>2.3.co;2

[11] Fan, X. M., Scaringi, G., Korup, O., et al., 2019. Earthquake-Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts. Reviews of Geophysics, 57(2): 421-503. https://doi.org/10.1029/2018rg000626

[12] Fan, X. M., Zhan, W. W., Dong, X. J., et al., 2018. Analyzing Successive Landslide Dam Formation by Different Triggering Mechanisms: The Case of the Tangjiawan Landslide, Sichuan, China. Engineering Geology, 243: 128-144. https://doi.org/10.1016/j.enggeo.2018.06.016

[13] Guan, M. F., Wright, N. G., Andrew Sleigh, P., 2014. Multimode Morphodynamic Model for Sediment-Laden Flows and Geomorphic Impacts. Journal of Hydraulic Engineering, 141(6): 04015006. https://doi.org/10.1061/(asce)hy.1943-7900.0000997

[14] Jiang, X. G., Wei, Y. W., 2020. Erosion Characteristics of Outburst Floods on Channel Beds under the Conditions of Different Natural Dam Downstream Slope Angles. Landslides, 17(8): 1823-1834. https://doi.org/10.1007/s10346-020-01381-y

[15] Kaurav, R., Mohapatra, P. K., 2019. Studying the Peak Discharge through a Planar Dam Breach. Journal of Hydraulic Engineering, 145(6): 06019010. https://doi.org/10.1061/(asce)hy.1943-7900.0001613

[16] Li, K., Cheng, Q.G., Lin, Q.W., et al., 2022. State of the Art on Rock Avalanche Dynamics from Granular Flow Mechanics. Earth Science, 47(3): 893-912 （in Chinese with English abstract).

[17] Liang, C. F., Abbasi, S., Pourshahbaz, H., et al., 2019. Investigation of Flow, Erosion, and Sedimentation Pattern around Varied Groynes under Different Hydraulic and Geometric Conditions: A Numerical Study. Water, 11(2): 235. https://doi.org/10.3390/w11020235

[18] Luo, J., Pei, X. J., Evans, S. G., et al., 2019. Mechanics of the Earthquake-Induced Hongshiyan Landslide in the 2014 Mw 6.2 Ludian Earthquake, Yunnan, China. Engineering Geology, 251: 197-213. https://doi.org/10.1016/j.enggeo.2018.11.011

[19] Marsooli, R., Wu, W. M., 2015. Three-Dimensional Numerical Modeling of Dam-Break Flows with Sediment Transport over Movable Beds. Journal of Hydraulic Engineering, 141(1): 04014066. https://doi.org/10.1061/(asce)hy.1943-7900.0000947

[20] Mastbergen, D. R., van den Berg, J. H., 2003. Breaching in Fine Sands and the Generation of Sustained Turbidity Currents in Submarine Canyons. Sedimentology, 50(4): 625-637. https://doi.org/10.1046/j.1365-3091.2003.00554.x

[21] Mei, S. Y., Chen, S. S., Zhong, Q. M., et al., 2021. Effects of Grain Size Distribution on Landslide Dam Breaching—Insights from Recent Cases in China. Frontiers in Earth Science, 9: 658578. https://doi.org/10.3389/feart.2021.658578

[22] Meyer-Peter, E., Muller, R., 1948. Formulas for Bed-Load Transport. Process of Congress IAHR, 6(2): 39-64.

[23] Movahedi, A., Kavianpour, M. R., Yamini, O. A., 2018. Evaluation and Modeling Scouring and Sedimentation around Downstream of Large Dams. Environmental Earth Sciences, 77(8): 320. https://doi.org/10.1007/s12665-018-7487-2

[24] Peng, M., Zhang, L. M., 2012. Breaching Parameters of Landslide Dams. Landslides, 9(1): 13-31. https://doi.org/10.1007/s10346-011-0271-y

[25] Qian, N., 1980. A Comparison of the Bed Load Formulas. Journal of Hydraulic Engineering, 11(4): 1-11 (in Chinese with English abstract).

[26] Roseberry, J. C., Schmeeckle, M. W., Furbish, D. J., 2012. A Probabilistic Description of the Bed Load Sediment Flux: 2. Particle Activity and Motions. Journal of Geophysical Research: Earth Surface, 117(F3): F03032. https://doi.org/10.1029/2012jf002353

[27] Samma, H., Khosrojerdi, A., Rostam-Abadi, M., et al., 2020. Numerical Simulation of Scour and Flow Field over Movable Bed Induced by a Submerged Wall Jet. Journal of Hydroinformatics, 22(2): 385-401. https://doi.org/10.2166/hydro.2020.091

[28] Shi, Z.M., Ma, X.L., Peng, M., et al., 2014. Statistical Analysis and Efficient Dam Burst Modelling of Landslide Dams Based on a Large-Scale Database. Chinese Journal of Rock Mechanics and Engineering, 33(9): 1780-1790 (in Chinese with English abstract).

[29] van Rijn, L. C., 1984. Sediment Transport, Part I: Bed Load Transport. Journal of Hydraulic Engineering, 110(10): 1431-1456. https://doi.org/10.1061/(asce)0733-9429(1984)110: 10(1431)

[30] van Rijn, L. C., 2020. Erodibility of Mud-Sand Bed Mixtures. Journal of Hydraulic Engineering, 146(1): 04019050. https://doi.org/10.1061/(asce)hy.1943-7900.0001677

[31] Walder, J. S., Iverson, R. M., Godt, J. W., et al., 2015. Controls on the Breach Geometry and Flood Hydrograph during Overtopping of Noncohesive Earthen Dams. Water Resources Research, 51(8): 6701-6724. https://doi.org/10.1002/2014wr016620

[32] Yakhot, V., Orszag, S. A., Thangam, S., et al., 1992. Development of Turbulence Models for Shear Flows by a Double Expansion Technique. Physics of Fluids A: Fluid Dynamics, 4(7): 1510-1520. https://doi.org/10.1063/1.858424

[33] Zhang, J. Y., Fan, G., Li, H. B., et al., 2021. Large-Scale Field Model Tests of Landslide Dam Breaching. Engineering Geology, 293: 106322. https://doi.org/10.1016/j.enggeo.2021.106322

[34] Zhang, L. M., Xiao, T., He, J., et al., 2019. Erosion-Based Analysis of Breaching of Baige Landslide Dams on the Jinsha River, China, in 2018. Landslides, 16(10): 1965-1979. https://doi.org/10.1007/s10346-019-01247-y

[35] Zhao, T.L., Chen, S.S., Wang, J.J., et al., 2016. Centrifugal Model Tests Overtopping Failure of Barrier Dams. Chinese Journal of Geotechnical Engineering, 38(11): 1965-1972 (in Chinese with English abstract).

[36] Zhong, Q. M., Chen, S. S., Wang, L., et al., 2020. Back Analysis of Breaching Process of Baige Landslide Dam. Landslides, 17(7): 1681-1692. https://doi.org/10.1007/s10346-020-01398-3

[37] Zhong, Q. M., Wang, L., Chen, S. S., et al., 2021. Breaches of Embankment and Landslide Dams-State of the Art Review. Earth-Science Reviews, 216: 103597. https://doi.org/10.1016/j.earscirev.2021.103597

[38] Zhu, X. H., Liu, B. X., Peng, J. B., et al., 2021. Experimental Study on the Longitudinal Evolution of the Overtopping Breaching of Noncohesive Landslide Dams. Engineering Geology, 288: 106137. https://doi.org/10.1016/j.enggeo.2021.106137

[39] 蔡耀军, 栾约生, 杨启贵, 等, 2019. 金沙江白格堰塞体结构形态与溃决特征研究. 人民长江, 50(3): 15-22.

[40] 曹鹏, 黎应书, 李宗亮, 等, 2021. 西藏昌都白格滑坡斜坡地质结构特征及成因机制. 地球科学, 46(9): 3397-3409.

[41] 陈生水, 陈祖煜, 钟启明, 2019. 土石坝和堰塞坝溃决机理与溃坝数学模型研究进展. 水利水电技术, 50(8): 27-36.

[42] 李坤, 程谦恭, 林棋文, 等, 2022. 高速远程滑坡颗粒流研究进展. 地球科学, 47(3): 893-912.

[43] 钱宁, 1980. 推移质公式的比较. 水利学报, 11(4): 1-11.

[44] 石振明, 马小龙, 彭铭, 等, 2014. 基于大型数据库的堰塞坝特征统计分析与溃决参数快速评估模型. 岩石力学与工程学报, 33(9): 1780-1790.

[45] 赵天龙, 陈生水, 王俊杰, 等, 2016. 堰塞坝漫顶溃坝离心模型试验研究. 岩土工程学报, 38(11): 1965-1972.

基金资助

国家自然科学基金联合基金重点项目(U2040221)

国家重点研发计划课题(2018YFC1508604)

中央级公益性科研院所基本科研业务费专项资金(Y721006)