The island of Taiwan is located between the boundaries of the Eurasia and the Philippines Plates and thus earthquakes occurred frequently. The excitation of earthquake affects the integrity of earth dams situated in the mountainous area of the island. A study was conducted to evaluate the dynamic response and safety of one of the earthquake dams. The computer program and soil model used were calibrated for their appropriate use for the subject dam against a well-instrumented centrifuge model. Numerical simulation was then conducted to examine the influence of upstream water storage level on the response of the earth dam. The numerical results identified three locations in the dam where attentions are required because these locations were found susceptible to liquefaction.

1. Introduction

Earthquakes due to active tectonic movements are frequent in Taiwan, which is located between the boundary of Eurasian and Philippines Sea Plates. The plate boundary tectonics are generally dominated by the subduction of the Philippines Sea Plate beneath Eurasia along the Ryukyu Trench that runs from southwest Japan to Taiwan [1].

A magnitude 7.3 earthquake on the Richter scale struck a small town, Chi-Chi (pronounced as Ji-Ji), in central Taiwan on September 21, 1999, as a result of the violent movement of the Chelongpu fault. On October 22 of the same year, another magnitude 6.4 earthquake on the Richter scale struck the town of Jiayi in southern Taiwan, and one of the accelerometers located on the right shoulder of the Renyitan Reservoir dam, which was located four kilometers away from the Jiayi town, detected a maximum ground acceleration of 992 gal. Cracks were observed at the crest of the dam. On March 4, 2010, a magnitude 6.4 earthquake on the Richter scale, with epicenter near the town of Jiaxian, jolted southern Taiwan especially the Gaoxioang county. According to the report summarized by Water-Watch Nongovernmental Organization [2], the earthquake resulted in a 3 cm wide, 75 cm deep, and 15 m long crack across the crest of the dam at the Hutoupi Reservoir, which was about 25 km away from Jiaxian. These are just some of the examples showing how frequent seismic activities occurred in the region and the consequences of the seismic activities on the integrity of earth dam. Earth dams may crumble under seismic loads as the result of soil liquefaction. The significance of conducting thorough dam investigation has been emphasized by Zhang [3].

Finite element (FE) analysis has been widely used for the evaluation of the seismic response and safety of earth dam. Clough and Chopra [4] first introduced the application of FE analysis in the dynamic analysis of earth dam. Seed et al. [5] and Seed [6] back analyzed the response of the 1971 San Fernando earthquake on the deformation of the Lower and Upper San Fernando dams, also using FE analysis. Zienkiewicz et al. [7] included material nonlinearity and liquefaction analysis to study the behavior of the dam materials. Zienkiewicz and Mroz [8] and Pastor et al. [9] proposed the use of the generalized plasticity theory in the response analysis of earth dams under seismic activities. Khoei et al. [10] compared the performance of Pastor-Zienkiewicz and cap plasticity models through the dynamic analysis of the failure of lower San Fernando dam under the 1971 earthquake and the Mahabad and Doroodzan dams under the 1978 Tabas earthquake. Wu et al. [11] employed the computer program LIQCA to simulate the dynamic responses of an earth dam. Gui and Chiu [12] simulated the dynamic response of Renyitan earth dam using the Finn model [13] provided in the two-dimensional finite difference program FLAC. All of these studies yielded convincing results.

When an earth dam is subjected to a seismic load pore-water pressure responded faster than the deformation at the top of dam. Thus, for safety management of earth dam, pore-water pressure response is an important indicator. This study aims at investigating the influence of the upstream water storage level on the pore-water pressure response in the body of an earth dam subjected to synthetic earthquake excitation. Experimental results from a dynamic centrifuge test were first used to calibrate the appropriate use of the chosen constitutive law for the following numerical simulation. Subsequently, effective stress analysis was conducted to evaluate the safety of the study earth dam under the synthetic seismic loading.

2. Governing Equations

Total stress analysis yields only the deformation of soil from the given stress-strain relationship but provides no information on the changes in excess pore-water pressure (EPWP) when the soil is subjected to excitation. Thus, the use of the total stress analysis for liquefaction study is insufficient. Biot [15] incorporated Darcy’s law for the representation of the fluids flow in the soil to the linear elastic solid mechanics framework and thus formulated the so-called two-phase mixing theory to deal with the pore-water pressure issues that the total stress analysis failed to handle. The equations of motion for elastic porous media saturated with a pore fluid consist of two variables’ fields: (i) pore pressure and (ii) skeleton displacements . A simplified numerical framework, known as the formulation, was formulated by assuming Darcy’s flow rule (laminar flow) and neglecting the fluid inertia [7, 1618]. By denoting was the flow rate of the fluid, we thus have the overall balance or momentum equilibrium equation of soil-liquid mixture as [19]where and are the total mass and total stiffness matrixes, is the discrete gradient operator coupling the soil and fluid phase, is the permeability matrix, and is the compressibility matrix. Vectors and represent the effect of body forces along with prescribed boundary conditions for the solid-fluid mixture and fluid phases.

The model used to simulate the dynamics behavior of the soil is the Pastor-Zienkiewicz mark III (P-Z III) model, which has been recognized as the model most suitable for modeling materials susceptible to liquefaction [19]. The model was proposed in 1985 by Pastor et al. for saturated sands subjected to dynamic excitation and liquefaction was imminent. The main advantage of the theory is that neither the yield surface nor the plastic potential surface needs to be explicitly defined; see, for example, Merodo et al. [20] and Khoei et al. [10]. Detailed description of the theory and constitutive relation for P-Z III can be found in, for example, Merodo et al. [20]; Khoei et al. [10]; Borowiec [19], Zienkiewicz and Mroz [8]; they are thus not repeated here.

3. Calibration of Soil Model and FE Program

To calibrate the appropriateness of the soil model chosen for this study, it is necessary to calibrate the result obtained from the finite element analysis with that of a well-instrumented experimental test.

3.1. Centrifuge Model

A dynamic centrifuge test, performed at the University of Colorado at Boulder, was selected for this purpose [11, 14, 21]. The centrifuge model test selected for this study was one of the three tests performed in the University of Colorado to obtain data for safety assessment of a dam in Taiwan that was rocked by a strong earthquake on 21 September, 1999. The soil used in the construction of the model dam was taken within the reservoir of the dam so that a better representation of the behavior of the prototype dam could be simulated in the centrifuge. The centrifuge has a maximum capacity of 400 g-ton and equipped with an electrohydraulic shaking table. Thus, the shock wave was generated by this in-flight electrohydraulic shaker table, which is a displacement-controlled system [21].

For all the three centrifuge models tested at the University of Colorado, the slope at the upstream and downstream sides was 1 : 3 and 1 : 3.5, respectively. Model I of the earth dam was built entirely using the prototype core material, which is the low-plasticity clay (CL), and it was founded directly on the impervious and rigid base of the aluminum model container [21]. Model II was analogous to Model I but with a 33 m thick low-plasticity silt (ML) foundation. Model III was built to mimic the actual section of the prototype dam, as shown in Figure 1. The materials used to construct the model dam included low-plasticity clay (CL) for the core, low-plasticity silt (ML) for the upstream and downstream slopes, and silty-sand (SM) for the 33 m thick foundation. All the models were scaled down 150 times from the actual dam; the models were then subjected to 150 g of gravitational acceleration to replicate the field stress state as in the actual dam [21]. Model III is the subject of this study.

To provide a better representation and understanding of the responses of the dam during excitation, a number of accelerometers, pore-water pressure transducers, and linear variable differential transformers were installed in the body of the model dam as shown in Figure 1. Specifically, eight accelerometers and six pore-water pressure transducers were installed to measure the time-history responses of horizontal and vertical accelerations and the EPWP, respectively, during the seismic excitation of the dam. In addition, four linear variable differential transformers (LVDTs) were installed at the dam crest and the shoulder of the downstream slope to monitor both the horizontal and vertical displacements of these locations during excitation.

3.2. Soil Parameters

The finite element program used for the following dynamic numerical analysis is called DIANA-SWANDYNE II. The Pastor-Zienkiewicz (P-Z) model was already built-in as its constitutive model. The parameters of the P-Z constitutive model were determined through a series of one-dimensional consolidation test, drained and undrained monotonic triaxial consolidated test (see, e.g., [22]), and undrained cyclic triaxial test and triaxial permeability test. The parameters of the constitutive model for the three types of the soil used in the simulation of the model dam are summarized in Table 1, in which is the slope of the critical state line, , , , , , , and are constitutive parameters of P-Z model [9, 16, 23], while is the mean effective confining stress.

3.3. Experimental and Numerical Results

The time-history accelerations and time-history pore-water pressure obtained from both the centrifuge test and the numerical analysis are compared and discussed in this section. The time-history of accelerometer ACC1, installed at the base of the centrifuge model, was used as the input of the earthquake acceleration in the numerical analysis.

3.3.1. Acceleration Comparison

The full 35 seconds of time-history horizontal accelerations recorded in the centrifuge and the numerical analysis are compared in Figure 2, while their corresponding maximum amplitudes are summarized in Table 2. The input shock wave was essentially the 1999 Chi-Chi earthquake ground motion but with all the frequencies higher than 400 Hz filtered out because of the capacity of the shaking table used in the University of Colorado [21]. The maximum amplitude recorded at the base of the centrifuge model (ACC1) was 0.099 g and as the shock wave was transmitted to the crest of dam the maximum amplitude at the crest (ACC2) became 0.125 g. Thus, the maximum amplitude of the acceleration was amplified by about 125% as the shock wave was transmitted from the base to the top of the dam. The numerical result overestimated the centrifuge result by about 25%, with an amplification of about 150%. This sort of topographic amplification is to be anticipated especially at the crest of a sloping ground; see, for example, Bouckovalas and Papadimitriou, 2005 [24] and Brennan and Madabhushi, 2009 [25]. The deviation between the two ACC4 (vertical accelerometer) profiles could be that it was placed near the dam crest of the dam and as such it might not be firmly hold in position, which resulted in its failure to capture the real excitation as compared to that of the numerical result.

3.3.2. Pore-Water Pressure Comparison

Pore-water pressure transducer PPT6 was installed at the upstream toe of the model dam to monitor whether the upstream water storage level was kept at the predetermined level during the in-flight test. Five more pore-water pressure transducers were installed in the body and the foundation of the dam, as shown in Figure 1, for monitoring the pore-water pressure response of the dam during in-flight earthquake excitation.

Figure 3 shows the time-history of EPWP recorded during the centrifuge test. The excitation of the EPWP was relatively small in the upstream silt shell (PPT1), clay core (PPT3), and downstream silt shell (PPT4) of the dam. The centrifuge data shows that most of the excitation was registered in the SM foundation layer of the model dam. However, as PPT5 and PPT2 were installed not very far apart, the discrepancy in the recorded EPWP of both transducers indicated that PPT5 could be malfunction at time of testing. The time-history of EPWP produced by the numerical analysis is shown in Figure 3. The numerical result shows that some 10~20 kPa of EPWP were observed in the SM foundation, CL core, and downstream ML shell. Both the centrifuge and numerical results confirmed that at the location of PPT1 in the upstream ML shell the EPWP was able to dissipate quickly such that no EPWP was recorded in this location. The centrifuge model registered some 20 kPa of EPWP at the end of the excitation while the numerical analysis underestimated the value by about 25%.

The contour of the initial hydrostatic PWP and the EPWP at = 5, 10, 15, 25, and 35 seconds of excitation are plotted in Figures 4(a)4(f). Figures 4(b)4(f) revealed that negative EPWP was generated between = 5 and 35 seconds near the intersection of the clay core, downstream ML shell, and the SM foundation of the dam and extended below the clay core, that is, along the interface between the clay core and the SM foundation layer, whereas positive EPWP was generated along the interface between the downstream ML shell and the SM foundation layer. The region of the positive EPWP was expanding with elapsed time and moving towards the downstream toe, as indicated in Figure 4(f). As shown in Figure 3, PPT4 revealed that the EPWP was slightly on the negative side throughout the earthquake excitation in the centrifuge but the numerical result shows otherwise with an EPWP of about 18 kPa towards the end of the simulation. One could argue that PPT4 could well be shifted during the centrifuge test from its original location to a location slightly above the water table in the dam.

4. Dynamic Finite Element Analysis

4.1. Model Setup

The 1,535 m long Renyitan dam is an off-stream dam with a height of 20.3 m and a 9 m wide crest; the dam has a water-storage area of about 3.66 km2 with a total water-storage capacity of about 29.11 million cubic meters, of which the effective water-storage capacity is 27.31 million cubic meters [12]. The construction of the dam at the upstream of the Bazhang River took seven years and was completed in 1987. A typical cross-section of the dam is shown in Figure 6. The core of the dam was built with low plasticity clay (CL), the outer shell was built with silty-sand (SM), the transition shells, sandwiched between the outer shell and the core, were built with low plasticity silt (ML), and the filter layer at the downstream side of the core was built with silty-sand (SM). The soil stratum of the foundation of the dam may be broadly divided into silty-sand (SM), low-plasticity silt (ML), and silty-sand (SM) layers.

The cross-section selected for the dynamic numerical analysis is shown in Figure 5. Figure 6 shows the layout of the pore-water pressure monitoring points in the study dam. In total, 28 pore-water pressure points were assigned to monitor the changes of the pore-water pressure at these locations when the dam was numerically subjected to earthquake excitation. Figure 6 also shows various water levels that would be used to study the influence of the upstream water storage level on the stability of the dam.

4.2. Synthetic Earthquake Simulation

The acceleration time history to be used in numerical analysis may be selected from the (i) past acceleration time history captured in real earthquakes; (ii) simulated synthetic ground acceleration from theoretical seismological models of seismic fault rupture; and (iii) simulated spectrum matched artificial acceleration using stochastic or random vibration theory [26]. This study adopted the simulated spectrum-matched artificial acceleration because it is generally recognized that it is able to produce results that present relatively lower dispersion, which is important for nonlinear analysis [27]. In this case, the design ground acceleration was obtained through the modification of the ground acceleration captured from the magnitude 7.3 earthquake on the Richter scale that rocked the island of Taiwan on September 21, 1999.

A trapezoidal spectrum was selected for deriving the synthetic earthquake. The trapezoid envelope rises to its highest point and then stays at the plateau before falling back down to zero. The trapezoid envelope thus represents the components of the rise time , the strong ground movement duration , and the total ground movement duration . The rise time signifies the time difference between the arrivals of P-wave and S-wave; the time for Taiwan is approximately 3 seconds according to Yeh and Tsai [28]. The duration of the strong ground movement is an important parameter for the energy spectrum density function; from the major earthquake records in Taiwan, the average in Taiwan is about 5.07 seconds with a standard deviation of 3.60 seconds [29]. The longest was assumed to be the sum of both, that is, 8.67 seconds. Based on the total duration relationship for various earthquake scales and ground vibrations greater than 0.05 g [30], the total ground movement duration corresponding to the scale 7.2 is 35 seconds. In addition, the maximum ground acceleration of 1 g and a magnitude of 7.2, determined from the past two largest earthquakes in 1941 and 1999 recorded nearby the study area, were also required information.

With the above information as input parameters and using the principle of natural random excitation, the synthetic earthquake program, SIMQKE-II [31], was then used to generate the synthetic earthquake that met the criteria of the design response spectrum of the Renyitan Reservoir. The time-history acceleration of the synthetic earthquake generated by SIMQKE-II is shown in Figure 7.

5. Results and Discussion

5.1. EPWP

Twelve (PPT1–PPT12) of the 28 pore-water pressure monitoring points were placed in the body of the earth dam and the rest (PPT13–PPT28) were installed in the foundation layers. PPT1, PPT6, and PPT10 were located in the outer shell of the dam at the upstream side, which is made of silty-sand (SM). PPT2, PPT7, and PPT11 were in the transition shell at the upstream side, which is made of low plasticity silt (ML). PPT3, PPT8, and PPT12 were in the low plasticity clay (CL) core of dam. PPT4 and PPT5 were in the transition shell and outer shell, respectively, of the downstream side. PPT9 was placed in the filter layer at the downstream side. PPT13 through PPT28 was in the foundation layers.

A filter layer was present on the right side of the core at the downstream side because of the need to reduce the phreatic surface running across the dam. The seepage and stress analyses provided the initial hydrostatic pore-water pressures and the initial vertical stresses of the dam for the dynamic analysis. The subsequent dynamic analysis provided the EPWP generated in the dam as a result of the synthetic earthquake excitation. Table 3 shows the initial stresses and EPWP measured at the pore-water pressure monitoring points when the upstream water storage level was at 17.3 m.

The dynamic numerical results are evaluated in terms of the normalized EPWP, defined as the ratio of the EPWP to the initial vertical effective stress. It is shown in Table 3 that the EPWP ratio was close to zero at PPT4, PPT5, and PPT9 at the downstream side of the dam. The EPWP excited was negative at PPT1, PPT7, and PPT11 at the upstream side, resulting in the negative EPWP ratios. Most of these ratios ranged from 10 to 80%. The EPWP excited at PPT3 located in the core reached as high as 215 kPa with a normalized EPWP ratio of nearly 90%. However, the EPWP ratio at PPT18 was greater than 100%, indicating that the EPWP was greater than the initial vertical effective stress in the shallow part of foundation at the downstream side, thus leading to instability in soil mass.

5.2. Influence of Upstream Water Storage Level

To examine the effect of upstream water storage level on the stability of the dam at the end of the synthetic earthquake excitation, five hypothetical upstream water storage levels at 7.6 m, 10.2 m, 12.7 m, 15.2 m, and 17.3 m were assumed to take place at the 20.3 m high Renyitan earth dam. The distribution of the EPWPs generated in the dam for water storage level of 0.37, 0.50, and 0.75 of the dam height is shown in Figure 8 while Table 4 shows the critical (maximum) EPWP generated and its coordinates in the dam.

Figure 8(a) shows the distribution of the EPWP when the water storage level was at 37% of the dam height. The critical EPWP of 274 kPa occurred near the bottom part of the upstream transition shell, albeit the affected area was very very small. In general, the whole dam suffered an EPWP of not more than 25 kPa. In the middle layer of the downstream foundation soil, mostly directly beneath the downstream transition shell, it was shown that a small zone was bounded by the 200 kPa EPWP.

Figure 8(b) shows the distribution of the EPWP when the water storage level was at 50% of the dam height. This time, the critical EPWP of 290 kPa occurred near the bottom of the core; in addition, the bottom part of the upstream transition shell also suffered from large EPWP, that is, greater than 200 kPa. Except for a small portion of the top left corner of the dam that has suction (negative pore-water pressure) generated during the excitation, in most parts of the dam their EPWP was less than 25 kPa. As for the foundation layer, the area bounded by the 200 kPa EPWP has expanded to a larger area compared to the case shown in Figure 8(a).

The distribution of the EPWP when the water storage level was at 75% of the dam height is shown in Figure 8(c). The critical EPWP also occurred near the bottom of the core of the dam and it was also shown that an even larger portion of the top of the dam, compared to Figure 8(b), has suction generated at the end of excitation. As for the foundation layer, the area bounded by the 200 kPa EPWP has expanded to an even larger area compared to the cases shown in Figures 8(a) and 8(b).

In general, it was seen that the area bounded by the 200 kPa EPWP in the foundation layer expanded with the increased of the upstream water storage level. The critical EPWP tended to occur near the bottom of the upstream transition shell and the core of the dam; this is because the filter layer placed on the right shoulder of the core and the bottom of the downstream transition shell was able to dissipate almost all the EPWPs generated at the end of the excitation.

The generation of excess pore-water pressure (EPWP) varied across the earth dam as the water storage level increased. The final EPWPs, normalized by their corresponding initial effective vertical stress, across the dam at the end of the earthquake excitation are presented in Figure 9. Figure 9(a) shows the normalized EPWP against the normalized water storage level for PPT1, PPT6, and PPT10 located in the upstream outer shell. PPT10, which was located near to the foundation layer, showed a linear increase between the EPWP and the increment of the water storage level. For PPT6, the normalized EPWP showed no significant variation with the increase of the water storage level. The normalized EPWP for PPT1 showed a negative relationship with the increment of the water storage level; thus, the EPWP near the top-third of the upstream outer shell was becoming more and more negative (suction) until the water level reached 75% of dam height where a reversal of trend is seen. In general, the maximum EPWP ratio in the upstream outer shell was 0.6 and, hence, no liquefaction is possible in this region.

Figure 9(b) shows the ratios of the EPWP and water storage level, recorded by PPT2, PPT7, and PPT11, in the upstream transition shell. PPT2, which is located near the midheight of the dam, showed a similar trend as PPT1 in Figure 9(a) albeit this time the decrement was not as significant as in PPT1. Thus, suction existed in the dam if the water storage level was less than 75% dam height.

The profiles of the EPWP ratio registered by PPT3, PPT8, and PPT12 in the dam core are shown in Figure 9(c). PPT3, which was located near the midheight of the CL core, showed a positive relationship with the increment of the water storage ratio. The EPWP reached some 90% of the initial effective stress when the water level ratio was at 0.85. Thus, the stability of the core requires attention at this instance. The EPWP ratio of PPT8 only increased slightly from 50% to 55% for water storage ratio that ranged between 0.37 and 0.85. The EPWP ratio of PPT12, which was located near the interface between the core and the foundation layer, remained reasonably constant at 0.2 until the water storage level reached 75% of the dam height where the EPWP ratio began to drop.

Figure 9(d) shows that no EPWP was generated in the transition shell, outer shell, and the filter layer at the downstream side of the dam regardless of the level of the upstream storage water. This was because these PPTs were placed in the vicinity of the highly permeable filtering layer.

The EPWP ratio profile registered by the PPTs in the upper layer of SM foundation is shown in Figure 9(e). As shown in Figure 6, PPT13, PPT14, PPT15, and PPT20 were at the upstream side while PPT16, PPT17, PPT18, and PPT21 were at the downstream side of the dam. PPT19 was placed directly below the core. It was found that the EPWP decreased slightly with the increase of the water storage level in the upstream upper layer of the foundation while the EPWP did not change with the water storage level in the downstream foundation layer. However, the EPWP ratio for PPT18 was greater than 1.0 for all values of upstream water storage level; this was due to the small overburden pressure at this point and hence suggested that the ground here would be unstable.

The EPWP generated in the middle layer of the ML foundation did not seem to approach the critical value, at which liquefaction would occur (Figure 9(f)). PPT22 and PPT23 at the upstream side have negative relationship with the increment of the water storage level, while PPT24 and PPT25 at the downstream side have positive relationship with the increment of the water storage level. The ratio for PPT25, located to the furthest right, increased from 60% to 80% of the initial effective stress for the range of the water storage level studied.

The relationship between the EPWP, normalized by the initial vertical stress, of all the 28 pore-water pressure monitoring points and the initial hydrostatic pore-water pressure is plotted in Figure 10. Some of the normalized EPWP ratios were greater than 1.0 at locations near to the ground surface with hydrostatic pore-water pressure of less than 30 kPa. For places with 100 to 250 kPa of hydrostatic pore-water pressure, in particular, those below the dam core and in the downstream foundation layer, their EPWP ratio was higher than 80% of the initial vertical stress. The soil mass here was susceptible to liquefaction and, therefore, their stability is a concern.

6. Conclusion

A series of dynamic finite element analysis under a synthetic earthquake was conducted using a two-dimensional program in this study. The computer program, DIANA-SWANDYNE II, and the soil model adopted in this study were first calibrated by simulating the response of a model earth dam under an earthquake excitation in the centrifuge. The following conclusions can be made.

The amplification of the earthquake wave as it was transmitting to the top of dam and the response of the EPWPs were captured and compared between the centrifuge test result and the numerical result. Good agreement has been achieved between both sets of result, indicating the reliability of the chosen program and soil model.

The silty-sand filter layer placed on the right shoulder of the core and the bottom of the downstream transition shell was effective in dissipating almost all the EPWPs generated at the end of the excitation.

The result of the analysis revealed that the EPWP near the top-third of the upstream outer shell was becoming more and more negative (suction) until the water level reached 75% of dam height indicating that this region will be in a stable state as no liquefaction is possible during excitation.

The results from dynamic analysis with five different water storage levels suggested that the zone just above the interfaces between the foundation soil and the upstream ML transition shell and the CL core was susceptible to liquefaction if the level of the upstream water storage was close to the dam height. In addition, it was observed that EPWP built up quickly in the middle ML layer of the downstream foundation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The authors are grateful to the financial support provided by Taiwan Water Corporation for conducting the above centrifuge tests.