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

Authors: Xiaoliang Xu; Jianlin Li; Jiawei Gong; Huafeng Deng; Liangpeng Wan
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

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 theg-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 over- or underestimating the slope reliability.


---

## Body

## 1. Introduction

Numerous problems regarding uncertainty in geotechnical engineering exist currently, such as soil and rock properties, environmental conditions, and theoretical models [1–3]. The application of reliability analysis for addressing uncertainty in geotechnical design is rather common [4–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–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–20], and coastal engineering [21, 22]. In geotechnical engineering, Li et al. [23] used copulas to simulate the bivariate probability distribution of two curve-fitting 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 copula-based sampling method. Huang et al. [27] used a copula-based 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.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 theg-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 copula-based slope reliability analysis. In Section 3, after finishing the detailed introductions of g-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 g-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.


## 2. Copulas for the Shear Strength Parameters

### 2.1. Copula Function

A copula function is a multivariate distribution function with uniform marginal distributions which are all over(0,1). For an n-dimensional uniform random vector U, a copula C is expressed as follows [30]:(1)Cu1,…,un=PU1≤u1,…,Un≤un,where Ui and ui are the ith uniform random variables on the interval (0,1) and the corresponding values and P represents the probability value.Using Sklar’s theorem [31], the joint CDF of an n-dimensional random variable can be established via the copula function as(2)Hx1,…,xn=CF1x1,…,Fnxn;θ1,…,θq,where C is a copula function, Fi(xi) is the marginal distribution function of the ith random variable and θ1,…,θq 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 parametersX1 and X2, corresponding to cohesion (c) and internal friction angle (φ), can be expressed as follows:(3)Hx1,x2=CF1x1,F2x2;θin which the θ parameter characterizes the dependence structure between the variables c and φ.The corresponding joint PDF is(4)fx1,x2=f1x1f2x2cF1x1,F2x2;θ,where f1(x1) and f2(x2) are the PDFs of X1 and X2 and c(F1(x1),F2(x2);θ) is the copula density function given by(5)cF1x1,F2x2;θ=cu1,u2θ=∂2Cu1,u2θ∂u1∂u2.


### 2.2. 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 first-line 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 (c 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 Dn, which can be expressed as(6)Dn=max⁡Fixi-F~ixi,where Fixi is the established marginal distribution functions and F~ixi is the cumulative frequency distribution; it is defined as(7)F~ixi=1n∑i=1nIxij≤xiin which Ixij≤xi is the indicator function, equal to 1 if xij≤xi and equal to 0 otherwise. Fixi will be accepted if Dn<Dn,α,where Dn,α is the critical value of a given significance level α. Dn,α can be obtained by using the established table [34] when n and α have been determined, and n 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 t copulas. The Archimedean class includes the Frank, Gumbel, and Clayton copulas.Most analyses performed on the shear strength parameters report a negative correlation coefficient betweenc 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(8)τ=2πarcsin⁡θ.For the Archimedean class of copulas, Kendall’s τ is given by(9)τ=1+4∫01φtφ′t,where φ(t) is the generator function. Kendall’s τ may be computed by (10)τ=N2-1∑i<jsign⁡x1i-x1jx2i-x2j,where N is the total number of samples, i,j=1,2,…,N, and sign⁡[·] is the sign function, which will be 1 for x1i-x1jx2i-x2j>0 or -1 in all other cases.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:(11)AIC=-2ln⁡maximized likelihood function for the model+2no. of fitted parameters.In this study, AICc [41] takes a correction into consideration for finite sample sizes, which can be determined by the AIC:(12)AICc=AIC+2kk+1n-k-1,where k is the number of estimated parameters and n denotes the number of observations (for a single parameter, k=1). The copula with the smallest AICc value is preferred. The Cramér-von Mises statistics, Sn, can be expressed as the sum of the squared differences between the true and empirical copulas for any considered point [39, 40]:(13)Sn=∑i=1nCnU1i,U2i-CθnU1i,U2i2,where Cθn is the copula with the estimated value of the copula parameter, θn, and Cn is the empirical copula given by(14)Cnu1,u2=1n∑i=1nIF1x1i≤u1IF2x2i≤u2in which F1x1i and F2x2i are empirical CDFs, I[·] is the indicator function, if F1x1i≤u1, I[F1x1i≤u1]=1, and I[F1x1i≤u1]=0 otherwise. A parametric bootstrapping approach is used to determine pc values [40]. The copula with the minimum pc and maximum Sn can be considered as optimal. The RMSE and bias can be determined by the following equations [38]:(15)RMSE=1n-1∑i=1nPei-Pi2(16)Bias=1N∑i=1nPei-Pi,where Pei is the ith observed value and Pi is the ith 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 (Pf) can be calculated by integrating the joint PDF established by copulas on the slope failure domain; this approach is called the copula-based 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.


## 2.1. Copula Function

A copula function is a multivariate distribution function with uniform marginal distributions which are all over(0,1). For an n-dimensional uniform random vector U, a copula C is expressed as follows [30]:(1)Cu1,…,un=PU1≤u1,…,Un≤un,where Ui and ui are the ith uniform random variables on the interval (0,1) and the corresponding values and P represents the probability value.Using Sklar’s theorem [31], the joint CDF of an n-dimensional random variable can be established via the copula function as(2)Hx1,…,xn=CF1x1,…,Fnxn;θ1,…,θq,where C is a copula function, Fi(xi) is the marginal distribution function of the ith random variable and θ1,…,θq 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 parametersX1 and X2, corresponding to cohesion (c) and internal friction angle (φ), can be expressed as follows:(3)Hx1,x2=CF1x1,F2x2;θin which the θ parameter characterizes the dependence structure between the variables c and φ.The corresponding joint PDF is(4)fx1,x2=f1x1f2x2cF1x1,F2x2;θ,where f1(x1) and f2(x2) are the PDFs of X1 and X2 and c(F1(x1),F2(x2);θ) is the copula density function given by(5)cF1x1,F2x2;θ=cu1,u2θ=∂2Cu1,u2θ∂u1∂u2.


## 2.2. 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 first-line 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 (c 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 Dn, which can be expressed as(6)Dn=max⁡Fixi-F~ixi,where Fixi is the established marginal distribution functions and F~ixi is the cumulative frequency distribution; it is defined as(7)F~ixi=1n∑i=1nIxij≤xiin which Ixij≤xi is the indicator function, equal to 1 if xij≤xi and equal to 0 otherwise. Fixi will be accepted if Dn<Dn,α,where Dn,α is the critical value of a given significance level α. Dn,α can be obtained by using the established table [34] when n and α have been determined, and n 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 t copulas. The Archimedean class includes the Frank, Gumbel, and Clayton copulas.Most analyses performed on the shear strength parameters report a negative correlation coefficient betweenc 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(8)τ=2πarcsin⁡θ.For the Archimedean class of copulas, Kendall’s τ is given by(9)τ=1+4∫01φtφ′t,where φ(t) is the generator function. Kendall’s τ may be computed by (10)τ=N2-1∑i<jsign⁡x1i-x1jx2i-x2j,where N is the total number of samples, i,j=1,2,…,N, and sign⁡[·] is the sign function, which will be 1 for x1i-x1jx2i-x2j>0 or -1 in all other cases.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:(11)AIC=-2ln⁡maximized likelihood function for the model+2no. of fitted parameters.In this study, AICc [41] takes a correction into consideration for finite sample sizes, which can be determined by the AIC:(12)AICc=AIC+2kk+1n-k-1,where k is the number of estimated parameters and n denotes the number of observations (for a single parameter, k=1). The copula with the smallest AICc value is preferred. The Cramér-von Mises statistics, Sn, can be expressed as the sum of the squared differences between the true and empirical copulas for any considered point [39, 40]:(13)Sn=∑i=1nCnU1i,U2i-CθnU1i,U2i2,where Cθn is the copula with the estimated value of the copula parameter, θn, and Cn is the empirical copula given by(14)Cnu1,u2=1n∑i=1nIF1x1i≤u1IF2x2i≤u2in which F1x1i and F2x2i are empirical CDFs, I[·] is the indicator function, if F1x1i≤u1, I[F1x1i≤u1]=1, and I[F1x1i≤u1]=0 otherwise. A parametric bootstrapping approach is used to determine pc values [40]. The copula with the minimum pc and maximum Sn can be considered as optimal. The RMSE and bias can be determined by the following equations [38]:(15)RMSE=1n-1∑i=1nPei-Pi2(16)Bias=1N∑i=1nPei-Pi,where Pei is the ith observed value and Pi is the ith 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 (Pf) can be calculated by integrating the joint PDF established by copulas on the slope failure domain; this approach is called the copula-based 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.


## 3.g-Line and Slope Failure Domain

### 3.1.g-Line

According to Klar et al. [29], the factor of safety (Fs) of a slope system can be expressed as follows based on the Mohr-Coulomb criterion:(17)Fs=τfτm=c+σtan⁡φcm+σtan⁡φm=tan⁡φtan⁡φm=ccm,where c and φ are the actual shear strength parameters of the soil slope, cm and φm are the mobilized parameters required for the limit equilibrium state, τf denotes the shear strength and τm denotes the mobilized shear strength with a factor of safety Fs, and σ is normal stress. Klar et al. [29] proposed the concept of the g-line based on (17) to determine the factor of safety geometrically. The g-line (Figure 1) represents the combinations of c and φ at which the slope will exist in limit equilibrium (Fs=1). The combinations of c and φ above the g-line imply a stable condition (Fs>1), whereas the combinations falling below the g-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 Fs can be written as(18)Fs=OAOB=xaxb=yaybin which |OA| denotes the distance between point A and the origin and the distance between the origin and point B is designated by |OB|. A is an arbitrary point in space c-φ, and B is the intersection of the g-line and OA.Figure 1
Schematics ofFs and Pf defined based on the g-line.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 (Fs=1) of a slope with given geometry and unit weight, and then the corresponding c and φ values in the limit equilibrium are obtained. The g-line is fitted by the combinations of c and φ in the limit equilibrium state.For any given shape of homogeneous slope, thec-φ values in the limit equilibrium state can always be used to determine the g-line in a unique manner when the slope geometry characteristic and unit weight are given. In other words, the g-line can be regarded as an inherent characteristic of the slope used to characterize the slope’s limit equilibrium state. The g-line makes it convenient to conduct the reliability analysis of the slope.


### 3.2.g-Line Fitting Shape

Theg-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 g-line. In this paper, several slopes with different slope heights (H), slope angles (β), and unit weights (γ) are used to calculate the c 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 g-line for general homogeneous slopes, and its coefficient of determination (R2) is generally in the range of 0.92 to 0.99. Considering the R2 value of the quadratic polynomial as an objective reference, the slope angle is more sensitive than the unit weight, for a same H 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.Table 1
Statistics ofR2 of the g-line fitting in the form of a quadratic polynomial under different conditions.
H (m) γ = 18 kN/m3 γ = 20 kN/m3 γ = 22 kN/m3 β = 30° β = 45° β = 6 0° β = 30° β = 45° β = 6 0° β = 30° β = 45° β = 6 0° 30 0.9868 0.9786 0.9441 0.9873 0.9790 0.9470 0.9877 0.9796 0.9483 60 0.9878 0.9788 0.9411 0.9882 0.9792 0.9436 0.9849 0.9735 0.9304 90 0.9844 0.9747 0.9118 0.9894 0.9857 0.9345 0.9862 0.9770 0.9148Figure 2
g-line curve-fitting forms. (a) g-line fitting by quadratic polynomial at different slope angles. ((b), (c), and (d)) g-line curve-fitting with different powers for changing slope angle β.
(a) (b) (c) (d)


### 3.3. Failure Domain Represented by the Shear Strength Parameters

Figure1 shows the g-line and failure domain Σ of a specific slope, where H=30 m, β = 45°, and γ = 20 kN/m3. The g-line can be expressed in terms of quadratic expressions, y=Ax2+Bx+C(A>0,B<0). Similarly, the failure domain Σ can be expressed as y≤Ax2+Bx+C. For x=0, the corresponding y value reaches a maximum (i.e., φmax=β). Similarly, the failure domain Σ is equivalent to x≤[-B-B2-4AC-y]/2A,  A>0 and B<0, if and only if y=0, the x value is the maximum, and the maximum value of the cohesive force c can be expressed as cmax=[-B-B2-4AC]/2A. The range of shear strength parameters in the failure domain must be determined to compute the probability of failure.


## 3.1.g-Line

According to Klar et al. [29], the factor of safety (Fs) of a slope system can be expressed as follows based on the Mohr-Coulomb criterion:(17)Fs=τfτm=c+σtan⁡φcm+σtan⁡φm=tan⁡φtan⁡φm=ccm,where c and φ are the actual shear strength parameters of the soil slope, cm and φm are the mobilized parameters required for the limit equilibrium state, τf denotes the shear strength and τm denotes the mobilized shear strength with a factor of safety Fs, and σ is normal stress. Klar et al. [29] proposed the concept of the g-line based on (17) to determine the factor of safety geometrically. The g-line (Figure 1) represents the combinations of c and φ at which the slope will exist in limit equilibrium (Fs=1). The combinations of c and φ above the g-line imply a stable condition (Fs>1), whereas the combinations falling below the g-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 Fs can be written as(18)Fs=OAOB=xaxb=yaybin which |OA| denotes the distance between point A and the origin and the distance between the origin and point B is designated by |OB|. A is an arbitrary point in space c-φ, and B is the intersection of the g-line and OA.Figure 1
Schematics ofFs and Pf defined based on the g-line.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 (Fs=1) of a slope with given geometry and unit weight, and then the corresponding c and φ values in the limit equilibrium are obtained. The g-line is fitted by the combinations of c and φ in the limit equilibrium state.For any given shape of homogeneous slope, thec-φ values in the limit equilibrium state can always be used to determine the g-line in a unique manner when the slope geometry characteristic and unit weight are given. In other words, the g-line can be regarded as an inherent characteristic of the slope used to characterize the slope’s limit equilibrium state. The g-line makes it convenient to conduct the reliability analysis of the slope.


## 3.2.g-Line Fitting Shape

Theg-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 g-line. In this paper, several slopes with different slope heights (H), slope angles (β), and unit weights (γ) are used to calculate the c 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 g-line for general homogeneous slopes, and its coefficient of determination (R2) is generally in the range of 0.92 to 0.99. Considering the R2 value of the quadratic polynomial as an objective reference, the slope angle is more sensitive than the unit weight, for a same H 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.Table 1
Statistics ofR2 of the g-line fitting in the form of a quadratic polynomial under different conditions.
H (m) γ = 18 kN/m3 γ = 20 kN/m3 γ = 22 kN/m3 β = 30° β = 45° β = 6 0° β = 30° β = 45° β = 6 0° β = 30° β = 45° β = 6 0° 30 0.9868 0.9786 0.9441 0.9873 0.9790 0.9470 0.9877 0.9796 0.9483 60 0.9878 0.9788 0.9411 0.9882 0.9792 0.9436 0.9849 0.9735 0.9304 90 0.9844 0.9747 0.9118 0.9894 0.9857 0.9345 0.9862 0.9770 0.9148Figure 2
g-line curve-fitting forms. (a) g-line fitting by quadratic polynomial at different slope angles. ((b), (c), and (d)) g-line curve-fitting with different powers for changing slope angle β.
(a) (b) (c) (d)


## 3.3. Failure Domain Represented by the Shear Strength Parameters

Figure1 shows the g-line and failure domain Σ of a specific slope, where H=30 m, β = 45°, and γ = 20 kN/m3. The g-line can be expressed in terms of quadratic expressions, y=Ax2+Bx+C(A>0,B<0). Similarly, the failure domain Σ can be expressed as y≤Ax2+Bx+C. For x=0, the corresponding y value reaches a maximum (i.e., φmax=β). Similarly, the failure domain Σ is equivalent to x≤[-B-B2-4AC-y]/2A,  A>0 and B<0, if and only if y=0, the x value is the maximum, and the maximum value of the cohesive force c can be expressed as cmax=[-B-B2-4AC]/2A. The range of shear strength parameters in the failure domain must be determined to compute the probability of failure.


## 4. 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 c-f joint PDF for the g-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 theg-line can be expressed as follows:(19)g=Ac2+Bc+C-tan⁡φin which A, B, and C are quadratic polynomial coefficients, c is the cohesion,and φ is the internal friction angle. The limit state equation of the g-line corresponds to g=0, and the failure domain Σ corresponds to g<0   (Fs<1); thus, the probability of failure can obtained as(20)Pf=∬Σfc,fdcdf.By substituting the copula-based joint probability density function given by (4) into (20), the expression of Pf can be transformed as(21)Pf=∬gc,tan⁡φ≤0f1cf2tan⁡φcF1c,F2tan⁡φ;θdcdtan⁡φ.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 Lu1,u2;θ=∂Cu1,u2;θ/∂u1 of the copula [9], and here u1 and u2 denote F1(c) and F2(tan⁡φ), respectively, such that(22)Pf=∫Σf1cLF1c,F2tan⁡φ;θdc.By substituting the range of shear strength parameters in theg-line failure domain into (22), the final expression of solving Pf is(23)Pf=∫0cmaxf1cLF1c,F2Ac2+Bc+C;θdc.


## 5. Illustrative Example

### 5.1. Slope Model and Itsg-Line Failure Domain

A hypothetical homogeneous slope shown in Figure3 is employed to investigate the probability of failure using the CDIM, with the slope height (H) of 30.0 m, the slope angle (β) of 30.0°, and the unit weight (γ) of 18 kN/m3. Figure 3 shows the g-line for the slope model and its expression of the quadratic polynomial fitting (R2 = 0.9788). From the analysis in Section 3.3, the parameter range of cohesion c characterizing the slope failure domain is 0<c<cmax and cmax = 139.8 kPa.Figure 3
Geometry of the illustrative slope and itsg-line.


### 5.2. 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, c and tan⁡φ coincide, respectively, with the normal distribution and lognormal distribution well. The gamma distribution is rejected. Figure 4(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).Table 2
Results of Kolmogorov-Smirnov test.
Variables n α D n , α D n Normal Lognormal Gamma c 63 0.05 0.1682 0.0656 0.1293 1 tan ⁡ φ 63 0.05 0.1682 0.1444 0.1224 0.9987Figure 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 forc, (c) bar chart and different PDFs for tan⁡φ, (d) different CDFs for c, and (e) different CDFs for tan⁡φ.
(a) (b) (c) (d) (e)From these 63 datasets, we can conclude that the mean and coefficient of variation ofc 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 c and tan⁡φ is −0.5; correspondingly, Kendall’s τ is −0.4035. In the analysis, the shear strength parameters, c 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 variabletan⁡φ 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 c and η will have the same magnitude but opposite in direction to that of c and tan⁡φ. The scatter diagram of c and η is shown in Figure 6. To compute the probability of failure, the g-line and its failure domain must be correspondingly transformed. The shaded part in Figure 5 shows the corresponding failure domain of the transformed g-line.Figure 5
Transformedg-line and failure domain.Figure 6
Contour plots of the joint PDFs for the transformed shear strength parameters (c and η) associated with different copulas.


### 5.3. 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 Sn and maximum pc. In this paper, Sn/pc is introduced to evaluate the goodness of fit of the copulas. The smallest value of Sn/pc indicates the optimum copula function. Using similar methods proposed by different researchers [38, 39, 41], the AICc, Sn/pc, RMSE, and bias of the copulas are computed by (11)–(16), respectively, and tabulated in Table 3. The Frank copula exhibited smaller AICc, Sn/pc, 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.Table 3
Values of the AICc,Sn/pc, RMSE, bias, Kendall’s τ, and copula parameter θ.
Family Gaussian Frank Clayton Gumbel Plackett Kendall’s coefficientτ Copula parameterθ 0.5922 4.2092 1.3528 1.6764 5.6180 0.4035 AICc - 41.0086 - 42.8635 - 36.9945 - 41.3524 - 41.7749 S n/pc 0.0437 0.0328 0.0502 0.2542 0.0367 RMSE 0.2203 0.2212 0.2298 0.2293 0.2282 Bias 0.1923 0.1867 0.2178 0.1940 0.1925Using 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 g-line, as shown in Table 4. In the analyses of MCS, 10000 random numbers of c 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.Table 4
Probabilities of failure and their relative error of the example slope.
Method Copula-based direct integration method (CDIM) FORM MCS10000 random numbers Gaussian Frank Clayton Gumbel Plackett P f 0.2606 0.2531 0.2668 0.2529 0.2571 0.2478 0.2674 Relative error (%) 5.17 2.14 7.67 2.06 3.75 0 / 2.54 5.35 0.22 5.42 3.85 / 0


## 5.1. Slope Model and Itsg-Line Failure Domain

A hypothetical homogeneous slope shown in Figure3 is employed to investigate the probability of failure using the CDIM, with the slope height (H) of 30.0 m, the slope angle (β) of 30.0°, and the unit weight (γ) of 18 kN/m3. Figure 3 shows the g-line for the slope model and its expression of the quadratic polynomial fitting (R2 = 0.9788). From the analysis in Section 3.3, the parameter range of cohesion c characterizing the slope failure domain is 0<c<cmax and cmax = 139.8 kPa.Figure 3
Geometry of the illustrative slope and itsg-line.


## 5.2. 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, c and tan⁡φ coincide, respectively, with the normal distribution and lognormal distribution well. The gamma distribution is rejected. Figure 4(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).Table 2
Results of Kolmogorov-Smirnov test.
Variables n α D n , α D n Normal Lognormal Gamma c 63 0.05 0.1682 0.0656 0.1293 1 tan ⁡ φ 63 0.05 0.1682 0.1444 0.1224 0.9987Figure 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 forc, (c) bar chart and different PDFs for tan⁡φ, (d) different CDFs for c, and (e) different CDFs for tan⁡φ.
(a) (b) (c) (d) (e)From these 63 datasets, we can conclude that the mean and coefficient of variation ofc 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 c and tan⁡φ is −0.5; correspondingly, Kendall’s τ is −0.4035. In the analysis, the shear strength parameters, c 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 variabletan⁡φ 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 c and η will have the same magnitude but opposite in direction to that of c and tan⁡φ. The scatter diagram of c and η is shown in Figure 6. To compute the probability of failure, the g-line and its failure domain must be correspondingly transformed. The shaded part in Figure 5 shows the corresponding failure domain of the transformed g-line.Figure 5
Transformedg-line and failure domain.Figure 6
Contour plots of the joint PDFs for the transformed shear strength parameters (c and η) associated with different copulas.


## 5.3. 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 Sn and maximum pc. In this paper, Sn/pc is introduced to evaluate the goodness of fit of the copulas. The smallest value of Sn/pc indicates the optimum copula function. Using similar methods proposed by different researchers [38, 39, 41], the AICc, Sn/pc, RMSE, and bias of the copulas are computed by (11)–(16), respectively, and tabulated in Table 3. The Frank copula exhibited smaller AICc, Sn/pc, 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.Table 3
Values of the AICc,Sn/pc, RMSE, bias, Kendall’s τ, and copula parameter θ.
Family Gaussian Frank Clayton Gumbel Plackett Kendall’s coefficientτ Copula parameterθ 0.5922 4.2092 1.3528 1.6764 5.6180 0.4035 AICc - 41.0086 - 42.8635 - 36.9945 - 41.3524 - 41.7749 S n/pc 0.0437 0.0328 0.0502 0.2542 0.0367 RMSE 0.2203 0.2212 0.2298 0.2293 0.2282 Bias 0.1923 0.1867 0.2178 0.1940 0.1925Using 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 g-line, as shown in Table 4. In the analyses of MCS, 10000 random numbers of c 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.Table 4
Probabilities of failure and their relative error of the example slope.
Method Copula-based direct integration method (CDIM) FORM MCS10000 random numbers Gaussian Frank Clayton Gumbel Plackett P f 0.2606 0.2531 0.2668 0.2529 0.2571 0.2478 0.2674 Relative error (%) 5.17 2.14 7.67 2.06 3.75 0 / 2.54 5.35 0.22 5.42 3.85 / 0


## 6. Discussion

As previously discussed, the optimum copula is not unique for different evaluation methods. Therefore, AICc,Sn/pc, 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 (c and tan⁡φ) can be expressed in terms of the weighted combination and density function as follows:(24)fEWc,φ=∑i=15λiciF1c,F2φ;θi,where λi, ci, and θi are the ith selected weight, density function, and parameter of the copula function, respectively, i=1,…,5. Similarly, the term fEW in (24) is integrated for the g-line failure domain. Consequently, the probability of failure of a slope system is determined for the considered EW copula.Table 5
P f and R values of the example slope calculated for different Fs values.
Number F s Gaussian Frank Clayton Gumbel Plackett EW copula P f/weight/R P f/weight/R P f/weight/R P f/weight/R P f/weight/R P f/R 1 1.083 0.2606/0.201/1.00 0.2531/0.211/0.97 0.2668/0.190/1.03 0.2529/0.195/0.97 0.2571/0.203/0.99 0.2599/1 2 1.218 0.0922/0.232/0.99 0.0927/0.169/1.01 0.0855/0.165/0.92 0.0956/0.252/1.03 0.0973/0.183/1.05 0.0930/1 3 1.357 0.0247/0.142/0.92 0.0269/0.111/1.10 0.0161/0.196/0.60 0.0300/0.364/1.12 0.0325/0.187/1.21 0.0269/1 4 1.486 0.0020/0.220/0.78 0.0095/0.127/1.28 0.0058/0.196/0.27 0.0090/0.228/1.22 0.0111/0.229/1.50 0.0073/1 5 1.615 0.0010/0.208/0.48 0.0028/0.207/1.33 0.0004/0.187/0.19 0.0024/0.190/1.14 0.0036/0.208/1.76 0.0021/1In addition, according to the example slope model, five-group factors of safety with differentc 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 (R) 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.Figure 7
P f and Fs of the illustrative slope calculated using different copulas.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 (whenFs=1.615); in addition, the value of R 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, the EW copula also can avoid over- or underestimating the slope reliability.


## 7. 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 theg-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, theg-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 g-line, and the coefficient of determination R2 ranges from 0.92 to 0.99. Higher-order polynomials exhibit a better fitting effect for high and steep slopes. The slope failure domain determined using the g-line can be represented by shear strength parameters (c 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.



---

*Source: 6141838-2016-09-29.xml*