EDDA 2.0: integrated simulation of debris flow initiation and dynamics considering two initiation mechanisms

Climate change is resulting in more frequent rainstorms and more rain-induced debris flows in mountainous areas. The prediction of likely hazard zones is important for debris flow risk assessment and management. Existing numerical methods for debris flow analysis often require the input of hydrographs at prescribed initiation locations, ignoring the initiation process and leading to large uncertainties in debris flow initiation locations, times, and volumes when applied to regional debris flow analysis. The evolution of the flowing mixture in time and space is also barely addressed. This paper presents a new integrated numerical model, EDDA 2.0, to simulate the whole process of debris flow initiation, motion, entrainment, deposition, and property changes. Two physical initiation mechanisms are modelled: transformation from slope failures and surface erosion. Three numerical tests and field application to a catastrophic debris flow event are conducted to verify the model components and evaluate the model performance. The results indicate that the integrated model is capable of simulating the initiation and subsequent flowing process of rain-induced debris flows, as well as the physical evolution of the flowing mixture. The integrated model provides a powerful tool for analysing multi-hazard processes, hazard interactions, and regional debris flow risk assessment in the future.

1 Introduction

Debris flows are one of the most catastrophic hazards in mountainous areas (e.g. Zhang et al., 2013; Raia et al., 2014), and can pose high risks to society (e.g. Tang et al., 2011; Gao et al., 2016). They are often triggered by heavy rainfall and sensitive to climate change (e.g. Wong, 2009; Lee et al., 2010). As extreme rainstorms become more frequent, coping with rain-induced debris flows becomes more critical in debris flow prone regions such as Italy, Japan, Hong Kong, and earthquake-affected areas in Sichuan, China.

During a storm, debris flows may be initiated by surface erosion, slope failures, or dam breaching (e.g. Takahashi, 2007), and enlarged during the subsequent flowing process (e.g. Iverson, 1997). The debris flow mixture finally deposits in a flatter area, while the interstice fluid still flows along the debris flow track without further material entrainment as rainfall continues. The evolution of the flowing mixture includes three phases in terms of sediment concentration: clear water flow, hyperconcentrated flow, and debris flow. The transition of the flowing mixture between any two phases occurs spatially and temporally during the whole rainfall process.

Many numerical programs have been successfully developed for debris flow analysis, such as DAMBRK (Boss Corporation, 1989), FLO-2D (O'Brien et al., 1993), DAN (Hungr, 1995), DMM (Kwan and Sun, 2006), Debris2D (Liu and Huang, 2006), FLATModel (Medina et al., 2008), MassMov2D (Beguería et al., 2009), DAN3D (Hungr and McDougall, 2009), PASTOR (Pastor et al., 2009), RAMMS (Bartelt et al., 2013), EDDA 1.0 (Chen and Zhang, 2015), DebrisInterMixing (von Boetticher et al., 2016), and AschFlow (Quan Luna et al., 2016). These programs can simulate the debris flow movement with either constant or varying properties of the flowing mixture. The entrainment and deposition processes can also be considered, such as in EDDA 1.0 (Chen and Zhang, 2015).

Until now, numerical simulation of the physical process of debris flow initiation has been largely avoided in the literature. Moreover, little attempt has been made to simulate the entire process from the initiation to the subsequent debris flow motion and deposition in an integrated manner. We address these two research gaps in this paper.

Experimental studies and field monitoring have been conducted to study the initiation mechanics of rain-induced debris flows (e.g. Johnson and Sitar, 1990; Cui, 1992; Cannon et al., 2001). A few physical models have been proposed (e.g. Takahashi, 1981; Iverson et al., 1997) to reveal the mechanisms of initiation using infinite slope stability models which are mathematically one-dimensional and statically determinate, leading to unambiguous quantitative results. However, these models do not simulate the debris flow initiation process, particularly the transformation from a slope failure to a debris flow. Statistical models have also been proposed to relate debris flow initiation to rainfall (e.g. Caine, 1980; Wieczorek, 1987; Chen et al., 2005; Godt et al., 2006; Cannon et al., 2008; Coe et al., 2008; Guzzetti et al., 2008; Baum and Godt, 2010; Berti et al., 2012; Staley et al., 2013; Zhou and Tang, 2014; De Luca and Versace, 2017a, b; Gao et al., 2017a) and other parameters such as surface runoff discharge (Berti and Simoni, 2005) or clay content (Chen et al., 2010). These models are not physically based.

Many of the existing computer programs do not simulate the initiation of debris flows. Instead, they require a predefined empirical hydrograph, created based on the estimated volumes of rainfall runoff and source materials, to initiate a debris flow, which is so-called “two-step” analysis (Fig. 1). Two-step analysis leads to large uncertainties in debris flow initiation locations, times, and volumes when applied to regional debris flow analysis. For instance, Shen et al. (2017) simulated hillslope debris flows initiated from surface erosion, in which the initiation location is artificially intervened (Fig. 1) and the slope failure mechanisms are not included. The integrated simulation of the whole process of debris flow (Fig. 1) remains an open challenge. In addition, the physical rainfall runoff and overland flow processes before the initiation of debris flows are overlooked. Currently, the study of the full evolution of the flowing mixture in time and space is limited.

Figure 1 Comparison between “two-step” simulation and integrated simulation of rain-induced debris flows.

Numerical tools have been developed for simulating single types of hazards (e.g. H. X. Chen et al., 2015; Shen et al., 2017). However, multiple types of hazards may be induced by a rainstorm (i.e. slope failures, debris flows, and flooding) (e.g. Zhang et al., 2014; Zhang and Zhang, 2017). One hazard can be the cause of another (e.g. rainfall triggers slope failures that in turn trigger debris flows). Different types of hazards can also interact with each other (e.g. several small debris flows from sub-channels can merge into a larger one). Therefore, hazard risk assessment requires hydrological, landslide, and debris flow analyses at a regional scale (e.g. Formetta et al., 2011; Archfield et al., 2013). The simulation of the complete processes of possible hazards and their interactions at a regional scale can be a powerful tool to help identify likely hazards, potentially affected areas, and elements at risk. However, the ability to numerically analyse hazard interactions is still limited (e.g. Kappes et al., 2012; Marzocchi et al., 2012). Using the existing “two-step” tools (Fig. 1) to analyse potential regional hazards can be challenging, as they involve tremendous uncertainties and it is time-consuming to conduct “two-step” analyses for each of all potential hazard locations (e.g. Chen and Zhang, 2015; Gao et al., 2016; Shen et al., 2017). Hence, the development of an integrated model for simulating multi-hazard processes and interactions (Fig. 1) is of great theoretical and practical importance.

The objectives of this paper are as follows: (1) to physically incorporate debris flow initiation into the debris flow motion simulation to enable the simulation of the whole process of rain-induced debris flows, (2) to study the full evolution of the flowing mixture in time and space during the whole rainfall process, and (3) to develop a tool to simulate multi-hazard processes and analyse hazard interactions.

2 Methodology

2.1 Strategy of modelling initiation, dynamics, and deposition of debris flows

Intense rainfall in mountainous regions could trigger debris flows from loose soil deposits on hillslopes or in channels. A conceptual model for rain-induced debris flows and likely initiation mechanisms are shown in Fig. 2. Debris flows can be initiated by three mechanisms: transformation from landslides, surface erosion, and dam breaching. Due to rainfall infiltration the hillslope gradually becomes saturated, and the soil loses its strength, causing shallow seated slope failures (Zhang et al., 2011). During a rainstorm, slope failures can occur at different times in space within a catchment. Some of the detached material may move into channels and form landslide dams, and some may directly transform into debris flows. As the surface runoff accumulates the landslide dam formed earlier in the channel may break, initiating a channelized debris flow (e.g. Liu et al., 2009; Chen et al., 2012; Peng and Zhang, 2012). At the same time, the surface runoff may cause bed erosion and initiate hillslope debris flows (e.g. Cannon et al., 2001). Some of the separate debris flows may merge in the main channel of the drainage basin, forming a larger catastrophic debris flow event (e.g. Iverson et al., 1997). The final magnitude of a debris flow could be many times that of its initial volume due to entrainment of materials along the path from additional slope failures, bed erosion, or bank collapses (e.g. Iverson et al., 2011; Chen et al., 2012; Ouyang et al., 2015). If the debris flow reaches a flat residential area downstream in the basin, it can cause severe loss of life and property.

Figure 2 Conceptual model of a rain-induced debris flow and three typical initiation mechanisms of debris flows: bed erosion, transformation from landslide, and dam breach.

Based on the conceptual model for the whole process of debris flow in Fig. 2 the strategy of the integrated model, including two debris flow initiation mechanisms (i.e. bed erosion and transformation from landslides), is shown in Fig. 3. The integrated model consists of a digital terrain module, a rainfall module, an infiltration module, an overland flow module, a slope stability module, a surface erosion module, a debris flow dynamics module, and a deposition module. The digital terrain module discretizes the study area into a grid system with geological, hydrological, and geotechnical information assigned for each cell. All the computations are based on the concept of cell. As the primary triggering factor, rainfall is simulated in the rainfall module. Water infiltration into the ground is then simulated to analyse the pore water pressure profile and compute the surface runoff. The slope stability and surface erosion are then evaluated in the slope stability module and surface erosion module, respectively. Once debris flows are initiated by the two physical mechanisms, the motion of the flowing mixture is analysed through the debris flow dynamics module. Material entrainment may occur along the flow path, incorporating solid materials from additional slope failures and surface erosion. Finally, the deposition process is assessed through the deposition module. The runout distance, inundation area, and deposition volume of the debris flows can all be assessed.

Figure 3 Framework of the integrated simulation of debris flows.

2.2 Debris flow dynamics

The core of the proposed integrated analysis is the debris flow dynamics simulation and constitutive modelling of the flowing mixture. The governing equations for debris flow dynamics describe the mixture movement and changes in debris flow properties, which are depth-integrated mass conservation equations (Eqs. 1 and 2) and momentum conservation equations (Eq. 3) (Chen and Zhang, 2015):

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(hv\right)}{\partial x}}={\displaystyle }& {\displaystyle }i[C_{\text{v* }}+(\mathrm{1}-C_{\text{v* }}\left)s_{\text{b}}\right]+A[C_{\text{vA}}+(\mathrm{1}-C_{\text{vA}}\left)s_A\right]\\ \text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial \left(C_{\text{v}}h\right)}{\partial t}}+{\displaystyle \frac{\partial \left(C_{\text{v}}hv\right)}{\partial x}}+=i{C}_{\text{v* }}+A{C}_{\text{vA}}\end{array}$$

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial v}{\partial t}}+v{\displaystyle \frac{\partial v}{\partial x}}={\displaystyle }& {\displaystyle }g\left[-\text{sgn}\left(v\right)S_f-{\displaystyle \frac{\partial (z_{\text{b}}+h)}{\partial x}}\right]\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{v\mathit{\left\{}i\right[C_{\text{v* }}+(\mathrm{1}-C_{\text{v* }})s_{\text{b}}]+A[C_{\text{vA}}+(\mathrm{1}-C_{\text{vA}})s_A\left]\mathit{\right\}}}{h}},\end{array}$$

where h is the flow depth; v is the depth-integrated flow velocity (m s⁻¹); i is the erosion rate (> 0) or deposition rate (< 0) (m s⁻¹); A is the rate of material entrainment from detached landslide materials (m s⁻¹); C_v is the volume fraction of solids in the flowing mixture; C_v*  and C_vA are the volume fraction in the erodible bed and in the entrained materials, respectively; s_b and s_A are the degree of saturation of solids in the erodible bed and in the entrained materials, respectively; S_f is the energy slope; z_b is the bed elevation (m); and the sgn (i.e. signum) function is used to ensure that the direction of the flow resistance is opposite to that of the flow direction.

One of the requirements of the integrated analysis is modelling different flowing mixtures simultaneously. The flowing mixture can be classified into three types: clear water flow, hyperconcentrated flow, and fully developed debris flow based on sediment concentration, combining grain-size distribution and particle densities (Pierson, 2005). In this study, the flowing types of mixtures are classified using the volumetric solid concentration C_v, following FLO-2D Software Inc. (2009):

1. If C_v < 0.2, the fluid mixture is deemed clear water flow which has a negligible yield stress and a dynamic viscosity similar to water.

2. If 0.2 < C_v < 0.45, a hyperconcentrated flow develops with a certain level of increased yield stress and dynamic viscosity.

3. If 0.45 < C_v < 0.6, the flowing mixture becomes a full debris flow with substantially increased yield stress and dynamic viscosity.

Therefore, a proper rheological model must involve C_v to account for the changing properties of the flowing mixture. We adopt different rheological models for different ranges of C_v to deal with this problem. For clear water flow, where C_v is less than 0.2, the energy slope S_f is based on the Manning equation. If C_v > 0.2, a quadratic rheological model developed by O'Brien et al. (1993) is used:

$$\begin{array}{}\text{(4)}& S_f={\displaystyle \frac{{\mathit{\tau }}_y}{\mathit{\rho }gh}}+{\displaystyle \frac{K\mathit{\mu }V}{\mathrm{8}\mathit{\rho }g{h}^{\mathrm{2}}}}+{\displaystyle \frac{n_{\text{td}}^{\mathrm{2}}{V}^{\mathrm{2}}}{h^{\mathrm{4}/\mathrm{3}}}},\end{array}$$

where ρ is the mass density of the flowing mixture (kg m⁻³); τ_y, μ, and n_td are the yield stress (Pa), dynamic viscosity (Pa s), and the equivalent Manning coefficient of the mixture, respectively; and K is the laminar flow resistance. n_td is expressed as follows (FLO-2D Software Inc., 2009):

$$\begin{array}{}\text{(5)}& n_{\text{td}}=\mathrm{0.0538}n{e}^{\mathrm{6.0896}{C}_{\text{v}}},\end{array}$$

where n is the Manning coefficient. The following empirical relationships are adopted to estimate τ_y and μ (O'Brien and Julien, 1988):

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_y={\mathit{\alpha }}_{\mathrm{1}}{e}^{{\mathit{\beta }}_{\mathrm{1}}{C}_{\text{v}}}\text{(7)}& {\displaystyle }& {\displaystyle }\mathit{\mu }={\mathit{\alpha }}_{\mathrm{2}}{e}^{{\mathit{\beta }}_{\mathrm{2}}{C}_{\text{v}}},\end{array}$$

where α₁, α₂, β₁, and β₂ are empirical coefficients.

2.3 Rainfall infiltration and convolution

Under heavy rainfall, excess rainwater becomes surface runoff when the rainfall intensity exceeds the infiltration capacity. In EDDA 2.0, the infiltration capacity is assumed to be the saturated permeability of the surface soil. The surface runoff process is simulated by solving the governing equations (Eqs. 1–3) and the Manning equation with i, A, and C_v equal to zero. The runoff water may cause surface erosion, or mix with landslide mass or flowing mixture, which will be described later.

Water infiltration will increase the subsurface pore water pressure, causing slope failures that are normally shallow-seated. The infiltration process is simulated in EDDA 2.0 by solving the Richards equation with a forward-time central-difference numerical solution. A non-uniform grid is created along the soil depth to enhance the accuracy of the solution near boundaries and interfaces. The integrated program calculates the instant pore water pressure profile to facilitate evaluating the slope stability of each cell at each time step.

2.4 Initiation of debris flows from slope failures

A debris flow may be initiated by the transformation of a mass flow of slope failure material at any location and at any time during a storm. The possible locations and approximate failing time can be identified in a cell-based slope stability analysis, if the topography, geology, soil properties etc. are properly defined. To consider this initiation mechanism, the slope instability evaluation must be performed over all the computational cells at each time step.

With the knowledge of real-time pore water pressure profiles provided by the infiltration module, a real-time slope instability analysis can follow. Considering that these rain-induced slope failures are shallow seated, the thickness of the failure mass is small compared to the large plan dimensions of these slopes. Therefore, an infinite slope model for two-layer soil slopes is a reasonable option to evaluate the factor of safety (F_s) (Wu et al., 2016). Following Chen and Zhang (2014), the search for the minimum F_s goes from the ground surface to the wetting front where the volumetric water content changes significantly. If the minimum F_s is smaller than one, slope failure will occur at the depth corresponding to the minimum F_s. The landslide mass is assumed to be a free-flowing mixture immediately after slope failure, with a predefined C_v value for the soil deposit and a flow depth the same as the failure depth.

2.5 Initiation of debris flows due to bed erosion

Intense rainfall can generate plentiful surface runoff, and the soil bed will erode in the runoff water. The initially clear overland flow can gradually develop into a hyperconcentrated flow and finally into a hillslope debris flow, as its C_v value increases through entrainment from bed erosion. To consider this initiation mechanism, the erosion process is analysed within each computational cell at each time step.

We consider the occurrence of erosion under the condition that the bed shear stress is equal to or larger than the critical erosive shear stress of the bed material and the volumetric sediment concentration is smaller than an equilibrium value. The equilibrium value proposed by Takahashi et al. (1992) is adopted in this study:

$$\begin{array}{}\text{(8)}& C_{\text{v}\mathrm{\infty }}={\displaystyle \frac{{\mathit{\rho }}_{\text{w}}\tan \mathit{\theta }}{({\mathit{\rho }}_{\text{s}}-{\mathit{\rho }}_{\text{w}})(\tan {\mathit{\varphi }}_{\text{bed}}-\tan \mathit{\theta })}},\end{array}$$

where ϕ_bed is the internal friction angle of the erodible bed; ρ_s is the density of soil particles (kg m⁻³); ρ_w is the density of water (kg m⁻³); and θ is the slope angle.

Many researchers have studied the relationship between the soil erosion rate and shear stress. A form of exponential expression has been used for bed erosion in the literature (e.g. Roberts et al., 1998; Z. Chen et al., 2015). More widely used is a linear function of shear stress (e.g. Graf, 1984; Hanson and Simon, 2001; Julian and Torres, 2006; Chang et al., 2011; Chen and Zhang, 2015):

$$\begin{array}{}\text{(9)}& i=K_{\text{e}}(\mathit{\tau }-{\mathit{\tau }}_{\text{c}}),\end{array}$$

where i is the erosion rate (m s⁻¹); τ is the shear stress at the soil–water interface (Pa); K_e is the coefficient of erodibility (m³ N-s⁻¹); and τ_c is the critical erosive shear stress at the initiation of bed erosion (Pa). The latter two parameters describe the erosion resistance of the bed soil and are related to soil index properties (e.g. Chang et al., 2011; Zhu and Zhang, 2016). The shear stress acting on the bed can be expressed as follows (e.g. Graf, 1984):

$$\begin{array}{}\text{(10)}& \mathit{\tau }=\mathit{\rho }gh{S}_f,\end{array}$$

where S_f is the energy slope.

2.6 Material exchange: entrainment and deposition

Material exchange occurs as a debris flow marches along its flowing path, including material entrainment (solid mass gain from outside of the flowing mixture) and deposition (solid mass loss from inside of the flowing mixture).

Entrainment from additional bed erosion or slope failure materials along its trajectory plays a significant role in debris flow volume amplification. The final volume of the debris flow deposit can be many times that of its initial volume. An excellent example of this is the 1990 Tsing Shan debris flow, which was the largest ever observed in Hong Kong. An originally small slip of 350 m³ developed into a final volume of 20 000 m³ by entraining colluvium along its flow path (King, 1996). In the integrated model, the landslide mass and surface erosion are considered as the sources of material entrainment. The slope stability and surface erosion evaluation module will be called for every computational cell at every time step; hence, the entrainment process is automatically considered once the two modules are called.

After flowing into a flatter area, deposition of some solid material will occur. Deposition is deemed to occur if the flow velocity is smaller than a critical value and C_v is larger than the equilibrium value described in Eq. 8. The deposition rate can be expressed as

$$\begin{array}{}\text{(11)}& i={\mathit{\delta }}_{\text{d}}\left(\mathrm{1}-{\displaystyle \frac{V}{p{V}_{\text{e}}}}\right){\displaystyle \frac{C_{\text{v}\mathrm{\infty }}-C_{\text{v}}}{C_{\text{v* }}}}V,\end{array}$$

where V_e is the critical flow velocity following Takahashi et al. (1992); δ_d is a coefficient of deposition rate; p (< 1) is a coefficient accounting for the location difference, and a value of 0.67 is recommended (Takahashi et al., 1992); V is the flow velocity; and C_v*  is the volume fraction of solids in the erodible bed. The deposition condition is also detailed in Chen and Zhang (2015).

2.7 Numerical scheme

The terrain is discretized into a grid of cells. Each cell is assigned with the input data, including topography, soil depth, geotechnical soil properties, rheological model parameters, and so on. There are eight flow directions in each cell: four compass directions and four diagonal directions. In each time step, the infiltration is first evaluated to compute the surface runoff and slope stability at each cell. Then changes in flow depth h and volumetric sediment concentration C_v within each cell are evaluated considering the surface runoff, slope failure mass entrainment, erosion, and deposition. This is followed by computing the flow velocity, discharge, and density along the eight flow directions of all the cells, with the averaged surface roughness and slope between two cells computed. The changes in h and C_v due to the flow exchange are finally evaluated at each cell.

After all the computations have been completed in each time step, numerical stability criteria are checked for each cell to limit the time step and avoid surging while allowing for large time steps. Three convergence criteria are adopted:

1. The Courant–Friedrichs–Lewy (CFL) condition, with the physical interpretation that a particle of fluid should not travel more than the cell size in one time step (Fletcher, 1990), is mostly used in explicit schemes. The time step is limited by

   $$\begin{array}{}\text{(12)}& \mathrm{\Delta }t\le C\mathrm{\Delta }x/(\mathit{\beta }V+c),\end{array}$$

   where C is the Courant number (C is not smaller than or equal to one); m is a coefficient (5∕3 for a wide channel); and c is the computed wave celerity.

2. The percent change of flow depth in one time step should not exceed a specified tolerant value, TOLP(h).

3. The change in flow depth in one time step should not exceed a specified tolerant value, TOL(h), which is applied when the flow moves to a cell with zero flow depth.

Adjusting these three criteria, the computational time and accuracy could reach a good balance. If all the numerical stability criteria are successfully satisfied, the time step can be increased for the next computational cycle. Otherwise the time step is reduced and the computation restarted. The volume conservation is computed at the end of each time step for the inflow, outflow, grid system storage, and infiltration loss.

Figure 4 A satellite image of the study area taken shortly after the Xiaojiagou debris flow on 14 August 2010.

3 Model verification

The previous version, EDDA 1.0 (Chen and Zhang, 2015), passed several verification tests including debris flow dynamics, erosion, and deposition. In this new version of integrated analysis, the new modules for surface runoff, coupled infiltration, and slope stability analysis, and the integrated program require further verification. The response of Xiaojiagou Ravine during a rainstorm in August 2010 is used to verify the new modules. The in situ conditions shortly after the 2010 Xiaojiagou debris flow event are shown in Fig. 4. The Xiaojiagou Ravine has an area of 7.84 km². The elevation of the ravine ranges between 1100 and 3200 m. The hillslopes within the ravine are very steep with an average slope angle of 46^∘. There is one main drainage channel and four branches within the Xiaojiagou Ravine. The loose soil deposits on the hillslopes and channels of the ravine before the debris flow event are identified based on field investigations and interpretation of a satellite image (e.g. Chen and Zhang, 2014). The rainstorm process triggering the catastrophic Xiaojiagou debris flow is presented in Fig. 5. The rainstorm lasted about 40 h with a total precipitation of 220 mm. In this study, the rainfall is assumed to be uniformly distributed. Spatially variable rainfall data can be used when a large area is considered, as spatial rainfall variation and the potential of triggering landslides are correlated (Gao et al., 2017b).

Figure 5 Rainfall process of the August 2010 rainstorm.

First the performance of the rainfall runoff module of the integrated program is compared with a commonly used program FLO-2D (FLO-2D Software Inc., 2009). The infiltration module is then checked against an analytical solution under steady rainfall. The slope stability analysis is verified by comparing it with the landslide satellite image and the computation results by Chen and Zhang (2014). Finally, the performance of the integrated model is checked against the 2010 Xiaojiagou debris flow event in Sect. 4.

Figure 6 Comparison of the maximum surface runoff flow depths and flow velocities simulated using FLO-2D (a, b) and EDDA 2.0 (c, d).

3.1 Verification test 1: rainfall runoff

The same input data are used in EDDA 2.0 and FLO-2D, including the digital elevation model, the Manning coefficient (n= 0.3), the limiting Froude number (L_f= 0.8), the saturated permeability of the surface soil (k_st= 3.6 mm h⁻¹ or 10⁻⁶ m s⁻¹), and the rainfall data (Fig. 5). Other hydrological parameters such as the soil porosities used in FlO-2D are adopted following Chen et al. (2013) and Shen et al. (2017).

The results from the two programs are compared in Fig. 6, including the distributions of the maximum flow depth and flow velocity. The result from FlO-2D (Fig. 6a and c) differ only slightly from those of EDDA 2.0 (Fig. 6b and d). During the rainstorm process, the maximum flow depth computed by FLO-2D is 3.2 m, while that computed by EDDA 2.0 is 3.4 m. The outflow hydrographs recorded at the mouth of the ravine of the two programs are shown in Fig. 7. The computed overall discharge processes from both programs are very close.

3.2 Verification test 2: infiltration process and resulting pore water pressure changes

Before applying the infiltration module to compute the pore water pressure profiles under the actual rainfall event, four cases of infiltration under steady rainfall are adopted to verify the infiltration module. The results are compared with those from an analytical solution by Srivastava and Yeh (1991) and Zhan et al. (2013). The scenario of two-layer soil is considered, which is also used in the field application. Table 1 presents the input parameters for the four cases. Four combinations are set up to represent likely in situ conditions. The results from the numerical infiltration module and the analytical solution are compared in Fig. 8. For all the four cases, the module performance is satisfactory.

Table 1 Parameters used in the infiltration module verification.

Notes: α is the  constitutive parameter; θ_s is the  saturated water content; θ_r is the  residual water content; k_s is the  saturated permeability; q_a is the  antecedent rainfall intensity; q_b is the  rainfall intensity for time greater than zero; and γ is the  slope angle. Parameters α, θ_s and θ_r are used in the constitutive relations between the hydraulic conductivity and moisture content and the pressure head (Srivastava and Yeh, 1991).

Table 2 Properties of four types of superficial materials.

Notes: c^′ is the  true cohesion of soil; ϕ^′ is the  friction angle of soil; γ_sat is the  unit weight of solid particles; and K_s is the  saturated permeability of soil.

Table 3 Soil properties for debris flow simulation.

Notes: d₅₀ is the  mean grain size; ρ_s is the  density of solid particles; C_v*  is the  volume fraction of solids in the erodible bed; s_b is the  degree of saturation of the erodible bed; τ_c is the  critical erosive shear stress; and K_e is the  coefficient of erodibility.

Table 4 Constitutive (rheological) parameters for debris flow simulation.

Notes: α₁, β₁ are the  empirical coefficients for calculating τ_y; α₂, β₂ are the  empirical coefficients for calculating μ; K is the  laminar flow resistance coefficient; δ_d is the  deposition coefficient; and n is the  Manning coefficient.

Figure 7 Comparison of the outflow hydrographs at the ravine mouth using FLO-2D and EDDA 2.0.

Figure 8 Pore water pressure profiles at various times: (a) case 1; (b) case 2; (c) case 3; and (d) case 4.

3.3 Verification test 3: slope stability analysis

The 2008 Wenchuan earthquake triggered over 50 000 landslides within the earthquake region, leaving a large amount of loose materials on hillslopes and in channels (Fig. 4). These materials became the source of numerous post-earthquake rain-induced landslides and debris flows. Presently, nearly 80 % of such material remains in the mountain regions, posing great potential threats (Zhang et al., 2016). EDDA 2.0 is used to reproduce the slope failures under the rainstorm conditions from August 2010 (Fig. 5) by Chen and Zhang (2014), who evaluated the slope stability of a 164.5 km² area near the epicentre. All the parameters in this research are the same as those in Chen and Zhang (2014), with the only difference being that the focus area in this study is Xiaojiagou Ravine (Fig. 4). The loose soil deposits are assumed to be two layers. Given the same parameters such as the topography, layer thicknesses, and soil properties, the unstable cells when rainfall terminates are computed using the slope failure module. Comparing the simulation results with the observation (Fig. 9), the computed unstable cells generally fall upon the landslide scars formed during the rainstorm event. Moreover, the results are compared with those from Chen and Zhang (2014), which have been verified using the confusing matrix method (e.g. Van Den Eeckhaut et al., 2006). It is found that the results of the two separate analyses are very similar. The computed total scar area is 4.42 × 10⁵ m², comparing well with 5.20 × 10⁵ m² from the satellite image; the difference is 15 %. It is concluded that the proposed slope stability module performs reasonably well.

4 Field application

4.1 Xiaojiagou debris flow on 14 August 2010

A heavy rainstorm swept the epicentre of the event, Yinxiu town, and its vicinity. The rainstorm lasted about 40 h from 12 to 14 August 2010, delivering a total of about 220 mm of precipitation (Fig. 5). A catastrophic debris flow was triggered by the storm in Xiaojiagou Ravine (Fig. 4). The debris flow was witnessed at the ravine mouth at approximately 05:00 LT on 14 August and lasted about 30 min. Roughly 1.17 × 10⁶ m³ of the soil deposit was brought out of the Xiaojiagou Ravine mouth in a form of a channelized debris flow. The runout material deposited in front of the mouth, burying 1100 m of Province Road 303 (PR303), blocking Yuzixi River, forming a debris flow barrier, and raising the river bed by at least 15 m.

Figure 9 Computed unstable cells vs. landslide scars on the satellite image.

Figure 10 Simulation results of the Xiaojiagou debris flow: (a) final shape and depth of the erosion zone; (b) maximum flow velocity.

Figure 11 Comparison of the simulated and observed deposition zones: (a) simulation result; (b) enlarged view of the observed deposition area (Chen and Zhang, 2015).

Figure 12 Outflow hydrograph and changes in C_v at the Xiaojiagou Ravine mouth during the simulation period.

Figure 13 Distributions of C_v at different times during the storm event: (a) clear water flow; (b) initiation of debris flows; (c) channelized debris flows; and (d) post hyperconcentrated/clear water flow.

4.2 Input information

In EDDA 1.0, the study area has to be divided into two domains for rainfall runoff simulation and debris flow runout simulation respectively. However, in the integrated simulation by EDDA 2.0, only one grid of 9500 30 × 30 m cells is created (Fig. 1). After the Xiaojiagou debris flow, detailed field investigations and laboratory tests were conducted (Chen et al., 2012), as well as numerical back analysis (Chen et al., 2013). The study area is divided into four zones by satellite interpretation: bare soil, vegetated soil, bedrock, and river bed (Chen and Zhang, 2014). The soil properties of each zone and the constitutive (or rheological) parameters used in the integrated simulation are determined following EDDA 1.0 (Chen and Zhang, 2015), shown in Tables 2–4. The erosion resistance parameters τ_c and K_e of the soils are determined using the empirical equations based on field tests in the Wenchuan earthquake zone (Chang et al., 2011):

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{c}}=\mathrm{6.8}{\text{PI}}^{\mathrm{1.68}}{P}^{-\mathrm{1.73}}{e}^{-\mathrm{0.97}}\text{(14)}& {\displaystyle }& {\displaystyle }{K}_{\text{e}}=\mathrm{0.020075}{e}^{\mathrm{4.77}}{C}_{\text{u}}^{-\mathrm{0.76}},\end{array}$$

where e is the void ratio; “PI” is the plasticity index; P is the fines content (< 0.063 mm); and C_u is the coefficient of uniformity. These four soil properties are determined to be 1.05, 18, 14, and 2000, respectively, according to Chang et al. (2011). Therefore, τ_c and K_e are estimated to be 8.7 Pa and 7.85 × 10⁻⁸ m³ N-s⁻¹, respectively. Uncertainties in the soil properties may also be included when considering soil spatial variability (Xiao et al., 2017).

4.3 Integrated simulation results

We examine the final output of the integrated simulation first. Erosion plays an important role in the volume magnification of debris flows. The final erosion depths in the eroded areas are shown in Fig. 10a. The most eroded areas during the Xiaojiagou debris flow event were in channels, where a huge amount of loose solid material was present (Chen et al., 2012). Loose deposits on the hillslopes also eroded after the landslide bodies detached from their original locations and slid down the slopes. The distribution of the eroded areas reflects that the debris flows were initiated from both slope failures and surface erosion, then developed along the channels by further erosion and entrainment of the slope failure materials; these are the two mechanisms considered in the integrated model. The distribution of the maximum flow velocity is shown in Fig. 10b, with the maximum value being 9.5 m s⁻¹. This value is very close to that from EDDA 1.0 (9.1 m s⁻¹). The slightly larger value of flow velocity from EDDA 2.0 is attributed to the consideration of the extra surface runoff within domain two created when using EDDA 1.0 (Fig. 1). The maximum velocity occurs in the ravine channels, indicating that the debris flow moves very rapidly.

The simulated and observed deposition areas are compared in Fig. 11. It is seen that the simulation results (Fig. 11a) match the observation (Fig. 11b) reasonably well. The simulated deposition depth is approximately 20 m, very close to that of the observed thickness of the deposit fan during the field investigations. The total volume of the observed deposition fan is about 1.17 × 10⁶ m³, while the simulated deposition volume of the debris flow is 0.9 × 10⁶ m³. The integrated model evaluates a smaller debris flow volume and the difference is about 23 %. The main uncertainty arises from the slope failure module and surface erosion module.

The changes in the volumetric sediment concentration C_v and the discharge hydrograph at leg 1-1 (see Fig. 4) are recorded during the simulation of the whole rainfall process, shown in Fig. 12. The integrated model simulates two peaks in the discharge process throughout the rainfall with a precursory boulder front arriving in advance. At around 12 h, the value of C_v increases very quickly to a peak value of 0.6, indicating the arrival of the debris flow. Afterwards, C_v decreases, which can be viewed as a hyperconcentrated flow or a clear water flow after the debris flow passes. Another large debris flow surge is simulated at around 32 h with the same pattern as the first one. The debris flow passes through leg 1-1 (see Fig. 4) first and continues to develop for some time. After most of the solid materials are carried away by the debris flow surge, the flow at leg 1-1 becomes a hyperconcentrated flow. The flowing mixture then gradually becomes a clear water flow as the rainwater continues to generate surface runoff without further material entrainment. The integrated simulation is capable of simulating multiple debris flow surges and the changes in the flowing mixture properties throughout a rainfall event.

To demonstrate the evolution of the flowing mixture within the drainage basin, the distributions of C_v at four snapshots during the storm are shown in Fig. 13. The recording times of these four figures span a complete evolution cycle, i.e. clear water flow (Fig. 13a), debris flow initiation (Fig. 13b), debris flow motion (Fig. 13c), and hyperconcentrated flow/clear water flow (Fig. 13d). This evolution cycle could occur within the basin several times in different branch channels, which can be captured by the integrated model.

5 Limitations of EDDA 2.0

We have successfully extended the “two-step” debris flow simulation to an integrated simulation of the whole process of rain-induced debris flows. However, there are still limitations in the underlying assumptions and simplifications:

1. EDDA 2.0 considers the initiation of debris flows from the transformation of slope failures and surface erosion. However, the initiation from dam breaching has not yet been tested.

2. The studies consider material entrainment from surface erosion and slope failure detachment, but the entrainment from bank failures can only be considered using an empirical rate, instead of via a three-dimensional physical model.

3. The governing equations are in a depth-integrated form; hence, particle segregation in the vertical direction cannot be considered.

4. The rheological models for the hyperconcentrated flow, fully developed debris flow, and slope failure mass flow need further study. The slope failure mass movement is particularly critical for estimating the transformation rate from a slope failure to a debris flow.

6 Summary and conclusions

A new integrated simulation model is developed for simulating rain-induced debris flow initiation, motion, entrainment, deposition, and property changes. The model is unique in that it simulates the whole process of rain-induced debris flow evolution and two physical initiation mechanisms (i.e. transformation from landslides and surface erosion). Previous “two-step analysis” with an assumed inflow hydrograph and an inflow location can now be conducted at once, scientifically, and without subjective assumptions.

Three numerical tests have been conducted to verify the performance of the newly added modules of the integrated model. The Xiaojiagou Ravine landslides and debris flows triggered by the rainstorm in August 2010 were used as a verification case. In test 1, the rainfall runoff simulation by EDDA 2.0 was compared to FLO-2D. The simulation results from the two models are very close, which indicates that EDDA 2.0 simulates rainfall runoff well. In test 2, an analytical solution for evaluating the pore water pressure profile under infiltration is adopted. Comparison between the model solution and the analytical solution indicates that the integrated model evaluates the infiltration process well. The regional slope stability within the study area under the same rainstorm was evaluated using the integrated model in test 3. The computed unstable cells compare well with the observations from the satellite image and the results from previous studies.

The new integrated model was finally applied to reproduce the Xiaojiagou debris flow event. The model can simulate the entire evolution process of rain-induced debris flows, and estimates the volume, inundated area, and runout distance of the debris flow reasonably well . It is concluded that the new integrated debris flow simulation model, EDDA 2.0, is capable of (1) simulating the whole process of rain-induced debris flow from debris flow initiation to post-initiation debris flow motion, entrainment and deposition, and (2) tracing the evolution of the flowing mixture in time and space during the whole process of rainfall. The integrated model will serve as a powerful tool for analysing multi-hazard processes and hazard interactions, and the assessment of regional debris flow risks in the future.