Simulation of landslide dam failure due to overtopping considering vertical soil particles distribution

Abstract

Accurate prediction of potential landslide dam breaches is crucial for timely and effective emergency responses to the resulting flood. Although many models for landslide dam overtopping have been proposed in past studies, most of these models did not consider the heterogeneity of dam materials, which is important in reflecting the erodibility of soil. In this study, a new mathematical model for landslide dam breaches due to overtopping was proposed based on a constructed vertical distribution model of dam particle size. Two real-world dams, the Tangjiashan landslide dam and Baige landslide dam, were used to evaluate the performance of the proposed model. In addition, the proposed model was compared with four existing dam failure models, and sensitivity analysis was conducted on the median particle size and particle size distribution. The results demonstrate that the computational outcomes of the proposed model align more closely with the actual results, indicating that the proposed model outperforms the four existing models in accurately simulating the landslide dam failure process. It also suggests that ignoring the variation of dam material particle size along the vertical direction inside the dam body would increase the prediction relative errors of key breach outcomes.

Keywords:

• Landslide dams
• overtopping failure
• particles distribution
• simulation

1. Introduction

Landslide dams are deposits of soil and rock that block valleys and rivers due to landslides caused by earthquakes, heavy precipitation and volcanic eruptions (Yan et al. 2009; Wu et al. 2024). There are great differences between landslide dams and artificial dams. For instance, the shape and height of landslide dams are irregular, and the accumulation range along the river is wider; the material gradation of the dam is wide and uneven and the internal structure within the dam is often loose (Zhong et al. 2018b). In addition, unlike artificial dams, landslide dams lack specific spillway arrangements. All these factors make landslide dams more prone to failure compared to artificial dams (Dong et al. 2009; Peng & Zhang 2012; Shen et al. 2020). Once a landslide dam breaks, floodwater will cause great damage in downstream areas. Therefore, it is important to estimate dam break flood of landslide dams accurately using appropriate mathematical models.

To forecast the parameters associated with landslide dam failures, including discharge, breach size and barrier lake water level, many researchers have studied the simulation technology of landslide dam breaks. Over the past few decades, a variety of parameter models have been developed as a result of these efforts (Zhong et al. 2018b; Chen et al. 2019). Costa (1985) proposed three regression equations to predict the peak discharge of the breach based on 10 cases of barrier dam break. Based on a Japanese dataset, Dong et al. (2011) utilized the logistic regression method to identify the important geomorphic variables and evaluated the failure probability of a landslide dam. Some scholars (Costa & Schuster 1988; Walder & O Connor1997) proposed several parameter models to predict peak flow. However, these early models only considered the hydrodynamic conditions of the barrier lake (Cencetti & Di Matteo 2022) and can only roughly predict the flow of the breach, while ignoring some important breach parameters such as breach size and reservoir water level. Recently, with increase in dam break data collected, the dam break parameter model has been developed from the single parameter models to multi-parameter models (Froehlich 1995; Xu & Zhang 2009; Froehlich 2016; Mei et al. 2018; Wu et al. 2023). The outputs of these models have expanded beyond peak discharge to encompass variables such as the final average width of the breach, duration of the dam break and others. Moreover, these models now incorporate factors such as dam shape, reservoir capacity and material characteristics of the dam (Duan et al. 2021).

Although there has been some significant progress in the simulation methods of landslide dam break, the current landslide dam failure models often neglect the significant heterogeneity characteristics of landslide dam materials. Typically, breaching outcomes of landslide dam are predicted assuming that parameters of the dam materials are constants (Zhong et al. 2020a). However, the formation of landslide dams inherently entails uneven spatial distribution of material composition within the dam body, with material composition varying significantly along the direction of gravity (Chang et al. 2011), resulting in changes in breaching parameters with depth, such as particle size, and their associated factors like incipient velocity and erodibility of soil, etc. Hence, it is difficult to accurately estimate the development of dam failure without considering these changes in parameters.

Therefore, in this study, a new landslide dam break model that incorporates a vertical distribution model of landslide dam material particle size is proposed to simulate the breaching process of landslide dams (Section 2), and the numerical computation procedure of the model is also established (Section 3). This model is evaluated using two case study dams, Tangjiashan landslide dam and Baige landslide dam (Sections 4 and 5), and against four commonly used existing models (Section 5). Finally, limitations of this model and future research are also discussed (Section 6).

2. Model development

For the calculation of landslide dam failure, previous physically based breach models usually include three parts: calculation of outflow, calculation of soil erosion and calculation of breach development (Chen et al. 2019). In order to consider the variation of dam material particle size along the depth, which has a significant impact on soil properties, a vertical distribution model of particle size is constructed as a new module of the breach mathematical model. Therefore, the new breach model is mainly composed of four modules. In addition, the sediment erosion rate and the associated equation for continuous breach expansion have also been re-derived. Each of these four modules is discussed as follows.

2.1. Module 1 – vertical distribution model of dam particles

An accurate mathematical description of the structure of a landslide dam can provide a reasonable initial input for the calculation of its failure process. In the process of landslide accumulation, especially for remote landslides, the sorting of soil and rock particles with different sizes is mainly reflected in the mixing and separation movement in the vertical direction. On the other hand, the autocorrelation distance of particle sizes in the horizontal direction is generally greater than that in the vertical direction, and it has a rough horizontal stratification feature (Chen et al. 2021). Therefore, this study considers landslide dams as horizontally isotropic heterogeneous materials and constructs a spatial particle size distribution model based on the vertical structure of reverse grading distribution.

Casagli et al. (2003) and Bu et al. (2017) discovered that the particle size and frequency of deposits resulting from geological disasters such as landslides and debris flows conform to a three-parameter Weibull distribution. The authors selected 12 landslide dams with clear grading data for curve fitting analysis, as depicted in Figure 1, and arrived at the same conclusion. In addition, the particle size of the landslide dam has self-similarity (Li et al. 2021), that is, for each soil layer of the landslide dam, the frequency distribution of the particle group is similar to the overall particle group frequency of the landslide dam. Therefore, a probability function of particle size distribution for landslide dams can be constructed as follows: (1) $$F\left(d\right)=\alpha +\gamma \left\{1-\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left[-{\left(\frac{d}{\lambda }\right)}^{\beta }\right]\right\}$$(1)

Figure 1. Analysis of particle grading curve for landslide dams.

Where d = particle size of landslide dam (mm), F (d) = accumulated content of soil particles smaller than particle size d, α = positional parameter, γ = scale parameter of particle size, λ = proportional parameter, β = shape parameter and α, γ, λ, β can be determined through regression analysis. According to Equation (1)(1) $$F\left(d\right)=\alpha +\gamma \left\{1-\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left[-{\left(\frac{d}{\lambda }\right)}^{\beta }\right]\right\}$$(1) , the calculation equation for any characteristic particle size (such as d₅₀) near the surface of the landslide dam crest can be derived as: (2) $$d_{\text{sur}}^{\text{*}}=\lambda \cdot {\left(\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\frac{\gamma }{\gamma -F(d)+\alpha }\right)}^{\frac{1}{\beta }}$$(2)

Where $d_{\text{sur}}^{\text{*}}$ = near surface characteristic particle size (mm). Then, a vertical distribution model for any characteristic particle size at any depth of the landslide dam can be constructed as: (3) $$d^{*}=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left[\frac{\left(H-H\prime \right)\cdot \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\left(d_{\text{sur}}^{\text{*}}-d_{\mathrm{\text{min}}}+1\right)}{H}\right]+d_{\mathrm{\text{min}}}-1$$(3)

Where $d^{*}$ = the characteristic particle size of a certain depth of soil layer (mm), H = height of landslide dam (m), H' = depth of breach (m), d_min = the minimum particle size of the landslide dam (mm). If the landslide dam is a clay dam, the minimum particle size is taken as the boundary particle size of the clay particles. In practical engineering, the surface gradation distribution is initially obtained through onsite investigation, such as the use of drones and image processing. Subsequently, the probability distribution of particle sizes is constructed, and the characteristic particle size value is determined. Then, a vertical distribution model of the characteristic particle size is established. The process for obtaining the spatial distribution of characteristic parameters of the dam body is shown in Figure 2.

Figure 2. Procedure for obtaining spatial particle size.

2.2. Module 2 – calculation of discharge

When a landslide dam overtops, the reservoir water flows through the lowest part of the dam crest and forms an initial breach. The development of the breach is closely related to the breach flow, whose development process also depends on the reservoir water and the erosion resistance of dam materials. Therefore, it is of great importance to accurately calculate the discharge at the dam break. During the process of landslide dams overtopping, the breach flow meets the water balance equation: (4) $$\frac{\mathit{\text{dV}}}{\mathit{\text{dt}}}=A_s\frac{d{z}_s}{\mathit{\text{dt}}}$$(4) (5) $$A_s\frac{d{z}_s}{\mathit{\text{dt}}}=Q_i-Q_b$$(5)

Where V = volume of barrier lake (m³), A_z = surface area of barrier lake (m²), z_s = water lever of barrier lake (m), t = breaching time (s), Q_i = inflow discharge (m³/s), Q_b = breach discharge (m³/s).

The initial breach of the dam crest is usually assumed to be trapezoidal or rectangular. According to the previous centrifugal model tests of landslide failure (Zhao et al. 2018; 2019), most of the breaches are developed with trapezoidal cross-sections. Therefore, in this study, the initial breach is assumed to be trapezoidal, and the breach flow could be calculated using the broad-crested weir equation: (6) $$Q_b=k_s{c}_d\sqrt{2g}\left[c_{1}b\Delta h^{1.5}+c_{2}m\Delta h^{2.5}\right]$$(6) (7) $$\Delta h=z_s-z_b$$(7)

Where Q_b = breach flow (m³/s), k_s = correction coefficient of submerged flow downstream, c_d = discharge coefficient, b = bottom width of the breach (m), $\Delta h$ = water depth of breach (m), z_b = elevation of breach bottom (m), g = gravitational acceleration (m/s²), m = breach slope coefficient, c₁ and c₂ are the discharge coefficients of rectangular and triangular parts of the breach respectively, and $c_1=1.7,$ $c_2=1.3$ according to the Singh (1996).

2.3. Module 3 – calculation of sediment transport

After the overtopping of landslide dam, the discharge flow usually erodes the breach vertically and horizontally, which makes the breach deepen and widen. In this study, the total sediment transport rate is obtained by calculating the bed load and suspended load sediment transport rates, respectively, based on the Bagnold Formula (1956). For bed load sediment transport rate: (8) $$\frac{{\rho }_s-\rho }{{\rho }_s}{g}_b=W_b{u}_b$$(8)

Where g_b = sediment transport rate of bed load (kg/m·s), W_b = weight of bed load sediment in unit area (kg/m²), u_b = bed load velocity (m/s), $\rho$ = water density (kg/m³), ${\rho }_s$ = sediment density (kg/m³).

In terms of energy, the energy consumed in carrying sediment particles in unit time is equal to the water flow power multiplied by its efficiency (Wu & Wang 2007), written as: (9) $$W_b{u}_b\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\alpha ={\tau }_{0}U{e}_b$$(9)

Where $W_b{u}_b\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\alpha$ = the work done by force in the direction of motion to transport sediment in unit time, α = internal friction angle of sediment (°), ${\tau }_0$ = drag force of bed surface (Pa), and ${\tau }_0\text{=}\rho \mathit{\text{gRJ}},$ with R = hydraulic radius (m), J = hydraulic gradient, and $J=\frac{\Delta h}{L+S},$ with L = breach length of dam crest (m), S = breach length of the downstream slope (m), ${\tau }_{0}U$ = potential energy loss of water flow in unit time, U = vertical average flow velocity (m/s), e_b = efficiency of bed load transport.

The sediment transport rate of the single-width bed load can be obtained by simultaneous Equations (5)(5) $$A_s\frac{d{z}_s}{\mathit{\text{dt}}}=Q_i-Q_b$$(5) and (6)(6) $$Q_b=k_s{c}_d\sqrt{2g}\left[c_{1}b\Delta h^{1.5}+c_{2}m\Delta h^{2.5}\right]$$(6) : (10) $$g_b=\frac{{\rho }_s}{{\rho }_s-\rho }\cdot \frac{e_b{\tau }_{0}U}{\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\alpha }$$(10)

Similarly, the sediment transport rate of suspended load can be written as (Wu et al. 2004, 2012): (11) $$\frac{{\rho }_s-\rho }{{\rho }_s}{g}_s=W_s\overline{u_s}$$(11)

Where g_s = sediment transport rate of suspended load (kg/m·s), $\overline{u_s}$= vertical average velocity in the flow direction (m/s).

As a whole, the centre of gravity of suspended sediment remains unchanged, that is, the sinking speed of suspended sediment is equal to the suspension one (Wu 2013). Then, the work done by suspended sediment is: (12) $$E_1=W_s\omega$$(12)

Where E₁ = work done by the suspended sediment (J), $\omega$ = settling velocity of sediment (m/s), which is determined by: (13) $$\omega =1.044\sqrt{\frac{{\gamma }_s-\gamma }{\gamma }g{d}_{50}}$$(13)

Where ${\gamma }_s$ = unit weight of soil (N/m³), $\gamma$ = unit weight of water (N/m³), g = gravitational acceleration (N/kg), $d_{50}$ = the average particle size of soil particles (mm). Based on the law of conservation of energy, the energy of water flow maintaining the movement of suspended load in unit time and bed surface is the potential energy of water flow loss minus that consumed by bed load movement and then multiplied by the efficiency coefficient. The energy of water flow maintaining the movement of suspended load is determined by: (14) $$E_2={\tau }_{0}U\left(1-e_b\right)e_s$$(14) Where E₂ = the energy of water flow maintaining the movement of the suspended load (J), e_s = efficiency of suspended load transport.

This part of the energy is used to maintain sediment suspension, that is, E₁ = E₂ (Wu & Li 2017). So, it has: (15) $$W_s\omega \text{=}{\tau }_{0}U\left(1-e_b\right)e_s$$(15)

According to Equations (11)(11) $$\frac{{\rho }_s-\rho }{{\rho }_s}{g}_s=W_s\overline{u_s}$$(11) , (12)(12) $$E_1=W_s\omega$$(12) and (15)(15) $$W_s\omega \text{=}{\tau }_{0}U\left(1-e_b\right)e_s$$(15) , there are: (16) $$g_s=\frac{{\rho }_s}{{\rho }_s-\rho }\frac{{\tau }_0\left(1-e_b\right)e_{s}U\overline{u_s}}{\omega }$$(16)

The vertical average velocity of the suspended load is directly proportional to the vertical average velocity of water flow. So, it can be written as: (17) $$\overline{u_s}\text{=}\mathit{\text{kU}}$$(17)

Where k = scale factor. Through the analysis of a large number of flume tests, Bagnold (1966) obtained $(1-e_b)e_{s}k=0.01.$ So, it has: (18) $$g_s=0.01\frac{{\gamma }_s}{{\gamma }_s-\gamma }\frac{{\tau }_0{U}^2}{g\omega }$$(18)

According to Einstein (1950), the combination of the suspended load and bed load sediment transport rates is known as the total sediment transport rate. Considering the high velocity and powerful scouring of a landslide dam breaching, a correction coefficient $\zeta$ (the value is 0.056 in this study) is introduced based on the total sediment transport rate, and the total sediment transport rate g_T may be calculated as: (19) $$g_T=\zeta \frac{{\gamma }_s}{{\gamma }_s-\gamma }\frac{{\tau }_{0}U}{g}\left(\frac{e_b}{\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\alpha }+0.01\frac{U}{\omega }\right)$$(19)

The total amount of sediment transfer is determined by: (20) $$Q_T=t{g}_T$$(20)

Where Q_T = total amount of sediment transfer (kg/m).

2.4. Module 4 – calculation of breach expansion

2.4.1. Continuous expansion of the breach

The main expansion types of the breach are surface erosion and traceable erosion. Due to the low contents of clay, the erosion mode of landslide dam breach is mainly surface erosion (Zhong et al. 2021), as shown in Figure 3.

Figure 3. Vertical and horizontal continuous expansion of the breach.

The increment of longitudinal cut depth of breach $\Delta z$ can be calculated by: (21) $$\Delta z=\frac{\Delta t{g}_T}{{\rho }_{s}l\left(1-n_0\right)}$$(21)

Where n₀ = porosity of dam material, l = generalized breach length (m), $\Delta z$ = time step (s). According to previous centrifugal model tests of landslide dam breaching (Zhao et al. 2018), the horizontal and vertical continuous development process of the breach roughly follows the following geometric relationship AB = BD. Therefore, based on geometric relationships in Figure 3, the increment of breach top width expansion $\Delta B$ and the bottom width expansion $\Delta b$ can be expressed as: (22) $$\{\begin{array}{c}\Delta B=\frac{2\Delta z}{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\beta }\\ \Delta b=2\Delta z\left(\frac{1-\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\beta }{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\beta }\right)\end{array}$$(22)

Where $\beta$ = breach slope angle (°), ΔB = the increment of the expansion of the breach top width (m), Δb = the increment of the expansion of the breach bottom width (m).

2.4.2. Collapse extension of the breach

According to Figure 3, with the continuous erosion of the breach, the angle of the breach slope will increase, resulting in the collapse of the breach slope at a certain angle. In this research, it is assumed that the sliding surface of the slope is a plane passing through the toe, that is, the collapse of the breach only causes the expansion of the top width (Figure 4). Therefore, it is important to determine the dip angle of the sliding surface for calculating the expansion increment of the breach top width.

Figure 4. Stability calculation of breach slope.

Assuming that the slope will collapse when the slope of the breach reaches the critical one, the mechanical state of the sliding surface can be expressed as: (23) $$W\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_k=\mathit{\text{cL}}+W\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}{\beta }_k\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\phi$$(23)

Where W = weight of the sliding body (kN), c = cohesion of soil mass (kPa), L = length of the sliding surface (m), $\phi$ = internal friction angle of soil (°), ${\beta }_k$ = dip angle of the sliding surface (°). The weight of the sliding body can be expressed as: (24) $$W=0.5{\gamma }_s\frac{H\prime }{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}}L\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_{k-1}-{\beta }_k\right)$$(24)

Where ${\beta }_{k-1}$ = critical dip angle of breach slope (°), H′ = breach depth (m). Combining Equations (23)(23) $$W\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_k=\mathit{\text{cL}}+W\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}{\beta }_k\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\phi$$(23) and (24)(24) $$W=0.5{\gamma }_s\frac{H\prime }{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}}L\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_{k-1}-{\beta }_k\right)$$(24) , it can be obtained: (25) $$c=0.5{\gamma }_{s}H\prime \frac{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_{k-1}-{\beta }_k\right)\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_k-\phi \right)}{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }$$(25)

According to Mohr-Coulomb Strength Theory, when the soil is destroyed in the limit equilibrium state, it meets ${\beta }_k\text{=}0.5({\beta }_{k-1}+\phi ).$ Combined with Equation (25)(25) $$c=0.5{\gamma }_{s}H\prime \frac{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_{k-1}-{\beta }_k\right)\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\left({\beta }_k-\phi \right)}{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }$$(25) , the following equation can be obtained. (26) $$H\prime =\frac{4c\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }{{\gamma }_s\left[1-\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\left({\beta }_{k-1}-\phi \right)\right]},\mathit{\text{k}}=1,2,3,\dots$$(26)

Whether the slope will collapse can be judged by the dip angle of the breach slope. According to Equation (26)(26) $$H\prime =\frac{4c\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}{\beta }_{k-1}\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }{{\gamma }_s\left[1-\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\left({\beta }_{k-1}-\phi \right)\right]},\mathit{\text{k}}=1,2,3,\dots$$(26) , the critical dip angle of breach slope can be expressed as: (27) $${\beta }_{k-1}=\mathrm{\text{arcsin}}\sqrt{1-\frac{{\left[2\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi -2\left(\frac{4c\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }{H\prime {\gamma }_s}+\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\phi \right)\sqrt{{\left(\frac{4c\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }{H\prime {\gamma }_s}+\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\phi \right)}^2-{\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}}^2\phi }\right]}^2}{4{\left[{\left(\frac{4c\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\phi }{H\prime {\gamma }_s}+\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\phi \right)}^2+{\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}}^2\phi \right]}^2}}$$(27)

When the breach slope collapses, the increment of the bottom width and depth of the breach is 0 according to Figure 4, and the increment of the top width can be calculated according to the following equation: (28) $$\Delta B\prime =H\prime \left(\frac{1}{\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}{\beta }_k}-\frac{1}{\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}{\beta }_{k-1}}\right)$$(28)

3. Numerical calculation process

The iterative method based on time step is used for running model. By inputting the initial conditions, boundary conditions and parameters of barrier dam and lake, the model can be used to calculate the development process of dam break due to overtopping. The calculation process is as follows.

• Step 1: input initial parameters, including dam geometry, crest breach, dam material, overtopping flow, reservoir capacity and downstream river.

• Step 2: set $\Delta t.$

• Step 3: make $t+\Delta t=t.$

• Step 4: construct the probability function of particle size distribution based on Equation (1)(1) $$F\left(d\right)=\alpha +\gamma \left\{1-\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left[-{\left(\frac{d}{\lambda }\right)}^{\beta }\right]\right\}$$(1) , calculate surface characteristic particle size $d_{\text{sur}}^{\text{*}}.$

• Step 5: construct vertical particle size distribution model based on Equation (2)(2) $$d_{\text{sur}}^{\text{*}}=\lambda \cdot {\left(\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\frac{\gamma }{\gamma -F(d)+\alpha }\right)}^{\frac{1}{\beta }}$$(2) , calculate characteristic particle size at a specific depth $d^{*}.$

• Step 6: calculate breach flow Q_b based on Equation (6)(6) $$Q_b=k_s{c}_d\sqrt{2g}\left[c_{1}b\Delta h^{1.5}+c_{2}m\Delta h^{2.5}\right]$$(6) .

• Step 7: calculate reservoir water level increment $\Delta z_s$ based on Equation (5)(5) $$A_s\frac{d{z}_s}{\mathit{\text{dt}}}=Q_i-Q_b$$(5) .

• Step 8: calculate total sediment transport rate g_T based on Equation (19)(19) $$g_T=\zeta \frac{{\gamma }_s}{{\gamma }_s-\gamma }\frac{{\tau }_{0}U}{g}\left(\frac{e_b}{\hspace{0.17em}\mathrm{\text{tan}}\hspace{0.17em}\alpha }+0.01\frac{U}{\omega }\right)$$(19) .

• Step 9: calculate the increment of longitudinal cut depth $\Delta z$ and the width expansion $\Delta B,$ $\Delta b$ of breach based on Equations (21)(21) $$\Delta z=\frac{\Delta t{g}_T}{{\rho }_{s}l\left(1-n_0\right)}$$(21) and (22). And calculate breach slope $\beta .$

• Step 10: judge whether $\beta <{\beta }_{k-1}.$ if yes, go to Step 7; if no, go to step 9.

• Step 11: make $B=B+\Delta B,$ $b=b+\Delta b,$ $z_b=z_b-\Delta z,$ $z_s=z_s-\Delta z_s.$

• Step 12: calculate the breach top width B, bottom width b, elevation of bottom z_b and reservoir water lever z_s according to step 11.

• Step 13: judge whether $z_s\le z_b.$ if yes, go to step 14; if no, go to step 5.

• Step 14: output Q, B, b, z_b and z_s.

The calculation process is shown in Figure 5.

Figure 5. The model calculation process of the mathematical model.

4. Case study and data

4.1. Tangjiashan landslide dam

Tangjiashan landslide dam was created during the 2008 Wenchuan earthquake in China’s Southwest Sichuan Province, with a height of approximately 100 m. It has a catchment area of 3550 km², a storage capacity of 145 million m³ and a volume of about 316 million m³ after full storage. The Tongkou River in the dam site is an asymmetric V-shaped valley, with a steep right bank slope of about 60° and a gentle left bank slope of about 30°. The geology of the site consists of Cambrian siliceous sandstone and phyllites with interbedded black slates and mylonites. During the earthquake, the hillslope on the right side of the river collapsed along the bedding planes that strike at N60°E and dip Northwest at an angle of 60° (Liu et al. 2010), as shown in Figure 6.

Figure 6. The (a) full extent, (b) cross-section, and (c) longitudinal section of the Tangjiashan landslide dam (Hu et al. 2010).

Due to the large dam height and storage capacity, once the Tangjiashan landslide dam breaks, it will greatly endanger the downstream people. To eliminate the danger as soon as possible, the earthquake relief headquarters of the Ministry of Water Resources, PRC decided to excavate the slugging channel for Tangjiashan landslide dam. As the similarity between the discharge process of the landslide dam and the overtopping, and the availability of detailed monitoring data, this study chose the Tangjiashan landslide dam discharge to validate the mathematical model.

According to the measured grading of the Tangjiashan landslide dam (Figure 7), the probability function of particle size distribution in the dam body can be determined. The particle size required for this model calculation is mainly the median particle size. Therefore, $F(d_{\text{sur,50}})=0.5$ can be incorporated into the probability distribution function and d₅₀ = 22.9 mm can be obtained. According to onsite investigation, the height of the Tangjiashan landslide dam is 103 m, and the minimum particle size d_min = 0.075 mm. Therefore, a vertical distribution model of d₅₀ can be constructed as: (29) $$d_{50}=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(3.17-0.0308H\prime \right)-0.925$$(29)

Figure 7. Grading analysis of Tangjiashan landslide dam.

The calculation parameters are shown in Table 1.

Table 1. Conditions of Tangjiashan landslide dam.

Index	Value	Index	Value
Dam height (m)	103	Initial breach depth (m)	13
Crest width (m)	300	Initial breach width (m)	8
Dam length (m)	612	ρs (kg/m^3)	2216
U/S slope (V/H)	0.36	d50 (m)	Equation (29)
D/S slope (V/H)	0.24	c (kPa)	25^a
Reservoir storage (m^3)	1.5 × 10^8	φ (º)	22
Initial reservoir lever (m)	92.5	n0	40%
Initial downstream water lever (m)	0^a	Time step (s)	1

^a Assumed values.

4.2. Baige landslide dam

At 22:06 on 10 October 2018, a large-scale landslide occurred in Baige Village (31°4′56.41′′N, 98°42′17.98′′E) at the junction of Sichuan Province, China and the Tibet Autonomous Region. The landslide blocked the mainstream of Jinsha River, forming a barrier lake, which naturally discharged on 12 October (Cai et al. 2020). On 3 November 2018, a second landslide occurred at the same location (Fan et al. 2019), with a volume of approximately 160 × 10⁴ m³ blocking the flow channel formed by the “10.10” landslide, causing a second river blockage.

The strata in the landslide area consist of the Permian Early Triassic Gangtuoyan Formation (PT₁g), distributed in a nearly north-south direction in the Jinsha River junction zone. The lithology predominantly comprises gray sericite quartz schist, with thin layers of quartz interbedded in the middle and lower parts. Certain local areas contain notable amounts of chlorite, charcoal and calcium. A regional reverse fault is developed at the rear of the landslide body (Huang et al. 2021). According to the exploration data, Figure 8 provides a schematic diagram of the cross-sectional morphology of the Baige landslide dam formed by these two landslides distinguished by two colours (Zhang et al. 2019).

Figure 8. The (a) cross-section, (b) longitudinal section of Baige landslide dam (Zhong et al. 2020b).

According to the height of the dam, the mode of collapse and the speed of upstream water level rise, the Yangtze River Defense Administration formulated a manual intervention plan on 4 November. The designed drainage groove is arranged on the left side of the “11.3” accumulation body, with a position slightly to the left compared to the “10.10” flow channel. The construction of the drainage channel began on 8 November and ended at 11:00 on the 11th. The excavation length was 220 m, with a bottom width of 3 m. The average depth was 13.48 m, and the maximum excavation depth was nearly 15 m. After the breaching of the barrier lake on 13 November, a new river channel was formed based on the drainage channel. This study selected the “11.3” Baige landslide dam (the second landslide) to further validate the breach model proposed.

According to the measured grading of the Tangjiashan landslide dam (Figure 9), the probability function of particle size distribution in the dam body can be determined. According to the probability distribution function, d₅₀ = 4.88 mm can be obtained. The height of Baige landslide dam is 96 m, and the minimum particle size d_min = 0.075 mm. Therefore, a vertical distribution model of d₅₀ can be constructed as: (30) $$d_{50}=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(1.759-0.018H\prime \right)-0.925$$(30)

Figure 9. Grading analysis of Baige landslide dam.

The calculation parameters are shown in Table 2.

Table 2. Conditions of Baige landslide dam.

Index	Value	Index	Value
Dam height (m)	96	Initial breach depth (m)	13.48
Crest width (m)	270	Initial breach width (m)	3
Dam length (m)	600	ρs (kg/m^3)	1854
U/S slope (V/H)	0.37	d50 (m)	Equation (30)
D/S slope (V/H)	0.18	c (kPa)	3
Reservoir storage (m^3)	7.57 × 10^8	φ (º)	32.8
Initial reservoir lever (m)	85.53	n0	37.5%
Initial downstream water lever (m)	0^a	Time step (s)	1

^a Assumed values.

4.3. Discharge coefficients

Although Singh (1996) provided an empirical value for the discharge coefficient in Equation (6)(6) $$Q_b=k_s{c}_d\sqrt{2g}\left[c_{1}b\Delta h^{1.5}+c_{2}m\Delta h^{2.5}\right]$$(6) , to determine its applicability in the two real cases discussed in this article, we conducted sensitivity analysis on the two coefficients using data of initial discharge from the Tangjiashan and Baige landslide dams. For the Tangjiashan landslide dam, we considered a bottom width of 8 m, an initial head of 3.3 m and a slope coefficient of 1.5; for the Baige landslide dam, a bottom width of 3 m, an initial water head of 3 m and a slope coefficient of 1.3 were used. The results of these calculations are shown in Figure 10.

Figure 10. Analysis of discharge coefficients.

From the sensitivity analysis of the discharge coefficient, it is evident that c₁ and c₂ are positively correlated with the breach flow rate. Furthermore, an increase of water head in the breach significantly amplifies the sensitivity of these two parameters to flow rate. Additionally, the calculated values of c₁ and c₂ were compared with the initial discharge rates of two real-world cases. The results indicate that for the initial discharge of the Baige landslide dam, the values of both parameters are relatively consistent. However, for the initial discharge of the Tangjiashan landslide dam, the calculated discharge based on the suggested values for these parameters is smaller than the measured initial discharge. Nevertheless, Singh’s (1996) suggested values for the two parameters generally fulfil the requirements of the two real-world case studies discussed in this research.

5. Results and discussion

5.1. Calculation results

Based on the parameters in Tables 1 and 2, the draining process of Tangjiashan landslide dam and Baige landslide dam is calculated using the breach model proposed in this study, including breach flow, width and depth development of breach and reservoir water lever, as shown in Figures 11 and 12.

Figure 11. The failure process of Tangjiashan landslide dam: (a) the discharge process; (b) the breach development; (c) the reservoir water level development.

Figure 12. The failure process of Baige landslide dam: (a) the discharge process; (b) the breach development; (c) the reservoir water level development.

According to Figures 11 and 12, the calculated the peak value and peak time of the breach flow and the breach width are consistent with the measured results, and there is a roughly linear relationship between the first derivative of break flow and the breach expansion rate. The breach slopes of Tangjiashan landslide dam and Baige landslide dam collapsed 2 and 7 times, respectively. The calculated reservoir water level of Tangjiashan landslide dam agree well with the measured results. However, the calculated reservoir water lever in the late stage of discharge decreased faster than the measured one. The possible reason is the error in the storage capacity parameters of the Baige barrier lake. The calculated and measured results of key breaching outcomes, i.e. peak breach flow (Q_p), time to peak discharge (T_p), final breach top width (B), final breach bottom width (b), final breach depth (D) and measured data, are shown in Tables 3 and 4.

Table 3. Comparison of breaching outcomes of Tangjiashan landslide dam.

	Qb (m^3/s)	B (m)	b (m)	Tp (h)	D (m)
Measured value	7500	186.2	83.8	15.1	35
Calculated value	7079	181.7	66.4	15.3	42.2
Relative error	−5.6%	−2.4%	−20.8%	1.3%	+20.6%

Table 4. Comparison of breaching outcomes of Baige landslide dam.

	Qb (m^3/s)	B (m)	b (m)	Tp (h)
Measured value	31,000.0	264.1	107.8	37.25
Calculated value	31,492.5	250.7	106.3	35.5
Relative error	+1.6%	−5.1%	−1.4%	−4.7%

From Tables 3 and 4, it shows an agreement between the calculated values and the measured ones in breach flow, such as peak value and time. The calculated breach sizes of the Baige landslide dam align well with reality, with a relative error of around 5%. While there are certain discrepancies in the calculated breach sizes of the Tangjiashan landslide dam, particularly in the breach bottom width and depth, which are around 20%. There is no process monitoring data for the size of the breach, making it difficult to compare the development process and identify the cause of the calculation result error.

5.2. Comparison with existing models

In this study, the commonly used mathematical model for landslide dam overtopping and collapse process was selected and compared with the model proposed in this research. Two models proposed by Zhong et al. (2018a) and Chang and Zhang (2010) were selected for the calculation of the draining process of the Tangjiashan landslide dam, and the other two models proposed by Zhang et al. (2019) and Chen et al. (2020), respectively, were selected for the calculation of the discharge process of the Baige landslide dam and compared with the calculation results of proposed model. Tables 5 and 6, respectively, provide a comparison of the failure process of Tangjiashan landslide dam and Baige landslide dam using different breach models.

Table 5. Comparison between calculated and measured data of Tangjiashan landslide dam breach.

Item	Qb (m^3/s)	B (m)	b (m)	Tp (h)	D (m)
Model by Zhong et al. (2018a)	6299.34	179.00	106.64	14.94	–
	−16.0%	−3.9%	+27.3%	−1.1%	–
Model by Chang and Zhang (2010)	6737	247	101	–	45
	−4.8%	+32.7%	+20.5%	–	+28.6%
Model proposed in this article	7079	181.7	66.4	15.3	42.2
	−5.6%	−2.4%	−20.8%	1.3%	+20.6%
Measured value	7500	186.2	83.8	15.1	35

Table 6. Comparison of measured and calculated data of Baige landslide dam breach.

Item	Qb (m^3/s)	B (m)	b (m)	Tp (h)
Model by Zhang et al. (2019)	34,348	337.4	139.8	38.32
	+10.8%	+27.8%	+29.7%	+2.9%
Model by Chen et al. (2020)	31,041.1	–	–	37.72
	+0.1%	–	–	+1.3%
Model proposed in this article	31,492.5	250.7	106.3	35.5
	+1.6%	−5.1%	−1.4%	−4.7%
Measured value	31,000.0	264.1	107.8	37.25

The comparison of breach models for the calculation of draining process is listed in Tables 5 and 6. When calculating the failure of the Tangjiashan landslide dam, Zhong et al. (2018a) set the final depth of the breach as a fixed value, that is, stop cutting when the breach erodes to that depth, and Chang and Zhang (2010) did not list the peak flow time in his research. Therefore, the relevant calculation data is not listed in Table 5. Furthermore, the result provided by Chen et al. (2020) is the water surface width at the breach, while the top and bottom widths of the breach can’t be determined. Therefore, the final size of the breach for the model by Chen et al. (2020) was not considered as listed in Table 6.

The comparison analysis based on Tables 5 and 6 shows that all five models can reflect the overtopping failure of landslide dams. The model proposed in this article considers the non-uniformity of particle size distribution in landslide dams, especially focusing on the vertical distribution pattern of particles. It controls the predicted peak flow velocity and the occurrence time within a smaller relative error range of approximately ±5%. However, it is undeniable that more cases will be needed to comprehensively assess the performance of each model, and the parameters of the dam material, especially the values of soil erodibility, may introduce uncertainty to the calculation data.

5.3. Influence of particle size

From the preceding discussion, it is evident that dam breach calculations involve numerous parameters. To further evaluate the model’s performance and understand the correlation between parameters and the dam breach process by analyzing the extent of their influence on the calculation results, certain factors were chosen for sensitivity analysis. Among the many input parameters for dam breach, the particle size of the dam material is particularly significant as it directly affects the initiation flow velocity of particles, thus serving as the most intuitive parameter reflecting soil erosion. Therefore, using the Tangjiashan landslide dam as an example, the dam breach process was calculated by separately adjusting the median particle size of the surface particles of the dam body by ±50%. The calculation results are presented in Figure 13.

Figure 13. The failure process of Tangjiashan landslide dam under different median particle sizes: (a) discharge of breach; (b) development of breach.

From the calculation results of dam failure with different median particle sizes, it can be observed that the median particle size of the dam material significantly influences the development process of dam failure. As shown in Figure 13, with an increase in the median particle size, the peak flow rate of the breach decreases, while the peak time extends, and the water rise and fall speeds decrease continuously. Furthermore, the width expansion speed of the breach decreases, and the final residual breach size also decreases progressively.

However, as seen in Figure 13(b), cases 2 and 3 show two collapses throughout the entire dam failure process, while in case 1, there is only one collapse. This discrepancy can be attributed to that, in case 1, the decreased particle size led to stronger erosion, resulting in rapid expansion of the top and bottom of the breach, with a small difference in expansion speed and no very steep slope formed. Consequently, the number of collapses at the breach was small. The key dam failure output results for the three cases are listed in Table 7.

Table 7. Influence of particle size on the dam breach process.

Cases	dsur,50 (mm)	Qp (m^3/s)	Tp (h)	B (m)	b (m)
1	11.45	9621.07	10.46	196.16	77.77
	−50%	+35.91%	−31.64%	+7.96%	+17.12%
2 (Prototype)	22.9	7079	15.3	181.7	66.4
	–	–	–	–	–
3	34.35	5009.10	22.00	171.23	59.18
	+50%	−29.24%	+43.78%	−5.76%	−10.87%

Table 7 quantitatively demonstrates the extent of the influence of dam material particle size on the flow rate and development of the breach. It is clear that the particle size of the material significantly affects the flow rate, including both the peak flow rate and the peak time, but has a weaker impact on the width of the breach. Overall, the output results help mitigate the bias of the input data. However, the median particle size of landslide dams plays a crucial role in dam failure scenarios. As the consequences of various outcomes of flood disasters like dam failures differ significantly, corresponding emergency plans will also vary greatly.

5.4. Influence of particle distribution

When the material of the landslide dam exhibits non-uniformity in the vertical direction, the grading of the surface soil of the breach will vary at different heights of the breach, resulting in different development and flow processes of the breach. Thus, the focus of this article is on the influence of particle size distribution on the process of dam failure. To investigate the impact of different particle size distributions on the dam failure process, the Tangjiashan landslide dam was selected as a case study, and the dam failure development process under both inverse and uniform grading conditions was calculated separately. The calculation results are presented in Figure 14.

Figure 14. The failure process of Tangjiashan landslide dam under different particle distributions: (a) discharge of breach; (b) development of breach.

From Figure 14, it is evident that different particle size distributions influence the dam failure process. Specifically, compared to the inverse grading dam body, the uniform grading dam body maintains coarse particles in the breach as it develops to the middle and lower parts of the dam body. This reduces the erosion of the breach soil by the overtopping water flow, resulting in a slower development of the breach and smaller flow peaks, with a corresponding delay in the peak time. However, in both cases, the breached slope exhibited two collapses, and the final width of the residual breach did not differ significantly.

The quantitative impact of particle distribution on the flow rate and size of the breach is listed in Table 8. It is evident that particle distribution significantly affects the flow rate at the breach, with peak time being more sensitive to changes in particle distribution than peak flow rate. Moreover, particle distribution also impacts the size of the breach, with the top width being more sensitive to particle distribution than the bottom width.

Table 8. Influence of particle distribution on the dam breach process.

Cases	Distribution	Qp (m^3/s)	Tp (h)	B (m)	b (m)
1 (Prototype)	Inverse grading	7079	15.3	181.7	66.4
		–	–	–	–
2	Uniform grading	5582.50	18.90	174.20	60.82
		−21.14%	+23.53%	−4.13%	−8.40%

This highlights the importance of considering the particle distribution of the dam body when calculating the actual flood caused by landslide dam failure. Emergency plans are mainly based on the evolution process of the dam failure flood, and the flow process at the breach serves as the initial input condition for such calculations. As the flow rate at the breach is sensitive to the particle distribution of the dam body, it is crucial to reasonably account for the non-uniformity of dam materials for calculation of dam breach.

6. Conclusion

In this study, we propose a new model for predicting the process of overtopping failure in landslide dams. The new model considers the vertical distribution of soil particles that is commonly ignored in existing dam failure models, therefore, directly addressing the shorting coming of current models and enables the consideration of the uneven spatial distribution of dam material properties. This will significantly improve the accuracy and reliability of landslide dam failure predictions.

The new model is mainly composed of four modules: (1) a vertical distribution model based on the three parameter Weibull distribution, for characteristic particle sizes of dam material; (2) the broad-crested weir equation based on water balance, for breach discharge; (3) a new derived sediment transport formular of total load that separately considers the transport of bed load and suspended load, for dam material erosion; (4) the breach deformation equations that separately computed for continuous erosion and intermittent collapse, for both horizontal and vertical expansion of the breach.

The performance of the new model has been evaluated using two case studies and compared to four existing models. It has been found that the model can describe the development of the breach and discharge, and the calculation results conform to the general law of overtopping-induced landslide dam break. The comparison of the calculated results with the measured data indicates that the new model proposed in this study can simulate the overtopping process of the landslide dam. Compared to the four existing models, the model proposed in this study can more accurately predict the peak flow rate and the occurrence time within a smaller error range. This improved performance is attributed to its more thorough consideration of heterogeneity that is present in the composition of landslide dam materials. To further verify this conclusion, sensitivity analysis was conducted on the distribution of particle size, and the results showed that the particle size distribution has a significant impact on the dam failure process, especially on the flow rate at the breach, including the peak flow rate and peak time. This aspect critically influences the rationality of the emergency plan for dam failure. Specifically, for inverse grading landslide dams, failure to consider the distribution of particle size could lead to overly conservative conclusions in the calculation results. This underscores the significance of the model established in this article.

It should be emphasized that the model proposed in this article is constructed for landslide dams characterized by particle size sorting, such as the dam formed by remote landslide. Dams formed by collapse typically consist of mixed sediments, often retaining the original strata sequence, resulting in insufficient particle sorting. As a result, the particle size distribution model and dam break simulation method established in this study may not be suitable for dam break calculations of this type of dam. For now, the method to simulate material composition for dams formed by collapses is the major limitation for the prediction of the dam break process in all available breach models. In addition, the soil erodibility of landslide dam fill varies greatly due to differences in soil type, compaction degree and moisture content, leading to a high level of uncertainty. To reduce these uncertainties, developing next generation or new physically based dam breach models and conducting uncertainty analyses by utilizing relevant simulation methods is indispensable.

Acknowledgement

The authors greatly acknowledge the constructive comments and suggestions of two anonymous reviewers.

Disclosure statement

The authors declare no conflict of interest.

Data availability statement

The input parameters and output results of the breach model used for real-world dam break calculation in this article are available from the digital data repository in Zhao (2024).

Additional information

Funding

This work is financially supported by the Chongqing Water Conservancy Science and Technology Project (Grant No. CQSLK-2023003), the National Natural Science Foundation of China (Grant Nos. 52279095 and 52109149), the Australian Research Council via the Discovery Early Career Researcher Award (Grant No. DE210100117), the Major Science and Technology Program of the Ministry of Water Resources (Grant No. SKS-2022076) and the Chongqing Water Conservancy Science and Technology Project (Grant No. CQSLK-2022001).