# Stability analysis of rock slopes subjected to seepage forces using the modified Hoek-Brown criterion

Zied Saada, Samir Maghous, Denis Garnier

► To cite this version:

Zied Saada, Samir Maghous, Denis Garnier. Stability analysis of rock slopes subjected to seepage forces using the modified Hoek-Brown criterion. International Journal of Rock Mechanics and Mining Sciences, 2012, 55, pp.45-54. 10.1016/j.ijrmms.2012.06.010. hal-01157363

HAL Id: hal-01157363

https://enpc.hal.science/hal-01157363

Submitted on 18 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stability analysis of rock slopes subjected to seepage forces using the modified Hoek-Brown criterion

Z. Saadaª, S. Maghousᵇ,*, D. Garnierᶜ

The paper deals with the stability analysis of a water-saturated rock slope by means of the kinematic approach of limit analysis theory. Particular emphasis is first given to the effects of pore water pressure on the global stability of geotechnical structures. The rock strength properties being formulated in terms of effective stresses, it is shown how the effect of seepage flow may be accounted for such an analysis through driving body forces derived from the gradient of excess pore pressure distribution is shown. The latter is obtained as the solution to a hydraulic boundary value problem, and then incorporated as external loading in the subsequent stability analysis. The rock strength properties are modeled by a modified Hoek-Brown failure criterion, for which closed-form expressions of the support functions have been derived in a previous paper. The approach is then applied to investigate the stability of rock slopes. Computational results are given, providing ample evidence of the destabilizing effects induced by the seepage forces. The influence of relevant geometrical, strength and loading parameters is also discussed.

## 1. Introduction

Assessment of the bearing capacity of rock structures, such as dam foundations, tunnels or rock slopes built in rock masses whose strength properties are described by traditional cohesive-frictional failure conditions, is a well-known classical problem in geotechnical engineering. Many works have devoted to this purpose with relatively satisfactory solutions in the absence of saturating fluid flow. However, such a problem has received less attention as soon as one deviates from this reference situation. This is notably the case when the rock structure is subjected to seepage forces which are likely to significantly reduce its stability, as has been early shown by several experimental or theoretical works. Pioneer works on this topic are reported for instance in Leca et al. [1].

Many of the theoretical investigations have been based on the assumption that the rock strength is governed by the linear Mohr-Coulomb failure criterion. It is however well-known that the strength envelopes of almost all geomaterials are non-linear, as it has been emphasized by several experimental works [2-8]. Among the non-linear failure criteria proposed in the literature,

the Hoek-Brown failure criteria are considered to reasonably well model the strength capacities of isotropic rocks. Nevertheless, few studies were devoted to the assessment of bearing capacity of structures in rocks modeled by the latter strength criteria. In this respect, one should quote the pioneer contributions of Baker and Frydman [9], Zhang and Chen [10], Drescher and Christopoulos [11] or Serrano and Olalla [12-14]. These fundamental contributions have inspired a series of works that were developed subsequently [15-19].

The problem addressed in this paper is related to the analysis of slopes stability when the rock materials exhibit non-linear failure envelope. In this context, several works have been investigated within the framework of a Hoek-Brown failure criterion: the problem of rock slopes stability and induced reduction in bearing capacity of foundations lying near rock slopes [20-22].

Less is known about the problem of reduction in slope stability due to the presence of seepage forces. More generally, stability analysis of geotechnical structure in the presence of seepage forces still remains a key challenge in soil and rock engineering. Most of the contributions to the field are in fact concerned with the specific problem of stability analysis of soil slopes. From a theoretical point of view, the pore water pressures distribution to be incorporated in the stability analysis is obtained from ground water conditions. In the context of slope stability, this distribution is classically simulated by a simplified flow net or a pore water pressure ratio. Pore water pressure may be considered as internal

<footnote>ªLaboratoire de Génie Civil, ENIT, Université de Tunis El Manar, Tunis, Tunisia</footnote>

<footnote>ᵇDepartment of Civil Engineering, Federal University of Rio Grande do Sul, Porto Alegre-RS, Brazil</footnote>

<footnote>ᶜLMSGC, Ecole Nationale des Ponts et Chaussées, Marne-la Vallée, France</footnote>

<footnote>ªCorresponding author. Tel.: +55 51 33083588; fax: +55 51 33083999.
E-mail addresses: samir.maghous@ufrgs.br, maghous@ppgec.ufrgs.br (S. Maghous).</footnote>
forces in the analysis of slope stability, as used for instance by Miller and Hamilton [23,24] which assumed the pore water distribution to be hydrostatic below a free water table. On the other hand, pore pressure can be included as external forces as developed in Michalowski [25] in the context of a Mohr-Coulomb strength criterion. The analysis developed in this work for slope stability implemented a log-spiral rotation failure mechanism, where the distribution of pore pressure along the failure surface was described by a pore pressure coefficient, originally introduced by Bishop and Morgenstern [26]. Later, Yang and Zou [27] investigated the same problem adopting the same concepts as the work of Michalowski [25] but the rock strength properties were modeled by means of a modified Hoek-Brown failure criterion. However, the main drawbacks, as clearly underlined in the paper, is its domain of applicability reduced to rigid-block failure mechanisms. It cannot thus be applied in the framework of “punching”-like failure mechanisms, which prove to be generally relevant for the problem of bearing capacity, since they involve deforming zones.

The present contribution is divided into two main sections. It is shown in Section 2 that, since the rock failure criterion is classically expressed as a function of the effective stress, the stability analysis of geotechnical structures subject to water flow involves a field of body forces derived from the gradient of excess pore pressure distribution. Provided certain conditions be fulfilled, the latter may be calculated as the solution to a simplified hydraulic problem in which a steady-state flow is assumed and the rock deformations are neglected. In the second part of the paper, the approach is applied to investigate the stability of rock slopes in the presence of seepage forces. The distribution of these seepage forces are numerically evaluated through a finite element method and incorporated into the stability analysis based on the implementation of a rotational log-spiral failure mechanisms.

## 2. Stability analysis in the framework of effective stresses

The objective of this section is to present a general framework, based on the validity of effective stress concept, for the stability analysis of geotechnical structures in the presence of seepage forces. Although the approach is illustrated on the case of a rock slope stability, the method that is proposed to account for the effect of seepage flow can be applied to more general situations, such as dam foundations, tunnels driven under the water table, etc.

### 2.1. Formulation of limit analysis problem

The plane strains stability analysis considered herein refers to a homogeneous and isotropic rock slope with angle $\beta$ and height

$H$, as shown in Fig. 1. In the sequel, the considered rock domain is defined by volume $\Omega$ and its boundary $\partial\Omega$.

In addition to mechanical loading (gravity forces and boundary conditions), the material system is also submitted to hydraulic force originating from pore water (or other saturating fluid) pressure and associated flow through the slope. It is observed that the particular hydraulic boundary conditions shown in Fig. 1 are only provided for illustrative purpose, without any loss of generality of the approach described below.

The starting point is the hydraulistic distribution for pore pressure $p_w = \gamma_w x_1$ prevailing along the slope in the absence of fluid flow, $\gamma_w$ referring to the water unit weight. Any water flow regime will induce through the rock slope a distribution of excess pore pressure defined as

$$u = p - p_w = p - \gamma_w x_1 \quad (1)$$

where $p$ is the pore water pressure at the considered point. The spatial variation of $u$ is a driven force of the hydraulic flow network defined by the filtration velocity field:

$$v = -k \cdot \text{grad} \frac{u}{\gamma_w} \quad (\text{Darcy's law}) \quad (2)$$

where $k$ is the permeability tensor (which has the dimension of a velocity $\overline{v}$). Introducing the strength capacity of the saturated rock material on the form $F(\sigma, p) \le 0$, the condition for the rock structure to remain safe under the prescribed loading can be formulated as follows:

$$\text{Stability} \Leftrightarrow \exists \sigma' \mid \left\{ \begin{array}{ll} \text{div } \sigma + \gamma e_1 = 0 & \text{on } \Omega \\ \frac{\sigma \cdot n}{\sigma} = \frac{T}{F(\sigma, p)} \le 0 & \text{on } \partial\Omega \quad (\text{stress-boundary conditions}) \\ \frac{\sigma'}{F(\sigma', p)} \le 0 & \text{on } \Omega \end{array} \right. \quad (3)$$

in which $\gamma$ denotes the rock unit weight. The prescribed stress boundary $T^d$ naturally depends on both the mechanical and the hydraulic loading conditions.

Let assume from now on that the strength capacities of the fluid saturated rock are controlled by the effective stress $\sigma' = \sigma + p1$

$$F(\sigma, p) = F(\sigma') \le 0 \quad (4)$$

The stability condition (4) is thus conveniently rewritten in terms of effective stresses as follows:

$$\text{Stability} \Leftrightarrow \exists \sigma' \mid \left\{ \begin{array}{ll} \text{div } \sigma' + \gamma' e_1 - \text{grad } u = 0 & \text{on } \Omega \\ \frac{\sigma'}{F(\sigma')} \cdot n = \frac{T}{F(\sigma')} \le 0 & \text{on } \partial\Omega \\ \frac{\sigma'}{F(\sigma')} \le 0 & \text{on } \Omega \end{array} \right. \quad (5)$$

where $\gamma' = \gamma - \gamma_w$ is the buoyant weight and $T'$ represents the distribution of effective stresses. Eq. (5) expresses in the context

<figure><img src="image_2.png" /><figcaption>Fig. 1. Problem geometry and loading mode.</figcaption></figure>
of effective stresses the stability conditions for the rock slope. It emphasizes that the seepage forces may be regarded as external body forces equal to -grad u, which is directly proportional to the filtration velocity -grad u = γ<sub>w</sub>k<sup>-1</sup>·v. The stability of the rock slope should thus be analyzed under its weight minus the buoyancy force, to the distribution of effective stress vectors T' on its boundary, and to the seepage body forces. From the stability analysis viewpoint, the loading mode associated with the framework of effective stresses can symbolically be represented by the prescribed set (γ', T'<sup>d</sup>, -grad u).

The dual form of the stability problem (5) is expressed by means of the kinematic approach of limit analysis. It is based upon the upper bound theorem [28,29] which implementation relies on the fundamental inequality

$$P_e(U) \leq P_{mr}(U) \quad \forall U \quad (6)$$

expressing a necessary condition for the structure to remain stable under the external loading. In the above inequality (6), U is any virtual kinematically admissible velocity field (referred to as failure mechanism in the sequel), P<sub>e</sub>(U) denotes the work done by the external forces

$$P_e(U) = \int_{\Omega} \gamma' e_1 \cdot \underline{U} \, d\Omega + \int_{\partial\Omega} T' \cdot \underline{U} \, dS + \int_{\Omega} -\text{grad} \, u \cdot \underline{U} \, d\Omega \quad (7)$$

accounted for through the concept of pore pressure coefficient, which is shown to be a convenient way to include pore pressure effects in the context of rigid-block failure mechanism. The general expression of the maximum resisting work P<sub>mr</sub>(U) developed in the failure mechanism U is

$$P_{mr}(U) = \int_{\Omega} \Pi[d(\underline{x})] d\Omega + \int_{\Sigma} \Pi[v(\underline{x}); [U(\underline{x})]] d\Sigma \quad (8)$$

where d is the strain rate field associated with U, [U] is the jump of U at a point x when crossing a possible velocity discontinuity surface Σ following its normal v(x). The Π-functions are the support functions and are defined by duality from the strength condition F(σ') ≤ 0:

$$\Pi[d] = \sup_{\underline{\sigma}'} \{\underline{\sigma}': d|F(\underline{\sigma}') \le 0\} \quad (9)$$

$$\Pi[v; [U]] = \sup_{\underline{\sigma}'} \{[U] \cdot \underline{\sigma}'; v|F(\underline{\sigma}') \le 0\} \quad (10)$$

We shall now specify the strength capacities of the constitutive rock through the definition accounted for through the concept of pore pressure coefficient, which is shown to be a convenient way to include pore pressure effects in the context of rigid-block failure mechanism of the yield condition F(σ') ≤ 0.

## 2.2. Modified Hoek-Brown failure condition

One of the main characteristics of rock masses is the existence at different scales of discontinuities of various sizes and orientations, usually referred to as joints. At a macroscopic scale, the rock material may be regarded as a homogeneous medium with anisotropic properties, whose strength capacities may be assessed experimentally or using upscaling approaches [30–33].

In the present analysis, the strength properties of the constitutive material are assumed to be isotropic. At the macroscopic scale, this assumption seems reasonable when the rock is intact or heavily jointed.

In the present contribution, the strength capacities of the rock material are modeled by a modified Hoek-Brown failure condition [34]

$$F(\underline{\sigma}') = \underline{\sigma}'_1 - \underline{\sigma}'_3 - \sigma_c \left( -m \frac{\underline{\sigma}'_1}{\sigma_c} + s \right)^n \le 0 \quad (11)$$

$\sigma'_1$ and $\sigma'_3$ are respectively the major and minor effective principal stresses (stresses are counted positive in tensile), and $\sigma_c \ge 0$ is the uniaxial compressive strength of rock. Parameters m, s and n depend on the geological strength index (GSI) [34] and take the following form:

$$\frac{m}{m_i} = \exp\left(\frac{GSI-100}{28-14D_0}\right) \quad (12)$$

$$s = \exp\left(\frac{GSI-100}{9-3D_0}\right) \quad (13)$$

$$n = \frac{1}{2} + \frac{1}{6} \left[ \exp\left(-\frac{GSI}{15}\right) - \exp\left(-\frac{20}{3}\right) \right] \quad (14)$$

where D<sub>0</sub> is a disturbance coefficient that varies from 0 for undisturbed in situ rock masses to 1 for very disturbed rock masses. The value of m<sub>i</sub> is obtained from compression tests on intact specimen of intact rock. Approximate values of this parameter are provided in Hoek [35] for some typical rocks.

The implementation of the kinematic approach requires the computation of the support functions defined by (9) and (10) for the modified Hoek-Brown criterion. These functions have been derived in [14]

$$\Pi[\underline{d}] = \frac{s\sigma_c}{m} \text{tr} \underline{d} + \sigma_c(n^{n/(1-n)} - n^{1/(1-n)}) \frac{m\mathcal{M}(\underline{d})}{\text{tr} \underline{d}} \quad \text{if } \text{tr} \underline{d} > 0 \quad (15)$$

where function M(d) is defined by the relation

$$\mathcal{M}(\underline{d}) = [\max(0, -d_1) + \max(0, -d_2) + \max(0, -d_3)]^{1/n} \quad (16)$$

in which d<sub>1</sub>, d<sub>2</sub> and d<sub>3</sub> represent the eigenvalues of d.

As regards the Π-function relative to a velocity jump defined by (10), it has been found that

$$\Pi[v; [U]] = \frac{s\sigma_c}{m} [U] \cdot \underline{v} + \sigma_c \left( n^{n/(1-n)} - n^{1/(1-n)} \right) \left( \frac{m\mathcal{M}(v; [U])}{[U] \cdot \underline{v}} \right)^{n/(1-n)} \quad \text{if } [U] \cdot \underline{v} > 0 \quad (17)$$

with

$$\mathcal{M}(v; [U]) = \frac{1}{2^{1/n}} (|[U]| - [U] \cdot v)^{1/n} \quad (18)$$

Conditions tr d ≥ 0 in (15) and [U] · v > 0 in (17) respectively express that Π[d] < +∞ and Π[v; [U]] < +∞. These conditions are necessary for the kinematic approach (6) results in non-trivial upper bound solutions.

# 3. Application to slope stability analysis

The purpose of the following section is to apply the kinematic approach of limit analysis described in the previous section to assess the stability of a rock slope (see Fig. 1) under loading of gravity and seepage forces. As a matter of fact, the latter ones appear as a component of the external loading of the rock structure.

In order to derive quantitative stability conditions in a given practical situation, the field -grad u representing the seepage forces must be first evaluated.

## 3.1. The hydraulic problem

Adopting the classical framework in which the coupling between the skeleton strains and pore pressure (hydromechanical
coupling) is disregarded, the fluid balance reads

$$ \mathrm{div} \underline{v} = 0 \qquad (19) $$

which can be rewritten by virtue of Darcy's law

$$ \mathrm{div} \left( -k \cdot \mathrm{grad} \frac{u}{\gamma_w} \right) = 0 \qquad (20) $$

The equations of the differential problem to be satisfied by the field of excess pore pressure $u$ are defined by (20) and the appropriate hydraulic boundary conditions, which naturally depend on the problem under consideration. The solution $u$ of the above mentioned problem may be sought numerically using commercial or ad hoc softwares specifically devised for this purpose [36].

For illustrative purposes, let us examine the particular seepage problem defined in Fig. 2. This case will be reconsidered for the rock slope stability analysis of Section 3.2. The pore water pressure is equal to zero (i.e., $p=0$) along the boundary surfaces $TT'$, $TS$ and $SS'$. The corresponding boundary conditions expressed in terms of excess pore pressure $u = p - p_w = p - \gamma_w x_1$ are

$$ \begin{cases} u = 0 & \text{along } TT' \text{ (top surface of the slope)} \\ u = -\gamma_w x_1 & \text{along } TS \\ u = -\gamma_w H & \text{along } SS' \text{ (bottom surface of the slope)} \end{cases} \qquad (21) $$

Due to the difference in hydraulic head between the water tables respectively located at the top and the toe of slope, a

ground water flow occurs through the slope. If the rock medium is homogeneous with principal values of permeability equals to $k_h = k_2 = k_3 = 3$ and $k_\nu = k_1 = 1$, the solution $u$ to the above hydraulic problem depends solely on the ratio $k_h/k_\nu$ between the horizontal and vertical permeabilities, and not on their individual magnitudes.

This above elliptic problem formed by Eqs. (20) and (21) is well-defined, and its solution can be sought numerically.

Fig. 3 shows an example of FE analysis of the problem: FE mesh used for the hydraulic problem (Fig. 3a), contours of excess pore pressure $u$ (Fig. 3b) and associated distribution of seepage forces -grad $u$ (Fig. 3c). These results have been obtained numerically using a FE freeware assuming isotropic permeability (i.e., $k = k_1$) for the rock.

$\overline{\overline{\mathrm{grad}}} \overline{\overline{u}}$ alternative way to evaluate the effects of pore water pressure is provided within the approximate framework described by Bishop and Morgenstern [26]. In the latter approach, the pore water pressure is conveniently given by expression

$$ p = r_p \gamma_h h \qquad (22) $$

where $h$ is the (vertical) depth of the point located on the failure surface below the rock surface, $\gamma$ denotes the unit weight of rock/ soil and $r_p$ is usually referred to as the pore pressure coefficient.

Most of the stability analyses involving pore water pressure effects are based on the concept of pore pressure coefficient [25-27]. Besides, the majority of these analyses assume a constant pore pressure coefficient $r_p$ throughout the cross-section,

<figure><img src="image_3.png" /><figcaption>Fig. 2. Hydraulic problem and associated boundary conditions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The hydraulic problem: example of finite element simulation. (a) FE mesh, (b) contours of pore pressure, (c) seepage forces.</figcaption></figure>
with values generally below 0.3. As one could naturally expect and as was originally emphasized by Bishop and Morgenstern [26], $r_p$ varies over the whole cross section involved within the failure mechanism of rock slope. These authors stated that, in most slope stability problems, an average value of $r_p$ can readily be calculated and used in the stability analysis with little loss in accuracy at the global level. Based on this idea, most of the traditional approaches to rock slope stability adopt the concept of pore pressure coefficient in the analysis. This simplified kind of approach is however not suitable in practice since the value of $r_p$ to be used for a particular problem is not known.

The whole stability analysis presented in the sequel is based on the rigorous approach of considering the pore water pressure distribution derived from a finite element code.

## 3.2. Rotational failure mechanism

Once the hydraulic problem is solved, the field of seepage forces is known and the stability analysis can be undertaken. Simple dimensional analysis arguments show that the stability of the rock slope is characterized by means of the following non-dimensional parameters: $\beta$, $\gamma H/\sigma_c$, $m_i$, GSI, $D_0$, $k_h/k_v$. They are related to slope inclination, gravity level, rock strength parameters and permeability anisotropy.

We apply herein the kinematic approach to address the classical problem of rock slope stability submitted to a loading mode defined by gravity forces and seepage forces. The kinematic approach is based on the rotational failure mechanism, shown in Fig. 4, which is usually employed for homogeneous Mohr-Coulomb soil or rock slope. In such a mechanism a volume of rock mass is rotating about a point $\Omega$ with an angular velocity $\omega$. The curve $I_1I_2$ separating this volume from the rest of the structure which is kept motionless is an arc of log-spiral of angle $\varphi$ and focus $\Omega$. It follows necessary that the velocity jump at any point of this line is inclined at angle $\varphi$ with respect to the tangent at the same point.

Two classes of failure mechanisms are considered in the analysis. Fig. 4a displays the situation of log-spiral surfaces

emerging above the slope toe (i.e., $I_2$ is located on the slope wall), whereas Fig. 4b shows the situation where the failure surface is extending below the slope toe (i.e., $I_2$ lies on the plane $x_1 = H$).

The velocity field within the volume in motion is orthonormal with expression on the form $U = \omega r e_0$. The distance $r_0 = \Omega I_1$ defines the radius of the log-spiral curve for $\theta = \theta_1$. Accordingly, this curve can be defined in polar coordinates $(\Omega, r, \theta)$ by equation $r = r_0 e^{(\theta - \theta_1) \tan \varphi}$. Such a failure mechanism involves four parameters: (i) three angles $(\theta_1, \theta_2, \varphi)$ and distance $r_0 = \Omega I_1$ in the case of log-spirales passing above the slope toe, or (ii) four angles $(\theta_1, \theta_2, \varphi, \beta')$ in the case of log-spirales passing below the slope toe.

In order to avoid heavy mathematical developments and for sake of clarity, the details of the approach will be presented in the sequel restricting ourself to the situation of log-spirales passing above the slope toe. The numerical results given in Section 4 are however obtained by considering both classes of log-spirales.

We shall recall that the loading is the conjunction of gravity (minus buoyancy) and seepage forces. Accordingly, the work done by the external forces in such failure mechanism is the sum of two contributions

$$P_e(U) = P_{\gamma'}(\underline{U}) + P_u(\underline{U}) \quad (23)$$

$$P_{\gamma'}(\underline{U}) = \int_{\Omega} \gamma' e_1 \cdot \underline{U} \, d\Omega \quad \text{and} \quad P_u(\underline{U}) = \int_{\Omega} -\text{grad } u \cdot \underline{U} \, d\Omega \quad (24)$$

It should be recalled that the seepage forces are in the limit analysis problem considered as given. As previously mentioned, the field of $u$ (or grad $u$) can be obtained either by solving the hydraulic problem via a finite element procedure or adopting the simplified formula (22) $p = r_p\gamma h$.

The first contribution $P_{\gamma'}$ to external work may be written as

$$P_{\gamma'} = r_0^3 \gamma' \omega (f_1 + f_2 + f_3) \quad (25)$$

in which the coefficients $f_1, f_2$ and $f_3$ are non-dimensional functions of parameters $\theta_1, \theta_2, \varphi$ and $r_0$. The corresponding

<figure><img src="image_5.png" /><figcaption>Fig. 4. Rotational failure mechanism for rock slope. (a) Log-spiral emerging on the slope wall, (b) log-spiral extending below the slope toe.</figcaption></figure>

<figure><img src="image_6.png" /></figure>
expressions are provided below:

$$f_1 = \frac{e^3 \tan \varphi (\theta_2 - \theta_1) (3 \tan \varphi \cos \theta_2 + \sin \theta_2) - (3 \tan \varphi \cos \theta_1 + \sin \theta_1)}{3(9 \tan^2 \varphi + 1)} \quad (26)$$

$$f_2 = -\frac{1}{6} \sin \theta_1 \left[ 2 \left( \frac{TI_1}{r_0} \right) \cos \theta_1 - \left( \frac{TI_1}{r_0} \right)^2 \right] \quad (27)$$

$$f_3 = -\frac{1}{3} \left(\frac{r_T}{r_0}\right)^3 \sin^3(\theta_T + \beta) \\ \times \left[ \frac{\cos \beta}{2} \frac{1}{\sin^2(\theta_T + \beta)} - \frac{1}{\sin^2(\theta_2 + \beta)} \right) \\ + \sin \beta \left( \frac{\cos(\theta_T + \beta)}{\sin(\theta_T + \beta)} - \frac{\cos(\theta_2 + \beta)}{\sin(\theta_2 + \beta)} \right] \quad (28)$$

where distance $TI_1$ is given by

$$TI_1 = r_0 \frac{\sin(\theta_1 + \beta) - e^{\tan(\varphi)(\theta_2 - \theta_1)} \sin(\theta_2 + \beta)}{\sin \beta} \quad (29)$$

and polar coordinates ($r_T, \theta_T$) of point T are

$$r_T = r_0 \sqrt{1 - 2 \frac{TI_1}{r_0} \cos \theta_1 + \left( \frac{TI_1}{r_0} \right)^2} \theta_T = \arctan \left( \frac{\sin \theta_1}{\cos \theta_1 - TI_1 / r_0} \right) \quad (30)$$

$\theta_c$ is the angle defining point C. It is the single solution within interval $[\theta_1, \theta_2]$ of the following equation:

$$\cos \theta_c e^{\theta_c \tan \varphi} + \left[ \frac{TI_1}{r_0} - \cos \theta_1 \right] e^{\theta_1 \tan \varphi} = 0 \quad (31)$$

which is solved numerically.

As regards the contribution of seepage forces $P_u$ to external work, it is first recalled that the distribution pore pressure $p$ as well as the density vector of seepage forces are computed at each point of the discretized domain of rock medium by means of the hydraulic software. The volume integral $P_u(U)$ defining the work of seepage forces is then computed by integrating over all elements or sub-elements of the mesh that are included within the zone in motion.

If the simplified form (22) for pore water pressure is adopted, i.e., $p = r_p \gamma h$ (with the assumption that $r_p$ is constant), the contribution $P_u$ can be evaluated analytically

$$P_u = r_0^3 \gamma_w \omega (f_1 + f_2 + f_3) + r_p r_0^3 \gamma \omega f_4 \quad (32)$$

leading to the following expression for the work of external forces:

$$P_e = r_0^3 \gamma \omega (f_1 + f_2 + f_3 + r_p f_4) \quad (33)$$

The expression of the non-dimensional parameter $f_4$ as a function of the geometrical parameters defining the failure mechanism is not given herein in order to avoid introducing additional heavy mathematical expressions.

On the other hand, the maximum resisting work results from the velocity jump along the log-spiral line $l_1 l_2$

$$P_{mr} = \sigma_c \omega r_0^2 f_5 \quad (34)$$

where expression of the non-dimensional function $f_5$ reads

$$f_5 = \frac{1}{2} \left[ \frac{s}{m} + \left( n^{n/(1-n)} - n^{1/(1-n)} \right) m^{n/(1-n)} \left( \frac{1 - \sin(\varphi)}{2 \sin(\varphi)} \right)^{1/(1-n)} \right] \left( e^{2 \tan(\varphi)(\theta_2 - \theta_1)} - 1 \right) \quad (35)$$

The fundamental kinematic inequality (6) implies that the stability of the rock structure can be characterized by the safety factor

$$\Gamma = \frac{P_{mr}}{P_e} = \frac{P_{mr}}{P_{\gamma'} + P_u} \quad (36)$$

which should comply with the necessary condition

$$\Gamma \ge 1 \quad (37)$$

The above necessary condition of stability implicitly requires that the work done by the external forces is positive $P_e(U) > 0$, otherwise the approach would be trivial.

The critical value of the safety factor is obtained by minimization with respect to the set of variables $(\theta_1, \theta_2, \varphi, r_0)$

$$\Gamma \le \min_{\theta_1, \theta_2, \varphi, r_0} \frac{P_{mr}}{P_e} \quad (38)$$

The minimization parameters are subjected to the following constraints:

$$\begin{cases} 0 < \varphi < \frac{\pi}{2}, & 0 < \theta_1 < \theta_2 < \pi - \beta \\ 0 < \sin(\theta_1 + \beta) - e^{(\theta_2 - \theta_1)\tan \varphi} \sin(\theta_2 + \beta) & \\ 0 < e^{(\theta_2 - \theta_1)\tan \varphi} \sin \theta_2 - \sin \theta_1 \le \frac{H}{r_0} & \end{cases} \quad (39)$$

### 3.3. Comments on the kinematic approach

Before the numerical implementation of the kinematic approach described previously, some comments deserve to be given:

(i) Actually, the kinematic approach presented in the previous section investigates a particular failure mechanism. It is based on the assumption that the inclination of velocity field with respect to the failure surface remains unchanged, equal to an arbitrary angle value $\varphi$, which obviously implies that the failure surface is in fact a log-spiral. Angle value $\varphi$ is therefore a purely geometric parameter that arises from the assumption made on the shape of the failure surface. In particular, it holds no connection with the concept of tangential friction angle of the rock material. Unlike the approaches based on the concept of tangential Mohr-Coulomb failure condition (see for instance [37]), one main specificity of the present approach is its ability to preserve the original modified Hoek-Brown failure criterion. This is made possible because closed-form expressions for the support functions $\pi$ have been derived, allowing the explicit calculation of the maximum resist in work in any virtual velocity field. Clearly enough, the kinematic approach can theoretically be improved by optimization with respect to the shape of failure surface, and not restrict to log-spiral curves. However, this would be a sophisticated task that lies beyond the scope of the present analysis.

(ii) The safety factor $\Gamma$ defined by (36) and characterizing the stability of the rock slope is a function of non-dimensional parameters introduced at the beginning of Section 3.2

$$\Gamma = \tilde{\Gamma}(\beta, \gamma H / \sigma_c, m_i, GSI, D_0, k_h / k_v) \quad (40)$$

(iii) For specific needs in soil mechanics applications, the safety factor is traditionally defined as the ratio between maximal and mobilized shear stresses. However, it should be emphasized that both the safety factor defined in such way and that introduced in the present approach by (36) can equivalently be used for the characterization of failure or stability of a given structure (see for instance [38,39]).

## 4. Computational results

This section provides numerical results for the safety factor of the rock slope under gravity and seepage forces derived from the
rotational failure mechanisms sketched in Fig. 4. For each selected model data, the numerical value of the upper bound solution has been obtained through constrained minimization procedures based on the non-linear sequential quadratic programming algorithm. As mentioned before, the seepage forces are computed by means of a finite element procedure, and then incorporated in the stability analysis as known external forces. The numerical analysis undertaken herein consists in selected examples of rocks and focuses on the effects of the strength, loading and geometry parameters on the safety factor. In particular, the effects of slope height H, slope inclination β, rock unit weight γ and parameters (GSI, mᵢ, σₑ) defining the generalized Hoek-Brown criterion are analyzed in the sequel.

## 4.1. First example: finite element approach versus simplified approach

The aim of this section is to compare the two modes of evaluating the contribution Pᵢ of seepage forces to the work done by external forces: the finite element computations and the simplified method based on the concept of pore pressure coefficient rₚ. For this purpose, we consider a particular problem of rock slope stability for which the kinematic method is implemented to evaluate the safety factor using both finite element and simplified approaches to compute the pore water pressures and seepage forces. The objective is two-fold: (1) to illustrate the crucial impact of seepage forces on the safety factor, and (2) to compare the two methods of computing the effects of pore water pressures. The model data for considered problem are β = 60°, GSI = 30, mᵢ = 17, σₑ = 10 MPa, γ = 20 kN/m³ and D₀ = 0. The numerical analysis is undertaken for several values of the rock slope height H. The results are displayed in Fig. 5.

As expected, Γ is a decreasing function of H. It could be observed that a mean value of 0.1 for the pore pressure coefficient rₚ leads to results similar to those obtained by finite element evaluation of seepage forces. This is consistent with the statement of Bishop and Morgenstern [26] that an average value of rₚ can readily be calculated and used in the stability analysis with little loss in accuracy at the global level.

Increasing the value of rₚ from 0.1 to 0.2, which expresses a higher intensity of water flow, results in a significant decrease of the safety factor for any value of the slope height. This provides clear evidence of the crucial role played by seepage forces in the rock stability.

However, the main limitation of the simplified method lies on the fact that the value of rₚ to be considered for a given problem is not a priori known.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Safety factor as a function of rock slope height: FE method versus simplified approach for seepage forces evaluation.</figcaption></figure>

**Remark.** In all what follows, only the contribution of seepage forces resulting from finite element calculations is used to evaluate the safety factor of rock slope.

## 4.2. Further illustrative examples

A series of numerical simulations are performed and corresponding results are presented in the situation of undisturbed in situ rock medium, that is D₀ = 0, and exhibiting isotropic permeability, that is kₕ/kᵥ = 1. These assumptions on strength and hydraulic condition properties are not the limitations for the applicability of the approach. Hence, the non-dimensional form of the safety factor will read under these assumptions

$$ \Gamma = \tilde{\Gamma}(\beta, \gamma H / \sigma_c, m_i, \text{GSI}) \quad (41) $$

It should be emphasized that the objective of the present study is not to provide charts for practical use in rock engineering, but only to illustrate the capabilities of the proposed method to deal with stability analysis problems involving seepage forces.

Still, the parametric simulations presented in the sequel are intended to give, in a particular configuration, some preliminary insights on the individual impact of slope inclination and strength properties on the safety factor of the rock slope.

More precisely, we consider the particular configuration defined by the fixed value of a dimensional parameter γH/σₑ = 1.2 × 10⁻². Fig. 6 shows the variations of the ratio Γ/Γ₀ as a function of the rock slope inclination β. For each value of β and strength parameters (mᵢ, GSI), Γ₀ refers to the safety factor of the rock slope in the absence of interstitial fluid (dry case). The whole calculations consisted in 120 evaluations of the rock slope safety factor, namely 60 configurations with seepage forces and 60 configurations under dry conditions. The results show that, even though Γ and Γ₀ strongly decrease as the inclination β increases, the variations of ratio Γ/Γ₀ with respect to β remain moderate. This may be a useful indication for future studies aiming at capturing quantitative effects of forces in rock slopes stability.

### 4.2.1. Comments

The parametric simulations undertaken and presented in the previous subsections should be completed by the following observations:

* Clearly enough, the comments related to the effects of any of the considered parameters on the safety factor correspond to preliminary results and should thus be considered with caution. They are only upper bound estimates of the safety factor, based on a particular failure mechanism, and must thus be interpreted as such. Further investigations regarding the dependence of the safety factor with respect to the geometry, loading and strength parameters are still to be made.

* Even it is achieved using an indirect procedure, comparisons have been performed between results derived from the present approach and upper bound solutions for the critical height of the rock slope obtained by Yang and Zou [27]. The latter approach consisted first in replacing the original Hoek-Brown strength criterion by an "optimal" tangential Mohr-Coulomb domain. Secondly, the kinematic method of limit analysis is implemented making use of a rotational log-spiral failure mechanism quite similar to that displayed in Fig. 4. Lastly, the generalized tangential technique is then used to evaluate the optimal tangential Mohr-Coulomb domain and the corresponding upper bound rock slope height. It should be, however, underlined that these authors modeled the hydraulic loading by means of prescribed values of the water
<figure><img src="image_8.png" /><figcaption>Fig. 6. Normalized safety factor versus rock slope inclination.</figcaption></figure>

pressure coefficient $r_p$. For a given set of problem data, including the value of $r_p$, the estimated value for critical rock slope height obtained in [27] is adopted as input for the present approach, and thus the upper bound estimate for the safety factor as defined by (38) is evaluated. In all considered configurations, it was found that $\Gamma \simeq 1$, which shows that the two approaches lead to similar results and consistent with the fact they are both based on a similar failure mechanism.

* In most general cases of hydraulic flow through rock slopes, the field of pore water pressure can theoretically be highly heterogeneous, and cannot thus be represented by a simple linear function of the vertical distance from the point on the failure surface to the slope surface $p = r_p \gamma h$. The distribution of pore pressure along the ground near the rock slope should therefore be properly evaluated, using for instance numerical procedures. However, following in this respect the pioneering

work due to Bishop and Morgenstern [26], and as confirmed by the preliminary analysis presented in Fig. 5, one could reasonably assume that for each hydraulic configuration of the problem, there always exists a constant (average) value of $r_p$ which can accurately represent the effect of seepage forces in the stability problem handled in the context of rotational failure mechanisms.

* A fundamental limitation of the stability analysis based on the concept of pore pressure coefficient lies on the fact that its applicability is reduced to piecewise rigid-block failure mechanisms (i.e., piecewise rigid body motion of elementary rock blocs), since in such approach the pore pressure is only prescribed along the failure surface. If a failure mechanism involving rock zone undergoing deformation, i.e., where the strain rate field associated with the considered velocity is not equal to zero, the pore pressure distribution should be necessarily known in order to
compute the term $P_u(U) = \int_{\Omega} -\mathrm{grad} \cdot u \cdot \underline{U} \, d\Omega$. This issue can be readily understood as follows.

Observing that along any continually deforming rock sub-domain: grad $u \cdot \underline{U} = \mathrm{div}(u \underline{U}) - u \mathrm{div} \underline{U}$, and using the Gauss theorem, the work of seepage forces for the saturated rock domain can be written as

$$P_u(U) = \int_{\Omega} u \mathrm{div} \underline{U} \, d\Omega + \int_{\Sigma} u[\underline{U}] \cdot \underline{n} \, d\Sigma \quad (42)$$

where $[\underline{U}]$ is the jump of $\underline{U}$ and $\underline{v}$ is the unit normal vector to discontinuity surfaces $\underline{\Sigma}$.

The above expression of $P_u$ emphasizes that if the considered velocity field complies everywhere with condition $\mathrm{div} \underline{U} = 0$, which is in particular the case for rigid body motions, the knowledge of the pore pressure field is required only along the failure surfaces. In this situation, the approach based on the concept of pore pressure coefficient seems appropriate and can be conveniently used in the stability analysis. In contrast, if the considered failure mechanism involve rock zones where $\mathrm{div} \underline{U} \neq 0$, which is the case for any velocity field with expansion zones (i.e., where $\mathrm{tr} \underline{d} \neq 0$), the pore pressure distribution is required in such zones. The approach based on the concept of pore pressure coefficient cannot therefore be used, and finite element solutions for $u$ are more appropriate.

## 5. Conclusion

The kinematic approach of limit analysis theory has been implemented to investigate the stability of rock slope subjected to pore water pressure. For this purpose, a comprehensive framework has been first set up for analyzing the stability of geotechnical structures in the presence of seepage force. The rock strength properties being formulated in terms of effective stress, how the effect of seepage flow may be accounted for by means of driving body forces derived from the gradient of excess pore pressure distribution is shown. The stability analysis requires being able to evaluate the latter. In the specific situation of steady state flow and the rock deformability can be neglected, the seepage forces distribution is obtained as the solution to a classical hydraulic problem with appropriate boundary conditions.

As regards the rock strength properties, a modified Hoek-Brown failure criterion is adopted. From a practical viewpoint, the availability of closed-form expressions of the $\Pi$-functions for such a criterion makes it possible to derive rigorous upper bound solutions for the stability problem of any structure involving such a failure criterion: rock slope stability, tunnel wall or face stability, bearing capacity of foundations, etc.

The method has been applied to analyze the stability of a rock slope subjected to gravity and seepage forces. The preliminary step consisted in evaluating the seepage forces distribution using a finite element approach to solve the hydraulic problem. These forces thus computed have been therefore incorporated as components of the external loading in the subsequent stability analysis. The latter is based on the block rotational log-spiral failure mechanism. A parametric study have been performed with the objective to give insight on the influence of relevant geometrical, strength or load parameters. The predicted upper bound estimates of the rock slope safety factor have been compared to those derived from the simplified approach in which the effects of seepage flow are accounted for through a pore pressure coefficient. Additionally, the limitations inherent to the classical approach that is based on the concept of pore pressure have been highlighted. The applicability of the latter is actually reduced to stability analyses based on rigid-block failure mechanisms.

The natural extensions of this work should undoubtedly consist in evaluating the bearing capacity of foundations located at the vicinity of rock slopes by implementation of both rotational and punching-like failure mechanisms, derived from the generalization of Prandtl or multi-wedge failure mechanisms. It is expected from this implementation an effective improvement of the upper bound estimate of the bearing capacity. However, as explained previously, the seepage forces distribution should be evaluated by means of a numerical procedure in the case of Prandtl mechanism.

Provided that the concept of effective stress is valid, the framework formulated for handling the stability analysis of geotechnical structures in the presence of seepage forces is rather general. The possible extension of the approach to other kind of structure in which water flow should be accounted for in the safety analysis may be considered.

## Acknowledgments

The preparation of the paper had received financial support from Ecole Nationale des Ponts et Chaussées (ENPC-France) and the French Institute of Cooperation in Tunisia (French Embassy-Tunisia). The authors have greatly appreciated the financial support.

## References

[1] Leca E, Garnier J, Atwa M, Chambon P, Skiker A, Dormieux L, et al. Theoretical and physical modelling of the face stability of shallow tunnels. In: Proceedings of 14th international conference on soil mechanics and foundation engineering, vol. 3. Hamburg; 1997. p. 1424-4.

[2] Hoek E, Brown ET. Empirical strength criterion for rock masses. J Geotech Eng, ASCE 1980;106(GT9):1013-36.

[3] Hoek E. Strength of jointed rock masses. Geotechnique 1983;33:187-223.

[4] Marinos P, Hoek E. Estimating the geotechnical properties of heterogeneous rock masses such as flush. Bull Eng Geol Environ 2001;60:84-92.

[5] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Mining Sci 1997;34(8):1165-86.

[6] Jiang JC, Baker R, Yamagami T. The effect of strength envelope nonlinearity on slope stability computations. Can Geotech J 2003;40:308-25.

[7] Agar JG, Morgenstern NR, Scott J. Shear strength and stress-strain behaviour of Athabasca oil sand at elevated temperatures and pressure. Can Geotech J 1985;24(1):1-10.

[8] Goodman RE. Introduction to rock mechanics. 2nd ed. New York: Wiley; 1989.

[9] Baker R, Frydman S. Upper bound limit analysis of soil with nonlinear failure criterion. Soil Found 1983;23(4):34-42.

[10] Zhang XJ, Chen WF. Stability analysis of slopes with general nonlinear failure criterion. Int J Numer Anal Methods Geomech 1987;11(1):33-50.

[11] Drescher A, Christopoulos C. Limit analysis slope stability analysis with nonlinear yield condition. Int J Numer Anal Methods Geomech 1988;12(3):341-5.

[12] Serrano A, Olalla C. Ultimate bearing capacity of rock masses. Int J Rock Mech Mining Sci 1994;31(2):93-106.

[13] Serrano A, Olalla C. Allowable bearing capacity of rock foundations using a non-linear failure criterion. Int J Rock Mech Mining Sci 1996;33(4):327-45.

[14] Serrano A, Olalla C, Gonzalez J. Ultimate bearing capacity of rock masses based on the modified Hoek-Brown criterion. Int J Rock Mech Mining Sci 2000;37:1013-8.

[15] Yang X, Yin JH, Li L. Influence of a nonlinear failure criterion on the bearing capacity of a strip footing resting on rock mass using a lower bound approach. Can Geotech J 2003;40:702-7.

[16] Yang XL, Yin JH. Upper bound solution for ultimate bearing capacity with a modified Hoek-Brown failure criterion. Int J Rock Mech Mining Sci 2005;42:550-60.

[17] Merifield RS, Lyamin AV, Sloan SW. Limit analysis solutions for the bearing capacity of rock masses using the generalised Hoek-Brown criterion. Int J Rock Mech Mining Sci 2006;43:920-37.

[18] Saada Z, Maghous S, Garnier D. Bearing capacity of shallow foundations on rocks obeying a modified Hoek-Brown failure criterion. Comput Geotech 2008;35:144-54.

[19] Yang XL. Seismic bearing capacity of a strip footing on rock slopes. Can Geotech J 2009;46:943-54.

[20] Collins IF, Gunn CIM, Pender MJ, Wang Y. Slope stability analyses for materials with a nonlinear failure envelope. Int J Numer Anal Methods Geomech 1988;12(5):533-50.
[21] Yang XL, Li L, Yin JH. Stability analysis of rock slopes with a modified Hoek-Brown failure criterion. Int J Numer Anal Methods Geomech 2004;28:181-90.

[22] Saada Z, Maghous S, Garnier D. Seismic bearing capacity of shallow foundations near rock slopes using the generalized Hoek-Brown failure criterion. Int J Numer Anal Methods Geomech 2011;35(6):724-48.

[23] Miller TW, Hamilton JH. A new analysis procedure to explain a slope failure at the Martin Lake mine. Géotechnique 1989;39:107-23.

[24] Miller TW, Hamilton JH. Discussion on a new analysis procedure to explain a slope failure at the Martin Lake mine. Géotechnique 1990;40:145-7.

[25] Michalowski RL. Slope stability analysis: a kinematical approach. Géotechnique 1995;45:283-93.

[26] Bishop AW, Morgenstern NR. Stability coefficients for earth slopes. Géotechnique 1960;10:129-50.

[27] Yang XL, Zou JF. Stability factors for rock slopes subjected to pore water pressure based on the Hoek-Brown failure criterion. Int J Rock Mech Mining Sci 2006;43:1146-52.

[28] Drucker DC, Prager W. Soil mechanics and plastic analysis or limit design. Q Appl Math 1952;10(2):157-65.

[29] Salençon J. An introduction to the yield design theory and its applications to soil mechanics. Eur J Mech, A/Solids 1990;9(5):477-500.

[30] Maghous S, Bernaud D, Fréard J, Garnier D. Elastoplastic behavior of jointed rock masses as homogenized media and finite element analysis. Int J Rock Mech Mining Sci 2008;45(8):1273-86.

[31] de Buhan P, Fréard J, Garnier D, Maghous S. Failure properties of fractured rock masses as anisotropic homogenized media. J Eng Mech, ASCE 2002;128(8):869-75.

[32] Maghous S, de Buhan P, Bekaert A. Failure design of jointed rock structures by means of a homogenization approach. Mech Cohesive-Frictional Mater 1998;3(3):207-28.

[33] de Buhan P, Maghous S. Comportement élastique non linéaire d'un matériau comportant un réseau de joints. C R de l'Acad Sci 1997;324(IIb): 209-18.

[34] Hoek E, Carranza-Torres C, Corkum B. Hoek-Brown failure criterion—2002 edition. In: Proceedings of the North American rock mechanics society meeting, Toronto; 8-10 July 2002. p. 267-73.

[35] Hoek E. Estimating Mohr-Coulomb friction and cohesion values from the Hoek-Brown failure criterion. Int J Rock Mech Mining Sci Geomech Abstr 1990;27(3):227-9.

[36] de Buhan P, Cuvillier A, Dormieux L, Maghous S. Face stability of shallow circular tunnels driven under the water table: a numerical analysis. Int J Numer Anal Methods Geomech 1999;23:79-95.

[37] Yang XL, Yin JH. Slope stability analysis with nonlinear failure criterion. J Eng Mech, ASCE 2004;130:267-73.

[38] Salençon J. Calcul à la rupture et analise limite. Presses de l'ENPC. Paris; 1983.

[39] Kovarik JB. Qu'est-ce qu'un coefficient de sécurité en génie civil? Rev Fr Génie Civil 2000;4(6):607-51.