Reliability analysis method for soil slopes permanent displacement under mainshock–aftershock sequences

After a primary seismic event, subsequent aftershocks frequently induce progressive damage to slopes. Evaluating the response of slopes to mainshock–aftershock sequences (MAS) from a probabilistic perspective is crucial for disaster prevention and mitigation. Current research primarily focuses on single mainshock events, and commonly adopts peak ground acceleration (PGA), with limited consideration of the cumulative aftershock effects. This study proposes a PDEM-based reliability framework for soil slopes subjected to MAS. First, the random input field of the MAS is first constructed by integrating theoretical models with real data. Then, considering the peak, cumulative, and spectral characteristics of the MAS, correlation analysis is conducted to identify cumulative absolute velocity (CAV) as the controlling parameter for the soil slope response among the 21 candidate parameters. Finally, based on the probability density evolution method (PDEM), a reliability assessment framework for soil slope behavior under MAS is constructed. Compared with existing methods, the proposed approach more effectively incorporates the effects of aftershocks and enables more accurate reliability assessment of slope permanent displacement under MAS. This study provides a probabilistic framework and methodological approaches for assessing slope stability under MAS.

1 Introduction

Major seismic events are commonly followed by a series of aftershocks, whose effects on slope stability should not be underestimated. Observations of earthquake-induced damage indicate that aftershocks often act as the final trigger for slope failure (Xu and Huang, 2008; Yin, 2008). As the safety requirements for engineering structures become increasingly stringent, slope stability under MAS has attracted considerable interest from both researchers and engineering practitioners. Li et al. (2009) developed and improved the PDEM framework for analyzing stochastic dynamic systems, providing a systematic approach to nonlinear stochastic dynamic analysis and reliability assessment of large-scale and complex engineering structures. The framework has been widely applied to stochastic dynamic analyses of slopes (Li and Chen, 2017). Building on this work, Pang et al. (2021, 2024) introduced an enhanced generalized-PDEM approach for reliability analysis of complex slopes. Their methodology considers multiple slope parameters and ground motion uncertainties. Additionally, accounting for the spatial variation in soil strength properties, a reliability assessment methodology integrating the Newmark method and PDEM is developed to evaluate the effects of aftershocks and spatial heterogeneity on slope dynamic reliability (Newmark, 1965; Zhou et al., 2023; Wang et al., 2022; Xu et al., 2025). Most existing studies use peak ground acceleration (PGA) as the primary measure of ground motion intensity. Nevertheless, intensity measures based only on peak values may not adequately capture the overall characteristics of MAS ground motions (Ruiz-Garcia and Negrete-Manriquez, 2011; Amiri et al., 2022).

For characterizing the ground motion random field of the MAS, the Monte Carlo method is a well-established and effective approach (Metropolis and Ulam, 1949; Hu et al., 2018; Nithin et al., 2020; Kim and Sitar, 2013). However, it requires a large number of samples and incomplete probability information of the sample set (Jiang et al., 2021). The random function-dimension reduction simulation technique generates MAS time histories with associated probabilities, thereby forming a complete probabilistic set. These time histories can then be coupled with the PDEM to conduct sophisticated dynamic response and reliability analyses for complex engineering structures subjected to MAS (Liu and Liu, 2017; Liu et al., 2019, 2022).

Current reliability analysis methods for seismic slopes are predominantly developed for single-mainshock scenarios and therefore may not adequately assess slope reliability under MAS. Moreover, existing studies on sequential ground motions predominantly use PGA as the intensity measure, which may not adequately capture the cumulative damage effect contributed by aftershocks. To address these limitations, this study adopts cumulative absolute velocity (CAV) as the intensity measure for MAS ground motions and develops a reliability analysis framework for the permanent displacement of soil slopes based on the PDEM. The proposed method accounts for the effects of aftershocks and provides a basis for stability assessment and disaster prevention of soil slopes subjected to MAS.

Figure 1 Geometric dimensions and mesh division of the slope model.

Table 1 Parameters of the slope.

2 Methods

2.1 Model

This investigation employs a two-layer slope model to investigate slope response under MAS. Figure 1 shows the double-layer slope model established by referring to Wang et al. (2021). The model is 150 m long and 30 m high. The grid size ΔL range from 0.5 to 1.8 m, with a finer resolution of approximately 0.5 m in the potential sliding zone. The soil consists of clay from southwestern China (Yang et al., 2022; Ma et al., 2023; Zhang et al., 2024), and its properties are listed in Table 1.

This investigation utilizes the finite-difference platform FLAC3D, with soil mechanical behavior characterized through a Mohr–Coulomb constitutive model with a tensile cutoff. Local damping is adopted in this study. Its implementation is given by Eq. (1) (Li and Yang, 2006). Because FLAC3D does not require the assembly of a global stiffness matrix, it is computationally efficient for nonlinear dynamic analysis. This feature makes FLAC3D suitable for simulating the nonlinear dynamic response of slopes. (Hu et al., 2017; Puthanpurayil et al., 2018; Yan et al., 2011).

$$\begin{array}{}\text{(1)}& \mathit{\xi }=\mathit{\pi }{D}_{\mathrm{s}}\end{array}$$

Within Eq. (1), the coefficient ξ quantifies localized damping intensity, D_s designates the critical damping fraction, while π maintains its conventional value of 3.14.

This investigation adopts a 5 % critical damping fraction to replicate energy dissipation characteristics during seismic wave transmission through soil media (Qu et al., 2015), yielding a local damping coefficient of 0.157 according to Eq. (1).

Figure 2 Correlation coefficient of MAS ground motion parameters (PGA: peak ground acceleration; PGV: peak ground velocity; PGD: peak ground displacement; a_sq: squared acceleration; v_sq: squared velocity; d_sq: squared displacement; a_rs: root square acceleration; v_rs: root square velocity; d_rs: root square displacement; I_A: Arias intensity; t_d: significant duration; CAV: cumulative absolute velocity; CAD: cumulative absolute displacement; CAI: cumulative absolute impulse; a_rms: root-mean-square acceleration; v_rms: root-mean-square velocity; d_rms: root-mean-square displacement; I_c: characteristic intensity; Sa_Ts: response spectral acceleration at natural vibration period T_s; Sa_1.5Ts: response spectral acceleration at 1.5 T_s; ASI: acceleration spectral intensity.).

2.2 Reliability analysis method of slope under MAS

Based on the principle of probability conservation, generalized density evolution equations can be derived for the analysis of stochastic dynamic systems. By combining these equations with the virtual stochastic process method, the probability characteristics of stochastic systems can be modeled. Furthermore, based on the PDEM, the reliability analysis of the permanent displacement of the slope under the MAS action is realized (Li and Chen, 2008; Liu and Liu, 2017; Jiang et al., 2021). The main steps are as follows:

• 1.

  The measured data of the MAS were statistically analyzed to obtain the frequency-domain energy distribution function curves for each ground motion. Employing optimal square approximation criteria, we determine the parameter vector λ defining each evolving power spectral density function. The MAS dataset in this research is modeled as zero-mean, fully non-stationary processes characterized by an evolutionary power spectral density (EPSD) model. The discrete representative points of each ground motion are expressed as ${\mathit{\theta }}_{\mathrm{q}}=({\mathit{\theta }}_{\mathrm{1},\mathrm{q}},{\mathit{\theta }}_{\mathrm{2},\mathrm{q}},\mathrm{\dots },{\mathit{\theta }}_{s,\mathrm{q}}$, q = $\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{96}$), with the corresponding probabilities given in Eq. (2):

  $$\begin{array}{}\text{(2)}& P_q={\int }_{Vq}p\mathrm{\Theta }\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\end{array}$$
• 2.

  Within the framework of reliability analysis based on PDEM, it is essential to identify a ground-motion parameter that exhibits the strongest correlation with slope displacement responses under MAS. We have demonstrated that among the 21 parameters considered in three categories – namely, peak types (PGA, PGV, PGD), spectral characteristics (Sa_Ts, Sa_1.5Ts, ASI), and cumulative types (a_sq, v_sq, d_sq, a_rs, v_rs, d_rs, I_A, CAV, CAD, CAI, a_rms, v_rms, d_rms, I_c, t_d) – CAV showed the highest correlation coefficient (ρ = 0.872) with the permanent displacement, and more effectively captured the characteristics of aftershocks, as shown in Fig. 2 (Zhang et al., 2024). Consequently, CAV is adopted as the primary intensity measure for slope displacement response under MAS and combined with the PDEM for reliability assessment.

  Based on trial calculations, 96 sets of MAS were scaled to three CAV levels: 12, 28, and 40 m s⁻¹. Using the FLAC3D, the slope response under MAS, denoted as ${\dot{X}}_j({\mathrm{\Theta }}_q,t)$; j = $\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{96}$.

• 3.

  Based on the above two steps to get assigned to probability P_q and slope of the response results ${\dot{X}}_j({\mathrm{\Theta }}_{\mathrm{q}},t)$ probability density evolution Eq. (3):

  $$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\displaystyle \frac{\partial p_{X\mathrm{\Theta }}(x,\mathit{\theta },t)}{\partial t}}+\dot{{\sum }_{j=\mathrm{1}}^m\dot{X_j}}({\mathrm{\Theta }}_{\mathrm{q}},t)\times {\displaystyle \frac{\partial p_{X\mathrm{\Theta }}(x,\mathit{\theta },t)}{\partial x}}=\mathrm{0},\\ & q=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }\mathrm{\dots },\mathrm{96}\end{array}\end{array}$$

  Within this formulation, X designates system response quantities, Θ encompasses the ensemble of random variables exclusive of initial conditions. $\mathrm{\Theta }=(X_{\mathrm{0}},\dot{\mathrm{\Theta }})$ represents the vector comprising all stochastic variables within the system; X=H(Θt) indicates the physical response of X; p_XΘ(xθt) signifies the joint probability density function (PDF) for X and Θ.

  When t=t₀, the following initial condition is as shown in Eq. (3):

  $$\begin{array}{}\text{(4)}& p_{X\mathrm{\Theta }}(x,\mathit{\theta },t){\mathrm{|}}_{t=t_{\mathrm{0}}}=\mathit{\delta }(x-x_{\mathrm{0}})P_{\mathrm{\Theta }}\left(\mathit{\theta }\right)\end{array}$$

  In this expression, δ denotes the Dirac delta function, while x₀ corresponds to the X value t=t₀. The PDF of X at any given time can be derived by summing the contributions from all discrete numerical solutions, as expressed in Eq. (4):

  $$\begin{array}{}\text{(5)}& p_X(x,t)=\int p_{X\mathrm{\Theta }}(x,\mathit{\theta },t)\mathrm{d}\mathit{\theta }\end{array}$$

  Probability density evolution equations require numerical methods for their solution. Deriving a direct analytical solution to Eq. (4) is mathematically challenging. Conventional practice uses finite-difference schemes to solve the equation.

  $$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial p(z,t)}{\partial t}}+a\left(t\right){\displaystyle \frac{\partial p(z,t)}{\partial z}}=\mathrm{0}\end{array}$$

  First-order Taylor series expansion of Eq. (5), yields approximate differential expressions:

  $$\begin{array}{}\text{(7)}& {\left[{\displaystyle \frac{\partial p}{\partial t}}\right]}_j\left(k\right)={\displaystyle \frac{p_j^{(k+\mathrm{1})}-p_j^{\left(k\right)}}{\mathrm{\Delta }t}}\end{array}$$

  For spatial coordinate z, analogous treatment produces:

  $$\begin{array}{}\text{(8)}& {\left[{\displaystyle \frac{\partial p}{\partial z}}\right]}_j\left(k\right)={\displaystyle \frac{p_j^{(k+\mathrm{1})}-p_j^{\left(k\right)}}{\mathrm{\Delta }z}}\end{array}$$

  Substitute Eqs. (6) and (7) into Eq. (5):

  $$\begin{array}{}\text{(9)}& \begin{array}{rl}{p}_j^{(k+\mathrm{1})}& =p_j^{\left(k\right)}-\mathit{\lambda }{a}^{\left(k\right)}\left[p_j^{\left(k\right)}-p_{j-\mathrm{1}}^{\left(k\right)}\right]\text{ or }{p}_j^{(k+\mathrm{1})}\\ & =(\mathrm{1}-\mathit{\lambda }{a}^{\left(k\right)})p_j^{\left(k\right)}+\mathit{\lambda }{a}^{\left(k\right)}{p}_{j-\mathrm{1}}^{\left(k\right)}\end{array}\end{array}$$

  Equation (9) is the one-sided difference scheme. Substitute ${\dot{X}}_j({\mathrm{\Theta }}_{\mathrm{q}},t)$ obtained from ${\dot{X}}_j$ into Eq. (3), and solve the partial differential equation using the L-W or TVD format of the finite difference method.

• 4.

  The PDF of slope displacement response p_X(xt) is obtained numerically through cumulative summation of the joint PDF p_XΘ(xθt); $q=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{96}$. Subsequently, both PDF and cumulative distribution function (CDF) graphs for slope displacement are developed. These PDF and CDF curves are then used to characterize the statistical distribution of slope displacement and evaluate slope reliability.

3 Stochastic simulations of MAS

Empirical MAS recordings provide more realistic representations of structural damage and cumulative degradation (Yang et al., 2022). This study uses a diverse set of real earthquake records to reduce the bias associated with specific ground motion types. Through iterative parameter identification based on evolutionary power spectral density modeling, this study combines theoretical models with data-driven methods to construct stochastic seismic inputs. These inputs provide a stochastic loading basis for investigating the effects of MAS on slope reliability.

This study selected 96 recorded mainshock–aftershock pairs from 15 seismic events in the NGA-West2 strong-motion database maintained by the Pacific Earthquake Engineering Research Center (PEER). We selected records based on five criteria (Bray and Macedo, 2019; Yeznabad et al., 2022; Yeznabad et al., 2026): (1) the mainshock and its aftershock must belong to the same seismic event, with the aftershock defined as the largest-magnitude event occurring within 12 months after the mainshock; (2) both events have a M_w ≥ 5.0 and a PGA ≥ 0.05 g; (3) both records are taken from the same component of the same station; (4) the station has an average shear-wave velocity in the upper 30 m (Vs30) between 100 and 700 m s⁻¹; and (5) the rupture distance of each event does not exceed 80 km. The 96 MAS records selected according to the above criteria encompass a broad range of sequence-type ground-motion characteristics. However, as these records are derived from the NGA-West2 database and primarily represent shallow crustal earthquakes, the applicability of the results is primarily limited to MAS ground motions from shallow crustal earthquakes.

Figure 3 Typical acceleration time history of a MAS.

Figure 4 The identification effect of the parameter vector λ. (a) The fitting results of the mainshock. (b) The fitting results of the aftershock.

Each MAS adopted a “mainshock + 20 s temporal gap + aftershock” configuration to preserve comprehensive non-stationary sequence characteristics (Wang et al., 2022). Figure 3 displays a typical mainshock–aftershock acceleration record from the Whittier Narrows event, with the mainshock and aftershock numbered 589 and 707, respectively, in the NGA-West2 database. In the sensitivity analysis of the 21 MAS ground-motion parameters, the original ground-motion records were used without amplitude scaling. Subsequently, when calculating the permanent displacement of the slope at different CAV levels, the ground-motion amplitudes were scaled to better represent practical engineering conditions.

Table 2 Identification results of MAS.

For the mainshock–aftershock sequences constructed above, frequency-domain energy distribution analysis is performed using an evolutionary power spectrum model for fully non-stationary ground motion processes (Priestley, 1965; Liu and Liu, 2017). The best square approximation criterion and the least square method are adopted to fit the energy distribution curve in the frequency domain. The parameter vector λ is further obtained through inversion and regression. In the process, considering the site damping ratio ξ_f excellence, site soil circular frequency ω_f and the influence of a frequency modulation function parameter. The identification results are shown in Fig. 4, and the corresponding data are listed in Table 2. The constructed ground motions can be regarded as zero-mean, real-valued, non-stationary ground motion processes generated by the evolutionary power spectral density function model.

Figure 5 Displacement response distribution of slopes under MAS. The horizontal axis (Ground motioN ID) represents the serial number of ground motion records (ID = 1–96).

4 Results and analysis

To assess slope reliability based on permanent displacement responses, it is necessary to define critical displacement thresholds. Jibson and Michael (2009) classified slope displacement into four ranges for seismic landslide risk evaluation: displacements of 0–1 cm indicate low risk, 1–5 cm indicate medium risk, 5–15 cm indicate high risk, and those exceeding 15 cm indicate very high risk. For permanent displacement, Ozkan (1998) suggests that the control standard for seismic sliding deformation is 1 m. In this study, based on the actual engineering requirements, three slope displacement thresholds of 0.05, 0.25, and 0.50 m were selected, corresponding to low-, medium-, and high-level failure states, respectively, to systematically evaluate slope reliability under MAS.

Figure 6 Upper bounds of permanent displacement at different CAV and PGA levels.

Figure 5 shows the distribution of slope permanent displacement under MAS. The distribution is widely dispersed, and the slope displacement response increases with increasing ground motion intensity. When the CAV of the MAS was 12, 28 and 40 m s⁻¹, respectively, the corresponding average slope permanent displacements were 0.067, 0.328 and 0.633 m, respectively, which are comparable to the slope displacement results when the PGA of the MAS was 0.4, 0.5 and 0.6 g in reference (Wang et al., 2022). Figure 6 compares the upper bounds of the permanent displacement ranges obtained at three corresponding CAV and PGA levels. It can be observed that, compared with PGA, the use of CAV leads to a more concentrated distribution of permanent displacement results, with the overall dispersion reduced by approximately 40 % on average. These results demonstrate that CAV provides a more concentrated representation of permanent displacement responses. Therefore, CAV is more suitable than the commonly used PGA for characterizing MAS ground motions in the reliability analysis of slope permanent displacement.

Nevertheless, even with the use of CAV as the MAS ground-motion intensity measure, the obtained slope displacement responses remain somewhat scattered. This suggests that additional ground-motion parameters – including frequency content and duration – also influence slope permanent displacement (Rathje and Saygili, 2008, 2009; Yeznabad et al.,2022). Given the current lack of a single comprehensive ground-motion parameter that can effectively predict slope displacement responses, the PDEM framework propagates this inherent variability into the reliability estimate, thus rendering a probabilistic treatment indispensable. Building on previous studies, this study develops a reliability probabilistic analysis method based on the PDEM.

Figure 7 Statistical analysis of slope permanent displacement probability under MAS. Probability Density Function (PDF), describing the probability density of permanent displacement, and Cumulative Distribution Function (CDF), representing the cumulative probability of permanent displacement.

Figure 7 presents the PDF and CDF curves of slope permanent displacement. The PDF curves show the distribution characteristics and main concentration ranges of slope permanent displacement, while the CDF curves provide the cumulative probability at different displacement levels for evaluating slope reliability based on the specified displacement thresholds. The PDF curves shown in Fig. 7a exhibit bimodal or multimodal characteristics, making them difficult to describe using standard distributions such as the normal or lognormal distribution. Under the MAS action of different intensities, the displacements of slopes have diverse distribution patterns, and their probabilities are also not the same. When the CAV is 12 m s⁻¹, slope permanent displacement ranges from 0 to 0.25 m, and the slope reliability is 48.5 % at the displacement threshold of 0.05 m corresponding to the low failure state. At a CAV of 28 m s⁻¹, the displacement range expands to 0–1.00 m, and the slope reliability is 45.5 % at the displacement threshold of 0.25 m corresponding to the medium failure state. At a CAV of 40 m s⁻¹, the displacement range further expands to 0–2.00 m, and the slope reliability decreases to 31.7 % at the displacement threshold of 0.50 m corresponding to the high failure state.

Figure 8 Slope permanent displacement at a CDF of 50 %.

Figures 8 and 9 compares the CDF of slope permanent displacement obtained using CAV and PGA. The CAV-based CDF curves show an overall leftward shift relative to the corresponding PGA-based curves, indicating smaller displacement values at the same cumulative probability. At a cumulative probability of 50 %, the slope permanent displacements are approximately 0.07, 0.32, and 0.55 m for CAV = 12, 28, and 40 m s⁻¹, respectively, compared with approximately 0.20, 0.50, and 0.70 m for PGA = 0.4, 0.5, and 0.6 g. The difference becomes more evident at higher cumulative probabilities. The PGA-based CDF exhibit more pronounced large-displacement tails and greater dispersion, whereas the CAV-based CDF show more concentrated displacement distributions. Overall, CAV better captures the cumulative effects of MAS on slope permanent displacement and therefore provides a more comprehensive intensity measure for slope reliability analysis under MAS.

Figure 9 Compared with the reliability results of slope displacement response in reference (Wang et al., 2022).

5 Conclusions

By identifying the intensity measure of MAS and integrating it with the PDEM, this study develops a reliability analysis method for slope permanent displacement under MAS. The main conclusions are as follows:

• 1.

  Among the 21 MAS ground-motion parameters considered, CAV showed the strongest correlation with slope permanent displacement (ρ = 0.872) and best represents the cumulative effects of MAS. Therefore, CAV can be used as the intensity measure to assess the response of soil slopes under MAS.

• 2.

  By incorporating the dynamic calculation results, the slope permanent displacement gradually increases as CAV increases from 12 to 40 m s⁻¹, while the displacement results still exhibit a degree of scatter. This variability highlights the need for probabilistic analysis when evaluating slope reliability under MAS.

• 3.

  By incorporating the characteristics of MAS into the PDEM framework and adopting CAV as the intensity measure, a methodology for evaluating the reliability of slope permanent displacement under MAS is developed. Compared with the conventional PGA-based method, the proposed method effectively reduces the dispersion of slope permanent displacement responses and improves the estimation of slope reliability at different displacement thresholds, demonstrating its theoretical and practical value.

This study focuses on clayey soil slopes in southwestern China and examines the feasibility of the proposed reliability analysis method for soil slope responses under shallow-crustal MAS. The effects of soil-parameter variability and slope geometric variability are not considered. The principal contribution of this study lies in the proposed methodological framework, which integrates CAV as the intensity measure for MAS with the PDEM for efficient reliability analysis. Further consideration of soil-parameter variability is expected to improve the applicability of the proposed method to practical engineering problems.