The first method for identifying RSSs was proposed by Zhang et al. (2011) based on the observation that FSs of many PSSs are somewhat correlated. With this idea, the contribution of different PSSs to the system failure probability of a slope can be different and represented by some important slip surfaces or RSSs. For simplicity purposes, the major steps of the method are summarized as follows:

Step 1

Calculate the reliability index β of each PSS using FORM with the following equation: β = min g θ − 1 = 0 α α T , (1) where g θ − 1 = 0 is the limit state function for a slip surface, θ is a vector of random variables considered, and α is a vector of uncorrelated reduced variables corresponding to θ .

Step 2

Find the slip surface corresponding to the smallest reliability index and take it as an RSS.

Step 3

Calculate the correlation coefficients between the FS of the RSS found in Step 2 and the FSs of other PSSs using the following equation: ρ i j = α α T β i β j , (2) where β i and β j are reliability indices for the ith and jth PSSs calculated in Step 1.

Step 4

Exclude the slip surfaces that have correlation coefficients with the RSS found in Step 2 larger than a threshold value ρ 0 from further consideration, as these PSSs are represented by the RSS identified in Step 2.

Step 5

Repeat Steps 2–4 until all PSSs are excluded and represented.

It should be noted that ρ 0 has a significant influence on the results of RSSs. A larger ρ 0 results in a larger number of RSSs, leading to more cost of computer resources, and vice versa. In the original work by Zhang et al. (2011), ρ 0 is determined by a parametric study where several values of ρ 0 are chosen to study its influence on the P f , although generally a value around 0.8 can reach accurate estimations of the system P f for the considered slope examples. In addition, the ISV of soil properties is not considered in this method.

This method is proposed based on Method I by Li et al. (2013) to overcome the problem of Method I in recursively selecting RSSs. It also uses the reliability index and correlation coefficient between two PSSs to identify the RSSs. The differences between the two methods lie in that Method II identifies all RSSs at one time, and a reference slip surface is used to calculate the correlation coefficients. The concrete procedures follow:

Step 1

Perform a deterministic stability analysis of a slope to obtain the critical deterministic slip surface (CDSS) and take it as a reference slip surface S 0 .

Step 2

Calculate the correlation coefficients between the FS of S 0 and all PSSs using the following equation (Chowdhury and Xu, 1995; Bhattacharya et al., 2003): ρ k l = ∑ j = 1 m h j k h j l σ x j 2 ∑ j = 1 m h j k 2 σ x j 2 ∑ j = 1 m h j l 2 σ x j 2 h j k = ∂ G k ∂ x j   h j l = ∂ G l ∂ x j ∂ G ∂ x j = G + − G − 2 σ x j , (3) where σ x j is the standard deviation of random variable x j , G k and G l are, respectively, the limit state function values corresponding to the kth and lth slip surface, G + and G − are limit state function values for the variable x j greater than and less than the mean value by σ x j , respectively, and ρ k l is the correlation coefficient between the FS of the k th and l th slip surface.

Step 3

Evaluate the equivalent reliability index β of each PSS using the following equation ( Ang and Tang, 2007 ): β k = G k μ X ∑ i = 1 m ∂ G k ∂ x i 2 σ x j 2 + ∑ i = 1 m ∑ j ≠ i m ∂ G k ∂ x i ∂ G k ∂ x j χ i j σ x i σ x j , (4) where the subscript k indicates the calculation is based on the kth slip surface, χ i j is the correlation coefficient between random variables x i and x j , and G μ X is the limit state function value when all random variables are set as their mean values.

Step 4

Sort the correlation coefficients obtained in Step 1 in decreasing order and then divide the slip surfaces corresponding to them into N w (e.g., 10) groups.

Step 5

Find the slip surface with the minimum equivalent reliability index in each group of slip surfaces and take that slip surface as the RSS of each group.

It is worthwhile to point out that Li et al. (2013) named the RSSs identified in Step 5 as candidate RSSs, and the final RSSs of slope are selected by the number of failures of each candidate RSS within an MCS analysis. The candidate RSSs, which have no contributions to the probability of slope failure, are not included in the final RSSs. In other words, slope reliability analysis can only be performed based on the candidate RSSs because the RSSs are the byproduct of reliability analysis. Therefore, to have a consistent comparison with other methods, herein the candidate RSSs identified in Step 5 are referred as the RSSs of the method. It is also noted that although the ISV of soil properties is properly considered in this method, the computation cost of modeling the ISV might be very high in the case of thousands of random field elements.

This method is conceptually simple. It is proposed by Jiang et al. (2015) based on the observation that the CSSs for different realizations of random variables or random fields might be the same PSS. Consider, for example, that an MCS with N p random field realizations would result in N p CSSs. However, most of the N p CSSs would be the same PSS (e.g., the CDSS obtained by considering soil properties as spatial constants), which finally leads to N r ( N p ≪ N r ) CSSs. Jiang et al. (2015) showed that reliability analysis performed based on N r CSSs produces very similar results as that performed on all PSSs. The N r CSSs, therefore, are considered as the RSSs and can be used for the subsequent slope reliability analysis. In general, the identification of RSSs of the method consists of the following three steps:

Step 1

Generate N p (e.g., about 1,000) realizations of random fields according to prescribed statistical information.

Step 2

Perform deterministic slope stability analysis with the N p realizations of random fields to obtain the corresponding N p CSSs.

Step 3

Find the duplicates of the N p CSSs and take the unique CSSs as the RSSs.

Compared with the aforementioned two methods, this method does not need to calculate the reliability index for each PSS and correlation coefficient between the FSs of the PSSs but requires some effort in performing MCS with a small number of MCS samples before conducting reliability analysis. Thus, the method is relatively simple and can easily deal with reliability problems with anisotropic ISV of soil properties. However, it is found that the RSSs identified by the method are very sensitive to the soil statistics, such as the auto correlation length of ISV. This is thus not convenient for slope reliability analysis involving many parametric studies that are often the case in practice.

Because the correlation coefficients between the FSs of two arbitrary PSSs in Methods II and III are calculated by an approximate approach that cannot consider the ISV of soil properties and by FORM that is computationally expensive, respectively, Jiang et al. (2017) proposed using the Pearson correlation coefficient (Benesty et al., 2009) as an improvement to measure the correlation between the FSs of PSSs. With the Pearson correlation coefficient, Jiang et al. (2017) proposed a method similar to Method I for identifying the RSSs of a slope that also involves dividing the correlation coefficients into different groups. However, the candidate RSS in each group of the method herein is taken as the slip surface with the smallest FS, rather than the smallest reliability index, among the PSSs in each group, which is computationally simpler and more efficient. For completeness, the major steps of the method are briefly described as follows:

Step 1

Perform deterministic slope stability analysis with N PSSs and select the CDSS as the first RSS.

Step 2

Calculate the correlation coefficients between the FSs of the CDSS and the other PSSs using the Pearson correlation coefficient that is estimated by the simulation method (Zheng et al., 2016): ρ 1 , i = ∑ j = 1 N s i m F S 1 X j − F S i ¯ F S i X j − F S i ¯ ∑ j = 1 N s i m F S 1 X j − F S 1 ¯ 2 ∑ j = 1 N s i m F S i X j − F S i ¯ 2 , (5) where N s i m is the number of MCS samples for the underlying random variables and/or fields, F S i X j is the FS of the ith slip surface for the jth samples of random variables and/or fields X j , and F S i ¯ is the mean value of F S i .

Step 3

Sort the N p correlation coefficients in ascending order and divide them into M groups based on the increment of Δ ρ = ρ max − ρ min / M , which indicates that the corresponding PSSs are sorted into M groups.

Step 4

Select the slip surface with the minimum FS within each group as the RSS for each group, resulting in N p RSSs.

Note that, similar to Method I, the RSSs identified in Step 4 are named as candidate RSSs in the work by Jiang et al. (2017), and the final RSSs of slope are selected based on the contributions of each candidate RSS to the system failure probability. However, slope reliability analysis is performed based on the candidate RSSs because the RSSs are the byproduct of reliability analysis. Again, to have a consistent comparison with other methods, herein the candidate RSSs identified in Step 4 are referred as the RSSs of the method. It is also noted that although the ISV of soil properties is properly considered in this method, the computation cost of modeling the ISV might be high in the case of thousands of random field elements. In addition, both the simulation times for evaluating correlation coefficients and the subdividing groups of the correlation coefficients need to be determined by a preliminary parametric study.

Although it is not necessary to search for the FS among all PSSs but only among the RSSs, direct calculation of the FS using LEM thousands of times is still not a trivial task. To improve the computation efficiency, a multiple response surface method that has been demonstrated to be effective and efficient is adopted here (Li et al., 2015; Li D.-Q. et al., 2016). The quadratic polynomial without cross terms is taken to construct the RSM as F S i X = a 1 i + ∑ j = 1 n b i j x j + ∑ j = 1 n c i j x j 2 , i = 1 , 2 , . . , N r , (7) where F S i X is the FS for the i t h RSS, N r is the number of RSSs, X is the random variable vector, and a 1 i , b i j , and c i j are unknown coefficients that can be calibrated using central composite design method.

Based on the established RSMs, MCS is then performed for reliability analysis. Given N s i m random field and/or variable samples, the P f of slope failure is calculated as P f = 1 N s i m ∑ j = 1 N s i m min i = 1,2 , ⋯ , N r F S i X j < 1 , (8) where I ⋅ is an indicator function that is equal to unity when min i = 1,2 , ⋯ , N r F S i X j < 1 and zero otherwise. The accuracy of P f is assessed by its coefficient of variation (COV) as C O V P f = 1 − P f P f N s i m . (9)

In this subsection, the AKCBR is applied to a layered slope without considering ISV. The slope has been widely studied in the literature (Chowdhury and Xu, 1995; Zhang et al., 2011). Figure 1 shows the geometry of the slope, which is a fill embankment resting on a clay layer. The statistics of soil properties for the two soil layers are tabulated in Table 1 and are the same as those used by Zhang et al. (2011). The shear strength parameters are subjected to normal distributions, whereas the unit weights are considered constants. The mean values of the cohesion c 1 and friction angle φ 1 for the embankment layer and are 10 kPa and 12°, respectively, whereas the mean value of the undrained strength of the clay layer c u 2 is 40 kPa. The COVs of c 1 and φ 1 of the fill embankment are, respectively, 0.2 and 0.25, whereas the COV of the c u 2 of the clay foundation is 0.2. The unit weights of the fill embankment and clay foundation are 18 kN/m³ and 20 kN/m³, respectively.

With the aforementioned mean parameters, the deterministic slope stability analysis is performed using the Bishop’s simplified method. A total of 4,250 PSSs are generated to cover the whole slope, and the FS of the CSS is 1.131, which is close to the value of 1.148 reported by Chowdhury and Xu (1995). Note that the FSs and the sliding volumes for the other PSSs are also calculated during the deterministic analysis, which can be conveniently used for the subsequent identification of RSSs. As mentioned before, a predefined range for the K value is a prerequisite for the adaptive K-means clustering method. From the previous study (Wang et al., 2020), an empirical range of [3, 200] is adopted herein. Then, the DUNN indices for different K values are obtained using the suggested procedure in Section 3 and Figure 2. The DUNN values are plotted against different K values. It can be seen from the figure that the K value has a significant effect on the DUUN index. As K is taken as 108, the DUNN index reaches the maximum, which means that clustering all PSSs into 108 sub-clusters is the best. Thereafter, the slip surface with the smallest FS in each of the 108 clusters is taken as the RSS of that cluster. Overall, 108 RSSs are finally identified and are also plotted in Figure 1. Note that the CSS signified by the dashed line in Figure 1 is also included in the final set of RSSs. To gain more insight into the effect of the clustering on PSSs, Figure 3 plots three typical clusters of the PSSs for the slope example. The figure shows that different clusters of slip surfaces exhibit different failure modes, and the slip surfaces in the same cluster present similar failure modes. For example, the slip surfaces of the 11th cluster plotted in Figure 3A are mainly shallow failure modes, whereas the slip surfaces in the 36th and 78th clusters tend to be much deeper.

Response surfaces are then constructed based on the previously identified RSSs using P f , followed by the MCS directly performed on the response surfaces to evaluate the P f of the slope. The P f is estimated as 0.398 and is very close to the value of 0.404 reported by Zhang et al. (2011), validating the accuracy of the method. The efficiency of the AKCBR can be approximately measured by a dimensionless nominal index N s s that indicates the total times of the FS evaluation of a single slip surface. Consider, for example, the deterministic slope stability analysis of the slope using 4,250 PSSs. The N s s is equal to 4,250 because 4,250 evaluations of the FS are required for the 4,250 PSSs. The computation cost for the reliability analysis of the AKCBR is mainly composed of 1) the time for identification of the RSSs and 2) the time for the calibration of the RSS-RSM. The first part can be equivalently evaluated using the actual physical time for the RSS identification divided by the time needed to calculate the FS of a single slip surface. For example, it takes about 0.001 s to calculate the FS of a single slip surface using a desktop computer with 16G RAM and an Intel(R) Core(TM) i9-9900x and 59 s to select the RSSs, so the N s s for the first part is approximately calculated as 59/0.001=59,000. In addition, the time of constructing the RSM and MCS with RSM is trivial and can be ignored. Therefore, the N s s for the AKCBR is about 63,250. In contrast, the time consumption of Method I (Zhang et al., 2011) is mainly for the evaluation of the reliability index of the 4,250 PSSs based on Eq. 1, and it takes about 80 s with the same computer. As such, the N s s for Method I is about 80,000, which is larger than that for the AKCBR.

This part applies the AKCBR for evaluating the reliability of an undrained cohesive slope to consider the ISV of the undrained shear strength c u in the vertical direction, which has also been studied in the literature (Wang et al., 2011; Li et al., 2013). The geometry of the slope is plotted in Figure 4, which has a slope height of 10 m and a slope angle of 26.6°. To enable a convenient and consistent comparison with previous studies, the mean of c u is 40 kPa, and the coefficient of variation is 0.25 (Wang et al., 2011; Li et al., 2013). The saturated unit weight γ s a t of the soil is 20 kN/m³, which is considered as a deterministic value in this analysis. The undrained shear strength c u , however, varies spatially in the vertical direction and has a mean μ c u and C O V c u of 40 kPa and 0.25, respectively. The ISV of c u is simulated by a stationary random field characterized by 40 random variables, which are the same as those used by Li et al. (2013) and described by the single exponential decaying correlation structure as R i j = exp ⁡ ⁡ − z i − z j 2 λ , (10) where R i j denotes the correlation between the soil properties at the depths z i and z j , and λ is the autocorrelation length. In theoretical engineering, λ can be generally determined from a large quantity of measurement data, which is, however, often not available in practice. Therefore, different values of λ varying from 0.5 to + ∞ are adopted to consider the effect of SOF on slope reliability assessment.

First, the deterministic slope stability analysis model was established using the aforementioned mean parameters, with 4,281 PSS covering the whole slope. The slope stability results are schematically shown in Figure 4, where the CSS passes through the bottom of the slope. The corresponding FS is 1.178, which is identical to the FS reported by Wang et al. (2011). Then, the 4,281 slip surfaces are divided into different clusters using the proposed adaptive K-means clustering method. The results show that when K is set as 154, the DUNN index reaches the maximum, signifying the optimal cluster number for the PSSs. Therefore, there are finally 154 RSSs identified for this slope example. All the RSSs are plotted in Figure 4, and the deterministic CSS is also included in the RSSs.

To illustrate the accuracy of the AKCBR, various parametric studies on the autocorrelation distances normalized by the slope height are conducted. As can be seen from Figure 5, the P f generally increases as the normalized autocorrelation distance increases. When λ H increases from 0.05 to 1 (or λ increases from 0.5 to 10 m for H=10 m), the P f significantly increases from 1.3% to 28%, showing that the ISV has a significant influence on slope reliability. When λ H > 0.5 , the P f becomes insensitive to the change of λ as the spatial variety is nearly ignored. Figure 5 also includes the results of Wang et al. (2011) (based on MCS) and Li et al. (2013) (i.e., Method II), which agree well with the results from the current study and validate the accuracy of AKCBR. It is, however, worth noting that although the number of RSSs identified in this study is larger than that selected by Method II (Li et al., 2013), the proposed method still has eminent efficiency in reliability analysis. To illustrate this point, the dimensionless nominal index N s s is adopted again to compare the computation efficiency between the two methods. A desktop with 16G RAM and the Intel (R) Core(TM) i7-9700 k CPU requires about 0.001 s to evaluate the FS of a single slip surface. The N s s here for the AKCBR is thus approximately calculated as the summation of 1) the time for the deterministic analysis (i.e., 4,281), 2) the time for the identification of the RSSs (i.e., 73/0.001=73,000), and 3) the time for RSM calibration (i.e., (40×2+1)×154=12,474), which is about 89,755. In contrast, most of the computation cost for Method II is consumed by the calculation of the equivalent reliability index for all PSSs, which requires about 81 times the deterministic analysis, thus resulting in a N s s of about 346,761. It is obvious that the computation cost for Method II, even ignoring the subsequent reliability analysis, is higher than AKCBR with full analysis.

This subsection takes a cohesive-frictional slope as an example to further illustrate the applicability and effectiveness of the AKCBR for reliability analysis considering 2-D ISV. The slope has been studied in the literature (Cho, 2010; Jiang et al., 2015; Liu et al., 2017; Liu et al., 2018), and the results from these available studies can be easily referred to. A similar slope study case was performed in Javankhoshdel et al. (2020). The cross-section of the slope is shown in Figure 6A, where the slope height and slope angle are 10 m and 45°, respectively. The total unit weight of the soil is 20 kN/m³. The cohesion c and friction angle φ are considered spatially varied and characterized by two cross-correlated 2-D lognormal stationary random fields. The cross-correlation coefficient between c and φ is −0.7. The mean values of c and φ are, respectively, 10 kPa and 30°. The COV values of c and φ are 0.3 and 0.2, respectively. The horizontal and vertical autocorrelation distances, θ h and θ v , are 20 m and 2 m, respectively. All these parameters are consistent with those used in the literature (Cho, 2010; Jiang et al., 2015; Liu et al., 2017; Liu et al., 2018).

Similar to the last slope example, the geometry of the slope is discretized into 1,210 random field elements with a side length of 0.5 m to consider the ISV of c and φ . The mid-point method is also adopted here to model the cross-correlated random fields of c and φ with the squared exponential autocorrelation function. Then, slope stability analysis is performed to obtain the FS of the slope using Bishop’s simplified method. The FS is calculated as 1.205, and the associated CSS among 7,436 predefined PSSs is plotted in Figure 6A. The results are the same as those reported by Jiang et al. (2015), suggesting the accuracy of the slope stability analysis model.

Then, with the proposed adaptive K-means clustering method, the 7,436 slip surfaces are divided into different clusters. The results show that when K is set as 110, the DUNN index reaches the maximum, signifying the optimal cluster number for the PSSs. Therefore, there are finally 110 RSSs identified for this slope example. All the RSSs are plotted in Figure 6B, and the deterministic CSS is also one of the RSSs. With the RSSs, reliability analysis is performed using the proposed RSS-RSM in Section 4 to validate the effectiveness and efficiency of the adaptive K-means clustering approach for RSS identification. The P f based on adaptive K-means-based RSS-RSM is 4.4 × 10⁻³, and the P f values corresponding to Method III and MCS are 4.9 × 10⁻³ and 3.9 × 10⁻³ (Cho, 2010; Jiang et al., 2015). This similarity validates the accuracy of the AKCBR as well as the applicability for reliability analysis considering 2-D ISV. Note that Wang et al. (2020) previously analyzed the reliability of this example with the same approach but using a different number of RSSs. The reason is two-fold: 1) the number and locations of the PSSs here are different from those used by Wang et al. (2020); 2) the initial clustering centers are determined by an optimized method in MATLAB 2018b; they are randomly selected in MATLAB 2014a. Nevertheless, the reliability results from the two methods are very similar, indicating that the two methods successfully identified the RSSs.

To gain more insight into the AKCBR, the effect of the horizontal and vertical spatial variability on the reliability of the slope is further checked by the AKCBR in order to illustrate its effectiveness and efficiency against the variations of θ h and θ v . The results are plotted in Figure 7. As a reference, the results from Liu et al. (2018) and Li et al. (2015) are also plotted in this figure. Note that, in the figure, only one parameter (i.e., θ h or θ v ) is changed for each case, while others are kept the same as previously defined. In general, the P f increases as the autocorrelation distances increase. The results obtained from this study show a good consistency with those from Liu et al. (2018) and Li et al. (2015). Such agreement suggests that the AKCBR is accurate and robust against the variation of ISV. Note that although Jiang et al. (2015) did not study the effect of ISV on the reliability of this slope example using their method, it is shown by a cohesive slope example application that their method is also robust against the variation of ISV. The AKCBR, therefore, is comparable with Method III in the accuracy of the reliability evaluation for a 2-D spatially varied slope.

Jiang et al. (2015) also conducted a sensitivity study to investigate the effect of cross-correlation on the reliability of the slope. To keep a consistent comparison with Method III, the proposed method is also applied to this example to examine the influence of the cross-correlation coefficient. The results obtained by this study and from Method III are plotted in Figure 8A. It is seen from the figure that the P f obtained by this study increases from 0.44% to 6.8% as ρ c φ varies from −0.7 to 0, suggesting a significant effect on P f . Meanwhile, the results from Method III by Jiang et al. (2015) almost overlap those from the current study by the proposed adaptive K-means clustering-based RSS-RSM method, showing a good agreement between the two methods.

The efficiency of the AKCBR can also be illustrated by the dimensionless index N s s used for the last two examples. The physical time for selecting the RSSs is about 69 s on a desktop computer with 16G RAM and Intel(R) Core(TM) i9-9900x, and 2,421 runs of the slope stability analysis model are required to calibrate the response surfaces of the 110 RSSs. This leads to the N s s of about 4 × 10 5 for the AKCBR. In contrast, it requires about an N s s of 7.4 × 10 6 to select the RSSs by Method III, where 1,000 LHS samples are generated to directly execute the slope stability analysis model with 7,436 PSSs. This comparison thus validates the efficiency of the AKCBR. In addition, note that the RSSs identified by Jiang et al. (2015) are sensitive to the variation of statistics of shear strengths (e.g., ISV and ρ c φ ), which, however, is bypassed by the current study and will be further illustrated in the later discussion section. For example, Figure 8B compares the efficiency of the two methods in identifying the RSSs. It is observed that the AKCBR is independent of the ρ c φ , whereas the number of RSSs identified by Method III shows an increasing trend toward the increase of ρ c φ . This further indicates the high efficiency of the AKCBR.

This part analyzes the reliability of a layered cohesive slope using the AKCBR. The slope was also widely studied in the literature (Li D. Q. et al., 2016; Jiang et al., 2017; Wang et al., 2020) for reliability analysis. The geometry of the slope is schematically shown in Figure 9A, which has a slope height of 24 m and a slope angle of 36.9°; the slope is composed of two clayey soil layers. The statistics of the undrained shear strengths are the same as those used by Jiang et al. (2017). The mean of c u of layer 1 is 80 kPa, and the corresponding COV is 0.3. The mean of c u of layer 2 is 120 kPa, and the corresponding COV is 0.3. The parameter distributions of two layers of soil are both lognormal distributions. For the two layers the θ h is 12 m, the θ v is 1.2 m, and the unit weight is 19 kN/m³. Following Jiang et al. (2017), the slope geometry is discretized into a finite number (856) of random field elements, and the size of the element is kept consistent with Jiang et al. (2017) for convenient comparison. Deterministic slope stability analysis is then performed with 6,231 PSSs using Bishop’s simplified method. The FS is calculated as 1.447, and the corresponding CSS passes through the slope toe. These results are comparable with those reported by Jiang et al. (2017).

With the AKCBR, 181 RSSs are identified, as schematically shown in Figure 9B, where the CSS is also included. Note that Wang et al. (2020) also previously analyzed the reliability of this example with the same approach but using a different number of RSSs. The reason is also two-fold and is the same as stated in the last example. The P f is subsequently estimated as 5.27×10⁻⁴ based on the RSS-RSM, which is comparable with the value of 5.13×10⁻⁴ by Method IV (Jiang et al., 2017) and those from other methods, as shown in Table 2. The consistency of the reliability results between different methods thus validates the accuracy of the AKCBR. Note that, however, the N s s of 2 × 10 5 for the AKCBR, estimated as that in the last example, is far less than that of 6 × 10 6 by Method IV. This consequently validates the computation efficiency of AKCBR.