Inclusion of Nature-Based Solution in the Evaluation of Slope Stability in Large Areas

Abstract

In areas affected by mining, which are undergoing reclamation, their geotechnical characteristics need to be monitored and the level of landslide risk should be assessed. This risk should preferably be reduced by nature-based solutions. This paper presents a KurZeS slope stability assessment technique based on areal data. This method is suitable for large areas. In addition, a procedure is presented for how to incorporate a prediction of the impact of nature-based solutions into this method, using the example of vegetation root reinforcement. The paper verifies the KurZeS method by comparing its results with the results of stability calculations by GEO5 software (version 5.2023.52.0) and validates the method by comparing its results with a map of closed areas in the area of the former open-cast mine Lohsa II in Lusatia, Germany. The original feature of the KurZeS method is the use of a pre-computed database. It allows the use of an original geometrical and geotechnical concept, where slope stability at each Test Point is evaluated not just along the fall line but also along different directions. This concept takes into account more slopes and assigns the Test Point the lowest safety factor in its vicinity. This could be important, especially in soil dumps with rugged terrain.

1. Introduction

Large-scale landslides can cause injuries, casualties, property damage and social as well as ecological problems [1]. Therefore, there is a strong need to prevent these events. This can be achieved and supported by an appropriate landslide susceptibility assessment, realized through a correctly selected evaluation method or software. The evaluation methods were intensively developed in the past decades. Nowadays, as mentioned in [2] and similarly in [3], they are usually divided into several classes: (1) qualitative methods that are based on expert experience, (2) semi-quantitative methods that are supplemented with scores/weights assigned to particularly assessed landslide susceptibility triggering aspects and (3) quantitative methods that can be further divided into geotechnical methods (physically based slope stability evaluation) [4,5,6,7,8], statistical analysis [2,9], neural network methods [10,11] or neuro-fuzzy logic methods [12,13].

In this paper, a physically based geotechnical model (a quantitative method) is considered. It is specifically a two- or three-layered geotechnical model with the use of the Bishop [14] (circular slip surface) and Sarma [15] methods (polygonal slip surface). In this paper, the method is referred to as KurZeS. Whereas many methods and software tools based on similar ideas use an assumption of infinite slope, e.g., [16,17,18,19,20,21,22], KurZeS does not, similarly to methods implemented in GRASS GIS [23,24,25], but it differs significantly in its geometrical concept described in Section 2 and also in its technical implementation using interpolation based on pre-computed data. The method provides an alternative procedure for slope stability assessment in large areas and an original implementation allowing calculations to be performed on desktop computers with limited computing power.

Shallow landslides, in particular, that can be assessed as described in [16,18,19,20,22,26,27,28,29] can be partially avoided through nature-based solutions (NBS). Root reinforcement [30,31] is the particular effect that stabilizes the slope in a natural way. Root reinforcement has not been investigated since at least 70 years ago. It is usually simulated either through additional cohesion of the rooted top soil layer [32] or through the Finite Element Method [33].

In KurZeS, the root reinforcement effect is simulated using the additional root cohesion concept. Evaluation of the additional root cohesion was described in [32,34,35,36,37]. KurZeS uses approximate values derived from values mentioned in [34,35,37].

For the visualization of the inputs and outputs of the model, we used the freely available geographic information system QGIS [38].

This article aims to analyze the influence of selected, nature-based solutions on the determination of the approximate slope stability on a large area. As typical examples, one can take areas intended for reclamation after open-pit coal mining. In the case of large areas, the prevailing reclamation method is afforestation or covering with other vegetation. Tree species recommended for reclamation of brown coal dumps are birch (Betula), larch (Larix), red oak (Quercus rubra) and pine (Pinus) [39]. In such a case, the main stabilization effect comes from the rooting of the surface layer.

The analysis is carried out in stages as follows:

• Evaluation of the top layer shear parameter increase due to the effect of rooting. This stage was based on the available literature data.
• The application of increased parameters in the large area surface stability model and realization of calculations using the original method created by the authors of the article.

The article is structured in five numbered sections. The first, the Introduction section, includes the background and context of the presented method. The second, the Materials and Methods section, is devoted to the description of the KurZeS method itself, its features and problems, the inclusion of nature-based solutions, the description of the testing problem and the description of the verification and validation approach. The third, the Results section, includes the results of the testing problem and results of the verification computations of specific details, and their comparison. The fourth, the Discussion section, includes specific observations concerning the KurZeS method and the inclusion of the NBS in it. The last, the Conclusions section, summarizes the main results and observations of the article.

2. Materials and Methods

2.1. Study Area

For this study, the area of the soil dump of the former open-cast mine Lohsa II in Lusatia, Germany, see Figure 1, was chosen. The map presented below was obtained from [40]. The area is interesting because several landslides have appeared there in recent times. Specifically, the slopes of the mining waste dump Scheibe at the northwest border of the closed coal mine Lohsa II moved in the spring of 2019. Even later (in the spring of 2021), landslides appeared on the slopes near the lake Knappensee.

2.2. Geometrical Concept

The used data describing the digital elevation model (DEM) of the terrain within the study area are organized in triplets of coordinates [X; Y; Z] with a regular step of 10 m on both the X and Y axes. The safety factor is evaluated in each point of the 10 × 10 m raster (this point is referred to as the TP—the Test Point). The safety factor corresponding to the TP is derived as the minimal safety factor of all relevant slopes in the TP’s neighborhood. The slopes relevant to the TP are defined by pairs of boundary points (${\mathrm{BP}}_1$ and ${\mathrm{BP}}_2$) satisfying the following conditions:

• The boundary points are different points located at the points of the 10 × 10 m raster not further than 100 m from the TP.
• The TP may coincide with one of the boundary points.
• The total length of the (possibly angled) line ${\mathrm{BP}}_1-\mathrm{TP}-{\mathrm{BP}}_2$ must be at least 20 m and at most 100 m.
• The convex obtuse angle of line segments ${\mathrm{BP}}_1-\mathrm{TP}$ and $\mathrm{theTP}-{\mathrm{BP}}_2$ shall be greater than or equal to 120^°.

The safety factor of the TP is then evaluated as the minimum of the safety factors of all slopes between ${\mathrm{BP}}_1$ and ${\mathrm{BP}}_2$ satisfying the above-specified conditions. The algorithm KurZeS includes an optimization trick based on the obvious equivalence of the slope ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ and ${\mathrm{BP}}_2-{\mathrm{BP}}_1$ reducing the number of safety factor estimations to a half (see Figure 2).

2.3. Geotechnical Concept

The safety factor of each slope ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ is evaluated with interpolation from a pre-calculated database based on the following geotechnical model constructed in GEO5 [41] (version 5.2023.52.0).

2.3.1. Model

The basic model is constructed as two-layered, with a constant gradient of the terrain in variants in the range from 1 in 2 to 1 in 10. The interface between the soil layers runs parallel to the terrain at a depth of 5 m. The parameters of the top layer are exactly those of the soil it contains. The underlying layer is supposed to consist of the same type of soil but firmer, which is modeled by higher stability parameters compared to the top layer. This increase reflects the usual improvement with depth (see Figure 3a and Table 1).

A more complex model is used for the case of root reinforcement. Here, the third layer of a constant thickness of 0.9 m is added to the surface (see Figure 3b and Table 1). The used soil parameters in it are also higher than in the original soil, this time due to root reinforcement. The specific values of increase are derived from the literature [34,35,37].

The groundwater level is considered parallel to the surface in variants from a depth of −20 to 0 m.

According to the method used, circular (Bishop method [14]) and polygonal (Sarma method [15]) slip surfaces with a fixed length from 20 to 100 m are considered.

2.3.2. Resulting Database

The result of an individual stability calculation in the given model is the search for the least favorable slip surface (either circular or polygonal, using the optimization procedure) and the determination of the safety factor ${\mathrm{F}}_{\mathrm{s}}$ by the limit equilibrium method (LEM). ${\mathrm{F}}_{\mathrm{s}}$ values were calculated for all options including 8 slope inclinations, 6 slope lengths, 7 groundwater levels and 2 cases with and without a rooted topsoil layer. For each selected soil type, the results database consists of 672 ${\mathrm{F}}_{\mathrm{s}}$ values.

The safety factor of a specific slope ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ is derived from this database so that the value of the corresponding slope inclination, slope length and groundwater level of the slope ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ is specified and the ${\mathrm{F}}_{\mathrm{s}}$ value is interpolated from neighboring values tabulated in the database.

It is typical for the particular two- or three-layer model with parallel construction described above that the type and shape of the slip surface changes as the least favorable slip surface lengthens. Short slip surfaces with the lowest ${\mathrm{F}}_{\mathrm{s}}$ are typically circular in shape (B—Bishop) confined in the uppermost layer, or in the central and rooted layer, not extending into the lowest firmer layer with an interface 5 m below the surface (Figure 4a). The intermediate slip surfaces are most commonly polygonal (S—Sarma) following the top of the interface at 5 m depth. In terms of type, they most closely resemble the infinite slope model (Figure 4b). In some cases, usually with the groundwater level between depths of 3 and 6 m below ground level, long slip surfaces of circular shape cross the aforementioned interface at a depth of 5 m below the surface and extend deep into the firmer layer (BD—Bishop, deep)—Figure 4c.

It is not clear in advance which of the slip surfaces will be the least favorable. The result is mainly influenced by the geometry of the model, soil parameters and the position of the water table or surface water level. It should be emphasized repeatedly that the results are based on the chosen model; for a differently designed model, the composition of the least favorable slip surfaces could be different.

2.4. Quantification of the Effect of Rooting

Figure 5 shows an example from the resulting database of ${\mathrm{F}}_{\mathrm{s}}$ for one selected soil type denoted as S5Y (sandy dump). The left half of the table applies to the model without a rooted topsoil layer, and the other to the model with a rooted layer. In the header of the table, the parameters used for effective strength (effective angle of internal friction (^°) and cohesion (kPa)) and bulk weight (kN/m³) are typed. The first column shows the depth of the water table, the second column shows the slope inclination and the first row under the heading shows the fixed length of the slip surface. The types of slip surfaces in the sense of the text above are distinguished by color for cells with the resulting safety factors ${\mathrm{F}}_{\mathrm{s}}$ (short circular slip surface in blue, medium-length polygonal slip surface in green and deep circular slip surface in red).

2.5. Data Sets

All data sets were obtained from freely available sources and they were transformed into the form of georeferenced rasters, where the centers of gravity of particular raster cells are the nodes of the network described in Section 2.2 dedicated to the geometrical concept of the model.

2.5.1. Digital Elevation Model

The main data source for performed simulations was the digital elevation model (DEM) DGM1 (German DEM [42]). This DEM (multiple shape files) was rasterized with the use of QGIS [38]. The selected cell edge length was 10 m. This selection of the cell edge length brought a compromise between acceptable computation times (huge number of elements/nodes) and sufficient accuracy of the slope description, see Figure 6.

2.5.2. Ground Water Table

Because of missing data, the groundwater table was derived from the DEM. For this purpose, raster calculator “gdal_calc”, which is a part of GDAL-extension [43] for QGIS, was used together with Equation (1), where $\mathrm{W}$ means water elevation, $\mathrm{A}$ represents a surface elevation in DEM-raster, 115 stands for the selected lake water level, $0.75$ defines the “speed” of water table increase and 20 and 0 are the maximal and minimal groundwater table depths, respectively.

$$W=\left\{\begin{array}{cc}115\hfill & \mathrm{if}\mathrm{A}\le 115\hfill \\ \mathrm{A}-20\hfill & \mathrm{if}(\mathrm{A}-115)·(1-0.75)>20\hfill \\ 115+(\mathrm{A}-115)·0.75\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(1)

The resulting water table raster is shown in Figure 7.

2.5.3. Type of Soils

A pedological map was obtained from LfULG (Bodenkarte BK50 [44]). Shapefiles from the original map were merged and divided into two or four (non-rooted soil and rooted soil, respectively) different groups, and they were rasterized in the same spatial extent and with the same cell size (10 × 10 m) as the DEM, see Figure 8a. For the case of rooted soils, two further soil types appeared, see Figure 8b. The distribution of rooted and non-rooted soil was based on the difference between the DEM and the selected lake water level (115 m MSL). Rooted soils were considered for positive values (above water) and non-rooted soils for negative (below water).

The strength parameters of soils can be found in Table 1.

Table 1. The strength parameters of used soil types. ${\phi }_{\mathrm{eff}}$ denotes effective angle of internal friction and ${\mathit{c}}_{\mathrm{eff}}$ stands for effective cohesion of soil. The thickness of the topsoil layer is 5 m and there is no rooted soil layer in case of no measures. In the case of 0.9 m rooting depth, the thickness of the rooted soil layer is 0.9 m and the thickness of topsoil layer is 4.1 m. Parameter values are based on expert judgment combined with information from the literature [34,35,37].

	Basal Soil Layer		Topsoil Layer		Rooted Soil Layer
	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$
S45 *	31°	16 kPa	27°	8 kPa	27°	24 kPa
S5Y *	29°	4 kPa	25°	2 kPa	25°	18 kPa

* The soil type notation comes from the Czech Technical Standard ČSN P 73 1005 (731005) Engineering geological survey [45]; the Y symbol denotes the relocated material in non-original placement.

Table 1. The strength parameters of used soil types. ${\phi }_{\mathrm{eff}}$ denotes effective angle of internal friction and ${\mathit{c}}_{\mathrm{eff}}$ stands for effective cohesion of soil. The thickness of the topsoil layer is 5 m and there is no rooted soil layer in case of no measures. In the case of 0.9 m rooting depth, the thickness of the rooted soil layer is 0.9 m and the thickness of topsoil layer is 4.1 m. Parameter values are based on expert judgment combined with information from the literature [34,35,37].

	Basal Soil Layer		Topsoil Layer		Rooted Soil Layer
	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$	${\mathit{\phi }}_{\mathbf{eff}}$	${\mathit{c}}_{\mathbf{eff}}$
S45 *	31°	16 kPa	27°	8 kPa	27°	24 kPa
S5Y *	29°	4 kPa	25°	2 kPa	25°	18 kPa

* The soil type notation comes from the Czech Technical Standard ČSN P 73 1005 (731005) Engineering geological survey [45]; the Y symbol denotes the relocated material in non-original placement.

2.6. Implementation

Slope stability was derived from interpolation tables (resulting database, see Section 2.3.2) pre-computed in the GEO5 software [41] (version 5.2023.52.0). derived the safety factor from the slope length (in range 20 to 100 m), slope inclination (in range 1:10 to 1:2), soil type (2 types of rooted soil and 2 types without reinforcement) and depth of the water table (in range 0 to 20 m).

Two variants of interpolation were tested: linear interpolation and polynomial interpolation. As they had comparable results and linear interpolation was noticeably faster, the presented computations were performed using multidimensional linear interpolation with the use of the linearNDInterpolator implemented in SciPy [46].

The calculations work with a circular or polygonal slip surface of general depth. Such a slip surface crosses the root-reinforced layer only along the possible landslide scarp. Thus, its area is the only place where the improvement in soil properties is accounted for by this method. In this way, the application of root reinforcement differs significantly from commonly used infinite slope stability models, e.g., [47].

For the presented computations, a computer with the following parameters was used:

• Intel^® Core^TM i7-3770K CPU; 3.50 GHz × 4;
• RAM 15.5 GiB;
• Ubuntu 22.04.4 LTS;
• Python 3.8.16.

Computation for such a large-scale area was not possible at once because of extensive memory demands. The results were obtained by dividing the area into 36 overlapping subareas (see Figure 9), which effectively decreased the memory requirements. Eventually, the computation could be realized in four parallel threads, i.e., in four subareas at the same time. A 100 m overlapping of subareas was necessary because of the geometrical concept of the model, see Section 2.2. One turn (4 subareas) took 14 min.

2.7. Verification

The mechanism for verification of the results can be found in Figure 10. Two safety factor computations were performed—one without root reinforcement and another with root reinforcement under consideration. Two hypotheses were tested as follows:

• Because all ${\mathrm{F}}_{\mathrm{s}}$ in the pre-calculated database are higher for root-reinforced soils than for the soils without vegetation, the safety factor should be higher or equal in each Test Point of prediction considering root reinforcement than in prediction without root reinforcement.
• The results of the prediction should be comparable with ${\mathrm{F}}_{\mathrm{s}}$ predicted by GEO5 at each slope in the area.

Several slopes in the area were taken, the surface and groundwater table cuts were taken from QGIS and included in GEO5 software and the following observations were made:

• Consideration of root reinforcement brought either the same or higher value of the predicted safety factor in every Test Point.
• For selected cuts, the predicted safety factor corresponded closely to the values obtained from GEO5.

For more details, see Section 3.

The two above-mentioned points verified the expected behaviour of the simulation code.

2.8. Validation

The approach for validation of the results is shown in Figure 11. The computed safety factor for all the cells in the study area with and without additional root cohesion under consideration was imported to QGIS [38] and compared with the shape file map of the area closed to the public because of safety reasons. The closed area correlates well with the slopes with a predicted ${\mathrm{F}}_{\mathrm{s}}$ very close to one, which validates the results with only independent accessible public data.

Another method of validation was the comparison of landslide risk maps with orthophoto maps of the study area and also with maps showing the slope gradient there. A noticeable correspondence was observed between these data, but due to the unclear, precise dating and spatial location of the landslides, these comparisons were not included in the study.

3. Results

3.1. Slope Stability Improvement Predicted by GEO5

To get a better idea about the expectable predicted safety factor improvement caused by additional root cohesion, selected data sets from the resulting database were visualized and analyzed. Safety factors of two selected soil types present in the study area, sandy clay (S45) and sandy dump (S5Y), each for one slope inclination, are depicted with and without vegetation cover (root reinforcement) in Figure 12a,b.

Both Figure 12a,b show the dependency of the absolute and relative improvement of the predicted safety factor on the depth of the water table as well as on the root reinforcement under consideration. Safety factor values (left-hand axis), as well as relative improvement (right-hand axis), decrease with growing slope length. The effect of root reinforcement is stronger in the case of sandy dump. However small the absolute improvement may seem, it can be crucial for the cases where driving and resisting forces are near to equilibrium (i.e., for ${\mathrm{F}}_{\mathrm{s}}$ around 1). The most significant relative improvement is predicted for slopes where groundwater reaches near the earth’s surface, i.e., the least stable slopes.

3.2. Slope Stability Evaluation with KurZeS

As a first step for the testing of KurZeS, the predicted safety factor ${\mathrm{F}}_{\mathrm{s}}$ without the root reinforcement (additional cohesion) under consideration was evaluated, see Figure 13. The most landslide-threatened slopes appeared on the lake shore (e.g., “AREA 79” on the Figure 1) and on the steep slopes in the northwestern part of the waste dump Scheibe (“AREA 46” on the Figure 1).

Figure 14 shows the comparison of the predicted safety factor of slopes for non-rooted soils with a shape file map of the area closed to the public. It shows good correspondence between both data sets and it shows the beneficial effect of the KurZeS approach for the identification of unstable areas suggested for further in situ investigation.

Another prediction was made with the addition of the root reinforcement under consideration. The resulting prediction of the safety factor is denoted as ${\mathrm{F}}_{\mathrm{sR}}$. The result can be seen in Figure 15. The values of ${\mathrm{F}}_{\mathrm{sR}}$ are generally higher (better) than the values of ${\mathrm{F}}_{\mathrm{s}}$ (without root reinforcement). It will be further analyzed in the following Section.

3.3. Improvement in Slope Stability Predicted by KurZeS Due to Root Reinforcement

From raster outputs of slope stability analyses, absolute and relative improvement in the predicted safety factor obtained with or without the roots under consideration was evaluated. The absolute improvement ${\mathrm{F}}_{\mathrm{sR}}-{\mathrm{F}}_{\mathrm{s}}$ reaches values in the range $\langle 0.0,0.25\rangle$. Relative improvement of the predicted safety factor ${\mathrm{F}}_{\mathrm{sR}}/{\mathrm{F}}_{\mathrm{s}}-1$ depicted in Figure 16 reaches values in the range $\langle 0.0,17.7\%\rangle$. Both the maximal values overreach the expectations constructed on the base from Section 3.1. The explanation coming from further analysis resides in the fact that the best improvement appears for extremely short slopes, which were not considered in Section 3.1.

3.4. Comparison of Results from KurZeS and GEO5

To verify the results of slope stability evaluation by the KurZeS approach, several further studies of slope stability, in particular mesh cells, were performed. The KurZeS results were compared with those obtained from GEO5 for the particular slopes. In this section, we present two representative examples as well as data sets prepared for the following research.

For nine Test Points, the boundary points ${\mathrm{BP}}_1$ and ${\mathrm{BP}}_2$ realizing the minimal associated ${\mathrm{F}}_{\mathrm{s}}$ were identified. The connection ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ was prolonged in QGIS [38] and on such a cross-section, the model input data (DEM, depth of water table, soil type) were extracted from intersected raster cells. These cross-section data served as inputs to GEO5 software for cross-sectional safety factor evaluations. The Test Points were selected from three areas that were somehow interesting. The areas were named “AREA 13” (it included Test Points 1 to 3 and corresponding cross-sections 1 to 3), “AREA 46” (it included Test Points 4 to 6 and corresponding cross-sections 4 to 6) and “AREA 79” (which included Test Points 7 to 9 and corresponding cross-sections 7 to 9). Here, we present the results of cross-sections 5 and 8 that represent all nine cross-sections well.

3.4.1. Cross-Section 5

In Figure 17, you can see the map from [40] of “AREA 46” that includes very steep slopes that are not at the lake shore with well-defined slope direction (almost planar slopes of a long straight valley). Test Point ${\mathrm{TP}}_5$ that is not a part of the cross-section was selected. In Figure 17, in red, 13 sampling points ($\mathrm{SP}\_\mathrm{i},\mathrm{j}$) formed the cross-section model in GEO5.

For the non-rooted soils, KurZeS evaluated the safety factor ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_5\right)=1.10$. Results for the whole “AREA 46” are shown in Figure 18.

For the scenario with non-rooted soil, GEO5 evaluated the safety factor ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_5\right)=$$1.08$. Figure 19 shows a 2D cross-section with visualization of the data extracted for complementary safety factor assessment.

For the model with rooted soil, the situation is very similar. KurZeS predicted a safety factor ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_5\right)=1.18$. Figure 20 depicts the results evaluated for the whole “AREA 46”. Safety factor ${\mathrm{F}}_{\mathrm{sR}}\left({\mathrm{TP}}_5\right)$ in Figure 20 evinces evident improvement in comparison with Figure 18.

GEO5 evaluated the safety factor ${\mathrm{F}}_{\mathrm{sR}}\left({\mathrm{TP}}_5\right)=1.16$ for the rooted soils. Visualization of the input data extracted along $\mathrm{cross}\_\mathrm{section}\_5$ is shown in Figure 21.

The difference between the resulting safety factors obtained from KurZeS and GEO5 can be denoted as negligible.

This comparison of two very close results from two different software tools verified the function of KurZeS for the case of a model without surface water.

3.4.2. Cross-Section 8

Another subarea selected for verification of the results takes place at the northeastern slope of the mining waste dump Scheibe and represents very steep slopes finishing in surface water, see Figure 22. Here, we present the comparison of safety factor predictions at Test Point ${\mathrm{TP}}_8$, which lies directly on cross-section 8.

The KurZeS approach predicted the safety factor ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_8\right)=1.24$ for slopes covered by non-rooted soil. You can see the results for the whole “AREA 79” in Figure 23. The direction of cross-section 8 is very close to fall lines (almost perpendicular to contour lines).

The safety factor of the slope covered by non-rooted soil and evaluated by GEO5 has a value ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_8\right)=1.25$. In Figure 24, you can see a 2D $\mathrm{cross}\_\mathrm{section}\_8$, which was extracted from input files for KurZeS and evaluated in GEO5.

The difference between safety factors predicted by KurZeS and GEO5 for slopes formed by non-rooted soils is even smaller than in the case evaluated in “AREA 46” and it can be denoted as negligible.

For the case of a slope formed by rooted soil, KurZeS predicted safety factor ${\mathrm{F}}_{\mathrm{sR}}\left({\mathrm{TP}}_8\right)=$$1.30$. The resulting raster for the whole “AREA 79” can be seen in Figure 25. Lighter colors than in Figure 23 indicate an improvement in slope stability due to root reinforcement.

Verification result obtained for the case of rooted soil from GEO5 has the value ${\mathrm{F}}_{\mathrm{s}}\left({\mathrm{TP}}_8\right)=1.27$ and the appropriate $\mathrm{cross}\_\mathrm{section}\_8$ with rooted soil is depicted in Figure 26.

For the case of a slope formed by rooted soil, the difference between safety factors obtained from KurZeS and GEO5 is slightly higher than in previously evaluated cases and it makes for approximately 2.3%, which we feel to be still very low.

This comparison of two very close results from two different software tools verified the function of KurZeS for the case of a model including surface water.

The above-mentioned results show that the results of KurZeS using the pre-calculated database are very comparable with the results of GEO5 for the same slopes. The analysis was conducted both for slopes finishing in surface water and slopes finishing in land, and both for cross-sections including their Test Point and cross-sections that passed it. The comparability of the results verifies the approach using the pre-calculated database as a good simplification of the problem. This approach efficiently allows many stability factor evaluations for each Test Point and finds the weakest slope in its vicinity.

4. Discussion

Commercial software tools for slope stability evaluation have offered the “Method of Slices” to compute safety factors in 2D slope cross-sections. The KurZeS approach uses the “Method of Slices” implemented in GEO5 for pre-computed interpolation tables for selected soil types. These tables are used to approximate the slope stability of large areas.

The use of a pre-computed database allows the KurZeS approach presented in this article to bring an original geometrical and geotechnical concept, where slope stability in each Test Point is not evaluated just along the fall line but also along rather different directions. This concept takes into account more slopes and assigns the TP the lowest safety factor in its vicinity. This “looking around” the TP approach could be important, especially in places where surface water or material interface appear. In such a place, the weakest slope may not be the one that runs exactly along the fall line passing the Test Point.

As in all applications, the accuracy must be balanced with simplicity. The price that should be paid for achievable large-scale models is a simplified two- (case without NBS) or three- (case with NBS) layer stratigraphy, which has been evaluated as sufficient detail for the soil dump slopes under consideration.

In most cases, the simplifications and inaccuracies of the method are related to the mesh discretization density and schematization of the model. In the geometric part, slip surfaces shorter than 20 m are not allowed.

The method allows the ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ line to pass the TP. The effect of this feature on the accuracy of the result is largely limited by the fact that the least favorable slip surfaces described by both pairs of line segments ${\mathrm{BP}}_1-\mathrm{TP}$ and $\mathrm{TP}-{\mathrm{BP}}_2$, or the line segment ${\mathrm{BP}}_1-{\mathrm{BP}}_2$, usually have the minimum angular deviation from the slope gradient and are close to each other in shape according to the experience analysis. This means that the angle between ${\mathrm{BP}}_1-\mathrm{TP}$ and $\mathrm{TP}-{\mathrm{BP}}_2$ is close to ${180}^{°}$, and therefore, the distance between $\mathrm{TP}$ and ${\mathrm{BP}}_1-{\mathrm{BP}}_2$ is small. Examples of cases with different probabilities of finding the least favorable slip surface with a comparison of the degree of inaccuracy are shown in Figure 27.

The inaccuracies of the method resulting from the influence of the mutual rotation of the 10 × 10 m mesh and the gradient direction increase with the decreasing length of the slip surfaces. Therefore they are analyzed in greater detail for the cases of the shortest slip surfaces with a changing angle between the contour lines and the orientation of the 10 × 10 m mesh. For this analysis, a completely homogeneous model is considered, without the influence of groundwater. Therefore, only the slope of the terrain on the line and its length have influence.

For $\mathrm{TP}\ne {\mathrm{BP}}_1$ and the ${120}^{°}$ angle, there are only two patterns with three lines (Figure 28, the first and second columns of three) with a matching ${\mathrm{BP}}_2$, that are repeated four times. In the case of $\mathrm{TP}\equiv {\mathrm{BP}}_1$, there is one pattern repeated four times (third column).

In the first triad, the contour lines have the same direction as the horizontal X axis and the slope as the vertical Y axis. In the second triplet, these lines form an angle of $26,6^{°}$ and in the third, ${45}^{°}$. One of the possible links for $\mathrm{TP}\ne {\mathrm{BP}}_1$ and $\mathrm{TP}\equiv {\mathrm{BP}}_1$ is always perpendicular to the contour lines, and its steepness is maximal under the given conditions, but its length changes by leaps from 20 m to 22.36 m ($10\times \sqrt{5}$) and 28.28 m ($10\times 2\times \sqrt{2}$)—Figure 29.

In favor of the adherence of the method to reality, the regularities in terrain modelling, especially in the case of reclamation, the “inverse proportion” of the slope gradient and length are at work. The longest slope cannot be the steepest one.

Further, the potential slope stability improvement achieved through root reinforcement represented by additional root cohesion was studied. Tree species recommended for reclamation of brown coal dumps are birch (Betula), larch (Larix), red oak (Quercus rubra) and pine (Pinus) [39]. In the presented simulations, a special type of slope with specific soil composition was considered. This kind of slope is common in soil dumps that are present in brown coal mining areas, e.g., in the Northern Bohemian Brown Coal Basin or Saxony and Lusatia, Germany.

The values of the stability parameter included in the presented model are rather conservative values coming from the literature and an expert’s experience. Such a choice of parameters can avoid damage caused by an overestimation of the predicted slope stability.

Predicted slope stability improvements may seem to reach small relative values; however, in very unstable areas where the safety factor approaches the value of 1, the improvement may mean an important difference between the risk and safety.

In this study, neither rainfall [16,18,48,49] nor groundwater flow [50] were considered. These effects could be taken into account by coupling the KurZeS method with different appropriate model tools.

A probable plant-induced water table decrease was not considered. The expected stabilization effect is in the order of a few percent. However, this small improvement may play a crucial role in the boundary cases where the slope is near to collapse.

Temporal dependency of the stabilizing effect, which may reach its optimum after a couple of decades, was also not considered. It could be implemented using time-dependent stability parameters of the soils.

For the sake of simplicity, we have also not considered in the model, the effect of earthquakes, as in [51,52].

In this study, the volume of soil removed by landslides was not predicted. The description of this phenomenon can be found in [8,47,53] and possible methods of its evaluation are described in [54,55].

Before the stabilization effects of the industrial forest fully appear, another NBS (fascines, grass, shrub) may be applied to avoid the risk of shallow landslides and erosion. On a small scale, at specific places, NBS may be complemented with standard technical solutions used for slope stabilization. The same holds for the cases where NBS are not able to stabilize the slope in sufficient measures and where the residual risk of landslides remains too high.

5. Conclusions

The presented method allows the geotechnical stability of large areas to be analyzed with generally accessible data on a usual computer in a reasonably short time, which was demonstrated on a specific area of soil dump of the former open-cast mine Lohsa II. It also allows evaluation of the influence of possible natural-based or technological measures, which was demonstrated in the case of root reinforcement.

The two most original aspects of the presented approach are the use of a pre-computed database of safety factors of slopes of various parameters and the testing of the stability of the whole neighborhood of each Test Point (the associated slopes do not need to cross their TPs). Both these aspects are somehow connected since the pre-computed database shortens computational time so taking more possible slopes into account is allowed.

A compromise between the detail and the large area of the study area is the choice of a square network with a node spacing (TP) of 10 m and a limitation of the shortest slip surface length to 20 m. However, the method includes, to a controlled extent, a rather detailed examination of the immediate vicinity of the TP by also considering links that do not directly pass through the TP. By taking into account slopes, directions, slip surface lengths and other parameters, the method is suitable for rugged terrain with variable soil composition or in contact with water bodies, despite the limitations.

The possibility of the inclusion of NBS was presented in the study of the influence of possible safety factor change due to vegetation (root reinforcement) that can be performed on a computer to help identify the slopes that can be stabilized by root reinforcement and do not need other (specifically technological) measures.

The presented approach was verified and validated in the area of soil dump of the former open-cast mine Lohsa II in Lusatia, Germany. Verification was performed by comparing several slope stability computations obtained from the KurZeS method and GEO5 software. The validation was performed by comparing the KurZeS results with the published map of areas closed to the public.

The maximum predicted relative increase in the stability factor due to the consideration of nature-based solutions was 17% and was observed for extremely short slopes with a length of 20 m or less. Shorter slopes are not accounted for by the model due to its geometrical concept.

However small this improvement may appear, it can be crucial in cases where the stability factor is close to 1.0, i.e., when the driving and motion-resisting forces are close to equilibrium.

The predicted maximum improvement exceeded expectations based on reference calculations from GEO5. This was due to the choice of the minimum tested slope length in the reference calculations being too high.

In further research, the authors will focus on the generalization of quantification of slope stabilization effects achieved through various NBS for various soil compositions. In the application, the focus will be directed towards the simulation of the stability of selected, large-scale mining waste dump slope areas. This will enable us not just to assess risks and threats but also to evaluate the potential to avoid shallow landslides using NBS.