r . randomwalk v 1 . 0 , a multi-functional conceptual tool for mass movement routing

r.randomwalk v1.0, a multi-functional conceptual tool for mass movement routing M. Mergili, J. Krenn, and H.-J. Chu Geomorphological Systems and Risk Research, Department of Geography and Regional Research, University of Vienna, Universitätsstraße 7, 1190 Vienna, Austria Institute of Applied Geology, University of Natural Resources and Life Sciences (BOKU), Peter-Jordan-Straße 70, 1190 Vienna, Austria Department of Geological Sciences, University of Canterbury, Private Bag 4800, Christchurch, New Zealand Department of Geomatics, National Cheng Kung University, 1 University Road, Tainan 701, Taiwan


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 andHungr, 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.
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
2.1 Computational implementation r.randomwalk is implemented as a raster module of the opensource 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 ) 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.mainwith 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.
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), andGruber andMergili (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).
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).During the pixel-to-pixel routing procedure, turns of > 90 • 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 ctrl defines the backward distance of each step over which the horizontal distance of motion has to increase (Fig. 2a).
The probability P px of any other neighbour pixel px to become the target pixel is where n is the total number of valid neighbour pixels, and β is the local slope between the current pixel and the considered neighbour pixel.f d and f β are weighting factors for the perpetuation of the flow direction and for the slope.
Control length Backward distance of each step of a random walk over which the horizontal distance of motion has to increase (see Fig. 2a) Segment length Length of segments used for computing L max (see Fig. 2b) Direction factor Factor for weighting the perpetuation of the movement direction during routing Slope factor Factor for weighting the slope during routing Applicable to lake outburst events (see Table 3) Parameters needed by the rules and relationships applied (see Table 3) Release probability Spatial probability that a mass movement is released from a given pixel (Mergili and Chu, 2015) RIS Published relationships (see Table 3) 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 seg .The travel distance L max is defined as the sum of the length of all segments (see Fig. 2b).Larger values of L 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 ω 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 semideterministic 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 I in the range 0-1 is computed from a user-defined sample of observed values of tan(ω 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(ω 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 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 ω T .The resulting CDF can be used as input to estimate P I .

Overlay of random walk results
The overlay of individual random walks operates at two levels: 1. 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 I -is derived

ID Equation Examples
Reference Process a b c from the random walk with the shortest travel distance (i.e., the straightest flow path and the highest value of ω) at the considered pixel.
2. 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 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 I or IIS can be automatically overlaid with a release probability (P R ; result: composite probability P I,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).
-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 ≥ 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 I : the true positive rate r TP (TP/OP) is plotted against the false positive rate r FP (FP/ON) for various levels of III or P I .The area under the curve connecting the resulting points, AUC ROC , is used as an indicator for the quality of the prediction (see Fig. 3).If the CDF for P 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 userdefined 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 ≥ 1 or P I > 0) travel distances are measured and compared along this profile.
3 Test cases and results

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 3 , is approximated from the reconstruction of the pre-failure topography and is lower than the value of V = 7.5 million m 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 walks , R max , L ctrl , L seg , f β , and f d on the model result.Preliminary tests have shown that r.randomwalk yields plausible results with the number of random walks: 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 ROC , and the size of the predicted impact area (Table 4).ω T = 11.62 • , 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 ω 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.

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 5 random walks: no substantial flat-  tening of the curves is observed.We conclude that (i) a very high value of n 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.
On the other hand, the quality of the prediction in terms of AUC ROC reaches a maximum at n walks ≈ 10 2 (pixel size 40 m) or n walks ≈ 10 3 (pixel size 20 m), decreasing with higher values of n walks .At a pixel size of 10 m, AUC ROC reaches a constant level at n walks ≈ 10 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 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 walks = 10 4 .
Figure 6c illustrates that, at L ctrl = 1000 m, the travel distance computed within the observed impact area decreases with increasing values of L 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 ≤ L seg ≤ 100 m.With L seg ≥ 3050 m, corre-  sponding 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 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 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 ctrl > 200 m.
On the other hand, the value of L 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 ROC is observed between L ctrl ≈ 2000 and L ctrl ≈ 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 ≈ 5-10 m (see Fig. 7b) and f β ≥ 5 (see Fig. 7c), whilst it reaches it maximum at f d ≈ 2-3 (see Fig. 7d).The predicted impact area in-creases with increasing R max and f d whilst it decreases with increasing f β .
Figures 6 and 7 indicate that the initial values of n walks , L ctrl , L seg R max , f β , f d , and the pixel size suggested in Sect.3.1.1and 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.1and 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 ω T yielded with the randomized values of b (see Fig. 5).Evaluation against the observed deposit yields a value of AUC ROC = 0.94 (see Fig. 8b).All values of AUC 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.

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).We consider a 61.5 km 2 subset of the Kao Ping Watershed for computing the landslide impact probability P 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 walks , R max , f β , f d , L ctrl , L seg resulting from the optimization procedure for the Acheron rock avalanche (see Sect. 3.1.1),and a pixel size of 20 m.P I is computed as follows: 1.A set of random walks (n walks = 10 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).
2. After completing all random walks for the study area, the statistical distribution of ω 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 ω 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 ω T .
3. We perform a forward analysis of P I by starting a set of random walks (n walks = 10 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.
4. Steps 2. and 3. are repeated for 100 randomly selected subsets (parallel processing is applied).The final map of P 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 3 random walks from each release pixel for computing P I .

Results
Starting sets of 10 4 random walks from the highest points of all landslides (test 1) results in a range of values of 16.0 ≤ ω T ≤ 43.5 • , an average of 30.4 • , and a standard deviation of 5.2 • (derived from n = 132 landslides, excluding those with L max < 100 m).Repeating the analysis with 10 4 random walks started from each pixel within the landslide release areas (test 2), we observe a range of values 16.4 ≤ ω T ≤ 44.1 • , an average of 26.9 • , and a standard deviation of 4.8 • (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 ω 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 ω T ; (ii) is most likely the consequence of a number of rather small landslides with high or low values of ω 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 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 I .We note that (i) test 2 predicts larger impact areas and higher values of P 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 I .In test 2, these times increased to 1719 and 9752 s, respectively.The relatively slight increase with regard to P I results from the reduced value of n walks in test 2.
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 ROC = 0.917 ± 0.038 for test 1 and 0.920 ± 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 demon- strated 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).
In either case the statistics (see Fig. 10) have to be derived with the same type of approach later used for producing the P I map.

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 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 p ) of 113 lakes in the area (Gruber and Mergili, 2013).
A set of random walks (n walks = 10 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 ω T = 11 • (test 1 -according to Haeberli, 1983and Huggel et al., 2003, 2004a, b

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 ω T = 11 • ).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 ≤ IIS ≤ 3) is predicted for much of the main valley.test 2 (rule 1 applied with ω T = 7 • ) 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.
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 ω T (see Figs. 2 and 6).Table 5. Empirical-statistical relationships and simple rules used for computing the IIS of GLOFs in the Gunt Valley (see Table 3).

IDtest Relationship
Reference Process L max = 1.9V 0.16 Z 0.83 Rickenmann (1999 1,2 ID(s) of test(s) where the rule or relationship is applied. 3A bulking factor of 5 is applied to (modified after Iverson, 1997).
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 ω 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 ω T ≥ 11 • 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), ω T = 9.3 • 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 ω T = 7 • as given by Zimmermann et al. (1997) for fine-grained debris flows might be more suitable as worstcase 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.

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.Re- lationships 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 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 I,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 highmountain 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 walks , R max , f β , f d , L ctrl , L 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 β , and f 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 seg (see Fig. 6c).Shorter travel distances at a certain pixel are associated with higher values of ω 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 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 broadscale 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.

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

Figure 2 .
Figure 2. Control length L ctrl and segment length L seg .(a) Application of L ctrl to avoid sharp bending of the flow.(b) Smoothing of the flow path by introducing segments with maximum length of L seg .
f d is governed by the input parameter d: f d = d 2 for the same flow direction as the previous one, f d = d for a 45 • turn and f d = 1 for a 90 • turn.

Figure 3 .
Figure3.Model validation with an ROC plot, relating the false positive rate r FP and the true positive rate r TP .This way of validation is suitable for predictor raster maps in the range 0-1, such as III or P I .It can also be used for binary predictor maps (0 or 1).In such a case AUC ROC is computed from two threshold levels only.

Figure 6 .
Figure 6.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 seg (in the legend, the corresponding value of L ctrl is given in parentheses).(d) Computed L max as a function of L ctrl .

Figure 7 .
Figure 7. Sensitivity of impact area and AUC ROC to selected input parameters.The numbers of the corresponding tests (see Table 4) are indicated in the yellow circles.(a) Control distance L ctrl ; (b) maximum run-up height R max ; (c) slope factor f β ; (d) direction factor f d .
Figure 7. Sensitivity of impact area and AUC ROC to selected input parameters.The numbers of the corresponding tests (see Table 4) are indicated in the yellow circles.(a) Control distance L ctrl ; (b) maximum run-up height R max ; (c) slope factor f β ; (d) direction factor f d .

Figure 8 .
Figure 8. 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 9 .
Figure 9. 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 10 .
Figure 10.Histograms, probability densities, and cumulative densities of ω T of mass movements in the test area in the Kao Ping Watershed.(a) Result for a set of 10 4 random walks started from the highest point of each landslide (test 1).(b) Result for a set of 10 4 random walks started from each pixel within the release areas of all landslides (test 2).

Figure 11 .
Figure 11.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 12 .
Figure 12.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).
for debris flows from glacieror moraine-dammed lakes, and Zimmermann et al., 1997 for coarse-and medium-grained debris flows) and with ω T = 7 • (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 ctrl , L seg , f β , and f d are set to the optimum values found for the Acheron rock avalanche, the pixel size is set to 60 m.

Figure 14 .
Figure 14.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.

Table 1 .
Summary of the key variables used in r.randomwalk.

Table 2 .
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 • is relevant for some applications.

Table 3 .
Types of rules and relationships supported by r.randomwalk.ω T is angle of reach, L max is travel distance, V is volume of motion, Z is elevation loss, Q p is peak discharge at release, and v T is velocity at termination.

Table 4 .
Tests of the parameters n walks , L ctrl , L seg , R max , f β , f 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.