r.randomwalk v1, a multi-functional conceptual tool for mass movement routing

Abstract

We introduce r.randomwalk, a flexible and multi-functional open-source tool for backward and forward analyses of mass movement propagation. r.randomwalk builds on GRASS GIS (Geographic Resources Analysis Support System – Geographic Information System), the R software for statistical computing and the programming languages Python and C. Using constrained random walks, mass points are routed from defined release pixels of one to many mass movements through a digital elevation model until a defined break criterion is reached. Compared to existing tools, the major innovative features of r.randomwalk are (i) multiple break criteria can be combined to compute an impact indicator score; (ii) the uncertainties of break criteria can be included by performing multiple parallel computations with randomized parameter sets, resulting in an impact indicator index in the range 0–1; (iii) built-in functions for validation and visualization of the results are provided; (iv) observed landslides can be back analysed to derive the density distribution of the observed angles of reach. This distribution can be employed to compute impact probabilities for each pixel. Further, impact indicator scores and probabilities can be combined with release indicator scores or probabilities, and with exposure indicator scores. We demonstrate the key functionalities of r.randomwalk for (i) a single event, the Acheron rock avalanche in New Zealand; (ii) landslides in a 61.5 km${}^{\mathrm{2}}$ study area in the Kao Ping Watershed, Taiwan; and (iii) lake outburst floods in a 2106 km${}^{\mathrm{2}}$ area in the Gunt Valley, Tajikistan.

Introduction

Mass movement processes such as landslides, debris flows, rock avalanches, or snow avalanches may lead to damages or even disasters when interacting with society. Computer models predicting travel distances, hazardous areas, impact energies, or travel times may help the society to mitigate the effects of such processes and, consequently, to reduce the risk and the losses (Hungr et al., 2005).

Physically based dynamic models are used for in-detailed analyses of specific events or situations (e.g., Savage and Hutter, 1989; Takahashi et al., 1992; Iverson, 1997; Pudasaini and Hutter, 2007; McDougall and Hungr, 2004, 2005; Pitman and Le, 2005; Christen et al., 2010a, b; Mergili et al., 2012b; Pudasaini, 2012; Hergarten and Robl, 2015; Mergili et al., 2015). Since the processes are complex in detail and the input parameters are uncertain, simplified conceptual models for the motion of mass flows are today used in combination with GIS (Geographic Information System). These models may be used for single events. However, they are particularly useful to indicate potential impact areas at broader scales. Hypothetic mass points are routed from a release pixel through a digital elevation model (DEM) until a defined break criterion is reached. Monte Carlo techniques (random walks, Pearson, 1905; Gamma, 2000) or multiple flow direction algorithms (Horton et al., 2013) are employed to simulate the lateral spreading of the flow.

The break criteria often consist in threshold values of the angle of reach (i.e., the average slope of the path) or horizontal and vertical distances (Lied and Bakkehøi, 1980; Vandre, 1985; McClung and Lied, 1987; Burton and Bathurst, 1998; Corominas et al., 2003; Haeberli, 1983; Zimmermann et al., 1997; Huggel et al., 2002, 2003, 2004a, b), sometimes related to volume (Rickenmann, 1999; Scheidl and Rickenmann, 2010). However, those relationships usually display a large degree of scatter. Further, key parameters for design issues, such as impact pressures, are not provided (Hungr et al., 2005).

Some approaches include simplified physically based models going back to the mass flow model of Voellmy (1955), relating the shear traction to the square of the velocity and assuming an additional Coulomb friction effect (Pudasaini and Hutter, 2007). They consider only the centre of the flowing mass, but not its deformation and the spatial distribution of the flow variables. This type of models is mainly used for snow avalanches and debris flows (Perla et al., 1980; Gamma, 2000; Wichmann and Becht, 2003; Mergili et al., 2012a; Horton et al., 2013).

Various – mostly open-source – software tools for conceptual modelling of mass movements (mainly flows) at medium or broad scales are available (e.g., Gamma, 2000; Wichmann and Becht, 2003; Mergili et al., 2012a; Horton et al., 2013). However, most of these tools lack substantial features: (i) they are limited to one single type of break criterion; (ii) they do not allow one to directly account for the uncertainty of the break criteria; (iii) they do not allow one to back calculate the statistics of a set of observed mass movements; and (iv) they do not offer built-in functionalities for evaluating the model results against observations. Consequently, the key objectives of the present study are

• to introduce r.randomwalk, a freely available, comprehensive and flexible tool for routing mass movements;

• to demonstrate the various functionalities of r.randomwalk, particularly in terms of overcoming the issues (i)–(iv);

• to discuss the potentials and limitations of this tool.

Next, we will describe the r.randomwalk software tool (Sect. 2). Furthermore, we will present the test areas and the results (Sect. 3). Finally, we will discuss the findings (Sect. 4) and conclude with some key messages of the work (Sect. 5).

The r.randomwalk application

Computational implementation

r.randomwalk is implemented as a raster module of the open-source software package GRASS (Geographic Resources Analysis Support System) GIS 7 (Neteler and Mitasova, 2007; GRASS Development Team, 2015). We use the Python programming language for data management, pre-processing and post-processing tasks (module r.randomwalk). The routing procedure (see Sect. 2.2–2.4) is written in the C programming language (sub-module r.randomwalk.main). The R software environment for statistical computing and graphics (R Core Team, 2015) is employed for built-in validation and visualization functions (see Sect. 2.5). Parallelization of multiple model runs is enabled. It allows for the exploitation of all computational cores available, speeding up analysis processes. The parallelization procedure is implemented at the Python level (analogous to the way described in Mergili et al., 2014): the module r.randomwalk produces a batch file for each model run. This batch file calls the Python-based sub-module r.randomwalk.mult, which is then used to launch r.randomwalk.main with the specific parameters for the associated model run. Thereby, the Python library “Threading”, a higher-level threading interface, and the Python module “Queue”, a class helping to block execution until all the items in the queue have been processed, are exploited. Parallel processing serves for reducing the computational time in the following contexts:

• Analyses with multiple random subsets of the release areas or coordinates. In each model run, one subset is used for back calculating the probability density function (PDF) of the angle of reach, the other subset is employed for validating the distribution of the impact probability derived with this PDF against the observed deposition areas.

• Analyses with multiple combinations of input parameters varied in a controlled or randomized way, enabling one to consider parameter uncertainties and to explore parameter sensitivity.

r.randomwalk was developed and tested with Ubuntu 12.04 LTS and is expected to also work on other UNIX systems. A simple user interface is available. However, the tool may be started more efficiently through command line parameters, enabling a straightforward batching on the shell script level. This feature facilitates model testing, the combination with other GRASS GIS modules and the consideration of process chains (i.e., using the output of one analysis as the input for the next one). The logical framework is illustrated in Fig. 1, the key variables used in r.randomwalk are summarized in Table 1.

Logical framework of r.randomwalk. Only those components covered in the present article are shown.

[Figure omitted. See PDF]

All tests (see Sect. 3) are performed on an Intel^® Core i7 975 with 3.33 GHz and 16 GB RAM (DDR3, PC3-1333 MHz), exploring a maximum of eight cores through hyperthreading.

Random walk routing

The term random walk refers to a Monte Carlo approach for routing an object through any type of space. The term was introduced by Pearson (1905). Constrained random walk approaches are used for routing mass movements such as debris flows through elevation maps (DEMs), e.g. by Gamma (2000), Wichmann and Becht (2003), Mergili et al. (2012a), and Gruber and Mergili (2013). Such methods enable a certain degree of spreading of the movement by also considering other routing directions than the steepest descent. It avoids the concentration of flows – or any other types of mass movements – to linear features, which would not be realistic for debris flows, snow avalanches, or other types of mass movements. However, the routing is constrained or weighted by factors such as the slope or the perpetuation of the flow direction. An alternative to constrained random walk routing would consist in a multiple flow direction algorithm (Horton et al., 2013).

Control length $L_{\mathrm{ctrl}}$ and segment length $L_{\mathrm{seg}}$. (a) Application of $L_{\mathrm{ctrl}}$ to avoid sharp bending of the flow. (b) Smoothing of the flow path by introducing segments with maximum length of $L_{\mathrm{seg}}$.

[Figure omitted. See PDF]

In the context of r.randomwalk, each random walk routes a hypothetic mass point from a release pixel through the DEM until a break criterion is reached (see Sect. 2.3). A large set of random walks is required for each mass point in order to achieve a satisfactory cover of the possible impact area. r.randomwalk is designed for

• one set of random walks for one mass movement, starting from a defined set of coordinates;

• multiple sets of random walks for one mass movement, one set starting from each pixel of the release area;

• sets of random walks for multiple mass movements in a study area (either starting from one set of coordinates per mass movement, or from all pixels defined as release areas);

• one set of random walks starting from each pixel in the study area.

Overlay rules for different random walks and sets of random walks are applied (see Sect. 2.4).

Summary of the key variables used in r.randomwalk.

Symbol	Unit	Name	Description/remarks
	[range]
Input
$n_{\mathrm{walks}}$	$-$[$\ge$ 0]	Size of a set of walks	Logarithm (base of 10) of the number of random walks
$L_{\mathrm{ctrl}}$	m [> 0]	Control length	Backward distance of each step of a random walk over which the horizontal distance of motion has to increase (see Fig. 2a)
$L_{\mathrm{seg}}$	m [> 0]	Segment length	Length of segments used for computing $L_{max}$ (see Fig. 2b)
$R_{max}$	m [$\ge$ 0]	Maximum vertical run-up height	–
$f_{\mathrm{d}}$	$-$[$\ge$ 1]	Direction factor	Factor for weighting the perpetuation of the movement direction during routing
$f_{\mathit{\beta }}$	$-$[$\ge$ 0]	Slope factor	Factor for weighting the slope during routing
$V$	m${}^{\mathrm{3}}$ [> 0]	Movement volume	–
$Q_{\mathrm{p}}$	m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ [> 0]	Peak discharge	Applicable to lake outburst events (see Table 3)
$a,b,c$	– [ ]	–	Parameters needed by the rules and relationships applied (see Table 3)
$P_{\mathrm{R}}$	$-$[0–1]	Release probability	Spatial probability that a mass movement is released from a given pixel (Mergili and Chu, 2015)
RIS	$-$[$\ge$ 0]	Release indicator score	Ordinal score denoting the tendency of a pixel to produce a mass movement
EIS	$-$[$\ge$ 0]	Exposure indicator score	Ordinal score denoting the exposed values at a given pixel
Output
$L_{max}$	m [$\ge$ 0]	Travel distance	Horizontal distance between the release pixel and the most distant pixel reached by a set of random walks, measured along the segments of the path (see Fig. 2b)
$Z$	m [$\ge$ 0]	Elevation loss	Vertical distance between the release pixel and the most distant pixel reached by a set of random walks
${\mathit{\omega }}_T$	${}^{\circ }$ [< |90|]	Angle of reach	Average slope angle measured between the release pixel and the most distant pixel reached by a set of random walks
IF	$-$[$\ge$ 0]	Impact frequency	Number of random walks impacting a given pixel
IIS	$-$[$\ge$ 0]	Impact indicator score	Number of rules and relationships predicting an impact on a given pixel
IHIS	$-$[$\ge$ 0]	Impact hazard indicator score	Ordinal score serving as a qualitative surrogate for the hazard of an impact on a given pixel
III	$-$[0–1]	Impact indicator index	Fraction of model runs impacting a given pixel out of all model runs
$P_{\mathrm{I}}$	$-$[0–1]	Impact probability	Spatial probability that a given pixel is impacted, building on user-defined release areas and a cumulative density function
$P_{\mathrm{I},\mathrm{C}}$	$-$[0–1]	Composite probability	Spatial probability that a given pixel is impacted, building on $P_{\mathrm{R}}$ and $P_{\mathrm{I}}$ (Mergili and Chu, 2015)
IRIS	$-$[$\ge$ 0]	Impact risk indicator score	Ordinal score denoting the expected/potential loss at a given pixel

During the pixel-to-pixel routing procedure, turns of > 90${}^{\circ }$ are not supported. Neighbour pixels are further considered invalid as target pixels in case they are out of the study area or conflict with at least one of the following limitations:

• In order to constrain upward movements, a user-defined maximum vertical run-up height $R_{max}$ is introduced. It takes the lowest elevation the random walk has passed through as reference.

• Certain types of mass flows (i.e., those with high viscosity) hardly change their flow direction sharply. The user-defined horizontal control distance $L_{\mathrm{ctrl}}$ defines the backward distance of each step over which the horizontal distance of motion has to increase (Fig. 2a).

Possibilities to define the break criteria. The flags provided through the command line or the user interface define the type of break criterion. RC is release coordinates (release from highest points of release areas), RP is release pixels (release from all pixels within release areas), $•$ is relevant for most applications, and $\circ$ is relevant for some applications.

Flag	Mode	Release		Output	Validation	Multiple cores
		RC	RP
q	Published relationships (see Table 3)	$•$		IIS	$\circ$
m	Published relationship, multiple runs	$•$		III	$\circ$	$•$
p	Impact probability	$\circ$	$\circ$	$P_{\mathrm{I}}$	$\circ$
p $+$ n	Impact probability, multiple runs	$\circ$	$\circ$	$P_{\mathrm{I}}$	$•$	$•$
b	Reconstruction of observation	$\circ$	$\circ$	CDF

The probability $P_{\mathrm{px}}$ of any other neighbour pixel px to become the target pixel is $$P_{\mathrm{px}}={\displaystyle \frac{p_{\mathrm{px}}}{\sum _{qx=\mathrm{1}}^{qx=n}{p}_{qx}}},p=f_{\mathrm{d}}{e}^{f_{\mathit{\beta }}\tan \mathit{\beta }},$$ where $n$ is the total number of valid neighbour pixels, and $\mathit{\beta }$ is the local slope between the current pixel and the considered neighbour pixel. $f_{\mathrm{d}}$ and $f_{\mathit{\beta }}$ are weighting factors for the perpetuation of the flow direction and for the slope. $f_{\mathrm{d}}$ is governed by the input parameter $d$: $f_{\mathrm{d}}=d^{\mathrm{2}}$ for the same flow direction as the previous one, $f_{\mathrm{d}}=d$ for a 45${}^{\circ }$ turn and $f_{\mathrm{d}}=\mathrm{1}$ for a 90${}^{\circ }$ turn.

The break criteria for the random walks (see Sect. 2.3) are directly or indirectly related to the travel distance $L_{max}$ i.e., the horizontal length between the release pixel and the terminal pixel measured along the flow path. Preliminary tests reveal that random walk routing through raster maps may result in quite uneven flow paths (see Fig. 2b). Consequently, the distance calculated by summing up all the pixel-to-pixel distances may be significantly longer than the more relevant distance along the observed main flow paths. Employing the sums of the pixel-to-pixel distances would lead to an underestimation of the angle of reach and, consequently, of the predicted travel distances and impact areas. We approach this problem by dividing the flow paths into straight segments with a user-defined maximum length of $L_{\mathrm{seg}}$. The travel distance $L_{max}$ is defined as the sum of the length of all segments (see Fig. 2b). Larger values of $L_{\mathrm{seg}}$ are expected to result in shorter travel distances due to the more pronounced smoothing of the path.

Break criteria

Each random walk continues until at least one neighbour pixel is outside the study area, or until the user-defined break criterion is fulfilled. The break criteria are the key parameters for estimating the mass flow impact areas and can be defined in various ways (Table 2):

• The angle of reach ${\mathit{\omega }}_T$ or the maximum travel distance $L_{max}$ is computed from empirical–statistical rules or relationships, based on the analysis of observed events (Table 3). They usually refer to the distance between the highest spot of the release area and the most distant spot of the impact area along the flow path (the Fahrböschung according to Heim, 1932). Consequently, random walks using this type of break criterion have to start from the set of coordinates defining the highest point of the observed or expected mass movement. Alternatively, also a semi-deterministic model (Perla et al., 1980) can be used.

• Empirical–statistical relationships or the semi-deterministic model may be applied in a large number of parallel computations with randomized values of the parameters $a$, $b$ and $c$ (see Fig. 1 and Table 3). This allows one to explore the effects of uncertainties in the relationships. Only one type of relationship is considered at once, and the output consists in a raster map of the impact indicator index III in the range 0–1, representing the fraction of tested parameter combinations predicting an impact on the pixel (i.e., where impact indicator score (IIS) $=$ 1). Further, the results of all model runs are stored in a way ready to be analysed with the parameter sensitivity and optimization tool AIMEC (Automated Indicator-based Model Evaluation and Comparison; Fischer, 2013).

• An impact probability raster map $P_{\mathrm{I}}$ in the range 0–1 is computed from a user-defined sample of observed values of tan(${\mathit{\omega }}_T$), which is employed to build a cumulative density function (CDF). The CDF represents the probability that the movement reaches the pixel associated with each value of tan(${\mathit{\omega }}_T$). The sample of observed values may be divided into one subset of mass movements for building the CDF, and another one for computing $P_{\mathrm{I}}$. This ensures a clear separation between parameter optimization and model validation (see Sect. 2.5). Parallel processing may be used to repeat the analysis for many random subsets in order to achieve a more robust result.

• If an inventory of events is available, the observed impact areas may be back calculated by routing each random walk until it leaves the observed impact area of the corresponding mass movement. This mode can be used to explore the statistical distribution of ${\mathit{\omega }}_T$. The resulting CDF can be used as input to estimate $P_{\mathrm{I}}$.

Types of rules and relationships supported by r.randomwalk. ${\mathit{\omega }}_T$ is angle of reach, $L_{max}$ is travel distance, $V$ is volume of motion, $Z$ is elevation loss, $Q_{\mathrm{p}}$ is peak discharge at release, and $v_{\mathrm{T}}$ is velocity at termination.

ID	Equation		Examples
			Reference	Process	$a$	$b$	$c$
1	${\mathit{\omega }}_T=a$	(2)	Haeberli (1983);	Debris flow from GLOF	11
			Huggel et al. (2002)
2	${\log }_{10}\tan {\mathit{\omega }}_T=a{\log }_{10}V+b$	(3)	Scheidegger (1973)	Rock avalanche	$-$0.15666	0.62419
3	$L_{\max }=a{V}^b{Z}^c$	(4)	Rickenmann (1999)	Debris flow	1.9	0.16	0.83
4	${\mathit{\omega }}_T=a{Q}_{\mathrm{p}}^b$	(5)	Huggel (2004)	GLOF	18	$-$0.07
5	$v_{\mathrm{T}}=\mathrm{0}$	(6)	Perla et al. (1980)

Overlay of random walk results

The overlay of individual random walks operates at two levels:

• Random walks of the same mass point: impact frequency (IF) is increased by 1 for each random walk predicting an impact. IIS is increased by 1 for each model where at least 1 random walk predicts an impact. The average angle of path – and therefore also $P_{\mathrm{I}}$ – is derived from the random walk with the shortest travel distance (i.e., the straightest flow path and the highest value of $\mathit{\omega }$) at the considered pixel.

• Sets of random walks for different mass points: the values of IF for all random walks impacting a pixel are just added up whilst the maximum of IIS is applied to each pixel. The issue gets more complex when it comes to $P_{\mathrm{I}}$: depending on the specific application, the maximum or the average out of all sets of random walks is more appropriate.

The resulting maps of $P_{\mathrm{I}}$ or IIS can be automatically overlaid with a release probability ($P_{\mathrm{R}}$; result: composite probability $P_{\mathrm{I},\mathrm{C}}$; Mergili and Chu, 2015) or a release indicator score (RIS; result: impact hazard indicator score – IHIS), and with an exposure indicator score (EIS) derived from the land cover (result: impact risk indicator score – IRIS; see Table 1). These steps are not further considered in this article and are therefore not shown in Fig. 1.

Validation

r.randomwalk includes three possibilities for validation of the model results. All three build on the availability of a raster map of the observed deposition area of the mass movement(s) under investigation. All parts of the observed impact areas outside of the observed deposition areas are set to no data (Fig. 3).

Model validation with an ROC plot, relating the false positive rate $r_{\mathrm{FP}}$ and the true positive rate $r_{\mathrm{TP}}$. This way of validation is suitable for predictor raster maps in the range 0–1, such as III or $P_{\mathrm{I}}$. It can also be used for binary predictor maps (0 or 1). In such a case AUC${}_{\mathrm{ROC}}$ is computed from two threshold levels only.

[Figure omitted. See PDF]

• For IIS, the true positive (TP), true negative (TN), false positive (FP), and false negative (FN) predictions are counted on the basis of pixels and put in relation. All pixels with IIS $\ge$ 1 are considered as observed positives (OP); all pixels with IIS $=$ 0 are considered as observed negatives (ON).

• ROC (receiver operating characteristics) plots are produced for III or $P_{\mathrm{I}}$: the true positive rate $r_{\mathrm{TP}}$ (TP/OP) is plotted against the false positive rate $r_{\mathrm{FP}}$ (FP/ON) for various levels of III or $P_{\mathrm{I}}$. The area under the curve connecting the resulting points, AUC${}_{\mathrm{ROC}}$, is used as an indicator for the quality of the prediction (see Fig. 3). If the CDF for $P_{\mathrm{I}}$ is derived from the same set of landslides, r.randomwalk includes the option to randomly split the set of observed landslides into a set for parameter optimization, and one for validation. This is done for a user-defined number of times, exploiting multiple processors (see Sect. 2.3 and Fig. 1). It results in an ROC plot with multiple curves. Note that two ROC plots are produced: one of them builds on the original number of TN pixels. For the other one, the number of ON pixels is set to 5 times the number of OP pixels. Whilst the number of FP pixels remains unchanged, the number of TN pixels is modified accordingly. This procedure aims at normalizing the ROC curves in order to enable a comparison of the prediction qualities yielded for different study areas.

• If only one mass movement is considered, a longitudinal profile may be defined by a set of coordinates of the profile vertices. The observed and predicted (IIS $\ge$ 1 or $P_{\mathrm{I}}$ > 0) travel distances are measured and compared along this profile.

Test cases and results

Acheron rock avalanche, New Zealand

Area description and model parameterization

The Acheron rock avalanche in Canterbury, New Zealand (Fig. 4), was triggered approx. 1100 years BP (Smith et al., 2006). Within the present study, the release volume, $V=6.4$ million m${}^{\mathrm{3}}$, is approximated from the reconstruction of the pre-failure topography and is lower than the value of $V=7.5$ million m${}^{\mathrm{3}}$ estimated by Smith et al. (2006). We use a 10 m resolution DEM derived by stereo-matching of aerial photographs. Impact, release and deposition areas are derived from field and imagery interpretation as well as from data published by Smith et al. (2006). All random walks start from the highest pixel of the release area.

We use this case study for demonstrating how to compute the impact indicator index III from an elevation map, the release area, and the release volume. Before doing so, we have to analyse the influence of the pixel size and the parameters $n_{\mathrm{walks}}$, $R_{max}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}$, $f_{\mathit{\beta }}$, and $f_{\mathrm{d}}$ on the model result. Preliminary tests have shown that r.randomwalk yields plausible results with the number of random walks: $n_{\mathrm{walks}}={10}^{\mathrm{4}}$, $R_{max}=$ 10 m, $L_{\mathrm{ctrl}}=$ 1000 m, $L_{\mathrm{seg}}=$ 100 m, $f_{\mathit{\beta }}=$ 5, $f_{\mathrm{d}}=\mathrm{2}$, and a pixel size of 20 m. These values are taken as a basis to explore the sensitivity of the model results to the variation of each parameter and the best fit of the parameters in terms of the travel distance, AUC${}_{\mathrm{ROC}}$, and the size of the predicted impact area (Table 4). ${\mathit{\omega }}_T=11.62$${}^{\circ }$ , the angle of reach observed for the Acheron rock avalanche, is applied as the break criterion for all tests. Some of the tests are run in the back-calculation mode (flag b; see Tables 2 and 4).

III is computed by executing r.randomwalk 100 times, with the parameter values optimized according to Table 4. We explore an empirical–statistical relationship for ${\mathit{\omega }}_T$ derived from a compilation of 127 case studies (Fig. 5). The offset of the equation ($b$ in Eq. 4 and Fig. 5) is randomly sampled between the lower and upper envelopes of the regression. The quality of the prediction is evaluated using the ROC plot (see Figs. 1 and 3). Note that the Acheron rock avalanche (not included in the relationship developed in Fig. 5) is found close to the lower envelope, meaning that it was very mobile compared to most of the other events.

Acheron rock avalanche. (a) Panoramic view; photo: M. Mergili, 28 February 2015. (b) Location and geometry.

[Figure omitted. See PDF]

Empirical–statistical relationship relating the angle of reach ${\mathit{\omega }}_T$ to the volume $V$ of avalanching flows of rock or debris. The data are compiled from Scheidegger (1973), Legros (2002), Jibson et al. (2006), Evans et al. (2009), Sosio et al. (2012), and Guo et al. (2014).

[Figure omitted. See PDF]

Results

Figure 6 summarizes the findings of the test s 1–3 (see Table 4). Test 1 leads to the expected result that the predicted impact area increases with the number of random walks. However, the predicted impact area is also a function of the pixel size: with larger pixels, less random walks are needed to cover an area of similar size than with smaller pixels. Figure 6a further indicates that the possible impact area is not fully covered even at 10${}^{\mathrm{5}}$ random walks: no substantial flattening of the curves is observed. We conclude that (i) a very high value of $n_{\mathrm{walks}}$ would be necessary to fully cover the possible impact area, and (ii) this would lead to a substantial overestimation of the observed impact area.

Tests of the parameters $n_{\mathrm{walks}}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}$, $R_{max}$, $f_{\mathit{\beta }}$, $f_{\mathrm{d}}$, and the pixel size. Where ranges of values are given in bold, the model is run with 100 random samples constrained by the minima and maxima indicated. Where values given in bold are separated by commas, in these cases exactly these values are tested.

Test	$n_{\mathrm{walks}}$	$L_{\mathrm{ctrl}}$ (m)	$L_{\mathrm{seg}}$ (m)	$R_{max}$ (m)	$f_{\mathit{\beta }}$	$f_{\mathrm{d}}$	Pixel size (m)
1${}^{\mathrm{1},\mathrm{3}}$	10${}^{\mathrm{0}}$–10${}^{\mathrm{6}}$	1000	100	10	5	2	10, 20, 40
2${}^{\mathrm{2}}$	10${}^{\mathrm{4}}$	50, 1000	10–150	1000	5	2	10, 20, 40
3${}^{\mathrm{1},\mathrm{2},\mathrm{3}}$	10${}^{\mathrm{4}}$	50–1000${}^{\mathrm{2}}$	100	1000${}^{\mathrm{2}}$	5	2	10, 20, 40${}^{\mathrm{2}}$
		1000–4000${}^{\mathrm{1},\mathrm{3}}$		10${}^{\mathrm{1},\mathrm{3}}$			20${}^{\mathrm{1},\mathrm{3}}$
4${}^{\mathrm{1},\mathrm{3}}$	10${}^{\mathrm{4}}$	1000	100	0–120	5	2	20
5${}^{\mathrm{1},\mathrm{3}}$	10${}^{\mathrm{4}}$	1000	100	10	0–10	2	20
6${}^{\mathrm{1},\mathrm{3}}$	10${}^{\mathrm{4}}$	1000	100	10	5	1–10	20

Test criteria: ${}^{\mathrm{1}}$ impact area; ${}^{\mathrm{2}}$ travel distance $L_{max}$ (flag b); ${}^{\mathrm{3}}$ AUC${}_{\mathrm{ROC}}$.

Results of the tests 1–3 (number of test indicated in the yellow circle). Number of random walks plotted against (a) the impact area and (b) the area under the ROC curve. (c) Computed travel distance $L_{max}$ as a function of $L_{\mathrm{seg}}$ (in the legend, the corresponding value of $L_{\mathrm{ctrl}}$ is given in parentheses). (d) Computed $L_{max}$ as a function of $L_{\mathrm{ctrl}}$.

[Figure omitted. See PDF]

On the other hand, the quality of the prediction in terms of AUC${}_{\mathrm{ROC}}$ reaches a maximum at $n_{\mathrm{walks}}\approx {10}^{\mathrm{2}}$ (pixel size 40 m) or $n_{\mathrm{walks}}\approx {10}^{\mathrm{3}}$ (pixel size 20 m), decreasing with higher values of $n_{\mathrm{walks}}$. At a pixel size of 10 m, AUC${}_{\mathrm{ROC}}$ reaches a constant level at $n_{\mathrm{walks}}\approx {10}^{\mathrm{4}}$ (see Fig. 6b). We may conclude that excessive numbers of random walks lead to an overestimation of the impact area rather than to a better prediction quality. Coarser pixel sizes allow one to achieve the same level of coverage and the same prediction quality at lower values of $n_{\mathrm{walks}}$. However, the pixel size has to be fine enough to account for the main geometric characteristics of the process under investigation (see Sect. 4). All further tests are performed with $n_{\mathrm{walks}}={10}^{\mathrm{4}}$.

Sensitivity of impact area and AUC${}_{\mathrm{ROC}}$ to selected input parameters. The numbers of the corresponding tests (see Table 4) are indicated in the yellow circles. (a) Control distance $L_{\mathrm{ctrl}}$; (b) maximum run-up height $R_{max}$; (c) slope factor $f_{\mathit{\beta }}$; (d) direction factor $f_{\mathrm{d}}$.

[Figure omitted. See PDF]

Figure 6c illustrates that, at $L_{\mathrm{ctrl}}=1000$ m, the travel distance computed within the observed impact area decreases with increasing values of $L_{\mathrm{seg}}$ (tests 2 and 3 in Table 4). This pattern is well explained by Fig. 2b. At short segment lengths, the effects of flow paths frequently changing their direction are particularly evident for pixel sizes of 10 m and 20 m. $L_{max}$ drops below the observed value of 3550 m (see Fig. 4b) at 75 $\le L_{\mathrm{seg}}\le$ 100 m. With $L_{\mathrm{seg}}\ge 3050$ m, corresponding to the Euclidean distance between the release point and the terminal point of the Acheron rock avalanche, $L_{max}$ would also take a value of 3050 m. At $L_{\mathrm{ctrl}}=50$ m (only shown for a pixel size of 20 m), r.randomwalk tends to predict too long travel distances, compared to the observation. This phenomenon occurs as flow directions are not well defined in the relatively plane deposition zone of the Acheron rock avalanche; therefore, flow paths may frequently change their direction or even go backwards or in a circular way if such a behaviour is not impeded by sufficiently high values of $L_{\mathrm{ctrl}}$ (see Fig. 2a). Figure 6d indicates that this undesired behaviour (visible in the area marked by the X in the gray circle) disappears at $L_{\mathrm{ctrl}}$ > 200 m.

On the other hand, the value of $L_{\mathrm{ctrl}}$ should not be chosen too high as this may negatively impact the model performance. In the case of the Acheron rock avalanche, a drop in AUC${}_{\mathrm{ROC}}$ is observed between $L_{\mathrm{ctrl}}\approx 2000$ and $L_{\mathrm{ctrl}}\approx 2500$ m (Fig. 7a). This drop is explained by an increasing number of false negative pixels in those areas, which cannot be reached by the random walks due to the strict constraint of flow direction.

Within the tested ranges of parameter values, the quality of the prediction is highest at values of $R_{max}\approx \mathrm{5}$–10 m (see Fig. 7b) and $f_{\mathit{\beta }}\ge \mathrm{5}$ (see Fig. 7c), whilst it reaches it maximum at $f_{\mathrm{d}}\approx \mathrm{2}$–3 (see Fig. 7d). The predicted impact area increases with increasing $R_{max}$ and $f_{\mathrm{d}}$ whilst it decreases with increasing $f_{\mathit{\beta }}$.

Impact indicator score for the Acheron rock avalanche. (a) Classified III map. (b) ROC plot, building on normalized ON area (see Sect. 2.5).

[Figure omitted. See PDF]

Figures 6 and 7 indicate that the initial values of $n_{\mathrm{walks}}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}{R}_{max}$, $f_{\mathit{\beta }}$, $f_{\mathrm{d}}$, and the pixel size suggested in Sect. 3.1.1 and Table 4 are within the optimum range of values (see Sect. 4). Therefore, they are used for computing the impact indicator index for the Acheron rock avalanche (Fig. 8a). Concerning the break criteria, this can be classified as a forward analysis. As expected from Fig. 5, where the Acheron rock avalanche falls in between the envelopes of the relationship employed, the upper part of the observed impact area displays a value of III $=$ 1, whilst the remaining part of the observed impact area displays values of 1 > III > 0, decreasing towards the terminus. As the event was comparatively mobile within the context of the relationship used (see Sect. 3.1.1 and Fig. 5), the values of III are close to zero in the terminal area, and the area with III > 0 does not reach far beyond the observed terminus. Note that the maximum value of III is 0.8, meaning that 20 % of all model runs did not even start due to very high values of ${\mathit{\omega }}_T$ yielded with the randomized values of $b$ (see Fig. 5). Evaluation against the observed deposit yields a value of AUC${}_{\mathrm{ROC}}=0.94$ (see Fig. 8b). All values of AUC${}_{\mathrm{ROC}}$ shown in Figs. 6 and 7 and the ROC plot of Fig. 8b build on normalized ON areas (see Sect. 2.5).

III was generated within a computational time of 188 s.

Kao Ping Watershed, Taiwan

Area description and model parameterization

Between 7 and 9 August 2009, Typhoon Morakot struck Taiwan and triggered enormous landslides, causing significant land cover change (Fig. 9). More than 22 000 landslides were recorded in southern Taiwan (Lin et al., 2011). One of the hot spots of mass wasting was the Kao Ping Watershed (Wu et al., 2011), where the extremely heavy rainfall (in total, more than 2000 mm depth and 90 h duration) triggered a catastrophic landslide in the Hsiaolin Village (Kuo et al., 2013).

Location, terrain and landslide inventory of the Kao Ping Watershed, Taiwan. Comparison of the satellite images illustrates the landslide-induced land cover changes associated with the Typhoon Morakot. The landslide inventory builds on the interpretation of the FORMOSAT-2 imagery.

[Figure omitted. See PDF]

We consider a 61.5 km${}^{\mathrm{2}}$ subset of the Kao Ping Watershed for computing the landslide impact probability $P_{\mathrm{I}}$, based on the observed landslide release areas. In all, 207 landslides are mapped in the shale, sandstone, and colluvium slopes (see Fig. 9). A 10 m DEM is used along with an inventory of the landslide impact areas. Release and deposition areas are extracted from the inventory. We employ the values of $n_{\mathrm{walks}}$, $R_{max}$, $f_{\mathit{\beta }}$, $f_{\mathrm{d}}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}$ resulting from the optimization procedure for the Acheron rock avalanche (see Sect. 3.1.1), and a pixel size of 20 m. $P_{\mathrm{I}}$ is computed as follows:

• A set of random walks ($n_{\mathrm{walks}}={10}^{\mathrm{4}}$) is started from each release point (i.e., the highest pixel of each landslide). Each random walk stops as soon as it would leave the impact area of the same landslide (back calculation, flag b).

• After completing all random walks for the study area, the statistical distribution of ${\mathit{\omega }}_T$ is analysed. All landslides with $L_{max}$ < 100 m are excluded. A fraction of 20 % out of all landslides (i.e., all values of ${\mathit{\omega }}_T$ associated with those landslides) is randomly selected and retained for validation. Using visual comparison, we have identified the log-normal distribution as the most suitable type of distribution for this purpose. Consequently, the log-normal CDF stands for the probability that a moving mass point leaves the observed impact area at or below the associated threshold of ${\mathit{\omega }}_T$.

• We perform a forward analysis of $P_{\mathrm{I}}$ by starting a set of random walks ($n_{\mathrm{walks}}={10}^{\mathrm{4}}$) from the release points of the retained landslides, and assigning the cumulative density associated with the average angle of path to each pixel. The result is validated against the observed deposition zones of the retained landslides by means of an ROC plot.

• Steps 2. and 3. are repeated for 100 randomly selected subsets (parallel processing is applied). The final map of $P_{\mathrm{I}}$ is generated by applying for each pixel the maximum of the values yielded by all the model runs.

We refer to this work flow as test 1 and repeat the analysis with starting random walks not only from the release points but also from all the pixels within the observed release areas (test 2). This means that the CDF is derived from a much larger sample of data than when considering only one point per landslide for starting random walks. We exclude all sets of random walks yielding $L_{max}$ < 100 m, use a log-normal CDF and start a set of only 10${}^{\mathrm{3}}$ random walks from each release pixel for computing $P_{\mathrm{I}}$.

Results

Starting sets of 10${}^{\mathrm{4}}$ random walks from the highest points of all landslides (test 1) results in a range of values of 16.0 $\le {\mathit{\omega }}_T\le 43.5$${}^{\circ }$, an average of 30.4${}^{\circ }$, and a standard deviation of 5.2${}^{\circ }$ (derived from $n=132$ landslides, excluding those with $L_{max}$ < 100 m). Repeating the analysis with 10${}^{\mathrm{4}}$ random walks started from each pixel within the landslide release areas (test 2), we observe a range of values 16.4 $\le {\mathit{\omega }}_T\le 44.1$${}^{\circ }$, an average of 26.9${}^{\circ }$, and a standard deviation of 4.8${}^{\circ }$ ($n=1563$). Figure 10 illustrates the histograms, probability density, and cumulative density functions derived from both analyses. Even though the ranges of values are similar in both tests, test 1 yields (i) a higher average of ${\mathit{\omega }}_T$ and (ii) a broader range of values than test 2; (i) is explained by the fact that those random walks starting from lower parts of the release areas are expected to leave the observed impact area at lower values of ${\mathit{\omega }}_T$; (ii) is most likely the consequence of a number of rather small landslides with high or low values of ${\mathit{\omega }}_T$ strongly reflected in the statistics. Such outliers are less prominent in the statistics of test 2 due to the much higher number of cases, most of them related to the larger landslides.

Each of the impact probabilities shown in Fig. 11 represents the overlay of 100 analyses where random sets of 80 % of the landslides are used for deriving the CDF and the remaining 20 % are used for computing the impact probabilities. The maps illustrate the maximum values of $P_{\mathrm{I}}$ out of the overlay of the 100 results. Each of the results is derived using a slightly different CDF. Both tests yield largely similar patterns of $P_{\mathrm{I}}$. We note that (i) test 2 predicts larger impact areas and higher values of $P_{\mathrm{I}}$ than test 1, and (ii) some random walks take the wrong direction in test 2 (indicated by “1” in the yellow circle in Fig. 11b), a phenomenon not observed for test 1; (i) is explained by the higher number – and the broader distribution – of release pixels in test 2, compared to test 1. The reason for (ii) is that random walks starting from the highest point of an observed landslide are forced to flow into the observed landslide area (test 1), a constraint not applicable when starting random walks from each release pixel (test 2). In this case it happens that pixels located at or near a crest produce random walks in both directions. In test 1, the computational time amounted to 63 s for deriving the CDF and 8613 s for calculating $P_{\mathrm{I}}$. In test 2, these times increased to 1719 and 9752 s, respectively. The relatively slight increase with regard to $P_{\mathrm{I}}$ results from the reduced value of $n_{\mathrm{walks}}$ in test 2.

Histograms, probability densities, and cumulative densities of ${\mathit{\omega }}_T$ of mass movements in the test area in the Kao Ping Watershed. (a) Result for a set of 10${}^{\mathrm{4}}$ random walks started from the highest point of each landslide (test 1). (b) Result for a set of 10${}^{\mathrm{4}}$ random walks started from each pixel within the release areas of all landslides (test 2).

[Figure omitted. See PDF]

The prediction quality is tested for each of the 100 model runs for the two tests, producing sets of 100 ROC curves (Fig. 12). AUC${}_{\mathrm{ROC}}=0.917\pm 0.038$ for test 1 and 0.920 $\pm$ 0.029 for test 2, both computed with the original number of TN pixels (see Sect. 2.5).

In contrast, the procedures demonstrated in the two tests vary strongly in their scope of applicability. We have demonstrated the methodologies by back calculating observed landslides. As soon as this is done, one may go one step further:

• The methodology shown in test 1 can be employed to make forward predictions for defined expected future landslides, given that a sufficient set of observed landslides of similar behaviour is available to derive the CDF.

• The methodology demonstrated in test 2 can be used in combination with maps of landslide release probability to explore the composite probability of a landslide impact (see Sects. 2.4 and 4).

Impact probability in the range 0–1. (a) Result of test 1 (random walks starting from the highest point of each landslide; cumulative density according to Fig. 10a). (b) Result of test 2 (random walks starting from all release pixels; cumulative density according to Fig. 10b).

[Figure omitted. See PDF]

ROC plots illustrating the prediction quality of (a) test 1 and (b) test 2, using the original number of TN pixels (see Sect. 2.5).

[Figure omitted. See PDF]

In either case the statistics (see Fig. 10) have to be derived with the same type of approach later used for producing the $P_{\mathrm{I}}$ map.

Gunt Valley, Tajikistan

Area description and model parameterization

As most mountain areas worldwide, the Pamir of Tajikistan experiences a significant retreat of the glaciers. One of the consequences thereof consists in the formation and growth of lakes, some of which are subject to glacial lake outburst floods (GLOFs), which may evolve into destructive debris flows (Mergili and Schneider, 2011; Mergili et al., 2013; Gruber and Mergili, 2013). No records of historic GLOFs in the test area are known to the authors. However, in August 2002 a GLOF in the nearby Shakhdara Valley evolved into a debris flow, which destroyed the village of Dasht, claiming dozens of lives (Mergili et al., 2011).

The frequency of such events is low and historical data are sparse. Consequently, possible travel distances of GLOFs may not be derived in a purely statistical way. Instead, we have to use published empirical–statistical relationships and simple rules to produce an impact indicator score (IIS) map.

We compute IIS with regard to GLOFs for a 2106 km${}^{\mathrm{2}}$ study area in the Gunt Valley (Fig. 13). The analysis builds on the ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer – Global Digital Elevation Model) V2 and the coordinates and characteristics (estimates of $V$ and $Q_{\mathrm{p}}$) of 113 lakes in the area (Gruber and Mergili, 2013).

A set of random walks ($n_{\mathrm{walks}}={10}^{\mathrm{4}}$) is routed from the outlet of each lake through the DEM. Six break criteria are combined to compute IIS, partly following Gruber and Mergili (2013). The relationships and rules employed as break criteria are summarized in Table 5. Rule 1 is applied with ${\mathit{\omega }}_T=11$${}^{\circ }$ (test 1 – according to Haeberli, 1983 and Huggel et al., 2003, 2004a, b for debris flows from glacier- or moraine-dammed lakes, and Zimmermann et al., 1997 for coarse- and medium-grained debris flows) and with ${\mathit{\omega }}_T=\mathrm{7}$${}^{\circ }$ (test 2 – Zimmermann et al., 1997 for fine-grained debris flows). All other rules and relationships are used for both tests. For each pixel, IIS consists in the number of relationships or rules predicting an impact (i.e., IIS takes values in the range 0–6).

$R_{max}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}$, $f_{\mathit{\beta }}$, and $f_{\mathrm{d}}$ are set to the optimum values found for the Acheron rock avalanche, the pixel size is set to 60 m.

Results

Figure 14 illustrates the possible impact areas of GLOFs in the Gunt Valley study area according to the relationships listed in Table 5.

Figure 14a shows the impact indicator score IIS i.e., the number of relationships predicting an impact, resulting from test 1 (rule 1 applied with ${\mathit{\omega }}_T=11$${}^{\circ }$). Except for one prominent exception, IIS > 3 (possible debris flow impact) only for the largely uninhabited upper portions of the tributaries to the Gunt Valley. In contrast, a possible flood impact (1 $\le$ IIS $\le \mathrm{3}$) is predicted for much of the main valley. test 2 (rule 1 applied with ${\mathit{\omega }}_T=\mathrm{7}$${}^{\circ }$) predicts a possible debris flow impact also for part of the main valleys (see Fig. 14b). The IF (per cent of random walks impacting each pixel) for test 1 is shown in Fig. 14c for a subsection of the test area, classified by quantiles. IF is strongly governed by the width of the movement, i.e. by the local topography, and may serve as a surrogate for the expected depth rather than as for the probability of an impact.

The test area in the Gunt Valley, Tajikistan. (a) Location, topography, glaciers and lakes. (b) Proglacial lake in the upper Varshedzdara Valley; photo: M. Mergili, 18 August 2011.

[Figure omitted. See PDF]

Empirical–statistical relationships and simple rules used for computing the IIS of GLOFs in the Gunt Valley (see Table 3).

IDtest	Relationship	Reference	Process
1${}^{\mathrm{1}}$	${\mathit{\omega }}_T=11$${}^{\circ }$	Haeberli (1983); Zimmermann et al. (1997);	Flood or debris flow
		Huggel et al. (2003, 2004a, b)
1${}^{\mathrm{2}}$	${\mathit{\omega }}_T=\mathrm{7}$${}^{\circ }$	Zimmermann et al. (1997)
2${}^{\mathrm{1},\mathrm{2}}$	${\mathit{\omega }}_T=18Q_{\mathrm{p}}^{-0.07}$	Huggel (2004)
3${}^{\mathrm{1},\mathrm{2}}$	$L_{\max }=1.9V^{0.16}{Z}^{0.83}$	Rickenmann (1999)${}^{\mathrm{3}}$
4${}^{\mathrm{1},\mathrm{2}}$	${\mathit{\omega }}_T=\mathrm{6}$${}^{\circ }$		Flood
5${}^{\mathrm{1},\mathrm{2}}$	${\mathit{\omega }}_T=\mathrm{4}$${}^{\circ }$
6${}^{\mathrm{1},\mathrm{2}}$	${\mathit{\omega }}_T=\mathrm{2}$${}^{\circ }$	Haeberli (1983); Huggel et al. (2004a)

${}^{\mathrm{1},\mathrm{2}}$ ID(s) of test(s) where the rule or relationship is applied. ${}^{\mathrm{3}}$ A bulking factor of 5 is applied to $V$ (modified after Iverson, 1997).

Note that Fig. 14 only indicates the tendency of an already released GLOF to impact certain pixels. It does not provide any information on the susceptibility of a certain lake to produce a GLOF at all. Earlier, Mergili and Schneider (2011) and Gruber and Mergili (2013) have attempted to combine GLOF release indicators with impact indicators and land cover maps to generate hazard and risk indicator maps. However, the results of their studies may underestimate the possible impact areas as the travel distance was computed on a pixel-to-pixel basis, possibly yielding too low values of ${\mathit{\omega }}_T$ (see Figs. 2 and 6).

The robustness and appropriateness of the rules and relationships for low-frequency events, such as GLOFs (see Table 5), is questionable. The rules building on a unique value of ${\mathit{\omega }}_T$ overpredict the possible impact areas for those lakes where not enough water is available to produce a flood in downstream valleys. Applying the rules and relationships for debris flows implies a blind assumption that enough entrainable sediment is available to produce a debris flow. Whilst ${\mathit{\omega }}_T\ge 11$${}^{\circ }$ is considered the worst case for debris flows of GLOFs from glacier- or moraine-dammed lakes in the European Alps according to Haeberli (1983) and Huggel et al. (2002), ${\mathit{\omega }}_T=9.3$${}^{\circ }$ was measured for the 2002 Dasht Event, the only well-documented GLOF near the test area (Mergili et al., 2011). Also the relationship proposed by Rickenmann (1999) severely underestimates the travel distance of this event, even when massive bulking is assumed. Applying ${\mathit{\omega }}_T=\mathrm{7}$${}^{\circ }$ as given by Zimmermann et al. (1997) for fine-grained debris flows might be more suitable as worst-case assumption for debris flows from GLOFs in the Pamir, even though this threshold leads to very conservative predictions.

We have measured computational times of 1520 s for test 1 and 1556 s for test 2.

Possible GLOF impact areas in the Gunt Valley, Tajikistan. (a) Impact indicator score derived with test 1. (b) Impact indicator score derived with test 2. (c) Impact frequency derived with test 1, classified by quantiles.

[Figure omitted. See PDF]

Discussion

Whilst conceptual tools are commonly applied for routing mass movements at medium and broad scales, most of them use single values or rules as break criteria, disregarding the high degree of uncertainty (e.g., Gamma, 2000; Wichmann and Becht, 2003; Huggel et al., 2002; Horton et al., 2013; Blahut et al., 2010). r.randomwalk introduces a set of tools to deal with uncertain break criteria in a flexible way, depending on the quality of rules or relationships available. In general, empirical–statistical relationships represent rough simplifications as mass movement processes may also stop when reaching valleys of higher order, run against opposite slopes or lose energy when bending sharply. However, relatively robust rules or relationships exist for the most common types of processes such as rock avalanches (Scheidegger, 1973; see Fig. 5) or debris flows (Rickenmann, 1999). They build on data sets large enough to derive meaningful envelopes and to compute impact indicator indices with r.randomwalk. Relationships for less frequent types of processes are less robust as it was illustrated for GLOFs (Haeberli, 1983; Zimmermann et al.; 1997; Huggel et al., 2002; Huggel, 2004; see Sect. 3.3.2). In such cases we recommend to compute impact indicator scores building on more than one model, as shown by Gruber and Mergili (2013) and in the present work. Impact indicator indices and scores are mainly useful for anticipating the possible impact area of expected single events (see Sect. 3.1.2), or for application at broader scales (see Sect. 3.3.2).

The impact probability is useful for predicting possible impact areas of mass movements in areas where many events are documented, but the volumes of possible future events are not known. Whilst in the present paper it was demonstrated how to compute impact probabilities related to observed release areas, r.randomwalk also includes the option to combine the impact probability with the release probability $P_{\mathrm{R}}$ (see Table 1). Landslide release probability (susceptibility) maps are often produced from a landslide inventory and a set of environmental layers (e.g., Guzzetti, 2006). Starting random walks from each single pixel of a study area, and combining the release probability of this pixel with the impact probability allows one to produce a composite probability $P_{\mathrm{I},\mathrm{C}}$ map. Doing this is non-trivial and requires specific strategies. It is therefore covered in a separate article (Mergili and Chu, 2015). Gruber and Mergili (2013) have combined release and impact indicator scores for various types of high-mountain hazards, and overlaid the results with a land cover data set to produce a risk indicator score.

The sensitivity of r.randomwalk to variations of the parameters $n_{\mathrm{walks}}$, $R_{max}$, $f_{\mathit{\beta }}$, $f_{\mathrm{d}}$, $L_{\mathrm{ctrl}}$, $L_{\mathrm{seg}}$ (see Sect. 2.2) and the pixel size were tested for the Acheron rock avalanche. Even though the optimized values are applied also to the other cases in the present work, this issue requires further investigation, also with regard to the scale of the processes. This is particularly true for the pixel size, which has to be fine enough not to lose the geometrical characteristics governing the motion (Blahut et al., 2010). Furthermore, coarser pixels and a larger number of random walks make results more conservative. $R_{max}$, $f_{\mathit{\beta }}$, and $f_{\mathrm{d}}$ control the degree of lateral spreading and therefore influence the conservativeness of the results. In the future we plan to compare the performance of r.randomwalk to software tools using multiple flow direction algorithms (e.g., Flow-R; Horton et al., 2013) in terms of computational times and prediction success.

Overestimating the travel distance at a certain pixel is avoided by choosing sufficiently high values of $L_{\mathrm{seg}}$ (see Fig. 6c). Shorter travel distances at a certain pixel are associated with higher values of $\mathit{\omega }$ and, consequently, larger predicted impact areas, i.e. more conservative results that are desirable for many applications. The values of $R_{max}$ leading to the best prediction quality are considerably lower than run-up height observed for the Acheron rock avalanche. This phenomenon is explained by the facts that (i) the observed maximum run-up height refers to a limited area, whilst r.randomwalk applies the run-up height defined by $R_{max}$ in any place; and (ii) not all random walks reach the bottom of the valley before running up.

We have demonstrated how to estimate the prediction quality of III and $P_{\mathrm{I}}$ maps. Where sufficient reference data are available to prove the validity of the model, the results may be applied for hazard zoning. Where data are not available, the outcomes of r.randomwalk are suitable for broad-scale overviews of possibly affected areas, which have to be considered as rough indicators only. A suitable level of spatial aggregation may be necessary in such cases (Gruber and Mergili, 2013).

r.randomwalk includes a break criterion building on the two-parameter friction model of Perla et al. (1980) (see Sect. 1 and Table 3), which can be used to compute flow velocities (e.g., Wichmann and Becht, 2013; Mergili et al., 2012a; Horton et al., 2013). Evaluating this functionality has to build on (i) specific strategies for the sensitivity analysis and optimization of multiple parameters and (ii) a sound comparison with the outcome of physically based models. This effort will be presented in a separate article (Krenn et al., 2015). Further, the parameter sensitivity and optimization code AIMEC (Fischer, 2013) can be directly coupled to r.randomwalk.

Conclusions

We have introduced the open-source GIS tool r.randomwalk, designed for conceptual modelling of the propagation of mass movements. r.randomwalk offers built-in functions for considering uncertainties and for validation. Employing a set of three contrasting test areas, we have demonstrated (i) the possibility to combine results yielded with various break criteria into one impact indicator score; (ii) the option to explore multiple computational cores for combining the results obtained with many randomized parameter combinations into an impact indicator index; (iii) the possibility to back calculate the CDF of the angles of reach of observed landslides, and to use this CDF to make forward predictions of the impact probability; and (iv) integrated functions for the validation and visualization of the results. This includes strategies to properly separate the data sets for parameter optimization and model validation.

We have further shown that controls for smoothing of the flow path and the avoidance of circular flows have to be introduced to avoid underestimating travel distances and impact areas. The number of random walks executed for each mass point and the pixel size influence the level of conservativeness of the results rather than the quality of the prediction. The scope of applicability of r.randomwalk strongly depends on the availability of robust break criteria and on the availability of reference data for evaluation.

Code availability

The model codes, a user manual, the scripts used for starting the tests presented in Sect. 3 and some of the test data are available at http://www.mergili.at/randomwalk.html.

Acknowledgements

The work was conducted as part of the international cooperation project “A GIS simulation model for avalanche and debris flows” funded by the Austrian Science Fund (FWF) and the German Research Foundation (DFG). Further, the support of Massimiliano Alvioli, Matthias Benedikt, Yi-Chin Chen, Ivan Marchesini, and Tim Davies is acknowledged. Edited by: T. Poulet