Numerical analysis of nonlinear dynamic behavior of earth dams

Abstract

A numerical study is conducted to investigate the dynamic behavior of earth dams. The numerical investigation employs a fully nonlinear dynamic finite difference analysis incorporating a simple elastic perfectly plastic constitutive model to describe the stress-strain response of the soil and the Rayleigh damping to increase the level of hysteretic damping. The extended Masing rules are implemented into the constitutive model to explain more accurately the soil response under general cyclic loading. The soil stiffness and hysteretic damping change with loading history. The procedures for calibrating the constructed numerical model with centrifuge test data and also a real case history are explained. For the latter, the Long Valley (LV) earth dam subjected to the 1980 Mammoth Lake earthquake as a real case-history is analyzed and the obtained numerical results are compared with the real measurements at the site in both the time and frequency domains. Relatively good agreement is observed between computed and measured quantities. It seems that the Masing rules combined with a simple elasto-plastic model gives reasonable numerical predictions. Afterwards, a comprehensive parametric study is carried out to identify the effects of dam height, input motion characteristics, soil behavior, strength of the shell materials and dam reservoir condition on the dynamic response of earth dams. Three real earthquake records with different levels and peak acceleration values (PGAs) are used as input motions. The results show that the crest acceleration decreases when the dam height increases and no amplification is observed. Further, more inelastic behavior and more earthquake energy absorption are observed in higher dams.

Keywords

earth dam / numerical / nonlinear response / dynamic analysis / earthquake / dam height

Introduction

The stability of earth dams subjected to dynamic action can be assessed through different approaches including the force-based pseudo-static methods, displacement-based sliding block methods and dynamic analysis [1,2]. The most common method, named the pseudo-static approach, is largely used in engineering practice to evaluate the seismic stability of earth fill dams. This approach takes into consideration many simplistic assumptions such as the modeling of the complex dynamic behavior in terms of static forces. Stability is expressed in terms of an overall factor of safety. The response of the dam subjected to an earthquake may be associated with different factors such as dam geometry, mechanical properties of construction soils, pre-seismic stress and pore water pressure distributions inside the dam, and input motion characteristics. Most of the factors are partially or totally ignored by the approaches traditionally adopted to assess earth dam seismic safety. The pseudo-static approach, for instance, does not take into account some dominant earthquake parameters such as frequency content and duration, known to significantly influence soil response [3]. On the other hand, the study of seismic response of earth dams is a complicated problem that commonly needs the use of dynamic analysis with different levels of sophistication in terms of proper problem formulation, characterization of material properties and modeling of stress-strain soil behavior. After the 1971 San Fernando earthquake in California, major advancements were attained in the understanding of earthquake action on dams [4]. Gazetas [5] discussed the historical developments of theoretical methods for estimating the dynamic response of earth dams to earthquake ground excitation and outlined their major features, their benefits and restrictions.

Generally, numerical techniques allow the most comprehensive analyses of the response of earth dams to dynamic loading. Progress in the area of computational geotechnics and numerical simulations offers powerful facilities to analyze the dam response and consider complex issues such as the soil nonlinearity, the development of the pore pressure during the dam construction procedure and real earthquake records.

Prevost and his co-workers [6] proposed two and three dimensional nonlinear dynamic finite element analyses of earth dams. Their investigation was based upon nonlinear hysteretic analyses employing a multi-surface plasticity theory. Afterwards, Lacy et al. [7] presented an efficient numerical methodology for analyzing the nonlinear dynamic response of earth dams. The appropriate coupled dynamic governing equations were outlined for the response of a soil system. Abouseeda et al. [8] investigated the nonlinear behavior of earth dam-foundation interaction during an earthquake using a coupled boundary element (BE) and finite element (FE) methods. Chen et al. [9] considered the stochastic response of the Santa Felicia earth dam in southern California with a spatially varying earthquake ground motion (SVEGM) model. They utilized a SVEGM model which takes into account both incoherence and propagation of seismic waves. Cascone et al. [10] studied the dynamic stability of earth dams using decoupled displacement analysis. Ming et al. [11] carried out a fully coupled analysis of earth dam failure and investigated remediation of the Lower San Fernando dam. To describe soil behavior over the wide range of loading conditions, a critical state model incorporating the concept of state dependent dilatancy was adopted. Adalier et al. [12] considered that the dynamic behavior of an embankment dam rested on a liquefiable foundation. Papalou et al. [13] investigated the nonlinear dynamic response of earth dams considering canyon interaction. They extended a finite element based method in which the dam was simulated as a shear beam and the surrounding medium as a half-space. The nonlinear behavior of the dam was studied using multi-yield surface plasticity theory. Ebrahimian et al. [14] studied the seismic response of earth dams during an earthquake using 2D fully nonlinear dynamic analysis. They used the finite difference method and adopted the Mohr-Coulomb elastic-perfectly plastic constitutive model to describe the stress-strain relation of the soil. They laid special attention on the nonlinear dynamic behavior of very high dams. Wang et al. [15] performed effective-stress-based dynamic analyses incorporating soil nonlinearity. They employed a bounding surface hypoplasticity model for sand and implemented the model into a finite difference code. The advantages of the proposed nonlinear approach were demonstrated by comparison with the results of the equivalent linear approach for a rockfill dam. The seismic performance of an earth dam was also evaluated to illustrate the model’s capability. Siyahi et al. [16] carried out a transient dynamic time history finite element simulation to investigate the performance of earth dams under seismic excitation. Then, they studied the different failure modes of earth dams as a consequence of earthquake shaking [17]. They presented a case study as an example. Sica et al. [18] considered the influence of past loading history on the seismic response of earth dams by interpreting the static and dynamic behavior of a real case-history. They utilized a coupled dynamic approach. The approach was solved numerically with the finite element method. Rampello et al. [19] studied the response of an earth dam to seismic loading through displacement-based analyses and finite element effective stress dynamic analyses. The finite element analyses were carried out using a constitutive model capable of reproducing soil nonlinearity, calibrated against laboratory measurements of the stiffness at small strains. The influence of the assumed input motion and bedrock depth on the seismic response of the dam was also investigated. Ebrahimian [20] presented a numerical modeling of the seismic behavior of an earth dam resting on liquefiable foundation. Numerical simulation was conducted using effective stress-based, fully coupled nonlinear dynamic analysis. The Finn-Byrne model with extended Masing rules was used to model pore pressure generation in the liquefied soils.

This paper presents a numerical study of the seismic behavior of earth dams overlaying bedrock subjected to real earthquake records using fully nonlinear dynamic analysis. The effect of nonlinear soil behavior is then accounted for the analyses from the early beginning of earthquake loading. The extended Masing rules [21] are implemented into the constitutive model to precisely explain the soil response under general cyclic loading. The soil stiffness and hysteretic damping change with loading history. The numerical analyses are carried out using a simple elastic perfectly plastic soil model with isotropic hardening, able to reproduce some aspects of soil behavior during cyclic loading. The procedure uses a nonlinear stress-strain relation incorporated into the finite difference computer program FLAC [22]. It will mainly focus on considering the effects of dam height, input motion characteristics, soil behavior, strength of the shell materials and dam reservoir condition on the dynamic response of earth dams. It is noted that this study does not consider the fluid-skeleton interaction, which may have a significant influence on the seismic response of the dam. At first, a numerical model is calibrated against the centrifuge model test and the Long Valley (LV) earth dam as a real case-history and some important features related to model calibration are discussed. Comparisons between the results obtained from the present study and the experimental observations and also the real measurements show that the current numerical procedure can capture fundamental aspects of the dynamic behavior of earth dams accurately.

The selected problem is a simplified representation of typical earth dam geometry. The dam section is a symmetric zone section with the central clay core rested on bedrock as shown in Fig. 1. Five earth dam cross sections with different heights, i.e., H = 40 m, 80 m, 120 m, 200 m and 280 m, are analyzed.

Numerical modeling procedure

Numerical analyses are carried out using the finite difference program FLAC based on a continuum finite difference discretization using the Langrangian approach [22]. Every derivative in the set of governing equations is replaced directly by an algebraic expression written in terms of the field variables (e.g. stress and displacement) at discrete points in space. For dynamic analysis, an explicit time marching scheme is utilized to solve the equation of motion using lumped grid point masses extracted from the real density surrounding zone. The calculation sequence first recalls the equations of motion to obtain new velocities and displacements from stresses and forces. Then, strain rates are derived from velocities, and new stresses from strain rates. Every cycle around the loop corresponds to one time step. Each box updates all of its grid variables from known values that remain fixed over the time step being executed (see Fig. 2).

Newton’s law of motion for the body relates the acceleration, {\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}t}}, of a mass, m, to the acting force, F. The equation is given as

(1)m{\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}t}}=F.

In a continuous solid body, Eq. (1) can be written as

(2)\rho {\displaystyle \frac{\partial {\dot{u}}_i}{\partial t}}={\displaystyle \frac{\partial {\sigma }_{ij}}{\partial x_j}}+\rho g_i,

where ρ = mass density, t = time, x_j = components of the coordinate vector, g_i = components of gravitational acceleration (body forces) and σ_ij = components of the stress tensor. In this equation, and those that follow, indices i denote components in a Cartesian coordinate frame, and summation is implied for repeated indices within an expression. To analyze a problem, the strain rate tensor and the spin tensor, having the velocity gradient, can be calculated from the following equations:

(3){\dot{e}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\dot{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\dot{u}}_j}{\partial x_i}}\right),

(4){\dot{\omega }}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\dot{u}}_i}{\partial x_j}}-{\displaystyle \frac{\partial {\dot{u}}_j}{\partial x_i}}\right),

where, {\dot{e}}_{ij} are the components of the strain rate and {\dot{u}}_iare the components of the velocity.

Constitutive model

Generally, the mechanical constitutive law has the following form:

(5){\sigma }_{ij}:=M_{ij}\left({\sigma }_{ij},{\dot{e}}_{ij},\kappa \right),

where M_ij is a tensor valued function depending on the stress, strain rate, and κ as a history parameter (based on the specific rules which may or may not exist). It is noted that := means “replaced by”.

In this study the Mohr-Coulomb constitutive relation is used to model the behavior of soil. The failure envelope for this model corresponds to a Mohr-Coulomb criterion (shear yield function) with tension cutoff (tensile yield function). The stress-strain relationship is linear elastic-perfectly plastic. The linear behavior is defined by elastic shear and bulk modulus. The plastic behavior is defined by the angle of internal friction and cohesion of the soil. The shear modulus of sandy soil as shell materials is calculated with the formula given [23]:

(6)G_{\max }=\mathrm{8400}{\displaystyle \frac{{(\mathrm{2.17}-e)}^{\mathrm{2}}}{\mathrm{1}+e}}{({\sigma }_{\mathrm{m}}^{\prime })}^{\mathrm{0.6}},

where G_max is the maximum (small strain) shear modulus in kPa, e is the void ratio and {\sigma }_{\mathrm{m}}^{\prime } is the mean effective confining stress in kPa. The Poisson’s ratio is taken as 0.3. The shear modulus of clayey soil as core materials is calculated with the following formula [24]:

(7)G_{\max }=\mathrm{3270}{\displaystyle \frac{{(\mathrm{2.973}-e)}^{\mathrm{2}}}{\mathrm{1}+e}}{({\sigma }_{\mathrm{m}}^{\prime })}^{\mathrm{0.5}}.

The Poisson’s ratio for the core is taken as 0.45.

To provide a constitutive model that can better fit with the curves of the shear modulus reduction and damping ratio derived from the experimental tests data, the basic elastic perfectly plastic model is modified to assess the potential for predicting dynamic behavior and associated deformations. To more accurately represent the nonlinear stress-strain behavior of soil that follows the actual stress-strain path during cyclic loading, the Masing behavior is implemented into FLAC via a FISH subroutine. The Masing model consists of a backbone curve and several rules that describe unload-reload behavior and cyclic modulus reduction. The backbone curve can be established from stiffness modulus degradation curves coupled with the small strain modulus (G_max). The rules related to unload-reload paths can similarly be formulated to reproduce hysteretic damping values expected from standard curves of damping ratio versus shear strain [25,26]. These formulations are described later.

In this study, the shear modulus reduction and damping ratio curves for sandy soils proposed by Seed et al. [25] and for clayey soils proposed by Vucetic et al. [26] are adopted as a reference. Geotechnical properties used in the analyses are presented in Table 1 for the earth dam materials.

Boundary condition

Many geotechnical problems can be idealized by assuming that the regions remote from the area of interest extend to infinity. To minimize the computation time and also to avoid the reflection of outwards propagating waves back into the model, the unbounded theoretical models have to be truncated to a manageable size by using artificial boundaries. The viscous boundary developed by Lysmer et al. [27] is used in the current calculations. In this case, independent dashpots are applied in the normal and shear directions at the lateral boundaries as shown in Fig. 3. During the static analysis, the bottom boundary is fixed in both horizontal and vertical directions and the lateral boundaries are just fixed in the horizontal direction. In dynamic analysis when the dam is laid on foundation (not bedrock), the lateral boundaries are changed to the free-field boundary available in the FLAC library, to eliminate wave reflections from the truncated boundaries.

Element size

Numerical distortion of the propagating wave can occur in a dynamic analysis as a function of the modeling conditions. Both the frequency content of the input excitation and the wave speed through the system will affect the numerical accuracy of wave transmission. Kuhlmeyer et al. [28] showed that for an accurate representation of the wave transmission through the soil model, the spatial element size must be smaller than approximately one-tenth to one-eighth of the wavelength associated with the highest frequency component of the input wave i.e.,

(8)\mathrm{\Delta }L={\displaystyle \frac{\lambda }{\mathrm{9}}},

where λ is the wave length associated with the highest frequency component that contains appreciable energy. Considering the above criteria, element size is defined small enough to allow seismic wave propagation throughout the analysis.

Damping

In general, material damping in a soil body is generated by the viscous properties, internal friction, and development of plasticity. Indeed, the role of damping in the numerical models is to cause energy losses in the natural system when it is subjected to a dynamic load. The dynamic damping in the model is provided by the Rayleigh damping option available in FLAC. Rayleigh damping is originally used in the analysis of structures and elastic continua, to damp the natural oscillation modes of the system. Rayleigh damping of R_d = 5% is adopted in the analyses to compensate for the energy dissipation through the medium [22]. In the dynamic analyses which incorporated plasticity constitutive models, a considerable amount of energy dissipation can occur during plastic flow. Thus, in this case the hysteretic damping is also considered in the calculations and only a minimal percentage of damping (e.g., 2%) may be required. The natural frequency of the dam as a Rayleigh damping parameter is determined by a Fourier analysis of the free response of the dam (see Fig. 4). This figure shows a fundamental frequency f₁ = 1.71 Hz for a dam with a height of 40 m. The fundamental frequencies of dams with different heights are tabulated in Table 2.

Time step

To complete the numerical solution, it is necessary to integrate the governing equations in time in an incremental manner. The time step used should be sufficiently small to accurately define the applied dynamic loads and to ensure stability and convergence of the solution. In the current calculations, the time step is approximately 10^-6 s.

Input excitation

One of the important decisions in the dynamic assessment process is the selection of the input motion. Usually in nonlinear dynamic analysis, the expected earthquake must be expressed as a set of ground motion time histories. For a safe design, it is important that among the acceptable acceleration time histories, the selected input motions generate an appropriate range of dam responses bounded by a reasonable conservatism in the assignment of seismic hazard parameters. Although this problem may seem intractable because of the number of time-history that could be associated with a given design/evaluation earthquake, rationally there exists a bound for the level of earthquake response likely to be achieved by the physical system. The quantification of this response requires a good understanding of the dynamic response of the system and the ground motion parameters that identify the damage potential of the input excitation [29]. Using acceleration time history as an expected earthquake ground motion, different parameters can be used to specify its intensity and damage potential. One of the noticeable parameters is the peak acceleration value (PGA). Using this factor is inherently natural since accelerations and the resulting inertial forces are directly related to Newton’s Second Law. However, there is not a direct correlation between the PGA and the structural response at the dominant natural frequencies of most typical dams. Furthermore, considering a large PGA is not sufficient to produce a critical response leading to significant damage. In spite of these shortcomings, PGA is one of the fundamental parameters used to evaluate the damage potential of a given acceleration time history. The frequency content of the earthquake ground motion has a significant influence on the dynamic response of the dam. Therefore, a better recognition of a given input motion can be achieved by using some form of spectral representation. In particular, the use of the Fourier amplitude spectrum (FAS) is at the core of earthquake engineering practice. However, these characterizations do not provide a direct description of the duration or the time varying features of a given input motion.

In this paper, to cover a wide rage of the dam responses, three different real acceleration time histories including Tabas (PGA= 0.93g in MCE level), Naghan (PGA= 0.72g in MDE level) and San Fernando (PGA= 0.21g in DBE level) are selected from a database of earthquake records. In the dynamic analyses of dams, the scaled records have been filtered to a maximum frequency of 10 Hz, transferred to the “inside” bedrock formation through a standard de-convolution analysis and applied at the base of the numerical model. Pertinent information on the earthquake records are summarized in Table 3, and the corresponding acceleration time histories and Fourier amplitude spectra are depicted in Figs. 5 and 6, respectively.

Fully nonlinear dynamic analysis

One of the common methods to determine the dynamic response of earth structures is the so-called equivalent linear analysis. In this approach, a linear analysis is first carried out using the initial values of the damping ratio and shear modulus. Then considering the maximum value of the shear strain, and using the curves obtained from laboratory test data, new values of damping ratio and shear modulus are calculated. These values are used to re-do the analysis. This procedure is repeated several times till no variation happens in the material properties. Thus, the method does not capture directly any nonlinear effects because it assumes linearity during the solution process; strain-dependent modulus and damping functions are only taken into account in an average sense, in order to approximate some effects of nonlinearity (damping and material softening).

In the nonlinear analysis which is employed in this study, the nonlinear stress-strain relationship is followed directly by each zone. The shear modulus and damping ratio of the materials are calculated automatically at different strain levels. The real behavior of soils is nonlinear and hysteretic under cyclic loading. This behavior can be simulated by the Masing model [21], which is able to reproduce the dynamic behavior of soils. In this model, the shear behavior of soil is expressed by a backbone curve in the form of

(9)F_{\mathrm{b}\mathrm{b}}(\gamma )={\displaystyle \frac{G_{\max }\gamma }{\mathrm{1}+(G_{\max }/{\tau }_{\max })\left|\gamma \right|}},

where F_bb(γ) is the backbone or skeleton function, γ is the shear strain amplitude, G_max is the initial shear modulus and τ_max is the maximum shear stress amplitude.

The stress-strain response follows the backbone curve in the first loading step (see Fig. 7(a)), but to explain the unloading-reloading paths, the aforementioned equation has to be revised. The next unloading and reloading curves have, necessarily, translation and dilatation in comparison with the backbone curve. In this case, if load reversal happens at a special point (τ_r, γ_r), the stress-strain curve is given by the following relation:

(10){\displaystyle \frac{\tau -{\tau }_{\mathrm{r}}}{\mathrm{2}}}=F_{\mathrm{b}\mathrm{b}}\left[{\displaystyle \frac{\gamma -{\gamma }_{\mathrm{r}}}{\mathrm{2}}}\right].

The shape of the unloading and reloading paths is identical to the shape of the initial backbone curve (with the origin shifted to the loading reversal point) but is enlarged by a factor of 2 (see Fig. 7(b)). This means that the dilation is described by a hysteresis scale factor. These first two rules describe Masing behavior [21].

These basic Masing rules are not sufficient to precisely describe the soil response under general cyclic loading. Finn et al. [30] presented additional rules to describe irregular loading conditions. They proposed that unloading and reloading curves follow two extra rules.

If the new unloading or reloading curve passes the last maximum strain and intersects the backbone curve, then it continues along the backbone curve till it reaches the next returning point (see Fig. 7(c)). If a new unloading or reloading curve exceeds the previous unloading or reloading curve, it will follow the former stress-strain curve (see Fig. 7(d)). According to the above rules, the tangent shear modulus at a point on the backbone curve can be expressed as below:

(11)G_{\mathrm{t}}={\displaystyle \frac{G_{\max }}{{\left(\mathrm{1}+{\displaystyle \frac{G_{\max }\left|\gamma \right|}{{\tau }_{\max }}}\right)}^{\mathrm{2}}}}.

The tangent shear modulus, at a point on the new reloading-unloading curve can be also described by the following equation:

(12)G_{\mathrm{t}}={\displaystyle \frac{G_{\max }}{{\left(\mathrm{1}+{\displaystyle \frac{G_{\max }}{\mathrm{2}{\tau }_{\max }}}\left|\gamma -{\gamma }_r\right|\right)}^{\mathrm{2}}}}.

Based on the research results, as the number of loading cycles increases, the shear stress decreases; this means that the shear stress-strain curves get more inclined. Actually, models that follow these four rules are often called extended Masing models. They inherently incorporate the hysteretic nature of damping and the strain-dependence of the shear modulus and damping ratio. To simulate the nonlinear shear stress-strain relationship in this study, the Masing rules are implemented into FLAC via a series of FISH functions.

Validation analysis

To validate the implementation of the Masing rules in the FLAC program, the simulation of a one-zone sample is conducted by using the unit cell as shown in Fig. 8. The one-zone sample consists of sandy soil with periodic motion exerted at its base. Vertical loading is established by gravity only. Equilibrium stresses are installed in the soil. The stress/strain loops of the one-zone sample for several cycles are shown in Fig. 9(a). It can be observed that shear modulus decreases with increasing shear strain. The hysteretic model seems to handle multiple nested loops in a reasonable manner. There is clearly energy dissipation and shear stiffness degradation during seismic loading as shown in Fig. 9(b). As it is clear, the shear modulus reduction curve obtained from the current study follows well the empirical relation proposed by Seed et al. [25] and test data.

To evaluate the applicability of the proposed model, the results obtained from numerical analysis are compared with the experimental counterparts. One of the centrifuge tests related to the embankment performed in the VELACS project (Verification of Liquefaction Analysis using Centrifuge Studies) [31,32] is chosen to calibrate the constructed numerical model in FLAC and also evaluate the ability of the constitutive model in predicting the dynamic response of the dam during seismic loading. It is attempted to create almost similar conditions between the laboratory model test and the numerical model. The model test configuration is depicted in Fig. 10(a). The numerical model is shown in Fig. 10(b). The numerical results are presented and compared to those obtained from the corresponding centrifuge test data. Comparisons between the computed and measured results (computed: numerical and measured: centrifuge test results) are made in Figs. 11 and 12. These comparisons show that the reference numerical model can predict the dynamic behavior of the earth dam in a rational way.

The last verification is related to a real well-documented case history to show the ability of the model to simulate dynamic behavior of an earth dam during a real earthquake. Thus, the results of nonlinear dynamic analyses of the LV earth dam in California subjected to the 1980 Mammoth Lake earthquake [33] are presented and compared with the real measurements recorded at the site and also the results presented by previous authors available in literature. The LV dam is located in the Mammoth Lake area (California) in close proximity to active faults. The dam is a rolled earth-fill with an impervious zone that forms the major portion of the embankment. The dam has a maximum height of 55 m, a length of 182 m at the crest, and upstream and downstream slopes of 3h/1v. The LV dam was instrumented in the 1970’s with a multiple-input-output array consisting of 3 accelerometer stations to monitor the boundary conditions and 5 stations to record the dam response (see Fig. 13). Thus, the array comprised a total of 22 accelerometers linked to a common triggering mechanism.

In May 1980, a series of 6 earthquakes occurred in the Mammoth Lake area. The magnitudes of these earthquakes ranged from M_L = 4.9 to M_L = 6.7, and the induced peak accelerations at the crest center reached 0.5g in the upstream-downstream direction (x direction, see Fig. 13) during the strongest event. An extensive array of 22 input-output (excitation-response) accelerations was recorded, providing a valuable source of information on the dam dynamic response over a wide range of deformation levels. In this study, the dam is subjected to the input motion recorded downstream at the outlet during Mammoth Lake earthquakes. The first 12 s of the recorded acceleration is used with data points at 0.02 s intervals and peak acceleration equal to 0.135g in the upstream-downstream direction and 0.084g in the vertical direction.

Detailed information on the LV dam can be found in Ref. [33]. The LV dam cross section is shown in Fig. 14. The numerical grid constructed in FLAC is presented in Fig. 15. The input accelerations are applied in the horizontal and vertical directions of the model base. A free field boundary contrition is exerted to the lateral boundaries of the model. Of particular interest is the computed acceleration at the crest which can be compared directly with measured values. Previous analyses of the LV dam have been reported by Lai et al. [34], Elgamal et al. [35], Griffiths et al. [33], Yiagos et al. [36], Zeghal et al. [37], Woodward et al. [38]. The first natural frequency obtained from the current study is presented in Table 4 and compared with other solutions available in the literature. The present study gives reasonably close agreement with the other numerical investigations. The crest acceleration responses of the LV dam are computed and compared with the motions recorded at the site in the both time and frequency domains.

Figure 16 shows the computed horizontal acceleration of the crest and indicates the amplification that has occurred between the base and the crest. The peak amplitude at the crest has a magnification factor of about 5.47 over the peak base amplitude. The computed response of the crest in the horizontal direction is compared with measured values in Fig. 17. The computed response is plotted with a dashed line. Excellent overall agreement is achieved, with the computed values giving somewhat higher amplitudes. The frequency content of the two time records is compared in the form of a Fourier amplitude spectrum in Fig. 18. The peaks are in close agreement although the computed values show rather more energy associated with the fundamental frequency around 1.8 Hz. The frequency content of the up/down stream motion presented in Fig. 18 shows that the energy is concentrated at a frequency of just under 2 Hz. The calculated acceleration in the vertical direction shows considerably less agreement with measured values. Part of the difficulty is illustrated in Fig. 19 which shows superimposed plots of vertical acceleration at the base and crest. This excitation is considerably noisier than in the horizontal direction and is less intensive. The maximum recorded vertical acceleration at the crest is 0.172g, compared with 0.64g in the horizontal direction. The computed accelerations in the vertical direction are compared with measured values in the crest of the LV dam in Fig. 20. The computed values show generally lower amplitudes than the measured values. The Fourier amplitude spectrum of these time histories is given in Fig. 21. The measured values indicate a broad band of frequencies with no particular frequency dominating the situation. The computed values also contain a broad band of frequencies, but with clear peaks in the ranges 2-3 Hz and 5-6 Hz. The frequency content of the vertical acceleration in the form of Fourier response spectra indicates that the computed values do not reproduce the higher frequencies present in the broad band of measured frequencies. It is noticed that the time and frequency-domain results give good agreement and high correlation in the horizontal direction but rather poor agreement in the vertical direction.

The results obtained from validation analysis of the LV dam in terms of crest acceleration are compared with the other numerical results presented by previous authors available in the literature and summarized in Table 5. These comparisons show that the current numerical procedure can capture fundamental aspects of the dynamic behavior of earth dams rather well. Due to the satisfactory modeling of the validation cases, the numerical model is then used to perform a parametric study on the hypothetical earth dam, as described earlier.

Numerical results

Numerous analyses are carried out to investigate the effects of dam height, input motion characteristics, soil behavior, strength of shell materials and dam reservoir condition on the dynamic behavior of earth dams. The effects of different types of earthquakes, mentioned before, on the horizontal permanent deformations, permanent shear strains and induced maximum accelerations at the crest of dams with different heights are shown in Fig. 22. Besides the displacements (see Fig. 22(a)), the relevant shear strains are also presented in Fig. 22(b). It is clear that the shear strain variation is similar to displacement. The horizontal displacements and shear strains in the dam body are increased with increasing dam height. The calculated quantities are much larger for the Tabas earthquake and failure occurs in the dam body. Fig. 22(a) shows that at the end of the Tabas earthquake, the computed maximum horizontal displacement at the crest of the dam is about 94 cm. It can be observed that the increase in the input motion energy leads to an important increase in the displacements and shear strains. Fig. 22(c) illustrates the coupled effect of the dam height and the type of earthquake on the induced maximum acceleration at the dam crest. It is noticed that the crest acceleration is reduced when the dam height increases and no amplification is seen. The reason may be attributed to more flexible behavior, larger damping and larger developed plastic zones observed in higher dams. Thus, because of these factors, there is more energy absorption in higher dams with respect to smaller dams. It can be seen that the reduction of accelerations in the dam crest is greater for higher dams than smaller ones. It is noteworthy that the PGA of the Naghan earthquake (0.72g) is much higher than the San Fernando earthquake (0.21g) but the created displacements and shear strains in the dam crest due to the Naghan earthquake is close to those obtained from the San Fernando input motion. It can be concluded that the PGA is not a sufficient parameter for gauging the potential of a particular earthquake time history to produce significant permanent deformations. Thus, the other earthquake parameters such as effective duration, magnitude, and frequency content should be considered in the analysis.

The pattern of failure mechanism with permanent shear strain contour in the dam body is shown in Fig. 23 for two different heights at the end of the Naghan earthquake. It is seen that failure occurs in the higher dam (280 m) compared with the smaller one (120 m). The slip surface in the dam with a 280 m height (see Fig. 23(b)) is much deeper and clearer than in the 120 m one (see Fig. 23(a)).

To investigate the effect of soil behavior on the dynamic response of a dam body, the dam with a height of 40 m subjected to the mentioned earthquakes is chosen as a reference with two different behaviors including elastic and elastic-plastic treatments. As it is expected, for elastic behavior, smaller displacements and shear strains are observed along the dam height (see Figs. 24(a) and (b)) but large amplification occurs especially for the strongest earthquake (see Fig. 24(c)). It means that plasticity reproduced more energy dissipation during dynamic loading. Thus, in this case, the accelerations are reduced across the dam height compared to the base acceleration, and larger displacement occurs in the dam. It is shown in Fig. 24(a) that for elastic behavior, maximum displacement occurs at about Z/H = 0.88 but for elastic-plastic behavior, maximum displacement happens at the crest of the dam. This figure shows that for elastic-plastic behavior a large increase happens for the dynamic induced residual (permanent) displacement in the upper part of the dam especially for the Tabas and Naghan earthquakes. In fact, previous research works have indicated this kind of behavior [39-41]. This is why in the design of embankment dams, due to the stronger shaking at the upper parts, special attention should be given to the crest to avoid undesirable deformations. The distribution of shear strain along the dam height is extremely nonlinear for stronger earthquakes (see Fig. 24(b)). The maximum acceleration in elastic dams occurs in the dam crest, as shown in Fig. 24(c). For the Tabas earthquake, there is a special increase in the acceleration profile along the dam centerline at Z/H = 0.38.

One of the important parameters which can significantly affect the dynamic response of the dam is the strength of dam materials. To clarify the effect of strength of shell materials, different friction angles are assumed in company with different dam heights subjected to the Naghan earthquake. Fig. 25(a) shows the quantities of horizontal displacement versus the dam height for different friction angles of shell materials. The variation of shear strain in the crest of dams with different heights for different friction angles of shell materials is indicated in Fig. 25(b). As expected, when the friction angle increases, the horizontal displacement and shear strain in the dam crest decrease. It can be seen that the variation in friction angle does not induce a significant change in the displacement and shear strain for Φ≤40°. However, the variation for Φ = 45° severely changes in comparison to the lower friction angles. The greatest displacement values are those developed at the crest of the dam with Φ = 30°. The horizontal displacement at the crest is computed to be about 16 cm and 13 cm for dams with heights of 40 m and 120 m, respectively. The computed shear strains at the crest are about 3.5e^-3 and 2.5e^-3 as well. The maximum induced acceleration at the top of the dam decreases as the friction angle decreases or the dam height increases (see Fig. 25(c)). Consideration of larger friction angles for the shell materials (e.g., Φ = 45°) leads to an increase of about 70% in the dynamic amplification. For Φ = 45°, the computed maximum crest acceleration of the dam with 40 m height is about 0.89g, compared to 0.52g for the dam with 120 m height. When the dam height decreases, the horizontal displacement and shear strain increase but the acceleration is reduced at the crest of the dam. The variation of all quantities for Φ = 45° is linear but for the other friction angles are slightly nonlinear, as shown in Fig. 25.

To clarify the influence of dam reservoir condition on the dynamic response of earth dams with various heights subjected to the Naghan earthquake, three different conditions are assumed for the dam reservoir including an empty reservoir (after dam construction), full reservoir and full reservoir with embedded flow net in the dam body (after seepage analysis). It is observed in Fig. 26(a) that the maximum horizontal displacements are related to a full reservoir with seepage in the dam body. It is noted that horizontal displacements for the empty reservoir are less than the case including the full reservoir without seepage flow in the dam body. As shown in Fig. 26(b), the same pattern is seen for shear strain in the crest of the dam. When the reservoir is full and also flow net is established in the dam body, the variation of displacement and shear strain changes so rigorously in contrast to the other cases (full and empty reservoir). Computed permanent horizontal displacements at the crest of the dam with full reservoir coupled with flow net is about 19 cm; computed shear strain at the crest of the dam is about 3.7e^-3. Maximum crest acceleration occurs for the full reservoir with flow net (see Fig. 26(c)). Similar to previous results, when the dam height increases, the crest acceleration decreases. The variation of displacement and shear strain versus dam height is nonlinear for a full reservoir with flow net in the dam body but linear for the other cases.

Conclusions

This paper presents nonlinear dynamic behavior of earth dams using an explicit finite difference method. A simple elastic perfectly plastic constitutive model with Mohr-Coulomb failure criterion is used to describe the stress-strain response of the soil. Rayleigh damping is utilized to increase the level of hysteretic damping during dynamic analysis. The Masing rules are implemented into the constitutive model to precisely explain the nonlinear soil response under general cyclic loading. The numerical model is then calibrated using centrifuge test data and also the field data obtained from the real measurements of the LV earth dam subjected to the 1980 Mammoth Lake earthquake. Comparisons between the dynamic analyses results obtained from the present study and the real measurements of the LV dam have been carried out in terms of accelerations computed at the crest of the dam in both time and frequency domains. The validation analyses confirm that the proposed numerical model reproduces the overall dynamic behavior of the earth dams under earthquake loading conditions well, qualitatively as well as quantitatively. After validation, a detailed parametric study has been performed to evaluate the influence of dam height, real earthquake loading, soil behavior, strength of shell materials and dam reservoir condition on the dynamic response of earth dams. Particular attention has been given to the influence of dam height on the nonlinear dynamic behavior of dams. The following conclusions can be drawn on the basis of the performed parametric study:

1) Nonlinear dynamic analysis shows that plasticity should be considered in the investigation of the seismic response of earth dams, because it leads to a decrease in the acceleration of the dam crest, an increase to displacements and shear strains of the dam body and an increase in energy dissipation, which all can significantly affect the dynamic response of earth dams.

2) The higher dams show more inelastic behavior than the smaller ones; as a result, this affects the shear strains which influence the shear modulus degradation and the attenuating coefficient, and all these effects tend to weaken the accelerations along the height.

3) If the dam materials keep their elastic behavior during dynamic loading, then the horizontal acceleration becomes larger along the dam height (from the base to the top). In this case, the higher dams show larger amplification, especially if the natural period of the dam body coincides with the periodical nature of earthquake waves.

4) When the dam body shows nonlinear treatment or materials go towards plastic behavior during a strong shaking, the attenuation of acceleration waves in the dam body becomes more effective, and consequently the earthquake accelerations descend in amplitude when passing from the base towards to the top. Large displacements and shear strains occur in this case and failure happens in the dam body.

5) Nonlinear elastic-plastic analyses show that the strongest dynamic loading (Tabas earthquake) induces plasticity in large parts of the dam body when the height of the dam increases. Consequently, the displacements and shear strains increase significantly in the dam.

6) The soils with less strength (suppose low friction angle) go towards yielding by small amounts of dynamic force, and this also causes the attenuation of acceleration along the dam height in the weaker materials compared to the stronger ones.

7) When the dam is subjected to the earthquakes with lower energy, the dam body behaves as an elastic material and thus the induced seismic accelerations inside the dam body become larger from the base to the top. In this case, small plasticity zones can be developed in the dam body and then the dam remains safe during dynamic loading.

8) When the dam reservoir is full and also the flow net is established in the dam body, displacements, shear strains and induced accelerations increase in the dam body in comparison with an empty reservoir and full reservoir without flow in the dam. It is noteworthy that the quantities for displacement, shear strain and acceleration for a full reservoir without flow in the dam are less than in the case of an empty reservoir (after dam construction).

References

[1] Sherard J L. Earthquake considerations in earth dam design. Journal of the Soil Mechanics and Foundations Division, 1967, 93: 377-401

[2] Seed H B, Lee K L, Idriss I M, Makdisi F I. The slides in the San Fernando dams during the earthquake of February 9, 1971. Journal of the Geotechnical Engineering Division, ASCE, 1975, 101(GT7): 651-688

[3] Ambraseys N N. The seismic stability of earth dams. In: Proceedings of the Second World Conference on Earthquake Engineering, Tokyo, Japan. 1960, 1345-1363

[4] USCOLD (US Committee on Large Dams). Observed Performance of Dams during Earthquakes. Committee on Earthquakes, Denver, CO, 1992

[5] Gazetas G. Seismic response of earth dams: some recent developments. Soil Dynamics and Earthquake Engineering, 1987, 6(1): 2-47

[6] Prevost J, Abdel-Ghaffar A M, Lacy S. Nonlinear dynamic analyses of an earth dam. Journal of Geotechnical Engineering, 1985, 111(7): 882-897

[7] Lacy S J, Prevost J H. Nonlinear seismic response analysis of earth dams. Soil Dynamics and Earthquake Engineering, 1987, 6(1): 48-63

[8] Abouseeda H, Dakoulas P. Non-linear dynamic earth dam-foundation interaction using a BE-FE method. Earthquake Engineering & Structural Dynamics. 1998, 27(9): 917-936

[9] Chen M, Harichandran R S. Response of an earth dam to spatially varying earthquake ground motion. Journal of Engineering Mechanics, 2001, 127(9): 932-939

[10] Cascone E, Rampello S. Decoupled seismic analysis of an earth dam. Soil Dynamics and Earthquake Engineering, 2003, 23(5): 349-365

[11] Ming H Y, Li X S. Fully coupled analysis of failure and remediation of lower San Fernando dam. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(4): 336-349

[12] Adalier K, Sharp M K. Embankment dam on liquefiable foundation-dynamic behavior and densification remediation. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(11): 1214-1224

[13] Papalou A, Bielak J. Nonlinear seismic response of earth dams with Canyon interaction. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(1): 103-110

[14] Ebrahimian B, Vafaeian M. Effects of dam height on the seismic response of earth dam. In: Proceedings of the 7nd International Conference on Civil Engineering, Tarbiat Modarres University, Tehran, Iran. 2005

[15] Wang Z L, Makdisi F I, Egan J. Practical applications of a nonlinear approach to analysis of earthquake-induced liquefaction and deformation of earth structures. Soil Dynamics and Earthquake Engineering, 2006, 26(2-4): 231-252

[16] Siyahi B, Arslan H. Nonlinear dynamic finite element simulation of Alibey earth dam. Environmental Geology (Berlin), 2008, 54(1): 77-85

[17] Siyahi B, Arslan H. Earthquake induced deformation of earth dams. Bulletin of Engineering Geology and the Environment, 2008, 67(3): 397-403

[18] Sica S, Pagano L, Modaressi A. Influence of past loading history on the seismic response of earth dams. Computers and Geotechnics, 2008, 35(1): 61-85

[19] Rampello S, Cascone E, Grosso N. Evaluation of the seismic response of a homogeneous earth dam. Soil Dynamics and Earthquake Engineering, 2009, 29(5): 782-798

[20] Ebrahimian B. Numerical modeling of the seismic response of an earth dam founded on liquefiable soils. In: Proceedings of the 2nd International Conference on Long Term Behavior of Dams, Graz University of Technology, Graz, Austria. 2009, 610-615

[21] Masing G. Eigenspannungen und Verfestigung Beim Messing. In: Proceedings of 2nd International Congress on Applied Mechanics, Zurich. 1926

[22] Itasca Consulting Group, Inc. FLAC (Fast Lagrangian Analysis of Continua). Version 4. Minneapolis, MN, 2004

[23] Kokusho T, Esashi Y. Cyclic tri-axial test on sands and coarse material. In: Proceedings of the 10th International Conference on Soil Mechanics and Foundation Engineering, 1981 (Quoted by Ishihara 1986)

[24] Hardin B O, Black W L. Vibration modulus of normally consolidated clay. Journal of the Soil Mechanics & Foundation Division, ASCE, 1968, 94(2): 353-369

[25] Seed H B, Wong R T, Idriss I M, Tokimatsu K. Moduli and damping factors for dynamic analyses of cohesionless soils. Journal of Geotechnical Engineering, ASCE, 1986, 112 (11): 1016-1032

[26] Vucetic M, Dobry R. Effects of the soil plasticity on cyclic response. Journal of Geotechnical Engineering, ASCE, 1991, 117(1): 89-107

[27] Lysmer J, Kuhlmeyer R L. Finite element method for infinite media. Journal of Engineering Mechanics, 1969, 95(EM4): 859-877

[28] Kuhlmeyer R L, Lysmer J. Finite element method accuracy for wave propagation problems. Journal of the Soil Mechanics and Foundations Division, 1973, 99(SM5): 421-427

[29] USCOLD (US Committee on Large Dams). Updated guidelines for selecting seismic parameters for dam projects. Committee on Large Dams, Denver, CO, 1999

[30] Finn W D L, Lee K W, Martin G R. An effective stress model for liquefaction. Journal of Geotechnical Engineering Division, ASCE, 1977, 103(6): 517-553

[31] Arulanandan K, Scott R F, eds. Verification of numerical procedures for the analysis of soil liquefaction problems. Vol. 1. In: Proceedings of the International Conference on the Verification of Numerical Procedures for the Analysis of Soil Liquefaction Problems, Davis, California, <day>17-20</day><month>October</month>1993. Rotterdam: Balkema A A, 1993

[32] Arulanandan K, Scott R F, eds. Verification of numerical procedures for the analysis of soil liquefaction problems. Vol. 2. In: Proceedings of the International Conference on the Verification of Numerical Procedures for the Analysis of Soil Liquefaction Problems, Davis, California, <day>17-20</day><month>October</month>1993. Rotterdam: Balkema A A, 1994

[33] Griffiths D V, Prevost J H. Two- and three-dimensional dynamic finite element analyses of the Long Valley dam. Geotechnique, 1988, 38(3): 367-388

[34] Lai S S, Seed H B. Dynamic Response of Long Valley dam in the Mammoth Lake Earthquake Series of May 25-27 1980. <patent>Report No. UCB/EERC-85/12</patent>. Earthquake Engineering Research Center, 1985

[35] Elgamal A M, Abdel-Ghaffar A M, Prevost J H. 2-D elastoplastic seismic shear response of earth dams: application. Journal of Engineering Mechanics, 1987, 113(5): 702-719

[36] Yiagos A N, Prevost J H. Tow-phase elasto-plastic seismic response of earth dams: applications. Soil Dynamics and Earthquake Engineering, 1991, 10(7): 371-381

[37] Zeghal M, Abdel-Ghaffar A M. Analysis of behavior of earth dam using strong-motion earthquake records. Journal of Geotechnical Engineering, 1992, 118(2): 266-277

[38] Woodward P K, Griffiths D V. Non-linear dynamic analysis of the Long Valley dam. Compute Methods and Advances in Geomechanics, 1996, 11(6): 635-644

[39] Ohmachi T, Kuwano J. Dynamic safety of earth and rock fill dams. Rotterdam: Balkema A A, 1994

[40] Ebrahimian B. Evaluating the seismic behavior (and considering the failure mechanism) of earth dams by using finite difference approach. Dissertation for the Master Degree. Isfahan: Department of Civil Engineering, Isfahan University of Technology, 2004

[41] Özkan M, Özyazicioglu M, Aksar U D. An evaluation of Güldürcek dam response during 6 June 2000 Orta earthquake. Soil Dynamics and Earthquake Engineering, 2006, 26(5): 405-419

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg