GIS-based two-dimensional numerical simulation of rainfall-induced debris ﬂow

. This paper aims to present a useful numerical method to simulate the propagation and deposition of debris ﬂow across the three dimensional complex terrain. A depth-averaged two-dimensional numerical model is developed, in which the debris and water mixture is assumed to be continuous, incompressible, unsteady ﬂow. The model is based on the continuity equations and Navier-Stokes equations. Raster grid networks of digital elevation model in GIS provide a uniform grid system to describe complex topography. As the raster grid can be used as the ﬁnite difference mesh, the continuity and momentum equations are solved numerically using the ﬁnite difference method. The numerical model is applied to simulate the rainfall-induced debris ﬂow occurred in 20 July 2003, in Minamata City of southern Kyushu, Japan. The simulation reproduces the propagation and deposition and the results are in good agreement with the ﬁeld investigation. The synthesis of numerical method and GIS makes possible the solution of debris ﬂow over a realistic terrain, and can be used to estimate the ﬂow range, and to deﬁne potentially hazardous areas for homes and road section.


Introduction
Debris flows are rapidly flowing mixtures of water, clay, and granular materials and often triggered by torrential rains in mountainous areas. There are three main possible initiations of debris flows: mobilization from rainfall-induced landslide (Mainali and Rajaratnam, 1994;Anderson, 1995;Bathurst et al., 1997;Lan et al., 2004;Fiorillo and Wilson, 2004;Fleming et al., 1989;Wen and Aydin, 2005;Dai et al., 1999); erosion of steep debris beds in gullies due to intense rainfall Correspondence to: C. Wang (chunxiangwang@hotmail.com) (Takahashi et al., 1991(Takahashi et al., , 1992; and the collapse of a natural dam (Han and Wang, 1996;Valiani and Caleffi, 2003). All the debris flows have at least four characters: rainfall or dam-break is the triggering factor; a debris flow is a gravitydriven flow with free upper surface that move across threedimensional terrain; the nature of the flow itself, which is rapid, transient, and includes a steep front mainly constituted of boulders (Laigle and Coussot, 1997); and debris flows have very strong destructive power and bring about extensive property damage and loss of life to the communities in their path (Takahashi, 1991;Hunt, 1994;Huang and Garcia, 1997;Lien and Tsai, 2003).
As debris-flows are mixtures of flowing sediment and water showing complicated flow behavior intermediate between clear-water flows and mass movements of solid material, a number of mathematical rheological models were developed to simulate the flow behavior. Many researchers have developed rheological models for mudflows and debris flows. Theses models can be classified as: Newtonian models (Johnson, 1970;Trunk et al., 1986;Hunt, 1994;Hungr, 1995;Rickenmann, 1999), Bingham model (Johnson, 1970;O'Brien and Julien, 1988;Liu and Mei, 1989;Jan, 1997;Whipple, 1997;Fraccarollo and Papa, 2000;Pastor et al., 2004), Herschel-Bulkley model García, 1997, 1998;Imran et al., 2001;Remaître et al., 2005;Rickenmann et al., 2006), generalized viscoplastic model (Chen, 1988), dilatant fluid models (Bagnold, 1954;Takahashi, 1978Takahashi, , 1991Mainali and Rajaratnam, 1994), dispersive or turbulent stress models (Arai and Takahashi, 1986;O'Brien and Julien, 1988;Hunt, 1994), biviscous modified Bingham model (Dent and Lang, 1983), and frictional models (Iverson, 1997;Chen and Lee, 1999;Arattano and Franzi, 2003;Pastor et al., 2004;Rickenmann et al., 2006;Naef et al., 2006). Takahashi and Tsujimoto (1984) presented a twodimensional finite difference model for debris flows based on a dilatant-fluid model coupled with coulomb flow resistance, and modified the model to include turbulence (Taka-hashi et al., 1991(Taka-hashi et al., , 1992. O'Brien et al. (1993) developed a two-dimensional flooding routing model that is a valuable tool for delineating flood hazards and simulating floodwave attenuation, mudflows, debris flows (FLO-2D). Iverson and Denlinger (2001) developed a generalization of the depth-averaged, two-dimensional grain-fluid mixture model that describes finite masses of variably fluidized grain-fluid mixtures that move unsteady across three-dimensional terrain. Egashira et al. (2003) presented a method of numerical simulation for 2-D debris flow on an erodible bed using the constitutive equations for sediment-water mixture when the equation of erosion rate is incorporated in the continuity equation. McDougall and Hungr (2003) developed a depth-averaged model for the simulation of rapid landslide motion across complex 3-D terrain. Pudasaini and Hutter (2003) presented a two-dimensional depth-integrated theory for the gravity-driven free-surface flow of a granular avalanche over an arbitrarily but gently curved and twisted topography which is an important extension of the original Savage and Hutter (1989) theory. Bouchut and Westdickenberg (2004) developed a multidimensional shallow water model for arbitrary topography. Pastor et al. (2004) presented a depth-integrated Bingham model which is discretized using a Taylor-Galerkin finite element algorithm. Pudasaini and Hutter (2006) provided a survey and discussion about the motion of avalanche-like flows from initiation to run out. Rickenmann et al. (2006) compared three two-dimensional debris-flow simulation models with field events, and these models are based on a Voellmy fluid rheology reflecting turbulent-like and basal frictional stresses, a quadratic rheologic formulation including Bingham, collisional and turbulent stresses, and a Herschel-Bulkley rheology representing a viscoplastic fluid.
In recent years, Geographic information systems (GIS) with their excellent data format and spatial data processing ability have attracted great attention in natural disaster assessment. This is because the collection, manipulation, visualization and analysis of the environmental data on landslide and debris flow hazard can be accomplished much more efficiently and cost effectively Guzzetti et al., 1999). Key requirements in the assessment of debris flow risk consist of the prediction of the flow trajectory over the 3-D complex topography, the potential runout distance and the extent of the hazard area in order to define a safety zone. Numerical simulation models by incorporating GIS are important prediction and analysis tools.
The objective of this paper is to develop a two-dimensional depth-averaged numerical model to simulate the rainfallinduced debris flow and to integrate GIS with the numerical model to analyze the hazard of debris flow. As raster grid networks of digital elevation model in GIS can be used as the finite difference mesh, the continuity and momentum equations are solved numerically using the finite difference method.

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 (Imran et al., 2001). It is wellknow 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 rheologies appear to vary with time, position, and feedbacks that depend on evolving debris-flow dynamics (Iverson, 2003). 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.
The continuity equation is and the Navier-stokes equations are as follows: in which u, v, w=velocity components in the x-, y-, and z-directions; ρ d =equivalent density of the debris and water mixture, and ρ d =ρ s v s +ρ w v w , ρ s and ρ w are the densities of solid grains and water, v s and v w are the volumetric concentrations of solids particles and water in the mixture; p=pressure; µ=dynamic viscosity; g=gravitational acceleration; and t=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 (cf. Vreugdenhil, 1994). As described by Savage and Hutter (1989), Iverson (1997), and Gray et al. (1999) and is deemed generally much smaller than 1. The timescale for order-unity changes is l g 1 2 , because the potential for free fall drives debris flow motion (Iverson, 1997). These time and length scales in turn lead to differing velocity scales for the flow direction, gl 1 2 , and z direction, ε gl 1 2 , which imply that u, 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: Using the notation and coordinate systems given in Fig. 1, Integration of (5) with respect to z from the elevation of the base z=η b to the upper surface of the flow at z=η yields, 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=η b to the free surface of the flow at z=η, using the Leibnitz rule to interchange the order of integration and differentiation. Herein, we define the depth-averaged velocities are as follows: in (7), and equations hereinafter, overbars denote depthaveraged quantities defined by integrals similar to that (7). The continuity Eq. (1) is integrated through the flow depth, using the kinematic boundary conditions at the base and at the free surface, we obtain, as h=h(x, y, t)=η−η b is the flow depth, we get, Inserting (7), (9), (10), (11) into (8) in which M=uh and N=vh are the x-and y-components of the flow flux respectively; Integrating the left-hand side of the Eq. (2), Two ad-hoc coefficients α 1 and α 2 have been introduced. As was noted by Savage and Hutter (1989), values of α 1 and α 2 in (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, α 1 =α 2 =1. For parabolic velocity profile and debris flows with no basal sliding, α 1 =α 2 =6/5, and for a stone-type debris flow on a rough inclined plane, α 1 = α 2 =1.25 (Takahashi et al., 1992). We may write for (13) using α=α 1 =α 2 , Furthermore, we integrate the right-hand side of the Eq.
(2), Using the Eq. (6), we get Where H =z=η b +h is 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.
Where the coefficients β 1 and β 2 represent the ration of the vertical normal stress to the horizontal one, for a rainfallinduced debris flow in which the material behaves more like a fluid, β 1 =β 2 =β=1. = τ bx ρ d , τ bx is the flow resistance. Naef et al. (2006) gave the comparison of flow resistance relations for debris flows using a onedimensional 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: where θ x is the angle of inclination at the bed along the xdirection; and ϕ is the dynamic friction angle, and tan ϕ is the dynamic friction coefficient. An analogous derivation must be performed for the righthand side of the second Navier-stokes Eq. (3).
The final form of the depth-averaged momentum equations are 3 Incorporation of numerical simulation with GIS

Digital elevation model in GIS
Recently, Geographic Information Systems (GIS), with their excellent data structures and spatial data processing capacity, have attracted great attention in landslides and debris flows disasters assessment. GIS is a computer-assisted system for the acquisition, storage, analysis and display of geographic data. And GIS provides strong functions in spatially distributed data processing and analysis. Digital Elevation Models (DEM) within GIS automatically extract topographic variables, such as basin geometry, stream networks, slope, aspect, flow direction, etc. from raster elevation data. A DEM is a numerical representation of landscape topography. Three schemes for structuring elevation data for DEMs are: triangulated irregular networks (TIN), grid networks, and vector or contour-based networks (Moore and Grayson, 1991). The most widely used data structures are square grid networks with rows and columns where each cell contains a value representing information, such as elevation. The gridbased discretization of the studied area is immensely useful for numerical solution of the partial differential equations governing the propagation of debris flows. Finite-difference method on rectangular grids are widely used in numerical models of environmental flows when using this method, the studied region is discretized into rectangular grids. Therefore, in this paper, we used grid nerworks in GIS as the rectangular grids of finite difference methods.

Numerical scheme
Numerical models are organized on a grid cell basis. In a raster-based DEM analysis, each cell has eight possible flow directions (left, right, up, down, plus the four diagonals), as show in Fig. 2. The flow direction of a cell is expressed in degrees: left=0 • , up=90 • , right=180 • , down=270 • ; and the diagonals: 45 • , 135 • , 225 • , 315 • . Within a cell overland flow is routed along one flow direction. The flow direction is the maximum downslope direction which is determined from the raster-based DEM (Fig. 3). The numerical solution is achieved using a finite difference formulation based on the DEM grid. The governing equations are approximated using leapfrog time-differencing. A staggered grid approach is followed to evaluate the spatial gradients of Eqs. (12), (20), and (21). A forward difference scheme is applied to discretize the linear terms, and a central difference scheme is applied to discretize the nonlinear terms. The thickness of the debris mass in each cell is computed at the midpoint of the cell and denoted in Fig. 3. The method of adjusting the time step and mesh size is used to prevent the appearance of numerical instability due to the use of too large a time interval and mesh size. The finite difference form of the continuity equation is: And the momentum equation of x-component A finite difference expression analogous to Eq. (23) represents the momentum equation of the y-component in Eq. (21). The model outlined above has been coded into a numerical model of general utility in ArcGIS(a GIS software developed by ESRI). The code is written in Visual Basic language and ArcObjects as a tool of ArcGIS application. ArcObjects is the development environment of the desktop ArcGIS applications (ArcMap, ArcCatalog and ArcScene). It is used to customize and extend ArcGIS using Visual Basic. 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.

Simulation of real debris flow event
To verify the model and illustrate its validation, the model is applied to a real debris flow occurred in Japan. On 19-20 July 2003, a short duration of high intensity rainfall event impacted on Minamata City in the Kumamoto prefecture, Japan. This rainstorm triggered many landslides and debris flows (Nakazawa et al., 2003;Iwao, 2003;Taniguchi, 2003;Mizuno et al., 2003;Hashimoto et al., 2003). The landslide and resultant debris flow at Hogawachi was the largest and most damaging of the 20 July disasters in Minamata area (Fig. 4). The debris flow occurred 4.3 h after the beginning of the rainstorm, at 4:20 a.m. on 20 July 2003, during the period of highest intensity (Fig. 5). The landslide that triggered the debris flow initiated in highly weathered andesite (An-7: mainly lava with subordinate tuff breccia) underlain by an almost impermeable layer (An-5: tuff breccia, brecciated lava, and lava) (Fig. 6). The maximum depth of the landslide is about 15 m. The debris flow began after the landslide entered the stream valley and traveled about 1.5km along the channel. The gradient of the slope that failed was 19 • near the top and 36 • in the lower section. The gradient of the channel is about 17 • ∼9 • , and the mean gradient of the fan is 5 • (Hashimoto et al., 2003). The volume of sediment discharge plunging into the village of Hougawachi-Atsumari district, Minamata City, was estimated to be about 68 000 m 3 , and the velocity of debris flow was estimated from about 2.9 m/s to 23.5 m/s (Mizuno et al., 2003;Taniguchi, 2003). This disaster killed 15 people and more than 14 houses were either damaged or destroyed. Based on a topographic map 1/2500 in scale, a vector contour line file is generated, with vertical spacing of 2 m. This file is converted to TIN (Triangulated Irregular Network), and subsequently DEM (Digital Elevation Model) within ArcGIS. The road and river are drawn as polylines, and the homes as polygons (Fig. 7). The DEM resolution is 2 m×2 m. The elevation of terrace and river are expressed in this DEM.

Debris avalanche initiation zone Region of fatalities
How to determine the input parameters is important because these parameters will affect the simulation results. The range of parameters for simulation was constrained using field observations because no direct measurements of the parameters are available . In this simulation, the square grid mesh in GIS is x= y=2 m and the depth-averaged velocities are considered as blunt, therefore, we have set α=1, β=1. As a rough approximation by field investigation, the sediment volume is about 68 000 m 3 , and the water volume is about 17 000 m 3 (Mizuno et al., 2003). The effective viscosity and the dynamic friction coefficient are key parameters that affect the flow behavior. Lien and Tsai (2003) gives the range of the dynamic friction coefficient from 0.32 to 0.75. Based on field investigation, the inundated area is 0.15 km 2 and the average thickness of deposits is assumed to be 3.5 m. Using these data and through iterative calculations, the best-fit pair of tan ϕ=0.4 and µ=0.11 are selected (Table 1). A time-lapse simulation of the dynamic progression and deposition of the debris flow over the 3-dimensional terrain is illustrated in Fig. 8. The simulation results show that it took about 200 seconds to travel 1500 m along the channel, and an average flow velocity is about 7.5 m/s. The affected region can be dynamically displayed again at different time. In order to know the flow depth and the deposition thickness in the simulation, a cross section BB and CC was given (Fig. 7). At the 60 s after the beginning of the event, the debris-flow arrived BB-cross section. Figure 9 shows the dynamic flow depth of BB-cross section at 60 s, 90 s, 120 s, 150 s, and 180 s. The deposition thickness and the dangerous homes in CC cross section are shown in Fig. 10. Figure 11a shows the simulation result of the area affected by debris flow. Comparing with the affected area tracked by aerial photograph and by field survey (Fig. 11b), the simulation result is in close agreement with filed observation. The homes affected by the simulation result are 15, which is the same with field observation also. This means that the numerical approach can be properly used to simulate the real debris-flow triggered by landslide and rainstorm in the study area.

Discussion
Debris flows are complex phenomena, due to spatial and temporal variability in material properties, and state of trace during their evolution. Their sudden occurrence and transient character make difficult to directly measure their rheological parameters and friction force. When using numerical simulation models for hazard assessment of debris flows, the selection of appropriate friction or flow resistance parameters is of great importance (Hungr, 1995;Ayotte and Hungr, 2000;Brufau et al., 2001;Arattano and Franzi, 2003;Revellino et al., 2004;Naef et al., 2006;Rickenmann et al., 2006). For the detailed reviews of different equations describing the flow resistance relations, represented by S f , see Hungr (1995), Chen (1988), Brufau et al. (2001) and Naef et al. (2006). The flow resistance term depends on the rheology of the material and is a function of several different known parameters of the flow (Hungr, 1995). On the other hand, a comprehensive understanding is still lacking, and a classification of debris flow into specific categories is often difficult, this is due to both to the wide range of water-sediment flows referred to as debris flows and to the flow behavior variability encountered (Sosio et al., 2007). In this study, the flow resistance is considered as the function of the flow depth, the effective viscosity, the gravity, the path slope angle, and the dynamic friction coefficient. The effective viscosity and the dynamic friction coefficient are the key parameters that affect the flow behavior.
The elevation values of Digital Elevation Model (DEM) are used in algorithms to calculate surface derivatives such as slope, aspect, flow direction and will affected the flow depth  in this model. The grid cell size of DEM within GIS will also affect the simulation results. Accurate representation of the channel and fan topography in the grid is especially important to achieve a good replication of the observed deposition pattern . To explore how different DEM resolutions influence the deposition thickness in the simulation, three other DEMs with the grid sizes 5 m×5 m, 8 m×8 m and 10 m×10 m were generated in ArcGIS. Figure 12 shows the depositions and inundated areas of the four different DEM resolutions (grid sizes of 2 m×2 m, 5 m×5 m, 8 m×8 m, 10 m×10 m). Using the same rheologic parameters, the ratio of the mean deposition thickness of grid sizes in 10 m×10 m to 2 m×2 m is 0.71, the ratio of the mean deposition thickness of grid sizes in 8 m×8 m to 2 m×2 m is 0.78, and the ratio of the mean deposition thickness of grid sizes in 5 m×5 m to 2 m×2 m is 0.92; the computation times of grid sizes 2 m×2 m, 5 m×5 m, 8 m×8 m, 10 m×10 m are 6.5 h, 5.2 h, 4.3 h, and 4.1 h, respectively. Comparing with the deposition and the inundated area of the field investigation, the results of the grid sizes of 2 m×2 m is better than the others. A coarse resolution DEM appears to be of poor quality results, but it can be used for rough estimation of the flow range and the deposition.

Conclusions
In the present study, we presented a two-dimensional depthaveraged numerical model to simulation the propagation and the inundated area of debris flow, and numerical simulation methods in combination with GIS-technology were applied. A GIS environment provides a good platform for coupling a numerical model of a debris flow. As raster grid networks of digital elevation model in GIS can be used as the finite difference mesh, the continuity and momentum equations are solved numerically using finite difference method. All the input and output data are processed in GIS. As a case study, we applied the model to a rainfall-induced real debris flow occurred in 20 July 2003, in Minamata City of southern Kyushu, Japan. The model achieved reasonable results in comparison with a field investigation. This simulation redisplays the propagation and deposition of the debris flow across the complex topography. How to prevent or mitigate disaster caused by landslides and debris flows is an urgent problem. Therefore, prediction of the characters of the landslide and debris flow, such as the possible inundated extent of the moving debris mass, the propagation progress across the real three-dimensional terrain, the area of deposition, and finding the dangerous homes   and other infrastructure is of great importance in hazard and risk assessment. The advantages of numerical method combining with GIS-technology are that the preprocessing routine in which the computation data are prepared, the postcomputation visualization, the results analysis, and providing an effective tool for risk analysis and hazard mapping.