Non-Linear Numerical Analysis of Earthquake-Induced Deformation of Earth-Fill Dams

is used to promote the level of hysteretic damping during dynamic analysis. Masing rules are implemented into the constitutive model to precisely explain the non-linear soil response under general cyclic loading. The numerical model is then calibrated using centrifuge test data as well as field data. The field data are obtained in real measuring of Long Valley earth-fill dam subjected to the 1980 Mammoth Lake earthquake. The results of dynamic analysis, obtained in this study, are compared with the real measurements of Long Valley dam in terms of accelerations computed at the crest of dam in both time and frequency domains. The proposed numerical model can properly reproduce the overall seismic behaviours of earth-fill dams, their qualities and quantities, under earthquake loading conditions, confirmed by validation analyses. After validation, the effects of dam height, real earthquake loading, soil behaviour and strength of shell materials on the seismic response of earth-fill dams are evaluated through a comprehensive This book sheds lights on recent advances in Geotechnical Earthquake Engineering with special emphasis on soil liquefaction, soil-structure interaction, seismic safety of dams and underground monuments, mitigation strategies against landslide and fire whirlwind resulting from earthquakes and vibration of a layered rotating plant and Bryan's effect. The book contains sixteen chapters covering several interesting research topics written by researchers and experts from several countries. The research reported in this book is useful to graduate students and researchers working in the fields of structural and earthquake engineering. The book will also be of considerable help to civil engineers working on construction and repair of engineering structures, such as buildings, roads, dams and monuments.

sophistication in terms of proper problem formulation, characterization of material properties and modelling of soil stress-strain behaviour. In dynamic analysis methods using numerical simulation techniques, comprehensive analysis of earth-fill dam responses to dynamic loading is allowed. The development of geotechnical computation and numerical modelling offers interesting facilities for dam response analysis, considering complex issues such as soil non-linearity, evolution of pore water pressures and real earthquake records. In this regard, Prevost et al. (1985) presented 2D and 3D non-linear dynamic finite element (FE) analyses of an earth-fill dam, based on non-linear hysteretic analysis using multi-surface plasticity theory. Lacy & Prevost (1987) proposed a general and efficient numerical procedure for analyzing the seismic response of earth-fill dams. In their procedure, the dams were considered as non-linear two-phase systems. They outlined appropriate coupled dynamic field equations for the response of two-phase soil system. Abouseeda & Dakoulas (1998) studied the non-linear seismic behaviour in earth-fill dam-foundation interaction using boundary element (BE) and finite element (FE) methods. Chen & Harichandran (2001) studied the stochastic response of Santa Felicia earth-fill dam, in southern California, to spatially varying earthquake ground motion (SVEGM). They used SVEGM model in which the propagation of seismic waves is taken into account. Cascone & Rampello (2003) investigated the seismic stability of an earth-fill dam using decoupled displacement analysis. Ming & Li (2003) conducted a full coupled analysis of failure and considered the remediation of Lower San Fernando Dam. They used a critical state model, incorporating the concept of state dependent dilatancy for describing the soil behaviour over full loading ranges. Adalier & Sharp (2004) studied the seismic behaviour and remediation of an embankment on a liquefiable foundation. Papalou & Bielak (2004) studied the non-linear seismic response of earth-fill dams with canyon interaction. In their developed FE-based method, the dam was idealized as a shear beam and the surrounding medium as a halfspace. The dam's non-linearity was considered using multi-yield surface plasticity theory. Ebrahimian & Vafaeian (2005) considered the seismic response of earth-fill dams during earthquake using 2D full non-linear dynamic analysis. They used finite difference (FD) method and adopted the Mohr-Coulomb elastic-perfectly plastic constitutive model to describe the stress-strain relation of the soil. They focused on the seismic behaviour of very high earth-fill dams. Wang et al. (2006) presented the dynamic analyses in which a nonlinear, effective-stress-based soil model is employed. They used bounding surface hypoplasticity model for sand and implemented the model into a 2D finite difference (FD) program. The advantages of the proposed non-linear approach, conducted for a rock-fill dam, were illustrated by comparing the obtained results with those of equivalent linear approach. The model's capability was demonstrated by evaluating the seismic performance of an earth-fill dam.  carried out the transient dynamic time history FE simulations to investigate the performance of earth-fill dams under seismic excitation. Then, they studied different failure modes of earth-fill dams as the earthquake aftermath. Sica et al. (2008) studied the effect of loading history on the seismic response of earth-fill dams. They considered the static and seismic behaviours of a real case-history using coupled dynamic approach. The approach was solved numerically by FE method. Rampello et al. (2009) studied the response of an earth-fill dam to seismic loading using the displacementbased analysis and the FE effective-stress dynamic analysis. The FE analysis was carried out using a constitutive model which was capable to reproduce the soil non-linearity and calibrated against laboratory measurements. They also investigated the effects of assumed input motion and bedrock depth on the seismic response of earth-fill dam. Ebrahimian (2009) presented a numerical modelling of seismic behaviour of an earth-fill dam rested on www.intechopen.com liquefiable foundation. The numerical simulation was carried out using effective-stressbased, full coupled non-linear dynamic analysis. In this regard, Finn-Byrne model with extended Masing rules was used to model the pore water pressure generation in the liquefied soils. Ebrahimian (2011) investigated the dynamic behavior of earth dams by using a full non-linear dynamic finite difference analysis. The effects of input motion characteristics and dam reservoir condition on the dynamic response of earth dams were identified in this study. For this purpose, three real earthquake records with different levels and PGAs were used as the input motions.
In many parts of the world, the repetition of medium-strong intensity earthquake ground motions at brief intervals of time has been observed. The design philosophies for dams in seismic regions are based on multi-level design approaches, which take into account more than a single damageability limit state. According to these approaches, a sequence of seismic actions may produce important consequences on the dam safety. In fact, dams have been among the first structures that have been designed systematically against different earthquake levels. Since 1989, the ICOLD guidelines have introduced several levels of seismic loading, namely the Operating Basis Earthquake (OBE), Maximum Credible Earthquake (MCE), Maximum Design Earthquake (MDE) and Safety Evaluation Earthquake (SEE). However, the terms MDE or SEE are used as substitutes for the MCE. In order to analyze the behaviour of dams for specified levels of seismic hazard, several requirements should be considered. The seismic input and performance levels associated with serviceability, damage control, and collapse prevention are also defined. A thorough review about the different earthquake levels for dam design has been given in Wieland (2008). Amadio et al. (2003) analyzed the effects of repeated earthquake ground motions on the response of single-degree-of-freedom systems (SDOF) with non-linear behaviour. Accordingly, a comparison study was performed to investigate the effect of a single seismic event on the originally non-damaged system for different hysteretic models in terms of pseudo-acceleration response spectra and damage parameters. They showed that the elastic-perfect plastic system is the most vulnerable one under repeated earthquake ground motions. Moustafa & Takewaki (2010) modelled ground motions of multiple sequences that produce the maximum damage in the structure. It was shown that critical repeated acceleration sequences produce larger structural damage compared to single critical earthquakes. Afterwards, Moustafa (2011) developed a new framework to model the design earthquake loads for inelastic structures. New measures of the structure performance that were based on energy concepts and damage indices were introduced in his paper.
Concerning the seismic-resistant design of dams, several international standards have been developed by scientific communities in the past 25 years. However, only few countries have their own guidelines and regulations for seismic design of dams. Therefore, the ICOLD Bulletins and the local seismic building codes (e.g., Eurocode 8) are used as references. In brief, other famous international codes are USCOLD (United States Committee on Large Dams), US Army Corps of Engineers (USACE), ANCOLD, (Australian National Committee o n L a r g e D a m s ) , I I T K -G S D M A G u i d e l i n e s f o r S e i s m i c D e s i g n o f E a r t h D a m s a n d Embankments and Canadian Dam Association Guidelines for Dam Safety.
In this chapter, a numerical study of seismic behaviour of earth-fill dams overlaying bedrock subjected to real earthquake records is presented. For this purpose, full non-linear dynamic finite difference (FD) analysis is employed incorporating a simple elastic perfectly plastic constitutive model and Rayleigh damping. The former is used to describe the stress-strain response of the soil and the latter to increase the hysteretic damping level. The effect of non-linear soil behaviour is then considered in the analysis from the very beginning of earthquake loading. In order to precisely explain the soil response under general cyclic loading, Masing rules (Masing, 1926) are implemented into the constitutive model. Soil stiffness and hysteretic damping change with loading history. Firstly, the procedures of calibrating the constructed numerical models with centrifuge test data as well as real case history are presented and explained. Moreover, some important aspects of model calibration are discussed. Long Valley earth-fill dam, subjected to the 1980 Mammoth Lake earthquake, is analyzed for the real case history and the obtained numerical results are compared with the real ones, measured at the site in both time and frequency domains. The computed values show relative good agreements with the measured ones. It is shown that the Masing rules, combined with the simple elastic-plastic model, offer reasonable numerical predictions. A comprehensive parametric study is also carried out to identify the effects of dam height, input motion characteristics, soil behaviour and strength of shell materials on the seismic response of earth-fill dams. It is demonstrated that the fundamental aspects of seismic behaviour of earth-fill dams can accurately be captured by the current numerical procedure. It should be mentioned that this study does not consider the fluid-skeleton interaction, which may have significant effects on the seismic response of earth-fill dams.

Numerical modeling procedure
Here, numerical analysis is conducted using FLAC program, based on a continuum finite difference discretization applying Lagrangian approach (Itasca, 2004). Every derivative in the set of governing equations is directly replaced by algebraic expression written in terms of field variables (e.g., stress or displacement) at discrete point in space. Regarding dynamic analysis, explicit finite difference scheme is applied to solve the full equation of motion using the lumped grid point masses derived from the real density of surrounding zone. The calculation sequence first invokes the equations of motion for deriving 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 grid variables from known values which are fixed over the time step being executed ( Fig. 1). The equation of motion, in the simplest form, relates the acceleration ( du dt  ) of a mass (m) to the applied force (F) which may vary with time. Newton's law of motion for the massspring system is: In a continuous solid body, Eq.
(1) is generalized as follows: where, ρ = mass density; t = time; x j = components of coordinate vector; g i = components of gravitational acceleration (body forces); ij = components of stress tensor; i = components in a Cartesian coordinate frame.
For problem analysis, the strain rate tensor and rotation rate tensor, having the velocity gradients, are calculated by the following equations: where, ij e = components of strain rate; ij  = components of rotation rate; i u  = components of velocity.
The specific mechanical relationship is used in order to obtain the stress tensor as below: where, M = specific rule of behaviour;  = history parameters (based on the specific rules which may or may not exist).
The problem selected here is the simplified representation of typical earth-fill dam geometry. The dam section is a symmetric zone section with central clay core rested on bedrock,

Constitutive model
Mohr-Coulomb constitutive relation is used to model the behaviour of soil. The failure envelope for this model corresponds to a Mohr-Coulomb criterion (shear yield function) with tension cutoff (tensile yield function). Stress-strain relationship is linear elasticperfectly plastic. Linear behaviour is defined by elastic shear and bulk modulus. While, plastic behaviour is determined by the angle of internal friction and cohesion of the soil. Shear modulus of sandy shell materials is calculated from the following formula (Kokusho & Esashi, 1986): where, G max = maximum (small strain) shear modulus in kPa; e = void ratio; m   = mean effective confining stress in kPa; Poisson's ratio is considered as 0.3 for the shell materials.
Shear modulus of clayey core materials is calculated by the below formula (Hardin & Black, 1968): Poisson's ratio for the core materials is taken as 0.45.
Here, the basic elastic-perfectly plastic model is modified in order to better fit with the curves of shear modulus and damping ratio derived from the experimental data. This modified model can predict the seismic behaviour and the associated deformations of earthfill dams. Masing behaviour is implemented into FLAC via FISH subroutine (Itasca, 2004) in order to represent more accurately the non-linear stress-strain behaviour of soil that follows the actual stress-strain path during cyclic loading. Masing model consists of a backbone curve and several rules that describe the unload-reload behaviour of soil and the cyclic modulus degradation. Backbone curve can be constructed by the modulus reduction curves coupled with the small strain modulus (G max ). Unload-reload rules can similarly be formulated to reproduce the hysteretic damping values expected from the standard curves of damping ratio versus shear strain (e.g., Seed et al., 1986;Vucetic & Dobry, 1991). These formulations are described later.
In this study, shear modulus and damping ratio curves, proposed by Seed et al. (1986) for sandy soils and Vucetic & Dobry (1991) for clayey soils, are adopted as the references. The geotechnical properties of earth-fill dam materials, used in the analyses, are presented in Table 1.

Boundary conditions
The geotechnical problems can be idealized by assuming that the regions far from the area of interest extend to infinity. The unbounded theoretical models should be truncated to the manageable size by using the artificial boundaries for minimizing the computation time as well as avoiding the outwards propagating waves form reflecting to the model. The viscous boundary, developed by Lysmer & Kuhlemeyer (1969), is used in the current calculations. In this case, independent dashpots are used in the normal and shear directions at the model boundaries, as shown in Fig. 3. During the static analysis, the bottom boundary is fixed in both the horizontal and the vertical directions and the lateral boundaries only in the horizontal direction. In dynamic analysis, when the dam is laid on the foundation (and not on the bedrock), lateral boundaries are changed into free-field boundaries, available in the FLAC library, in order to eliminate the wave reflections from the truncated boundaries.

Element size
Numerical distortion of propagating wave can occur in dynamic analysis as a function of modelling condition. The numerical accuracy of wave transmission is affected by both the frequency content of input wave and the wave speed characteristics of system. Kuhlemeyer & Lysmer (1973) showed that for an accurate representation of wave transmission through the soil model, the spatial element size should be smaller than 1/10 to 1/8 of the wavelength associated with the highest frequency component of input wave i.e., where, λ = wave length associated with the highest frequency component that contains significant energy. Considering the above mentioned criteria, the element size is defined small enough to allow the seismic wave propagation throughout the analysis.

Damping
Material damping in soil is generally because of its viscosity, friction and plasticity development. Indeed, the role of damping in the numerical models is the reproduction of energy losses in the natural systems subjected to dynamic loads. The dynamic damping is provided in the model by Rayleigh damping option available in FLAC. Rayleigh damping was primarily used in analyzing the structures and elastic continua in order to damp the natural oscillation modes of system. Rayleigh damping R d =5% is used to compensate the energy dissipation in the media (Itasca, 2004). In the dynamic analysis incorporated plasticity constitutive models, considerable amount of energy dissipation can occur during the plastic flow. In the calculations of such cases, minimal percentage of hysteretic damping (e.g. 2%) is considered as well. The dam's natural frequency is determined as the Rayleigh damping parameter by Fourier analysis of its free response, as shown in Fig. 4. The fundamental frequency (f 1 = 1.71 Hz) of the dam with 40 m height is shown in this figure and those of dams with different heights are tabulated in Table 2.

Time step
The governing equations in time should be integrated incrementally for completing the numerical solution. The solution time step should be small enough in order to accurately define the applied dynamic loads and guarantee the stability and convergence of solution.
In this regard, time step is about 10 -6 second in the current FLAC model.

Input excitations
Selecting dynamic input motion is an important task in the seismic evaluation processes. In non-linear dynamic analysis, the expected earthquakes should be expressed as a set of ground motion time histories. For more correct evaluation, the input motions which offer an appropriate range of dam responses should be selected in the adaptable time history realizations. This procedure may be intractable due to the number of time-history realizations. However, in reality, the level of earthquake responses, probably achieved by physical system, is limited. Quantifying such responses demands good understanding of the seismic response of the system as well as the ground motion parameters that characterize the damage potential of seismic input (USCOLD, 1999). Different parameters can be employed to identify the severity and damage potential of a certain acceleration time history, assumed as the representative of expected earthquake ground motion; peak ground acceleration value (PGA) is of such kind. The use of this descriptor is intuitively natural since accelerations and resulting inertial forces are directly related by Newton's second law. However, there is no direct relation between PGA and structural response at the dominant natural frequencies of most typical dams. Moreover, large PGA values are not sufficient for generating response conditions which lead to significant damage. Despite these limitations, PGA is still the fundamental parameter used to judge the damage potential of certain acceleration histories. On the other hand, the seismic response of system is strongly affected by the frequency content of earthquake. Therefore, the better characterization of a given input motion can be achieved by using some forms of spectral representations. In particular, using Fourier amplitude spectrum is at the core of earthquake engineering practice. However, such characterizations do not provide direct description of the duration or time variation features of a given input motion.
Based on above, in this chapter, three different real acceleration time histories are selected from a database of earthquake records: Tabas, PGA=0.93g in MCE level; Naghan, PGA=0.72g in MDE level; San Fernando, PGA=0.21g in DBE level. In the dynamic analysis of dams, the scaled records are filtered to the maximum frequency of 10 Hz, transferred to the "inside" bedrock formation by standard de-convolution analysis and applied at the base of numerical model. The information of earthquake records are summarized in Table 3

Full non-linear dynamic analysis
Equivalent linear analysis is the common method used for evaluating the seismic behaviour of earth structures. In this approach, first, the responses are linearly analyzed using the initial values of damping ratio and shear modulus. Then, the new values of damping ratio and shear modulus are estimated, using maximum value of shear strain and laboratory curves. These values are used for redoing the analysis. This procedure is repeated several times until the material properties show no variation. Therefore, no non-linear effect is directly captured by this method as it assumes linearity during the solution process. Straindependent modulus and damping functions are considered roughly in order to approximate some effects of non-linearity (damping and material softening).
In the non-linear analysis, employed in this study, the non-linear stress-strain relationship is directly followed by each zone. Damping ratio and shear modulus of the materials are calculated automatically at different strain levels. The real behaviour of soil, under cyclic loading, is non-linear and hysteretic. Such behaviour can be simulated by Masing model (Masing, 1926), which can model the dynamic behaviour of soil. In this model, the shear behaviour of soil may be explained by a backbone curve as: where, F bb(γ) = backbone or skeleton function;  = shear strain amplitude; G max = initial shear modulus; max = maximum shear stress amplitude.
Stress-strain curve follows the backbone curve in the first loading, as shown in Fig. 7(a); however, for explaining the unload-reload process, the above equation should be modified. If load reversal occurs at the point ( r , γ r ), stress-strain curve follows the path given by the below formula: In other words, the shapes of unload-reload curves are similar to that of backbone curve (with the origin shifted to the loading reversal point) except they are enlarged by a factor of 2, as shown in Fig. 7(b). The Eqs. (9) & (10) describe the Masing behaviour (Masing, 1926).
Masing rules seem not to be enough for precise explanation of soil response under general cyclic loading. Finn et al. (1977) developed modified rules to describe the irregular loading. They suggested that unloading and reloading curves follow the concerning two rules. If the new unloading or reloading curve exceeds the last maximum strain and cut the backbone curve, it will follow the backbone curve up to meeting the next returning point, as shown in Fig. 7(c). If a new unloading or reloading curve passes through the previous one, it will follow the former stress-strain curve, as shown in Fig. 7(d). According to this model, the tangent shear modulus can be defined at the points on the backbone and new reloadingunloading curves by the Formulas (11) & (12), respectively, as: Based on the results, obtained in this research, the shear stress decreases as the number of load cycles increases; it means that shear stress-strain curves are more inclined. In this study, Masing rules are implemented into FLAC via a series of FISH functions in order to simulate the non-linear stress-strain relationships.

Validation analysis
In this research, one-zone sample is simulated using the unit cell as shown in Fig. 8, in order to validate the implementation of Masing rules in FLAC program. The one-zone sample consists of sandy soil and a periodic motion is exerted at its base. Vertical loading is established only by gravity and then, the Equilibrium stresses are installed in the soil. The stress/strain loops of the sample are shown in Fig. 9(a) for several cycles. According to the figure, shear modulus decreases as shear strain increases. It seems that the hysteretic model can reasonably handle the multiple nested loops. Energy dissipation and shear stiffness degradation are clearly observed during seismic loading, as shown in Fig. 9(b). Shear modulus reduction curve, obtained in this study, follows well the empirical relation proposed by Seed et al. (1986) and the test data. The results obtained from the numerical analyses are compared with those of experimental ones in order to evaluate the capability of proposed model. One of the centrifuge tests related to the embankment which was performed in VELACS project (VErification of Liquefaction Analysis using Centrifuge Studies Scott, 1993, 1994)) is chosen as a reference. It is attempted to create almost similar conditions for both laboratory model test and numerical model, as shown in Figs Finally, the model's ability in simulating the seismic behaviour of earth-fill dam during a real earthquake is verified by a real well-documented case history. In this regard, the results of non-linear dynamic analysis of Long Valley (LV) earth-fill dam in California subjected to the 1980 Mammoth Lake earthquake (Griffiths and Prevost, 1988) are presented. Then, the results are compared with the real measurements recorded at the site and also with the results presented by previous authors. LV dam is located in Mammoth Lake area www.intechopen.com (California) in close proximity of active faults. The dam is a rolled earth-fill one formed mainly with the impervious zone. The dam has maximum height of 55 m, 182 m length 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; it has 3 accelerometer stations to monitor the boundary conditions, and 5 stations to record the dam response ( Fig. 12(a)). 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 Lakes area. The magnitudes of these earthquakes were M L = 4.9 -6.7, and the induced peak accelerations at the crest centre was 0.5g in the upstream-downstream direction (x direction, as shown in Fig. 13(a)) during the strongest event. Extensive arrays of 22 input-output (excitationresponse) accelerations were recorded, providing a valuable information source of the dam seismic responses 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 seconds of the recorded acceleration is used with data point at 0.02 second intervals and the peak acceleration is 0.135g in the upstream-downstream (x) direction and 0.084g in the vertical direction (y).
The cross section of LV dam is shown in Fig. 12(b) and its detailed information is found in (Griffith and Prevost, 1988). The numerical grid constructed in FLAC is presented in Fig.  12(c). The input accelerations are applied in the horizontal and vertical directions of the model base. Free Field boundary conditions are exerted to the lateral boundaries of numerical model. This research focuses on the computed acceleration at the crest which can be compared directly with the measured values. Previously, LV dam has been analyzed by : Lai & Seed, 1985;Elgamal et al., 1987;Griffiths & Prevost, 1988;Yiagos & Prevost, 1991;Zeghal & Abdel-Ghaffar, 1992;Woodward & Griffiths, 1996 Fig. 13(a) shows the computed horizontal acceleration of the crest; it indicates that the amplification occurs between the base and the crest. The magnification factor of peak amplitude at the crest is about 5.47 over the peak base amplitude. The crest response, computed in the horizontal direction, is compared with the measured values, as shown in Fig. 13(b); the dashed line corresponds to the computed response there. Excellent overall agreement is achieved between the computed and measured values; however, the computed values show higher amplitudes. The frequency contents of two time records are compared in the form of Fourier amplitude spectra (FAS), as shown in Fig. 13(c). Their 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. 13(c), shows that the energy is concentrated just at the frequencies below 2 Hz.
In the vertical direction, the calculated acceleration shows low agreement with the measured values, as shown in Fig. 14(a). According to this figure, the plots of vertical acceleration are superimposed at the base and crest. This excitation is considerably noisier and less intensive in the vertical direction in compared with that of horizontal one. The maximum accelerations at the crest, recorded in the vertical and horizontal directions are 0.172g and 0.64g, respectively. The computed accelerations in the vertical direction are compared with the measured values in the crest of LV dam, as shown in Fig. 14(b). According to this figure, the computed values have generally lower amplitudes in compared with those of measured values. The Fourier amplitude spectra of these time histories are given in Fig. 14 The results, obtained in the validation analysis of LV dam, in term of crest acceleration are given in Table 5 and compared with the other numerical results presented by previous authors. According to the comparisons, the numerical procedure, presented in this study, can properly capture the fundamental aspects of seismic behaviours of earth-fill dams. As mentioned earlier, the numerical model is then used for parametric studying of hypothetical earth-fill dams due to the satisfactory modelling of validation cases.

Parametric study
Here, the analyses are carried out to investigate the effects of dam height, input motion characteristics, soil behaviour and strength of shell materials on the seismic behaviour of earth-fill dams. The effects of different earthquakes are studied on the horizontal permanent deformations, permanent shear strains and maximum accelerations, as shown in Fig. 15. The values have been induced at the crests of dams with different heights. The displacements are shown in Fig. 15(a) and the relevant shear strains in Fig. 15(b). It is clear that the shear strain variation is similar to displacement. The horizontal displacements and shear strains in the dam body increase with dam height increasing. The calculated values are much higher in Tabas earthquake and the failure occurs in the dam body. According to Fig. 15(a), the maximum horizontal displacement computed at the crest of dam is about 94 cm at the end of Tabas earthquake. It can be observed in Figs. 15 (a) & (b), that increasing in the input motion energy leads to significant increase of displacements and shear strains. Fig. 15(c) illustrates the coupled effects of dam height and earthquake type on the maximum acceleration induced at the dam crest. According to the figure, the crest acceleration decreases as the dam height increases and no amplification is seen maybe due to more flexible behaviour, larger damping and larger developed plastic zones, observed in higher dams. Therefore, because of these factors, more energy is absorbed in higher dams in compared with that in the shorter ones. It can be seen in Fig. 15 (c) that the accelerations in the dam crest are more reduced in higher dams comparing with the smaller ones. It should be mentioned that PGA of Naghan earthquake (0.72g) is much higher than that of San Fernando earthquake (0.21g). However, the created displacements and shear strains in the dam crest caused by Naghan earthquake are close to those of San Fernando input motion. It can be concluded that using just PGA parameter is not sufficient for evaluating the effect of a certain earthquake time history on the dam response. Therefore, other earthquake parameters such as effective duration, magnitude and frequency content should be considered in the analysis. Failure mechanism with permanent shear strain contour in the dam body is shown in Fig.  16, regarding two different heights at the end of Naghan earthquake. The slip surface is much deeper and more obvious in the dam with 280 m height ( Fig. 16(b)) in compared with that of 120 m height ( Fig. 16(a)). A dam with 40 m height is subjected to the mentioned earthquakes and chosen as a reference with two different behaviours, elastic and elastic-perfectly plastic, in order to investigate the effect of soil behaviour on the seismic response of dam body. As it is expected, regarding linear elastic behaviour, smaller displacements and shear strains are observed along the dam height, Figs. 17(a) & (b). However, large amplification occurs especially for the strongest earthquake, Fig. 17(c). It means that plasticity reproduce more energy dissipation during dynamic loading. In such cases, the accelerations are reduced across the dam height and therefore become lower than the base acceleration. According to Fig. 17(a), maximum displacement occurs at about Z/H = 0.88 in linear elastic behaviour, while it happens at the crest of dam in linear elastic-perfectly plastic behaviour. Furthermore, in elastic-perfectly plastic behaviour, the dynamic induced residual (permanent) displacement increases largely in the upper part of dam, especially for Tabas and Naghan earthquakes, confirmed in the previous research works (Ohmachi and Kuwano, 1994;Ozkan et al., 2006). That is why the crest should especially be considered in designing the embankment dams, due to the stronger shaking at the upper parts, for avoiding undesirable deformations. The distribution of shear strain is extremely non-linear along the dam height in the stronger earthquakes, as shown in Fig. 17(b). In the elastic dams, maximum acceleration occurs in the dam crest, as shown in Fig. 17(c). In the acceleration profile of Tabas earthquake, a special increase is seen along the dam centerline at Z/H = 0.38. Strength of dam materials is an important parameter which can significantly affect the seismic response of dam. In this regard, different friction angles are assumed accompanying with different dam heights, subjected to Naghan earthquake, in order to clarify the above mentioned effect. Fig. 18(a) shows horizontal displacement values versus dam height for different friction angles of shell materials. The variations of shear strains in the crests of dams, with different heights, are shown in Fig. 18(b) for different friction angles of shell materials. As it is expected when the friction angle increases, the horizontal displacement and shear strain in the dam crest decrease. It can be seen in the above figure that the variation of friction angle causes no significant change in the displacement and shear strain regarding ≤ 40˚. However, the variation is more significant in = 45˚, compared with the lower friction angles; the highest displacement values correspond to = 30˚. The horizontal displacements computed at the crests of dams with 40 and 120 m heights are about 16 and 13 cm, respectively, and their shear strains are about ‫³-5.3‬ and ‫,³-5.2‬ respectively. The maximum acceleration induced at the top of dam decreases as the friction angle decreases or the dam height increases, as shown in Fig. 18(c). Considering larger friction angles for the shell materials (e.g., = 45˚) leads to about 70% increase in the dynamic amplification. The computed maximum crest acceleration of dam with 40 m height is about 0.89g for = 45˚, while that of 120 m height is 0.52g. When the dam height deceases, the horizontal displacement and shear strain increase but the acceleration decreases at the crest of dam. All variations are linear for = 45˚, but for the other friction angles are slightly non-linear, as shown in Fig. 18.

Lessons learned
The author experienced several interesting points and noteworthy items during numerical model calibration and numerous dynamic analyses which are listed below:  In full non-linear dynamic analysis, soil stiffness degradation is automatically taken into account upon constitutive model of soil and just the initial shear modulus is needed as an input parameter. Therefore, it is important to be sure that in the numerical model, the trend of shear modulus decrease and damping ratio increase are in agreement with those of laboratory test results during dynamic loading.  The poisson's ratio of about 0.5 should not be used in the calculating of bulk modulus of undrained soil layers (such as clay) in the analysis. Otherwise, bulk modulus increases irrationally and the time step of analysis decreases rigorously and consequently the calculation time increases excessively. Therefore, the poisson's ratio should not be more than 0.45 in such cases.  If a "raw" acceleration record from a site is used as a time history, then FLAC model may exhibit residual displacements once the motion is finished. This arises from the fact that the integral of complete time history may not be zero. Therefore, the process of baseline drift correction should be performed in such cases.  The input motion should be filtered before being applied to the FLAC grid in order to eliminate all high frequency components form it.  The stages of construction should be considered in the numerical simulation of earth-fill dams. In the present study, the stages are: initial state of foundation (if any); layer by layer dam replacement; applying the hydrostatic water pressure due to the replacement of dam reservoir; seepage analysis in the dam body; mechanical adjustment to new flow field; and finally dynamic analysis. Regarding the mentioned stages, one is run to equilibrium and then the next stage is started. However, construction sequences have much greater effects on static results than dynamic results.

Conclusions
This chapter presents the non-linear seismic behaviour of earth-fill dams using explicit finite difference method. In this regard, a simple elastic-perfectly plastic constitutive model with Mohr-Coulomb failure criterion is used to describe the stress-strain response of the soil.
Here, Rayleigh damping is used to promote the level of hysteretic damping during dynamic analysis. Masing rules are implemented into the constitutive model to precisely explain the non-linear soil response under general cyclic loading. The numerical model is then calibrated using centrifuge test data as well as field data. The field data are obtained in real measuring of Long Valley earth-fill dam subjected to the 1980 Mammoth Lake earthquake. The results of dynamic analysis, obtained in this study, are compared with the real measurements of Long Valley dam in terms of accelerations computed at the crest of dam in both time and frequency domains. The proposed numerical model can properly reproduce the overall seismic behaviours of earth-fill dams, their qualities and quantities, under earthquake loading conditions, confirmed by validation analyses. After validation, the effects of dam height, real earthquake loading, soil behaviour and strength of shell materials on the seismic response of earth-fill dams are evaluated through a comprehensive parametric study. The effect of dam height on the non-linear seismic behaviour is particularly focused in this research. The following conclusions are obtained based on the performed parametric study:  If the dam materials keep their elastic behaviours during dynamic loading, then the horizontal acceleration increases along the dam height (from the base to the top). In this case, the higher dams show larger amplifications, especially if the natural periods of their bodies coincide with the periodical nature of earthquake waves.  When the dam body shows non-linearity or the materials go towards plastic behaviour during a strong shaking, the attenuation of acceleration waves in the dam body becomes more effective. Consequently, the amplitudes of earthquake accelerations decrease when moving from the base towards the top.  According to the non-linear elastic-plastic analyses, when the height of the dam increases, then the strongest dynamic loading (Tabas earthquake) induces plasticity in large parts of the dam body. In fact, strong earthquakes are more effective in changing the material behaviour from elastic to plastic condition in comparison with weak earthquakes.  The higher dams are more flexible than the smaller ones. Consequently, the flexibility affects the shear strains which influence the shear modulus degradation and attenuating coefficient. All these effects are on the trend of weakening the accelerations along the height.  Soils with less strength (suppose low friction angle) go towards yielding by small amount of dynamic force which cause the attenuation of acceleration along the dam height in the weak materials compared with the strong ones.  Regarding a dam subjected to the earthquake with lower energy, the dam body behaves as an elastic material. Therefore, the induced seismic accelerations inside the dam body become larger from the base of dam to its top. In this case, small plasticity zones are developed in the dam body and the dam remains safe during dynamic loading.  Finally, non-linear dynamic analysis shows that plasticity should be considered in the investigation of seismic response of earth-fill dams, because of which the acceleration of the dam crest decreases and the displacements and shear strains of dam body as well as the energy dissipation increase. All these can significantly affect the seismic response of earth-fill dams.