GIS-based numerical simulation of Amamioshima debris flow in Japan

Abstract

Debris flow is a rapid flow which could lead to severe flooding with catastrophic consequences such as damage to properties and loss of human lives. It is important to study the movement of debris flow. Since during a debris flow process, the erosion and deposition processes are important, the no entrainment assumption is not acceptable. In this study, first we considered the debris flow as equivalent fluid and adopted the depth-averaged govern equations to simulate the movements and evolution of river bed. Secondly, the set of partial differential equations was solved numerically by means of explicit staggered leap-frog scheme that is accurate in space and time. The grid of difference scheme was derived from GIS raster data. Then the simulation results can be displayed by GIS and easily used to form the hazard maps. Finally, the numerical model coupled with GIS is applied to simulate the debris flow occurred on Oct. 20th, 2010, in Amamioshima City, Japan. The simulation reproduces the movement, erosion and deposition. The results are shown to be consistent with the field investigation.

Keywords

debris flow / numerical simulation / GIS / movement / erosion / deposition

Introduction

Debris flow is a rapid flow which could lead to severe flooding with catastrophic consequences such as damage to properties and loss of human lives. For example, a giant debris flow burst on Aug. 8th, 2010 in Zhouqu City, China, 1765 people were killed. More than 5500 houses were inundated and the total economic losses reached to 212 million CNY (Chinese Yuan) [1]. In the “Minamarta debris flow disaster” in 2003, 21 people were killed by the debris flow. 240 houses, many office buildings and the Hakata-eki subway station were inundated with muddy water [2]. Therefore, debris flow disasters have been recognized as a critical problem facing the world today and increasing attention has been focused on the study of debris flow.

Debris flow can be considered of sediment and water mixture in a manner as if it was a flow of continuous fluid driven by gravity, and it attains large mobility from the enlarged void space saturated with water or slurry. Associated with debris flow movement there is often significant erosion and deposition that can dramatically change the channel bed. At the end of an alluvial fan, where the slope of the bed decreases significantly, debris flow slows significantly, depositing large amounts of sediments [3].

Due to the complexity of the debris flow process, a number of models were developed to simulate the flow behavior. These models can be classified as: single-phase models [4-8], and two-phase models [9-17]. One of these processes which is still difficult to include inside a numerical model and that plays an important role in debris flow is the erosion and deposition. The erosion and deposition mechanisms are able to change significantly the movement of debris flow, through rapid changes of the volume and its rheological behavior [18,19]. Models using a constant volume cannot yield accurate forecasts of debris flow characteristics [20,21]. After the failure at the source zone, the erosion materials may accumulate 10-50 times in volume with respect to initially mobilized mass [22]. Some effects have already been made to quantify the erosion processes and entrained volumes, trying to propose a physical explanation for the extreme bulking rates [14,19,23,24]. However, the introduction of erosion and deposition in the models requires additional parameters, which will certainly complicate calculations even further.

For hazard mapping and risk assessment, the geographic information system (GIS) has been recognized as a useful tool to process spatial data and to display results. The GIS-based approaches in assessing debris flow hazards were reported in recent years [25-28]. Numerical simulation by incorporating GIS is important prediction and analysis tools. One significant advantage of numerical simulation coupled with GIS is that the grid networks for simulation can be extracted from GIS raster data and all the calculated results can be displayed in GIS and used to the hazard mapping directly.

In this study, following a derivation of the debris flow model considering both erosion and deposition processes, we developed a debris flow simulation program incorporating with GIS by deriving the computation networks from GIS raster data. This coupled model was used to simulate a well-documented Amamioshima debris flow in Japan. Comparing the results with the detailed field investigation, it shows that the present coupled model is rational and effective.

Description of numerical model

The derivation of governing equations

Modeling debris flows require rheological models (or constitutive equations) for solid-liquid mixtures. Identification of an appropriate rheology has long been regarded as the key to interpretation, modeling, and prediction of debris-flow behavior, and debates about the most suitable rheological formula have persisted for several decades. The rheological property of a debris flow depends on a variety of factors, such as the water concentration, the solid concentration, cohesive properties of the fine material, particle size distribution, particle shape and grain friction [29]. It is well know that water is the main contributor to rainfall-induced debris flow initiation and the role played by the water in such flows will affect the rheological property. At the same time, field observations and video recordings of debris flows have provided clear evidence that no unique rheology is likely to describe the range of mechanical behaviors exhibited by poorly sorted, water-saturated debris. Instead, apparent rheology appear to vary with time, position, and feedbacks that depend on evolving debris-flow dynamics [30]. Therefore, in this paper, the debris and water mixture is assumed to be uniform continuous, incompressible, unsteady flow. The flow is governed by the following forms of the continuity and Navier-stokes equations.

(1){\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial w}{\partial z}}=\mathrm{0},

(2){\rho }_{\mathrm{m}}\left({\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}+w{\displaystyle \frac{\partial u}{\partial z}}\right)=-{\displaystyle \frac{\partial P}{\partial x}}+\mu \left({\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial z^{\mathrm{2}}}}\right),

(3){\rho }_{\mathrm{m}}\left({\displaystyle \frac{\partial v}{\partial t}}+u{\displaystyle \frac{\partial v}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}+w{\displaystyle \frac{\partial v}{\partial z}}\right)=-{\displaystyle \frac{\partial P}{\partial y}}+\mu \left({\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial z^{\mathrm{2}}}}\right),

(4){\rho }_{\mathrm{m}}\left({\displaystyle \frac{\partial w}{\partial t}}+u{\displaystyle \frac{\partial w}{\partial x}}+v{\displaystyle \frac{\partial w}{\partial y}}+w{\displaystyle \frac{\partial w}{\partial z}}\right)={\rho }_{\mathrm{m}}g-{\displaystyle \frac{\partial P}{\partial z}}+\mu \left({\displaystyle \frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}w}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}w}{\partial z^{\mathrm{2}}}}\right),

in which, u, v, w are velocity components in the x, y, z directions; {\rho }_{\mathrm{m}} is equivalent density of the debris and water mixture, Pis pressure; \muis dynamic viscosity; gis gravitational acceleration; and t is time.

A key step in simplifying Eqs. (2)-(4) involves scaling that is similar but not identical to the well-known shallow water or Saint-Venant scaling [31]. As described by Savage and Iverson (1997) [18], Savage and Hutter (1989) [32], and Gray et al. (1999) [33], which imply thatu,v>>w. In Eq. (4), all terms are small relative to the gravitational acceleration, only the pressure gradient remains to balance it. Therefore Eq. (4) can be rewritten as follows:

(5){\rho }_{\mathrm{m}}g-{\displaystyle \frac{\partial P}{\partial z}}=\mathrm{0.}

Using the notation and coordinate systems given in Fig. 1, Integration of (5) with respect to z from the elevation of the bottom basez=Z_b to the free surface of the flow atz=Hyields,

(6)P={\rho }_{\mathrm{m}}g(H-Z_b).

The final step in further simplifying the continuity and momentum equations of motion is to adapt depth averaging to eliminate explicit dependence on the coordinate normal to the bed, z. Depth averaging requires decomposing the Eqs. (1), (2) and (3) into component equations in locally defined x-y-z-orthogonal directions, then integrating each component equation from the base of the flow at z=Z_b to the free surface of the flow at z=H, using the Leibnitz rule to interchange the order of integration and differentiation. Herein, we define the depth-averaged velocities are as follows:

(7)\overline{u}={\displaystyle \frac{\mathrm{1}}{h}}{\displaystyle {\int }_{Z_b}^{H}u\mathrm{d}z},\overline{v}={\displaystyle \frac{\mathrm{1}}{h}}{\displaystyle {\int }_{Z_b}^{H}v\mathrm{d}z}.

In Eq. (7), and equations hereinafter, overbars denote depth-averaged quantities defined by integrals similar to (7).

The continuity Eq. (1) is integrated through the flow depth,

(8){\displaystyle {\int }_{Z_b}^H\left({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial w}{\partial z}}\right)\mathrm{d}z}={\displaystyle \frac{\partial }{\partial x}}{\displaystyle {\int }_{Z_b}^{H}u\mathrm{d}z}+{\displaystyle \frac{\partial }{\partial y}}{\displaystyle {\int }_{Z_b}^{H}v\mathrm{d}z}-{\left(u{\displaystyle \frac{\partial H}{\partial x}}+v{\displaystyle \frac{\partial H}{\partial y}}-w\right)|}_{z=H}+{\left(u{\displaystyle \frac{\partial H}{\partial x}}+v{\displaystyle \frac{\partial H}{\partial y}}-w\right)|}_{z=Z_b}.

Using the kinematic boundary conditions at the base and at the free surface, we obtain,

(9)u{\displaystyle \frac{\partial z}{\partial x}}+v{\displaystyle \frac{\partial z}{\partial y}}-w=E\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}z=Z_b,

(10){\displaystyle \frac{\partial z}{\partial t}}+u{\displaystyle \frac{\partial z}{\partial x}}+v{\displaystyle \frac{\partial z}{\partial y}}-w=\mathrm{0}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}z=H,

as h=h(x,y,t)=H-Z_b is the flow depth, we get:

(11){\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial (H-Z_b)}{\partial t}}.

Inserting (7), (9), (10), (11) into (8), we get the continuity equation of debris flow:

(12){\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial M}{\partial x}}+{\displaystyle \frac{\partial N}{\partial x}}=E.

Integrating the left-hand side of the Eq. (2):

(13)$$\begin{gathered} {\displaystyle {\int }_{Z_b}^H\left({\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}+w{\displaystyle \frac{\partial u}{\partial z}}\right)\mathrm{d}z}={\displaystyle {\int }_{Z_b}^H\left({\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial (uu)}{\partial x}}+{\displaystyle \frac{\partial (uv)}{\partial y}}+{\displaystyle \frac{\partial (uw)}{\partial z}}\right)\mathrm{d}z} \\ ={\displaystyle \frac{\partial }{\partial t}}{\displaystyle {\int }_{Z_b}^{H}u\mathrm{d}z}+{\displaystyle \frac{\partial }{\partial x}}{\displaystyle {\int }_{Z_b}^H{u}^{\mathrm{2}}\mathrm{d}z}+{\displaystyle \frac{\partial }{\partial y}}{\displaystyle {\int }_{Z_b}^{H}uv\mathrm{d}z}-u{\left({\displaystyle \frac{\partial z}{\partial t}}+u{\displaystyle \frac{\partial z}{\partial x}}+v{\displaystyle \frac{\partial z}{\partial y}}-w\right)|}_{z=Z_b}^{z=H}. \end{gathered}$$

With

(14)\begin{array}{l}{\displaystyle {\int }_{Z_b}^H{u}^{\mathrm{2}}dz}={\alpha }_{\mathrm{1}}h{u}^{\mathrm{2}}={\alpha }_{\mathrm{1}}uM,\hfill \\ {\displaystyle {\int }_{Z_b}^{H}uvdz}={\alpha }_{\mathrm{2}}huv={\alpha }_{\mathrm{2}}vM\hfill \end{array}.

Two ad hoc coefficients {\alpha }_{\mathrm{1}} and {\alpha }_{\mathrm{2}} have been introduced. As was noted by Savage and Hutter (1989) [32], values of {\alpha }_{\mathrm{1}} and {\alpha }_{\mathrm{2}} in Eq. (14) that deviate from unity provide information about the deviation of the vertical velocity profile from uniformity. If a debris flow velocity profile is reasonably blunt, {\alpha }_{\mathrm{1}}={\alpha }_{\mathrm{2}}=\mathrm{1}. For parabolic velocity profile and debris flows with no basal sliding, {\alpha }_{\mathrm{1}}={\alpha }_{\mathrm{2}}=\mathrm{1.2}, and for a stone-type debris flow on a rough inclined plane, {\alpha }_{\mathrm{1}}={\alpha }_{\mathrm{2}}=\mathrm{1.25} [14].

We may write for (13) using\alpha ={\alpha }_{\mathrm{1}}={\alpha }_{\mathrm{2}},

(15){\displaystyle {\int }_{Z_b}^H\left({\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}+w{\displaystyle \frac{\partial u}{\partial z}}\right)\mathrm{d}z}={\displaystyle \frac{\partial M}{\partial t}}+\alpha {\displaystyle \frac{\partial \left(\overline{u}M\right)}{\partial x}}+\alpha {\displaystyle \frac{\partial \left(\overline{v}M\right)}{\partial y}}.

Furthermore, we integrate the right-hand side of the Eq. (2), Using the Eq. (6), we get:

(16){\displaystyle {\int }_{Z_b}^H{\displaystyle \frac{\partial P}{\partial x}}\mathrm{d}z}={\displaystyle \frac{\partial }{\partial x}}{\displaystyle {\int }_{Z_b}^H{\rho }_{m}g(H-Z_b)\mathrm{d}z}={\rho }_{m}hg{\displaystyle \frac{\partial H}{\partial x}},

where H=z=Z_b+his the height of the free surface.

Similarly, we can integrate term by term the velocity derivatives on the right-hand side of Eq. (2) to establish the relationship between viscous stress gradients and their depth averages.

(17){\displaystyle \frac{\mu }{{\rho }_m}}{\displaystyle {\int }_{Z_b}^H{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}dz}={\displaystyle \frac{\mu {\beta }_{\mathrm{1}}}{{\rho }_m}}{\displaystyle \frac{{\partial }^{\mathrm{2}}hu}{\partial x^{\mathrm{2}}}},{\displaystyle \frac{\mu }{{\rho }_m}}{\displaystyle {\int }_{Z_b}^H{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}dz}={\displaystyle \frac{\mu {\beta }_{\mathrm{2}}}{{\rho }_m}}{\displaystyle \frac{{\partial }^{\mathrm{2}}hu}{\partial y^{\mathrm{2}}}},

where the coefficients {\beta }_{\mathrm{1}} and {\beta }_{\mathrm{2}} represent the ration of the vertical normal stress to the horizontal one, for a debris flow in which the material behaves more like a fluid,

{\beta }_{\mathrm{1}}={\beta }_{\mathrm{2}}=\beta =\mathrm{1},

(18){\displaystyle \frac{\mu }{{\rho }_{\mathrm{m}}}}{\displaystyle {\int }_{Z_b}^H{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}dz}={\displaystyle \frac{\mu }{{\rho }_{\mathrm{m}}}}{{\displaystyle \frac{\partial u}{\partial z}}|}_{z=H}-{\displaystyle \frac{\mu }{{\rho }_{\mathrm{m}}}}{{\displaystyle \frac{\partial u}{\partial z}}|}_{z=Z_b}.

Using the shear stress boundary conditions, at the free surface, {{\displaystyle \frac{\mu }{{\rho }_{\mathrm{m}}}}{\displaystyle \frac{\partial u}{\partial z}}|}_{z=H}={\displaystyle \frac{{\tau }_{sx}}{{\rho }_{\mathrm{m}}}}=\mathrm{0}, at the base, {{\displaystyle \frac{\mu }{{\rho }_{\mathrm{m}}}}{\displaystyle \frac{\partial u}{\partial z}}|}_{z=Z_b}={\displaystyle \frac{{\tau }_{bx}}{{\rho }_{\mathrm{m}}}}=\mathrm{0}, {\tau }_{bx} is the flow resistance. Naef et al. (2006) [34] gave the comparison of flow resistance relations for debris flows using a one dimensional finite element simulation model. For the two dimensional numerical simulation in this study, a combination of a viscous and Coulomb friction flow resistance is used:

(19){\tau }_{bx}=\mu {\rho }_{\mathrm{m}}\sqrt{gh}\cos {\theta }_x\tan \varphi ,

where {\theta }_x is the angle of inclination at the bed along the x-direction; and \varphi is the dynamic friction angle.

An analogous derivation must be performed for the right-hand side of second Navier-stokes Eq. (3).

The final form of the depth-averaged momentum equations are:

(20){\displaystyle \frac{\partial M}{\partial t}}+\alpha {\displaystyle \frac{\partial (uM)}{\partial x}}+\alpha {\displaystyle \frac{\partial (vM)}{\partial y}}=-{\displaystyle \frac{\partial H}{\partial x}}gh+{\displaystyle \frac{\mu \beta }{{\rho }_{\mathrm{m}}}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}M}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}M}{\partial y^{\mathrm{2}}}}\right)-\mu \sqrt{gh}\cos {\theta }_x\tan \varphi ,

(21){\displaystyle \frac{\partial N}{\partial t}}+\alpha {\displaystyle \frac{\partial (uN)}{\partial x}}+\alpha {\displaystyle \frac{\partial (vN)}{\partial y}}=-{\displaystyle \frac{\partial H}{\partial y}}gh+{\displaystyle \frac{\mu \beta }{{\rho }_{\mathrm{m}}}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial y^{\mathrm{2}}}}\right)-\mu \sqrt{gh}\cos {\theta }_y\tan \varphi .

Implementation of erosion and deposition

An estimation of the initial and final avalanche mass is of crucial importance in many debris flow problems, especially in mitigation studies where flow heights and volumes are required. Therefore, debris flow cover erosion and deposition must be modeled. Erosion and deposition is using a rate-controlled approach [35] which allows us to regulate both the mass uptake and the time delay required to accelerate the mass to the debris flow velocity. The effective rate E is parameterized by the dimensionless entrainment coefficient, K:

(22)E={\displaystyle \frac{{\rho }_{\mathrm{c}}}{{\rho }_{\mathrm{m}}}}K\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}.

Where{\rho }_{\mathrm{c}} is the soil density in channel. If K>\mathrm{0}, the erosion will be happened; On the contrary, the deposition will be happened.

The equation for the change of river bed elevation is:

(23){\displaystyle \frac{\partial Z_b}{\partial t}}=-E.

Numerical model coupled with GIS

Generation of DEM for simulation

To generate the topographical data required for the simulation using this model, a GIS-based digital is first converted to a Digital Elevation Model (DEM) (shown in Fig. 2(b)).The resolution of the DEM is 2.5 m × 2.5 m and is saved as raster in GIS. This raster-based DEM is then converted to a grid data file format readable by the T-model for the computer simulation. Each data point contains elevation and coordinates (Fig. 2(b)), and this data set forms our topographical model.

Source identification and upstream boundary setting

Next step is to identify the potential source zone. After field trip work in the area and after studying the available aerial photographs, potential source zone is identified in Fig. 2(b) in shaded area. According to field observation, the total volume of soil and source zone boulder that may be available in GIS. So the average thickness of slide discharge in source zone is estimated from volume divide area. This average thickness is used as the upstream boundary of simulation.

Numerical solution and results display

Numerical models are organized on a grid cell basis. Each cell has eight possible flow directions (left, right, up, down and the four diagonals) (see Fig. 2 (a)), but in fact that a cell overland flow is only routed along one flow direction which is the maximum down-slope direction (Fig. 2 (b)). The numerical solution of the above Eq. (1) to (5) is based on a leap-frog difference scheme that is accurate in space and time, the linear terms use forward difference scheme, and the nonlinear terms use central difference scheme. As denoted in Fig. 3, the flow depth, concentration or elevation of the debris flow mass in each grid is arranged at the midpoint, and the flux M and N are arranged at the boundary central point of the grid. The specific difference equations were derived by means of explicit staggered leap-frog scheme proposed by Yoshiaki (1988). Finally, these five equations have been coded by C++ language, the simple flow chart of the numerical simulation coupled with GIS is shown in Fig. 4. The yellow area is input area, the red area is computation area, and green area is output area. All the input and output data are processed in ArcGIS. The tool is embedded in a GIS to simplify the specification of input and to help the interpretation of numerical simulation results. Through all the calculated results displayed in GIS we can get the hazard maps for debris flow.

Application to Amamioshima debris flow

In this section some applications of the model are presented, the model was applied to a real debris flow occurred in Japan. On Oct. 20th, 2010, a debris flow induced by landslide happened in Yohutagawa (28^.24^’ N, 129^.31^’ E) in Amamioshima city. Although this debris flow only destroyed one building and one house, however there are kindergarten and a primary school in its path (Fig. 5). According to the detailed investigation of Kokusai Kogyo Group (KKG),the soil discharge and water depth are about 1.5 and 1.9 m. In this simulation, the total average thickness 3.4 m is used in the initial simulation. In addition, other material properties and rheological parameters well-documented by investigation of KKG used for simulation are listed in Table 1. The actual movement of soil and sand is shown in Fig. 6, we can see that the total volume in initial and erosion area is approximately 8697 m³, the volume in intercept area by dam is about 5621 m³, and the volume rush out of the dam in deposition area is about 3076 m³. The specific location of each area is also shown in Fig. 6.

A time-lapse simulation of the movements and affected regions of debris flow over the three-dimensional complexity terrain is illustrated in Fig. 7. The simulation results show that it took about 220 s to travel 1100 m along the channel, and an average flow velocity is about 5.0 m/s. The affected region can be dynamically displayed again at different times. Figure 8 shows the simulated hazard area and bed variation of this debris flow, we can see that the maximum erosion depth and maximum deposition depth are 1.51 and 1.96 m, respectively. Comparing with field investigation (Fig. 6), we can see that each area has good agreement in terms of both location and volume. It shows that the present model based on GIS which considered erosion and deposition is effective and rational. The max-flow depth and max-velocity distributions are shown in Figs. 9 and 10, respectively. We can see that the max-flow depth and max-velocity are 4.70 m and 23.08 m/s, respectively.

Conclusions

We have presented an approach to estimate the hazard of debris flow by incorporating the results of numerical simulation and GIS technology. A GIS environment provides a good platform for coupling a numerical model of a debris flow. As rater grid networks of digital elevation model in GIS can be used as the finite difference mesh, the governing equations are solved numerically using Leap-frog difference scheme. All the input and output data are processed in GIS. As a real case study, the model achieved reasonable results in comparison with field investigation. The conclusions and advantages of this method are listed as follows:

1) The numerical model describes debris flow covering both erosion and deposition processes. Comparing with the simulation results with field investigation, it shows that this model can well simulate not only movement but also erosion and deposition. Also that, some important parameters for evaluating debris flow disaster like max flow depth and max velocity can be obtained.

2) The pre-processing and post-processing of numerical simulation are easily realized by using GIS technology. The grid networks for computation can extracted from GIS, and the simulation results can be converted to GIS to form the hazard map for debris flow.

References

[1] Tang C, Rengers N, Yang Y H, Wang G F. Triggering conditions and depositional characteristics of a disastrous debris flow event in Zhouqu city, Gansu Province, northwestern China. Natural Hazards and Earth System Sciences, 2011, 11(11): 2903–2912

[2] Roy C S, Masahiro C. Landslides and debris flows strike Kyushu, Japan. Eos, Transactions, American Geophysical Union, 2004, 85: 145–156

[3] Armanini A, Fraccarollo L, Rosatti G. Two-dimensional simulation of debris flows in erodible channels. Computers & Geosciences, 2009, 35: 993–1006

[4] Hungr O. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Canadian Geotechnical Journal, 1995, 32: 610–623

[5] Hungr O, Evans S G. A dynamic model for landslides with changing mass. In: Marinos, Koukis, Tsiambaos and Stournaras, eds. Engineering Geology and the Environment. Athens, Greece, 1997, 719–724

[6] Rickenmann D, Koch T. Comparison of debris flow modelling approaches. In: Chen C L, ed. Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment. Proceedings of the first International DFHM Conference. ASCE, New York, USA, 1997, 576–585

[7] Coussot P. Some considerations on debris flow rheology. In: Ergenzinger P, Schmidt K H, eds. Dynamics and Geomorphology of Mountain Rivers, Dynamics and Geomorphology of Mountain Rivers. Lectures Notes in Earth Sciences 52, Berlin, 1994, 315– 326

[8] Whipple K X. Open channel flow of Bingham Fluids: Applications on debris flow research. Journal of Geology, 1997, 105: 243–262

[9] Brufau P, Garcia-Navarro P, Ghilardi P, Natale L, Savi F. 1D mathematical modelling of debris flow. Journal of Hydraulic Research, 2000, 38: 435–446

[10] Nakagawa H, Takahashi T. Estimation of a debris flow hydrograph and hazard area. In: Chen C L, eds. Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment. Proceedings of the first International DFHM Conference. San Francisco, CA, USA, 1997, 64–73

[11] Nakagawa H, Takahashi T, Satofuka Y. A debris-flow disaster on the fan of the Harihara River, Japan. In: Wieczorek G F, Naeser N D, eds. Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment. Proceedings of Second International DFHM Conference. Taipei, China, 2000, 193–201

[12] Shieh C L, Jan C D, Tsai Y F. A numerical simulation of debris flow and its applications. Natural Hazards, 1996, 13: 39–54

[13] Takahashi T. Debris Flow, IAHR Monograph Series. Rotterdam: Balkema, 1991

[14] Takahashi T, Nakagawa H, Harada T, Yamashiki Y. Routing debris flows with particle segregation. Journal of Hydraulic Engineering, 1992, 118(11): 1490–1507

[15] Zanre D D L, Needham D J. On simple waves and weak shock theory for the equations of alluvial river hydraulics. Philosophical Transactions of the Royal Society A, 1996, 354: 2993–3054

[16] Lai C. Modelling alluvial channel flow by multimode characteristic method. Journal of Engineering Mechanics, 1991, 117: 32–53

[17] Morris P H, Williams D J. Relative celerities of mobile bed flows with finite solid concentration. Journal of Hydraulic Engineering, 1996, 122: 311–315

[18] Iverson R M. The physics of debris flows. Reviews of Geophysics, 1997, 35: 245–296

[19] McDougall S, Hungr O. Dynamic modeling of entrainment in rapid landslides. Canadian Geotechnical Journal, 2005, 42: 1437–1448

[20] Chen H, Crosta G B, Lee C F. Erosional effects on runout of fast landslides debris flows and avalanches: A numerical investigation. Geotechnique, 2006, 56(5): 305–322

[21] Remaitre A. Mophologie et dynamique des laves torrentielles: Applications aux torrents des Terres Noires du basin de Barcelonnette (Alpes du Sud). Dissertation for the Doctoral Degree. Caen: Université de Caen, 2006

[22] Vandine D F, Bovis M. History and goals of Canadian debris flow research, a review. Natural Hazards, 2002, 26: 69–82

[23] Sovilla B, Margreth S, Bartelt P. On snow entrainment in avalanches dynamics calculations. Cold Regions Science and Technology, 2007, 47: 69–79

[24] Mangeney A, Roche O, Hungr O, Mangold N, Faccanoni G, Lucas A. Erosion and mobility in granular collapse over sloping beds. Journal of Geophysical Research, 2010, 115: 20–25

[25] Lin M L, Jeng S S. A preliminary study of risk assessment of debris flow using GIS. Journal of the Chinese Institute of Civil and Hydraulic Engineering, 1995, 7: 475–486 (in Chinese)

[26] Lin J Y, Yang M D, Lin P S. An evaluation on the applications of remote sensing and GIS in the estimation system on potential debris flow. Annual Conference of Chinese Geographic Information Society (CD edition), 1998

[27] Lin P S, Lin J Y, Hung J C, Yang M D. Risk assessment of potential debris flows using GIS. In: Proceedings of the Second International Conference on Debris-Flow. Balkema, 2000, 431–440

[28] Cheng K Y, Lin L K, Chang S Y. The field investigation and GIS application in a potential hazardous area of debris flow. In: Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, Proceedings of the First International Conference. ASCE, New York. 1997, 83–92

[29] Imran J, Harff P, Parker G. A number model of submarine debris-flow with a graphical user interface. Computers & Geosciences, 2001, 27: 717–729

[30] Iverson R M. The debris-flow rheology myth. In: Rickenmann D, Chen C L, eds. Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment. Rotterdam: Millpress, 2003, 303–314

[31] Vreugdenhil C B. Numerical Methods for Shallow-Water Flow. Norwell: Kluwer Acadamic, 1994

[32] Savage S B, Hutter K. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics, 1989, 199: 177–215

[33] Gray J M N T, Wieland M, Hutter K. Gravity-driven free surface flow of granular avalanches over complex basal topography. Proceedings of the Royal Society of London, 1999, 455: 1841– 1874

[34] Naef D, Rickenmann D, Rutschmann P, McArdell B W. Comparision of flow resistance relations for debris flows using a one-dimensional finite element simulation model. Natural Hazards and Earth System Sciences, 2006, 6: 155–165

[35] Christen M, Kowalski J, Bartelt P. RAMMS: Numerical simulation of snow avalanches in three-dimensional terrain. Cold Regions Science and Technology, 2010, 63: 1–14

[36] Yoshiaki I, Akihide T. Numerical simulations of flows and dispersive behaviours in southern part of lake Biwa. Disaster Prevention Research Institute Annuals, 1978, 31(B-2): 293–305 (in Japanese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg