The Ordinary method ignores the shear forces and the interslice normal forces while satisfying the moment equilibrium for a circular slip surface. This method’s benefit is that it does not require an iteration procedure, making it simple to solve the FoS. The FoS is based on moment equilibrium and computed by Abramson et al., (2002) and Nash (1987) shown in Equations 2, 3: F o S = ∑ c ′ l + N ′ ⁡ tan ⁡ ϕ ′ ∑ W Sin α (2) N ′ = W ⁡ cos ⁡ α − u l (3) where, c′ = cohesion, u = pore pressure, l = slice base length and α = inclination of slip surface at the middle of slice.

It is assumed that the soil mass fails due to rotation on a circular slip circle centered on a point (Zhang and Li, 2021; Bishop, 1955). The normal force is taken to apply at the center of the base of each slice, and the shear stress between the slices is ignored because the forces on the sides of the slices are taken to be horizontal. The simplified Bishop’s approach produces reasonably precise FoS values, but it does not satisfy the entire static equilibrium (Verma et al., 2013). The Bishop method suggests that interslice shear forces can be neglected. The Bishop approach proposes that interslice shear forces can be ignored. The FoS can be computed using Bishop’s approach using Equation 4. F o S = ∑ c ′ l ⁡ cos α + P − r u l ⁡ cos α tan ϕ ′ cos α + sin α tan ϕ ′ / F o S ∑ W ⁡ sin α (4) where c′ = cohesion, l = width of slice, α = angle at the base of sliding slice, ϕ′= angle of internal friction, r _u = Pore water pressure, W = effective weight of the slice, p = Normal force acting at the base of the slice, and FoS = Factor of safety.

Figure 1 depicts a simplified schematic diagram of circular failure of slopes using Bishop’s method, as well as the forces acting on a single slice.

Janbu’s approach is suitable to non-circular soil mass failure. Similar to Bishop’s method, shear stress acting between the slices is ignored because horizontal forces should exist on the sides of each slice, and the normal force is assumed to act at the center of each slice’s base. Janbu’s corrected approach, which took into account the inter-slice shear force, was applied by Verma et al. (2013). The Janbu approach suggests that the interslice forces are normal. Equation 5 provides the Janbu expression for FoS. F o S = c ′ l + P − r u l tan ϕ ′ sec α ∑ W ⁡ sin α (5)

Figure 2 depicts the slice approach for analyzing non-circular failures on slopes using Janbu’s method and the forces acting on a single slice. Researchers employed non-circular failure analysis as well as Janbu’s corrected approach for circular failure in slopes.

The primary purpose of Lowe–Karafiath’s Method (Lowe, 1960) is to analyze and assess the safety of slopes, embankments, and earth structures by calculating the factor of safety against potential sliding along a slip surface. In FoS computation, Lowe-Karafiath’s technique (Lowe, 1960) satisfies only force equilibrium. Lowe-Karafiath’s technique, like other methods, makes the assumption that the inclination of the interslice force is equal to the average of the inclinations of the slope surface (β) and the slice base (α), that is, θ = 1 2 α + β where θ is the inclination of the interslice resultant force. Hence, the interslice forces can be expressed in Equation 6 as follows: T = E ⁡ tan ⁡ θ (6)

The Corps of Engineers approach is similar to Lowe-Karafiath’s method, with the exception of the assumption of interslice force inclination. The Corps of Engineers method is used to determine the FoS for slope stability. There are two possible approaches to assume the angle of the interslice resultant force using this method. Initially, it can be presumed that it is parallel to the ground, meaning that θ = β, where β represents the slope angle. Secondly, it can be equivalent to the average angle of slope between the critical shear surface’s entry and exit points.

Sarma’s Method Sarma (1973) is a well-established technique for determining the FoS for slope stability, especially in complex scenarios involving non-circular slip surfaces. Both of the equilibrium requirements are met by this approach (Abramson et al., 2002). Additionally, a linear Mohr-Coulomb expression (Equation 7) is assumed for the interslice force interaction. T = c h + E ⁡ tan ⁡ ϕ (7) where, c and ϕ are the shear strength parameters, and h is the slice height.

The Morgenstern-Price method, which assumes the interslice force function, also fulfills both force and moment equilibriums. The interslice force inclination can vary with an arbitrary function (f(x)) as per the Morgenstern-Price approach (Morgenstern and Price, 1965) and presented in Equation 8: T = f x ⋅ λ ⋅ E (8)

Where, λ = the scaling factor of the assumed function and f(x) = the interslice force function, which varies continuously through the slip surface. The approach proposes assuming any kind of force function, such as user-defined, half-sine, or trapezoidal. The base normal force (N) and the interslice forces (E, T) have the same relationships as those found in Janbu’s generalized technique. The interslice forces for a particular force function are calculated iteratively until F _f equals F _m in the Equations 9 and 10 (Nash, 1987). F f = ∑ c ′ l + N − u l tan ϕ ′ sec α ∑ W − T 2 − T 1 tan α + ∑ E 2 − E 1 (9) F m = ∑ c ′ l + N − u l tan ϕ ′ ∑ W ⁡ sin α (10)

Spencer’s approach is identical to the Morgenstern-Price method, with the exception of the assumption for interslice forces. For interslice forces, a constant inclination is assumed, and the FoS is calculated for both equilibrium (Spencer, 1967). Using this approach, the interslice shear force is associated with Equation 11: T = E ⁡ tan ⁡ θ (11)

The friction circle method is part of the broader category of limit equilibrium methods used for analyzing slope stability. This method divides the slip surface into vertical slices and evaluates forces and moments for each slice to determine equilibrium (Terzaghi et al., 1996). It calculates the FoS by comparing resisting forces, derived from soil friction along the slip surface, to driving forces. The method is particularly useful for preliminary assessments of slopes with cohesive soils, where circular failure surfaces are likely. While it simplifies calculations, it assumes uniform soil properties and a circular failure surface, which may not always reflect real conditions.

Taylor’s method is applicable for circular failure of the soil mass (Taylor, 1937). The assumption that soils are homogeneous and isotropic is incorrect in practice (Sakellariou and Ferentinou, 2005). Researchers previously used Taylor’s slope stability chart to estimate the factor of safety of slopes having simple geometry with homogenous and isotropic soil properties in clays under single-value undrained shear strength (Javankhoshdel and Bathurst, 2014; Liu et al., 2022). Taylor (1937) developed charts for calculating the factor of safety for simple slopes of cohesive-frictional (c-ϕ) shear strength soils. The disadvantages of Taylor’s chart include the requirement for repeated processes to determine the FoS. Equation 12 can be used to compute the FoS. Taylor’s slope stability chart is shown in Figure 3. F o S = S u γ H N s (12)

where the depth factor (D) and slope angle (α) determine the stability number (N _s ). In Figure 3, “H” represents the slope’s height and “DH” the hard stratum’s depth from the slope crest. It is generally accepted that “H” and “α” are deterministic. Additionally, S _u , a lognormally distributed or random variable, and γ, a constant, or both, can be used to compute the probability of failure for lognormal distributions with uncorrelated random variables.

The Linear Surface Failure Method, used for analyzing infinite slopes in granular soils, assumes a planar failure surface. It is particularly effective for evaluating stability in granular materials like sand where cohesion is minimal (Duncan et al., 2014). The method involves calculating the factor of safety (FS) by comparing the resisting forces (due to friction) to the driving forces (due to gravity). This approach simplifies analysis by assuming an infinite slope and a linear failure surface, making it suitable for preliminary assessments of steep slopes.

Table 2 contains detail comparison of various limit equilibrium methods used in slope stability analysis, including the Ordinary Method, Bishop’s method, Jambu’s method, Corps of Engineers method, Sarma’s method, Lowe and Karafiath’s method, Morgenstern and Price’s method, Spencer’s method, Finite element method, Taylor’s method, and Linear surface failure method.

In the ANN mathematical model, neurons are regarded as processing elements or nodes. A network with input vectors x ₁ , x ₂ , … , x _m and output vectors y ₁ , y ₂ , … , y _p as shown in Figure 5, can be expressed in Equation 13 as follows: y j = f o u t p u t ∑ h = 1 k w j h f h i d d e n ∑ i = 1 m w h i x i + w h 0 + w j 0 (13) where f _hidden and f _output represent the respective transfer functions for the hidden and output layers. w _hi and w _jh represent the neuron weights of input neuron i to hidden neuron h and hidden neuron h to output neuron j, respectively. x _i is the i ^th input unit; k is the number of neurons in the hidden layer; m is the number of input units; p is the number of output units; w _ho is the threshold (or bias) for neuron h; w _jo is the threshold for neuron j.

Neural network models can be classified as feed forward or feedback based on their topology. The back propagation (BP) and radial basis function (RBF) algorithm models are commonly used in soft computing prediction models.

Zhou et al. (2019) employed ANN to predict slope stability, with a 0.827 accuracy. Similarly, Qi and Tang (2018) used ANN to predict slope stability and found an accuracy of 0.84.

The optimal classification hyper-plane problem with linear separability is the basis for the SVM method (Samui, 2008; Zhu and Zhang, 2004). Maximizing the training set interval is the algorithm’s objective. Nonnegative relaxation variables ξ n , n = 1 , 2 , . . . , k are necessary for generalized optimum classifications under conditions that are not linearly separable. Then, the Equations 14, 15 present classification hyper-plane optimization problem: min w , b = w 2 2 + C ∑ i = 1 n ξ i 2 (14)

which subject to y i w × x i + b ≥ 1 − ξ i ξ i ≥ 0 , i = 1 , 2 , . . . , n (15)

where C ∈ 0 , + ∞ is the parameter for penetration, w is the normal vector, and b is the hyper-plane’s bias. Through nonlinear transformations, the SVM algorithm can convert the input space into a high-dimensional space. The kernel function K x i , x j = Φ x i × Φ x j can be used to define the nonlinear transformation. There are four kernel functions, with the most popular being the radial basis function (RBF) kernel. The function of the radial basis function kernel is given in Equation 16: K x i , x j = ⁡ exp − σ x i − x j 2 (16) where σ is the essential kernel parameter that specifies the kernel’s width, and σ ∈ 0 , 1 .

The SVM for slope stability was applied by Lin et al. (Lin et al., 2018) with six typical slope parameters: pore water ratio, unit weight, cohesion, internal friction angle, slope angle, and slope height and found an accuracy of 0.6667.

An ensemble learning technique called RF is used to address regression and classification issues. Breiman in 2001 developed a combinatorial classification method called the RF (Breiman, 2001). It is comprised of many decision-tree classification models Θ h X , k , k = 1 , 2 , . . . , where the input vector and the parameter set is a random vector that is independently and identically distributed. The inconsistent results are decided by the voting process. The following is the main idea of RF: The original sample set is resampled to create numerous sample sets, each of which reflects the entirety of the training data for a given category tree. When each sample set evolves into a tree, m _try properties are selected at random from the M properties at each node. Once the purity of every node reaches its lowest point, one of the m _try attributes is chosen to mature. The generated n _tree trees can then be used to build the RF classifier, which can then be used to classify the newly data set. The number of votes by tree classifiers influences the classification result and presented in Equation 17 is: c = a r g m a x c ∑ k = 1 n t r e e I h X , Θ k = c n t r e e (17) where c is the result of the classification and n _tree is the number of trees. Figure 6 depicts the RF functional architecture.

Zhou et al. (2019) found that the accuracy value was 0.808 after applying RF to 221 datasets for slope stability. In order to determine slope stability, Lin et al. (2018) employed RF on 132 datasets, yielding an accuracy of 0.83. In the same way, Qi and Tang (2018) used RF to analyze the soil slope stability of 148 slope instances and achieved an accuracy of 0.93.

Rasmussen’s hypothesis—that neighboring observations should to exchange information—forms the foundation of the GPR model (Rasmussen and Williams, 2006). A joint multivariate Gaussian distribution can be found for any finite number of the random variables in a gaussian process. Let a×b stand for the input and output domains, respectively, from which n uniformly and randomly distributed pairs (a _i , b _i ) are taken. For regression, let b ⊆ R ; then, a Gaussian process on a is distinct by the mean function μ : a → R and a covariance function k : a × a → R . The main supposition of GP regression is that y is given as b = f a + ξ , where ξ ∼ N 0 , σ 2 . In GPR, for every input x, there is a corresponding random variable f(a), which is the value of the stochastic function f at that location. That is, In GPR, there is a corresponding random variable f(a)that represents the value of the stochastic function f at that specific position for each input x. In this work, it is assumed that the observational error n is normal independent and identically distributed, with a mean value of zero μ a = 0 , a variance of σ 2 and f(a) drawn from the Gaussian process on a specified by k. This is given in Equation 18, B = b 1 , . . . , b n ∼ N 0 , K + σ 2 I (18) where K _ij = k (ai, aj) and I is the identity matrix. As B / A ∼ N 0 , K + σ 2 I is normal, so is the conditional distribution of test labels given the training and test data of B / B , A , A * . Then, one has B * / B , A , A * ∼ N μ , ∑ where μ and ∑ are given in Equations 19 and 20: μ = K A * , A K A , A + σ 2 I − 1 B (19) ∑ = K A * , A * − σ 2 I − K A * , A K A , A + σ 2 I − 1 K A , A * (20)

This is also true for the other values of K(A, A), K(A*, A), and K(A*, A*); here, A and B are the vectors of the training data and training data labels b _i , whereas A* is the vector of the test data. If there are n training data and n* test data, then K(A, A*) represents the n × n* matrix of covariance, which is evaluated at all pairs of training and test datasets. Tien Bui et al. (2019) used GPR for predicting slope stability FoS and found that the R² value is 0.9467.

DT addresses both regression and classification issues using a structure like a tree (Kardani et al., 2021; Pekel et al., 2020). Nodes in a DT model that have outgoing edges are called internal nodes, while nodes without any edges are called leaf nodes (leaves). DT models also use branches to connect the nodes. Classification trees and regression trees are the two forms of DTs. Classification trees divide data into two subsets according to class labels, and they keep doing this until a stopping criteria is reached. Another kind of DT used in machine learning to address regression issues is regression trees. They are used for predicting continuous output variables, as compared to classification trees, which are used to predict a discrete set of values. Regression trees split data into two groups and continue doing so until a stopping condition is met. DTs have numerous advantages, including ease of understanding and explanation (Witten and James, 2013). But DTs can be extremely sensitive. Minor adjustments to the input data can have considerable effects on the trees and the results (Witten and James, 2013). Because of their reliance on the greedy algorithm, DTs may occasionally fail in achieving a globally optimal outcome (Ben-Gal et al., 2014). A tree structure may effectively show the relationship between binary dependent variables and related independent variables, which is one of the advantages of DT analysis. Decision trees are hence a popular and effective data mining technique (Duch et al., 2004). Amirkiyaei and Ghasemi (2022) used J48 to assess the stability of slopes with 92.3% accuracy.

The k-NN approach is a well-known soft computing technique that is used to solve regression and classification issues. In the feature space, it is assumed that similar samples are frequently found near to one another (Peterson, 2009). The k-NN technique determines the k closest samples within the dataset for each sample point and computes their distance from the previous sample points. The samples that satisfy k with the smallest distances comprise the input data. The outcome of classification problems is class membership (labels). The majority class among the k closest samples is used to classify a sample. The sample is assigned to the closest class if k = 1. The average value of the closest samples in a regression issue is the output value (Mahmoodzadeh et al., 2022). The k-NN method’s quick and easy training process is one of its benefits. Furthermore, the k-NN method can reduce noise to some degree (Bhatia, 2010). But as k is a sensitive parameter, a smaller value could result in overfitting, while a bigger value could result in underfitting. Furthermore, processing costs might be high, particularly for large datasets (Bhatia, 2010). To ensure performance, hyperparameters must be properly chosen. Image recognition is one of the k-NN method’s application scenarios (Zhang et al., 2006; Homaeinezhad et al., 2012), text classification (Trstenjak et al., 2014; Tan, 2006), and recommendation system (Adeniyi et al., 2016). Huang et al. (2020) used k-NN algorithm for slope stability to achieve prediction accuracy up to 92.30%.

The MLR, often known as multiple regression, is a simple machine learning technique that explains the relationship between two or more predictor variables and a response variable. The model is shown in Equation 21 (Pandey et al., 2022): y = α 0 + α 1 × x 1 + α 2 × x 2 + . . . + α n × x n + ξ (21) where y is the dependent variable (output data), α ₀ is the intercept or constant term, and α ₁ to α _n are the coefficients of the independent variables (input data) x ₁ to x _n , respectively. The MLR aims to minimize the sum of squared errors between the dependent variable’s actual values and its predicted values by estimating the values of the α coefficients. The fact that MLR can only take into account a linear relationship between the input and output data is one of its drawbacks. Chakraborty and Goswami (2017) used MLR to predict FoS of slope stability using 200 cases and found that R² is 0.847.

MARS is a nonlinear regression approach that is used for classification and regression applications. MARS uses a collection of piecewise linear functions known as basis functions to model the relationship between several inputs and one output (Falae et al., 2021; Zhang and Goh, 2016). There are similarities between other basis functions, using the linear basis function as an example. The linear basis function Bi(x) can be defined in Equation 22 as three types:

(1) Constant values T (the intercept).

MARS model M can be expressed as a sum of basis functions given in Equation 23: M = ∑ i = 1 t t i B i x (23) where t _i is the constant coefficient. The model iteratively adds new functions after beginning with a single basis function. The algorithm chooses each hinge’s ideal location in a stage-by-stage, greedy fashion (Friedman and Roosen, 1995). Using a limited number of basis functions, MARS has the advantage of developing a straightforward model that captures intricate interactions between inputs and outputs. On the other hand, MARS might be affected by the initial basis function and stopping criteria selected. MARS is often used in credit scoring (Lee et al., 2006; Lee and Chen, 2005), and species distribution models (Leathwick et al., 2006) and also has applications in time series analysis (Lewis and Stevens, 1991). Liao and Liao (Liao et al., 2020) used MARS to predict FoS of slope stability with R² value is 0.8629.

To increase prediction accuracy, GBM is an ensemble of weak prediction models, like DTs. GBM constructs the model step-by-step and, because it can maximize any differentiable loss function, has far wider applicability than previous boosting techniques. GBM is an accurate and effective technique that works well for both regression and classification issues. Web search ranking and ecology are two areas of research where GBM modeling has been widely applied. The benefits of GBM include its ability to handle mixed-type data naturally, its strong prediction ability, and its proficiency in handling output space outliers. In Civil Engineering, GBM has been employed for predicting hanging wall stability (Qi et al., 2018). Zhou et al. (2019) found that 0.865 was the accuracy after using GBM to 221 datasets for slope stability.

The ANFISs are a soft computing technique that combines ANNs with Takagi-Sugeno fuzzy inference systems (Jang, 1993; Jang, 1991; Mu’azu, 2023). The essential component is a collection of fuzzy rules—conditional statements that read, “if x is A, then y is B”—that are used to simulate a system’s input–output relationship. An ANFIS uses backpropagation and gradient descent methods to train a model in about five stages (Jang, 1993):

(1) Fuzzy sets are created from the input data by applying membership functions.

ANFIS’s main benefit is its ability to express structured knowledge and nonlinearity (Jang, 1993). To train the model, ANFIS needs an adequate dataset; yet, choosing the right input-output data is far more important. Additionally, the number of input data increases exponentially with the number of fuzzy rules, thus raising the cost of computation and may effect a model’s performance (Al-Mahasneh et al., 2016). Mohamed et al. (2012) used ANFIS and obtained R² of 0.9997 to predict FoS with high accuracy compared with MLR.

Some other soft computing methods are sparse polynomial chaos expansion is based on polynomial chaos expansion and hybrid soft computing models for predicting slope stability.

The use of soft computing models for slope stability analysis offers distinct advantages in terms of accuracy, flexibility, and efficiency, particularly in addressing the complexities and uncertainties associated with real-world slope conditions. On the other hand, numerous studies have been undertaken in recent years to develop soft computing models for predicting slope stability status (see Table 3). For example, Zhou et al. (2019) used the gradient boosting machine (GBM) method to analyze slope stability and compared to the artificial neural network (ANN), random forest (RF) and support vector machine (SVM). It was found that the GBM has highest accuracy i.e., 0.865. Sari (2018) used binary logistic regression analysis as an alternative technique to predict stability condition of slopes with 90.9% accuracy. Accordingly, Qi and Tang (2018) used basic ensemble classifier consisting of six individual classifiers and majority voting as a combination method to improve analysis of slope stability. It is worthwhile to note that the coefficient of determination (R²)/accuracy in Table 4 ranges from 0.556 to 0.96. Although, these models have greatly improved our understanding; however, the performance and accuracy of the predictive models are not well interpreted, and the problem of slope stability is still far from being fully solved. Therefore, the issue of slope failure still posed considerable challenge for geotechnical professionals. The existing slope stability issues and damages caused by landslides, it is necessary to put forward the idea that more systematic and in-depth research should be carried out on predicting the slope stability. Consequently, there is a need for more accurate and reliable methods to predict slope stability. Soft computing methods provide powerful tools for tackling complex, uncertain, and imprecise problems. Their flexibility, adaptability, and robustness make them suitable for a wide range of applications. Soft computing techniques are known for their proficiency in non-linear modeling, and there is evidence in the literature from a number of technical and scientific fields that these techniques can establish correlations between desired outcomes and a variety of influencing parameters, whether those parameters have direct or indirect impacts (Asteris et al., 2022; Koopialipoor et al., 2019). However, choosing a suitable soft computing model is challenging for the reasons listed below: (a) Inadequate modeling and validation; (b) models in use not being able to identify the precise global optimum; (c) problems with overfitting, etc. However, the challenges related to computational cost, interpretability, and data dependence must be carefully managed to fully leverage their potential.