Slope stability of deep surface coal mines in the presence of a weak zone

Coal and lignite (brown coal) are geo-resources that govern energy production for decades. Geomechanical challenges, particularly slope stability, related to surface coal and lignite mines are critical during operation and determine the post-coal era. Several failure incidents have been reported in mining areas, typically associated with a sub-horizontal failure surface on a weak—clay or marl—layer or an interface of low strength. This weak zone controls the soil profile in terms of stability and is common in several mines globally. In this work, the finite element method with the shear strength reduction technique is primarily employed to evaluate the slope stability of this profile. Three geotechnical software are initially compared, and results are identified as practically identical. Moreover, slopes with benches, as typically in mines, and without benches, as typically in slope stability analysis, are compared, with the no-benches analysis being consistently more conservative. The crucial parameters' effect is then examined: the height and the inclination of the slope, the inclination, thickness, and strength of the weak zone, and the strength of the overburden soil. Their effect on slope stability is quantified by combining the probabilistic point estimate method with the finite element method. It is concluded that the inclination and strength of the weak zone and the water conditions are the most critical parameters and control the stability. This work can support a preliminary slope stability analysis and expands the knowledge and understanding of slope stability of a weak zone profile. Slope stability is performed on a soil profile with a weak zone, describing a common type of coal mines' failure in several countries. Numerical issues of practical importance are addressed for the simulation and analysis of the problem. The inclination and shear strength of the weak zone together with the water conditions control the stability of the mines' slopes. Slope stability is performed on a soil profile with a weak zone, describing a common type of coal mines' failure in several countries. Numerical issues of practical importance are addressed for the simulation and analysis of the problem. The inclination and shear strength of the weak zone together with the water conditions control the stability of the mines' slopes.


Introduction
Coal and lignite mining are geo-resources that have significantly contributed to the global energy needs during several decades, accounting for 25% of the global energy production in 2000, 30% in 2010and 27% in 2019(BP 2020. When the coal deposits are located relatively near the surface, their cost-efficient extraction is usually achieved by surface (also called open-pit or opencast) mining methods. The sustainability of such mines is significantly affected by the stability of the corresponding excavations. Several slope failure events in open-pit lignite and coal mines have been reported in the past (e.g. Prountzopoulos Solak et al. 2017;Zevgolis et al. 2019;Zhigang et al. 2020); such incidents set at risk human lives, infrastructure, and mining operation. Additionally, landslides pose a severe risk to the reclamation of mining areas; thus, even though several countries worldwide have pledged to cease any lignite or coalrelated activities soon (e.g. Canada, Germany, France, the United Kingdom and Greece), mining areas' slope stability will still be of importance.
Stability failure mechanisms are often associated with a horizontal or sub-horizontal failure surface developing on a low strength zone, usually located near the bottom of the excavation. This zone can be a layer or an interface between layers of low shear strength, hereafter named as the weak zone. As a result, a composite failure mechanism (rather than a circular surface) represents the failure surface (e.g. Kavvadas et al. 2020;Leonardos 2004), with sliding occurring along the weak zone and reaching the top of the slope with a planar or curvilinear transition. This mechanism is dominated by the weak zone's residual strength due to the immense shear strains developed on the relatively thin interface during ground movements. This type of failure is encountered in open-pit coal and lignite mining in several countries, such as Greece (Kavvadas et al. 2020;Leonardos 2004;Zevgolis et al. 2019), Turkey (e.g. Ulusay et al. 2014Ulusay et al. , 2001Ural and Yuksel 2004), Poland (e.g. Bednarczyk 2017), Czech Republic (e.g. Mencl 1977) and Australia (e.g. Ghadrdan et al. 2020;McQuillan et al. 2018). The weak zone mechanism is widespread in several landslide phenomena and several lignite and coal mines' slope failure events globally.
Failure incidents related to a weak zone (also named a weak layer) have been discussed and studied by several researchers outside the context of mining (e.g. Bromhead 1992;Urciuoli et al. 2007). Landslides governed by such failure mechanisms have been globally documented, e.g. in Canada, Hong Kong, Scandinavia, Spain and the United Kingdom (L' Heureux et al. 2012;Lastras et al. 2004;Puzrin et al. 2017;Quinn 2009;Zhang and Wang 2020). In these cases, various reasons may have triggered the reported incidents, such as the removal of the downslope support due to soil degradation, or a local increase of the pore water pressures due to seismic loading, or gas hydrate disassociation. Nevertheless, the failure surface's non-circular shape and the weak zone's importance are two features that remain alike for all these events.
Different approaches have been investigated towards the design and prevention of such stability issues. Huang (1983) has discussed this failure mechanism and proposed an analytical approach based on limit equilibrium analysis to assess the safety factor (SF). Kavvadas et al. (2020) also developed a simplified analytical model for Greek lignite mines and investigated the importance of critical parameters on the SF. However, the limitations of the analytical models and the increase of the computational power have made numerical analysis the primary tool (e.g. Cała et al. 2020;Ozbay and Cabalar 2015). Nevertheless, in most cases, numerical modelling is performed for back analysis of slope failures, assessing the weak zone's strength parameters (e.g. Tutluoglu et al. 2011). Overall, the weak zone's importance is discussed in many studies, but its effect on stability has not been systematically examined.
The present study fills the gap of systematic numerical analysis on the slope stability of deep surface excavations in the presence of a weak zone, typically encountered in open-pit lignite and coal mines. It contributes to the understanding and quantifying the effect of various parameters on the failure mechanism controlled by a weak zone.
The finite element method is primarily employed, combined with the shear strength reduction technique to calculate the SF. Additionally, the point estimate method was used to estimate the statistical moments of the SF and the probability of slope failure. Several researchers have combined probabilistic methods with finite elements to evaluate slopes' reliability (e.g. Ghadrdan et al. 2021;Griffiths and Fenton 2004). Through this framework, the critical geotechnical parameters' uncertainty can be quantitatively considered. Moreover, the probability of failure provides a more rational assessment of the risks associated with slope stability than the SF. It also quantifies better the slope stability than the SF: differences in the probability of slope failure quantify actual different risk levels. On the contrary, the safety factor works in principle qualitatively in terms of risk evaluation: low SFs (''relatively'' close to 1) denote possible failure, while high SFs denote that the slope is stable.
Initially, the critical parameters that govern the failure mechanism are identified based on the literature and experience from Greek lignite surface mines. Two finite element (Plaxis2D and RS2) and one limit equilibrium (Slide2) geotechnical software are compared for the SF calculation. Additionally, slopes with benches, as in engineering practice, and without benches, as in slope stability analysis and design, are compared.
The central part of this work includes evaluating the effect of the geometrical and geotechnical parameters and pore water pressures on slope stability. Overall, the present analysis on the influence of the critical parameters on the stability allows for the preliminary stability evaluation of a lignite/coal mine and the appropriate quantification of each parameter's effect through the probability of slope failure.
2 Development of the reference model 2.1 Soil stratigraphy and properties Figure 1 illustrates a characteristic geometry with the weak zone crossing at the bottom of the excavation, lying above the bedrock and below the overburden soil. In practice, the prevalent weak zone can be found close to the toe, but sometimes also below or above the bottom of the excavation. Herein, the weak zone is assumed to cross the slope's surface at the toe of the slope (i.e. at the bottom of the excavation) as this is most critical for stability analysis. Suppose the weak zone crosses the stratigraphy beneath the excavation's bottom; in that case, the failure surface will go deeper to pass through the weak zone and have a larger part in the overburden soil of higher strength, leading to a more stable condition. In contrast, a failure surface above the excavation could be interpreted as an equivalent slope of smaller height; thus, a higher SF is expected.
The overburden soil consists of layers of sterile materials (for instance, in Greek lignite mines, these are mostly stiff to hard clays and marls) and lignite seams. As shown in the sequel, the overburden soil's properties are less crucial than the properties of the weak zone. Additionally, these layers could be characterised for slope stability in terms of their shear strength only, and they frequently have similar shear strength. As a result, the overburden soil is simulated as a homogeneous material ( Fig. 1). At the bottom of the excavation, a bedrock formation with high strength and stiffness is commonly encountered. This region (underneath the weak zone) does not directly contribute to the slope stability analysis if a weak zone is present.
Geometry and soil properties critically affect the stability of the slope. Geometrical and geotechnical parameters that are expected to be important are the height (H) and the inclination (b) of the slope, the inclination (b z ) and thickness (d z ) of the weak zone, and the soil's strength for both the overburden soil (/, c) and the weak zone (/ z , c z ). Slope stability literature recurrently discusses the effect of the parameters H, b, / and c (e.g. Rahardjo et al. 2007). Leonardos (2004) identified the parameters of the weak zone (b z , / z and c z ) as critical, and Kavvadas et al. (2020) tentatively analysed them for Greek lignite surface mines, based on a simple analytical model.
Typical values for lignite mines are considered in this work to specify these geometrical and Fig. 1 Simplified stratigraphy of a surface coal mine slope with a weak zone at the bottom of the excavation geotechnical parameters. Deep slopes characterise excavations as they have become exceptionally large, reaching depths of 200 m. Moreover, inclinations vary between 1:7-1:4 (8°-14°) for the mines' sustainable operation as large bucket-wheel excavators are usually employed to exploit such deep excavations (Kavvadas et al. 2020). According to Kavvadas et al. (2020), the weak zones' thickness varies between a few millimetres and several tens of centimetres. At the same time, its inclination may differ from the horizontal up to 6°, based on reported failure incidents (e.g. Kavvadas et al. 2013;Steiakakis et al. 2018).
In a recent research study, Theocharis et al. (2021) organised and presented a large geotechnical database collected from various boreholes in several Greek lignite mines. Extensive laboratory results, including index and classification properties, triaxial, ring shear and direct shear tests, were statistically analysed. The authors reported that fine-grained soils dominate Greek lignite mines, and shear strength parameters (friction angle and cohesion) of the overburden soil formations often vary in a narrow range. In that vein, considering the overburden soil formations as a single material for slope stability analyses is a reasonable assumption. However, its strength parameters may vary significantly between different mines.
The mean values of friction angle and cohesion were equal to 28°and 185 kPa, while the standard deviations were 7°and 145 kPa, respectively. In many cases, results indicated specimens with highly diminished residual strength, indicating the weak zone's presence. This zone is characterised by a mean residual friction angle of 10°and a standard deviation of 2.6°, and its residual cohesion lies between 0 kPa and 5 kPa. The mean value of soils' unit weight was 17 kN/m 3 , and its standard deviation 1.6 kN/m 3 .

Background information
The Finite Element Method (FEM) is primarily used in this work, being the most versatile and widely used numerical method in geotechnical engineering. FEM need to be combined with the Shear Strength Reduction (SSR) technique to perform ultimate limit state analysis. In the SSR type of analysis, herein based on the Mohr-Coulomb failure criterion, the safety factor is calculated by progressively reducing the material's shear strength to bring the slope to the verge of failure.
A reference model was established with a 200 m slope height, the maximum reported in the literature, and an inclination of 10°, a mild inclination for a high slope. The weak zone was simulated as a 5 m thick horizontal soil layer. For the overburden soil, the mean values of Theocharis et al. (2021) for the friction angle (/ = 28°) and cohesion (c = 185 kPa) were used. The weak zone was assigned significantly reduced strength parameters with respect to the overburden material, presumed to be the zone's residual strength; as mentioned above, the weak zone's strength is considered to be at its residual state due to the immense shear strains developed on this thin interface during ground movements. Four residual friction angles of 5°, 8°, 12°a nd 15°and a residual cohesion of 5 kPa were implemented; these values were obtained from the laboratory data presented in Theocharis et al. (2021). Young's modulus for the overburden soil and the weak zone was 50 MPa; however, this parameter has a minor impact on slope stability. Additionally, the strength and Young's modulus of the bedrock do not affect the results.
For the moist unit weight, the mean value of Theocharis et al. (2021) c = 17kN/m 3 , typical for the materials analysed, was used, while Poisson's ratio and dilation angle were m = 0.3 and w = 0°, respectively. These values are typical for normally consolidated fine-grained soils. Additionally, the dilation angle has a minor but systematic effect on the SF and the failure mechanism; less than 5% difference on SF between w = / and w = 0° (Liu et al. 2015;Tutluoglu et al. 2011). However, the flow rule has a significant impact on steep slopes (inclination angle more than 40°) of high friction angle materials (e.g. / [ 40°) (Tschuchnigg et al. 2015a, b). In the present work, a non-associated flow rule with w = 0°has been employed as mild inclinations are considered, the friction angle is less than 35°, and the dilation's angle impact on the stability calculations would be minor. Last, the Mohr-Coulomb elastic-perfectly plastic constitutive model was employed. In summary, the geometrical and geotechnical parameters of the reference model are shown in Table 1.
Two-dimensional (2D) plane strain and drained conditions were used. Vertical boundaries were located at a 2H distance from the slope's crest and toe, respectively, and the horizontal bottom boundary at a distance H from the slope's toe to ensure the minimization of the boundary conditions effect, where H stands for the slope's overall height (Fig. 1). A very fine discretisation was used for FEM, and mesh density was increased in some critical areas (Fig. 2).
Two widely available software, Plaxis2D (Plaxis 2020) and RS2 (Rocscience 2020), were used and compared. In Plaxis2D, 6-node or 15-node triangular elements are available, while in RS2, the selection is restricted to 3 and 6 nodes for triangular, and 4 and 8 nodes for rectangular elements. The limited number of Weak zone's thickness d z (m) 5 0.5-5  nodes in RS2 may lead to inaccuracies, especially on initial stresses, the model boundaries, and high stress gradient areas. Furthermore, RS2 results varied more with mesh density and demanded caution in the discretisation and meshing procedures. Analyses in both FEM software were conducted for 6-noded triangular elements and similar mesh densities, and they were further validated with 15-node triangular elements. All the reference model analyses presented identical failure mechanisms and SF values (differences less than 1%). Figure 3 presents the shear strains that denote the failure surface position for two indicative cases conducted with RS2 and Plaxis2D.
Finally, it is noteworthy that in FEM, the weak zone can be simulated in two distinct ways: as a thin soil layer or as a ''zero-thickness'' interface element. For the scope of the present work, as the layer's thickness is considered a parameter under investigation, only the layer-case approach is implemented. However, preliminary analyses applying an interface element with similar strength parameters led to identical results.

Finite element versus limit equilibrium method
The limit equilibrium method (LEM) and the finite element method (FEM) are two of the most commonly used methods in evaluating the stability of excavated slopes (e.g. Carranza-Torres et al. 2017;Khandelwal et al. 2015;Wang et al. 2019). In this section, FEM-SSR results are compared with LEM results as a reference comparison for the weak zone soil profile. LEM leads the slope stability analysis over the last decades providing simplicity and speed in the calculation process. However, LEM does not consider the stress-strain relationship of the materials. Simultaneously, the shape and location of the failure surface should be defined in advance, mainly for non-circular failures as those of lignite mines.
On the contrary, FEM does not require any assumption on the failure surface but identifies it based on soil's stress-strain response. Nonetheless, FEM demands higher computational time and a more detailed description of soil behaviour. The two FEM (Plaxis2D and RS2) software used in the previous section were compared with the LEM software Slide2 (Rocscience 2019) based on the stratigraphy of Fig. 1. Comparison is based primarily on the SF and the failure surface.
For LEM analyses, Spencer's method was implemented, combined with the Mohr-Coulomb failure criterion. All LEM software demand the assumption of the shape (circular or non-circular) and, in some cases, also the location of the failure surface to provide reliable results, depending on the search method. Figure 4a and b illustrate two LEM analyses (Slide2) for the reference case and / z = 5°, assuming a noncircular and a circular type of failure surface, respectively. The two analyses provided different failure surfaces and SFs, even though the geometry and the geotechnical parameters of the materials are the same.
To simulate the complex failure mechanism and the non-circular failure surface, the block search method (Slide2) was employed. In that method, various noncircular failure surfaces are employed that necessarily pass through a predefined area of the model. LEM results were practically identical to FEM results (differences less than 2%).

The role of benches in the overall stability
A typical slope on surface mining consists of several benches, depending on the excavation method and equipment used. The exact number and geometry of benches vary in a wide range, making it difficult to study the slope's stability systematically. As a result, benches are often neglected without justification, even though the consideration of a slope consisting of several benches is more realistic. In that vein, it is essential to quantify the error for such an assumption.
Different slope geometries and benches' configurations were implemented for comparison in Plaxis2D. The slope's geometrical parameters varied within their reported range (Table 1). The reference parameters (Table 1, Reference values) describe the strength and stiffness of the soils. Moreover, there is a need to assign a linear elastic material to the benches to avoid minor regional bench failures, insignificant on the overall safety analysis. Under this perspective, the weak zone was 1 m thick to prevent regional failure surfaces at the benches closest to the toe. As for the weak zone's strength, the comparison is conducted for residual friction angles between 5°and 15°. Different configurations were created with 1 to 9 benches and heights between 20 and 100 m; the bench face inclination varied from 12°to 45°, and the bench width from 100 to 200 m. The benches' number for each configuration depended on the geometry and the selected range of the benches' parameters. The configurations' characteristics are summarised in Table 2. Overall, 21 different configurations were analysed, and Fig. 5 shows two indicative ones.
The failure surface was similar for all models, regardless of the benches' number and geometry. Moreover, the models with benches provided higher SF values in all cases, and the analysis without benches was always more conservative by 3%-11%. Similarly to the present results, Verma et al. (2013) conducted FEM analyses for an internal dump in a coalfield in India and concluded that omitting benches was always more conservative with an SF error of 3-6%.
This work aims to systematically examine the influence of the weak zone on slope stability of lignite mines, and one particular configuration with benches is not of interest. Thus, for the analysis that follows, benches are ignored, noticing that the results are on the conservative side and the error is less than 11%.   3 Effect of geometry and geotechnical properties on slope stability In this section, the effect of the geometry and the geotechnical parameters on the safety factor (SF) was investigated through the FEM (Plaxis2D). The critical parameters, as identified in Sect. 2, were investigated within their typical range in Greek lignite mines (Table 1). Nonetheless, this particular range is rather broad, and lignite mines in other countries have similar properties (e.g. Bednarczyk 2017; Ulusay et al. 2014). Each set of analysis was conducted by varying one critical parameter, while the rest of the parameters kept their reference values (see Table 1); the thickness of the weak zone was 1 m. FEM results indicate that the weak zone's shear strength, particularly its friction angle, is a key parameter to slope stability. The range of the weak zone's cohesion is narrow (0-5 kPa), and when cohesion decreases to near zero, numerical issues may arise; thus, it is not investigated in this work. As shown in Fig. 6, the weak zone's friction angle impact on the slope's stability is profound; decreasing / z from 15°to 5°reduces SF between 45% and 58%, regardless of the weak zone's inclination. The maximum reduction was encountered for the case of the steepest unfavourable inclination of the weak zone (b z = 6°), where SF drops from 2.00 to 0.85. Furthermore, notice that even a small alteration of / z by 3°can change SF up to 30%; for instance, for b z = 6°SF drops from 1.22 (/ z = 8°) to 0.85 (/ z = 5°).
The tremendous effect of the weak zone on the SF is even more highlighted if the weak zone is assigned strength parameters equal to the overburden soil's, i.e. / z = 28°and c z = 185 kPa. Then, the slope is composed of one soil material and a bedrock formation, and for the reference parameters, SF increases to 4.18. At the same time, the failure surface is almost circular, with a small fraction being affected by the bedrock formation (Fig. 7b). If the bedrock is also considered to have the same properties as the soil, the SF drops slightly to 3.99, and the failure surface is circular (Fig. 7a).
However, the weak zone's presence dictates a composite failure surface and a considerable decrease of SF, from 4.18 to 2.47 (/ z = 15°), and 1.35 (/ z = 5°). Sliding occurs along the horizontal/subhorizontal weak zone and reaches the slope's top with a planar transition. For lower weak zone's friction angles, the inclination of the planar transition gets steeper, and the transition point gets closer to the crest ( Fig. 7c and d). Notice that the shear strains inside the weak zone are very high but do not appear in Fig. 7c or d due to the very small thickness of the zone respecting the overall geometry (e.g. see the detail of Fig. 7d).
Both the strength and the inclination of the weak zone are critical parameters (Fig. 6). On the contrary, the effect of the weak zone's thickness was found insignificant (Fig. 8). Although d z varies in a wide range of 0.5-5 m (d z /H = 0.25-2.5%, for H = 200 m), the difference in SF is less than 1.5%. Notice that for d z B 1 m, the required computational time was remarkably higher.
Design of slopes in surface mining is usually performed based on empirical safety factors (deterministic approach). Nonetheless, safety factors do not accurately reflect the probability of failure of a slope. Regardless of the type of geotechnical structure, safety factors not only tell little about the possibility that a failure may occur, but they can even be misleading because of the ambiguity between them and the underlying level of risk (Whitman 1984). For instance, large safety factors do not necessarily reflect a lower level of risk because the presence of considerable uncertainties can negate their effect. Probabilistic slope stability analysis offers the framework for Fig. 6 SF versus weak zone's friction angle / z for different weak zone inclinations b z addressing the above shortcoming by explicitly accounting for uncertainty sources in the computation of the probability of failure and its complement reliability. This type of analysis greatly facilitates the use of reliability as a criterion of design optimization that plays a fundamental role in mining (Zevgolis et al. 2018).
In this work, Rosenblueth's point estimate method is employed to calculate the probability of failure (P F ) and assess each parameter's importance on slope stability (Rosenblueth 1975). This method is widely used in geotechnical practice due to its simplicity and accuracy (Baecher and Christian 2003). The low-order moments of the dependent random variable (herein the safety factor) can be approximated based on the independent random variable's low-order moments. Since the weak zone strength largely governs the slope's stability (see also Fig. 6), its friction angle / z is considered as a random variable with mean value l[/ z ] = 10°. Two scenarios are considered in terms of the uncertainty of / z through two coefficient of variations (COVs), i.e. COV[/ z ] = 20% and 30%. Both values are typical for friction angle (e.g. see table S1 in Zevgolis et al. 2021). Therefore, SF is a function of one independent variable; for comparison purposes, SF is assumed to be both normally and lognormally distributed. Once the mean value and standard deviation of SF are computed, the reliability index may also be computed. Ultimately, the probability of failure is assessed through the normally distributed cumulative function of the reliability index. Figure 9 shows the relation between P F and SF for b z varying between -6°and 6°. The Figure clearly illustrates that SF, by itself, cannot reflect the risk level, given its highly non-linear relation with P F . This non-linear relation indicates that even a small decrease in SF might significantly increase P F and the Fig. 7 Failure surface for a a homogeneous slope, b a slope with a bedrock formation and no presence of a weak zone, c a slope with a weak zone of / z = 15°, and d a slope with a weak zone of / z = 5°considering the reference case corresponding risk level. For example, for the given range of b z , SF varies in a wide range, and results indicate that the slope is relatively stable in all cases (minimum SF equals 1.45). Nonetheless, for COV[/ z ] = 30%, a normally distributed SF with a mean value equal to 1.45 corresponds to P F = 9.9%, while an SF = 1.68 corresponds to P F = 2.3%. For higher SFs, P F is diminished, as for SF = 2.42, P F equals 0.00005%. It is noteworthy that this trend is more significant for a lognormally distributed SF and for COV[/ z ] = 20%.
In the following sections, the significance of parameters /, c, H, b and b z on the probability of failure (based on treating / z as a random variable) is parametrically evaluated. The equivalent SFs for deterministic analysis (not presented) are always higher than 1. Results for COV[/ z ] = 30% are always more conservative than for COV[/ z ] = 20%. Besides, the normal distribution typically provides more conservative results than the lognormal distribution, as expected for low probabilities of failure, approximately less than 15% (Baecher and Christian 2003). Because the two distributions lead to similar conclusions, the authors interpret the following Figures, focusing only on the normal distribution results.
Notice that the vertical axes are in logarithmic scale, extending from a maximum value of 1 (100%) up to a minimum value of 10 -7 (10 -5 %). This range is wide, and small P F values close to the minimum are typically insignificant in engineering practice. However, this range allows for a better comparison of the effect of different parameters, noticing that the trends are as important in this analysis as P F values.

Overburden soil's strength
The overburden soil's friction angle took the values 21°, 28°and 35°(i.e. the mean value plus and minus one standard deviation), while the cohesion 40 kPa and 185 kPa (i.e. the mean value minus one standard deviation) (Fig. 10). In Fig. 10a, for COV[/ z ] = 30%, a friction angle increase from 21°to 35°reduces P F from 2.4% to 0.3%; for COV[/ z ] = 20%, reduces P F from 0.1% to 0.002%. In Fig. 10b, cohesion increases to 185 kPa, and P F varies between 0.7% and 0.06% for COV[/ z ] = 30%, and between 0.009% and 0.00007% for COV[/ z ] = 20%. The comparison between Figs. 10a and b shows the influence of cohesion. Increasing cohesion from 40 to 185 kPa leads to a maximum P F decrease for / = 21°, with P F decreasing from 2.4% to 0.7% (for COV[/ z ] = 30%), and from 0.1% to 0.009% (for COV[/ z ] = 20%).  Figure 11 illustrates the effect of the slope's geometry on P F ; Fig. 11a presents the influence of parameter H on P F, and Fig. 11b of parameter b on P F . As expected, increasing the slope height increases P F (Fig. 11a). The increase of the slope height causes an increase of the sliding mass' weight but also contributes to the increase of the shearing resistance along the base of the sliding mass. This double contribution of H on slope stability leads to the non-linear relation of height with P F (and SF). Moreover, as the slope height is increased from 50 to 200 m, the planar transition of the weak zone to the ground surface becomes larger; thus, the total sliding surface has a more significant part inside the overburden soil of higher strength (than the weak zone). As a result, the overall stability is enhanced compared to the lower slope heights, and the plotted curves around 100 m-200 m are milder than  for lower heights 50 m-100 m. It is noteworthy that the analysis considering COV[/ z ] = 20% and a lognormal distribution for P F is not presented herein, as all P F s were very low (\ 0.00001%). Figure 11b presents an almost linear relation between P F in log scale and b, with a 2°increase in the slope's inclination resulting in a P F increase of approximately 50%, for COV[/ z ] = 30%, and 80%, for COV[/ z ] = 20%. Comparing the reference value (b = 10°) and a steeper slope of b = 14°, P F rises from 0.2 to 0.9%, for COV[/ z ] = 30%, and from 0.001% to 0.02%, for COV[/ z ] = 20%. Figure 12 illustrates the variance of P F with b z . Indicatively, for COV[/ z ] = 30% and b z = -6°(a dip of strata in favour of stability), P F equals 0.00005%, while the respective value for b z = 6°is 9.9%. The P F range is broader for the lognormal distribution, with P F ranging from 7.5% to a value less than 0.00001%. In any case, the results indicate the critical role of the dip of the weak zone. Figure 13 illustrates the shear strains and the failure surface for four indicative configurations with different geometric parameters; Fig. 13a shows the reference model, and  (Fig. 13a). Considering slope height reduced by 50 m, SF is increased by 4% (SF = 2.01), but the failure surface remains practically identical (Fig. 13b). The respective P F is reduced to 0.13% (41% decrease). The decrease of the slope's inclination b by 2°leads to an SF increase of 15% (SF = 2.24) and a P F decrease of 50% (0.11%), with the failure surface moving slightly to the left (Fig. 13c). Last, in Fig. 13d, a slight deviation of the weak zone's inclination (b z = -3°) results in an SF rise of 13% (SF = 2.19) but also in a huge decrease of P F by 96% (0.008%). In this specific case, even if the effect of b on SF is marginally greater than the effect of b z , their vast P F difference indicates that b z is crucial.

Weak zone's inclination
All analyses present similar failure surfaces; they have a part in the weak zone and reach the ground surface with a planar transition of approximately 45°. Moreover, a reflection of the shear band that identifies the transition to the ground surface is always present. It is concluded that, for dry conditions, the most critical parameters are those of the weak zone, i.e. the weak zone's inclination and friction angle. Specifically, the uncertainty of the weak zone's strength, considered a random variable and investigated through two different COVs, tremendously affects the P F . For all simulations probability of failure changed almost one order of magnitude for large probabilities and more than that for lower P F s when this COV increased from 20 to 30%. Finally, it is noteworthy that the weak zone characteristics are also challenging to determine in practice, making them even more important for stability analysis.

The effect of water on slope stability
In this section, groundwater's influence is considered through pore water pressure using the simplified approach of r u coefficient (Bishop and Morgenstern 1960). This method is frequently used for slope stability in engineering practice, as the groundwater regime is often difficult or impossible to identify. The pore water pressure is calculated as: where r u is the coefficient of pore pressure (ranging from 0, zero pore pressure, to 1, zero effective stress), c is the soil unit weight, and z is the depth below the soil surface. Furthermore, if the saturated soil unit weight is 20 kN/m 3 , and the water unit weight is 10 kN/m 3 , r u = 0.5 is equivalent to a phreatic water table lying at the soil surface. The r u method cannot be implemented in Plaxis2D, and thus, the analyses that follow were carried out with RS2. Note that RS2 results for dry conditions (r u = 0) are practically identical with the values derived from Plaxis2D. Different pore water pressures were simulated, using r u equal to 0 (dry conditions), 0.2 and 0.4, with the soil's moist unit weight presumed to be 17 kN/m 3 . It should be noted that all FEM analyses performed for r u = 0.4 and / z B 8°led to SF values less than or very close to 1. Figures 14, 15, 16 illustrate the variance of P F with the parameters /, c, H, b and b z, and for three different r u (considering l[/ z ] = 10°and COV[/ z ] = 20/30%). The probability of failure was computed considering SF as normally distributed. Increasing r u results in a substantial P F increase as expected, with r u = 0.4 leading to immense probabilities of failure. For example, in Fig. 14a, for the reference model (/ = 28°) and COV[/ z ] = 30%, increasing r u from 0 to 0.4 increases P F from 0.2 to 33.2%.
Extreme values for the P F are encountered for the highest values of b and b z . In Fig. 15b, the combination of b = 14°and r u = 0.4 results in P F = 89.6% (COV[/ z ] = 20%), while in Fig. 16, P F = 89.6% for b z = 6°. It is noteworthy that in these analyses, the mean values of the SF distributions are less than 1 and almost identical to each other. In such cases, smaller COVs will provide larger probabilities of failure, as a greater part of the SF distribution is less than 1. For instance, for COV[/ z ] = 30%, the respective results of P F are smaller, as they equal 79.8% and 79.6% for b = 14°and b z = 6°, respectively.
Overall, it is concluded that the water pressure's effect on stability is vital. For the reference model, r u values of 0, 0.2 and 0.4 correspond to P F 0.2%, 2.8% and 33.2%, respectively. Although the value of r u = 0.4 might represent high pore pressures (r u = 0.5 is equivalent to a phreatic water table lying at the soil surface), it further supports the groundwater's crucial role. Immense probabilities of failure were encountered, even for the case that r u = 0.2. As a result, groundwater conditions are identified as one of the most critical parameters on slope stability, along with b z and / z .

Conclusions
In the present work, the slope stability of surface mines has been numerically investigated within a probabilistic framework. The focus was on the weak zone often encountered near the bottom of open-pit lignite and coal mines excavations. Several numerical issues on analysis and design were addressed, and the critical parameters of stability were identified and discussed. Their effect on SF was numerically quantified and evaluated, primarily through the FEM combined with SSR. The point estimate method was implemented to evaluate their significance through the probability of slope failure.
Initially, two FEM (Plaxis2D and RS2) and one LEM (Slide2) geotechnical software were compared, with slope stability results being practically identical. An excellent agreement was demonstrated on the failure mechanism, and the difference between the SFs was less than 2%.
Moreover, slopes with benches, as typically in mines, and without benches, as typically in slope stability analysis, were compared. Those without benches provided similar failure mechanisms and a marginally lower SF, compared to a slope with benches and equivalent geometrical characteristics; the difference in SF was always less than 11%. Thus, benches' influence can be ignored, noticing that the results are always on the safe side.
Furthermore, the effect of geometrical and geotechnical parameters on slope stability was examined. FEM analysis, literature and experience all agree that the most critical parameter is the weak zone's strength; decreasing friction angle / z from 15°to 5°leads to an SF decrease between 45% and 58%, for the weak zone's inclination ranging between 0°-6°. The point estimate method was employed to calculate the probability of failure and assess each parameter's importance; / z was considered a random variable since the weak zone strength largely governs the slope's stability. The uncertainty of the weak zone's strength, investigated through two different COVs, affects the P F dramatically. For all simulations probability of failure changed almost one order of magnitude for large P F s and more than that for lower P F s when COV[/ z ] increased from 20 to 30%. The inclination of the weak zone was treated parametrically and is also crucial, as changing b z from 0°to 6°r esulted in a P F increase from 0.2% to 9.3%. Finally, the layer's thickness is insignificant for the overall stability.
The soil strength of the overburden material and the height and the inclination of the corresponding slope are essential in slope stability. The influence of /, c, H and b on the P F for this case is notable but not as crucial as the effect of b z . Moreover, it was shown that small variations of their values might cause noticeable deviations to SF but do not crucially affect P F . As a result, it is concluded that the most critical parameters are those of the weak zone, i.e. the weak zone's inclination and strength. It is noteworthy that the weak zone characteristics are also challenging to determine, making them even more critical for stability analysis.
Groundwater was considered through the pore water pressure using the simplified approach of r u coefficieny. Water pressure's effect is vital as P F vastly increases with the increase of r u . The highest P F s are encountered when b and b z increased from their initial values, 10°and 0°, to 14°and 6°r espectively, as P F equals 89.6% for r u = 0.4. The vast increase of the P F indicates that even very stable slopes in dry conditions can reach failure if the pore water pressure increases. As a result, groundwater conditions are identified as one of the most critical parameters, along with b z and / z . Overall, the present work has employed numerical analysis within a probabilistic framework to systematically investigate the critical parameters' effect on slope stability. The crucial role of the weak zone's strength and inclination was identified, and their effect on slope stability was quantified through the probability of slope's failure. At the same time, the present work included the groundwater pore pressures in a practical and standard way, quantifying the crucial role of water on stability.