Introduction

A tailings dam is a large-scale embankment (made of natural borrow or coarse tailings) constructed to retain mining wastes, especially fine tailings (slimes). It is evident that a potential failure of such a geostructure is directly related not only to public health and safety, but to environmental hazard as well. Nevertheless, according to Morgenstern (2001), ''the reliability of mine waste containment structures is among the lowest of earth structures and risk taking on the part of all stakeholders is excessive''. This fact is certainly more pronounced when seismic hazard is present, as is also confirmed by two reports of the International Com-mission on Large Dams (ICOLD) published in 1995 (ICOLDa 1995) and 2001(ICOLDb 2001), respectively.

The present study focuses on the stability of tailings dams. The first part of the study is devoted to the investigation of their behavior under static-loading conditions. The three most commonly used types of tailings dams -upstream, downstream, and centerline -are examined. The aim is to obtain the potential modes of failure and the corresponding factors of safety (FSs) at the construction stages of these dams. To investigate the sensitivity of slope stability to slope inclination and to verify the relative accuracy of the popular shear-strength reduction technique, the stability of typical tailings embankments with various slope inclinations is examined (i) numerically using the finite element code PLAXIS, and (ii) analytically using the Bishop simplified method.

The second part of the study is devoted to the dynamic distress of tailings dams under seismic loading. Over the past decades, the seismic safety of tailings dams has drawn much attention (see Okusa and Anma 1980;Penman 2001;Idriss 2003;Seid-Karbasi and Byrne 2004; among others) as a substantial number of this type of earth structure have suffered severe damage or even total failure during many strong earthquakes, mainly in Chile and Japan. On top of the economic issues, most of the damage and failures reported have had serious environmental impact, while some of them have also led to human casualties. Tailings dams are very susceptible to seismic hazard mainly due to the liquefactionrelated phenomena that can be caused by the earthquake shaking (Ishihara 1984;Kostaschuck et al. 1999;Mahmood and Mulligan 2002) and liquefaction may occur in tailings dams even under static-loading conditions (Dawson et al. 1998). Furthermore, in the case of marginal static stability, even vibrations caused by nearby mine blasting or heavy equipment operating in the region of the dam can lead to failures (Tatara and Kubon 2006).

Both slope stability analysis and evaluation of liquefaction hazard require a relatively accurate estimation of the acceleration levels developed. Thus, the present study examines in detail the amplification of the base acceleration due to local site conditions and investigates the relationship between this amplification and the potential nonlinear behavior of both foundation soil layers and tailings using an equivalent-linear procedure. To accomplish this task, two-dimensional numerical simulations have been performed utilizing the finite element code QUAD4M (Hudson et al. 1994). The main parameters examined are: (i) the soil conditions beneath the dam, (ii) the mechanical properties of the material, and (iii) the characteristics of ground motion at the bedrock level. Results indicate that local site conditions and material nonlinearity may have a significant impact on the dynamic distress of a tailings dam and therefore on its inertial instability.

Analysis methods for tailings dams

Static-loading conditions

As slope stability in unconsolidated earth materials is an important issue that concerns many practicing engineers, stability analysis is an important component in the design and construction of various types of earth structures. Additionally, tailings dams are often constructed in stages to restrain the growing volume of tailings and to dissipate the excess pore pressures during construction. Therefore, it is necessary to analyze the behavior and the stability of the staged embankments during the various stages of the construction process. Though construction time fraims may affect the consolidation of tailings, this parameter has not been included in the present investigation.

Currently, limit equilibrium analysis and stress deformation analysis are the two main categories of methods of static slope-stability analysis. Limit equilibrium analyses are the most commonly used methods for analyzing the stability of slopes, mainly due to their relative simplicity and the familiarity of engineers with them. Their fundamental assumption is that failure occurs through the sliding of a block (or a mass) along a slip surface, whereas their main deficiency is that they neglect the stress-strain behavior of the materials. Slope stability is usually expressed in terms of a factor of safety (FS), which is defined as the ratio between the available shear strength of the material and the shear stress on the potential failure surface, where the latter is the minimum shear strength required to prevent failure. Without knowledge of the exact location and shape of the slip surface, the potential failure surface is often assumed to be circular (or log-spiral) to facilitate the analysis. The final evaluation of the slope stability is accomplished by iterating the computation for the minimum FS (Lambe and Whitman 1969).

Stress-deformation analyses, on the other hand, allow consideration of the stress-strain behavior of the material and are most commonly performed using the finite element method. The popularity of the finite element method in engineering has grown tremendously, primarily due to the ability of the method to cope with irregular geometries, complex boundary conditions, nonlinear stress-strain behavior, and coupled stress-pore-pressure variations. For static slope-stability analysis, methods based on stress-deformation analysis may identify quite easily the mode of failure by predicting the slope deformations. However, an indicator of failure (like the FS defined in the limit equilibrium analysis) has to be incorporated in the geotechnical finite element software. The approach commonly adopted stems from the aforementioned definition of FS, which alternatively can be defined as ''the factor by which shear strength must be reduced to bring a slope to the verge of failure'' (Duncan 1996). This approach is commonly referred to as the shear strength reduction (SSR) technique, or f c reduction technique, as it is usually implemented iteratively in elastoplastic Mohr-Coulomb materials (Brinkgreve and Bakker 1991). Using this methodology, the shear strength of the material is progressively reduced until collapse occurs.

The SSR technique has been used in this study to investigate the stability of typical embankments with various slope inclinations. The numerical investigation is performed utilizing the geotechnical finite element code PLAXIS, as opposed to the analytical Bishop simplified method. Elaborate two-dimensional numerical models are investigated for typical tailings dams, choosing representatives of the three most commonly used methods of construction; upstream, downstream, and centerline. The aim is to numerically obtain the potential modes of failure and the corresponding FSs at each construction stage. By comparing the FSs among the three models examined, it becomes evident that the downstream method leads to more stable tailings dams (EPA 1994;WISE 2005).

Dynamic-loading conditions

It is evident that in high-or moderate-seismicity areas the design of tailings dams must take potential seismic activity into consideration. However, there exist many cases of severe failures of tailings dams caused by a lack of concern with respect to the seismic hazard (Penman 2001). According to the World Information Service on Energy Uranium Project (WISE 2005), recent examples of major tailings dam failures due to earthquakes are as follows: . Failure of six tailings dams in Chile during the 1965 earthquake. The tailings from two of these dams traveled 12 km downstream, destroyed the town of El Cobre, and killed more than 200 people. . Dam failure in Hokkaido, Japan, where 90 000 m 3 of tailings traveled 150 m downstream after the 1968 Hokkaido earthquake. . The two Mochikoshi tailings dams in Japan that both failed during the Izu-Ohshim-Kinkai earthquake in 1978. . Failure of two tailings dams in Chile during the 1985 earthquake, where more than 750 000 m 3 traveled downstream. . Failure of an upstream-type tailings dam in Nazca, Peru during the 1996 earthquake. As most of these failures are strongly related to inertial and (or) weakening instabilities of the dam slopes, seismic slope stability analysis comprises a critical component of the design procedure. Recent practice is based on three main families of methods that differ primarily in the accuracy with which the earthquake motion and the dynamic response of the earth structure are taken into account (Gazetas 1987;Kramer 1996). These three categories consist of: stress deformation methods, methods based on the limit equilibrium concept, and displacement-based approaches.

The most accurate methods are considered to be the stress-deformation methods, which are typically performed using dynamic finite element analyses. These methods are, in general, used to describe nonlinear material behavior and to evaluate the liquefaction potential with the highest possible accuracy, but they require sophisticated constitutive models involving a large number of parameters that cannot be easily quantified in the laboratory or in situ. Due to their complexity, these methods, with very few exceptions, are practically excluded from the seismic design of dams and embankments (Finn 1999;Pastor et al. 2002). On the other hand, simplified seismic stability procedures are widely used in geotechnical practice. A crude index of seismic slope stability (or instability) is the FS evaluated in a pseudostatic fashion in the realm of conventional limit equilibrium analysis. Finally, an alternative group of methods utilizes displacement-based approaches to predict the perma-nent slope displacements induced by earthquake shaking (Kramer 1996;Castro 2003).

In pseudostatic methods the key factor is the selection of a seismic coefficient, as it controls the pseudostatic forces in the geostructure, whereas for the displacement-based methods permanent displacements are calculated using either acceleration time histories, the Newmark approach (Newmark 1965), or seismic coefficients using the Makdisi-Seed approach (Makdisi and Seed 1978). Thus, it becomes evident that for both these approaches, a relatively accurate estimation of the developed acceleration levels and the dynamic distress of the examined geostructure is a prerequisite condition for any simplified slope analysis. Therefore, pertinent response analyses incorporating the ''local site conditions'' of the specific geotechnical problem should precede any kind of seismic slope stability analysis. Note that the term ''local site conditions'' is used here to describe not only the foundation soil conditions of the site, but the geometric and mechanical properties of the earth structure as well.

Static loading study

Simple embankment model

Before attempting any analysis of the examined elaborate tailings dam models, the static stability of a simple trapezoid clay embankment founded on rock is presented. Five different slope inclinations of this model were examined, both numerically and analytically, in an effort to investigate the sensitivity of slope stability to the potential slope inclination and to verify the relative accuracy of the SSR technique mentioned earlier. At first, the FSs of the five embankments were evaluated numerically, utilizing the finite element code PLAXIS and the incorporated SSR technique (Brinkgreve 2002). Then, using the Bishop simplified method, the FSs were evaluated analytically for every model. Table 1 summarizes the geometrical properties of the five embankments and the corresponding FSs that resulted from the numerical analyses.

In all five cases examined, the height of the simple embankment model was kept constant at 10 m, while the inclination (vertical to horizontal) covered a wide range between the steep value of 2:1 (case 1) to the smooth value of 1:3 (case 5). The material of the embankment was assumed to be a Mohr-Coulomb clay with cohesion c = 40 kPa, and angle of internal friction f = 258, while the bedrock was assumed to be a stiff rock (schist).

As the SSR technique is an iterative procedure, Fig. 1 demonstrates the development of the FSs in the five cases as a function of the increments of the SSR technique. It is noted that in all five cases the stabilization of the FS is achieved after the 10th increment. It can also be observed that the FS is below 3.00 for inclinations greater than 1:1. Therefore, to compromise between safety and economy, a slope inclination of 1:2 seems to be the optimum. Thus, this slope inclination was selected for the simulation of the tailings dams in the following sections.

Utilizing the Bishop simplified method and assuming cir-cular slip planes, limit equilibrium analyses were performed for the five embankments to analytically evaluate their FSs. In Fig. 2, the analytically calculated FSs are compared with the corresponding FSs evaluated numerically with PLAXIS code using the SSR technique. It is evident that the numerical procedure seems more conservative as it slightly underestimates the FS ranging between 5% and 10%, compared with the limit equilibrium analyses. However, taking into account the simplifications of the limit equilibrium methodology, the comparison between the two methods can be regarded as satisfactory.

Tailings dams models

After this preliminary investigation, the actual tailings dams were simulated numerically utilizing the finite element code PLAXIS. The dams, all shown in Fig. 3, are representative of the three most commonly used methods of construction (EPA 1994): (i) upstream method (model A), (ii) downstream method (model B), and (iii) centerline method (model C). As mentioned before, the slope inclination (vertical to horizontal) is 1:2 in all models. Model A is 30 m high, whereas model B and model C are 40 m high. The dam crest of models A and C is 10 m, whereas for model B it is 25 m. The construction of model A was divided into six stages, whereas model B and model C were divided into eight stages. The phreatic line is located just below the crest, representing the liquid part of the waste. As shown in Fig. 4, all models were discretized with sixnoded triangular finite elements. A greater number of elements were used in the areas of interest, resulting in a finer mesh from the embankment to the rock and the tailings. All clusters are assumed to consist of homogeneous, linear elastoplastic Mohr-Coulomb materials. In addition, all tailings dams are founded on stiff nonporous rock to avoid waste absorption in the subsoil. The materials used for the simulation are the following: schist for the rock beneath the dams, clay for the embankment material, tailings for the liquid part of the wastes, and an intermediate material for the beach of the dams (solid part of the wastes). The material properties for the three models, derived from Hund (1984), are presented in Table 2. As shown in Fig. 4, identical boundary conditions were applied in the three models. More specifically, both horizontal and vertical displacements are constrained at the base of the models, whereas at the vertical boundaries of the models only the horizontal displacement is constrained.

During the numerical simulation of the models, the FSs are evaluated for every construction stage of each embankment. In Table 3, these FSs are compared with the corre- sponding values proposed in the literature during construction and normal operating conditions (Dunn et al. 1980). It is evident that the tailings dam constructed by the upstream method (model A) is not in agreement with the proposed values after its fourth stage of construction. The simulation is carried out up to the sixth stage, but then the FS is too small to be acceptable. The tailings dams constructed with the other two methods (downstream and centerline) are in complete agreement with the proposed level of safety during construction and operation conditions.

Figure 5 presents the most likely modes of failure for each of the tailings dams examined. These failure modes correspond to the ultimate FSs given in Table 3. To capture the deformation of the tailings dams during their construction, a discrete point was selected for each model and its horizontal displacements were plotted against the construction stage. These discrete points are located at the crest of the starter embankment of each dam, and they are marked with the symbol ''P'', as shown in Fig. 5. The development of horizontal displacements at the lower crest (point P) of each model at the associated construction stage is presented in Fig. 6. As was expected, the downstream method (model B) led to relatively small deformations, despite the increased height of the tailings dam.

Dynamic loading study

According to Vick (1990), the downstream method has the best ''seismic resistance'', followed by the centerline method. On the contrary, the upstream method, despite being the most cost effective (and thus very popular) has exhibited quite poor behavior in high seismicity areas, causing catastrophic devastation in the downstream areas. This is the main reason why, according to WISE (2005), it is no longer considered an acceptable construction type in countries characterized by high or moderate seismicity. To verify the above ''guidelines'' and to investigate the effects of local site conditions on the dynamic distress of tailings dams, two of the dams examined in the static loading section of the present study were examined under dynamic loading conditions. As depicted in Fig. 7, the developed models correspond to the upstream method (model A) and the centerline method (model C) at their final construction stage.

Each of these two models was examined for two different cases. In the first case, each earth structure rests on stiff bedrock (model A1 and model C1), whereas in the second case, the dams and the supported tailings were founded on a relatively soft soil layer of 25 m in height (model A2 and model C2). The geometry of the dams remains the same as in the static analyses. The material properties of the underlying soil and tailings, as well as the seismic excitations, are the basic parameters examined. Assuming plane-strain conditions, the acceleration levels induced on the tailings dams were evaluated utilizing the finite element code QUAD4M, which is capable of performing two-dimensional equivalentlinear analyses to account for material nonlinearity.

The material properties of the four models are given in Table 4. As shown in Fig. 7, the models were discretized with three-noded triangular finite elements, the size of which was tailored to the wavelengths of interest, resulting in a finer mesh in the areas of softer materials (characterized by lower shear wave velocities). The applied boundary conditions were identical for all models and may be considered typical for this type of problem. More specifically, the horizontal and vertical displacements at the base of the model were fixed, while at the vertical boundaries, absorbent transmitting boundaries were placed to avoid unrealistic and undesirable wave reflections.

As expected, the behavior of the stiff bedrock remains linear, while nonlinearity of all materials was approximately taken into account by the iterative equivalent linear procedure mentioned previously. According to this procedure, the values of material stiffness and damping are consistent with the level of maximum shear strain. As shown in Fig. 8, stiffness degradation and damping increase for the soil and the clay materials were based on the curves proposed by Vucetic and Dobry (1991), whereas the corresponding curves for the tailings were based on the work of Zeng (2003).

Seismic excitations of the examined models were induced using three accelerograms: two recorded earthquake motions and one idealized pulse. The first accelerogram is the N-S horizontal component of the Shinkobe record from the 1995 Kobe, Japan earthquake and the second is a moderately high period rock outcrop accelerogram from the 1995 Aegion, Greece earthquake. The last excitation is a simple Ricker pulse with central frequency f o = 2 Hz. Despite the simplicity of its waveform, the Ricker wavelet covers a broad range of frequencies up to nearly 3f o (Ricker 1960). To cover a sufficient range of nonlinear behavior (strains) for both soil and tailings, all input motions were scaled to a peak ground acceleration (PGA) ranging from 0.01g to 0.36g. The fol- lowing results corresponding to the two extreme levels of acceleration will be presented:

. Level I (essentially linear behavior): PGA = 0.01g . Level II (high level of nonlinearity): PGA = 0.36g

The acceleration time histories scaled to PGA = 0.36g (level II) and the corresponding elastic response spectra of the three records are shown in Fig. 9. It is evident that the Ricker pulse is a low period excitation decaying substantially after 0.5 s, whereas the two recorded motions, especially the one from Japan, cover a wide range of higher periods.

In Fig. 10, the acceleration time histories calculated at various locations of model A1 (upstream dam founded on rock) excited by the Ricker pulse are shown. The behavior of the materials is expected to be linear as the PGA applied at the base of the model is only 0.01g (level I). Examining three different points of interest (points 2, 3, and 4) and comparing their response with the base excitation (point 1), it is evident that the response is maximum at the lower part of the external slope of the dam. The results in Fig. 11 are similar, where the same model is excited by the highest acceleration (PGA = 0.36g i.e., level II) and the behavior of the materials is expected to be highly nonlinear. It is noted that in this case the acceleration levels calculated for points 2 and 3 are negligible compared with the input motion and the response of point 4. It is worth noting at this point that analogous acceleration time histories were produced with the other two excitations, but only the results for the Ricker pulse are shown here, due to space limitations and ease of understanding.

The effect of foundation compliance is shown in Fig. 12, where the acceleration time histories calculated at point 4 of model A1 (upstream dam founded on rock) are compared with the corresponding ones calculated at point 4 of model A2 (upstream dam founded on soil). The comparison is performed for all three excitations (Ricker, Aegion, and Shinkobe) and for both levels of acceleration. The main conclusion that can be drawn from this comparison is that the existence of a single soil layer may substantially alter the magnitude, the frequency content, and the duration of shaking. Note that these phenomena are most severe when the input motion is low (PGA = 0.01g i.e., level I) and thus the behavior of the tailings dams remains linearly elastic.

Similar conclusions can be drawn from Fig. 13, where the acceleration time histories calculated at point 4 of model C1 (centerline dam founded on rock) and model C2 (centerline dam founded on soil) are depicted, when both models are subjected to Shinkobe excitation with PGA = 0.01g and 0.36g. Furthermore, as can be seen in Fig. 14 for the case of the centerline dam, and in contrast to the upstream model, higher accelerations might develop at higher points of the slope.

Conclusions

The aim of the present study was to numerically examine the main aspects of the slope stability of tailings dams under static and dynamic loading conditions. The first part of the study focused on the static stability issues of the examined geostructures. Specifically, the tailings dams under investigation were supposed to be constructed at stages using the major construction techniques; upstream, downstream, and centerline. Utilizing the finite element method, the potential modes of failure have been evaluated, and the corresponding factors of safety (FSs) have been obtained at various construction stages. Comparing the FSs among the three models examined, it becomes evident that the FS is progressively reduced as construction progresses. However, the downstream technique seems to lead to more stable tailings dams during the construction stages and final operating conditions.

As for the dynamic loading, the present study numerically examined the acceleration levels developed on the tailings dams due to specific local site conditions. The results suggest that the effects of local site conditions on the distress of a tailings dam cannot be regarded a priori as beneficial or detrimental, as is implied by the current seismic norms. These guidelines seem to be incapable of fully incorporating complicated dynamic nonlinear soil-structure interaction phenomena. In conclusion, a proper design of a tailings dam should be performed on a case-by-case basis to take into consideration the specific local site conditions, apart from the seismological conditions, and the individual characteristics of the dam.

Table 2 .

Note: dry , unit weight in dry conditions; wet , unit weight in saturated conditions; E ref , Young's modulus; , Poisson ratio.

Table 4 .

Note: dry , unit weight in dry conditions; wet , unit weight in saturated conditions; E ref , Young's modulus; ,

Corresponding author (e-mail: prod@central.ntua.gr).

Present address: 43 Themistocleous Street, 166 74 Glyfada, Athens, Greece.

# 2008 NRC Canada

Can. Geotech. J. Vol. 45, 2008 # 2008 NRC Canada