© 2008 Science Publications Geo-Seismic Environmental Aspects Affecting Tailings Dams Failures

Seismic performance evaluations of tailings dams are essential for characterizing the geo- environmental risks posed by these earthen structures, which should include the geotechnical hazards implied by slope instability failure, free board loss and the potential release of contaminants. The observed damage is more important when liquefaction occurs on the dam body and foundation, which often leads to cracking, settlements, tilting and general distortion of dam geometry. Analyses based on limit equilibrium are generally sufficient to establish hazard zones. However, numerical models with solution schemes formulated in the time domain, which are capable of taking into account the kinematics of soil movement more realistically, are needed to quantify the geotechnical risk. This paper reviews the main geotechnical earthquake engineering aspects to account for when designing tailings dams and describes the application of a practice-oriented simplified constitutive model, which implemented in a lagragian finite difference platform, is capable of predicting the accumulation of pore pressure in fine-grained saturated materials due to earthquake loading, the reduction of shear strength and the corresponding permanent displacements. The model uses the Mohr-Coulomb failure criterion coupled with an incremental pore pressure generation scheme. Pore pressure is accumulated as a function of the number of stress cycles. The secant soil stiffness and hysteretic damping change with loading history. The numerical simulation is able to properly capture the kinematics of dam failure and provides the parameters to assess potential environmental impacts on the nearby areas of the dam.


INTRODUCTION
Mining activity and other industrial processes generate large quantities of solid waste that annually need to be disposed of and stored in the so-called tailings dams in different parts of the world. The search for an economically, technically sound and environmentally admissible alternative to deal with this type of materials is a multidisciplinary activity where the geotechnical engineer plays a key role. Among the different aspects to consider in tailings dam design, seismic dam stability analysis is the most relevant in projects located in earthquake prone areas to ensure that catastrophic failures will not occur. In practice, these evaluations very often become challenging, due to presence of very fine nonplastic uniform granular soils, both in the tailings dam body and foundation. In addition, these materials can be totally or partially saturated and exhibit a very loose structure if the hydraulic fill method is used to transport and dispose of them. Hydraulic filling is the preferred method by miners for construction of tailings dams due to its low cost, especially considering that tailings dams do not have a direct economic value. Thus, detailed geotechnical investigations are not carried out in most cases and tailings dams are built based on simplified empirical recommendations, usually given by mine owners that increase dam height as a function of the storage volume required.
In order to reduce the environmental hazard associated with contaminant release, dam seismic stability analyses are necessary to ensure that slope instability, freeboard loss, or a larger catastrophic failure, such as those reported in technical literature [1,-4] will not occur. Seismic analysis of tailings dam should account for the degree of internal drainage and the characteristics of the design earthquake, to properly estimate the amount of pore pressure built-up during the seismic event, as well as the pertinent impact on the reduction of shear strength and the modification of dynamic soil properties (i.e., shear stiffness and damping). Based on these evaluations, proper assessment of the volume of material mobilized during a potential failure and the pattern and speed of deformation can be achieved, thereby providing guidance to estimate the extension of the affected zones. Thus, maps of geotechnical hazard and quantification of specific risks can be developed.
Pseudo-static methods, based on limit equilibrium and others formulated considering the sliding block mechanism are sufficient to define geotechnical hazard zones but insufficient to quantify the seismic risk, which requires a complete understanding of the kinematics of the slope failure process [5] . A complete representation of the failure mechanism usually requires using constitutive models to simulate the dynamic response of the soil, coupled with numerical techniques, such as finite differences or finite element methods. Such models have evolved over time, from the classical hyperbolic Masing-type models [6] to more sophisticated models based on bounding surface hypoplasticity theory [7] .
A fully nonlinear dynamic analysis of tailings dams requires using a numerical model able to account for the effect of generation and dissipation of pore pressure within the dam body and foundation and its impact on the variation of shear strength which, in turn, will lead to permanent displacements. Perhaps, the best representation of this phenomenon can be achieved with a fully coupled constitutive model, where the equation of motion is solved simultaneously with the diffusion equation during an effective stress analysis [7] . However, in general, the more sophisticated a constitutive model is, the more cumbersome it becomes to use in engineering practice.
This study reviews some of the key geo-seismic environmental aspects to consider during seismic evaluations of tailings dams and presents a practiceoriented numerical scheme which, implemented on a lagragian platform, appears to capture the overall behavior of ground response in geotechnical problems well in which large deformations are likely to occur, such as liquefaction or cyclic mobility. This approach enhances the formulation proposed by Dawson et al. [8] . The method allows for incorporating the degree of internal distortion of the mobilized soil mass, drainage conditions, localization of failure planes near the slope surface and the plastic yield at depth. The methodology is illustrated through its application to the analysis of a case history and the model predictions are compared with field measurements, both qualitatively and quantitatively. This approach, along with a risk analysis framework [9] , can be used in dam safety evaluation assessments.

ANALYTICAL FRAMEWORK
Seismic parameters: The parameters required to characterize the seismic loads to use in the evaluation of the seismic response of an earth dam, such as a tailings dam, depend on the method of analysis considered. Nevertheless, in general, the following parameters should be accounted for: the earthquake magnitude, the acceleration response spectrum for design and the maximum displacement expected in the dam foundation. In particular for tailings dams, acceleration time histories compatible with the aforementioned parameters are required as well. The geotechnical subsoil conditions and its spatial variability (e.g., from rock to deep soil), at a particular site, which may significantly affect ground motions, should be fully identified. Based on this information, it is easier to establish possible dam failure scenarios to be considered in risk analyses such as the rupture of the dam due to movement in its foundation, free board loss due to differential tectonic movements, slope stability failures associated to both foundation movement and ground shaking, slide of a section of the dam, cracking and tilting, among other things.

Liquefaction:
The problem of liquefaction during strong ground shaking, which can be understood as a shear strength drop exhibited by loose granular poorlygraded saturated soils, is due to the development of excess pore pressure during cyclic loading and has been traditionally studied by following two approaches: cyclic stresses and cyclic strains. Historically, methodologies based on stress cycles have been the most popular in practice, mainly because they can be estimated directly from an acceleration time history measured or computed in a given soil profile [11] , whereas the estimation of strain cycles requires making assumptions regarding the stress-strain relationship of the soil. The cyclic stress approach considers that the liquefaction potential of a soil stratum is a function of the number and magnitude of shear stresses applied to the soil during the dynamic event. Cyclic shear stresses are expressed in terms of the Cyclic Stress Ratio (CSR), defined as the ratio between the cyclic shear stress τ cy acting on the failure plane, divided by the initial vertical effective stress σ´v o .
This relationship depends mainly on the initial state of stresses existing in the ground and seismic load, expressed in terms of the number of equivalent cycles, which in turn is defined by the duration, intensity and frequency content of the earthquake.
Soil resistance to liquefaction: It is common practice to describe soil resistance to liquefaction through a curve of cyclic strength, which relates the number of cycles required to have a pore pressure ratio of one,   [10] ) with its corresponding cyclic stress ratio (Fig. 1). The pore pressure ratio, r u , is defined as r u = ∆ ug /σ vo´, where ∆ ug   is the generated pore pressure. These curves can be obtained from laboratory testing, or derived from CSR plots developed from field measurements and back analysis of case histories. Thus, the cyclic strength of a soil is mainly a function of its relative density, including grain gradation and shape and fines content.
Cyclic stress approach: An analysis based on the cyclic stress approach should at least include the following steps: C Evaluation of shear stress time histories at each soil layer conducting numerical analyses of wave propagation or soil-structure interaction. C Approximation of the shear stress histories to an equivalent number of uniform stress cycles, N eq C Comparison of the number of equivalent cycles with the liquefaction resistance curve of the soil to determine if there will be liquefaction and how close the soil is to the condition r u = 1. C Evaluation of the effect that the computed number of equivalent stress cycles has in the soil, in terms of change of volume and increase of pore pressure.

NUMERICAL SCHEME
The numerical scheme utilized in this study is based on the so-called simplified procedure, proposed    Shear strength is assumed to follow a modified bilinear Mohr-Coulomb law (Fig. 2). In this manner, soil strength is a function of the relative density of the material, fines content and residual strength. Even though this type of formulations considerably simplifies many aspects of the real soil behavior because it considers that pore pressure generated is due directly to cyclic shear stress instead of a contraction of the soil skeleton, this type of models seems to provide adequate estimations to geotechnical problems where plastic deformations are significant. The validity of similar approaches to predict permanent deformations due to earthquake loading has been demonstrated by Roth et al. [12] through comparisons with centrifuge test data for dry sands and by Dawson et al. [8] , for saturated sands, by comparing numerical predictions with actual measurements obtained directly from case histories.
Mechanism of pore pressure generation: The mechanism of pore pressure generation is an enhancement of an existing formulation proposed by Dawson et al. [8] and is illustrated in Fig. 3. The history of horizontal shear stresses is monitored and each time the shear stress passes through zero twice, a semi-cycle of that particular cyclic stress ratio, CSR, obtained as τ cy /σ´v o is accumulated (Fig. 3), developing and amount of excess pore pressure, ∆ ug . This ∆ ug is described in terms of the pore pressure ratio. Each semi-cycle will generate an increment in pore pressure proportionate to the total number of cycles which are necessary to reach the condition r u = 1 for a given CSR. The increment of pore pressure, ∆u gi , is computed as ∆u gi = ∆r ui σ´v. When pore pressure goes up, effective stress goes down and, in turn, shear strength and soil stiffness decrease. Tangent soil stiffness is modified using the following expression G = G max (σ´v/σ´v o ), where σ´v is the current value of vertical effective stress for that semi-cycle. Thus, the secant modulus will decrease not only due to the loss of shear strength, which decreases as a function of σ´v, but also directly due to changes in the tangent modulus of plastic stress-strain hysteretic loops.

Fig. 5: Model response for irregular cycles
Residual strength after liquefaction: After liquefaction, the shear strength of soil is mostly provided by its residual strength, s r . Thus, s r , is a key parameter to assess levels of permanent displacement levels and the overall dam stability at the end of the earthquake. The model presented herein incorporates the residual strength using a bilinear failure envelope. This envelope consists of an initial cohesion and friction value of zero, which extends until it intersects the Mohr-Coulomb envelope (Fig. 2). Residual strength can be estimated from Standard Penetration test, SPT, corrected blow counts for clean sand (N 1 ) 60-cs [13] . Thus, residual strength is the lower bound of shear strength and will essentially define soil resistance when the pore pressure ratio, r u , becomes 1.
Model Performance: To illustrate the model performance, a parametric study was conducted using three different uniform stress histories, with amplitudes and frequencies varying from 9.81 kPa to 49.03 kPa and from 0.10 to 0.02 Hz, respectively, as shown in Fig. 4. The material was assumed as saturated clean sand with a corrected SPT (N 1 ) 60-cs value of 17, a shear wave velocity of 250 m secG 1 and a residual strength, s r , of approximately 33 kPa.
The initial vertical stress was 101.3 kPa (i.e., 1 atm) and the initial pore pressure was assumed to be 20 kPa. As can be observed, pore pressure increases linearly each semi-cycle, starting from initial pore water conditions. For tests 2 and 3, the pore increases until it reaches total vertical stress after approximately 12 sec and thus, r u becomes one. For test 1 no appreciable increase in pore pressure is observed, due to the relative low shear stress amplitude.
Similarly, Fig. 5 shows the model response for a nonuniform stress history. The maximum stress amplitude was about 49 kPa and the mean frequency was 0.16 Hz. The material properties were the same as those considered in the previous case. Initial vertical stress was 101.325 kPa. Again, the initial pore pressure was assumed to be 20 kPa. The material reaches initial liquefaction after approximately 31 sec.

MODEL APLICATION TO THE EVALUATION OF A CASE HYSTORY
The model described above was implemented in a finite difference platform to be used for analyzing the dynamic response of potentially liquefiable earth dams, such as tailings dams. A case study, which has been analyzed in the past by numerous researchers, namely, the performance of the Lower San Fernando Dam during the Northridge Earthquake, was considered to illustrate the methodology.
Background: Located in Los Angeles County, the Lower San Fernando dam, along with the Los Angeles and Upper San Fernando dams forms the Van Norman hydraulic complex (Fig. 6), which controls from 50% to 75% of the total water supply of the city of Los Angeles [4] . Due to its strategic importance, its seismic behavior observed during several seismic events has been studied by numerous researches and practitioners [2,4,8,[14][15][16] in the past. During the 1994 Northridge earthquake, the Van Norman Complex was severely shaken, leading to generation of sand boils, lateral spreading, settlement and cracking along the complex and particularly in the Upper and Lower San Fernando hydraulic fill dams. Similarly, the Los Angeles reservoir underwent measurable movement.
Lower San Fernando Dam: Although the basin was almost empty at the time of the earthquake, the water level was usually found from 3 to 10 m below the crest of the dam. A typical cross section of the dam is presented in Fig. 7. The upstream slope of the dam was reconstructed after the 1971 San Fernando Earthquake. The dam is approximately 40 m high and was originally built using the hydraulic fill method. However, on several occasions, compacted fill was added to increase its storage capacity [8] . Material properties are summarized in Table 1. The core clay shear modulus was considered for the analysis of 244600 kPa. For the other materials, it was obtained as a function of the state of stresses according to the expression: G =22 K 2 (σ´m) 0.5 , where σ´m is the mean effective stress (in kPa). The values of K 2 are included in Table 1. The resistance to liquefaction was obtained using the applicable blow counts, corrected by energy, obtained during SPT test. The seismic performance evaluation of the dam was conducted considering the 1995 Northridge earthquake. CPT measurements taken after the earthquake indicate that the slide debris has an equivalent clean-sand blow count (N 1 ) 60 of approximately 15. Figure 8 shows the acceleration time history utilized in the analysis, which corresponds to what was recorded in a seismological station located in the LA Dam west abutment record and has a Peak Ground Acceleration (PGA) of 0.43 g, a fundamental period of 0.3 sec and an approximate duration of the severe part of the earthquake of 15 sec. Although a much higher PGA (0.86 g) was recorded at the Rinaldi Receiving station, located near the downstream toe of the lower San Fernando dam, the record apparently shows some site amplification due to its alluvial nature.

Input ground motion:
Analysis results: The finite difference model of the dam utilized in the dynamic analysis is shown in Fig. 9. Figure 10 shows the pore pressure distributions before and after the earthquake and Fig. 11   of maximum vertical displacements. Significant amounts of pore pressure developed after the earthquake can be seen in the majority of the granular material located in the upstream face of the dam. A comparison between the computed permanent horizontal displacement and field measurements are presented in Fig. 12. It can be seen that the model predicts the measured displacements favourably. The failure mechanism of the slopes can be observed in the plot of displacement vectors presented in Fig. 13. It is clear that according to the simulation, when the soil liquefies, the failure does not occur as a rotation, but preferably as a progressive movement in horizontally oriented soil layers, which concurs with the observation of cracks and general ground movement observed along the crest of the dam, during the damage survey conducted after the earthquake [14,15] .   Figure 14 shows the horizontal and vertical displacement computed at the tip of the upstream slope of the dam. It can be seen that the predicted displacement (0.15 m) concurs well with the measured response (0.16 m). Figure 15 shows the evolution of the pore pressure ratio, r u = ∆u g /σ v´, during the earthquake, for points A, B and C identified in Fig. 10. How r u reaches a value of one can be observed at points A, B Fig. 15: Pore pressure in points A, B and C of the dam and C. Pore pressure increases significantly from the second 3.5, which as can be seen in the acceleration time history, corresponds to the beginning of the intense part of the earthquake. Furthermore, a pulse can be observed in t = 3.5 sec approximately, clearly identified in the velocity time histories computed in the aforementioned locations (Fig. 16), which induced quite large shear stresses in the soil in a relatively short period of time, thereby leading to a significant

CONCLUSIONS
Quantifying seismic risk zones around dams built with fine-grained saturated materials, such as tailings dams, requires performing analyses that take into account the degree of internal distortion of the soil mass generated during the event, drainage conditions, the localization of strength planes close to the surface of the soil and the plastic yielding to depth. Simple models, such as the one presented in this paper, together with lagragian finite difference formulations allow for capturing the overall response of the dam during the earthquake qualitatively as well as quantitatively. They are also able to predict the evolution of pore pressure, shear strength during the earthquake and the permanent displacement associated, as shown in the example discussed of the lower San Fernando Dam. The aforementioned analysis methodology, coupled with a risk analysis framework, can be used to assess dam safety and potential environmental risks.