The Setting of Strength Parameters in Stability Analysis of Open-Pit Slope Using the Random Set Method in the Bełchat ó w Lignite Mine, Central Poland

: The slopes of open-pit mines are often at risk of failure. To identify this hazard, stability analyses are performed. An important element of these stability analyses is the reliable selection of input parameter values for the calculations. This selection is difﬁcult because the slopes of the open pit are disturbed by mining activities. In such conditions, rheological processes, intensiﬁed by weathering, develop in open-pit slopes. This study is aimed at setting the strength parameters for the stability analysis of open-pit slopes with a developed slide process, using the random set method. The study was performed on the example of the open pit of the Bełchat ó w lignite mine, central Poland. A four-stage methodology, according to the random set method, was proposed. The methodology covered the following: site investigation, sensitivity analyses, shear strength reduction (SSR) analyses using numerical calculations, and probability analyses of the factor of safety (FoS) calculation results. The setting of the input parameters took into account the peak and residual strength parameters for each lithological unit in the physical model of the open-pit slope. Samples for laboratory tests were taken from the cores of nine test boreholes. The sensitivity analysis included all peak and residual strength parameters for each lithological unit in the body. As a result of the sensitivity analysis, speciﬁc strength parameters were adopted that would have a great impact upon the results of the calculations. Selected sets of parameter values were then used for the FoS calculations. The resultant FoS values revealed the probable slide planes. The positions of the slide planes were consistent with the interpreted slide surfaces based on the control boreholes and terrain observations. Knowledge of the slide planes positions and the values of the strength parameters enabled the designing of a securing approach for this landslide, and the taking of preventive measures to reduce this risk.


Introduction
The slopes of exploited open-pit mines of different mineral resources are often at risk of failure (e.g., [1][2][3][4][5][6]). To more accurately recognise this risk, various types of stability analyses are performed, using, in the broadest sense, deterministic and probabilistic methods (e.g., [4][5][6][7]). An important element of the deterministic stability analysis is a reliable assessment of the strength parameters to be used in the factor of safety (FoS) calculations.
The strength parameters are important in calculating the actual FoS, representing the slope's stability. They are also used to forecast the critical FoS when the slope loses its stability. The strength parameters can be obtained from test results, either directly or through correlation, theory or empiricism, and other data (EN 1997(EN -1:2004 [8].

calculations.
The strength parameters are important in calculating the actual FoS, representing the slope's stability. They are also used to forecast the critical FoS when the slope loses its stability. The strength parameters can be obtained from test results, either directly or through correlation, theory or empiricism, and other data (EN 1997(EN -1:2004 [8].
The proper setting of the strength parameters is extremely difficult in an open-pit slope, due to ground degradation resulting from mining activities. In such conditions, a rheological process accelerated by the weathering process develops more intensively. Most frequently, complex groundwater conditions exist. As a result, a large variation in the values of ground parameters develops. These reasons are the cause of practical difficulties in setting the characteristic strength parameters to analyse the open-pit slope's current stability.
It is advantageous to describe the great variability of the strength parameters, using peak and residual parameters. This requires knowledge of a stress-strain relationship [9,10]. The typical stress-strain relationship presenting a reduction in shear strength from peak to residual values in the case of clay is shown in Figure 1. The stress-displacement curves for over-consolidated (OC) and normal consolidated (NC) clay are shown on the left-hand side. On the right-hand side, the shear strength-effective normal stress lines with failure lines are presented. Failure lines of shear strength and effective friction angles for peak parameters of over-consolidated (φ'OC) and normal consolidated (φ'NC) clay and residual parameters (φ'r) describe the variability of the strength parameter values. Peak strength is a condition that occurs at the beginning of local failure along the weakest surface. High plasticity soils indicate a low value of shearing displacements at the peak value ( Figure 1). This failure occurs in the stress-softening stage when the shear surface is forming, e.g., [9]. When the state of particle clusters in the soil structure exceeds the peak strength, the destructing process starts. Unbalanced stresses may affect neighbouring particles in the destructing process, eventually leading to a residual strength condition. Usually, bonds are destroyed, and the strong reorientation of soil particles occurs parallel to the direction of shearing when the residual strength is reached [11,12]. In this new condition, both shear stress and changes in the displacement become constant ( Figure 1). Romero et al. [13] stated that the residual shear strength in soils depends on the level of normal stress, soil grading, bulky and platy-like particle mineralogy, dolomite and/or calcite content, rate of shearing, and the chemical interaction between water chemistry and soil mineralogy.
There are two approaches used to calculate the FoS: the limit equilibrium method (LEM) and the shear strength reduction (SSR) method [14]. The LEM is used to define the FoS as a ratio of nominal capacity to the system's demand, and it is a force and/or moment equilibrium calculation. The SSR is used to determine the strength reduction factor Peak strength is a condition that occurs at the beginning of local failure along the weakest surface. High plasticity soils indicate a low value of shearing displacements at the peak value ( Figure 1). This failure occurs in the stress-softening stage when the shear surface is forming, e.g., [9]. When the state of particle clusters in the soil structure exceeds the peak strength, the destructing process starts. Unbalanced stresses may affect neighbouring particles in the destructing process, eventually leading to a residual strength condition. Usually, bonds are destroyed, and the strong reorientation of soil particles occurs parallel to the direction of shearing when the residual strength is reached [11,12]. In this new condition, both shear stress and changes in the displacement become constant (Figure 1). Romero et al. [13] stated that the residual shear strength in soils depends on the level of normal stress, soil grading, bulky and platy-like particle mineralogy, dolomite and/or calcite content, rate of shearing, and the chemical interaction between water chemistry and soil mineralogy.
There are two approaches used to calculate the FoS: the limit equilibrium method (LEM) and the shear strength reduction (SSR) method [14]. The LEM is used to define the FoS as a ratio of nominal capacity to the system's demand, and it is a force and/or moment equilibrium calculation. The SSR is used to determine the strength reduction factor (SRF), using a stress/strain analysis in numerical calculations. In the open-pit slope stability analysis, both LEM and SSR methods can be applied.
Many studies have been performed to analyse the stability of slopes in open-pit mines. Ural and Yuksel [15] analysed slope instability, using the Spencer-Wright limit equilibrium method of SIROQUANTM in the Kislakoy open-pit mine, Turkey. They deter-mined the most representative shear strength parameters from typical shear stress-shear strain obtained from the direct shear tests. They also performed a back analysis to obtain critical strength parameters from known, or assumed, failure surfaces independent of a priori knowledge of the shear strength. They noted that the contact zone presents a shear strength at, or approaching, the residual condition. The back-calculated parameter values were in good agreement with the laboratory-derived residual shear strength parameter of the contact zone.
Soren et al. [16] used the SSR method in the FLAC/Slope numerical calculations to analyse the instability of a theoretical slope in a typical open-pit mine. They showed how to find the factor of safety of an open-pit slope by initiating a systematic reduction sequence for the available shear strength parameters c and ϕ to the slope failure. The reduction values of the shear strength parameters of cohesion ctrial and internal friction angle ϕtrial are defined as the following: SRF is a strength reduction factor, which is equal to the factor of safety [16]. Tao et al. [17] performed a stability analysis in the Niukutou open-pit mine in Qinghai Province, China, using the simplified Bishop limit equilibrium method of Slide 6.0 and the shear strength reduction method of Phase2 v.7.0. They assumed characteristic strength parameters of the rock mass based on the Chinese standard recommended methods, using the Hoek-Brown strength criterion. The two methods both indicated that the factor of safety was consistent and much smaller than 1.0.
Bednarczyk [6] analysed the slope stability of a lignite open-pit mine in central Poland. He performed SSR analysis of FLAC v. 7.0. Due to the lack of appropriate strength parameter values, he assumed characteristic and design values of effective strength parameters from the so-called analogue method. Therefore, the values of the parameters used for the numerical calculations were assumed from lignite open-pit mines located nearby. This allowed the prediction of the slide surface location, the greatest deformation zones, and the velocity of soil displacement within the analysed slopes.
Fang et al. [18] studied an unstable slope in the Buzhaoba open-pit mine, China. This slope was at critical stability mainly due to mining operations and rheological processes accelerated by weathering. They implemented an analytical SSR method to calculate the FoS. They assumed characteristic strength parameters of the ground obtained from standard tests. In the calculation scheme of critical FoS, they decreased strength parameters to residual values. They concluded that the greater the residual strength of the tested material, the more slowly the FoS decreases.
One of the approaches suitable for the selection of strength parameters for both LEM and SSR methods can be the random set method (RSM) [19]. Specifically, it is a theory of setvalued stochastic processes. As coarse data represent a case of imprecise observation, the concept of random sets can be extended to random fuzzy sets to model perception-based information [20]. Random sets theory serves as a bridge between several useful measures of uncertainty in decision analysis. RSM is useful in conditions of a large number and great variability of parameters [21][22][23][24]. It allows analysis of the impact of changes in the value of strength parameters on FoS with a specific probability measure assignment [22]. The advantage of this method with regard to the stability analysis is the recognition of shallow and deep-seated slide planes assigned to specifically the most influential parameter sets.
This study aimed to set the strength parameters for the stability analysis of an open-pit slope with a developed slide process, using the random set method. The SSR approach was proposed in this study to enhance the prediction of FoS by numerical analysis, using FLAC v. 7.0 [25]. We analysed a large landslide that occurred on open-pit slopes to identify its mechanism, particularly the location of possible slide surfaces. Knowing the slide surface location and values of strength parameters was necessary for designing a protection system for landslides and taking other preventive measures to reduce this risk. The study was performed on the example of the open pit of the Bełchatów lignite mine, central Poland ( Figure 2). The stability analysis was calculated for an 18S landslide that occurred in the northern slope of the Bełchatów open pit (Figures 2 and 3). A landslide developed in a Miocene formation on the bedrock of Jurassic formation. Highly different values of strength parameters characterised the landslide colluvium. The soil strength parameters were tested by the laboratory method, using the direct shear apparatus with the Polish standard (PN-B-04481:1988). Peak and residual parameters were determined, based on the stress-strain relationship [12,26,27].
This study aimed to set the strength parameters for the stability analysis of an open-pit slope with a developed slide process, using the random set method. The SSR approach was proposed in this study to enhance the prediction of FoS by numerical analysis, using FLAC v. 7.0 [25]. We analysed a large landslide that occurred on open-pit slopes to identify its mechanism, particularly the location of possible slide surfaces. Knowing the slide surface location and values of strength parameters was necessary for designing a protection system for landslides and taking other preventive measures to reduce this risk.
The study was performed on the example of the open pit of the Bełchatów lignite mine, central Poland ( Figure 2). The stability analysis was calculated for an 18S landslide that occurred in the northern slope of the Bełchatów open pit (Figures 2 and 3). A landslide developed in a Miocene formation on the bedrock of Jurassic formation. Highly different values of strength parameters characterised the landslide colluvium. The soil strength parameters were tested by the laboratory method, using the direct shear apparatus with the Polish standard (PN-B-04481:1988). Peak and residual parameters were determined, based on the stress-strain relationship [12,26,27].

Geological and Mining Setting
The research area is located within the Bełchatów lignite mine, approximately 180 km southwest of Warsaw ( Figure 2). The lignite deposit is located in the Kleszczów trough, the deepest neotectonic subsidence zone in the Polish Lowlands. Faults heavily

Geological and Mining Setting
The research area is located within the Bełchatów lignite mine, approximately 180 km southwest of Warsaw ( Figure 2). The lignite deposit is located in the Kleszczów trough, the deepest neotectonic subsidence zone in the Polish Lowlands. Faults heavily cut the Bełchatów mine deposit in the NW-SE and WSW-ENE directions and subordinate WNW-ESE and SW-NE (Figures 2 and 3). The geological units of up to 180 m in depth ( Figure 4) include the Mesozoic bedrock and the following four Miocene formations in the overburden: sand-lignite, clay-lignite-silt, clay, and clay-sand.  Exploitation in the Bełchatów mine is conducted using the opencast method with a side and frontal block system. The open-pit mining system requires the extraction and relocation of very large volumes of soil. The removal of the overburden changes the stress and groundwater conditions in the open pit. In such a situation, ground relaxation develops, its volume increases and the shear strength decreases. These changes are particularly developed on the main structural borders.
In the geological and mining conditions of the Bełchatów mine, landslides occur on the slopes of the open pit. These threats result from the presence of numerous lithological and tectonic contact zones, layers of weak soil, the presence of residual waters, and paleo-landslide structures. About 90% of the slide surface develops in the contact zones under the conditions of the Bełchatów mine. In the history of the open-pit mining operations of the Bełchatów mine, about forty-eight large landslides were recorded, covering a volume greater than 2000 m 3 [28]. The largest number of landslides in the Bełchatów mine is associated with lithological surfaces of geological layers, such as on the clay-coal contact (51%) and, to a lesser extent, with discontinuities of tectonic origin (44%) and other causes, e.g., erosion surfaces, or with palaeo-landslides (5%).
One  Exploitation in the Bełchatów mine is conducted using the opencast method with a side and frontal block system. The open-pit mining system requires the extraction and relocation of very large volumes of soil. The removal of the overburden changes the stress and groundwater conditions in the open pit. In such a situation, ground relaxation develops, its volume increases and the shear strength decreases. These changes are particularly developed on the main structural borders.
In the geological and mining conditions of the Bełchatów mine, landslides occur on the slopes of the open pit. These threats result from the presence of numerous lithological and tectonic contact zones, layers of weak soil, the presence of residual waters, and paleolandslide structures. About 90% of the slide surface develops in the contact zones under the conditions of the Bełchatów mine. In the history of the open-pit mining operations of the Bełchatów mine, about forty-eight large landslides were recorded, covering a volume greater than 2000 m 3 [28]. The largest number of landslides in the Bełchatów mine is associated with lithological surfaces of geological layers, such as on the clay-coal contact (51%) and, to a lesser extent, with discontinuities of tectonic origin (44%) and other causes, e.g., erosion surfaces, or with palaeo-landslides (5%).
One of the largest is the 18S landslide, which gradually developed in the hanging wall of the southern fault on the southern slope of the Bełchatów open pit (Figures 2 and 3). August 1993, and then on 4 March 1996, and 30 April 1998. The volume of soil affected by the landslide process was estimated at 2 million m 3 . The landslide length was about 750 m and covered the levels of the open-pit slope at the ordinates of +205 to +78 m above sea level ( Figure 3). The maximum width of the landslide was 230 m. Characteristic geological profiles from the area of the 18S landslide are shown on the example of boreholes B-3 and B-7 ( Figure 4). Figure 4 shows an example of the slide surface that was found on the basis of geological engineering studies on the cores from the test boreholes.

Methods and Data
The scheme of the research methodology is shown in Figure 5. The methodology consists of four main stages: site investigation, sensitivity analysis, numerical analysis, and probability analysis.

Methods and Data
The scheme of the research methodology is shown in Figure 5. The methodology consists of four main stages: site investigation, sensitivity analysis, numerical analysis, and probability analysis.

Site Investigation
Based on the analysis of archival documentation and the geological and engineering mapping of landslides, it was found that the properties and structure in the research area

Site Investigation
Based on the analysis of archival documentation and the geological and engineering mapping of landslides, it was found that the properties and structure in the research area are complicated. There are different types of zones of deformed layers with the inclusion of different soils (Figure 6a), deformed contact zones (Figure 6b) that are often saturated with water (Figure 6c), or zones of strongly weathered soils. Based on geological and engineering mapping and the archival data results, nine cored boreholes with lengths of 50 m to 235 m were drilled (Figure 3a).  The soil and rock samples for laboratory tests were collected from core fragments from investigation boreholes and outcrops from the exposed geological layers. Firstly, the determination of soil parameters: grain composition, volumetric density, natural moisture, soil condition, and consistency limits was performed (Table 1). Soil shear strength tests were then conducted, using a direct shear apparatus. All tests were performed following the Polish standard (PN-B-04481:1988) [29]. For rocks, compressive strength and tensile strength tests were executed, following the Polish standard (PN-G-04303:1997) [30]. Shear strength parameters were calculated using the empirical formula for cohesion c and the angle of internal friction φ: where: σT -tensile strength, σC -compressive strength. The soil and rock samples for laboratory tests were collected from core fragments from investigation boreholes and outcrops from the exposed geological layers. Firstly, the determination of soil parameters: grain composition, volumetric density, natural moisture, soil condition, and consistency limits was performed (Table 1). Soil shear strength tests were then conducted, using a direct shear apparatus. All tests were performed following the Polish standard (PN-B-04481:1988) [29]. For rocks, compressive strength and tensile strength tests were executed, following the Polish standard (PN-G-04303:1997) [30]. Shear strength parameters were calculated using the empirical formula for cohesion c and the angle of internal friction ϕ: where: σ T -tensile strength, σ C -compressive strength.   The parameters of the lithological units are summarised in Table 1. This table presents the basic physical parameters (grain size distribution, moisture content, Atterberg limits, and bulk density) and strength parameters (cohesion and internal friction angle). For the lithological units from I to IV, relationships of shear stress-shear strain differentiated themselves (Figure 7). The greatest diversity is shown in units III and IV, which include lignite seams with carbonised clays and silts.
• LI-a clay and sand formation consisting of various-grained, brown-beige quartz sands and clay. In the top part of the formation, there are variegated clay and silt sediments with a sand interbed; • LII-clay formation with lignite lenses; • LIII-a clay, lignite, and silt formation consisting of the main layer of lignite A, black, carbonised clays, and silts. In the lower part of the formation, there are light-grey silt and sandy loams; • LIV-a formation consisting of a multi-layer lignite seam separated by clay and silt; • LV-a formation consisting of fine and medium sands with inserts of sandy silt, silted sands and silt; • LVI-Mesozoic bedrock formation consisting of limestones and marls at a depth of approx. 90 to 190 m.
The parameters of the lithological units are summarised in Table 1. This table presents the basic physical parameters (grain size distribution, moisture content, Atterberg limits, and bulk density) and strength parameters (cohesion and internal friction angle). For the lithological units from I to IV, relationships of shear stress-shear strain differentiated themselves (Figure 7). The greatest diversity is shown in units III and IV, which include lignite seams with carbonised clays and silts. The L1 unit is approximately homogeneous and shows little variation in parameter values (Figure 7a). The arithmetic mean of the peak shear stress is 150 kPa, and the residual stress is 120 kPa.
Unit LII has more variable shear strength values than unit LI. For the LII unit, the arithmetic mean of the peak is 80 kPa, and the residual value is 70 kPa.
The LIII unit differs significantly in terms of the shear strength, due to the presence of lignite. The shear surface was revealed on the contacts of layers in the samples. It was found that the shear failure developed from 3% of the shear strain for all samples from unit III. These samples show brittle behaviour. The maximum peak shear strength varies from 88 to 272 kPa. The arithmetic mean of the peak is 183 kPa, and the residual is 65 kPa. The L1 unit is approximately homogeneous and shows little variation in parameter values (Figure 7a). The arithmetic mean of the peak shear stress is 150 kPa, and the residual stress is 120 kPa.
Unit LII has more variable shear strength values than unit LI. For the LII unit, the arithmetic mean of the peak is 80 kPa, and the residual value is 70 kPa.
The LIII unit differs significantly in terms of the shear strength, due to the presence of lignite. The shear surface was revealed on the contacts of layers in the samples. It was found that the shear failure developed from 3% of the shear strain for all samples from unit III. These samples show brittle behaviour. The maximum peak shear strength varies from 88 to 272 kPa. The arithmetic mean of the peak is 183 kPa, and the residual is 65 kPa.
The LIV unit containing the lignite seams has less varied shear strength values than the LIII unit. The arithmetic mean of the peak is 130 kPa, and the residual is 75 kPa. For unit IV, shear failure developed from 6 to 9% of the shear strain for all samples.
LV and LVI units are approximately homogeneous and are characterised by low variability of the determined parameters. The peak and residual parameters were not considered since these units were not involved in the development of the landslide process.
Based on the distinguished lithological units, a cross section was made along the longitudinal axis of the landslide, which was used for numerical calculations (Figure 8).
The LIV unit containing the lignite seams has less varied shear strength values than the LIII unit. The arithmetic mean of the peak is 130 kPa, and the residual is 75 kPa. Fo unit IV, shear failure developed from 6 to 9% of the shear strain for all samples.
LV and LVI units are approximately homogeneous and are characterised by low variability of the determined parameters. The peak and residual parameters were no considered since these units were not involved in the development of the landslide pro cess.
Based on the distinguished lithological units, a cross section was made along the longitudinal axis of the landslide, which was used for numerical calculations (Figure 8). .

Sensitivity Analysis
The sensitivity analysis consisted of determining the input parameters that have the greatest impact on the stability of the open-pit slope. We have determined these param eters as being the most influential. The others' parameters are determined as the less in fluential parameters.
The measure of the impact of the input parameter on the resulting quantity is the sensitivity ratio. It is defined as the ratio of the relative change in the resulting quantity to the relative change in the input parameter (see [31,32]).
The sensitivity analysis consists of three steps, as shown in Figure 5: 1. Selection of input parameters; 2. Calculation of sensitivity ratio based on numerical calculations; 3. Determination of the most influential parameters based on the quantity of the nor malized sensitivity ratio.
The sensitivity analysis includes choosing parameters of lithological units as inpu parameters for the FoS numerical calculation, such as bulk density, cohesion, and angle of internal friction. All identified parameters concerning the stability of the open-pi slope with their estimated ranges can be considered sets of input parameters in the RSM calculation scheme ( Figure 5).
Based on the borehole data and field observations, we found that the landslide process developed within four lithological units (I to IV). This was assumed on the basi of the slide surfaces identified in test boreholes (Figure 8). For each of these units, the minimum and maximum ranges of the possible values of the parameters were deter

Sensitivity Analysis
The sensitivity analysis consisted of determining the input parameters that have the greatest impact on the stability of the open-pit slope. We have determined these parameters as being the most influential. The others' parameters are determined as the less influential parameters.
The measure of the impact of the input parameter on the resulting quantity is the sensitivity ratio. It is defined as the ratio of the relative change in the resulting quantity to the relative change in the input parameter (see [31,32]).
The sensitivity analysis consists of three steps, as shown in Figure 5: 1. Selection of input parameters; 2.
Calculation of sensitivity ratio based on numerical calculations; 3.
Determination of the most influential parameters based on the quantity of the normalized sensitivity ratio.
The sensitivity analysis includes choosing parameters of lithological units as input parameters for the FoS numerical calculation, such as bulk density, cohesion, and angle of internal friction. All identified parameters concerning the stability of the open-pit slope with their estimated ranges can be considered sets of input parameters in the RSM calculation scheme ( Figure 5).
Based on the borehole data and field observations, we found that the landslide process developed within four lithological units (I to IV). This was assumed on the basis of the slide surfaces identified in test boreholes (Figure 8). For each of these units, the minimum and maximum ranges of the possible values of the parameters were determined, specifically for bulk density, cohesion, and angle of internal friction ( Table 2). The values of the parameters of the V and VI bedrock units ( Figure 8) were taken in the form of the arithmetic mean value.  RSM requires determining the number of sources from which the parameter values are derived. These sources can be related to various field and laboratory tests or arbitrarily determined by an expert. In the study, we assumed that the parameter values were derived from two sources: "peak" stress-strain relations (source A) and "residual" stress-strain relations (source B) ( Table 2). We assumed the occurrence of a given parameter value from both sources as equally probable at 0.5 each.
When constructing the sets of input parameters, combinations of the maximum and minimum values of the parameters from different sources are made ( Table 2). In the analysed case, three strength parameters were determined for each of the four lithological units. Assuming that all parameters from units I to IV could be the most influential, the number of input sets n c for numerical calculations is 16,777,216, according to the following formula [33]: where: N-number of the most influential parameters (N = 12), n-number of sources of information about parameters (n = 2).
To limit the number of computational implementations, a normalised sensitivity ratio for a given parameter was calculated, according to Shen and Abbas (2013), using the numerical calculations in FLAC 2D v. 7.0.
Having data on the value ranges of the parameters of units I to IV and their probability of occurrence (0.5) can be presented in a random set. In the RSM, the sensitivity analysis is recommended to be performed over both small and large changes in input parameters, called local and global intervals [34]. The graphical representation of the calculated random sets for the parameters of units I to IV is shown in Figure 9. In this step, it is recommended to investigate the correlation coefficient of the analysed parameters [34].
To determine the sensitivity ratio for a given parameter, five computational implementations are required for the extreme values of the local and global range of the parameter value. If the number of all parameters is defined as N S , then the number of all L calculations is given by the formula: The calculated sensitivity ratio after normalization can be expressed as a percentage ( Figure 10). If parameters whose impact on the required result are less influential, their average value is taken for further calculations. To determine the sensitivity ratio for a given parameter, five computational implementations are required for the extreme values of the local and global range of the parameter value. If the number of all parameters is defined as NS, then the number of all L calculations is given by the formula: The calculated sensitivity ratio after normalization can be expressed as a percentage ( Figure 10). If parameters whose impact on the required result are less influential, their average value is taken for further calculations.

Numerical Stability Analysis
Stage III consists of the following steps ( Figure 5): 1. Selection of sets of input parameters;  To determine the sensitivity ratio for a given parameter, five computational implementations are required for the extreme values of the local and global range of the parameter value. If the number of all parameters is defined as NS, then the number of all L calculations is given by the formula: The calculated sensitivity ratio after normalization can be expressed as a percentage ( Figure 10). If parameters whose impact on the required result are less influential, their average value is taken for further calculations.

Numerical Stability Analysis
Stage III consists of the following steps ( Figure 5): 1. Selection of sets of input parameters;

Numerical Stability Analysis
Stage III consists of the following steps ( Figure 5):

1.
Selection of sets of input parameters; 2.
Firstly, sets of input parameters for the calculations are developed. The input data set contains a combination of the most influential parameters from various sources (Table 3). For example, with four of the most influential parameters and two data sources (A and B), the result is a combination of 16 four-element subsets. Within each subset, combinations of each parameter's minimum and maximum values are then determined, providing the next 16 subsets. Finally, 256 sets of the most influential parameter values are obtained, supplemented with average values of parameters for which the impact on the required result is less influential but required for the calculations. Table 3. Combinations of the most influential parameters from different sources (A and B).

Combinations of the Most Influential Parameters from Different Sources A and B
1 The main part of stage III is the numerical analysis of the stability of the open-pit slope. In the case under consideration, the analysis was performed using the finite difference method of FLAC 2D v.7.0. The Coulomb-Mohr strength criterion was adopted, using the associated flow rule. The vertical stress in the model increased linearly with depth. We used the SSR method to calculate the FoS. As a result of the shear-strain velocity calculations, the locations of the slide planes were obtained. On this basis, the locations of the shallow slide planes with local impact and deep-seated slide planes disturbing the greater parts of the landslide body were graphically determined.

Probability Analysis
This analysis was aimed at assigning the probability for the FoS value for the analysed sets of input parameters. The stage IV was conducted in the following steps:

1.
Development of a P-box graph of the lower and upper limits of the cumulative probability; 2.
Verification of the calculation results.
To construct a P-box graph of the upper and lower probability of the occurrence of the FoS value, the following calculations were made (e.g., Peschl, 2004):

•
FoS values for each of the source combinations were sorted in ascending order; • The FoS values within the variability range were rejected, and the minimum and maximum values were adopted for further calculations; • Each FOS value was assigned a probability of occurrence; • P-box graphs of FoS values on the horizontal axis were constructed along with the assigned probabilities on the vertical axis. As a result, the upper and lower bounds of the cumulative probability of FoS occurrence were obtained, approximated by the continuous normal distribution function.
The P-box graph is standardly used in the random set method [21,32,35]. Based on the graph, it is possible to determine the range of the most probable values of the FoS for the analysed set of input parameters.

Sensitivity Analysis
Assuming a 10% threshold [32,35], four of the most influential parameters were distinguished: the cohesion of the lithological unit I cLI, the cohesion of unit II cLII, the cohesion of unit III cLIII, and the internal friction angle of unit III ϕLIII (Figure 10). The cohesion and internal friction angle were assumed to be uncorrelated. The impact of the other six parameters on the required result was less influential. In the four most influential parameters and six less influential parameters, the number of combinations was 16 (Table 3). In accordance with formula 5, we obtained 256 sets of input parameters for the FoS calculations.

Numerical Analysis
The numerical model was developed on the basis of the I-I' section along the longitudinal axis of the landslide (Figure 11). This model has six lithological units. The length of the model was 750 m. The left frame was 283 m high, while the right frame was 166 m. The boundary conditions were introduced in such a way as to block the possibility of horizontal displacements on the right and left frames and vertical displacements on the bottom one. The continuous line reflecting the terrain's morphology was a free surface in which displacement was allowed in any direction. The position of the groundwater table was introduced at various depths up to 20 m. The model was digitised with a regular 1 × 1 m grid. The calculation model is illustrated in Figure 11.
The P-box graph is standardly used in the random set method [21,32,35]. Based on the graph, it is possible to determine the range of the most probable values of the FoS for the analysed set of input parameters.

Sensitivity Analysis
Assuming a 10% threshold [32,35], four of the most influential parameters were distinguished: the cohesion of the lithological unit I cLI, the cohesion of unit II cLII, the cohesion of unit III cLIII, and the internal friction angle of unit III φLIII ( Figure 10).
The cohesion and internal friction angle were assumed to be uncorrelated. The impact of the other six parameters on the required result was less influential. In the four most influential parameters and six less influential parameters, the number of combinations was 16 (Table 3). In accordance with formula 5, we obtained 256 sets of input parameters for the FoS calculations.

Numerical Analysis
The numerical model was developed on the basis of the I-I' section along the longitudinal axis of the landslide (Figure 11). This model has six lithological units. The length of the model was 750 m. The left frame was 283 m high, while the right frame was 166 m. The boundary conditions were introduced in such a way as to block the possibility of horizontal displacements on the right and left frames and vertical displacements on the bottom one. The continuous line reflecting the terrain's morphology was a free surface in which displacement was allowed in any direction. The position of the groundwater table was introduced at various depths up to 20 m. The model was digitised with a regular 1 × 1 m grid. The calculation model is illustrated in Figure 11. We performed 256 variants of numerical calculations for sets of input parameters determined as a result of the sensitivity analysis. In each case, the FoS values were calculated, which varied from 0.28 to 2.00. For further analysis, based on experiments, a FoS We performed 256 variants of numerical calculations for sets of input parameters determined as a result of the sensitivity analysis. In each case, the FoS values were calculated, which varied from 0.28 to 2.00. For further analysis, based on experiments, a FoS of less than 1.1 was assumed, which meant the loss of stability or the state of critical equilibrium of the analysed landslide on the open-pit slope of the Bełchatów mine. In total, 166 calculation variants were within the range of FoS changes, from 0.28 to 1.1. Many of these calculation variants had similar slide planes positions. Figure 12 shows the generalised slide planes grouped, according to the similarity of their course. The shallowest slide planes with a local range are marked with a purple dotted line. The remaining deepseated slide planes are marked with a blue line. The red line marks the location of the slide surface determined from borehole data and field observations.
It should be noted that most of the slide planes developed in the lower part of the landslide on the border of the II and III units and inside the III unit close to the border with the V unit. Presumably, this is the effect of a thinning of the III unit. Both units II and III contain lignite seams separated by clay. Thus, the course of the slide plane is very diverse in these lithological units. In particular, unit III shows strong heterogeneity, expressed by a large variability of parameter values.
The calculated slide planes located closest to the slide surface identified in the boreholes have an FoS value between 0.69 and 0.70. There are five such slide planes denoted by numbers of calculation variants: 105, 161, 169, 225, and 233. They are generalised in Figure 12 in the form of a single slide plane denoted with a green dotted line. This generalised slide plane develops in unit I, at the bottom of unit II, and at the top of unit III. In the upper part of the landslide body, in unit I, there is quite a large inaccuracy in the position of the generalised calculated slide plane, compared to the observed slide surface in the boreholes. This inaccuracy amounts to a maximum of around 30 m. However, in the lower part of unit II and on the border with unit III, both surfaces are consistent with an inaccuracy of several dozen centimetres. This results from a significant diversification of the strength parameter between the values of both units II and III and in the weathered unit I ( Table 2).  Figure 12 shows the generalised slide planes grouped, according to the similarity of their course. The shallowest slide planes with a local range are marked with a purple dotted line. The remaining deep-seated slide planes are marked with a blue line. The red line marks the location of the slide surface determined from borehole data and field observations. It should be noted that most of the slide planes developed in the lower part of the landslide on the border of the II and III units and inside the III unit close to the border with the V unit. Presumably, this is the effect of a thinning of the III unit. Both units II and III contain lignite seams separated by clay. Thus, the course of the slide plane is very diverse in these lithological units. In particular, unit III shows strong heterogeneity, expressed by a large variability of parameter values.
The calculated slide planes located closest to the slide surface identified in the boreholes have an FoS value between 0.69 and 0.70. There are five such slide planes denoted by numbers of calculation variants: 105, 161, 169, 225, and 233. They are generalised in Figure 12 in the form of a single slide plane denoted with a green dotted line. This generalised slide plane develops in unit I, at the bottom of unit II, and at the top of unit III. In the upper part of the landslide body, in unit I, there is quite a large inaccuracy in the position of the generalised calculated slide plane, compared to the observed slide surface in the boreholes. This inaccuracy amounts to a maximum of around 30 m. However, in the lower part of unit II and on the border with unit III, both surfaces are consistent with an inaccuracy of several dozen centimetres. This results from a significant diversification of the strength parameter between the values of both units II and III and in the weathered unit I ( Table 2). Figure 13 shows the graphs of the lower and upper bounds of the cumulative probability of the occurrence of FoS values. Within the lower bound of the discrete cumulative probability for FoS values smaller than 1.1, FoS = 0.28 had the highest probability (25%). However, the slide plane course was not entirely consistent with the borehole  Figure 13 shows the graphs of the lower and upper bounds of the cumulative probability of the occurrence of FoS values. Within the lower bound of the discrete cumulative probability for FoS values smaller than 1.1, FoS = 0.28 had the highest probability (25%). However, the slide plane course was not entirely consistent with the borehole data and field observations ( Figure 12). The next highest probability of FoS (12.5%) occurred in the range from 0.69 to 0.70. This group of slide planes is located closest to the slide surface identified in borehole data and field observations. Energies 2021, 14, x FOR PEER REVIEW 17 of 20 data and field observations ( Figure 12). The next highest probability of FoS (12.5%) occurred in the range from 0.69 to 0.70. This group of slide planes is located closest to the slide surface identified in borehole data and field observations. The probability of the occurrence for all calculated slide planes were approximated using a normal distribution, separately for lower and upper bounds. The matching parameters are presented in Table 4. The presented data show that the probability of the occurrence of FoS values smaller than 1.1 is from 10 to 97%.  The probability of the occurrence for all calculated slide planes were approximated using a normal distribution, separately for lower and upper bounds. The matching parameters are presented in Table 4. The presented data show that the probability of the occurrence of FoS values smaller than 1.1 is from 10 to 97%.

Setting the Critical Strength Parameters
Each of the slide planes for FoS = 0.69 − 0.70 has an assigned set of input parameters. Table 5 shows the most influential input parameters for these slide planes assigned to the generalised slide plane in Figure 12 (denoted with green dotted line). The strength parameters of units II and III, shown in Table 5, particularly the residual cohesion values, had the greatest impact on landslide formation. In Table 6, all critical parameters for each lithological unit required to perform numerical stability analysis are summarised, including those that are the most and least influential. Due to the variability of the strength parameters listed in Table 5, they were averaged to determine the median. The numerical stability analysis performed for the set of parameters in Table 6 confirmed the consistent position of the resultant slide plane with the generalised slide plane shown in Figure 12. The new FoS calculated for this resulted slide plane is 0.70. In the study, the method of setting the critical parameters for the landslide of the Bełchatów open-pit mine is certainly a more labour-intensive exercise than the methods presented in the works mentioned in the introduction section (e.g., [6,15,16]). However, our method enables the consideration of a great variability of parameter values, due to the severe impact of mining works in an open-pit slope. Thus, the execution of reliable stability analysis is a challenging task. Based on the proposed methodology, we reconstructed the failure plane that caused the landslide in the Bełchatów mine. The great advantage of the SSR method is the possibility of determining many sliding planes. In conjunction with the RSM, it can identify the very probable slide planes and set of critical strength parameters. In the study, we chose the slide plane located closest to the slide surface identified in the boreholes and field observations.
Knowledge of the critical values of the strength parameters enables the design of a protection system for such a landslide and the taking of preventive measures to reduce this risk, such as decreasing the inclination of the open-pit slope, changing the geometry of open-pit benches, or designing a hydrogeological barrier.

Conclusions
In the study, we developed an approach for setting strength parameters for numerical stability analysis of landslides in an open-pit slope, using the random set method, on the example of the Bełchatów lignite mine. Based on the study results, the following conclusions can be drawn: 1.
The methodology of calculations allowed determination of the slide surface with a factor of safety, most closely to the slide surface identified in boreholes and field observations. In the considered model of the landslide, the slide plane of the landslide developed in lithological units of a complex structure and properties with strongly variable values of strength parameters.

2.
Numerical calculations using the random set method allowed determination of the set of the most influential input parameters for assessing the stability of the open-pit slope model. The strength parameters of units II and III, particularly the residual cohesion values, had the greatest impact on landslide formation. 3.
The stability analysis, using the shear strength reduction method and the random set method provided information about different courses of the slide surfaces with the assigned probability of the factor of safety. The proposed methodology can be used in stability analyses of other types of slopes with strongly different properties, such as activated colluviums and slopes of embankments of communication routes. It is useful in conditions where the determination of the average input parameters is subject to high uncertainty, and residual parameters give unrealistic solutions. A significant limitation of this method is the fact that it is highly labour intensive. It is related to many calculations, increasing with the complexity of the structure of the calculation model and with the number of information sources. In more complex calculations, it is possible to perform partial calculations sequentially in the FLAC software, making the process more user-friendly.