Physically based modeling of coseismic landslide, debris flow, and flood cascade
 ^{1}Faculty of GeoInformation Science and Earth Observation (ITC), University of Twente, Enschede, the Netherlands
 ^{2}Institute of Mountain Hazards and Environment, Chinese Academy of Sciences & Ministry of Water Conservancy, Chengdu, China
 ^{1}Faculty of GeoInformation Science and Earth Observation (ITC), University of Twente, Enschede, the Netherlands
 ^{2}Institute of Mountain Hazards and Environment, Chinese Academy of Sciences & Ministry of Water Conservancy, Chengdu, China
Correspondence: Bastian van den Bout (b.vandenbout@utwente.nl)
Hide author detailsCorrespondence: Bastian van den Bout (b.vandenbout@utwente.nl)
The 2008 Wenchuan earthquake lead to various complex multihazard chains that included seismically triggered landslide initiation, landslide runout, river damming, dam breaching, and flooding. The modeling of the interactions between such hazardous processes is challenging due to the complexity and uncertainty. Here we present an eventbased physically based model that is able to simulate multihazard land surface process chains within a single unified simulation. The final model is used to simulate a multihazard chain event in the Hongchun watershed, where coseismic landslides led to a landslide dam and, 2 years later, a debris flow that breached the landslide dam. While most aspects of the multihazard chain are predicted well, the correct prediction of slope failures remains the biggest challenge. Although the results should be treated carefully, the development of such a model provides a significant progress in the applicability of multihazard chain simulations.
Many of the damaging events that involve land surface processes are not caused by individual but rather multiple interacting hazardous processes. Such combinations can take place because the same triggering events (e.g., extreme rainfall) trigger various hazardous processes (e.g., flash floods, landslides, and debris flows) that interact and that may impact the same elements at risk. For multihazard assessment, the hazard intensities and impact can differ significantly when compared to the individual hazardous processes (Gill and Malamud, 2014; Van Westen and Greiving, 2017). Hazardous events may also occur in sequence as cascading or domino events whereby one hazardous process triggers another either directly or later in time. One particular hazardous example of such cascading hazard events is the natural damming of rivers by landslides (Costa and Schuster, 1988; Walder and O'Conner, 1997). These landslide dams occur in mountainous landslideprone areas around the world (Swanson et al., 1986; Chai et al., 2000; Dai et al., 2005; Korup, 2005; Harp and Crone, 2006; Nash et al., 2008; Liu et al., 2009; Fan et al., 2012). This sequence of hazardous events starts with a catastrophic slope failure, often caused by intense precipitation, a seismic trigger, or a combination of these. The released material moves downslope and can lose momentum through friction or by entering a sharply incised channel or river. When the volume of the solid materials is large enough, or when the moved landslide body retains mostly its original strength characteristics, it will form a barrier for water flow which will accumulate to form a dammed lake. Depending on the strength of the materials composing the dam, the barrier may breach when the barrier lake level exceeds the height of the dam due to either accelerated erosion, piping, or barrier collapse (Ermini and Casagli, 2003). During the breach, extreme discharges and solidladen floods with high velocity may devastate the downstream areas (Schuster, 1993; Walder and O'Connor, 1997).
Natural dynamical systems such as the one described above can be complicated, containing many interactions and numerous fundamental processes related to hydrology and sediments (Walder and O'Conner, 1997). Slope failures, mass movement runout, and flooding are influenced by catchmentscale hydrology (Bout et al., 2018). Furthermore, interhazard interactions exist in many varieties, a review of which can be found in Kappes et al. (2012). In the case of landslide dam formation and breaching, many interactions exist between processes that are typically approached individually in modeling. Landslide dambreak floods have been analyzed using both empirical and physically based models (Evans, 1986; Costa and Schuster, 1988; Peng and Zhang, 2012). Empirical models are simpler to apply but provide less comprehensive results than physically based models. On the other hand, physically based models require detailed physical parameters as input and can be computationally costly. Whereas the individual components of the hazard chain – landslide initiation, landslide runout, dam breach, and flooding – have been modeled using physically based models, the interactions between these processes are generally not simulated within a single model due to their high complexity.
The prediction of landslide volumes resulting from earthquakes is a complex problem and requires specialized numerical models. Several physically based simulation tools for slope failure volume modeling have been developed, such as CLARA (Hungr et al., 1989), TSLOPE3 (Pyke, 1991), 3DSLOPE (Lam and Fredlund, 1993), 3DSLOPEGIS (Xie et al., 2003), r.slope.stability (Mergili et al., 2014a), Scoops3D (Reid et al., 2015), EDDA (Chen and Zhang, 2015), and OpenLISEM (Bout et al., 2018). CLARA, TSLOPE3, and 3DSLOPE can only be applied on individual slopes, while r.slope.stability, Scoops3D, and OpenLISEM are spatially distributed models, which are based on geographic information systems (GISs). These can be applied for landslide volume estimation over a large area up to several hundred square kilometers. Numerical modeling of mass movement runout using twodimensional approaches has been implemented in a variety of models (Malet et al., 2004; Rickenmann et al., 2006; Van Asch et al., 2007; Hürlimann et al., 2007; Domenech et al., 2019). They need detailed information on initial volume, rheology, and entrainment, as well as an accurate and detailed digital elevation model (DEM) (Hürlimann et al., 2007). Erosion, the waterdriven uptake of sediment, and entrainment, the graindriven uptake of sediment, have been used in understanding mass flow soil interactions. Erosion models, while traditionally focused on agricultural processes, come in great variety and provide insight into the flow–surface interactions. Examples are WEPP (Flanagan et al., 2001), EUROSEM (Morgan et al., 1998), and Delft3D Sediment (Roelvink and van Banning, 1995).
In the case of multihazard chains including landslide dams, integrated simulations face several critical issues. When a mass movement enters a channel with a certain water level, the landslide material mixes with the water in a dynamic manner, which is generally ignored in existing flow models such as Flo2D (O'Brien et al., 1993). Similarly, the volumetric sediment content of water increases when a landslide dam is breached, and the material of the dam is entrained by the water flow. The entrainment of bed material is simulated in a limited number of spatial twodimensional mass movement models but is rarely initiated from lowconcentration water flow (Hungr et al., 2005). Several models have shown this functionality but lack the capability of modeling breaching behavior, and they ignore the resulting changes in the digital elevation model (Chen and Zhang, 2015; Hu et al.,2016; Baggio et al., 2021). Ignoring these changes makes the simulation of any breaching behavior impossible since increasing outflow must be the result of the entrainment of a flow path on the landslide dam.
Simplified coupled model approaches have been tested using separate and nonspatially distributed models for water flow and dam breaching. Empirical equations for dambreach discharge have been developed and implemented by Singh and Snorasson (1984), Wang et al. (2008), and in the BREACH model (Fread, 1988). In these models, mathematical expressions for the dynamics of the outgoing discharge during a dam breach are derived from simplified landslide dam examples. Typically, a feedback loop between outflowing discharge and the amount of material entrained from the landslide dam determines the dynamics of the hydrograph. While these provide a useful estimation of the relevant physical processes during a dam breach, only outgoing discharge is simulated, and downstream processes are unknown. Moreover, the accuracy is generally low for more complicated cases (Zhu, 2006). The BREACH model simulates the increasing breach depth in a landslide dam using an iterative numerical solution. At the sides of the entrained channel, a limiting angle determines the additional collapse of material. Valiani et al. (2002) improved this by simulating dambreach discharge using a twodimensional finite element method.
Fan et al. (2012) provided an insightful step towards integrated modeling by linking a onedimensional breach outflow model with a hydraulic twodimensional flood simulation. The outflow from the BREACH model determines the boundary condition for the flood model (SOBEK). With this combined setup, it was possible to predict the dynamic dam breaching and the resulting largescale flood behavior with significant accuracy. However, this integrated setup was still limited by the assumptions in the model. The breach model is onedimensional and uses a simplified shape for the estimation of breaching dynamics. The setup ignores catchmentscale hydrological processes that could influence the surface flow. Furthermore, breach outflow can typically contain large amounts of solid material, altering the dynamics of the mixture. Fan et al. (2012) implemented a flood model in which flow is calculated using the Saint Venant equations for shallow flow, which ignores forces such as viscosity and implements a fluidbased frictional model (Dhondia and Stelling, 2002). Li et al. (2011a) provided a different approach to integrated simulations of landslide dams by linking the BREACH model with both a regional rainfallrunoff model and the Flo2D debris flow model for modeling the runout of the breach material. Despite their improved method, the landslide dambreach modeling depends on assumptions such as constant flow material properties, a limiting region for entrainment, and landslide initiation coming from predefined boundary conditions. More recently, Mergili et al. (2017) have shown the application of diluting mass flows to modeling glacial lake outburst floods. They developed the r.avaflow model, based on twophase mixture equations by Pudasaini (2012), that allows for mixed flow of water and solids and implements a simplified entrainment process. However, this work does not take into account catchmentscale processes or actual landslide dam formation.
In this research, we aim to simulate a complex multihazard chain using a physically based integrated model. We aim to present an updated version of the OpenLISEM hazard model and analyze the application of this model to complex coseismic hazard chains. The model includes a complete spatial simulation of a landslide dam process chain, including initial slope failure, landslide runout, deposition, lake formation, dam breaching, and flooding. To test the behavior of the developed model, simulations and validation will be shown for a case study of a dambreak flood event in the Hongchun watershed, which occurred in 2010 near the town of Yinxiu, close to the epicenter of the 2008 Wenchuan earthquake, in Sichuan Province, China. Finally, we investigate the predictive capabilities of complex multihazard multistage simulations by analyzing the sensitivity of the model to changes in input parameters.
The investigation of the Hongchun watershed builds on previous works in the literature. Tang et al. (2011, 2015) describe the coseismic and postseismic landslide events between 2008 and 2011 in this area. Several other studies have simulated the event using a variety of modeling techniques. Ouyang et al. (2015) applied shallowflow depthaveraged debris flow equations in order to understand the event as a simplified runout process. Zhang et al. (2018) utilized a depthaveraged smooth particle hydrodynamics model to simulate both runout from landslides and the later debris flow. The authors show a novel application of such methods to a multistage event, however, without an integration in catchmentscale hydrology and a physical implementation of entrainment and breaching of the landslide dam. Domenech et al. (2019) used a multievent debris flow model including entrainment to study the effect of material depletion on debris flow initiation in the Hongchun watershed. Using results of modeling studies, Chen et al. (2016) performed a cost–benefit analysis for the mitigation measures to protect the touristic town of Yinxiu, located directly opposite to the outlet of the Hongchun watershed on the other side of the Min River.
In this research, there is a focus on the integrated nature of the developed multihazard model. Multihazard multistage behavior is simulated for a series of interacting earthquake, landslide, debris flow, and flood processes in the Hongchun watershed. In order to simulate the behavior of this complex event, we develop an extended and improved version of OpenLISEM Hazard. To analyze the uncertainties in modeling such process chains, we employ ensemble simulations and analyze spatial hazard probabilities to estimate reliability. Finally, we discuss the benefits, downsides, and potential application of modeling methods that involve integrated multihazard process chains.
This section presents the theoretical background and governing equations for the components of our multihazard model. One of the important cornerstones is the work by Pudasaini (2012), who developed a set of physically based twophase mass movement equations that can adapt the internal forces in the flow based on the local volumetric solid concentration. This allows us to simulate the behavior of landslides, the flow of water, and the interactions between mass movements and water flow (Mergili et al., 2018a). Using these equations, Bout et al. (2018) developed an integrated model for slope hydrology, slope failure, mass movements, and runout.
We have implemented a series of new processes within the physically based OpenLISEM Hazard model. In this section, we describe the relevant theory of existing functionality, as well as the addition of terrainaltering entrainment by mixture flows and the simulation of coseismic shallow landslides. An overview of the processes is provided in Fig. 1.
2.1 Hydrology
The principle setup of the model starts with the surface aspects of the hydrological cycle. Spatially and temporally distributed rainfall forms the main input of water. We implement a set of equations for interception and microsurface storage based on the original version of OpenLISEM (Jetten and de Roo, 2001; Baartman et al., 2012). The Green and Ampt infiltration model is implemented, which assumes a wetting front moving down into the soil due to infiltrating rainfall (Green and Ampt, 1911). The resulting potential infiltration is subtracted from the available surface water (Eq. 1).
where f_{pot} is the potential infiltration rate (m s^{−1}), F is the cumulative infiltrated water (m), θ_{s} is the porosity (m^{3} m^{−3}), θ_{i} is the initial soil moisture content (m^{3} m^{−3}), ψ is the matric pressure at the wetting front ($h=\mathit{\psi}+Z$) (m), and K_{s} is the saturated conductivity (m s^{−1}).
A separate modeling component is used to simulate several months of groundwater dynamics and generate the initial ground water conditions for the OpenLISEM Hazard simulation, which simulates the effects of a single rainstorm. Groundwater flow is estimated using a depthaveraged Darcy equation (Eq. 2), as is common in regional hydrological models (Van Beek, 2002).
where h is the hydraulic head (m), and V_{d} is the Darcy flow velocity (m s^{−1}).
2.2 Slope stability and failure
Hydrology influences slope stability and can eventually lead to failure. Slope failure is based on the iterative failure method (Bout et al., 2018). This technique reverses the factor of safety (Eq. 3) to solve for the remaining depth of material at which the local situation becomes stable. The locally altered terrain then results in changed forces in the surrounding cells. Through iteration, the method keeps removing material until no unstable cell is left and the minimum required material for a stable terrain has been removed. We added the seismic forcing in the factor of safety calculation following Morgenstern and Sangrey (1978).
where FOS is the factor of safety (−), β is the slope of the soil section (−), c^{′} is the apparent cohesion of the soil (kPa), α is the peak horizontal earthquake acceleration (m s^{−2}), ϕ is the internal friction angle of the soil (−), γ is the density of the slope material (kg m^{−3}), γ_{w} is the density of water (kg m^{−3}), m is the fraction of the soil depth that is saturated from the basal boundary (−), and h_{s} is the depth of the failure plane (m).
The apparent cohesion consists of additional root cohesion and a matric suction term. The acceleration is assumed to be, as in the most critical situation, the estimated peak horizontal ground acceleration. To account for subsurface force propagation, we include an iterative force solution. This is an extension of similar approaches in Cheng and Zhou (2013) and Chen et al., (2010). In the proposed implementation, forces are iteratively solved throughout the entire terrain description (Fig. 2).
In the threedimensional case, where the x and y components of the seismic forcing and slope steepness influence the propagation, this can be expressed as Eq. (4).
with C the force capacity (numerator of Eq. 3), D the force demand (denominator of Eq. 3), S the normalized slope vector (−), and F_{lat} the vector of laterally acting forces (kg m s^{−2}).
We assume that excess force is transferred downslope, but the fractured top material is unable to transfer force resistance upslope.
We assume in the setup of the model that a failure plane can develop in the soil materials at any depth. Equation (3) can be inverted to find the value of h, for which the safety factor is 1. To solve this equation for h, we first express the slope based on the local elevation differences and create a shortened equation for the factor of safety (Eqs. 5 and 6).
where the simplified constants are given by equations C_{1}=c, ${C}_{\mathrm{2}}=\left(\left(\mathit{\gamma}m{\mathit{\gamma}}_{\mathrm{w}}\right){h}_{\mathrm{s}}+m{\mathit{\gamma}}_{\mathrm{w}}{h}_{\mathrm{s}}\right)$, and ${C}_{\mathrm{3}}=\left(\left(\mathit{\gamma}m{\mathit{\gamma}}_{\mathrm{w}}\right){h}_{\mathrm{s}}\right)$, β is the slope angle (−), z is the elevation above the failure surface (m), and z_{0} is the lowest neighboring elevation (m).
Solving this equation can be done using the trigonometric identities (Eqs. 7 and 8).
Finally, we find the lowest real root to the secondorder polynomial equation of the form
where
Using Eq. (9) the critical depth can be found where the material on the slope is in equilibrium. By multiplying the area of a pixel with the slope material height above this critical depth, a failure volume can be calculated. This volume consists of solids and water, depending on the soil saturation level, and is then added to the flow equations that simulate Mohr–Coulomb mixture flow.
2.3 Flow dynamics
Once a volume of flowing material is introduced (either through hydrology or slope failure), its movement must be described. The dynamics of mass and momentum are calculated using a discretization of the continuity equation for mass and momentum of both solids and fluids (Eqs. 13, 14, and 15). Equation (13) described the conservation of mass, and Eq. (14) describes the conservation of momentum. On the lefthand side the storage term and advective terms are found. External forces and mass sources are added on the righthand side.
with h being the flow height (m), u the flow velocity (m s^{−1}), R the rainfall (m), I the infiltration (m), g the gravitational acceleration (m s^{−2}), S the friction term (m s^{−2}), and S_{f} the momentum source term (m s^{−2}).
The momentum source terms predominantly determine the behavior of the flow. Here we implement pressure and friction forces for fluids as is typical for shallow water assumptions (O'Brien et al., 1993).
The friction slope terms, which are the friction forces divided by the water height and the gravitational acceleration, can be calculated using the Darcy–Weisbach friction law (Chow, 1959) (Eq. 16).
with n being the Manning's n friction coefficient (s m^{−⅓}).
The momentum source terms within OpenLISEM Hazard are based on the work by Pudasaini (2012). This set of equations contains a physically based twophase momentum balance. Besides pressure and gravitational forces, it includes viscous forces, nonNewtonian viscosity, twophase drag, and a Mohr–Coulomb type friction force for the solid phase (Eqs. 19, 20, 21, 22). Based on the current and local state of flow, forces change in magnitude. This approach allows for a smooth transition between nonviscous flow, hyperconcentrated streamflow, debris flows, or landslide runout by automatically scaling and solving the interactions between solids and fluids.
S_{s} is the momentum source term for solids (m s^{−2}), S_{f} is the momentum source term for fluids (m s^{−2}), α_{s} and α_{f} are the volume fraction for solid and fluid phases (−), P_{b} is the pressure at the base surface (Kg m^{−1}s^{−2}), b is the basal surface of the flow (m), N_{Re} is the Reynolds number (−), ${N}_{R{e}_{A}}$ is the quasiReynolds number (−), C_{DG} is the drag coefficient (−), ρ_{f} is the density of the fluid (kg m^{−3}), ρ_{s} is the density of the solids (kg m^{−3}), γ is the density ratio between the fluid and solid phase (−), χ is the vertical shearing of fluid velocity (m s^{−1}), ε is the aspect ratio of the model (−), and ξ is the vertical distribution of α_{s} (m^{−1}).
Our model follows the equations by Pudasaini (2012) for the definitions of the regular and interface Reynolds number, which represents the ratio between inertial and viscous forces within the fluid (Eqs. 23 and 24).
where L is the length scale of the flow (m), H is the height of the flow (m), η is the viscosity (kg s^{−1} m^{−1}), and A is the mobility of the interface, which is a scaling parameter for nonNewtonian viscousfluid stresses ( $\approx \mathrm{1}{\mathit{\alpha}}_{\mathrm{s}}$) (−).
To apply these twophase equations successfully in a catchmentbased model, we replace the frictional force for the water phase with the Darcy–Weisbach equation for water flow friction. To complete the set of equations that govern debris flow dynamics, several flow properties were estimated based on the volumetric sediment content. Viscosity is based on an empirical relation by O'Brien and Julien (1985) (Eq. 25).
where α_{s} is the volumetric solid content of the flow (−), α is the first viscosity parameter (−), and β is the second viscosity parameter (−).
The drag coefficient is based on the relation provided by Pudasaini (2012) (Eqs. 26, 27, 28, and 29).
where F and G are the scalar functions describing the flow velocity of solids and fluids respectively (−), Re is the particle Reynolds number (−), d is the median grain diameter (−), U_{T} is the settling velocity (m s^{−1}), and M is an empirical parameter depending on the Reynolds number (≈0.2) (−).
Finally, the settling velocity of small (d<100 µm) grains is estimated by Stokes equations for a homogeneous sphere in water (Stokes, 1850). For larger grains (>1 mm), the equation by Zanke (1977) is used (Eq. 30).
in which U_{T} is the settling (or terminal) velocity of a solid grain (m s^{−1}), η is the dynamic viscosity of the fluid (Pa s), ρ_{f} is the density of the fluid (kg m^{−3}), ρ_{s} is the density of the solids (kg m^{−3}), and d is the grain diameter (m).
2.4 Deposition and dam formation
Flows such as the ones described by the equations presented above have complex interactions with the basal surface. In particular, material in the flow can be deposited, or material from the terrain can be entrained. Deposition occurs when the flow velocities of a solid–fluid mixture have sufficiently low velocities and when drag forces are insignificant compared to shear resistance. Then, water and solids are subtracted from flow volumes to form a saturated deposits layer. Further deposition occurs as an active process, during flow, based on the deposition equations of Takahashi (1992). These equations use local stability analysis and the ratio between the flow velocity and the critical velocity to estimate deposition of solids (Eqs. 31, 32, 33, and 34). Using these generalized deposition equations, a variety of more specific depositionbased processes can be potentially simulated. For example, deposition of landslide material in a river would equate to landslide dam formation in rivers.
where D is the deposition rate (m s^{−1}), u_{cr} is the critical velocity for deposition (m s^{−1}), p is the calibration factor for the critical velocity for deposition (−), α_{eq} is the equilibrium volumetric solid concentration (−), α^{b} is the volumetric solid concentration of the bed material (−), and d_{50} is the median grain size (m).
2.5 Entrainment equations
In order to estimate entrainment, we implement the equations by Takahashi et al. (1992) in a similar manner as was done in the Edda model (Chen and Zhang, 2015). The expressions for the entrainment rate are provided by Eqs. (35), (36), and (37).
where E is the rate of change in the basal topography (erosion rate) (m s^{−1}), τ is the shear stress (Pa), τ_{c} is the critical shear stress (Pa), S_{f} is the surface friction term (−), τ_{y} is the yield stress (Pa), K is the resistance parameter for laminar flow (−), n_{td} is the turbulent dispersive coefficient (m^{½} s^{−1}), c^{b} is the cohesion of the bed material (Pa), and C_{s} is the coefficient of suspension (−).
The surface shear term is calculated from the momentum conservation equation as the sum of all surface frictional terms. We adapt the expressions to conform to the momentum jump boundary condition, as provided by Iverson and Ouyang (2015) (Eq. 38).
where ρ_{eff} is the total effective density of the flow (kg m^{3}), and u_{b} is the basal velocity (m s^{−1}).
Iverson and Ouyang (2015) indicate that the equations by Takahashi et al. (1992) do not include a distinction between basal and mean velocities and do not conserve momentum. The basal velocity differs from the mean velocity of the flow according to a vertical velocity profile. Our model description does not include an estimate of the basal flow velocity. Therefore, we use the entrainment coefficient as a scaling device to convert from the mean velocity to a basal velocity (${u}_{\mathrm{b}}\approx K\overline{u}$), with $\overline{u}$ the depthaveraged velocity of the flow. This approach was similarly taken by Pudasaini and Fischer (2016). Therefore,
Furthermore, we ignore vertical velocities in the model setup, a both common and necessary assumption that has given good results in other models (Iverson and Ouyang, 2015). Finally, we alter the mass and momentum source terms to include the produced mass and momentum from entrainment and thereby hold to conservation of these quantities. We introduce into the flow, based on the entrainment rates, a new mixture of solids and fluids. The momentum of the flowing material is additionally increased by adding the momentum of the entrained mass, which gains a velocity equal to the estimated basal velocity. Lateral entrainment includes collapse due to slope failure which is included by solving for stable depth. This is automatically solved by the iterative failure method as described previously. Slope failure is based on a local safety factor estimation which is inverted as is done for slope stability.
2.6 Flowchart and numerical implementation
A flowchart of the full model is given in Fig. 3. Hydrology forms the basis of the simulated cycle. From this, other sediment and solidrelated processes are linked. At each time step, properties of the flow are determined based on the actual water and solid content. Furthermore, solid properties such as density, friction angle, and size are advected with the solids. The flowcharts indicate the data required for each step of the simulation. This highlights the downside of integrated modeling approaches: the increase in required input data.
Similar catchmentintegrated simulations that include hydrology, entrainment, and debris flow initiation have been used successfully by the EDDA model (Chen and Zhang, 2015; Hu et al. 2016) and by the STREP TRAMM model (Lehmann et al., 2019). These models predict entrainment of slope material or landslide deposits but do not alter the elevation model. By using recent developments in numerical schemes, elevation model changes can be simulated without erroneous feedbacks between terrain and flow. The numerical scheme within the model is based on the monotonic upstream cellcentered scheme (MUSCL) (Van Leer, 1979). This scheme uses piecewise linear reconstructions of both the terrain surface and the flow properties. This terrain reconstruction estimates local slope and, in rough terrain, also includes cell boundary elevation differences. Besides shear stress from the bottom surface in a cell, cell boundary barriers can also provide shear resistance to the flow. Our implementation of entrainment therefore takes into account direct and lateral entrainment.
The integrated OpenLISEM multihazard model was applied in the Hongchun watershed, located near the town of Yingxiu in the epicentral area of the 2008 Wenchuan earthquake, in Sichuan province, China. This watershed experienced numerous coseismic landslides and a catastrophic debris flow event in 2010. The debris flow deposits from this watershed dammed the main Min River and flooded the town of Yinxiu (Fig. 4).
This watershed has an area of 5.3 km^{2}, and elevation ranges between 900 and 1700 m. The very steep slopes of this watershed (>30^{∘} slope for more than 86 % of the catchment) were covered by dense forests before the earthquake in 2008. The steep terrain and sharply incised channels provide few locations for human settlements. The touristic town of Yinxiu is located adjacent to the Min River and opposite of the outlet of the Hongchun watershed. The lithology of the area consists mainly of highly fractured granitic rocks, with some pyroclastic rocks, limestones, and sandstones (Tang et al., 2011). The texture of the weathered material is predominantly clay loam with large amounts of gravel. The Beichuan thrust fault runs straight through the watershed (Fig. 4; Mahodja et al., 2016). In 2008, the fault was ruptured in the Wenchuan earthquake (M_{w}7.9). After the earthquake, numerous landslides occurred in the area, leaving large volumes of deposits in the streams and channels and removing the vegetation from 50 % of the area. A detailed coseismic landslide inventory from Tang et al. (2016) is shown in Fig. 4. The northern part of the catchment, which is part of the hanging wall, was more impacted than the southern part. One of the largest landslides in the watershed blocked the main channel of the Hongchun catchment (Fig. 4), with the deposits exceeding 8 m in depth and 70 m in width.
3.1 The multihazard event of 2010
In 2010, 2 years after the earthquake the area experienced several rainy weeks followed by a highintensity rainfall event. This event consisted of two peaks of several hours of rainfall, with intensities up to 33 and 31 mm h^{−1} and a total event duration of 15 h. In total, 220 mm of rainfall fell in the two days surrounding the event.
During this event, a debris flow was generated by entrainment of loose landslide deposits, and the landslide dam located in the central channel was breached. Due to the dam breaching, the volume of the flow increased substantially, also due to more entrainment downstream. Upon leaving the Hongchun watershed, the debris flow material deposited in the Min River, with a thickness of up to 15 m high and an estimated total volume of 7.11×10^{5} m^{3} (Tang et al., 2011). The Min River, which experienced a high discharge at this moment, was diverted laterally and flooded parts of the nearby newly reconstructed the town of Yinxiu (Fig. 6). For a more detailed description of the event, see also Tang et al. (2011) and Ouyang et al. (2015). A schematic overview of the event is provided in Fig. 5.
In summary, the study area experienced two different multihazard chains that will be simulated.

Coseismic chain. The first one was experienced during the earthquake, when ground shaking and topographic amplification triggered a series of coseismic landslides, some of which blocked the Hongchun stream.

Postearthquake chain. Heavy rainfall triggered debris flows due to entrainment, which broke the coseismic landslide dam. The resulting debris flow dammed the Min River, which was diverted into the town of Yingxiu.
All simulations include entrainment of estimated soil depth, hydrology, and the related surface processes as described in the theory section.
3.2 Model input and parameters
The input data are based on a combination of laboratory or field measurements and drone or satellitebased spatial data products (Table 1).
A preearthquake digital terrain model with 20 m spatial resolution was obtained from the local government. Unfortunately, we were not able to obtain betterquality preearthquake elevation data. Because of this a comparison of the preearthquake contour lines with postearthquake elevation does not yield useful data. Postearthquake surface elevation data were available at 2 m spatial resolution, acquired from fixedwing drone flights (see Fig. 6) and filtered using the Pix4D digital terrain model (DTM) filter to remove vegetation (Pix4D, 2017). Flow simulations typically are most dependent on the accuracy and spatial resolution of the elevation model. In order to have an effective compromise between detail and computation time, we resampled all input base maps to the final resolution of 10 m, at which the simulations were carried out.
Preearthquake normalized difference vegetation index (NDVI) values were calculated at 30 m resolution Landsat8 images from 2008. A postearthquake NDVI map was obtained from 4 m spatial resolution SPOT images acquired in 2009. NDVI values were used to estimate leaf area index (LAI) by applying an empirical relation obtained from tropical forest data (Kalacska et al., 2004). We used it to estimate fractional vegetation cover using a similar empirical relationships (Jiang et al., 2006).
Spatial seismic acceleration data were obtained from the USGS ShakeMap product (Worden et al., 2010), developed using combinations of measurements and intensity prediction equations, and peak ground acceleration (PGA) values were around 1.5 g for this earthquake in the Hongchun area. The PGA value reflects estimated surface shaking but does not take into account the complex variability in site effects such as amplification due to soil or topography. Landslide material strength parameters were obtained from triaxial strength tests performed for engineering reports ordered by the local government (Table 2; Yang, 2010; Hao et al., 2011; Li et al., 2011b). The consolidated drained triaxial tests were performed on samples acquired from both the deposition and source areas, although exact locations are unknown. The resulting average values for cohesion of 7.3 kPa and internal friction angle of 27^{∘} were calculated from 16 samples. While the cohesion and density match the values found in other studies (Ouyang et al., 2015; Domenech et al., 2019), the internal friction angle was estimated to be between 30 and 40^{∘} in a previous study (Ouyang et al., 2015), based on measurements of internal friction angle on similar igneous rocks in the nearby Xiaojigou ravine. The difference between these values can be explained by the material used in the triaxial tests. In our case, material from landslide depositional areas or areas close to a landslide source was used. These samples contain more fractured and loosened material. Thus, these measurements will provide a more appropriate value in postearthquake entrainment simulations. The root cohesion is not used for slope stability in the model as roots do not reach the depth of the typical failure plane. Instead, the root cohesion is used to increase shear strength for the entrainment prediction. Here, vegetation cover influences the predicted erosion. However, for preearthquake slope stability estimations, a value of 35 ^{∘} was used, as was done in other studies (Ouyang et al., 2015; Domenech et al., 2019). Textures were found to be clay loam with large amounts of gravel. Saturated infiltration rates were measured in the field by using field ring infiltration tests. Because of the large amounts of gravel and macropores within the materials structure, infiltration values were relatively high. Values were obtained for a total of four locations, which gave an average value of 65 mm h^{−1} for saturated conductivity. Other soilrelated parameters such as porosity, density, and matric suction were obtained using the pedotransfer functions from Saxton et al. (2006).
Soil depth values were obtained by applying the spatial soil depth model from Ruette et al. (2013), which uses steadystate assumptions to balance material production with soil movement according to an empirical formulation of erosion and lateral transport. The production of weathered material depends on the weathering rate of bedrock. Since weathering rates are generally difficult to obtain, this parameter becomes the main calibration parameter. Validation was done using soil depth values obtained from landslide scarps. In this study, we used the pre and postearthquake elevation model differences provided by Tang et al. (2019). From their data, we sampled elevation model differences in landslide source areas to obtain failure depths, assuming that the failures were at least as deep as the top layer of weathered material. In the case of shallow landslides, we took the maximum landslide depth as soil depth. For the larger, deepseated landslides, no samples were taken since it was not possible to obtain a good estimate of the depth of the top layer. Additionally, we simulated a second layer within the slope stability calculations with an additional depth of 20 m to capture any deeper rotational slides within the bedrock material. We assumed that this second layer did not contain groundwater, had an internal friction angle of 35^{∘}, was subject to subsurface lateral forcing, and was influenced by the weight of the upper layer. Calibration results and cumulative distributions of depth values are shown in Fig. 7.
The effect of vegetation on slope stability and entrainment rates were taken into account by adding root cohesion. The recovery of vegetation in the landslideaffected areas developed through distinct stages (Huang et al., 2017). We estimated the average root cohesion for the various vegetation classes using a combination of literature data and measurements (Schmidt et al., 2001; De Baets et al., 2008; Chok et al., 2015). We used the method as described by De Baets et al. (2008) which requires measuring root critical force and root diameters for all significant roots within a specific surface area. Then, after converting force to tensile strength, the total effect of the roots per unit area can be estimated using Eq. (40).
where c_{root} is the added apparent root cohesive strength (kPa), i is the root diameter class, T_{i} is the root tensile strength (kPa), n_{i} number of roots within the diameter class, a_{i} is the crosssectional area of the root, A is the area of the soil occupied by roots (m^{2}), θ is the angle of shear distortion in the shear zone (^{∘}), and ϕ is the internal friction angle (^{∘}). Note that the term related to the shearing angle (sin (θ)+cos (θ)tan (ϕ)) is usually assumed to be approximately 1.2 (Baets et al., 2008). Here, we used a spatial calculation of this variable.
A total of 12 measurements of root cohesion averaged over a 0.1 m^{2} area were done. On average, 108 roots were found per test site of 0.1 m^{2}, with an average diameter of 3.6 mm. Average values for young, postlandslide vegetation were found to be 4.8 kPa, and averages for mediumsized vegetation were found to be 6.2 kPa. For mixed forest, literature values were used, and root cohesion was estimated to have a value of 8 kPa (Chok et al., 2015; Huang et al., 2017). For each land cover class, the root cohesion from the closest measured vegetation type was used. Finally, within each land cover class, we linearly scaled the root cohesion to fractional vegetation cover values that were derived from the NDVI both before and after the earthquake.
3.3 Calibration and validation
The predicted coseismic failure areas were calibrated using the landslide inventory and mapped deposition based on aerial imagery and pre and postevent elevation data differences. Within the inventory, no separation between source and deposit areas was provided. Therefore we assumed that, based on field visits, on average the 25 % highest part of each landslide polygon represented the source area. In order to compare the outputs of our model with others, we parameterized two other regional slope stability models. The first of these is Scoops3D, which uses random spheroid sampling to find the landslides with the lowest failure volumes and lowest factors of safety for each pixel (Reid et al., 2015). This model can also use seismic shake maps as input. The second model is r.slope.stability, which uses random ellipsoid sampling to find the lowest factor of safety and failure depth for each pixel (Mergili et al., 2014). The version of the model that was used did not incorporate seismic acceleration, so the output is calibrated without addition of seismic forcing. Recent developments of this software have added these effects. The Newmark displacement method (Newmark, 1965) was not used in our comparison. The primary reasons were the lack of predicted failure depths and the high similarity to the infinite slope model. When the critical acceleration is taken as threshold and failures occur when accelerations go beyond this value, results are identical to an infinite slope model that incorporates seismic acceleration. For both Scoops3D and r.slope.stability, the subsurface description is identical to the input for OpenLISEM Hazard. As a measure of model fit, we used Cohen's kappa value. This metric has benefits over simple accuracy, especially for modeling landslide occurrence, since it compensates for the large amounts of true negative predictions (no landslide in model and inventory) that usually dominate landslide study sites (Bout et al., 2018).
The second phase of the modeling first simulates hydrology and a debris flow initiated by entrainment. Then, within the same simulation, the Min River is blocked by debris flow material, and the town of Yinxiu is flooded. Calibration for this part of the modeling is based on mapped deposition extent in the Min River, the flood extent in the town of Yinxiu, and the estimated depositional volume (approximately 7.11×10^{5} m^{3}).
Calibration of the entrainment and debris flow runout was based on the final spatial deposit extent at the catchment outlet. The calibration parameters for each part of the simulation are shown in Table 3. An exception to this was the entrainment constant, for which no clear guideline was known to determine the value for specific terrain types. Based on earlier simulations and flume tests from Takahashi et al. (1992), a starting value was chosen of 0.05. For initial soil moisture content, the value is cut off at full saturation. Initially, each parameter was varied by choosing values between 60 % and 140 % of the original value. After initial calibration was done, the parameters were adjusted according to the steepest descent principle in order to find the best set of parameter values. To have a level of validation in the simulations, the parameters resulting from calibration of the first chain were used as input for the second chain.
Simulation results for the first chain, coseismic slope failure and runout, are shown in Figs. 8 and 9. Both the accuracies and Cohen's kappa values for each of the slope stability models used are shown in Table 8.
The general spatial patterns are predicted reasonably by all models, but the models differ considerably in details. For the ellipsoid and spheroid sampling done by r.slope.stability and Scoops3D, the failures are larger than those actually mapped (Table 4). For the iterative method, sizes are mixed but much more similar. Major landslides, in particular in the north, are predicted with similar size, although not in the exact location, by all models. The highest accuracy (81 %) is obtained using r.slope.stability and the highest Cohen's kappa (0.232) value with OpenLISEM Hazard. The OpenLISEM Hazard iterative method shows the best reproduction of the general pattern, in particular since both Scoops3D and r.slope.stability lack slope failures in the southern half of the catchment. Of the three models, r.slope.stabiltiy and OpenLISEM Hazard show similarity in total failure area when compared to the inventory. Despite the achievements of the model, errors remain large, and Cohen's kappa values indicate all models have low to moderate reliability.
4.1 Runout and the blocking of the Hongchun stream
When slope failures are simulated, OpenLISEM Hazard automatically introduces landslide runout by transferring the failed volume and its properties to the Mohr–Coulomb solid–fluid mixture flow equations. The depth of slope failures determine directly the amount of solids and fluids introduced (Fig. 9D). The landslide material moved down the slopes into the main channel of the Hongchun watershed, blocking it in at least one location (Fig. 9A). The accuracy for the calibrated simulation was 64 % with a Cohen's kappa value of 0.28 (Table 5), which is mainly attributed to the accuracy of the failures that started this process (Fig. 9B). Runout distances are similar to those mapped, indicated by the landslides reaching the main channels within the Hongchun catchment but not reaching the Min River. In one location the channel was blocked by a deposit 12 to 18 m deep which was simulated with high accuracy as compared to the mapped blockage (Fig. 9B). Engineering reports indicate a similar depth (16 m) for the landslide dam (Yang et al., 2010; Hao et al., 2011; Li et al., 2011b). We considered using the landslide initiation polygons from the inventory to increase the accuracy of the runout simulation. However, since the aim of this research is to provide a true multistage modeling setup, we used the integrated prediction of slope failures as input in the runout modeling.
4.2 Validation of failure and runout for the major central landslide using elevation model differences
Two large landslides occurred in the northern part of Hongchun catchment. One of those, near the center of the area, blocked the main channel of the catchment. For this landslide, pre and postearthquake elevation models from lidar data are available (Tang et al., 2019), which were used to calculate the landslide volume. For this particular area, we compared the predicted elevation model differences based on slope failure and deposition from our modeling chain with the results from the lidar DEMs (Fig. 10).
Both the slope failure and runout thicknesses are predicted with high accuracy for this landslide. The landslide deposits of around 20 m thickness remained in the main channel without spreading significantly. These deposits have later been mined as materials for local construction. As can be seen in Fig. 10 the failure, deposition, and elevation differences are highly similar with matching spatial patterns.
4.3 Simulation of the second multihazard chain
The results from the modeling of the second stage show a full physically based simulation that reproduces the behavior and the impact of the event (Fig. 11). This rainfall event with two distinct peaks resulted in a rainfall amount of 220 mm in 2 d, which was modeled in time steps of 0.5 s. Due to the large rainfall volume, runoff increased rapidly leading to large amounts of sediments that were entrained from the coseismic landslide deposits. Within the simulation, entrainment takes place after runoff has been converted into streams. There, higher water pressures and velocities provide shear stress on the surface that is sufficient to overcome the materials internal stability. Near the outlet of the watershed, entrainment decreased due to a decrease in slope steepness. This decreased flow velocities but also increased the internal stability of the available material. Because of this, entrainment stopped, the flowing material lost momentum, and it was finally deposited in the Min River.
The simulated debris flow reaches the Min River, and the deposits accumulate in the river. Most of the momentum has been lost before the material leaves the Hongchun watershed, and even though the Min River has a velocity of 7 m s^{−1}, the deposited volume is too large to be eroded by the Min River. Both the total deposit volume in the river and the extent of the deposits compare well to the measured and mapped values. The deposits were modeled with an accuracy of 92 % and a Cohen's kappa value of 0.84 (Table 6). The modeled deposited volume is 5.82×10^{5} m^{3} (the observed estimation of the deposited volume is 7.11×10^{5} m^{3}).
As a final stage of the event, the model predicted the flood behavior in the town of Yinxiu by the Min River as its main course was blocked by the debris flow deposits. The modeled damming area correlates well with the observed deposit volumes and flood extents, obtained from the interpretation of aerial images and field photos from the event (Fig. 14). The flood accuracy is 91 % with a Cohen's kappa value of 0.81 (Table 6). Since the elevation data were upscaled from 2 m resolution to 10 m, the multistory buildings merged with small streets to become a joined obstacle for the flow, which corresponds with the visual observations.
While validation using comparison to mapped extent is a useful tool, it does not guarantee all aspects of the hazard are correctly simulated (e.g., timing, velocities, heights). While a model might have correctly predicted the extent of impact, it might predict the occurrence too early or late, with incorrect velocities or flow heights. While we do not have data on flow heights or velocities, we confirmed with local authorities the debris flow was reported at 03:30 LT (Sichuan local time, GMT + 8 h). In our simulation, the first peak of solid material at the outlet occurs at this time (Fig. 12). The major peak of debris material occurs later, at 06:00 LT, but we could not confirm in the town of Yinxiu for how long the debris kept arriving at the Min River.
5.1 Uncertainties in modeling the multihazard chains
Several major obstacles in the simulation of coseismic landslide occurrence can be identified, which are related to the assumptions and techniques used in the models. The iterative method that is implemented in OpenLISEM Hazard assumes that, at least initially, failure surfaces are parallel to the terrain surface. This assumption can lead to a variety of issues when the terrain has smallscale variations that do not represent the overall topography. This assumption is not present in the random ellipsoid sampling methods. The iterative method shows similar failure depths and locations but generally separates slope failures more, whereas the other models show larger joined failures. The predominant cause of this behavior most likely lies in the subsurface force iteration, where the potential failure surface is also initially assumed to be parallel to the surface. Then, when an obstacle exists that alters the slope around several pixels, the subsurface forces are blocked, and potentially larger failures become disconnected (Fig. 13).
Secondly, structural geological input data are not available and cannot be used in the iterative method. For threedimensional analysis using random ellipsoid sampling, structural weaknesses could be implemented by altering strength properties of specific layers. When such data are available, applications of the detailed random ellipsoid sampling method can allow for a greater predictive value (Mergili et al., 2014; Cance et al., 2017; Tun et al., 2019).
Furthermore, seismic acceleration of material is a complex and dynamic process, driven by seismic waves compressing and stretching the sloping materials. Such waves reflect and refract based on material properties and boundary conditions, which leads to a variety of topography or materialbased amplifications. These amplifications are generally important in landslide hazard modeling (Jafarzadeh et al., 2015). However, shake maps produced by the USGS utilize empirical predictive functions to spatially extend ground accelerations and ignore these crucial local amplification effects. As a result there is high spatial uncertainty in the peak ground acceleration values we have used as input. Such uncertainties can be overcome by linking slope stability approaches to seismic wave modeling, a computationally heavy task that has of yet not been performed for the Wenchuan earthquake.
The final calibrated values of the input parameters were within a reasonable range (between 50 % and 150 %) with respect to the original estimated ones. Two parameters stood out significantly in the calibration process. The first of these, soil depth, is a crucial parameter as it determines the amount of material that is released, which is of direct influence in the runout modeling. This can indirectly influence the amount of material available for entrainment by flow at a later stage. In an active landscape, soil depth patterns are determined by rock weathering rates and mass transport/wasting. Soil depth increases until there is sufficient material to induce slope failures or significant erosion. As a result, at many locations, combinations of soil depth and slope can be near a critical state. The spatial component of soil depth is therefore of high importance. Despite its importance, soil depth is very difficult to measure over larger areas and has a high uncertainty in most model applications (Kuriakose et al., 2009). Because information on soil depth is difficult to obtain through direct observation or remote sensing, it is generally obtained through modeling approaches (Kuriakose et al., 2009; Ruette et al., 2013). Another reason for the uncertainty in soil depth information is that a soil layer is a theoretical concept that does not always translate well into reality, whereas weathering is a more gradual process. A second variable which was difficult to estimate is the entrainment coefficient. Currently, there are several types of entrainment equations available in the literature, but most of these lack significant guidelines for selecting practical entrainment parameters (Iverson and Ouyang, 2015). Besides this coefficient, other parameters used in the entrainment equations, such as soil cohesion, internal friction angle, and moisture content, were measured in the laboratory and are common values for geotechnical research. However, the number of samples that are tested for these geotechnical parameters are always limited, and their spatial variability is generally high.
5.2 Ensemble simulations
To analyze the uncertainties within the multihazard multistage modeling setup, we extended the calibration process to an ensemble analysis. In total, for each of the six calibration parameters, three equalinterval values were used (calibrated values), as well as the values plus or minus a given range of 10 %–50 %, as indicated in Table 7. Thus, the first simulation was repeated 27 times. The best estimate was used as input for the second simulation, which was similarly repeated 27 times. In total, 56 simulations were thus performed, out of which normalized frequencies were generated for hazard occurrence. We defined a threshold of 0.25 m above which a flow is counted as an actual hazard occurrence to avoid that the results would include runoff with insignificant depths.
Figure 14 shows an ensemble plot for each stage of the simulation. For each location, these maps show the normalized probability of hazard occurrence within the ensemble of simulations.
The four maps in Fig. 14 illustrate how much the variation of the input parameters influences the simulated hazard for the four stages of the multihazard chain. The prediction of coseismic slope failures shows the highest spread and therefore the highest influence on the total event variability. The uncertainty in the input parameters influences the slope failure equations most since we have simulated the multistage event as an integrated sequence in which the uncertainties in one process influence various other processes. Despite the uncertainties, there is a substantial certainty that flooding of the town of Yinxiu will occur, independent of the input parameters. For all performed simulations, the deposition volume in the Min River is never below approximately a third (2.1×10^{5} m^{3}) of the estimated volume (7.11×10^{5} m^{3}). For at least 80 % of the simulations, there is at least some flooding experienced with a depth above 50 cm in Yinxiu. In comparison to other multistage hazard studies, these numbers are relatively high (Mergili et al., 2018b). Mergili et al. (2018b) reported that threshold behavior and nonlinear effects dominated the alterations in model behavior with changing parameters. While the event they describe is quite different from the Hongchun disaster, there is a similar influence from threshold effects on the interpretation and certainty of model outcomes. In this current study, threshold effects most likely exist related to the landslide dam behavior in both the Hongchun watershed and the Min River. With altered input, there could be insufficient material for landslide dam creation or not enough flow and entrainment for breaching. However, we do not encounter scenarios that fall below the threshold of breaching in our ensemble.
There are several factors that drive the high likelihood of flooding in the town of Yinxiu. The event is initiated by several triggering conditions, or forcings, such as the seismic acceleration and rainfall, which occur in a very steep watershed and determine the behavior of the hazardous processes. Seismic accelerations during the 2008 earthquake reached the 10th and highest USGS category for seismic acceleration (>139 % g). Similarly, the 2010 rainfall event resulted in 220 mm of accumulated rainfall within 48 h, with significant rainfall in the weeks before. For events with a smaller return period and lower intensity of the triggering process, variability and uncertainty in the potential outcomes might be larger. An example is the occurrence of landslides' long duration through lowintensity rainfall, which is a much larger challenge for physically based modeling (van Beek, 2002). In our case, material strength parameters must deviate significantly from their measured and estimated values in order not to result in slope failure under such extreme seismic acceleration. Thus, the severity of the triggering events causes the simulations to predominantly show major slope failures, debris flows, deposition in the Min River, and finally flooding of the town of Yinxiu.
5.3 Usability
A major challenge for the application of such a type of multihazard modeling is the step from backanalyzing past events to forward prediction of future events. In the case of backanalysis, calibration is possible on data from the event itself, whereas predicting future events must be done without knowledge of the hazard that will occur. By calibrating on data of multiple events and validating on other events (checking for accuracy using calibrated input from other events), a more robust model can be established. When validated on similar events, a higher level of confidence can be generated in the abilities of the model to predict future events. In the case of the Hongchun watershed, the accuracy of the presented numerical experiments depended on thorough calibration. This is caused by several factors, such as the uncertainties in the input data and the limited accuracy of the assumptions made within the model equations. In studies where multistage multihazard events are considered, this becomes even more difficult due to threshold effects and stacking uncertainties (Mergili et al., 2018a). Despite these issues, calibration is possible by performing model runs with a large set of varying input parameters and comparing the output with measured data on the event (e.g., discharge, extent of landslides, debris flows, or floods). Validation options, on the other hand, may be very limited as similar events are very rare. Debris flows occur more frequently within the watershed and can therefore be both calibrated and validated (Domenech et al., 2019). However, the multihazard chain as described in this research is a rather unique sequence of events that does not easily allow for validation except when looking at individual components. Even then the coseismic slope failure and runout stages do not occur frequently within the same region. Thus, we run into two primary issues in the predictive use of multistage and multihazard modeling. Firstly, the complex modeling has internal uncertainties that propagate, stack and provide a challenge in calibration. Secondly, such events occur with very low frequency (Kappes et al., 2012) and can therefore rarely be further calibrated and validated. Despite the difficulties in the calibration and validation process for multihazard multistage modeling, we have shown that for the Hongchun watershed, many aspects of the first sequence had a high likelihood despite the uncertainties. This indicates the possible use in predicting future events, especially when ensemble analysis is performed to show the relative spread of model results. However, this is computationally intensive and cannot be carried out over large areas yet.
We have investigated the potential for understanding and predicting complex multihazard coseismic process chains. A physically based model is used that allows for the simulation of multistage and multihazard chain events in mountainous terrain. Catchmentbased hydrology was combined with slope stability analysis under seismic acceleration. Entrainment of loose material was implemented to alter both topography and the composition of the flow. The modeling setup is able to fully simulate the behavior and impact of multistage and multihazard events. The developed model code is available as part of the ongoing development of the opensource OpenLISEM Hazard model (https://lisemmodel.com/, last access: 23 August 2022).
The model was tested in the Hongchun watershed, located near the epicenter of the 2008 Wenchuan earthquake, in two sequences of hazard interactions. The first sequences was the generation of earthquakeinduced landslides that resulted in runout of landslide materials and the blockage of a stream channel. The second sequence was the triggering of debris flows by extreme rainfall, resulting in the breaking of the existing landslide dam, and the formation of a debris flow fan in the main river, forcing the river to flood the city on its banks. The modeled slope failures show a reasonable accuracy when compared with the actual landslide inventory. The failure simulation is the least accurate process within the simulation (accuracy 64 %). Improvements in this process might be possible using a higherquality ground shaking input map that includes topographic amplification and other nonlinear terraininduced effects. Furthermore, a better spatial estimate of the subsurface structure and the threedimensional variability of strength parameters could improve the accuracy. However, current techniques for simulation of regional failure volumes are still often not achieving the desired accuracy (Wasowski et al., 2011). Over time, such a type of modeling should be able to translate the modeling of seismic acceleration for many earthquake scenarios into a probabilistic earthquakeinduced landslide hazard map.
The runout of the coseismic landslides and the debris flow that was triggered 2 years later were predicted with high accuracy. Pre and postearthquake elevation model differences were used to validate the modeling outcomes spatially for the largest central landslide, where highaccuracy elevations models were available for both pre and postearthquake periods. In particular, the blocking of the Hongchun stream by coseismic landslide runout was simulated with good agreement (accuracy 91 %). During the simulation of the 2010 rainfall event, a debris flow was initiated by entrainment, and the landslide dam was breached. Eventually, the model simulated with relatively high accuracy the flood behavior in the second stage of the event (accuracy 92 %).
The sensitivity of the model was further analyzed by looking at ensemble simulation results. Uncertainty propagated through the different stages of the simulations. Despite this uncertainty, the ensemble models showed a high likelihood of flooding in Yinxiu. The model results show high certainty due to both the magnitude of the seismic accelerations, rainfall intensities, and the equations guiding the gravitational flows. This does not mean that the model is able to predict future events with equal reliability. However, it does indicate that, when the inherent uncertainties are taken into account, predictions from advanced spatially distributed and physically based models can provide insight and have some predictive value. Simulation results should be treated with caution. However, with validation and advancements in the understanding of the physical processes, multistage multihazard modeling might become a useful and trusted tool for decision makers in multihazard risk assessment.
The LISEM model, which was used in this research, is an opensource multihazard modeling tool. The software and source code are available from https://github.com/bastianvandenbout/LISEM (last access: 23 August 2022; bastianvandenbout, 2022), together with documentation and example datasets. The source is written in C$++$ and compiles on Windows and Linuxbased operating systems. Visit the Github repository to get started (https://github.com/bastianvandenbout/LISEM, last access: 23 August 2022; https://doi.org/10.5281/zenodo.7016476; van den Bout, 2022).
A dataset for Hongchun gully, which was studied in this work, is provided on the LISEM website (https://lisemmodel.com/, last access: 23 August 2022, or https://doi.org/10.17026/dansz6jeepm; van den Bout, 2019). For more information on how to download and use this datasets, visit the downloads page.
All authors contributed significantly to this final work. CT carried out fieldwork and laboratory tests for the model, as well as helped develop the model theory. VJ and CvW provided assistance with the hydrological modeling, interpretation and discussion of results, and writing. BvdB developed the model code, carried out the simulations, and wrote the methodology.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This research was supported by National Key Research and Development Program of China (2017YFC1501004) and National Natural Science Foundation of China (41672299). We would like to acknowledge the help of State Key Laboratory of Geohazard Prevention and Geoenvironment. In particular, the assistance and advice from Chenxiao Tang, Theo W. J. van Ash, and Xuanmei Fan were used in improving the quality of this research.
This research has been supported by the National Key Research and Development Program of China Stem Cell and Translational Research (grant no. 41672299) and the National Key Research and Development Program of China (grant no. 2017YFC1501004).
This paper was edited by Andreas Günther and reviewed by Martin Mergili and Tommaso Baggio.
Baartman, J. E., Jetten, V. G., Ritsema, C. J., and Vente, J.: Exploring effects of rainfall intensity and duration on soil erosion at the catchment scale using openLISEM: Prado catchment, SE Spain, Hydrol. Process., 26, 1034–1049, https://doi.org/10.1002/hyp.8196, 2012.
Baggio, T., Mergili, M., and D'Agostino, V.: Advances in the simulation of debris flow erosion: The case study of the Rio Gere (Italy) event of the 4th August 2017, Geomorphology, 381, 107664, https://doi.org/10.1016/j.geomorph.2021.107664, 2021.
bastianvandenbout: LISEM Project, GitHub [software], https://github.com/bastianvandenbout/LISEM, last access: 23 August 2022.
Bout, B., Van Westen, C. J., Lombardo, L., and Jetten, V.: Integration of twophase solid fluid equations in a catchment model for flashfloods, debris flows and shallow slope failures, Environ. Model. Software, 105, 1–16, https://doi.org/10.1016/j.envsoft.2018.03.017, 2018.
Cance, A., Thiery, Y., Fressard, M., Vandromme, R., Maquaire, O., Reid, M. E., and Mergili, M.: Influence of complex slope hydrogeology on landslide susceptibility assessment: the case of the Pays d'Auge, Normandy, France, In Journées Aléas Gravitaires, https://www.researchgate.net/ (last access: 23 August 2022), 2017.
Chai, H., Liu, H., Zhang, Z., and Xu, Z.: The distribution, causes and effects of damming landslides in China, Journal of the Chengdu Institute of Technology, 27, 302–307, 2000.
Chen, N. S., Zhou, W., Yang, C. L., Hu, G. Sh., Gao, Y. Ch., and Han, D.: The processes and mechanism of failure and debris flow initiation for gravel soil with different clay content, Geomorphology, 121, 222–230, https://doi.org/10.1016/j.geomorph.2010.04.017, 2010.
Chen, H. X. and Zhang, L. M.: EDDA 1.0: integrated simulation of debris flow erosion, deposition and property changes, Geosci. Model Dev., 8, 829–844, https://doi.org/10.5194/gmd88292015, 2015.
Chen, Y., Duan, H. Y., and Zhou, C. T.: The numerical calculation of highway slope stability under the influence of rainfall, in: Applied Mechanics and Materials, Trans Tech Publications Ltd, 260, 907–911, https://doi.org/10.4028/www.scientific.net/AMM.260261.907, 2013.
Chen, M. L., Hu, G. S., Chen, N. S., Zhao, C. Y., Zhao, S. J., and Han, D. W.: Valuation of debris flow mitigation measures in tourist towns: a case study on Hongchun gully in southwest China, J. Mt. Sci., 13, 1867–1879, 2016.
Chok, Y. H., Jaksa, M. B., Kaggwa, W. S., and Griffiths, D. V.: Assessing the influence of root reinforcement on slope stability by finite elements. International Journal of GeoEngineering, 6, 12, https://doi.org/10.1186/s4070301500125 2015.
Chow, V.: Open channel hydraulics, McGrawHill Book Company, Inc., New York, USA, 314–321, ISBN 070107769, 1959.
Costa, J. E. Schuster, R. L.: The formation and failure of natural dams, Geol. Soc. Am. Bull., 100, 1054–1068, 1988.
Dai, F. C., Lee, C. F., Deng, J. H., and Tham, L. G.: The 1786 earthquaketriggered landslide dam and subsequent dambreak flood on the Dadu River, southwestern China, Geomorphology, 65, 205–221, 2005.
De Baets, S., Poesen, J., Reubens, B., Wemans, K., De Baerdemaeker, J., and Muys, B.: Root tensile strength and root distribution of typical Mediterranean plant species and their contribution to soil shear strength, Plant Soil, 305, 207–226, 2008.
Dhondia, J. F. and Stelling, G. S.: Application of the one dimensionaltwo dimensional integrated hydraulic model for flood simulation and damage assessment, Hydroinformatics 2002: Proceedings of the Fifth International Conference on Hydroinformatics, Cardiff, UK, ISBN 9812387870, 2002.
Domènech, G., Fan, X., Scaringi, G., van Asch, T. W., Xu, Q., Huang, R., and Hales, T. C.: Modelling the role of material depletion, grain coarsening and revegetation in debris flow occurrences after the 2008 Wenchuan earthquake, Engineering Geology, 250, 34–44, 2019.
Ermini, L. and Casagli, N.: Prediction of the behavior of landslide dams using a geomorphological dimensionless index, Earth Surf. Proc. Land., 28, 31–47, 2003.
Evans, S. G.: The maximum discharge of outburst floods caused by the breaching of manmade and natural dams, Can. Geotch. J., 23, 385–387, 1986.
Fan, X., Tang, C. X., van Westen, C. J., and Alkema, D.: Simulating dambreach flood scenarios of the Tangjiashan landslide dam induced by the Wenchuan Earthquake, Nat. Hazards Earth Syst. Sci., 12, 3031–3044, https://doi.org/10.5194/nhess1230312012, 2012.
Flanagan, D. C., Ascough, J. C., Nearing, M. A., and Laflen, J. M.: The water erosion prediction project (WEPP) model, in: Landscape erosion and evolution modeling, Springer, Boston, MA, 145–199, https://doi.org/10.1007/9781461505754_7, 2001.
Fread, D. L.: BREACH, an erosion model for earthen dam failures, Hydrologic Research Laboratory, National Weather Service, NOAA, https://www.rivermechanics.net/images/models/breach.pdf (last access: 23 August 2022), 1988.
Gill, J. C. and Malamud, B. D.: Reviewing and visualizing the interactions of natural hazards, Rev. Geophys., 52, 680–722, 2014.
Green, W. H. and Ampt, G. A.: Studies in soil physics, J. Agric. Sci., 4, 1–24, 1911.
Hao, H., Yang, X., Huang, Y., Liu, K., and Zeng, Z.: Emergency Investigation Report of Niujuan Catchment Debris Flows, Unpublished Report, Sichuan Huadi Construction Engineering Co., Ltd., 2011, in Chinese, unpublished.
Harp, E. L. and Crone, A. J.: Landslides Triggered by the 8 October 2005, Pakistan Earthquake and Associated LandslideDammed Reservoirs, U.S. Geological Survey Openfile Report, 13, 2006–1052, 2006.
Hu, W., Dong, X. J., Xu, Q., Wang, G. H., Van Asch, T. W. J., and Hicher, P. Y.: Initiation processes for runoff generated debris flows in the Wenchuan earthquake area of China, Geomorphology, 253, 468–477, 2016.
Huang, Y. Y., Han, H., Tang, C., and Liu, S. J.: Plant community composition and interspecific relationships among dominant species on a postseismic landslide in Hongchun Gully, China, J. Mt. Sci., 14, 1985–1994, 2017.
Hungr, O., Salgado, F. M., and Byrne, P. M.: Evaluation of a threedimensional method of slope stability analysis, Can. Geotech. J., 26, 679–686, 1989.
Hungr, O., McDougall, S., and Bovis, M.: Entrainment of material by debris flows. In Debrisflow hazards and related phenomena, Springer, Berlin, Heidelberg, 135–158, https://doi.org/10.1007/3540271295_7, 2005.
Hürlimann, M., Medina, V., Bateman A., Copons, R., Altimir, J.: Comparison of different techniques to analyse the mobility of debris flows during hazard assessmentCase study in La Comella catchment, Andorra, edited by: Chen I. and Majors, DebrisFlow Hazard MitigationMechanics, Prediction and Assessment, Millpress, Netherlands, 411–422, ISBN 9789059660595, 2007.
Iverson, R. M. and Ouyang, C.: Entrainment of bed material by Earthsurface mass flows: Review and reformulation of depthintegrated theory, Rev. Geophys., 53, 27–58, 2015.
Jafarzadeh, F., Shahrabi, M. M., and Jahromi, H. F.: On the role of topographic amplification in seismic slope instabilities, Journal of Rock Mechanics and Geotechnical Engineering, 7, 163–170, 2015.
Jetten, V. G. and de Roo, A. P.: Spatial analysis of erosion conservation measures with LISEM. In Landscape erosion and evolution modeling, Springer, Boston, MA, 429–445, https://doi.org/10.1007/9781461505754_14, 2001.
Jiang, Z., Huete, A. R., Chen, J., Chen, Y., Li, J., Yan, G., and Zhang, X.: Analysis of NDVI and scaled difference vegetation index retrievals of vegetation fraction, Remote Sens. Environ., 101, 366–378, 2006.
Kappes, M. S., Keiler, M., von Elverfeldt, K., and Glade, T.: Challenges of analyzing multihazard risk: a review, Natural Hazards, 64, 1925–1958, https://doi.org/10.1007/s1106901202942, 2012.
Kalacska, M., SánchezAzofeifa, G. A., Rivard, B., CalvoAlvarado, J. C., Journet, A. R. P., ArroyoMora, J. P., and OrtizOrtiz, D.: Leaf area index measurements in a tropical moist forest: A case study from Costa Rica, Remote Sens. Environ., 91, 134–152, 2004.
Korup, O.: Geomorphic hazard assessment of landslide dams in South Westland, New Zealand: Fundamental problems and approaches, Geomorphology, 66, 167–188, 2005.
Kuriakose, S. L., Devkota, S., Rossiter, D. G., and Jetten, V. G.: Prediction of soil depth using environmental variables in an anthropogenic landscape, a case study in the Western Ghats of Kerala, India, Catena, 79, 27–38, 2009.
Lam, L. and Fredlund, D. G.: A general limit equilibrium model for threedimensional slope stability analysis. Canadian geotechnical journal, 30, 905–919, 1993.
Lehmann, P., von Ruette, J., and Or, D.: Deforestation effects on rainfall‐induced shallow landslides: Remote sensing and physically‐based modelling, Water Resour. Res., 55, 9962–9976, 2019.
Li, D., Hao, H., Ma, j., Wu, X., Yan, Z., Li, G., Wang, H., Gao, J., Liu, H., Huang, Y., Yang, X., Zeng, Z., Liu, J., Gao, L., Shen, T., Cao, N., Zhang, Y., Li, Z., Liu, K., Li, D., and Xian, Z.: Emergency Mitigation Engineering Design on the Catastrohic Hongchun Catchment Debris Flow, Yingxiu, Wenchuan County, Unpublished report, Guanghan Institute of Geological Engineering Investigation, 2011a, in Chinese, unpublished.
Li, M. H., Sung, R. T., Dong, J. J., Lee, C. T., and Chen, C. C.: The formation and breaching of a shortlived landslide dam at Hsiaolin Village, Taiwan – Part II: Simulation of debris flow with landslide dam breach, Engin. Geol., 123, 60–71, 2011b.
Liu, N., Zhang, J., Lin, W., Cheng, W., and Chen, Z.: Draining Tangjiashan Barrier Lake after Wenchuan Earthquake and the flood propagation after the dam break, Sci. China Ser. E., 52, 801–809, 2009.
Malet, J.P., Remaître, A., Maquaire, O.: Runout modeling and extension of the threatened area associated with muddy debris flows, Geomorphologie: relief, processus, Environnement, 3, 195–210, 2004.
Mergili, M., Marchesini, I., Alvioli, M., Metz, M., SchneiderMuntau, B., Rossi, M., and Guzzetti, F.: A strategy for GISbased 3D slope stability modelling over large areas, Geosci. Model Dev., 7, 2969–2982, https://doi.org/10.5194/gmd729692014, 2014a.
Mergili, M., Marchesini, I., Rossi, M., Guzzetti, F., and Fellin, W.: Spatially distributed threedimensional slope stability modelling in a raster GIS, Geomorphology, 206, 178–195, 2014b.
Mergili, M., Fischer, J.T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced opensource computational framework for the propagation and interaction of twophase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd105532017, 2017.
Mergili, M., Emmer, A., Juřicová, A., Cochachin, A., Fischer, J. T., Huggel, C., and Pudasaini, S. P.: How well can we simulate complex hydrogeomorphic process chains? The 2012 multilake outburst flood in the Santa Cruz Valley, Cordillera Blanca, Perú, Earth Surf. Proc. Land., 43, 1373–1389, https://doi.org/10.1002/esp.4318, 2018a.
Mergili, M., Frank, B., Fischer, J. T., Huggel, C., and Pudasaini, S. P.: Computational experiments on the 1962 and 1970 landslide events at Huascarán (Peru) with r. avaflow: Lessons learned for predictive mass flow simulations, Geomorphology, 322, 15–28, 2018b.
Mohadjer, S., Ehlers, T. A., Bendick, R., Stübner, K., and Strube, T.: A Quaternary fault database for central Asia, Nat. Hazards Earth Syst. Sci., 16, 529–542, https://doi.org/10.5194/nhess165292016, 2016.
Morgan, R. P. C., Quinton, J. N., Smith, R. E., Govers, G., Poesen, J. W. A., Auerswald, K., Chisci, G., Torri, D., and Styczen, M. E.: The European Soil Erosion Model (EUROSEM): a dynamic approach for predicting sediment transport from fields and small catchments, The Journal of the British Geomorphological Group, Earth Surf. Proc. Land., 23, 527–544, 1998.
Morgenstern, N. R., and Sangrey, D. A.: Methods of stability analysis, Transportation Research Board Special Report, 176, 155–171, https://trid.trb.org/view/86172 (last access: 23 August 2022), 1978.
Nash, T., Bell, D., Davies, T., and Nathan, S.: Analysis of the formation and failure of Ram Creek landslide dam, South Island, New Zealand, New Zealand, J. Geol. Geophys., 51, 187–193, 2008.
Newmark, N. M.: Effects of earthquakes on dams and embankments, Geotechnique, 15, 139–160, 1965.
O'brien, J. S., Julien, P. Y., and Fullerton, W. T.: Twodimensional water flood and mudflow simulation, J. Hydraul. Eng., 119, 244–261, 1993.
Ouyang, C., He, S., and Tang, C.: Numerical analysis of dynamics of debris flow over erodible beds in Wenchuan earthquakeinduced area, Eng. Geol., 194, 62–72, 2015.
Peng, M. and Zhang, L. M.: Breaching parameters of landslide dams, Landslides, 9, 13–31, 2012.
Pix4D, S. A.: Pix4dmapper 3.2 user manual, https://support.pix4d.com/hc/enus/articles/204272989OfflineGettingStartedandManualpdf (last access: 23 August 2022), 2017.
Pudasaini, S. P.: A general twophase debris flow model, J. Geophys. Res.Earth Surf., 117, F03010, https://doi.org/10.1029/2011JF002186, 2012.
Pudasaini, S. P. and Fischer, J. T.: A mechanical erosion model for twophase mass flows, preprint, arXiv:1610.01806, 2016.
Pyke, R., TSLOPE3: User Guide, Taga Engineering System and Software, Lafayette, https://tagasoft.com/wpcontent/uploads/2016/10/TSLOPE_Getting_started_Guide.pdf (last access: 23 August 2022), 1991.
Reid, M. E., Christian, S. B., Brien, D. L., and Henderson, S.: Scoops3D – Software to Analyze ThreeDimensional Slope Stability Throughout a Digital Landscape, In U.S. Geological Survey Techniques and Methods, 14, 218, https://doi.org/10.3133/tm14A1, 2015.
Rickenmann, D., Laigle, D., McArdell, B. W., and Hübl, J.: Comparison of 2D debrisflow simulation models with field events, Computat Geosci, 10, 241–264, 2006.
Roelvink, J. A. and Van Banning, G. K. F. M.: Design and development of DELFT3D and application to coastal morphodynamics, Oceanographic Literature Review, 11, 925 p., 1995.
Ruette, J. V., Lehmann, P., and Or, D.: Rainfalltriggered shallow landslides at catchment scale: Threshold mechanicsbased modeling for abruptness and localization, Water Resour. Res., 49, 6266–6285, 2013.
Saxton, K. E. and Rawls, W. J.: Soil water characteristic estimates by texture and organic matter for hydrologic solutions, Soil Sci. Soc. Am. J., 70, 1569–1578, 2006.
Schuster, R. L.: Landslide dams – a worldwide phenomenon, Proceedings Annual Symposium of The Japanese Landslide Society, Kansai Branch, 27 April, Osaka, 31, 38–49, 1993.
Schmidt, K. M., Roering, J. J., Stock, J. D., Dietrich, W. E., Montgomery, D. R., and Schaub, T.: The variability of root cohesion as an influence on shallow landslide susceptibility in the Oregon Coast Range, Can. Geotech. J., 38, 995–1024, 2001.
Singh, K. P. and Snorrason, A.: Sensitivity of outflow peaks and flood stages to the selection of dam breach parameters and simulation models, J. Hydrol., 68, 295–310, 1984.
Stokes, G. G.: On the effect of internal friction of fluids on the motion of pendulums, Transactions of the Cambridge Philosophical Society, 9, p. 8, https://www.biodiversitylibrary.org/item/19878#page/7/mode/1up (last access: 23 August 2022), 1850.
Swanson, F. J., Oyagi, N., and Tominaga, M.: Landslide dams in Japan, in: Landslide dams: process, risk, and mitigation New York: American Society of Civil Engineers Special Publication, edited by: Schuster, R. L., 3, 273–378, ISBN 9780872625242, 1986.
Tang, C., Zhu, J., Ding, J., Cui, X. F., Chen, L., and Zhang, J. S.: Catastrophic debris flows triggered by a 14 August 2010 rainfall at the epicenter of the Wenchuan earthquake, Landslides, 8, 485–497, 2011.
Tang, C., Jiang, Z., and Li, W.: Seismic landslide evolution and debris flow development: a case study in the Hongchun Catchment, Wenchuan area of China, in: Engin. Geol. for Society and Territory, Springer, Cham, 2, 445–449, https://doi.org/10.1007/9783319090573_72, 2015.
Tang, C., Van Westen, C. J., Tanyas, H., and Jetten, V. G.: Analysing postearthquake landslide activity using multitemporal landslide inventories near the epicentral area of the 2008 Wenchuan earthquake, Nat. Hazards Earth Syst. Sci., 16, 2641–2655, https://doi.org/10.5194/nhess1626412016, 2016.
Tang, C., Tanyas, H., van Westen, C. J., Tang, C., Fan, X., and Jetten, V. G.: Analysing postearthquake mass movement volume dynamics with multisource DEMs, Engin. Geol., 248, 89–101, 2019.
Takahashi, T., Nakagawa, H., Harada, T., and Yamashiki, Y.: Routing debris flows with particle segregation, J. Hydraul. Eng., 118, 1490–1507, 1992.
Tun, Y. W., LlanoSerna, M. A., Pedroso, D. M., and Scheuermann, A.: Multimodal reliability analysis of 3D slopes with a genetic algorithm, Acta Geotechnica, 14, 207–223, 2019.
Valiani, A., Caleffi, V., and Zanni, A.: Case study: Malpasset dambreak simulation using a twodimensional finite volume method, J. Hydraul. Eng., 128, 460–472, 2002.
Van Asch, T. W. J., Malet, J.P., van Beek, L. P. H., Amitrano, D.: Techniques, issues and advances in numerical modelling of landslide hazard. Bull. Soc. Géol. Fr., 2007, 178, 65–88, 2007.
van Beek, L. P. H.: Assessment of the influence of changes in land use and climate on landslide activity in a Mediterranean environment, Nederlandse Geografische Studies, 294, 360 p., ISBN 9068093290, 2002.
van den Bout, B.: OpenLisem MultiHazard Model and Datasets, Faculty of GeoInformation and Earth Sciences, ITC, Twente University, DANS [data set], https://doi.org/10.17026/dansz6jeepm, 2019.
van den Bout, B.: bastianvandenbout/LISEM: Beta 0.2.4 (source), Zenodo [code], https://doi.org/10.5281/zenodo.7016476, 2022.
Van Leer, B.: Towards the ultimate conservative difference scheme, a secondorder sequel to Godunov's method, J. Comput. Phys., 32, 101–136, 1979.
van Westen, C. J. and Greiving, S.: Risk Assessment for Decision Making, Environmental hazards Methodologies for Risk Assessment and Management, Dalezios, N. R., 31–94, ISBN 9781780407128, 2017.
Walder, J. S. and O'Connor, J. E.: Methods for predicting peak discharge of floods caused by failure of natural and constructed earthen dams, Water Resour. Res., 33, 2337–48, 1997.
Wang, G., Liu, F., Fu, X., and Li, T.: Simulation of dam breach development for emergency treatment of the Tangjiashan Quake Lake in China, Science in China Series E: Technological Sciences, 51, 82–94, 2008.
Wasowski, J., Keefer, D. K., and Lee, C. T.: Toward the next generation of research on earthquakeinduced landslides: current issues and future challenges, Engin. Geol., 122, 1–8, 2011.
Worden, C. B., Wald, D. J., Allen, T. I., Lin, K., Garcia, D., and Cua, G.: A revised groundmotion and intensity interpolation scheme for ShakeMap, Bull. Seismol. Soc. Am., 100, 3083–3096, https://doi.org/10.1785/0120100101, 2010.
Xie, M. W., Esaki, T., Zhou, G. Y., and Mitani, Y.: Threedimensional stability evaluation of landslides and a sliding process simulation using a new geographic information systems component, Environ. Geol., 43, 503–512, https://doi.org/10.1007/s0025400206553, 2003.
Yang, J.: Emergency Investigation Report of Shaofang Catchment Debris Flows, Unpublished report, Sichuan Institute of Geological Engineering Investigation, 2010, in Chinese, unpublished.
Zhang, N., Matsushima, T., and Peng, N.: Numerical investigation of postseismic debris flows in the epicentral area of the Wenchuan earthquake, Bull. Engin. Geol. Environ., 78, 3253–3268, https://doi.org/10.1007/s1006401813596, 2018.
Zanke, U.: Berechnung Der Sinkgeschwindigkeiten Von Sedimenten, Mitteilungen des FranziusInstituts für Wasserbau und Küsteningenieurwesen der Technischen Universität Hannover, Hannover, Germany, 46, 234–245, ISSN 03400077, 1977.
Zhu, Y.: Breach growth in claydikes, Doctoral dissertation, TU Delft, Delft University of Technology, ISBN 9090209646, 2006.