## ABSTRACT

Landslide dams are formed by the blockage of rivers with landslides or any debris flows. Such dams are prone to failure and therefore, simulation of their failure and induced morphological changes is of great importance for hazard mitigation. Multi-layer and non-uniform sediments are usually observed in landslide dam materials, but these are rarely considered in simulations. In this paper, a two-dimensional numerical model considering multi-layered and non-uniform sediments is proposed. The model solves the shallow water equations, the Exner equation considering the different size fractions and Hirano's active layer equation based on a finite volume scheme in a weakly coupled approach. Two experimental tests are used to assess the applicability and accuracy of the model. The results show that it can well simulate the flow and morphological changes, as well as the surface coarsening and fining phenomena related to multi-layer and non-uniform sediments. The numerical results of dam failure process and morphological changes are in overall good agreement with the experiments, but the peak discharge and bed elevation tend to be underestimated when finer material is considered. Finally, the influence of sediments fractions and active layer depth on numerical results, as well as the limitations of the model are discussed.

## HIGHLIGHTS

A numerical model considering multi-layered and non-uniform sediments is developed.

The applicability and accuracy of the model in simulating dam failure and dam-break flow with multi-layered and non-uniform sediments are verified.

The influence of sediments fractions and the active layer depth on numerical results is discussed.

## INTRODUCTION

Landslide dams are natural dams that are formed by the blockage of river channels with landslides, rockfalls, or debris flows (Costa & Schuster 1988). The failure of landslide dams is one of the essential issues related to disasters in mountainous areas. The flood induced by the failure of such a dam may bring catastrophic disasters in the downstream area, and significantly alter the local topography and environment (Korup 2004). For instance, the flood induced by the failure of Yigong landslide has caused hundreds of deaths in the downstream area in India (Delaney & Evans 2015). The failure of Tangjiashan landslide dam, which was formed by the 2008 Wenchuan earthquake in China, has induced several meters of deposits and dramatically altered the local morphology (Cui *et al.* 2010). More recently, the flood induced by the failure of two Baige landslide dams, which occurred in 2018 at the boundary of Tibet and Sichuan Province, China, has inundated quantities of farmland and resident area over 600 km in the downstream area (Yang *et al.* 2022). Therefore, analysis of dam failure-induced flooding and local morphological changes is of great importance for hazard management.

Research on landslide dam failure and induced morphological change are carried out based on field investigation, experimental and numerical studies (Coleman *et al.* 2002; Casagli *et al.* 2003; Cao *et al.* 2011), the latter being widely used. The numerical study of landslide dam failure can be further divided into the following two categories: simplified models and physics-based models. Simplified models are usually used to predict the discharge hydrograph during dam breaching and the results can be further used in flood routing analysis in combination with hydrodynamic software, such as HEC-RAS. These models use empirical erosion rate equations and the assumption of a broad crested weir to calculate the discharge and breach parameters, such as the DABA model and the DB-IWHR model (Zhong *et al.* 2018). Although these models are efficient in simulating landslide dam failure, several input parameters and also the empirical breach evolution pattern have to be determined before using the model. As the results are highly sensitive to these choices, this may lead to a series of uncertainties during the simulation. In addition, as the flood routing in the downstream area is simulated separately, the morphological changes induced by the dam failure are not considered.

Physics-based models, in which the hydrodynamic and morpho-dynamic equations are considered, have been proposed to simulate landslide dam breaching and induced morphological changes (Balmforth *et al.* 2008; Cao *et al.* 2011). These models can well represent the interaction between flow and sediments and are widely used in dam breaching simulations (Soares-Frazão & Zech 2010; Chen & Zhang 2015; Juez *et al.* 2014; Meurice & Soares-Frazão 2020). However, in these models, the sediments are treated as a uniform material with a unique median diameter (d50), thus, leading to inaccurate simulation results. In recent years, several numerical models considering non-uniform sediments have been developed. These models can be divided into the following two categories: the two-phase models or the quasi-single-phase models (Li *et al.* 2018; Xia *et al.* 2018; Martínez-Aranda *et al.* 2022) and the active layer models (Stecca *et al.* 2014; Juez *et al.* 2016; Chavarrías *et al.* 2019, 2022).

In both types of numerical models, the non-uniform sediments are divided into several fractions with different median diameters (*d*_{50}), while different algorithms are considered to calculate the transport of these sediments. In the two-phase model or the quasi-single-phase model, the flow and sediments are assumed to be well mixed in the flow, and the models solve the solid and fluid phases using distinct equations. The mass conservation equation is then considered in each of the sediment fractions associated with the non-uniform sediments (Pudasaini 2012; Johnson *et al.* 2012). These models have been widely used to simulate the propagation of sediment laden flow such as mudflow or debris flow, while ignoring the variation of the bed material composition during flow propagation. In terms of the active layer model, two layers of sediments are considered in the erodible bed, including the upper active layer and the bottom substrate layer. The interactions between the flow and the sediments are mainly considered in the upper active layer. Thus, the sediment transport rate, the evolution of morphology and the change of surface grain size distribution can be calculated within the active layer based on the sediment transport rate of each fraction and its percentage in the total sediment (Hirano 1971).

The active layer model has been widely used to simulate the flow propagation and morphological changes related to a non-uniform erodible bed, but not under severe transient flow conditions. However, when considering dam failure and the induced morphological changes, especially in the case of landslide dams, the material composition is one of the main characteristics that significantly influences the dam failure and the sediment transport (Mei *et al.* 2021). Current numerical models rarely consider the material characteristics of landslide dams, especially their non-uniform material composition and multi-layered characteristics, leading to the inaccurate results associated with the peak discharge, breach evolution duration during dam failure.

In this paper, a two-dimensional numerical model is developed with non-uniform and multi-layer sediments to simulate dam failure and induced morphological changes. The model is solved in a weakly coupled way based on a finite volume scheme by Soares-Frazão & Zech (2010), and two flume tests are used to verify the applicability and accuracy in simulating dam failure and morphological changes. The paper is organized as follows: In Section 2, the three parts of the numerical model, including the hydrodynamic module, the morpho-dynamic module and the active layer module are introduced. Then, the numerical scheme is explained in Section 3. In Section 4, the proposed model is verified by simulating two laboratory tests: (i) dam-break flow test over movable bed and (ii) dam failure test with different non-uniform sediments. In Section 5, the influence of the sediment fraction and active layer depth on numerical results, and also the limitations of the model are discussed. Finally, the conclusions arising from the current work are described in Section 6.

## MODEL CONCEPTION

The numerical model in which the non-uniform and multi-layer sediments are considered is shown in Figure 1. In this model, the sediments in the movable bed are divided into several layers numbered from #1 to #n, with the upper active layer included in layer #1. In each layer, sediments are further divided into fractions with different median diameter and percentage content. The numerical model is composed of three parts, including the hydrodynamic module, the morpho-dynamic module and the module for the variation of surface grain size distribution. The governing equations and characteristics of each module are introduced in the following section.

### Governing equations

#### Hydrodynamic module

*h*is the water depth,

*u*and

*v*are the flow velocity in the

*x*and

*y*directions, respectively,

*g*is the gravitational acceleration, and are the bed slopes in

*x*and

*y*directions, with being the bed level, andare the friction slopes in the

*x*and

*y*directions, with

*n*being the Manning coefficient.

#### Morpho-dynamic module

*N*is the total number of sediments fractions, is the percentage content of the

*i*th sediment fraction (Figure 1), and are the unit sediment transport rates of the

*i*th sediment fraction in the

*x*and

*y*directions. These transport rates are calculated using a closure equation as described in the following.

#### Surface grain size variation module

*i*from the

*N*fractions:where is the active layer depth. The left-hand side of Equation (5) refers to the update of the percentage content of each sediment fraction in the active layer. In the right-hand side of the equation, the first term accounts for the interactions between the bed load and the active layer, while the second term refers to the sediment exchange between the active layer and the substrate layer. The coefficient in the second term is an empirical parameter that can be determined by considering whether the bed is under degradation or aggradation condition. If the bed degrades, the active layer moves downwards, then the exchange between the active and substrate layers is determined by the substrate layer, while if aggrades, the sediment exchange is controlled by the bed load and active layers. Therefore, can be expressed as follow:where

*a*is an empirical parameter ranging from 0 to 1 and is the bedload transport rate of the

*i*th sediment fraction among the total bedload sediment transport rate in the bedload transport layer and can be calculated as follows:

*i*th sediment fraction. The function linking the median diameter of the total mixture and that of each sediment fraction can either be determined by a simplified linear relationship or using a regression analysis.

### Closure equations

The numerical model consists of *N**+**4* equations, i.e. Equations (1)–(4) and the *N* Equation (5) for the *N* fractions, with *N**+**4* unknow parameters. Closure equations are required for the sediment transport rate *q _{s}*, the active layer depth

*L*and the Manning coefficient

_{a}*n*of the bed that changes according to bed surface composition.

#### Sediment transport rate

*n*is the Manning friction coefficient. Taking the hiding-exposure effect into consideration, the formula proposed by Ashida & Michiue (1972) is used to calculate the dimensionless critical shear stress in each fraction. The formula can be expressed as follow:where is the dimensionless critical shear stress of the whole sediment mixture, which is directly determined by the median diameter of the sediment mixture and can be described as follow (Wu

*et al.*2000):withbeing the non-dimensional particle size of the mixture with median diameter

*d*

_{50}.

#### Active layer depth

*k*is an empirical coefficient ranging from 1 to 3 (Juez

*et al.*, 2016). In addition, since multiple sediment layers are considered in the numerical model, the variation of the active layer depth is described according to three scenarios associated with the change of layers, as shown in Figure 2. In Scenario (a), the active layer is totally located in the upper layer, then it can be directly calculated by Equation (13). In Scenario (b), the active layer overlaps two layers, then is the maximum value according to Equation (13) considering the

*d*

_{90}of the two layers. In Scenario (c), the active layer is located directly above the fixed bed, in this way, equals to the depth of the movable sediments.

#### Manning coefficient

*U*is the flow velocity, is the energy slope, , , and are the hydraulic radiuses corresponding to the sediment mixture, to the first and the th fraction of sediment, respectively. With the Einstein assumption, the hydraulic radius can be further expressed as:

## NUMERICAL SCHEME

The model equations are solved based on a first order finite volume scheme over an unstructured triangular mesh (Soares-Frazão & Zech 2010). As the solution of the active layer model coupled with the shallow water equations may be ill-posed due to the change of mathematical character under specific conditions (Chavarrías *et al.* 2019), a weakly coupled approach proposed by Juez *et al.* (2014, 2016) is used here to solve the governing equations. In this approach, the shallow water equations, Exner equation and Hirano's active layer equation are solved separately. The details of the numerical scheme of each module are given in the following.

### Hydrodynamic numerical solution

*n*)

_{x}, n_{y}*the unit vector normal to the interface, Equations (1)–(3) then can be written in a matrix form as (Soares-Frazão & Zech 2010):where is conservation term in the local co-ordinate system, and are the fluxes term perpendicular and along the interface, respectively, and*

^{T}*S*is the source term. These four vectors are expressed as:

with being the matrix related to the co-ordinate system transformation.

*m*, refers to the time step, is the cell area,

*k*is the interface number ( in a triangular mesh), is the th-interface length, is the relevant flux term through the

*j*th interface and can be calculated using the lateralized HLLC scheme as follow:where

with , and being the hydrodynamic characteristics related to the Jacobian matrix of shallow water equations (Soares-Frazão & Zech 2010).

### Morpho-dynamic numerical solution

*et al.*(2014, 2016) and Meurice & Soares-Frazão (2020), a morpho-dynamic numerical characteristic is proposed to solve the Exner equation using a finite volume with HLLC-type flux approach. Equation (4) can be written in the local co-ordinate system attached to the considered interface with as:wherewith being the unit sediment transport rate of the

*i*th sediment fraction in the direction perpendicular to the interface.

### Active layer model solution

### Numerical stability

*C*is the CFL number usually taken as , is a reference length in the computational cell (e.g., the radius of the incircle for a triangular cell), is the maximum value of numerical characteristics related to the shallow water equations, the modified Exner equation and the Hirano's active layer equation, and can be calculated as follow:

## MODEL VALIDATION

In this section, a synthetic test and two laboratory tests are considered to assess the accuracy and suitability of the proposed numerical model, and also the clarify the active layer assumption. The synthetic test is a sediment feed test and is used to clarify if the model can adjust to the equilibrium state under initial aggradation and degradation conditions considering non-uniform sediments. The first experimental test is used to investigate a dam-break flow and its induced morphological changes over movable bed with uniform sediments (Soares-Frazão *et al.* 2012). Based on this test, three numerical scenarios with both uniform and non-uniform sediments are analyzed. The second experimental test aims to investigate landslide dam failure and the induced morphological changes, in which three scenarios with different non-uniform sediments are considered. Furthermore, the numerical results related to the variation of water level and discharge during dam failure, the dam breaching process and the downstream morphological change are compared with the experimental data.

### Synthetic test: sediment feed test under aggradation and degradation conditions

The synthetic test was carried out in a 4 m long and 10 m wide rectangular channel with a triangular mesh of 0.5 m (average edge length). The bedrock level at the bottom of the channel remained at 0, and the channel was filled with sediments with a specific bed slope. The following two numerical tests were considered: first a test considering uniform sediments and then a test considering non-uniform sediments with two fractions. The sediments used to fill the channel in both tests had a median diameter of 1.7 mm, a Manning coefficient of 0.00167 s/m^{1/3} a porosity of 0.44 and a relative specific gravity of 2.65. A constant unit water inflow rate and unit sediment feed rate of 0.05 and 0.00098 m^{3}/s/m, respectively, were imposed at the upstream channel boundary, while a transmissive boundary was considered at the downstream boundary of the channel. Based on the transport equation for calculating the sediment transport rate (Equation (9)) in the present model, the equilibrium bed slope was estimated to be 5%, with an initial water depth of 0.035 m. Furthermore, two additional bed slopes of 6 and 4%, which represent the aggradation and degradation conditions, were considered as initial conditions to test the accuracy and stability of the present model. It should also be noted that in the numerical test considering non-uniform sediments, the median diameter of each fraction was 1 and 2.4 mm, respectively, and the percentage content of each fraction was equal to 50%. Therefore, the corresponding unit sediment feed rate of each fraction was estimated to be 0.0005 and 0.00048 m^{3}/s/m, respectively.

### Test 1: Dam-break flow over movable bed

*et al.*(2012). In this test, the upstream and downstream water levels were 0.470 and 0.085 m, respectively, and a movable bed extending over 1.5 m upstream and 9.0 m downstream of the gate was considered. The movable bed had a depth of 8.5 cm, which consisted of uniform sand with a median diameter of 1.61 mm, a relative specific gravity of 2.63 and a porosity of 0.42. In the present simulations, this case is referred to as Scenario 1. Furthermore, another two scenarios considering three layers of uniform and non-uniform sediments (Scenarios 2 and 3) were designed based on Scenario 1 to verify the ability of the model to simulate the flow and the entrainment of non-uniform and multi-layer sediment induced by a dam-break flow. The details of the three numerical scenarios are shown in Table 1. In Scenario 2, three layers of uniform sediments are considered, and the sediments in each layer is as same as in Scenario 1. In Scenario 3, three layers of non-uniform sediments with depth of 0.028, 0.028 and 0.029 m are considered. The sediments in each layer consist of three fractions (with a median diameter of 0.5, 2 and 9 mm, respectively). The median diameter of the three-layer sediment ranges from 1.85 to 5.84 mm by varying the percentage content of each fraction. The manning coefficient of Scenarios 1 and 2 is the same as the experimental value (0.0165 s/m

^{1/3}), while for Scenario 3, the coefficient is calculated with Equation (16).

. | . | Parameters . | . | . | . | . | . | . | . | |||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Scenario . | L
. | d_{50}/mm
. | s
. | n/sm^{−1/3}
. | . | F
. | P
. | d_{50}/mm
. _{,i} | ||||

S1 | 1 | 1.61 | 2.63 | 0.0165 | 0.42 | 1 | 1 | 1.61 | ||||

S2 | 3 | 1.61 | 2.63 | 0.0165 | 0.42 | 1 | 1 | 1.61 | ||||

1 | 1 | 1.61 | ||||||||||

1 | 1 | 1.61 | ||||||||||

S3 | 3 | 1.85 | 2.63 | 0.0174 | 0.4 | 3 | 0.57 | 0.33 | 0.10 | 0.5 | 2 | 9 |

3.89 | 2.63 | 0.0185 | 0.3 | 3 | 0.33 | 0.33 | 0.34 | |||||

5.84 | 2.63 | 0.023 | 0.25 | 3 | 0.10 | 0.33 | 0.57 |

. | . | Parameters . | . | . | . | . | . | . | . | |||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Scenario . | L
. | d_{50}/mm
. | s
. | n/sm^{−1/3}
. | . | F
. | P
. | d_{50}/mm
. _{,i} | ||||

S1 | 1 | 1.61 | 2.63 | 0.0165 | 0.42 | 1 | 1 | 1.61 | ||||

S2 | 3 | 1.61 | 2.63 | 0.0165 | 0.42 | 1 | 1 | 1.61 | ||||

1 | 1 | 1.61 | ||||||||||

1 | 1 | 1.61 | ||||||||||

S3 | 3 | 1.85 | 2.63 | 0.0174 | 0.4 | 3 | 0.57 | 0.33 | 0.10 | 0.5 | 2 | 9 |

3.89 | 2.63 | 0.0185 | 0.3 | 3 | 0.33 | 0.33 | 0.34 | |||||

5.84 | 2.63 | 0.023 | 0.25 | 3 | 0.10 | 0.33 | 0.57 |

*Note*: *L* is the number of layers and *F* is the number of sediments fractions in each layer, *P* is the percentage content of each sediments fraction, *d*_{50} is the median diameter of each layer, *s* is the relative specific gravity, is the porosity and *n* is the Manning coefficient, *d*_{50,i} is the median diameter of each sediment fraction.

Numerical simulations for the three scenarios were conducted on an unstructured triangular mesh with a spatial resolution of 0.1 m. Two monitor points located at (0.64, −0.5) for Point 1 and (1.94, −0.33) for Point 2 are selected to compare the water level variation that is obtained from Scenario 1, Scenario 2 and the experimental results. The final topography and the bed elevation along four profiles (*y* = 0, *y* = 0.2 m, *y* = 1 m and *y* = 1.45 m) are compared between the numerical and experimental results. Furthermore, the surface grain size distribution and sediment transport rate related to the non-uniform and multi-layer sediments are analyzed based on the results of Scenario 3.

#### Numerical results of Scenarios 1 and 2

Based on the analysis of Scenarios 1 and 2, it can be concluded that the proposed numerical model can well simulate the dam-break flow and the induced morphological change. In addition, since the results of Scenario 1 are identical to those of Scenario 2, the accuracy of the multi-layer numerical scheme applied to uniform sediments in the proposed numerical model can be confirmed.

#### Numerical results of Scenario 3

*d*

_{50,i}= 0.5 mm). As shown in the first column of Figure 9, the movable bed sediments in the first layer are mainly eroded at

*t*= 5 s. The erosion and deposition of finer sediments in the first layer lead to a slight decrease of the percentage content of the finest sediments on the main channel, further leading to the increase of median diameter on the surface of the movable bed. Meanwhile, at the corner where the flume widens abruptly, the erosion depth is larger than the depth of the first layer, therefore the erosion of the following sediment layers leads to the increase of the surface median diameter at this area. With the increase of erosion depth and area (

*t*= 10 s), the percentage content of finest sediments at the central downstream of the gate decreases from 0.57 to 0.33, indicating that the second layer is eroded in this area. In the meantime, the surface median diameter at the corner where flume suddenly widens decreases to zero due to the total erosion of the movable bed (i.e., the fixed bed has been reached), while around the corner, it greatly increase due to the deposition of the sediments in the third layer. At the final state of the simulation (

*t*= 20 s), sediments in different layers are eroded at different locations on the movable bed. The surface median diameter mainly increases on the movable bed, and the increase rate decreases along the downstream of the flume. During the whole simulation, the increase of surface median diameters induced by the change of sediment layer is much larger than that induced by the variation of percentage content of each sediment fraction in each layer.

*y*= 0) are analyzed, as shown in Figure 10. The sediment transport rate in each fraction first increases and then decreases along the flow direction. In addition, the transport rate in the sediment fraction with finer material is much larger, and it decreases with the increase of the median diameter in each sediment fraction. The difference of sediment transport rate in each fraction leads to the variation of surface median diameter. As the larger transport rate is observed in the finer sediment fraction, the percentage content of coarsest sediment fraction increases gradually, and finally leading to the increase of the median diameter on the surface of the movable bed. In addition, in Figures 10(b)–(d), a sudden decrease of sediment transport rate occurs at

*x*= 2.0 m, where the change of sediment layers and also the mixing of sediments with different diameters are observed, as shown in Figure 9. With these two phenomena, the percentage content of finest sediment decreases around the boundary of sediment layers, while it increases downstream of the boundary, further leading to the variation of sediment transport rate in each sediment fraction.

Based on the analysis of Scenario 3, it can be concluded that the proposed numerical model can well simulate the variation of surface median diameter relate to multi-layered and non-uniform sediment.

### Test 2: Dam failure with different non-uniform sediments

*et al.*2024). The set-up of the test is shown in Figure 11. The bottom of the flume consisted of 1.5 m fixed bed and 3.5 m movable bed with a depth of 0.08 m, and the modeled landslide dam was located at 2 m downstream of the flume. The dam was in a trapezoidal shape, with a height of 0.25 m, a top width of 0.12 m and a bottom width of 1.08 m. In addition, a diversion channel with a depth of 0.05 m and a width of 0.06 m was considered at the dam crest. The slope gradient of the flume was 0.016 and the inflow rate remained constant with 0.001 m

^{3}/s during the test. For the material composition of dam and movable bed, three materials with sand content (

*d*< 2 mm) ranging from 30% - 70% were considered, as shown in Figure 12. Based on these three materials, three numerical scenarios are designed as described in Table 2. The material in each scenario is divided into two fractions with same percentage content. A relationship between the median diameter of the sediment mixture and the median diameter and percentage content of each sediment fraction is determined based on a polynomial regression from the grain size distribution, as shown in Figure 12(b). The Manning coefficient and the active layer depth are calculated based on Equations (16) and (12), respectively.

Scenario . | Parameters . | ||||||||
---|---|---|---|---|---|---|---|---|---|

N
. | f
. _{i} | d_{50}/mm
. _{,i} | s
. | n/sm^{−1/3}
. | . | L/m
. _{a} | |||

S2-1 | 2 | 0.5 | 0.5 | 14.74 | 01.27 | 2.65 | 0.026 | 0.18 | 0.03 |

S2-2 | 9.00 | 0.62 | 0.024 | 0.25 | |||||

S2-3 | 4.02 | 0.36 | 0.022 | 0.31 |

Scenario . | Parameters . | ||||||||
---|---|---|---|---|---|---|---|---|---|

N
. | f
. _{i} | d_{50}/mm
. _{,i} | s
. | n/sm^{−1/3}
. | . | L/m
. _{a} | |||

S2-1 | 2 | 0.5 | 0.5 | 14.74 | 01.27 | 2.65 | 0.026 | 0.18 | 0.03 |

S2-2 | 9.00 | 0.62 | 0.024 | 0.25 | |||||

S2-3 | 4.02 | 0.36 | 0.022 | 0.31 |

*Note*: *N* is the number of sediment fractions, *f _{i}* is the percentage content of each fraction,

*d*

_{50,i}is the percentage content of each fraction,

*s*is the relative specific gravity,

*n*is the Manning coefficient, is porosity and

*L*is the active layer depth.

_{a}Numerical simulations of this dam failure test were conducted on an unstructured triangular mesh with a spatial resolution of 0.02 m. In addition, in order to simulate the lateral slope collapse during the test, a bank failure module is applied in our simulation. The details of the bank failure module are given by Swartenbroekx *et al.* (2010) and the relevant parameters used in our simulation are shown in Table 3. In the simulation, the critical and residual angles are the same, the values below the water level is determined by the repose angle of the material, while the values above the water level are derived from experimental observation as steeper slopes can be achieved because of apparent cohesion due to moisture of the material. Based on the three scenarios in Table 2, the parameters change, dam failure process and also the downstream morphological changes are analyzed and compared with the experimental data.

Scenario . | Critical angle(°) . | Residual angle(°) . | ||
---|---|---|---|---|

Above water . | Under water . | Above water . | Under water . | |

S2-1 | 87 | 37.8 | 87 | 37.8 |

S2-2 | 75 | 36.6 | 75 | 36.6 |

S2-3 | 50 | 33.1 | 50 | 33.1 |

Scenario . | Critical angle(°) . | Residual angle(°) . | ||
---|---|---|---|---|

Above water . | Under water . | Above water . | Under water . | |

S2-1 | 87 | 37.8 | 87 | 37.8 |

S2-2 | 75 | 36.6 | 75 | 36.6 |

S2-3 | 50 | 33.1 | 50 | 33.1 |

*Note*: The critical and residual angle are same and are determined by experimental data and the rest angle of the sediments.

#### Parameters change and dam failure process

^{−3}m

^{3}/s, while the related value obtained from experimental data is 5.41 × 10

^{−3}m

^{3}/s.

*t*= 80 s, the dam profile obtained from the numerical results is much higher than that of the experimental results, and the slope of the erosion surface is also much smaller. However, this discrepancy has little influence on the final profile after dam failure.

The difference between the numerical and experimental results when the finer material is considered can be attributed to the dam failure characteristics. The erosion of the dam mainly occurs as surface erosion when the coarser material is considered. However, for the dam with the finer material, the erosion of the dam may be transformed form surface erosion to backward erosion with mass failures of blocks of material as it happens for cohesive material, further leading to a steeper erosion surface during dam failure. This steeper erosion surface cannot be properly simulated with the present model, thus, leading to the difference between the numerical and experimental results.

#### Morphological changes

*y*= 0.1 m and

*x*= 3.5 m (Figure 15) are selected to analyze the bed elevation. The comparison between numerical and experimental results of different scenarios is shown in Figure 16. In the different scenarios, the numerical bed elevation along the profiles is always little bit lower than that of the experimental results, this may be attributed to the fact that only two sediment fractions are considered in the simulation. Therefore, the influence of the largest sediments (20–40 mm) may be ignored, leading to a longer transport distance and thus less deposition. However, the difference between numerical and experimental bed elevation is relatively small, and the present model is considered able to simulate the bed variation induced by dam failure. Still, the numerical and experimental bed elevations on Profile 2 (Figures 16(b), (d) and (f)) show larger differences when the finer material is considered, as was already observed in the bed topography simulation: the difference in dam failure simulation leads to a difference in the creation of the erosion channel in the middle part of the movable bed, further leading to the observed difference in bed elevation along Profile 2.

## DISCUSSION

### Influence of sediment fractions and active layer depth on numerical results

In this section, we consider S2-2 as the reference case, and four additional numerical scenarios are considered to analyze the influence of sediment fractions and active layer depth on the numerical results. The four additional numerical scenarios based on S2-2 are described in Table 4. In SA-1 and SA-2, one and four sediment fractions are considered respectively and the others parameters remain as those in S2-2. In addition, the four sediment fractions in SA-2 are obtained by sub-dividing the two fractions considered in S2-2 around the median diameter of each fraction. In SA-3 and SA-4, only the active layer depth is changed from 0.03 m to 0.02 and 0.04 m. The relationship between the median diameter and the fine content in SA-2 is the same as that in Figure 12(b), but the fine content here is the sum of the last two fractions shown in Table 4.

Scenario . | N
. | f
. _{i} | d_{50}/mm
. _{,i} | ||||||
---|---|---|---|---|---|---|---|---|---|

SA-1 | 1 | 1 | 2.00 | ||||||

SA-2 | 4 | 0.25 | 0.25 | 0.28 | 0.22 | 19.80 | 6.03 | 1.00 | 0.18 |

SA-3 | 2 | 0.5 | 0.5 | 9.00 | 0.62 | ||||

SA-4 | 2 |

Scenario . | N
. | f
. _{i} | d_{50}/mm
. _{,i} | ||||||
---|---|---|---|---|---|---|---|---|---|

SA-1 | 1 | 1 | 2.00 | ||||||

SA-2 | 4 | 0.25 | 0.25 | 0.28 | 0.22 | 19.80 | 6.03 | 1.00 | 0.18 |

SA-3 | 2 | 0.5 | 0.5 | 9.00 | 0.62 | ||||

SA-4 | 2 |

Note: *N* is the number of sediment fractions, *f _{i}* and

*d*

_{50,i}are the percentage content and median diameter of each sediment fraction, respectively.

### Limitations

The numerical results in this paper are in good agreement with the experimental data, but there are still some limitations. First, the accuracy in simulating dam failure with fine sediments by using the present numerical model could be improved. Indeed, the model cannot well capture the steep slope induced by backward erosion and mass failures when the fine content is relatively large in the material composition. This would require further assumptions and the implementation of geomechanical failure mechanisms adapted to cohesive sediments in the present numerical model. Secondly, the model is unable to properly consider the remobilization and redeposition of the sediments. Indeed, in the numerical scheme related to the change of layer, only the information related to the sediments in the active layer is recorded, which is sufficient for dam failure simulation. However, when the sediments are remobilized or redeposited, this may cause inaccuracies in the simulation results related to morphological changes and surface grain size distribution. In addition, the deposition of sediments with different diameters during the simulation should be further optimized. Third, up to now, only two experimental tests are simulated to validate the numerical model. Further applications should be considered in the future, especially related to real dam failure cases.

## CONCLUSION

In this paper, a two-dimensional numerical model considering multi-layered and non-uniform sediments is proposed to simulate landslide dam failures and the induced morphological changes. The model solves the shallow water equations, the Exner equation and Hirano's active layer equation based on a weakly coupled finite volume scheme. Two experimental tests are used to testify the applicability and accuracy of the present numerical model.

The first test considers a dam-break flow case with movable bed composed of uniform sediments. Three numerical scenarios were tested, with one-layer uniform sediments, multi-layer uniform sediments and multi-layer non-uniform sediments. The model can well simulate the water level and morphological changes related to uniform sediments, and the numerical results are in good agreement with the experimental data. In addition, the numerical model is able to simulate the surface coarsening and fining phenomena related to multi-layer and non-uniform sediments. The variation of the median diameter on the surface of the movable bed can be attributed to the change of sediment layers and the percentage content of each sediment fraction in each layer, among which, the change of sediment layers shows more obvious influence on surface grain size distribution.

The second series of tests aims at replicating a landslide dam failure with a wide grain size distribution and the induced morphological changes, and three scenarios with different sand content in the materials are analyzed. The numerical results of parameters change, dam failure process and the morphological changes are in good agreement with those of the experimental data when the material with smaller sand content are considered. With the increase of sand content in the material, the numerical model tends to inaccurately simulate the dam failure process, further leading to an underestimation of the peak discharge and of the decrease rate of the upstream water level, as well as to the inaccurate results related to the final bed elevation.

The influence of sediment fractions and active layer depth on numerical results is further analyzed based on the landslide dam failure test. Compared to the scenario with one sediment fraction, scenarios with two and four sediment fractions can better simulate the dam failure and the induced morphological changes. However, the difference in scenarios with two and four sediment fractions are relatively small. The active layer depth shows no influence on the dam failure and the morphological changes. The numerical surface grain size distribution is in good agreement with the experimental result. In addition, a larger median diameter in the zones where the median diameter increases and a smaller one in the zone where the median diameter decreases are observed with the increase of sediment fractions and the decrease of the active layer depth.

Future research work will focus on improving the simulations for finer sediments, the erosion and deposition of sediments and the real dam failure cases.

## ACKNOWLEDGEMENTS

The authors acknowledge the support from the Natural Science Foundation of China (No. 41731283 and No. 42307196) and the China Postdoctoral Science Foundation (2023M733029).

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

## CONFLICT OF INTEREST

The authors declare there is no conflict.