Reliability analysis of excavated slopes in undrained clay

Abstract

A novel approach based on the upper bound theory is proposed to assess the stability of excavated slopes with spatially variable clay in undrained conditions. The proposed random limit analysis is a combination of the deterministic slope stability limit analysis together with random field theory and Monte Carlo simulation. A series of analyses is conducted to verify the potential application of the proposed method and to investigate the effects of the soil undrained shear strength gradient and the spatial correlation length on slope stability. Three groups of potential sliding surfaces are identified and the occurrence probability of each sort of failure mechanism is quantified for various slope ratios. The proposed method is found to be feasible for evaluating slope reliability. The obtained results are helpful in understanding the slope failure mechanism from a quantitative point of view. The paper could provide guidance for slope targeted local reinforcement.

Graphical abstract

Keywords

slope stability / spatial variability / limit analysis / random field / clay

1 Introduction

Slope stability analysis is one of the fundamental problems in geotechnical engineering and has received extensive literature attention. Researchers have gained new perspectives for slope stability analysis over the past decades by developing various methods. Conventional methods, including the limit equilibrium method (LEM), limit analysis method, and finite element method have been extensively adopted to identify the critical sliding surface and the associated minimum factor of safety (FS) [1–5]. While these methods may present quantitative differences between the FS under some circumstances, these are typically not large in practical terms.

As presented by Dasaka and Zhang [6], natural soils exhibit strongly spatial variability in properties due to the complicated interactions between geology, topography and climate during the formation. It has been well recognized that in slope stability analysis, similar apparent conditions may in fact produce very different failure mechanisms and probability of failure, mainly due to the spatial variability of soil properties.

There is a growing body of literature that recognizes the importance of considering the inherent soil variability in various geotechnical problems [7–15]. With respect specifically to the slope reliability analysis, most works have employed the random limit equilibrium method (RLEM) and random finite element method. Cho [16] combined the spatial variability of the soil properties with Spencer’s LEM to calculate the slope reliability index, and, further, in Ref. [17], a modified RLEM without assuming the critical sliding surface was proposed. Javankhoshdel and Bathurst [18] estimated the FS for simple slopes with cohesive and c−φ spatially variable soils. They provided a series of design charts using conventional limit equilibrium-based circular sliding analyses. Jiang et al. [19] quantified the risk of slope failure considering the spatial variability of soil parameters within the LEM framework, and it was found that the proposed method could reasonably identify the critical failure modes efficiently for two-layered slopes. Griffiths and Fenton [20] pioneered investigation of the probability of failure for undrained slopes using a self-developed intrusive finite element code. Research on similar issues, considering different random variables and various conditions, can be found in studies such as Griffiths and Marquez [21], Griffiths et al. [22], and Zhu et al. [23]. Another technique is the random finite element limit analysis (RFELA) based on the lower and upper bound theorems of classical plasticity theory. Taking advantage of the finite element approach, it can seek out the weakest region and find the critical failure mechanism for slopes with arbitrary geometry, constitutive model, loading, and boundary conditions. Meanwhile, rigorous upper and lower solutions can be obtained by constructing allowable stress and velocity fields in random soils, respectively. The RFELA method has been developed and applied to assess the reliability for slopes in purely cohesive soils [24] and c–φ soils [25].

Though RFELA has shown its feasibility in the stability evaluation of slope, it is necessary to conduct meshing and search optimization similar to the finite element method, which is inefficient at computing the stochastic slope responses. Several approaches under the framework of the limit analysis have been proposed, based upon the idea of discretization, in order to overcome the limits of applicability in the slope reliability analysis. For instance, Zhang et al. [26] carried out a reliability analysis for soil and rock slopes, based on the upper bound theorem in plasticity theory and on Winner’s polynomial chaos expansion, to solve the random FS and the critical failure surface. Zhao et al. [27] back-calculated the shear strength parameters for a 3D slope by upper bound limit analysis theory, considering the randomness and uncertainty of geotechnical parameters. However, that research ignored the spatial correlation of soil properties, which essentially amounts to an assumption of the infinite scale of fluctuation.

In the present study, a novel approach based on the upper-bound limit analysis theorem is proposed to evaluate the stability of slopes with spatially variable soils. This approach has the potential to inherit the merit of the simplicity of RLEM and provide accurate stochastic responses with moderate computational costs. By this method, the effects of the spatial correlation structure of soil properties and the soil strength gradients on the slope reliability index for various slope ratios can be highlighted. The slope failure mechanism is classified into three modes and the proportions of failure mode are quantified under each combination of random inputs.

2 Methodology

2.1 Random field modeling

The spatial variability of soil properties has been widely reported in the literature. The inherent soil heterogeneity can be simulated as spatially correlated random variables using random fields. c_u is modeled as random fields and the probabilistic distribution of c_u is confirmed to follow a log-normal distribution [28]. The log-normal distributed random fields of c_u are utterly defined by the mean {\mu }_{{\mathrm{c}}_{\mathrm{u}}}, the coefficient of variation CO{V}_{{\mathrm{c}}_{\mathrm{u}}}, and the spatial correlation length l. The trend of {\mu }_{{\mathrm{c}}_{\mathrm{u}}} is assumed to be increased linearly with depth and is expressed by the equation below:

(1){\mu }_{{\mathrm{c}}_{\mathrm{u}}}=k\mathit{\text{z}}+{\mu }_{{\mathrm{c}}_{\mathrm{u}0}},

where {\mu }_{{\mathrm{c}}_{\mathrm{u}}} is the mean undrained shear strength at the ground surface, and k is the strength gradient which varies with depth z. The coefficient of variation CO{V}_{{\mathrm{c}}_{\mathrm{u}}} represents the dispersion of undrained shear strength, which is defined as:

(2)CO{V}_{{\mathrm{c}}_{\mathrm{u}}}=\frac{{\sigma }_{{\mathrm{c}}_{\mathrm{u}}}}{{\mu }_{{\mathrm{c}}_{\mathrm{u}}}},

where {\sigma }_{{\mathrm{c}}_{\mathrm{u}}} is the standard deviation of undrained shear strength. An exponential autocorrelation function is adopted in this study to describe the spatial correlation of soil strength between two arbitrary locations and expressed as:

(3)\rho =\exp \left(-\sqrt{{\left(\frac{2{\tau }_x}{l_x}\right)}^2+{\left(\frac{2{\tau }_y}{l_y}\right)}^2}\right),

where \rho is the autocorrelation function; {\tau }_x and {\tau }_y are the horizontal lag distance and vertical lag distance, respectively; l_x and l_y are the spatial correlation lengths in the horizontal and vertical directions, respectively. Anisotropic spatially variable random fields correspond to the case where l_x and l_y are not equal. For the isotropic case, and Eq. (3) can be simplified to:

(4)\rho =\exp \left(-\frac{2\tau }{l}\right).

Several methods have been developed to generate random fields, such as the local average method [29], the K−L expansion method [30], the midpoint method [19] and the spectral representation method [31]. The spectral representation method is adopted in this study due to high accuracy in terms of the autocorrelation function and lower-order moments. Detailed procedure for generation of stationary and non-stationary random fields via the spectral representation method can be found in Shu et al. [31].

2.2 Failure mechanism

As suggested by Gao et al. [32] and Shu et al. [33], three groups of sliding surfaces are identified for different ranges of sliding surface depth D (see in Fig.1), namely.

1) Toe circle (TC): the failure occurs along a circle arc passing through the toe and the depth of the sliding surface D is less than or equal to the slope height, H, as shown in Fig.1(a).

2) Deep toe circle (DTC): this is similar to the TC case, but D is greater than H, as shown in Fig.1(b).

3) Deep circle (DC): the critical sliding surface reaches the region beneath the toe with the exit point on the slope base, as shown in Fig.1(c). The horizontal contained angle of the connection between the exit point and the crest, β′, is introduced to distinguish the failure type from the others, i.e., β′ is inferior to β while D is greater than H.

2.3 Random limit analysis (RLA)

Drucker and Prager [34] first applied the limit analysis based on the principle of plastic mechanics to estimate the slope stability. The limit theory was put forward by combining the static field with the velocity field, and the upper- and lower-bound limit theorems were established. Within the framework of upper-bound limit analysis, the admissible velocity field is anticipated via virtual work principle, the following work-energy balance equation holds:

(5){\int }_{V^{\ast }}{\sigma }^{\ast }{\dot{\epsilon }}^{\ast }\text{d}V+{\int }_{{\mathrm{\Gamma }}^{\ast }}\text{d}{D}_{\mathrm{s}}={\int }_L{T}^{\ast }{\dot{v}}^{\ast }\text{d}s+{\int }_{V^{\ast }}{W}_{\mathrm{s}}{\dot{v}}^{\ast }\text{d}V.

In Eq. (5), {\dot{\epsilon }}^{\ast } and {\dot{v}}^{\ast } are the rate of plastic admissible strain and displacement, respectively, {\sigma }^{\ast } is the stress in the plastic zone, D_s is the energy dissipation along the sliding surface, and W_s is the self-weight of the sliding body.

Li and Lumb [35] pointed out that slope stability was controlled by the average soil strength rather than the soil strength at a specific location on the sliding surface, since soils behave plastically. To capture the influence of variation of the average soil strength on the slope stability, the proposed algorithm requires both the soil domain and the sliding surface to be discretized. The sketch of the discrete failure mechanism of a slope is shown in Fig.2. The soil domain is discretized by a family of uniformly-spaced points. The sliding surface is evenly divided into a series of segments P_nP_n+1 (n = 1, 2, 3,…), where O_n is the midpoint of each segment. For a random field realization, the generated random values of c_u are assigned to each point according to the coordinates S(x_n, y_m).

The circle-shaped region DCAB rotates as a rigid body about the center of rotation O with the angular velocity ω. The soil mass below the velocity discontinuity surface (or sliding surface) BD remains at rest. The assumed mechanism can be specified by three variables θ₀, θ_h, and r, where θ₀ and θ_h are the slope angles of the chords OB and OD, respectively, and r is the radius of rotation. From the geometrical relations, it is shown that the radius of rotation r can be expressed in terms of θ₀ and θ_h in the form:

(6)r=\frac{H}{\sin {\theta }_{\mathrm{h}}-\sin {\theta }_0}.

Further, the distance between the entry point and the crest, L, can be expressed by:

(7)L=r\frac{\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\text{h}}-{\theta }_0)}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\text{h}}}-r\frac{\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\text{h}}+{\beta }^{\prime })}{(\sin {\theta }_{\text{h}}-\sin {\theta }_0)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\text{h}}\mathrm{s}\mathrm{i}\mathrm{n}{\beta }^{\prime }},

for the DC failure mechanism or by:

(8)L=r\frac{\sin {\theta }_{\mathrm{h}}-\sin {\theta }_0}{\tan \beta }+r\cos {\theta }_{\mathrm{h}},

for the TC and DTC failure mechanisms.

The coordinates (x_c, y_c) of the center of rotation O can be determined from the following relations:

(9)\begin{array}{c}{x}_{\mathrm{c}}=H\cot \beta +L-r\cos {\theta }_0,\\ y_{\mathrm{c}}=H+r\sin {\theta }_0.\end{array}

The central angle of each segment P_n{P}_{n+1}, denoted as Δθ, and the angle between OO_n and the positive direction of the x-axis is found as:

(10){\theta }_n={\theta }_0+(n-1)\mathrm{\Delta }\theta +\frac{\mathrm{d}\theta }{2}.

The coordinates of the midpoint of each segment, O_n(x_n, y_n), can be obtained by substituting Eqs. (9) and (10) into the following equations:

(11)\begin{array}{c}{x}_n=x_{\mathrm{c}}+r\cos {\theta }_n,\\ y_n=y_{\mathrm{c}}-r\sin {\theta }_n.\end{array}

It is noteworthy that the midpoint of each segment O_n(x_n, y_n) might not always fall on the discrete point in the discretized slope domain. Therefore, the algorithm automatically searches out the discrete point nearest to the midpoint O_n, and then assigns the c_u value of this point to O_n in the calculation of energy dissipation. As illustrated in Fig.2, if the length of O_nS(x_n, y_m) is smaller than that of O_nS(x_n+1, y_m), O_nS(x_n, y_m+1), and O_nS(x_n+1, y_m+1), then the c_u value of the point S(x_n, y_m) is assigned to the midpoint O_n to represent the random strength property for this segment.

Hence the work rate of energy dissipation for an arbitrary segment can be solved by:

(12)\mathrm{\Delta }{D}_n=c_{\mathrm{u}}{r}^2({\theta }_{\mathrm{h}}-{\theta }_0)\omega .

Consequently, the total dissipation on the sliding surface can be calculated from the following expression:

(13)D=\sum _{n=1}^i\mathrm{\Delta }{D}_n.

The rate of work of the external force provided by the soil weight is calculated as the dot product of the total weight of the block ABDC and the angular velocity. The detailed derivation process can be found in Ref. [36].

Finally, the FS can be obtained via the energy balance equation. Although the concept of FS is well established and widely employed in geotechnical problems, there is no universal definition. The FS in this study is defined as the ratio of the undrained shear strength to that needs to maintain limit equilibrium, as:

(14)FS=\frac{c_{\mathrm{u}}}{c_{\mathrm{u}\mathrm{d}}},

where c_ud denotes the undrained shear strength necessary only to maintain the slope in limit-equilibrium state. The least upper-bound solution to the safety factor was found using the optimization procedure based on the random search technique adopted by Chen [37]. Changing the variables θ₀, θ_h, and β′ can be used to calculate the minimum value of FS until the variables become less than predefined values (i.e., 0.01° for angles).

2.4 Reliability index

The most effective applications of probability theory, to the slope stability analysis, state the uncertainties in the form of a reliability index, as done by Christian et al. [38]. If the computed FS is normally distributed, the reliability index β_S can be defined as:

(15){\beta }_{\mathrm{S}}=\frac{{\mu }_{\mathrm{F}\mathrm{S}}-1}{{\sigma }_{\mathrm{F}\mathrm{S}}},

where {\mu }_{\mathrm{F}\mathrm{S}} and {\sigma }_{\mathrm{F}\mathrm{S}} are the mean and the standard deviation of FS, respectively. For a log-normally distributed FS, the estimation of the reliability index is given by:

(16){\beta }_{\mathrm{S}}=\frac{\text{ln}{\displaystyle \frac{{\mu }_{\mathrm{F}\mathrm{S}}}{1+{\left({\sigma }_{\mathrm{F}\mathrm{S}}/{\mu }_{\mathrm{F}\mathrm{S}}\right)}^2}}}{\sqrt{\mathrm{l}\mathrm{n}\left[1+{\left({\sigma }_{\mathrm{F}\mathrm{S}}/{\mu }_{\mathrm{F}\mathrm{S}}\right)}^2\right]}}.

2.5 Implementation procedure

The implementation procedure for the proposed RLA method is summarized as follows.

Step 1: definition of random soil properties, i.e., probabilistic distribution, mean value, coefficient of variation, and spatial autocorrelation structure.

Step 2: establishment of the slope stability model and discretization of the slope domain and the assumed failure mechanism.

Step 3: generation of random fields of c_u and assignment to each discretized point in the slope domain.

Step 4: identification of the appropriate c_u according to coordinates to represent the property of each segment of sliding surface.

Step 5: calculation of the work rates of external force and energy dissipation and of the slope FS.

3 Results and discussion

3.1 Case description and model validation

The slope geometry considered in this study is shown in Fig.3, with a slope height H = 10 m and a varying slope ratios Rs of 1:1, 1:2, and 1:3, corresponding to slope inclinations β = 45°, 26.6°, and 18.4°, respectively. The soil unit weight \gamma is assumed to be 16 kN/m³. The soil under undrained loading conditions is modeled by a perfectly plastic Tresca yield criterion with an associated flow rule. The slope is analyzed under undrained conditions with the undrained shear strength characterized statistically by Eq. (1). {\mu }_{{\mathrm{c}}_{\mathrm{u}0}} is set to be 30 kPa and k values are set to be 0, 0.5, 1, 1.5, 2, 2.5, and 3 kPa/m. The selected parameters fall in the typical ranges of normally consolidated clay, as Koppula [39] and Javankhoshdel et al. [40] reported. It has been suggested that typical values of the coefficient of variation for undrained shear strength fall in the range from 0.1 to 0.5 [41]. Therefore, CO{V}_{{\mathrm{c}}_{\mathrm{u}}} is chosen as 0.3 in this study. Both the isotropic and anisotropic spatial correlations lengths are considered. For the isotropic case, the spatial correlation lengths are set to be l = l_x = l_y = 5, 10, 15, 20, 25, 30, 35, and 40 m, while l_y = 5 m and l_x = 10, 15, 20, 25, 30, 35, and 40 m are adopted from the anisotropic case.

The accuracy of the proposed RLA depends on the number of segments into which the sliding surface and the soil domain are divided. In this study, the sliding surface is divided into 100 segments and the distance between each discretized point of the soil domain is 0.5 m. The present model is verified with the deterministic results provided in Ref. [42]. Given a uniform 1:1 slope with H = 10 m, c_u = 30 kPa, and \gamma = 16 kN/m³, Michalowski [42] provides the stability number N = c_u/H\gammaFS = 0.179. The result of this study gives N = 0.178, which compares well with the theoretically exact solution.

Fig.4 presents the deterministic sliding surfaces for different slope ratios. It can be seen that the parameter Rs and k have significant influences on the critical sliding surface. For Rs = 1:3, the sliding surfaces are of DC type, and the depth of sliding surfaces decreases gradually with the increase of k (Fig.4(a)). The DC mechanism dominates the sliding surfaces when k is between 0 to 1.5 kPa/m and then translates to DTC as k increases from 2 to 3 kPa/m, as shown in Fig.4(b). In Fig.4(c), the failure mechanism becomes independent of k when k is greater than 0, which falls within the group of DC, while the sliding surface belongs to DTC for k = 0.

Monte Carlo (MC) simulations are performed for each random field realization. The accuracy of MC simulation is highly dependent on the number of simulations. Fig.5 plots the results of MC simulation for the case of Rs = 1:2, k = 1 kPa/m, and l = 10 m. As indicated in Fig.5(a), convergences are observed for both the mean and standard deviation of FS (μ _FS and σ_FS) at approximately the 400th simulation. Hence, 1000 MC simulations are carried out for each combination of random parameter inputs throughout the analysis. Fig.5(b) shows a 10-bin histogram of FS together with a fitted normal distribution. The normal distribution can fit the probability distribution function (PDF) columns well, characterizing FS at a 5% significance level based in the chi-square tests. Hence, Eq. (15) can be adopted to quantify the reliability of slopes in the following sections. Another exact comparison is made with the solutions provided by OptumG2 software for the same stochastic case. In the OptumG2 analysis, a triangle upper bound element type is used and the number of elements is set to be 2000. After 1000 MC simulations, the reliability index performed by OptumG2 software is β_S = −0.90 and the proposed RLA method gives β_S = −0.85. The two results are similar with a difference being less than 6%, which means that the proposed RLA method can perform well in a probabilistic framework for the slope reliability analysis. Moreover, the computational efforts for the OptumG2 software and the proposed method to conduct 1000 MC simulations are approximately 70 and 10 min on the authors’ computer (3.40GHz CPU and 64G RAM), respectively, indicating that the proposed method can obviously reduce the computational cost with an identical accuracy level.

3.2 Effect of isotropic spatial correlation lengths

Fig.6 shows the variation of β_S as a function of l for different values of Rs. In most cases, β_S decreases with the increase of l, except in the case of k = 0. The reduction in β_S is mainly because the spatial variability of c_u is averaged over a region or sliding surface as l becomes large. Fig.6(b) and Fig.6(c) clearly indicate two branches, with β_S less or greater than zero for stationary and non-stationary random fields, respectively. However, the stochastic results for the cases, such as k = 0 and l = 5 m, give β_S being less than 0. The findings emphasize that the deterministic slope stability analysis might lead to a conservative result and underestimate the risk of failure.

The influence of k on β_S is further plotted in Fig.7. For all the cases considered, β_S increases with increasing k, which is expected. Comparison of Fig.7(a) with Fig.7(b) and Fig.7(c) highlights the importance of slope ratio in the relationship between l and k to the slope reliability: under a fixed combination of k and l, the greater the slope ratio, the smaller the reliability index is; a large slope ratio results in similar β_S, especially when l is large and k is small.

The reliability index method of slope stability evaluation is effective and apparent. However, the stability analysis based on the reliability index is not fully informative because identical slopes with the same reliability index may exhibit different failure modes. As an illustration of the results, 100 sliding surfaces randomly selected from 1000 realizations for the cases of k = 0 and 3 kPa/m, and l = 10 and 30 m, respectively, with Rs = 1:1 are shown in Fig.8. The corresponding deterministic sliding surfaces are also depicted in red. The soil spatial variability caused various failure mechanisms, which are not manifested in the deterministic analysis. The following section will evaluate the variation of the distance between the entry point and the crest L and the occurrence of probability of each sliding surface type from a qualitative perspective.

The computed results of the mean distance from the entry point to the crest, μ _L, are listed in Tab.1, and some of the computed results are reproduced in Fig.9 and Fig.10 for illustrative purposes. The corresponding deterministic results of L are tabulated in Tab.2. A direct comparison between Tab.1 and Tab.2 leads to some interesting points: (i) results derived from the stationary random field (i.e., k = 0) are always smaller than the corresponding deterministic results; (ii) for 1:2 and 1:3 slopes in non-stationary random fields, the stochastic results are smaller than the corresponding deterministic results when k is small and the effect of l seems to be negligible and (iii) for steep slopes with Rs = 1:1, the stochastic results in non-stationary random fields are greater than the corresponding deterministic results. In addition, σ_L generally diminishes with k and l, indicating the dispersion degree of L reduces, and the large slope ratio corresponds to the great value of σ_L. Changes in the failure mechanism caused by soil spatial variability might contribute to these notable observations.

To further demonstrate the effect of k on μ _L, an example under l = 10 m is selected as plotted in Fig.9. It shows that μ _L decreases with the increase of k, with apparently small values for Rs = 1:1 compared with Rs = 1:2 and 1:3. Besides, the decreasing rate of μ _L increases with the rise of Rs. Note that μ _L of the slope with Rs = 1:2 is largest at k = 0 irrespective of l, as can be seen in Tab.1 and Fig.9.

As shown in Fig.10 for the cases of k = 1 kPa/m, a consistent pattern of variation in μ _L is found where it first increases and then decreases as l increases with a critical point observed at around l = 15 m. As l extends over 25 m, the stochastic response of L almost remains constant. The slope ratio effect gives a similar trend to that observed previously in Fig.9, namely the position of the entry point moves toward the crest when the slope becomes steep.

The occurrence probabilities of each failure mode are summarized in Tab.3 for a total of 168 combinations for l, k, and Rs. The slope failure mechanisms are classified into three types, as previously mentioned and shown in Fig.1. For gentle slopes (i.e., Rs = 1:3 and 1:2), the behaviors of the failure mechanism are similar. There is no occurrence of TC with the variation of k and l, and the DC failure mechanism is significantly dominant. As the slope becomes steep (i.e., Rs = 1:1), the failure mode of TC emerges whereas its occurrence probability is consistently lower than that of DTC and DC. It is notable that the proportion of DTC becomes prevalent for large values of k and l, such as the cases with k ≥ 2 kPa/m and l ≥ 35 m. The observations emphasize the importance of the slope ratio on the failure mechanism in spatially variable soils.

To further explore the influence of k and l on the occurrence of probability for each failure mechanism, the results under a fixed k (i.e., 1 kPa/m) and a fixed l (i.e., 10 m) are reproduced in Fig.11 and Fig.12 in which the probability of occurrence for each failure mode is counted in histogram charts. Fig.11 indicates that the influence of k on the occurrence probability of various failure modes is pronounced. A rise in k parallels an increase in the percentage of DTC but a reduction in the percentage of DC for the 1:3 and 1:2 slopes (Fig.11(a) and 11(b)). As plotted in Fig.11(c), there is a positive correlation between k and the occurrence probabilities of TC and DTC; meanwhile, the percentage of DTC is always higher than that of TC. Conversely, the percentage of DC failure mode decreases gradually as k increases. With the increase of k, more soils with large values of c_u are present in the deeper zone, inducing the sliding surface depth to become shallow and the exit point to move toward the toe.

The influence of l on the variation of failure modes is plotted in Fig.12 for the cases with k = 1 kPa/m. Fig.12 demonstrates that the increase of l from 5 to 20 m causes the proportion of the DC mechanism to diminish first and then to augment for the slopes with Rs = 1:3 and 1:2, and the reverse pattern occurs for the DTC mechanism. The critical point is found at l = 20 m, which is equal to the restrained depth. In the 1:1 slopes, there is a gradual reduction of the occurrence probability of DC as l increases, and in turn, also of the occurrence probability of DTC. The variation of TC is similar compared with that of μ _L shown in Fig.10, so it could be speculated that l has no impact on the failure mode as l is greater than 20 m. Fig.12 also shows that, under a given combination of slope ratio and soil strength gradient, the fluctuation in the proportions of each failure mode is not obvious for gentle slopes (Rs = 1:2 and 1:3) and the effect is opposite for steep slopes (Rs = 1:1).

3.3 Effect of anisotropic spatial correlation lengths

The spatial correlation length of soil properties in the horizontal direction is likely to be greater than in the vertical direction. To observe the effect of anisotropic spatial variability of undrained shear strength on the slope reliability index and failure mechanism, l_y is fixed at 10 m l_x varies from 10 to 40 m and k is equal to two values, 0 and 1 kPa/m.

The results of reliability index are shown in Fig.13. As for the isotropic cases (Fig.6), β_S increases or decreases monotonically with the increase of l. However, a pronounced critical point can be observed at l_x = 30 m for the isotropic spatial variability cases. The results highlight that the spatial correlation length should be carefully selected in the slope design and stability evaluation if the anisotropic spatial variability needs to be considered.

The variation of μ _L with l_x under k = 0 and 1 kPa/m is illustrated in Fig.14. In each figure, the variation trend of μ _L is consistent while it is difficult to identify a clear relationship between the stationary and the non-stationary cases. Another observation is that, as is apparent in Fig.14(a), the values of μ _L with Rs = 1:3 and 1:1 are very similar and smaller than those with Rs = 1:2. However, as depicted in Fig.14(b), a steep slope corresponds to a small μ _L at a fixed l_x.

The occurrence probabilities of each failure mode for the anisotropic case are summarized in Tab.4 and reproduced in Fig.15 and Fig.16.

In general, the failure mode for Rs = 1:2 and 1:3 is composed of DC and DTC modes and the TC mode only occurs when Rs = 1:1, which is a similar pattern to that of the results with isotropic spatial variability. It seems that l_x has negligible effect on the variation of failure mechanism. For example, in Fig.15(a), the difference between the maximum and minimum occurrence probability of the DC failure mode is less than 1% when l_x extends from 10 to 40 m.

4 Conclusions

Probabilistic slope stability assessments are performed within the framework of upper-bound limit analysis, random field theory and MC simulation. The influences of the clay undrained shear strength gradient and of both isotropic and anisotropic spatial correlation lengths on slope stability are investigated. In addition, the failure mechanism under various random fields is classified and discussed. The following conclusions can be drawn.

1) The proposed method, which is a combination of random field and upper bound limit analysis, is shown to be capable of efficiently evaluating the reliability index and the failure mechanism of slopes caused by spatial variability of soil properties.

2) The combination of slope ratio, strength gradient and spatial correlation length significantly impacts slope stability. There is a gradual reduction of the slope reliability index with increase of the slope ratio. For a given slope ratio, the slope reliability index is positively correlated with the strength gradient and negatively correlated with the isotropic spatial correlation length. When the anisotropic spatial correlation length is considered, the slope reliability index first increases and then decreases, with a critical point found at 30 m.

3) The failure mechanism of the DC, which is reached at an occurrence probability of 97%, is predominant in all the cases considered. For slope ratios equaling 1:3 and 1:2, no TC failure mechanism is observed.

4) With the increase of soil strength gradient, the proportion of the DC failure mode reduces and the proportion of the circle or DTC failure modes increases accordingly. The effect of the isotropic spatial correlation length on the slope failure mechanism is complicated. It is found that, under a given combination of slope ratio and soil strength gradient, the fluctuation in the proportions of each failure mode is not obvious for gentle slopes and the effect is opposite for steep slopes. When the spatial correlation length is greater than 20 m, the spatial correlation length effect becomes vanishingly small. In addition, the horizontal spatial correlation length has negligible effect on the variation of failure mechanism when the anistropic spatial variability is considered.

References

[1] Griffiths D V, Yu X. Another look at the stability of slopes with linearly increasing undrained strength. Geotechnique, 2015, 65(10): 824–830

[2] Lim K, Li A J, Lyamin A V. Three-dimensional slope stability assessment of two-layered undrained clay. Computers and Geotechnics, 2015, 70: 1–17

[3] Michalowski R L. Slope stability analysis: A kinematical approach. Geotechnique, 1995, 45(2): 283–293

[4] Renani R H, Martin D C. Factor of safety of strain-softening slopes. Journal of Rock Mechanics and Geotechnical Engineering, 2020, 12(3): 473–483

[5] Yu H S, Salgado R, Sloan S W, Kim J M. Limit analysis versus limit equilibrium for slope stability. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(1): 1–11

[6] Dasaka S M, Zhang L M. Spatial variability of in situ weathered soil. Geotechnique, 2012, 62(5): 375–384

[7] Griffiths D V, Fenton G A. Bearing capacity of spatially random soil: The undrained clay Prandtl problems revisited. Geotechnique, 2001, 51(4): 351–359

[8] Jiang S H, Huang J, Griffiths D V, Deng Z P. Advances in reliability and risk analyses of slopes in spatially variable soils: A state-of-the-art review. Computers and Geotechnics, 2022, 141: 104498

[9] Li J, Hu P, Uzielli M, Cassidy M J. Bayesian predication of peak resistance of a spudcan penetrating sand-over-clay. Geotechnique, 2018, 68(10): 905–917

[10] Montgomery J, Boulanger R W. Effects of spatial variability on liquefaction-induced settlement and lateral spreading. Journal of Geotechnical and Geoenvironmental Engineering, 2017, 143(1): 04016086

[11] Roberts L A, Misra A. Reliability-based design of deep foundations based on differential settlement criterion. Canadian Geotechnical Journal, 2009, 46(2): 168–176

[12] Soubra A H, Al-Bittar T, Thajeel J, Ahmed A. Probabilistic analysis of strip footings resting on spatially varying soils using kriging metamodeling and importance sampling. Computers and Geotechnics, 2019, 114: 103107

[13] Wang Y, Cao Z J, Li D. Bayesian perspective on geotechnical variability and site characterization. Engineering Geology, 2016, 203: 117–125

[14] Zhang T T, Guo X F, Baroth J, Dias D. Metamodel-based slope reliability analysis—Case of spatially variable soils considering a rotated anisotropy. Geosciences, 2021, 11(11): 465

[15] Zhang T T, Guo X F, Dias D, Sun Z B. Dynamic probabilistic analysis of non-homogeneous slopes based on a simplified deterministic model. Soil Dynamics and Earthquake Engineering, 2021, 142: 106563

[16] Cho S E. Effects of spatial variability of soil properties on slope stability. Engineering Geology, 2007, 92(3−4): 97–109

[17] Cho S E. Probabilistic assessment of slope stability that considers the spatial variability of soil properties. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(7): 975–984

[18] Javankhoshdel S, Bathurst R J. Simplified probabilistic slope stability design charts for cohesive and c−φ soils. Canadian Geotechnical Journal, 2014, 51(9): 1033–1045

[19] Jiang S H, Liu X, Huang J, Zhou C B. Efficient reliability-based design of slope angles in spatially variable soils with field data. International Journal for Numerical and Analytical Methods in Geomechanics, 2022, 46(13): 2461–2490

[20] Griffiths D V, Fenton G A. Probabilistic slope stability analysis by finite elements. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(5): 507–518

[21] Griffiths D V, Marquez R M. Three-dimensional slope stability analysis by elasto-plastic finite elements. Geotechnique, 2007, 57(6): 537–546

[22] Griffiths D V, Huang J, Fenton G A. Probabilistic infinite slope analysis. Computers and Geotechnics, 2011, 38(4): 577–584

[23] Zhu D, Griffiths D V, Huang J, Fenton G A. Probabilistic stability analyses of undrained slopes with linearly increasing mean strength. Geotechnique, 2017, 67(8): 733–746

[24] Huang J, Lyamin A V, Griffiths D V, Krabbenhoft K, Sloan S W. Quantitative risk assessment of landslide by limit analysis and random fields. Computers and Geotechnics, 2013, 53: 60–67

[25] Chwała M. Upper-bound approach based on failure mechanisms in slope stability analysis of spatially variable c−φ soils. Computers and Geotechnics, 2021, 135: 104170

[26] ZhangHZhaoYLiC. Probabilistic slope stability analysis based on the upper bound theorem. In: Proceedings of 2010 International Conference on E-Product E-Service and E-Entertainment. Henan: IEEE, 2010, 1–4

[27] Zhao L H, Zuo S, Lin Y L, Li L, Zhang Y B. Reliability back analysis of shear strength parameters of landslide with three-dimensional upper bound limit analysis theory. Landslides, 2016, 13(4): 711–724

[28] Phoon K K, Kulhawy F H, Grigoriu M D. Reliability-based design for transmission line structure foundations. Computers and Geotechnics, 2000, 26(3−4): 169–185

[29] Fenton G A, Vanmarcke E H. Simulation of random fields via local average subdivision. Journal of Engineering Mechanics, 1990, 116(8): 1733–1749

[30] Li D, Qi X, Cao Z, Tang X, Zhou W, Phoon K K, Zhou C. Reliability analysis of strip footing considering spatially variable undrained shear strength that linearly increases with depth. Soil and Foundation, 2015, 55(4): 866–880

[31] Shu S, Gao Y, Wu Y. Probabilistic bearing capacity analysis of spudcan foundation in soil with linearly increasing mean undrained shear strength. Ocean Engineering, 2020, 204: 106800

[32] Gao Y, Zhang F, Lei G, Li D, Wu Y, Zhang N. Stability charts for 3D failures of homogeneous slopes. Journal of Geotechnical and Geoenvironmental Engineering, 2013, 139(9): 1528–1538

[33] Shu S, Ge B, Wu Y, Zhang F. Probabilistic assessment on 3D stability and failure mechanism of undrained slopes based on the kinematic approach of limit analysis. International Journal of Geomechanics, 2023, 23(1): 06022037

[34] Drucker D C, Prager W. Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics, 1952, 10(2): 157–165

[35] Li K S, Lumb P. Probabilistic design of slopes. Canadian Geotechnical Journal, 1987, 24(4): 520–535

[36] ChenW F. Limit Analysis and Soil Plasticity. Amsterdam: Elsevier, 1975

[37] Chen Z Y. Random trials used in determining global minimum factors of safety of slopes. Canadian Geotechnical Journal, 1992, 29(2): 225–233

[38] Christian J T, Ladd C C, Baecher G B. Reliability applied to slope stability analysis. Journal of Geotechnical Engineering, 1994, 120(12): 2180–2207

[39] Koppula S D. On stability of slopes in clays with linearly increasing strength. Canadian Geotechnical Journal, 1984, 21(3): 577–581

[40] Javankhoshdel S, Cami B, Chenari R J, Dastpak P. Probabilistic analysis of slopes with linearly increasing undrained shear strength using RLEM approach. Transportation Infrastructure Geotechnology, 2021, 8(1): 114–141

[41] Jiang S H, Huang J, Qi X H, Zhou C B. Efficient probabilistic back analysis of spatially varying soil parameters for slope reliability assessment. Engineering Geology, 2020, 271: 105597

[42] Michalowski R L. Stability charts for uniform slopes. Journal of Geotechnical and Geoenvironmental Engineering, 2002, 128(4): 351–355

RIGHTS & PERMISSIONS

Higher Education Press