A GIS-based numerical model for simulating the kinematics of mud and debris flows over complex terrain

This article presents MassMov2D, a twodimensional model of mud and debris flow dynamics over complex topography, based on a numerical integration of the depth-averaged motion equations using a shallow water approximation. The core part of the model was implemented using the GIS scripting language PCRaster. This environment provides visualization of the results through map animations and time series, and a user-friendly interface. The constitutive equations and the numerical solution adopted in MassMov2D are presented in this article. The model was applied to two field case studies of mud flows on torrential alluvial fans, one in the Austrian Tyrol (Wartschenbach torrent) and the other in the French Alps (Faucon torrent). Existing data on the debris flow volume, input discharge and deposits were used to back-analyze those events and estimate the values of the leading parameters. The results were compared with modeling codes used by other authors for the same case studies. The results obtained with MassMov2D matched well with the observed debris flow deposits, and are in agreement with those obtained using alternative codes.


Introduction
Mud and debris flows can be defined as gravity-driven flows of a highly saturated mixture of debris in a steep channel, which exhibit properties of viscous and turbulent flows. According to Hungr et al. (2001), the difference between debris and mud flows is related to the significantly greater water content of the latter, which causes a plastic behavior (plasticity index higher than 5%). Mud and debris flows can origi-Correspondence to: S. Beguería (sbegueria@eead.csic.es) nate either at a single source, typically from the fluidization of a failed mass (landslide), and/or by the re-entrainment of sediment accumulated in a torrential catchment. Fluidization of the mass usually requires a large amount of water, but rather dry clastic avalanches can also show flow-like behavior due to grain collisions, and flow on surfaces less inclined than the angle of repose of the material. When debris flows are confined on a torrent, they can propagate over very large distances before final spreading over an alluvial fan. Debris flows are a common and important factor in erosion and sediment transfer in mountainous areas, and constitute an important risk to the population. The threats to human life and property from mud and debris flows are greater than those of some other landslide types, due to their higher velocity (sometimes in the order of tens of m s −1 ) and capacity to propagate even on very gentle slopes (below 5 • ). Issues identified in recent reviews of the physics of debris flows include the phenomenon of dilatancy; the roles of grain collisions, internal friction and cohesion; the mechanisms of fluidization; and the phenomena of particle segregation and deposition (e.g. Iverson and Denlinger, 1987;Takahashi, 1991;Hutter et al., 1996;Iverson, 1997).
In addition to theoretical interest, the development of robust numerical models of debris flow has important applications in hazard evaluation and the design of mitigation measures, as the models enable quantitative estimates of debris flow characteristics including flow depth and velocity, the force exerted against obstacles, and evaluation of preventive structures. Different simulation scenarios can be easily constructed by varying the governing parameters and boundary conditions (e.g. the input discharge), which allows assessment of the likelihood that a given area will be affected under differing circumstances. This technique is especially useful if information is available about variability in the parameters controlling the flow. A number of runout models 1898 S. Beguería et al.: Kinematics of mud and debris flows over complex terrain have been developed for both one-dimensional (Hungr, 1995;Fraccarollo and Papa, 2000;Arattano and Franzi, 2003;Naef et al., 2006) and two-dimensional cases (Savage and Hutter, 1989;O'Brien et al., 1993;Shieh et al., 1996;Laigle and Coussot, 1997;Chen and Lee, 2000;Iverson and Denlinger, 2001;Pudasaini et al., 2005;Hungr and McDougall, 2008). Two-dimensional (2-D) models are becoming common tools for the analysis and prediction of mass movements including rock avalanches, earth and debris flows, and mud flows. On 11-12 December 2007, a benchmarking exercise on landslide runout modeling was organized in Hong Kong as part of the 2007 International Forum on Landslide Disaster Management, comparing various 2-D codes (Ho and Li, 2007).
The main objective of this article is to describe the development and implementation of a numerical code, Mass-Mov2D (a copy of which can be downloaded from http: //hdl.handle.net/10261/11804), for simulating mud and debris flows over complex topography, implemented in the PCRaster GIS environment (Wesseling et al., 1996;Karssenberg et al., 2001). The model is based on a 2-D finite difference solution of a depth-averaged form of the fluid dynamics equations. The flow is thus treated as a one-phase material, whose behavior is controlled by rheology (i.e. by a functional relationship between strain and stress). Several flow laws can be implemented within a common numerical scheme in the model, which can accept a detailed description of a complex topography through a digital elevation model (DEM), thus facilitating simulation of field case studies.
In contrast to other codes, the purpose of the model described here is to serve as a generic framework within which various modeling concepts can be tested. The model is flexible enough to work with various configurations of the initial and boundary conditions, a characteristic that will allow researchers to adapt the simulation to a variety of realistic situations, ranging from alluvial fan dispersion to channel overflow problems. The implementation of several rheological models allows comparison of various flow types. The model was implemented in a geographical information system (GIS) package, thus benefiting from a set of tools that facilitate the preparation of input data and evaluation of the results. However, implementation in a GIS platform is beneficial beyond being the mere pre-and post-preprocessing tools available, as the core equations of the model (mass balance, equation of motion, rheology) are accessible to the user in an easy-to-learn scripting language. This enables easy testing of new modeling concepts through modification of the original code.
We first describe the governing equations of the model, the various flow laws, the numerical implementation and the specification of the initial and boundary conditions, and provide examples of the model operation. We then present a case study showing the results of application of the model to two debris flow events on torrential alluvial fans, one in the Austrian Tyrol (Wartschenbach torrent) and the other in the French Alps (Faucon torrent).

Governing equations
The approach was based on the classical Savage-Hutter theory (Savage and Hutter, 1989), which assumes a one-phase homogeneous material with rheological properties. The flow was modeled as a 2-D continuum medium using a depthintegrated approximation of the flow dynamics equations; this has become a classical approach to debris flow modeling (Savage and Hutter, 1989;Laigle and Coussot, 1997;Fraccarollo and Papa, 2000;Denlinger and Iverson, 2001;Mc-Dougall and Hungr, 2004;Mangeney-Castelnau et al., 2005;Pirulli et al., 2007). Depth integration is based on the shallow water assumption (Saint Venant equations), which applies where the horizontal length scale (length of the flowing mass) is much greater than the vertical length scale (thickness of the flowing mass). In these conditions the vertical velocity of the fluid is small, so that the vertical pressure gradient is nearly hydrostatic. Integration over the flow depth allows a 2-D description of the flow, avoiding the much more complex 3-D formulations. There is experimental evidence of agreement with the experimental data (Savage and Hutter, 1989;Greve and Hutter, 1993).
Contrary to depth-integrated models that use a local reference frame linked to the topography, the equations governing MassMov2D are referenced in a 2-D Euclidean space with Cartesian coordinates x, y (Fig. 1). Very similar formulations have recently been proposed . The mass (1) and momentum (2) balance equations were developed as: ∂u ∂t + c x u ∂u ∂x + c y v ∂u ∂y = −c x g S x + k ∂(c x h) ∂x + S f q x ∂v ∂t + c y v ∂v ∂x + c x u ∂v ∂y = −c y g S y + k ∂(c y h) ∂y + S f q y , (2) respectively, where h is the flow thickness in the direction normal to the bed (L); (u,v) are the x and y components of the velocity vector along the bed (L T −1 ); the coefficients c x =cosα x and c y =cosα y are the direction cosine of the bed (geometry factors to correct from local to global reference systems); and α x and α y are the values of the angle between the bed and the horizontal plane in the x and y directions, which take a negative value if the down slope direction is towards positive x and y, respectively. The momentum equation (2) is expressed in terms of acceleration (L T −2 ), where g is the acceleration due to gravity. The second and third terms on the left side of Eq.
(2) represent the convective acceleration, i.e. the time rate of change due to change in position in the spatial field. The right side of the equation represents the local or time acceleration, expressing the time rate of change at a fixed position. The first term between the brackets represents the acceleration due to gravity, and S x =tanα x and S y =tanα y are the bed slope gradient in the x and y directions (L L −1 ), respectively. The spatial derivative in the second term is the pressure acceleration, i.e. the time rate of change due to pressure differences within the flow. S f is the flow resistance gradient (L L −1 ), which accounts for momentum dissipation within the flow due to frictional stress with the bed (see next section). The term k in Eq.
(2) is the earth pressure coefficient, i.e. the ratio between the tangential and normal stresses. It has a value of 1 for a perfect fluid, but can vary greatly for plastic materials (Savage and Hutter, 1989;Hungr, 1995), and ranges between two extreme values corresponding to the active and passive states in the Rankine theory, i.e. k a ≤1≤k p . These values depend on the internal friction angle of the mixture, δ: The active state k a applies to regions of the flow that undergo expansion, The passive state k p corresponds to regions experiencing compression, Finally, q x and q y are coefficients: where the minus sign before u and v ensures that S f opposes the direction of the velocity.

Flow laws
The flow term S f in Eq.
(2) represents the bed shear stress of the flow, which is responsible for energy dissipation. This variable describes the rheological properties of the flow, which control flow behavior. The model defined by Eqs. (1) and (2) is therefore valid for different types of flows, depending on the formulation of S f . One-phase, depth-integrated models characterize flow behavior by a rheological relationship, and typically assume homogeneous and constant flow properties. These simplified models are preferred for applied field studies because of their parsimony, which makes them more suitable for calibration and back-analysis. A number of rheological models for the flow term have been reviewed by Naef et al. (2006). Mud and debris flows have often been modeled as viscoplastic materials with a laminar flow regime, i.e. as Bingham fluids with constant yield strength and viscosity (Coussot, 1994). A linear stress-strain rate relationship is assumed once the yield strength is exceeded. However, a convex relationship defining a shear thinning behavior (i.e. decreasing viscosity with applied shear) is in better agreement with most experiments, and led to definition of the Herschel-Bulkley model. Both models can be described by the following relationship: where τ (z) is the resisting shear stress at a given depth z (M L −1 T −2 ), τ C is a constant yield strength due to cohesion, ∂v/∂z is the shear rate, and µ (M L −1 T −1 )-beta is a consistency index acting as a viscosity parameter. The viscosity of the flow is closely related to the percent concentration of solids (O'Brien and Julien, 1988;Parsons et al., 2002), and can range from 10 −2 to 10 3 . The exponent β is an empirical parameter of value 1 in the Bingham model, and in the Herschel-Bulkley model has to be determined with the aid of experimental data. A number of debris flow models are based on rheological formulations for a Bingham viscous fluid, or the Herschel-Bulkley model with different values of β (O' Brien et al., 1993;Coussot, 1994;Han and Wang, 1996;Laigle, 1997;Laigle and Coussot, 1997;Fraccarollo and Papa, 2000;Zanuttigh and Lamberti, 2004;. The Bingham/Herschel-Bulkley model is valid for materials where the fine fraction is large enough to lubricate contacts between grains, i.e. approximately 10% clays size fraction (Ancey, 2001;Rickenmann et al., 2006). A more general alternative model is obtained if the constant yield strength in Eq. (5) is replaced by a combined cohesive-frictional component, leading to the Coulomb-viscous model (Johnson and Rodine, 1984): which incorporates a frictional component depending on the is the internal pore fluid pressure, and ϕ ( • ) is the basal friction angle of the flow. The internal pore fluid pressure is a transient property that is coupled to the normal stress and can dissipate during motion, making it extremely difficult to model. Thus, in most practical applications the pore pressure ratio (u/σ ) is assumed to be constant, and hence the frictional dissipation can be lumped into a single parameter (tanϕ ), where ϕ is an apparent friction angle. This approach was implemented in MassMov2D, based on the hypothesis that the pore pressure is conserved and remains constant during the simulation. Mixed frictional-viscous models similar to that in Eq. (6) are suited to describing the shear behavior of granular solid materials (Naef et al., 2006). The Coulomb frictional-viscous model has also been used to model debris flows (O'Brien et al., 1993;Coussot, 1994;Han and Wang, 1996;Laigle, 1997;Laigle and Coussot, 1997;Fraccarollo and Papa, 2000;Zanuttigh and Lamberti, 2004;Van Asch et al., 2007) and has also been applied to very large landslides in which the formation of a melted rock layer adds extra lubrication to the flow (De Blasio and Elverhøi, 2008). In applying Eqs. (5) and (6) there is a need to integrate the stress over the flow depth, to enable computation of the total stress, S f . In depth-integrated models this is solved by assuming that the stress distribution in the direction normal to the flow is triangular. Following Coussot (1997), we have used a simplified form of the third-order Bingham expression for the shear stress to derive S f : Equation (7) is a good approximation to the original formulation for values of the stress ratio (τ C /τ ) smaller than 0.5 (Naef et al., 2006). For the Coulomb-viscous rheology we have used: Hence, for particular simulations the rheology of the flow is defined by a parameter vector (ρ,τ C ,µ) or ρ,ϕ ,τ C ,µ for the Bingham and the Coulomb-viscous laws, respectively.

Numerical implementation
The numerical implementation of hyperbolic partial differential equations, such as Eqs. (1) and (2) is a challenging task, especially for systems exhibiting an advancing front (see the review by McDougall, 2006). We adopted a reasonably simple compromise solution that achieved a desired level of stability, accuracy and controlled diffusivity.
It is convenient to write Eqs. (1) and (2) in more compact vector notation, in order to describe the numerical solution used in MassMov2D: The model was implemented in an explicit finite difference (Eulerian) mesh, i.e. the flow was described by variation in the conservative variables at points of fixed coordinates (i, j ) as a function of time (n). The mesh is defined as a regular grid with size s= x= y, in accordance with the grid format common to most GIS platforms. Equation (9) is solved numerically using a central difference forward scheme: where t is the time step duration (s), and the pressure gradient term in q is computed by central differences. A common problem with such simple methods is the introduction of dispersive effects that lead to unphysical oscillations, especially in the presence of large gradients. A solution to this problem is the addition of a certain amount of numerical regularization, as introduced in Eq. (10) by the function W w n i,j : This function performs a weighted spatial averaging over w, and the amount of numerical regularization is controlled by the value of the Courant-Levy-Friedrichs condition (CFL; see Eq. 12, below) at each point, so it is applied with preference to the areas of the flow that are experiencing sudden changes and have values of CFL typically in the range 0.5-1. This solution proved to provide only the required amount of regularization needed to avoid numerical instability, while adding just a minimum amount of artificial dispersion. Reflecting boundary conditions are applied to Eqs. (10) and (11), so if any of w n i−1,j , w n i+1,j , w n i,j −1 or w n i,j +1 lie outside the boundaries of the flow, then they take the value of w n i+1,j , w n i−1,j , w n i,j +1 or w n i,j −1 , respectively. Another problem with first-order time solutions is the over-and underestimation of the flow resistance term, which typically happens in accelerating and decelerating flows. To deal with this problem a two-step solution was adopted. Hence, the source term s n+1/2 i,j in Eq. (10) was evaluated at interleaved time steps to reduce over-and undershoots. The velocity components of w were estimated at times n+1/2 by applying Eq. (10) to w n with t = t/2. The resulting solution is thus similar to the Lax-Wendroff second-order scheme (Press et al., 2002). To maintain the stability of the solution throughout the simulation, the Courant-Friedrichs-Levy condition CFL<1 must be satisfied. CFL was calculated at each point of the mesh as: The time step duration t was adapted at each iteration to the highest velocities found in the spatial domain, to ensure that the Courant-Friedrichs-Levy condition was met over the whole computation domain. Despite this variable internal time step, reporting is made at a regular time interval set by the user. An underlying assumption of the finite difference solution is that of continuity and smoothness of the spatial domain. Although this generally holds true, in some cases the basal topography presents horizontal discontinuities or singularities such as channels (artificial or natural), which are very relevant to the correct simulation of the flow. Such structures can not be adequately represented in a finite differences mesh, as the central difference scheme requires a continuous surface. Specification of topographical discontinuities was incorporated into the model such that they are treated as no-flow boundaries unless the flow thickness becomes larger than the difference in height at both sides of the discontinuity. Whenever this occurs, flow is automatically allowed between the two sides of the topographical discontinuity by applying Eq. (10) to the fraction of h that exceeds the height of the obstacle.
The model was implemented in the PCRaster environmental modeling language (Wesseling et al., 1996;Karssenberg et al., 2001). PCRaster is a scripting language; hence the computation time involved is greater than it would be using a compiled language such as C or Fortran. For this reason the basic vector operators (gradient, divergence) and some functions that are used repeatedly in the code, including the W () function in Eq. (10) and other parts of the numerical scheme, were programmed in C and compiled as new PCRaster functions. However, the core functions of the model were implemented in the PCRaster scripting language, promoting easy understanding and modification of the original code. This makes MassMov2D a convenient tool for researchers interested in testing new modeling concepts.

Initial and boundary conditions
The model requires three input files comprising the topography of the bed, and the locations of the inlet and outlet cells. The topography of the bed is a basic input to the model and is provided in the form of a DEM. The DEM defines not only the basal boundary of the flow, but also the spatial computation domain and the mesh size. No flow is allowed outside the spatial limits of the DEM (closed boundary) except at two specific areas: i) the inlet, a set of contiguous cells where input flow is allowed at a rate specified by the user (see below); and ii) the outlet, a set of cells where flow is allowed to exit the computation domain (open boundary). Additionally, a map with the location of topographical discontinuities, such as channels or dikes, can also be set. This map identifies the domain cells that lie inside the channel, as well as the height of the discontinuity in the direction normal to the basal surface (m). These files can easily be created as raster maps using standard PCRaster commands.
In addition to these maps, the model requires specification of the input discharge at the cells defined in the inlet map during the simulation time, w t inlt (boundary condition). This can be set either as a constant inflow rate for a specified time interval, or as a transient flow by creating a text file containing three columns that indicate the time variation of the state variables (h, u, v) at the inlet.
An example of the configuration of a simulation is provided in Fig. 2. A rectilinear channel 150 m long and 3 m deep with a constant gradient α=10.5 • was defined. The channel borders were treated as spatial discontinuities, meaning that only in the case of the flow depth in the channel exceeding 3 m would there be overflow and connection with the flow on the banks. The lateral limits of the spatial domain were defined as impermeable boundaries, and inlet and outlet boundaries were defined at the upper (upstream) and the lower (downstream) limits. An input discharge file was defined with a triangular shape, rising from h=0.65 m at 0 s to h=0.85 m at 25 s, and then falling to h=0 m at 30 s. Bingham rheology was used, with parameters: ρ=1800 kg m −3 , ϕ =5.5 • and µ=50 Pa s. Deceleration of the flow along the channel generated an increase of the input wave, and resulted in overflow a few seconds after the peak discharge. After that moment the flow over the lateral banks remained disconnected from the flow in the channel.

Case studies
There are many difficulties involved in the simulation of field events, particularly related to the availability of data of sufficient quality, and the simplifying assumptions made in most models. Despite these difficulties, simple simulation models assuming constant flow properties have been successful in predicting the major characteristics of mud and debris slurries (e.g. Shieh et al., 1996;Chen and Lee, 2000;Laigle, 2000;Laigle et al, 2003;Bertolo and Wieczorek, 2005;Rickenmann et al., 2006;Wang et al., 2008;Medina et al., 2008). In most cases the extent and thickness of the deposited material were the only known properties, and rarely was there information about the flow velocity (e.g. Revellino et al., 2004). The total volume of the event was generally estimated from the final deposit, which only includes the solid part of the debris flow mixture, requiring that the liquid part be estimated.
The simulation of observed events is not only a procedure for validating models, but can also deliver results relevant to case studies. For example, there are no ways to measure directly the rheological properties of the flow, so these need to be calculated by back-analysis using simulation tools.
We used MassMov2D to simulate field events of mud and debris flows in torrential alluvial fans, one in Austria and the other in France. In both cases the simulations were performed using the Bingham and Coulomb-viscous rheologies, and the most appropriate values of the rheological parameters for each case were established using a heuristic process, comparing the outcomes from different parameter sets with the field deposition patterns. These two events were selected because both had previously been analyzed using other, wellestablished codes, thus enabling comparison of these results with those from MassMov2D.

The event of 16 August 1997 in the Wartschenbach torrent
The Wartschenbach torrent is located near Lienz, eastern Tyrol, Austria. It is a small mountain torrent with a catchment area of 2.5 km 2 and an altitude ranging from 600 to 2500 m a.s.l. At the bottom the torrent forms an alluvial fan with a mean gradient of 16% (around 9 • ). A debris flow event occurred in the Wartschenbach torrent on 16 August 1997 after an intense storm. A large amount of sediment was trapped in an upstream debris retention dam. Between 25 000 and 30 000 m 3 of debris and water reached the fan apex. This caused overflow of the little channel, and spread of the debris flow over the fan affected 15 buildings (Hübl and Steinwendtner, 2000). The same debris flow event was modeled by Rickenmann et al. (2006), who compared the results of three different codes for debris flow simulation: DFEM-2D (Naef et al., 2006), FLO-2D (O'Brien et al., 1993) and HB (Laigle, 1997;Laigle and Coussot, 1997). A pre-event DEM with a grid resolution of 2 m was used in the simulation. The inlet area was located at the fan apex, in the upper part of the DEM. The density of the flow was fixed at ρ=2000 kg m −3 . In the absence of accurate information about the input discharge we used a triangular shape with a total duration of 3240 s (54 min) and a total volume of 23 576 m 3 . Discharge increased from 0 to 15.86 m 3 s −1 at t=1080 s; it remained at this value until t=1440 s and then decreased progressively until t=3240 s. The total simulation time was extended until t=4000 s to allow for the flow to stop. These conditions mimic those used by Rickenmann et al. (2006), so the results of the different models could be compared.
Two simulation experiments were undertaken, using Bingham and Coulomb-viscous rheologies. A set of parameter values varying within a reasonable range were chosen, and successive model runs were performed. The results were compared with the field extension and thickness of the deposit in order to find the best set of parameters. For the Bingham rheology the best parameter set was τ C =2500 Pa and µ=525 Pa s. These results are close to those of Rickenmann et al. (2006), who used the Herschel-Bulkley rheology and found that the best results were obtained with a yield stress value within the range 0.8<τ C /ρ<1.35. In our simulation the value of this ratio was 1.25. Rickenmann et al. (2006) kept the µ/τ C ratio constant at 1/3, following Coussot et al. (1998). We found that the optimum value for this ratio was slightly lower, at approximately 1/5. In the case of the Coulomb-viscous rheology, the best parameter set was τ C =2500 Pa, µ=1300 Pa s and ϕ =4.8 • .
The shape of the final deposit was reasonably well predicted by both rheologies, although the results using the Bingham rheology were superior and almost perfectly matched the observed debris flow deposition boundaries (Fig. 3). However, the Coulomb-viscous model provided the best simulation of the thickness of the deposit. The thickest area was located in the upper central part of the fan, coinciding with the coarsest deposits observed. In the Bingham simulation the deposit thickness was more regular, with the thickest area located slightly more to the south of the fan. In both simulations the most problematic area was the zone around the upper part of the fan, where the simulated flow was thicker than the observed deposit, and there was some exceeding overflowing. Compared with the results of Rickenmann et al. (2006), the final shape of the debris flow in the present study provided a marginally better fit to the observed deposition. Another difference is that with respect to the Coulomb-viscous rheology, the thickest deposition was located further upstream, coinciding more closely to observations. However, these differences can be considered marginal and do not allow preferring one code over the other. In addition to the final deposition, map sequences of the flow thickness and velocity at different times (Fig. 4) were produced. The event was characterized by relatively slow velocity, which coincided with eyewitness reports and helped reduce the damage caused by the overflow.
The simulation proved to be more sensitive to changes in τ C and ϕ than to µ. As found by , it was difficult to make an independent calibration of the model, as the final results were relatively insensitive to small variations in the parameters and because of interactions among the effects of the parameters. For example, results close to the optimum were obtained using Mass-Mov2D with τ C =2200 Pa and µ=550 Pa s for the Bingham rheology, and with τ C =2500 Pa, µ=1400 Pa s, and ϕ =5.0 • for the Coulomb-viscous rheology.

The event of 5 August 2003 in the Faucon torrent
The Faucon torrent is located on the south-facing slope of the Barcelonnette Basin. It has a catchment area of 8 km 2 and an altitude ranging from 1170 to 2982 m a.s.l. At the bottom the torrent forms an alluvial fan with a slope of 0.12-0.08% (5-7 • ). The main outcrops are very finely laminated Jurassic black marls, which are mostly covered by moraines and/or slope deposits. The land cover is mainly forest (60%), abandoned agricultural lands, and bare lands with intense gullying. The Barcelonnette Basin has a dry and mountainous Mediterranean climate with strong inter-annual rainfall variability (733±412 mm over the period , strong storm intensities (over 50 mm h −1 ), and more than 100 days of freezing temperatures per year.
The Barcelonnette Basin is known to be prone to mass movements of various types including slow and fast landslides, and torrential debris flows . In recent times the Faucon torrent has been subject to debris flows; of particular note are the two most recent, in 1996 and 2003, for which good data are available. We chose the 2003 event for the modeling experiment, as it was the only one producing significant overflowing in the alluvial fan area. This debris flow event occurred after a summer storm, which generated a series of landslides in the upper catchment. According to eyewitnesses (Remaître, 2006;Remaître et al., 2009) the event consisted of four surge waves lasting about 3-5 min, separated by normal water flow. The first two surges were more liquid and stayed within the channel, with almost no overflow on the banks. The third surge carried a lot of debris, including wood, which blocked a bridge situated immediately upstream of the built sector. When the fourth surge arrived it hit the obstacle formed by the bridge and the accumulated debris. This caused major overflow on the banks, and destruction of the bridge. The total volume of the mudflow was estimated to be 25 000-30 000 m 3 , based on sediment budget analysis immediately after the event, and from geomorphologic fieldwork along the torrent (Remaître, 2006).
A pre-event DEM with a grid resolution of 1 m was used in the simulation. The DEM was generated from a 5 m contour lines map transformed to points data, enriched with GPS points. Since the area is an artificial channel with concrete banks, it was possible to use a thin plate interpolator to generate a smooth surface fitting the gaps between the points. A map defining the channel depth and its borders was used to represent topographical discontinuities. The location of the houses was determined from aerial photographs, and these were treated as obstacles. The density of the flow was set at ρ=1850 kg m −3 , following Remaître (2006). We simulated the fourth surge, which was associated with most of the overflow. The input discharge file consisted of one surge wave rapidly increasing from 0 to 16 m 3 s −1 in 60 s, and lasting for another 180 s for a total volume of 11 000 m 3 . Flow velocity at the inlet increased from 7 to 10 m s −1 during the surge, The Bingham and Coulomb-viscous rheologies were tested. As in the previous case, a calibration procedure was followed to determine the best set of parameters. Results for the shape of the deposits were good with both rheologies, although the Bingham rheology tended to under-predict the extension of the deposit area (Fig. 5), and the Coulomb-viscous rheology tended to over-predict it. Most of the deposit area was predicted to be close to the first bridge, where the obstacle was located. From that point a debris flow body advanced downstream on the left bank, where the houses were located (Fig. 6). On the lower right bank there was some flooding directly from the torrent; the shape of the deposition in this area was well reproduced. There was more overflow from the channel on the right bank between the two bridges, coinciding with the observed deposition. The thickness of the deposit was best predicted with the Coulomb-viscous rheology, as deposits more than 2 m thick were removed from the torrent banks near the bridge after the event finished. As in the previous case, the Bingham rheology resulted in a more homogenous and thin deposit. The best parameter sets were τ C =400 Pa and µ=67 Pa s (Bingham) and τ C =200, ϕ =3.8 • and µ=10 Pa s (Coulomb-viscous).
There are few reports of the modeling of canalized mud flows, but the Faucon 2003 event was modeled by Remaître (2006) using the BING code (Imran et al., 2001). The best parameter set for this simulation was τ C =404 Pa and µ=122 Pa s, which are very close to the best Bingham rheology parameters in this study. The results of both simulations are similar, although the final deposition maps of Remaître (2006) more accurately reproduced the deposition shape around the bridge. Our model correctly reproduced the observed deposition thickness, with most of the material remaining inside the channel, but in the results of Remaître (2006) the thickest deposition was concentrated around the houses on the left bank of the torrent. However, as these models were not run under equal initial and boundary conditions, comparison of the results is not definitive.

Discussion
Following a common approach for one-phase models, Mass-Mov2D assumes a homogeneous fluid with constant rheological properties. Although this assumption may be appropriate for mud flows, it is more critical in the case of debris flows where the separation of the solid and liquid phases is fundamental during the deposition phase. The inability to simulate the differing behaviors of the two phases is a limitation of one-phase models. Bingham-type flow models like the ones used in MassMov2D artificially simulate deposition of the debris flow by introduction of a yield stress value. Other authors have proposed to change the flow regime during the simulation, following pre-specified rules or rules imposed as external conditions (e.g. McDougall et al., 2006). Other alternative flow laws are rheological models based on the collisional dilatant model following the Bagnold theory (e.g. see Takahashi, 1991).
Another important assumption in the MassMov2D model is that the initial pore pressure is maintained during the simulation. This approximation is valid for flows where the consolidation time is much greater than the propagation time (e.g. mud flows), but this may not be the case during the deposition phase of debris flows, where the pore pressure dissipation becomes the main controlling process (Major and Iverson, 1999). The model could be further extended by incorporating a mechanism for pore pressure dissipation during the flow. This has proved to be difficult under the assumption of a one-phase flow, although some authors have proposed parametric models to simulate the hysteretic behavior of granular flows (Pouliquen and Forterre, 2002;Mangeney et al., 2007). An alternative solution may be to include the coupling between the solid and the fluid fractions in order to model transient pore pressures, based on the grainfluid mixture theory (Iverson, 1997;Iverson and Denlinger, 2001;Denlinger and Iverson, 2004). Such approaches, however, have not yet generally been implemented in engineering practice, and are beyond the original intentions of Mass-Mov2D.
The assumption of a fixed bed may also not be suitable for debris flows, for which erosion by scouring has been described as an important source of material. Sediment entrainment can also affect the properties of the debris flow, and thus greatly affect its behavior. A new version of MassMov2D, which will allow for a non-fixed bed, is currently under development.

Conclusions
MassMov2D is a numerical simulation model of the dynamics (runout and deposition) of highly concentrated mass movements over complex topography. Due to the difficulties involved in obtaining direct information about the characteristics of debris and mud flows, numerical models are useful tools for hazard evaluation and mitigation, as they allow estimation of relevant characteristics including the flow depth and the impact force.
The model proposed here is based on a depth-averaged form of the equations of motion for a fluid continuum, solved in a 2-D finite differences mesh. It can incorporate various constitutive relations for the flow resistance term, which allows the model to adapt to materials having different characteristics. The main drawback of the model is the assumption of a constant rheology. Further iterations of the model should address important issues, including the decrease in pore pressure during the simulation, allowing for a much more realistic simulation of debris flow behavior. However, this issue will remain a challenge for most one-phase codes.
In contrast with other codes for debris flow modeling, MassMov2D is implemented in a GIS environment. This facilitates the tasks of preparing the input files and analyzing the model outputs, resulting in a user-friendly interface that facilitates use by non-experts in numerical computing. Moreover, as the core equations of the model are implemented in an easy-to-learn scripting language, the model is especially suited to testing new modeling concepts.
The key inputs to the model can be created and manipulated using common GIS tools, making it is easy to develop simulation scenarios to fit real case studies. In this study we considered torrent debris flows in Wartschenbach (Austria) and Faucon (France). Through a back-analysis process to establish the best set of parameters, the model was able to simulate the final spread and thickness of the debris flow deposits. These results provide evidence that the model can accurately simulate real events, and justify its use as a tool in hazard analysis.