Copula-Based Slope Reliability Analysis Using the Failure Domain Defined by the g-Line

The estimation of the cross-correlation of shear strength parameters (i.e., cohesion and internal friction angle) and the subsequent determination of the probability of failure have long been challenges in slope reliability analysis. Here, a copula-based approach is proposed to calculate the probability of failure by integrating the copula-based joint probability density function (PDF) on the slope failure domain delimited with the g-line. Here, copulas are used to construct the joint PDF of shear strength parameters with specific marginal distributions and correlation structure. In the paper a failure (limit state) function approach is applied to investigate a system characterized by a homogeneous slope. The results show that the values obtained by using the failure function approach are similar to those calculated by means of conventional methods, such as the first-order reliability method (FORM) and Monte Carlo simulations (MC). In addition, an entropy weight (EW) copula is proposed to address the discrepancies of the results calculated by different copulas to avoid overor underestimating the slope reliability.


Introduction
Numerous problems regarding uncertainty in geotechnical engineering exist currently, such as soil and rock properties, environmental conditions, and theoretical models [1][2][3].The application of reliability analysis for addressing uncertainty in geotechnical design is rather common [4][5][6].To determine the reliability of a slope with high accuracy, it is imperative to investigate the relevance of parameters and establish their joint cumulative distribution function (CDF) or joint probability density function (PDF).However, establishing an accurate joint CDF or PDF remains a practical challenge in geotechnical engineering due to the limited available datasets [7][8][9].Copulas may be used to address this problem [10,11].The use of copulas has the distinct advantage of enabling the construction of a joint CDF of random variables with any form of marginal distribution, overcoming the limitation of conventional approaches [12,13] in which nonnormal variables must be converted into normal variables [9,14].
In recent years, copulas have been applied to different disciplines, including finance [10,15], biomedicine [16,17], hydrology [18][19][20], and coastal engineering [21,22].In geotechnical engineering, Li et al. [23] used copulas to simulate the bivariate probability distribution of two curvefitting parameters underlying the load-displacement models of piles.Tang et al. [9,14] calculated the failure probabilities for an infinite slope and a retaining wall using the direct integration method and investigated the impact of different types of copulas on reliability.Marchant et al. [24] introduced a more general multivariate function for the spatial prediction of soil properties based on copulas.Wu [25,26] determined the failure probability of a slope using a copulabased sampling method.Huang et al. [27] used a copulabased method to estimate the shear strength parameters of rock mass.Motamedi and Liang [28] proposed a probabilistic methodology for landslide hazard assessment at regional scale based on the copula modelling technique considering the possible dependence between the hazard elements.

Mathematical Problems in Engineering
Although copulas have been well developed and widely applied to determine the relevance of material properties and to establish the CDF or PDF, it is difficult to determine the performance function of a slope system in the critical state.Analysing the reliability of a slope remains a significant challenge in practice.
With the help of the -line [29], we propose a reliability analysis method based on copulas, which can be applied to any given homogeneous slopes.The remainder of this paper is organized as follows.Section 2 briefly introduces the copula functions and presents the procedure used in the copulabased slope reliability analysis.In Section 3, after finishing the detailed introductions of -line, we present an approach to identify the failure domain of a slope of any shape.The range of shear strength parameters, which can represent the slope failure domain, is also deduced.Section 4 explains how to calculate the failure probability by directly integrating the joint PDF of the considered variables on the failure domain identified by the -line.Section 5 presents a demonstrative application with evaluation of the failure probability for an example slope.Finally, the results calculated from different copulas are compared and discussed for an entropy weight (EW) copula in Section 6.

Copula Function.
A copula function is a multivariate distribution function with uniform marginal distributions which are all over (0, 1).For an -dimensional uniform random vector , a copula  is expressed as follows [30]: where   and   are the th uniform random variables on the interval (0, 1) and the corresponding values and  represents the probability value.Using Sklar's theorem [31], the joint CDF of an dimensional random variable can be established via the copula function as where  is a copula function,   (  ) is the marginal distribution function of the th random variable and  1 , . . .,   are the copula parameters.Copulas offer the possibility of separately modelling the marginal distributions and the dependence structure among the random variables.
Using a one parameter copula function, the joint CDF of the shear strength parameters  1 and  2 , corresponding to cohesion () and internal friction angle (), can be expressed as follows: in which the  parameter characterizes the dependence structure between the variables  and .

Slope Reliability Assessment Using
Copulas.The establishment of the joint CDF or PDF in terms of copulas for the implementation of the slope reliability analysis involves the following 4 steps.
(1) Estimation of the Marginal Distributions of Each Random Variable.As for fitting copulas for shear strength, the firstline procedure is to determine the appropriate marginal distribution of each random variable.In various previous studies the normal, lognormal, and gamma distributions were selected as the best-fitting marginal distributions for the shear strength parameters ( and ) [25,32,33].The lognormal distribution is defined only for nonnegative values [32].Kolmogorov-Smirnov (K-S) test commonly used to check whether a set of data follows a certain distribution is adopted to examine marginal distribution [27].The key procedure of K-S test is computing the test statistic   , which can be expressed as where   (  ) is the established marginal distribution functions and F (  ) is the cumulative frequency distribution; it is defined as in which    ≤  is the indicator function, equal to 1 if   ≤   and equal to 0 otherwise.  (  ) will be accepted if   <  , , where  , is the critical value of a given significance level . , can be obtained by using the established table [34] when  and  have been determined, and  is the number of observed datasets.
(2) Copula Parameter Estimation.Two-dimensional copulas commonly used in geotechnical engineering belong to two different classes.The elliptical class includes the Gaussian and  copulas.The Archimedean class includes the Frank, Gumbel, and Clayton copulas.Most analyses performed on the shear strength parameters report a negative correlation coefficient between  and , with value between −0.24 and −0.70.However, a positive correlation was also found in several cases [9,25].The shear strength parameters of rocks and soils typically have symmetric structures [9] and asymmetric structures [35,36].Considering the correlation characteristics and symmetry between shear strength parameters, the capability of the Gaussian, Frank, Clayton, Gumbel, and Plackett copulas is tested to model the dependence structure of the data [9,25].The considered copulas are symmetric copulas, except for the Clayton and Gumbel copulas.In addition, negative correlations may not be allowed in some of those copulas; nonetheless, a simple conversion which negates the values of one variable can achieve a positive value for the correlation.
Copula parameters can be obtained by using maximum likelihood estimation (MLE) [30], the Pearson linear correlation coefficient , and the Kendall or Spearman rank correlation coefficient  [11].In this paper, the Kendall rank correlation coefficient  is chosen to estimate the copula parameter, .The parameter is independent of the marginal distribution.For the elliptical class of copulas, the relationship between  and  is given by For the Archimedean class of copulas, Kendall's  is given by where () is the generator function.Kendall's  may be computed by where  is the total number of samples, ,  = 1, 2, . . ., , and sign[⋅] is the sign function, which will be 1 for Then, the corresponding  value is determined by using (8) or (9) depending on the copula class.The parameter of the Plackett copula was determined with the maximum likelihood method [30].
(3) Identification of the Best-Fitting Copula.There are various methods to identify the best-fitting copula model, among these candidate copulas, such as the Akaike information criterion (AIC) [37], root mean square error (RMSE), bias [38], and Cramér-von Mises statistical method [39,40].The AIC is the method most commonly used in engineering practice, and it can be defined as follows: AIC = −2 ⋅ ln (maximized likelihood function for the model) + 2 (no. of fitted parameters) .(11) In this study, AICc [41] takes a correction into consideration for finite sample sizes, which can be determined by the AIC: where  is the number of estimated parameters and  denotes the number of observations (for a single parameter,  = 1).
The copula with the smallest AICc value is preferred.The Cramér-von Mises statistics,   , can be expressed as the sum of the squared differences between the true and empirical copulas for any considered point [39,40]: where    is the copula with the estimated value of the copula parameter,   , and   is the empirical copula given by and  [ 1 ( 1 )≤ 1 ] = 0 otherwise.A parametric bootstrapping approach is used to determine   values [40].The copula with the minimum   and maximum   can be considered as optimal.The RMSE and bias can be determined by the following equations [38]: where    is the th observed value and   is the th estimated value.The best-fitting model is the one that has the smallest RMSE or bias.
(4) Slope Reliability Analysis.Following three successive steps, the joint CDF or PDF of the random variable can be obtained.Subsequently, the probability of failure (  ) can be calculated by integrating the joint PDF established by copulas on the slope failure domain; this approach is called the copulabased direct integration method (CDIM).The determination of the slope failure domain and the integral solution of the probability of failure will be discussed in Sections 3 and 4, respectively.

𝑔-Line and Slope Failure Domain
3.1.-Line.According to Klar et al. [29], the factor of safety (  ) of a slope system can be expressed as follows based on the Mohr-Coulomb criterion: where  and  are the actual shear strength parameters of the soil slope,   and   are the mobilized parameters required for the limit equilibrium state,   denotes the shear strength and   denotes the mobilized shear strength with a factor of safety   , and  is normal stress.Klar et al. [29] proposed the concept of the -line based on (17) to determine the factor of safety geometrically.The -line (Figure 1) represents the combinations of  and  at which the slope will exist in limit equilibrium (  = 1).The combinations of  and   above the -line imply a stable condition (  > 1), whereas the combinations falling below the -line (namely, in the failure domain symbolized by Σ) are impossible and represent nonphysical conditions since the system cannot mobilize more than its capacity.The geometrical expression of   can be written as in which || denotes the distance between point  and the origin and the distance between the origin and point  is designated by ||. is an arbitrary point in space -, and  is the intersection of the -line and .
In the present paper, the simplified Bishop method is used to calculate the factor of safety.Using the software SLOPE/W, hundreds of runs are made to search the critical state (  = 1) of a slope with given geometry and unit weight, and then the corresponding  and  values in the limit equilibrium are obtained.The -line is fitted by the combinations of  and  in the limit equilibrium state.
For any given shape of homogeneous slope, the - values in the limit equilibrium state can always be used to determine the -line in a unique manner when the slope geometry characteristic and unit weight are given.In other words, the -line can be regarded as an inherent characteristic of the slope used to characterize the slope's limit equilibrium state.The -line makes it convenient to conduct the reliability analysis of the slope.

𝑔-Line
Fitting Shape.The -line determines the failure domain of the slope, and its shape also directly affects the results of the slope reliability analysis.However, there are few studies regarding the shape of the -line.In this paper, several slopes with different slope heights (), slope angles (), and unit weights () are used to calculate the  and tan  values in the limit equilibrium, as shown in Table 1. Figure 2(a) and Table 1 clearly indicate that the quadratic polynomial is effective for fitting the -line for general homogeneous slopes, and its coefficient of determination ( 2 ) is generally in the range of 0.92 to 0.99.Considering the  2 value of the quadratic polynomial as an objective reference, the slope angle is more sensitive than the unit weight, for a same  value.In contrast, the sensitivity of the slope height is relatively low.Quadratic polynomial fitting has a relatively poor performance when considering high and steep slopes, as shown in Figure 2; in contrast, higher-degree polynomials exhibited a better fitting effect.

Failure Domain Represented by the Shear Strength Parameters.
Figure 1 shows the -line and failure domain Σ of a specific slope, where  = 30 m,  = 45 ∘ , and  = 20 kN/m 3 .The -line can be expressed in terms of quadratic expressions,  =  2 +  +  ( > 0,  < 0).Similarly, the failure domain Σ can be expressed as  ≤  2 +  + .For  = 0, the corresponding  value reaches a maximum (i.e.,  max = ).Similarly, the failure domain Σ is equivalent to  ≤ [− − √ 2 − 4( − )]/2,  > 0 and  < 0, if and only if  = 0, the  value is the maximum, and the maximum value of the cohesive force  can be expressed as  max = [− − √  2 − 4]/2.The range of shear strength parameters in the failure domain must be determined to compute the probability of failure.

Application of the Direct Integration Method
The probability of failure can be accurately determined through the direct integration approach [9], which can mitigate the errors inherent in many reliability methods, such as the first-order reliability method (FORM), second-order reliability method (SORM), and Monte Carlo simulation (MCS).In this study, the direct integration approach is adopted to integrate the copula-based - joint PDF for the -line failure domain (CDIM).Thus, the probability of failure of a slope with any complex shape can be calculated by CDIM.
Based on the limit state, the performance function of slope for the -line can be expressed as follows: in which , , and  are quadratic polynomial coefficients,  is the cohesion, and  is the internal friction angle.The limit state equation of the -line corresponds to  = 0, and the failure domain Σ corresponds to  < 0 (  < 1); thus, the probability of failure can obtained as By substituting the copula-based joint probability density function given by ( 4) into (20), the expression of   can be transformed as Because of the complexity of the copula function, the above double-integral calculation is often difficult to perform, (21) can be converted into a single integral by defining first derivative ( 1 ,  2 ; ) = C( 1 ,  2 ; )/ 1 of the copula [9], and here  1 and  2 denote  1 () and  2 (tan ), respectively, such that By substituting the range of shear strength parameters in the -line failure domain into (22), the final expression of solving   is

Slope Model and Its 𝑔-Line Failure Domain.
A hypothetical homogeneous slope shown in Figure 3 is employed to investigate the probability of failure using the CDIM, with the slope height () of 30.0 m, the slope angle () of 30.0 ∘ , and the unit weight () of 18 kN/m 3 .Figure 3 shows the -line for the slope model and its expression of  the quadratic polynomial fitting ( 2 = 0.9788).From the analysis in Section 3.3, the parameter range of cohesion  characterizing the slope failure domain is 0 <  <  max and  max = 139.8kPa.

Datasets.
In this study, 63 datasets of measured shear strength parameters collected from the Xiaolangdi hydropower station [42] are used to examine the marginal distributions and to fit the joint behaviour of the soil shear strength.The K-S test results are listed in Table 2. Obviously,  and tan  coincide, respectively, with the normal distribution and lognormal distribution well.The gamma distribution is rejected.Figure 4 From these 63 datasets, we can conclude that the mean and coefficient of variation of  are 65.97 kPa and 0.3, respectively.Similarly, the average value of tan  is 0.42, with the coefficient of variation of 0.15.The Pearson correlation coefficient between  and tan  is −0.5; correspondingly, Kendall's  is −0.4035.In the analysis, the shear strength parameters,  and tan , are selected as random variables.The slope shape parameters and unit weight are assumed as deterministic.
Considering some copulas that do not exhibit negative parameter correlations, the variable tan  can be transformed into  = 2mean(tan ) − tan , where mean(⋅) is the average value.It is worth keeping in mind that both  and tan  have the same mean and standard deviation.Also the correlation coefficient between  and  will have the same magnitude but opposite in direction to that of  and tan .The scatter diagram of  and  is shown in Figure 6.To compute the probability of failure, the -line and its failure domain must be correspondingly transformed.The shaded part in Figure 5 shows the corresponding failure domain of the transformed -line.

Probability of Failure Using the CDIM.
Based on those measured datasets, (10) is employed for the aforementioned copulas to obtain Kendall's rank correlation coefficient , and then the parameters of the copulas can be calculated using ( 8), (9), and the maximum likelihood method (MLE) [30], as shown in Table 3.The setting for the optimum copula is the minimum   and maximum   .In this paper,   /  is introduced to evaluate the goodness of fit of the copulas.The smallest value of   /  indicates the optimum copula function.Using similar methods proposed by different researchers [38,39,41], the AICc,   /  , RMSE, and bias of the copulas are computed by ( 11)-( 16), respectively, and tabulated in Table 3.The Frank copula exhibited smaller AICc,   /  , and bias values, whereas the Gaussian copula exhibited smaller RMSE values.The optimum copula is not unique for different evaluation approaches.Figure 6 shows the contour plots of the joint PDFs for the transformed shear strength parameters associated with different copulas.The following points can be made following a careful investigation of Table 3 and Figure 6.Frank copulas showed the best fits, and the Gaussian, Gumbel, and Plackett copulas demonstrated better fits than the Clayton copulas.
Using the CDIM, one can determine the probability of slope failure by substituting the correlation parameter  (from Table 3) into (23).The corresponding probability of failure is calculated by using FORM and MCS, considering the performance function represented by the -line, as shown in Table 4.In the analyses of MCS, 10000 random numbers of  and tan  are generated by using the method proposed by Ilich [43].The probabilities of failure computed using the CDIM, FORM, and MCS yield similar results.The maximum relative error values for the Clayton copula and Gumbel copula are 7.67% and 5.42%, respectively.

Discussion
As previously discussed, the optimum copula is not unique for different evaluation methods.Therefore, AICc,   /  , RMSE, and bias are considered as the evaluation index, and the weights of the selected five copulas are calculated and tabulated in Table 5 using the entropy analysis method [44].The bivariate EW copula for the shear strength parameters ( and tan ) can be expressed in terms of the weighted combination and density function as follows: where   ,   , and   are the th selected weight, density function, and parameter of the copula function, respectively,  = 1, . . ., 5. Similarly, the term  EW in ( 24) is integrated for the -line failure domain.Consequently, the probability of failure of a slope system is determined for the considered EW copula.In addition, according to the example slope model, fivegroup factors of safety with different  and tan  values under the same coefficient of variation are given in Table 5 to further investigate the variation characteristics of the probability of failure with the factor of safety under different copulas.The probabilities of failure of the example slope under different copulas, including the EW copula, and the ratios () of the probabilities of failure calculated by the five copulas compared to that calculated by the EW copula are also shown in Table 5.For a more intuitive comparison, Figure 7 shows a bar graph of the probability of failure of the slope corresponding to each copula function under different factors of safety.
The probability of failure clearly decreases with increases in the factor of safety, but the probability of failure differs notably between each copula function with different factors of safety.The computed results of the different copulas show insignificant discrepancies for a relatively small factor of safety.Nonetheless, the variation increased rapidly for higher factors of safety.For example, the probability of failure of the Plackett copula is 9-fold higher than that of the Clayton copula (when   = 1.615); in addition, the value of  gradually deviated from "1" and showed progressive discrepancies.For a higher factor of safety, the Gaussian and Clayton copulas demonstrated relatively small probabilities of failure.In contrast, the Frank and Plackett copulas exhibited relatively higher probabilities of failure.Therefore, the Frank and Plackett copulas are relatively conservative, whereas the Gaussian and Clayton copulas overestimate the reliability of the slope.
The different structures of the copulas lead to different probabilities of failure.In particular, for a low probability of failure (i.e., high factor of safety), the result of the reliability analysis shows greater sensitivity to the type of copula function.The EW copula combines the results of different copula functions and effectively controls the discrepancies between the probabilities of failure of the slopes.Thus, Factor of safety (F s ) Figure 7:   and   of the illustrative slope calculated using different copulas.
Mathematical Problems in Engineering the EW copula also can avoid over-or underestimating the slope reliability.

Conclusions
In this paper, a CDIM was used to compute the probability of failure of a slope system.The failure domain of a homogeneous slope was determined by using the -line.In addition, the bivariate joint CDFs or PDFs for shear strength parameters were established based on the copula function.
Based on the main findings of this paper, the following conclusions are noted: (1) For a certain homogeneous slope, the -line characterizes its limit equilibrium state.This phenomenon can be regarded as the inherent property of the slope.A quadratic polynomial can accurately fit the homogeneous slope -line, and the coefficient of determination  2 ranges from 0.92 to 0.99.Higherorder polynomials exhibit a better fitting effect for high and steep slopes.The slope failure domain determined using the -line can be represented by shear strength parameters ( or ).
(2) The computed probabilities of failure using the CDIM, FORM, and MCS yield highly similar results.The maximum relative error values for the Clayton copula and Gumbel copula are 7.67% and 5.42%, respectively.This result validated the rationality of CDIM method.The CDIM provides a further approach to slope reliability analysis.However, further research is required to compute the probability of failure for heterogeneous slopes.
(3) The computed results of the copula functions showed minor discrepancies for relatively small factors of safety.Nonetheless, the variation increased rapidly for higher values of factor of safety.The EW copula combines the results of different copula functions and effectively controls the discrepancies between the probability of failure of slopes.The EW copula can also eliminate bias in slope reliability estimates.

Figure 1 :
Figure 1: Schematics of   and   defined based on the -line.

Figure 3 :
Figure 3: Geometry of the illustrative slope and its -line.
(a) shows the scatter plot of datasets, and the PDFs and CDFs of different distribution forms are given in Figures 4(b), 4(c), 4(d), and 4(e).

Figure 4 :
Figure 4: Scatter plot, PDFs, and CDFs of 63 measured datasets.(a) Scatter plot of the measured shear strength parameters, (b) bar chart and different PDFs for , (c) bar chart and different PDFs for tan , (d) different CDFs for , and (e) different CDFs for tan .

Figure 6 :
Figure 6: Contour plots of the joint PDFs for the transformed shear strength parameters ( and ) associated with different copulas.

Table 1 :
Statistics of  2 of the -line fitting in the form of a quadratic polynomial under different conditions.

Table 4 :
Probabilities of failure and their relative error of the example slope.

Table 5 :
and  values of the example slope calculated for different   values.