# Application of Probabilistic Method to Stability Analysis of Gravity Dam Foundation over Multiple Sliding Planes

Authors: Gang Wang; Zhenyue Ma
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The current challenge to the engineering profession is to carry out probabilistic methods in practice. The design point method in generalized random space (DPG method) associated with the method of divided difference can be utilized to deal with the complex problem of probability calculation of implicit performance function with nonnormal and correlated variables. For a practical concrete gravity dam, the suggested method is performed to calculate the instability probability of the dam foundation over multiple sliding places. The general conclusions drawn in the paper are identical to those in other research and the method is proved to be feasible, accurate, and efficient. As the same analysis principle, the method can also be used in other similar fields, such as in fields of slopes, earth-rock dams, levees, and embankments.


---

## Body

## 1. Introduction

It is a trend to use the probabilistic method for evaluation of the risk of failure in almost all engineering fields, including some in geotechnical or structural engineering [1–3]. The current challenge to the geotechnical and structural engineering profession is to carry out probabilistic methods in practice [4]. Many researchers have focused on the research topics of reliability-based design and risk analysis and made progress in resolving the problem about the instability failure of slopes, dams, levees, embankments, and other geotechnical or structural engineering fields in recent years [1, 4, 5].The Monte Carlo method (MC) is usually used to estimate the reliability index,β. However this method is rarely adopted due to its huge calculation time [5, 6]. Besides MC, many other methods have been proposed for reliability analysis, such as the first-order reliability method (FORM) [1, 7], second-order reliability method (SORM) [8, 9], and some improved methods. In order to obtain β, the partial derivatives of performance function, g
(
X
), are needed in these methods. But, in geotechnical or structural engineering, g
(
X
) is usually implicit, and its partial derivatives are complex or difficult to be derived from implicit to explicit. Therefore, these conventional reliability methods only can be used to analyze small structures.Always, the response surface methods (RSMs) are utilized to obtain the solution ofβ for reliability problem with implicit g
(
X
) for complex structures. Wong applied RSM to evaluate the reliability of a homogeneous slope [10]. Moore and Sa constructed confidence intervals about the difference in mean responses at the stationary point and alternate points based on the proposed delta method and F-projection method and compared coverage probabilities and interval widths [11]. Zheng and Das proposed an improved response surface method and applied that to the reliability analysis of a stiffened plate structure [12]. Guan and Melchers evaluated the effect of response surface parameter variation on structural reliability [13]. Gupta and Manohar used the response surface method to study the extremes of Von Mises stress in nonlinear structures under Gaussian excitations [14]. Wong et al. proposed an adaptive design approach to overcome the problem, which was that the solution of the reliability analysis initially diverged when the loading was applied in sequence in the nonlinear finite element (NLFE) analysis, and made several suggestions to improve the robustness of RSM [15]. Xu and Low used RSM to approximate the performance function of slope stability in slope reliability analysis, in which the response surface is taken as a bridge between standalone numerical packages and spreadsheet-based reliability analysis [16]. Cheng et al. presented a new artificial neural network (ANN) based response surface method in conjunction with the uniform design method for predicting failure probability of structures [17]. Gavin and Yau described the use of higher order polynomials in order to approximate the true limit state more accurately in contrast to recently proposed algorithms which focused on the positions of sample points to improve the accuracy of the quadratic stochastic response surface method (SRSM) [18]. Zou et al. presented an accurate and efficient MC for limit state-based reliability analysis at both component and system levels, using a response surface approximation of the failure indicator function [19]. Nguyen et al. proposed an adaptive construction of the numerical design, in which the response surface was fitted by the weighted regression technique, which allowed the fitting points to be weighted according to their distance from the true failure surface and their distance from the estimated design point [20]. Similar to support vector machine- (SVM-) based RSM, Samui et al. adopted relevance vector machine- (RVM-) based first-order second-moment method (FOSM) to build a RVM model to predict the implicit performance function and evaluate the partial derivatives with sufficient accuracy [21]. Tan et al. discussed similarities and differences between radial basis function networks (RBFN) based RSMs and SVM-based RSMs, which indicated that there is no significant difference between them, and then proposed two new sampling methods and a hybrid RSM to reduce the number of evaluations of the actual performance function [22, 23].However, RSM and its improved methods mentioned above are relatively complicated and need much more computing cost, since they have a lot of iterating calculations at different design points associated with numerical method (i.e., finite element method) to fit limit state curved surface, which is represented byg
(
X
)
=
0. In this paper, a novel design point approach in generalized random space (DPG method) associated with the method of divided difference is proposed and applied to analyze the probability of gravity dam foundation instability over multiple sliding planes. In the method, implicit performance function with nonnormal and correlated variables is considered, and iterating calculation for β is performed in generalized random space directly. The whole procedure for reliability analysis does not need much preparation mathematically and is relatively simple for the complex engineering problems.


## 2. Basic Principle of Computing Instability Probability of Gravity Dam Foundation

Due to the complexity of foundation stability against sliding of gravity dam, for simplification, the instability probability is commonly evaluated using models of single or dual sliding plane(s) for the dam foundation [24]. However, instability of foundation over multiple sliding planes is a general case occurring in bedrock under most engineering geological conditions for a gravity dam. Therefore, it is necessary to apply a model of multiple sliding planes to analysis on instability probability of gravity dam foundation, but there is limited existing research regarding this case. In traditional deterministic analysis, equal safety coefficient method (namely, equal-K method) is used to analyze and assess the stability failure of dam foundation. In the method, factor of safety, K, is regarded as an unknown number in nonlinear equations and needs to be solved by iterative calculating numerically, similar to the application of the methods of Bishop, Spencer, and Janbu as well, which are widely used for stability analysis on slope. The reliability index β or failure probability p
f is calculated by implicit performance function, g
(
X
)
=
K
-
1, associated with the factor of safety K. Accounting for this, the conventional FORM and its improved methods (such as JC method which combines the method of equivalent normalization transformation with FORM and is recommended by the Joint Committee of Structural Safety) cannot be directly applied since the partial derivatives of g
(
X
) with respect to the variable X cannot be evaluated directly [1, 25]. Hence, if a partial derivative of g
(
X
) is calculated by divided difference mathematically, the value of the derivative at any design point x
∗ will be worked out easily, and the limit state curved surface will not be fitted completely as RSM [1, 25]. Usually, the JC method associated with the method of divided difference could calculate design points by iteration on the actual limit state curved surface so that it could avoid fitting the overall curved surface. The method can provide sufficient accuracy and efficiency to calculate probability of failure as prone to dealing with complex problems for large projects, which has been verified by many examples in some literatures. So it is a good idea to use the similar method to calculate instability probability of gravity dam foundation over multiple sliding planes.


### 2.1. Uncertainty in the Analysis

Probabilistic methods can reveal the contributions of different components to the uncertainty in the analysis on the probability of failure of dam foundation for gravity dam. The uncertainty related to instability of dam foundation could be classified as follows: (1) the uncertainty of hydrographic and hydraulic parameters; (2) the uncertainty of mechanical parameters of rocks and soils; (3) the uncertainty of seismic excitation; (4) other uncertainty, underlying the process of dam design, construction, operation, management, and so on.The second type of uncertainty above is studied mainly in the paper, and others are ignored by deterministic engineering applications just like in conventional analysis for dam stability against sliding. The uncertainty of mechanical parameters of rocks and soils might be enslaved to analysis methods, due to subjective or objective conditions, derived from the lack of knowledge, and so on. Always, the laboratory andin situ measured data of the properties of rocks or soils are inadequate. So there is demonstrable distinction between the measured values from the statistics of the small sampling and the true values of the actual physical and mechanical properties existing in real world.Many of studies have shown that friction coefficientf
′ and cohesion force c
′ against shear fracture of rocks or discontinuities are two important uncertainty parameters influencing the foundation stability of gravity dam. Firstly, by a great number of data of field measurements and their statistical back-analyses, there are obvious scatter in the spatial variation of the two shear strength parameters not only for rock mass and faults of different projects but also for different types of rock mass and faults at different locations in the same project. Chen et al. indicated that the variability of f
′ is nearly at the same level, but the variability of c
′ is not so by the study on the variability of the shear strength parameters of the foundation rock mass and faults at several dam sites in China [26, 27]. Secondly, the probability distributions and negative correlation of the two parameters cannot be ignored. The random parameters might be normally distributed or ln-normally distributed.In our probabilistic analysis, the two parameters against shear fracture,f
′ and c
′, are to be regarded as random variables following some distributions, such as normal or ln-normal distribution, and considering their variability and correlation.


### 2.2. Performance Function for Computing Instability Probability

#### 2.2.1. Equivalence Safety Factor Method for Foundation Instability over Multiple Sliding Planes

The sliding patterns of gravity dam foundation could be divided by slip paths over single, dual, and multiple plane(s). The last one is the general case for analysis of foundation instability of gravity dam by assuming the potential failure surface as shown in Figure1. There are n slip planes of n sliding wedges from upstream to downstream in the foundation. The ith wedge is subjected to (a) the vertical normal stress σ and the horizontal shear stress τ induced by action of dam or other upside loads along the topside of the wedge, (b) dead load R
i, (c) force N
i and uplift U
i normal to the ith slip plane, (d) resistance Q
i provided by the (
i
+
1
)th wedge, (e) angle of inclination α
i of the interface between the ith and (
i
+
1
)th wedges, and (f) angle of inclination β
i, the length l
i, and the shear strength parameters f
i
′ and c
i
′ of the ith slip plane.Figure 1
Illustration of multiple sliding planes for stability calculation of dam foundation.Based on the limit equilibrium method, the normal forceN
i and applied shear T
i of the ith slip plane, respectively, are expressed as(1)
N
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
Q
i
cos
⁡
α
i
+
β
i
-
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
+
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
(2)
T
i
=
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i
.The stability factor of safety for theith wedge is defined as(3)
K
i
=
N
i
f
i
′
+
c
i
′
l
i
T
i
.So substituting the expressions forN
i and T
i into (3) for the typical wedge and according to the definition of equivalence safety factor method, K
1
=
K
2
=
⋯
=
K
n
=
K, the general wedge and wedge interaction equation can be written as(4)
K
=
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
f
i
′
-
Q
i
cos
⁡
α
i
+
β
i
f
i
′
+
p
i
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
q
i
,where p
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
f
i
′
+
c
i
′
l
i, and q
i
=
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i.The nonlinear equations (4) are in terms of K and Q
i (i
=
0,1
,
…
,
n), and Q
0
=
Q
n
=
0. So there are eventually n unknown numbers for n sliding wedges of dam foundation. Assuming α
i
=
90° according to conventional vertical slice method [26, 28–30], K and Q
i (i
=
1
,
…
,
n
-
1) can be worked out by iterative calculation numerically.


#### 2.2.2. Performance Function for Probability Calculation

The basic expression of the performance function for reliability analysis on dam foundation stability could be given as(5)
Z
=
K
-
1
,where K is the factor of safety mentioned in (4), which is the function in terms of variables X
1
,
X
2
,
…
,
X
l. So (5) can be rewritten as(6)
Z
=
K
-
1
=
g
X
1
,
X
2
,
…
,
X
l
,where g
X
1
,
X
2
,
…
,
X
l is the function associated with the random variables X
1
,
X
2
,
…
,
X
l.In the analysis of instability of foundation over multiple sliding planes,Z
=
g
(
X
1
,
X
2
,
…
,
X
l
) is implicit function, and limit state equation, Z
=
0, can be represented by nonlinear equation (4) above. Correspondingly, the random variables X
1
,
X
2
,
…
,
X
l (l
=
2
n) are actually the strength parameters against shear fracture, f
i
′ and c
i
′, where i
=
1,2
,
…
,
n.


### 2.3. Computing Method of Failure Probability with Nonnormal and Correlative Random Variables

It is a problem to calculate instability probability with nonnormal and correlated random variables of dam rock base over multiple sliding planes. What is more, the random characteristics of the variables will influence the accuracy of the value of instability probability.A well known method for calculating reliability index or probability of failure of a structure, MC simulation, is derived from sampling principle. In the method, each continuous variable is replaced by a large number of discrete values generated from the underlying distribution. These values are used to compute a large number of values of functionZ and its distribution. Finally, the probability of failure can be estimated to be equal to the number of failures which failed divided by the total number of calculations statistically. MC has been widely applied for reliability calculation due to its clear definition in spite of the low calculating efficiency, and under some circumstances it is always the unique method to verify the accuracy of any other probabilistic method.Besides MC method, JC method is also a conventional method for calculating reliability index or probability of failure of a structure. When the random variables are nonnormally distributed, the treatment of equivalent normalization transformation [31] for the variables is performed rationally; however, when the random variables are correlative for each other, the problem becomes complicated because the calculation with nonnormal and correlative variables involves complex combining processes of equivalent normalization transformation and orthogonal transformation. Calculation of the eigenvalues of matrices must be conducted in the procedure in particular. Other approaches, such as the method of Rosenblatt or Nataf transformation, are also utilized to deal with the calculation with nonnormal and correlated variables [32, 33]. However, in the former method the joint cumulative distribution function of variables must be known and in the latter the ratio of coefficients of correlation between pre- and posttransformation of variables needs to be determined previously. The preconditions in the two methods are hardly ensured usually.So a new reliability-based method, which is called design point method in generalized random space (DPG method), can be applied, since iterative calculation forβ is performed in generalized random space directly by assuming that the value of coefficient of correlation of any two variables is not changed before and after variable transformation in generalized space based on the principle of isoprobabilistic transformation [34–37]. Obviously, the DPG method can be used in rather more general cases than JC method. When the variables are independent of each other, the former becomes the latter.A set of random variablesX
1
,
X
2
,
…
,
X
n, describing the randomness in the geometry, material, loading, and so on, can form an n-dimensional generalized random space, in which there is a performance function Z
=
g
(
X
1
,
X
2
,
…
,
X
n
) for the structure. Similar to the Euclidean random space in conventional reliability theory, the limit state of the structure is represented by Z
=
0. Z
>
0; then the structure is on the reliable state (“safe”). Z
<
0; then the structure is on the failure state (“fail”). Figure 2 shows the principle of the DPG method.Figure 2
A conceptual illustration of the design point and reliability index in 2-dimensional generalized random space.To perform DPG method, firstly the angle of any two axes is defined as(7)
θ
X
i
X
j
=
π
-
arccos
⁡
ρ
X
i
X
j
,where ρ
X
i
X
j is the coefficient of correlation for variables X
i and X
j. The coefficient of correlation for the two variables after transformation in generalized random space based on the principle of isoprobabilistic transformation is(8)
ρ
X
i
′
X
j
′
≈
ρ
X
i
X
j
,where X
i
′ and X
j
′ are arbitrary two variables by transformation in generalized space from X
i and X
j, respectively.The mean and standard deviation of any variableX
i
′ are defined as(9)
μ
X
i
′
=
x
i
∗
-
σ
X
i
′
Φ
-
1
F
i
x
i
∗
,
σ
X
i
′
=
φ
Φ
-
1
F
i
x
i
∗
f
i
x
i
∗
,where x
i
∗ is any component of the vector x
∗
=
(
x
1
∗
,
x
2
∗
,
…
,
x
i
∗
,
…
,
x
n
∗
), which is the value of the vector X
=
(
X
1
,
X
2
,
…
,
X
n) at design point p
∗ shown in Figure 2; f
i
(
·
) denotes the probability density function (PDF) of the random variables; and F
i
(
·
) denotes the cumulative distribution function (CDF) of the random variables.The reliability index atx
∗ is(10)
β
=
∑
i
=
1
n
∂
g
/
∂
X
i
p
∗
μ
X
i
′
-
x
i
∗
∑
i
=
1
n
∑
j
=
1
n
ρ
X
i
′
X
j
′
∂
g
/
∂
X
i
∂
g
/
∂
X
j
p
∗
σ
X
i
′
σ
X
j
′
1
/
2and the sensitivity coefficient of any variable X
i is calculated as(11)
α
i
′
=
∑
j
=
1
n
ρ
X
i
′
X
j
′
∂
g
/
∂
X
j
p
∗
σ
X
j
′
∑
j
=
1
n
∑
k
=
1
n
ρ
X
j
′
X
k
′
∂
g
/
∂
X
j
∂
g
/
∂
X
k
p
∗
σ
X
j
′
σ
X
k
′
.x
i
∗ for next iterative calculation is(12)
x
i
∗
=
μ
X
i
′
+
β
α
i
′
σ
X
i
′
.The initialx
∗ are always the mean of random variable. When x
∗ converges by repeated calculation from (9) to (12), the resultant reliability index β and probability of failure p
f
=
1
-
Φ
(
β
) are obtained finally.


### 2.4. Method of Divided Difference to Computing Partial Derivative of Implicit Performance Function

According to (10) and (11), the value of partial derivative of performance function at design point p
∗, ∂
g
/
∂
X
i
p
∗, needs to be calculated. When the performance function, g
(
X
1
,
X
2
,
…
,
X
n
), is implicit, the method of one order divided difference can be utilized as(13)
∂
Z
∂
X
i
p
∗
=
∂
g
∂
X
i
p
∗
≈
Δ
g
Δ
x
i
=
g
x
1
,
x
2
,
…
,
x
i
+
Δ
x
i
,
…
,
x
n
Δ
x
i
-
g
x
1
,
x
2
,
…
,
x
i
,
…
,
x
n
Δ
x
i
.From the practical experience of calculation, the computing efficiency and accuracy can be ensured whileΔ
x
i
=
(
0.01

~

0.05
)
x
i.


### 2.5. Computing Steps and Procedures

DPG method associated with method of divided difference to calculate the stability reliability indexβ or instability probability p
f of gravity dam foundation over multiple sliding planes has been proposed above and can be implemented easily within Matlab. There are six steps of the solution process:(1)
Regard all of the strength parameters against shear fracture,f
i
′ and c
i
′ (i
=
1,2
,
…
,
n), of the n slip planes under n sliding wedges as the random variables X that affect the foundation stability of gravity dam; generally assume the means of f
i
′ and c
i
′ as the initial value x
∗ at current design point p
∗.
(2)
Use the currentx
∗ to calculate the means and standard deviations of the variables by (9).
(3)
Calculateβ at the current design point p
∗ by (10), and then perform sensitivity analyses according to (11). The critical calculation of partial derivative of performance function, Z
=
g
(
X
)
=
K
(
X
)
-
1 (see (6)), is performed by divided difference calculation according to (13), in which the factor of safety K solved by nonlinear equations (4) at the limit state (Z
=
0) is considered.
(4)
Use (12) to obtain the new x
∗ for next iterative calculation.
(5)
Perform repeated calculation from Step(
2
) to Step (
4
) until the values of vector x
∗ are converged.
(6)
Obtain the ultimateβ and calculate the instability probability p
f
=
1
-
Φ
(
β
) finally.


## 2.1. Uncertainty in the Analysis

Probabilistic methods can reveal the contributions of different components to the uncertainty in the analysis on the probability of failure of dam foundation for gravity dam. The uncertainty related to instability of dam foundation could be classified as follows: (1) the uncertainty of hydrographic and hydraulic parameters; (2) the uncertainty of mechanical parameters of rocks and soils; (3) the uncertainty of seismic excitation; (4) other uncertainty, underlying the process of dam design, construction, operation, management, and so on.The second type of uncertainty above is studied mainly in the paper, and others are ignored by deterministic engineering applications just like in conventional analysis for dam stability against sliding. The uncertainty of mechanical parameters of rocks and soils might be enslaved to analysis methods, due to subjective or objective conditions, derived from the lack of knowledge, and so on. Always, the laboratory andin situ measured data of the properties of rocks or soils are inadequate. So there is demonstrable distinction between the measured values from the statistics of the small sampling and the true values of the actual physical and mechanical properties existing in real world.Many of studies have shown that friction coefficientf
′ and cohesion force c
′ against shear fracture of rocks or discontinuities are two important uncertainty parameters influencing the foundation stability of gravity dam. Firstly, by a great number of data of field measurements and their statistical back-analyses, there are obvious scatter in the spatial variation of the two shear strength parameters not only for rock mass and faults of different projects but also for different types of rock mass and faults at different locations in the same project. Chen et al. indicated that the variability of f
′ is nearly at the same level, but the variability of c
′ is not so by the study on the variability of the shear strength parameters of the foundation rock mass and faults at several dam sites in China [26, 27]. Secondly, the probability distributions and negative correlation of the two parameters cannot be ignored. The random parameters might be normally distributed or ln-normally distributed.In our probabilistic analysis, the two parameters against shear fracture,f
′ and c
′, are to be regarded as random variables following some distributions, such as normal or ln-normal distribution, and considering their variability and correlation.


## 2.2. Performance Function for Computing Instability Probability

### 2.2.1. Equivalence Safety Factor Method for Foundation Instability over Multiple Sliding Planes

The sliding patterns of gravity dam foundation could be divided by slip paths over single, dual, and multiple plane(s). The last one is the general case for analysis of foundation instability of gravity dam by assuming the potential failure surface as shown in Figure1. There are n slip planes of n sliding wedges from upstream to downstream in the foundation. The ith wedge is subjected to (a) the vertical normal stress σ and the horizontal shear stress τ induced by action of dam or other upside loads along the topside of the wedge, (b) dead load R
i, (c) force N
i and uplift U
i normal to the ith slip plane, (d) resistance Q
i provided by the (
i
+
1
)th wedge, (e) angle of inclination α
i of the interface between the ith and (
i
+
1
)th wedges, and (f) angle of inclination β
i, the length l
i, and the shear strength parameters f
i
′ and c
i
′ of the ith slip plane.Figure 1
Illustration of multiple sliding planes for stability calculation of dam foundation.Based on the limit equilibrium method, the normal forceN
i and applied shear T
i of the ith slip plane, respectively, are expressed as(1)
N
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
Q
i
cos
⁡
α
i
+
β
i
-
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
+
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
(2)
T
i
=
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i
.The stability factor of safety for theith wedge is defined as(3)
K
i
=
N
i
f
i
′
+
c
i
′
l
i
T
i
.So substituting the expressions forN
i and T
i into (3) for the typical wedge and according to the definition of equivalence safety factor method, K
1
=
K
2
=
⋯
=
K
n
=
K, the general wedge and wedge interaction equation can be written as(4)
K
=
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
f
i
′
-
Q
i
cos
⁡
α
i
+
β
i
f
i
′
+
p
i
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
q
i
,where p
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
f
i
′
+
c
i
′
l
i, and q
i
=
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i.The nonlinear equations (4) are in terms of K and Q
i (i
=
0,1
,
…
,
n), and Q
0
=
Q
n
=
0. So there are eventually n unknown numbers for n sliding wedges of dam foundation. Assuming α
i
=
90° according to conventional vertical slice method [26, 28–30], K and Q
i (i
=
1
,
…
,
n
-
1) can be worked out by iterative calculation numerically.


### 2.2.2. Performance Function for Probability Calculation

The basic expression of the performance function for reliability analysis on dam foundation stability could be given as(5)
Z
=
K
-
1
,where K is the factor of safety mentioned in (4), which is the function in terms of variables X
1
,
X
2
,
…
,
X
l. So (5) can be rewritten as(6)
Z
=
K
-
1
=
g
X
1
,
X
2
,
…
,
X
l
,where g
X
1
,
X
2
,
…
,
X
l is the function associated with the random variables X
1
,
X
2
,
…
,
X
l.In the analysis of instability of foundation over multiple sliding planes,Z
=
g
(
X
1
,
X
2
,
…
,
X
l
) is implicit function, and limit state equation, Z
=
0, can be represented by nonlinear equation (4) above. Correspondingly, the random variables X
1
,
X
2
,
…
,
X
l (l
=
2
n) are actually the strength parameters against shear fracture, f
i
′ and c
i
′, where i
=
1,2
,
…
,
n.


## 2.2.1. Equivalence Safety Factor Method for Foundation Instability over Multiple Sliding Planes

The sliding patterns of gravity dam foundation could be divided by slip paths over single, dual, and multiple plane(s). The last one is the general case for analysis of foundation instability of gravity dam by assuming the potential failure surface as shown in Figure1. There are n slip planes of n sliding wedges from upstream to downstream in the foundation. The ith wedge is subjected to (a) the vertical normal stress σ and the horizontal shear stress τ induced by action of dam or other upside loads along the topside of the wedge, (b) dead load R
i, (c) force N
i and uplift U
i normal to the ith slip plane, (d) resistance Q
i provided by the (
i
+
1
)th wedge, (e) angle of inclination α
i of the interface between the ith and (
i
+
1
)th wedges, and (f) angle of inclination β
i, the length l
i, and the shear strength parameters f
i
′ and c
i
′ of the ith slip plane.Figure 1
Illustration of multiple sliding planes for stability calculation of dam foundation.Based on the limit equilibrium method, the normal forceN
i and applied shear T
i of the ith slip plane, respectively, are expressed as(1)
N
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
Q
i
cos
⁡
α
i
+
β
i
-
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
+
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
(2)
T
i
=
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i
.The stability factor of safety for theith wedge is defined as(3)
K
i
=
N
i
f
i
′
+
c
i
′
l
i
T
i
.So substituting the expressions forN
i and T
i into (3) for the typical wedge and according to the definition of equivalence safety factor method, K
1
=
K
2
=
⋯
=
K
n
=
K, the general wedge and wedge interaction equation can be written as(4)
K
=
Q
i
-
1
cos
⁡
α
i
-
1
+
β
i
f
i
′
-
Q
i
cos
⁡
α
i
+
β
i
f
i
′
+
p
i
Q
i
-
1
sin
⁡
α
i
-
1
+
β
i
-
Q
i
sin
⁡
α
i
+
β
i
+
q
i
,where p
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
f
i
′
+
c
i
′
l
i, and q
i
=
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i.The nonlinear equations (4) are in terms of K and Q
i (i
=
0,1
,
…
,
n), and Q
0
=
Q
n
=
0. So there are eventually n unknown numbers for n sliding wedges of dam foundation. Assuming α
i
=
90° according to conventional vertical slice method [26, 28–30], K and Q
i (i
=
1
,
…
,
n
-
1) can be worked out by iterative calculation numerically.


## 2.2.2. Performance Function for Probability Calculation

The basic expression of the performance function for reliability analysis on dam foundation stability could be given as(5)
Z
=
K
-
1
,where K is the factor of safety mentioned in (4), which is the function in terms of variables X
1
,
X
2
,
…
,
X
l. So (5) can be rewritten as(6)
Z
=
K
-
1
=
g
X
1
,
X
2
,
…
,
X
l
,where g
X
1
,
X
2
,
…
,
X
l is the function associated with the random variables X
1
,
X
2
,
…
,
X
l.In the analysis of instability of foundation over multiple sliding planes,Z
=
g
(
X
1
,
X
2
,
…
,
X
l
) is implicit function, and limit state equation, Z
=
0, can be represented by nonlinear equation (4) above. Correspondingly, the random variables X
1
,
X
2
,
…
,
X
l (l
=
2
n) are actually the strength parameters against shear fracture, f
i
′ and c
i
′, where i
=
1,2
,
…
,
n.


## 2.3. Computing Method of Failure Probability with Nonnormal and Correlative Random Variables

It is a problem to calculate instability probability with nonnormal and correlated random variables of dam rock base over multiple sliding planes. What is more, the random characteristics of the variables will influence the accuracy of the value of instability probability.A well known method for calculating reliability index or probability of failure of a structure, MC simulation, is derived from sampling principle. In the method, each continuous variable is replaced by a large number of discrete values generated from the underlying distribution. These values are used to compute a large number of values of functionZ and its distribution. Finally, the probability of failure can be estimated to be equal to the number of failures which failed divided by the total number of calculations statistically. MC has been widely applied for reliability calculation due to its clear definition in spite of the low calculating efficiency, and under some circumstances it is always the unique method to verify the accuracy of any other probabilistic method.Besides MC method, JC method is also a conventional method for calculating reliability index or probability of failure of a structure. When the random variables are nonnormally distributed, the treatment of equivalent normalization transformation [31] for the variables is performed rationally; however, when the random variables are correlative for each other, the problem becomes complicated because the calculation with nonnormal and correlative variables involves complex combining processes of equivalent normalization transformation and orthogonal transformation. Calculation of the eigenvalues of matrices must be conducted in the procedure in particular. Other approaches, such as the method of Rosenblatt or Nataf transformation, are also utilized to deal with the calculation with nonnormal and correlated variables [32, 33]. However, in the former method the joint cumulative distribution function of variables must be known and in the latter the ratio of coefficients of correlation between pre- and posttransformation of variables needs to be determined previously. The preconditions in the two methods are hardly ensured usually.So a new reliability-based method, which is called design point method in generalized random space (DPG method), can be applied, since iterative calculation forβ is performed in generalized random space directly by assuming that the value of coefficient of correlation of any two variables is not changed before and after variable transformation in generalized space based on the principle of isoprobabilistic transformation [34–37]. Obviously, the DPG method can be used in rather more general cases than JC method. When the variables are independent of each other, the former becomes the latter.A set of random variablesX
1
,
X
2
,
…
,
X
n, describing the randomness in the geometry, material, loading, and so on, can form an n-dimensional generalized random space, in which there is a performance function Z
=
g
(
X
1
,
X
2
,
…
,
X
n
) for the structure. Similar to the Euclidean random space in conventional reliability theory, the limit state of the structure is represented by Z
=
0. Z
>
0; then the structure is on the reliable state (“safe”). Z
<
0; then the structure is on the failure state (“fail”). Figure 2 shows the principle of the DPG method.Figure 2
A conceptual illustration of the design point and reliability index in 2-dimensional generalized random space.To perform DPG method, firstly the angle of any two axes is defined as(7)
θ
X
i
X
j
=
π
-
arccos
⁡
ρ
X
i
X
j
,where ρ
X
i
X
j is the coefficient of correlation for variables X
i and X
j. The coefficient of correlation for the two variables after transformation in generalized random space based on the principle of isoprobabilistic transformation is(8)
ρ
X
i
′
X
j
′
≈
ρ
X
i
X
j
,where X
i
′ and X
j
′ are arbitrary two variables by transformation in generalized space from X
i and X
j, respectively.The mean and standard deviation of any variableX
i
′ are defined as(9)
μ
X
i
′
=
x
i
∗
-
σ
X
i
′
Φ
-
1
F
i
x
i
∗
,
σ
X
i
′
=
φ
Φ
-
1
F
i
x
i
∗
f
i
x
i
∗
,where x
i
∗ is any component of the vector x
∗
=
(
x
1
∗
,
x
2
∗
,
…
,
x
i
∗
,
…
,
x
n
∗
), which is the value of the vector X
=
(
X
1
,
X
2
,
…
,
X
n) at design point p
∗ shown in Figure 2; f
i
(
·
) denotes the probability density function (PDF) of the random variables; and F
i
(
·
) denotes the cumulative distribution function (CDF) of the random variables.The reliability index atx
∗ is(10)
β
=
∑
i
=
1
n
∂
g
/
∂
X
i
p
∗
μ
X
i
′
-
x
i
∗
∑
i
=
1
n
∑
j
=
1
n
ρ
X
i
′
X
j
′
∂
g
/
∂
X
i
∂
g
/
∂
X
j
p
∗
σ
X
i
′
σ
X
j
′
1
/
2and the sensitivity coefficient of any variable X
i is calculated as(11)
α
i
′
=
∑
j
=
1
n
ρ
X
i
′
X
j
′
∂
g
/
∂
X
j
p
∗
σ
X
j
′
∑
j
=
1
n
∑
k
=
1
n
ρ
X
j
′
X
k
′
∂
g
/
∂
X
j
∂
g
/
∂
X
k
p
∗
σ
X
j
′
σ
X
k
′
.x
i
∗ for next iterative calculation is(12)
x
i
∗
=
μ
X
i
′
+
β
α
i
′
σ
X
i
′
.The initialx
∗ are always the mean of random variable. When x
∗ converges by repeated calculation from (9) to (12), the resultant reliability index β and probability of failure p
f
=
1
-
Φ
(
β
) are obtained finally.


## 2.4. Method of Divided Difference to Computing Partial Derivative of Implicit Performance Function

According to (10) and (11), the value of partial derivative of performance function at design point p
∗, ∂
g
/
∂
X
i
p
∗, needs to be calculated. When the performance function, g
(
X
1
,
X
2
,
…
,
X
n
), is implicit, the method of one order divided difference can be utilized as(13)
∂
Z
∂
X
i
p
∗
=
∂
g
∂
X
i
p
∗
≈
Δ
g
Δ
x
i
=
g
x
1
,
x
2
,
…
,
x
i
+
Δ
x
i
,
…
,
x
n
Δ
x
i
-
g
x
1
,
x
2
,
…
,
x
i
,
…
,
x
n
Δ
x
i
.From the practical experience of calculation, the computing efficiency and accuracy can be ensured whileΔ
x
i
=
(
0.01

~

0.05
)
x
i.


## 2.5. Computing Steps and Procedures

DPG method associated with method of divided difference to calculate the stability reliability indexβ or instability probability p
f of gravity dam foundation over multiple sliding planes has been proposed above and can be implemented easily within Matlab. There are six steps of the solution process:(1)
Regard all of the strength parameters against shear fracture,f
i
′ and c
i
′ (i
=
1,2
,
…
,
n), of the n slip planes under n sliding wedges as the random variables X that affect the foundation stability of gravity dam; generally assume the means of f
i
′ and c
i
′ as the initial value x
∗ at current design point p
∗.
(2)
Use the currentx
∗ to calculate the means and standard deviations of the variables by (9).
(3)
Calculateβ at the current design point p
∗ by (10), and then perform sensitivity analyses according to (11). The critical calculation of partial derivative of performance function, Z
=
g
(
X
)
=
K
(
X
)
-
1 (see (6)), is performed by divided difference calculation according to (13), in which the factor of safety K solved by nonlinear equations (4) at the limit state (Z
=
0) is considered.
(4)
Use (12) to obtain the new x
∗ for next iterative calculation.
(5)
Perform repeated calculation from Step(
2
) to Step (
4
) until the values of vector x
∗ are converged.
(6)
Obtain the ultimateβ and calculate the instability probability p
f
=
1
-
Φ
(
β
) finally.


## 3. Assessment Standard

Dam failure is different from general safety accident. On the one hand, the consequence is very serious; sometimes it even causes destructive disasters; on the other hand, dam failure is affected by many natural factors and cannot be controlled by human. Therefore, risk standard of dam failure is a comprehensive index involving many factors, such as politics, economy, engineering technology, society, natural environment, and cultural background. If the risk standard is selected relatively low, the safety reliability resisting dam failure will be insufficient, and once failure takes place, the consequence will be serious; on the contrary, if the standard is relatively high, the project cannot bring anticipated economic benefit and may cause water waste and investment loss. From this point, many experts and researchers have studied the establishment of appropriate design standards. Based on the related data of different projects in the worldwide, the allowable probability of failure,P
y, is 1.1
×
10
-
4
/
a for large and middle dams and the failure probability is unacceptable when it is larger than 1.0
×
10
-
3
/
a. Among all the failure modes, stability failure accounts for about 10% of the total of dam failure according to ICOLD (the International Commission on Large Dams) Bulletin 99 (1995), “Dam FailuresStatistical Analysis,” and Chinese “unified standard for reliability design of hydraulic engineering structures.” Therefore, it is considered that the acceptable probability of stability failure is less than 1.1 × 10−5/a, and the instability probability is unacceptable when it is larger than 1.0
×
10
-
4
/
a.In contrast, the allowable safety factor [K] is 3.0 under normal loading condition according to deterministic limit equilibrium method in Chinese Design Specification of Concrete Gravity Dam (DL 5108-1999), and the target reliability index β
T is 4.2 for the hydraulic structures at Grade 1 and 3.7 for the hydraulic structures at Grade 2 in the relative reliability-based design codes.


## 4. Application

### 4.1. Basic Data

A hydropower station located in southeast China is for power generation singly without other comprehensive benefits. Its main body is roller compacted concrete gravity dam which has a maximum height of 119.0 m, the crest length of 284.0 m, and the elevation of 2564.00 m. This dam consists of 14 blocks and is numbered from number 1 to number 14 from left bank to right bank. Number 5 to number 8 are overflow structures; others are water retaining structures. The regulating storage of the reservoir is 0.254 billion m3. The usual pool elevation under normal operation is 2560.0 m, and the corresponding tailwater elevation is 2471.61 m. The silt elevation is 2497.70 m. The uplift intensity factor is 0.25. The structure grade of the dam is Grade 2, so the design reference period N
d is 50 a according to the design specification. The allowable instability probability is [
P
]
=
P
y
×
N
d
=
5.5
×
10
-
4, and the probability is unacceptable when it is larger than 5
×
10
-
3.


#### 4.1.1. Engineering Geological Conditions

The geological exploration in the dam site has shown that there are a great number of lowly inclined weak structural planes and several fault-fractured zones within the foundation which are adverse to the dam stability. The geologic parameters of the bed base rock and the structural planes are listed in Tables1 and 2, respectively.Table 1
Geologic parameters of bed base rock in the dam site.
Class of rock mass
Dry densityρ/g/cm3
Saturated compressive strengthR
b/MPa
Bedrock/bedrock
Concrete/bedrock
Deformation moduleE
0/GPa
Poisson’s ratioμ
Strength against shear fracture
Strength against shear fracture
f
′
c
′/MPa
f
′
c
′/MPa
II
2.65–2.70
80–100
15~21
0.20–0.25
1.20~1.35
1.50~1.60
1.10~1.20
1.10~1.15
III1
2.60–2.65
60–80
8~12
0.25–0.30
1.10~1.20
1.00~1.30
1.00~1.07
0.90~1.00Table 2
Geologic parameters of structural planes in the dam site.
Classification
Strength against shear fracture
Category
Subcategory
f
′
c
′/MPa
Hard structural plane
0.60~0.70
0.10~0.15
Weak structural plane or fault-fractured zones
Clastic and fragment
0.45~0.50
0.10~0.15
Clastic with mud
0.35~0.40
0.05~0.10
Mud with silty grain clastic
0.25~0.30
0.01~0.05
Clayey mud
0.18~0.25
0.002~0.010Figure3 shows the geologic section of typical sliding modes of number 6 overflow section. There are lowly inclined structural planes (mainly fh09, fh10, and fh11) and joints (fh152-2, etc.) below the dam base and steeply inclined faults (F16, F17, etc.) near to the dam heel. The geological characteristics of lowly inclined weak structural planes are listed in Table 3. Table 4 lists the slip paths and the values of the material properties of the sliding planes used for stability calculation.Table 3
Geologic characteristics of lowly inclined weak structural planes.
Name of structural planes
Occurrence
Width of fractured zone for weak structural planes/cm
Extension length/m
Type of structural planes
fh09
N40–60°W/NE∠ 20–25°
2–30
210
Clastic with mud
fh10
N40–60°W/NE∠ 20–25°
2–20
200
Clastic with mud
fh11
N40–60°W/NE∠ 20–25°
2–8
125
Clastic and fragment
fh152-5
N80°E/NW∠ 20°
2~5
/
Mud with silty grain clasticTable 4
Description of the slip paths and the value of the material properties of the sliding planes.
Number of slip paths
Name of sliding planes comprising the paths
Angle of apparent dip/°
Characteristics
Strength against shear fracture
f
′
c
′/MPa
A
A1-A2
90
Tensile fractured plane at dam heal
0.000
0.000
A2-A3
−16
Clastic with mud
0.375
0.075
A3-A4
−16
Clastic and fragment
0.475
0.125
A4-A5
−16
Clastic with mud
0.375
0.075
A5-A6
0
Dam base
1.150
1.125
B
B1-B2
90
Tensile fractured plane at dam heal
0.000
0.000
B2-B3
−14
Clastic with mud
0.375
0.075
B3-B4
−16
Clastic and fragment
0.475
0.125
B4-B5
−14
Clastic with mud
0.375
0.075
C
C1-C2
90
Tensile fractured plane at dam heal
0.000
0.000
C2-C3
−15
Clastic and fragment
0.475
0.125
C3-C4
−19
Rock bridge (II)
1.275
1.550
C4-C5
−16
Mud with silty grain clastic
0.275
0.030
C5-C6
−16
Rock bridge (III1)
1.050
0.850Figure 3
Geological profile and assuming slip paths of the dam foundation.


#### 4.1.2. Analyses of Slip Paths and Wedges

The slip paths of dam foundation are consisted of the faults, weak structural planes, joints, and rock bridges by their probable combinations (shown in Figure3). For stability analysis of the dam foundation, the bed base rock will be divided in several wedges according to any different slip path. Figure 4 gives the divided wedges on Slip Path B, for example. There are one tensile fractured plane (B1-B2) below dam heal and three sliding planes (B2-B3, B3-B4, and B4-B5) at the bottom of the three wedges (Wedges I, II, and III).Figure 4
Wedges on Slip Path B.V: vertical load due to dam action. H: horizontal load due to dam action.


#### 4.1.3. Results from Deterministic Analyses

Under the normal loading condition, the safety factors are, respectively, 2.56, 2.45, and 10.71 for the three slip paths mentioned above by deterministic equivalence safety factor method without taking the variation of parameters into consideration. The first two factors of safety are less than the allowance safety factor,[
K
]
=
3.0. So it is necessary to meet the stability of dam foundation against sliding by some effective engineering treatments.In the next sections, the influence of randomness of strength parameters against shear fracture of sliding planes on foundation stability will be researched by DPG method associated with the method of divided difference proposed by the paper. That is, the foundation instability probability affected by the probability distribution, variability, and negative correlation of the two parameters,f
′ and c
′, is considered in detail.


### 4.2. Limit State Equations and Random Parameters

For the practical dam project, the limit state equations for the most adverse slip path, Slip Path B, of number 6 overflow section are expressed by (4) as follows:(14)
K
Q
0
sin
⁡
α
0
+
β
1
-
K
Q
1
sin
⁡
α
1
+
β
1
+
K
q
1
-
Q
0
cos
⁡
α
0
+
β
1
f
1
′
+
Q
1
cos
⁡
α
1
+
β
1
f
1
′
-
p
1
=
0
,
K
Q
1
s
i
n
α
1
+
β
2
-
K
Q
2
s
i
n
α
2
+
β
2
+
K
q
2
-
Q
1
c
o
s
α
1
+
β
2
f
2
′
+
Q
2
c
o
s
α
2
+
β
2
f
2
′
-
p
2
=
0
,
K
Q
2
s
i
n
α
2
+
β
3
-
K
Q
3
s
i
n
α
3
+
β
3
+
K
q
3
-
Q
2
c
o
s
α
2
+
β
i
f
3
′
+
Q
3
c
o
s
α
3
+
β
3
f
3
′
-
p
3
=
0
,where(15)
α
0
=
α
1
=
α
2
=
9
0
∘
,
β
1
=
-
1
4
∘
,
β
2
=
-
1
6
∘
,
β
3
=
-
1
4
∘
,
p
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
f
i
′
+
c
i
′
l
i
,
i
=
1,2
,
3
,
q
i
=
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i
,
i
=
1,2
,
3
,
∫
σ
d
A
i
=
V
i
,
i
=
1,2
,
3
,
∫
τ
d
A
i
=
H
i
,
i
=
1,2
,
3
,
l
1
=
49.9
m
,
l
2
=
33.5
m
,
l
3
=
117.4
m
,
Q
0
=
Q
3
=
0
.According to the practical engineering experience and the special geological and mechanical parameters of the bed base rock of the example dam in the paper, the degree of scatter of the variation coefficients of the strength parameters against shear fracture is obvious, so the coefficients of variation off
′ and c
′ are assumed to change in certain ranges. The statistical characteristics of f
′ and c
′ of Slip Path B used for calculation are listed in Table 5. According to Chinese “unified standard for reliability design of hydraulic engineering structures” (1994), the parameter distribution type of slip planes can be assumed in normal or ln-normal distributions for reliability analysis when the statistical data are not sufficient.Table 5
Statistical characteristics of random variables of slip path B.
Slip plane
Random variable
Mean value
Variation coefficient
Probability distribution
B2-B3
f
1
′
0.375
0.14~0.26
Normal/ln-normal distribution
c
1
′
/MPa
0.075
0.15~0.45
Normal/ln-normal distribution
B3-B4
f
2
′
0.475
0.18~0.33
Normal/ln-normal distribution
c
2
′/MPa
0.125
0.25~0.75
Normal/ln-normal distribution
B4-B5
f
3
′
0.375
0.14~0.26
Normal/ln-normal distribution
c
3
′/MPa
0.075
0.15~0.45
Normal/ln-normal distribution


### 4.3. Resultant Analyses

#### 4.3.1. Feasibility of DPG Method

Two methods, DPG and MC methods, are used to calculate the instability probabilities of Slip Path B during design reference period. Table6 gives the calculated results obtained by the two methods; when f
′ and c
′ are both normally distributed, their variation coefficients, δ
f
′ and δ
c
′, are 0.2 and 0.3, respectively, and their coefficient of correlation, ρ
f
′
c
′, is 0.Table 6
Comparison of resultant instability probabilities calculated by MC and DPG method.
Method
Design points (sliding plane B2-B3)
β
p
f
Notes
f
1
′
∗
c
1
′
∗/MPa
MC
/
/
3.650
1.310 × 10−4
The numbers of calculating performance function,n
=
1
,
000
,
000, and the variation coefficient for p
f, V
p
f
=
0.0873, when 95% confidence level is specified
DPG
0.107
0.039
3.664
1.244 × 10−4
The numbers of calculating performance function,n
=
36, with the step Δ
x
i
=
0.03
x
i, for divided difference calculationThe instability probabilities of Slip Path B are, respectively,1.310
×
10
-
4 and 1.244
×
10
-
4 calculated by MC and DPG methods, and the absolute and relative errors between the two methods are, respectively, 0.0007% and 5.3%. Thus, the two methods have nearly the same accuracy, but the latter requires less iterating calculation for convergence. Meanwhile, by DPG method the design points can be obtained easily, and the probability distributions, variability, and negative correlation of f
′ and c
′ can be taken into consideration for special research. So the DPG method suggested in the paper has much merit and is rational and feasible for the major research goal of the paper.


#### 4.3.2. Analyses of the Effects of Correlation off
′ and c
′ on Instability Probability

Table7 lists the calculated results of the instability probabilities of Slip Path B by DPG method; when f
′ and c
′ are both normally distributed, δ
f
′ and δ
c
′ are 0.2 and 0.3, respectively, and the coefficient of correlation, ρ
f
′
c
′, is within the range of −0.7~0.0. The negative value of ρ
f
′
c
′ indicates the degree of negative correlation of f
′ and c
′. It is clear from Table 7 that the resultant instability probability calculated decreases with the increasing of the absolute value of negatively correlated coefficient. p
f
=
1.244
×
10
-
4 when ρ
f
′
c
′
=
0.0 (independent case), and p
f = 3.106 × 10−7 when ρ
f
′
c
′
=
-
0.7 (strongly negatively correlated case). By comparison, the values of the instability probability under the two cases are different in 3 orders of magnitude. Correspondingly, the values of the reliability index under the two cases are 3.664 and 4.984, respectively. Obviously the independent case is on the conservative side for simulation.Table 7
Calculated results under the conditions of differentρ
f
′
c
′ by DPG method.
Number
ρ
f
′
c
′
Design points (sliding plane B2-B3)
β
p
f
f
1
′
∗
c
1
′
∗/MPa
1
0.0
0.107
0.039
3.664
1.244 × 10−4
2
−0.1
0.100
0.044
3.789
7.554 × 10−5
3
−0.2
0.092
0.050
3.929
4.262 × 10−5
4
−0.3
0.082
0.057
4.086
2.194 × 10−5
5
−0.4
0.071
0.065
4.264
1.005 × 10−5
6
−0.5
0.058
0.075
4.467
3.961 × 10−6
7
−0.6
0.042
0.087
4.704
1.275 × 10−6
8
−0.7
0.022
0.102
4.984
3.106 × 10−7
Note: ifρ
f
′
c
′ is less than −0.8, f
1
′
∗ will be less than 0, and the simulation will be not converged at design points.


#### 4.3.3. Analyses of the Effects in Different Distributions off
′ and c
′ on Instability Probability

Figure5 illustrates the resultant instability probability under four different cases of combinations of normal or ln-normal distributions of f
′ and c
′ for δ
f
′
=
0.2 and δ
c
′
=
0
.
3 by DPG method. In all cases, p
f tends to increase as ρ
f
′
c
′ increases from −0.7 (negatively correlated) to 0 (independent). When f
′ is normally distributed, the instability probabilities are approximately identical whether c
′ is normally or ln-normally distributed. When f
′ is ln-normally distributed, the instability probabilities are smaller and the values calculated under the case that c
′ is in ln-normal distribution are less than in normal distribution at the same level of correlation of the two parameters. From the figure, p
f = 1.24 × 10−4 if f
′ and c
′ are both normally distributed, and p
f
=
4.64
×
10
-
9 if the two parameters are both ln-normally distributed, for δ
f
′
=
0.2, δ
c
′
=
0.3, and ρ
f
′
c
′
=
0.0.Figure 5
Relation ofp
f and ρ
f
′
c
′ assuming that f
′ and c
′ follow different distributions.


#### 4.3.4. Analyses of the Effects of Variability off
′ and c
′ on Instability Probability

On the basis of the above analyses, it is reasonable, simple, and conservative to assume thatf
′ and c
′ are normally distributed and independent of each other unless otherwise demonstrated. Considering the range of the values of coefficients of variation, δ
f
′
= 0.14~0.26, δ
c
′
= 0.15~0.45, the calculations are made, and the curves about p
f varying with δ
f
′ and δ
c
′ are displayed in Figures 6 and 7. According to the resultant curves, p
f starts to apparently climb as δ
f
′ increases for different δ
c
′ (shown in Figure 6), and, in addition, p
f also tends to increase as δ
c
′ increases for different δ
f
′ (shown in Figure 7). However, p
f
=
1.92
×
10
-
3 if δ
f
′
=
0.26, and p
f
=
4.04
×
10
-
7 if δ
f
′
=
0.14 for δ
c
′
=
0.30; meanwhile, p
f
=
3.18
×
10
-
4 if δ
c
′
=
0.45, and p
f
=
6.08
×
10
-
5 if δ
c
′
=
0.15 for δ
f
′
=
0.20. There is the difference of p
f in 4 orders of magnitude between δ
f
′
=
0.14 and 0.26 for δ
c
′
=
0.30; in contrast, under the latter case, there is the deference of p
f in one order of magnitude between δ
c
′
=
0.15 and 0.45 for δ
f
′
=
0.20. So the influence of the variability of f
′ is more apparent than that of c
′ on p
f despite the fact that the varying range of coefficients of variation of f
′ is less than that of c
′.Figure 6
Relation ofp
f and δ
f
′.Figure 7
Relation ofp
f and δ
c
′.


#### 4.3.5. Stability Safety Analyses and Assessments of the Practical Dam

The instability probability of the most adverse slip path of the example dam section (number 6 overflow section) mentioned above is1.244
×
10
-
4 for δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0. It is less than the allowable value of the instability probability, [P] = 5.5 × 10−4, but they are in the same order of magnitude. Meanwhile, it is 2.861 × 10−3 for δ
f
′ = 0.26, δ
c
′
=
0.45, and ρ
f
′
c
′
=
0.0. It is greater than [P] but less than 5 × 10−3, which is defined as unacceptable probability. In contrast, the factor of safety is 2.45 by deterministic method without taking the variation of parameters into consideration. The factor of safety is less than the allowable safety factor, [
K
]
=
3.0; the reliability index is 3.66 for δ
f
′ = 0.20, δ
c
′ = 0.30, and ρ
f
′
c
′
=
0.0, which nearly equals the target reliability index, β
T
=
3.7. Obviously, the conclusions are more or less different following three different safety standards. Generally, it has been verified that the deterministic method is most conservative in the three methods by the analyses of the practical dam in the paper and the conclusions from other literatures as well. From above analyses, the margins of safety of the foundation stability are not adequate for the example dam section and the enhancing treatments should be carried out especially for the fault-fractured zones and weak structural planes.Similarly, the instability probabilities of the other two slip paths can be calculated by the DPG method. Assume thatf
′ and c
′ of each sliding plane of each slip path are normally distributed and independent, and δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0 for any slip plane. The resultant instability probabilities are listed in Table 8. The three slip paths can be seen as three main failure modes for the rock foundation of the example dam, and each failure mode is independent due to the independence of each shear strength parameter of any sliding plane. Three failure modes comprise a series failure system of dam foundation. So the overall probability of failure for the overall dam foundation system can be calculated by summing over the instability probabilities of the three slip paths. The final instability probability of the overall rock foundation is 2.053 × 10−4 for a conservative consideration.Table 8
Resultant instability probabilities of three slip paths calculated by DPG method.
Slip paths
Reliability method (DPG)
Deterministic method
β
p
f
K
A
3.772
8.094 × 10−5
2.56
B
3.664
1.244 × 10−4
2.45
C
5.612
9.993 × 10−9
10.71
Note: the deterministic method is equivalence safety factor method without considering the variation of the parameters against shear fracture.It should also be noted that for the example dam the coefficients of variation and negative correlation of the two key mechanical parameters of sliding planes were selected in the relatively wider ranges, which are summarized from the field and experimental data. In the calculation for Slip Path B, four combinations of distributions of the two parameters were detailedly taken into consideration. So the final results obtained for the example dam can reveal the performance about foundation stability of any similar gravity dam in general. We confirm that the randomness off
′ is more significant than c
′, and, what is more, it is reasonable to assume a normal distribution and independence of each parameter in the absence of further information. The assumption will probably overestimate the probability of foundation instability of gravity dam in some degree.


## 4.1. Basic Data

A hydropower station located in southeast China is for power generation singly without other comprehensive benefits. Its main body is roller compacted concrete gravity dam which has a maximum height of 119.0 m, the crest length of 284.0 m, and the elevation of 2564.00 m. This dam consists of 14 blocks and is numbered from number 1 to number 14 from left bank to right bank. Number 5 to number 8 are overflow structures; others are water retaining structures. The regulating storage of the reservoir is 0.254 billion m3. The usual pool elevation under normal operation is 2560.0 m, and the corresponding tailwater elevation is 2471.61 m. The silt elevation is 2497.70 m. The uplift intensity factor is 0.25. The structure grade of the dam is Grade 2, so the design reference period N
d is 50 a according to the design specification. The allowable instability probability is [
P
]
=
P
y
×
N
d
=
5.5
×
10
-
4, and the probability is unacceptable when it is larger than 5
×
10
-
3.


### 4.1.1. Engineering Geological Conditions

The geological exploration in the dam site has shown that there are a great number of lowly inclined weak structural planes and several fault-fractured zones within the foundation which are adverse to the dam stability. The geologic parameters of the bed base rock and the structural planes are listed in Tables1 and 2, respectively.Table 1
Geologic parameters of bed base rock in the dam site.
Class of rock mass
Dry densityρ/g/cm3
Saturated compressive strengthR
b/MPa
Bedrock/bedrock
Concrete/bedrock
Deformation moduleE
0/GPa
Poisson’s ratioμ
Strength against shear fracture
Strength against shear fracture
f
′
c
′/MPa
f
′
c
′/MPa
II
2.65–2.70
80–100
15~21
0.20–0.25
1.20~1.35
1.50~1.60
1.10~1.20
1.10~1.15
III1
2.60–2.65
60–80
8~12
0.25–0.30
1.10~1.20
1.00~1.30
1.00~1.07
0.90~1.00Table 2
Geologic parameters of structural planes in the dam site.
Classification
Strength against shear fracture
Category
Subcategory
f
′
c
′/MPa
Hard structural plane
0.60~0.70
0.10~0.15
Weak structural plane or fault-fractured zones
Clastic and fragment
0.45~0.50
0.10~0.15
Clastic with mud
0.35~0.40
0.05~0.10
Mud with silty grain clastic
0.25~0.30
0.01~0.05
Clayey mud
0.18~0.25
0.002~0.010Figure3 shows the geologic section of typical sliding modes of number 6 overflow section. There are lowly inclined structural planes (mainly fh09, fh10, and fh11) and joints (fh152-2, etc.) below the dam base and steeply inclined faults (F16, F17, etc.) near to the dam heel. The geological characteristics of lowly inclined weak structural planes are listed in Table 3. Table 4 lists the slip paths and the values of the material properties of the sliding planes used for stability calculation.Table 3
Geologic characteristics of lowly inclined weak structural planes.
Name of structural planes
Occurrence
Width of fractured zone for weak structural planes/cm
Extension length/m
Type of structural planes
fh09
N40–60°W/NE∠ 20–25°
2–30
210
Clastic with mud
fh10
N40–60°W/NE∠ 20–25°
2–20
200
Clastic with mud
fh11
N40–60°W/NE∠ 20–25°
2–8
125
Clastic and fragment
fh152-5
N80°E/NW∠ 20°
2~5
/
Mud with silty grain clasticTable 4
Description of the slip paths and the value of the material properties of the sliding planes.
Number of slip paths
Name of sliding planes comprising the paths
Angle of apparent dip/°
Characteristics
Strength against shear fracture
f
′
c
′/MPa
A
A1-A2
90
Tensile fractured plane at dam heal
0.000
0.000
A2-A3
−16
Clastic with mud
0.375
0.075
A3-A4
−16
Clastic and fragment
0.475
0.125
A4-A5
−16
Clastic with mud
0.375
0.075
A5-A6
0
Dam base
1.150
1.125
B
B1-B2
90
Tensile fractured plane at dam heal
0.000
0.000
B2-B3
−14
Clastic with mud
0.375
0.075
B3-B4
−16
Clastic and fragment
0.475
0.125
B4-B5
−14
Clastic with mud
0.375
0.075
C
C1-C2
90
Tensile fractured plane at dam heal
0.000
0.000
C2-C3
−15
Clastic and fragment
0.475
0.125
C3-C4
−19
Rock bridge (II)
1.275
1.550
C4-C5
−16
Mud with silty grain clastic
0.275
0.030
C5-C6
−16
Rock bridge (III1)
1.050
0.850Figure 3
Geological profile and assuming slip paths of the dam foundation.


### 4.1.2. Analyses of Slip Paths and Wedges

The slip paths of dam foundation are consisted of the faults, weak structural planes, joints, and rock bridges by their probable combinations (shown in Figure3). For stability analysis of the dam foundation, the bed base rock will be divided in several wedges according to any different slip path. Figure 4 gives the divided wedges on Slip Path B, for example. There are one tensile fractured plane (B1-B2) below dam heal and three sliding planes (B2-B3, B3-B4, and B4-B5) at the bottom of the three wedges (Wedges I, II, and III).Figure 4
Wedges on Slip Path B.V: vertical load due to dam action. H: horizontal load due to dam action.


### 4.1.3. Results from Deterministic Analyses

Under the normal loading condition, the safety factors are, respectively, 2.56, 2.45, and 10.71 for the three slip paths mentioned above by deterministic equivalence safety factor method without taking the variation of parameters into consideration. The first two factors of safety are less than the allowance safety factor,[
K
]
=
3.0. So it is necessary to meet the stability of dam foundation against sliding by some effective engineering treatments.In the next sections, the influence of randomness of strength parameters against shear fracture of sliding planes on foundation stability will be researched by DPG method associated with the method of divided difference proposed by the paper. That is, the foundation instability probability affected by the probability distribution, variability, and negative correlation of the two parameters,f
′ and c
′, is considered in detail.


## 4.1.1. Engineering Geological Conditions

The geological exploration in the dam site has shown that there are a great number of lowly inclined weak structural planes and several fault-fractured zones within the foundation which are adverse to the dam stability. The geologic parameters of the bed base rock and the structural planes are listed in Tables1 and 2, respectively.Table 1
Geologic parameters of bed base rock in the dam site.
Class of rock mass
Dry densityρ/g/cm3
Saturated compressive strengthR
b/MPa
Bedrock/bedrock
Concrete/bedrock
Deformation moduleE
0/GPa
Poisson’s ratioμ
Strength against shear fracture
Strength against shear fracture
f
′
c
′/MPa
f
′
c
′/MPa
II
2.65–2.70
80–100
15~21
0.20–0.25
1.20~1.35
1.50~1.60
1.10~1.20
1.10~1.15
III1
2.60–2.65
60–80
8~12
0.25–0.30
1.10~1.20
1.00~1.30
1.00~1.07
0.90~1.00Table 2
Geologic parameters of structural planes in the dam site.
Classification
Strength against shear fracture
Category
Subcategory
f
′
c
′/MPa
Hard structural plane
0.60~0.70
0.10~0.15
Weak structural plane or fault-fractured zones
Clastic and fragment
0.45~0.50
0.10~0.15
Clastic with mud
0.35~0.40
0.05~0.10
Mud with silty grain clastic
0.25~0.30
0.01~0.05
Clayey mud
0.18~0.25
0.002~0.010Figure3 shows the geologic section of typical sliding modes of number 6 overflow section. There are lowly inclined structural planes (mainly fh09, fh10, and fh11) and joints (fh152-2, etc.) below the dam base and steeply inclined faults (F16, F17, etc.) near to the dam heel. The geological characteristics of lowly inclined weak structural planes are listed in Table 3. Table 4 lists the slip paths and the values of the material properties of the sliding planes used for stability calculation.Table 3
Geologic characteristics of lowly inclined weak structural planes.
Name of structural planes
Occurrence
Width of fractured zone for weak structural planes/cm
Extension length/m
Type of structural planes
fh09
N40–60°W/NE∠ 20–25°
2–30
210
Clastic with mud
fh10
N40–60°W/NE∠ 20–25°
2–20
200
Clastic with mud
fh11
N40–60°W/NE∠ 20–25°
2–8
125
Clastic and fragment
fh152-5
N80°E/NW∠ 20°
2~5
/
Mud with silty grain clasticTable 4
Description of the slip paths and the value of the material properties of the sliding planes.
Number of slip paths
Name of sliding planes comprising the paths
Angle of apparent dip/°
Characteristics
Strength against shear fracture
f
′
c
′/MPa
A
A1-A2
90
Tensile fractured plane at dam heal
0.000
0.000
A2-A3
−16
Clastic with mud
0.375
0.075
A3-A4
−16
Clastic and fragment
0.475
0.125
A4-A5
−16
Clastic with mud
0.375
0.075
A5-A6
0
Dam base
1.150
1.125
B
B1-B2
90
Tensile fractured plane at dam heal
0.000
0.000
B2-B3
−14
Clastic with mud
0.375
0.075
B3-B4
−16
Clastic and fragment
0.475
0.125
B4-B5
−14
Clastic with mud
0.375
0.075
C
C1-C2
90
Tensile fractured plane at dam heal
0.000
0.000
C2-C3
−15
Clastic and fragment
0.475
0.125
C3-C4
−19
Rock bridge (II)
1.275
1.550
C4-C5
−16
Mud with silty grain clastic
0.275
0.030
C5-C6
−16
Rock bridge (III1)
1.050
0.850Figure 3
Geological profile and assuming slip paths of the dam foundation.


## 4.1.2. Analyses of Slip Paths and Wedges

The slip paths of dam foundation are consisted of the faults, weak structural planes, joints, and rock bridges by their probable combinations (shown in Figure3). For stability analysis of the dam foundation, the bed base rock will be divided in several wedges according to any different slip path. Figure 4 gives the divided wedges on Slip Path B, for example. There are one tensile fractured plane (B1-B2) below dam heal and three sliding planes (B2-B3, B3-B4, and B4-B5) at the bottom of the three wedges (Wedges I, II, and III).Figure 4
Wedges on Slip Path B.V: vertical load due to dam action. H: horizontal load due to dam action.


## 4.1.3. Results from Deterministic Analyses

Under the normal loading condition, the safety factors are, respectively, 2.56, 2.45, and 10.71 for the three slip paths mentioned above by deterministic equivalence safety factor method without taking the variation of parameters into consideration. The first two factors of safety are less than the allowance safety factor,[
K
]
=
3.0. So it is necessary to meet the stability of dam foundation against sliding by some effective engineering treatments.In the next sections, the influence of randomness of strength parameters against shear fracture of sliding planes on foundation stability will be researched by DPG method associated with the method of divided difference proposed by the paper. That is, the foundation instability probability affected by the probability distribution, variability, and negative correlation of the two parameters,f
′ and c
′, is considered in detail.


## 4.2. Limit State Equations and Random Parameters

For the practical dam project, the limit state equations for the most adverse slip path, Slip Path B, of number 6 overflow section are expressed by (4) as follows:(14)
K
Q
0
sin
⁡
α
0
+
β
1
-
K
Q
1
sin
⁡
α
1
+
β
1
+
K
q
1
-
Q
0
cos
⁡
α
0
+
β
1
f
1
′
+
Q
1
cos
⁡
α
1
+
β
1
f
1
′
-
p
1
=
0
,
K
Q
1
s
i
n
α
1
+
β
2
-
K
Q
2
s
i
n
α
2
+
β
2
+
K
q
2
-
Q
1
c
o
s
α
1
+
β
2
f
2
′
+
Q
2
c
o
s
α
2
+
β
2
f
2
′
-
p
2
=
0
,
K
Q
2
s
i
n
α
2
+
β
3
-
K
Q
3
s
i
n
α
3
+
β
3
+
K
q
3
-
Q
2
c
o
s
α
2
+
β
i
f
3
′
+
Q
3
c
o
s
α
3
+
β
3
f
3
′
-
p
3
=
0
,where(15)
α
0
=
α
1
=
α
2
=
9
0
∘
,
β
1
=
-
1
4
∘
,
β
2
=
-
1
6
∘
,
β
3
=
-
1
4
∘
,
p
i
=
R
i
+
∫
σ
d
A
i
cos
⁡
β
i
-
∫
τ
d
A
i
sin
⁡
β
i
-
U
i
f
i
′
+
c
i
′
l
i
,
i
=
1,2
,
3
,
q
i
=
R
i
+
∫
σ
d
A
i
sin
⁡
β
i
+
∫
τ
d
A
i
cos
⁡
β
i
,
i
=
1,2
,
3
,
∫
σ
d
A
i
=
V
i
,
i
=
1,2
,
3
,
∫
τ
d
A
i
=
H
i
,
i
=
1,2
,
3
,
l
1
=
49.9
m
,
l
2
=
33.5
m
,
l
3
=
117.4
m
,
Q
0
=
Q
3
=
0
.According to the practical engineering experience and the special geological and mechanical parameters of the bed base rock of the example dam in the paper, the degree of scatter of the variation coefficients of the strength parameters against shear fracture is obvious, so the coefficients of variation off
′ and c
′ are assumed to change in certain ranges. The statistical characteristics of f
′ and c
′ of Slip Path B used for calculation are listed in Table 5. According to Chinese “unified standard for reliability design of hydraulic engineering structures” (1994), the parameter distribution type of slip planes can be assumed in normal or ln-normal distributions for reliability analysis when the statistical data are not sufficient.Table 5
Statistical characteristics of random variables of slip path B.
Slip plane
Random variable
Mean value
Variation coefficient
Probability distribution
B2-B3
f
1
′
0.375
0.14~0.26
Normal/ln-normal distribution
c
1
′
/MPa
0.075
0.15~0.45
Normal/ln-normal distribution
B3-B4
f
2
′
0.475
0.18~0.33
Normal/ln-normal distribution
c
2
′/MPa
0.125
0.25~0.75
Normal/ln-normal distribution
B4-B5
f
3
′
0.375
0.14~0.26
Normal/ln-normal distribution
c
3
′/MPa
0.075
0.15~0.45
Normal/ln-normal distribution


## 4.3. Resultant Analyses

### 4.3.1. Feasibility of DPG Method

Two methods, DPG and MC methods, are used to calculate the instability probabilities of Slip Path B during design reference period. Table6 gives the calculated results obtained by the two methods; when f
′ and c
′ are both normally distributed, their variation coefficients, δ
f
′ and δ
c
′, are 0.2 and 0.3, respectively, and their coefficient of correlation, ρ
f
′
c
′, is 0.Table 6
Comparison of resultant instability probabilities calculated by MC and DPG method.
Method
Design points (sliding plane B2-B3)
β
p
f
Notes
f
1
′
∗
c
1
′
∗/MPa
MC
/
/
3.650
1.310 × 10−4
The numbers of calculating performance function,n
=
1
,
000
,
000, and the variation coefficient for p
f, V
p
f
=
0.0873, when 95% confidence level is specified
DPG
0.107
0.039
3.664
1.244 × 10−4
The numbers of calculating performance function,n
=
36, with the step Δ
x
i
=
0.03
x
i, for divided difference calculationThe instability probabilities of Slip Path B are, respectively,1.310
×
10
-
4 and 1.244
×
10
-
4 calculated by MC and DPG methods, and the absolute and relative errors between the two methods are, respectively, 0.0007% and 5.3%. Thus, the two methods have nearly the same accuracy, but the latter requires less iterating calculation for convergence. Meanwhile, by DPG method the design points can be obtained easily, and the probability distributions, variability, and negative correlation of f
′ and c
′ can be taken into consideration for special research. So the DPG method suggested in the paper has much merit and is rational and feasible for the major research goal of the paper.


### 4.3.2. Analyses of the Effects of Correlation off
′ and c
′ on Instability Probability

Table7 lists the calculated results of the instability probabilities of Slip Path B by DPG method; when f
′ and c
′ are both normally distributed, δ
f
′ and δ
c
′ are 0.2 and 0.3, respectively, and the coefficient of correlation, ρ
f
′
c
′, is within the range of −0.7~0.0. The negative value of ρ
f
′
c
′ indicates the degree of negative correlation of f
′ and c
′. It is clear from Table 7 that the resultant instability probability calculated decreases with the increasing of the absolute value of negatively correlated coefficient. p
f
=
1.244
×
10
-
4 when ρ
f
′
c
′
=
0.0 (independent case), and p
f = 3.106 × 10−7 when ρ
f
′
c
′
=
-
0.7 (strongly negatively correlated case). By comparison, the values of the instability probability under the two cases are different in 3 orders of magnitude. Correspondingly, the values of the reliability index under the two cases are 3.664 and 4.984, respectively. Obviously the independent case is on the conservative side for simulation.Table 7
Calculated results under the conditions of differentρ
f
′
c
′ by DPG method.
Number
ρ
f
′
c
′
Design points (sliding plane B2-B3)
β
p
f
f
1
′
∗
c
1
′
∗/MPa
1
0.0
0.107
0.039
3.664
1.244 × 10−4
2
−0.1
0.100
0.044
3.789
7.554 × 10−5
3
−0.2
0.092
0.050
3.929
4.262 × 10−5
4
−0.3
0.082
0.057
4.086
2.194 × 10−5
5
−0.4
0.071
0.065
4.264
1.005 × 10−5
6
−0.5
0.058
0.075
4.467
3.961 × 10−6
7
−0.6
0.042
0.087
4.704
1.275 × 10−6
8
−0.7
0.022
0.102
4.984
3.106 × 10−7
Note: ifρ
f
′
c
′ is less than −0.8, f
1
′
∗ will be less than 0, and the simulation will be not converged at design points.


### 4.3.3. Analyses of the Effects in Different Distributions off
′ and c
′ on Instability Probability

Figure5 illustrates the resultant instability probability under four different cases of combinations of normal or ln-normal distributions of f
′ and c
′ for δ
f
′
=
0.2 and δ
c
′
=
0
.
3 by DPG method. In all cases, p
f tends to increase as ρ
f
′
c
′ increases from −0.7 (negatively correlated) to 0 (independent). When f
′ is normally distributed, the instability probabilities are approximately identical whether c
′ is normally or ln-normally distributed. When f
′ is ln-normally distributed, the instability probabilities are smaller and the values calculated under the case that c
′ is in ln-normal distribution are less than in normal distribution at the same level of correlation of the two parameters. From the figure, p
f = 1.24 × 10−4 if f
′ and c
′ are both normally distributed, and p
f
=
4.64
×
10
-
9 if the two parameters are both ln-normally distributed, for δ
f
′
=
0.2, δ
c
′
=
0.3, and ρ
f
′
c
′
=
0.0.Figure 5
Relation ofp
f and ρ
f
′
c
′ assuming that f
′ and c
′ follow different distributions.


### 4.3.4. Analyses of the Effects of Variability off
′ and c
′ on Instability Probability

On the basis of the above analyses, it is reasonable, simple, and conservative to assume thatf
′ and c
′ are normally distributed and independent of each other unless otherwise demonstrated. Considering the range of the values of coefficients of variation, δ
f
′
= 0.14~0.26, δ
c
′
= 0.15~0.45, the calculations are made, and the curves about p
f varying with δ
f
′ and δ
c
′ are displayed in Figures 6 and 7. According to the resultant curves, p
f starts to apparently climb as δ
f
′ increases for different δ
c
′ (shown in Figure 6), and, in addition, p
f also tends to increase as δ
c
′ increases for different δ
f
′ (shown in Figure 7). However, p
f
=
1.92
×
10
-
3 if δ
f
′
=
0.26, and p
f
=
4.04
×
10
-
7 if δ
f
′
=
0.14 for δ
c
′
=
0.30; meanwhile, p
f
=
3.18
×
10
-
4 if δ
c
′
=
0.45, and p
f
=
6.08
×
10
-
5 if δ
c
′
=
0.15 for δ
f
′
=
0.20. There is the difference of p
f in 4 orders of magnitude between δ
f
′
=
0.14 and 0.26 for δ
c
′
=
0.30; in contrast, under the latter case, there is the deference of p
f in one order of magnitude between δ
c
′
=
0.15 and 0.45 for δ
f
′
=
0.20. So the influence of the variability of f
′ is more apparent than that of c
′ on p
f despite the fact that the varying range of coefficients of variation of f
′ is less than that of c
′.Figure 6
Relation ofp
f and δ
f
′.Figure 7
Relation ofp
f and δ
c
′.


### 4.3.5. Stability Safety Analyses and Assessments of the Practical Dam

The instability probability of the most adverse slip path of the example dam section (number 6 overflow section) mentioned above is1.244
×
10
-
4 for δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0. It is less than the allowable value of the instability probability, [P] = 5.5 × 10−4, but they are in the same order of magnitude. Meanwhile, it is 2.861 × 10−3 for δ
f
′ = 0.26, δ
c
′
=
0.45, and ρ
f
′
c
′
=
0.0. It is greater than [P] but less than 5 × 10−3, which is defined as unacceptable probability. In contrast, the factor of safety is 2.45 by deterministic method without taking the variation of parameters into consideration. The factor of safety is less than the allowable safety factor, [
K
]
=
3.0; the reliability index is 3.66 for δ
f
′ = 0.20, δ
c
′ = 0.30, and ρ
f
′
c
′
=
0.0, which nearly equals the target reliability index, β
T
=
3.7. Obviously, the conclusions are more or less different following three different safety standards. Generally, it has been verified that the deterministic method is most conservative in the three methods by the analyses of the practical dam in the paper and the conclusions from other literatures as well. From above analyses, the margins of safety of the foundation stability are not adequate for the example dam section and the enhancing treatments should be carried out especially for the fault-fractured zones and weak structural planes.Similarly, the instability probabilities of the other two slip paths can be calculated by the DPG method. Assume thatf
′ and c
′ of each sliding plane of each slip path are normally distributed and independent, and δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0 for any slip plane. The resultant instability probabilities are listed in Table 8. The three slip paths can be seen as three main failure modes for the rock foundation of the example dam, and each failure mode is independent due to the independence of each shear strength parameter of any sliding plane. Three failure modes comprise a series failure system of dam foundation. So the overall probability of failure for the overall dam foundation system can be calculated by summing over the instability probabilities of the three slip paths. The final instability probability of the overall rock foundation is 2.053 × 10−4 for a conservative consideration.Table 8
Resultant instability probabilities of three slip paths calculated by DPG method.
Slip paths
Reliability method (DPG)
Deterministic method
β
p
f
K
A
3.772
8.094 × 10−5
2.56
B
3.664
1.244 × 10−4
2.45
C
5.612
9.993 × 10−9
10.71
Note: the deterministic method is equivalence safety factor method without considering the variation of the parameters against shear fracture.It should also be noted that for the example dam the coefficients of variation and negative correlation of the two key mechanical parameters of sliding planes were selected in the relatively wider ranges, which are summarized from the field and experimental data. In the calculation for Slip Path B, four combinations of distributions of the two parameters were detailedly taken into consideration. So the final results obtained for the example dam can reveal the performance about foundation stability of any similar gravity dam in general. We confirm that the randomness off
′ is more significant than c
′, and, what is more, it is reasonable to assume a normal distribution and independence of each parameter in the absence of further information. The assumption will probably overestimate the probability of foundation instability of gravity dam in some degree.


## 4.3.1. Feasibility of DPG Method

Two methods, DPG and MC methods, are used to calculate the instability probabilities of Slip Path B during design reference period. Table6 gives the calculated results obtained by the two methods; when f
′ and c
′ are both normally distributed, their variation coefficients, δ
f
′ and δ
c
′, are 0.2 and 0.3, respectively, and their coefficient of correlation, ρ
f
′
c
′, is 0.Table 6
Comparison of resultant instability probabilities calculated by MC and DPG method.
Method
Design points (sliding plane B2-B3)
β
p
f
Notes
f
1
′
∗
c
1
′
∗/MPa
MC
/
/
3.650
1.310 × 10−4
The numbers of calculating performance function,n
=
1
,
000
,
000, and the variation coefficient for p
f, V
p
f
=
0.0873, when 95% confidence level is specified
DPG
0.107
0.039
3.664
1.244 × 10−4
The numbers of calculating performance function,n
=
36, with the step Δ
x
i
=
0.03
x
i, for divided difference calculationThe instability probabilities of Slip Path B are, respectively,1.310
×
10
-
4 and 1.244
×
10
-
4 calculated by MC and DPG methods, and the absolute and relative errors between the two methods are, respectively, 0.0007% and 5.3%. Thus, the two methods have nearly the same accuracy, but the latter requires less iterating calculation for convergence. Meanwhile, by DPG method the design points can be obtained easily, and the probability distributions, variability, and negative correlation of f
′ and c
′ can be taken into consideration for special research. So the DPG method suggested in the paper has much merit and is rational and feasible for the major research goal of the paper.


## 4.3.2. Analyses of the Effects of Correlation off
′ and c
′ on Instability Probability

Table7 lists the calculated results of the instability probabilities of Slip Path B by DPG method; when f
′ and c
′ are both normally distributed, δ
f
′ and δ
c
′ are 0.2 and 0.3, respectively, and the coefficient of correlation, ρ
f
′
c
′, is within the range of −0.7~0.0. The negative value of ρ
f
′
c
′ indicates the degree of negative correlation of f
′ and c
′. It is clear from Table 7 that the resultant instability probability calculated decreases with the increasing of the absolute value of negatively correlated coefficient. p
f
=
1.244
×
10
-
4 when ρ
f
′
c
′
=
0.0 (independent case), and p
f = 3.106 × 10−7 when ρ
f
′
c
′
=
-
0.7 (strongly negatively correlated case). By comparison, the values of the instability probability under the two cases are different in 3 orders of magnitude. Correspondingly, the values of the reliability index under the two cases are 3.664 and 4.984, respectively. Obviously the independent case is on the conservative side for simulation.Table 7
Calculated results under the conditions of differentρ
f
′
c
′ by DPG method.
Number
ρ
f
′
c
′
Design points (sliding plane B2-B3)
β
p
f
f
1
′
∗
c
1
′
∗/MPa
1
0.0
0.107
0.039
3.664
1.244 × 10−4
2
−0.1
0.100
0.044
3.789
7.554 × 10−5
3
−0.2
0.092
0.050
3.929
4.262 × 10−5
4
−0.3
0.082
0.057
4.086
2.194 × 10−5
5
−0.4
0.071
0.065
4.264
1.005 × 10−5
6
−0.5
0.058
0.075
4.467
3.961 × 10−6
7
−0.6
0.042
0.087
4.704
1.275 × 10−6
8
−0.7
0.022
0.102
4.984
3.106 × 10−7
Note: ifρ
f
′
c
′ is less than −0.8, f
1
′
∗ will be less than 0, and the simulation will be not converged at design points.


## 4.3.3. Analyses of the Effects in Different Distributions off
′ and c
′ on Instability Probability

Figure5 illustrates the resultant instability probability under four different cases of combinations of normal or ln-normal distributions of f
′ and c
′ for δ
f
′
=
0.2 and δ
c
′
=
0
.
3 by DPG method. In all cases, p
f tends to increase as ρ
f
′
c
′ increases from −0.7 (negatively correlated) to 0 (independent). When f
′ is normally distributed, the instability probabilities are approximately identical whether c
′ is normally or ln-normally distributed. When f
′ is ln-normally distributed, the instability probabilities are smaller and the values calculated under the case that c
′ is in ln-normal distribution are less than in normal distribution at the same level of correlation of the two parameters. From the figure, p
f = 1.24 × 10−4 if f
′ and c
′ are both normally distributed, and p
f
=
4.64
×
10
-
9 if the two parameters are both ln-normally distributed, for δ
f
′
=
0.2, δ
c
′
=
0.3, and ρ
f
′
c
′
=
0.0.Figure 5
Relation ofp
f and ρ
f
′
c
′ assuming that f
′ and c
′ follow different distributions.


## 4.3.4. Analyses of the Effects of Variability off
′ and c
′ on Instability Probability

On the basis of the above analyses, it is reasonable, simple, and conservative to assume thatf
′ and c
′ are normally distributed and independent of each other unless otherwise demonstrated. Considering the range of the values of coefficients of variation, δ
f
′
= 0.14~0.26, δ
c
′
= 0.15~0.45, the calculations are made, and the curves about p
f varying with δ
f
′ and δ
c
′ are displayed in Figures 6 and 7. According to the resultant curves, p
f starts to apparently climb as δ
f
′ increases for different δ
c
′ (shown in Figure 6), and, in addition, p
f also tends to increase as δ
c
′ increases for different δ
f
′ (shown in Figure 7). However, p
f
=
1.92
×
10
-
3 if δ
f
′
=
0.26, and p
f
=
4.04
×
10
-
7 if δ
f
′
=
0.14 for δ
c
′
=
0.30; meanwhile, p
f
=
3.18
×
10
-
4 if δ
c
′
=
0.45, and p
f
=
6.08
×
10
-
5 if δ
c
′
=
0.15 for δ
f
′
=
0.20. There is the difference of p
f in 4 orders of magnitude between δ
f
′
=
0.14 and 0.26 for δ
c
′
=
0.30; in contrast, under the latter case, there is the deference of p
f in one order of magnitude between δ
c
′
=
0.15 and 0.45 for δ
f
′
=
0.20. So the influence of the variability of f
′ is more apparent than that of c
′ on p
f despite the fact that the varying range of coefficients of variation of f
′ is less than that of c
′.Figure 6
Relation ofp
f and δ
f
′.Figure 7
Relation ofp
f and δ
c
′.


## 4.3.5. Stability Safety Analyses and Assessments of the Practical Dam

The instability probability of the most adverse slip path of the example dam section (number 6 overflow section) mentioned above is1.244
×
10
-
4 for δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0. It is less than the allowable value of the instability probability, [P] = 5.5 × 10−4, but they are in the same order of magnitude. Meanwhile, it is 2.861 × 10−3 for δ
f
′ = 0.26, δ
c
′
=
0.45, and ρ
f
′
c
′
=
0.0. It is greater than [P] but less than 5 × 10−3, which is defined as unacceptable probability. In contrast, the factor of safety is 2.45 by deterministic method without taking the variation of parameters into consideration. The factor of safety is less than the allowable safety factor, [
K
]
=
3.0; the reliability index is 3.66 for δ
f
′ = 0.20, δ
c
′ = 0.30, and ρ
f
′
c
′
=
0.0, which nearly equals the target reliability index, β
T
=
3.7. Obviously, the conclusions are more or less different following three different safety standards. Generally, it has been verified that the deterministic method is most conservative in the three methods by the analyses of the practical dam in the paper and the conclusions from other literatures as well. From above analyses, the margins of safety of the foundation stability are not adequate for the example dam section and the enhancing treatments should be carried out especially for the fault-fractured zones and weak structural planes.Similarly, the instability probabilities of the other two slip paths can be calculated by the DPG method. Assume thatf
′ and c
′ of each sliding plane of each slip path are normally distributed and independent, and δ
f
′
=
0.20, δ
c
′
=
0.30, and ρ
f
′
c
′
=
0.0 for any slip plane. The resultant instability probabilities are listed in Table 8. The three slip paths can be seen as three main failure modes for the rock foundation of the example dam, and each failure mode is independent due to the independence of each shear strength parameter of any sliding plane. Three failure modes comprise a series failure system of dam foundation. So the overall probability of failure for the overall dam foundation system can be calculated by summing over the instability probabilities of the three slip paths. The final instability probability of the overall rock foundation is 2.053 × 10−4 for a conservative consideration.Table 8
Resultant instability probabilities of three slip paths calculated by DPG method.
Slip paths
Reliability method (DPG)
Deterministic method
β
p
f
K
A
3.772
8.094 × 10−5
2.56
B
3.664
1.244 × 10−4
2.45
C
5.612
9.993 × 10−9
10.71
Note: the deterministic method is equivalence safety factor method without considering the variation of the parameters against shear fracture.It should also be noted that for the example dam the coefficients of variation and negative correlation of the two key mechanical parameters of sliding planes were selected in the relatively wider ranges, which are summarized from the field and experimental data. In the calculation for Slip Path B, four combinations of distributions of the two parameters were detailedly taken into consideration. So the final results obtained for the example dam can reveal the performance about foundation stability of any similar gravity dam in general. We confirm that the randomness off
′ is more significant than c
′, and, what is more, it is reasonable to assume a normal distribution and independence of each parameter in the absence of further information. The assumption will probably overestimate the probability of foundation instability of gravity dam in some degree.


## 5. Conclusions

In this paper the main goal has been to try applying probabilistic approach to calculate the instability probability of gravity dam foundation over multiple sliding planes rather than single or dual sliding plane(s) in traditional analysis. The simple and efficient probabilistic method, DPG method combined with a divided difference approximation to the partial derivatives, is proposed. The effects of randomness of the shear strength parameters of sliding planes for a practical typical gravity dam are analyzed in detail. The following conclusions can be obtained:(1)
By comparison between MC method and DPG method, the two methods have nearly the same accuracy, but the latter needs less cost for calculation. Meanwhile, the DPG method has much other merit. For example, the design points can be obtained easily during the calculating process, and probability distributions, variability, and correlation of variables can be taken into consideration rationally. Hence, the DPG method associated with the method of divided difference is feasible for the calculation of probability of failure with implicit performance function just like the case of instability of dam foundation over multiple sliding planes.
(2)
The randomness of material properties of sliding planes is a major contributor to the uncertainty of the foundation stability of a gravity dam. The randomness of the two parameters against shear fracture,f
′ and c
′, of the sliding planes may be described by assumption of their distributions, variability, and correlation. So the effects of the randomness of the parameters on resultant instability probability of gravity dam foundation can be analyzed based on these aspects. The study of the practical dam shows that (1) when both of the two parameters follow normal distribution, the resultant value of instability probability is greater than that which is calculated when one or both of the parameters do not follow normal distribution, and the distributed type of f
′ is dominant; (2) the stronger the negative correlation of the two parameters is, the smaller the instability probability calculated is; and (3) the instability probability tends to increase as the coefficients of variability of f
′ and c
′ increases, and the influence of the variability of f
′ is more apparent than that of c
′ on the instability probability. So it is reasonable to assume a normal distribution and independence of each parameter in the absence of further information. The assumption will probably overestimate the instability probability of gravity dam foundation. The spatial variation of f
′ is more significant than c
′ on the instability probability. All of the conclusions drawn in the paper are identical to those of other research.
(3)
Currently, there are three kinds of safety indices—factor of safety, reliability index, and probability of failure—for engineering stability problem. The factor of safety is used in deterministic method. It is so convenient for engineers to utilize and have so much experience in designs of many existing gravity dam projects that has become a major index to assess foundation stability of gravity dam for a long time. However, the assessment method based on this index is relatively more conservative and should be improved upon undoubtedly. The trend is application of probabilistic approach, in which the reliability index or probability of failure will be utilized. Now, a fundamental problem is the lack of data. Field exploration is expensive and laboratory tests are always based on small numbers of samples, so there are seldom enough data to support meaningful statistical analyses and broad conclusions about the statistics of rock or soil properties. In a word, the current challenges to the profession are to make use of appropriate probabilistic methods in practice and to sharpen the investigations and analyses so that each additional data point provides maximal information.
(4)
In the paper, we applied a simple, accurate, and efficient method to calculate the instability probability of concrete gravity dam foundation over multiple sliding places. It can also be utilized in other similar analyses of slopes, levees, embankments, earth-rock dams, and so on.



---

*Source: 4264627-2016-06-29.xml*