Stochastic Dynamic Response Analysis of the 3D Slopes of Rockﬁll Dams Based on the Coupling Randomness of Strength Parameters and Seismic Ground Motion

: Because rockﬁll strength and seismic ground motion are dominant factors affecting the slope stability of rockﬁll dams, it is very important to accurately characterize the distribution of rockﬁll strength parameters, develop a stochastic ground motion model suitable for rockﬁll dam engineering, and effectively couple strength parameters and seismic ground motion to precisely evaluate the dynamic reliability of the three-dimensional (3D) slope stability of rockﬁll dams. In this study, a joint probability distribution model for rockﬁll strength based on the copula function and a stochastic ground motion model based on the improved Clough-Penzien spectral model were built; the strength parameters and the seismic ground motion were coupled using the GF-discrepancy method, a method for the analysis of dynamic reliability of the 3D slope stability of rockﬁll dams was proposed based on the generalized probability density evolution method (GPDEM), and the effectiveness of the proposed method was veriﬁed. Moreover, the effect of different joint distribution models on the dynamic reliability of the slope stability of rockﬁll dams was revealed, the effect of the copula function type on the dynamic reliability of the slope stability was analysed, and the differences in the dynamic reliability of the slope stability under parameter randomness, seismic ground motion randomness, and coupling randomness of parameters and seismic ground motion were systematically determined. The results were as follows: the traditional joint distribution models ignored related nonnormal distribution characteristics of rockﬁll strength parameters, which led to excessively low calculated failure probabilities and overestimations of the reliability of the slope stability; in practice, we found that the optimal copula function should be selected to build the joint probability distribution model, and seismic ground motion randomness must be addressed in addition to parameter randomness.


Introduction
The uncertainty of structural parameters and seismic ground motion is a dominant factor affecting the engineering safety of rockfill dams.Therefore, rational representations of the uncertainties of structural parameters and seismic ground motion have become key points of aseismic designs for rockfill dams.With the continuous development and improvement of reliability analysis theories and methods, reliability analysis has become a leading method for analysing and evaluating uncertainty factors in civil engineering.However, though traditional deterministic analysis methods still predominate in rockfill dam engineering, the use of reliability analysis is just beginning.
The current principal issues in the reliability-theory-based seismic safety analysis of the slope stability of rockfill dams are as follows.(1) The rational characterization of the uncertainty factors of rockfill strength parameters and seismic ground motion is required.
The amount of test data for rockfill dams is limited; with limited data, it is only possible to acquire a marginal distribution function and correlation coefficients.Most traditional methods employ rockfill strength parameters as a standalone normally or lognormally distributed variables without considering their correlations in the study of the reliability of the slope stability of rockfill dams.Additionally, the frequency and amplitude of seismic ground motion, which are highly random, exhibit nonstationarity over time; such nonstationarity has an important impact on the dynamic response and reliability of engineering structures.Accordingly, it is extremely important for studies on the dynamic reliability of the slope stability of rockfill dams to rationally characterize rockfill strength parameter-related nonnormal distribution characteristics from limited data and to build stochastic ground motion models suitable for the nonstationary random process.(2) The rockfill strength parameters are coupled with the uncertainty of the seismic ground motion.Traditional methods normally deal with the dynamic reliability of the slope stability of rockfill dams by separating the rockfill strength parameters from the randomness of seismic ground motion.For instance, Lizarraga and Lai [1], Sukkarak et al. [2,3], and Wu et al. [4] studied the reliability of the slope stability of rockfill dams based on the uncertainty of the rockfill strength parameters, and Xu et al. [5], Pang et al. [6][7][8][9], Kemal [10], Huang and Xiong [11], Ji et al. [12,13], Johari A et al. [14,15], and Wang et al. [16] studied the reliability of the slope stability of rockfill dams based on seismic ground motion randomness.When using the random combination [17,18] of sample strength parameters and seismic ground motions, problems such as an insufficient sample size and a heavy computing load may arise, so the combination is not coupled.Therefore, it is critical to carry out in-depth studies on the seismic safety analysis of the slope stability of rockfill dams based on the coupling randomness of rockfill strength parameters and seismic ground motion.(3) The dynamic reliability analysis method is based on the coupling randomness of the strength parameters and seismic ground motion.Current dynamic reliability analysis methods generally simplify high-dimensional nonlinear stochastic dynamic analysis by separately analysing the random structure problem and the random vibration problem.For example, the stochastic simulation method [19][20][21][22], the stochastic perturbation method [23][24][25], and the orthogonal polynomial expansion method [26,27] were developed based on the stochastic structure problem.The statistical linearization method [28,29], the Fokker-Planck-Kolmogorov (FPK) method [30][31][32], the stochastic averaging method [33][34][35], the path integral method [36][37][38], and the non-Gaussian closure method [39] were developed based on the random vibration problem.However, for engineering structures with a large number of degrees of freedom (DOFs) and nonlinear and random coupling effects, all the aforementioned methods face tremendous challenges.The use of Monte Carlo simulation (MCS) is an excellent solution to such problems [40,41], but it involves an immense computing load that restricts its practical engineering applications.In recent years, Li and Chen [42][43][44][45] developed a probability density evolution concept for stochastic dynamic system analysis called the generalized probability density evolution method (GPDEM).The basic idea of the GPDEM is to obtain the decoupled generalized probability density evolution equation through random events of the probability conservation principle, to analyse the probability of the reliability of complex structures, and to obtain the probability density function and reliability of the structural dynamic response at any time.This method not only analyses all the random factors in the system but also solves the problem of a high degree of freedom.The calculated results can contain all the probability information of a structural response.The GPDEM has become an effective method for analysing the random vibration and dynamic reliability of structures.For engineering design, the method has the following problems.(1) Compared with the traditional deterministic calculation method, the GPDEM requires a large number of deterministic calculations.Therefore, finding ways to reduce the number of samples of random variables and ensure that the extracted sample points can reasonably characterize the distribution characteristics of random variables is particularly critical to improving the calculation efficiency.(2) Most engineering designs use a single safety index to evaluate the safety performance of structures.The evaluation method is simple, but this index cannot truly reflect the impact of material parameters and the uncertainty of the actual loads on the safety performance of the structures.On the basis of deterministic calculations, the GPDEM can quantitatively consider a variety of uncertain factors and comprehensively evaluate the safety performance of structures based on probability analysis.However, this method increases the use of the reliability analysis process, and the analysis theory and evaluation method are relatively complex.Hence, developing a stochastic dynamic analysis method suitable for high dimensions and nonlinearity is key to the analysis of the dynamic reliability of the slope stability of rockfill dams.
In conclusion, in this study, a joint probability distribution model of rockfill strength based on the copula function was built by fully considering the coupling randomness of rockfill strength parameters and seismic ground motion, a stochastic ground motion model was built based on the improved Clough-Penzien spectral model, the strength parameters were coupled with seismic ground motion by the GF-discrepancy method, and an analysis method for the dynamic reliability of the three-dimensional (3D) slope stability of rockfill dams by using the GPDEM was proposed.With a concrete-faced rockfill dam as an example, the calculation results from MCS were compared to demonstrate the effectiveness of this method.On this basis, the effect of different joint distribution models on the dynamic reliability of the slope stability was revealed, the effect of the copula function type on the dynamic reliability of the slope stability was analysed, and the differences in the dynamic reliability of the slope stability under parameter randomness, seismic ground motion randomness, and the coupling randomness of parameters and seismic ground motion were systematically determined.

3D Slope Stability Analysis Method
Currently, the reliability analysis of rockfill dam slope stability is mostly performed using a two-dimensional (2D) model based on a clear concept and mature theory; this model has convenient calculations and wide engineering applications [46,47].It is simple to implement this 2D model in rockfill dam engineering projects in wide river valleys and low dams, and the results have been conservative.Advances in construction technology have enabled the heights of rockfill dams to increase from 200 to 300 m, making 3D effects more significant.Thus, it is no longer reasonable to use 2D model analysis, and the effective reinforcement range is difficult to determine.To avoid the shortcomings of 2D model analysis and obtain more practical results, it is essential to carry out a threedimensional (3D) slope stability analysis.At present, 3D slope dynamic stability analysis mainly includes the limit equilibrium, limit analysis, and numerical analysis methods.The limit equilibrium method [48][49][50] has a simple principle and a clear concept, and it is widely used in 2D slope stability analysis.It plays an important role in promoting slope stability analysis and research.However, some assumptions should be made in its analysis process, and there will be differences between predicted and practical results.If the limit equilibrium method is adopted in 3D slope stability analysis, more assumptions should be introduced to make the equilibrium problem statically determinable and solvable.The existing methods mainly focus on precision and complexity.Additionally, limited by its basic theoretical framework, the limit equilibrium method cannot consider the nonlinear constitutive nature of soil, the action of water, and the propagation of ground motion.The limit analysis method [51,52] is based on plastic mechanics theory, and the energy balance equation is solved via the numerical analysis method, which avoids the complex stress-strain incremental iteration process and directly analyses the final limit state.However, this method is similar to the limit equilibrium method, and certain assumptions are still needed in the calculation process.Restricted by the sliding surface search method, it is difficult to determine the most dangerous failure mechanism of a 3D slope.Previous research has primarily focused on the theoretical method, which has limited software development and engineering applications.With rapid improvements in computer capabilities, numerical calculation methods have quickly developed and now mainly include the finite element method [53,54], the finite element strength reduction method [55][56][57], and the finite difference method [58,59].Numerical calculation methods can solve slope problems under complex conditions without assumptions of the conditions and obtain stress, strain, displacement, and other information for any position.However, the determination criteria of the critical state of the slope and the position of the critical sliding surface remain unsolved.
Therefore, this research group analysed the advantages and disadvantages of various slope stability analysis methods using the large finite element software GEODYNA in conjunction with the limit equilibrium and finite element methods.A 3D slope stability analysis method and a calculation program (FEMSTABLE-3D) were developed, and the stability analysis method was improved by considering dangerous surface pre-storage.A highly efficient and feasible 3D slope stability analysis method for slip surface stresses was thus created.The basic theory and verification process of the program are detailed elsewhere [60].For the finite element calculation and analysis, we adopted the analysis software GEODYNA [61] developed by the Institute of Engineering Earthquake Resistance of the Dalian University of Technology.A flow chart of the slope stability analysis program is shown in Figure 1.method, which avoids the complex stress-strain incremental iteration process and directly analyses the final limit state.However, this method is similar to the limit equilibrium method, and certain assumptions are still needed in the calculation process.Restricted by the sliding surface search method, it is difficult to determine the most dangerous failure mechanism of a 3D slope.Previous research has primarily focused on the theoretical method, which has limited software development and engineering applications.With rapid improvements in computer capabilities, numerical calculation methods have quickly developed and now mainly include the finite element method [53,54], the finite element strength reduction method [55][56][57], and the finite difference method [58,59].Numerical calculation methods can solve slope problems under complex conditions without assumptions of the conditions and obtain stress, strain, displacement, and other information for any position.However, the determination criteria of the critical state of the slope and the position of the critical sliding surface remain unsolved.
Therefore, this research group analysed the advantages and disadvantages of various slope stability analysis methods using the large finite element software GEODYNA in conjunction with the limit equilibrium and finite element methods.A 3D slope stability analysis method and a calculation program (FEMSTABLE-3D) were developed, and the stability analysis method was improved by considering dangerous surface pre-storage.A highly efficient and feasible 3D slope stability analysis method for slip surface stresses was thus created.The basic theory and verification process of the program are detailed elsewhere [60].For the finite element calculation and analysis, we adopted the analysis software GEODYNA [61] developed by the Institute of Engineering Earthquake Resistance of the Dalian University of Technology.A flow chart of the slope stability analysis program is shown in Figure 1.

Building the Joint Distribution Model of Rockfill Strength Parameters Based on the Copula Function
Many types of copula functions exist, but most can only simulate a limited range of correlation coefficients.Table 1 shows six commonly used copula functions that can effectively characterize the significant positive correlations between the rockfill strength parameters.
According to Sklar's theorem for a binary distribution [62], a joint distribution function is composed of the distribution function of chosen variables and the copula function between those variables.
The joint distribution function of the strength parameters ϕ 0 and ∆ϕ is F(φ 0 , ∆φ) is: where u 1 = F 1 (ϕ 0 ) and u 2 = F 2 (∆ϕ) are the distribution functions of the intensity parame- ters and θ is a parameter of the copula function.
The parameter θ must be determined to define the copula function.This parameter can be obtained using the Pearson correlation coefficient ρ and the Kendall rank correlation coefficient [63].The definition of the correlation coefficient can be used to relate θ and ρ (Nelsen, 2006) as follows: Thus, θ can be obtained.However, except for the Gaussian copula function, most copula functions are difficult to solve for use in integration.In the literature [64,65], θ of the Gaussian copula function has been obtained from the Pearson correlation coefficient as follows: where are the variables of the standard normal distribution.The θ obtained for the Gaussian copula function can be used to determine the Kendall rank correlation coefficient τ as follows: Finally, θ for different copula functions can be determined using the formula given below.

Nonstationary Stochastic Ground Motion Model
The scientific community has proposed and developed a series of stable random ground motion models, such as K-T [66] and C-P [67], to provide reasonable descriptions of the random process of ground motion.Most of these models produce relatively reasonable results for linear structures.
However, stationary and nonstationary ground motion models produce quite different results for nonlinear structures.The steady ground motion model underestimates the cumulative damage and strength degradation or structural rigidity of a structure.As the intensity and frequency of ground motion are characterized by strong randomness, ground motion should be a strictly nonstationary excitation process.Therefore, it is very important to develop a reasonable nonstationary ground motion model for use in the dynamic reliability analysis of nonlinear structures.Some advances have been made in recent studies on frequency nonstationary power spectrum models, such as a time-frequency modulation function for nonstationary ground motion [68] and the Clough-Penzien power spectrum model [69].In 2011, Cacciola and Deodatis proposed an improved Clough-Penzien bilateral evolution power spectrum model [70], which laid a foundation for the nonstationary random excitation analysis of nonlinear structures.
(1) Improved Clough-Penzien power spectrum model The power spectrum density function S is defined within the improved Clough-Penzien bilateral evolution power spectrum model as follows: •S 0 (t) (7) where A(t) is an intensity modulation function, for which the recommended form is [71]: where d is a parameter that controls the shape of the function and c is the moment at which the peak acceleration occurs.Frequency nonstationarity is incorporated into the power evolution spectrum using the following parameters: where ω 0 is the initial circular frequency and ξ 0 is the initial damping ratio, both of which are determined by the site characteristics; a and b are parameters determined from the site characteristics and ground categories; and T is the duration of the ground motion acceleration time course and depends on the site category.The spectral parameter S 0 (t) is expressed as: where γ is the equivalent peak factor and a max is the average peak acceleration.
(2) Using the spectral expression-random function method to generate random ground motions The generalized Clough-Penzien power spectrum model was used in this study.Considering the hydraulic seismic code, a random function was selected to generate nonstationary random ground motions, and GPDEM theory was used to analyse the dynamical response and reliability of a rock fill dam.
The spectral expression of the random process of nonstationary ground motion acceleration is simulated using the following formula [8,72]: ..
where ω k = k∆ω and S .. X g denotes the power spectral density function for bilateral evolution, which satisfies the following condition at a frequency ω = 0.
A random function expression is constructed using the orthonormal random variables {X k , Y k }.Let X n and Y n (n = 1, 2, . . ., N) be any two sets of standard orthogonal random variables.X n and Y n are functions of two independent basic random variables Θ 1 and Θ 2 .The random function is: where cas(x) = cos(x) + sin(x) is the Hartley function and the basic random variables Θ 1 and Θ 2 are independent of each other and follow a uniform distribution over [0, 2π].After certain deterministic mapping, X n , Y n can become the orthonormal random variable X k , Y k required by Equation ( 12) and be uniquely determined.

Generalized Probability Density Evolution Method (GPDEM)
The GPDEM was proposed by Li and Chen [42][43][44][45].The principle of conservation of probability is applied to a dynamic equation of state to produce a decoupled generalized probability density evolution equation.A stochastic dynamic reliability analysis is carried out by combining the spatial point selection technique, the deterministic analysis method, and the finite difference technique.The GPDEM can be used to consider a variety of random factors and evaluate the dynamic reliability of engineering structures with several or even dozens of variables.Thus, the GPDEM constitutes a novel approach that can be used to investigate the stochastic dynamic response, reliability, and optimal control of nonlinear engineering structure systems.
In general, the differential equation of motion of a dynamic system with n degrees of freedom can be expressed as:
where M, C, and K denote the effective mass, damping, and stiffness matrix, respectively, of a structure, for which the basic parameters may be random; ..

X(t),
. X(t), and X(t) denote the acceleration, speed, and displacement vector of the structural response, respectively; .. X(Θ, t) denotes the random dynamic excitation process; and Θ is the random vector for the entire system.
For general engineering dynamic systems, a unique and continuous solution exists under the given initial conditions and can be written as: where The speed can be expressed as: .
where h = (h 1 , h 2 , . . ., h n ) T .In engineering practice, there may be other physical quantities of interest, such as a safety factor, the displacement, the internal force of the section, and the stress or strain at key points.The physical quantity analysed in this study is denoted as Z = (Z 1 , Z 2 , . . . ,Z n ) T , where n is the number of physical quantities of interest, and the corresponding displacement and velocity are expressed as follows: Any physical parameter Θ can be used as a random variable in the GPDEM evolution equation.In this study, the safety factor of the dynamic stability of a rockfill dam slope was used as a random variable, and its dynamic reliability was analysed.Based on the principle of conservation of probability [42][43][44][45], the generalized probability density evolution equation is expressed as: The initial conditions are: The joint PDF is:

Method for Selecting Representative Discrete Points in Probability Space
The selection of discrete representative points plays a key role in solving for the probability density evolution and has an important influence on the accuracy and efficiency of the results.In the GF-discrepancy method [73,74], a probability distribution replaces equal weights to incorporate the influence of the global characteristics of the point set on the molecular domain of the spatial section, thereby reflecting the connection between the points, spatial structural information, and global characteristics of the probability distribution.This method uses few selected points, has strong representativeness, and exhibits high precision.Thus, the complexity incurred from an exponential increase in the computational workload with the dimension is prevented, which increases the computational efficiency.
The GF-discrepancy method consists of two main steps: First, an initial point set is generated from a Sobol sequence [75] and is rearranged to minimize the GF-deviation.The initial point set x q = (x q,1 , x q,2 , . . ., x q,i ) consists of the Sobol point set u q = (u q,1 , u q,2 , . . ., u q,i ) (q = 1, 2, . . ., n); where i represents the i-th random variable) and is given as follows: x where is the inverse cumulative distribution function of the i-th random variable.Second, a probability is assigned to bring a set of n point sets closer together.

Analysis Procedure for the Dynamic Reliability of the 3D Slope Stability of Rockfill Dams Based on Coupling Randomness
A joint distribution model of rockfill strength parameters based on the copula function was constructed in MATLAB, a stochastic ground motion model was built based on the improved Clough-Penzien spectral model, and the GF-deviation method with the GPDEM was used to achieve an efficient and feasible analysis method for the dynamic reliability of the 3D slope stability of rockfill dams.The principal computational procedure of this method consisted of the following four steps: (1) Representative points were chosen by the GF-discrepancy method.A total of 144 sets of four-dimensional (4D) coupled, standard, and uniformly distributed random variables (i.e., 144 sets of 4D rockfill strength parameters-seismic ground motion parameter coupled sample data, with uniform distribution) were simultaneously generated within unit space [0,1] 4 .
(2) The statistical analysis of the data was performed for the rockfill strength parameters, the second-order statistics and correlation coefficients of the strength parameters were calculated, the optimal marginal distribution function and the optimal copula function were identified, and the joint distribution model of the rockfill strength was built based on the optimal copula function.
(3) The irrelevant standard uniformly distributed random variables were converted into related nonnormal variables to suit the strength distribution characteristics based on the joint distribution model of the rockfill strength, 144 sets of random nonstationary time-history samples of seismic ground motion acceleration were generated based on the improved Clough-Penzien spectral model, and finally, 144 sets of random sample data were generated through the coupling of rockfill strength parameters and seismic ground motion.
(4) A deterministic dynamic time-history analysis was performed based on the resulting 144 sets of coupled random sample data, and the dynamic reliability of the 3D slope stability was analysed with the GPDEM.Figure 2 shows the details of the calculation process.
improved Clough-Penzien spectral model, and finally, 144 sets of random sample data were generated through the coupling of rockfill strength parameters and seismic ground motion.
(4) A deterministic dynamic time-history analysis was performed based on the resulting 144 sets of coupled random sample data, and the dynamic reliability of the 3D slope stability was analysed with the GPDEM.Figure 2 shows the details of the calculation process.

Computation Models and Parameters
Figure 3 shows a typical CFRD model with a dam height of 200 m.The dam crest was set to 16 m wide, and the upstream and downstream slopes were set to 1.4 and 1.6, respectively.The panel thickness was set to 0.3 + 0.0035H, where H is the dam height.A contact surface unit was arranged between the face plate and the dam rockfill.The 3D finite element model of the dam consisted of the dam body, bedrock, and mountains on both sides of the dam, with 80,573 nodes and 75,980 elements.The grid size was set to 4 × 4 × 8 m.The element grid size was applied to the dam body, the mountains on both sides, and the bedrock.The foundation was set to extend a certain distance in depth and the horizontal direction.Viscoelastic artificial boundaries were added at the bottom and around the bedrock.The dam was filled in layers.The water storage height was set to 175 m.There was no water downstream.The stability of the dam slope under dynamic loads

Computation Models and Parameters
Figure 3 shows a typical CFRD model with a dam height of 200 m.The dam crest was set to 16 m wide, and the upstream and downstream slopes were set to 1.4 and 1.6, respectively.The panel thickness was set to 0.3 + 0.0035H, where H is the dam height.A contact surface unit was arranged between the face plate and the dam rockfill.The 3D finite element model of the dam consisted of the dam body, bedrock, and mountains on both sides of the dam, with 80,573 nodes and 75,980 elements.The grid size was set to 4 × 4 × 8 m.The element grid size was applied to the dam body, the mountains on both sides, and the bedrock.The foundation was set to extend a certain distance in depth and the horizontal direction.Viscoelastic artificial boundaries were added at the bottom and around the bedrock.The dam was filled in layers.The water storage height was set to 175 m.There was no water downstream.The stability of the dam slope under dynamic loads was analysed by imposing constraints on the bottom of the bedrock and surrounding mountains.
To accurately characterize the probability distribution types and correlation characteristics of the rockfill strength parameters and establish a reasonable joint probability distribution model, detailed information was compiled from 124 rockfill dam projects to yield 767 sets of rockfill strength parameter data (shown in Figure 4).Outliers in the statistical data were eliminated using 3σ rules to yield 755 sets of valid data.
A static calculation was performed using the E-B model with the parameters shown in Table 2.A dynamic calculation was performed using an equivalent linear viscoelastic model with the parameters shown in Table 3.The panel, surrounding bedrock, and mountain were modelled using a linear elastic model with the panel parameters shown in Table 4.
was analysed by imposing constraints on the bottom of the bedrock and surrounding mountains.To accurately characterize the probability distribution types and correlation characteristics of the rockfill strength parameters and establish a reasonable joint probability distribution model, detailed information was compiled from 124 rockfill dam projects to yield 767 sets of rockfill strength parameter data (shown in Figure 4).Outliers in the statistical data were eliminated using 3σ rules to yield 755 sets of valid data.A static calculation was performed using the E-B model with the parameters shown in Table 2.A dynamic calculation was performed using an equivalent linear viscoelastic model with the parameters shown in Table 3.The panel, surrounding bedrock, and mountain were modelled using a linear elastic model with the panel parameters shown in Table 4.To accurately characterize the probability distribution types and correlation characteristics of the rockfill strength parameters and establish a reasonable joint probability distribution model, detailed information was compiled from 124 rockfill dam projects to yield 767 sets of rockfill strength parameter data (shown in Figure 4).Outliers in the statistical data were eliminated using 3σ rules to yield 755 sets of valid data.A static calculation was performed using the E-B model with the parameters shown in Table 2.A dynamic calculation was performed using an equivalent linear viscoelastic model with the parameters shown in Table 3.The panel, surrounding bedrock, and mountain were modelled using a linear elastic model with the panel parameters shown in Table 4.  Five marginal distribution functions were selected to reasonably and comprehensively determine the optimal probability distribution type of the rockfill strength parameters and are shown in Table 5 (where µ denotes the mean and σ denotes the standard deviation).The Akaike information criterion (AIC) [76] was used to determine the optimal marginal distribution function.The AIC is a standard to measure the goodness of statistical model fitting.Its basic idea is to evaluate the goodness of model fitting data through the concept of entropy.It has become the main method to evaluate the goodness of data fitting.The calculation results are shown in Table 6.
Table 5.Five alternative marginal distribution functions.

Distribution Type Probability Distribution Function Probability Density Function Remarks
Truncated Normal Φ x−p q The identification of the optimal marginal distribution function is: In the formula, x i (i = 1, 2, . . ., N) is the test data of the nonlinear strength parameters, N is the number of samples, f (x i ; p, q) is the probability density function of the alterna- tive distribution type, p and q are the distribution parameters, and k 1 is the number of distribution parameters of the alternative edge distribution function.
The test data of the geophysical parameters, a variety of edge distribution functions, are known.The AIC values of various edge distribution functions were calculated by Formula (25).The edge distribution with the minimum AIC values was found to be the best edge distribution function for fitting the rock and soil mass parameter data.
(2) Optimal Copula Function The AIC was used to determine the optimal copula function.The calculation results are shown in Table 7.The calculation formula used to identify the optimal copula function is: where (u i , v i ), i = 1, 2, . . ., N represents the test data of the strength parameters of the dam material, N represents the total number of data, D(u i , v i ; θ) is the probability density function of the alternative copula function, and k 1 is the number of parameters of the alternative copula function.
The test data of the geophysical parameters, a variety of copula probability density functions, are known.The AIC values of various copula functions were calculated using Formula (26).The copula functions with the minimum AIC values were found to be the best copula functions for fitting the rock and soil mass parameter data.

Seismic Response of 3D Concrete-Faced Rockfill Dam Slope Stability
The dam slope dynamic stability reliability under coupled random conditions could be calculated following the steps presented in Section 2.6.To identify the differences in the computational results from traditional joint distribution models, 144 sets of coupled random sample data were obtained based on traditional joint distribution models using the aforementioned method.The accuracy of this method was verified by the MCS method (2000 sets of data were randomly selected).The MCS sampling method was conducted similarly to the flow shown in Figure 1, as 2000 groups of 4D independent standard uniformly distributed random variables were randomly selected.Based on the established joint distribution model of rockfill parameters and the standard normal joint distribution model, 2000 groups of 2D data were converted into relevant nonnormal distribution variables consistent with the strength parameters of the rockfill.Another 2000 groups of 2D data generated nonstationary motion time-history samples and then formed 2000 groups of random samples while considering the coupling of the rockfill strength and ground motion.
The safety factor time-history of the stability of the dam slope was determined through a time-history analysis of over 4288 sets of deterministic dynamic finite elements that were substituted into the GPDEM equation to identify the seismic response of the 3D slope stability of concrete-faced rockfill dams.
Figures 5 and 6 show a comparison of the safety factor mean and standard deviation (SD) calculated by the GPDEM and MCS and a comparison of the computational results from different joint distribution models, respectively.As shown in the figures, the results of the GPDEM and MCS methods were in good agreement, which indicates that this method is accurate and reliable, with a greatly improved calculation efficiency.The computational results from both joint distribution models exhibited a similar distribution trend, but there were still significant differences in the local areas.Compared to traditional joint distribution models, the joint distribution model built based on the Plackett copula function yielded smaller results.The cause is stated as follows: the joint distribution model based on the optimal copula function could reasonably characterize the distribution characteristics of parameters; hence, the samples were more concentrated, and the discrete range was narrower compared to traditional joint distribution models.
computational results from both joint distribution models exhibited a similar distribution trend, but there were still significant differences in the local areas.Compared to traditional joint distribution models, the joint distribution model built based on the Plackett copula function yielded smaller results.The cause is stated as follows: the joint distribution model based on the optimal copula function could reasonably characterize the distribution characteristics of parameters; hence, the samples were more concentrated, and the discrete range was narrower compared to traditional joint distribution models.Figures 7 and 8 show the evolution and contour map, respectively, of the safety factor PDF calculated with both joint distribution models over time within the time interval of [9,10] s.As shown in the figures, under the coupling randomness of the strength parameters and seismic ground motion, the evolution of the seismic response of the concrete-faced rockfill dam was found to be complex and irregular, which indicates that the dam slope stability safety factor significantly changed over time.Compared to the joint distribution model based on the optimal copula function, the sample data points extracted by the traditional joint distribution model had a wider distribution, larger values, a larger probability density, and a wider threshold range.Figures 9 and 10 show the PDF and cumulative distribution function (CDF) curves, respectively, calculated by the two joint distribution models at typical times (4.1, 8.2, and 10 s).As shown in the figures, the shape of the PDF curve was completely different from that of the CDF curve at different times, and the curves fluctuated and significantly changed.Compared to the joint distribution model based on the optimal copula function, the results calculated by the traditional distribution model were larger and the distribution was wider.The GPDEM and MCS calculation results exhibited considerable consistency, which further demonstrates that the method is accurate.Figures 9 and 10 show the PDF and cumulative distribution function (CDF) curves, respectively, calculated by the two joint distribution models at typical times (4.1, 8.2, and 10 s).As shown in the figures, the shape of the PDF curve was completely different from that of the CDF curve at different times, and the curves fluctuated and significantly changed.Compared to the joint distribution model based on the optimal copula function, the results calculated by the traditional distribution model were larger and the distribution was wider.The GPDEM and MCS calculation results exhibited considerable consistency, which further demonstrates that the method is accurate.Compared to the joint distribution model based on the optimal copula function, the results calculated by the traditional distribution model were larger and the distribution was wider.The GPDEM and MCS calculation results exhibited considerable consistency, which further demonstrates that the method is accurate.10 s).As shown in the figures, the shape of the PDF curve was completely different from that of the CDF curve at different times, and the curves fluctuated and significantly changed.Compared to the joint distribution model based on the optimal copula function, the results calculated by the traditional distribution model were larger and the distribution was wider.The GPDEM and MCS calculation results exhibited considerable consistency, which further demonstrates that the method is accurate.

Analysis of the Dynamic Reliability of the 3D Slope Stability of Rockfill Dams
In this study, equivalent extremum events were constructed in the safety factor time-history, a "virtual" time parameter stochastic process was created, and the dynamic stability reliability (i.e., the PDF and CDF curves of the equivalent extremum safety factor (Figure 11) based on the GPDEM process) were determined.

Analysis of the Dynamic Reliability of the 3D Slope Stability of Rockfill Dams
In this study, equivalent extremum events were constructed in the safety factor timehistory, a "virtual" time parameter stochastic process was created, and the dynamic stability reliability (i.e., the PDF and CDF curves of the equivalent extremum safety factor (Figure 11) based on the GPDEM process) were determined.
As shown in Figure 11, the traditional joint distribution models ignored the related nonnormal distribution characteristics of the rockfill strength parameters, which led to excessively low calculated failure probabilities and an overestimation of the dynamic reliability of the slope stability.Therefore, for the analysis of the reliability of the stability of a slope based on the coupling randomness of strength parameters and seismic ground motion, the related nonnormal distribution characteristics among the rockfill strength parameters should be fully considered, which is of great significance for the determination of the overall seismic performance, seismic protection, and disaster mitigation of rockfill dams.

Analysis of the Dynamic Reliability of the 3D Slope Stability of Rockfill Dams
In this study, equivalent extremum events were constructed in the safety time-history, a "virtual" time parameter stochastic process was created, and the dy stability reliability (i.e., the PDF and CDF curves of the equivalent extremum safety (Figure 11) based on the GPDEM process) were determined.As shown in Figure 11, the traditional joint distribution models ignored the r nonnormal distribution characteristics of the rockfill strength parameters, which

Discussion
The type of copula function determines the dependency structure of the joint distribution model for rockfill strength, directly affects the selection of sample points, and indirectly leads to differences in calculated reliability.Stochastic factors also affect the computational results of the dynamic reliability of the slope stability of rockfill dams.Hence, the authors of this study analysed and compared the dynamic reliability of the slope stability of rockfill dams from two perspectives, i.e., the copula function type and the different stochastic factors.

Effect of Copula Function Type on the Dynamic Reliability of the Slope Stability of Dams
Based on the 755 sets of strength parameter data presented in Section 3.1 and the six types of copula functions in Table 6, 864 sets of sample data with coupling randomness of strength parameters and seismic ground motion were generated following the steps given in Section 2.6.Then, the slope stability was analysed.Finally, the effect of the type of copula function on the dynamic reliability of the slope stability was analysed in detail.
(1) Effect on Seismic Response of Slope Stability Figures 12 and 13, respectively, show the evolution and isoline distribution contour of the safety factor PDF calculated with both joint distribution models over time within the time interval [8,9] s.As shown in the figures, the evolution of the PDF over time calculated with the different models exhibited a similar overall trend, but there were still significant differences in the peaks and troughs in local areas and obvious differences can be observed in the calculation range and threshold range.The primary cause of these differences was that the type of dependency structure of the copula function was different, which caused the joint distribution model to be different and led to the different distribution characteristics of the sample points, eventually resulting in the dramatic differences in the computational results.(2) Effect on the Dynamic Reliability of the Slope Stability The PDF and CDF curves of the equivalent extremum safety factor (Figure 14) calculated with different models were obtained by constructing equivalent extremum events.
As shown in the figures, the extremums and distribution ranges significantly varied depending on the PDF and CDF calculated by the different models.The t copula model produced the largest result, while the Clayton copula produced the smallest result, and the computational results from the Frank and Plackett copulas were close to each other.
(2) Effect on the Dynamic Reliability of the Slope Stability The PDF and CDF curves of the equivalent extremum safety factor (Figure 14) calculated with different models were obtained by constructing equivalent extremum events.As shown in the figures, the extremums and distribution ranges significantly varied depending on the PDF and CDF calculated by the different models.The t copula model produced the largest result, while the Clayton copula produced the smallest result, and the computational results from the Frank and Plackett copulas were close to each other.
To further compare the effect of copula function type on reliability, the authors of this paper used the joint distribution model built by the Plackett copula function as a reference.The failure probabilities and relative errors (safety factor = 1.0) calculated by the different functions are shown in Table 8, which shows a quantitative analysis of the effect of copula function type on reliability.In conclusion, constructing a joint distribution model of rockfill strength by rationally selecting the optimal copula function is particularly important.The selection of the optimal copula function directly affects the dependency structure of a joint distribution model and may eventually result in differences in the computational results of the dynamic reliability of the slope stability.

Effect of Stochastic Factors on the Dynamic Reliability of the Slope Stability
The dynamic reliability of the slope stability was analysed and compared from three perspectives (i.e., parameter randomness, seismic ground motion randomness, and the coupling randomness of the parameters and seismic ground motion), and the effect of the stochastic factors on the dynamic reliability of the slope stability was systematically studied.To further compare the effect of copula function type on reliability, the authors of this paper used the joint distribution model built by the Plackett copula function as a reference.The failure probabilities and relative errors (safety factor = 1.0) calculated by the different functions are shown in Table 8, which shows a quantitative analysis of the effect of copula function type on reliability.In conclusion, constructing a joint distribution model of rockfill strength by rationally selecting the optimal copula function is particularly important.The selection of the optimal copula function directly affects the dependency structure of a joint distribution model and may eventually result in differences in the computational results of the dynamic reliability of the slope stability.

Effect of Stochastic Factors on the Dynamic Reliability of the Slope Stability
The dynamic reliability of the slope stability was analysed and compared from three perspectives (i.e., parameter randomness, seismic ground motion randomness, and the coupling randomness of the parameters and seismic ground motion), and the effect of the stochastic factors on the dynamic reliability of the slope stability was systematically studied.
Among the structural parameters, only the rockfill strength parameters were considered.A total of 144 sets of sample data were separately generated for each strength parameter, seismic ground motion, and coupling randomness of parameters (432 sets in total); then, the slope stability was analysed; and finally, the effect of different stochastic factor types on the dynamic reliability of the slope stability was analysed in detail.
(1) Effect on the Seismic Response of Slope Stability Figure 15 shows the second-order statistics of the safety factors calculated with different stochastic factors.According to the mean time-history curve, the safety factor calculated by parameter randomness obviously fluctuated over time, and the time-history trends of the safety factors calculated by the coupling randomness and seismic ground motion randomness were close, except that certain differences observed in the intensified stage (3-8 s) of the seismic ground motion.According to the time-history curve of SD, the computational results of the three stochastic factors were quite different from each other.The SD calculated by the seismic ground motion randomness gradually increased from zero and then gradually decreased after reaching the peak, and the initial SD of the parameter randomness was affected by changes in the parameters.If the initial value was large, the fluctuation did not show the development process of gradually increasing from zero to the peak and then decreasing.The trend of SD calculated by coupling randomness was between the computational results calculated by the seismic ground motion randomness and the parameter randomness.
(1) Effect on the Seismic Response of Slope Stability Figure 15 shows the second-order statistics of the safety factors calculated with different stochastic factors.According to the mean time-history curve, the safety factor calculated by parameter randomness obviously fluctuated over time, and the timehistory trends of the safety factors calculated by the coupling randomness and seismic ground motion randomness were close, except that certain differences observed in the intensified stage (3-8 s) of the seismic ground motion.According to the time-history curve of SD, the computational results of the three stochastic factors were quite different from each other.The SD calculated by the seismic ground motion randomness gradually increased from zero and then gradually decreased after reaching the peak, and the initial SD of the parameter randomness was affected by changes in the parameters.If the initial value was large, the fluctuation did not show the development process of gradually increasing from zero to the peak and then decreasing.The trend of SD calculated by coupling randomness was between the computational results calculated by the seismic ground motion randomness and the parameter randomness.Figures 16 and 17 show the evolution and contour map, respectively, of the safety factor PDF calculated with different stochastic factors over time within the time interval [7,8] s.As shown in the figures, the evolution of the PDFs calculated with different stochastic factors greatly varied.The PDF calculated with parameter randomness exhibited a slower evolution, a smaller variation liquidity, a higher probability density, and a more concentrated distribution of thresholds, while the PDFs calculated with seismic ground motion randomness and coupling randomness exhibited a faster evolution and a higher variation liquidity.Compared to seismic ground motion randomness, the coupling randomness caused fewer peak points, a smaller probability density, and a wider threshold distribution range.Figures 16 and 17 show the evolution and contour map, respectively, of the safety factor PDF calculated with different stochastic factors over time within the time interval [7,8] s.As shown in the figures, the evolution of the PDFs calculated with different stochastic factors greatly varied.The PDF calculated with parameter randomness exhibited a slower evolution, a smaller variation liquidity, a higher probability density, and a more concentrated distribution of thresholds, while the PDFs calculated with seismic ground motion randomness and coupling randomness exhibited a faster evolution and a higher variation liquidity.Compared to seismic ground motion randomness, the coupling randomness caused fewer peak points, a smaller probability density, and a wider threshold distribution range.Figure 18 shows the PDF and CDF curves of the safety factor at a typical time (5.2 s) calculated with different stochastic factors.As shown in the figures, for each stochastic factor, the calculated PDF was totally different from the calculated CDF.The PDF calculated with parameter randomness had a high peak and a narrow threshold distribution range; in comparison, the PDFs calculated with seismic ground motion randomness and coupling randomness had a lower peak and a wider threshold distribution range.The trends of seismic ground motion randomness and coupling randomness were found to be close, and differences only existed in the peaks and troughs in local areas.Figure 18 shows the PDF and CDF curves of the safety factor at a typical time (5.2 s) calculated with different stochastic factors.As shown in the figures, for each stochastic factor, the calculated PDF was totally different from the calculated CDF.The PDF calculated with parameter randomness had a high peak and a narrow threshold distribution range; in comparison, the PDFs calculated with seismic ground motion randomness and coupling randomness had a lower peak and a wider threshold distribution range.The trends of seismic ground motion randomness and coupling randomness were found to be close, and differences only existed in the peaks and troughs in local areas.Figure 18 shows the PDF and CDF curves of the safety factor at a typical time (5.2 s) calculated with different stochastic factors.As shown in the figures, for each stochastic factor, the calculated PDF was totally different from the calculated CDF.The PDF calculated with parameter randomness had a high peak and a narrow threshold distribution range; in comparison, the PDFs calculated with seismic ground motion randomness and coupling randomness had a lower peak and a wider threshold distribution range.The trends of seismic ground motion randomness and coupling randomness were found to be close, and differences only existed in the peaks and troughs in local areas.(2) Effect on the Dynamic Reliability of the Slope Stability (2) Effect on the Dynamic Reliability of the Slope Stability The PDF and CDF curves of the equivalent extremum safety factor (Figure 19) calculated with three different stochastic factors were obtained by constructing equivalent extremum events.
As shown in the figures, the PDF and CDF curves were significantly different for the three different stochastic factors.The minimum safety factor calculated with seismic ground motion randomness was the smallest, and its distribution was more concentrated; the safety factor calculated with parameter randomness was the largest, and its threshold distribution range was the largest; and the computational result of coupling randomness was between the other two results.
The PDF and CDF curves of the equivalent extremum safety factor (Figure 19) calculated with three different stochastic factors were obtained by constructing equivalent extremum events.As shown in the figures, the PDF and CDF curves were significantly different for the three different stochastic factors.The minimum safety factor calculated with seismic ground motion randomness was the smallest, and its distribution was more concentrated; the safety factor calculated with parameter randomness was the largest, and its threshold distribution range was the largest; and the computational result of coupling randomness was between the other two results.
In conclusion, all three stochastic factors were found to have an important impact on the safety of the slope stability.Compared to the strength parameter randomness, the seismic ground motion randomness had a more obvious effect; however, the parameter randomness was not negligible, which may have resulted in significant changes in the probability density of the dam slope stability safety factor.

Conclusions
The uncertainty of structural parameters and ground motion is an important factor affecting the stability of a dam slope.The reasonable description of the uncertainty of the structural parameters and ground motion has become a key link in the seismic design of CFRDs.The authors of this study considered the coupling randomness of rockfill strength parameters and seismic ground motion, and they developed a method to analyse the dynamic reliability of the 3D slope stability of rockfill dams.With a concrete-faced rockfill dam as an example, the authors of this paper analysed the differences in the dynamic reliability of dam stability calculated with different joint distribution models and systematically studied the effects of the copula function type and different stochastic factors on the dynamic reliability of dam stability.The following major conclusions were drawn: (1) The copula theory can rationally characterize the nonnormal distribution characteristics related to rockfill strength parameters, and the GF-discrepancy method and the GPDEM can be used to accurately acquire second-order statistics and probability information, which improves the computational efficiency and offers a new approach to analyse the structural dynamic reliability with the coupling randomness of parameters and seismic ground motion.
(2) Traditional joint distribution models ignore the nonnormal distribution characteristics related to rockfill strength parameters, which leads to excessively low calculated failure probabilities and overestimations of the dam slope stability reliability.Therefore, In conclusion, all three stochastic factors were found to have an important impact on the safety of the slope stability.Compared to the strength parameter randomness, the seismic ground motion randomness had a more obvious effect; however, the parameter randomness was not negligible, which may have resulted in significant changes in the probability density of the dam slope stability safety factor.

Conclusions
The uncertainty of structural parameters and ground motion is an important factor affecting the stability of a dam slope.The reasonable description of the uncertainty of the structural parameters and ground motion has become a key link in the seismic design of CFRDs.The authors of this study considered the coupling randomness of rockfill strength parameters and seismic ground motion, and they developed a method to analyse the dynamic reliability of the 3D slope stability of rockfill dams.With a concrete-faced rockfill dam as an example, the authors of this paper analysed the differences in the dynamic reliability of dam stability calculated with different joint distribution models and systematically studied the effects of the copula function type and different stochastic factors on the dynamic reliability of dam stability.The following major conclusions were drawn: (1) The copula theory can rationally characterize the nonnormal distribution characteristics related to rockfill strength parameters, and the GF-discrepancy method and the GPDEM can be used to accurately acquire second-order statistics and probability information, which improves the computational efficiency and offers a new approach to analyse the structural dynamic reliability with the coupling randomness of parameters and seismic ground motion.
(2) Traditional joint distribution models ignore the nonnormal distribution characteristics related to rockfill strength parameters, which leads to excessively low calculated failure probabilities and overestimations of the dam slope stability reliability.Therefore, based on the coupling randomness of rockfill strength parameters and ground motion, the relevant nonnormal distribution characteristics between rockfill strength parameters should be fully considered in the dynamic stability reliability analysis of dam slopes, which has a certain reference significance for the seismic stability safety evaluation and ultimate seismic capacity analysis of concrete-faced dam slopes.
(3) The type of copula function has an important impact on the dynamic reliability of the slope stability of rockfill dams.Therefore, the optimal copula function that characterizes the rockfill strength parameters should be selected as a normal engineering practice to build joint probability distribution models.
(4) The parameter randomness, the seismic ground motion randomness, and the coupling randomness of parameters and seismic ground motion have important impacts on the dynamic reliability of slope stability.The effect of seismic ground motion randomness is more prominent; however, parameter randomness is not negligible, which may cause significant changes in the probability density of the rockfill dam slope stability safety factor.Geotechnical engineering parameters not only have correlations but also strong spatial variability.Subsequent work will combine the random field model to analyse the random dynamic response analysis of rockfill dams with parameter correlation and spatial variability.

Figure 2 .
Figure 2. Analysis process of the dynamic reliability of the 3D slope stability of a rockfill dam based on the coupling randomness of parameters and seismic ground motion.

Figure 2 .
Figure 2. Analysis process of the dynamic reliability of the 3D slope stability of a rockfill dam based on the coupling randomness of parameters and seismic ground motion.
(a) Distribution of rockfill dam projects (b) Strength parameters

Figure 4 .
Figure 4.A total of 124 rockfill dam projects and 767 sets of strength parameters.
(a) Distribution of rockfill dam projects (b) Strength parameters

Figure 4 .
Figure 4.A total of 124 rockfill dam projects and 767 sets of strength parameters.

Figure 4 .
Figure 4.A total of 124 rockfill dam projects and 767 sets of strength parameters.(a) Distribution of rockfill dam projects.(b) Strength parameters.

Figure 5 .
Figure 5. Curve of dam slope safety factor mean under different types of distributions.

Figure 7 .
Figure 7. Evolution of the PDF of the dam slope safety factor.

Figure 8 .
Figure 8. Iso-probability density lines of the dam slope safety factor.

Figure 7 .Figure 7 .
Figure 7. Evolution of the PDF of the dam slope safety factor.(a) Independent normal.(b) Plackett copula.

Figure 8 .
Figure 8. Iso-probability density lines of the dam slope safety factor.

Figures 9 and 10
Figures 9 and 10 show the PDF and cumulative distribution function (CDF) curves, respectively, calculated by the two joint distribution models at typical times (4.1, 8.2, and 10 s).As shown in the figures, the shape of the PDF curve was completely different from that of the CDF curve at different times, and the curves fluctuated and significantly changed.Compared to the joint distribution model based on the optimal copula function, the results calculated by the traditional distribution model were larger and the distribution was wider.The GPDEM and MCS calculation results exhibited considerable consistency, which further demonstrates that the method is accurate.

27 Figure 9 .
Figure 9. PDF curves of the dam slope safety factor at typical times.

Figure 9 .
Figure 9. PDF curves of the dam slope safety factor at typical times.

Figure 9 . 27 Figure 10 .
Figure 9. PDF curves of the dam slope safety factor at typical times.

Figure 10 .
Figure 10.CDF curves of the dam slope safety factor at typical times.

Figure 10 .
Figure 10.CDF curves of the dam slope safety factor at typical times.

Figure 11 .
Figure 11.Equivalent extremum PDF and CDF curves of 3D dam slope safety factors for CFRDs.

Figure 11 .
Figure 11.Equivalent extremum PDF and CDF curves of 3D dam slope safety factors for CFRDs.(a) PDF curve.(b) CDF curve.

Figure 12 .
Figure 12.Evolution of dam slope safety factor PDF.

Figure 13 .
Figure 13.Iso-probability density lines of the dam slope safety factor.

Figure 13 .
Figure 13.Iso-probability density lines of the dam slope safety factor.Figure 13.Iso-probability density lines of the dam slope safety factor.

Figure 13 .
Figure 13.Iso-probability density lines of the dam slope safety factor.Figure 13.Iso-probability density lines of the dam slope safety factor.

Figure 14 .
Figure 14.Equivalent extremum PDF and CDF curves of the dam slope safety factor calculated by different copula functions.

Figure 14 .
Figure 14.Equivalent extremum PDF and CDF curves of the dam slope safety factor calculated by different copula functions.(a) PDF curves.(b) CDF curves.

Figure 15 .
Figure 15.Second-order statistics of dam slope safety factors under different influencing factors.

Figure 15 .
Figure 15.Second-order statistics of dam slope safety factors under different influencing factors.(a) Safety factor mean.(b) Safety factor SD.

Figure 16 .
Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.(a) Seismic ground motion randomness.(b) Coupling randomness.(c) Parameter randomness.

Figure 16 .
Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.(a) Seismic ground motion randomness.(b) Coupling randomness.(c) Parameter randomness.

Figure 16 .Figure 17 .
Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.

Figure 18 .( 2 )
Figure 18.PDF and CDF curves of the safety factor at a typical time calculated with different stochastic factors.

Figure 17 .
Figure 17.Iso-probability density lines of dam slope safety factor under different stochastic factors.(a) Seismic ground motion randomness.(b) Coupling randomness.(c) Parameter randomness.

Figure 16 .Figure 17 .
Figure 16.Evolution of dam slope safety factor PDF under different stochastic factors.

Figure 18 .
Figure 18.PDF and CDF curves of the safety factor at a typical time calculated with different stochastic factors.

Figure 18 .
Figure 18.PDF and CDF curves of the safety factor at a typical time calculated with different stochastic factors.(a) PDF curve.(b) CDF curve.

Figure 19 .
Figure 19.Equivalent extremum PDF and CDF curves of the safety factor calculated with different stochastic.

Figure 19 .
Figure 19.Equivalent extremum PDF and CDF curves of the safety factor calculated with different stochastic.(a) PDF curve.(b) CDF curve.

Table 1 .
Six types of 2D copula functions.

Table 2 .
Static E-B model parameters.

Table 6 .
Calculation results for the optimal marginal distribution function of the rockfill strength parameters.

Table 7 .
Determination of optimal copula function.

Table 8 .
Failure probabilities calculated using different copula functions.

Table 8 .
Failure probabilities calculated using different copula functions.