Numerical analysis of earthquake response of an ultra-high earth-rockfill dam

Abstract. Failure of high earth dams under earthquake may cause disastrous economic damage and loss of lives. It is necessary to conduct seismic safety assessment, and numerical analysis is an effective way. Solid-fluid interaction has a significant influence on the dynamic responses of geotechnical materials, which should be considered in the seismic analysis of earth dams. The initial stress field needed for dynamic computation is often obtained from postulation, without considering the effects of early construction and reservoir impounding. In this study, coupled static analyses are conducted to simulate the construction and impounding of an ultra-high earth rockfill dam in China. Then based on the initial static stress field, dynamic response of the dam is studied with fully coupled nonlinear method. Results show that excess pore water pressure accumulates gradually with earthquake and the maximum value occurs at the bottom of core. Acceleration amplification reaches the maximum at the crest as a result of whiplash effect. Horizontal and vertical permanent displacements both reach the maximum values at the dam crest.


Introduction
With the development of construction technology and a better understanding of the earth-rockfill dams' mechanism (Yasuda and Matsumoto, 1993;Hunter and Fell, 2003;Soroush and Jannatiaghdam, 2012), high earth-rockfill dams have become a more competitive type of dams.A number of high core earth dams and concrete faced rockfill dams (CFRDs) with heights of more than 200 ∼ 300 m are being under construction or to be built in Southwest China located in the Himalayas seismic belt with heavy potentially seismic activity.
Failures of earth dams, especially for high earth dams with a height more than 200 m, may cause catastrophic economic damage and loss of lives.Among these failures, earthquake attack is one of the major causes.In general, the principal failures induced by earthquake attack can be classified as sliding failure, liquefaction failure, longitudinal cracks, transverse cracks, and piping failure (Chen, 1995).Numerical computation based on constitutive model of geotechnical materials is one of the effective ways to conduct seismic analyses and seismic safety assessment.
Earth dams are complicated nonhomogeneous geotechnical structures composed of different materials such as rockfill and clayey soils.This results in great difficulty to analyse and predict the responses of earth dams under various conditions, such as construction, reservoir impounding and earthquake attacks (Costa and Alonso, 2009;Seo and Ha, 2009).Therefore, existing theories and numerical analyses of earth-rockfill dams include many simplifying assumptions (Prevost et al., 1985).
Pseudo static method (Seed and Martin, 1966) is largely used in the numerical analyses to assess the seismic response of earth dams.In this approach, seismic effect on dam body is considered with an equivalent static force, which is the product of soil mass multiplied by a seismic coefficient, and soil is assumed as rigid-perfectly plastic behaving.
Newmark (Newmark, 1965) presented the important concept that the seismic stability of slopes is more important rather than the traditional factor of safety, and he proposed a double integration method to compute the permanent displacement under cyclic loading, which is based on the concept of analyses of a rigid block sliding on a plane.Once the acceleration exceeds the "yield acceleration", permanent displacement can be obtained by double integrating the acceleration.
Up to now, most numerical earthquake response analyses of earth-rockfill dams were conducted using equivalent linear models, in which shear modulus and damping ratio are obtained by iterative procedure (Abdel-Ghaffar and Scott, 1979;Mejia et al., 1982).Cascone and Rampello (2003)  adopted to describe the behaviour of both shells and core of the earth dam.Their analyses showed that plasticity should be considered in the analysis of the seismic response of the dam, as it leads to a decrease in the natural frequency of the dam, which could significantly affect the seismic response of the dam.However, this study did not consider the fluid-skeleton interaction, which should have a significant influence on the seismic response of earth dams.Siyahi and Arslan (2008) studied the dynamic behaviour and earthquake resistance of Alibey earth dam located in Turkey by adopting the OpenSees software (Open System for Earthquake Engineering Simulation) (2000).Elastoplastic constitutive models of soils such as pressure dependent and independent multi-yield materials were implemented in the analyses.They obtained the permanent displacements and concluded that the liquefaction effect is not significant.
However, so far, application of solid-fluid coupled elasto-plastic dynamic analyses in safety assessment of earth-rockfill dams is not much.Furthermore, the initial stress field needed for dynamic computation is often provided without considering the effects of early construction and reservoir impounding.
In this work, the static and dynamic responses of Nuozhadu high earth-rockfill dam have been studied by solid-fluid coupled nonlinear method, since solid-fluid interaction has a significant influence on the behaviour of earth dams.First, the dam site and detailed design data are introduced.The adopted numerical analysis method and soil constitutive model parameters are presented.
In the static analysis, the practical construction was simulated.Then, reservoir impounding was simulated by conducting solid-fluid coupled analysis.
At last, the coupled nonlinear dynamic analyses were conducted based on the initial stress field obtained from the preceding computation.The acceleration component of the El Centro 1940 earthquake was chosen as seismic motion.The excess pore water pressure, acceleration responding, and permanent displacements induced by seismic move are presented and analysed.

Nuozhadu dam
Nuozhadu hydropower station is located in the Lancang River which is also named Mekong River in the downstream, in Yunnan Province, China, as shown in Fig. 1.The installed capacity of the power-station is 5850 MW.The most important part of Nuozhadu project is an ultra-high earth-rockfill dam with a height of 261.5 m, which is the highest earth-rockfill dam in China and the fourth highest in the world.The reservoir has a storage capacity of 237.0×10 8 m 3 , with the normal storage water level of 812.5 m (Fig. 2b).
Figure 2 shows the material zoning and construction stages of the maximum cross section.The elevation of the earth core bottom and the crest of the dam are 562.6 m and 824.1 m, respectively.The dam crest has a longitudinal length of 630 m with a width of 18 m.The upstream and downstream slopes are at 1.9 : 1 and 1.8 : 1, respectively.The dam body is composed of several different types of materials.The shells of upstream and downstream are composed of decomposed rock materials.Anti-seepage material in the earth core is clay mixed with gravel.Adding gravel to the clay can improve the strength of core and reduce the arching effect between shells and earth core.
The gravel material consists of fresh crushed stone of breccia and granite with a maximum diameter of 150 mm.In addition to these, fine rockfill and filter materials are filled against the earth core to prevent the fine particle from being washed away.
The dam construction was started in 2008, and was completed at the end of 2012. Figure 1c shows the construction site on the dam crest.Figure 2b demonstrates the practical construction process.

Numerical analyses method
Since solid-fluid interaction plays a significant role in the behaviour of earth dams during impounding and earthquake, coupled method was adopted in the work.The solid-fluid interaction consists of two mechanical effects.First, changes in pore pressure will affect the response of the solid according to the principle of effective stress.Second, changes in mechanical volume will also affect the pore pressure.
Governing differential equations to describe the interaction between soil and water include constitutive laws, balance law, transport law, and compatibility equations below (Itasca Consulting Group, 2005).
The constitutive response for the porous solid has the form where σ i j is the co-rotational stress rate, p is the pore pressure, H is the functional form of the constitutive law, κ is a history parameter, δ i j is the Kronecher delta, and ξ i j is the strain rate.
The constitutive law for fluid content is related to change in pore pressure, p, saturation, s, and mechanical volumetric strains, ε, as where M is Biot modulus, n is the porosity, α is Biot coefficient and ζ is the variation of fluid content.The fluid mass balance maybe expressed as where q v is the volumetric fluid source intensity.The fluid transport is described by Darcy's law.For a homogeneous, isotropic solid and constant fluid density, this law is given in the form where q i is the specific discharge vector, k is the tensor of absolute mobility coefficient of the medium, k(s) is the relative mobility coefficient, ρ f is the fluid density, and g j is the component of the gravity vector.
And the relation between strain rate and velocity gradient is (5) The Mohr-Coulomb model is the conventional model used to represent shear failure in soils and rocks.The failure criterion used in the constitutive model is a composite Mohr-Coulomb criterion with tension cut-off, which can be represented in the σ 1 − σ 3 plane, as illustrated in Fig. 3.The failure envelope f (σ 1 , σ 3 ) = 0 is defined from point A to B by the Mohr-Coulomb failure criterion f s = 0 which can be derived from the classic form τ f = σ tan Φ + c and from B to C by a tension failure criterion f t = 0 where φ is the friction angle, c, the cohesion, σ t , the tensile strength, and The potential function is described by means of two functions, g s and g t , used to define shear plastic flow and tensile plastic flow, respectively.The function g s corresponds to a non-associated law and has the form where ψ is the dilation angle and The function g t corresponds to an associated flow rule and is written as

Numerical model and parameters
The numerical analyses were performed by the effective stress analysis to simulate the performance of the dam during construction stage and under earthquake action.
Figure 4 shows the mesh of the maximum cross section of the dam according to the material zoning and construction stage (see Fig. 2).There are 648 quadrilateral and triangular elements and 682 nodes.
Geotechnical properties used in the analyses are presented in Table 1 for earth dam materials.The friction angle and cohesion were determined through conventional triaxial test.The Young's moduli were chosen more close to reality.
Numerical analyses presented in this study were performed by using the finite differential code FLAC3D, which is based on a continuum finite difference discretization using the Fast Lagrangian Approach.Seismic loading was input at the model base as a stress history excitation.In the study, El Centro earthquake was chosen as the horizontal acceleration and vertical acceleration by multiplying a factor of 2/3.The maximum acceleration occurs at 2.3 s with a value of 0.36 g.Velocity and displacement histories were obtained by integrating the acceleration history (see Fig. 5).The local damping was applied in the dynamic simulation with a value of 0.05.
The velocity history is transformed into stress history using the formulas as where ρ is soil mass density, σ n and σ s are applied normal stress and shear stress, v n and v s are input normal and shear velocities, C p and C s are the speeds of p wave and s wave propagating through medium given by Viscous lateral boundary was applied to the model in order to prevent spurious reflections of outward propagating waves from artificial boundaries.Free field boundary was set in parallel with the main-grid.
In order to conduct dynamic analysis, an initial stress field is needed.In this study, Mohr-Coulomb model was used to provide initial stress field.Since the real initial stress field has a significant influence on the seismic behaviour of the dam, practical construction process was simulated according to the design (as shown in Fig. 2).

Construction stage
In the static analysis, staged construction was simulated by applying the null model in FLAC3D.Forty five stages were divided according to the practical construction from 2008 to 2012. Figure 6 shows the evolution of vertical displacement during construction stage, from which we can see that the maximum vertical displacement occurred in the middle of earth core.Figure 7 presents the displacement and stress distributions when the dam construction is completed.The largest settlement occurred in the middle of earth core with a value of about 3.72 m, which is slightly larger than the monitoring value of 3.55 m due to the three dimensional effect.The monitoring value was obtained from the electromagnetic settlement rings buried in the middle of earth core during construction.Large horizontal displacements concentrated in the upstream and downstream shells with a maximum value of about 0.35 m.Because of the big differences of moduli between rockfill materials and clay, there exists obvious arching effect in the core wall, as shown in Fig. 7c.

Reservoir impounding
Then, reservoir impounding was simulated by conducting solid-fluid coupled analyses.The normal storage level of 812.5 m was applied on the upstream slope, as shown in Fig. 8. Earth core and dam shells with filter materials were set with different permeability coefficients, respectively (as shown in Table 1).The result shows that contour line of pore pressure dropped sharply in the earth core, which proves the effectiveness of the earth core in anti-seepage.
Figure 9 gives the displacements distribution after reservoir impounding.Compared with the analyses results of construction completion, horizontal displacement developed toward downstream due to the water pressure applied on the upstream slope.
As a result of reservoir storage, soils suffered an upward buoyant force, which resulted in the decrease of settlement without considering the wetting deformation of rockfill materials.

Dynamic analysis
Coupled solid-fluid dynamic analysis was conducted based on the initial stress field obtained from the preceding computation.Without empirical formulas for computing permanent displacements and accumulated excess pore water pressure like traditional equivalent linear method (Seed, 1966;Finn et al., 1977), nonlinear coupled dynamic analysis can give those results directly.

Seismic wave input validation
The N-S component of the El Centro 1940 earthquake was chosen as ground motion with a 20 s duration time, as given in Fig. 5.The earthquake wave was applied to the dam in the upstream-downstream direction, and in the vertical direction multiplied by 2/3.To verify the correctness of seismic wave input, acceleration time history was recorded at the model base, as shown in Fig. 10a.Furthermore, acceleration response spectral were obtained and compared in Fig. 10b.As shown in the figure, we can see that the computed acceleration was quite close to the input acceleration.

Result and analyses
Permeability coefficient of the earth core was so small that pore water pressure cannot dissipate during the earthquake.Solid-fluid coupled analysis can give the accumulation process of excess pore water pressure.Figures 11 and 12 present the recorded excess pore pressure in the core center at different elevations (see monitoring points in Fig. 4) and different times, respectively.
As shown in Fig. 11, the excess pore water pressure changed with the earthquake action and accumulated gradually.The maximum excess pore pressure occurred at the bottom of core, as excess pore pressure in the middle of the core bottom was the most difficult to dissipate.
The distribution of excess pore pressure varied with time, as presented in Fig. 12.As we can see, the distribution of excess pore pressure varied with earthquake action and the major excess pore pressure concentrates in the bottom of the core.
The acceleration amplification along the dam axis is presented in Fig. 13, where the acceleration amplification was defined as the ratio of the maximum recorded acceleration to the maximum input acceleration.As the figures shown, the amplification factors increased with elevation with a maximum value of 1.87 in the horizontal direction and 2.91 for vertical direction.And it can be observed that the amplification factor increased sharply near the crest of dam due to the whiplash effect (Gazetas, 1987).
Time histories of acceleration at different elevations were monitored (see monitoring points in Fig. 4), as shown in Figs. 14 and 15.The maximum acceleration increased with elevations, and showed a delayed effect.The maximum accelerations occurred at 2.1 s, 5.6 s, and 6.8 s, respectively.
Figure 16 shows the Fourier amplitude spectrum of accelerations at different elevations.It is shown that the amplifier of acceleration components mainly concentrated in the lower frequency ranging from 0.4 Hz to 2.0 Hz and 0.4 Hz to 3.4 Hz for horizontal and vertical acceleration, respectively.This indicates that the components of input The permanent horizontal and vertical displacements induced by seismic move are presented in Figs. 17 and 18.As seen from the figures, no matter horizontal or vertical permanent displacements, they reached the maximum value at the crest with 0.83 m and 0.84 m, respectively.As a result, special attention should be paid to the design and construction of the dam crest, especially for high earth dams such as Nuozhadu high earth-rockfill dam.
In summary, solid-fluid coupled dynamic analysis can give a reasonable simulation result of earth dams.Strong earthquake can bring seriously damage to the dam crest, which may lead to overtopping resulting in huge floods.Excess pore pressure generated in the core can cause hydraulic fracture, which should be considered in the design of earth dams.

Conclusions
This paper includes analyses of the static and dynamic behaviour of an ultra-high earthrockfill dam in the Southwest China.Mohr-Coulomb nonlinear elasto-plastic model was adopted to describe the behaviour of dam materials.
In the static analysis, staged construction and reservoir impounding were simulated.Then based on the initial static stress field, dynamic response of the dam was studied with fully-coupled nonlinear method.The dynamic behaviour of the dam during the earthquake was analysed, and the seismic characteristics of response were investigated.
The excess pore water pressure changed with the earthquake action and accumulated gradually.The maximum excess pore pressure occurred at the bottom of core.Acceleration amplification factors increased sharply in the upper part of the dam.As a result, special attention should be paid to the design and construction of the dam crest, especially for high earth dams such as Nuozhadu high earth-rockfill dam.
evaluated the seismic stability of an earth dam via the decoupled displacement analysis to compute the earthquake-induced displacements.Parish et al. (2009) performed a study of the seismic response of a simplified typical earth dam using the FLAC3D software.Non-associated Mohr-Coulomb criterion was 2321 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2323
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

) 2325 Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2327
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2329
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | earthquake whose dominant frequency was close to the natural frequency of the dam are amplified significantly due to resonance.

Fig. 2 .
Fig. 2. Maximum cross section with material zoning and construction stage. 17

Table 1 .
Material properties of earth dam soils.