Numerical Simulation on the Longitudinal Breach Process of Landslide Dams Using an Improved Coupled DEM-CFD Method

An accurate investigation of the landslide dam breach process is crucial for the understanding the breach mechanism and disaster prediction. However, the numerical research on the landslide dam breach process to date is rarely reported, especially regarding the soil-water flow coupling effect incorporated in the erosion process. This paper presents a numerical investigation on the longitudinal breach process of landslide dams via a coupled discrete element method (DEM) and computational fluid dynamics (CFD) with the volume of fluid (VOF). Moreover, a virtual sphere model is proposed to overcome the computational instability caused by the particle size approaching the mesh size. The accuracy and validity of the improved coupled method are verified using a series of single particle sedimentation cases. By employing this method, the longitudinal breach process of landslide dams featuring different materials and hydrodynamic conditions has been simulated. It is found to satisfactorily reproduce the longitudinal breach process of landslide dams including surface flow erosion, backward erosion, head-cut erosion, and water and sediment rebalance or complete breach. The effects of the inflow discharges and dam materials on the erosion process are systematically resolved. The breach flow can cause the rotation trend of particles and lead to the increase of tangential contact force at the initial stage of the dam breaching. During the breach process, both the strength and density of the force chain continue to attenuate. The results obtained from the improved coupled DEM-CFD simulations can reasonably explain the particle-fluid interaction mechanisms, physical and morphological evolution and breach process at both macroscopic and mesoscopic scales.


INTRODUCTION
Landslide dams, presented as natural dams formed by landslides blocking the river channels, are widely distributed in mountainous regions all over the world (Schuster and Costa, 1986;Korup, 2004;Nian et al., 2018). The natural lake formed upstream of the landslide dam is called the barrier lake. As a natural dam without a spillway to control the reservoir capacity of the barrier lake, it has a high risk of failure that may cause outburst floods with potentially catastrophic consequences in downstream areas (Casagli et al., 2003;Chang et al., 2011;Wu et al., 2020). For example, the landslide dam on the Dadu River in 1786 broke 10 days after forming and caused severe floods that killed more than 100,000 people (Dai et al., 2006). Recently, a huge landslide dam (called Baige landslide dam) was formed on the Jinsha River in China on October 11, 2018, and created a barrier lake with a maximum reservoir capacity of more than 200 million m 3 (Ouyang et al., 2019). Two days later, the breach of the Baige landslide dam occurred, affecting more than 100,000 people and causing up to 15 billion economic losses. Therefore, there is an urgent need to deepen the understanding the landslide dam breach mechanism that is essential for breaching prediction and disaster assessment .
The research on the breach mechanism of the landslide dam has received considerable attention in recent years. Model experiments in the laboratory and field are currently the most common approaches for investigating the breach process of landslide dams. Many researchers carried out a series of experiments to explore the breach process of landslide dams under different geomorphological conditions, hydrodynamic conditions and dam geometries mainly focusing on the breach flow discharge, evolution of breach morphology, and erosion characteristics (Davies et al., 2007;Cao et al., 2011;Chen et al., 2015;Xiangang et al., 2018;Wang et al., 2018;Jiang, 2019;Zhou et al., 2019). The breach mechanism of landslide dams was discussed based on macroscopic observation, and the primary factors affecting the landslide dam breach process were identified (such as inflow discharge, barrier lake volume, dam material and downstream slope angle). While experiment investigations of the landslide dam breach mechanism are still challenging due to the limitations of the experimental conditions and measurement methods, and thus may not be able to fully reveal the microscopic characteristics (Cao et al., 2011;Wang et al., 2018).
The mathematical models are currently another kind of methods to investigate the landslide dam breach issue. In general, they can be categorized as the parametric model and physically based model (ASCE/EWRI Task Committee on Dam/Levee Breaching, 2011). The parametric models are mainly based on the statistical data of landslide dam failure events, and breach parameters are estimated from empirical formulas without any breach process calculations (Zhong et al., 2016). Besides, these models can only provide discrete characteristic values such as failure time and peak breach flow (Xu et al., 2009;Peng and Zhang, 2012;Sattar, 2014). The physically based models are developed based on the hydrodynamic and erosion process during the landslide dam breach that is widely used in practice (Chang and Zhang, 2010;Zhong et al., 2018Zhong et al., , 2020Zhang et al., 2019). However, these models bear a considerable degree of uncertainty as the longitudinal evolution mode and erosion development process are simplified and pre-defined (Wang et al., 2015;Zhong et al., 2020). Particularly, the landslide dam breach involves complicated interactive processes between breach flow, breach morphology and the multi-scale particle transport, the mathematical models can not enough describe the soil-water flow coupling effect incorporated in the erosion process. To date, the understanding of the breach mechanism has still remained poor (Bovis and Jakob, 2000;Cao et al., 2011;Jiang et al., 2020).
This article aims to provide new insights of the breach mechanism of landslide dams through numerical simulation. To achieve this goal, a coupled discrete element method (DEM) and computational fluid dynamics (CFD) with the volume of fluid (VOF) was used to model the longitudinal breach process of landslide dams, and a virtual sphere was proposed to overcome the computational convergence caused by the limitation of particle and mesh cell size. Then, the dam morphological evolution, force evolution and the influence of hydrodynamic conditions and dam materials on the breach process are investigated systematically, and the longitudinal breach mechanism of landslide dams is discussed at both macroscopic and mesoscopic scales. The erosion process at the dam crest as well as the capability of the improved coupled DEM-CFD modeling of landslide dam breach is also discussed.

IMPROVED COUPLED DEM-CFD METHOD FOR LANDSLIDE DAM BREACH SIMULATION
Currently, the coupled methods of continuous and discontinuous seem to facilitate modeling the interaction of soil and water flow. In this regard, the coupled DEM and CFD numerical method has already been used to investigate the underwater sedimentation of the landslide dam and seepage failure mechanism of dam material (Zhao et al., 2017;Shi et al., 2018;Nian et al., 2021). However, it is noted that the traditional DEM-CFD model is difficult to be directly applied in the simulation of the landslide dam breach process due to the shortcoming in describing the evolution of the free river surface. Besides, the critical size ratio of particle to mesh cell also limits the modeling of coarse particles in dam material and the accurate calculation of the fluid domain.
In this section, an improved coupled DEM-CFD method is developed to simulate the longitudinal breach process of landslide dams, and a virtual sphere model is proposed as a new local porosity calculation method to overcome the limitation of the particle and mesh cell size. The granular landslide dam material is modeled by the DEM, the fluid domain is described by the Navier-Stokes equation to be solved by the CFD, and the free fluid surface is simulated by the VOF model by taking into account the presence of the interface between fluid phases. The derivation of the governing equations, calculation of coupling forces, and realization of the virtual sphere model are described in detail in the following sections.
Governing Equations for the Improved Coupled DEM-CFD Method For the particle in the fluid, the governing equations for the translational and rotational motions of a single particle are determined according to Newton's second law and can be expressed as (Cundall and Strack, 1980): Where m i is the mass of particle i, v i , and ω i are the translation and angular velocities of particle i, respectively. F ij is the contact force between particle i and particle j, respectively. F f is the particle-fluid interaction force acting on particle i. r c is the distance from the center of the particle to the contact point and M ij is the rolling resistant moment. The Hertz-Mindlin contact model is adopted to describe the interaction behavior in normal and tangential directions (Utili et al., 2015;Li et al., 2020): Where F n and F d n are the normal force and the normal damping force, while F t and F d t are the tangential force and the tangential damping force, δ n and δ t are the normal overlap and the tangential overlap, v − → rel n and v − → rel t are the normal component and tangential component of the relative velocity, E * , R * , and m * are the equivalent Young's Modulus, the equivalent radius and the equivalent mass, respectively. S n is the normal stiffness and S t is the tangential stiffness. µ r is the coefficient of rolling friction and R i is the distance of the contact point from the center of mass.
The governing equation of fluid flow considering the presence of particles is as follows (Anderson and Jackson, 1967): Where u is the fluid velocity, ρ f is the density of the fluid, µ f is the fluid viscosity, p is the fluid pressure, g is the gravity acceleration. The terms on the right side of Equation (9) represent the pressure gradient, stress, gravity, and particlefluid interaction force, respectively. f pf in Equation (9) can be related to F f in Equation (1) withf pf = − n i=1 F f /V cell , with V cell being the volume of a CFD mesh cell. ε is the local porosity that is used to represent the influence of particles on the fluid calculations. F s = σκ∇α 1 is the surface tension force, with σ being the surface tension and κ being the local curvature at the interface.
In this research, based on the multiphase fluid theory in the CFD, the VOF model is introduced to describe the dynamic behaviors of the barrier lake and breach flow with free fluid surface evolution (Hirt and Nichols, 1981). In each mesh cell, the volume fraction of all the fluid phases is ε, for which, the volume fraction of the primary phase (e.g., water) is α 1 , while the volume fraction of the secondary phase (e.g., air) is α 2 (α 2 =1 − α 1 ). Hence, for the fluid phase is fully occupied by water, α 1 =1, for the fluid phase is full of air, α 1 = 0. The case of 0 < α 1 <1 normally refers to the presence of a free fluid surface. The free surface motion between the two phases can be tracked by solving the continuity equation of the secondary phase volume fraction: At the interface of air and water in the VOF model, the characteristics of fluid density ρ f and viscosity µ f in Equation (9) are derived from the characteristic of each phase by linear interpolation with the volume fraction:

Particle-Fluid Interaction
The interaction force between particles and fluid mainly includes the buoyant force, drag force, lift force, and virtual mass force (Zhu et al., 2007). The buoyant force F b , drag force F d , and lift force F l are generally considered significant in the sediment entrainment problem Ibrahim and Meguid, 2020) and are introduced to describe the particle-fluid interactions in this research. The buoyant force F b acting on the particle under consideration can be calculated by: Where V p is the volume of the particle. The drag force is caused by the viscous shearing effect of fluid, and is induced by the relative motion between particle and fluid. The Di Felice drag force model (Di Felice, 1994) is used herein to define the drag force: Where C d = 0.63 + 4.8 √ Re p is the drag force coefficient, d is the diameter of the particle, v is the velocity of the particle, is the corrective coefficient and Re p = εdρ f |u−v| µ f is the particle Reynolds number. The lift force accounts for the rotational movement of particles, and the Saffman lift force F lS (Saffman, 1965) and Magnus lift force F lM (Rubinow and Keller, 1961) are both considered in this research. Equations to calculate these forces are given as follow:

Calculation of Local Porosity for Coarse Dam Material
In the locally averaged DEM-CFD method, the local porosity ε plays a crucial role in the solution of fluid governing equations and the calculation of particle-fluid interaction. The calculation of porosity in previous researches is to directly calculate the volume occupied by particles in a single fluid mesh cell [in the central model, the particle volume is assigned to the fluid grid where the center of mass is located, while in the divided model, the particle volume is divided according to the actual position of the particle (Kloss et al., 2012;Zhao and Shan, 2013)]. In either way, the fluid domain would be discontinuous and cause numerically incorrect results when the particle size approaches the minimum fluid mesh cell size. In 3D coupled models, the size ratio of the fluid cell to particle diameter was suggested keeping not less than 4 to keep the accuracy of the interaction calculation (Zhao et al., 2014). For the simulation of landslide dam material with non-negligible coarse particles, a larger mesh cell size would force to be used that is contrary to the requirement of accurate fluid domain calculation.
To overcome the limitation of the size ratio of the fluid cell to dam material, Shi et al. (2018) adopted a method of cluster element generation to replace coarse particles with fine particle clusters. But this would lead to volume errors of the replaced large particles. In this research, based on the porous cube method (Link et al., 2005), an alternative way to calculate local porosity, the virtual sphere model, is presented. As shown in Figure 1, for coarse particles with a diameter greater than a quarter of the mesh cell size (represented by a solid circle), the local porosity calculation is based on the larger mesh cell range (represented by colored cells) defined by the virtual sphere (represented by dotted circles). The virtual sphere has the same centroid as the real particle, and the diameter is 4 times the real particle diameter. The volume of the real particles is equally divided into the related mesh cells involved in the local porosity calculation: Where ϕ VS,i is the volume contributed by particle i to a related mesh cell participating in the local porosity calculation, V p is the actual particle volume and V VS is the volume of all related mesh cells involved in the local porosity calculation. The mesh cells related to different virtual spheres are allowed to overlap, so the porosity of one certain mesh cell involved in the calculation should be:

Benchmarking Example
The sedimentation of a single particle from air to water has been simulated to verify the accuracy and validity of the present improved DEM-CFD program. A fluid domain of 0.1 m × 0.1 m × 0.2 m is set, with the upper and lower parts filled with air and water, respectively. The motionless particle with a diameter of 1    mm is placed at a distance of 0.05 m above the water surface. The main simulation parameters are shown in Table 1. The movement of the particle can be expressed as: The fluid domain in the simulation is meshed by equilateral hexahedrons with an equal size of 4 times the particle diameter.
In Figure 2, the time series of particle sedimentation velocity from numerical results are compared with the analytical solutions. The numerical results are in good agreement with the analytical solutions, and the particle velocities are both stabilized at 0.134 m/s. Furthermore, the single-particle sedimentation simulations with different size ratios of the fluid cell to particle diameter have been carried out to verify the validity of the coupling force calculation by the virtual sphere model. The particle diameter is kept constant, and the mesh size varies between 1 and 5 mm. The simulation is carried out simultaneously in the central model and virtual sphere model. The simulation results are shown in Figure 3. In the central model, as the size ratio continues to decrease, the error of the terminal velocity keeps increasing. This result reflects the critical size effect of the locally averaged DEM-CFD method, that is, keeping the size ratio above 4 would obtain more accurate results. On the contrary, the simulation results of the virtual sphere model under different mesh sizes are relatively stable. Compared with the central model, the accuracy of terminal velocity is improved by 40% maximum with the virtual sphere model. It indicates that the virtual sphere model effectively eliminates the limitation on the size ratio of the mesh cell to particle diameter in the local porosity calculation.

MODEL CONFIGURATIONS
Both the geometric configuration of the landslide dam and the reservoir capacity of the barrier lake play an important role in the evolution of the dam breaching, which should be considered when designs the landslide dam model configuration . Peng and Zhang (2012) proposed three dimensionless FIGURE 5 | Grain size distribution of the simulation materials. The light shaded region is the cumulative grain size distribution of 14 landslide dams in Southwest China , the gray shaded region is the cumulative grain size distribution of Yangjiagou landslide dam (Li et al., 2021), and the solid lines represents the grain size distribution of the modeled landslide dams.   The schematic diagram of model configurations is illustrated in Figure 4. The simulation was conducted in a flume measuring 7 m long, 0.5 m wide, and 0.5 m deep, and the inclination angle of the river bed was 5 • . The landslide dam was placed 3 m away from the inlet. The height of the landslide dam was 0.3 m, the top width was 0.3 m, the angle between the upstream slope and the river bed was 35 • , and the angle between the downstream slope and the river bed was 25 • . An "idealized" geometric configuration of the landslide dam across the river channel was used in this research, which would still very meaningful to advance the understanding of the longitudinal breach mechanism of landslide dams (Cao et al., 2011 were 0.18, 2.2, and 2.0, respectively, and were all conformed to the distribution range of dimensionless parameters given by Zhou et al. (2019). The upstream and downstream boundaries were set as inlet and outlet, and the top boundary was set as an open-air condition. In this research, the upstream inlet flow discharge was changing to study the influence of different hydrodynamic conditions on the dam breaching process. The ratio of the cube root of the inflow discharge per unit time to the dam height (Q in ·T u ) 1/3 H d was defined as the hydrodynamic coefficient, where Qin is the inflow discharge and Tu is the time scale. The dimensionless coefficient was evaluated with the data of 70 landslide dams reported by Shen et al. (2020), and it was found that the inflow discharge distributes in a relatively wide range from 0.007 to 1.176, T u = 1s). According to the evaluation results of the dimensionless coefficient, it was determined that the inflow discharges in this research are 0.005 m 3 /s, 0.01 m 3 /s, 0.015 m 3 /s, and 0.02 m 3 /s ( (Q in ·T u ) 1/3 H d range from 0.57 to 0.9) to represent the inflow from low to high. Figure 5 shows the grain size distribution of Yangjiagou landslide dam (Li et al., 2021) and 14 landslide dams in Southwest China .Considering the grain size effect and the computational limitations, the particles with diameters larger than 50 mm and less than 1 mm were removed. Two types of non-cohesive landslide dam materials were set and shown in Figure 5, one for unsize coarse-grained granular material (consisting of gravels, represented by C) and one for wide size range fine-grained granular material (consisting of sands and gravels, represented by F). The serial numbers for the unsize coarse-grained granular material specimens with these four inflow discharges were C1, C2, C3, and C4, and those for the wide size range fine-grained granular material specimens were F1, F2, F3, and F4, respectively, as summarized in Table 2. The fluid domain was discretized by fixed-grid parallelepiped cells with a size of 50 mm. The parameters used in the model are given in Table 3. The values of these DEM parameters have been calibrated by Li et al. (2020). A total of 8 sets of simulations were carried out (see Table 2), and the dam breaching process as well as the influence of different characteristics of dam materials and inflow discharges were discussed in the following sections.

NUMERICAL SIMULATION OF THE LANDSLIDE DAM BREACH PROCESS General Description
Although each scenario has its own characteristics, they all share important universal characteristics. In general, with the continuous inflow from upstream, the water level of the barrier lake gradually increases. During this process, due to the fragmented dam materials with large pores inside, the seepage flow would first appear at the foot of the downstream slope. Once the water level overtops the dam crest, dam breaching commences through overtopping erosion . Figure 6 shows the typical temporal sequences of the landslide dam breach process (scenario C2). The numerical results indicate that the longitudinal breach of landslide dams experiences four phases: surface flow erosion, backward erosion, head-cut erosion, and water and sediment rebalance or complete breach. The characteristics of different breaching phases are considered to describe the dam breaching process, as follows: (a) Surface flow erosion. The dam breaching initially begins at the downstream slope of the dam, and the scouring of the surface flow is the key factor for the particle activation on the downstream surface . During this phase, it can be observed that the water flow on the slope is accelerated by gravity, the flow velocity at the foot of the downstream slope increases significantly, and the erosion ability of the water flow is enhanced ( Figure 6A). The failure of the landslide dam starts from the toe of the downstream slope, and the activation of fine particles is earlier than that of coarse particles. (b) Backward erosion. With the further development of erosion, the breach of the dam body started to transform from the failure of the particles at the downstream slope to the activation of a certain scale dam material. In this phase, the hydraulic-gravity coupling [consisting of the scouring and seepage of the breach flow as well as the gravity potential energy of the dam material (Dang et al., 2008)] is the main factor for the activation of the dam. The backward erosion of the landslide dam without initial breach develops from the toe to the crest of the downstream slope, accompanied by the continuous slowing of the downstream dam slope angle (Figures 6B,C), which is in agreement with the experiment results from Cao et al. (2011). (c) Head-cut erosion. As the backward erosion continues to develop and reaches the upstream dam crest, the headcut erosion begins and causes the dam height to decrease significantly. In this phase, the breach flow discharge FIGURE 7 | Time series of normal contact force, tangential contact force and particle-fluid interaction force. The black axis, red axis, and blue axis on the y-axis represent the average normal contact force Fn, average tangential contact force Ft and total particle-fluid interaction force Ff, respectively. increases rapidly due to the decrease of the dam height, and the increase in the flow velocity at the dam crest could be significantly observed in Figure 6D. The increase of the flow velocity at the dam crest would further strengthen the erosion capacity of the breach flow and in turn accelerates the development of head-cut erosion. (d) Water and sediment rebalance or complete breach. As the water level of the barrier lake continues to drop, the outflow discharge is gradually decreased, and the erosion capacity of the breach flow continues to weaken. The coarse particles on the downstream slope accumulate and form a coarse layer as the dam breaching progresses, which hinders the further erosion of the lower particles. In the case of low inflow discharge, the inflow and the breach flow could reach a new balance, and the dam no longer undergoes significant erosion damage. In the case of high inflow discharge, the dam would fail completely.

Force Evolution During Landslide Dam Breach
The time series of normal contact force, tangential contact force, and particle-fluid interaction force are shown in Figure 7. As shown in Figure 7, on the one hand, the breach flow reduces the normal contact force between the particles through the buoyancy force, on the other hand, the breach flow causes the rotation trend of particles through the seepage and scouring, resulting in an increase of tangential contact force between particles. About 50 s, the backward erosion has progressed to the upstream dam crest, and this marks the beginning of the head-cut erosion phase. The erosion capacity of the breach flow is enhanced due to the loss of the dam height. This would cause the normal contact force and the tangential contact force to decrease rapidly. For the particlefluid interaction force, there will be a certain increase in the initial breaching stage due to the appearance of the overtopping flow. As the breach progresses and the storage capacity decreases, the particle-fluid interaction force continues to decrease. The total energy E of the particle system is defined as the sum of the kinetic energy, rotational kinetic energy and potential energy, which can comprehensively reflect the energy evolution in the dam breaching process. The dimensionless total energy [E] is normalized by the total energy at t divided by the initial total energy at rest state: The evolution process of the dimensionless total energy of each scenario is shown in Figure 8. As the breach progresses, the dam material is activated by the breach flow and propagates downstream, accompanied by the transformation of potential energy to kinetic energy. The activated particles are continuously lost due to the carrying of the breach flow, resulting in the continuous reduction of the dam volume and the continuous attenuation of the total energy. With the formation of the coarse layer and the continuous decrease of the storage capacity, the breach of the dam tends to stagnate, and the attenuation of the total energy gradually stabilizes. In this process, the increase in the inflow discharge can significantly increase the initial decay rate of the total energy and reduce the residual value of the total energy, which also marks the intensification of the breach and the continuous decrease of the residual dam volume. With the same inflow discharge, the residual total energy of the wide size range fine-grained granular dam is always lower than that of the unsize coarse-grained granular dam. This is because fine particles are easier to activate and aggravate the dam breach. Figure 9 shows the time series of breach morphology evolution processes of landslide dams along the longitudinal direction with different inflow discharge and dam materials. The morphology curve gradually changes from sparse to dense with time, and finally basically no longer changes, indicating the end of the failure process. In the case of 0.005 m 3 /s inflow discharge and unsize coarse-grained dam material, the overtopping flow does not appear. It can be seen in Figure 9A that when the deformation progresses to the middle of the dam crest, the dam failure stops. In the case of low inflow discharge and coarse granular landslide dam with high porosity, the seepage effect is responsible for the main drainage function and may lead to substantial dam subsidence . This kind of subsidence is similar to backward erosion and the failure point develops from the downstream slope to the dam crest. The dam material composed of fine particles can reduce the void ratio and improve the impermeability of the dam (Dhungana and Wang, 2020). Compared with Figure 9A, the dam in Figure 9E has experienced a complete overtopping breach process. As shown in Figures 9B-H, for scenarios with overtoppinginduced dam breaching, the head-cut erosion is the main phase of the dam breaching process. With the increase of inflow discharge, the duration of the surface flow erosion phase and the backward erosion phase continues to decrease. The increase in the inflow can significantly increase the erosion capacity of the breach flow so that the erosion point quickly develops to the upstream dam crest. With the evolution of the failure process, the dam morphology gradually transforms into a triangle, and the downstream slope gradually slows down. It can be observed from Figure 9 that the inflow discharge and the dam material can significantly affect the downstream slope angle. With the increase of the inflow discharge, the carrying capacity of the breach flow increases, so that the dam material can be carried to a farther place and the downstream slope is slowed down. Also, the erosion is directly related to the particle size, that is, the erosion resistance of fine particles is much lower than that of coarse particles. Therefore, under the same inflow conditions, the head-cut effect of the unsize coarse-grained dam material is more significant than that of the wide size range fine-grained dam material, and the downstream dam slope angle is also significantly slower than that of the wide size range finegrained dam material. Figure 10 shows the evolution process of the landslide dam height and the residual dam height of each scenario. It can be seen that one of the typical effects of inflow discharge on the landslide dam breach is to reduce the residual height of the dam. With the increase of the inflow discharge, the erosion capacity of the breach flow is continuously enhanced and causes more dam body failure. Under the same inflow conditions, the residual dam height of the wide size range fine-grained dam material is always lower than that of the unsize coarse-grained dam material, which indicates that the existence of fine particles reduces the erosion resistance of the landslide dam. On the one hand, this phenomenon is due to the existence of fine particles that reduces the porosity inside the dam body and reduces the drainage capacity of seepage; on the other hand, the erosion resistance of fine particles is lower than that of coarse particles. This also explains the phenomenon that there are still residues in scenario C4 and the dam body has completely failed in scenario F4.

Erosion Process at the Dam Crest
To further research the breach process of the landslide dam, the dam height reduction rate ε is defined as the loss of dam height per unit time, and is used to quantify the erosion rate at the dam crest: The time series of dam height reduction rates in each scenario is shown in Figure 11. It is shown that the dam height reduction rate fluctuates with time. In fact, the experiment proposed by Jiang et al. (2016) also confirmed the multi-peak phenomenon of the erosion rate, which supports the reliability of our simulation. The incipient velocity of fine particles is much lower than that of coarse particles. Therefore, the fine particles are taken away first, while the coarse particles are more difficult to be eroded. The failure at the dam crest is dominated by the small-scale collapse . This is mainly caused by the loss of fine particles and the erosion of the breach flow and would lead to the rapid increase of the dam height reduction rate in a short time.
With the increase of the inflow discharge, the local activation of the dam crest becomes more frequent and more intense, which makes the dam height reduction rate fluctuate more intensely. Besides, with the same inflow discharge, the height reduction rate of wide size range fine-grained granular dams fluctuates more frequently than unsize coarse-grained granular dams, but the peak point of height reduction rate of unsize coarse-grained granular dams is higher than that of wide size range fine-grained granular material dams. This indicates that the activation scale of the unsize coarse-grained granular dam is larger than that of the wide size range fine-grained granular dam, and the number of activation is less than that of the wide size range fine-grained granular dam. The fluctuations in each scenario continue to decay over time. This is due to the continuous decrease of the storage capacity of the barrier lake and the gradual formation of the coarse layer. It is increasingly difficult for the breach flow to cause activation of the dam. For scenario F4, a significant rise in the dam height reduction rate is observed at the end of the breach, which is due to the overall failure of the residual dam.

Application and Limitation
Apparently, the landslide dam breach is a coupling sequence of breach flow and landslide dam (Zhong et al., 2018). In this study, the proposed improved coupled DEM-CFD method is innovatively used to simulate the longitudinal breach process of landslide dams. The advantage of DEM is allowed to model loose landslide dam material through discontinuous particles. The improved coupled DEM-CFD method cannot only fully consider the soil-water coupling mechanism, but also reproduce the development process of breach flow with free water surface evolution, which is conducive to understanding the landslide dam breach mechanism. Although dimensionless parameters have been used to guide the design of the model configuration, the landslide dam breach mechanism from this research is inevitably constrained by the spatial and time scales. This is mainly due to the limitation of current computational efficiency, which also leads to the landslide dam materials used in this research are still distributed in a relatively narrow range. A landslide dam on-site can contain a wide range of sediment sizes, from clay, sand to gravel (Costa and Schuster, 1988;Korup, 2002), and would lead to a dramatic increase of the particle number in the DEM modeling. In this research, the influence of cohesive clay on the dam breaching process is ignored, so the results inevitably bear uncertainty. Despite these uncertainties and limitations, the landslide dam breach mechanism as illustrated in this research still offers valuable insights for the understanding of disaster evolution. Further research on this topic will be the improvement of calculation speed by Message Passing Interface (MPI) parallelization and comprehensive analyses of the landslide dam breach process based on real cases.

CONCLUSION
In this work, the improved coupled DEM-CFD method is developed to simulate the landslide dam breach process. Besides, a virtual sphere model is proposed to overcome the computational instability caused by the particle size close to the mesh size. The longitudinal breach process of the landslide dam is simulated with different dam materials and inflow discharges, the breaching phases of landslide dams are proposed according to the dam morphological evolution, and the breach characteristics are revealed from both macroscopic and microcosmic scales.
The simulation shows that the landslide dam breach process can be identified as four phases: surface flow erosion, backward erosion, head-cut erosion, and water and sediment rebalance or complete breach. This is different from the classification based on hydrological characteristics , but it can better reflect the characteristics of the longitudinal morphological evolution of the landslide dam. The normal contact force, tangential contact force and particle-fluid interaction force during the breach process have been evaluated. During the surface flow erosion, the particle-fluid interaction force increases due to the appearance of the overtopping flow. In addition, the rotation trend of particles caused by the breach flow leads to an increase in the normal contact force. The head-cut erosion is the main breach phase of the landslide dam, and the obvious increases in the breach flow velocity at the dam crest is captured in this stage. This can accelerate the dam breach and cause a rapid decrease in the normal contact force, tangential contact force and particle-fluid interaction force. The failure at the dam crest in the head-cut erosion phase is dominated by small-scale collapse and leads to the dam height reduction rate fluctuates with time. From the view of energy, the total energy of the landslide dam continues to attenuate with the continuous loss of dam material during the breach process. The inflow discharge can significantly increase the erosion capacity of the breach flow, thereby accelerating the dam breaching and significantly reducing the residual landslide dam height. The presence of fine particles can slow down the downstream slope during the breach process and reduce the residual dam height.

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/s.