Probability Methods for Stability Design of Open Pit Rock Slopes: An Overview

: The rock slope stability analysis can be performed using deterministic and probabilistic approaches. The deterministic analysis based on the safety concept factor uses ﬁxed representative values for each input parameter involved without considering the variability and uncertainty of the rock mass properties. Probabilistic analysis with the calculation of probability of failure instead of the factor of safety against failure is emerging in practice. Such analyses offer a more rational approach to quantify risk by incorporating uncertainty in the input variables and evaluating the probability of the failure of a system. In rock slope engineering, uncertainty and variability involve a large scatter of geo-structural data and varied geomechanical test results. There has been extensive reliability analysis of rock slope stability in the literature, and different methods of reliability are being employed for assessment of the probability of failure and the reliability of a slope. Probabilistic approaches include Monte Carlo simulation (MCS), the point estimate method (PEM), the response surface method (RSM), ﬁrst- and second-order reliability methods (FORMs and SORMs), and the ﬁrst-order second-moment method (FOSM). Although these methods may be complicated, they provide a more complete deﬁnition of risk. Probabilistic slope stability analysis is an option in most commercial software; however, the use of this method is not common in practice. This paper provides an overview of the literature on some of the main probabilistic reliability-based methods available for the design of the rock slope in open pit mining. To demonstrate its applicability, the paper investigates the stability of a rock slope in an open pit mine in the Goldﬁelds region, Western Australia. Two different approaches were adopted: deterministic stability analysis using two-dimensional limit equilibrium and ﬁnite element shear strength reduction methods using SLIDE and RS2 software, respectively, and probabilistic analysis by applying the MCS and RSM methods in the limit equilibrium method. In this example, the slope stability analysis was performed using the Spencer method with Cuckoo search optimization to locate the critical slip surface. The results obtained were compared and commented on. procedures. A brief description of


Introduction
Probability analysis of rock slopes has gained considerable attention in the design of open pit mines (e.g., [1][2][3][4][5][6]). The unavoidable uncertainties involved in geotechnical design parameters has, however, attracted significant research in the use of reliability analysis of slope stability over the past few decades [7][8][9][10][11][12]. In most cases, the stability of the open pit slope is expressed in terms of the factor of safety. The factor of safety is determined by the deterministic methods, which include limit equilibrium (LEM) and the shear strength reduction method (SSRM). For complex cases of variable slope geometry and geological settings, deterministic methods cannot be applied. The natural variability Table 1. Relationship between the reliability index (β) and probability of failure (P f ) [22].

Reliability Index β
Probability of Failure P f =Φ(−β) Expected Performance Level Both LEM and SSRM, although originally deterministic, can easily be adapted to suit probabilistic models [23]. The probability of failure and the reliability index are used to quantify risks and hence evaluate the consequence of failure [24]. Various probability methods have been proposed to estimate P f and/or β (e.g., [14,[25][26][27][28]). The most widely used methods are the first-order reliability method (FORM) (e.g., [29][30][31]), first-order second-moment (FOSM) method (e.g., [7,26]), and second-order reliability method (SORM) (e.g., [32,33]). The reliability-based design approach has been developed and utilized to address the shortcomings of the deterministic approach and considers uncertainties explicitly [14,34]. Several applications of solution methods using Monte Carlo simulation (e.g., [35][36][37]), the point estimate method (e.g., [19,38,39]), and the response surface method (RSM) have been reported in the literature (e.g., [32,40,41]). [7] explored the mean firstorder reliability method, which is a simplification of the more general first-order reliability method. [21] used FORM to model rock plane failure. [42] used Monte Carlo simulation as sensitivity tool for calculating the probability of failure [43]. The reliability methods and the response surface methods have been used for slope reliability problems with an implicit performance function (e.g., [40]).
The early applications of reliability methods (e.g., [44][45][46][47][48]) were limited to theoretical case studies. Currently, the probability and reliability methods for slope stability constitute an already well-developed research area and most of the methods are currently implemented in numerical code for the design of rock slopes in practice. The most important literature in this area includes journal articles and books (e.g., [14,[49][50][51]). Although the system reliability methods have been used extensively in theoretical cases and in civil engineering applications, these approximation methods may not be available for practical open pit slopes with complicated geometry and/or multiple slope bench due to the perceived mathematical complexity. More importantly, with the advent of computers in data handling, Monte Carlo simulation, Latin hypercube, and the response surface method (RSM) have been adopted in LEM software programs, such as SLIDE. The SSRM methods, however, have been adopted in several well-known finite element (RS2/3) or finite difference (FLAC) programs. In addition, the FOSM, for instance, is time consuming, which is quite common for the implicit performance for LEM analysis [52]; however, researchers have tried to overcome these problems using methods like RSM [53]. Recently, [53] used Gaussian process regression (GPR) and genetic programming (GP) to resolve the problems of the implicit function performance function. This GPR is a nonparametric kernel-based probabilistic model.
However, while LEM is the most widely used method for evaluating slope stability in practice, extensive research has been undertaken to improve its performance, particularly in finding the global critical slip surface [23]. Again, while the determination of the critical slip surface may be affected by the experience of the researcher or engineer, several search algorithms, such as the cuckoo, particle swam optimization, simulated annealing, etc., have been proposed to optimize the search for the critical slip surface. On the other hand, to carry out a probabilistic analysis using SRM, a spatially correlated field is typically developed in random field theory and then solved using finite element or finite difference methods [54].
This paper aims to review some probability methods that have been used in the design of rock slope reliability problems. A brief discussion of the uncertainty and variability associated with rock mass and some methods of quantification is presented. Then, the cuckoo and particle swam optimization search methods are briefly discussed. Based on existing knowledge of the uncertainty associated with rock slopes, and the limitation of the implicit reliability functions, the paper adopts the cuckoo techniques coded in the Rocscience SLIDE computer software to analyze the stability of a case pit slope located in the Goldfields, Western Australia. Deterministic and probabilistic analyses were carried out using LEM using the Rocscience SLIDE computer software to determine the factor of safety, probability of failure, and reliability index. The reliability of the case pit slope stability was evaluated using Monte Carlo simulation and the response surface method and the results were compared.

Uncertainty in Rock Slope Engineering and Methods of Quantification
In general, sources of uncertainty in slope engineering include the natural variability of the rock material, limited availability of information about the ground condition, and errors made during measurements and testing of samples [55]. There are various classifications of uncertainty in the literature; however, one classification that proves adequate for geotechnical engineers places uncertainty into three groups ( Figure 1): geological uncertainty, model uncertainty, and parameter uncertainty [14,56,57]. Geological uncertainty comprises uncertainties associated with the geometry of geological structures and their relationships between lithologies and those associated with the boundaries of lithologies [56,58]. Model uncertainty exists if there is a possibility of obtaining an incorrect result even if the exact values are available for all the model parameters [56]. It reflects the inability of a model or design technique to represent the true physical behavior of a system under consideration [14,58,59]. Parameter uncertainty reflects the inability to account for the various attributes of the geotechnical model [56]. It includes uncertainties associated with the values adopted for the rock mass and hydrogeological model parameters, which can stem from data scatter, such as the spatial variability and random testing error and systematic error, which involves statistical error and bias in measurements. under consideration [14,58,59]. Parameter uncertainty reflects the inability to account for the various attributes of the geotechnical model [56]. It includes uncertainties associated with the values adopted for the rock mass and hydrogeological model parameters, which can stem from data scatter, such as the spatial variability and random testing error and systematic error, which involves statistical error and bias in measurements. As uncertainty is present in the in situ material, it means uncertainty will also exist in the expected performance function of the design [60]. If a value of the performance function like cohesion or the friction angle varies, and this variation over space and time cannot be predicted, one cannot say with certainty what the value of the variable will be but only the likelihood or probability that it will be within some specific range of values. Unfortunately, the overall uncertainty is rarely quantified in rock engineering. Instead, conservative designs are generally adopted through deterministic analysis. However, as single values are assigned in the deterministic analyses, there is no guarantee that the design will perform as expected. Therefore, there is the need to have an approach that can consider varying rock variabilities and uncertainties in rock slope conditions during the design of rock slopes and analysis of rock slope stability.
Various mathematical frameworks have been developed for the assessment of uncertainty and variability in slope stability analysis (e.g., [7,8,11,12,14,19,25,32,34,35,37,55,[61][62][63][64][65]). Two main methodologies have been proposed to deal with uncertainties in the rock properties in the assessment of slope stability, i.e., the reliability method and non-deterministic methods [66]. As shown in Figure 2, the non-deterministic methods consist of probabilistic methods and non-probabilistic methods (also called imprecise methods). In the probabilistic methods, the rock properties affecting the stability of the slope are considered as random variables that have a certain probability distribution. Some of the common non-probabilistic methods are evidence theory, fuzzy set theory, the interval ap- As uncertainty is present in the in situ material, it means uncertainty will also exist in the expected performance function of the design [60]. If a value of the performance function like cohesion or the friction angle varies, and this variation over space and time cannot be predicted, one cannot say with certainty what the value of the variable will be but only the likelihood or probability that it will be within some specific range of values. Unfortunately, the overall uncertainty is rarely quantified in rock engineering. Instead, conservative designs are generally adopted through deterministic analysis. However, as single values are assigned in the deterministic analyses, there is no guarantee that the design will perform as expected. Therefore, there is the need to have an approach that can consider varying rock variabilities and uncertainties in rock slope conditions during the design of rock slopes and analysis of rock slope stability.
Various mathematical frameworks have been developed for the assessment of uncertainty and variability in slope stability analysis (e.g., [7,8,11,12,14,19,25,32,34,35,37,55,[61][62][63][64][65]). Two main methodologies have been proposed to deal with uncertainties in the rock properties in the assessment of slope stability, i.e., the reliability method and non-deterministic methods [66]. As shown in Figure 2, the non-deterministic methods consist of probabilistic methods and non-probabilistic methods (also called imprecise methods). In the probabilistic methods, the rock properties affecting the stability of the slope are considered as random variables that have a certain probability distribution. Some of the common non-probabilistic methods are evidence theory, fuzzy set theory, the interval approach, possibility theory, random set theory, etc. [58,66]. For instance, uncertainty characterized by fuzziness is treated with a branch of methodologies based on a fuzzy representation of uncertain variables. A complete description of each method is outside the scope of this paper and the reader is referred to the documents cited for more information on the mathematical formulations and procedures. A brief description of FORM, FOSM, SORM, PEM, MCS, and RSM is given in the next section to provide some insight on the mathematical meaning. Reliability analyses provide a more rational approach to quantity slope design risk than a deterministic method by incorporating uncertainty in the input parameters in the analysis. By doing so, the probability of failure can be established for a specific failure mode.
proach, possibility theory, random set theory, etc. [58,66]. For instance, uncertainty characterized by fuzziness is treated with a branch of methodologies based on a fuzzy representation of uncertain variables. A complete description of each method is outside the scope of this paper and the reader is referred to the documents cited for more information on the mathematical formulations and procedures. A brief description of FORM, FOSM, SORM, PEM, MCS, and RSM is given in the next section to provide some insight on the mathematical meaning. Reliability analyses provide a more rational approach to quantity slope design risk than a deterministic method by incorporating uncertainty in the input parameters in the analysis. By doing so, the probability of failure can be established for a specific failure mode.

Reliability Index and Probability of Failure
The main task of reliability analysis is to calculate the probability of failure and reliability index. To perform a reliability analysis, a performance function g(X) must be defined that relates to the resistance R(X) and disturbance S(X) acting on the system. This is mathematically expressed as: where X represents the collection of random input variables. By this definition, failure occurs if g(X) < 0, while g(X) > 0 denotes stable conditions ( Figure 3). The surface defined by g(X) = 0 is referred to as the critical limit state as it defines the boundary between these two conditions [60]. When considering the critical limit state for a rock slope under planar failure, the performance function can be expressed as the shear strength that resists sliding minus the shear forces that initiate sliding. Such equations can be evaluated analytically with little additional effort. For more complex problems, such as an analysis of rock slope

Reliability Index and Probability of Failure
The main task of reliability analysis is to calculate the probability of failure and reliability index. To perform a reliability analysis, a performance function g(X) must be defined that relates to the resistance R(X) and disturbance S(X) acting on the system. This is mathematically expressed as: where X represents the collection of random input variables. By this definition, failure occurs if g(X) < 0, while g(X) > 0 denotes stable conditions ( Figure 3). The surface defined by g(X) = 0 is referred to as the critical limit state as it defines the boundary between these two conditions [60]. When considering the critical limit state for a rock slope under planar failure, the performance function can be expressed as the shear strength that resists sliding minus the shear forces that initiate sliding. Such equations can be evaluated analytically with little additional effort. For more complex problems, such as an analysis of rock slope deformation under seismic conditions, it is difficult to define the loads and resistances explicitly. Approximate methods of evaluation are therefore required. deformation under seismic conditions, it is difficult to define the loads and resistances explicitly. Approximate methods of evaluation are therefore required. Figure 3. Illustration of the linear and nonlinear limit state function [58].
To calculate the reliability of the system, the distance between the mean value of the performance function and the critical limit state at g(X) = 0 must be determined. When the distance between these two points is normalized with respect to the standard deviation of the performance function, this is referred to as the reliability index β for the system ( Figure  3). The reliability index is defined as [60]: where μg and σg are the mean and standard deviation of the performance function, respectively. However, to solve for the value of β in Equation (2), the exact shape of the performance function must be known, which is not always the case. Based on this, a more versatile measurement of reliability is the Hasofer-Lind reliability index βHL. This method, also known as the first-order reliability method (FORM), calculates the minimum distance in units of the directional standard deviation from the mean value point of the multivariate distribution of the random variables to the boundary of the critical limit state ( Figure  3). This provides a more consistent and invariant measure of reliability for the system and can also be easily calculated for correlated or uncorrelated variables using the approach outlined in [29]. The matrix equation for βHL can be defined as [60]: or: where X is the vector of random variables of the set random variable, μ is the vector of mean values of random variables, σ is the standard deviation, and F defines the failure region of g(X) < 0. The variable C defines the covariance matrix and R is the correlation matrix, which allows the user to establish either a positive or negative relationship between random variables. Equation (4) is preferred to Equation (3) because the correlation To calculate the reliability of the system, the distance between the mean value of the performance function and the critical limit state at g(X) = 0 must be determined. When the distance between these two points is normalized with respect to the standard deviation of the performance function, this is referred to as the reliability index β for the system ( Figure 3). The reliability index is defined as [60]: where µ g and σ g are the mean and standard deviation of the performance function, respectively. However, to solve for the value of β in Equation (2), the exact shape of the performance function must be known, which is not always the case. Based on this, a more versatile measurement of reliability is the Hasofer-Lind reliability index β HL . This method, also known as the first-order reliability method (FORM), calculates the minimum distance in units of the directional standard deviation from the mean value point of the multivariate distribution of the random variables to the boundary of the critical limit state ( Figure 3). This provides a more consistent and invariant measure of reliability for the system and can also be easily calculated for correlated or uncorrelated variables using the approach outlined in [29]. The matrix equation for β HL can be defined as [60]: or: where X is the vector of random variables of the set random variable, µ is the vector of mean values of random variables, σ is the standard deviation, and F defines the failure region of g(X) < 0. The variable C defines the covariance matrix and R is the correlation matrix, which allows the user to establish either a positive or negative relationship between random variables. Equation (4) is preferred to Equation (3) because the correlation matric R is easier to set up and conveys the correlation structure more explicitly than the covariance matrix. For uncorrelated variables, the matrix R simplifies to a symmetric unit matrix. The matrix algebra may be unfamiliar to some; however, programs, such as Microsoft Excel or MATLAB, can be used to easily complete these calculations (e.g., [29,31,67]).
After the reliability index has been determined, the probability of failure P f for the system can be found by calculating the probability of g(X) < 0 ( Figure 3). This is related to the reliability index using the following equation: where Φ is the cumulative distribution function (CDF) for the performance function evaluated at 0 with a unit standard deviation and a mean β. As mentioned earlier, the shape of the distribution is rarely known and therefore must be assumed. In most cases, a normal distribution is reasonable; however, a truncated or lognormal distribution may be more appropriate when the performance function depends on positive functions, such as the factor of safety, extent of yield in a material, or displacements [60]. Again, from Equation (2), a constant mean value suggests that the reliability index increases, and the uncertainty in the estimate of the performance function decreases. This results in a narrower distribution for the performance function and a decrease in the probability of failure for the system. FORM is widely used because of its efficiency. However significant errors may arise when the nonlinearity of the failure/performance function increases. This nonlinearity is due to the nonlinear relationship between random variables, the consideration of non-normal random variables, and/or the transformation from a correlated to uncorrelated random variable [24]. As many of the reliability methods are based on approximations, it is likely that different methods will produce different results. Different methods must be compared to obtain a more accurate understanding of a system. There are two variations of FORM, i.e., FOSM and AFOSM. In FOSM, the statistical distribution of the random variables is ignored whereas in AFOSM, the distribution of the random variables is considered.

Reliability Methods
As mentioned earlier, for a more a complex problem, the performance function cannot be stated explicitly. The reliability methods are normally coupled with the finite element method to evaluate the performance function. The number of evaluations and what input parameters are selected depend on the reliability methods used [60]. The following section briefly describes five methods that can be used to approximate the statistical moments of the performance function.

First Order Second Moment
FOSM has long been applied to assess the reliability of slopes [68]. It consists of the Taylor series of approximations of the mean and variance of the performance function. Where the performance function is smooth and regular, the mean and variance can be calculated using the first terms of the Taylor series expansion method to expand g(X) [69,70]. It requires a linearized form of the performance function at the mean values of the random variables [24,70]. This method assumes that the expected value of the performance function is approximately equal to the value of the function calculated with the mean values of all variables [60,70]. The variance is determined by calculating the partial derivatives of the performance function with respect to each uncertain variable. For uncorrelated input random variables, the variance of the function is given as: where X i denotes the random variables and n is the number of random variables. As stated earlier, since the performance function cannot be stated explicitly in most geotechnical engineering applications, a linear approximation for the partial derivatives is required [70]. To achieve this, each of the variables is changed by a small (∆X i ) amount while all other variables are kept at their mean values. The change in the performance (∆G) that results is then divided by the difference in the input. Therefore, to maintain a consistent level of uncertainty, the input variables are chosen at the mean plus and minus one standard deviation, and Equation (6) can be revised as [60]: Knowing the moments of the performance function, the reliability index and probability of failure can be calculated using Equations (2) and (5) assuming a normal distribution [69,70].

Second-Order Reliability Method
SORM was initially studied by [71,72]. The method was developed to enhance the accuracy of the estimated probability of failure. Since the performance function is approximated by a linear function, the accuracy of FORM deteriorates when the nonlinearity of a limit state function increases. In other words, SORM overcomes the disadvantage of FORM. SORM is more computationally expensive than FORM since a second derivative is required. Once the second-order surface is obtained, according to the asymptotic formula of Breitung, the probability of failure can be calculated. Breitung's formulation [72] for SORM is given by [73]: where v i (i = 1, . . . , n − 1) are the principal curvatures of the limit state function at the most probable point (MPP) and β is the reliability index obtained from FORM. MPP is the point that has the highest probability density on the performance g(X) = 0 as shown in Figure 4. The principal curvatures of the limit state surface are obtained as the eigenvalues of the rotational transformed second-order derivatives matrix called the Hessian matrix of the performance function in the standard normal space [74]. In attempt to address this problem, a number of studies have been performed aiming to eliminate the calculation of the Hessian matrix (e.g., [73,[75][76][77]). The other popular formulation is given by [78], which is considered more accurate than Breitungs's formulation [73,76,77]. Ref [73] modified SORM, called the second-order reliability method with first-order efficiency (SORM-FOE). According to [73], if the derivatives of the MPP search are evaluated numerically, the number of functions required for FORM will be linearly proportional to the number of random variables n as [73,74]: where k represents the number of iterations of the MPP search. However, if the finite difference formula is used for the derivative evaluation, the number of functions required by SORM is: Clearly, the SORM is second-order efficient because N SORM is quadratic in terms of n. SORM-FOE improves the accuracy of FORM while maintaining a similar level of efficiency [73].
1 Figure 4. Illustration of the probability density function (PDF) showing the reliability index β and probability of failure pf for a performance function of a system.

Point Estimate Method
The point estimate method (PEM) originally proposed by [79,80] is one of the most popular numerical procedures that approximates the expected value and the variance of a performance function by evaluating it as a series of specifically chosen discrete points. Evaluation of the points is chosen at the mean plus and mean minus one standard deviation for each variable, resulting in 2 n evaluations for n random variables. A weighting value ρ is used at each evaluation point to ensure the expected (mean) value and standard deviation of the input parameters are recovered [19,39,59,60,81,82]. The method accounts for a maximum of three statistical moments, namely mean, variance or standard deviation, and skewness. It does not require prior knowledge of the shape of any probability density function of the input variables and the spatial correlation; however, this approximate method may lead to incorrect interpretations of the reliability if the var function g(X) is highly nonlinear or the random variables are asymmetric. According to [60], if all the evaluation points are weighted equally, this value is 1/n for each variable. The statistical mean and variance are given by the following equations: where x i is the input variable, g(x i ) is the function of the input variable x i , and n is the number of variables. When the coefficients of variation for the input parameters are small, PEM is found be robust; however, the number of evaluations can be significantly high when a large number of variables are considered. [83], including several other authors, have developed methods to reduce the number of evaluations, and users need to be mindful of the assumptions [60,81,82].

Monte Carlo Simulation
The Monte Carlo (MC) simulation method is considered a very powerful tool. It was developed in 1949 by John von Neumann and Stanislav Ulam when they published a paper "the Monte Carlo method". Due to its robustness and concept simplicity, it has been widely used in reliability analyses [14,34,37,84,85]; it provides efficiency for engineers with basic working knowledge of probability and statistics for risk evaluation of risk or reliability of complex engineering systems. Thus, when the behavior of the performance function is difficult to evaluate, the probability of failure can be calculated directly by using Monte Carlo simulation. In this method, large sets of randomly selected input variables are generated according to their probability distribution function in the analytical model to determine the behavior of the system [60,82]. The method generates a random number for each of the input random variables in the problem and it makes combinations amongst all these random variables to perform several deterministic computations [5,86]. The accuracy of the method depends on the number of simulations performed and increases with an increasing number of simulations. Though the method has some advantages, it can be computationally intensive and time consuming. However, the lack of approximations makes Monte Carlo an ideal standard to compare to other reliability methods.

Response Surface Method
The response surface method (RSM) has been presented as an efficient tool to identify the likelihood of the failure behavior of rock slopes [87]. The method has been applied in using the central composition method design (CDD), which is one of the statistical design methods used to implement experiments to examine the effects of the main interactions of different levels of independent variables on the resulting response (dependent variable). As a result, an equation for the response, i.e., factor of safety is established as a function of the design variable (independent variables) of the response surface. The resulting mathematical model or response equation is used to estimate the probability of an unsatisfactory performance in the rock slope [88]. The RSM can be used to approximate the performance function by relating the input and output parameters for a system by a simple mathematical expression. It uses a small number of strategically selected computations to create a response surface of factor of safety (FS) values for various combinations of input parameters. It then predicts the factor of safety values for any combination of samples and provides an estimated probability of failure. It has been shown that for a potential slip surface of a slope, the relationship between the factor of safety and the input parameters can be approximated by a quadratic polynomial function [41,89,90]. Thus, in reliability analyses, the exact limit state function g(X) can be approximated by the polynomial function g'(X) [60,81]: where X = (x 1 , ..., x i, ..., x n ) is the vector of the input random variables, n is the number of input random variables or number of random field elements; and a = (a 0 , b 1 , . . . , b n , c 1 , . . . ., c n ) T is the vector of unknown coefficients that must be determined [91]. According to [60], to properly evaluate the number of unknowns in the quadratic equation, 2n + 1 evaluations are required [91][92][93]. A regression-based approach is used to compute the unknown coefficients (e.g., [91]).
From [94], the method works by: (1) converting all random variables to standard normal random variables (0, 1); (2) representing the resulting factor of safety in polynomial chaos expansion form; and (3) using a small number of computations to determine the coefficients of the polynomial in step 2 and (4) generate Latin hypercube samples and plugging them into the polynomial to estimate the factor of safety. According to [92], the initial random variables are converted to standard normal random variables using transformation equations [92,94]. Hence, the Hermite polynomial chaos expansion of the factor of safety for a given failure surface looks like this: where F is a random output of the model; a i 1 i 2,... i n are deterministic coefficients in the expansion to be estimated; n is the number of variables used to represent the uncertainty in the model inputs; U = U i 1 , U i 2 , . . . , U i n is a vector of independent standard normal variables; and Γ n U i 1 , U i 2 , . . . , U i n is the polynomial chaos of order n [92]. A complete description of the above mathematical formulation is outside the scope of this paper and readers are referred to articles and monograph (e.g., [92,95]).
Once the approximate limit state function has been established, FORM Equation (3) is used to determine the reliability index directly. This is more accurate than the FOSM method as it uses geometric interpretations to determine the reliability index rather than determining statistical moments through a linear extrapolation of the mean input values. The advantage of the combined RSM/FORM method is that it can be used for correlated and non-normal input variables and is suitable for any linear limit state surface. One disadvantage is the assumption that the inputs and outputs are related through a quadratic equation, which may not be valid in all situations.

Cuckoo Search Optimization Method
The cuckoo search (CS) is a metaheuristic optimization algorithm that is inspired by [96], based on cuckoos' breeding behavior. This method has been used in welded beam and spring designs, data fusion in wireless network sensors, and recently in slope stability designs. [97,98] used this method in 2-D slope stability analysis [99]. The application of the cuckoo search method in 3-D slope stability analysis has also been reported (e.g., [99,100]). From [99,101], the following rules applies to the CS algorithm [96]: a.
Each cuckoo lays one egg at a time and dumps it in a randomly chosen nest. b.
The best nests with high-quality eggs will carry over to the next generations. c.
The number of available host nests is fixed, and hosts can discover an alien egg with a probability pa ∈ (0, 1). Based on these rules, the host bird can either throw the egg away or abandon the nest and build a new nest [101]. By the last rule, the fraction p a can be used to determine the worst solutions of n nest that will be replaced with a new nest randomly [99,101]. To solve the problem, the illustration is that every egg in a nest represents one new solution. The aim is to use the new and better solution to replace the current solution in the nest. In some cases, the nest may have two eggs (solution), but the problem can be simplified so one nest has only one solution [99,102]. From [97], the CS algorithm begins by initializing a fixed number of n valid solution vectors {P 0 , . . . ,P i , . . . ,P N−1 |F(P i ) > 0 exist ∀i ∈ [0, N − 1]}. A number of iterations I max is defined; both N and I max depend on the dimensionality of the problem. The solution vectors are sorted from worst fitness (i.e., highest factor of safety) to the best fitness [94].
Several other authors have developed methods to refine the solution and readers are referred to the documents cited. The method has been compared with other algorithms, such as PSO and the genetic algorithm, and the results show that CS has a higher success rate [96,99].

Particle Swarm Optimization Method
The particle swarm optimization (PSO) was initialized by Kennedy and Eberhart in 1995 [103,104]. The method simulates bird flock activities when they randomly search for food in their path. Each solution is considered a particle in the search space and each particle has a fitness value. During movement, each particle adjusts its position by changing its velocity according to its own experience and the group's experience, finally moving to the optimal position [105][106][107]. PSO has gained popularity in the field of structural engineering (e.g., [108][109][110]), hydrogeological (e.g., [111,112]), and geotechnical engineering (e.g., [113][114][115]). Examples of PSO applications to slope stability problems include [107,[116][117][118][119][120][121][122]. In the literature, PSO's dynamics is governed by five principles of swarm intelligence: (a) proximity, i.e., ability to perform simple and time computation; (b) quality, i.e., ability to respond to quality factors in the environment; (c) diverse response, i.e., the method should not commit its activities along excessive narrow channels; (d) stability, i.e., the behavior of the method must not change with small changes in the environment; and (e) adaptability, i.e., the method must be able to alter its behavior when the computational cost is not excessive [104,122,123].
PSO reaches its goal if it meets the termination criteria. The commonly used termination criteria are set as follows: (i) reaching a maximum number of iterations; (ii) finding a satisfactory solution; and (iii) achieving constant fitness for a certain number of iterations. These criteria are set to guarantee the completion of the iterative search process [118,119]. In essence, PSO uses a population of search points to probe the search space. The population is called the swarm and the search points are called particles. Each particle moves in the search space with adaptable velocity, recording the best position it has ever visited in the search space, i.e., the position with the lowest function value. The adaption of the velocity is based on the information coming from the particle itself, as well as from the rest of the particles. As each particle has a neighboring prescribed particle, the best position attained by any neighbor is communicated to the particle and influences the movements. For the mathematical formulations of PSO, readers are referred to [107,[116][117][118][119]122].

Adopted Methods
In this study, the adopted methods can be grouped into deterministic factors of safety evaluation, probability of failure by means of MCS and RSM methods, and comparison of the results. MCS is an ideal standard to replace other reliability methods. RSM determines the optimum condition of the model's input variable that leads to the maximum or minimum response within the region of interest. These methods were chosen because they are coded in most slope stability commercial software, and they are readily available for use in practice.

Limit Equilibrium Analysis with Monte Carlo Simulation and the Response Surface Method
There are various alternative methods that are available in this group. The main difference between different limit equilibrium methods is in the assumptions made about the shape of the slip surface and the equilibrium equation that can be satisfied. As part of this effort, the factor of safety and probability of failure were evaluated by adopting the cuckoo search (CS) optimized slip surface in the Spencer method in Slide2. The CS is a very fast and efficient global optimization method, which is used in Slide2 for locating critical non-circular slip surfaces and hence locates a lower factor of safety than other methods, such as PSO. It requires no user input of trial surfaces or search objects.

Finite Element Shear Strength Reduction Method
The finite element elasto-plastic analyses assess the magnitude of deformation. The RS2 two-dimensional mode yields a deterministic factor of safety by means of the shear strength reduction (SSR) technique, during which the cohesion and friction angle of linear materials and the shear strength envelope of nonlinear materials are simultaneously reduced by a reduction factor until numerical convergence within the specified tolerance is no longer possible. The greatest SSR factor that allows convergence is considered the factor of safety against slope instability. The finite element method more realistically models actual failure mechanisms by allowing the failure surface to implicitly emerge as strain occurs within the continuum during the shear strength reduction process.

Application to Case Study
To demonstrate the application of probability methods, a case study was examined for an open pit mine. The case mine is a gold mining operation located in the Goldfields region of Western Australia. Slope stability analyses for a large open pit mine are relevant, such as increasing the slope angles of the existing design. A comprehensive slope stability project was conducted to determine the engineering properties of the rock mass to assess the failure mechanism and investigate alternatives for improving the overall stability of the slopes. The slope stability analysis was conducted for the eastern slope of the pit. The material in this location comprises three types, namely volcaniclastic sediment, porphyry, and basalt.

Results
The mean values of the material properties for the rock types are presented in Table 2.
The expected values and standard deviations were determined for each parameter by analyzing the geotechnical test results ( Table 2). All variables of the slope were found to be normally distributed, so all the analyses were based on normal distribution data. The mechanical properties of the rock mass were estimated from Hoek-Brown investigation. The results of the laboratory test show that the rock strength values, as obtained from the uniaxial compression strength test for the volcaniclastic sediment, basalt, and porphyry, are 140, 168, and 215 MPa, respectively. The unit weight of the rock mass is 27.4 kN for volcaniclastic, 28.8 kN/m 3 for basalt, and 26.3 kN/m 3 for the porphyry. The geological strength index (GSI) for the rock mass is estimated to be 58 for volcaniclastic, 68 for basalt, and 61 for porphyry. The mi is a constant estimation that is related to rock types and is based on considerations from RocData. In addition, the disturbance factor was taken as 1 for conventional production (poor) blast and the effect of the slope height was also considered in the calculation of the rock mass where the overall height used was 280 m. Given the high quality and strength of the three rock types, and the structural conditions, rock mass failure was considered the most likely failure mechanism. The traditional Hoek-Brown failure criterion was considered for the rock mass.

Discussion
To reduce the number of random variables, UCS and GSI were assumed to have engineering significance and were treated as random variables. All other parameters were treated deterministically, and their mean values used. From field observation, and based on the geological cross-section, no anisotropy matching fault/shear was included in the model. The location of the groundwater table in the model was based on water level readings from standpipe piezometers and the pressure head from a vibrating wire piezometer as obtained from the site.
To ensure the slip surface is global, the maximum iteration in CS was 500, the number of nests was 50 with an initial number of surface vertices of 8, and analyzed for the noncircular mode of failure. The analysis was performed to verify the performance of CS in probabilistic Monte Carlo simulation and the response surface search method. Table 3 shows the minimum factor of safety and probability of failure with the two iterations for the slope. This approach finds the optimal solution. It is clear from Table 3 that CS with RSM could be used in analyzing slope stability. The respective factor of safety (FS) of the overall slope is 2.59 and 2.90 using Monte Carlo simulation and the response surface search, and 0% probability of failure (PF). This suggests the slope is stable ( Figure 5). model. The location of the groundwater table in the model was based on water lev ings from standpipe piezometers and the pressure head from a vibrating wire pie as obtained from the site. To ensure the slip surface is global, the maximum iteration in CS was 500, the of nests was 50 with an initial number of surface vertices of 8, and analyzed for circular mode of failure. The analysis was performed to verify the performance probabilistic Monte Carlo simulation and the response surface search method. Table 3 shows the minimum factor of safety and probability of failure with iterations for the slope. This approach finds the optimal solution. It is clear from that CS with RSM could be used in analyzing slope stability. The respective factor (FS) of the overall slope is 2.59 and 2.90 using Monte Carlo simulation and the r surface search, and 0% probability of failure (PF). This suggests the slope is stable 5).

Conclusions
In this paper, an overview of probabilistic reliability methods for slope stability calculation considering material uncertainty was presented. The study highlights the need for a conscientious understanding of uncertainty and variability for effective representation of rock slope design. Uncertainty is a common fact in geological engineering problems and three types of uncertainty are normally identified for slope design. These are geological uncertainty related to natural variability in the rock properties, which cannot be reduced no matter our knowledge and expertise displayed in estimating them; model uncertainty related to the possibility of obtaining an incorrect result even if all the exact values are available for all model parameters; and parameter uncertainty related to the inability to account for the various characteristics of the geotechnical model that stem from either data scatter and/or systematic error, such as statistical and bias in measurements. These uncertainties can have a significant impact on the design performance of a slope if not properly accounted for.
There are two main approaches to deal with uncertainties. These are reliability and non-deterministic methods. There are two categories of non-deterministic probabilistic methods and non-probabilistic methods. Probabilistic methods are commonly used to represent and quantify uncertainty in the slope design process. Likewise, the reliabilitybased design approach quantifies and provides a consistent measure of safety by determining the probability of failure for a slope system. Even though the reliability-based methods may appear complex, they have been logically applied to a variety of rock slope engineering problems. For probability methods, the techniques have been coded in limit equilibrium and finite element computer software. However, as results may differ depending on the reliability (or probability) method chosen, two or more methods should be used to gain an understanding of the errors involved.
A case example was presented to demonstrate the value of probabilistic-based analyses in the slope design process. A rock slope from an open pit in a gold mining operation was studied. Geotechnical parameters are considered significant factors and any inconsistencies in these parameters and difficulty in the selection of appropriate data are paramount in rock slope design. The uniaxial compression strength (UCS), geological strength index (GSI), disturbance factor (D), rock mass constant (mi), and unit weight were determined to be key parameters. UCS and GSI were treated as random variables and a normal distribution was chosen. The Monte Carlo simulation and response surface methods implemented in SLIDE software and solved with cuckoo Spencer limit equilibrium methodology were used for this purpose. The study highlights the need for verification with deterministic finite element methodology in RS2 software. The results

Conclusions
In this paper, an overview of probabilistic reliability methods for slope stability calculation considering material uncertainty was presented. The study highlights the need for a conscientious understanding of uncertainty and variability for effective representation of rock slope design. Uncertainty is a common fact in geological engineering problems and three types of uncertainty are normally identified for slope design. These are geological uncertainty related to natural variability in the rock properties, which cannot be reduced no matter our knowledge and expertise displayed in estimating them; model uncertainty related to the possibility of obtaining an incorrect result even if all the exact values are available for all model parameters; and parameter uncertainty related to the inability to account for the various characteristics of the geotechnical model that stem from either data scatter and/or systematic error, such as statistical and bias in measurements. These uncertainties can have a significant impact on the design performance of a slope if not properly accounted for.
There are two main approaches to deal with uncertainties. These are reliability and non-deterministic methods. There are two categories of non-deterministic probabilistic methods and non-probabilistic methods. Probabilistic methods are commonly used to represent and quantify uncertainty in the slope design process. Likewise, the reliability-based design approach quantifies and provides a consistent measure of safety by determining the probability of failure for a slope system. Even though the reliability-based methods may appear complex, they have been logically applied to a variety of rock slope engineering problems. For probability methods, the techniques have been coded in limit equilibrium and finite element computer software. However, as results may differ depending on the reliability (or probability) method chosen, two or more methods should be used to gain an understanding of the errors involved.
A case example was presented to demonstrate the value of probabilistic-based analyses in the slope design process. A rock slope from an open pit in a gold mining operation was studied. Geotechnical parameters are considered significant factors and any inconsistencies in these parameters and difficulty in the selection of appropriate data are paramount in rock slope design. The uniaxial compression strength (UCS), geological strength index (GSI), disturbance factor (D), rock mass constant (mi), and unit weight were determined to be key parameters. UCS and GSI were treated as random variables and a normal distribution was chosen. The Monte Carlo simulation and response surface methods implemented in SLIDE software and solved with cuckoo Spencer limit equilibrium methodology were used for this purpose. The study highlights the need for verification with deterministic finite element methodology in RS2 software. The results were compared and were useful to highlight the benefit of probabilistic analysis of slope designs over the deterministic method.
Geotechnical engineers must be aware that probabilistic slope design methods are evolving and that the deterministic methods are not always consistent slope design methods. They must be aware that when using the deterministic method, significant design errors can occur if the design technique fails to represent the true physical behavior of the slope or rock mass conditions, such as calculating the factor of safety in deterministic limit equilibrium methods or finite element methods. Based on these conditions and for the purposes of expedience, it is necessary to integrate deterministic analysis with probabilistic design analysis to give a credible economic slope design. Slope design methods that use probabilistic analysis give a wider review of failure probability and risk-based decisionmaking, which is warranted in complex mining situations. Despite their mathematical complexity, the probabilistic-reliability methods discussed are in use today and due to progress in computation, they provide an invaluable reference to decision-making.