Stability of Riverbanks Under Unsteady Flow Conditions - J-Stage

International Journal of Erosion Control Engineering Vol. 7, No. 2, 2014

Original Article

Stability of Riverbanks Under Unsteady Flow Conditions Ahmed ALY EL-DIEN1, Hiroshi TAKEBAYASHI2, and Masaharu FUJITA2 1 Dept. of Civil and Earth Resources Engineering, Graduate School of Engineering, Kyoto University (Kyoto daigaku-katsura, Nishikyo-ku, Kyoto 615-8530, Japan) E-mail: [email protected] 2 Disaster Prevention Research Institute, Kyoto University, Ujigawa Open Laboratory (Higashinoguchi, Shimomisu, Yokoouji Fushimi, Kyoto, 612-8235, Japan)

Stability of riverbanks under three unsteady river flow conditions is simulated. First is the case of constant water level in the river, second is the case of filling the river to the riverbank top and finally the case of drawdown. Two modes of drawdown are analyzed; the rapid drawdown and the slow drawdown. A finite element model of saturated unsaturated seepage flow was de-coupled with a plain strain elastic-plastic finite element model using strength reduction technique to calculate the stability factor of safety (FOS) of bank material. The influences of location of phreatic surface, pore-water pressure, drawdown rate and ratio on the riverbank stability are discussed in details for each case. Results showed that safety factor for saturated riverbank is nearly 83% of its dry value (for the studied case). The case of filling is most likely to stabilize the riverbank. For the case of slow drawdown, the minimum FOS occurs when the water depth is about 0.25~0.3 the bank height. The case of rapid drawdown is the most critical case. Key words: Riverbank stability, drawdown, seepage, Strength reduction method, Finite element

Montaro in Peru [Duncan et al., 1990], the Three Gorges Reservoir slope [Zheng et al., 2005], the Pilarcitos Dam south of San Francisco [Berilgen, 2007]. Riverbank stability, like with the slope stability analysis, may be performed using either the traditional limit equilibrium methods (LEM) or by using the numerical deformation finite element methods (FEM). Duncan [1996] reviewed a large number of limit equilibrium methods such as the Bishop's Modified Method, Force Equilibrium Methods, Janbu's Generalized Procedure of Slices, and Spencer's Method. He, also, outlined the advantages and limitations of the finite-element method for use in practical engineering problems. Although the limit equilibrium methods allow engineers to evaluate the stability factors quickly and simply, the finite element approach provides solutions to problems which cannot be solved by conventional methods, such as complex geometry, material anisotropy, non-linear behavior. The advantages of finite element methods over limit equilibrium methods are as follows [Griffiths and Lane, 1999]; (1) No assumption needs to be made in advance about the shape or location of the failure

1. INTRODUCTION Stability of riverbank depends on its geometry, its material properties, and the forces to which it is subjected. These forces include the effects of water both internally, in the form of pore-water pressures and seepage forces, and externally, in terms of hydrostatic and hydrodynamic forces [Lane and Griffiths, 2000]. Many riverbanks and channel slopes are frequently subjected to sudden changes of water level (rise or fall), like the riverbanks during flood events, reservoir banks during dam operation, and irrigation channels that are operated by rotation system. Such fluctuation in water level changes the flow conditions inside the bank and produces changes in hydraulic boundary conditions which, in turn, induce a change in the effective stresses, which finally leads to riverbank failure. Scientific research on the riverbanks and slope failure under different degrees of submergence and drawdown was first quantified by Morgernstern [1963]. Detailed research on water drawdown in the natural and constructed banks are investigated by many researchers, such as the Mississippi earth banks [Desai, 1971, 1972, 1977], river bank along the Rio 48


surface. (2) There is no need for assumptions about slice side forces since there is no concept of slices in the finite element approach. (3) The finite element solutions give information about deformations at working stress levels. (4) The finite element method is able to monitor progressive failure up to and including overall shear failure. Many methods for slope stability analysis using finite elements have been proposed during the last two decades. Among those methods, gravity increase method [Swan and Seo, 1999] and strength reduction method [Matsui and San, 1992; Griffiths and Lane, 1999] are considered the most widely used methods. In the gravity increase method, gravity forces are increased gradually until the slope fails then the factor of safety is defined as the ratio between the gravitational acceleration at failure (gf) and the actual gravitational acceleration (g). In the strength reduction method, soil strength parameters , are reduced until the slope becomes unstable, therefore, the factor of safety is defined as the ratio between the initial strength parameter and the critical strength parameter. The gravity increase method is used to study the stability of embankments during construction since it gives more reliable results while the strength reduction method is used to study the stability of existing slopes [Albataineh N., 2006]. Riverbank stability during rapid drawdown can be analyzed in two approaches; namely the effective stress approach and the total stress approach [Abramson Lee W. 2002]. Some researchers had investigated the problem of riverbanks or slope stability under various flow conditions. For example, Lane and Griffiths [2000] , Berilgen [2007], Alonso and Pinyol [2011] assessed and investigated the stability of slopes under drawdown conditions. Besides Nian et al. [2011], Fredlund et al. [2011] combined seepage and slope stability models to study the stability of bank slopes subjected to transient unsaturated seepage. In the current paper the riverbank stability is studied by using the finite element method with strength reduction technique and by considering the effective stress approach. Two finite element models are decoupled to predict the stability factor of safety (FOS) of riverbanks, namely, seepage model, and slope stability model. Decoupled solution means the governing equations of mathematical models are solved sequentially and not at the same time. The objectives of our paper are (1) understand the response of natural riverbank slopes to both the variation of water surface level in the river and the internal unsteady seepage conditions ,(2) determine the minimum factor of safety under unsteady flow

conditions,(3) investigate the most critical hydraulic condition that affect the riverbank stability.

2. SEEPAGE MODEL The finite element (FE) model introduced by [Istok, 1989] is employed to simulate 2-D transient saturated-unsaturated seepage flow through the riverbanks. The governing equation for water flow through saturated and unsaturated soil can be obtained by introducing Darcy’s law into the mass continuity equation. The general governing differential equation for two-dimensional seepage is given by Richards equation as: (

( )

)

(

( )(

))

(1)

( )

( ) where is the pressure head, [m]; ( ) are the unsaturated hydraulic conductivity in x, z directions respectively, [m/sec]; ( ) is the specific water capacity, [1/m]; is the volumetric water content [m3/m3] , ; n is the soil porosity; S is the degree of saturation (S ranges from 0.0 in dry soil to 1.0 in fully saturated soil); n, S are dimensionless, x, z are the cartesian coordinates,[m]; and t is time, [sec]. In an unsaturated soil, both the volume of water stored within the voids, and hydraulic conductivity will vary depending on the negative pore-water pressure. A conductivity function ( ) and a storage function ( ) should be well defined to complete the seepage model. These two functions could be obtained, for every soil type, by direct measurement in the laboratory and then inserted into the seepage model to get ( ) and ( ). Since it can sometimes be difficult or time consuming to obtain the storage and conductivity functions by laboratory measurements, the functions could be estimated by using closed-form solutions as proposed by [Van Genuchten, 1980] or [Fredlund and Xing, 1994]. In the current study, the air phase is assumed to be in an atmospheric state which means that the flow of air phase is negligible. Applying the Galerkin method to Eq.(1) results in systems of equations that can be written in a matrix form as: [ ( )]{ }

{ } ; where:

[ ( )] [ ( )] [ ( )] { } { } { } ([ ( )] ( ) [ ( )]){ } ){ } { } (( )

(2) (3)

where [ ( )] is the global capacitance matrix, [ ( )] is the global conductance matrix , is the 49


matrix (or the constitutive matrix) which, plane strain problems, is given by:

relaxation factor, is the time step, and {F} are the specified rates of seepage flow representing sources and sinks at Neumann nodes. {F} is known at all time steps. { } is known from the initial conditions or from the solution for the previous time step. { } contains the unknown values of pressure head at time . Eq.(2) is a system of nonlinear equations because the entries of the element capacitance and conductance matrices contain the specific water capacity and the hydraulic conductivity for the element, and for unsaturated flow these are functions of pressure head. The matrix formulation for the transient unsaturated flow equation can be solved using Picard iteration, [Istok, 1989]. The output from the seepage model is the pore-water pressure distribution inside the riverbank and the locations of phreatic surface with time.

[ ]

[ [ ]

Where constant)

∫ ∫[ ] [ ][ ]

]

(6)

(7)

]

[

]

3.2 Mohr-Coulomb Failure criteria and Strength Reduction Factor Mohr circle gives the relation between the strength parameters of soil ( (effective cohesion), (effective friction angle)), and the acting principal stresses. After obtaining the stresses and principal stresses inside the riverbank, Mohr-Coulomb's criterion, Eq.(8), will judge whether failure takes place at a certain node or not. are the major and minor principal and ( effective stresses). (8) If the failure function F < 0, then the stresses are inside failure envelope (elastic). If F = 0, then the stresses are on failure envelope (yielding). Finally if F > 0, then the stresses are outside failure envelope which means that yielding and failure took place at that location and stresses must be redistributed. The Factor of Safety (FOS) of a riverbank slope is defined as the factor by which the original shear strength parameters , must be divided in order to bring the slope to the point of failure [Griffiths and Lane, 1999]. The factored shear strength parameters , and are therefore given by:

(4)

where{ }and { }are the nodal displacements and force components for the element (e).The stiffness matrix [ ] could be expressed for an element as: [ ]

[

N1,N2,N3, ….Nn are the interpolation functions for a certain element with (n) nodes. After solving Eq.(4) for the nodal displacements { }, we can obtain strains { } and stresses { } for all elements by using:{ } [ ]{ },and{ } [ ]{ }.

3.1 Plain strain Model Stresses within the soil material of embankments and slopes could be estimated by assuming the soil as an elastic-perfect plastic material. Riverbanks can be treated as plane strain problem. In the plane strain riverbank problem, the strain normal to the cross section plane of riverbank and the shear strain are assumed to be zero. In the finite element formulation for an elastic soil element with (n) nodes, the displacements are related to the applied loads through the relation: { }

)

, and

[ ]

The stability model is based on the programs explained by [Smith and Griffiths, 2005]. It is a 2-D plane strain slope stability analysis of an elastic-plastic soil with a Mohr-Coulomb failure criterion and utilizing the strength reduction method. The model consists of three sub-models, (1) plane strain sub-model to calculate the effective stress distribution through the riverbank, (2) sub-model to compare the effective stresses with the Mohr-Coulomb failure envelop and determine which nodes have been yielded and overstressed, (3) Visco-plastic algorithm to redistribute the stresses of yielding nodes through the mesh.

}

)(

in which, E is the modulus of elasticity, and  is Poisson’s ratio, and [ ] [ ][ ]which are given by:

3. STABILITY MODEL

[ ]{

(

for

(

(5)

)

(9)

where SRF is a "Strength Reduction Factor", , are the original shear strength parameters. In

is the element thickness (assumed , and [D] is called the stress-strain 50


order to obtain a true factor of safety, the strength reduction factor needs to be gradually increased until the reduced strength parameters, and bring the slope to the failure state (numerical non convergence occurs within specified maximum number of iterations), at that time, the factor of safety for the slope is equal to the strength reduction factor and FOS=SRF. In this paper, the effective cohesion and effective friction angle are reduced by the same factor.

Step. 6

Step. 7

3.3 Elastic-visco-plastic model For the nodes that have a failure function F > 0, this node has been yielded and cannot resist loads anymore and the stress state can only be sustained momentarily and visco-plastic straining occurs. The magnitude of the visco-plastic strain rate is determined by the value of the failure function F. The load at the yielded nodes should be redistributed among the neighboring nodes which are not yielded yet. The elastic-plastic model involves five input parameters; elasticity modulus (E), Poisson’s ratio (), effective cohesion (c’), effective friction angle (’), and angle of dilation (*). Refer to Potts et al. [2001] for more details.

Step. 8

failure envelope (F≥0), then that location is assumed to be yielding. Yielding stresses are redistributed throughout the mesh utilizing the visco-plastic algorithm. The shear strength reduction factor SRF is increased incrementally until the global failure of the slope occurs, which means that the finite element calculation diverges. The lowest factor of safety of the embankment slope lies between the shear strength reduction factor SRF at which the iteration limit is reached, and the immediately previous one. The FOS for the slope is defined by division of the original shear strength parameters, c’, ’ by the values at failure, c’f , ’f.

5. NUMERICAL SIMULATION OF RIVERBANK STABILITY UNDER UNSTEADY FLOW Stability of riverbanks under three unsteady river flow conditions is simulated. These cases start with the case of constant water level in the river, then the case of filling the river to the riverbank top level and finally the case of drawdown. Two modes of drawdown were analyzed; the rapid drawdown and the slow drawdown, as shown in Table 1.

4. PROCEDURE Calculations of safety factor for riverbank can be summarized in the following steps: Step. 1 Define the geometry of riverbank, soil properties for each layer and water level in the river. Step. 2 Calculate the pore water pressure and water table variation with time by the transient-unsaturated seepage model. Step. 3 The soil is initially assumed to be elastic, calculate stresses by the plane-strain finite element model. Then obtain major and minor principal effective stresses ’1, ’3. Step. 4 Calculate the factored shear strength parameters cf’, f’ using Eq.(9) by selecting initially a small value of shear strength reduction factor SRF, for example 0.01, which makes the shear strength large enough to keep the slope in elastic stage. Step. 5 Substitute the principal stresses and the factored shear strength parameters (from steps 3,4) into the Mohr-Coulomb failure criterion in Eq.(8). If the stresses at a particular Gauss-point lie within the Mohr-Coulomb failure envelope (F1) when it is dry or fully submerged but became unstable (FOS