An Efficient Solution for Probabilistic Slope Seismic Stability Analysis Based on Polynomial Chao Kriging Metamodel

Abstract

Slope stability analysis plays a crucial role in geotechnical engineering, particularly in regions susceptible to seismic activity. The inherent non-homogeneity and uncertainty of soil properties pose significant challenges in assessing slope stability under seismic conditions. To address these complexities, a novel and efficient methodology named DUBLA-PDM-PCK is proposed. In this methodology, the effects of soil non-homogeneity and uncertainty, along with the time and spatial variations of seismic loading, are systematically considered. The deterministic framework integrates discretized upper bound limit analysis (DUBLA) to accommodate soil non-homogeneous characteristics, and the pseudo-dynamic method (PDM) to model seismic loading variability. Then, a robust and efficient probabilistic analysis method, PCK-MA, is implemented utilizing adaptive Polynomial Chaos Kriging metamodeling, Monte Carlo Simulation, and Analysis of Covariance to investigate the uncertainty of the parameters. This approach treats nine key parameters, including soil cohesion, friction angle, non-homogeneous coefficients, horizontal and vertical seismic coefficients, period, and amplification factor, as random variables to assess their uncertainty effects on failure probability (stability level) and sensitivity indices. The DUBLA-PDM-PCK methodology offers a streamlined and reliable tool tailored for assessing slope stability in seismic environments, demonstrating notable efficiency in addressing soil variability and seismic loading uncertainties. Its application holds promise for guiding engineering practices and enhancing understanding of slope behavior in regions prone to seismic hazards.

1. Introduction

Landslides, characterized by the downward movement of rock, soil, and debris on a slope, represent prevalent geological disasters worldwide. These events occur as a result of various factors, such as heavy rainfall, human-induced modifications, or earthquakes, with earthquakes being particularly significant triggers of instability in susceptible regions. Earthquake-induced landslides pose substantial threats to natural environments and human infrastructure, causing widespread damage, altering landscapes, and resulting in loss of life. Moreover, inherent soil variability due to natural geological processes persists, and this variability can influence the probability of slope failure. Therefore, gaining a deeper understanding of slope seismic stability while considering uncertainties is crucial for mitigating the impacts of geohazards.

Numerous studies have been conducted offering valuable insights into the analysis of slope seismic stability [1,2,3,4,5,6]. However, certain limitations and shortcomings persist, which will be discussed in the following sections.

(1)

The pseudo-static method is widely adopted for its simplicity in considering seismic effects [5,7]. However, this method treats seismic effects as uniform inertial forces with constant accelerations, neglecting the dynamic aspects of earthquakes, including acceleration amplification, shaking duration, and intensity–frequency characteristics [1,8,9]. To address these limitations, the pseudo-dynamic approach has been introduced to account for time and space variations in ground shaking [10]. This method has been applied in various geotechnical contexts, including retaining walls [11], slopes [12], and tunnels [13]. However, it is predominantly used for deterministic analyses, often overlooking parameter uncertainties such as variations in amplification factor or seismic period, aspects that will be explored in this study.

(2)

Deterministic methods commonly employed for slope seismic analysis include the limit equilibrium method (LEM) [14,15], numerical simulation [2,6], and upper bound limit analysis (UBLA) [1,16,17]. However, the LEM lacks the capability to consider the soil stress–strain relationship, which can lead to biased outcomes due to assumptions regarding failure surfaces and inter-slice forces [18]. While numerical simulations can yield reasonable results, they often come with heavy computational burdens. The upper bound limit analysis offers a good balance of accuracy and efficiency; however, traditional approaches like the log-spiral mechanism may involve tedious integral calculations, especially for non-homogeneous cases [1]. To address these challenges, the discretized upper bound limit analysis (DUBLA) was introduced based on a “point-to-point” technique. DUBLA allows for the creation of a discretized mechanism with spatially varying soil parameters while ensuring computational efficiency by avoiding mass discretization [19]. In this study, DUBLA is utilized for slope seismic stability analysis to overcome the limitations of traditional methods and provide accurate and efficient results.

(3)

The probabilistic analysis of slope seismic stability employs various methods, including Monte Carlo Simulation (MCS) [2,3,4], Subset Simulation (SS) [12], and the First Order Reliability Method (FORM) [20]. However, sampling methods like MCS and SS are time-consuming as they require multiple simulations. Approximation methods such as FORM may not provide accurate estimates for complex cases due to assumptions of linear limit state surfaces. The active learning method has received much attention with the development of surrogate models (such as the metamodel and Artificial Neural Network) as cheap proxies of computational models [21]. Pan et al. [22] utilized Polynomial Chaos Expansion (PCE) metamodels to enhance the efficiency of slope seismic stability analyses, demonstrating the effectiveness of PCE metamodels in improving computational efficiency and accuracy in these analyses.

To enhance the performance of the active learning method in slope stability analysis, this study implements an advanced approach based on Polynomial Chaos Kriging (PCK) metamodeling. PCK combines the benefits of Polynomial Chaos Expansion (PCE) and Kriging interpolation, offering a more efficient solution compared to using PCE and Kriging independently. Additionally, traditional probabilistic methods such as Monte Carlo Simulation (MCS) are applied to generate comprehensive results, including Failure Probability (P_f), Probability Density Function (PDF), Cumulative Distribution Function (CDF), and statistical moments of the system response. This integrated approach leverages the strengths of PCK metamodeling and traditional probabilistic methods to enhance the accuracy and efficiency of slope stability analysis.

Furthermore, in probabilistic analysis, several random variables are considered, and it is valuable to investigate the impact of these input variables on the model response through sensitivity analyses. Global Sensitivity Analysis (GSA), which comprehensively assesses the entire space of input variables, is a widely utilized technique [23,24]. In this study, Analysis of Covariance (ANCOVA), a type of GSA, will be employed to generate informative sensitivity indices, particularly for correlated input variables. This approach aims to provide insights into the relative importance and interactions of different parameters on the outcomes of the slope stability analysis.

This study introduces an efficient methodology, DUBLA-PDM-PCK, which combines DUBLA, the PDM, and PCK to facilitate a rapid and reliable assessment of slope stability in seismic-prone regions with consideration of the soil non-homogeneity and parameter uncertainty effects. The methodology is initially presented, followed by its implementation. The effectiveness of DUBLA-PDM-PCK is then validated through comparison with existing studies. Discussions are conducted on topics related to slope seismic stability and uncertainties in parameters. Additionally, sensitivity analyses are performed to determine the importance of different parameters in the analysis.

2. Deterministic Seismic Analysis of the Slope Stability

The analysis focuses on the seismic stability of a non-homogeneous slope. This section aims to outline the fundamental principles and the optimization process for determining the safety factor.

2.1. Pseudo-Dynamic Method

The pseudo-static method simplifies seismic loadings by assuming uniform accelerations applied to structures, where both horizontal and vertical seismic accelerations remain constant. This approach overlooks the dynamic nature of earthquake motions. In contrast, the pseudo-dynamic approach accounts for amplitude and phase variations in seismic vibrations, as well as the period of ground shaking. This method offers a more realistic representation of dynamic effects compared to the pseudo-static method, while retaining computational simplicity relative to fully dynamic analyses [1,12].

The propagation of seismic waves in soils is influenced by wave velocities, specifically shear wave velocity $V_{\mathrm{s}}$ and primary wave velocity $V_{\mathrm{p}}$. These velocities can be expressed as:

$$\{\begin{array}{l}{V}_{\mathrm{s}}=\sqrt{G/\rho }\hfill \\ V_{\mathrm{p}}=\sqrt{[2G(1-\nu )]/[\rho (1-2\nu )]}\hfill \end{array}$$

(1)

where $G$, $\rho$, and $\nu$ are the soil shear modulus, density, and Poisson’s ratio.

The vibrational effects are amplified as seismic waves propagate toward the ground surface. The PDM assumes the magnitudes of the horizontal and vertical seismic accelerations ($a_h$ and $a_v$) increase linearly from the slope base to the ground surface with an amplification factor $f$, and the seismic accelerations can be represented by $f{a}_h$ and $f{a}_v$.

Additionally, a sinusoidal function is used to model the seismic wave, which simplifies analysis and enhances comprehension. Both horizontal and vertical shakings have the same period, denoted as T. The horizontal acceleration and vertical acceleration for any time t and depth h, can be determined by

$$\{\begin{array}{c}{a}_h=[1+\frac{h}{H}(f-1)]\cdot k_{h}g\cdot \sin {\omega }_e(t-\frac{h}{V_s})\\ a_v=[1+\frac{h}{H}(f-1)]\cdot k_{v}g\cdot \sin {\omega }_e(t+t_0-\frac{h}{V_{\mathrm{p}}})\end{array}$$

(2)

where $k_h$ and $k_v$ are the horizontal and vertical seismic coefficients, respectively, and $g$ is the gravity acceleration. ${\omega }_e$ denotes the angular frequency of the seismic shaking and can be determined by ${\omega }_e=2\pi /T$ for the sinusoidal wave [1,12,16].

2.2. Discretized Upper Bound Limit Analysis

Figure 1a illustrates the mechanism of the discretized upper bound limit analysis (DUBLA). H and β, respectively, represent the slope depth and angle. A coordinate system is established with the slope toe C as the origin point. The mechanism rotates around the rotational center O with an angular velocity ω. The mechanism is related to the determination of two parameters, including the $r_0$ (length of OC) and ${\theta }_0$ (the angle between x axis direction and line OC). The slip surface AC is discretized into a series of straight segments $P_i{P}_{i+1}$. It is important to note that each segment $P_i{P}_{i+1}$ should form an angle ${\phi }_i$ with the velocity vector $v_i$ to satisfy the associated flow rule. Therefore, the generation process aims to define all the discretization points $P_i$.

A “point-to-point” technique is implemented to determine $P_i$ along the critical slip surface, where each point is derived sequentially from the previous one. $\mathsf{\delta }\theta$ is also introduced to control the accuracy of the failure mechanism and it is often set as 0.1°, which is recommended as a great combination of accuracy and time cost [25]. The formulation of the $P_i$ generation can be expressed as

$$\{\begin{array}{c}{x}_{i+1}=x_i+\frac{\sqrt{{(x_i-x_O)}^2+{(y_i-y_O)}^2}\times \sin \delta \theta }{\sin (\pi /2+\phi -\delta \theta )}\times \cos ({\theta }_i-\pi /2+{\phi }_i)\\ y_{i+1}=y_i+\frac{\sqrt{{(x_i-x_O)}^2+{(y_i-y_O)}^2}\times \sin \delta \theta }{\sin (\pi /2+\phi -\delta \theta )}\times \sin ({\theta }_i-\pi /2+{\phi }_i)\end{array}$$

(3)

where $x_O$ and $y_O$ are point $O$ coordinates in horizontal and vertical directions, and the polar angle for line $l_i$ is ${\theta }_i={\theta }_0+\left(i-1\right)\delta \theta$. A linear interpolation method is subsequently implemented to adjust the last generated point on the slope crest (point A).

2.3. Stability Analysis

The kinematic approach uses the work rate balance equation [26], which states that the external work rate W is equal to the energy dissipation D, to evaluate the sliding event. Herein, the energy dissipation D occurred along the slip surface, which is governed by an associated flow rule and is calculated by summing the energy dissipation of discretized elements due to the sliding surface discretization:

$$D=\omega {\displaystyle \sum c_i{L}_i{R}_i}\cos {\phi }_i$$

(4)

where $c_i$ and ${\phi }_i$ denote the cohesion and friction angle at point $P_i$, $L_i$ is the length of segment $P_i{P}_{i+1}$, and $R_i$ is the distance between point O and point $P_i$.

The external work rate W is contributed to by the weight of the collapse block (W_G) and the forces induced by earthquakes (W_H and W_V, respectively, for the horizontal and vertical directions), which can be defined as

$$W_G=-\gamma \cdot \omega {\displaystyle \sum S_i}\cdot R_{Gi}\cdot \cos {\theta }_{Gi}$$

(5)

$$W_H=\gamma \cdot \omega {\displaystyle \sum S_i\cdot [1+\frac{h}{H}(f-1)]\cdot k_h\cdot \sin {\omega }_e(t-\frac{h}{V_s})\cdot (y_O-y_{{\mathrm{G}}_i})}$$

(6)

$$W_V=\gamma \cdot \omega {\displaystyle \sum S_i\cdot [1+\frac{h}{H}(f-1)]\cdot k_v\cdot \sin {\omega }_e(t+t_0-\frac{h}{V_{\mathrm{p}}})\cdot (x_{{\mathrm{G}}_i}-x_O)}$$

(7)

where $\gamma$ is the unit weight and $S_i$ is the area of the discretized element. An infinitesimal trapezoidal element $P_i{P}_{i+1}{Q}_{i+1}{Q}_i$ as depicted in Figure 1b is utilized given that seismic accelerations are uniform at the same depth. $R_{Gi}$ represents the distance between the rotating center O and the $P_i{P}_{i+1}{Q}_{i+1}{Q}_i$ centroid $P_{Gi}$ and ${\theta }_{Gi}$ is the angle between $O{P}_{Gi}$ and the horizontal direction; h is the ordinate of $P_{Gi}$.

Following the determination of work rates, the shear strength is implemented to calculate safety factor $F_s$ for a given slip surface. An optimization work is given with respect to three parameters, including $r_0$, ${\theta }_0$ and t.

$$F_s=f\left(r_0,{\theta }_0,t\right)$$

(8)

An optimization strategy utilizing both a bisection approach and genetic algorithm methods is employed to determine the critical value of $F_s$. Further details can be found in Zhang et al. [12].

3. Probabilistic Analysis

This section provides a concise overview of the probabilistic methods utilized in this study, namely the metamodel PCK, MCS, and ANCOVA. Additionally, a procedural outline of the probabilistic analysis is presented.

3.1. Presentation of the Used Probabilistic Methods

3.1.1. Metamodel: PCK

The Polynomial Chaos Kriging (PCK) method integrates the strengths of Polynomial Chaos Expansion (PCE) and Kriging. It employs a sparse set of orthonormal polynomials to capture the global behavior of the computational model and utilizes Kriging to address local variability in model outputs [27]. This approach facilitates the development of an efficient and accurate metamodel compared to using PCE and Kriging separately [28]. The basic function can be expressed as

$$\mathit{Y}\approx M^{PCK}(\mathit{x})={\displaystyle \sum _{i\in A}{\alpha }_i{\mathsf{\Phi }}_i(\mathit{x})}+{\sigma }^{2}Z(\mathit{x})$$

(9)

where x is a vector with input random variables, $M^{PCK}(\mathit{x})$ is the model output approximation using PCK, and A is the index set of polynomials. ${\displaystyle \sum _{i\in A}{\alpha }_i{\mathsf{\Phi }}_i(\mathit{x})}$ is the sum of orthonormal polynomials, which is used to capture the trend within the universal Kriging framework. To ensure practical feasibility, it is necessary to truncate this sum to a finite number of terms. Standard truncation methods can be inefficient, especially for high-dimensional problems. An improved hyperbolic truncation scheme proposed by Blatman and Sudret [29] is used in this study for more efficient model approximation.

The model construction involves determining a multivariate polynomial basis ${\mathsf{\Phi }}_i(\mathit{x})$, which is the tensor product of univariate orthonormal polynomials, and the corresponding unknown coefficients ${\alpha }_i$, which are estimated by the least-square minimization method [30]. Various families of univariate orthonormal polynomials are considered, with Hermite polynomials corresponding to standard normal random variables used in this study. Additionally, Sparse Polynomial Chaos Expansion (SPCE) is implemented using the Least Angle Regression (LAR) algorithm to select important candidate polynomial bases. Insignificant PCE coefficients are then discarded to reduce the number of polynomials.

${\sigma }^2$ and $Z(\mathit{x})$ denote, respectively, the variance and the zero mean. The unit variance stationary Gaussian process is defined by an autocorrelation function between two sample points $R(\left|x-x^{\prime }\right|;\theta )$, where $\theta$ is the hyper-parameter to be estimated. Several autocorrelation functions have been proposed, i.e., linear, exponential, Gaussian, and Matern functions. The Matern-5/2 model is adopted as the autocorrelation function in this study [31].

3.1.2. Reliability and Sensitivity Methods

The limit state function for the slope stability can be expressed as

$$G(x)=F_{sc}-F_s$$

(10)

where $F_{sc}$ and $F_s$ are, respectively, the limiting safety factor and the calculated safety factor. The slope is considered a failure when $F_s$ exceeds $F_{sc}$.

The failure probability can then be determined directly based on the crude sampling method MCS [32], and can be expressed as

$$P_f=\frac{1}{N_{MCS}}{\displaystyle \sum _{k=1}^{N_{MCS}}{I}_k}$$

(11)

where $N_{MCS}$ is the size of the Monte Carlo population, and the indicator function $I_k$ is equal to 1 for G(x) < 0 and to 0 for G(x) ≥ 0. While the increased samples enhance accuracy, more computational effort is necessary. To balance computational accuracy and efficiency, the coefficient of variation for $P_f$ is introduced and can be expressed as

$$CO{V}_{P_f}=\sqrt{\frac{1-P_f}{N_{MCS}{P}_f}}$$

(12)

The value of $CO{V}_{P_f}$ should be below one specified threshold value, which is often determined as 0.05.

For sensitivity analysis, the Analysis of Covariance (ANCOVA)-based GSA, which can consider the dependence between the considered random variables, is implemented in this study to analyse the impacts of parameters on the model response. Three parts are considered for the sensitivity indices calculation, including $S_i^{\mathrm{U}}$, $S_i^{\mathrm{I}}$ and $S_i^{\mathrm{C}}$, which denote the un-correlative index of $x_i$, and the interactive and correlative index of $x_i$, respectively. $S_i^{\mathrm{U}}$ represents the independent contributions, while $S_i^{\mathrm{I}}$ and $S_i^{\mathrm{C}}$ denote the dependent contributions, which can be defined by

$$\begin{array}{c}{S}_i^{\mathrm{U}}=\frac{\mathrm{Var}[M_i^{PC}(x_i)]}{\mathrm{Var}[Y]}\\ S_i^{\mathrm{I}}=\frac{\mathrm{Cov}[M_i^{PC}(x_i),{\displaystyle {\sum }_{v\subset \left\{1,\dots ,M\right\}}{M}_v^{PC}(x_v)}]}{Var[Y]}\\ S_i^{\mathrm{C}}=\frac{\mathrm{Cov}[M_i^{PC}(x_i),{\displaystyle {\sum }_{w\subset \left\{1,\dots ,M\right\}}{M}_w^{PC}(x_w)}]}{Var[Y]}\end{array}\}$$

(13)

where $M_i^{PC}(x_i)$ denotes the univariate component function(s), and Y is the full PCE.

The sensitivity index can then be obtained by the sum of these three parts [33].

3.2. Automation of Probabilistic Analysis

The PCK metamodel-based probabilistic analysis is conducted to reduce computational burden, leveraging the efficiency of metamodel construction. Subsequent probabilistic methods (such as MCS and ANCOVA) are performed using the metamodel to derive probabilistic results. It is crucial to ensure the accuracy of the metamodel for reliable analysis. However, using a single set of experimental design (ED) may not suffice for constructing a highly accurate metamodel. Therefore, enriching the ED is necessary to enhance its representativeness and ability to capture the complexities of computationally expensive models. Detailed information about the proposed DUBLA-PDM-PCK procedure can be found in Figure 2 and is provided below for reference.

• Step 1: Preparation.
  • Construct a deterministic model for slope seismic stability;
  • Determine the statistical information of considered parameters (i.e., the distribution, mean value, coefficient of variation, and the cross-correlation coefficient);
• Step 2: Input–output sets’ determination.
  • Generate an initial experimental design N_ED based on the given statistical information and a large number of samples S_Pf used for the failure probability P_f calculation. N_ED is set to be 2N_RV [34], and N_RV is the number of considered random variables.
  • Compile the batch commands of MATLAB and map the generated samples on the deterministic models;
  • Export the model response (i.e., the safety factor) based on the DUBLA and save the results automatically.
• Step 3: PCK model construction.
  • Construct a PCK metamodel based on the input–output sets;
  • Evaluate the P_f based on the metamodel PCK.
• Step 4: Accuracy control (two stop** criteria need to be satisfied).
  • The first one is related to the leave-one-out (LOO) error of the PCK metamodel as presented by

$$Er{r}_{\mathrm{LOO}}=\frac{1}{N}\left[\frac{{\displaystyle \sum _{i=1}^N{({\mathit{Y}}^{(i)}-{\mu }_{y,(-i)}(x^{(i)}))}^2}}{\mathrm{Var}(\mathit{Y})}\right]<Er{r}_{\mathrm{LOO}\_\mathrm{c}}$$

(14)

where ${\mu }_{y,(-i)}(x^{(i)})$ is the values predicted by the PCK metamodel based on the ED except $x^{(i)}$, ${\mathit{Y}}^{(i)}$ is the DUBLA output, and $\mathrm{Var}(\mathit{Y})$ is the corresponding estimated variance [33]. The $Er{r}_{\mathrm{LOO}\_\mathrm{c}}$ is the predefined threshold value that is used for the accuracy control, which means that only the LOO value based on the existing input–output sets is smaller than the $Er{r}_{\mathrm{LOO}\_\mathrm{c}}$, and the next step can proceed.

The second one corresponds to the failure probability convergence, which can be determined by

$$Er{r}_{P_f}=\max (E_{P_f})<Er{r}_{P_f\_c}$$

(15)

where $E_{P_f}=\mathrm{abs}[(P_{fi}-P_{fi-1})/P_{fi-1}]$, and $P_{fi-1}$ and $P_{fi}$ are the P_f values for the (i − 1)th and ith iterations. i is in the range of $[N_{\mathrm{en}}-N_{\mathrm{tg}}+1,N_{\mathrm{en}}]$, where $N_{\mathrm{en}}$ is the enrichment samples number, and $N_{\mathrm{tg}}$ is the required comparative failure probabilities. $Er{r}_{P_f\_c}$ is the threshold value for the failure probability error. Both threshold values (i.e., the $Er{r}_{P_f\_c}$ and $Er{r}_{P_f\_c}$) are set as 0.01 to ensure the metamodel’s accuracy.

Metamodel construction cannot be finished until the criteria are met. Conversely, the informative samples from the S_Pf need to be selected to enrich the ED for the metamodel’s construction. The learning function used in this study is U function and can be defined as

$$U(x)=\frac{\left|\mu (x)\right|}{\sigma (x)}$$

(16)

The newly added sample is chosen by $S_n=\arg \min U(x)$ [34,35].

• Step 5: Post-processing of MCS results.
  • Perform MCS and ANCOVA based on the well-built PCK metamodel;
  • Perform post-processing of the probabilistic results, including the failure probability, the statistical moments, the PDF/CDF of the model response, and the sensitivity indices.

4. Comparison and Validation

This section presents the probabilistic analysis of slope seismic stability based on the proposed DUBLA-PDM-PCK methodology to indicate its efficiency and accuracy. A slope exhibiting linear variation of shear strength parameters with depth is considered, followed by probabilistic analyses and comparison of results.

4.1. Problem Statement

A non-homogeneity case presented in Figure 3, which assumes that the soil friction angle and cohesion vary linearly with depth (referred to as the depth-dependent case), is discussed.

${\phi }_1$ and $c_1$ represent the friction angle and cohesion at h₁. ${\phi }_2$ and $c_2$ denote the friction angle and cohesion at a depth of h₂. In this context, h₁ and h₂ are set to 0 and H, which means that ${\phi }_1(c_1)$ and ${\phi }_2(c_2)$ are the friction angle (cohesion) at the slope crest and H below the slope crest, respectively. ${\lambda }_{\phi }$ and ${\lambda }_c$ are the dimensionless non-homogeneity coefficients and can be defined as

$$\{\begin{array}{l}{\lambda }_{\phi }={\phi }_2/{\phi }_1\\ {\lambda }_c=c_2/c_1\end{array}$$

(17)

The slope soil is homogeneous when ${\lambda }_{\phi }=1$ and ${\lambda }_c=1$.

In probabilistic analysis, the dimensionless coefficients of cohesion and friction angle, the cohesion and friction angle at slope crest, and the unit weight are considered as random variables. In addition, the uncertainties of the horizontal and vertical seismic coefficients, seismic period, and the amplification factor are accounted for.

Table 1. Parameters of the reference case.

Input Parameters		Notation	Statistics of Parameter		Correlation Coefficient	Unit
			Mean	COV
Pseudo-dynamic	Horizontal seismic coefficient	kh	0.2	0.25	-	-
	Vertical seismic coefficient	kv	0.1	0.25	-	-
	Period	T	0.3	0.2	-	s
	Amplification factor	f	1.2	0.2	-	-
Soil	Cohesion of slope crest	c2	20	0.15	−0.5	kPa
	Friction angle of slope crest	φ2	30	0.15		°
	Dimensionless coefficient of cohesion	λc	1.3	0.15	-	-
	Dimensionless coefficient of friction angle	λφ	1.3	0.15	-	-
	Unit weight	γ	20	0.1	-	kN/m3
Slope geometry	Angle	β	45	-	-	°
	Height	H	10	-	-	m

summarizes the detailed parameter information for the reference case [1,36,37].

4.2. Probabilistic Results for the Reference Case

In this subsection, a comprehensive probabilistic analysis for the reference case outlined in Section 4.1 is provided using the DUBLA-PDM-PCK procedure.

Initially, an ED comprising 18 samples is generated, which is followed by the initial metamodel PCK. Subsequently, an iterative process is initiated involving the enrichment of the ED and refinement of the PCK metamodel until specific criteria (the LOO and P_f convergence) are met. The iterative process, illustrated in Figure 4, demonstrates that the leave-one-out error $Er{r}_{\mathrm{LOO}}$ is satisfied (as presented in Equation (14)) after the enrichment of five samples. The failure probability is converged from the 62nd sample and the following 10 $E_{P_f}$ satisfy the requirements of Equation (15). Therefore, a total of 90 evaluations, including 18 initial EDs and 72 enrichments, are employed to construct the PCK metamodel. Subsequently, the MCS and ANCOVA are conducted to derive probabilistic results.

Following the MCS analysis, the failure probability, probability density functions (PDFs), and cumulative distribution functions (CDFs) of the calculated safety factors, and corresponding statistical moments, are determined. The results are illustrated in Figure 5. The analysis reveals a failure probability of 1.65 × 10⁻² with safety factors predominantly distributed around 1.4. The CDF provides a clear representation of the failure probability based on the safety factor distribution.

The sensitivity analysis (as presented in Figure 6) reveals that the ANCOVA indices of the friction angle and its non-homogeneity coefficient, as well as the horizontal seismic coefficient, are relatively significant compared to the cohesion and its non-homogeneity coefficient. Additionally, for the seismic parameters, the amplification factor is identified as important, whereas the vertical seismic coefficient and seismic period exhibit minor influences on slope seismic stability. These findings underscore the importance of considering soil non-homogeneity and spatial variations in seismic loading for accurate slope stability analysis.

4.3. Comparison and Validation

To validate the accuracy and precision of the developed DUBLA-PDM-PCK, some comparisons with the existing studies are carried out and can be found below.

4.3.1. With MCS

The accuracy of the metamodel is validated by comparing safety factors (from 1000 samples generated by LHS) evaluated using the PCK metamodel against those obtained from the analytical method DUBLA. The results are illustrated in Figure 7. The analysis reveals that nearly all scatter points fall within ±5% of the diagonal line, indicating a close correspondence between safety factors derived from the PCK metamodel and those from the DUBLA method. This validation confirms the accuracy of the PCK metamodel.

However, it should be noted that conducting MCS calculations requires at least 24,000 samples when the coefficient of variation of the failure probability COV_Pf is less than 5%, which can be time-consuming. Additionally, sensitivity analysis typically necessitates around 9000 deterministic simulations. In contrast, the proposed DUBLA-PDM-PCK procedure achieves probabilistic results using just 90 samples, significantly reducing computational effort.

4.3.2. With Kriging and SPCE

Metamodel-based probabilistic methods like Kriging and SPCE-based MCS are also employed in this study for slope seismic probabilistic analysis, aiming to assess the performance of the proposed DUBLA-PDM-PCK methodology. A detailed comparison is provided in

Table 2. Comparison of the failure probability based on three methods.

Metamodel	Pf	MeanFS	Std.FS	Ntotal
PCK	1.65 × 10−2	1.405	0.189	90
Kriging	1.72 × 10−2	1.384	0.186	110
SPCE	1.83 × 10−2	1.398	0.194	113

and illustrated in Figure 8 to elucidate the relative effectiveness of different approaches.

These three methods demonstrate similar probabilistic results, including failure probability, mean, and standard deviation of the safety factor (as presented in Table 2), along with sensitivity indices (depicted in Figure 8). This consistency validates the accuracy of the DUBLA-PDM-PCK methodology. However, it is noteworthy that the PCK method only requires 90 simulations for metamodel construction, whereas the Kriging and SPCE methods necessitate more (110 and 113 simulations, respectively), indicating greater computational effort with the latter approaches.

Furthermore, the performance of metamodels in predicting safety factors is discussed and summarized using a graphical Taylor diagram (Figure 9), which considers standard deviation (represented by black circles), root mean square error (RMSE, indicated by green circles), and correlation coefficient (depicted by blue lines). Using 1000 samples, as mentioned in Figure 7, the Kriging and PCE metamodels predict the corresponding values for evaluation. The RMSE and correlation coefficient of observed values (DUBLA) are 0 and 1, respectively. The metamodel showing the highest similarity to the observed results is considered the most effective. The PCK method demonstrates superior performance, being closest to the observed point with corresponding indicator values (correlation coefficient, RMSE, standard deviation) of 0.986, 0.024, and 0.189, respectively.

In summary, the proposed DUBLA-PDM-PCK procedure offers cost-effective probabilistic analysis with accurate results, making it highly recommended for slope seismic stability assessments.

5. Discussion

This section delves into the effects of soil non-homogeneity, correlation coefficients, and seismic forces on slope seismic stability. Various results are presented for discussion, encompassing failure probability, distributions, statistical moments of system response, and sensitivity indices. Unless otherwise specified, parameters align with those detailed in the reference case outlined in Table 1. Further details are provided below.

5.1. Influence of the Soil Non-Homogeneity

The non-homogeneity coefficient is often assumed to be a fixed value, yet its actual variability due to spatial variations and soil sedimentation necessitates consideration. Therefore, it is crucial to explore the effects of uncertainty associated with the non-homogeneity coefficient. Different coefficients of variation of the non-homogeneity coefficient for the friction angle ${\lambda }_{\phi }$ and the cohesion ${\lambda }_c$ on the probabilistic results are discussed. Figure 10 illustrates these results, showing that as the COV values of the non-homogeneity coefficient increase, the failure probability also increases, with corresponding sensitivity indices exhibiting significant growth. For instance, the sensitivity index can reach up to 0.4 when the COV value is 0.25. Therefore the soil non-homogeneity cannot be ignored for the slope stability analysis, particularly given the substantial uncertainty associated with the non-homogeneity coefficient.

5.2. Influence of Cross-Correlation Coefficients

In this section, the correlation coefficient is examined to assess its impact on slope seismic stability; specifically, how the correlation between cohesion and the friction angle influences outcomes. The results are detailed in Figure 11. The analysis reveals that as the correlation coefficient increases, the failure probability also rises. This relationship is attributed to lower cohesion and higher friction angle combinations associated with smaller correlation coefficients. Conversely, larger correlation coefficients result in more extreme scenarios where cohesion and friction angle values are simultaneously low or high, leading to notably varied safety factors.

This relationship is further illustrated in Figure 12, which shows the PDF of the safety factor under different correlation coefficients. The wider PDF curves observed with increasing correlation coefficients contribute to higher failure probabilities.

Additionally, Figure 11 presents the sensitivity indices of strength-related parameters (cohesion and friction angle at the slope crest) and their corresponding non-homogeneity coefficients. Observations reveal that with an increase in the correlation coefficient, the sensitivity indices of cohesion-related parameters increase, while those related to the friction angle decrease. When the correlation coefficient exceeds 0.5, cohesion sensitivity indices can surpass those of friction angles. Therefore, the correlation coefficient plays a significant role in the slope probabilistic stability analysis and should be determined with attention. In addition, the sensitivity indices of the non-homogeneity coefficients are decreased with the increase in the correlation coefficient, while the strength parameters at the slope crest are more significant.

5.3. Influence of the Seismic Effects

The seismic coefficient represents earthquake intensity, and a range of [0.1, 0.3] with an interval of 0.05 is examined, as illustrated in Figure 13. There is no doubt that with the increase in the horizontal seismic coefficient, the failure probability increases significantly. The observed failure probabilities range from 3.71 × 10⁻⁴ to 1.31 × 10⁻¹, corresponding to “good” and “unsatisfactory” performance levels in practical terms. [38]. Therefore, determining the horizontal seismic coefficient is crucial for slope stability analysis, given its substantial impact on failure probability.

Furthermore, Figure 13 presents the sensitivity indices of seismic-related parameters. The analysis reveals that the sensitivity index of the horizontal seismic coefficient increases with higher k_h values. Conversely, the vertical seismic coefficient shows only a slight effect on slope stability. Notably, the amplification factor demonstrates greater significance compared to seismic duration, underscoring the importance of spatial variations in seismic loading for slope stability.

6. Conclusions

An efficient solution for probabilistic slope seismic stability analysis is provided in this study, and the main contributions and conclusion can be summarized as follows:

(1)

A comprehensive method, DUBLA-PDM-PCK, is proposed to investigate slope seismic stability by combining the DUBLA analytical method with PCK-metamodel-based MCS and ANCOVA (PCK-MA). This approach accounts for soil non-homogeneity, time and space variations in seismic loading, and uncertainties in parameters such as soil strength and seismic characteristics.

(2)

DUBLA-PDM-PCK enables affordable probabilistic analysis with high accuracy and fewer deterministic simulations compared to other methods like crude Monte Carlo Simulation, Kriging, and SPCE-based probabilistic approaches. Additionally, sensitivity indices are provided to assess the importance of considered random variables, offering valuable insights for designers.

(3)

The effects of soil non-homogeneity, cross-correlation coefficient, and seismic factors on slope probabilistic stability are discussed. The results highlight the importance of considering soil non-homogeneity and spatial variations in seismic loading for slope seismic stability analysis, which means that the parameters’ determination should be given attention.

(4)

As the cross-correlation coefficient of cohesion and friction angle increases, the slope stability decreases and the effects of the friction angle-related parameters on the safety factor decreases. Moreover, the effects of the non-homogeneity coefficient uncertainty decrease with the increase in the cross-correlation coefficient.

(5)

The horizontal seismic coefficient and amplification factor significantly impact results, with the former showing pronounced effects as it increases. In contrast, the vertical seismic coefficient and seismic period have less influence on outcomes.

This paper provides valuable insights into deterministic and probabilistic analyses of slope stability under seismic loads. However, certain limitations persist, including the treatment of shear strength parameters (cohesion and friction angle) as random variables for simplicity, despite their spatial variability in practice. The spatial variability will be addressed in future research endeavors.