#### Abstract

This paper aims to develop an effective tool for quantifying the risk of slope failure and identifying the sources of failure risk by combining the limit equilibrium method and the assumption of maximum area of sliding mass with factor of safety = 1. The assumption adopted in this study is firstly validated through the results from the homogeneous slope model, the laboratory experiment, and the smoothed particle hydrodynamics (SPH) program, respectively. Secondly, the proposed method is implemented through the quantification of slope failure risk and the identification of failure sources for a homogeneous slope and a cohesive slope. The conventional method which quantifies the failure risk based on the slip surface with minimum factor of safety (FS) is also performed to enable the comparison with the proposed method. The comparative study has demonstrated that the conventional method tends to underestimate the failure risk due to the negligence of the whole failure process as compared with the proposed method. The failure risk has a tendency to increase as vertical spatial variability of friction angle and *S*_{u} grow less significant for both proposed method and conventional method. However, the failure sources identified by the conventional method are more likely to decrease as the vertical spatial variability of *S*_{u} becomes less significant for cohesive slope, whereas the proposed method is able to find a nearly constant number of failure sources by considering the whole process of slope failure. As a result, it is worthwhile to point out that attention is highly recommended to be focused on the failure sources when the spatial variability is less significant, even if it is not considered during the risk mitigation and reinforcing works.

#### 1. Introduction

Risk assessment of slope failure accounting for both failure probability and the consequences plays a vital role in slope design and risk mitigation [1, 2]. Much attention has been focused on the issue of quantitative risk assessment of slope failure [3–11]. The failure probability and the consequence were evaluated based on the critical slip surface with the minimum factor of safety (FS) in most of the previous research [12–29]. Although the area (volume for 3D case) of sliding mass corresponding to the critical slip surface with the minimum FS can serve as an index to quantify the consequence of a slope failure, current studies have shown that the critical slip surface with minimum FS is only an initial location of a slope failure, and it cannot represent the whole process of a slope failure, including not only the initiation stage but also the propagation and the evolution stages [8]. Therefore, the use of critical slip surface to determine the consequence is a simple tool, and further studies should be conducted. To properly address this issue, numerical simulations of slope failure using appropriate algorithms have been conducted using material point method (MPM) [30], smoothed particle hydrodynamics (SPH) [31–34], and discrete element method (DEM) [35].

Although the numerical simulations can model the whole process of a slope failure and determine the subsequent consequence, the demanding computational effort involved makes it an unacceptable approach, especially when the Monte Carlo simulation (MCS) is adopted to perform the risk assessment of slope failure where a large number of numerical simulations are inevitable.

To achieve a balance between the computational efficiency and accuracy of calculated results, an alternative tool combining the simplicity of limit equilibrium method (LEM) and the assumption of maximum area of sliding mass (MASM) for slope failure is developed in the current study. As will be demonstrated in the later sections, the proposed method is an effective and efficient tool for the risk assessment of slope failure. The paper starts with the validation of MASM assumption based on respective experimental and homogeneous slope models, followed by the brief review of framework for risk assessment of slope failure by MCS. Then, the proposed methodology is illustrated through three slope stability examples. Finally, conclusions are drawn, and discussions are made to provide insights into geotechnical reliability analysis.

#### 2. Validation of Assumption of MASM and Its Use in Determination of Consequence

##### 2.1. Validation Based on Homogeneous Slope Model

A simplified homogeneous slope model with slope height *H* = 5 m and slope angle *α* = 45° is adopted to validate the assumption of MASM. The soil of slope has a cohesion *c* = 0 and friction angle *φ* = 30°. The critical slip surface with the minimum FS for this slope model is determined by SLOPE/W (LEM) and FLAC (FDM), respectively. The minimum respective FS values calculated by LEM (based on Morgenstern–Price method) and FDM are 0.58, and it is identical to that calculated by using equation (1), that is, FS = tan 30°/tan 45° = 0.58:where is the factor of safety, is the friction angle, and is the slope angle of the simplified slope model. Figure 1 plots the critical slip surface obtained by LEM and the contour of shear strain rate obtained by FDM. It can be noticed from Figure 1 that the critical slip surface by LEM is nearly parallel to the slope surface line. The area of sliding mass encompassed by the critical slip surface and the slope surface line is 0.34 m^{2}. Although the same FS is obtained, FDM indicates a different critical slip surface which is located at the upper part of the slope surface compared with that by LEM, as shown in Figure 1(b). The difference in the critical slip surface may be attributed to inherent assumptions between LEM and FDM, and detailed explanations are not provided herein. As noted in [33], the critical slip surface indicates only the initial location of slope failure and provides limited information for the whole process of slope failure. It can be intuitively appreciated that the slope is still instable after the removal of the critical slip surface obtained by LEM (i.e., the extreme case when the 0.34 m^{2} sliding mass slides away). A new critical slip surface will be obtained based on the updated slope geometry. This process will be repeated several times until the updated geometry has a slope angle of 30° (at that time, the FS calculated using equation (1) is equal to 1.0). The simplified final slip surface with FS = 1 as shown in Figure 1 by dashed blue line has an area of sliding mass = 9.15 m^{2}. Significant difference in the area of sliding mass between the critical slip surface and the simplified final slip surface is observed, and the choice of the slip surface must be properly dealt with before the consequence of a slope failure is rationally quantified.

**(a)**

**(b)**

**(c)**

To further validate the simplified final slip surface, the whole failure process of the simplified slope model is simulated using an in-house FORTRAN-based SPH program which has been validated and used to identify the failure slip surfaces for landslide risk assessment [33]. The fundamentals of SPH are not presented herein, and readers can refer to the literatures including but not limited to Lucy [36]; Gingold and Monaghan [37]; Monaghan, and Lattanzio [38]; Liu and Liu [39]; and Bui et al. [32, 40]. The sliding process of the simplified slope model is simulated using the SPH algorithm, and the slope deposited surface at the steady state (i.e., the landslide is stopped) is plotted in Figure 1(c). To facilitate the comparison of failure slip surface between the SPH and analytical method, the original slope surface is sketched by red line and the simplified final slip surface is shown by dashed blue line in Figure 1(c). The simplified final slip surface predicts fairly well with the boundary line between sliding particles group and stable particles group, especially at the upper part of the slope. At the lower part, the actual boundary line between sliding particles group and stable particles group diverges from the simplified final slip surface. The reason for this may be that the previous failed masses deposit at the slope toe contributing to the stability of the residual slope, which has not been considered in the determination of the simplified final slip surface. Finally, the assumption of MASM is validated, and it is more likely to be used in the landslide risk assessment, especially in the determination of area of sliding mass and further in the quantification of failure consequence.

##### 2.2. Validation Based on Experimental Model Results

To verify the assumption of MASM of slope failure, the experimental model results from Jia et al. [41] were adopted and revisited. A large scale slope model which has a length of 15 m, height of 6 m, and width of 5 m was prepared and adopted to simulate the performance of a sandy silt soil slope subjected to fluctuations in the water level [41]. After the construction of the slope model, the cohesion and internal friction angle of the sandy silt are 1 kPa and 30°, respectively, based on the results from isotropic consolidated, undrained triaxial compression tests. After the slope model was submerged for 72 h, slope geometry changed from an initial slope angle of 45° to an intermediate slope angle of 33°. Figure 2 shows the intermediate slope geometry after 72 h submersion. The matric suction of the sandy silt has decreased to zero after 72 h submersion and therefore resulted in a reduction in the friction angle [41]. The previous research outputs have shown that the sandy silt exhibits a reduction in the friction angle by 10%, and the cohesion value tends to be zero [42]. Hence, the reduced shear strength parameters of sandy silt are 0 kPa and 27°, respectively. The drawdown condition has been conducted by Jia et al. [41], and three failure slip surfaces are observed, as shown in Figure 2. The FS values for these three slip surfaces are recalculated by the Spencer method based on the reduced shear strength parameters. The FS for S1 is 0.89, for S2 is 0.95, and for S3 is 1.01. The calculations indicate that the slip surface with minimum FS is the initiating sliding surface and differentiates from the true sliding surface (refer to S3). The real sliding surface in a landslide has an approximate FS value of 1.0. Therefore, the proposed new method for locating failure slip surface with FS ≈ 1 instead of minimum FS has been validated using experimental results.

##### 2.3. Determination of Consequence

The assumption of MASM is validated through respective experimental and analytical results in the previous sections, and it is used in the determination of consequence of a slope failure. The main steps are demonstrated in Figure 3. As shown in Figure 3, after the soil parameters are determined, the conventional LEM methods such as Bishop method and Morgenstern–Price method can be adopted to calculate the FS value for a priori slip surfaces. If the FS for a given slip surface is lower than 1.0, this slip surface is classified into failure slip surface group. By calculating the FS values for a sufficiently large number of slip surfaces, the failure slip surface group is determined. The consequence is determined based on that failure slip surface with FS ≈ 1 and having MASM within the failure slip surface group. Assume *M* out of *N* slip surfaces are within failure slip surface group and their respective area of sliding mass values and FS are denoted by *A*_{n1}, *A*_{n2}, …, *A*_{n,M−1}, *A*_{n,M}, and FS_{n1}, FS_{n2}, FS_{n3}, and FS_{n,M}. Assume that *Q* out of *M* FSs are approaching 1 and their respective area of sliding mass are denoted by *A*_{m1}, *A*_{m2}, *A*_{m3}, …, and *A*_{m,Q}. The maximum one among *Q* area of sliding mass values is denoted by *A*_{max} and is used to quantify the failure consequence. For a given set of input parameters, the FS for any of *N* priori slip surfaces is larger than 1.0, indicating that no slope failure occurs (i.e., *A*_{max} = 0). The proposed method for calculating the consequence of a slope failure is used in the risk assessment of slope failure by MCS in the next section.

##### 2.4. Revisit of the Homogeneous Slope

The proposed method is used to locate the failure consequence for homogeneous slope discussed in Section 2.1. A large number of priori circular slip surfaces based on (*N* = 20000 herein) are firstly generated to cover the potential failure domain. The circular instead of noncircular slip surface is selected considering the convergence issue of FS for noncircular slip surface using the Morgenstern–Price method [43–45]. The Bishop method is used to calculate the FS for a given priori slip surface, considering that the factor of safety from the Bishop method is usually close to that from the Morgenstern–Price method. Finally, the failure slip surface group is determined and is composed of *M* = 8209 priori slip surfaces. If the FS for a slip surface within the failure slip surface group is larger than 0.99, the area of sliding mass corresponding to this slip surface is calculated. The FSs for a total of *Q* = 235 slip surfaces are larger than 0.99, and therefore, their respective areas of sliding mass are compared, and the maximum one is adopted to quantify the failure consequence. The slip surface corresponding to the maximum area of sliding mass is plotted in Figure 1(c) by a green line. The boundary line between stable particles and sliding particles in the SPH program indicates the slip surface for a slope failure, and it is compared with that determined using the proposed method. It is noticed from Figure 1(c) that the slip surface determined by the proposed method agrees well with that sketched by the SPH program. The proposed method is extended to the risk assessment of slope failure by MCS.

#### 3. Risk Assessment of Slope Failure by MCS

After the failure consequence has been determined using the proposed method described in Section 2.3, the risk assessment of slope failure can be performed by MCS. As shown in Figure 4, before the risk assessment of slope failure can be conducted, the necessary input information for risk assessment of slope failure must be characterized. The required information includes, but not limited to, slope geometry and statistical information of geotechnical properties, such as mean values, standard deviations, marginal distributions, autocorrelation functions, and scale of fluctuations. Secondly, the random field realizations of geotechnical properties are simulated using the covariance matrix decomposition method. Thirdly, each set of random samples is taken as the deterministic inputs into the slope stability analysis model, and the area of sliding mass for the slope failure is determined based on the methodology described in Figure 3, which is adopted to quantify the failure consequence. Finally, the mean value of the area of sliding mass for the total sets of random samples is used to evaluate the risk. The details for the generation of random fields using the covariance matrix decomposition method and the steps following the calculation of FS using the limit equilibrium method when considering the spatial variability of properties can be found by Li et al. [22].

#### 4. Case Studies

##### 4.1. Homogeneous Slope

The homogeneous slope with geometry shown in Figure 1 is analyzed to assess the risk of slope failure. The friction angle of the soil has a mean value *μ* = 45° and standard deviation *σ* = 5° following a lognormal distribution. The scale of fluctuation *λ* and the autocorrelation function are used to consider the vertical variability of the friction angle. As shown in Figure 5, the 7-m-thick soil layer is divided into fourteen 0.5-m-thick sublayers, and *φ* at each sublayer is represented by an entry in a *φ* vector with a length of 14. Let and denote the value of at location *i* and *j*, respectively; *d*_{ij} denote the vertical distance between locations *i* and *j*. The simple exponential correlation function is adopted herein to simulate the correlation coefficient *ρ*_{ij} between and [22, 28]:

The covariance matrix decomposition method [22, 28, 46, 47] is used to generate the realizations of the vertical random field of the friction angle. The area of sliding mass for each realization (i.e., a *φ* vector with length = 14) can be determined using the proposed methodology, and the final risk value can be obtained following the steps described in Figure 4. Different values for *λ* (=2, 4, 5, 10, 20, 100, and 1000 m) are assumed to investigate the effect of *λ* on risk of slope failure, *R*. A large number (*N* = 20000) of priori circular slip surfaces used in Section 2.4 are still adopted to cover the potential failure domain. The number of samples in MCS, *T*, is equal to 1000.

###### 4.1.1. Effect of *λ* on Failure Risk for Homogeneous Slope

Table 1 summarizes the simulation results for different *λ* values using the proposed method. It is noted from Table 1 that, as *λ* increases from 2 m to 1000 m, the failure probability increases from 0.22 to 0.35, and this trend has been validated by many researchers [22–24, 28]. The risk of the slope failure represented by the area of sliding mass is calculated by the proposed method and the conventional method based on the slip surface with minimum FS, respectively. The computed results are shown in Figure 6. It is noted from Figure 6 that, as *λ* increases from 2 m to 100 m, the risk value, *R*, increases from 1.17 to 2.3 m^{2} for the proposed method, while the risk value, *R*, calculated by the conventional method increases from 0.74 to 1.13 m^{2}. The comparison of the simulation results has shown that the risk value determined by the proposed method is larger than that determined by the convention method at the same *λ* value. This is due to the consideration of the whole failure process in the proposed method, whereas only the failure initiation is considered in the conventional method.

###### 4.1.2. Risk Deaggregation of Slope Failure Risk for Homogeneous Slope

The areas of sliding mass with different values indicate various sources of failure consequence, and in order to further deaggregate the risk to the separate sources, the histogram of area of sliding mass is plotted in Figures 7 and 8, respectively, for the conventional and proposed method. The area of sliding mass for a given set of MCS samples is classified into fourteen area intervals increasing from [2, 3] to [14, 15] with increment of 1. By counting the number of MCS samples within specified interval, the frequency is obtained and is plotted in *y*-axis. It is noticed from Figure 7 that the area of sliding mass for 225 samples is within the interval [3, 4], and only one area of sliding mass is within the interval [4, 5] when the vertical scale of fluctuation *λ* is 2 m. As *λ* increases from 2 m to 1000 m, the frequency for the interval [3, 4] increases from 225 to 347. The comparison of the results indicates that the area of sliding mass related with the slope failure at different *λ* values remains unchanged although the frequency of slope failure increases as the spatial variability grows more insignificant (i.e., *λ* has a larger value). The conventional method deaggregates the failure risk to only one source, as shown in Figure 7. The one source means that the slope failure occurs along one same slip surface, no matter how the friction angle of soil changes.

As for the proposed method, the failure risk is attributed to a finite number of sources rather than only one source owing to the consideration of whole process of slope failure. When *λ* is equal to 2 m, the respective frequency corresponding to area of sliding mass interval [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 15], and [15, 16] is 0, 43, 72, 61, 26, 19, 4, 0, 1, 0, 0, 0, 0, and 0. When it is compared with the results from the conventional method, the failure risk is deaggregated to 7 sources by the proposed method. Besides the area of sliding mass interval [3, 4], the slope failure may occur leading to the area of sliding mass within other intervals such as [4, 5], [5, 6], [6, 7], and [7, 8]. Accordingly, the risk obtained by the proposed method is 1.17 m^{2} compared with the risk of 0.74 m^{2} calculated by the conventional method. It can be concluded that the conventional method tends to underestimate the failure risk compared with the proposed method. When *λ* increases to 4 m, the respective frequency corresponding to the same area of sliding mass intervals as in *λ* = 2 m is 0, 40, 60, 65, 37, 28, 15, 6, 2, 0, 1, 0, 0, and 0. The respective number of sources at *λ* = 4 m, 5 m, 10 m, 20 m, 100 m, and 1000 m is 9, 10, 11, 12, 12, and 12. It can be seen that, as the vertical scale of fluctuation *λ* increases, the number of sources contributing to the failure risk varies slightly about 10. Unlike the conventional method where only one source is deaggregated, the proposed method provides an effective tool to quantify the failure risk and deaggregate the risk to separate sources for landslide risk mitigation treatment. The proposed method is extended to the risk assessment of slope failure for cohesive slopes in the next subsection.

##### 4.2. Cohesive Slope

The proposed methodology is extended to assess the failure risk for a cohesive soil slope example that has been analyzed by Wang et al. [28] to study the effect of spatial variability of soil properties and uncertainty in the critical slip surface. As shown in Figure 9, the cohesive soil slope has a height *H* = 10 m and slope angle of 26.6°, corresponding to an inclination ratio of 1 : 2. The cohesive soil is underlain by a firm stratum at 20 m below top of the slope. Short-term shear strength of the cohesive soil is characterized by undrained shear strength *S*_{u}, and the saturated unit weight of soil is *γ*_{sat}. *γ*_{sat} is taken as a deterministic value of 20 kN/m^{3}, and the *S*_{u} is modeled by a one-dimensional random field spatially varying along the vertical direction. The value of at the same depth is assumed to be fully correlated. The spatial variability with depth is modeled by a homogeneous lognormal random field with an exponentially decaying correlation structure. As shown in Figure 9, the 20-m-thick cohesive soil layer is divided into forty 0.5-m-thick sublayers, and at each sublayer is represented by an entry in a *S*_{u} vector with a length of 40. The mean and standard deviation of are approximately equal to 40 kPa and 10 kPa (i.e., 25% coefficient of variation (COV)), respectively. Different values for *λ* (=2, 4, 5, 10, 20, 100, and 1000 m) are assumed to investigate the effect of *λ* on *R*. A large number (*N* = 20000) of priori circular slip surfaces are adopted to cover the potential failure domain. The number of samples in MCS, *T*, is equal to 1000.

###### 4.2.1. Effect of *λ* on Failure Risk for Cohesive Slope

The failure probability calculated by the proposed method at different *λ* values is summarized in Table 2. It can be seen from Table 2 that, as *λ* increases from 2 m to 5 m, the failure probability increases from 0.18 to 0.28, and the failure probability varies slightly about the value of 0.3 after *λ* is larger than 5 m. This variation trend is similar to that observed for homogeneous slope, and it also agrees well with that reported by Wang et al. [28]. The risk of the slope failure represented by the area of sliding mass is calculated by the proposed method and the conventional method. Figure 10 summarizes the results in a plot of risk values versus the vertical scale of fluctuation *λ*. For the proposed method, as *λ* increases from 2 to 10 m, the value of *R* increases significantly from about 116 to about 260 m^{2}. When *λ* is larger than 10 m, the effect of *λ* on *R* begins to diminish, and *R* varies slightly as *λ* further increases. Similar trend has been observed for the conventional method. Owing to the negligence of the whole process for slope failure, the conventional method leads to a smaller risk value compared with the proposed method at the same different *λ* values. It is concluded that the conventional method has a tendency to underestimate the failure risk.

###### 4.2.2. Risk Deaggregation of Slope Failure Risk for Cohesive Slope

The potential range of area of sliding mass is equally divided into 27 intervals growing from [100, 150] m^{2} to [1400, 1450] m^{2} with 50 increments. If the obtained area of sliding mass for a specified set of MCS samples is within one interval, the frequency of this interval will increase by one. 1000 area of sliding mass values for *T* = 1000 sets of MCS samples can be obtained both by the proposed method and the conventional method, and the frequency of each interval is counted. By doing this, the failure risk is deaggregated into separate sources, which can provide much more insight into risk mitigation and reinforcing works. Figure 11 compares the frequency corresponding to each area interval between the proposed method and the conventional method under different *λ* values. Figure 11(a) summarizes the histogram of area of sliding mass for the case of *λ* = 2 m, where the vertical spatial variability of *S*_{u} is significant. The failure probability of 0.18 indicates that 180 out of 1000 MCS samples leads to slope failure. Therefore, there are 180 areas of sliding mass values, and they are classified into the 27 intervals. If the frequency corresponding to one interval is zero, this interval is called zero interval; otherwise, it is called nonzero interval. The number of nonzero intervals for the conventional method and the proposed method are 12 and 17, respectively. For the conventional method, the most frequent source is within the area of sliding mass interval [650, 700], and the next seven equal frequent sources are [300, 350], [500, 550], [400, 450], [550, 600], [450, 500], [350, 400], and [600, 650]. Finally, there are four least frequent sources at [150, 200], [200, 250], [250, 300], and [700, 750]. For the proposed method, there is no most frequent source compared with the conventional method. The eight dominant sources are within intervals of [450, 500], [500, 550], [550, 600], [600, 650], [650, 700], [700, 750], [750, 800], and [800, 850]. The nine less-dominant sources are within intervals of [200, 250], [250, 300], [300, 350], [350, 400], [400, 450], [850, 900], [900, 950], [950, 1000], and [1050, 1100]. The details in Figures 11(b)–11(f) are not described. As the vertical spatial variability of *S*_{u} grows less significant (i.e., *λ* increases from 4 m to 1000 m), the number of nonzero intervals varies from 13, 14, 12, 12, 10, and 4 for the conventional method and varies from 21, 20, 23, 25, 23, and 15 for the proposed method, as can be seen from Figure 12. It is inferred from Figure 12 that, as the vertical spatial variability of *S*_{u} becomes less significant, the number of nonzero intervals tends to decreases for the results obtained by the conventional method. This trend agrees well with that reported by Wang et al. [28]. Since the area of sliding mass in the conventional method represents the failure initiation, the observed trend indicates that the initiation source of slope failure tends to widen as the vertical spatial variability of *S*_{u} grows significant. Whereas for the proposed method, the number of nonzero intervals varies slightly about a fixed number (i.e., 20 in the cohesive slope model), and no significant decrease tendency has been found. Furthermore, as *λ* value increases from 2 m to 1000 m, the sources of failure risk move from the left part towards the right part of the *x*-axis. The comparative study has shown that the mere identification of failure initiation can not fully address the issue of risk deaggregation using the conventional method. The proposed method provides an efficient tool to identify the full sources of failure risk using the assumption of maximum area of sliding mass and to rationally quantify the failure risk.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

#### 5. Discussion and Conclusions

The area of sliding mass of slope failure must be properly determined before the corresponding failure consequence can be quantified in risk assessment of slope failure. The area of sliding mass calculated using the slip surface with minimum FS tends to underestimate the “true” area of sliding mass for the slope failure, because the slip surface with minimum FS is only the initial slip surface for a slope failure without considering the propagation and evolution of the slope failure process. By combining the assumption of maximum area of sliding mass and the limit equilibrium method, an effective approach is developed in this study, and it is used to quantify the slope failure risk and to deaggregate the failure risk to separate sources for a homogeneous slope and a cohesive slope. The failure probability and the failure risk have been calculated using the proposed method and the conventional method which bases the failure risk on the slip surface with minimum FS. The results have demonstrated that the conventional method tends to underestimate the slope failure risk as compared with the proposed method. The number of failure sources decreases as the vertical spatial variability of *S*_{u} becomes less significant, whereas it remains almost unchanged as vertical spatial variability of friction angle varies when the conventional method is used. However, the number of failure sources varies slightly about a specific value as the vertical spatial variability of *S*_{u} and/or friction angle changes when the proposed method is used. Therefore, it is worthwhile to point out that same attention should be paid to focus on the failure sources when the spatial variability is less significant, even if it is not considered during the risk mitigation and reinforcing works.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by grants from the National Natural Science Foundation of China (Grant no. 51778313) and Cooperative Innovation Center of Engineering Construction and Safety in Shandong Blue Economic Zone. The financial supports are gratefully acknowledged.