Interactive comment on “ EDDA : integrated simulation of debris flow erosion , deposition and property changes ” by

EDDA is written in FORTRAN, which can be compiled by Intel FORTRAN Compilers. The source code is enclosed as supplement files. The main subroutine is dfs.F90, which contains the numerical solution algorithm for solving the governing equations. Two input files are needed. One is edda_in.txt, which is the file for inputting material parameters and setting controlling options. EDDA is designed as the debris-flow simulation part of a cell-based model for analysing regional slope failures and debris flows, so the edda_in.txt file also includes the material parameters and controlling options for slope stability analysis. The other is inflow.txt, which is the inflow hydrograph file. Digital terrain data (e.g. surface elevation, slope gradient, erodible layer thickness) are included in separate ASCII grid files and enclosed in the data folder. Output files are stored in the result folder. Investigated variables at selected points are stored in EDDALog.txt.


Introduction
Debris flow is a flow of a sediment-water mixture driven by gravity.The mechanical triggers of debris flows can be classified into three types, namely erosion by surface runoff, transformation from landslides, and collapse of debris dams (Takahashi, 2007).Basal erosion, side erosion, and any other surficial material entrainment during the marching process entrain additional material into the flow; the final volume can be several or dozens of times of the initial volume (e.g.Hungr et al., 2005;Chen et al., 2006Chen et al., , 2012Chen et al., , 2014;;Berger et al., 2010).When the debris flow moves to a flatter area, the coarse materials can deposit gradually.During the entire movement process, not only the debris flow volume, flow velocity and flow depth change significantly, the properties of the debris flow mixture also change substantially, which in turn influence the runout characteristics.
The mechanisms of changes in debris flow mass are important and have attracted the attention of many researchers (e.g.Cannon and Savage, 1988;Takahashi et al., 1992;Hungr, 1995;Egashira et al., 2001;Iverson, 2012).Cannon and Savage (1988) and Hungr (1995) proposed one-dimensional lumped-mass models based on momentum conservation to describe the entrainment or loss of material during the movement of a debris flow.Takahashi et al. (1992) proposed a model to describe erosion and deposition based on volumetric sediment concentration and flow velocity.Researchers have also described the erosion process from a stress point of view (e.g.Medina et al., 2008;Iverson, 2012;Quan Luna et al., 2012): erosion occurs when the basal shear stress exceeds the critical erosive shear stress of the bedding material.
During the entire process of a debris flow, the debris flow properties can change significantly.The volumetric sedi-H.X. Chen and L. M. Zhang: EDDA 1.0: erosion deposition debris-flow analysis ment concentration (i.e.ratio of the solid volume to the total volume of the debris flow mixture) can increase substantially due to entrainment of solid materials (e.g.Takahashi et al., 1992;Egashira et al., 2001) or decrease due to deposition (e.g.Takahashi et al., 1992) and dilution (e.g.Pierson and Scott, 1985).Accordingly, the rheological characteristics of debris flows (e.g.yield stress and dynamic viscosity) will change with the volumetric sediment concentration, which has been observed in a large number of experiments (e.g.O'Brien and Julien, 1988;Rickenmann, 1991;Major and Pierson, 1992;Sosio and Crosta, 2009;Bisantino et al., 2010).Various rheological models have been adopted to describe debris flows, such as the laminar flow model (e.g.Takahashi, 2007), the Bingham fluid model (e.g.Fraccarollo and Papa, 2000), the Voellmy model (e.g.Medina et al., 2008), and the quadratic rheological model (Julien and Lan, 1991).
Depth-integrated models have been widely adopted to describe erosion and deposition (e.g.Takahashi et al., 1992;McDougall and Hungr, 2005;Armanini et al., 2009;Hungr and McDougall, 2009;Iverson et al., 2011;Quan Luna et al., 2012;Ouyang et al., 2015).The Mohr-Coulomb failure process was adopted to simulate bed erosion (e.g.Medina et al., 2008;Quan Luna et al., 2012).Ouyang et al. (2014) further combined the Mohr-Coulomb model and the Voellmy model to overcome the flaws of each of these two models.The changes in flow depth, flow velocity and debris mass have been accounted for in the literature.Limited attempts have also been made to consider the evolution of volumetric sediment concentration (Takahashi et al., 1992;Denlinger and Iverson, 2001;Ghilardi et al., 2001).Several key problems, however, still remain.How can one describe the various phases of a debris flow (e.g.clear water flow, hyperconcentrated flow, and fully developed debris flow) using a general rheological model?How do the properties of debris flows (e.g.volumetric sediment concentration, yield stress, viscosity) change in the erosion and deposition processes?How do these changes affect the runout characteristics of the debris flow?These problems are very important for the risk assessment of debris flows.
The objective of this paper is to develop a numerical model to consider the erosion and deposition processes and debris flow property changes during these processes.The paper is organized as follows.The methodology is introduced in Sect.2, including the problem description, governing differential equations, constitutive models, initiation of erosion and deposition, numerical solution algorithm, time stepping and numerical stability.The model is tested and verified in Sect. 3 using an analytical solution, numerical solutions, and experimental tests.A large-scale debris flow event in the Wenchuan earthquake zone is simulated as a field application in Sect. 4. The limitations of the model are indicated in Sect. 5.

Problem description
The volume of a debris flow can increase due to erosion or entrainment and decrease due to deposition.Due to changes in sediment concentration, a debris flow triggered by surface runoff may experience several flow regimes.The debris flow can evolve from a clear water flow to a hyper-concentrated flow, a fully developed debris flow, and finally a deposit on the debris fan.The erosion and deposition processes and property changes in debris flow are illustrated in Fig. 1.The debris flow entrains and incorporates materials from the channel bed if the volumetric sediment concentration, C v , is smaller than an equilibrium value, C v∞ , for the channel gradient and the shear stress is sufficiently large.Some material separates from the debris flow mixture and deposits on the channel bed when the volumetric sediment concentration is larger than an equilibrium value for the channel gradient and the flow velocity is not sufficient to take all the material.Due to erosion and deposition, debris flow properties change significantly.The changes in the volumetric sediment concentration, yield stress and viscosity are shown in Figs. 2 and 3. When it is a clear water flow, only a small amount of solid particles moves with the flow.The yield stress is negligible, and the dynamic viscosity is close to that of water.When solid materials are entrained into the flow due to erosion, the flow may evolve into a hyper-concentrated flow.The flow develops a significant yield stress, and the dynamic viscosity increases to a certain level.A debris flow can fully develop after sufficient solid materials are entrained into the flow.The yield stress and dynamic viscosity increase to relatively high levels.The debris flow decelerates when moving to a flatter area, and deposition occurs along the flow path.The volumetric sediment concentration hence decreases in the deposition process.

Governing differential equations
In this study, an integrated numerical model is developed to simulate debris flow erosion, deposition, and the induced property changes.The model is named EDDA, which stands for Erosion-Deposition Debris flow Analysis.The reference frame is defined in Fig. 1.Depth-integrated mass conservation equations (Eqs. 1, 2) and momentum conservation equations (Eqs.3, 4) are adopted to describe the movement of a debris flow:

∂(C
where h is the flow depth; t is time; v x and v y are the depthaveraged flow velocity in x and y directions, respectively; i is the erosion rate (> 0) or deposition rate (< 0); A is the rate of surficial material entrainment from collapse of bank materials or detached landslide materials; C v * and C vA are the volume fraction of solids in the erodible bed and the entrained surficial materials, respectively; s b and s A are the degree of saturation of the erodible bed and entrained surficial materials, respectively; C v is the volumetric sediment concentration of the debris flow mixture; g is the gravitational acceleration; S fx and S fy are the flow resistance slopes in x and y directions, respectively; z b is the bed elevation; and the sgn (i.e.signum) function is used to make sure the direction of the flow resistance is opposite to that of the flow direction.Similar to the two-dimensional model proposed by O'Brien et al. (1993), the governing equations above use a global coordinate system, which has been proven to simulate well flows in channels and alluvial fans (Akan and Yen, 1981;O'Brien et al., 1993).The difference is that EDDA considers changes in debris flow properties due to material entrainment and the induced momentum exchange.In Eqs.(3) and (4), the flow velocity gradient in the orthogonal direction is neglected, since very little accuracy is sacrificed by neglecting this term (Akan and Yen, 1981;O'Brien et al., 1993).In this study, erosion and deposition are investigated while surficial material entrainment is not.The bed elevation changes in the 832 H. X. Chen and L. M. Zhang: EDDA 1.0: erosion deposition debris-flow analysis erosion and deposition processes and can be expressed as (5)

Constitutive models
Various forms of rheological models can be implemented in the momentum conservation equation, which allows for the simulation of various types of flows.Several of the most widely used rheological models are introduced below to compute S f , namely the laminar flow model, the turbulent flow model, the Bingham fluid model, the Voellmy model, and the quadratic rheological model.The laminar flow model is useful to describe the movement of a fully liquefied flow, which is governed by viscous behaviour.The flow resistance slope is expressed as where µ is the dynamic viscosity, V is the absolute value of depth-averaged flow velocity, ρ is the debris flow density, and g is the gravitational acceleration.Turbulent flows with low volumetric sediment concentration are often analysed using the Manning equation: where n is the Manning coefficient.The Bingham fluid model considers both plastic and viscous behaviours.A Bingham fluid does not move if the shear stress is smaller than a threshold yield stress, but behaves as a viscous material when the shear stress exceeds the threshold.The model is expressed as where where τ 0 is the basal shear stress, and τ y is the yield stress of debris flow.The Voellmy model (Voellmy, 1955) combines the effects of frictional and turbulent behaviours: where θ is the bed slope, φ is the friction angle of the solid particles contacting the bed, and ξ is a turbulence parameter.
The quadratic rheological model proposed by Julien and Lan (1991) considers the effects of frictional behaviour, viscous behaviour, and turbulent behaviour plus the resistance arising from solid-particle contacts, which are represented by three terms as follows: where n td is the equivalent Manning coefficient, which accounts for both the turbulent behaviour and the resistance arising from solid-particle contacts and is expressed as (FLO-2D Software Inc., 2009) Since the quadratic rheological model accounts for the most comprehensive flow behaviour, it is adopted into the governing differential equations in this paper.
O'Brien and Julien (1988) proposed the following empirical relationships to estimate the yield stress, τ y , and the dynamic viscosity, µ, based on laboratory tests: where α 1 , α 2 , β 1 and β 2 are empirical coefficients.The equations describe well the changes of τ y and µ with C v when the C v value is sufficiently large.But, for very small C v values (e.g.0 for water), Eqs. ( 13) and ( 14) give α 1 and α 2 , respectively, while in reality τ y and µ are 0 and 0.001 Pa•s at 20 • , respectively, when C v is 0. In this study, a new equation is derived to estimate τ y .Assuming a hydrostatic pressure distribution within the debris flow, the effective normal stress on the inclined channel bed can be expressed as where ρ s is the density of solid particles, ρ w is the density of water, and θ is the bed slope.If suspension of particles is not considered, the yield stress at limit equilibrium can be calculated using the Mohr-Coulomb equation: The effective cohesion of the debris flow material is taken as zero in the above equation.If only the particles in contact are considered, the yield stress can be calculated by incorporating a coefficient of suspension of solid particles as follows: where C s is the coefficient of suspension of solid particles, and (1-C s ) represents the portion of solid particles that are in contact.
Three typical suspension scenarios are shown in Fig. 4: partial suspension, 0 < C s < 1; full suspension, C s = 1; and no suspension, C s = 0.The three scenarios have the same volumetric sediment concentration but different C s values.C s is in the range between 0 and 1 in most cases, which means that the solid particles are partly suspended (Fig. 4a).C s is 1 when all the solid particles are suspended and do not make contact with the bed (Fig. 4b); C s is 0 when all the solid particles are retained on the bed and no suspension occurs (Fig. 4c).In reality, some of the solid particles are in suspension due to buoyant forces, collision between solid particles, and turbulent fluid forces.The value of C s is related to particle size and flow discharge.A smaller particle size gives a larger C s value since smaller particles are more likely to suspend in water.Larger flow discharges also give larger C s values based on field observations (Alexandrov et al., 2003).
Equation ( 17) is suitable for calculating the changing yield stress especially at low solid concentrations.Equation ( 13) is suitable for calculating the changing yield stress especially at high solid concentrations and performs well on alluvial fans (O'Brien et al., 1993).Therefore, Eqs. ( 17) and ( 13) are adopted to calculate the yield stress in a confined channel with erodible materials and an unconfined flat area, respectively.The combination of the two equations overcomes the drawbacks of each of the two equations.
The values of C v in most laboratory tests range between 0.1 and 0.8 (e.g.O'Brien and Julien, 1988;Coussot et al., 1998;Schatzmann et al., 2009;Sosio and Crosta, 2009;Bisantino et al., 2010) and Eq. ( 14) can be adopted to estimate µ when C v is greater than or equal to 0.1.When C v is smaller than 0.1, Eq. ( 14) is not valid and µ is assumed to increase linearly from 0.001 Pa•s for water to α 2 e 0.1β2 .

Initiation of erosion and deposition
Erosion occurs when the bed shear stress is sufficiently large 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: where C v∞ is the equilibrium volumetric sediment concentration, and φ bed is the internal friction angle of the erodible bed.The computed C v∞ value is larger than 1 when θ approaches φ bed , and smaller than 0 when θ is larger than φ bed , both indicating an unstable or quasi-stable bed.Solid materials are difficult to be retained on such a steep slope.Hence, no erosion is expected to occur on such a slope.The erosion rate can be described approximately by the following equation: where i is the erosion rate; K e is the coefficient of erodibility, which is a soil property that describes the erosion speed; τ is the shear stress; and τ c is the critical erosive shear stress.The  2012) consider bed erosion as a Mohr-Coulomb failure process.In this study, the critical erosive shear stress can be calculated by considering the partly suspended particles at limit equilibrium using the Mohr-Coulomb equation: where c is the effective cohesion of the bed material.Based on Eqs. ( 18)-( 21), erosion occurs if τ is larger than τ c and C v is smaller than C v∞ (Fig. 1).K e and τ c can also be measured in situ using a jet index method (e.g.Chang et al., 2011).
When the debris flow moves to a flatter place, deposition occurs if C v is larger than the respective C v∞ value (Fig. 1) and the flow velocity is smaller than a critical value.Takahashi et al. (1992) proposed the following critical flow velocity for the initiation of deposition: where ρ is the density of debris flow and can be computed as follows: where d 50 is the mean particle size of the debris flow mixture.
The deposition rate can be expressed as follows (Takahashi et al., 1992): where δ d is a coefficient that describes the deposition rate; and p (< 1) is a coefficient introduced to account for the difference between the locations where actual deposition takes place in the experiment and where the velocity becomes less than V e (Takahashi et al., 1992), a value of 0.67 is recommended by Takahashi et al. (1992).

Numerical solution algorithm
The analysis domain is discretized into a grid first, with properties of each cell assigned, including the initial flow depth, the thickness and properties of the erodible soil layer, the elevation of the non-erodible layer, Manning's coefficient and so on.As shown in Eqs. ( 1) and ( 2), the changes in h and C v are governed by two effects, namely erosion or deposition, and the flow exchange among cells.The change in flow velocity is governed by four effects, namely convective acceleration, flow resistance, total head, and momentum exchange due to erosion or deposition.A volume conservation algorithm is developed to describe the changes in debris flow properties and the movement of debris flow.
As shown in Fig. 5, each cell has eight flow directions, namely four compass directions (i.e.north, east, south and west) and four diagonal directions (i.e.northeast, southeast, southwest and northwest).In each time step, the changes in h and C v at each cell due to erosion or deposition are first evaluated.After that, the flow velocity, the flow discharge, and the density of the exchange flow across each flow boundary (i.e.N , E , S , W , NE , SE , SW , NW ) of all the cells are computed; and the changes in h and C v at each cell due to the flow exchange among the cells are then evaluated.The computation of the flow velocity in each of the eight directions is independent.Therefore, Eqs.(3) and ( 4) are reduced to one equation.This type of method has been proven to be sufficient and efficient for simulating overland flows (FLO-2D Software Inc., 2009).The numerical solution algorithm is introduced step by step as follows.
1.At the beginning of each time step, the erosion rate or deposition rate is computed for each cell.The flow depth and volumetric sediment concentration of each cell are updated using Eqs.( 1) and (2) as follows: where superscript n notes the sequence of time stepping.The spatial differences of hv and C v hv in Eqs. ( 1) and ( 2) are not considered at this stage but will be accounted for in step 4. The updated bed elevation and density of flow, ρ predi , can be computed using Eqs.( 5) and ( 25), respectively.2. At each flow boundary, the average flow depth, flow density, volumetric sediment concentration and roughness of the two cells bounded at the boundary are computed.The bed slope between the two cells is defined using the gradient between the centres of the cells.
3. The new flow velocity across each flow boundary is obtained by solving the momentum conservation equation: 4. The discharge, q, across the flow boundaries is then computed, and the flow depth and density at a cell are updated as follows: where h new and ρ new are the updated flow depth and density, respectively; q b and ρ b are the discharge and density of the exchange flow across a boundary, respectively; nb is the number of flow boundaries of the cell (i.e.eight); and A cell is the area of the cell.By considering the changes of flow and density due to the flow between any two cells, the influence of the spatial differences of hv and C v hv that are not considered in step 1 is considered in this step.
5. To make the solution more robust, the average values of v n and v predi are computed and steps (1)-( 4) are repeated until the value of v predi converges.Once this is achieved, the values of v predi , h new and ρ new are assigned to v n+1 , h n+1 and ρ n+1 , respectively, and the time step moves forward.

Time stepping and numerical stability
On the one hand, the time step should be sufficiently small to ensure the numerical stability.On the other hand, the time step should be large enough to attain reasonable computational efficiency.An adaptive time stepping scheme is adopted in this research to ensure both the numerical stability and the computational efficiency, especially for cases which involve a large number of cells so that the simulation time is likely long.The algorithm for the adaptive time stepping scheme is shown in Fig. 6.Three convergence criteria are adopted in this study.The first criterion is the Courant-Friedrichs-Lewy (CFL) condition; namely, a particle of fluid should not travel more than the cell size in one time step, t.The second criterion states that the percent change of flow depth in one time step should not exceed a specified tolerant value, TOLP (h) (e.g. 10 %), which ensures that the flow depth at one cell will not change from a positive value to a negative value within one time step.If the flow moves to a cell with a zero flow depth, the second criterion cannot be satisfied and the third criterion is needed.The third criterion states that the change of flow depth in one time step should not exceed a specified tolerant value, TOL(h) (e.g.0.1 m), which makes the time step move forward even though the second criterion cannot be satisfied.The values of TOLP (h) and TOL(h) depend on required accuracy and the maximum flow depth.Larger values of TOLP (h) and TOL(h) lead to higher computation efficiency but lower accuracy.Hence, if the first criterion and either one of the last two criteria are satisfied for all cells, the time step can move forward successfully and the time step can be enlarged.Otherwise, the computation for that time step must be abandoned and the time step should be shortened until the required criteria are satisfied.

Model verification
In this section, four numerical tests are conducted to verify the performance of the proposed model.In Test 1, an analytical solution to one-dimensional debris flow is adopted to validate the performance of the model in simulating the movement of a debris flow with constant material properties.In Test 2, a two-dimensional dam-break flow problem is adopted to validate the performance of the model in sim-  ulating two-dimensional problems.In Test 3, a flume test is adopted to validate the performance of the model in describing the erosion process and material property changes.In Test 4, another flume test is adopted to validate the performance of the model in describing the movement of a debris flow considering the material property changes due to both erosion and deposition.

Test 1: analytical solution to one-dimensional debris flow with constant properties
The problem described by Liu and Mei (1989) is adopted in this test.The materials are initially retained as a triangular pile by a board and the initial profile is shown in Fig. 7.The materials start moving after the board is removed and cease moving finally due to the presence of yield stress.The final profile for the one-dimensional flow is (Liu and Mei, 1989) where a is a coefficient, and ρ and µ are 1200 kg m −3 and 5 Pa•s, respectively (Liu and Mei, 1989).Given an a value of 0.5 m, τ y is 475 Pa (Liu and Huang, 2006).The bottom surface is assumed to be smooth and the Manning coefficient is hence 0. The resistance parameter for laminar flow, K, is set to be 2500 after several trial computations.The cell size is set to be 0.1 m following Liu and Huang (2006).Since the problem is one-dimensional, there are only two flow directions and two flow boundaries for each cell.The width of the flow boundary equals the cell size.The final profiles from the analytical solution and the numerical solution are compared in Fig. 7.The maximum error is only 1 %.The model simulates the movement of the debris flow exceptionally well.

Test 2: two-dimensional dam-break water flow
A two-dimensional partial dam-breach problem reported by Fennema and Hanif Chaudhry (1987) is adopted.A sketch of the problem is shown in Fig. 8a.The computation domain is a channel 200 m in length and 200 m in width.The depth of the reservoir water is 10 m, and the depth of the tail water is 5 m.The boundary is assumed to be frictionless.The dam is assumed to fail instantaneously and the breach width is 75 m.The computation domain is discretized into a grid with cell dimensions of 2.5 × 2.5 m.The time step is kept at 0.01 s.The flow resistance slope, S f , is taken as 0 in this test as the channel is assumed frictionless.The water depth at 7.1 s after the dam breaches is shown in Fig. 8a, which agrees with the result of Fennema and Hanif Chaudhry (1987).Two points in Fig. 8a are selected for investigating the variation of water depth with time.The results from the numerical solution using EDDA and two numerical solutions by Fennema and Hanif Chaudhry (1987) are compared in Fig. 8b, which again agree reasonably well.     1 and 2, respectively.Some parameters are adopted from Takahashi et al. (1992), namely c , d 50 , ρ w , ρ s , φ, φ bed , s b , and C v * .The values of α 2 and β 2 for computing the dynamic viscosity are adopted from Bisantino et al. (2010).The mean particle size of the samples in the tests of Bisantino et al. (2010) ranged from 0.6 to 0.9 mm, which was close to the value of the materials used in the flume tests by Takahashi et al. (1992).The resistance parameter for laminar flow, K, is determined as 500, which is within the recommended range for erodible sand surfaces (FLO-2D Software Inc., 2009).The Manning coefficient is estimated to be 0.10, which is within the recommended range for open grounds with debris (FLO-2D Software Inc., 2009).The values of C s and K e are determined based on several trial calculations.The time step is kept at 0.002 s.The calculated C v values when the debris flow front moves 1, 2, 3 and 4 m are 0.23, 0.39, 0.46, 0.46, respectively.The errors when compared with the measured values are 8, 2.5, 17.9, and 15 %, respectively.The model reproduces the erosion process reasonably well.

Test 3: flume test considering changes in debris flow properties due to erosion
A sensitivity analysis is conducted to investigate the influence of K e on the erosion process.The results are summarized in Table 3.Four values of K e are adopted, namely 1 × 10 −5 , 3.5 × 10 −5 , 7.5 × 10 −5 , 3.5 × 10 −4 m 3 (Ns) −1 .The calculated C v values when the debris flow front moves 1, 2, 3 and 4 m are recorded.With the increase of K e , the erosion process becomes more intensive.For example, C v reaches 0.46 when the flow marches by only 1 m if K e is 3.5 × 10 −4 m 3 (Ns) −1 ; while C v is only 0.16 when the flow marches by 4 m if K e is 1 × 10 −5 m 3 (Ns) −1 .

Test 4: flume test considering changes in debris flow properties due to erosion and deposition
Another series of flume tests conducted by Takahashi et al. (1992) is simulated in Test 4. The experiment setup is shown in Fig. 10.In the test series, the flume width was also 1 × 10 −5 0.04 0.1 0.12 0.16 3.5 × 10 −5 0.13 0.24 0.32 0.39 7 × 10 −5 0.23 0.39 0.46 0.46 3.5 × 10 −4 0.46 0.46 0.46 0.46 10 cm.The bed layer had a length of 3.0 m and a thickness of 10 cm, which was located 5.5 m from the outlet of the flume.
A partition with a height of 10 cm was used to retain the sediment in the experiment, and the partition is assumed to have no influence on the flow.A board inclined at 5 • longitudinally (Fig. 10b) was connected to the outlet to observe the temporal variations of the shape and elevations of the debris fan.The mean particle size, d 50 , of the bed material was 3.08 mm, and the C v * value was 0.65.The surface of the bed was glued with the same material to increase the roughness.
Water was later introduced from the upstream end at a constant discharge of 600 cm 3 s −1 for 20 s to produce a debris flow.To consider the uncertainties of the sample properties, the tests were conducted six times repeatedly.In each run, the discharge at the outlet of the flume was measured.The first two runs are treated as trial runs in this study, and the results of the last four runs are adopted for comparison.
In the simulation, the flume is discretized into a grid with cell dimensions of 0.02 × 0.02 m.No flow is allowed across the flume walls.Adaptive time steps are adopted in this test following the algorithm in Fig. 6, with a minimum time step of 0.0001 s, a maximum converging time step of 0.001 s, t I of 0.0001 s, t D of 0.0001 s, TOLP (h) of 20 %, and TOL(h) of 5 cm.The soil properties and hydrological parameters used to simulate the erosion and deposition processes are summarized in Table 1 and Table 2  which are obtained following the same methods in Test 3. The value of β 1 for computing the yield stress is adopted from Bisantino et al. (2010).The value of α 1 is back calculated following the method proposed by Chen et al. (2013).Takahashi et al. (1992) found that the erosion rate was inversely proportional to the mean particle size.Since the mean particle size in Test 4 was nearly twice of that in Test 3, the coefficient of erodibility, K e , is therefore taken as one-half of that in Test 3. The Manning coefficient is determined as 0.1 for the flume covered by the saturated sand, and 0.05 for the other parts which is within the recommended range for open grounds without debris (FLO-2D Software Inc., 2009).The coefficient of deposition rate, δ d , is determined as 0.03 after trial calculations.
The computed discharges at the outlet and the measured results from the last four experiments are compared in Fig. 11.Time t = 0 in Fig. 11 denotes the time when the debris flow front reaches the outlet.As shown in Fig. 11, the model describes very well the movement of the debris flow in the confined channel.
When the debris flow moves to the flood board (Fig. 10), it decelerates and deposits gradually.The flow depth, deposit thickness, volumetric sediment concentration, and flow velocity can be monitored for all cells.If deposition occurs somewhere, the deposit thickness there will be larger than zero.The thickness of the debris fan is the sum of the flow depth and the deposit thickness.The debris fans in the experimental tests and numerical solution are compared in Fig. 12, with contours of the thickness of the debris fans.At t = 5 s, after the debris flow runs out of the outlet, the calculated runout distance of the debris fan (55 cm) is slightly smaller than the experimental result (70 cm) while the calculated width of the debris fan (70 cm) is slightly larger than the experimental result (50 cm).At t = 10 s, after the debris flow runs out, the calculated runout distance of the debris fan (75 cm) is slightly smaller than the experimental result (85 cm).The calculated width of the debris fan (85 cm) is slightly larger than the experimental result (65 cm).At the final stage, the calculated runout distance (102 cm) is slightly  position process influences the runout characteristics of the debris flow significantly.
4 Field application

Xiaojiagou debris flow event on 14 August 2010
Rainfall-induced landslides are one of the most catastrophic hazards in mountainous areas (e.g.Chen et al., 2012;Chen and Zhang, 2014;Raia et al., 2014;Zhang et al., 2013Zhang et al., , 2014)).Decisions for effective risk mitigation require hydrological and landslide analyses at the regional scale (e.g.O'Brien et al., 1993;Formetta et al., 2011;Archfield et al., 2013;Chen et al., 2013).The 2008 Wenchuan earthquake triggered numerous landslides, leaving a large amount of loose materials on the hill slopes or channels.From 12 to 14 August 2010 a storm swept the epicentre, Yingxiu, and its vicinity, triggering a catastrophic debris flow in Xiaojiagou Ravine (Fig. 14).About 1.01 × 10 6 m 3 of deposit was brought out in the form of a channelized debris flow.According to the rainfall record at Yingxiu that is 5 km from Xiaojiagou Ravine, the total rainfall amount was 220 mm in a period of 40 h (Chen et al., 2012).The debris flow was witnessed to occur at the ravine mouth at about 05:00 time) on 14 August (i.e.36 h after the storm started) and lasted about 30 min.The cumulative rainfall from the beginning of the storm to the occurrence of the debris flow was 188 mm.The runout materials of the debris flow buried 1100 m of road, blocked the Yuzixi River, formed a debris flow barrier lake and raised the river bed by at least 15 m.Interpretation of the satellite images taken before and after the debris flow reveals that the source material of this debris flow was mainly the channel colluvium (Chen et al., 2012).The deposits in the main channel marked "location of the main source material" in Figs. 14 and 15 had a volume of approximately 0.74 × 10 6 m 3 before the debris flow event and much of it had been washed away (Chen et al., 2012).Based on interpretation of satellite images and field investigations, the observed deposition zone is determined and shown in Fig. 14.

Determination of input information
The study area is divided into two domains, one for rainfallrunoff simulation and the other for debris flow runout simulation (Fig. 15).Grid systems are created within the two domains with grid sizes of 30 × 30 m for domain one and 15 × 15 m for domain two.After the Xiaojiagou debris flow, detailed field investigations and laboratory tests were conducted.The hydrological parameters for rainfall runoff simulation and debris flow runout simulation have been proposed by Chen et al. (2013) and are summarized in Table 4.
Rainfall runoff simulation is conducted first in domain one using FLO-2D (FLO-2D Software Inc., 2009).The rainfall data at Yingxiu is adopted.The runoff water would be retained by the colluvium and accumulate behind the landslide deposits, forming landslide barrier ponds.The cumulative runoff water at Section 1-1 (Figs.14, 15) can be com-puted, which is applied at Section 1-1 as the inflow hydrograph for debris flow runout simulation in domain two.Debris flow would occur when the barrier ponds breach.The source materials are assumed to be saturated before the occurrence of the debris flow.As water flows over the source material, erosion occurs if the conditions in Eqs. ( 18)-( 21) are met.Since the debris flow was witnessed to occur at the ravine mouth about 36 h after the storm started and lasted about 30 min, the cumulative runoff water at Section 1-1 in Fig. 14 at 36 h after the storm started is adopted to create the inflow hydrograph, and the surface runoff is determined to be 0.5 × 10 6 m 3 .The observed outflow hydrographs for landslide dams in the Wenchuan earthquake zone (e.g.Tangjiashan and Xiaogangjian) show a peak during the dam collapse process (Chang and Zhang, 2010).The two landslide dams breached due to overtopping erosion.Therefore, the inflow hydrograph for the Xiaojiagou debris flow is assumed to be an isosceles triangle for simplicity (Fig. 16).The duration and the peak discharge are 0.5 h and 556 m 3 s −1 , respectively.The area of the inflow hydrograph is equal to 0.5 × 10 6 m 3 .
In domain two, the source material is distributed into 329 cells with a thickness of 10 m.The internal friction angle, φ bed , for the source material is determined as 37 • according to the test results of Zhao et al. (2013).Since the source material has a very low content of silt and clay (< 2 %) according to sieving tests (Chen et al., 2012), the effective cohesion, c , is assumed to be 0 in the debris flow runout simulation.The values of d 50 , ρ s , and C v * are determined based on field and laboratory tests.The method and testing results have been reported in detail by Chen et al. (2012).Since the source material is assumed to be saturated, s b is 1.K e is determined using an empirical equation developed based on field tests in the Wenchuan earthquake zone by Chang et al. (2011): where e is the void ratio and C u is the coefficient of uniformity.The values of e and C u are 0.54 and 39, respectively, based on the sieving tests.Hence, K e is 6.6 × 10 −5 m 3 (Ns) −1 .The soil properties are summarized in  adopted in this test following the algorithm in Fig. 6, with a minimum time step of 0.01 s, a maximum converging time step of 1.0 s, t I of 0.01 s, t D of 0.02 s, TOLP (h) of 20 %, and TOL(h) of 0.5 m.

Changes in debris flow properties during the flow process
The values of the volumetric sediment concentration, C v , when the debris flow reaches Sections 2-2, 3-3, and 4-4 in Fig. 14 are computed, which are 0.13, 0.23, and 0.49, respectively.With the increase of C v , the yield stress and dynamic viscosity also increase significantly as Fig. 3 shows.The change of C v with time at Section 3-3 is shown in Fig. 17.C v is 0.23 at T = 0 when the forefront reaches the section, which can be viewed as the precursory surge in Fig. 1.Afterwards, C v increases very quickly to a peak value of about 0.5, which can be viewed as the boulder front in Fig. 1.After the boulder front passes, C v sustains at the peak value for some time (about 60 s), which can be viewed as the fully developed debris flow in Fig. 1.After that, C v decreases gradually to a lower level, which can be viewed as the hyper-concentrated flow in Fig. 1.This flow region is erosive and the bedding solid materials can be entrained.The erosion process upstream Section 3-3 lasts about 300 s.

Comparison between simulated and observed results
The volume of the simulated debris flow is about 1.0 × 10 6 m 3 .Since the sum of the inflow water (0.5 × 10 6 m 3 ) and the saturated source material (0.74 × 10 6 m 3 ) is 1.24 × 10 6 m 3 , approximately 0.24 × 10 6 m 3 of the source material is not entrained into the debris flow in the simulation.Hence, not all the source material erodes.The simulated and observed deposition zones are shown in Fig. 18a.The simulated inundated area and runout distance match the observed results reasonably well.It is noted that the debris fan front has very small depths; the fan front is the precursory surge in front of the boulder front in Fig. 1.The distribution of the maximum flow velocity is shown in Fig. 18b, which indicates that the debris flow moves very rapidly, especially in the ravine channel.Taking into account the large volume, the debris flow is very destructive, which has been observed by Chen et al. (2012).
The comparison between the simulation results and the observations suggests that the model evaluates well the debris flow volume considering the erosion process.The inundated area and the runout distance can also be predicted reasonably well.

Limitations of EDDA
The mathematical model proposed in this study has limitations due to the simplifying assumptions and approximations in the underlying theory.The main limitations are as follows.
1.The model is suitable for describing the initiation and movement of debris flows originated from runoff-driven channel bed failure or breaching of landslide dams by overtopping erosion, which has been tested in this study.
The model is also able to consider surficial material entrainment from collapses of bank material or detached landslide material as shown in the governing differential equations.But the latter capability has not been tested and further work is needed.
2. The new model is suitable for channels and flat alluvial fans, but may not be ideal for steep terrains.Such influence is not significant in a confined channel since the orthogonal gradient is small.In an unconfined flat area, the eight flow directions account for the influence to certain extent, but further work is needed to test the performance of the model.
4. Further work is needed to test the performance of some empirical equations adopted in the model.
5. The governing equations are in a depth-integrated form; hence, the particle segregation in the vertical direction cannot be considered.
6.A hydrostatic pressure distribution is assumed along the vertical direction, which affects the computed yield stress.
7. The suspension coefficient, C s , can vary with mean particle size and discharge.In Tests 3 and 4, C s is assumed to be a constant for simplicity.Further work is desired to properly determine this parameter.

Summary and conclusions
A new depth-integrated numerical model for simulating debris flow erosion, deposition, and property changes (EDDA) is developed in this study.The model considers the changes in debris flow density, yield stress, and dynamic viscosity, as well as the influences of such changes on the runout characteristics of the debris flow.
The model is unique in that it considers erosion and deposition processes, changes in debris flow mass, debris flow properties and topography due to erosion and deposition.Considering the partly suspended solid particles at limit equilibrium, the yield stress of the debris flow mixture esti-mated using the Mohr-Coulomb equation is suitable for water flows, hyper-concentrated flows, and fully developed debris flows.An adaptive time stepping algorithm is developed to assure both numerical stability and computational efficiency.
Four numerical tests have been conducted to verify the performance of the model.In Test 1, an analytical solution to one-dimensional debris flow with constant properties is adopted.Comparison between the numerical solution and the analytical solution indicates that the model simulates exceptionally well the one-dimensional movement of debris flow with constant properties.In Test 2, a two-dimensional dambreak water flow problem is adopted.The model simulates very well the two-dimensional dam-break water flow.Flume tests are simulated in Tests 3 and 4. The calculated volumetric sediment concentration at the debris flow front agrees with the experimental results reasonably well in Test 3. In Test 4, the model simulates reasonably well the erosion and deposition processes, and the movement of the multidirectional debris flows in the confined channel and the unconfined flat area in terms of the discharge hydrographs at the outlet and the time-varying geometry and elevations of the debris fan.Sensitivity analyses in Tests 3 and 4 indicate that erosion and deposition processes influence the property changes and runout characteristics significantly.
The model is also applied to simulate a large-scale debris flow in Xiaojiagou Ravine to test the performance of the model in catchment-scale simulations.The model describes well the changes in debris flow properties and estimates the volume of debris flow.Considering the deposition process, the inundated area and the runout distance are predicted properly.The model is shown to be a powerful tool for debris flow risk assessment in a large area and intended for use as a module in a real-time warning system for debris flows.

Code availability
EDDA is written in FORTRAN, which can be compiled by Intel FORTRAN Compilers.The source code is enclosed as supplement files.The main subroutine is dfs.F90, which contains the numerical solution algorithm for solving the governing equations.Two input files are needed.One is edda_in.txt,which is the file for inputting material properties and hydrological parameters and setting controlling options.EDDA is designed as the debris flow simulation part of a cell-based model for analysing regional slope failures and debris flows, so the edda_in.txtfile also includes the material properties and controlling options for slope stability analysis.The other is inflow.txt,which is the inflow hydrograph file.Digital terrain data (e.g.surface elevation, slope gradient, erodible layer thickness) are included in separate ASCII grid files and enclosed in the data folder.Output files are stored in the results folder.Investigated variables at selected points are stored in EDDALog.txt.
The Supplement related to this article is available online at doi:10.5194/gmd-8-829-2015-supplement.

Figure 1 .
Figure 1.Erosion, deposition and property changes in debris flow.

Figure 3 .
Figure 3. Changes in dynamic viscosity and yield stress with volumetric sediment concentration: (a) dynamic viscosity, (b) yield stress.

Figure 5 .
Figure 5. Eight flow directions and flow boundaries of each cell.

Figure 6 .
Figure 6.Algorithm for adaptive time stepping.I : cell I ; h(I ) i : change of flow depth of cell I during time step i; h(I ) i : flow depth of cell I during time step i; TOLP (h): tolerable value of percent change of flow depth during a time step; TOL(h): tolerable value of change of flow depth during a time step; t I : increment of time step; t D : decrement of time step.

Figure 7 .
Figure 7.Comparison of the final debris flow depth profiles from the analytical solution and the numerical solution in Test 1.

AFigure 8 .
Figure 8. Numerical solutions in Test 2: (a) water depth at 7.1 s after the dam breaches computed by EDDA; (b) comparison of the computed water depths at selected points using EDDA developed in this study and two numerical schemes reported byFennema and Hanif Chaudhry (1987).

Figure 11 .
Figure 11.Comparison of discharge hydrographs at the downstream end of the flume in Test 4.

Figure 12 .Figure 13 .
Figure 12.Comparison of the time-varying geometry and elevations of the debris fan in Test 4 from the numerical solution and the experimental tests.

Figure 14 .
Figure 14.Location of the study area and a satellite image shortly after the Xiaojiagou debris flow.

Figure 15 .
Figure 15.Grid system for rainfall runoff simulation and debris flow runout simulation.

Figure 16 .
Figure 16.Surface water inflow hydrograph for debris flow runout simulation.

3.Figure 18 .
Figure 18.Simulation results of the Xiaojiagou debris flow: (a) final shape and depth of the deposition zone, (b) the maximum flow velocity.

Table 1 .
Soil properties in Test 3, Test 4 and field application.

Table 2 .
Hydrological parameters for simulating the erosion and deposition processes in Test 3, Test 4 and field application.

Table 3 .
Results of erosion sensitivity analysis in Test 3.

Table 4 .
Hydrological parameters for rainfall-runoff and debris flow runout simulations in field application.

Table 1
Chen et al. (2013)1 , β 1 , α 2 , β 2 , K, and n are determined followingChen et al. (2013).The value of C s and δ d are determined after several trial computations.The hydrological parameters for simulating the erosion and deposition processes are summarized in Table2.Adaptive time steps are