Flow-R , a model for susceptibility mapping of debris flows and other gravitational hazards at a regional scale

The development of susceptibility maps for debris flows is of primary importance due to population pressure in hazardous zones. However, hazard assessment by processbased modelling at a regional scale is difficult due to the complex nature of the phenomenon, the variability of local controlling factors, and the uncertainty in modelling parameters. A regional assessment must consider a simplified approach that is not highly parameter dependant and that can provide zonation with minimum data requirements. A distributed empirical model has thus been developed for regional susceptibility assessments using essentially a digital elevation model (DEM). The model is called Flow-R for Flow path assessment of gravitational hazards at a Regional scale (available free of charge under www.flow-r.org) and has been successfully applied to different case studies in various countries with variable data quality. It provides a substantial basis for a preliminary susceptibility assessment at a regional scale. The model was also found relevant to assess other natural hazards such as rockfall, snow avalanches and floods. The model allows for automatic source area delineation, given user criteria, and for the assessment of the propagation extent based on various spreading algorithms and simple frictional laws. We developed a new spreading algorithm, an improved version of Holmgren’s direction algorithm, that is less sensitive to small variations of the DEM and that is avoiding over-channelization, and so produces more realistic extents. The choices of the datasets and the algorithms are open to the user, which makes it compliant for various applications and dataset availability. Amongst the possible datasets, the DEM is the only one that is really needed for both the source area delineation and the propagation assessment; its quality is of major importance for the results accuracy. We consider a 10 m DEM resolution as a good compromise between processing time and quality of results. However, valuable results have still been obtained on the basis of lower quality DEMs with 25 m resolution.


Introduction
Landslide susceptibility has been defined by Fell et al. (2008) as "a quantitative or qualitative assessment of the classification, volume (or area), and spatial distribution of landslides which exist or potentially may occur in an area".Debris flows are expected to be more frequent in the most susceptible areas, but as time frame is not taken into account, the susceptibility does not directly represent temporal occurrence probability.Regional susceptibility mapping is used for delineation of potentially threatened areas based on minimum data requirements (Aleotti and Chowdhury, 1999;Glade, 2005;Jakob, 2005;van Westen et al., 2006).Although several models for debris flows propagation assessment exist (O'Brien et al., 1993;Hungr, 1995;Crosta et al., 2003;Iovine et al., 2005;Rickenmann, 2005;Pirulli and Mangeney, 2008;Beguería et al., 2009;Bregoli et al., 2011), very few can be used at the regional scale (Iverson et al, 1998;van Westen et al., 2006;Berti and Simoni, 2007).
Pure mechanical modelling of debris flows is difficult because of the complex nature of the phenomenon (Coussot, 1993;Hungr, 1995;Iverson, 1997), the variability of controlling factors, and the uncertainty in modelling parameters (He et al., 2003).Debris flow models rely on physical variables that cannot be acquired for a whole region at reasonable costs (Carrara et al., 2008).Iverson et al. (1998) suggest the use of simplified spatially distributed models for Published by Copernicus Publications on behalf of the European Geosciences Union.P. Horton et al.: Flow-R, a model for susceptibility mapping of gravitational hazards regional case studies based on empirical or semi-empirical approaches (Costa, 1984;Hungr et al., 1984;Johnson, 1984;Rickenmann, 1999).Empirical parameters can usually be transposed to similar environments.It allows us to calibrate the model for small areas where inventory of past events exists, and to apply the optimal parameters set to a whole region.The weakness of this approach is that site-specific characteristics cannot be accounted for (Hürlimann et al., 2008;Kappes et al., 2011).
The Flow-R model was built to process GIS-based regional susceptibility assessments of debris flows (see also Carrara et al., 1995;Chung and Fabbri, 1999;Hofmeister and Miller, 2003;Melelli and Taramelli, 2004;Guinau et al., 2007).It allows the identification of potential source areas and corresponding propagation extent.A similar approach was used for modelling either the propagation of glacier lake outbursts (Huggel et al., 2003) or debris flows with a procedure called random walk based on multiple modelling of stochastic propagations with limited spreading (Gamma, 2000).Flow-R is available free of charge for Windows and Linux at www.flow-r.org.It is based on Matlab, but is a standalone application thanks to the Matlab Compiler Runtime.
Flow-R has been used in different countries to produce regional debris flow susceptibility maps with satisfying accuracy.It was first used in Switzerland for the Canton de Vaud (Horton et al., 2008) and the Val de Bagnes (Jaboyedoff et al., 2012).Later, it was used also by other universities and national geological services in France (Kappes et al., 2011), Italy (Blahut et al., 2010, Lari et al., 2011a, b), Norway (Fischer et al., 2012), Argentina (Baumann et al., 2011) and Pakistan (Horton et al., 2011).Flow-R is also part of training materials (van Westen et al., 2010;Kappes et al., 2012).
The goal of this paper is to give a full overview of the model and its recent improvements, and to provide relevant examples.The model and its algorithms are explained, along with details of its implementation.Sensitivity analyses of the various parameters are then performed to depict the influence of the DEM resolution, the processing precision, and the parameterization of the propagation.

The model concepts
Flow-R is a spatially distributed empirical model developed under Matlab ® .Application of the model requires two distinctive steps based on a digital elevation model (DEM): (1) the source areas are first identified by means of morphological and user-defined criteria, and then (2) debris flows are propagated from these sources on the basis of frictional laws and flow direction algorithms.The debris flow volume and mass are not taken into account, as accurate values cannot be easily assessed over a large region and due to the significant mass changes occurring through erosion and deposition (Iverson and Denlinger, 2001), which are excessively diffi-cult to estimate.The software has a graphical user interface (Fig. 1) which allows the user to specify criteria for the delineation of the source areas and to choose the algorithms and parameters for the assessment of the propagation.
Flow-R is a tool for susceptibility assessment, but it is not suitable for individual event modelling because the calculated propagation provides a range of possible events.It is, however, advisable to compare the assessed susceptible zone with specific events in order to evaluate the accuracy of the results and to adjust the model parameters.Parameterization is empirical and the user has to specify the parameters according to the hazard type.Different hazards, such as debris flow and mudflow, have to be modelled separately and can later be combined using a GIS software.
The model sets normalized or absolute values for the assessed initial source areas.These values are then propagated downslope and represent a notion of weight in a relative hazard scale, which can be adjusted according to the user's need to account, as an example, for frequency values (see Michoud et al., 2012).When normalized, the values never exceed 1, and thus get close to a notion of spatial probability.For the sake of simplicity, we will further on use the term susceptibility value to denote this weight.

Assessment of the source areas
The source area delineation uses an index-based approach.Input datasets can represent different types of spatial information, and are handled with user-defined parameters.Accordingly, grid cells of each input dataset are classified as (1) favourable, when initiation is possible, (2) excluded when initiation is unlikely, or (3) ignored when no decision can be taken on this parameter (Fig. 2).Datasets are combined according to the following rule: a cell is a source area if it was at least once selected as favourable, but never excluded (Fig. 2).Alternatively, the user can directly import source areas which have been generated by another (GIS-based) approach.
The susceptibility value of the source cells is 1 by default.However, the user can assign a spatial grid of susceptibility or frequency values to characterize the source cells (see Michoud et al., 2012).
According to Rickenmann and Zimmermann (1993) and Takahashi (1981), three criteria in a critical combination are relevant for the initiation of a debris flow: terrain slope, water input, and sediment availability.
Most debris flows occur from terrain with a slope gradient higher than 15 • (Rickenmann and Zimmermann, 1993;Takahashi, 1981).We normally consider this value as the lower initiation threshold.Plan curvature can be used to identify the gullies, provided DEM resolution is sufficient (Delmonaco et al., 2003;Wieczorek et al., 1997), and thus allows refining the delineation of the source areas when the user is interested in debris flows found in channels.
Water input can be represented by the upslope contributing area (flow accumulation) as frequently done in distributed  hydrological models (Tarboton, 1997;Erskine et al., 2006).Literature diverges in regards to the minimal contributing area for initiation.We observed on various case studies that a threshold of 0.01 km 2 seems to be a reasonable value for the central Alps, but Fisher et al. ( 2012) found lower values for Norway.By analysing past events, a relationship was found between the upslope contributing area and the terrain slope.Combining the work of Rickenmann and Zimmermann (1993) on the extraordinary 1987 event in Switzerland, and of Heinimann (1998), we can define two curves (Fig. 3) identifying the lower limit for a debris flow source initiation, one for rare events and the other for extreme events (Horton et al., 2008).The threshold for rare events is given in Eq. ( 1) and the threshold for extreme events in Eq. ( 2):  Heinimann (1998), andRickenmann andZimmermann (1993). tan where tan β thres is the slope threshold, and S uca the surface of the upslope contributing area.These values were established for the central Alps area.Potential users can change these relationships to customize the method to other locations.
Although the plan curvature is often used to identify gullies, there is no recognized threshold in the literature.For a 10 m DEM in western Switzerland, we found an optimum at −2/100 m −1 .Fisher et al. ( 2012) considered values of −1.5/100 m −1 to −0.5/100 m −1 as more accurate for Norway.This threshold is expected to vary with the location and the type of debris flows, and strongly depends on the DEM resolution and accuracy.Gullies may not be identifiable on a Any other dataset can also be added, such as geological maps to account for sediment availability, or land use maps to help removing inaccurate source areas.

Assessment of the propagation
The propagation routine provides a full assessment of possible events in one run starting from the previously assessed source areas (see also Huggel et al., 2003).Two types of algorithms are involved in the propagation assessment: 1. Spreading algorithms controlling the path and the spreading of the debris flows; 2. Friction laws determining the runout distance.

Algorithms for the spreading assessment
The spreading is controlled by flow direction algorithms and persistence functions.

Flow direction algorithms
Several flow direction algorithms are implemented in the software, but not all are relevant for debris flow modelling: -The Holmgren (1994) algorithm adds a parameter to the multiple flow direction algorithm as an exponent x allowing control of the divergence: where i, j are the flow directions, p fd i the susceptibility proportion in direction i, tan β i the slope gradient between the central cell and the cell in direction i, and x the variable exponent.For x = 1 the spreading is similar to the multiple flow direction.When x increases, the divergence is reduced up to resulting into the single flow direction when x → ∞.This parameter gives us control over the spreading and thus allows us to reproduce a wide range of other flow accumulations.On the basis of field and laboratory measurements, Claessens et al. (2005) suggested a value of the exponent equal to 4 for debris flows.
-We developed a modified version of Holmgren's algorithm by changing the height of the central cell by a factor dh, which will change the gradients values.This allows smoothing of DEM roughness and production of more consistent spreading, particularly in the case of high-resolution data.Indeed, a cell that is 20 cm higher than the central pixel will certainly be reached by the real debris flow, whereas this cannot be reproduced by any flow direction algorithm as none integrates a notion of height of the flow.The interest of this new algorithm is illustrated in Sects.4.1.1 and 4.1.3.
-The D8 (O'Callaghan and Mark, 1984;Jensen and Domingue, 1988) algorithm provides convergent but unrealistically straight flow, and too limited spreading (Desmet and Govers, 1996;Endreny and Wood, 2003;Erskine et al., 2006;Tarboton, 1997).It can be, however, useful for a rapid assessment of the parameters related to the runout distance.
-The D∞ algorithm (Tarboton, 1997) allows for small spreading, but is still insufficient in certain conditions such as on the alluvial fans.
-The Rho8 algorithm (Fairfield and Leymarie, 1991) introduces a stochastic approach to improve the paths, but is extremely convergent and induces a lack of deterministic results (Erskine et al., 2006).
-The multiple flow direction approach (Quinn et al., 1991) often results in much too large spreading (Huggel et al., 2003).
-The Freeman (1991) algorithm is similar to the multiple flow direction, with a slightly lower spreading.
For the susceptibility assessment of debris flows, Holmgren's algorithm or its modified version should be chosen for any location as it allows reproducing most of the other flow direction algorithms.Moreover, it is the only algorithm that allows parameterizing of the spreading.

Inertial parameter (persistence function)
The persistence function (Eq.5) aims at reproducing the behaviour of inertia, and weights the flow direction based on the change in direction with respect to the previous direction (see Fig. 4) (Gamma, 2000).
where p p i is the flow proportion in direction i according to the persistence, and α(i) the angle between the previous direction and the direction from the central cell to cell i.
Three implementations of the persistence were chosen (Table 1): the first is called proportional, the second uses a cosine, and the third is based on Gamma (2000).In every persistence distribution, the cell opposed to the flow direction is nulled (w 180 = 0) to avoid eventual backward propagation, and thus to save computing time.

Overall susceptibility
The values given by the flow direction algorithm and the weighting of the persistence are combined according to Eq. ( 6) where i, j are the flow directions, p i is the susceptibility value in direction i, p fd i the flow proportion according to the flow direction algorithm, p p i the flow proportion according to the persistence, and p 0 the previously determined susceptibility value of the central cell.
The result of Eq. ( 6) is a 3 × 3 array with assigned susceptibility values.An additional step is to check the energy values (Sect.2.2.2) and to null the susceptibility of the cells that cannot be reached due to lack of energy.A normalization stage is then required again to bring the sum of the cells to the value of p 0 .This aims at avoiding loss of susceptibility.

Implemented algorithms for the runout distance assessment
The runout distance assessment is based on simple frictional laws; as the source mass is unknown, the energy balance is unitary (Eq.7).The processing takes place at the cell level and controls which other cells the flow would be able to reach.Thus, these algorithms control the runout distance and, in addition, may reduce lateral spreading (when a cell on the border of the spreading cannot be reached because of insufficient energy).
where E i kin is the kinetic energy of the cell in direction i, E 0 kin is the kinetic energy of the central cell, E i pot is the change  Gamma (2000) 1.5 1 1 1 0 in potential energy to the cell in direction i, and E i f is the energy lost in friction to the cell in direction i.
The friction loss can be assessed by two types of algorithms: the two parameters friction model by Perla et al. (1980) and a simplified friction-limited model (SFLM).Both methods can result in similar propagation areas, depending on the parameters choice (Jaboyedoff et al., 2011).

Perla's two parameters friction model
The friction model from Perla et al. (1980) (similar to the one by Voellmy, 1955) was developed for avalanches, but has also been used for debris flows (Zimmermann et al., 1997).It is based on a non-linear friction law, which is the solution of the equation of movement, leading to the velocity V i of the flow at the end of the segment i (Perla et al., 1980): with where µ is the friction parameter, ω is the mass-to-drag ratio, originally expressed as M/D (see Perla et al., 1980), β i is the slope angle of the segment, V 0 is the velocity at the beginning of the segment, L i is the length of the segment, and g the acceleration due to gravity.When the terrain slope decreases rapidly, a correction factor based on the conservation of linear momentum can be applied (Eq.9): This correction factor requires DEM data outside a 3 × 3 array, and thus a larger array must be given to the algorithm.

Simplified friction-limited model
The simplified friction-limited model (SFLM hereafter) is based on the maximum possible runout distance, which is characterized by a minimum travel angle, also named angle of reach (Corominas, 1996) or fahrböschung angle (Heim, 1932).It is the angle of the line connecting the source area to the most distant point reached by the debris flow, along its path: where E f i is the energy lost in friction from the central cell to the cell in direction i, x the increment of horizontal displacement, tan ϕ the gradient of the energy line, and g the acceleration due to gravity.
The minimum angle amongst a set of observations in the Swiss Alps is about 11 • for coarse-and medium-grained debris flows, and 7 • for fine-grained flows (Zimmermann et al., 1997).The value of 11 • is mentioned in various other studies (Rickenmann and Zimmermann, 1993;Bathurst et al., 1997;Huggel et al., 2002).
This approach may result in improbable runout distances in steep catchments due to unrealistic energy amounts reached during the propagation.To keep the energy within reasonable values, a maximum limit can be introduced to ensure not to exceed realistic velocities (Fig. 5).Putting it together with Eqs. ( 7) and ( 10), we can express the velocity as in Eq. ( 11): where h is the difference in elevation between the central cell and the cell in direction i, V max is the given velocity limit.With the observed maximum velocity of debris flows in Switzerland being 13 to 14 m s −1 (Rickenmann and Zimmermann, 1993), we often choose a limit of 15 m s −1 .

Handling the processing precision
A susceptibility threshold has to be applied under which the processing is stopped.It is a necessary limitation otherwise infinitesimal susceptibility values take an important part of the processing time.Thus, the cells with susceptibility below that threshold are nulled and the removed values are redistributed to active cells.Redistribution of the residuals is a step required to ensure conservation of the total susceptibility value.Indeed, if these residual values were only removed, the propagation may stop due to loss of susceptibility values, disregarding its energy.Runout distances may then vary according to the grid size and the spreading algorithm, which is undesirable.
As it is best to ensure local conservation of the susceptibility value, and to avoid a drift of these values to another part of the propagation, the residuals are redistributed locally to directly adjacent cells whenever possible.

Computational approach
The processes involved in the source areas and the propagation assessment are illustrated in Fig. 6.After assessment of the source areas, post-processing may be performed before propagation: -According to the desired type of run (described at the end of the present section), a procedure of identification of the superior source areas may be needed.
-If the option to trigger connected areas is enabled, these have to be identified first.This option is interesting when the connected areas are discretized in a meaningful way (e.g. from a catchment to another), and are not too wide (for performance purposes).It is typically of interest for avalanche modelling.
The propagation routine considers one source area (a cell or a connected group of these) at a time and transfers it into the active cells list (see Fig. 6).From this list the routine selects one cell at a time and spreads it according to the algorithms chosen by the user.It results in new propagation cells that are added to the active cells list.When the list is empty, the current propagation is over and the next source area is selected until all of them are processed.The processing of the propagation is based on 4 levels of matrices, as illustrated in a simplified way in Fig. 7: -The processing level is the small domain (e.g. 3 × 3 pixels) where the propagation is processed for one given cell.Once it has been spread, the new pixels are added to the active cells list.
-The active cells list is an array of cells that are considered as active, meaning their susceptibility and energy are not null, and thus will continue to propagate.Once an active cell has been propagated, its former value is stored in the current propagation level.
-The current propagation level contains the values of the cells that are no longer active.It stores the evolution of the propagation from only one source area.When some resulting cells overlay, which is common, these are merged by saving the sum of susceptibility values and the maximum of energy values.When propagation of one source area is over, it is stored in the results layer.
-Finally, the results layer contains the outcomes of all propagations that are merged by summing or maximizing the susceptibility values, according to the user's choice.The result of the model is the total area that can potentially be reached by debris flows, with an associated susceptibility value.The different overlapping propagations can be combined according to two different approaches: either we sum the susceptibility or we select the maximum value.For susceptibility mapping, the maximum is sufficient as we are mainly interested in the total extent.The sum option may be interesting when the user aims at assessing a notion of probability for analysis purpose, eventually with frequency values assigned to the source areas (see e.g.Michoud et al., 2012).
Three different types of run were implemented that require very different processing time: 1.An overview of the propagation area that is processed by triggering only the sources located at the top of every catchment.This involves a routine to identify the uppermost source cells.
2. A quick assessment that triggers first all the uppermost source cells (as in the overview type), but which also propagates the remaining sources in a subsequent stage.
A test is then performed during the processing to identify previously assessed propagations that resulted in higher energy and susceptibility values.If such a former propagation is found at the same place, the current processing is stopped.
3. The complete option triggers every source area, with no check on previously processed propagations.This option can be chosen in order to control that the quick assessment did not miss any propagation area, which can marginally happen on a few cells due to differences in the propagation direction.Another and more important application is when the user is interested in the sum of the propagations, e.g. when frequency values are assigned to the source areas.
The energy map, which is created by considering the maximum value of overlapping propagations, is a sub-product of the model.It has to be considered with care as the mass of the debris flows is not taken into account.

Case study sites
Three case study sites have been selected for parameters sensitivity analyses.These case studies are located in Switzerland, in the Rhone River catchment (Fig. 8).The source material of debris flows originates either from glacial deposit, colluvium, slope mass movement or/and fluvial deposits.The volumes of significant events are given for documentation, but we would like to remind that these cannot be assessed by Flow-R.

Fully
The Fully case study is a well-documented debris flow that occurred in October 2000 after heavy rainfall on   saturated Quaternary deposits due to hydroelectric pipe leaking (FOWG, 2002).It travelled 3 km downslope and reached the Rhone Valley floor, next to the village of Fully.The event, transporting approximately 350 000 m 3 (Petrascheck, 2003), only affected vineyards and access roads, with no human casualties.The propagation profile is fairly regular, but presents a brutal slope transition to the valley floor, which makes it a non-regular debris fan.

Saint-Barthélémy
Saint-Barthélémy creek is situated 10 km west of Fully.It is historically known for massive debris flows reaching the debris fan and disrupting human activities.The last flows took place between 1926 and 1930, with a total transported volume of about 1 500 000 m 3 , following a rock collapse in the upper part of the catchment (Virieux, 1931).The debris flows, originating from a steep gorge (˜1600 m a.s.l.), repeatedly reached and blocked the Rhone River, 6.5 km downstream of the initiation zone.The longitudinal profile of the creek is more complex than in the Fully case study, with changes in slope at multiple locations between the initiation zone and the debris fan's apex.However, the slope angle of the debris fan is very constant, at 5 • .

Solalex
Solalex is a small village located south of the Diablerets Range, at 1470 m a.s.l.Debris flow events are very common on the Diablerets Range south side, regularly blocking a private road to access the Anzeinde pastures, located 400 m higher.We will focus on a fan on which recent debris flows follow a complex path, on the very edge of the fan (Fig. 9).
These fans are very active because the upper part is made of folded limestone with included marl layers belonging to the Diablerets nappes (Badoux and Gabus, 1990).This creates steep, small, impervious catchments highly productive in rock fragments, having a very short time of concentration.

Parameters sensitivity analyses
The parameters sensitivity analyses illustrate the effect of some choices in the algorithms and parameters.One has to keep in mind that even though limited interactions may exist between the spreading algorithms and the runout distance parameters, they can be considered independently during calibration of the model.As one can see on Fig. 10, the propagation extent is satisfying at a resolution of 10 m for the results with the standard version of Holmgren's algorithm (dh = 0).Finer resolutions are missing significant areas due to the roughness of the DEM and the strong effect of channelization at such scale, while rougher resolutions artificially enlarge the extent.It is consistent with the observations made by Zhang and Montgomery (1994), who concluded that "for many landscapes, 10 m grid size presents a rational compromise between increasing resolution and data volume for simulating geomorphic and hydrological processes".Quinn et al. (1995) also recommend a 10 m resolution for flow direction algorithms to capture the variability of the topographic form for hillslopes, and point out that lower cell sizes do not bring significant information.
The chosen case study in Solalex presents a main channel that is located on the west edge of the fan.One can see in Fig. 10 that most of the susceptibility value follows this channel for the 2 m and 5 m grid, whereas from a 10 m resolution, the susceptibility value distribution begins to shift eastward.This involves that, for people interested in the pattern of the susceptibility values, a higher resolution may be needed.However, as it was said before, the standard Holmgren's algorithm results in an unsatisfying extent with many omissions.This is where our modified version shows its biggest relevance -by allowing the spreading to be guided by the general topography, and not by the DEM details.One can see that the extent is significantly more consistent in between the 2 m, 5 m, and 10 m resolutions with a central cell elevated by 2 m in the modified Holmgren's algorithm (Fig. 10).For the 25 m and 50 m DEM, the effect is negligible.

Susceptibility threshold
Figure 11 illustrates the impact of the choice of the susceptibility threshold on propagations with similar parameters on two DEM resolutions (2 m and 10 m).The first threshold (10 −3 ) does not produce satisfying extents for both resolutions.On the 2 m resolution, many small straight patterns (oriented east-west and north-south) appear.The next two thresholds (10 −4 and 10 −5 ) show pretty similar extents, particularly at 10 m resolution.
The 2 m DEM is very sensitive to the susceptibility threshold, which should be set to 10 −4 .Indeed, as the processing time exponentially increases with the decreasing of the threshold, the value of 10 −4 is considered as a good compromise between accuracy and required calculation resources.However, processing time at that resolution remains very important (approximately 50 times more than on the 10 m DEM).
For the 10 m resolution, very similar results were found for a precision of 10 −4 or above, up to approximately 5 × 10 −4 .We would recommend a threshold of 3 × 10 −4 for that resolution.

Holmgren's exponent
The exponent in the classic Holmgren's algorithm has a strong effect on the spreading on a 2 m resolution grid, as one can see on Fig. 12.As the exponent increases, the flow quickly converges.The impact on the extent here is to be interpreted by the use of a susceptibility threshold as all extents should look similar in case of an infinite precision.The propagation extent on a 10 m resolution DEM is less sensitive to this parameter.
Figure 12 also depicts the extent of the modified Holmgren's version with exponents of 4 and 6, which is a commonly chosen range.This version does not present the very strong channelization related to the DEM fine-scale features.

Sensitivity analysis
A parametric analysis of both the Perla and the SFLM models was performed on the Fully and Saint-Barthélémy case studies at a 10 ṁ resolution, and on an artificial topography (smooth polynomial longitudinal profile, confined spreading).Longitudinal profiles are illustrated in Fig. 13.Parameter sets were processed by simultaneously varying µ and ω in Perla's model and tan ϕ and the velocity limit in the SFLM model.The compared output is the runout distance according to the most probable path (identified by means of a D8 algorithm).The parameters controlling the spreading are fixed and have no influence on the runout distance.
Figure 14 shows that the results of Perla's model applied to the virtual topography are affected mainly by the µ coefficient (expressed as ϕ in the figure).The ω parameter plays an increasing role for high values of µ.The SFLM model is shown to be more sensitive to the friction angle than the velocity limit for parameters in the usual range.These observations can also be found in the real-life case studies, but with significant influence of the actual topography.The Perla's model runout on Fully and Saint-Barthélémy is affected by both parameters for values of ω inferior to 40 or 50.Above that limit, the runout distance seems also to be controlled by the µ coefficient.For the SFLM model, we thus recommend to set the velocity limit to accurate values and to calibrate the friction angle.

Parameterization using expert knowledge
It is not trivial to choose the parameters of both the Perla and the SFLM models.Nevertheless, some simple considerations may be used to define the range of values for the friction coefficient (µ), the friction angle (ϕ), and the mass-to-drag ratio (ω).Experts have the capacity to get information from a simple field trip by observing the runout distance of debris flows, the topography, and the type of material.In addition, videos of events are more and more available.Some simple considerations presented hereafter show how estimates of the parameters can be performed along profiles.
In the case of the SFLM model applied to a debris flow that reaches a fan and stops on its slope (Fig. 15), it is possible to estimate the friction coefficient µ starting from the relationship (Jaboyedoff et al., 2011): where x is the horizontal travel distance over the fan, β f the fan slope angle and ϕ the friction angle.Based on direct observations or by visiting the apex of the fan, experts can assess the velocity v a at the apex, which permits approximation of the value of the energy line at the apex, using v a = √ 2g h el .Knowing x and v a , an estimate of the friction angle can be performed using (Jaboyedoff et al., 2011) Equation ( 13) lets us assume that the maximum velocity for the slope is reached at the apex.Thus, the mass-to-drag ratio (ω) can be assessed from the maximum velocity for an infinite slope in Perla's model (Perla et al., 1980;Jaboyedoff et al., 2011): This permits us to get a first guess of the parameters.It has been shown for Perla's model that the results of the simulations are very close to the SFLM model when the velocity limit of the latter is equal to the velocity for an infinite slope in Perla's model, using the average slope (β a ) of the studied profile (Jaboyedoff et al., 2011).As a consequence, using the horizontal travel distance L, the total time of travel t can be obtained from These equations permit evaluation of whether the chosen parameters are relevant.
5 Published case studies Jaboyedoff et al. (2012) used Flow-R to create, among other studies, a susceptibility map for debris flows at a scale of 1 : 25 000 (see also Horton et al., 2011).The model was calibrated on four test sites within the catchment.The fieldwork showed that most of the estimated source areas (86%) were relevant.Horton et al. (2009Horton et al. ( , 2011) illustrated a successful application of the model despite the lack of data and the coarse resolution of the DEM in Pakistan.It is also, as far as the authors know, the first regional susceptibility mapping of hyperconcentrated flows.Blahut et al. (2010) proceeded to a hazard zonation of debris flows in part of the Valtellina Valley (Lombardy Region, Italian central Alps), which "provided results for consequent hazard and risk studies".The specificity of the study is the use of spatial probability and of a probability of occurrence (hazard frequency) for classification of the source areas.Lari et al. (2011a) used Flow-R in the context of a study on societal and economic quantitative risk assessment related to rockfall and debris flow events in lower Valtellina (Lombardy Region, Italy).Lari et al. (2011b) went further in this study and provided more details of the method to assess risk based on vulnerability estimation and source frequency spatialization with Flow-R.Two sets of parameters were calibrated for open-slope or channelled debris flows, and the model was found to "allow a good modelling of both".Kappes et al. (2011Kappes et al. ( , 2012) ) assessed the model in the Barcelonnette Basin (France) for three classes of events: high, medium, and low frequency events.The study concludes that "a comparison with the footprints of a few mapped events indicates reasonable results but suggests a high dependency on the quality of the digital elevation model".The results accuracy was assessed by means of methods proposed by Beguería (2006).Fischer et al. (2012) presented a case study that aimed at assessing the feasibility of a national debris flows susceptibility map for Norway.The model was found to be very robust in identifying the source areas in channelled topography, whereas the results for open-slope topography are sensitive to fine parameterization.Different sub-regions homogenous in terms of geological settings were identified to allow for distinct parameterization.Finally, after propagation modelling, the authors concluded that "good correlation between the modelling results and the field observations was found".

Application of Flow-R to other phenomena
The model contains various algorithms that were found to be also accurate for other natural hazards such as rockfall, snow avalanche and flooding.A successful application to rockfall is presented in Michoud et al. (2012) by using the SFLM approach.This algorithm has similarities with the CONEFALL method (Jaboyedoff and Labiouse, 2011), and resulted in quite similar zonation, sometimes even more realistic due to better integration of the topography.Block release frequencies were first assessed and then propagated with Flow-R to result, after calibration, in a regional hazard map for rockfall.Values of frequency of block releases from the source areas were thus integrated within the propagation calculation.
Flood assessment was conducted by Jaboyedoff et al. (2010) for the Bagnes Valley (Switzerland).Flow-R is not able to determine the location of an overflow of river banks.This part was assessed in the field, and factors favouring flooding, such as bridges or low banks, were noted on a map.Hereafter, Flow-R was used to simulate the flow of water outside the river.The goal here was not to model a frictional behaviour, but to spread water on the topography as a static flood inundation.
Finally, snow avalanches were also assessed for the Bagnes Valley (Switzerland) (Jaboyedoff et al., 2010(Jaboyedoff et al., , 2012)).The sources were estimated by means of geomorphological indices.Perla's friction model was obviously chosen to process the propagation.As this algorithm was first developed for the modelling of snow avalanche propagation, it is no surprise that the model outputs fit accurately the inventory of past events (results are shown for the Verbier area in Fig. 16).Avalanche susceptibility mapping has also been performed for both studied districts in Pakistan (Horton et al., 2009).

Discussion and conclusions
The Flow-R model relies on an empirically distributed approach for regional susceptibility mapping, and thus cannot integrate local controlling factors and actual physical behaviours.However, numerous successful case studies demonstrated its suitability for debris flow susceptibility mapping.One of its main advantages is its low data requirement; a DEM only may be sufficient to assess potential source areas -if not known a priori -and to process the propagation.Another highlight is its opening to the user in terms of inputs and algorithms, which makes it versatile for very different case studies and even hazard types.
Resulting propagation areas from the model are in general larger than observed events in the field, which is intentional in the framework of susceptibility mapping.Indeed, the map should contain every possible event, even the worst case.Susceptibility maps offer a good overview to point out where field investigation should be conducted to establish a detailed hazard map.
Amongst the data input the DEM is the most important one.It allows source areas delineation by deriving morphological criteria that were shown to be useful, and is the basis of the propagation.Both the source areas identification and the propagation strongly depend on the DEM quality.Its resolution and accuracy are key elements which condition the quality of outputs.Artefacts, such as smoothed gullies, important errors or a hidden flow, e.g.under a bridge or in deep gorges, can be misleading.However, small irregularities do not affect the outputs if the algorithms and parameters are chosen accordingly.If modelling with very coarse DEMs is performed, the user has to keep in mind that results are likely to represent major torrents, but it is less probable that every small debris flow in minor gullies will be shown.
The analysis of the parameters influencing the propagation may help Flow-R users to choose their data and algorithms according to their needs.These choices mainly depend on the study objectives.Over-precision is time consuming, as processing time increases exponentially with the increasing DEM resolution.If one aims at creating a susceptibility map over a large region, a 10 m resolution is sufficient.Fischer et al. (2012) also noticed that a DEM with a 10 m resolution shows good agreement between the assessed source areas and the known flow tracks.A 25 m resolution produces results of lower quality, but still usable, whereas the 50 m resolution is too coarse and the resulting extent should be interpreted with care.If one aims at studying some torrents, a finer resolution can be recommended, but with the use of the modified version of Holmgren's algorithms exclusively.
The modified version of Holmgren's algorithm improves the spreading extent by making it less sensitive to DEM small features, and so less dependent on the DEM resolution.It provides a coverage that is more realistic, and allows for more accurate spreading in flat areas.We recommend using it instead of the classic algorithms for any type of hazard modelling, with a dh consistent with the studied phenomena.
The analysis of the parameters acting on the runout distance gives a useful overview of their respective influence in both the Perla and the SFLM models.The user can thus have an idea of the governing parameters and their sensitivity in the framework of real case studies.For Perla's model, the linear friction parameter is globally more controlling the runout distance than ω.For the SFLM the velocity limit is not a factor on which the user should play, but it should be fixed to realistic values according to observations of the modelled hazard.Jaboyedoff et al. (2011) show that both the Perla and the SFLM models reproduce well the non-linearity of the relationship between the velocity and the terrain slope and that they can be used both in a similar way.The SFLM model presents the advantage that parameters like the velocity limit and the angle of reach can be estimated through simple observation of debris flow events.
Flow-R is available to download for free at www.flow-r.org.As the software utility has been proven and users' interest is growing, a new version will be developed in the near future.It will be written in C++ for performance purposes, will be cross platform, and will enclose a GIS engine.

Fig. 1 .
Fig. 1.Screen capture of the main frame of the Flow-R software.

Fig. 2 .
Fig. 2. Illustration of the combination of various datasets for the assessment of the source areas.

P.
Horton et al.: Flow-R, a model for susceptibility mapping of gravitational hazards rough DEM, and thus curvature may become useless, if not misleading.

Fig. 4 .
Fig. 4. Illustration of the spreading of susceptibility value to the neighbouring cells.

Fig. 5 .
Fig. 5. Illustration of the travel angle and the velocity limitation of the simplified friction-limited model (SFLM).

Fig. 7 .
Fig. 7. Illustration of the 4 different layers of processing of the propagation and their associated data flow.

Fig. 8 .
Fig. 8. Location of the case studies in the Swiss Alps.The catchments are displayed in red for Saint-Barthélémy and Fully, along with the location of the longitudinal profiles (yellow) given in Fig. 13.The extent of the Solalex case study (Fig.9) and the Verbier area (Fig.16) are depicted in blue (Geodata ©swisstopo -DV084371).

Fig. 9 .
Fig. 9.A fan on the Diablerets Range near Solalex (see location in Fig.8).The red shape shows the extent of the whole susceptibility zone modeled with Flow-R from a chosen source.Propagation parameters are the same as the case with precision 10 −4 and 10 m DEM in Fig.11(Geodata ©swisstopo -DV084371).
Propagation was processed for the Solalex case study (Sect.3.3) by means of the modified version of Holmgren's algorithm.The DEM, originally with a 1 m resolution, was degraded to 2 m, 5 m, 10 m, 25 m, and 50 m to illustrate the effect of the data resolution and the role of the modification in Holmgren's algorithm.

Fig. 10 .
Fig. 10.Effect of the DEM resolution and the modification of Holmgren's spreading algorithm on a debris flow fan in Solalex, Switzerland (Geodata ©swisstopo -DV084371).

Fig. 13 .
Fig. 13.Longitudinal profiles of the three case studies (see location in Fig. 8).

Fig. 14 .
Fig. 14.Effect of the parameters of both algorithms for energy assessment on the runout distance for the virtual topography, Fully and Saint-Barthélémy case studies, Switzerland.

Table 1 .
Implemented weightings of the persistence function in the assessment of the spreading.