Interactive comment on “ A strategy for GIS-based 3-D slope stability modelling over large areas ”

Abstract. GIS-based deterministic models may be used for landslide susceptibility mapping over large areas. However, such efforts require specific strategies to (i) keep computing time at an acceptable level, and (ii) parameterize the geotechnical data. We test and optimize the performance of the GIS-based, 3-D slope stability model r.slope.stability in terms of computing time and model results. The model was developed as a C- and Python-based raster module of the open source software GRASS GIS and considers the 3-D geometry of the sliding surface. It calculates the factor of safety (FoS) and the probability of slope failure (Pf) for a number of randomly selected potential slip surfaces, ellipsoidal or truncated in shape. Model input consists of a digital elevation model (DEM), ranges of geotechnical parameter values derived from laboratory tests, and a range of possible soil depths estimated in the field. Probability density functions are exploited to assign Pf to each ellipsoid. The model calculates for each pixel multiple values of FoS and Pf corresponding to different sliding surfaces. The minimum value of FoS and the maximum value of Pf for each pixel give an estimate of the landslide susceptibility in the study area. Optionally, r.slope.stability is able to split the study area into a defined number of tiles, allowing parallel processing of the model on the given area. Focusing on shallow landslides, we show how multi-core processing makes it possible to reduce computing times by a factor larger than 20 in the study area. We further demonstrate how the number of random slip surfaces and the sampling of parameters influence the average value of Pf and the capacity of r.slope.stability to predict the observed patterns of shallow landslides in the 89.5 km2 Collazzone area in Umbria, central Italy.


Introduction
Landslide susceptibility is the spatial probability of landslide occurrence, based on local terrain conditions (Brabb, 1984;Guzzetti et al., 1999).The susceptibility to landslides can be determined using statistical and physically based models (Guzzetti et al., 1999;Van Westen, 2000;Guzzetti, 2006;Van Westen et al., 2006).Most commonly, modelling of the spatial probability of shallow landslides for small catchments relies upon the use of physically based ("deterministic") models (Van Westen et al., 2006).These models build on the limit equilibrium concept, resulting in a factor of safety (FoS; Carson and Kirkby, 1972;Crozier, 1986;Duncan and Wright, 2005) of the failure surface.FoS is given by the dimensionless ratio between the resisting (stabilizing) forces and the driving (destabilizing) forces.
In combination with a raster GIS, FoS calculations most commonly build on the assumption of a planar slope of infinite length, with the potential failure surface parallel to the topographic surface (infinite slope stability model; Van Westen and Terlien, 1996;Burton and Bathurst, 1998;Xie et al., 2004a;Baum et al., 2008;Godt et al., 2008;Mergili et al., 2014;Raia et al., 2014).However, the infinite slope stability model is well suited only for shallow slope stability in frictional materials, and is less appropriate for cohesive materials (Mergili et al., 2014).It is suitable only for sufficiently large length-to-depth ratios of landslide release areas (Griffiths et al., 2011a;Milledge et al., 2012;Tiwari et al., 2014).
To evaluate the stability/instability conditions of slopes susceptible to deep-seated landslides, the zone above a known, inferred or hypothetical failure surface is divided into vertical slices of equal or different sizes.The resisting (stabilizing) forces R and the driving (destabilizing) forces T are computed for each slice, and summed up linearly to obtain a single value of FoS for the entire slope.First applied to 2-D cross-sections (Duncan and Wright, 2005), this type of model was extended to 3-D topographies and failure surfaces (e.g.Hovland, 1977;Hungr, 1987;Hungr et al., 1989).More details are given in Mergili et al. (2014).Few attempts were made to develop sliding-surface models applicable at the regional scale, coupled to GIS (Reid et al., 2000;Xie et al., 2003Xie et al., , 2004bXie et al., , c, 2006;;Marchesini et al., 2009;Jia et al., 2012).A recent study of Mergili et al. (2014) indicates that, also for shallow landslides, more complex slip surface models might perform slightly better in reproducing the observed landslide areas than the infinite slope stability model.
The broad-scale GIS implementation of sliding surface models faces two major challenges: (i) the spatial distribution of the geotechnical and geometric parameters is uncertain; and (ii) a very large number of possible slip surfaces has to be tested using a reasonably fine pixel spacing, requiring strategies to keep computational times at an acceptable level.
Geotechnical uncertainty is due to inherent spatial and temporal variability of terrain materials (Hicks and Spencer, 2010;Suchomel and Mašin, 2010;Griffiths et al., 2011b) but also due to limitations of laboratory analysis (Di Matteo et al., 2013).It can affect also the geometric uncertainty, i.e. the depth and shape of the failure surface.Predefining the failure surface in a deterministic way most likely leads to an overestimation of FoS or to an underestimation of the slope failure probability (Griffiths et al., 2011b).
A considerable amount of literature has been published on the uncertainty of geotechnical parameters and slip surface geometry.Recently, Johari et al. (2013) introduced five classes of probabilistic methods that have been used for the analysis of the stability of slopes.He defined these classes as (i) random sampling (Monte Carlo simulation), (ii) analytic methods, (iii) approximate methods, (iv) response surface method and (v) stochastic finite element method (see Johari et al., 2013, and references therein for detailed descriptions of each method).Even though a number of researchers have tried to demonstrate that the methods pertaining to the last four classes can perform reasonably well and in a faster and more efficient way, there is broad agreement that the Monte Carlo approach remains the best method (e.g.Low, 2007;Husein Malkawi et al., 2000;Wong, 1985;Ishii and Suzuki, 1987;Tan et al., 2013).
Existing examples and practices of exploiting parallel computing in GRASS GIS show relevant progress in the topic and are worth being considered for a vast number of applications (Sorokine, 2007;Liu et al., 2009;Alvioli et al., 2013).Such parallel codes are implemented with a variety of techniques and ideas.Given the design of the present work (C code wrapped into a GRASS GIS module written in Python, see Sect.2.2) they may act on three different levels of the overall implementation, namely (i) at the global GIS level, (ii) at the script level and (iii) at the level of the main code.
At the global GIS level, a few examples of integration of GRASS GIS and other GIS software tools into Grid, cloud or similar distributed computing environments are known (e.g.Huang et al., 2011;Aji et al., 2013;Agarwal et al., 2012;Wang et al., 2013).Nevertheless, the mentioned examples suffer from the problem that they are not of sufficient widespread use and the underlying implementation is not always reproducible since the corresponding middleware is not disclosed to the community.In other cases, pursuing such a strategy would require a higher-level expertise on the end user side, preventing potentially interested researchers from using the code.
On the other hand, GRASS GIS allows the use of script, a list of calls to the numerous software modules (written in Bash or Python) to be executed sequentially.Scripting can be exploited for simple and effective parallelization by properly adopting a strategy which is usually dictated by the particular problem one is trying to solve (see e.g.GRASS Wiki, 2014).
In the case of parallelization at the main code level, one can in principle choose to work within message-passing interface or OpenMP schemes, known to be extremely effective when applicable.As a matter of fact, no or very few working examples of such GRASS GIS parallel modules exist, leaving room for unpredictable inconveniences both at the stage of developing the code and at run time.
In the present work, we propose an implementation of a 3-D sliding surface model where the study area is partitioned into overlapping tiles, processed in parallel on a multi-core computer.For our experiment, we use the model r.slope.stability,a further development of r.rotstab (Mergili et al., 2014), to demonstrate strategies to optimize model performance in terms of computing time and of the quality of the model results.r.slope.stabilitywas implemented as a raster module of the open source software package GRASS GIS (Neteler and Mitasova, 2007;GRASS Development Team, 2014).GRASS GIS is well suitable for the task due to its open structure, modular design, and the compatibility with various programming languages.Further, GRASS GIS is frequently used as the basis for GIS-based models related to mass movements (Mergili et al., 2012a(Mergili et al., , b, 2014;;Gruber and Mergili, 2013).Our parallel implementation is performed at the Python level, whereas the core of the model is written in C. The model code and a user manual can be obtained from the model's web site http: //www.slopestability.org.The present work largely builds on r.slope.stability20140520 (20 May 2014), GRASS 6.4 and R 3.0.2.
In the following sections, we first introduce the 3-D (strictly speaking, 2.5-D, as the vertical dimension is represented by attributes, not by coordinates) slope stability model r.slope.stability(Sect.2).We then present the study area and the data used in the experiment (Sect.3), and we define the framework for testing the performance of the software in terms of quality of the model results and computing time (Sect.4).Next, we present (Sect.5) and discuss (Sect.6) the results before concluding with the key messages of the work (Sect.7).

Modelling approach
The r.slope.stabilitymodel evaluates the slope stability conditions for a large number of randomly selected ellipsoidal or truncated slip surfaces (Fig. 1).The ellipsoidal slip surfaces are defined by the geographic coordinates of the centre, the length of the three half-axes a e , b e and c e , the aspect α, the inclination β, and the offset of the ellipsoid centre above the terrain z b .The a e half-axis follows the steepest slope, and c e is aligned perpendicular to the terrain surface.a e , b e and c e are derived from landslide length L, landslide width W , maximum depth of the bottom of the ellipsoid D, β and z b (see Fig. 1b).Simple pseudo-random numbers, generated separately for each parameter of each ellipsoid, are used to define the centre coordinates as well as the values of L, W , D and z b , constrained by user-defined minima and maxima for each parameter.Testing a sufficiently large number of ellipsoids ensures a proper repartition of the ellipsoids over the study area, and the consideration of a large variety of possible ellipsoid dimensions.The tested slip surfaces correspond well to ideal ellipsoids only for reasonably small pixels in relation to the ellipsoid size (Mergili et al., 2014).When using larger pixels, the shapes of the tested slip surfaces represent systems of discrete plane surfaces strongly depending on the discretization of the pixels.For the modelling of more realistic shallow failure surfaces, r.slope.stabilitycan use truncated ellipsoids to consider the bottom of soil, shallow weak layers, or shallow discontinuities bounded by hard bedrock as possible failure surfaces.As a consequence, more than one slip surface may be associated with each ellipsoid (Mergili et al., 2014).To compute FoS, r.slope.stabilityadopts a modified version of the 3-D sliding surface model of Hovland (1977), revised and extended by Xie et al. (2003Xie et al. ( , 2004bXie et al. ( , c, 2006)): In Eq. (1), the upper term corresponds to the resisting forces R, and the lower term corresponds to the driving forces T .R and T are summed over all columns C of the slip surface.c (N m −2 ) is the effective cohesion, A (m 2 ) is the 3-D area of the slip surface of the considered pixel, G (N) is the weight of the moist soil, β c is the inclination of the slip surface at the considered column, ϕ is the effective internal friction angle, and β m is the apparent dip of the slip surface at the considered column in the direction of α.N s and T s (N) are the contributions of the seepage force to the normal force and the shear force.No inter-column forces or external forces, such as seismic loading, are considered by the model.The geotechnical, hydraulic, and geometric principles of the FoS calculation are discussed in detail by Mergili et al. (2014).Upon completion of the slope stability calculation for all the slip surfaces, each pixel in the modelling domain is intersected by various slip surfaces, and each slip surface is associated with a value of FoS.For each pixel, the lowest value of FoS of all the intersecting slip surfaces is taken as the representative FoS.
Compared to the r.rotstab model (Mergili et al., 2014), r.slope.stabilityintroduces the following innovations: 1.An improved data management strategy to meet the standards of GRASS GIS, including built-in functions for model validation and graphic presentation.
2. The ability to fully exploit multi-core computers.
3. The ability to consider complex systems of geological layers, relevant for the modelling of deep-seated landslides.
4. The ability to compute the slope failure probability P f in addition to FoS, based on the statistical distribution of c , ϕ and, for truncated ellipsoids, the truncated depth d.Four types of PDFs can be used: rectangular, normal, log-normal or exponential (see Sect. 4).
Points ( 1) and ( 2) are explained in Sect.2.2, and point ( 3) is not exploited in this work.The rationale of point ( 4) relies on the high natural variability of the geotechnical and geometric parameters, resulting in an uncertain definition of the horizontal and vertical distributions of c and ϕ (see Sect. 3).
A map of FoS building on data from a single site, or a limited number of sites, may fail to account for the details of the landscape.To overcome this limitation, we adopt an approach to compute the slope failure probability P f .This approach allows considering the full range of measured values of c and ϕ .The statistical properties of the parameters are assumed constant in space (see Sect. 6).A range of values of the truncated depth d can be considered, which is particularly useful for modelling shallow landslides in soils of uncertain depth.This approach is implemented in the following three steps: 1. Computing the arithmetic mean µ, standard deviation σ , minima and maxima of c , ϕ and d.The number of statistic samples n of parameter combinations to be considered is defined by the user.
2. c , ϕ and d are varied as a function of the defined minima, maxima and intervals in order to exploit the full range of possible parameter values.The variation of d builds on truncating the ellipsoid at various depths.FoS is computed for each combination using Eq. ( 1), building the ratio of the sums of the shear resistances and the shear forces over all columns of the ellipsoid (see Fig. 1).
3. The slope failure probability P f for the ellipsoid is computed as a function of the fraction of parameter combinations where FoS < 1, related to all the tested parameter combinations: where f i = 1 for FoS i < 1, f i = 0 for FoS i ≥ 1, and w i is the weight assigned to the parameter combination i (see below).The sum of w i over all parameter combinations n is 1.
At the end, the largest value of P f out of all intersecting slip surfaces is taken as the value representative for each pixel.
The sample of parameters to be tested has to represent the probability of occurrence of the parameter combination.In order to explore the influence of choosing different sampling strategies on the model results, we test three contrasting strategies.Figure 2 illustrates the strategies, using a sample size of n = 100 for two normally distributed, arbitrary parameters.Strategy (a) shows random sampling of parameter combinations: 100 parameter combinations are randomly sampled, where the probability of a parameter combination to be sampled directly relates to the product of the probability densities of the parameter values.Strategy (b) shows random sampling of parameters.Here, 10 values of each parameter are randomly sampled, where the probability of a parameter value to be sampled directly relates to its probability density.All possible pairs of sampled parameter values are then considered, resulting in 100 tests.Strategy (c) shows equal density sampling.Ten values of each parameter are sampled, equally distributed along the cumulative density function associated to each parameter (see Fig. 2d).This ensures that the distribution of samples reflects the PDF.All possible pairs of sampled parameter values are then considered, resulting in a total number of 100 tests.
Panels (a) and (b) represent the two possible ways to use a Monte Carlo approach -the most established strategywhilst (c) represents a deterministic approach.As the sample is fixed in the deterministic approach, it has to be determined only once at the beginning of the entire computation.The w i (see Eq. 2) represents the product of the cumulative density intervals associated with the values of the combined parameters (see Fig. 2d).This means that the edge samples are down-weighted as they only represent half of the area under the PDF, compared to the other samples.For the Monte Carlo approaches, the sample is determined separately for each ellipsoid tested, and w i = 1/n.
In the present work, this approach is applied to three (c , ϕ , d) instead of two (c , ϕ ) parameters.

Computational implementation
r.slope.stability is a raster module of the open source software package GRASS GIS 6.4 (Neteler and Mitasova, 2007;GRASS Development Team, 2014).The software exploits the Python programming language for data management, pre-processing and post-processing tasks.The slope stability model itself (see Sect. 2.1) is implemented as a C code (sub-module r.slope.stability.main).r.slope.stabilityalso includes a built-in validation and presentation module.Output maps and plots are produced using R, a free software environment for statistical computing and graphics (R Core Team, 2014).The logical framework of r.slope.stability is illustrated in Fig. 3.The numerical implementation presented in this work extends the applicability of the slope stability model to large study areas.This requires a very large number of ellipsoids to be tested.Assuming a test site with an area of A s = 100 km 2 , average ellipsoids of length L avg = 100 m and width W avg = 80 m, and an average number of ellipsoids per pixel (the "density" of ellipsoids) d e = 1000, the total number of ellipsoids n e to be tested sums up to roughly 16 million, The pixel spacing used for the slope stability model has to be small enough to capture the geometry of the assumed slope failure, which may fall into a very broad range of sizes (see e.g.Alvioli et al., 2014, and references therein).Given a study area of 100 km 2 and a pixel size of 5 m, four million pixels need to be processed.The potentially large number of pixels in combination with the large number of ellipsoids, and the complex processing of each ellipsoid, pose challenges in terms of (i) computer memory and (ii) computing time.We combine two strategies to overcome these computational challenges: 1.In the C programming environment, raster data sets are commonly held in memory as arrays.This allows a fast and efficient access to each pixel.If the data sets become too large, or if too many large arrays are held in memory at the same time, the available memory may be exceeded, causing the model execution to fail.We use the GRASS GIS Segment Library (GRASS Development Team, 2014) to avoid this problem.The library enables storage and use of very large raster data sets independently from the available computer memory, however at the expense of computing time.r.slope.stability.mainuses the GRASS Segment Library for data input, preparation, and output.For ellipsoid-specific computations, where a lot of data covering a smaller number of pixels have to be accessed frequently, it uses arrays by default.
In this study, we apply a segment size of 16 × 16 pixels to all computations, maintaining 16 segments in memory.As the most time-consuming operations of r.slope.stabilitymake only limited use of the GRASS Segment Library, preliminary studies have shown that, within a certain range, the computing time displays a weak dependence on changes in those settings.
2. To reduce the computing time when modelling the slope stability of large areas, r.slope.stabilityprovides the option to divide the study area into a user-defined number of tiles processed in parallel, if the code is run on an ordinary multi-processor or multi-core machine (see Fig. 3).In this case, r.slope.stability.main is run separately for each tile.The final result is obtained by collecting and combining the results for the single tiles.
To ensure a full coverage of the study area, an overlap between the tiles of at least the maximum ellipsoid dimension is required.Each tile is sent to a free computing core as soon as one is available, and until all the tiles are processed.This procedure is implemented in the way that the r.slope.stability.pymodule produces a batch file for each tile.The batch file calls the submodule r.slope.stability.multicore,which is then used to launch r.slope.stability.main with the tile-specific parameters (see Fig. 3); the actual parallel processing is performed in the Python part of the module, exploiting the "Threading" Python library (a higher-level threading interface) and the "Queue" Python module (a class for managing the "producer-consumer" problem able to block execution until all the items in the queue have been processed).
We note that neither the use of the GRASS Segment Library nor multi-core processing affect the model results in terms of FoS or P f .For our experiments we use a 48-core (AMD Opteron, frequency of 2.2 GHz and cache of 512 KB) computer with 140 GB of RAM and running a 12.04 LTS Ubuntu GNU/Linux OS with the 3.5.0-26-generickernel image.

Study area and data
We test the r.slope.stabilitycode in the Collazzone area, Umbria, central Italy (Fig. 4).Covering an area of 89.5 km 2 , this hilly area ranges from 145 m a.s.l.along the Tiber River flood plain, to 634 m a.s.l. at Monte di Grutti.Various types of continental sediments, Pliocene to Pleistocene in age, cover the area.Landslides are frequent and abundant in the Collazzone area, and a detailed landslide inventory (Fig. 4) is available along with geologic and morphologic information and maps (Guzzetti et al., 2006a(Guzzetti et al., , b, 2009;;Ardizzone et al., 2007;Galli et al., 2008;Rossi et al., 2010;Fiorucci et al., 2011).Intense or prolonged rainfall periods are the primary natural triggers of landslides in the area (Ardizzone et al., 2013), followed by rapid snow melt (Cardinali et al., 2000).Recent landslides are most frequent in cultivated areas, indicating a relationship with agricultural practices.
In the present work we focus on shallow landslides, considering an inventory of 2381 landslides (Fig. 4) for model evaluation (see Sect. 5).The 5th and 95th percentiles of landslide length L, width W , and of the L / W ratio for selected shallow landslides are used for constraining the randomization of possible slip ellipsoids (Table 1).As most landslides in the area are not very mobile, we use the entire landslide areas instead of the scarps only.
Most commonly, the sliding surface of shallow landslides coincides with the lower boundary of the soil which, in cultivated areas, we define as the layer disturbed by agricultural practices.Statistics of soil depth d s in the continental sediments of the Collazzone area were obtained from a set of 90 measurements, considering the lower boundary of the C v horizon, where present.Analysis of the measurements resulted in an arithmetic mean of the soil depth µ = 0.60 m, with a standard deviation σ = 0.27 m.The minimum soil depth measured in the area was zero, and the maximum soil depth was 1.22 m.
Figure 4 illustrates that landslides are rare where hard bedrock crops out (hatched areas), and abundant in the continental sediments (all other areas).In the present work, we consider all areas with hard bedrock outcrops as uncondi- 15.8 6.9 29.9 σ 1.0 7.2 * 7.5 * For the exponential distribution applied to cohesion, the standard deviation is set to the mean instead of using the value given in the table.
tionally stable, and concentrate to the areas where continental sediments crop out.The geotechnical characteristics of the continental sediments the study area were estimated using direct shear tests on 13 samples taken from a variety of lithological conditions (Table 2, see Fig. 4).The variation of geotechnical parameters within each class is considerable, partly exceeding the variation between the classes.For this reason, we decide not to consider separate sets of geotechnical parameters for the different lithological classes present in the study area.Instead, we explore the statistics of the parameters for the entire area with continental sediments.The same approach is used for the parameterization of the soil depth.
In addition to the landslide inventory, soil depths, and the geotechnical data, we use a 5 m × 5 m digital elevation model (DEM) derived by the automatic interpolation of 10 and 5 m contour lines, obtained from 1 : 10 000 scale topographic base maps.

Model parameterization
In this work, we consider only shallow slope stability, truncating the ellipsoids at the depth of the soil.We set the dry specific weight of the soil γ d = 15.8 kN m −3 (see Table 2), and the saturated water content s = 40 vol.%.Within a reasonable range of values, both parameters are not decisive for the outcome of the slope stability computation.Instead, FoS and P f are most sensitive to the effective cohesion c , the effective angle of internal friction ϕ , the depth of the potential  2).

As commonly observed for shallow landslides in the
Collazzone area, the maximum slip surface depth -at which all ellipsoids are truncated -is set to the soil depth.A log-normal PDF is used to model the variability of truncated depth.
4. We further assume the hydraulically most unfavourable case of fully saturated soil with slope-parallel seepage.
A separate map of FoS is computed, considering the most probable values (modes) of c , ϕ and truncated depth d, deriving those from their respective PDF: c = 0 kN m −2 ; ϕ = 27.3 • ; d = 0.46 m.Table 2 lists a range of c = 0-24.5kN m −2 and ϕ = 18.1-42.4• .These values are used to constrain the variation of the parameters during r.slope.stabilityruns.For c and ϕ , this range is justified by the rule-of-thumb values given by Prinz and Strauss (2011) for the possible range of geotechnical parameters for various soil types: for the continental sediments in the Collazzone area the relevant ranges would be c = 0-25 kN m 2 and ϕ = 15-45 • .For the soil depth d, the maximum of 1.22 m corresponds reasonably well to the approximately 1.3 m maximum depth of disturbance by agricultural practices observed in the Collazzone area (Mergili et al., 2014).Table 3 summarizes the parameters' minima, maxima, and assumed statistical distributions used for the computation of P f .The geotechnical parameterization is kept constant for all the tests.
Table 4 lists the parameters tested, and the settings applied in our numerical experiments.The ellipsoid size is constrained according to Table 1.The maximum depth of the bottom of the ellipsoid is constrained with D = 2.5-10 m.Considering all the combinations of the parameter values listed in Table 4 would result in a very large number of model runs, with excessive computing times.We therefore divide the task into two parts: 1. Multi-core processing: influence of multi-core processing on the computing time of r.slope.stabilityfor the entire Collazzone area (Fig. 4).A few combinations of the number of tiles t, and the number of processors p given in Table 4 are tested for d e = 100 and 2500, and dx = dy = 5 m.
2. Factor of safety and slope failure probability: influence of d e and -in the case of P f -sample size n (number of tested values of c , ϕ and d) and sampling strategy (see Fig. 2) on the model results (average value of P f and correspondence with observed shallow landslides).Part of this test is performed for a subset of the Collazzone area (see Fig. 4).All possible values of d e and n given in Table 4 are considered.

Test 1: multi-core processing
The gain in computing time due to parallel processing is most easily summarized by the speedup S p : where p is the number of processes, T 0 is the execution time of the sequential algorithm, T p is the execution time of the parallel algorithm with p processes, f s is the sequential fraction, summarizing the overhead, or irreducible serial part, of the code, and f p is the parallel fraction (f s + f p = 1).S p = p or f s = 1 − f p = 0 would indicate a linear (or ideal) speedup.
In such a case, the efficiency would be 1, with f s > 0 and E p < 1 in the case of r.slope.stabilitydue to (i) shared use of the RAM by multiple cores; (ii) non-optimized sequential use of cores; (iii) operations such as creating tiles and combining the results from each single tile.Further, the total area to be processed increases with t due to the overlap between the tiles.We now show the patterns of f s , S p and E p when using r.slope.stability to compute FoS for the entire Collazzone area at a pixel size of 5 m × 5 m, constraining the ellipsoid size according to Table 1. Figure 5 clearly illustrates that the values of f s , S p and E p depend on t, p and d e .Figure 5ac   of t (42).f s is much lower with d e = 2500, resulting in optimum values of S p and E p using a large number of tiles (see Fig. 5d-f).These observations are easily explained by the fact that speedup and efficiency are turned down with large values of t by the high cost of combining the results from the different tiles into one set of raster maps for the entire study area.The relative impact of this effect -and therefore also f s -decreases with increasing values of d e .With d e = 2500, the optimum speedup and efficiency are observed with 182 tiles.S p does not follow a linear increase with p, reflected in decreasing values of E p with p (see Fig. 5c and f).This phenomenon is most likely explained by the shared use of the RAM by multiple cores.
There is, of course, no gain in terms of speedup at p > t (not shown in Fig. 5).However, for the lower values of t, speedup becomes constant with increasing p already at p < t.This observation reflects a non-optimized sequential use of cores.Particularly with low values of p or t, and varying numbers of null cells among the tiles, it likely happens that one core is assigned much more work load than another.This type of effect, illustrated by the irregular patterns of E p For d e = 2500, the absolute values of T p reduce from 110 000 s for t = 1 and p = 1 to 4700 s with t = 182 and p = 42.However, the effects of considering other study areas, different pixel sizes or different ellipsoid dimensions on T p have to be noted.In principle, we expect a nearlinear dependency of T p on the number of pixels to be processed.However, increasing the pixel size results in an underproportional gain of T p .Areas of null cells due to the irregular shape of the study area cause computations on ellipsoids or entire tiles to break in an early stage of processing.This leads to a relative increase of operations not depending on the number of pixels.
We further expect that the computing time does not depend on the dimensions of the ellipsoids: a given value of d e means that all pixels of the study area have to be considered for approx.d e times (see Eq. 3).If larger ellipsoid dimensions are chosen, fewer ellipsoids need to be processed.However, larger ellipsoids have a higher chance to be cancelled as they touch areas with null cells rather than smaller ellipsoids.As a consequence, T p decreases with larger ellipsoids.In the specific setting considered here, doubling the constraints of L and W given in Table 3, resulting in a fourfold increase in size of an average ellipsoid, decreases the computing time by 11 % whilst executing the model with halved values of L and W , leading to a quarter of the original average ellipsoid size, increases the computing time by 21 %.

Factor of safety and slope failure probability
Next, we compute FoS for shallow landslides in the study area with values of d e = 100, 500, 2500 and 12 500.We eval-uate the modelling results against observed shallow landslide areas (see Fig. 4).Larger values of d e result in a more conservative prediction in terms of FoS -if more ellipsoids are tested, the chance is higher for each pixel that at least one slip surface is associated with FoS < 1 (Fig. 6).All tests result in a successful rather than unsuccessful prediction, even though the false prediction rates are significant.There is no optimum value for d e , per se.Strictly speaking, d e ∼ ∞ would be needed -as the rate of positive predictions may increase also at very high values of d e , there will always be a tradeoff between the computing time and the quality of the results.However, we note that, in this example, the overall quality of the prediction does not increase with larger values of d e (i.e. the polygon does not significantly shift towards a successful prediction), indicating that most areas with FoS < 1 were detected at earlier stages of the computation, and the additional areas with FoS < 1, detected at later stages of the computation, consist equally in true positive and false positive predictions.For the purpose of the present study, we consider d e = 2500 a sufficiently reasonable approximation.
We compute slope failure probability for a subset of the Collazzone area (see Fig. 4) with five different sample sizes, applying each of the sampling strategies (a), (b) and (c) introduced in Fig. 2. The c , ϕ and d are sampled.We assume that the accuracy of the results increases with increasing values of d e and n.However, so does the computing time.Therefore, we attempt to identify those values where the results converge -i.e. the ideal values in terms of accuracy and time efficiency.We take the average value of P f over the study area as reference.Figure 7a illustrates how the average value of P f increases with increasing d e .It further indicates the sample size n needed for convergence, i.e. the value of n where the average P f remains constant when n is further increased.Equal density sampling (c) performs best whilst random sam- pling of parameters (b) is not a valid alternative: with a very high number of tested ellipsoids, it is likely that at least one of the random samples is biased towards low values of c and ϕ .Therefore, on the logarithmic scale used in Fig. 7a, average P f steadily increases with increasing d e .This effect is less pronounced for larger values of n, but it could only be diminished by testing excessively large samples, i.e. at the cost of a very long computing time.Sampling strategy (a), with randomly sampled parameter combinations, is less susceptible to these effects as the samples are better distributed within their range (see Fig. 2).Still, with the assumptions tested, the curves converge at a higher average of P f and flatten out more slowly than the curves for equal density sampling.Further, strategy (a) is highly inefficient.With similar values of d e and n, the computing time is roughly 20-25 times longer than for strategy (c).The reason for this phenomenon is that the number of truncated depths to be tested is n with strategy (a) and the cubic root of n with the other strategies.Hence, with (a), the geometry of a much larger number of slip surfaces has to be built than for (b) and (c), which is costly in terms of computing time, given the current implementation of r.slope.stability.
Independently of the sampling strategy, the average slope failure probability decreases with the number of samples.For the strategies (a) and (b), this is a result of the lower tendency of outliers with larger sample sizes, which is more pronounced with (b) than with (a).With sampling strategy (c), it is a result of the fact that an exponential PDF is assumed for c .With lower sample sizes, the relative weight of the minimum value c = 0 kN m −2 is higher than with higher values of n, resulting in higher values of P f .Among all the tests shown, we expect (c) equal density sampling with n = 15 3 to perform best in terms of accuracy.All results shown in Figs.7b and  Figure 8a illustrates the modelled distribution of FoS in the study area, and Fig. 8b portrays the spatial patterns of P f .Table 5 summarizes the evaluation outcomes and computation times for FoS and P f for all the parameter combinations considered in Figs. 6 and 7b.

Discussion
Exploiting multi-processor computing environments enables the execution of complex slope stability models for reason-   to the same slope units as those used by Rossi et al. (2010), the value of A ROC increases to 0.77, indicating that, despite the uncertain geotechnical and geometric parameters, the physically based model slightly outperforms the statistical ones.However, (i) the results are not fully comparable due to different sets of landslides used for validation and (ii) a better understanding of the spatial distribution of the geotechnical and geometric parameterization is needed to improve the reliability of landslide susceptibility and hazard maps.
In the present work, we assume constant statistical properties (µ, σ , minimum, maximum) of the geotechnical parameters c and ϕ and of the soil depth d over the relevant part of the study area.Even though, in this specific case, we can well justify this generalization, it may be too simplistic in other cases.The ability to better constrain the geotechnical parameters would possibly also allow to reduce the size of the statistical sample and therefore the computing time.
We further assume independent statistical properties of c and ϕ .However, this is a rough simplification as these parameters -representing the offset and inclination of the linear regression in the Mohr-Coulomb diagram -are often negatively correlated.A future challenge will consist in finding an appropriate way to build PDFs considering the interdependency of the two parameters.
A further limitation consists in the assumption of fully saturated soil with slope-parallel seepage.The computed values of P f are therefore only valid for this worst-case assumption in terms of slope hydraulics.Partial saturation is more difficult to treat from a geotechnical point of view and shall be the subject of future studies.
Finally, the PDFs that were used in the study may be improved.Whilst the density functions for d and ϕ are reasonably well supported by the empirical observations, the exponential PDF used for c was derived for soils with a high content of sand and silt (El-Ramly et al., 2005;Petrovic, 2008).For clay, a log-normal function seems to better describe the observations.A joint, two-variable PDF depending on both c and ϕ may by hypothesized.Such a function is expected to yield significantly less conservative results.Given a sufficiently large data set, we suggest using the PDF for ϕ and couple the function for c to the tested value of ϕ .An appropriate geotechnical parameterization requires a detailed knowledge of the area under investigation.As an example, if deep-seated slope stability is considered, this understanding should include the strike and dip directions of bedding surfaces (Santangelo et al., 2014).

Conclusions
We have described and tested r.slope.stability,a multi-core numerical GRASS GIS implementation of a 3-D slope stability model for large areas, highlighting (i) the gain in computing time, and the consequent applicability to large areas, and (ii) the possibility of modelling the spatial probability of slope failures, based on the natural variability of geotechnical characteristics of the soils.Using commonly available multi-core hardware, the use of parallel processing may reduce running times by a factor larger than 20.Our parallel implementation is transparent to the r.slope.stabilityuser in GRASS GIS, since it is based on the automatic partitioning of the study area in tiles, processed in parallel.The modelling results are presented for the entire area, and validated against observed landslides.
We conclude that parallel processing enables the application of complex slope stability models for large areas in a reasonable amount of time.A remaining challenge for this type of task is the geotechnical parameterization of the area under investigation.In the present paper, we have demonstrated a simple approach to compute slope failure probabilities by using PDFs of c , ϕ , and d.This approach is considered sufficient for the purpose of the present work.The model results reasonably correspond to the distribution of observed shallow landslides in the Collazzone area.However, we have identified a considerable potential for improvement with regard to (i) regionalization of the parameters, (ii) consideration of the interrelation of c and ϕ and (iii) optimization of the PDFs used.

Figure 1 .
Figure 1.Typical ellipsoid used as slip surface in r.slope.stability:(a) ground plot, (b) longitudinal section, (c) forces acting at each column.The factor of safety is computed for the ellipsoid bottom and (as shown in the figure) for the combination of the ellipsoid bottom and each intersecting layer bottom.The geotechnical, hydraulic and geometric details are outlined by Mergili et al. (2014).

Figure 2 .
Figure 2. Sampling of parameters for computing slope failure probability.We test the sampling strategies (a), (b), and (c).For clarity, sampling of two normally distributed arbitrary parameters is shown.In reality, we sample three parameters (c , ϕ and truncated depth) according to different types of statistical distributions (see Sect. 4 for details).Panel (d) illustrates how the cumulative density function is employed for equal density sampling.

Figure 3 .
Figure 3. Logical framework of r.slope.stability.Plain text denotes steps directly implemented in the module r.slope.stability,and text in boxes denotes sub-modules.Italic letters indicate the programming environment used for the modules.

Figure 4 .
Figure 4. Collazzone study area, Umbria, central Italy.IDs of sample points correspond to IDs listed in Table2.

Figure 5 .
Figure 5.The serial fraction of code f s , speedup S p and efficiency E p , plotted against the number of processes p for different values of the number of tiles t and average number of ellipsoids per pixel, d e .See text for further explanations.

Figure 6 .
Figure 6.Influence of d e on the performance of r.slope.stability in terms of prediction rates, building on FoS.Pixels representing observed shallow landslide areas (observed positives, OP) with a modelled value of FoS < 1 represent true positive predictions (TP).OP pixels with FoS ≥ 1 represent false negative predictions (FN).Pixels not representing observed shallow landslide areas (observed negatives, ON) with a modelled value of FoS < 1 represent false positive predictions (FP).Finally, ON pixels with FoS ≥ 1 represent true negative predictions (TN).

Figure 7 .
Figure 7. Slope failure probability computed with r.slope.stability.(a) Evolution of P f for a subset of the Collazzone study area (see Fig. 4) with increasing value of d e .The outcomes of the sampling strategies (a), (b) and (c) introduced in Fig. 2 are compared.(b) ROC plot relating P f obtained with equal density sampling (strategy c) to the observed shallow landslide areas in the entire Collazzone study area for different values of n and d e .
8 therefore build on sampling strategy (c).With d e = 12 500, n = 15 3 yields an average P f = 0.094.n = 12 3 and d e = 12 500 yields an average P f = 0.095.In a certain range, reducing n and d e affects moderately the model results, but improves significantly the computational efficiency.Setting n = 9 3 and d e = 12 500 gives an average P f = 0.097, saving 75 % of the computing time, and setting n = 9 3 and d e = 2500 gives an average P f = 0.090, saving 95 % of the computing time.Given the level of uncertainty in the geotechnical parameterization, reducing the values of n and d e can be a strategy for the computation of very large areas, keeping the computing time within reasonable limits.Even though we do not recommend using values of n < 9 3 and d e < 2500, Fig. 7b shows that, within a certain range, changes of n and d e do not affect significantly the capability of the model to reproduce the patterns of observed landslide/non-landslide areas in terms of the area under the ROC curve A ROC .This indicates that changes of the results for larger values of n and d e affect equally areas with low and high values of P f .

Figure 8 .
Figure 8. Spatial patterns of shallow slope stability in the Collazzone study area, computed with r.slope.stabilityapplying equal density sampling (strategy c).(a) FoS for d e = 2500.(b) P f for d e = 2500 and n = 9 3 .

Table 1 .
5th and 95th percentiles of length L, width W , and the L / W ratio for selected shallow landslides mapped in the Collazzone area.L is measured in the direction of the steepest slope.
Percentile L (m) W (m) L / W

Table 2 .
Geotechnical key parameters derived for 13 samples from the Collazzone study area (see Fig.4for location of the sites).γ d = dry specific weight (kN m −3 ), c = effective cohesion (kN m −2 ), ϕ = effective angle of internal friction (degree).Arithmetic mean µ and standard deviation σ are listed.

Table 3 .
Constraints and assumed statistical distribution of geotechnical parameters and soil depth for the generation of a slope failure probability P f map.c = effective cohesion (kN m −2 ), ϕ = effective angle of internal friction (degree), d = soil depth (m).
illustrate f s , S p and E p for d e = 100.The graphs clearly reflect high values of f s for high values of t.Consequently, speedup and efficiency are highest with relatively low values

Table 4 .
Parameters tested for their influence on model performance.Subscripts x and y refer to the x and y directions.
at lower values of t, is smoothed out at high values of t, where load balance is roughly done automatically.This phenomenon also results in increasing values of S p for t > p.

Table 5 .
Characterization of model results and computing times.Prediction rates: TP = true positive, TN = true negative, FP = false positive, FN = false negative.