Stability Analysis of Earth Dam Slopes Based on the Second-Order Work Criterion in Finite Element Modeling

Abstract

Slope instability issues could cause severe damage and endanger the population, especially when dams are concerned. Over the past years, more and more refined numerical approaches have been proposed to evaluate the slope stability. However, a criterion based on the lack of numerical convergence is generally considered to compute the safety factor with this kind of approaches, which may result in a dramatic increase in the overall computation time for the probabilistic assessment of the safety factor (e.g., using Monte Carlo simulations that require the automation of a large number of simulations). This paper proposes an original approach coupling the shear strength reduction technique with the second-order work criterion. This approach is implemented in the open finite element software Cast3M, version 22.0. The relevance and efficiency of the developed approach is illustrated with two case studies: a heuristic slope and an existing earth dam. Safety factors similar to those calculated by classical approaches are obtained, but without pushing the calculation to the point of non-convergence. Among other advantages described in the paper, the proposed approach makes it possible to compute safety factors using a rational and physically based criterion, while facilitating probabilistic calculations.

1. Introduction

The stability of slopes has been intensively discussed in many contributions over the past decades. The likelihood of a slope to fail is usually accounted for through the definition of a safety factor. This safety factor is estimated from two main approaches: with traditional limit equilibrium methods (LEMs), or with numerical methods like the finite element method (FEM). With the substantial development of numerical tools and computing efficiency [1], the latter is becoming increasingly used in research works as an alternative to the LEM approach [2,3,4,5,6,7,8].

When applying the LEM approach, one has to keep in mind the shortcomings concerning assumptions on the inter-slice interactions and the failure surface that is specified beforehand. Other shortcomings are listed in [9,10]. Compared with LEMs, the FEM enables us to model the soil behavior in terms of the stress/strain relation (by different constitutive laws) including an elasto-plastic flow rule (that is usually non-associated as far as geomaterials are concerned) [3,9]. However, the LEM is still commonly used for engineering purposes because it remains the easiest and quickest way to assess slope stability with a satisfactory accuracy and few input data for many common applications.

In the review carried out by Duncan [11], attention was mostly focused on the use of numerical methods to assess the role of deformations on the slope stability [9]. This work introduced the basic framework for the computation of the slope stability assessment with elasto-plastic models implemented in numerical methods. The most common numerical method used in the literature to assess the slope stability remains the strength reduction technique coupled with the FEM. First described by Zienkiewicz [12], this technique consists of progressively decreasing the shear resistance of the soil until failure occurs. The strength reduction technique is extensively used in slope stability analyses using the FEM. These investigations concern principally relatively simple slopes [2,13], but some authors also applied this technique to analyze the stability of more complex geotechnical structures such as earth dams [9,14]. In many contributions, the safety factors obtained with the LEM and the strength reduction technique are very similar, except in some specific situations pointed out by Cheng et al. (2007) [3], i.e., when the cohesion of the soil is small and the friction angle is high, and vice versa. Cheng et al. (2007) also notice that most of the studies giving similar safety factors with the LEM and the FEM are limited to homogeneous soil slopes with a relatively regular geometry involving no special features (e.g., the presence of a thin layer of soft material or special geometry).

Although the FEM brings several advantages for a slope stability analysis, the lack of numerical convergence is frequently used as a failure criterion in the shear strength reduction method, which raises the issue of its physical relevance in predicting soil failure. Indeed, the notion of failure in geomechanics is clearly defined as a mechanical issue related to material properties [15,16], but it has yet to be integrated at the scale of a structure in order to control the strength reduction routine in a slope stability analysis by the FEM (cf. Section 2.1). In other words, there is a need for the definition of an explicit criterion at the scale of the overall structure (and not only at the material point scale) to attest the effective failure of the slope or the structure. Even though the material failure criterion is well-known, failure at the engineering scale depends on other factors, such as the type of boundary conditions and the specific geometry of the structure. As a result, there is no unique criterion for a global failure [10] and several failure criteria can be proposed [9,17]. The failure criteria used in a slope stability analysis by the FEM are based either on purely numerical considerations (e.g., divergence of the computation), or on further assumptions that bring numerical methods closer to the LEM (cf. Section 2.2).

Classically, failure in soils is driven by a plastic failure criterion at the material point scale [10]. This criterion is represented by a single plastic limit failure envelope in the stress space (the Mohr–Coulomb criterion is widely selected), and all failure states are assumed to belong to the stress points constituting this envelope [15,18]. Nevertheless, it has been widely observed in practice that failure can occur before this failure surface is reached, (e.g., for mobilized friction angles smaller than the friction angle if the material follows a Mohr–Coulomb plastic criterion). An illustrative example is the liquefaction of loose soils in undrained triaxial tests [19], which stems from the fact that geomaterials are known to be non-associated (the plastic flow rule is not normal to the yield surface) [18]. In other words, volume changes (dilatancy or contractancy) can occur under pure shear loading. This motivated the definition of a more general stability criterion at the material point scale. Following the pioneering work of Hill [20], one rational candidate is the second-order work criterion. Within the framework of small strain continuum mechanics, this criterion can be formulated as follows: “for a given equilibrium (σ, ε) reached after a given loading history, if there exists at least one stress increment dσ associated with a strain response dε such that W₂ = dσ∶dε < 0, the material point is unstable”. The physical meaning behind this corresponds to a situation in which the deformation of the mechanical system can be pursued without any input of energy from the observer (on the contrary, the system is spontaneously releasing some stored energy). Over the past years, this criterion was proved to be powerful, when dedicated to the characterization of failure in geomaterials from both numerical [16,21,22] and experimental considerations [15,18,23]. By integrating the second-order work over the entire structure domain, the vanishing of this structural second-order corresponds to the moment when the entire system switches from quasi-static to dynamic conditions (prior to the onset of any potential discontinuity in the structure integrity). This structural criterion was employed to investigate slope stability issues, mainly in the context of heuristic slopes without hydraulic conditions [10,24]. Moreover, Hu et al. (2020) [25] show that this criterion is independent of the mesh size (on the example of a triaxial compression test modelled using the finite difference method).

The present manuscript proposes a combination of the strength reduction technique and the second-order work criterion in order to have a rational and physically based criterion to detect the structure failure and, thus, define a related safety factor for earth dam stability analysis. Taken alone, the strength reduction technique and the second-order work are not original concepts as they are widely used in the literature. The novelty here is to combine these two aspects in order to give a robust criterion able to compute a safety factor without relying on purely numerical aspects (non-convergence of the computations). In this work, the second-order work criterion is integrated at the structure scale and used to anticipate the effective failure of the structure. This allows us to stop the strength reduction procedure and save computation time. The open-access FEM code Cast3M, working as a command language, is used to code the proposed approach. The methodology is then tested on two different types of embankments: a heuristic slope without considering seepage and an existing earth dam in a normal operating situation. These two examples allow us to show that the proposed approach can be used to obtain relevant safety factors for cases of varying complexity.

2. Materials and Methods

2.1. Shear Strength Reduction Technique

As discussed in introduction, slope stability analyses are traditionally performed using LEM. The principal methods using limit-equilibrium are discussed in Duncan’s review [11]. In the LEM approach, the safety factor is given by an analytical expression as follows [10,11]:

$$\mathrm{FoS} =\frac{{\mathsf{\tau }}_{\max }}{{\mathsf{\tau }}_{\mathrm{eq}}}=\frac{\mathrm{maximum}\mathrm{shear}\mathrm{strength}\mathrm{of}\mathrm{the}\mathrm{soil}}{\mathrm{shear}\mathrm{strength}\mathrm{required}\mathrm{for}\mathrm{equilibrium}}$$

(1)

${\mathsf{\tau }}_{\max }$ is the maximal resulting force of the shear efforts along the potential slip surface. A slope is considered stable as soon as the safety factor is greater than one, and, inversely, a slope will fail if the safety factor is beneath the unit. The safety factor can be regarded as the factor by which the shear strength has to be reduced to bring the slope to the failure state [10,11].

The numerical approach, using FE (or finite difference) method, makes it possible to include the elasto-plastic behavior of the soil material by implementing the stress–strain relation at the material point scale.

Basically, three classes of approaches to assess slope stability using FEM can be distinguished [2]: (i) the enhanced LEM where the stress field is stemming from an FEM computation; (ii) the gravity-increasing approach where the gravity-driven loading is increased until the failure occurs; and (iii) the strength reduction technique which has been introduced by Zienkiewicz [12], and which remains the most popular technique for FE slope stability analyses. The latter approach will be considered thereafter in the manuscript.

The general principle of the strength reduction technique is to decrease gradually the shear strength properties until stability limit is reached. Considering a Mohr–Coulomb yield surface, the initial shear strength parameters $\mathrm{c}\prime$ (cohesion) and $\tan {\mathsf{\phi }}^{\prime }$ (friction coefficient) are reduced by a factor ${\mathrm{F}}_{\mathrm{red}}$ to give the parameters ${\mathrm{c}\prime }_{\mathrm{red}}$ and $\tan {\mathsf{\phi }\prime }_{\mathrm{red}}$ used in the non-linear calculations as follows:

$$\begin{gathered} {\mathrm{c}\prime }_{\mathrm{red}}=\raisebox{1ex}{$\mathrm{c}\prime $}\!\left/ \!\raisebox{-1ex}{${{\mathrm{F}}_{\mathrm{red}}}_{}$}\right. \\ {\mathsf{\phi }\prime }_{\mathrm{red}}={\tan }^{-1}\left(\frac{\tan {\mathsf{\phi }}^{\prime }}{{{\mathrm{F}}_{\mathrm{red}}}_{}}\right) \end{gathered}$$

(2)

The safety factor ${\mathrm{F}}_{\mathrm{s}}$ is assumed to be equal to the reduction factor ${\mathrm{F}}_{\mathrm{red}}$ at failure. This leads to a similar definition of the safety factor as in LEM, as the safety factor is the number by which the original shear strength parameters must be divided in order to bring the slope to the point of failure [9,14].

The principal shortcoming of the strength reduction technique is that this technique requires an iterative calculation scheme, which is time-consuming. This calculation scheme also has to include a failure criterion, indicating when the slope failure has occurred. The classical criteria are discussed thereafter in Section 2.2. In addition, in some cases, some mesh dependency issues may be observed. Finally, Cheng et al. (2007) demonstrate that the strength reduction technique is sensitive to non-linear solutions algorithms (i.e., the numerical solver) in some specific cases [3]. They also illustrate that the elastic modulus $\mathrm{E}\prime$ has no influence on the strength reduction technique as the computed safety factor is not sensitive to changes in this parameters, confirming that only the total unit weight $\mathsf{\gamma }$, the shear strength parameters $\mathrm{c}\prime$ and $\mathsf{\phi }\prime$, and the geometry of the problem are important in a FE slope stability analysis [9].

2.2. Limit-State Criteria in Strength Reduction Technique

According to [9,26], several failure criteria are available but they are not necessarily equivalent which may result in different safety factor values. Three main definitions are as follows [19]: (i) bulging of the slope line; (ii) limit shear; and (iii) non-convergence of the solution.

The first criterion is based on the analysis of the horizontal displacement of a specific point (referred to as the bulging point) close to the failure surface (e.g., the nodal point immediately above the slope toe). Failure is thought to have occurred when the displacement of this point exceeds a maximum limit. This agrees with the observation reported in [9]. When examining the evolution of the maximal displacement of the bulging point (or a dimensionless expression of this displacement, $\mathrm{E}\prime {\mathrm{d}}_{\max }/{\mathsf{\gamma }\mathrm{H}}^2$, where ${\mathrm{d}}_{\max }$ is the maximum nodal displacement at convergence and $\mathrm{H}$ is the height of the slope) versus the reduction factor ${\mathrm{F}}_{\mathrm{red}}$, the authors pointed out that there is a sudden acceleration of maximal displacement close to the LEM value of safety factor (see Figure 1 based on the results from [9]). This method introduces two sources of arbitrariness in the choice of the bulging point and on the maximum displacement allowed before failure.

The second criterion consists of using the computed FEM stresses along the same critical slip surface from which the LEM solution is calculated. This critical slip surface is determined by LEM and corresponds to the surface where the safety factor is the smallest upon a set of investigated surfaces. The safety factor corresponds to the ratio between the average mobilized shear strength along the slip surface and the prescribed value. The main advantage here is that the safety factor is computed with more realistic stresses obtained with FEM but the method cannot help identify a failure surface of a different shape.

The third and most popular criterion is the non-convergence of the numerical computation. This criterion makes sense insofar as the force equilibrium equations are solved by the non-linear iterative procedure. When the algorithm is unable to converge within a given number of iterations, no stress distribution can be found to satisfy both the failure criterion and global equilibrium: the failure is said to have occurred [9]. The safety factor is obtained in progressively reducing the shear strength parameters in a strength reduction technique procedure until non-convergence of the algorithm occurs. As shown by Nicot et al. [18], the failure corresponds to a bifurcation from a quasi-static regime to a dynamical one. As the governing equations do not take inertial mechanisms into account, the problem becomes ill-posed when failure has occurred, which also corresponds to a loss of uniqueness.

These three criteria rely on the use of arbitrary thresholds that lack a physical basis. The non-convergence is broadly used, certainly because it can be conveniently implemented within most existing FE codes. However, this failure criterion is sensitive to both the number of iterations and the algorithm of the FE code. Moreover, the non-convergence will generally cause fatal errors in FE codes, which can be a true issue for the user, for example, in the case where parametric studies are launched for probabilistic analyses [27]. The results of the strength reduction technique are, thus, dependent on the interpretations of the user, including, to some extent, his personal expertise in geotechnical and numerical modeling.

As an alternative to stopping the strength reduction procedure in earth dam stability analysis, this manuscript investigates the relevance of the second-order work criterion. This criterion has already been used in slope stability analyses by several authors [19,28]. The novelty here is to couple the second-order work criterion and the strength reduction technique.

2.3. Second-Order Work Instability Criterion

2.3.1. Introduction

Within the framework of elasto-plasticity, a first class of material failure can be defined through the use of a plastic limit surface. Any stress state outside this envelope cannot be applied to the considered material in the static regime. This definition is, however, not general enough, especially for geomaterials, because of the non-associated character of the flow rule (the yield surface and the plastic potential do not coincide, which means that the dilatancy angle does not match the mobilized friction angle for a material following Mohr–Coulomb criterion) [22]. Indeed, some experimental observations have shown that instability can be encountered strictly inside the failure surface [18]. One classical example of such contradiction is the static liquefaction of geomaterials under undrained triaxial conditions [18,23].

Therefore, there is a need for a more general criterion to define soil instability. Hill [20] proposed such a criterion based on the second-order work, which is defined as the scalar product between the stress increment $\underset{\_}{\mathsf{\delta }\mathsf{\sigma }}$ and the associated strain increment $\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}=\underset{\_}{\underset{\_}{\mathrm{M}}}:\underset{\_}{\mathsf{\delta }\mathsf{\sigma }}$, where $\underset{\_}{\underset{\_}{\mathrm{M}}}$ is the constitutive operator (non-symmetric for non-associated materials):

$${\mathrm{w}}_2=\underset{\_}{\mathsf{\delta }\mathsf{\sigma }}:\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}$$

(3)

According to Hill’s approach, a material system is considered unstable if, for a small stress perturbation, the deformation can develop without requiring external energy from the system [20,22]. On the other hand, the Hill criterion states that a stress–strain state is stable if, for any couple $(\underset{\_}{\mathsf{\delta }\mathsf{\sigma }},\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }})$, the second-order work ${\mathrm{w}}_2$ is strictly positive:

$$\forall \left(\underset{\_}{\mathsf{\delta }\mathsf{\sigma }},\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}\right)\in {\mathrm{R}}^{2\mathrm{n}}\backslash \left\{0\right\}{, \mathrm{w}}_2=\underset{\_}{\mathsf{\delta }\mathsf{\sigma }}:\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}>0$$

(4)

where $\mathrm{n}$ is the dimensions of the stress space.

When dealing with the hydro-mechanical FE modeling of the initiation of failure in soils, Kakogiannou et al. present three different expressions of the second-order work according to the definition of the stress considered [22]: one expressed in terms of effective stress, one expressed with total stress, and, finally, the third one taking into account the hydraulic contribution in a partially saturated medium. A simulation of rapid desaturation of a soil sample demonstrates that all of these three expressions give very close numerical results and are all able to capture the instability mechanism induced by cavitation of the pore water [22]. It can be underlined that the second-order work expressed with the effective stress usually vanishes a little earlier than the other two formulations.

In order to implement the second-order work criterion, Khoa [19] proposed a normalized form of the second-order work, as follows:

$${\mathrm{w}}_{2\mathrm{n}}=\frac{\underset{\_}{\mathsf{\delta }{\mathsf{\sigma }}^{\prime } }: \underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}}{‖\underset{\_}{\mathsf{\delta }{\mathsf{\sigma }}^{\prime }}‖‖\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}‖}\in \left[-1, 1\right]$$

(5)

In a FE model, the calculation of this normalized term is performed on the integration points of the continuous medium between two successive equilibrium states. This is, therefore, referred to as a local (or material) second-order work. By integrating the quantity ${\mathrm{w}}_2$ over the whole domain $\mathrm{V}$, the global (or structural) second-order work ${\mathrm{W}}_{2\mathrm{n}}$ can be determined:

$${\mathrm{W}}_{2\mathrm{n}}=\frac{{\int }_{\mathrm{V}}^{}\underset{\_}{\mathsf{\delta }{\mathsf{\sigma }}^{\prime }}:\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }} \mathrm{dV}}{{\int }_{\mathrm{V}}^{}‖\underset{\_}{\mathsf{\delta }{\mathsf{\sigma }}^{\prime }}‖‖\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}‖ \mathrm{dV}}$$

(6)

In this manuscript, Hill’s criterion is used as an alternative way to simulate failure occurrence. It consists of analyzing the sign of both local (${\mathrm{w}}_{2\mathrm{n}}$) and global (${\mathrm{W}}_{2\mathrm{n}}$) second-order works to detect the apparition of instability on both local and global scales.

At the structure scale, the vanishing of the second-order work corresponds to the moment when the entire system switches from quasi-static to dynamic conditions (with the potential creation of discontinuities in the structure). It should be acknowledged that we could imagine situations where the second-order work vanishes locally while it remains strictly positive at the structure scale. However, even in this case, the onset of discontinuities in the structure takes time to develop and the global second-order work usually vanishes prior to it.

2.3.2. Implementation Procedure of the Second-Order Work Criterion for Slope Stability in a FE Code

This section focuses on the technical implementation of Hill’s criterion in an FEM code to assess the stability of slopes from the strength reduction technique. The global implementation procedure is schematically described in Figure 2. It consists of three main parts: the construction of the FE model, the strength reduction technique loop, and the second-order work computation.

In this article, the FEM code Cast3M is used [29]. (Unlike others FEM code, Cast3M does not work like a ‘black box’. Cast3M has a command language consisting of a series of operators that allow the user to manipulate data and results in the form of objects, giving them names with a specific language, which is both a programming and a modeling language [29]. A complete code has been written for this work using Cast3M tools).

A detailed description of the different steps of the implementation procedure is given as follows:

Step 1—Initialization (pre-processing) of the FE model: The model geometry and the mesh are created and the boundary conditions are set. For the material properties, an elasto-plastic constitutive behavior following Mohr–Coulomb yield surface, with no hardening/softening and a non-associated flow rule (ψ ≠ φ ≠ 0), is selected and the initial parameters are set.

Step 2—FEM processing n°0 (Initial equilibrium state): An initial equilibrium state is reached by FEM processing after the incremental application of the different loads up to their nominal values. This step constitutes a classical FE routine computation using an elasto-plastic model based on Mohr–Coulomb criteria. The obtained initial equilibrium state is the starting point for the strength reduction technique procedure. Strain, stress, and displacement fields are then extracted for further analysis.

Step 3—Shear parameters reduction n°i: The effective shear strength parameters $\mathrm{c}\prime$ and $\mathsf{\phi }\prime$ are reduced by a factor ${\mathrm{F}}_{\mathrm{red}}$ using Equation (2). This reduction in shear parameters is carried out gradually (loop iterations) by increasing the factor ${\mathrm{F}}_{\mathrm{red}}$ by an increment $∆{\mathrm{F}}_{\mathrm{red}}$: ${\mathrm{F}}_{\mathrm{red},\mathrm{i}}={\mathrm{F}}_{\mathrm{red},\mathrm{i}-1}+∆{\mathrm{F}}_{\mathrm{red}}$. The value of the increment $∆{\mathrm{F}}_{\mathrm{red}}$ is generally chosen between 0.01 and 0.05, possibly adopting a higher value for the first iterations i (when failure is unlikely to occur) and a lower value for the last iterations. For reduction n°1, if $∆{\mathrm{F}}_{\mathrm{red}}=0.05$, the factor ${\mathrm{F}}_{\mathrm{red},\mathrm{i}}=1.05$.

Step 4—FEM processing n°i: The reduced shear strength parameters ${\mathrm{c}\prime }_{\mathrm{red}}$ and ${\mathsf{\phi }\prime }_{\mathrm{red}}$ are then used to compute a new equilibrium state (n°i) from the previous one. Strain and stress fields, and maximal displacement are extracted. To reduce the number of FEM computations, the reduction increment $∆{\mathrm{F}}_{\mathrm{red}}$ can be modified as the loop progresses as a function of the growth rate of the “maximal displacement vs. reduction percentage” curve. A sharp increase in the slope of the “maximal displacement vs. reduction percentage” curve indicates that the failure is approaching, which suggests that a new calculation loop has to be performed with a lower reduction increment $∆{\mathrm{F}}_{\mathrm{red}}$.

Step 5—Increments calculation: The stress and strain increments $\underset{\_}{\mathsf{\delta }\mathsf{\sigma }}$ and $\underset{\_}{\mathsf{\delta }\mathsf{\epsilon }}$ are calculated from stress and strain fields in states $\mathrm{i}$ and $\mathrm{i}-1$.

Step 6—Second-order work calculation: Local and global normalized second-order works are computed from these increments using Equations (5) and (6).

Step 7—Second-order work stop condition: The sign of the global second-order work ${\mathrm{W}}_{2\mathrm{n}}$ is analyzed. If the global second-order work is positive (i.e., “stable”), the strength reduction loop is pursued: the index i of the loop is incremented ($\mathrm{i}+1$) and the procedure restarts at Step 3. If the global second-order work ${\mathrm{W}}_{2\mathrm{n}}$ is negative (i.e., “unstable”), the strength reduction loop is stopped and the procedure continues at Step 8. However, a very small negative value ζ (e.g., absolute value lower than ζ = 10⁻⁴) should not stop the strength reduction technique loop because the global second-order work ${\mathrm{W}}_{2\mathrm{n}}$ could fluctuate slightly around the nil value during the computation due to numerical noise.

Step 8—Safety Factor: The safety factor SF is equal to the reduction factor ${\mathrm{F}}_{\mathrm{red}}$ at the end of iterations of the strength reduction technique loop.

3. Results

The applicability of the methodology exposed in the previous section is illustrated by two examples: a simple homogeneous slope without seepage and a real dam where seepage and soil properties stemming from field data are considered.

3.1. Heuristic Case: Application for a Homogeneous Slope

The first example corresponds to a heuristic slope in cohesive and frictional soil with a height of 10 m and a 1:1 slope inclination. The geometry considered also includes a foundation layer 5 m thick, as shown in Figure 3. This slope has already been analyzed in several studies [4,5,6]. The objective here is to validate the proposed approach with a simple case for which the standard approaches is known to be reliable. A homogeneous material for both the foundation and slope is considered, with values of $\mathrm{c}\prime$, $\mathsf{\phi }\prime$, and $\mathsf{\gamma }$, respectively, equal to 10 kPa, 30°, and 20 kN/m³. These values are the same as those used in the previous studies, in which no indications were reported concerning the elastic parameters (Young modulus $\mathrm{E}\prime$ and Poisson’s ratio $\nu \prime$) and the plastic potential (characterized by the dilatancy angle $\mathsf{\psi }\prime$). Standard values have been selected for these parameters, which are believed to not have much influence on the safety factor value [9]. Note that $\mathsf{\psi }\prime$ is chosen to be different from $\mathsf{\phi }\prime$ which corresponds to the non-associated plasticity. No seepage is supposed to occur in the soil, so that the total stresses are considered in the computation. The values of the different parameters are reported in

Table 1. Soil properties in example 1.

Properties	Cohesion–Frictional Soil
Dry unit weight ${\mathsf{\gamma }}_{\mathrm{d}}$ (kN·m−3)	20.0
Cohesion $\mathrm{c}\prime$ (kPa)	10.0
Friction angle $\mathsf{\phi }\prime$ (°)	30.0
Dilatancy angle $\mathsf{\psi }\prime$ (°)	10.0
Elastic modulus $\mathrm{E}\prime$ (MPa)	50.0
Poisson ratio $\nu \prime$	0.33

.

The stability of this slope has been preliminary analyzed using the LEM with Bishop’s simplified method. This approach provides a safety factor equal to 1.21, corresponding to those computed with similar methods (i.e., 1.20, 1.21, and 1.21) by Cho (2010) [4], Li et al. (2015) [5], and Liu et al. (2017) [6], respectively. The obtained failure surface is represented by the red dashed line in Figure 3.

To apply the methodology developed in the last section, the homogeneous slope is modeled using the FE code Cast3M [29]. The mesh of the slope is similar to that used in previous studies and is composed of four-node quadrilateral elements. These elements are triangular, close to the slope surface (see Figure 3). The FE model includes a non-associated elasto-plastic behavior for the soil. In this case study, the Drucker–Prager criterion is selected as an approximation of the Mohr–Coulomb criterion. It smoothes out the singularities of the hexagonal pyramidal yield surface of the Mohr–Coulomb criterion. With the plane strain hypothesis, relationships exist to link the Drucker–Prager parameters with the shear strength parameters $\mathrm{c}\prime$ and $\tan {\mathsf{\phi }}^{\prime }$ [30].

Concerning the boundary conditions, the bottom of the mesh is considered fixed and only vertical displacement is allowed on the left boundary. A gravity load is increasingly applied up to its nominal value to obtain the initial equilibrium state. The strength reduction technique routine described in Section 2 is then performed until a reduction factor ${\mathrm{F}}_{\mathrm{red}}$ corresponding to the first occurrence of a negative value of the global second-order work is reached. For the homogeneous slope, the parameter reduction procedure results in a reduction factor ${\mathrm{F}}_{\mathrm{red}}=1.23$, corresponding to a reduction of 18.79% in the parameter values. Figure 4 gives an overview of the results obtained in this example.

The convergence of the algorithm can be visually followed with the curve of Figure 5 showing the normalized second-order work’s progression. The notations A, B, C, and D refer to Figure 4. The lack of convergence of the computation occurs for a reduction of 21.75% in the parameter values, corresponding to a safety factor SF = 1.28 (materialized by the green vertical line on Figure 5).

The four graphs Figure 4a–d represent the negative local second-order work field at different reduction levels (3%, 9%, 15%, and 18.79%, respectively). On these graphs, the black dashed line represents the LEM failure surface, corresponding to a safety factor of 1.21 (i.e., a 17.35% reduction).

The global second-order work ${\mathrm{W}}_{2\mathrm{n}}$ (red curve in Figure 5) is fluctuating very slightly around the nil value with absolute values ranging within [10⁻⁹–10⁻⁵] until the 15th iteration. The initial percentage of the reduction step is equal to 1.5% and this step decreases progressively and continuously when approaching the failure (as seen on Figure 5). Then, for a reduction of 18.79% at the 16th iteration, ${\mathrm{W}}_{2\mathrm{n}}$ drops to a significantly negative value, inferior to 10⁻⁴, as explained in Section 2.3.2. This stops the numerical strength reduction technique loop and provides a safety factor equal to 1.23, close to the value obtained with the LEM. The LEM can be assumed to give an accurate safety factor for this classic example, considering a Mohr–Coulomb behavior and no pore water pressures. This validates the proposed approach.

Concerning the local second-order fields, it can be seen that very low negative values are visible after a 3% reduction of the shear strength parameters (Figure 4a). These negative values are clearly observed for larger reduction levels of 9% and 15% (Figure 4b,c, respectively). It can be noticed that these negative values are localized along the failure surface obtained with the LEM. For the last reduction step (Figure 4d), the negative values become significantly larger and make it possible to visually identify the failure surface without making a priori assumptions about it (which is almost identical to that obtained with the LEM).

3.2. Application for an Existing Earth Dam

The previous example was investigated to illustrate the effectiveness of the second-order work criterion combined with the strength reduction technique in the case of a heuristic slope issue in order to validate the proposed approach. The example considered in this section explores the real case of an embankment dam, with a more complex geotechnical situation including a variety of materials and the occurrence of a seepage mechanism.

The dam considered is a 23 m-high pseudo-zoned dam located in western France. The dam body includes three zones: a core zone (COR) composed of sandy silts, and two shoulders (UPS: upstream shoulder and DOS: downstream shoulder) made up of coarse sands. The downstream shoulder is composed of a material slightly coarser than that of the upstream shoulder. All materials composing the dam stem from the degradation of claystones that are present in the foundation. The geometry of the main cross-section of the structure is presented in Figure 6a.

The FE model, also developed with Cast3M [29], considers the different zones of the dam, as well as its drainage system. The modeled geometry is visible in Figure 6b. The seepage mechanism is accounted for in the embankment. A hydraulic FE seepage analysis is performed beforehand of the mechanical calculation involving the strength reduction technique. For the hydraulic calculation, the upstream boundary condition corresponds to the normal operating water level (50.0 m). Zero pressure is considered in the drainage system as the downstream boundary.

The computed pore water pressures act as a part of the loading considered in the mechanical computation, which also includes the hydrostatic pressure of the reservoir at the normal water level and the weight of the structure. The bottom of the mesh is considered fixed and only the vertical displacements are allowed on the edges. As in the previous heuristic case, the Drucker–Prager criterion is implemented within the elasto-plastic constitutive model of the material in hand. The parameters considered in this example are obtained from available in-field (dam) data. The values of these parameters are reported in

Table 2. Soil properties in example 2.

Properties	Coarse Sands (UPS and DOS)	Sandy Silts (COR)
Dry unit weight ${\mathsf{\gamma }}_{\mathrm{d}}$ (kN·m−3)	19.8	17.9
Cohesion $\mathrm{c}\prime$ (kPa)	8.9	13.4
Friction angle $\mathsf{\phi }\prime$ (°)	34.8	34.1
Dilatancy angle $\mathsf{\psi }\prime$ (°)	12.0	12.0
Elastic modulus $\mathrm{E}\prime$ (MPa)	45.0	40.0
Poisson ratio $\nu \prime$	0.33	0.33

.

The loads (gravity, pore water pressures, and hydrostatic pressure) are increasingly applied through FEM processing, leading to an initial equilibrium state, which is the starting point of the strength reduction technique routine. The shear strength parameters of the different materials are then reduced until the global second-order work becomes negative (according to the iterative procedure described in Section 2). The results are shown in Figure 7. The four graphs Figure 7a–d represent the negative local second-order work field at different reduction levels, 10%, 30%, 50%, and 60.03%, respectively. On these graphs, the white dashed line represents the LEM failure surface, corresponding to a safety factor of 2.70 (i.e., a reduction of 62.3%).

In Figure 8, the curve of the normalized second-order work’s progression shows the convergence of the algorithm. The notations A, B, C, and D refer to Figure 7. The lack of convergence of the computation occurs for a reduction of 61.25% in the parameter values, corresponding to a safety factor SF = 2.58 (materialized by the green vertical line in Figure 8 and Figure 9). As in the previous example, the value of ${\mathrm{W}}_{2\mathrm{n}}$ (red curve in Figure 8 and Figure 9) fluctuates around zero with very small absolute values until the 27th iteration. Then, a drop is clearly observed at the 28th iteration for a reduction of 60.03%. This corresponds to a safety factor of 2.50.

Figure 9 compares the evolution of the second-order work and the maximum displacement as a function of the parameter reduction. The initial reduction step used in this case was equal to 2.5% and the reduction step is reduced since the 21st iteration (reduction of 52.5%). With a traditional approach, the evolution of the maximal displacement (grey curve on Figure 9) displays a classical trend, with a sharp increase in the displacement as the failure approaches. It can be seen that the asymptote of this curve is the line corresponding to the non-convergence (reduction of 61.25%).

Overall, the safety factor obtained with the second-order work approach (2.50) is lower than the one obtained with the non-convergence criterion (2.58) that is also lower than the safety factor obtained from the LEM (2.70). For the present case, the LEM assumptions and the non-physical failure criterion lead to an overestimation of the safety factor, which is not a conservative approach.

Table 3. Comparison of the results of the different approaches.

	LEM	FEM with Non-Convergence	FEM with Second-Order Work
Safety factor	2.70	2.58	2.50
Failure surface	White dashed line in Figure 7 and Figure 10	Maximal shear strain field on Figure 10	Visible on Figure 7d
Computational efficiency	≈2 s	≈3 min	≈3 min

compares the different safety factors obtained with the three different approaches: the LEM, the FEM with the non-convergence criterion, and the FEM with the second-order work criterion.

Concerning the failure surface, it is worth noting that the evolution of the local second-order field is different from that observed in the case of the heuristic slope. The Figure 7a shows that a zone in which the local second-order work is negative develops early. This zone is located at the interface between the core and the upstream shoulder of the dam. As a result, the plastic strains are the more important in this zone of the dam at the beginning of the strength reduction loop. The plastic strains, and the negative local second-order work field as well, evolve as the shear parameters are reduced. Figure 7b,c illustrate the evolution of the negative local second-order work from the zone at the interface between the core and the upstream shoulder of the structure. Another zone of negative local second-order work develops during the reduction process at the downstream toe of the dam. As the factor of reduction increases, the two zones connect (see Figure 7c). Then, failure occurs and the ultimate field of the negative values of the local second-order work (see Figure 7d) clearly shows the development of a shear band representative of the failure surface.

It is worth noting that this band is coherent with the maximal shear strain (MSS) band obtained in the computational step before the lack of numerical convergence (Figure 10). On Figure 10, it could be seen that the MSS band is not continuous. The area beneath the drainage system has few deviatoric deformations but it is here where the local second-order work is clearly negative, as seen on Figure 7d.

Table 3 gives a comparison between the location of the failure surface and the computational efficiency of these approaches. It could be seen in Figure 7 and Figure 10 that the unstable volume is quite similar for the different approaches. However, the computational efficiency is lower in the case of the LEM than for the FEM approaches.

Thereby, the failure surface obtained from the FEM computation (local second-order work field in Figure 7) is deeper than the LEM critical circle (that remains a good approximation of it). The developed approach, combining FEM calculations with the second-order work criterion, is able to assess the slope stability of complex soil structures in determining a safety factor without making any assumption on the shape and location of the failure surface.

4. Discussion

This section discusses the advantages and limits of the proposed approach, as well as possible improvements.

In this work, an original use of the combination between the strength reduction technique and the second-order work criterion has been used for assessing the stability analysis of slopes and earth dams. It highlights some features concerning the computation of the safety factor from FEM simulations.

Although the FEM approach enables us to avoid LEM assumptions concerning the shape and position of the failure surface and allows us to consider a more realistic behavior of the soil using a constitutive behavior, this approach needs more input parameters to be accurate. The plentiness of these parameters and the uncertainty on the exact values of each parameter certainly introduce a bias in the calculation of the safety factor. In addition, these constitutive parameters need to be calibrated, often from limited field data.

The study cases considered in this paper correspond to a simple case of a homogeneous slope and a more complex case corresponding to a zoned earth dam subject to internal flows. The sliding surfaces are relatively deep, which is favourable for identifying a significant reduction in the second-order work. Future studies could test the approach to detect or discriminate between superficial slip surfaces that lead to non-convergence with the classic parameter reduction method but that do not constitute a safety issue for the overall stability of the slope.

The case studies considered involve different types of soil. Further studies could also consider the case of slopes containing materials with different levels of rigidity, such as slopes with masonry walls or reinforcements with anchors.

The examples considered in this paper show that, during the strength reduction loop, the global second-order work could remain almost equal to zero before dropping to significantly negative values. This allows a clear identification of the occurrence of failure but it could also cause computational issues and it needs the arbitrary definition of the parameter ζ (see Section 2.3.2).

In the same way, improvements can be made in terms of implementing the proposed approach in an FE code. For example, an optimization can be sought for the reduction steps considered in the strength reduction technique routine, by adopting higher values for the first iterations (for a faster result) and adopting lower values for the last iterations (for a better safety factor accuracy).

In the developed approach, the second-order work is used as a criterion to stop the strength reduction method. A practical advantage is that the loop for reducing strength parameters is stopped before a non-convergence occurs. This can be crucial for analyses where non-convergence is not desirable, for example, in parametric studies related to sensitivity analyses or probabilistic analyses using Monte-Carlo simulations that require the automation of a large number of simulations.

Finally, it is well-known that the LEM and FEM approaches give similar results in terms of the safety factor for simple cases. In the case of a complex dam geometry or soil behavior, the developed approach takes on its full meaning. This can be observed in the case of the real dam studied in this article: the failure surface and the safety factor are not exactly the same as those obtained with the LEM. The prospect of applying the developed approach to more complex cases, with the use of more sophisticated constitutive behaviors of soils, can be the purpose of a future work.

5. Conclusions

This paper proposes to use a criterion based on the second-order work to stop the strength reduction procedure. Although these two concepts are well-known in the literature, the novelty stems from their combination to compute a safety factor to assess the slope stability. This approach offers several advantages:

• The failure of the geomaterial is described with a criterion based on a physical consideration rather than on the non-convergence of the numerical computation.
• Obtaining the safety factor before the occurrence of the non-convergence allows us to overcome the numerical crash of the FEM code. This can be useful for parametric studies related to sensitivity analyses or probabilistic analyses that require the automation of a large number of simulations.
• The failure surface emerges naturally without any a priori shape assumptions by computing the local second-order work.

This primary work brings two main prospects: First, the proposed approach could be applied to more complex slope stability analyses, involving unusual geometries or complex soil behaviors. Secondly, the characterization of mechanical stability using the second-order work criterion is versatile and can pave the way for the definition of new concepts for stability analysis such as the definition of a bifurcation domain at the structure scale [31].