Introduction

As a slope-forming process, landslides depict a natural part of landscape evolution that can be triggered by various agents, mainly rainfall and earthquakes (Glade et al. 2005). The size, shape, intensity and predisposition of the triggered landslides, however, are mainly determined by (un)certain environmental factors, such as slope geometry, composition of the soil or regolith and the vegetation cover (Guzzetti et al. 1999; Glade et al. 2005; Reichenbach et al. 2014). Particularly forest cover—and thus the distribution and density of trees—underlies a continuous spatial change in populated mountainous woodland areas all around the world. Hence, the susceptibility of a landscape to be affected by landslides is highly associated to the spatio-temporal change of the forest cover (Papathoma-Köhle et al. 2013; Promper et al. 2014). Landslide susceptibility is defined as the spatial probability of landsliding and can be approximated using various approaches, e.g. statistical assessments or physically based models (Guzzetti et al. 1999; Van Westen 2000; Van Westen et al. 2006). Latter are based on the limit equilibrium concept and presume that slopes composed of rigid soils can fail along a single failure plane (slip surface) according to the Mohr-Coulomb failure criterion. A dimensionless Factor of Safety (FoS) is calculated that measures the relation between the resisting forces (R) and the driving forces (T) that stabilise or destabilise the soil package on a slope:

$$ \mathrm{FoS}=\frac{\mathrm{R}}{\mathrm{T}}. $$
(1)

Herein, FoS > 1 depicts stable conditions, whereas FoS ≤ 1 indicates unstable slope conditions.

Slope stability can be computed in one, two or three dimensions (Xie et al. 2003). In general, a 1D slope stability model only considers soil thickness as a linear, metric parameter to compute FoS for an individual pixel (Van Westen et al. 1997). In geotechnics, 2D models are commonly used to assess the mechanical stability along a pre-defined failure plane of the vertical cross-section of a slope. However, both 1D and 2D models are not able to represent the three-dimensionality of slip surfaces. GIS environments are generally defined as 2.5D environments, since only a single z value per GIS raster layer (i.e. a set of discrete z values in the case that more than one layer is used) is assigned to a plane coordinate (x- and y-tuple) of a raster pixel in a Digital Terrain Model (DTM). This distinguishes 2.5D from 3D environments, in which a location is represented by an isotropic volume element (voxel). Thus, a location with the plane coordinates x and y might contain multiple discrete z values in a 2.5D environment, whereas a continuous representation of the patterns in z direction requires a 3D environment. The expression ‘3D’ or ‘three-dimensionality’ is often misleading in literature, due to the discriminative practices of the expression in different communities (e.g. geotechnical engineering and geomorphology) or simply because of inconsistent (or negligent) use. Hence, it is not always discernible whether a 2.5D or 3D was used.

Embedded in a GIS environment (2.5D), the infinite slope model evaluates the stability of a terrain for individual raster cells without taking into account forces that are apparent in adjacent raster cells (e.g. Montgomery and Dietrich 1994; Pack et al. 1998; Baum et al. 2002; Malet et al. 2005). However, this model is not or only conditionally suitable for slip surface length-to-depth ratios L/D < 16–25 (Griffiths et al. 2011; Milledge et al. 2012), or for soils with discontinuities such as root systems. In such conditions, slip surface models are more appropriate, applying more or less complex limit equilibrium relations to calculate the FoS along pre-defined or randomly determined failure planes (e.g. Hovland 1979; Hungr 1987, 1988; Hungr et al. 1989; Lam and Fredlund 1993; Xie et al. 2004). Several software applications employing 2.5D or 3D slip surface models are available, such as CLARA (Hungr 1988), TSLOPE3 (Pyke 1991), 3D-SLOPE (Lam and Fredlund 1993) and the r.slope.stability (Mergili et al. 2014a). 2.5D and 3D slip surface models have been applied either for model testing (e.g. Hovland 1979; Hungr 1988; Lam and Fredlund 1993; Xie et al. 2003; Mergili et al. 2014a) or for particular case studies (Seed et al. 1990; Chen et al. 2001; Chen et al. 2003; Jia et al. 2012; Mergili et al. 2014b).

The reinforcement potential of plant roots on the shallow subsurface and its importance as a model input parameter in physically based slope stability models was addressed in many previous studies (Sidle et al. 1985; Greenway 1987; Sidle 1991; Sidle and Ochiai 2006; Stokes et al. 2008; Fan and Chen 2010; Ghestem et al. 2011; Schwarz et al. 2012). Tree roots increase the soil strength during shear due to certain root characteristics, such as areal density, root distortion and root tensile strength (Wu et al. 1979). Moreover, the additional cohesion provided by roots is strongly connected to the distribution of the roots, determining the root system architecture (Sidle et al. 1985; Greenway 1987; Sidle and Ochiai 2006; Stokes et al. 2008; Ghestem et al. 2011). Multiple investigations focused on the distribution and architecture of root systems (Schmidt et al. 2001; Schmid and Kazda 2002; Puhe 2003; Bischetti et al. 2005; Danjon et al. 2008; Bischetti et al. 2009) and the differences in root tensile strength and cohesion (Schmidt et al. 2001; Preti and Giadrossich 2009). Cohesion forces of forest root biomass that provide additional stability to the soil are often modelled in 2D (Temgoua et al. 2016). Several modelling approaches were developed to appropriately integrate root cohesion in deterministic models in 1D or 2D environments (e.g. Sidle 1991; Hammond et al. 1992; Van Beek et al. 2007; Kokutse et al. 2016). However, tree roots usually do not form a continuous network throughout their maximum rooting depth and might penetrate differently composed soil layers or slip surfaces. Thus, in order to simulate the additional stability to hillslopes by roots in an adequate way, the 3D spatial distribution of root systems has to be considered. Various approaches were developed to model the 3D root system distribution and 3D root cohesion, such as fibre or root bundle models (e.g. Thomas and Pollen-Bankhead 2010; Schwarz et al. 2010), 3D finite element models (e.g. Dupuy et al. 2007; Temgoua et al. 2016) or within a infinite slope environment (e.g. Bischetti et al. 2005; Danjon et al. 2008; Bischetti et al. 2009). So far, however, information about forest cover and stand density in general, and mechanical root reinforcement by trees in particular, have been disregarded in GIS-based slip surface models. Appropriate strategies to parameterise the spatial distribution of the root systems in relation to the tested slip surfaces still have to be developed. In the present study, we explore this gap. For this purpose, we construct a generic hillslope topography, structured into eight different plots. On each of these plots, a set of scenarios are simulated that comprise different forest stand conditions by means of root system composition (uniform or mixed root systems with specific depth and radius), stand density (sparsely or densely distributed trees) and stand distribution (patchiness of forest stands over the plot). The assumptions made upon the spatial extension of the approximated root systems in these scenarios are derived from characteristic root systems found in central European mountain forests (Kutschera and Lichtenegger 2002). The spatial patterns of slope stability within each of the zones are analysed with the 2.5D slip surface model of the tool r.slope.stability, developed by Mergili et al. (2014a), in order:

  • to allow for quantitative statements on the influence of different types and densities of root systems on the slope stability conditions;

  • to build a basis for developing strategies towards appropriately parameterising root reinforcement in real-world slope stability analyses in a GIS environment.

Our study is based on three main assumptions:

  1. i.

    the stabilising potential of root systems are diverse according to both root architecture and spatial position along a hillslope;

  2. ii.

    potential slip surface size influences whether a distinct root system type (or a tree stand with a diverse distribution of various root systems) is able to contribute to slope stabilisation.

  3. iii.

    the approximation of root systems by paraboloid solids allows a straight-forward parameterisation of the spatial root reinforcement distribution in 2.5D slip surface models.

Next, we explain the generic landscape and root systems used for all the analyses. We then introduce the relevant functionalities of the software r.slope.stability and the computational experiments performed (“Materials and methods”). We present (“Results”) and discuss the results and conclude with the key messages and an outlook (“Discussion”).

Materials and methods

Generic landscape and root system morphology

We create a generic hillslope landscape that ensures similar topographic conditions for all root system assumptions. The characteristics of the generic landscape are summarised as follows:

  • Three identically shaped octagonal pyramidoids (Fig. 1);

  • Each pyramidoid has an extent of 2500 m × 2500 m, 2 m × 2 m raster cell grid and is divided into eight equal plots of triangular shape, so that 24 identical plots are available in total. The surface area of each plot is 51.75 ha.

  • All plots have an equal slope inclination of 25°.

  • The slope stability analyses are performed only for the lower 75% of the pyramidoid since the convergence of the margins in the upper 25% would not allow a proper analysis of the r.slope.stability results.

  • Only the mechanical reinforcing effect of roots is considered; hydrological effects of roots are not considered.

Fig. 1
figure 1

Pyramidoid of the generic landscape. Only the shaded area is considered for modelling. The entire landscape consists of three pyramidoids identical in shape

Full size image

We further define three sets of root system scenarios, each associated to one of the pyramidoids: set I assumes uniform tree distribution; set II assumes separated tree stands (patches); and set III assumes mixed patches (Fig. 2). Each scenario consists of eight plots that vary among themselves in stand density (S) and type of the root system. We distinguished two root system types: shallow root systems (SRS) and taproot systems (TRS). Mixed stands are composed by a mixture of both (MRS). Two plots of scenario I are kept without any tree coverage as a reference. Hence, 23 different types of root systems are tested (Table 1; see Fig. 2).

Fig. 2
figure 2

Root system scenarios imposed on the generic landscape. The detail figures D1–D3 exemplify the peak area of the landscape. We note that the patterns shown within each plot are representative for the entire plot. Shaded raster cells depict the simulated position of a tree stem for a distinct root system (blue: SRS; black: TRS)

Full size image
Table 1 Land cover conditions and properties of distinguished plots
Full size table

We simulate sparse and dense stand conditions (S = 50, 150 and 2000 trees ha−1, respectively) with different maximum distances between the tree trunks (cf. Schmid and Kazda 2002; Puhe 2003). The overall distribution of trees on each plot are categorised in ‘uniform’ coverage (scenario I), ‘separated patches’ (scenario II), ‘uniform patches’ and ‘mixed patches’ (both scenario III). Each raster cell of a plot (2 m × 2 m) can only accommodate a single tree, which fixes minimum tree distance to 2 m (raster cell length). For a uniform coverage, we assume a full forest cover of the plot and tree distances depending on the simulated stand density condition (S). We create randomly generated values for each raster cell ranging from 0 to 1 and define those raster cells as tree locations, where the random raster value does not exceed the percentile of the tree cells that can occupy the plot. For example, for a uniform tree density of 50 trees ha−1, all raster cells are assigned as tree locations, Ntree (Ntree ≔ N0), where the random cell value ≤ 50 trees ha−1 × 51.75 ha (plot size). For mixed uniform conditions (e.g. S1.7 and S1.8), S halves for each species. Thus, stand density S for both species remains equal, whereas only 50% of trees of each species can occupy a stand. In case of scenario II and III, where forest is simulated to not cover the whole plot, we define patches (P) first that function as tree stands. Therefore, 25–30 or 35–40 random points over the whole plot are selected from which radially ~ 2500 adjacent raster cells define the patch area, where trees can be located. All cells that are defined as tree cells in a uniform coverage and are located within the defined patches are then assigned as tree cells. Raster cells that are located outside the patches are assigned as non-tree cells. Herein, the spaces between patches create voids that simulate forest glades in which potential reinforcement from the tree root systems is remarkably reduced or non-existent. In cases where patches within a single plot shall be separately occupied by SRS and TRS (scenario III, plots 1–4), the number of patches P is again halved for each species, whereas S remains constant. For cases of mixed conditions within the patches (scenario III, plots 5–8), P remains constant but S halves (cf. description above for S1.7 and S1.8).

We approximate the vertical cross section of the root system of any given tree by a parabola, which is determined by a pre-defined radius rmax and depth zmax that spreads from the stem:

$$ {z}_{ij}=\frac{z_{\mathrm{max}}}{{r_{\mathrm{max}}}^2{r}_{ij}}, $$
(2)

where zij is any value at vertical position i and horizontal position j, and rij is the lateral distance between the stem and rmax at z0,j. Considering that rmax is deemed to be equal all around the stem, the spread of the parabola forms a paraboloid, approximating the basal surface of the root system in the 2.5D environment (Fig. 3).

Fig. 3
figure 3

Root system approximation with a paraboloid. a Conceptual illustration of approximated rmax and zmax. b Modelled root system paraboloid in a 2.5D mesh (each node represents the edge point of a raster cell)

Full size image

The dimensions of zmax and rmax are determined by the root system type and define the shape of the paraboloid. Environmental factors that normally influence the distribution of the roots (Sidle and Ochiai 2006; Ghestem et al. 2011; Schwarz et al. 2012) are not considered. According to findings from in-field root excavations (e.g. Smith 1964; Stone and Kalisz 1991; Kutschera and Lichtenegger 2002; Bolte et al. 2003, Kalliokoski et al. 2008) zmax and rmax are kept constant for the two root system types, assuming adult trees with the same age (Table 2). Shallow (SRS) and taproot systems (TRS) are approximated by the paraboloids. Mixed root systems (MRS) represent a combination of SRS and TRS. Root cohesion was then defined to be uniform for the whole root system and similar for all root system types. A reduction of root cohesion by inter tree competition was not considered, since we intended to focus on the effects of the root system composition and the location of the respective root systems in respect to the ellipsoid.

Table 2 Constant values for vertical and horizontal distribution of the two considered types of root systems
Full size table

The tool r.slope.stability

The 2.5D slip surface model r.slope.stability, developed by Mergili et al. (2014a,b), evaluates the slope stability conditions in a given area based on predefined or randomly distributed ellipsoid-shaped or truncated slip surfaces. The slip surface dimensions in terms of length (L), width (W) and maximum depth (D) (all in m) of the tested slip surfaces can be defined either as fixed values or determined randomly. Not only the ellipsoid bottom, but optionally also user-defined layer interfaces intersecting the ellipsoid surface are considered as potential slip surfaces. Relying on a DTM and a set of geotechnical parameters, the r.slope.stability computes a FoS that is based on the relation of Hovland (1979), revised by Xie et al. (2003, 2004, 2006):

$$ \mathrm{FoS}=\frac{\sum \left({c}^{\hbox{'}}A+\left({G}^{\hbox{'}}\cos {\beta}_{\mathrm{c}}+{N}_{\mathrm{s}}\right)\tan {\varphi}^{\hbox{'}}\right)\cos {\beta}_{\mathrm{m}}}{\sum \left({G}^{\hbox{'}}\cos {\beta}_{\mathrm{m}}+{T}_{\mathrm{s}}\right)\cos {\beta}_{\mathrm{m}}}, $$
(3)

where c′ (N m−2) is the apparent cohesion (soil plus roots), A (m2) is the slip surface area at the considered raster cell, G′ (N) is the weight of the soil, βc (°) is the inclination of the slip surface, βm (°) is the dip of the slip surface in the column with the aspect α, and φ′ is the effective internal friction angle. Ns and Ts (both in N) are contributors of the seepage force to the normal and the shear force, respectively (Fig. 4).

Fig. 4
figure 4

Example of a slip surface approximated as ellipsoid by the r.slope.stability (modified after Mergili et al. 2014b). Depth and length are given by D and L. a Ground plot of the ellipsoid in a raster-based GIS environment. The boundary of the ellipsoid is indicated by the green dashed line, the aspect of the ellipsoid by α. b Longitudinal section of the ellipsoid. The lower boundary (green dashed line) indicates the slip surface with the inclination β. Computed forces are given in the lower right box

Full size image

To take into account the uncertainty of the input parameters, the r.slope.stability includes the option to consider a user-defined space of c′ and φ′ values instead of fixed values (Mergili et al. 2014a). Following a probability density function (PDF), FoS is computed for a large number of parameter combinations within the given space. Assuming a uniform PDF, the slope failure probability Pf for a given slip surface essentially corresponds to the fraction of parameter combinations yielding FoS < 1. Consequently, Pf can assume values in the range 0–1.

The value of FoS or Pf for a given raster cell results from the overlay of a large number of slip surfaces: the lowest value of FoS or the highest value of Pf out of all intersecting slip surfaces is considered most relevant and therefore applied as the final result. We refer to Mergili et al. (2014a,b) for the technical and mathematical details of the r.slope.stability.

Computational experiments

We apply the r.slope.stability to the generic landscape in order to compute the spatial distribution of Pf for the different root system scenarios. Thereby, three computational experiments are performed in terms of varying the spatial dimensions of the slip surfaces (Table 3). Within the experiments E1 and E2, the slip surface dimensions are fixed, whereas the size of each individual slip surface is determined randomly within E3. In all of the experiments E1–E3, W is set to two thirds of L. An L/W ratio of 1.5 depicts an acceptable compromise for landslides with rotational character, as reported in previous case studies (Parise and Jibson 2000; Dewitte and Demoulin 2005; Cislaghi et al. 2017). D is set to 2 m in order to represent the shallow subsurface of the soil where vegetation has its major biomechanical impact on the soil (Greenway 1987; Sidle and Ochiai 2006; Stokes et al. 2008). This also means that SRS only penetrate the edges of a given potential slip surface, and not in the central area where D > zmax(SRS). The slip surface density (Mergili et al. 2014b) is set to 200 for the experiments E1 and E2, and to 1000 for E3. Moreover, the critical length of the ellipsoid or slip surface (Lcrit), respectively, was computed for all scenarios. Lcrit is the length of the ellipsoid for which the lowest value of FoS or the highest value of Pf was calculated for a certain pixel. That is, Lcrit is calculated from the results for all the ellipsoids in a plot that intersect the respective pixel under consideration.

Table 3 Spatial dimensions of the slip surfaces in each computational experiment
Full size table

The geotechnical parameterisation remains unchanged throughout all computational experiments and for all root system scenarios (Table 4). The parameter space applied is supposed to cover a broad variety of conditions typical for central European mountain forests and the associated soils. c′ is defined as the sum of soil cohesion c′s and root cohesion c′r. The range between minimum and maximum of c′r was set according to reported values (4.7–94.3 kN m−2 for coniferous and 5–100 + kN m−2 for deciduous trees, depending on location, stand composition, distance from trunk and rooting depth) in former studies (cf. Schmidt et al. 2001; Roering et al. 2003; Sakals and Sidle 2004; Bischetti et al. 2009). Pf is computed for each slip surface from 36 combinations of c′ and φ′. Herein, we assume a rectangular (uniform) probability density function for both parameters. More precisely, the values of each parameter are uniformly distributed within the parameter range defined in Table 4 and analysed for six parameter values each, covering the range between the respective minimum and maximum of the parameter at equal intervals. Further, Pf is computed as a fraction for all tested parameter combinations where FoS < 1 at the boundary of the ellipsoid (congruent to the slip surface). Finally, the intersections of all slip surfaces with the largest value of Pf are taken as representative Pf value for the respective raster pixel. By selecting a set of samples within a parameter value range and computing Pf with different parameter combinations, we account for a probabilistic parameter representation, which enables to cope with uncertainties of the required geotechnical input parameters in the modelling process. However, different parameter sampling strategies can be applied, e.g. when some ground-truth information are available. Therefore, we refer the reader to Mergili et al. (2014a) for a concise synopsis of parameter sampling strategies in the r.slope.stability. The specific weight of the dry soil γ is kept constant. We further assume fully saturated conditions with a saturated water content θs = 40% per volume.

Table 4 Values of geotechnical parameters used in the r.slope.stability
Full size table

Results

Slip surface size

The spatial patterns of Lcrit for all scenarios are illustrated in Fig. 5. Note that the values of Lcrit are considerably higher along the edges as (i) they are influenced by the steeper slope of the edge, compared to the sides of the solid and (ii) they may be influenced by the pattern in the adjacent plots. This effect increases with slip surface size. Therefore, all analyses of the results disregard those areas adjacent to the edges. With roots, an increase in the slip surface size generally results (i) in blurred patterns (Fig. 5) and (ii) in lower values of Pf. (i) is explained by the increasing area covered by each individual—and therefore also the most critical—slip surface. (ii) is related to the fact that smaller slip surfaces are more likely to squeeze in the spaces between the individual root systems, and therefore escape the influence of root cohesion. However, this only applies with low values of S or with separated patches: otherwise the spaces within the individual root systems are too small even for E1 slip surfaces, so that the results yielded with E1 and E2 are almost identical.

Fig. 5
figure 5

Distribution of the critical slip surface length (m) (Lcrit) for all scenarios in E3

Full size image

Considering the densities of the uniform stands (scenarios S1.3–S1.6), TRS stands show a remarkably higher decrease of Lcrit than SRS when stands are sparse (Fig. 6). Allowing for a variable slip, surface size (E3) combines the patterns observed with E1 and E2 (Fig. 5). Whereas local minima of Pf derived for E1 are smoothed out (e.g. scenarios S1.5 and S1.7 in Fig. 7) by large slip surfaces, local maxima are associated to small slip surfaces. In general, E3 is supposed to yield higher values of Pf for all raster cells.

Fig. 6
figure 6

Critical ellipsoid length (Lcrit) at which the lowest FoS values (or highest Pf values, respectively) for raster cells with high and low stand densities of SRS and TRS are computed in E3. Lcrit decreases with sparse stand conditions—particularly for stands with TRS

Full size image
Fig. 7
figure 7

Distribution of Pf for all scenarios and computational experiments

Full size image

Scenario I: uniform tree distributions

We now analyse the variations in Pf for scenario I (Fig. 8), considering E3. Pf is generally highest in the plots without roots, and lowest in the plots with mixed roots or taproots, whilst it is intermediate in the plots with shallow roots. However, the influence of the stand density is more pronounced than the influence of the type of the root system. With E1 and S = 50 trees ha−1, Pf decreases from the initial value of 0.84 (no trees) to 0.78 (SRS), 0.74 (MRS) and 0.70 (TRS). With S = 2000 trees ha−1, Pf decreases to 0.13 (SRS), and even 0.00 (MRS and TRS). Also, the standard deviations (SDs) are higher for the scenarios with low stand densities, reflecting the more pronounced spatial pattern in Pf associated with a decrease in stand densities (ability of small slip surfaces to squeeze in; see “Scenario I: uniform tree distributions”). Larger slip surfaces (E2) lead to equal or slightly higher averages (M) of Pf than E1 for all scenarios except for those with S = 50 trees ha−1 (see “Slip surface size” and Fig. 8). Due to smoothing of the patterns, SDs are much lower than those associated with E1. Variable slip surfaces result in higher M and SD of Pf.

Fig. 8
figure 8

Pf for all plots of scenario I. The colours of the dots indicate the Pf value. The circles around the dots indicate the standard deviation (SD)

Full size image

Scenario II: separated patches

Considering separated patches, the stability conditions are determined by an interplay between the type of the root system (SRS vs. TRS), the stand density and the patch density. Those slopes with TRS are more stable than those slopes with SRS. Equally, increasing stand and patch densities lead to decreasing values of Pf (Fig. 9).

Fig. 9
figure 9

Pf for all plots of scenario II (description of colours, dots and circles cf. Fig. 8)

Full size image

Whilst E1 and E2 yield comparable M of Pf for S = 150 and P = 25–30, we note that for TRS, increasing the patch density has a comparatively larger effect on the average Pf with E2 (ΔPf = − 0.30) than with E1 (ΔPf = − 0.20): With the higher patch densities, large slip surfaces cannot any more squeeze between the patches (the same principle as observed for E1 with stand density). In general, Pf is lower for E2 than for E1 with TRS. In contrast, the average values of Pf are similar or higher with E2 than with E1 when assuming SRS. The reason for this phenomenon consists in the fact that the slip surface bottom strongly interacts with TRS, whilst only the edges (< 1.5 m) interact with SRS, so that the influence of SRS is much smaller.

Similarly to scenario I, SD of Pf is larger for E1 than for E2 due to a reduced amount of smoothing. However, due to the additional variation in the spatial distribution of the root systems, SD is generally higher than in scenario I.

E3 generally yields higher M of Pf than E1 and E2. This phenomenon is particularly pronounced for S2.7: Whilst E1 leads to a distinctive spatial pattern with several raster cells of low Pf, E2 results in a smoothed pattern of intermediate values of Pf without notable peaks. E3 represents a combination of both, imposing the peaks (small slip surfaces) upon the smoothed background (large slip surfaces). As in scenario I, the M of Pf are highest for E3 compared to E1 and E2. SD of Pf is intermediate between the values yielded for E1 and E2.

Scenario III: mixed patches

Surprisingly, both SRS/TRS and MRS stands with S = 2000 and P = 25–30 (S3.3 and S3.6) show higher M of Pf compared to those SRS/TRS and MRS stands (S3.2 and S3.7) where P (35–40 patches per plot) is higher but S (150) is much smaller (Fig. 10). This is also observed for the respective counterparts in scenario II. Moreover, SD is considerably higher for plots of scenario III with higher P. Totally mixed stands (S3.5–S3.8) show a generally higher SD compared to separated mixed patches (S3.5–S3.8). Except of S3.3, all plots show M of Pf for E2 than for E1. Interestingly, two main patterns are recognisable: (i) M of Pf in SRS/TRS plots is generally higher compared to MRS plots when P is low; (ii) Pf is generally lower for E1 at SRS/TRS plots, whereas Pf is lower for E2 at MRS plots. This means that dense and sparse MRS stands with low values of P (S3.6 and S3.8) show lower Pf compared to their respective counterparts (S3.3 and S3.1).

Fig. 10
figure 10

Pf for all plots of scenario III (description of colours, dots and circles cf. Fig. 8)

Full size image

Considering 50% around the median of Pf, MRS plots show a tendency of lower Pf than SRS/TRS plots (Fig. 11). This is particularly the case for low values of P (Fig. 11A, E, I and C, G, K). MRS stands seem to have a higher variation than SRS/TRS stands; however, generally, MRS stands show lower M of Pf.

Fig. 11
figure 11

Computed Pf values for scenario III for increasing numbers of trees per plot with varying stand and patch densities (indicated in blueish shades). The experiments (E1–E3) applied on each plot combination of scenario III are represented by the stacked plots (indicated in reddish shades) with details about the respective experiment setup. Circles depict separated stands (SRS/TRS), diamonds mixed stand conditions (MRS). The presented patches of each column are illustrated at the bottom of the figure

Full size image

Discussion

Impact of the slip surface size

Root systems are only able to contribute to reinforcement when the potential slip surface intersects but not entirely contains the paraboloid, since (Greenwood et al. 2004; Cammeraat et al. 2005; Van Beek et al. 2007). In the context of the present work, this concerns SRS, whereas TRS are always able to add root cohesion to the slope because rmax(TRS) > D. We construe that the reinforcement potential of TRS decreases with slip surface size for low stand densities due to the decreasing share of the potential slip surface coinciding with the root systems. The paraboloids seem to influence stability more effectively when they are rather located at the edges of the slip surface. This appears to be plausible since the intersected area of paraboloid and slip surface boundary of both SRS and TRS is larger when the paraboloid is located at the edge. However, since the angle that determines the curvature of the SRS paraboloid is more pointed than the one of TRS paraboloid, the added cohesion of SRS is larger at the slip surface edge compared to TRS.

This finding is in accordance with the issue highlighted in the research of Danjon et al. (2008), who stated that trees with many thick vertical roots fix a soil package better when located in the middle of a slope. In contrast, the contribution to slope reinforcement of root systems with oblique roots is rather higher, when located at the edge of a potential slip surface. Therefore, we assert that the r.slope.stability can reconstruct these relations, which appear consistent with the findings of Danjon et al. (2008).

Effect of species-related root system type and stand density

The results indicate that TRS contribute more to stability than SRS for all observed slope sections, expressed by a lower Pf throughout all scenarios (see Figs. 7, 8 and 9). TRS and MRS generally show lower values of Pf than SRS, indicating that paraboloids representing deep-rooting have a positive effect on slope stability. This is a consequence of the spatial interaction of the ellipsoid-shaped slip surfaces and the paraboloid root systems, where SRS do not penetrate the potential slip surface where D ≥ 1.5 m. The spatial extent of the root system itself—and thus its ability to penetrate a potential slip surface—strongly depends on the soil depth and inclination (Nilaweera and Nutalaya 1999). Kusakabe (1984) gave indications on the stabilising potential of different root systems in respect to specific slope conditions. In this regard, the relation between soil thickness (that strongly determines the depth of a potential slip surface) and the root system morphology should be under particular consideration in further studies.

The findings of scenario I have shown that stand density strongly influences Pf (see Figs. 7a and 8). This is in line with field observations and tree mapping of Roering et al. (2003) in a landslide prone area on the Oregon Coast range, which showed the effect of thinning (e.g. due to fires or timber harvesting) on slope failure occurrences. The varying extent of root systems of different tree species in landslide scarp adjacencies has thus a considerable effect on slope stability. The results of scenario I reproduce the patterns reported by Roering et al. (2003), since Pf is strongly associated with the planar extent of the root systems. The effects on slope stability caused by thinning and tree removal that lead to a decrease of stand density (and thus to a decreased spatial extent of root systems that can contribute to stability) were also highlighted by Genet et al. (2008).

Our results further indicate that TRS tends to react more sensitively to sparse stand densities than SRS (Figs. 10 and 11). This suggests that the r.slope.stability is sensitive to spatially delimited root cohesion values, depicted by differently distributed root system morphologies. Although the intention of the present study was to distinguish two different root system morphologies (SRS and TRS), these results might enable to implement more sophisticated information about the forest stand (e.g. species or tree type, growth stage or development stage of the root system). This can be important when the spatio-temporal assessment of species-related slope reinforcement by various root system morphologies is intended. Watson et al. (1999) compared two tree species with different spatio-temporal growth patterns of the root systems according to their respective slope reinforcement potential for different growth stages. They highlighted that distinct root systems would be the favoured to rapidly reinforce a slope since the initial lateral root growth, and thus, the development of the planar root system coverage is faster in an early stage of growth. The fact that Lcrit of sparse TRS stands (S1.5) is remarkably smaller (ΔLcrit = − 8.6 m) compared to dense TRS stands, indicates that the planar extent of the simulated root systems is highly associated with the stand density of the plot. In contrast, SRS stands seem to behave more robust (ΔLcrit = − 0.4 m) on stand density decrease than TRS (Fig. 6). This might be caused by the larger horizontal extent of SRS (rmax = 7 m) compared to TRS (rmax = 5 m) and thus leads to a higher potential contribution to stabilisation when trees are sparsely distributed. Consequently, a sparse density of trees might lead to the effect that a landslide (slip surface) is located right between the stems. The reinforcement potential of root systems with smaller planar extent (TRS) decreases more than the reinforcement potential of SRS when stand density is sparse (Fig. 12).

Fig. 12
figure 12

Effect of planar root system extent and stand density. a SRS with tree position Si. b TRS with tree position Ti. Slip surface size corresponds to the computational experiment E1 with the ellipsoid centre P. The overlapping areas between root systems and slip surfaces indicate the areas where the root systems contribute to stability

Full size image

Separated vs. mixed stands

The results of scenario III show that mixed stands have a beneficial effect on slope reinforcement when P is lower (Fig. 10). We assume that a consistent portioning of SRS and TRS considerably reduces Pf. In this regard, P, the distance between the patches and the size of the patches themselves might be the determining factors whether separated mixed stands (SRS/TRS) or totally mixed stands (MRS) can reduce slope failure probability. Thus, MRS appear to be preferable when P is low because then, the chance is higher that a root system occupies a position where it can deploy its maximal reinforcement potential (e.g. SRS close to the landslide margins and TRS in the middle of the landslide). Results presented in Fig. 11 indicate that completely mixed stands tend to decrease Pf when stand density is sparse, compared to plots where species are separated (SRS/TRS). Moreover, the findings suggest that in general, a higher amount of trees within a plot contribute more to slope stability and thus to a reduction of Pf. This is observable in scenario III (cf. Figs. 6 and 11); when the patch density is higher, overall, more trees occupy the whole plot, and voids between patches are smaller. This indicates that particularly a higher patch density leads to a reduction of Pf in environments with voids (forest glades), because more trees are able to reinforce the slope and number/size of voids are reduced (compare plots AEI and CGK with plots BFJ and DHL in Fig. 11). Herein, particularly plots with mixed stands (MRS) tend to be more stable. We suggest that the stand mixture (MRS) leads to this remarkable decrease of Pf. This is because SRS and TRS trees are then located at positions where they might be able to deploy a maximum of their position-dependent stabilising potential (e.g. SRS at ellipsoid edges, TRS in the middle of ellipsoids), compared to stand compositions where SRS and TRS are separated. Genet et al. (2010) reported that stand diversity does not primarily affect slope stability, however, rather the tree position on the slope and the architecture of the respective root system that crosses the potential slip surface (also suggested by results of Danjon et al. 2008). The model outputs indicate that the r.slope.stability is able to reproduce the positive effect of favourable positions of the different root systems, considering their distinct root system morphologies.

Implications for model parameterisation

Temgoua et al. (2016) recently used a partly comparable approach, where root system morphologies were approximated by geometric solids and implemented in a finite element environment. Referring to these findings, we hypothesise that the application of paraboloids depict a further step in representing root system morphologies in 2.5D and 3D slip surface models. We assume that paraboloids are flexible in reflecting the rooting peculiarities of different tree species, since root systems tend to form tropisms and diverse distribution of lateral and oblique roots when topographic conditions change (e.g. increased slope inclination). We expect that the estimation of root systems by paraboloids facilitates the adaption of geometric solids to several types of root architecture, e.g. SRS or TRS. Considering the results of Danjon et al. (2008) it can be anticipated that the implementation of these approximated morphologies into slip surface models might yield a better representation of various root systems and their spatial impact (Fig. 13).

Fig. 13
figure 13

Differences of root system positions along the slope. Reticles depict the position of the tree stems (Si or Ti). Sx depict the intersected area between root systems and slip surface, β is the slope angle. rmax and zmax depict the maximum radius and maximum depth that determines the spatial spread of the paraboloid

Full size image

It was reported by many former studies that soil characteristics highly influence the rooting behaviour of plants and thus the development of the root system (e.g. Greenway 1987; Bischetti et al. 2009; Ghestem et al. 2011; Stokes et al. 2014). However, we did not consider any spatial changes of the approximated root systems due to changing geotechnical soil parameters. We reason that a connection between the distribution of the root system architecture and apparent soil properties appears to be challenging for implementation in 2.5D slip surface models and applications in real-world case studies. Further investigations could tackle this issue and consider the variability of root systems in different soil environments and incorporate hydrological impacts of root systems on the soil, e.g. by root water uptake (Zhu and Zhang 2015). Based on the results of this study, we highlight that the spatial distribution of cohesion values within the approximated root system should be favoured compared to uniform cohesion values. For example, Bischetti et al. (2005), Bischetti et al. (2009), Thomas and Pollen-Bankhead (2010), Schwarz et al. (2010) and Schwarz et al. (2012) addressed the spatial distribution of roots, root tensile strength and thus dispersal of root-soil cohesion values. However, the spatial diversity of cohesion values within root systems remains to be a vague physical input which is hard to parameterise. In terms of 2.5D slip surface modelling, the consideration of the spatial distribution of root cohesion should be considered in further studies. Moreover, the complexity that arises, when estimating root systems properly (and considering soil characteristics or external environmental influences), identifies the general challenge of an accurate implementation of root properties in physically based models. Many deterministic approaches use the FoS to state whether a hillslope becomes unstable according to distinct physical parameters. Mergili et al. (2014a) and De Lima Neves Seefelder et al. (2016) emphasised that FoS—derived with a fixed set of geotechnical parameters—might fail to capture the details of a landscape. The wide range of root cohesion forces associated to different root systems would rather promote misinterpretations of the FoS due to the uncertainties associated with the dissimilarities of root distribution. Therefore, we suggest to use Pf rather than FoS as slope stability indicator particularly in those cases where root properties are used as input. Regarding a practical application of our approach, we highlight that studies applied on real-world conditions require ground-truth data about the geotechnical parameters, information about landslide dynamics (e.g. provided by a multi-temporal landslide inventory) and knowledge about forest stand conditions to avoid the induction of uncertainties.

The findings presented shall be employed as a basis to better parameterise root system morphology in real-world slope stability modelling. However, finding strategies to transform the results to real-world conditions remains a challenge:

Compared to traditional catchment-scale slope stability modelling, a finer spatial resolution if the GIS raster cells is needed to appropriately capture the root morphology. Whilst this would not be a problem with the infinite slope stability model, computationally, more complex approaches such as employed by the r.slope.stability could run into difficulties with the computational systems usually available, particularly for larger areas. Our study is performed on a 2 m × 2 m raster grid with a relatively manageable extent. However, it should be considered to elaborate in more detail on the ‘optimum’ cell size for real-world case studies that allows a reliable representation of the root system morphology whilst ensuring a justifiable computational time.

Alternatively, we propose the derivation of a set of rules to implement the findings gained on the generic topography in real-world case studies. However, we note that the results obtained in the present study build on particular parameter assumptions. Even though we postulate that the general patterns obtained have a broader range of validity, more computational experiments with a broad range of conditions of root system morphology, slope, water status etc. have to be performed. We further note that, in the present study, all non-tree roots are disregarded.

The findings reveal that the implementation of even simplified approximated root systems in a 2.5D slip surface model yield highly nonlinear effects on the model output. Further studies that apply root system approximations in 2.5D slip surface models should therefore particularly focus on an accurate representation of sensitive parameters such as stand density, stand arrangement (mixed or separated) and tree species.

Conclusions

The added cohesion from root systems and its implementation in 2.5D or 3D slip surface models has not been performed so far. To explore the potential of root system morphology implementation in a 2.5D slip surface model and its impact on the model performance, we used a set of idealised root systems as input for the computational tool r.slope.stability. The model was tested for 23 plots depicting different types and configurations of root systems, summarised in three scenarios (uniform tree distribution, separated patches and mixed conditions). Moreover, the scenarios were tested for three different sizes of potential shallow slip surfaces. We computed Pf to determine the added soil reinforcement by the assumed root systems. Our results show that the differently approximated root system morphologies exert distinct effects on the slope stability, where stand density, approximate position on the simulated slip surface and the ability of the root system to cross the potential slip surface are the driving factors for the model output. These findings are in accordance with the results of studies that revealed similar findings in field. Further work will focus on the implementation of more sophisticated ways of root system approximation, considering the spatial distribution of root cohesion values within the root system. In this regard, a better representation of the biomechanical interactions of the plant-soil-continuum is expected. We will further attempt to employ the insights gained for better parameterising real-world slope stability modelling campaigns.