A Numerical Study of the Critical Threshold for Landslide Dam Formation Considering Landslide and River Dynamics

Liu, Wei; Hu, Yu-xiang; He, Si-ming; Zhou, Jia-wen; Chen, Kun-Ting

doi:10.3389/feart.2021.651887

ORIGINAL RESEARCH article

Front. Earth Sci., 28 May 2021

Sec. Geohazards and Georisks

Volume 9 - 2021 | https://doi.org/10.3389/feart.2021.651887

This article is part of the Research Topic Water-Related Natural Disasters in Mountainous Area View all 39 articles

A Numerical Study of the Critical Threshold for Landslide Dam Formation Considering Landslide and River Dynamics

Wei Liu^1,2,3

Yu-xiang Hu³

Si-ming He^1,2

Jia-wen Zhou³

Kun-Ting Chen^4,5*

¹Key Laboratory of Mountain Hazards and Surface Process, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China
²University of Chinese Academy of Sciences, Beijing, China
³State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu, China
⁴Compound Disaster Prevention Research Center, General Research Service Center, National Pingtung University of Science and Technology, Pingtung, Taiwan
⁵Department of Soil and Water Conservation, National Pingtung University of Science and Technology, Pingtung, Taiwan

Landslide dam formation can be influenced by the erosive capacity of river flow and the dynamic characteristics of the landslide. When the deposition rate of a landslide that reaches a river is higher than the erosion rate of river flow, the landslide can form a dam by blocking the channel. Hence, in this paper, a dimensionless discharge threshold for landslide dam formation considering landslide and river dynamics is established and studied numerically. A two-layer depth-averaged model coupled with an erosion term is presented to simulate river and landslide movements and their interactions. Several numerical cases are simulated to study the influence of landslide and river dynamics on the critical threshold for dam formation by considering some key factors, such as landslide velocity and the angle between the river and landslide transport directions. Through the simulations, three types of landslide intrusion into river can be reflected: a dam forms quickly, a dam forms or does not form close to a critical state, and no dam forms. The results show that these factors together affect the process of dam formation if the difference between the landslide and river discharges is relatively small. All results are helpful to further clarify the formation of such dams for natural hazard prevention under future climate change conditions.

Introduction

Landsides occurring in river valleys have the potential to block river channels by forming dams, and a cascade of negative consequences, such as dam-break floods and debris flows, can be induced if a dam forms (Romeo et al., 2017; Liang et al., 2019). Such a chain effect can greatly enlarge the broad scope and destructive power of disasters, resulting in serious economic loss and high casualties. Recent examples include the 2014 Bujumbura floods resulting from the failure of a landslide dam in Burundi, which caused 64 casualties and destroyed more than 940 houses Nibigira et al. (2018), and the 2018 Baige landslide dam in Southwest China, which caused economic losses of approximately RMB 74.3 billion (Deng et al., 2019). Finding the critical condition for forming a landslide dam has become a key issue that needs solution in landslide dam disaster prevention.

Numerous studies have demonstrated that geomorphological features Korup (2004); Fan et al. (2012); Chen and Chang (2016); Chen et al. (2021) and hydrological conditions Zhao et al. (2017); Chen et al. (2019a); Liao et al. (2019) jointly determine whether a landslide dam can be formed. Based on these studies, the following three conditions are essentially mandatory to form a landslide dam: first, the landslide needs to cross the river channel; second, the erosive rate of river flow must be smaller than the depositional rate of the landslide; and three, the thickness of landslide deposits in the river channel must be greater than the water depth (Fan et al., 2014; Chen and Chang, 2016; Chen et al., 2019b). To find the threshold, (e.g. landslide runout distance and dam height) for satisfying these conditions, empirical database analysis Fan et al. (2012); Chen and Chang (2016) and experimental measurements Okada and Uchida (2014); Chen and Orense (2020) are widely used, and several dimensionless critical indexes that are composed of variables characterizing the different elements involved, (e.g. landslide and river) have been proposed; these indexes have high significance (Ermini and Casagli, 2003; Fan et al., 2012; Dal Sasso et al., 2014). However, although the indexes can be used to forecast and discriminate between possible dam evolutions, they cannot quantitatively describe the formation process of landslide dam. From this point of view, numerical modeling in recent years has been devoted to studying landslide dam formation by using different physical models (Liu and He, 2016; Zhao et al., 2017; Chen and Wu, 2018; Zhao et al., 2019; Li et al., 2020). In summary, these models apply two distinct equations with corresponding rheological properties to describe the dynamics of the landslide and the river. Furthermore, some behaviours of the landslides, such as the high mobility of landslide Pastor et al. (2014), the entrainment induced by landslide Liu and He (2016), and the interactions between river flow and submerged landslide Zhao et al. (2017), are also investigated by these models since these behaviours may have appreciable impacts on the process of landslide dam formation. However, most of the existing studies focus on the first and three conditions mentioned above, and works related to the second condition are still rare.

The erosive capacity of the river and the deposition rate of the landslide determine whether a dam can be formed when a landslide reaches a river. The former factor is significantly influenced by river conditions, including flow depth, flow velocity and river slope (Whipple et al., 1998; Carroll et al., 2004; Zhang et al., 2009; Chen et al., 2019a). The latter factor depends mostly on the characteristics, (e.g. mass volume, velocity and material composition) of the landslide (Fan et al., 2014; Pastor et al., 2014; Liu et al., 2016). Although some studies have investigated the mechanism of granular deposition in fluid and its influence on dam formation Shan and Zhao (2014); Zhao et al. (2019); Li et al. (2020), they do not consider the impacts of the erosive capacity of river flow on dam formation. Recently, Chen et al. (2019a) suggested a critical threshold that reflects the influence of river erosive capacity on dam formation in a quantitative way. However, one drawback of this threshold is that it does not consider the dynamic characteristics of the landslide. This factor determines the coverage area and deposition rate of a landslide in a river. In addition, in most field cases, the landslides enter river in an orthogonal (or oblique) direction, which may have an influence on dam formation and needs to be considered.

Because only average river erosive capacity can be observed in the experiment performed by Chen et al. (2019a), there are no way to quantify dynamic processes and dynamic change of different enter river direction (between the river channel and landslide movement direction). Therefore, we need to use numerical simulation to clarify the dynamic process of landslide movement, landslide dam formation and river erosive capacity. At the same time we also observe the dynamic effects of different river entry directions and thus more accurately capture the block point of landslide dam based on river erosive capacity.

In this study, focusing on the second condition, we attempt to determine a critical threshold that incorporates the dynamic characteristics of both landslide and river, in order to reflect the erosive capacity of river flow and its effect on dam formation more reasonably. To describe the dynamics of a landslide and a river during the process of landslide dam formation, a two-layer model based on the depth-averaged theory is used, which incorporates the erosion term between the landslide and the water flow. By analyzing the existing laboratory experimental data and numerical simulation results in combination, a critical threshold value is determined. Finally, the variation in the critical threshold value is discussed by simulating several numerical cases that consider different landslide dynamic characteristics.

Critical Threshold for Landslide Dam Formation

Several works have demonstrated that the formation of a landslide dam is determined by the erosive capacity of the river and the deposition rate of the landslide, both of which are related to the landslide and river discharges, particle diameter and river slope (Yan et al., 2009; Zhao et al., 2017; Chen et al., 2019a). Based on experiments, Chen et al. (2019a) suggest that landslide-generated dams form once the ratio of the erosion rate of river flow to the deposition rate of the landslide exceeds a threshold value. As a matter of fact, these rates represent the discharges of eroded mass and intruding landslide per unit time and therefore, the dimensionless critical discharge can be written in the form of discharge as

q_{*} = \frac{q_{e}}{q_{s} \sqrt{\frac{ρ_{s} - ρ_{f}}{ρ_{f}} g D_{50}^{3}}}, (1)

where q_e is the discharge of mass eroded by river flow with a discharge q_w, q_s is the discharge of landslide intrusion into river, ρ_f is the flow density, ρ_s is the landslide density, D₅₀ is the median grain size of the landslide, and g represents the gravitational acceleration. As mentioned above, the experiments by Chen et al. (2019a) were performed with a constant discharge of sediment into the flow channel. In general, during the process of landslide intrusion into a river, the landslide discharge changes over time in practical cases, as reflected by the variations in the depth and velocity of the sliding mass. This in turn affects the flow erosive capacity, which depends on the interactions between the landslide and the river. Thus, the dimensionless discharge also changes with time and a time-averaged value of this variable may be more suitable for predicting landslide dam formation. To achieve this, instantaneous discharges of the landslide and the river are required, which are relative to their dynamic characteristics.

Physical Model Framework

Governing Equations

The dynamics of the landslide and river are influenced by many different factors, including initial and boundary conditions, material properties, and topography (Chen and Chang 2016). This means that a reliable method of predicting both landslide and river dynamics is needed. From this point of view, a two-layer model that describes the landslide and river dynamics simultaneously has been developed and widely used (Capart and Young, 2002; Chen and Peng, 2006; Liu and He, 2016; Li et al., 2020). Thus, a two-layer model incorporating the erosion term between the landslide and river flow is presented here, following Adduce et al. (2012) and Liu and He (2016). A detailed derivation of the model equations is presented in Supplmentary Appendix A. Since we focus on studying the erosive capacity of river flow and its effect on dam formation, bed entrainment is not considered. By assuming that both layers are incompressible, the mass and momentum equations in a Cartesian coordinate system for each layer can be written as

\begin{array}{l} \frac{\partial h_{1}}{\partial t} + \frac{\partial (h_{1} u_{1})}{\partial x} + \frac{\partial (h_{1} v_{1})}{\partial y} = E, \\ \frac{\partial (h_{1} u_{1})}{\partial t} + \frac{\partial}{\partial x} (h_{1} u_{1}^{2} + \frac{1}{2} g h_{1}^{2}) + \frac{\partial (h_{1} u_{1} v_{1})}{\partial y} = u_{1 m} E - g h_{1} \frac{\partial (z_{b} + h_{2})}{\partial x} - C_{f s} (u_{1} - u_{2}) | u_{1} - u_{2} |, \\ \frac{\partial (h_{1} v_{1})}{\partial t} + \frac{\partial (h_{1} u_{1} v_{1})}{\partial x} + \frac{\partial}{\partial y} (h_{1} v_{1}^{2} + \frac{1}{2} g h_{1}^{2}) = v_{1 m} E - g h_{1} \frac{\partial (z_{b} + h_{2})}{\partial y} - C_{f s} (v_{1} - v_{2}) | u_{1} - u_{2} |, \end{array} (2)

\begin{array}{l} \frac{\partial h_{2}}{\partial t} + \frac{\partial (h_{2} u_{2})}{\partial x} + \frac{\partial (h_{2} v_{2})}{\partial y} = - E \\ \frac{\partial (h_{2} u_{2})}{\partial t} + \frac{\partial}{\partial x} (h_{2} u_{2}^{2} + \frac{1}{2} k_{a p} g h_{2}^{2}) + \frac{\partial (h_{2} u_{2} v_{2})}{\partial y} = - u_{2 m} E - k_{a p} γ g h_{2} \frac{\partial h_{1}}{\partial x} - k_{a p} g h_{2} \frac{\partial z_{b}}{\partial x} + γ C_{f s} (u_{1} - u_{2}) | u_{1} - u_{2} | - \frac{u_{2}}{| u_{2} |} g h_{2} (1 - γ) \tan φ_{b e d}, \\ \frac{\partial (h_{2} v_{2})}{\partial t} + \frac{\partial (h_{2} u_{2} v_{2})}{\partial x} + \frac{\partial}{\partial y} (h_{2} v_{2}^{2} + \frac{1}{2} k_{a p} g h_{2}^{2}) = - v_{2 m} E - k_{a p} γ g h_{2} \frac{\partial h_{1}}{\partial y} - k_{a p} g h_{2} \frac{\partial z_{b}}{\partial y} + γ C_{f s} (v_{1} - v_{2}) | u_{1} - u_{2} | - \frac{v_{2}}{| u_{2} |} g h_{2} (1 - γ) \tan φ_{b e d}, \end{array} (3)

where t is the time; h₁ is the river flow depth; and h₂ is the landslide depth. The two flowing layers, river and landslide, are assumed to have distinct densities ρ_f and ρ_s, with corresponding velocities u₁ = (u₁, v₁) and u₂ = (u₂, v₂), respectively. γ = ρ_f/ρ_s is the density ratio; z_b is the fixed bed surface; E is the erosion rate; k_ap is the earth-pressure coefficient, which reflects the state of stress when a material element deforms Gray et al. (1999); and φ_bed is the basal frictional angle. The term C_fs(u₁–u₂)|u₁–u₂| represents the shear stress at the interface when the landslide moves underneath the water flow, where C_fs = gn²/h₁^1/3Li et al. (2020) and n is the Manning roughness coefficient. In other cases, this term for water flow can be reduced to C_fsu₁|u₁|. u_1m = (u_1m, v_1m) and u_2m = (u_2m, v_2m) are the velocities for the landslide and water flow at the interface boundary, respectively. In simplified situations, as suggested by Adduce et al. (2012), u_(1m,_2m) is considered simply as u_{(1, 2)}.

To close the model, the quantity E must be expressed in terms of variables such as river flow depth, flow velocity and solid density. Currently, researchers have gradually reached an agreement that the erosion rate results from the inequality between the shear stress imparted by water flow and the shear resistance by sediment material (Fraccarollo and Capart, 2002; Spinewine, 2005; Zech et al., 2008). Thus, hydraulic erosion rates are quantified by using a Meyer-Peter and Muller (1948)-style equation that can be empirically fitted to each shear stress-bed load relation (Vericat et al., 2008; Darby et al., 2010).

E = a {(τ_{s} - τ_{b})}^{b} (4)

where a is an erodibility coefficient and b is an empirically derived exponent; the shear stress τ_s can be expressed as τ_s = ρ_fC_fs|u₁–u₂|²; the critical shear stress τ_b is calculated using Shields (1936) equation, τ_b = τ_cD₅₀g (ρ_s–ρ_f) Vericat et al. (2008), in which τ_c is the dimensionless shear stress or Shield’s number modified for sediment materials.

In summary, Eqs. 2, 3 control the state of river flow and landslide, respectively. The first Eqs. 2, 3 represent mass conservation. The second and third Eq 2 represent the momentum conservation in the x and y directions, and the terms on the right-hand side represent the effects of momentum production due to erosion, the gradient induced by the river bed and landslide, and the interface shear stress. Similarly, the terms on the right-hand side of the momentum conservation Eqs. 3 represent the effects of momentum production generated by erosion, buoyancy-related force, gradient induced by the river bed, interface shear stress and friction loss. By coupling (2)–(4), the process of dam formation can be quantitatively described while considering the dynamic characteristics of the river and the landslide. To verify the feasibility of the presented model, the numerical case proposed by Fernández-Nieto et al. (2008) and further used by Kurganov and Miller (2014) is calculated (see Supplmentary Appendix B).

Computational Scheme

In this paper, the Godunov-type scheme based on the finite volume method is adopted to solve the presented model equations. The Godunov-type scheme is a conservative numerical scheme which solves exact or approximate Riemann problems at each inter-cell boundary (Brufau et al., 2004). Here, the Riemann problem at the cell interface is solved by using Harten–Lax–van Leer contact (HLLC) approximation as a robust and efficient solver (Benkhaldoun et al., 2012). In convenience, the model equations can be written in vector format as following:

\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = S + T, (5)

where,

U = (\begin{array}{l} h_{1} \\ h_{1} u_{1} \\ h_{1} v_{1} \\ h_{2} \\ h_{2} u_{2} \\ h_{2} v_{2} \end{array}), F = (\begin{array}{l} h_{1} u_{1} \\ h_{1} u_{1}^{2} + \frac{1}{2} g h_{1}^{2} \\ h_{1} u_{1} v_{1} \\ h_{2} u_{2} \\ h_{2} u_{2}^{2} + \frac{1}{2} k_{a p} g h_{2}^{2} \\ h_{2} u_{2} v_{2} \end{array}), G = (\begin{array}{l} h_{1} v_{1} \\ h_{1} u_{1} v_{1} \\ h_{1} v_{1}^{2} + \frac{1}{2} g h_{1}^{2} \\ h_{2} v_{2} \\ h_{2} u_{2} v_{2} \\ h_{2} v_{2}^{2} + \frac{1}{2} k_{a p} g h_{2}^{2} \end{array}),

S = (\begin{array}{l} E \\ u_{1 m} E - g h_{1} \frac{\partial (z_{b} + h_{2})}{\partial x} - C_{f s} (u_{1} - u_{2}) | u_{1} - u_{2} | \\ 0 \\ - E \\ - u_{2 m} E - k_{a p} γ g h_{2} \frac{\partial h_{1}}{\partial x} - k_{a p} g h_{2} \frac{\partial z_{b}}{\partial x} + γ C_{f s} (u_{1} - u_{2}) | u_{1} - u_{2} | - \frac{u_{2}}{| u_{2} |} g h_{2} (1 - γ) \tan φ_{b e d} \\ 0 \end{array}),

T = (\begin{array}{l} 0 \\ 0 \\ v_{1 m} E - g h_{1} \frac{\partial (z_{b} + h_{2})}{\partial y} - C_{f s} (v_{1} - v_{2}) | u_{1} - u_{2} | \\ 0 \\ 0 \\ - v_{2 m} E - k_{a p} γ g h_{2} \frac{\partial h_{1}}{\partial y} - k_{a p} g h_{2} \frac{\partial z_{b}}{\partial y} + γ C_{f s} (v_{1} - v_{2}) | u_{1} - u_{2} | - \frac{v_{2}}{| u_{2} |} g h_{2} (1 - γ) \tan φ_{b e d} \end{array}) .

A simplest space-splitting type has also been used for dividing the model equations into two 1-D problems as following (Liang and Marche, 2009; Liu and He, 2016).

{\begin{cases} \frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = S \\ \frac{\partial U}{\partial t} + \frac{\partial G}{\partial y} = T \end{cases} (6)

After that, the solution at next time step can be obtained by an efficient step as following:

U^{n + 1} = L_{x} (\frac{d t}{2}) L_{y} (\frac{d t}{2}) L_{y} (\frac{d t}{2}) L_{x} (\frac{d t}{2}) U^{n} (7)

where n represents the time level; L_x and L_y represent the operator in x and y directions, respectively. For L_x, the internal flux, e.g., F_w, is computed as follows:

F_{w} = {\begin{matrix} F_{l} i f 0 \leq S_{l} \\ F_{* l} i f S_{l} \leq 0 \leq S_{m} \\ F_{* r} i f S_{m} \leq 0 \leq S_{r} \\ F_{r} i f S_{r} \leq 0 \end{matrix} (8)

where F_l and F_r are the interface fluxes on both sides of a cell interface; F_*l and F_*r represent the left and right sides of the contact wave, respectively. Both of them are calculated from the left and right Riemann states U_l and U_r. S_l, S_m, and S_r represent the speeds of the left, middle, and right waves, respectively, for a local Riemann problem. The fluxes F_* in the middle region are needed to calculate F_*l and F_*r, which is obtained from the Harten–Lax–van Leer (HLL) formula (Harten et al., 1983).

F_{*} = \frac{S_{r} F_{l} - S_{l} F_{r} + S_{r} S_{l} (U_{r} - U_{l})}{S_{r} - S_{l}} (9)

Considering the dry-bed condition from the two-rarefaction approximate Riemann solver, the wave speeds are calculated as follows (Fraccarollo and Toro, 1995; Soares-Frazão and Zech, 2011).

\begin{array}{l} S_{l} = {\begin{matrix} u_{r} - 2 c_{r} i f h_{l} = 0 \\ \min (u_{l} - c_{l}, u_{*} - c_{*}) i f h_{l} > 0 \end{matrix} \\ S_{r} = {\begin{matrix} u_{l} + 2 c_{l} i f h_{r} = 0 \\ \max (u_{r} + c_{r}, u_{*} + c_{*}) i f h_{r} > 0 \end{matrix} \\ S_{m} = \frac{S_{l} h_{r} (u_{r} - S_{r}) - S_{r} h_{l} (u_{l} - S_{l})}{h_{r} (u_{r} - S_{r}) - h_{l} (u_{l} - S_{l})} \end{array} (10)

where c is the speed of gravity waves; u_l, u_r, h_l, h_r are the components of the left and right Riemann states for a local Riemann problem; u_* and h_* are the components of the middle Riemann states, which are calculated as follows:

u_{*} = \frac{1}{2} (u_{l} + u_{r}) + c_{l} - c_{r}; h_{*} = \frac{1}{g_{z}} {(\frac{1}{2} (c_{l} + c_{r}) + \frac{1}{4} (u_{l} - u_{r}))}^{2}, (11)

In order to obtain high-order of accuracy and avoid spurious oscillations, we couple the monotonic upstream-centered scheme for conservation laws (MUSCL) with HLLC scheme to reconstruct the interface data. The reconstruction form can be expressed as

U_{l} = U_{i}^{n} + \frac{1}{2} (U_{i}^{n} - U_{i - 1}^{n}) \cdot M (q_{i}^{-}); U_{r} = U_{i + 1}^{n} - \frac{1}{2} (U_{i + 1}^{n} - U_{i}^{n}) \cdot M (q_{i}^{+}) (12)

where,

{\begin{cases} q_{i}^{-} = \frac{U_{i + 1}^{n} - U_{i}^{n}}{U_{i}^{n} - U_{i - 1}^{n}} \\ q_{i}^{+} = \frac{U_{i + 2}^{n} - U_{i + 1}^{n}}{U_{i + 1}^{n} - U_{i}^{n}} \end{cases}

The function M is a Roe’ Superbee flux limiter and can be written as

M (x) = m a x (0, m i n (1,2 x), m i n (2, x)) (13)

The time step ∆t that satisfies the demand of two layers dynamic computing simultaneously can be calculated by the stability criterion (Simpson and Castelltort, 2006).

Δ t = \min (\begin{matrix} Δ t_{l a y e r 1}, & Δ t_{l a y e r 2} \end{matrix}); which Δ t_{l a y e r 1 / 2} \leq \min (\frac{c f l \cdot η}{\max (| u_{(1 / 2)} | + \sqrt{g h_{(1,2)}})}) (14)

where cfl is the Courant number and its value should be less than one; η is the ratio of the area of the grid to its perimeter.

Computing Dimensionless Discharge

Landslide and river dynamics and mass exchange between two layers are simulated for each of the cases under different conditions, which allows us to compute the dimensionless discharge as a function of time for every channel location in the basin. Since we are interested in the value of q_* within channel areas at the times when a landslide intrusion enters the river, we compute time-averaged values of q_e and q_s (which allows us to calculate q_*) with the channel network over a length of time from the moment when the landslide reaches the river to the moment when the landslide dam forms. The reasons why we compute time-averaged values of q_* rather than instantaneous values or final values are that the formation of a landslide dam is a gradual process and the erosive capacity of river flow changes over time.

Results

The present model is first applied to simulate laboratory experiments on landslide dam formation over a fixed bed. Then, based on numerical case studies, the value changes in the dimensionless critical discharge for landslide dam formation are presented by considering different dynamic conditions for the landslide.

The setup of the laboratory experiments by Chen et al. (2019a) consisted of two acrylic flumes: one for transporting water and the other for transporting sediment. The water channel was rectangular, 200 cm long, 15.5 cm wide and 20 cm high, with an adjustable slope and a flow valve. The sediment flume was located above the water channel with a longitudinally adjustable gate and was used to supply sediment to simulate a landslide mass entering a river channel. A water reservoir at the upstream end of the flume was used to provide water. The flow valve was attached to the water supply line behind the reservoir. The slope of the water channel was adjusted by attaching a shaft to the upstream end of the flume and a height-adjustable cross-bar at the downstream end. The sediment flume was 180 cm long, 14 cm wide, and 20 cm high, positioned in parallel above the water channel and inclined at a 40° angle, (i.e. the slope angle was always greater than the internal friction angle of the sediment). This arrangement ensured that the deposited sediments were evenly distributed along the width of the water channel. The rate of sediment supply discharge was controlled through an adjustable gate and an acrylic panel at the upstream end of the sediment flume (see Figure 1). Based on experimental results of Chen et al. (2019a), a dimensionless velocity index $v_{d e}$ equation was proposed. If $v_{d e} < 54$ , the landslide mass would block the river channel and form a landslide dam; if $v_{d e} <47$ , the landslide mass would not be able to block the river flow; and if $v_{d e}$ was between 47 and 54, the formation of a landslide dam would be inconclusive. The experimental results showed an 89% accuracy when the dimensionless velocity index ( $v_{d e}$ ) was used to evaluate conditions for which a landslide dam forms. As suggested by Wu and Chou (2003), the empirical value of τ_c for pure sand materials was found to be 0.004 for grain sizes larger than 0.5 mm. The value of q_s for each case was in the Chen et al. (2019a) experiments. Values of other parameters were adjusted in a trial-and-error procedure until empirical adequacy was reached. An overview of the required parameters for the model is shown in Table 1.

FIGURE 1

FIGURE 1. Diagram of the experimental setup on front view (A) and top view (B), modified from Chen et al. (2019a).

TABLE 1

TABLE 1. Model parameter values used in numerical simulations of landslide dam formation.

Figure 2 shows the computed erosion discharge q_e and dimensionless discharge q_* obtained by simulating the processes of 46 experiments in which landslide dams formed in 20 experimental runs and did not form in 26 runs. The dynamic processes of the landslide and river for some cases (1, 3, and 40) are provided in Supplementary Figures S1–S3 (see Supplementary Material). The results show that the mass of the landslide obstructs flow in the river channel and forms a landslide dam when the ratio of q_e and q_s is larger than a certain value (Figure 2B). This trend is similar to the results from Chen et al. (2019a) (Figure 2A), which verifies the feasibility of using our numerical approach to investigate the variations in the dam formation process under different conditions. Simulations indicate that for some cases (especially for the cases that form a dam), there is a large difference between values of q_e calculated by Chen et al. (2019a) and values of q_e calculated by our approach. On the one hand, this contrast is due to the different formula for the erosion rate that we employed. The erosion rate derived by Chen et al. (2019a) is an average value, and that obtained in this study is derived by a process-based method. On the other hand, this difference may result from the fact that the erosion rate gradually decreases to zero during the process of dam formation; thus, a small average value of q_e is obtained. Moreover, the critical threshold of dimensionless discharge q_* as a function of river slope θ is also presented. This threshold clearly shows that q_* is lower for cases where a landslide dam forms than for cases that form no dam (Figures 2C,D). With steeper river slopes, the discharge of the landslide has to be relatively high to form a dam since the erosive capacity of river flow is enhanced by slope.

FIGURE 2

FIGURE 2. Relationship between depositional capacity of landslides q_s and erosive capacity of river flows q_e obtained from Chen et al. (2019a)(A) and numerical simulation (B); critical threshold of dimensionless discharge q_* derived from Chen et al. (2019a)(C) and experimental simulations (D) as a function of river slope θ. Red squares represent the experimental runs with no dam, and gray dots represent the experimental runs with dams. The values of these variables are obtained as an average from the time that the landslide starts entering the river until the time that the landslide dam forms.

Next, the processes of dam formation under different dynamic conditions for the river and landslide are simulated. In general, the dynamic condition of surface flow can be reflected by flow velocity and flow depth; thus, these two variables are considered in our simulation. In addition, field surveys indicate that landslide debris always enters river in an oblique (or orthogonal) direction Tang et al. (2011); Xu et al. (2018), so the angle β between the river channel and landslide movement direction is considered another variable. Based on this consideration, the simulation setups for all the numerical cases are set as h₁ = 1.03 cm, u₁ = 62.8 cm/s, h₂ = 3–25 cm (an interval of one between two adjacent cases), u₂ = 5–35 cm/s (an interval of 15 between two adjacent cases), and θ = 2–9° (an interval of one between two adjacent cases). The sizes of river and landslide channels in the simulation are the same as those in Chen et al. (2019a), except that the two channels are linked at an angle (see Figure 1B). Three values of β are chosen as β = 60, 90, and 120°. For the sake of simplicity, we ran these simulations by using the same model parameter values applied in the prior runs. It should be noted that the value of q_s is not constant in this simulation and calculated at the contact surface between two channels at each time step as well as q_e. The dynamic processes of the landslide and river with different values of β are also provided in Supplementary Material (see Supplementary Figures S4–S6). Using the threshold based on dimensionless discharges (Figure 3), we find that the dynamic characteristics of the landslide have significant impacts on dam formation. In general, typical values of q_* for dam formation are high compared with q_* obtained from flume experiment simulations. With the same angle β, the values of q_* decrease when the landslide velocity increases, and the gap between the values of q_* under different landslide velocity conditions increases as the river slope increases. This means that the erosive capacity of river flow plays a key role in dam formation when the discharge of landslide intrusion into the river is not very large. If the ratio of landslide discharges to river discharges is large, the river channel is blocked quickly, and the effects of flow erosion and river slope on dam formation become less visible. In addition, there is a gap between the values of q_* obtained with different intrusion directions of landslides into rivers. The interesting point is that the value of q_* is largest when β = 90°. One potential explanation for this trend is that the river erosive capacity is enhanced due to the intensified interaction between landslide and river flows when the landslide moves against the river (β = 120°). A larger discharge of landslide intrusion into the river is needed to form a dam; thus, a smaller value of q_* is generated.

FIGURE 3

FIGURE 3. Critical threshold of dimensionless discharge q_* derived from numerical simulations as a function of river slope θ, by considering different dynamic conditions for the river and the landslide. Values of q_* are obtained as an average from the time that the landslide starts entering the river until the time that the landslide dam forms.

Discussion

The results demonstrate that the thresholds for dimensionless discharge obtained in our experimental simulations are smaller than the thresholds for dimensionless discharge suggested by Chen et al. (2019a) (Figures 2C,D). Moreover, there are large differences between the erosion rates derived here and the values suggested by Chen et al. (2019a). As a consequence, the threshold proposed by Chen et al. (2019a) may lead to underestimates of the erosive capacity of river flow required to produce landslide dam under a given hydrodynamic condition. The apparent discrepancy between thresholds derived here and those derived previously may arise from the fact that the previous study focused on a threshold over the whole time period (from landslide entry into river to the end of landslide movement). Whether landslides will continue to be eroded by river flow after dam formation is unclear since they require adequate hydrodynamic conditions. In addition, both our results and those of Chen et al. (2019a) show transitional conditions under which several dam formed and no dams formed (see Figures 2A,B). This is principally because the cases that form no dam have small values of q_e and the cases that form a dam have large values of q_e. Experimental observations indicate that these cases are in (or close to) the critical state of dam formation (see Supplementary Videos S3–S6, which refer to groups 7, 19, 23, and 40, respectively, in Chen et al., 2019a). This means that the erosion rate caused by river flow experiences major changes over the whole time period. Variation histories of erosion rates in these cases calculated by our approach are shown in Figure 4A, which illustrates that the rates of erosion increase first and then decrease (close to zero) and finally increase again. This trend is in line with that observed from the experiments. The changing trend of q_* calculated by our approach is similar with that of q_e (Figure 4B). By comparing the time averaged and Chen et al. (2019a) values of q_* (refer to q_*sa and q_*e), it can be found that q_*sa can reflect the state of dam more accurately. For example, the dam in experiment case seven is more close to the critical state of dam formation than that in experiment case 23, so the value of q_* for experiment case seven should be less than that for experiment case 23. However, the value of q_*e for experiment case seven is greater than that for experiment case 23, which do not reflect this trend. Thus, the dimensionless critical discharge calculated by our approach can reflect the characteristics of dam evolution that is determined by the dynamics of landslide and river, and make it possible to more accurately capture the block point of landslide dam. However, some behaviours, (e.g. local collapses of debris and hydraulic jumps) cannot be considered here due to the limitations of the applied model, and these behaviours may influence the value of the erosion rate. Nevertheless, the proposed approach for generating dimensionless discharge thresholds based on model simulations provides a promising alternative to empirical methods for assessing the potential for dam formation.

FIGURE 4

FIGURE 4. Variation histories of (A) erosion rates and (B) dimensionless critical discharge for the experiment simulation cases 7, 19, 23, and 40. The arrow in the figure represents the critical point of dam formation. q_*sa and q_*e represent the average value of q_* calculated with the data from the simulation and Chen et al. (2019a), respectively.

The choice to average values of dimensionless discharge over a length of time from when the landslide starts entering the river until a dam forms is based on past observations in laboratory experiments that no landslide mass is transported after dam formation (see Supplementary Videos S1). Simulations indicate similar trends between the results obtained by Chen et al. (2019a) and the results obtained from our approach (see Figure 2). Based on the physical model, averaging q_e and q_* over time periods appears reasonable and produces similar results. An additional consideration is whether including landslide and river dynamics in the simulations is necessary to obtain dimensionless thresholds. Our results show that dimensionless thresholds derived from the physical model are larger than those estimated using an empirical formula that cannot consider landslide and river dynamics. The primary reason is that the landslide velocity and the direction of landslide intrusion into river can change the erosion capacity of river flow due to the enhanced or weakened interaction between the landslide and the flowing river. The dimensionless critical threshold in conjunction with numerical simulations potentially reflects the dynamic characteristics of both landslide and river flow. This result is encouraging since landslide and river dynamics may be poorly considered in many applications. An additional benefit of deriving the dimensionless discharge based on model simulation for landslide dam formation is that the erosion rate varies due to the instantaneous states of the landslide and the river. For example, landslide velocity, river flow height, and riverbed elevation change with time following sediment transport (Liu and He, 2016; Zhao et al., 2017; Li et al., 2020). The erosive capacity of river flow therefore changes with time, but thresholds based on dimensionless discharge remain constant. Thus, given data that constrain the initial conditions of the landslide and the river, the variables (q_e and q_s) derived here can be applied to determine how dimensionless discharge changes with time and can be used to identify the landslide-induced terrain changes that have great impacts on dam formation potential. The applicability of the dimensionless discharge threshold, however, assumes that the difference between the discharges of the landslide and the river is not too large.

On the other hand, the presented model assumes that landslide materials are uniform, which is simple if more complex scenarios are considered e.g., different grain size distributions of landslide material. Some complex behaviours that may influence the value of the erosion rate are also not considered. Thus, further research is needed to improve the physical model for providing more accurate results that are closer to reality.

Conclusion

In this study, we derive critical thresholds for the formation of a landslide dam based on slope-dependent values of dimensionless discharge. Furthermore, we present a method for estimating dimensionless discharge thresholds using a process-based two-layer model and the proposed physically based thresholds. The erosion rate and dimensionless discharge derived from the present method indicate trends similar to those estimated by the empirical formula for the experimental cases. The results establish a new method to estimate the thresholds for dam formation focusing on the relationship between the erosive rate of river flow and the deposition rate of a landslide. Several dynamic conditions for the landslide are considered to study their influences on the dimensionless discharge threshold. The physically based dam formation thresholds derived here also make it possible to incorporate the effects of changes in dynamic conditions on the landslide and the river, which could be particularly valuable in addressing landslide dam hazards when landslides and river flow have small differences in discharge.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

WL and YH did a lot of work in numerical data analysis; SH and JZ helped perform the analysis with constructive discussions; WL and KC performed the data analyses, wrote the manuscript, made all figures and approved the final version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos. 41907241, U19A2049, and U20A20111), Open Project of State Key Laboratory of Hydraulics and Mountain River Engineering (Grant No. SKHL 2025), the Soil and Water Conservation Innovation Research Projects in 2021 (SWCB-110–026), the CAS “Light of West China” Program and the Foundation for Young Scientists of the Institute of Mountain Hazards and Environment, CAS (Grant No. SDS-QN-1912, No. SDS-QN-1901).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors wish to thank the Editor (Mark Bebbington) and two anonymous reviewers for their constructive comments which helped to improve the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.651887/full#supplementary-material

References

Adduce, C., Sciortino, G., and Proietti, S. (2012). Gravity Currents Produced by Lock Exchanges: Experiments and Simulations with a Two-Layer Shallow-Water Model with Entrainment. J. Hydraul. Eng. 138 (2), 111–121. doi:10.1061/(asce)hy.1943-7900.0000484