Physically-based modelling of granular flows with Open Source GIS

Computer models, in combination with Geographic Information Sciences (GIS), play an important role in up-to-date studies of travel distance, impact area, velocity or energy of granular flows (e.g. snow or rock avalanches, flows of debris or mud). Simple empirical-statistical relationships or mass point models are frequently applied in GISbased modelling environments. However, they are only appropriate for rough overviews at the regional scale. In detail, granular flows are highly complex processes and physicallybased, distributed models are required for detailed studies of travel distance, velocity, and energy of such phenomena. One of the most advanced theories for understanding and modelling granular flows is the Savage-Hutter type model, a system of differential equations based on the conservation of mass and momentum. The equations have been solved for a number of idealized topographies, but only few attempts to find a solution for arbitrary topography or to integrate the model with GIS are known up to now. The work presented is understood as an initiative to integrate a fully physicallybased model for the motion of granular flows, based on the extended Savage-Hutter theory, with GRASS, an Open Source GIS software package. The potentials of the model are highlighted, employing the Val Pola Rock Avalanche (Northern Italy, 1987) as the test event, and the limitations as well as the most urging needs for further research are discussed.


Introduction
Granular flows -including avalanches of snow, mud, debris or even rocks -are highly destructive phenomena putting people, buildings and infrastructures at risk.Delineation of possible impact areas as well as flow velocities is an essential precondition for efficient action towards risk reduction, e.g.designation of hazard zones for land use planning, dimensioning of technical structures etc. (Hungr et al., 2005;Pudasaini and Hutter, 2007).
Since the 1990s, Geographic Information Sciences (GIS) have played an increasing role in the mapping and prediction of landslide, debris flow, and avalanche hazard and risk.They enable an efficient management of spatial data at all scales, usually in raster or vector format.On a regional or even national level, GIS have a high popularity for identifying areas with high landslide hazard and risk and for generating so-called landslide susceptibility maps.A large array of statistical methods has been applied for this purpose (e.g.Corominas et al., 2003;Lee, 2004;Guzzetti et al., 2006).Also regarding physicallybased methods, which require more detailed topographic and geotechnical information and are therefore mainly applicable at the small catchment scale, various approaches have been worked out and several GIS-based studies on shallow slope failures have been published (e.g.Godt et al., 2008).Programs like SINMAP (Pack et al., 1998) or SHALSTAB (Dietrich and Montgomery, 1998) work as extensions to proprietary GIS software packages.Mergili and Fellin (2011) implemented a model for rotational failures with the Open Source package GRASS GIS (GRASS Development Team, 2011).However, GIS-based model development has focused to a lesser extent on the movement itself, but rather on the onset of the movement.
Models for the motion of rapid mass movements such as granular flows may be classified as shown in Table 1.Most of the existing approaches are potentially suitable for more than one type of flow (snow avalanche, debris flow, rock avalanche etc.).Many of them were developed independently from GIS, which are just a mode of implementation offering efficient spatial analysis tools.The main types of models are shortly outlined below, partly following Pudasaini and Hutter (2007).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Table 1.Types of models for granular flows, partly based on Pudasaini and Hutter (2007).
Mass point models (sometimes referred to as semideterministic or conceptual models) for snow avalanches or debris flows go back to the work of Voellmy (1955), who related the shear traction at the base to the square of the velocity and assumed an additional Coulomb friction effect (Pudasaini and Hutter, 2007).This approach considers only the centre of the flowing mass but not its deformation and the spatial distribution of the flow parameters.The Voellmy (1955) model was further developed and applied e.g. by Perla et al. (1980); Gamma (2000); Wichmann and Becht (2003); and Mergili et al. (2011).With statistical and mass point models, random walk procedures based on Monte Carlo techniques (Gamma, 2000) are used for routing the flow.As Van Westen et al. (2005) showed, the influence of the routing algorithm on the modelled areas of deposition is considerable.
Physically-based distributed models require more specialized input information and are therefore suitable for localscale studies.Movements are computed based on physical laws, assuming specific flow rheologies (Savage and Hutter, 1989;Takahashi et al., 1992;Hungr, 1995;Iverson, 1997;Pudasaini and Hutter, 2003;McDougall andHungr, 2004, 2005;Pudasaini et al., 2005a, b).GIS applications of such models are still rather scarce.Some examples exist where numerical models can be coupled with proprietary GIS software, e.g.FLO-2D (O'Brien, 2003), SAMOS (Sampl and Zwinger, 2004) or an application of the Takahashi et al. (1992) model (Chau and Lo, 2004).Some other models have been used in combination with GRASS GIS, like DAN (Hungr, 1995;Revellino et al., 2008) or TI-TAN2D (Pitman et al., 2003).The RAMMS model (Christen et al., 2010a, b), primarily applied to snow avalanches, uses Voellmy (1955) viscous drag and runs with GRASS GIS in the background.However, a full GIS implementation of more complex models (e.g.Savage and Hutter, 1989;Iverson 1997;Pudasaini and Hutter, 2003; see Pudasaini and Hutter, 2007 for a comprehensive list of references) still remains a challenge.Recently, there has been rapid progress in theoretical modelling, numerical simulation and validation of the models with experimental and field data for granular flows (Gray et al., 1999;Pudasaini et al., 2005a;Pudasaini et al., 2008;Luca et al., 2009a, b, c, d;Pudasaini, 2011b).Zwinger et al. (2003) adapted the Savage-Hutter model for arbitrarily shaped mountain terrain, introducing some geometric simplifications.Based on Nessyahu and Tadmor (1990) and Tai et al. (2002), Wang et al. (2004) devised a numerical scheme applicable to flows over complex shallow terrain along a horizontally straight flow path (Gray et al., 1999).Pudasaini and Hutter (2003) proposed an avalanche model for generally curved and twisted channels which was further developed by Pudasaini et al. (2005aPudasaini et al. ( , b, 2008)).The simulation results were verified with experimental data.The Gray et al. (1999) and the Pudasaini and Hutter (2003) model are both based on curvilinear coordinate systems.Whilst such general models are adequate for the problems related to arbitrarily shaped mountain terrain, substantial efforts are still needed to apply these equations to real-world mountain topographies in connection with conventional raster-based GIS environments.This requires some adaptations which will be explained later.
The present paper describes an attempt to implement a fully physically-based model for the motion of granular flows with GIS.The software GRASS GIS was selected for this purpose.As an Open Source GIS package with focus on raster processing, GRASS facilitates model development, distribution and evaluation: the program code is freely accessible and new modules may be added by anybody.Therefore, the entire scientific community has the chance to contribute to the evaluation of the model and the further development of the program code.Examples of GRASS implementations of mass movement models are provided e.g. by Cannata and Molinari (2008) and Mergili et al. (2011).
The GRASS module developed was named r.avalanche.It builds on a numerical scheme devised by Gray et al. (1999), Pudasaini (2003), andWang et al. (2004) for simple concave topographies with an only vertically curved flow line (Fig. 1).Pudasaini et al. (2005aPudasaini et al. ( , b, 2008) ) further applied these numerical schemes to more complex curved and twisted topographies, and compared the model simulations with laboratory and field data for dry avalanches and mixture debris flows.Future work should be directed towards extending the scope of the physical and numerical model to arbitrary mountain terrain.

Model basics
The extended Savage-Hutter Model (Savage and Hutter, 1989;Gray et al., 1999;Pudasaini and Hutter, 2003) is expressed as a system of partial differential equations describing the conservation of mass and momentum for the rapid motion of dry granular avalanches: where h is the avalanche thickness, and u and v are the depthaveraged downslope and cross-slope velocities.s x and s y are the net driving accelerations: where ζ is the downslope inclination angle of the reference surface, κ is the local curvature of the reference surface, b is the elevation above the reference surface, and ε and λ are some scaling factors (see Sect. 2.3.2). β x and β y are defined as β y = εcosζ K y .
(7) K x and K y are the earth pressure coefficients in downslope and cross-slope directions (as computed at the basal surface): φ is the angle of internal friction, and δ is the bed friction angle.Active stress rates (subscript act) are connected to local elongation of the flowing mass, passive stress rates (subscript pass) are connected to local compacting -it depends on acceleration or deceleration of the flow whether active or passive stress rates are applied (Gray et al., 1999).The system of equations described above is only valid for the shallow flow of cohesionless and incompressible granular materials which can be considered as continuum.The model is often applied to cases with non-shallow onset areas.This is acceptable as long as the movement quickly transforms into a shallow flow, since the geometry during the flow and the deposition does not much depend on the geometry of the onset mass in such cases.The curvature of the terrain has to be relatively small.It has to be emphasized that all variables are dimensionless.The model is scale-invariant and small-scale laboratory tests can be used as reference for large-scale problems in nature.

Numerical solution of the equations
The differential equations Eqs.
(1) to (3) were solved using a high resolution Total Variation Diminishing Non-Oscillatory Central Differencing (TVD-NOC) Scheme, a numerical scheme useful to avoid unphysical numerical oscillations (Nessyahu and Tadmor, 1990).Cell averages of h, u and v are computed using a staggered grid: the system value real-time period after which the simulation stops x path1 , y path1 , x path2 , y path2 (m) values coordinates of two points defining the main flow line is moved half of the cell size with every time step, the values at the corners of the cells and in the middle of the cells are computed alternatively.The TVD-NOC scheme with the minmod limiter has been successfully applied to a large number of avalanche and debris flow problems (Tai et al., 2002;Wang et al., 2004;Pudasaini et al., 2005a, b, 2008, Pudasaini and Kroener, 2008).Here, this scheme is implemented in the r.avalanche computational tool.
The simulation is run for a user-defined, real-time period; in the future it is planned to implement a criterion based on the absolute value of the average velocity of the flow, with a user-defined sufficiently low threshold value.The time steps have to be kept short enough to fulfil the Courant-Friedrichs-Levy (CFL) condition: where x is the cell size in x direction, c max is the global maximum wave speed and i and j are the coordinates of the cells in x and y direction, respectively.t is determined dynamically during the run-time of r.avalanche, based on the CFL condition from the previous time step.

Layout
r.avalanche was developed as a raster module for GRASS GIS, using the C programming language.Raster maps and a set of parameters compiled in a text file serve as input (see Table 2).The implementation of the model and the numerical scheme described above into GIS should consider the following: -Eqs.(1) to (3) provide dimensionless values.However, in real flow simulations, for better physical intuition, it is desirable to use dimensional values; -in general the solution is obtained for a curvilinear reference system.However, this must be converted into a rectangular GIS system.
In reality the avalanche valleys are not only curved in the downslope direction, but also in the cross-slope direction.
Here, for simplicity, the idealized situation that the slope is only curved in the downslope direction but laterally flat is assumed (Gray et al., 1999).However, as in Pudasaini et al. (2005aPudasaini et al. ( , b, 2008)), it would be more realistic to construct the mountain topography based on a curved and twisted surface or even better, on arbitrary topography.This could be an interesting future direction, but is not within the scope of this paper.

Non-dimensionalization of the variables
The first issue concerns the non-dimensionalization of the governing variables and parameters, which is explained e.g. in Pudasaini (2003).With the typical avalanche length L, the typical avalanche depth H and the typical radius of curvature R, the dimensional variables are derived by using x, y = (Lx,Ly) where the variables denoted with a cap are the dimensional counterparts of the variables used in Eqs.
(1) to (3) and g is gravitational acceleration.The implemented model computes the dimensionless variables and converts them to dimensional values according to Eq. ( 12) for output.The scaling parameters L, H , and R are set to 1 in this paper.The factors ε = H /L and λ = L/R (Pudasaini, 2003) are therefore 1 as well.However, the choice of the scaling parameters does not influence the numerical results since Eqs.(1) to ( 9) are scale invariant.

Adaptation of the coordinate system
In general the extended Savage-Hutter model and its solution are formulated by using a curvilinear coordinate system based on a talweg which can be generally curved and twisted (Wang et al., 2004;Pudasaini et al., 2005a, b;2008).The talweg is considered to be the predominant flow direction.Here, as in Gray et al. (1999) and Wang et al. ( 2004), it is assumed that the talweg is only curved in the downslope direction so that its projection to a horizontal plane is a straight line (see Fig. 1).Three steps are required for converting the original rectangular coordinate system of the input raster maps into the coordinate system for the simulation: -the coordinate system is rotated around the z axis so that the direction of the main flow line is aligned with the new x (downslope) direction.The main flow line is based on two user-defined pairs of coordinates (see Table 2); -a reference surface is created, defined by the given talweg and an inclination of zero in y (cross-slope) direction.As the talweg follows the natural terrain, its longitudinal section (Fig. 2) does not represent a straight or regularly curved line like in the idealized topography shown in Fig. 1; -based on this reference surface, the cell size x for each x parallel to the reference surface is computed.The offset b (m) -defined as the distance between terrain surface and reference surface perpendicular to the reference surface -is derived.Since the terrain surface is not given analytically and z 0,terrain in Fig. 2 depends on the corresponding x 0 , b has to be determined iteratively.This is done by varying the horizontal shift until the tested value of z 0,test converges to the terrain surface.
The terrain height is interpolated between the centres of the two closest raster cells in the direction of x 0 .Initial avalanche thickness h is derived in an analogous way.As a final step, all raster maps are resampled in order to set x = y = const.(a precondition for applying the numerical scheme used by Wang et al., 2004).
After completing the simulation, the entire system is reconverted into the rectangular coordinate system used in the  GRASS GIS mapset in order to enable a proper display of the results.

The 1987 Val Pola Rock Avalanche
The Val Pola Rock Avalanche (Valtellina, Lombardy Region, Northern Italy, Fig. 3a), like most large mass movements, has been subject of numerous detailed studies.A review with the relevant references is provided by Govi et al. (2002).The good state of knowledge on the geometric characteristics as well as the mechanics and kinematics of the event makes it suitable for model evaluation.Granular flow modelling has already been applied there (e.g. by Crosta et al., 2003Crosta et al., , 2004)).
The event (also referred to as Val Pola Landslide) occurred on 28 July 1987 after a period of heavy and persistent rainfall.During this period, erosion was intense, and numerous debris flows were reported from the tributaries to the Valtellina.The Val Pola Creek was deeply eroded, undercutting the northern edge of an old landslide body.Preceded by the opening of a prominent crack, a block of highly fractured and faulted igneous and metamorphic rocks detached and rushed into the valley (Crosta et al., 2003).Figure 3b provides a comparison of satellite views of the area before and after the event.
The volume of the released mass was estimated to 34-43 million m 3 (Crosta et al., 2003) and 35 million m 3 (Govi et al., 2002), the entrainment by the resulting rock avalanche to further 5-8 million m 3 .After detaching, the mass first moved northwards until rushing against the Sassavin-Motta Ridge, and then proceeded as rock avalanche eastwards down to the main valley.Govi et al. (2002) distinguished six phases of the movement, with a duration of the main avalanching phase of 8-12.5 s.The velocity of the flow was estimated to peak at 76-108 m s −1 (Crosta et al., 2004), indicating that the event was one of the most rapid documented mass movements in history.
The mass moved up 300 m on the opposite slope.In the main valley, it continued 1.5 km downstream and 1.5 km upstream, with a maximum thickness of 90 m (Crosta et al., 2003; see Fig. 3b).The Adda river was dammed, and a lake started to form, the level rising several metres per day.An artificial drainage was constructed rapidly in order to avoid an uncontrolled sudden drainage of the lake.The event claimed 27 human lives and caused high economic costs (Crosta et al., 2003).

Data
The following data sets were obtained for the Val Pola Rock Avalanche: -Digital Elevation Models (DEMs) before and after the event; -geotechnical parameters of the sliding mass and the sliding surface; -reference information on the distribution of the deposit.
A DEM provided by the Geoportal of the Region Lombardy (Italy; www.cartografia.regione.lombardia.it) at a cell size of 20 m was used.It represents the situation before the 1987 Val Pola event.An SRTM-4 DEM (Jarvis et al., 2008) from February 2000 was obtained in order to capture the situation after the event.It has a cell size of 3 arc seconds and was resampled to 20 m, too.However, it was mainly used for defining the depth of the failure plane whilst the pre-failure DEM was applied for the simulation.
The distribution of the released mass was estimated from the difference between the two DEMs.A volume of approx.40 million m 3 was derived, being well within the range reported by Crosta et al. (2003Crosta et al. ( , 2004)).The higher value compared to that one reported by Govi et al. (2002) can be explained by the fact that material detached from the Sassavin-Motta Ridge -considered as entrained material by those authors -is included here (Fig. 4).Furthermore, the defined onset area is larger than shown by Crosta et al. (2003).This counterbalances the fact that entrainment is disregarded by r.avalanche.In general, it is hard to clearly delineate the zones of onset and entrainment by the mass flow.The accuracy of the DEMs did not allow for a reliable estimate of entrained and deposited volume.However, the maximum thickness of deposited material derived by the DEM overlay was 83 m and therefore close to the value reported by Crosta et al. (2003) and Govi et al. (2002).However, this overlay has to be considered as a rough approximation since significant anthropogenic reshaping of the terrain took place after the event.Only areas with a derived thickness of the deposit ≥20 m are shown in Fig. 4. The affected area as reference data was mapped from Landsat TM imagery and verified with the Italian Landslide Inventory (IFFI Project) and a map published by Crosta et al. (2003).Possible inaccuracies are caused by secondary processes, e.g.mud flows downstream blurring the delineation of the impact area downvalley.
r.avalanche was run with different combinations of parameter sets (Table 3), partly following Crosta et al. (2003Crosta et al. ( , 2004)).Crosta et al. (2003) back-calculated an angle of internal friction φ = 45 • for the released mass, and a bed friction angle δ = 18 • for the sliding surface.However, they used  2004) back-calculated the event in a three-dimensional way with a Lagrangian Finite Element Model (LFEM), assuming a basal friction angle of 22 • and a dynamic internal friction angle of 35 • .These values were used as a base also for the present study.However, it has to be emphasized that the parameters were not directly derived from laboratory tests but indirectly by fitting the model results to the observations.The standard cell size for the computation was set to 20 m, according to the DEM used.The model was also tested with 40 m and 12 m cell size, respectively.Furthermore, two different assumptions of the talweg (see Fig. 4) were tested: talweg 1 starting farther N than the centre of the affected area (accounting for the fact that the detached mass moved northwards before converting into a rock avalanche), and talweg 2 along the centre of the affected area.All simulations were run for a real-time duration of 300 s, following Crosta et al. (2004) and the observation that no substantial temporal changes of the height and geometry of the modelled deposit were observed at that time.
The single-phase model used does not allow for considering landslide-river interactions that may have occurred when the rock avalanche plunged into the Adda river.However, future model development will go in this direction (see Conclusions).

Test with simple artificial topography
Before being applied to the Val Pola Rock Avalanche, r.avalanche was tested for a simple plane slope with an inclination of 35 • , running out into a horizontal plane.Assuming a hemispherical starting mass, the simulation was run with four combinations of φ (30 • and 37 • ) and δ (23 • and 28 • ).The model responded much less sensitively to variations of

Simulation with published friction parameters
The model was first run at a cell size of 20 m (the same as the DEM resolution) with φ = 35 • and δ = 22 • , assuming talweg 1 (Simulation 1; see Table 3 and Fig. 4).This simulation can be considered as a validation with the parameters used by Crosta et al. (2004).Figure 6 illustrates the distribution of flow depth at selected time intervals from the onset (t = 0) to t = 300 s, when the flow is assumed having stopped.86 % of the modelled deposit coincide with the observed area of deposition, and 79 % of the observed deposition area are occupied by the modelled deposit (considering only those areas with depth of deposition >1 m; Table 4).The travel distance modelled for the central part of the deposit corresponds very well to the observation.
However, some issues become obvious from Fig. 6: 1. some portions of the lateral (northern and southern) parts of the modelled deposit deviate from the observed deposit.Particularly in the north, the travel distance is underestimated;  2. the depth of the deposit is underestimated: whilst the measured maximum depth is approx.90 m (Govi et al., 2002), the model yields a maximum depth of 75 m only; 3. already shortly after the onset of the flow (t = 10 s), the simulation predicts a degree of lateral spreading of the flow beyond the observed extent.As a consequence, part of the material crosses the delineation of the catchment and follows gullies north and south of the  main flow path, a behaviour not at all observed in reality.
The lake dammed by the rock avalanche deposit was modelled by filling the sink behind the simulated deposit (see t = 300 in Fig. 6).The extent of the modelled lake corresponds well to the observation on 20 September 1987 (see Fig. 3b), which represents the lake close to its maximum extent.
Besides flow depth, flow velocity is another important characteristics of granular flows (Pudasaini et al., 2005b(Pudasaini et al., , 2007;;Pudasaini and Kroener, 2008;Pudasaini, 2011b).Flow velocity is more difficult to verify than flow depth distribution since direct reference measurements are hardly applicable.Figure 7 illustrates a time series of the distribution of depth-averaged flow velocity.The maximum simulated velocity is 94 m s −1 at t = 34 s.This value is well within the estimated range of 76-108 m s −1 (Crosta et al., 2004), but appears very high, anyway.Very rapid mass movements may reach 100 m s −1 and more (e.g.Scheidegger, 1973;Highland and Bobrowsky, 2008).Evans et al. (2009b) modelled a maximum velocity of approx.80 m s −1 for the Khait Rock Avalanche in Tajikistan, characterized by a much lower average slope angle than the Val Pola event, with DAN.However, velocities of rock avalanches comparable to the Val Pola event typically range from 30-50 m s −1 .A two-dimensional model yielded a maximum of 50 m s −1 for the Val Pola event (Crosta et al., 2003).Usually only rock-ice avalanches and related processes move at much higher velocities (e.g.Evans et al., 2009a).
The values used for φ and δ -and their spatial distribution in particular -are uncertain.Therefore, the sensitivity of the model to variations of these governing parameters has to be evaluated.r.avalanche is supposed to be sensitive also to changes of the assumed talweg and to the cell size used for the simulation.Table 4 summarizes the key output parameters for each simulation and the correspondence of modelled and observed deposit.

Sensitivity to cell size
The simulation was repeated with 40 m (Simulation 2) and 12 m cell size (Simulation 3), leaving the remaining parameters unchanged (see Table 3).Longitudinal profiles comparing the results to those yielded with 20 m cell size as well as the maps of the simulated deposits are shown in Fig. 8.It was observed that the variations of the cell size do not lead to a substantial change in the behaviour of the simulated rock avalanche.However, with coarser resolution, the spreading of the flow is more pronounced and the simulated flow mass is smoothed, leading to a larger impact area but smaller depth of flow and deposit (see Table 4).With 12 m resolution, the predicted maximum depth of the final deposit (85 m) comes close to the observation (approx.90 m), whilst with 40 m resolution, it is clearly underestimated (39 m).
This effect is particularly visible when comparing the results for 40 m and 20 m and much less pronounced between 20 m and 12 m.A maximum flow velocity of 86 m s −1 was predicted with 40 m cell size, 96 m s −1 with 12 m cell size.
These findings are not surprising as the cell size governs the distance of spreading during each time step of the simulation.Tests with channelized debris flows have shown that r.avalancheonly works well if the cell size is much smaller than the width of the flow, otherwise lateral spreading is overestimated.A cell size of 40 m is definitely too coarse for the simulation of the Val Pola Rock Avalanche.The less pronounced difference between the results yielded with 20 m and 12 m suggests that these values are close to the cell size ideal for the simulation of this specific event (choosing a too fine resolution would unnecessarily increase computing time).

Sensitivity to assumption of talweg
Whilst the talweg is clearly defined in the idealized topography suggested by Gray et al. (1999) and Wang et al. (2004), the identification of the talweg is non-trivial for the real terrain.A talweg for the simulation has to be defined manually, based on the geometry and the mechanics of the flow.Assumption 1, used for most of the simulations presented here (see Table 3), builds on the fact that the detached mass first moved towards north before rushing down as a rock avalanche: therefore, the starting point of the talweg  is shifted north from the centre line of the affected area (see Fig. 4).Assumption 2 works with a talweg following the centre line of the affected area and was tested against assumption 1 (Simulation 4).The result (comparison of time series of flow depths) is shown in Fig. 9.
As prescribed by the scheme used, the change in the assumption of the talweg leads to a rotation of the dominant flow direction, adopting to the direction of the talweg (x-axis in Figs. 1 and 2).Even though the flow follows the talweg assumption 1 in its initial and intermediate stages, this effect becomes particularly significant at the valley bottom.The counter-clockwise rotation of the system, compared to Simulation 1, is best visible in the southernmost part of the deposit.In this particular case, the modelled travel distance is only slightly shorter with talweg assumption 2, and a depth of the final deposit of 72 m is predicted instead of 75 m (see Table 4).Despite the comparable outcome regarding these key parameters, the different shapes of the deposits indicate that there are some cross-slope curvature effects and these effects should have been included (Pudasaini et al., 2005a(Pudasaini et al., , b, 2008)).However, only the curvature in the downslope direction is considered here.

Sensitivity to friction parameters
The simulation was repeated with δ = 18 • (Simulation 5) and 26 • (Simulation 6; Fig. 10; see Table 4).The expected effects of reduced flow velocity and resulting delayed motion with increasing values of δ were observed (Fig. 11).With δ = 18 • , the predicted flow velocity peaks at 107 m s −1 , with δ = 26 • it reaches 80 m s −1 only.Whilst the simulated maximum travel distance is longer when assuming a lower bed friction (t = 40 s and t = 60 s in Fig. 10a), the final travel distance at t = 300 s is almost the same with δ = 18 • and with δ = 22 • (see Fig. 10b).This phenomenon is explained by the higher tendency of the material to flow back into the main valley with lower bed friction.Also with δ = 26 • , the modelled final deposit reaches approx.as far east as in the other simulations.
An overall look at the deposits shown in Fig. 10b shows that varying the bed friction angle may result in a highly nonlinear response governed by several factors (see Fig. 10): with δ = 18 • , the maximum thickness of the deposit decreases to 67 m (see Table 4), compared to δ = 22 • (75 m), an effect to be attributed primarily to increased lateral spreading (north-south direction).With δ = 26 • , part of the flow  material would remain in the transit and even onset area.The deposit in general would therefore assume a much more stretched shape in east-west direction, with a maximum depth of 69 m.However, with δ = 26 • , the modelled deposit extends farther south than even with δ = 18 • .This counterintuitive effect illustrates that the numerical model is not designed to appropriately describe flows on curved flow paths.In this specific case, this affects only the lateral parts of the flow (flow height ≤5 m) in the downstream direction, larger flow heights are modelled in a plausible way.
In general, lower values of φ should lead to a thinned deposit due to an increased tendency to spread out, higher values of φ should have the reverse effect.Also here, the situation is more complex in natural terrain, particularly in narrow valleys (Fig. 12): internal friction angles of φ = 40 • (Simulation 7) and 45 • (Simulation 8) were tested against the value of φ = 35 • used by Crosta et al. (2004).The lowest value of φ caused the model to predict the largest maximum depth of the deposit (75 m).This actually unexpected outcome is a consequence of the stronger tendency of material with reduced internal friction to level out and therefore fill up the valley bottom.When assuming higher values of φ, there is less tendency of the modelled rock avalanche to level out the valley bottom, and the deposit is spread over a larger area, but with a lower maximum depth (67 m for φ = 40 • and 60 m for φ = 45 • ).For laboratory-scale flows, the effects of the basal and internal friction angles are analyzed in detail in Pudasaini and Kroener (2008).

Conclusions
r.avalanche represents a first attempt to implement a physically-based, distributed granular flow model with Open Source GIS.With respect to travel distance and impact areas, the physical model considered and computational tool developed are potentially suitable for real avalanche flows.With its current layout, the model is applicable to large-scale mass movements rather than to channelized debris flows.If the width of the flow would be much smaller than its length, a very small cell size (much smaller than the flow width) would be required in order to avoid an unrealistically high degree of lateral spreading.Such a small cell size, however, would lead to unacceptably long computing times.The travel distance and partly also the impact area of the Val Pola Rock Avalanche were well reconstructed without re-calibration of the input parameters.In contrast, the maximum depth of the deposit was rather underestimated.The friction parameters used by Crosta et al. (2004) in combination with the talweg assumption 1 lead to the most realistic model results regarding travel distance along the talweg and maximum depth of the deposit.However, also the other simulations yield parameters of comparable or, for specific parameters, even better quality (see Table 4).The impact area is best predicted when assuming φ = 45 • (Simulation 8).Most simulations failed to reconstruct the northern distal part of the impact area but the travel distance in the central part was predicted very well.The same phenomenon is visible in the simulation results of Crosta et al. (2004), where the underestimation of the northern portion of the impact area is even more pronounced than in the present work.
Whilst a cell size of 40 m is definitely too coarse, the results are only moderately sensitive to a variation between 12 m and 20 m, so that 20 m cell size is considered acceptable.The modelled flow velocities of all simulations are in the range of published estimates (Crosta et al., 2004) which, however, appear high when compared to those specified for other large and rapid mass movements (see Scheidegger, 1973, for examples).
The quality of the results may be limited by uncertainties in the governing parameters and their spatial patterns.However, as it stands now, there are some limitations of the physical model and the numerical model used: -r.avalanche is primarily applicable to granular flows with approximately horizontally projected straight flow lines, as this was a basic assumption in the derivation of the model equations for flows over super-imposed complex terrain by Gray et al. (1999).However, this concept can further be extended to more general topographies.
Whilst the requirement of a straight flow line is fulfilled by many large mass movements, this assumption fails for most channelized debris or mud flows where the travel distance is much larger than the flow width.For such processes, the Pudasaini and Hutter (2003) model, developed for curved and twisted channels, could be considered (Pudasaini et al., 2005a(Pudasaini et al., , b, 2008).An alternative model for arbitrary topography can also be considered.However, it still needs considerable efforts to develop such a model and to apply it to raster-based GIS; -entrainment of regolith during the flow is disregardedin the case of the Val Pola Rock Avalanche, entrainment was not extremely significant.However, in general, entrainment plays a major role for the travel distance of granular flows (e.g.McDougall and Hungr, 2005;Sovilla et al., 2006Sovilla et al., , 2007;;Quan Luna et al., 2012).Therefore, high priority should be given to this aspect when pushing the model development further; -the role of pore water for the motion of the flow is neglected, but many geophysical mass flows in nature are a mixture of solid and fluid components and should be considered as such (Iverson and Denlinger, 2001;Pudasaini et al., 2005a).
Attacking the limitations in a comprehensive way will require at least the following steps: 1. selecting and adapting a sound method for modelling rapid granular flows over arbitrary topography, using and extending the existing theories.Such an approach would have to build on the latest extensions of the Savage-Hutter model (Luca et al., 2009a, b, c, d), incorporating particle entrainment and the role of pore fluid (besides water, also air may play an important role).As explained in Pudasaini (2011a), the rheological model should take into account different dominant physical aspects as observed in two-phase geophysical mass flows.
Such aspects include the solid volume fraction gradient enhanced non-Newtonian fluid extra stress; the generalized interfacial momentum transfer that takes into account the viscous drag, buoyancy, and the virtual mass; and also the generalized drag which covers both the solid-like and fluid-like contributions.Such a general model is essential to describe e.g.complex landslideriver/lake interactions.
2. devising an appropriate high resolution shock capturing numerical scheme for solving the differential equations derived in (1)-(3).Numerical solutions of the avalanche model equations for arbitrary topography would have to be elaborated (e.g.Bouchut and Westdickenberg, 2004;Pudasaini et al., 2005aPudasaini et al., , b, 2008)); 3. performing more simulation tests and calibrating the model with well-studied granular flow events.
Since no user-friendly Open Source software for the motion of granular flows (fully incorporating the relevant physical and geometrical processes) is available at present, the further development of r.avalanche in terms of the above points would be highly relevant for reliable process modelling and the delineation of hazard zones.

Fig. 1 .
Fig. 1.Idealized topography for the solution of the extended Savage-Hutter Model by Wang et al. (2004), also used for r.avalanche.

Fig. 2 : 30 Fig. 2 .
Fig. 2: Iterative determination of b based on the coordinate system defined by the reference plane.
Fig. 3: Val Pola Rock Avalanche: (a) location; (b) comparison of Landsat ima after the event.

Fig. 3 .
Fig. 3. Val Pola Rock Avalanche: (a) location; (b) comparison of Landsat imagery before and after the event.

Fig. 4 .
Fig. 4. Val Pola Rock Avalanche.Distribution of detached and deposited volumes as derived from the DEM change detection, the flow trimline was mapped from Landsat TM imagery (thick black line).
Fig. 8: (a) Longitudinal profiles of the flow along the talweg at different time s φ = 35° and δ = 22°: comparison of the results for 12 m, 20 m and 40 m cell sizefor 20 m is also shown proportionally to the terrain; (b) flow depth distribution after computed at three different cell sizes.

Fig. 8 .
Fig. 8. (a) Longitudinal profiles of the flow along the talweg at different time steps with φ = 35 • and δ = 22 • : comparison of the results for 12 m, 20 m and 40 m cell size -the result for 20 m is also shown proportionally to the terrain; (b) flow depth distribution after t = 300 s, computed at three different cell sizes.
Fig. 10: (a) Longitudinal profiles of the flow at different time steps with φ = 35°: co of the results with δ = 18, 22 and 26° with a cell size of 20 m -the flow depth for also shown proportionally to the terrain; (b) flow depth distribution at t = 300 s, with three different assumptions of δ.

Fig. 10 .
Fig. 10.(a) Longitudinal profiles of the flow at different time steps with φ = 35 • : comparison of the results with δ = 18, 22 and 26 • with a cell size of 20 m -the flow depth for δ = 22 • is also shown proportionally to the terrain; (b) flow depth distribution at t = 300 s, computed with three different assumptions of δ.

Fig. 11 :Fig. 11 .
Fig. 11: Velocity profiles along the assumed talweg with different assumptions of δ for selected time steps and for the maximum over the entire flow.The values are positive in downslope direction (from left to right).
Fig. 12: (a) Longitudinal profiles of the flow at different time steps with δ = 22°: co of the results with φ = 35, 40 and 45° with a cell size of 20 m -the flow depth for φ also shown proportionally to the terrain; (b) flow depth distribution after t = 300 s different assumptions of φ.

Fig. 12 .
Fig. 12.(a) Longitudinal profiles of the flow at different time steps with δ = 22 • : comparison of the results with φ = 35, 40 and 45 • with a cell size of 20 m -the flow depth for φ = 35 • is also shown proportionally to the terrain; (b) flow depth distribution after t = 300 s for three different assumptions of φ.

Table 2 .
Input parameters for r.avalanche.

Table 3 .
Parameter combinations applied for the simulations.
Crosta et al. (2003)el, whilst the Savage-Hutter Model is based on a Mohr-Coulomb model.Therefore, and sinceCrosta et al. (2003)used a two dimensional model with plane strain condition, these values are not directly applicable to the present study.Crosta et al. (