Abstract
This paper aims to investigate the differences in factor of safety (FS) and failure mechanism (FM) for spatially variable undrained soil slope between using finite element method (FEM) , finite difference method (FDM), and limit equilibrium method (LEM). The undrained shear strength of cohesive soil slope is modeled by a one-dimensional random field in the vertical direction. The FS and FM for a specific realization of random field are determined by SRT embedded in FEM- and FDM-based software (e.g., Phase2 6.0 and FLAC) and LEM, respectively. The comparative study has demonstrated that the bishop method (with circular failure surface) exhibits performance as fairly good as that of SRT both in FS and FM for the undrained slope cases where no preferable controlling surfaces such as hydraulic tension crack and inclined weak seams dominate the failure mechanism. It is, however, worthwhile to point out that unconservative FM is provided by the Bishop method from the aspect of failure consequence (i.e., the failure consequence indicated by the FM from the Bishop method is smaller than that from SRT). The rigorous LEM (e.g., M-P and Spencer method with noncircular failure surface) is not recommended in the stability analysis of spatially variable soil slopes before the local minima and failure to converge issues are fully addressed. The SRT in combination with FEM and/or FDM provides a rigorous and powerful tool and is highly preferable for slope reliability of spatially variable undrained slope.
1. Introduction
Both the failure probability and the respective failure consequence are involved in the landslide risk assessment [1–5]. The shallow slide and deep slide as shown in Figure 1 clearly demonstrate the difference in failure consequence. It is presumed that deep slide leads to a more severe consequence than the shallow one with other things being equal. Although the failure consequences are in general site-specific (e.g., loss of life and property), the volume of sliding mass (i.e., failure mechanism) can be a preliminary measure of failure consequence [4–6].

(a)

(b)
The failure probability and failure mechanism can be calculated in a probabilistic manner combing the finite element method (FEM) [3, 5, 7–15], finite difference method (FDM) [16] with the strength reduction technique (SRT), and limit equilibrium method (LEM) [17–27]. On the one hand, FEM and FDM are preferably adopted for deterministic slope stability evaluation owing to their abilities of automatically searching the critical failure mechanism, modeling complicated constitutive law, and considering sophisticated external loads. On the other hand, however, the extensive computation effort within the implementation of FEM and FDM in slope reliability analysis always enforces the geotechnical practitioners to the use of simplified LEM. It is, therefore, of practical interest and importance to know the difference in factor of safety (FS) and failure mechanism (FM) among the aforementioned three methods especially for the probabilistic slope stability problem taking into account the spatial variability of soil properties.
The objective of this paper is to investigate the FS and FM obtained from different numerical procedures of FEM, FDM, and LEM as well under specific realizations of random field for undrained shear strength in cohesive soil slope (i.e., undrained soil slope). The comparison is expected to provide some insights into the applicability of abovementioned numerical procedures in reliability analysis for spatially variable undrained soil slope. The paper starts with a brief introduction of probabilistic characteristics of spatially random soil, followed by the description of an illustrative example slope. Then, the results from FEM, FDM, and LEM are presented and compared. Finally, the conclusions and discussion are given, and some insights into the applicability of FEM, FDM, and LEM in reliability analysis of spatially variable undrained soil slope are provided.
2. Probabilistic Characteristics of Spatially Random Soil
Soil properties exhibit spatial variability even in a homogeneous soil layer [28]. A straightforward and well-recognized approach to modeling the inherent spatial variability is stationary random field. Within the context of the stationary random field, the spatially variable soil property (e.g., undrained shear strength, Su) is characterized by the mean value μ, standard deviation σ, scale of fluctuation λ, and autocorrelation function. In this study, Su is assumed to follow a lognormal distribution, which guarantees that all the realizations of sampling are nonnegative and benefits from a simple transformation to the classical normal (Gaussian) distribution [2]. Based on this assumption, the lognormal random field of Su can be easily transformed to a normal random field of lnSu. The mean value and standard deviation of lnSu, denoted as μln and σln, are determined, respectively, in [2]:
Scale of fluctuation λ indicates a measure of distance within which soil properties at different locations are highly correlated. Let lnSu(i) and lnSu(j) denote the value of lnSu at location i and j, respectively; dij denotes the vertical distance between locations i and j. The simple exponential correlation function is adopted herein to simulate the correlation coefficient ρij between lnSu(i) and lnSu(j) [22, 29, 30]:
There are several available random field generation techniques, such as the local average subdivision method [31], Fourier series method [32], and covariance matrix decomposition method [22, 30, 33–38]. Consider, for example, the covariance matrix decomposition method is used to generate the random field of Su represented by a vector Su that comprised a set of component random variables Su(1), Su(2), …, Su(n). n is the properly selected number to discretize the continuous random field into a discrete one for practical use. Su can be written as follows:where I is the identity vector with n components all equal to unity, L is the a n-by-n lower triangular matrix determined from Cholesky decomposition of the correlation matrix ρ = [ρij]n×n = LLT, and ξ is the standard normal (Gaussian) vector with n independent elements.
A set of values for Su(1), Su(2), ..., Su(n) is a specific realization random file of Su. These values are taken as the input parameters in deterministic slope stability analysis procedures and FS, and the corresponding FM can be obtained under specific realization.
3. Illustrative Example
3.1. Input Parameters
The cohesive soil slope as shown in Figure 2 has a height of 10 m and a slope angle of 45°. The short-term shear strength of undrained soil slope is characterized by undrained shear strength Su, and the saturated unit weight of soil is γsat. In this study, γsat has a deterministic value of 20 kN/m3. Su is modeled by a lognormal random field varying in the vertical direction. The mean value μ and standard deviation σ are 40 kPa and 12 kPa, respectively. In order to discretize the continuous random field of Su, the homogeneous soil layer of 12 m thickness as shown in Figure 2 is equally divided into 12 sublayers. The Su at the top sublayer is Su(1), and similarly, that one of the bottom sublayer is Su(12). It is noted that the variance reduction (i.e., σln in equation (3) will be reduced in accordance with the discretization size) due to the local averaging effect should be properly accounted for within the discretization of random field [29]. Five values of scale of fluctuation are considered; that is, λ = 1 m, 2 m, 5 m, 10m, and +∞ (a sufficiently large value, e.g., 10,000,000 is used herein). For each of four λ values, 10 realizations of Su are generated using equation (3). Table 1 summarizes the 10 realizations of Su at each λ value. The first realization of Su (i.e., 12 Su values) at λ = 1m is shown in Figure 2 for illustration. For each of 10 realizations at each of λ values, the Su values at all the sublayers and γsat are fundamental input parameters for slope stability analysis. Besides, Young’s modulus E and Poisson’s ratio ν are required in the numerical procedure such as FEM- and FDM-based SRT to conduct deterministic slope stability analysis and to evaluate the FS and corresponding FM. Identical parameters of Young’s modulus E = 100 MPa and Poisson ratio ν = 0.3 as indicated in [39] are used in the ensuing studies.

3.2. Three Numerical Procedures for the Calculation of FS
The stability of undrained soil slope can be assessed by numerical procedures such as FEM and FDM in combination with SRT. For a stable slope, the SRT is implemented by gradually reducing the original strength parameters until the failure limit state is reached and the FS is defined as the factor by which the original strength parameters have been reduced. This technique has been embedded in several commercial numerical programs such as Rocscience Phase2 6.0, FLAC, Abaqus, and so on. Rocscience Phase2 6.0 (http://www.rocscience.com/rocscience/products/rs2) is a powerful 2D elastoplastic finite element stress analysis program for underground or surface excavations in rock or soil. It is employed in the present study to calculate the FS with the choice of FEM-based SRT. The left and right boundaries of the slope model are rollers, and the bottom boundary is a hinge in Phase2 6.0. Besides, as a mini-version of FLAC that is designed specifically to perform FS calculations for slope stability analysis, FLAC/Slope is used to obtain the FS based on FDM.
In addition to FEM and FDM approaches, the traditional method for FS calculation is LEM including, but not limited to, the Bishop method [40], Morgenstern–Price (M-P) method [41], and Spencer method [42]. The LEM-based commercial software Rocscience Slide 5.0 (http://www.rocscience.com/products/8/Slide) is adopted for deterministic slope stability analysis. The Bishop method is selected for circular failure surface, and the minimum FS is optimized using the auto refine search method, whereas M-P (with half sine interslice force function) and Spencer methods are selected for noncircular failure surfaces and the minimum FS is obtained by a path search with 10,000 potential trial surfaces in Slide 5.0.
It is worthwhile to note that the numerical models built in Slide 5.0 can be directly imported into Phase2 6.0 and there is no need to rebuild them, thereby improving the efficiency. This feature is one of reasons for the choice of Phase2 6.0 instead of Abaqus to conduct FEM-based deterministic slope stability analysis.
The minimum FS and the related FM (failure surface) in Slide 5.0 are compared with those obtained by Phase2 6.0 and FLAC/Slope. For the convenient comparison of FM, herein a toe failure mechanism refers to that goes through the toe but not tangent to the firm base, while the one that goes deep tangent to the firm base and outcrops at the foundation layer is named a deep failure mechanism. The contour of maximum shear strain in Phase2 6.0 and that of maximum shear strain rate in FLAC/Slope provide an approximate view of the FM, while the failure surfaces in Slide 5.0 are added to the contour for comparison.
3.3. Convergence Study on the Mesh Size in Phase2 6.0 and FLAC/Slope
The convergence study on the effect of mesh size on the FS and especially on the FM should be performed before the comparative study on the FS and FM between FEM, FDM, and LEM. Consider, for example, the first realization of Su is at λ = 1 m, the coarse mesh is first adopted to conduct deterministic slope stability analysis, and thereafter, the finer mesh is used. The negligible variations in FS and FM are achieved by gradually improving the mesh precision. As illustrated in Figure 3, the number of 6-node triangular elements increases from 423 to 21860 representing the coarsest mesh and the finest mesh, respectively. The FS corresponding to the coarsest mesh is 1.19, and after the number of elements is greater than 2718, FS = 1.18 does not change with the increasing mesh size, which indicates the trivial influence of mesh size on the FS. For the FM, as the mesh size increases from 423 to 11535, the failure region indicated by the contour of maximum shear strain tends to grow narrower, and after the number of elements is greater than 11535, slight difference in the failure zone is found. Therefore, if both the FS and FM are of interest, the number of elements discretizing the slope profile has to be sufficiently large (e.g., 11535 or 21860). In the ensuing comparative studies, 21860 triangular elements are used to determine the FS and the FM in Phase2 6.0. As shown in Figure 4, the number of zones used in FLAC/Slope is determined in a similar manner, and finally, 19675 zones are selected for the following numerical calculations. It is noted that since the main aim of Figures 3 and 4 is to demonstrate the variation of failure zones, the legends indicating the detailed magnitude of shear strain (for FEM) or strain rate (for FDM) in contours is omitted. For the similar reason, the legends in Figures 5–9 are not shown.

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(a)

(b)
4. Results
4.1. Comparison of FS
Table 2 summarizes the FSs from FEM (e.g., Phase2 6.0), FDM (e.g., FLAC/Slope), and LEM (e.g., Slide 5.0, with the Bishop method for circular failure surface, M-P method, and Spencer method for noncircular failure surface). Figure 10(a) compares the FS by FEM with that presented by FDM. It is evident that the FSs determined by the FEM plot closely to that obtained by FDM and the majority of data pairs lie on the solid 45-degree line. The correlation coefficient for the linear regression analysis of the data pairs is named r and is included in Figure 10. It is 0.999 for Figure 10(a) case, and it implies a trivial dispersion in FS between FEM and FDM. The similar trend and correlation coefficient r have been noticed from Figure 10(b) for the data pairs of FS between FEM and Bishop method. It is demonstrated that the Bishop method provides almost identical FSs as compared with those by FEM and FDM. However, significant differences in FS between FEM and M-P and between FEM and Spencer are noticed as shown in Figures 10(c) and 10(d). The correlation coefficient r is 0.85 and 0.935, respectively, for Figures 10(c) and 10(d). The reasons for this dispersion in Figures 10(c) and 10(d) are rather complicated. On the one hand, as Tabarroki et al. [39] pointed out, the choice of noncircular failure surface leads to a significant increase in complexity of the optimization problem as compared to the selection of the circular failure surface where only three optimized variables are of interest. With the increasing complexity, the search algorithm tends to locate the local minimum instead of global minimum (i.e., local minima issue). To properly address the local minima issue, several heuristic optimization algorithms such as improved Harmony search algorithm [43, 44] and discontinuous flying particle swarm optimization algorithm [45] have been proposed and used in deterministic slope stability analysis. The local minima issue arises in Slide 5.0 when the path search method with 10,000 potential failure surfaces cannot find the global minimum FS. It is therefore necessary to try several search methods and choose the minimum FS as the final solution. In this study, the block search method has been performed for the realizations where the FSs from LEM deviate significantly from those by FEM and/or FDM. However, despite the further effort, there still exists the cases where the FSs from M-P and Spencer methods are greater than those of Bishop (e.g., 4th, 5th, 9th, and 10th realization at λ = 1 m), which may be attributed to the issue of failure to converge while calculating the FS for a specified noncircular failure surface. As pointed out by Cheng et al. [46], the interslice force function should be varied according to the variation of the failure surface and input parameters for the proper calculation of FS. However, fixed interslice force functions like half sine (f(x) = sin(x)) or constant (f(x) = 1) are assumed, respectively, in M-P and Spencer method for all the potential failure surfaces. This assumption may provide a biased FS (unrealistically high or small) and, furthermore, can create prohibited search zones misleading the search direction. These two issues must be adequately addressed before the LEM can provide reliable solutions in slope stability analysis especially in spatially variable slope stability analysis. The Bishop method is manifested to yield comparable results with those from SRT because the aforementioned two issues are fully addressed. It is therefore expected similar failure probability results can be obtained by Bishop method and SRT for spatially variable undrained soil slope.

(a)

(b)

(c)

(d)
4.2. Comparison of FM
It is not straightforward for the comparison of FM as compared with that of FS because failure regions are indirectly indicated by the contour of maximum shear strain or shear strain rate in FEM and FDM, whereas the circular and/or noncircular failure surface is provided in LEM. A simple criterion is specified herein for comparison between failure region and failure surface. It is taken for granted that the more parts of the failure surface are within the failure region, to a greater extent the failure surface is comparatively well with the failure region. If the majority (e.g., more than 80%) of the failure surface is within the failure zone, the term well is applicable, whereas if more than 50% and lower than 80% of the failure surface is in the failure zone, the term moderate are pertinent to this case; otherwise, the term bad is applicable to this case.
Figures 5–9 plot the failure regions obtained by FEM and FDM, and the failure surfaces in LEM as well for each realization at λ = 1 m, 2 m, 5 m, 10 m, and +∞, respectively. The failure surfaces from LEM are separately incorporated into the contours from FEM and FDM for comparison for each of 10 realizations. Therefore, each of Figures 5–8 comprises 20 comparisons. Unlike the other four cases, λ = +∞ case represents no consideration of the spatial variability of Su yielding the same Su value at each of 12 sublayers. Although the FS with respect to each realization at λ = +∞ case is different, the FM remains unchanged for all the 10 realizations. Hence, only the comparison at one realization is plotted in Figure 9. Note that the FS and the calculation method are also included in Figures 5–9 with the form of the“FS-calculation method.” For example, “1.21-Bishop” means the FS of the plotted FM is 1.21 and it is calculated by the Bishop method. In addition, if Spencer is absent, this means both M-P and Spencer method have the same FS and FM. Based on the above criterion, the extents to which the failure surfaces by the LEM agree well with those contours from FEM and FDM are subjectively and qualitatively evaluated by the authors. Table 3 summarizes the qualitative results of the comparison in FM between LEM and FEM/FDM. To quantitatively assess the performances of different LEM methods, the qualitative results are specified to different values. For example, if a result is well and it will be replaced with a number of 80, and the moderate and bad results are replaced with 50 and 25, respectively. The average scores for Bishop, M-P, and Spencer methods are 75, 39, and 46, respectively, for 40 realizations at λ = 1 m, 2 m, 5 m, and 10 m. It is apparent that the FMs (circular failure surface) obtained by the Bishop method agree favorably well with those in FEM and FDM. M-P and Spencer methods do not perform as satisfactorily as Bishop for the case where spatial variability of Su is involved, although they display as good as performance for the case where spatial variability of Su is ignored. Although the performance of the Bishop method in the determination of FM is fairly acceptable, it is worthwhile to note that slightly considerable difference in the entry location of failure mechanism on the crest is noticed (e.g., 6th, 7th, and 10th realization at λ = 1 m; 5th, 8th, and 10th realization at λ = 2 m; 6th, 7th, and 8th realization at λ = = 5 m; and 5th∼10th realization at λ = 10 m). The Bishop method tends to provide a smaller volume of sliding mass than FEM- and FDM-based SRT, and therefore, it yields an unconservative FM from the aspect of failure consequence.
As for the FM type, either toe or deep or both is defined for each of 10 realizations at λ = 1 m, 2 m, 5 m, 10 m, and +∞, respectively. It is noted that both toe and deep failure mechanisms occur, and this phenomenon has been observed in many previous researches such as in [47–49]. The results are summarized in Table 2 (see columns 3 and 4). It is found that six out of ten realizations show deep failure mechanism at 1 m, 2 m, and 5 m. As λ increases, the number of deep failure mechanism increases. Eight out of ten realizations follow deep mechanism at λ = 10 m, and all the ten realizations exhibit deep failure mechanism for the λ = +∞ case. For the λ = +∞ case, the 12 sublayers in Figure 2 reduce to a homogeneous soil layer and either deep or toe failure mechanisms depend on whether the slope angle is lower than or greater than 53° for undrained soil slope. Therefore, the failure mechanism for undrained soil slope shown in Figure 2 (slope angle = 45°<53°) at λ = +∞ case falls into the category of deep failure mechanism. Zhu et al. [3] studied the effect of spatial variability on the failure mechanism for undrained soil slope, and they concluded that, for certain combinations of random field properties, relatively flat slopes (slope angles <53°) may display a significant number of toe mechanisms. This conclusion has also been verified in this study. Five out of ten realizations show toe failure mechanisms at λ = 1 m, 2 m, and 5 m, respectively. Only three out of ten realizations display toe failure mechanisms at λ = 1 m, and no toe failure mechanisms are present for the λ = +∞ case. The number of realizations showing toe failure mechanisms decreases with increasing λ values. The recognition of uncertainty of failure mechanisms is important because the consequences may be more serious in a deep failure as it involves a greater volume of soil.
5. Discussion and Conclusions
As a geotechnical routine, the slope stability evaluation is always performed by SRT in combination with FEM and/or FDM numerical procedures and by traditional LEM (e.g., Bishop method, M-P method, and Spencer method). The FEM- and FDM-based SRT is superior to LEM in that the former can easily account for the stress-strain curve of soils and be able to automatically locate the failure regions, thereby not requiring the prior assumption on the failure surface as compared to LEM. The key limitation of FEM- and FDM-based SRT when extended to probabilistic slope stability problem is computationally demanding as compared with LEM. The use of LEM in probabilistic slope stability analysis is therefore found in literature [21–23, 30]. Accordingly, it is of practical interest to investigate the performance of LEM in slope probabilistic slope stability analysis especially for the case considering the spatial variability of soil properties.
The performance of LEM (e.g., Bishop, M-P, and Spencer methods) in the determination of FS and FM was compared against that of FEM- and FDM-based SRT for spatially variable undrained soil slope. The spatially variable soil properties were generated using one-dimensional lognormal random field in the vertical direction. The comparative investigations have demonstrated that the Bishop method performs as satisfactorily as FEM- and FDM-based SRT both in FS and FM for the cases where no preferable surfaces (e.g., the existence of inclined soft band and hydraulic tension crack) dominate the failure mechanism. Although M-P and Spencer methods are rigorous approaches satisfying both force and moment equilibriums, they cannot provide a comparable FS and an equivalent FM as compared with FEM- and FDM-based SRT for spatially variable undrained soil slope due to the biased calculation of FS ignoring the variation of interslice functions. Therefore, for the reliability analysis of spatially variable undrained soil slopes, Bishop method can be used as a preliminary approach for the evaluation of failure probability. However, it is worthwhile to note that the FM obtained by the Bishop method tends to be unconservative from the aspect of failure consequence. A precautionary note is given that rigorous LEM (e.g., M-P and Spencer methods) may be carefully employed in the stability analysis of spatially variable soil slopes if the local minima and failure to converge issues are fully addressed. It is highly recommended that FEM- and/or FDM-based SRT be adopted for the reliability analysis of spatially variable undrained soil slope.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant no. 51778313). The financial support is gratefully acknowledged.