The Effects of Dam–Reservoir Interaction on the Nonlinear Seismic Response of Earth Dams

ABSTRACT The objective of this study is to investigate the effects of dam–reservoir interaction (DRI) on the nonlinear seismic response of earth dams. Although DRI effects have for long been considered as insignificant for earth dams, that conclusion was mainly based on linear elastic investigations which focused only on the acceleration response of the crest without examining the seismic shear stresses and strains within the dam body. The present study explores further the impact of DRI focusing on the nonlinear behavior of earth dams. The effects of reservoir hydrodynamic pressures are investigated in terms of both seismic dam accelerations and nonlinear dynamic soil behavior (seismic shear stresses and strains). It is shown that although dam crest accelerations are indeed insensitive to DRI, the stress and strain development within the dam body can be significantly underestimated if DRI is ignored.


Introduction
The dynamic behavior of a dam with a full reservoir is known to be different from that of a dam with an empty reservoir. The developed hydrodynamic pressures [Chopra, 1968;Chopra, 1995] affect the motion of the dam, while the dam response influences in turn the dynamic response of the reservoir. This phenomenon is termed dynamic dam-reservoir interaction (DRI) and it could be catastrophic in cases of resonance, that is when the two domains (dam and reservoir) are vibrating in phase.
There are numerous studies related to the numerical modeling of DRI in which the reservoir domain is discretized and included in the analysis. Several methods exist, including finite element (FE) [Demirel, 2015] and boundary element (BE) approaches or coupled BE-FE [Antes and Von-Estorff, 1987;Yazdchi et al., 1999]. The reservoir can be discretized using elastic "solid finite elements" [Wilson, 1975;Zienkiewicz et al., 1986;Wilson, 1995;Dakoulas and Gazetas, 2008;Pelecanos et al., 2013] in which the bulk modulus of the water is assigned to the elements with a very small shear modulus, or "fluid elements", following Eulerian [Kucukarslan et al., 2005;Gogoi and Maity, 2007;Fan and Li, 2008] or Lagrangian [Akkose et al., 2008;Bilici et al., 2009;Kartal and Bayraktar, 2013] or MPM [Kularathna and Soga, 2017a; La Villita earth dam Description of the dam La Villita is a 60 m high zoned earth dam in Mexico with a crest about 420 m long, founded on a 70 m thick alluvium layer. The dam cross-section is composed of a central clay core of very low permeability, with sand filters and rockfill shells. The alluvial layer beneath the clay core was grouted below the dam, while there is also a 0.6 m thick concrete cut-off wall to control seepage through the alluvium below the dam. Fig. 1(a) shows a schematic representation of the transverse cross section of the dam, whereas (b) shows a longitudinal section.
A summary of known material properties is given in Table 1. According to Elgamal [1992], the maximum shear stiffness, G max , for all the materials in the dam embankment varies and is about 140-260 MPa from top to bottom, whereas the foundation alluvium has a constant value of around 200 MPa.  The dam was built in 1967, the reservoir was impounded over 6 months in 1969 and it operated safely until the major seismic event in 1975. It subsequently experienced six more major seismic events between 1975 and 1985 (Table 2), which resulted in some permanent deformations, though not detrimental for the stability of the dam. The earthquake motions were recorded by three accelerometers installed on the dam (at the crest shown as C in Fig. 1) and the downstream berm, that is the small embankment in the downstream side of the dam (shown as B in Fig. 1) and the right rock bank (shown as R in Fig. 1). However, due to the instrument malfunction only EQ2 and EQ5 (Fig. 2) are of sufficient quality to be useful for numerical analysis [Elgamal, 1992] (see Pelecanos et al. [2015] for more details about the instrumentation). The seismic response of the dam was investigated by previous researchers, who were particularly interested in dynamic dam behavior [Elgamal, 1992;Pelecanos, 2013;Pelecanos et al., 2015], permanent displacements [Elgamal et al., 1990;Succarieh et al., 1993;Gazetas and Uddin, 1994;Uddin, 1997], and dam-canyon interaction [Papalou and Bielak, 2001;Papalou and Bielak, 2004].

Numerical model and calibration
In this study two-dimensional (2D) plane-strain static and dynamic-in-the-time-domain coupled-consolidation FE analyses, employing the Imperial College Finite Element Program (ICFEP) [Potts and Zdravković, 1999;Potts and Zdravković, 2001], were performed. The dynamic coupled hydro-mechanical (u-p) formulation of ICFEP is presented in Hardy [2003], Kontoe [2006], Kontoe et al. [2008c], Kontoe et al. [2008a], and Han  [2015]. It should be noted however that Pelecanos [2013] modeled the seismic response of La Villita dam using both coupled and uncoupled analyses and found very minor differences in the dam response. The full stress history of the dam prior to the earthquake events (including layered embankment construction, reservoir impoundment and consolidation) is modeled to establish a realistic starting point for the subsequent time-domain dynamic analyses. La Villita dam was previously analyzed by Pelecanos et al. [2015], incorporating both the static (layered construction, reservoir impounding, and operation-consolidation) and dynamic (seismic) behavior of the dam. That numerical model was compared with available field recorded data and focussed in particular on the dynamic stiffening effect of the narrow canyon geometry. That study did not consider DRI effects; instead it completely ignored the hydrodynamic pressures and only considered the hydrostatic reservoir pressure as an externally applied boundary stress (BS) on the upstream side of the dam. This current paper builds upon the already verified (against available field data) model from the previous study [Pelecanos et al., 2015] and focuses on the nonlinear dynamic DRI effects on earth dams during earthquakes, looking in particular at the acceleration response of the dam structure, together with the mobilized shear stresses and strains. The material constitutive model employed is a cyclic nonlinear elastic (CNL) model, which adopts a logarithmic function to describe the backbone curve of soil's monotonic response [Puzrin and Burland, 2000;Taborda, 2011], coupled with a Mohr-Coulomb yield criterion to capture plasticity [Kontoe et al., 2011]. The logarithmic relation dictates the degradation of shear stiffness, G, and the increase of damping, ξ, with cyclic shear strain, γ. The CNL model is able to reproduce hysteretic cyclic soil behavior, it depends on G max and requires 3 further distinct model parameters (E dL , J L , c-see Table 3). The backbone curve is given by the logarithmic relation in Eq. 1.
where J and E d are the deviatoric stress and strain invariants respectively defined by Eq. 2 and Eq. 3, whereas α and R are auxiliary constants defined by Eq. 4 and Eq. 5.
And χ L is defined by Eq. 6.
More details about this logarithmic model and its implementation in ICFEP may be found in Taborda [2011]. This model has been successfully used in static and dynamic-seismic analyses of soil-structure interaction problems [Kontoe et al., 2012;Kontoe et al., 2014;Han, 2014;Han et al., 2016;Han et al., 2017]. Due to the lack of experimental data, the CNL model is calibrated on empirical relations. The curves of Vucetic and Dobry [1991], Seed et al. [1986], and Rollins et al. [1998] have been used for the clay core, sand filters and rockfill-alluvium materials respectively. The parameters derived from the calibration are listed in Table 3 and the resulting stiffness degradation and damping development curves are shown in Fig. 3. The latter figures show that the employed model is able to reproduce the empirical stiffness and damping curves up to shear strains of 0.2-0.3%. More details about the calibration may be found in Pelecanos et al. [2015].
The FE mesh employed is shown in Fig. 4. The dam embankment is constructed in ten successive layers to simulate incremental construction. Similarly, reservoir impounding is simulated by constructing these elements (i.e. activating their weight) over ten increments.
In order to examine the effects of DRI, two approaches are followed to simulate the seismic reservoir hydrodynamic pressures. In the first approach, the hydrodynamic pressures are ignored and instead a constant hydrostatic variation of the water pressure, in the form of a BS, is applied along the upstream slope of the dam to model the reservoir pressures and at the same time, a hydrostatic variation of pore water pressures is specified in the upstream rockfill. This hydrostatic BS is kept constant for the entire static and dynamic analysis, thus ignoring DRI effects.
In the second approach, the reservoir domain is discretized with displacement-based finite elements . Its behavior is assumed to be linear elastic (using 8-noded iso-parametric displacement-based solid elements) and interface elements are placed along the interface between the reservoir and the dam and between the reservoir and the foundation alluvium. The material properties assigned for the reservoir are the bulk modulus of the water, K w = 2.2 10 6 kPa and a nominal value of shear modulus, G w = 100 kPa. The values for the shear and normal stiffnesses of the interface elements are K s = 1 kN/m 3 and K N = 10 8 kN/m 3 , respectively, the former prescribed as negligible in order to simulate the lack of shear transfer in water. Damping of the Rayleigh type is specified in the reservoir domain corresponding to a target damping, ξ t = 1% (ξ t < 2% was found to be an appropriate value by Pelecanos [2013] to avoid unrealistic reductions in the simulated hydrodynamic pressures). The two Rayleigh damping circular frequencies, ω 1 and ω 2 , are taken as the first and the third natural circular frequencies of the reservoir and they are listed in Table 4 along with the Rayleigh damping coefficients A and B (Eq. 7) [Bathe, 1996].
The boundary conditions at the vertical upstream reservoir boundary during the static analysis are the same as adopted on the lateral boundaries of the dam foundation, that is zero displacement in the horizontal and zero force in the vertical direction. To simulate the hydrostatic distribution of water pressure, an equal (hydrostatic) horizontal stress in applied as BC. For the dynamic analysis, the standard viscous boundary condition [Lysmer and Kuhlemeyer, 1969;Pelecanos et al., 2013] is applied normal to the vertical upstream reservoir boundary. Finally, vertical and horizontal displacement tied-degreesof-freedom (TDOF) are applied on the lateral foundation boundaries, which means that  nodes of the same elevation at the two lateral foundation sides are prescribed to have the same horizontal and vertical displacement. A summary of the static and dynamic BCs are shown graphically in Fig. 4 with black and grey colour respectively.

Numerical analysis of DRI
Dynamic response of the dam Figure 5 shows the acceleration time-histories at the crest of the dam for EQ2 and EQ5 respectively, whereas Fig. 6 shows the associated response spectra for the analyses considering the presence of the reservoir, shown as "Reservoir". The results of the analysis of Pelecanos et al. [2015] that ignored DRI and considered the reservoir water pressure simply as a constant hydrostatic "boundary stress", BS, are also included in the same graphs. Clearly, very minor differences exist between the analyses which either consider or ignore the DRI, with the "Reservoir" case yielding slightly smaller accelerations than the case with the hydrostatic BS (i.e. which ignores DRI). This is in agreement with earlier observations [Hall and Chopra, 1982a;Zhao et al., 1993;Guan and Moore, 1997;Pelecanos et al., 2016] who stated that seismic dam accelerations are not significantly  affected by the presence of the upstream dam reservoir. Additionally, the response spectra of the filtered acceleration records monitored during the two seismic events (EQ2 & EQ5) on the crest of La Villita dam are included in Fig. 6, showing a good comparison between the calculated and monitored acceleration values and thus confirming the validity of the adopted numerical model. Likewise, Fig. 7 presents a comparison of the calculated (for the BS approach) and recorded time-histories of horizontal acceleration and displacement at the crest of the dam during EQ5, which shows that a very good match is obtained between the calculated response from the numerical model and the actual response recorded in the field. Moreover, the vertical profiles of the maximum horizontal displacement and the corresponding maximum horizontal acceleration in the core for EQ2 and EQ5 are shown in Fig. 8 and Fig. 9, respectively. It should be clarified that the figures plot the "envelope" of the maximum horizontal displacement and acceleration at every point along the central vertical section within the core of the dam and not snapshots of displacement or acceleration. It may be observed that minor differences exist between the two analyses. For both EQ2 and EQ5, the analyses considering DRI resulted in slightly smaller values of acceleration at the dam crest, as also shown for the acceleration time-histories and the associated response spectra. Moreover, smaller values were calculated for the horizontal accelerations and displacements in the dam core (up to 10%) for the analyses considering the reservoir. Finally, it should be noted that the effects of DRI were more pronounced for the higher intensity ground motion, that is EQ5, which resulted in higher values of accelerations in the dam structure. It may be concluded here that the presence of the reservoir "damped" the response of the dam and resulted in slightly smaller values of accelerations. The effects of DRI were theoretically investigated by Pelecanos et al. [2016] for linear elastic dams founded on a rigid foundation. A direct comparison of this study with those results cannot be made as (a) Pelecanos et al. [2016] considered dams on a rigid base, whereas La Villita dam is founded on a compliant soil layer and (b) the fundamental period of the dam cannot be accurately estimated here, as nonlinear behavior was taken into consideration. However, a qualitative agreement between the two studies may be noted, as both studies show that (a) DRI results in some additional "damping" of the dam response (except in the rare case of dam-reservoir resonance) and (b) that DRI effects on dam accelerations are not significant in earth dams (difference in accelerations and displacements less than 10%). Fig. 10 shows the time-histories of the total hydrodynamic force, F dyn , on the upstream face of the dam during EQ2 and EQ5, whereas Fig. 11 shows the distribution of the peak hydrodynamic pressure, p dyn , on the upstream face of the dam for the two earthquakes. The values of F dyn were obtained by integrating the hydrodynamic pressures on the upstream face with respect to the vertical length of the face. It may be observed that significantly higher values of the hydrodynamic pressure (and hence the total hydrodynamic force on the upstream face) occur for EQ5, which is larger than EQ2 and therefore induces higher values of acceleration in the dam structure.

Hydrodynamic pressures
The absolute values of the maximum total hydrodynamic force on La Villita dam are listed in Table 5 along with the fraction with respect to the hydrostatic value, F st , calculated using Eq. 8 (where, γ w and h w are the unit weight of water and height of the reservoir respectively).
It is shown that the maximum total hydrodynamic force for EQ5 was found to be larger than that for EQ2. This was expected, as EQ5 induced higher values of acceleration in the dam structure and the reservoir domain and this is in agreement to earlier theoretical studies [Westergaard, 1933;Zangar, 1952;Chopra, 1967a;Liu, 1986]    maximum value of the induced accelerations. Moreover, the maximum value of the hydrodynamic pressure does not occur at the base of the dam, but rather at some elevation. This is again in agreement with the earlier observation of Zangar [Zangar, 1952;Zangar and Haefeli, 1952] who suggested that the maximum hydrodynamic pressure for dams with an inclined upstream face occurs at a vertical distance from the base approximately equal to one third of the reservoir height. The validity of the adopted reservoir modeling approach with solid finite elements has already been extensively verified  against available closed form solutions [Chopra, 1967a] and other numerical studies [Kucukarslan et al., 2005]. Therefore, since it is considered an appropriate approach to model hydrodynamic pressures, the numerically predicted hydrodynamic pressures are used herein as a benchmark against which the Zangar [1952] predictions are assessed. The relationship of Zangar [1952] has been widely used in engineering practice and research to provide an estimate of the "added mass" of the water in cases where the dam reservoir is not discretized [Dowdell and Fan, 2004;Yamaguchi et al., 2004;Leung et al., 2008;Huang, 2011;Dakoulas, 2012;Hellgren and Gasch, 2015]. Therefore, such a comparison allows assessment of the impact of the assumptions of Zangar [1952] on the accuracy of the simplified predictions. Consequently, Fig. 11 includes also a comparison of the computed peak pressure distribution (grey line) with the commonly employed analytical relationship of Zangar [1952] for the vertical profile of the hydrodynamic pressure, p dyn (z) (Eq. 9).
where, β is the angle of the upstream dam slope with the vertical direction, h w is the height of the reservoir, y is the vertical elevation from the bottom of the reservoir, and a c is the peak earthquake acceleration. The Zangar [1952] relationship (Eq. 9) estimates the peak hydrodynamic pressure distribution on a dam with a sloped upstream face assuming an incompressible reservoir. It is shown that the hydrodynamic pressures calculated in this study are larger than those obtained from Zangar [1952] for both earthquake events. This difference may be attributed to the following: • The analytical results of Zangar [1952] referred to rigid incompressible dams and therefore neglected any effects of reservoir-dam interaction, whereas La Villita dam is (and was analyzed as) a deformable earth dam. The hydrodynamic pressures on deformable dams were investigated analytically by Lee and Tsai [1991], who showed that, although deformable dams experience generally smaller values of hydrodynamic pressures, there might be DRI cases where large values of such hydrodynamic pressures may be generated, for example when dam-reservoir resonance occurs. However, a direct comparison with the work of Lee and Tsai [1991] is not possible, as they considered only very simplified geometries of elastic cantilever (modeled as thin uniform Euler-Bernoulli beam) dams. • The analytical results of Zangar [1952] were based on dams built on a rigid base, whereas La Villita dam is built on a compliant soil layer foundation. Dam-foundation interaction has been studied by previous researchers [Chopra and Perumalswami, 1971;Dakoulas and Gazetas, 1987;Papalou and Bielak, 2001;Papalou and Bielak, 2004] showing that in most cases it results in smaller dam accelerations, although for certain cases where resonance occurs it may result in amplification of accelerations. Therefore, in cases of dam-foundation resonance, larger dam accelerations will be expected and therefore larger hydrodynamic pressures. However, it is not clear whether there is resonance between dam and foundation in this case due to the nonlinear response of the soil and the complicated problem geometry. • The analytical results of Zangar [1952] assumed an incompressible reservoir (i.e. infinitely stiff), whereas the upstream reservoir was modeled as compressible (adopting a finite value of compressibility corresponding to the actual compressibility of water, see Numerical model and calibration). The effects of reservoir compressibility were initially investigated by Chopra [1967a], who showed that the assumption of incompressible reservoir in certain cases may considerably underestimate the hydrodynamic pressures.
Finally, it should be noted that although Zangar's [1952] simple analytical relation (Eq. 9) predicts the location of the maximum hydrodynamic pressure to be at the bottom of the dam, the experimental results in the same study showed that the maximum value occurs at a distance from the base of about one third of the reservoir height. This is in agreement with the numerically predicted profiles of hydrodynamic pressure (Fig. 11), which exhibit the location of the maximum hydrodynamic pressure to be at about h w /3 and h w /2 for EQ2 and EQ5 respectively. The experimental observations of Zangar [1952] were also later proved mathematically too [Chwang and Housner, 1978;Chwang, 1978;Liu, 1986]. Despite the existence of a plethora of subsequent and more advanced studies of hydrodynamic pressures on dams, there is no analytical expression that can precisely predict the hydrodynamic pressures on La Villita dam (with this geometry of the dam and the underlying soft soil foundation along with the dynamic seismic excitation) and therefore the fundamental results of Zangar [1952] are considered sufficient here for comparison. Although it has already been proven that hydrodynamic pressures do not seem to affect significantly the crest accelerations of earth dams [Hall and Chopra, 1982a;Pelecanos et al., 2016] (due to the small inertia of the reservoir compared to the inertia of a large embankment dam), it is argued that the observed high values of hydrodynamic pressures from severe earthquakes (e.g. here, for EQ5, F dyn /F st = 0.11) may induce localized effects on the upstream dam slope increasing the soil stresses which could potentially lead to localized soil yielding. Therefore, it is useful to investigate in more detail the dynamic soil behavior of the dam materials.

Dynamic soil behavior
As mentioned earlier, most studies of DRI for earth dams focused on evaluating the dam crest accelerations only. However, seismic wave propagation and reservoir hydrodynamic pressures can affect the dynamic behavior of the earth dam materials and therefore a more comprehensive investigation of the dynamic response within the dam body is needed in order to properly assess the effects of DRI on earth dams.
Considering EQ5 as the largest seismic event, Fig. 12 and Fig. 13 show the timehistories of shear strain at the first integration point of each of the four elements within the rockfill, both in the upstream (elements URL and URH, shown in Fig. 4) and downstream sides (elements DRL and DRH, shown in Fig. 4) for the two analyzed cases. The results show no major differences between the two analyses in the calculated shear strains in the downstream rockfill (Fig. 13). This was expected as the hydrodynamic pressures do not affect the downstream side of the dam. In contrast, considerably higher values of shear strain are reported from the analysis considering DRI for the upstream rockfill (Fig. 12), which also show some clear step-changes resulting in permanent plastic strains. This is more pronounced in element URH, suggesting that this element has experienced significant yield due to the additional lateral hydrodynamic loads from the reservoir. This observation is supported by the shear stress-shear strain plots for the upstream rockfill elements, shown in Fig. 14. The figure shows that although both elements experience shear stresses of comparable magnitude in both analyses, they exhibit higher shear strains in the DRI analysis case which results in permanent residual plastic strains. Again, of particular interest is element URH, which seems to engage plasticity a number of times (nine times of significant strain development can be clearly observed) as evidenced by the large development of strains. Plots of stress paths J-p' (where p' is the mean effective stress), for the two upstream rockfill elements considered, are shown in Fig. 15 and Fig. 16 for both analysis cases. These figures include the stress paths for the whole analysis, namely, the static and dynamic (both EQ2 and EQ5) parts. On the same figures, the Mohr-Coulomb yield surface is also plotted; the yield surface plotted here corresponds to the point that the yield surface was first engaged; that is, the first time that plastic strain developed at that point. As it may be observed from the figures, and especially for element URH, the stress paths for the DRI analysis case engage and travel along the yield surface a number of times during the seismic analysis. This results in significant plastic strain development as observed earlier in Fig. 12  upstream rockfill (it was reported that during EQ5 the reservoir load was increased by 11.22% due to the additional hydrodynamic component, see Table 5).

Comments
The results of this study, which considers DRI effects on the nonlinear response of earth dams, confirm the observations of previous researchers who suggested that dam crest accelerations are rather insensitive to reservoir hydrodynamic pressures. Moreover, when examining the peak horizontal accelerations and displacements within the dam body, it was found that DRI "damps" the dynamic response of the earth dam resulting in slightly smaller values of horizontal accelerations and displacements (generally less than 10%), perhaps due to the asynchronous vibration of the dam and reservoir domains. Nevertheless, significantly high values of reservoir hydrodynamic pressure may be experienced on the upstream dam face during earthquakes of large intensity. These hydrodynamic pressures may induce large stresses on the upstream dam slope, which in turn result in substantial strains due to localized yielding of the soil materials. The current study showed that some significant shear strains develop in the upstream dam slope potentially due to these additional reservoir hydrodynamic pressures. This also shows that examining the magnitude of dam crest accelerations and displacements only may not be representative of the actual soil behavior within the entire dam body and therefore might significantly underestimate the effects of DRI.
It is therefore highlighted that careful consideration of hydrodynamic pressures from the upstream reservoir is required when analyzing the seismic behavior of earth dams. The role of DRI may be more detrimental for dams prone to liquefaction. In those cases the additional load from the upstream reservoir may contribute to the generation of higher stresses leading potentially to flow liquefaction failures of dams. Finally, further parametric studies would be needed to extend this investigation, to examine, for example, the DRI effects for different dam material properties and to identify specific cases in which DRI effects are critical.
Comparing the results from the two seismic events considered (EQ2 and EQ5) it is found that larger values of earthquake accelerations result in larger effects on the dam embankment.  To obtain a more complete assessment and quantify the effects of DRI on the nonlinear response of earth dams, it is suggested that further and systematic parametric studies are performed which investigate factors like the influence of dam slope angle, dam material permeability and magnitude and frequency content of earthquake accelerations.

Conclusions
This paper describes an investigation carried out to assess the effects of the dam-reservoir interaction (DRI) on the nonlinear seismic behavior of dams. This study considers further the effects of DRI by investigating the response of a well-documented earth dam employing nonlinear finite element analysis. Two cases are considered for the upstream reservoir: (a) modeled simply as an external constant hydrostatic BS (i.e. ignoring hydrodynamic pressures) and (b) discretized with solid finite elements representing water (i.e. taking into account hydrodynamic pressures).
The main conclusions of this study may be summarized as follows: • The dam crest accelerations and associated response spectra are found to be insensitive to the presence of the upstream reservoir and the related hydrodynamic pressures. This is in agreement with previous studies which have observed essentially no major difference in crest accelerations when studying visco-elastic dams under harmonic and random seismic loads. • The profiles of horizontal displacements and accelerations in the dam core are also found not to be greatly sensitive to the reservoir hydrodynamic pressures. However, it was still found that the latter pressures result in smaller values of horizontal displacements and accelerations (in some cases up to 10% smaller). It is believed that the presence of the reservoir "damps" slightly the dynamic response of the dam, resulting in dam vibrations of smaller intensity due to the asynchronous vibration of dam and reservoir domains. • Although the predicted hydrodynamic pressures on the upstream dam face were found to be larger than those calculated from relevant analytical closed-form solutions, their magnitude was still comparable. This discrepancy can be attributed to the various simplifying assumptions adopted by the simple analytical solutions. There are no analytical solutions for complicated geometries, material behavior and loading, as in the examined case. It was also noticed that high values of hydrodynamic pressures (more than 10% of the hydrostatic values) may be observed for severe seismic events, which could potentially result in soil yielding on the upstream dam slope. • In contrast to the insensitive dam crest accelerations, significant differences were observed in the mobilized shear stresses and shear strains in the upstream rockfill close to the dam slope. Shear strains of considerable magnitude were predicted for the DRI case, which resulted from soil yielding. On the contrary, only minor differences were observed in the downstream rockfill, as this is not directly affected by the upstream reservoir hydrodynamic pressures. Therefore, ignoring DRI has severely underestimated the induced seismic strains.
In conclusion, it is believed that although the dam crest accelerations seem to be insensitive to the modelling of the upstream reservoir, the hydrodynamic pressures may impose some significant additional localized stresses on the upstream dam slope. These may result in soil yielding and associated large plastic strains. It is consequently recommended that the effects of the upstream reservoir on embankment dams should not be totally ignored, but instead, assessed in more detail taking account of realistic nonlinear soil behavior of earth dams.

Nomenclature
A Rayleigh damping coefficient a c earthquake acceleration B Rayleigh damping coefficient F dyn total reservoir hydrodynamic force F st total reservoir hydrostatic force G shear modulus g acceleration of gravity h w height of reservoir J deviatoric stress K N , K s normal and shear stiffness of interface elements K w bulk modulus of water p' mean effective stress S a spectral acceleration t time T fundamental period of vibration z elevation β slope of upstream dam slope with the vertical direction γ shear strain γ w unit weight of water θ Lode's angle ξ damping σ' 1 , σ' 2 , σ' 3 principal effective stresses τ shear stress ω 1 , ω 2 circular frequencies of vibration for Rayleigh damping

Acknowledgments
The work presented in this paper forms part of a PhD research of the first author, conducted at Imperial College London (2009)(2010)(2011)(2012)(2013) and funded by EPSRC. Their support is greatly appreciated.

Funding
This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC).