Dynamic Reliability Analysis of Layered Slope Considering Soil Spatial Variability Subjected to Mainshock–Aftershock Sequence

: The slope instability brought on by earthquakes frequently results in signiﬁcant property damage and casualties. At present, the research on displacement response of a slope under earthquake has mainly emphasized the action of the mainshock, without accounting for the impact of an aftershock, and the spatial variability of material parameters is often neglected. The spatial variability of parameters is fully accounted for in this paper, and dynamic reliability of permanent displacement ( D P ) of a slope produced by the mainshock–aftershock sequence (MAS) is studied. A slope reliability analysis method is proposed based on the Newmark displacement method and the generalized probability density evolution method (GPDEM) to quantify the effect of the spatial variability of materials parameters on dynamic reliability. Firstly, the parameter random ﬁeld is generated based on the spectral representation method, and the randomly generated parameters are assigned to the ﬁnite element model (FEM). In addition, the random simulation method of MAS considering the correlation between aftershock and mainshock is adopted based on the Copula function to generate the MAS. Then, the D P of slopes caused by the MAS considering the spatial variability is calculated based on the Newmark method. The impacts of the coefﬁcient of variation ( COV ) and aftershock on the D P of slope is analyzed by means of mean values. Finally, the effect of COV and aftershock on the reliability of D P is explained from a probabilistic point of view based on the GPDEM. The results revealed that with the increase in the COV , the mean of the D P of the slope shows a trend of increasing gradually. The D P of slope is more sensitive to the coefﬁcient of variation of friction angle ( COV F ). The mean D P of the slope induced by the MAS is larger compared to the single mainshock, and the PGA has a signiﬁcant impact on the D P .


Introduction
The instability of slopes brought on by earthquakes is a significant geological risk. Strong earthquakes have a significant impact on large-scale geological disasters, such as landslides and debris flows brought on by slope instability, which frequently result in catastrophic losses and negative social repercussions [1]. It is reported that the Chi-Chi earthquake initiated in excess of 10,000 landslides and slope instability in 11,000 m 2 in Central Taiwan [2]. Around 20,000 people perished in the 2008 Wenchuan earthquake as a result of a large number of landslides and slope instability issues, which made up nearly half of all earthquake fatalities [3,4]. In the Yushu earthquake in 2010, the earthquake created more than 2000 landslides, resulting in a direct economic loss of about CNY 600,000, 8 deaths, and 14 injuries [5]. These aforementioned disaster consequences show that reasonable consideration of dynamic response and sliding displacement of slopes induced by strong earthquakes is very essential for predicting the potential damage possibility of ground motion and conducting rapid seismic risk assessment.
A large amount of historical seismic data has shown the occurrence of strong earthquakes is frequently complemented by multiple aftershocks [6,7]. In the two months after the Wenchuan earthquake in 2008, more than 20,000 aftershocks were triggered, including dozens of strong aftershocks of magnitude 5 or greater. In the three days after the 2013 Lushan earthquake, there were more than 3000 aftershocks and about 4 aftershocks with magnitude stronger than 5 [8]. The spatiotemporal distribution characteristics of aftershocks and mainshocks also play an essential part in influencing the dynamic response of building structures [9]. Hence, it is worthwhile to effectively examine the dynamic response of a slope induced by the MAS. At present, researchers have paid attention to the damage of building structures subjected to the MAS, and more attention has been paid to the dynamic response caused by the combined action of the mainshock and the aftershock with the maximum magnitude. Pang et al. [10] discussed the susceptibility of a CFRD with a height of over 200 m produced by the MAS based on the multiple analysis. Zhou et al. [11] investigated the association between the intensity parameters of ground motion and the structural destruction under the MAS, and established the damage prediction model of the MAS based on the optimal parameters. Based on the Copula theory, Shen et al. [12] established a sequential random model of ground motion that can better represent the spatial correlation of sequential earthquakes. However, there are few studies on the response characteristics of slope caused by mainshock-aftershock sequence at present. In addition, the repeated method is mostly used to construct artificial mainshock-aftershock sequences in the above research, which can neither truly reflect the characteristics of the real mainshock-aftershock sequence nor properly consider the association between the aftershock and mainshock intensity. Therefore, the effect of an aftershock on the structural dynamic response cannot be reasonably responded to. The researchers cannot really grasp the safety of the slope when it is further subjected to an aftershock after the mainshock's initial damage due to the lack of research on the response of the slope induced by the MAS. Therefore, the advanced stochastic simulation method of the MAS requires of further investigation.
In order to consider the seismic slip danger of slopes, many methods have been developed, such as statistical analysis, the permanent displacement method, the pseudo-static method, the safety factor method, and the stress-deformation method [13][14][15][16]. Compared with the safety factor method, the permanent displacement caused by seismic action can assess the damage condition and seismic performance of the slope more reasonably [17][18][19]. In addition, traditional slope stability analysis methods are in general primarily based on deterministic analysis, which considers a series of elements affecting the slope stability as definitive factors. However, a many disaster results and geotechnical tests have revealed the apparent stochastic nature of variables affecting slope stability, such as external load, performance of materials, and model geometry [20][21][22]. Calculating and analyzing the stability of slopes by using deterministic methods can create many errors. The adoption of reliability theory has provided the opportunity of quantitative consideration of uncertainties in recent years [23,24]. By establishing extreme state equations, metrics such as probability of failure and reliability are employed to describe the safety of systems, thus providing a more complete guide for engineering design. Some traditional probabilistic methods, such as the First Order Second Moment method [25], the Monte Carlo method [26], the response surface method [27], and their improved forms, have been used to analyze the reliability of results and proven to be effective [28]. However, these methods have the characteristics of having difficulty obtaining the random dynamic information of the structure, a huge calculation scale, and a need to be coupled with the structural response analysis and continuous sample training and iteration. Therefore, it is demanding and challenging to apply the seismic random dynamic response and probability analysis of slopes with strong nonlinearity, complex problems, and a large calculation scale. The GPDEM is a newly developed method for probabilistic analysis [29,30]. At present, the GPDEM has been implemented for dynamic reliability assessment of bridges, slopes, earth-rock dams, and other engineering structures, and its efficiency has been confirmed by comparing it with the MCS method [30][31][32][33][34]. However, the feasibility of GPDEM in the evaluation of slope D P caused by a mainshock-aftershock sequence needs further verification.
Soil parameters have significant spatial variability due to differences in depositional conditions, loading history, and other geological processes [35][36][37]. Moreover, the soil parameters at different spatial locations have a certain relevance and are not completely independent, which makes the slope stability research more complex. In slope reliability analysis, two methods are generally adopted to imitate the variability of soil parameters. Assuming that the parameters are spatially homogeneous, the probability distribution model is applied to describe the inherent variability, which is called random variable model [38,39]. The random variable model assumes that the parameters in the study area are perfectly correlated and that differences in the physical and mechanical properties of the local and overall geotechnical properties at different points in space cannot be considered. This obviously does not conform to the actual situation of geotechnical engineering and cannot meet the needs of objective analysis and evaluation of the spatial variation characteristics of geotechnical parameters. The random field theory was first proposed and adopted by Cornel to describe the random characteristics of parameters in 1972 [40]. On this basis, the theory was gradually developed and refined by Vanmarcke [41]. The main idea is to treat the soil parameters at a certain location in the space as random variables subject to certain statistical laws and describe the spatial variability of soil parameters at different locations through variance reduction function, correlation distance, correlation function, etc. In contrast, the spatial variability of soil material parameters would be better depicted by the random field theory.
In this paper, a slope reliability analysis method based on the GPDEM and the Newmark displacement method is proposed to quantify the effect of spatial variability of soil parameters on dynamic reliability. The MAS and random field are generated by the random simulation method of MAS and the spectral representation method (SRM), and the D P of the slope is obtained by a nonintrusive analysis. Firstly, the impact of COV on a slope's dynamic stability is investigated from the mean value of D P , and then the influence of COV and PGA on a slope's dynamic reliability is explained from a probabilistic point of view by combining the GPDEM. The flowchart of the evaluation framework is depicted in Figure 1.

Generalized Probability Density Evolution Method
This paper considers the spatial variability of layered slope material parameters. Therefore, the slope dynamic equation under the action of MAS may be represented as:

Generalized Probability Density Evolution Method
This paper considers the spatial variability of layered slope material parameters. Therefore, the slope dynamic equation under the action of MAS may be represented as: where, K, C, and M represent the stiffness matrices, damping, and effective mass of the structure, respectively, and their basic parameters may be random. For convenience, the solution of Equation (1) can be formulated as: where H = (H 1 , H 2 , . . . , H n ) T ; n is the number of degrees of freedom. In addition, it is worth noting that it can be regarded as a variable Θ and t. Accordingly, variables such as acceleration, velocity, strain, and stress could also be expressed in the form similar to Equation (2). Therefore, in order to be more general, we can uniformly express the physical quantities of interest in the following form: where Hz = (H z,1 , H z,2 , . . . , H z , m ) T . Under the framework of probability conservation, the generalized probability density evolution equation of the stochastic process may be represented as: where p zΘ (z, θ, t) is the joint probability density function of the system (z, Θ), and m is the dimension of the equation, independent of the number of degrees of freedom of the system n. ∂p zΘ (z, θ, t) ∂t If only a single physical quantity is considered, Equation (4) is simplified to the form of Equation (5). Combining the initial conditions represented by Equation (6) and the boundary conditions represented by Equation (7), the structural reliability is finally obtained through mathematical processing. For simple problems, the analytical solution can be obtained in this way. For complex problems, such as large complex nonlinear systems, the numerical solution can only be obtained by mathematical methods, which could be achieved by the procedure as follows: (1) Point selection and probability assignment in probability space.
Discrete representative points are selected by some means in Ω Θ of random variables (such as the number-theory method, the point-selecting method by cutting the ball, the quasi-rotational symmetry point method, and the GF-deviation method).
The physical Equations (1) and (3) are solved, and the velocity of the required physical quantity is found for each given Θ = θ q .
(3) Solving probability density evolution equation. After the representative points are selected and the probabilities are assigned in step (1), the Equation (4) is transformed into: The corresponding initial conditions are transformed into: The result of the partial differential equation can be obtained by substituting • Z θ q , t m , as obtained in step (2), into Equations (8) and (9).
The result of p z (z, t) is acquired by the summation of all the above single results p zΘ z, θ q , t .

Spectral Representation Method
Many methods are currently employed to decompose random fields, for example, the midpoint method [42], spectral representation method [43], the spatial averaging method [44], the K-L decomposition method [45], and other methods. The spectral representation method has gradually become a widely used random field simulation method because of its good speed and accuracy in convergence to the objective function, and the generated sample function has ergodicity in all states.
By using the spectral representation method [43], the establishment of one-dimensional stationary random field could be expressed as: where Φ n is the independent phase angle uniformly distributed within the region of [0, 2π]; A n is the amplitude; and κ n is the frequency.
In Equation (12), S ff (κ n ) is the power spectrum function. The relationship between S ff (κ n ) and autocorrelation function can be acquired by Fourier transform, as shown below: Note that: The two-dimensional stationary random field can be characterized as: A n 1 n 2 cos κ 1n 1 x 1 + κ 2n 2 x 2 + Φ where Φ (1) n 1 n 2 and Φ (2) n 1 n 2 are individual random phase angles uniformly distributed within region of [0, 2π]; A n 1 n 2 and A n 1 n 2 are amplitude; κ 1n 1 and κ 1n 2 are frequency.
where κ 1u and κ 2u denote the number of truncated waves and meet the following relationship: In Equation (17), S f 0 f 0 (κ 1 , κ 2 ) is the power spectrum function. The relationship between S f 0 f 0 (κ 1 , κ 2 ) and autocorrelation function can be obtained through Fourier transform, as shown below:

Generation of Parametric Random Fields Based on Spectral Representation Method
In the random field simulation of slope strength parameters, since the value of strength parameters is usually positive, lognormal random field is employed to simulate the spatial difference and correlation of material parameters. The logarithmic stationary random field of slope strength parameters is established based on Equation (16).
where the V ij (θ) and W ij (θ) are mutually independent and obey the standard normal distribution; ω 1i and ω 2j are frequency coordinate values. x and z are the horizontal and vertical coordinate values of space. ξ ln and λ ln are the logarithmic standard deviation and logarithmic mean of parameters. σ ij is the standard deviation of i * M + j + 1.
where S ωω is the power spectral density function corresponding to the correlation function, which can be obtained by two-dimensional Fourier transform of the autocorrelation function.
where ρ(x, z) is the autocorrelation function. The Gaussian autocorrelation function with good stability and continuity is used for calculation. ∆ω 1 and ∆ω 2 are the discrete intervals of the frequency coordinate axes ω 1 and ω 2 , respectively.

Random Simulation of Mainshock-Aftershock Sequence (MAS)
Due to the limited number of measured records, it is necessary to generate the MAS ground motion through artificial simulation for seismic analysis of engineering structures. The existing method for constructing the MAS is to develop the magnitude relationship between the mainshock and aftershock and then separate and adjust the actual ground motion records (or artificial ground motion) to obtain the time histories of the MAS. However, in addition to the magnitude, the mainshock and aftershock are intimately associated in respect to source, propagation path, and local site impact, i.e., they are highly correlated in terms of spectrum characteristics, ground motion intensity, and duration. Obviously, a single magnitude parameter cannot accurately reflect the characteristics of the MAS. In addition, by adjusting the recorded ground motion or adopting the ground motion model of single shock, the changes of the amplitude, duration, and frequency spectrum of ground motion in the process of seismic wave propagation cannot be well reflected. Therefore, a random MAS simulation method accounting for the relevance between aftershock and mainshock based on Copula function is adopted to generate the MAS. This approach is characterized in greater depth in previous studies [46], and the primary steps of the approach can be simplified as follows: (1) Establishment of a physical random function model of the MAS.
(2) The real MASs are collected from the PEER to determine the physical parameters in the physical random function model of the mainshock-aftershock sequence. (3) Select a representative set of points of seismic parameters according to the GF difference. Then, establish the relevance between the aftershock and mainshock parameters based on the Copula theory. (4) Generate of a series of random MASs by using the narrowband harmonic superposition method.

Nonintrusive Analysis of Slope Dynamic Reliability
The biggest advantage of noninvasive randomness analysis is that the process of deterministic analysis and randomness analysis are independent of each other. The FE method is adopted to perform deterministic analysis without modifying the finite element kernel; therefore, the integration of deterministic analysis and stochastic analysis is realized, which significantly improves the reliability of the stochastic analysis results. By combining dynamic reliability analysis with finite element batch processing, this paper proposes a nonintrusive analysis frame of slope reliability considering spatial variability subjected to the MAS and compiles the interface program between dynamic reliability analysis and GeoStudio finite element software.
(1) Establish the slope of the FE model, divide the model mesh, set the boundary conditions, define the load loading method, define the material properties, and assign the elements in the SIGMA/W module with the parameter averages. Then, establish the corresponding relationship between the elements, groups, and material properties. Additionally, establish the stability analysis model in SLOPE/W, and save the FEM as a file with the extension name of ".xml". (2) The slope strength parameters are simulated by the spectral representation method.
N groups of data of parameters will be generated, and the parameters in the "xml" file will be replaced in batches with the newly generated n groups of data through MATLAB programming to obtain n new ".xml" files. (3) Use the UE text editor software to directly use GeoStudio to batch calculate the stability of the n new ".xml" files obtained in step (2). Then, output the calculation result files corresponding to each group of parameters. (4) The calculation results corresponding to all parameter groups are extracted in batch, and the D P is statistically analyzed.

Finite Element Model
In this study, the FEM adopted is a simplified layered soil slope based on geological data and field survey along the G317 Sichuan-Tibet Highway. The two-dimensional FEM was adopted to carry out the dynamic reliability analysis of layered soil slope. As demonstrated in Figure 2, the layered soil slope model is 40 m long and 24 m high. The FEM has been used in other studies to research the failure mode of a slope through numerical simulation and model tests [47]. In this study, two different types of layered soil slopes are used to research the dynamic reliability of the slope considering the spatial variability subjected to the MAS. The size of the grid is chosen to be 0.5 m, which ensures both computational efficiency and accuracy. Other information about the finite element model is introduced in detail in previous studies [6]. According to different soil layer distribution types, the two layered soil slopes are clayey soil-gravel soil-sandy soil-foundation soil and clayey soil-sandy soil-gravel soil-foundation soil. (The soil mass is arranged from top to bottom.) Water 2023, 15, x FOR PEER REVIEW 9 of 22

Finite Element Model
In this study, the FEM adopted is a simplified layered soil slope based on geological data and field survey along the G317 Sichuan-Tibet Highway. The two-dimensional FEM was adopted to carry out the dynamic reliability analysis of layered soil slope. As demonstrated in Figure 2, the layered soil slope model is 40 m long and 24 m high. The FEM has been used in other studies to research the failure mode of a slope through numerical simulation and model tests [47]. In this study, two different types of layered soil slopes are used to research the dynamic reliability of the slope considering the spatial variability subjected to the MAS. The size of the grid is chosen to be 0.5 m, which ensures both computational efficiency and accuracy. Other information about the finite element model is introduced in detail in previous studies [6]. According to different soil layer distribution types, the two layered soil slopes are clayey soil-gravel soil-sandy soil-foundation soil and clayey soil-sandy soil-gravel soil-foundation soil. (The soil mass is arranged from top to bottom.) During the initial static analysis and dynamic response analysis of the layered soil slope, the bottom of the FEM is constrained both horizontally and vertically. In addition, the right and left boundaries of the FEM are restrained horizontally during initial static analysis of the slope but not during the dynamic analysis.

Calculation Parameters
Various constitutive models are adopted to characterize the mechanical properties of different soil materials of the layered soil slope. The three layers of soil above the foundation soil are described by the equivalent linear model. The equivalent linear model is applied to the foundation soil because it is compacted. The correlation between damping ratio, shear modulus, and shear strain is presented in Figure 3. The material calculation parameters are the same as those employed by Huang [37], as shown in Table 1.  During the initial static analysis and dynamic response analysis of the layered soil slope, the bottom of the FEM is constrained both horizontally and vertically. In addition, the right and left boundaries of the FEM are restrained horizontally during initial static analysis of the slope but not during the dynamic analysis.

Calculation Parameters
Various constitutive models are adopted to characterize the mechanical properties of different soil materials of the layered soil slope. The three layers of soil above the foundation soil are described by the equivalent linear model. The equivalent linear model is applied to the foundation soil because it is compacted. The correlation between damping ratio, shear modulus, and shear strain is presented in Figure 3. The material calculation parameters are the same as those employed by Huang [37], as shown in Table 1.  The spatial variability of material is accounted for, and the parameter random field is generated by the above SRM. Then, the parameters are assigned to the well-constructed FEM. In this study, the parameters of each soil layer are presumed as independent of each other. The respectively, which are greater than the accuracy requirements (5.7~7.6) given by Ching and Phoon [48].

Input of the Mainshock-Aftershock Sequence
One of the randomly generated MAS is employed as the deterministic seismic wave input. The acceleration curve of MAS is presented in Figure 4.

Effect of Coefficient of Variation on Dynamic Reliability of Layered Soil Slope
The DP is adopted to assess the dynamic stability of slope, so it is necessary to define the critical DP of the soil slope. According to previous research [6], three DP thresholds (0.05 m, 0.5 m, and 1 m) were used to assess the dynamic reliability of the layered slope. The spatial variability of material is accounted for, and the parameter random field is generated by the above SRM. Then, the parameters are assigned to the well-constructed FEM. In this study, the parameters of each soil layer are presumed as independent of each other. The COV of c and ϕ (COV C and COV F ) are set as 0.1, 0.2, and 0.3. The vertical and horizontal autocorrelation distances (l h and l v ) are 20 m and 2 m, respectively. The horizontal dimension of the random field unit is 2 m, and the vertical dimension is 0.5 m. The ratios of the vertical and horizontal fluctuation ranges to the vertical and horizontal dimensions of the random field are δ h /l x = 20 √ π/2 = 17.7 and δ v /l y = 2 √ π/0.5 = 7.08, respectively, which are greater than the accuracy requirements (5.7~7.6) given by Ching and Phoon [48].

Input of the Mainshock-Aftershock Sequence
One of the randomly generated MAS is employed as the deterministic seismic wave input. The acceleration curve of MAS is presented in Figure 4.  The spatial variability of material is accounted for, and the parameter random field is generated by the above SRM. Then, the parameters are assigned to the well-constructed FEM. In this study, the parameters of each soil layer are presumed as independent of each other. The

Input of the Mainshock-Aftershock Sequence
One of the randomly generated MAS is employed as the deterministic seismic wave input. The acceleration curve of MAS is presented in Figure 4.

Effect of Coefficient of Variation on Dynamic Reliability of Layered Soil Slope
The DP is adopted to assess the dynamic stability of slope, so it is necessary to define the critical DP of the soil slope. According to previous research [6], three DP thresholds (0.05 m, 0.5 m, and 1 m) were used to assess the dynamic reliability of the layered slope.

Effect of Coefficient of Variation on Dynamic Reliability of Layered Soil Slope
The D P is adopted to assess the dynamic stability of slope, so it is necessary to define the critical D P of the soil slope. According to previous research [6], three D P thresholds (0.05 m, 0.5 m, and 1 m) were used to assess the dynamic reliability of the layered slope. A total of 86 sets of random material parameters were generated based on spectral representation to explore the impact of spatial variability on slope dynamic reliability. Figure 5 presents the distribution of D P discrete points for Case 1 when the COV C values are 0.1, 0.2, and 0.3. It is significant that when the COV C is small, the D P is small and relatively concentrated. With the increase in COV C , the range of variation of soil cohesion increases, the distribution of discrete points of D P becomes more discrete, and the mean of D P gradually increases. When the COV C is 0.1, the mean of D P caused by the MAS is 0.63 m, while the mean D P for the slope subjected to the single mainshock is 0.339 m. It is obvious that the mean D P of the slope induced by MAS is wider than the D P caused by the single mainshock. Moreover, the mean D P values of the slope under the MAS are 0.668 m and 0.725 m, respectively, when the COV C is 0.2 and 0.3. At this time, the mean D P of the slope due to the single mainshock is 0.368 m and 0.42 m. The mean value of D P increases continuously along with the increment of COV C , and the discrepancy of D P also shows a gradual tendency to increase. A total of 86 sets of random material parameters were generated based on spectral representation to explore the impact of spatial variability on slope dynamic reliability. Figure 5 presents the distribution of DP discrete points for Case 1 when the COVC values are 0.1, 0.2, and 0.3. It is significant that when the COVC is small, the DP is small and relatively concentrated. With the increase in COVC, the range of variation of soil cohesion increases, the distribution of discrete points of DP becomes more discrete, and the mean of DP gradually increases. When the COVC is 0.1, the mean of DP caused by the MAS is 0.63 m, while the mean DP for the slope subjected to the single mainshock is 0.339 m. It is obvious that the mean DP of the slope induced by MAS is wider than the DP caused by the single mainshock. Moreover, the mean DP values of the slope under the MAS are 0.668 m and 0.725 m, respectively, when the COVC is 0.2 and 0.3. At this time, the mean DP of the slope due to the single mainshock is 0.368 m and 0.42 m. The mean value of DP increases continuously along with the increment of COVC, and the discrepancy of DP also shows a gradual tendency to increase.     A total of 86 sets of random material parameters were generated based on spectral representation to explore the impact of spatial variability on slope dynamic reliability. Figure 5 presents the distribution of DP discrete points for Case 1 when the COVC values are 0.1, 0.2, and 0.3. It is significant that when the COVC is small, the DP is small and relatively concentrated. With the increase in COVC, the range of variation of soil cohesion increases, the distribution of discrete points of DP becomes more discrete, and the mean of DP gradually increases. When the COVC is 0.1, the mean of DP caused by the MAS is 0.63 m, while the mean DP for the slope subjected to the single mainshock is 0.339 m. It is obvious that the mean DP of the slope induced by MAS is wider than the DP caused by the single mainshock. Moreover, the mean DP values of the slope under the MAS are 0.668 m and 0.725 m, respectively, when the COVC is 0.2 and 0.3. At this time, the mean DP of the slope due to the single mainshock is 0.368 m and 0.42 m. The mean value of DP increases continuously along with the increment of COVC, and the discrepancy of DP also shows a gradual tendency to increase.   When the COV C is 0.1, the maximum value of PDF is 3.98, the fluctuation region of D P is 0.4-1.0 m, and the D P is primarily focused around 0.45 m. When the COV C is 0.3, the maximum value of PDF is 1.2, the fluctuation region of D P of slope is 0-1.5 m, and the D P is relatively centralized around 0.8 m. With the growth of the COV C , the PDF value gradually decreases, the curve gradually shifts to the right, the D P distribution range is wider, and the failure probability of the slope is higher. Table 2 shows the reliability information of the slope when the cumulative slips are 0.05 m, 0.5 m, and 1 m under different COV C and PGA. The dynamic reliability of the slope caused by the MAS decreases by 13% with the COV C increasing from 0.1 to 0.3 when the PGA is 0.5 g, and the displacement threshold is 1 m. When the COV C is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 13% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the dynamic reliability of the slope induced by MAS decreased by 35%.  Figures 7 and 8 show the probability information of DP of slope under different COVC. When the COVC is 0.1, the maximum value of PDF is 3.98, the fluctuation region of DP is 0.4-1.0 m, and the DP is primarily focused around 0.45 m. When the COVC is 0.3, the maximum value of PDF is 1.2, the fluctuation region of DP of slope is 0-1.5 m, and the DP is relatively centralized around 0.8 m. With the growth of the COVC, the PDF value gradually decreases, the curve gradually shifts to the right, the DP distribution range is wider, and the failure probability of the slope is higher. Table 2 shows the reliability information of the slope when the cumulative slips are 0.05 m, 0.5 m, and 1 m under different COVC and PGA. The dynamic reliability of the slope caused by the MAS decreases by 13% with the COVC increasing from 0.1 to 0.3 when the PGA is 0.5 g, and the displacement threshold is 1 m. When the COVC is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 13% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the dynamic reliability of the slope induced by MAS decreased by 35%.    Figures 7 and 8 show the probability information of DP of slope under different COVC. When the COVC is 0.1, the maximum value of PDF is 3.98, the fluctuation region of DP is 0.4-1.0 m, and the DP is primarily focused around 0.45 m. When the COVC is 0.3, the maximum value of PDF is 1.2, the fluctuation region of DP of slope is 0-1.5 m, and the DP is relatively centralized around 0.8 m. With the growth of the COVC, the PDF value gradually decreases, the curve gradually shifts to the right, the DP distribution range is wider, and the failure probability of the slope is higher. Table 2 shows the reliability information of the slope when the cumulative slips are 0.05 m, 0.5 m, and 1 m under different COVC and PGA. The dynamic reliability of the slope caused by the MAS decreases by 13% with the COVC increasing from 0.1 to 0.3 when the PGA is 0.5 g, and the displacement threshold is 1 m. When the COVC is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 13% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the dynamic reliability of the slope induced by MAS decreased by 35%.   Figure 9 provides the distribution of D P discrete points for Case 1 when the COV F values are 0.1, 0.2, and 0.3. As the COV F is small, the D P for the slope is low and concentrated. With the increase in COV F , the fluctuating region increases and the mean value of D P of slope gradually increases and becomes more discrete. When the COV F is 0.1, the mean D P of the slope under the MAS is 0.674 m and the mean D P of the slope under the single mainshock is 0.372 m. The mean D P of the slope induced by the MAS is larger than that under the single mainshock. When the COV F values are 0.2 and 0.3, the mean D P values of the slope caused by the MAS are 0.743 m and 0.795 m, respectively. At this time, the mean D P values of the slope due to the single mainshock are 0.409 m and 0.432 m. With increasing COV F , the mean value of D P of the slope continuously increases.   Figure 9 provides the distribution of DP discrete points for Case 1 when the COVF values are 0.1, 0.2, and 0.3. As the COVF is small, the DP for the slope is low and concentrated. With the increase in COVF, the fluctuating region increases and the mean value of DP of slope gradually increases and becomes more discrete. When the COVF is 0.     Figures 11 and 12 show the probability information of the DP of the slope under different COVF values. The maximum PDF value is 1.89 when the COVF is 0.1, and the DP of the slope is principally around 0.7 m. With the growth of the COVF, the PDF value gradu-  Figures 11 and 12 show the probability information of the D P of the slope under different COV F values. The maximum PDF value is 1.89 when the COV F is 0.1, and the D P of the slope is principally around 0.7 m. With the growth of the COV F , the PDF value gradually decreases, and the D P of the slope is more widely distributed. Table 3 shows the information of dynamic reliability for the slope when the cumulative slip is 0.05 m, 0.5 m, and 1 m under different COV F and PGA. When the PGA is 0.5 g and the threshold is 1 m, the dynamic reliability of the slope induced by the MAS decreases by 25% with the COV F increasing from 0.1 to 0.3. When the COV F is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 31% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the dynamic reliability of the slope under the MAS decreased by 17%.  Figures 11 and 12 show the probability information of the DP of the slope under different COVF values. The maximum PDF value is 1.89 when the COVF is 0.1, and the DP of the slope is principally around 0.7 m. With the growth of the COVF, the PDF value gradually decreases, and the DP of the slope is more widely distributed. Table 3 shows the information of dynamic reliability for the slope when the cumulative slip is 0.05 m, 0.5 m, and 1 m under different COVF and PGA. When the PGA is 0.5 g and the threshold is 1 m, the dynamic reliability of the slope induced by the MAS decreases by 25% with the COVF increasing from 0.1 to 0.3. When the COVF is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 31% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the dynamic reliability of the slope under the MAS decreased by 17%.    For Case 1, with the increase in COV C and COV F , the dynamic reliability gradually decreases, and the failure probability gradually increases under different displacement thresholds. In contrast, the dynamic reliability of slopes is more sensitive to COV F . Additionally, the dynamic reliability of slopes is more sensitive to the COV F . Similar conclusions were also obtained by Huang et al. [1], but the impact of aftershocks was not considered in their research.     Figures 15 and 16 describe the probability information of D P of slope under different COV C . The maximum PDF value is 3.2 when the COV C is 0.1, and the D P is mainly concentrated around 0.2-0.4 m. As the COV C improves, the PDF value gradually decreases, the curve gradually shifts to the right, and the D P of the slope is more widely distributed. Table 4 shows the information of dynamic reliability of D P for the slope when the cumulative slip is 0.05 m, 0.5 m, and 1 m under different COV C and PGA values. When the PGA is 0.5 g and the displacement threshold is 1 m, the reliability of the slope subjected to the MAS decreases by 7% with the COV C increasing from 0.1 to 0.3. When the COV C is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 7% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the reliability of D P of the slope produced by the MAS decreased by 8%. Figure 14 summarizes the distribution of DP discrete points of the slope under different PGA when the COVC is 0.3. When the PGA is 0.4 g and 0.6 g, the mean of DP values under the MAS are 0.261 m and 0.881 m, respectively. However, the mean of DP values under a single mainshock are 0.173 m and 0.577 m. The DP of the slope is raised step by step with the increase in PGA. The maximum PDF value is 3.2 when the COVC is 0.1, and the DP is mainly concentrated around 0.2-0.4 m. As the COVC improves, the PDF value gradually decreases, the curve gradually shifts to the right, and the DP of the slope is more widely distributed. Table 4 shows the information of dynamic reliability of DP for the slope when the cumulative slip is 0.05 m, 0.5 m, and 1 m under different COVC and PGA values. When the PGA is 0.5 g and the displacement threshold is 1 m, the reliability of the slope subjected to the MAS decreases by 7% with the COVC increasing from 0.1 to 0.3. When the COVC is 0.3, the dynamic reliability of the slope under the action of the MAS is also reduced by 7% compared with the single mainshock. In addition, with the PGA increasing from 0.4 g to 0.6 g, the reliability of DP of the slope produced by the MAS decreased by 8%.       Figure 17 presents the distribution of D P discrete points for Case 2 when the COV F values are 0.1, 0.2, and 0.3. At a small COV F value, the D P of the slope is small and concentrated. The variation range of soil increases with the increase in COV F , so that the mean of D P gradually increases and becomes more discrete. When the COV F is 0.1, the mean of D P subjected to the MAS is 0.299 m, and the mean of D P caused by the single mainshock is 0.247 m. The mean of D P under the MAS is larger than that caused by a single mainshock. When the COV F values are 0.2 and 0.3, the mean of D P produced by MAS is 0.404 m and 0.553 m, respectively. At this time, the mean of D P values due to a single mainshock are 0.293 m and 0.364 m. The mean D P of the slope keeps increasing with the growth of COV F , and the discrepancy of D P also displays a gradual increasing tendency.    Figures 19 and 20 show the probability information of D P of the slope under different COV F . The maximum PDF value is 3.2 when the COV F is 0.1, and the D P of the slope is mainly focused around 0-0.5 m. The PDF value gradually decreases, the curve gradually shifts to the right with the increase in the COV F , and the D P of the slope is more widely distributed. Table 5 lists the dynamic reliability of the slope when the cumulative slips are 0.05 m, 0.5 m, and 1 m under COV F and PGA. When the PGA is 0.5 g and the threshold is 1 m, the dynamic reliability of the slope induced by the MAS decreases by 17% with the COV F increasing from 0.1 to 0.3. When the COV F is 0.3, the dynamic reliability of the slope subjected to the MAS is also reduced by 14% compared with the single mainshock. In addition, the dynamic reliability of the slope under the MAS decreased by 14% with the PGA increasing from 0.4 g to 0.6 g. Figure 18 shows the distribution of DP discrete points of slope under the action of different PGA when the COVF is 0.3. When the PGA values are 0.4 g and 0.6 g, the mean of DP values of the slope produced by the MAS are 0.261 m and 0.881 m, respectively. However, the mean DP values of the slope due to the single mainshock are only 0.173 m and 0.577 m. The DP of the slope increases continuously with the increase in PGA. The maximum PDF value is 3.2 when the COVF is 0.1, and the DP of the slope is mainly focused around 0-0.5 m. The PDF value gradually decreases, the curve gradually shifts to the right with the increase in the COVF, and the DP of the slope is more widely distributed. Table 5 lists the dynamic reliability of the slope when the cumulative slips are 0.05 m, 0.5 m, and 1 m under COVF and PGA. When the PGA is 0.5 g and the threshold is 1 m, the dynamic reliability of the slope induced by the MAS decreases by 17% with the COVF increasing from 0.1 to 0.3. When the COVF is 0.3, the dynamic reliability of the slope subjected to the MAS is also reduced by 14% compared with the single mainshock. In addition, the dynamic reliability of the slope under the MAS decreased by 14% with the PGA increasing from 0.4 g to 0.6 g.      On the basis of the above reliability information, it can be seen that the conclusion of Case 1 is comparable to that of Case 2, i.e., the COV F has a significantly greater influence on the dynamic reliability of the slope. However, due to the different distribution of soil layers in the layered slopes, the influence of the COV on the failure probability of a slope is different. Compared with Case 2, the lower layer located on the slope is a sandy soil with poorer properties and its thickness is relatively deep. Therefore, the dynamic reliability of Case 1 is more significantly affected by the COV. The influence of different soil layers on the dynamic reliability and sliding surface position of slopes has been discussed in our previous research and detailed content can be found in [6].

Conclusions
A slope reliability analysis method based on GPDEM and the Newmark displacement method is proposed to quantify the impact of spatial variability of soil strength parameters on the dynamic reliability. The MAS and parameter random field are generated by the random simulation method of MAS and spectral representation method. Based on the Newmark method, the D P of layered soil slope is calculated by nonintrusive reliability analysis, and the influence of the COV C and COV F on the dynamic reliability of slope is compared. The conclusions of this study are as follows: (1) A reliability analysis method for D P of the slope is established based on the GPDEM and Newmark methods. Combined with the noninvasive stochastic analysis method, the failure probability of a slope can be quickly obtained. (2) According to the stochastic dynamic calculation results of the layered soil slope, COV C and COV F have a significant impact on the D P of the slope induced by the MAS. The mean of D P of the slope also presents a trend of increasing gradually with an increase in the COV C and COV F values. In contrast, the D P of slope is more sensitive to the COV F . (3) Affected by the randomness and nonlinearity of the materials, the PDF curve has nonuniform single or double peaks. As the COV increases, the PDF curve becomes lower and wider, and the failure probability of the layered soil slope increases. When the D P threshold is 1 m and PGA is 0.5 g, the dynamic reliability of the soil slope is continuously reduced, and the failure probability is even increased by about 20% with the COV increasing from 0.1 to 0.3. (4) The impact of aftershocks on the D P of the soil slope cannot be ignored. The mean of D P of the slope induced by the MSA is larger than that under a single mainshock. The dynamic reliability of the slope caused by the MAS can even be reduced by 7-30% compared with a single mainshock when the displacement threshold is 1 m and the COV C is 0.3. Additionally, the impact of aftershocks on the D P of slope increases with an increase in PGA.  Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.