Numerical Analyses of Earthquake Induced Liquefaction and Deformation Behaviour of an Upstream Tailings Dam

Abstract

Much of the seismic activity of northern Sweden consists of micro-earthquakes occurring near postglacial faults. However, larger magnitude earthquakes do occur in Sweden, and earthquake statistics indicate that a magnitude 5 event is likely to occur once every century. This paper presents dynamic analyses of the effects of larger earthquakes on an upstream tailings dam at the Aitik copper mine in northern Sweden. The analyses were performed to evaluate the potential for liquefaction and to assess stability of the dam under two specific earthquakes: a commonly occurring magnitude 3.6 event and a more extreme earthquake of magnitude 5.8. The dynamic analyses were carried out with the finite element program PLAXIS using a recently implemented constitutive model called UBCSAND. The results indicate that the magnitude 5.8 earthquake would likely induce liquefaction in a limited zone located below the ground surface near the embankment dikes. It is interpreted that stability of the dam may not be affected due to the limited extent of the liquefied zone. Both types of earthquakes are predicted to induce tolerable magnitudes of displacements. The results of the postseismic slope stability analysis, performed for a state after a seismic event, suggest that the dam is stable during both the earthquakes.

Muhammad Auchar Zardari 1 and Hans Mattsson 2 and Sven Knutsson 2 and Muhammad Shehzad Khalid 3 and Maria V. S. Ask 2 and Björn Lund 4

Academic Editor:Rui Wang

1, Department of Civil Engineering, Quaid-e-Awam University of Engineering Science and Technology, Nawabshah, Sindh, Pakistan

2, Department of Civil, Environmental and Natural Resources Engineering, Luleå University of Technology, 97187 Luleå, Sweden

3, Department of Urban Management, Kyoto University, Kyoto, Japan

4, Department of Earth Sciences, Uppsala University, Uppsala, Sweden

Received 13 October 2016; Accepted 10 January 2017; 15 February 2017

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

This paper presents finite element analyses of dynamic behaviour of the Aitik tailings dam in northern Sweden in response to the impact of an earthquake. It is one of the first attempts to investigate the stability of a tailings dam in Sweden, using numerical methods to study dynamic loading. The purpose of this study is to evaluate liquefaction potential, permanent deformations, and the stability of the dam under both a commonly occurring Swedish earthquake and a more extreme event.

In Sweden, approximately one earthquake with magnitude above three occurs every year [1]. Using the available database of earthquakes that have occurred in Sweden since the 19th century, there are indications of approximately one magnitude 5 event every century and one magnitude 6 event every millennium [2]. Guidelines on tailings dams in Sweden [3] suggest that a tailings dam needs to be stable for at least 1000 years. In such a long time perspective, it is therefore relevant to consider the possibility of a large magnitude earthquake [3].

The Aitik tailings dam is mainly raised using the upstream construction method [4]. In this raising method, dikes are sequentially constructed on the previously deposited tailings that could be in a loose and saturated state. In this type of dam structure, pore water in saturated tailings may not drain as rapidly as seismic shaking occurs. As a result, excess pore pressures increase cumulatively with each cycle of loading. If the excess pore pressures increase to such an extent that effective stresses are reduced to near zero, it could lead to liquefaction (see, e.g., [5]). In the liquefied state, large deformations could occur due to low stiffness of the tailings. The deformations progressively increase with each cycle of loading and may lead to failure of the dam. It is relevant to mention here that several tailings dams, located in different parts of the world, have failed due to seismic liquefaction (see, e.g., [6]).

Seismic liquefaction in a tailings dam is a complex phenomenon. Therefore, it is appropriate to investigate seismic behaviour of the tailings dam with advanced numerical software, based on, for example, the finite element method. To the authors' knowledge, there are just a few case studies of tailings dams published, where such advanced numerical tools have been used (see, e.g., [7-10]). Some similar case studies regarding seismic behaviour of geotechnical structures (mainly embankment dams) could be found in the literature (see, e.g., [11-17]). More case studies of this type are needed in order to gain more experience of numerical modelling of seismic liquefaction in tailings dams. In this study, numerical modelling of seismic liquefaction of the Aitik tailings dam is performed with the finite element program PLAXIS 2D [18], using a recently implemented advanced constitutive model called UBCSAND (see [19, 20]).

The above information on seismicity in Sweden suggests that it is important to analyse liquefaction potential and stability of the Aitik dam by taking into account both normal and extreme seismic conditions in Sweden. The modelling technique presented in this paper is not unique for the Aitik tailings dam and can be used for other tailings dam facilities.

1.1. Seismicity of Northern Sweden

Intraplate seismicity is characterized by few, dispersed, and irregularly occurring earthquakes. In general, these earthquakes tend to be relatively small in magnitude, although there are many examples of devastating intraplate earthquakes. England and Jackson [21] report that ~100 out of the ~130 earthquakes that have resulted in a thousand or more casualties over the last 120 years worldwide took place in continental interiors. Therefore, they conclude that research efforts should focus on the threat of unanticipated earthquakes located in intraplate regions.

Sweden is located within the Baltic Shield, far away from plate boundaries. Figure 1 shows that the micro-earthquake activity in northern Sweden is surprisingly high [22-24] and focused along postglacial faults (e.g., [25]). It has been proposed that these faults hosted earthquakes with magnitudes between 7 and 8 when the ice disappeared after the latest glaciation, some 10000 years ago (e.g., [26]). The roles of the ice, in situ effective stresses, and the nature of the crust have been analysed in terms of fault instability (e.g., [27, 28]). The Swedish National Seismic Network [1] has recorded 9 earthquakes larger than magnitude 3 in northern Sweden since 2000, and in 1983, a magnitude 4.1 earthquake occurred in Solberga in the southern part of the region [29]. The largest earthquake with waveform data, from south of Skellefteå, was used in this study.

Figure 1: Observed earthquakes from year 1398 and onwards from the joint Nordic catalogue at Helsinki University FENCAT (2011) and from 2000-2011 by the Swedish National Seismic Network (SNSN) at Uppsala University (e.g., [2]). The location of the Aitik tailings dam is indicated by the circle. Red dots show the epicentral location of the earthquakes. Blue rectangles show the location of permanent Swedish seismic stations. Postglacial faults are highlighted by numbers: (1) the Pärvie fault, (2) the Merasjärvi fault, (3) the Lainio-Suijavaara fault, (4) the Stuoragurra fault, (5) the Suasselka fault, (6) the Pasmajärvi/Venejärvi fault, (7) the Lansjärv fault, and (8) the Burträsk fault. Modified from [22].

[figure omitted; refer to PDF]

Near the Pärvie fault, the longest known postglacial fault in the world (cf. (1) in Figure 1), a number of events have been reported by people living near the fault, and the instrumental records contain earthquake data of up to magnitude 4 [30]. Further to the southeast, historical data from the Gulf of Bothnia area derived from newspaper reports suggest that nine earthquakes occurred during the period from 1757 to 1909 (e.g., [31, 32]). The majority of these events have estimated magnitudes of 3.5 to 3.9, with a few events having magnitudes above 4, and the largest earthquake with magnitude 4.7 recorded in the Gulf of Bothnia in 1898 (e.g., [31, 32]). Mäntyniemi [32] proposed that earthquake activity in this region may occur in clusters over time. The largest known historical earthquake near the study area had a magnitude of 5.3 and is proposed to have occurred in 1819 in Lurøy, Norway [33].

The Aitik tailings dam is located within about 20 km and 100 km from the large underground mines of LKAB in Malmberget and Kiruna, respectively. As the mining progress to deeper levels, mining-induced seismicity has been increased. LKAB has established local seismic networks to monitor mining-induced seismicity, and the larger events are also recorded by the SNSN [1]. Two magnitude 2.9 events were recorded in 2007 and 2008 in Malmberget and Kiruna, respectively.

2. The Aitik Tailings Dam

Aitik is the largest open pit copper mine in Sweden. This mine is situated near Gällivare in northern Sweden. The annual production of the mine was expanded to 34.3 million tonnes of ore yearly in 2012.

Tailings are the by-products which are generated in a mine processing plant when minerals are extracted from the ore. The tailings, produced at Aitik, are transported in slurry form to the tailings impoundment. The tailings slurry is then discharged through spigots in the pipeline along the dam embankments. The tailings which are deposited close to the embankments can be classified as silty fine sands according to the European standard [34] and as silty sands in accordance with the unified soil classification system [35].

Figure 2 shows Aitik tailings impoundment which occupies an area of 13 square kilometres. The impoundment comprises four dams named A-B, C-D, E-F (including E-F2 extension), and G-H. The clarification pond is situated downstream of dam E-F.

Figure 2: Aerial view of Aitik tailings dam and impoundment [36].

[figure omitted; refer to PDF]

Stability of the dam E-F, as well as the corner between the dam parts E-F and G-H, has already been analysed under static loading conditions for gradual raising [37, 38]. Among the various dams at the Aitik impoundment, the dynamic finite element analyses have been performed on the dam E-F. The consequences of a failure of the dam E-F are expected to be more serious compared to failures of other dam parts at the Aitik impoundment. The reason is that a failure of the dam E-F may also lead to a failure of the dam I-J which is located downstream of the clarification pond (Figure 2). Such a failure would lead to discharge of tailings and waste water which might affect the surrounding environment.

3. Finite Element Model of the Dam

The construction of a finite element model basically consists of the following steps before calculation: specify the type of analysis and the type of elements, input of geometry, input of boundary conditions, input of material properties, mesh generation, input of initial conditions, and definition of calculation phases.

Figure 3 presents the finite element model of the Aitik tailings dam section named E-F. The different material zones are illustrated in the model. The tailings material zones 2, 3, 5, 6, 7, and 8 are categorized as silty sands. The bottom of the dam is at an elevation of 340 m. It has been assumed, based on mine production plans, that the tailings are deposited in the impoundment at a rate of three meters per year from a level of 376 m. The elevation of the dam in 2012 was about 390 m and it is planned to raise the dam to the level of 409 m in 2018 (cf. Figure 3). The gradual raising of the dam, from level 376 m to 409 m, was simulated in eleven stages. Each stage consisted of a raising phase over 10 days and a consolidation phase over 355 days.

Figure 3: Finite element model of dam E-F. Dimensions are in meters.

[figure omitted; refer to PDF]

The phreatic level of the dam was evaluated using Piezometers. The phreatic level for first raising is shown in Figure 4. It can be seen that the phreatic level is on the surface of tailings and then it is parallel to the downstream slope of the dam. The same trend of the phreatic level was utilized in further raising of the dam (see Figure 5).

Figure 4: Location of phreatic level at first raising of the dam.

[figure omitted; refer to PDF]

Figure 5: Location of phreatic level at eleventh raising of the dam.

[figure omitted; refer to PDF]

A finite element program for two-dimensional models PLAXIS 2D [18] was utilized for numerical analyses of the dam. A plane strain condition was adopted which is a suitable assumption for both static and dynamic analyses of the long and straight dam section E-F.

The numerical analyses were conducted in two steps. In the first step, static analyses were performed to simulate the sequential raising of the dam. For this purpose, coupled deformation and consolidation analyses [18, 39, 40] were carried out to estimate the development and dissipation of excess pore pressures, as a function of time, in the dam.

In the second step, dynamic analyses were performed directly after the 11th raising of the dam (elevation 409 m, cf. Figure 3) to evaluate the seismic behaviour of the dam for two different earthquake events. In both steps, slope stability analyses were conducted to compute safety factors and to observe the potential failures zones in the dam. Slope stability analyses for static loading (i.e., gradual raising) of the dam are described in more detail in Ormann et al. [37].

Figure 6 illustrates the finite element mesh adopted in the analyses. Fifteen noded triangular elements were utilized in the mesh. These elements give a fourth-order (quartic) interpolation for displacements [18]. The elements used in the mesh were 14077 in number. Average size of an element was about 1.88 m. The finite element mesh was sufficiently refined, as suggested by Kuhlemeyer and Lysmer [41], to ensure that seismic waves are realistically transmitted through the numerical model. The horizontal width of the numerical model was large enough in order to minimize the influence of seismic waves that might reflect back from the left vertical boundary.

Figure 6: Finite element mesh of dam E-F.

[figure omitted; refer to PDF]

In both the static and the dynamic analyses, groundwater flow and/or consolidation can occur through all the boundaries except the left vertical boundary and the base. In the static analyses, displacements along the left vertical boundary were restrained in the horizontal direction, while displacements at the base were restricted in both the horizontal and the vertical directions. In the dynamic analyses, (i) an absorbent boundary was utilized on the left vertical boundary to absorb the seismic waves, which otherwise would reflect back into the dam body [18], and (ii) a prescribed horizontal displacement factor of 0.01 m was applied at the bottom of the model, and vertical displacements are set to zero there. The prescribed horizontal displacement factor of 0.01 m at the bottom of the model is scaled to the actual magnitudes of the displacements that vary with respect to time during an earthquake. In the finite element program PLAXIS, an input for a dynamic load can be given with either displacements or velocities or accelerations produced due to a seismic shaking [18]. The velocities or accelerations are then converted into displacements in the computation process.

The foundation of the dam lies on a dense moraine layer. There is bedrock below the moraine layer. The average depth of the moraine deposit is about 8 m. It was assumed that the foundation of the dam on the moraine deposit is rigid and nonliquefiable; so the moraine layer was not included in the analyses.

It is to be noted, in the static analyses, rockfill berms were placed on the downstream side of the dam to increase slope stability during various raising (see [37]). The rockfill berms P, Q, R, S, T, U, V, and W (Figure 7) were placed at the start of the second, fourth, fifth, sixth, seventh, eighth, ninth, and tenth raising, respectively.

Figure 7: Placement of the rockfill berms P, Q, R, S, T, U, V, and W.

[figure omitted; refer to PDF]

3.1. Constitutive Models and Material Properties

In the static analyses, the Mohr Coulomb model was used to represent the constitutive behaviour of all the material zones in the dam. For the dynamic analyses, (i) the UBCSAND model [19, 20] was chosen to simulate the constitutive behaviour of the tailings material zones, and (ii) the Mohr Coulomb model was applied to the filter, moraine, and rockfill zones (cf. Figure 3) which are assumed to be nonliquefiable. A detailed description of the two selected constitutive models is presented in the following sections.

3.2. Mohr Coulomb Model

The Mohr Coulomb (MC) model is a simple linear elastic perfectly plastic model which requires five input parameters. These input parameters are Young's modulus, Poisson's ratio, friction angle, cohesion, and dilatancy angle. These parameters can be obtained from basic soil tests [18]. In the finite element program PLAXIS, some additional material parameters such as unit weights and hydraulic conductivities are also needed.

The values of Young's modulus, friction angle, cohesion, unit weights, and hydraulic conductivities for tailings material zones were determined from field and laboratory tests. Appropriate values of these parameters for moraine, filter, and rockfill zones were obtained from the available literature [3, 42-45]. The abovementioned parameters for all the material zones in the dam are presented in Table 1. It is to be remembered that for the dynamic analyses the MC model is just utilized for the nonliquefiable zones 1, 4, 9, 10, and 11 (cf. Figure 3).

Table 1: Parameters of the Mohr Coulomb model [3, 42-45].

Material zone	Material type	γ u n s a t kN/m3	γ s a t kN/m3	k x m/s	k y m/s	E kN/m2	c [variant prime] kN/m2	[straight phi] [variant prime] °
1	Moraine (initial dike)	20	22	1 × 10-7	5 × 10-8	20000	1	35
2	Soft sand tailings	18	18	1 × 10-7	1 × 10-8	9800	6	18
3	Layered sand tailings	17	19	5 × 10-7	5 × 10-8	9300	10	22
4	Moraine (dikes)	20	22	5 × 10-8	1 × 10-8	20000	1	37
5	Compacted sand tailings	16	19	1 × 10-6	1 × 10-7	8800	13	26
6	Soft sand tailings	18	18	1 × 10-7	1 × 10-8	3000	6	18
7	Compacted sand tailings (dikes)	16	19	1 × 10-6	1 × 10-7	7200	13	26
8	Layered sand tailings	17	19	5 × 10-7	5 × 10-8	3900	10	22
9	Filter	18	20	1 × 10-3	1 × 10-3	20000	1	32
10, 11	Rockfill (berms + external erosion protection)	18	20	1 × 10-1	1 × 10-1	40000	1	42

Note: _γunsat is the unit weight above phreatic level, _γsat is the unit weight below phreatic level, _kx is the hydraulic conductivity in horizontal direction, _ky is the hydraulic conductivity in vertical direction, E is Young's modulus, ^{c[variant prime]} is the effective cohesion, and ^{[straight phi][variant prime]} is the effective friction angle.

The value of Poisson's ratio is assumed to be 0.33 for all the materials in the dam. This value of Poisson's ratio is considered to be suitable for this type of analysis [18]. For all material zones in the dam, the value of the dilatancy angle is assumed to be zero. This is a convenient assumption. The reason is that unrealistically high negative pore pressures may occur due to the use of a positive value of the dilatancy angle, whereas unreasonably large positive pore pressures may develop if the value of the dilatancy angle is negative [18].

3.3. UBCSAND Model

UBCSAND is an effective stress elastic-plastic model which is capable of simulating the liquefaction behaviour of sands and silty sands under seismic loading [20]. The name UBCSAND implies that this model was developed at the University of British Columbia for prediction of liquefaction behaviour of sand. An earlier version of the UBCSAND model [46, 47] was used in a case study of dynamic analyses of Mochikoshi tailings dam, in Japan, and the results of these analyses were consistent with the observed failure pattern of the dam induced due to seismic liquefaction (see [9, 10]). The UBCSAND model [20], with some modifications, has been implemented as a user defined soil model in the finite element program PLAXIS [19]. The PLAXIS version of the UBCSAND model is utilized in this study.

The material parameters required for the UBCSAND model are constant volume friction angle _[varphi]cv , peak friction angle _[varphi]p , cohesion c, elastic shear modulus number ^KGe , plastic shear modulus number ^KGp , elastic bulk modulus number ^KBe , elastic bulk modulus index me, elastic shear modulus index ne, plastic shear modulus index np, failure ratio _Rf , atmospheric pressure _PA (100 kPa), tension cut-off _σt (0 kPa), densification factor _fachard , corrected standard penetration test value of the soil _{(N1 )60} , and postliquefaction factor _fpost .

The constant volume friction angle, the peak friction angle, and cohesion were evaluated from direct shear tests on tailings material. The values of the constant volume friction angle were assumed to be the same as the values of the effective friction angle evaluated for the MC model (cf. Table 1). The values of the parameter _{(N1 )60} for tailings were obtained by using correlations (see [48]) with the results of cone penetration tests which were performed on the Aitik deposit. The elastic shear modulus number, the plastic shear modulus number, and the failure ratio were obtained by curve fitting with the direct shear test results. The elastic bulk modulus number was related to the elastic shear modulus number using Poisson's ratio [19]. For geotechnical problems which involve unloading and reloading (e.g., dynamic loading), a low value of Poisson's ratio is recommended to be used in advanced soil models [18]. Poisson's ratio of 0.12 was utilized to evaluate stiffness parameters of the UBCSAND model. This value of Poisson's ratio is considered to be suitable for this type of analysis [20, 49]. For Poisson's ratio of 0.12, the elastic bulk modulus number is approximately equal to the elastic shear modulus number. Therefore, the same value is used for the elastic shear modulus number and the elastic bulk modulus number. The elastic shear modulus index, elastic bulk modulus index, and plastic shear modulus index were assigned as 0.5, 0.5, and 0.4, respectively (see [20]). Appropriate values of the densification factor and the postliquefaction factor were taken as 0.3 and 0.01, respectively (see [19, 50]). Table 2 presents the following input parameters: the peak friction angle, the elastic shear modulus number, the plastic shear modulus number, the failure ratio, and the corrected standard penetration test value of the soil.

Table 2: Some input parameters of the UBCSAND model.

Material zone	Material type	[straight phi] p [variant prime]	K G e	K G p	R f	( N 1 ) 60
2, 6	Soft sand tailings	18.6°	120	36	0.95	6
3, 8	Layered sand tailings	22.6°	120	36	0.90	6
5, 7	Compacted sand tailings	26.8°	120	72	0.86	8

Note: ^{[straight phi]p[variant prime]} is the peak friction angle, ^KGe is the elastic shear modulus number, ^KGp is the plastic shear modulus number, _Rf is the failure ratio, and _{(N1 )60} is the corrected standard penetration test value of the soil.

It is to be noted that, due to unavailability of cyclic direct shear and/or cyclic triaxial test results, suitable values of the damping ratio, densification factor, and postliquefaction factor were assumed based on results of a case study of an upstream tailings dam presented in the literature [9, 10]. It would of course be better if these parameters had been evaluated from cyclic soil tests on Aitik tailings material. However, it is believed that the results of the analyses give reasonable estimates of the liquefaction potential and the stability of the dam under seismic loading conditions.

3.4. Damping

The amplitude of elastic waves in soils may decrease with time. This property of elastic waves is known as material damping (see, e.g., [51]). In dynamic computations performed with the finite element program PLAXIS 2D, material damping occurs to some extent due to viscous properties of soil, friction, and development of plastic strains [18]. However, the hysteretic damping, exhibited by PLAXIS 2D, is generally lower than the actual damping observed in soils. Therefore, an extra damping is required to simulate damping effects in soils in a realistic manner [18]. This additional damping is provided by means of Rayleigh damping which is a function of mass, and stiffness of the soil, and is defined as (see [52]) [figure omitted; refer to PDF] where C, M, and K are the damping matrix, the mass matrix, and the stiffness matrix, respectively; α and β are the Rayleigh coefficients.

The Rayleigh damping coefficients α and β can be calculated as a function of two representative frequencies of the waves which propagate through the soil [18] and an associated damping ratio [53]. The two representative frequencies, 0.3 Hz and 1.8 Hz, were obtained by taking into consideration both the natural frequencies of the dam and the input frequencies of the seismic waves (Figure 8). The natural frequencies were evaluated from a free vibration analysis of the Aitik dam. It can be seen from the figure that the two dominant frequencies of the input motions are 0.3 Hz and 1.8 Hz. The average of these two frequencies is close to the dominant natural frequency (i.e., 1.2 Hz).

Figure 8: Frequency spectrum for free vibration analysis of the dam and two specific earthquakes.

[figure omitted; refer to PDF]

The UBCSAND model is capable of producing damping. However, for nonsymmetric cyclic loading conditions, the UBCSAND model predicts a less amount of damping compared to the damping of soils observed in laboratory tests [20]. As additional damping, 1% of the critical damping (i.e., damping ratio 0.01) was assigned to the tailings material zones whose constitutive behaviour was represented with the UBCSAND model (see [54, 55]). On the other hand, 3% of the critical damping was assigned to the moraine, rockfill, and filter zones where the MC model was utilized (see [54, 55]).

3.5. Input for Ground Motion

The seismic behaviour of the dam was evaluated for one somewhat larger than normal and one extreme earthquake in northern Sweden. For the normal case, an earthquake of local magnitude 3.6 was selected. This earthquake occurred in 2010 outside Skellefteå, in northern Sweden. The input data for this earthquake was obtained from the Swedish National Seismic Network (SNSN). For the extreme case, an earthquake of moment magnitude 5.8 was selected. This earthquake occurred in 2011 in Virginia, eastern United States. It was chosen because the geology of eastern United States is similar to that of northern Sweden and no records of such a large magnitude earthquake exist in Sweden. The data for the Virginia earthquake was downloaded from the U.S. Geological Survey. As the Swedish local magnitude scale is moment based, a comparison between the magnitude scales is straightforward and in this case the Virginia earthquake radiated approximately 14,000 times more energy than the Skellefteå earthquake. From now on, the Skellefteå earthquake will be referred to as the normal earthquake/case and the Virginia earthquake as the extreme earthquake/case.

The input data for both the normal earthquake and the extreme earthquake are shown in Figures 9 and 10, respectively. Both the Skellefteå earthquake and the Virginia earthquake were recorded on the bedrock. The instrument response has been deconvolved such that Figures 9 and 10 show true ground velocity and acceleration for the respective earthquakes. It can be seen that the horizontal ground motions of the normal earthquake and the extreme earthquake are expressed in velocity and acceleration, respectively. The normal earthquake and the extreme earthquake were recorded at a distance of 11 km and 53 km, respectively.

Figure 9: Input ground motion of Skellefteå earthquake (normal case).

[figure omitted; refer to PDF]

Figure 10: Input ground motion of Virginia earthquake (extreme case).

[figure omitted; refer to PDF]

In this study, only the horizontal ground motions (accelerations and/or velocities) were applied and the vertical ground motions were not considered. The reason is (see [56]) that (i) vertical motions cause much less shear stresses compared to those due to horizontal motions and (ii) pore pressures that develop due to vertical motions are smaller than those caused by horizontal motions.

Remember, in the finite element program PLAXIS, a seismic loading can be applied with any of the following quantities: velocity, acceleration, or displacement. In order to reduce computational time, low values of acceleration and/or velocity are filtered out. It is to be noted that the results of the dynamic analyses are not affected by this filtering process.

4. Results

Dynamic analyses were performed for the state directly after the 11th raising of the Aitik dam. This state was chosen because it is considered as critical since excess pore pressures have developed due to the gradual raising of the dam. In the following sections, the seismic behaviour of the dam is analysed in terms of excess pore pressures, liquefaction, earthquake induced permanent deformations, and postseismic stability.

4.1. Excess Pore Pressures

Figure 11 shows excess pore pressures, after 11th raising of the dam, prior to the earthquake shakings. Figures 12 and 13 illustrate the development of excess pore pressures when the dam is shaken with the normal and the extreme earthquake, respectively. It can be observed that the largest initial excess pore pressures, due to the static loading, are of magnitude 95 kPa (Figure 11). As expected, the largest magnitudes of excess pore pressures, due to the seismic loadings, have increased to (i) about 220 kPa for the normal earthquake (Figure 12) and (ii) about 320 kPa for the extreme earthquake (Figure 13). As anticipated, the magnitudes of the excess pore pressures were higher for the extreme earthquake compared to those generated from the normal earthquake. There was an increase in excess pore pressures throughout the dam body. The maximum increase in excess pore pressures was observed at the lower parts of the dam below the embankment dikes (Figures 12 and 13). As expected due to dynamic loading, positive excess pore water pressures developed in contractant zones (soft sand tailings and layered sand tailings) and negative excess pore water pressures generated in dilatant zones (moraine dikes, compacted sand tailings, and rockfill) of the dam.

Figure 11: Excess pore pressures in the dam after 11th raising. Negative sign indicates suction.

[figure omitted; refer to PDF]

Figure 12: Excess pore pressures after the dam is shaken with normal earthquake. Negative sign indicates suction.

[figure omitted; refer to PDF]

Figure 13: Excess pore pressures after the dam is shaken with extreme earthquake. Negative sign indicates suction.

[figure omitted; refer to PDF]

4.2. Liquefaction

Liquefied zones in the dam were identified by evaluating the excess pore pressure ratio _ru which is defined as [20] [figure omitted; refer to PDF] where f stands for function, _ue is the excess pore pressure, and ^{σv[variant prime]} is the vertical effective stress in the dam after seismic shaking; _u0 is the initial excess pore pressure and ^{σv0[variant prime]} is the initial vertical effective stress in the dam prior to earthquake loading.

It is assumed that liquefaction is likely to occur in the dam at such zones where the values of the excess pore pressure ratio are equal to or greater than 0.7 [54]. For the normal earthquake, no liquefaction was observed in the dam with the selected input parameters. For the extreme earthquake, potential liquefied zones (_ru ≥0.7) are illustrated in Figure 14. It can be seen that liquefaction is likely to occur in the deposited tailings below the ground surface. The potential liquefied zone is located about 4 m below the ground surface near the embankment dikes. It can be observed (Figure 14) that the extent of the liquefied zone is limited. It is, therefore, interpreted that this liquefied zone could not affect the stability of the dam. The occurrence of a liquefied zone below the surface is consistent with general observations that liquefaction initiates at the surface and then extends downwards if the seismic shakings generate high excess pore pressures in the deep layers (see, e.g., [57]).

Figure 14: Potential liquefied zones after shaking with extreme earthquake. Dimensions are in meters. For visibility, scale of vertical axis is larger than the horizontal axis.

[figure omitted; refer to PDF]

4.3. Earthquake Induced Permanent Deformations

The horizontal displacements in the dam induced due to the normal earthquake and the extreme earthquake are shown in Figures 15 and 16, respectively. The maximum horizontal displacements of about 0.3 m and 0.6 m are predicted, in the vicinity of the downstream slope, for the normal earthquake, and the extreme earthquake, respectively. As expected, the horizontal displacements are controlled by the rockfill berms (see Figures 15 and 16).

Figure 15: Horizontal displacements in the dam after shaking with normal earthquake. The positive direction is to the right.

[figure omitted; refer to PDF]

Figure 16: Horizontal displacements in the dam after shaking with extreme earthquake. The positive direction is to the right.

[figure omitted; refer to PDF]

Figures 17 and 18 show the predicted vertical displacements after the dam is shaken with the normal earthquake and the extreme earthquake, respectively. It can be seen that the maximum vertical displacements of about 0.3 m and 0.4 m mainly occurred on the crest of the embankment, for the normal earthquake, and the extreme earthquake, respectively. The vertical displacements are tolerable because they are much smaller than the freeboard (vertical distance between the surface of tailings and crest) which is 2 m. In general, total displacements of these magnitudes in a dam (see [58]) are considered to be within acceptable limits.

Figure 17: Vertical displacements in the dam after shaking with normal earthquake. The positive direction is upward.

[figure omitted; refer to PDF]

Figure 18: Vertical displacements in the dam after shaking with extreme earthquake. The positive direction is upward.

[figure omitted; refer to PDF]

4.4. Slope Stability

Slope stability of the dam, under static and seismic loading conditions, was assessed on the basis of a safety factor. In the finite element program PLAXIS, a strength reduction technique is employed to compute safety factors [18]. In this technique, the tangent of the friction angle and the cohesion of the soil are gradually reduced in the same proportion until the geotechnical structure fails [18]. In the strength reduction method, the safety factor is defined in the same way as described in the limit equilibrium methods [59]. It is to be noted that the strength reduction technique is applicable to the MC model only. Therefore, the MC model was utilized for computation of postseismic safety factors after performing dynamic analysis with the UBCSAND model.

A safety factor of about 1.45 was achieved directly after the 11th raising (static loading). Postseismic safety factors of about 1.32 and 1.22 were obtained from the slope stability analyses for the normal earthquake and the extreme earthquake, respectively. As expected, these postseismic safety factors are less than the safety factor obtained under static loading. This is because shear strength is reduced due to an increase in excess pore pressures caused by the earthquakes.

Postseismic stability of the dam was evaluated with the safety factor criterion suggested by ANCOLD [60] which states that a safety factor of at least 1.2 is required for a tailings dam to be safe enough against seismic shakings. The stability of the dam is considered to be satisfactory for both the earthquakes according to the guidelines provided by ANCOLD [60]. It is to be noted, for the extreme earthquake, the value of the postseismic safety factor of the dam is just on the limit recommended by ANCOLD [60]. This implies that the dam may not be safe enough for earthquakes stronger than the Virginia earthquake presented in this paper.

In this type of analysis, in addition to the value of the safety factor, it is important to identify possible failure zones that would have developed if the dam had failed. Figures 19 and 20 show such failure zones that were obtained from the postseismic slope stability analyses for the normal earthquake and the extreme earthquake, respectively. It can be seen that the failure zone caused by the extreme earthquake is less wide and deep compared to the one that developed due to the normal earthquake. This is because the extreme earthquake has led to larger negative pore pressures and higher strength in the moraine dikes and the rockfill berms (cf. Figure 3); therefore, the failure zone did not extend to these zones.

Figure 19: Failure zone of the dam after shaking with normal earthquake.

[figure omitted; refer to PDF]

Figure 20: Failure zone of the dam after shaking with extreme earthquake.

[figure omitted; refer to PDF]

5. Concluding Remarks

In this paper, seismic behaviour of the Aitik tailings dam was evaluated, based on results of dynamic finite element analyses, in terms of liquefaction potential, permanent deformations, and postseismic stability. The seismic behaviour of the dam was analysed for two cases: (i) a normal case (an earthquake of Swedish local magnitude 3.6) and (ii) an extreme case (an earthquake of moment magnitude 5.8). The results, for the extreme case, indicate that liquefaction could occur in a limited zone located below the surface near the embankments. Because of the limited extent of this liquefied zone, it is interpreted that the stability of the dam could not be affected. The predicted displacements induced by both the earthquakes are considered to be within acceptable limits. Postearthquake slope stability analyses suggest that the dam could resist the abovementioned earthquakes. However, for the Virginia earthquake, the postseismic safety factor is just on the limit suggested by the available guidelines. Since no vertical ground accelerations were used in the numerical modelling, the results may slightly differ if the vertical ground accelerations are applied in addition to the horizontal ground accelerations.

The stability of the Aitik tailings dam has been studied due to the occurrence of one earthquake at a time. The possibility of subsequent weakening of the structure by multiple earthquakes of various magnitudes has not been examined.

Due to unavailability of cyclic laboratory test results on the Aitik tailings material, dynamic material properties were evaluated with engineering judgement based on results of a case study relevant to this dam. If cyclic laboratory tests are conducted in the future, they could provide valuable information in further understanding of the seismic behaviour of the dam.

Acknowledgments

The authors would like to express their sincere thanks to Boliden AB, Sweden, for giving them an opportunity to carry out the presented study regarding the tailings dams at the Aitik mine and its disposal facilities and for giving access to site information. Kerstin Pousette at Luleå University of Technology, Sweden, and Fredrik Jonasson at Sweco Infrastructure AB, Sweden, are to be acknowledged for performing laboratory tests and providing information from field tests and for evaluation of material parameters. The first author would like to thank Lars Erik Lundbergs Foundation, Sweden, for granting him a scholarship to conduct some part of his Ph.D. research. The "IRIS" project within the EU Fp7 framework together with Luleå University of Technology, Boliden AB, and "Swedish Hydropower Centre-SVC" are acknowledged for financial support, which made the work possible. SVC has been established by the Swedish Energy Agency, Elforsk, and Svenska Kraftnät together with Luleå University of Technology, The Royal Institute of Technology, Chalmers University of Technology, and Uppsala University. Participating hydropower companies are Andritz Hydro Inepar Sweden, Andritz Waplans, E.ON Vattenkraft Sverige, Fortum Generation, Holmen Energi, Jämtkraft, Karlstads Energi, Linde Energi, Mälarenergi, Skellefteå Kraft, Sollefteåforsens, Statkraft Sverige, Statoil Lubricants, Sweco Infrastructure, Sweco Energuide, SveMin, Umeå Energi, Vattenfall Research and Development, Vattenfall Vattenkraft, VG Power, and WSP.

References

[1] SNSN, "Svenska nationella seismiska nätet," 2013 http://snsn.geofys.uu.se

[2] R. Bödvarsson, B. Lund, R. Roberts, R. Slunga, "Earthquake activity in Sweden: study in connection with a proposed nuclear waste repository in Forsmark or Oskarshamn," http://www.skb.se/upload/publications/pdf/R-06-67.pdf , SKB Rapport , no. R-06-67, 2006.

[3] GruvRIDAS, Gruvindustrins Riktlinjer för Dammsäkerhet , Svensk Energi AB/SveMin, Stockholm, Sweden, 2007.

[4] S. G. Vick, Planning, Design and Analysis of Tailings Dams , BiTech Publishers, Richmond, Canada, 1990.

[5] P. K. Robertson, C. E. Wride, "Evaluating cyclic liquefaction potential using the cone penetration test,", Canadian Geotechnical Journal , vol. 35, no. 3, pp. 442-459, 1998.

[6] WISE, "Chronology of major tailings dam failures," 2013, http://www.wise-uranium.org/mdaf.html

[7] M. Seid-Karbasi, H. Hawson, U. Atukorala, "Seismic stability of a Peruvian tailings earth-rockfill dam with liquefiable foundation,", Proceedings of International Symposium on Dams and Reservoirs under Changing Challenges--79 Annual Meeting of ICOLD, Swiss Committee on Dams, Lucerne, Switzerland, 1 June 2011 , pp. 613-620, CRC Press/Balkema, Leiden, The Netherlands, 2011.

[8] M. Seid-Karbasi, U. Atukorala, "Deformations of a zoned rockfill dam from a liquefiable thin foundation layer subjected to earthquake shaking," in Proceedings of the 31st Annual USSD Conference on 21st Century Dam Design-Advances and Adaptations, pp. 1351-1367, United States Society on Dams, San Diego, Calif, USA, April 2011.

[9] M. Seid-Karbasi, P. M. Byrne, "Embankment dams and earthquakes,", International Journal on Hydropower and Dams , vol. 11, no. 2, pp. 96-102, 2004.

[10] P. M. Byrne, M. Seid-Karbasi, "Seismic stability of impoundments," in Proceedings of the 17th Annual Symposium, Vancouver Geotechnical Society, Vancouver, Canada, 2003.

[11] G. Elia, M. Rouainia, "Seismic performance of earth embankment using simple and advanced numerical approaches,", Journal of Geotechnical and Geoenvironmental Engineering , vol. 139, no. 7, pp. 1115-1129, 2013.

[12] G. Elia, A. Amorosi, A. H. C. Chan, M. J. Kavvadas, "Fully coupled dynamic analysis of an Earth dam,", Geotechnique , vol. 61, no. 7, pp. 549-563, 2011.

[13] M. Alyami, M. Rouainia, S. M. Wilkinson, "Numerical analysis of deformation behaviour of quay walls under earthquake loading,", Soil Dynamics and Earthquake Engineering , vol. 29, no. 3, pp. 525-536, 2009.

[14] S. Sica, L. Pagano, A. Modaressi, "Influence of past loading history on the seismic response of earth dams,", Computers and Geotechnics , vol. 35, no. 1, pp. 61-85, 2008.

[15] K. K. Muraleetharan, S. Deshpande, K. Adalier, "Dynamic deformations in sand embankments: centrifuge modeling and blind, fully coupled analyses,", Canadian Geotechnical Journal , vol. 41, no. 1, pp. 48-69, 2004.

[16] O. Aydingun, K. Adalier, "Numerical analysis of seismically induced liquefaction in earth embankment foundations. Part I. Benchmark model,", Canadian Geotechnical Journal , vol. 40, no. 4, pp. 753-765, 2003.

[17] M. M. Dewoolkar, H.-Y. Ko, R. Y. S. Pak, A. Engrg, "Seismic behavior of cantilever retaining walls with liquefiable backfills,", Journal of Geotechnical and Geoenvironmental Engineering , vol. 127, no. 5, pp. 424-435, 2001.

[18] R. B. J. Brinkgreve, W. M. Swolfs, E. Engine, PLAXIS User's Manual , PLAXIS, Delft, The Netherlands, 2011.

[19] A. Petalas, V. Galavi, "Plaxis liquefaction model UBC3D-PLM: documentation report," January 2013, https://www.plaxis.com/content/uploads/import/kb/kb-downloads/UBC3D-PLM%20%28REPORT%29.Jan2013.pdf

[20] M. H. Beaty, P. M. Byrne http://www.itasca-udm.com/media/download/UBCSand/UBCSAND_ UDM_ Documentation.pdf , UBCSAND Constitutive Model Version 904aR: Documentation Report , 2011.

[21] P. England, J. Jackson, "Uncharted seismic risk,", Nature Geoscience , vol. 4, no. 6, pp. 348-349, 2011.

[22] E. Lindblom, Microearthquake study of end-glacial faults in northern Sweden [Licentiate thesis] , Uppsala University, Uppsala, Sweden, 2011.

[23] E. Lindblom, B. Lund, E. Lindblom, "Focal mechanisms and the state of stress along the Pärvie end-glacial fault, Northern Sweden,", Microearthquake Study of End-Glacial Faults in Northern Sweden. Licentiate thesis , Uppsala University, Uppsala, Sweden, 2011.

[24] E. Lindblom, B. Lund, A. Tryggvason, M. Uski, C. Juhlin, R. Bödvarsson, R. Roberts, E. Lindblom, "Microearthquake activity on the Pärvie end-glacial fault system, northern Sweden,", Microearthquake Study of End-Glacial Faults in Northern Sweden , of Licentiate thesis, Uppsala University, Uppsala, Sweden, 2011.

[25] R. Lagerbäck, M. Sundh, "Early Holocene faulting and paleoseismicity in northern Sweden," 2008, http://resource.sgu.se/produkter/c/c836-rapport.pdf, Research Paper , no. C836, Geological Survey of Sweden

[26] R. Arvidsson, "Fennoscandian earthquakes: whole crustal rupturing related to postglacial rebound,", Science , vol. 274, no. 5288, pp. 744-746, 1996.

[27] R. Muir Wood, S. Gregersen, P. W. Basham, "Extraordinary deglaciation reverse faulting in northern Fennoscandia,", Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound , pp. 141-173, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1989.

[28] B. Lund, P. Schmidt, C. Hieronymus, "Stress evolution and fault instability during the Weichselian glacial cycle," http://www.skb.se/upload/publications/pdf/TR-09-15.pdf , no. TR-09-15, Swedish Nuclear Fuel and Waste Management, 2009.

[29] W. Y. Kim, O. Kulhánek, T. van Eck, R. Wahlström, "The Solberg, Sweden, Earthquake of September 29, 1983,", Report No. , no. 1-85, Seismological Department, University of Uppsala, Uppsala, Sweden, 1985.

[30] FENCAT http://www.seismo.helsinki.fi/bulletin/list/catalog_format.html , Fennoscandian Earthquake Catalog , 2011.

[31] P. Mäntyniemi, "Earthquake of 4 November 1898 in northern Europe: new insights,", Journal of Geophysical Research: Solid Earth , vol. 113, no. 11, 2008.

[32] P. Mäntyniemi, S. D'Amico, "Intraplate seismicity and seismic hazard: the Gulf of Bothnia Area in Northern Europe Revisited," 2012 http://www.intechopen.com/books/howtoreference/earthquake-research-and-analysis-new-frontiers-in-seismology/intraplate-seismicity-and-seismic-hazard-the-gulf-of-bothnia-in-northern-europe-revisited, Earthquake Research and Analysis-New Frontiers in Seismology , pp. 283-298, InTech, Rijeka, Croatia

[33] H. Bungum, O. Olesen, "The 31st of August 1819 Lurøy earthquake revisited,", Norwegian Journal of Geology , vol. 85, no. 3, pp. 245-252, 2005.

[34] European Standard, "Geotechnical investigation and testing--identification and classification of soil--part 1: identification and description,", EN ISO , no. 14688-1, 2002.

[35] ASTM, "Standard practice for classification of soils for engineering purposes,", ASTM , no. D 2487-10, ASTM International, West Conshohocken, Pa, USA, 2010.

[36] Google Maps, "Satellite view of Aitik tailings dam near Gällivare, Sweden," 2011, http://maps.google.com

[37] L. Ormann, M. A. Zardari, H. Mattsson, A. Bjelkevik, S. Knutsson, "Numerical analysis of strengthening by rockfill embankments on an upstream tailings dam,", Canadian Geotechnical Journal , vol. 50, no. 4, pp. 391-399, 2013.

[38] L. Ormann, M. A. Zardari, H. Mattsson, A. Bjelkevik, S. Knutsson, "Numerical analysis of curved embankment of an upstream tailings dam,", Electronic Journal of Geotechnical Engineering , vol. 16, pp. 931-944, 2011.

[39] A. M. Britto, M. J. Gunn, Critical State Soil Mechanics via Finite Elements , Ellis Horwood, West Sussex, UK, 1987.

[40] D. M. Potts, L. Zdravkovic, Finite Element Analysis in Geotechnical Engineering. Theory , Thomas Telford Ltd, London, UK, 1999.

[41] R. L. Kuhlemeyer, J. Lysmer, "Finite element method accuracy for wave propagation problems,", Journal of Soil Mechanics and Foundations Division , vol. 99, pp. 421-427, 1973.

[42] T. Stål, P. Wedel, S. Avén 1984 (Swedish), Handboken Bygg Geoteknik , Liber, Stockholm, Sweden

[43] Vattenfall, Jord-och Stenfyllningsdammar , Vattenfall, Stockholm, Sweden, 1988.

[44] U. Bergdahl, E. Ottosson, B. S. Malmborg 1993 (Swedish), Plattgrundläggning , AB Svensk Byggtjänst och Statens Geotekniska Institut

[45] Bro. Vägverket 2004 (Swedish), Vägverkets Allmänna Tekniska Beskrivning för Nybyggande och Förbättring av Broar , Borlänge

[46] H. Puebla, P. M. Byrne, R. Phillips, "Analysis of CANLEX liquefaction embankments: prototype and centrifuge models,", Canadian Geotechnical Journal , vol. 34, no. 5, pp. 641-657, 1997.

[47] M. H. Beaty, P. M. Byrne, "An effective stress model for predicting liquefaction behaviour of sand,", Geotechnical Special Publication , vol. 75, pp. 766-777, 1998.

[48] P. K. Robertson, K. L. Cabal, Guide to Cone Penetration Testing for Geotechnical Engineering , Gregg Drilling & Testing, Inc. Signal Hill, Calif, USA, 2010.

[49] B. O. Hardin, "The nature of stress strain behavior of soils," in Proceedings of the ASCE Geotechnical Engineering Division Specialty Conference Earthquake Engineering and Soil Dynamics, pp. 3-90, American Society of Civil Engineers, Pasadena, Calif, USA, June 1978.

[50] E. Naesgaard, A hybrid effective stress--total stress procedure for analyzing soil embankments subjected to potential liquefaction and flow [Ph.D. thesis] , The University of British Coulombia, 2011.

[51] B. M. Das, G. V. Ramana, Principles of Soil Dynamics , Cengage Learning, Stamford, Conn, USA, 2011.

[52] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, Finite Element Method--Its Basis and Fundamentals , Elsevier Butterworth-Heinemann, 2005.

[53] ITASCA, Fast Lagrangian Analysis of Continua (FLAC) User's Guide , Itasca Consulting Group, Inc. Minneapolis, Minn, USA, 2011.

[54] M. H. Beaty, V. G. Perlea, "Several observations on advanced analyses with liquefiable materials," in Proceedings of the 31st Annual USSD Conference and 21st Conference on Century Dam Design-Advances and Adaptations, pp. 1369-1397, US Society on Dams, San Diego, Calif, USA, April 2011.

[55] T. D. Stark, M. H. Beaty, P. M. Byrne, G. Castro, F. C. Walberg, V. G. Perlea, P. J. Axtell, J. C. Dillon, W. B. Empson, D. L. Mathews, "Seismic deformation analysis of Tuttle Creek Dam,", Canadian Geotechnical Journal , vol. 49, no. 3, pp. 323-343, 2012.

[56] H. B. Seed, K. L. Lee, I. M. Idriss, F. Makdisi, "Analysis of the slides in the San Fernando dams during the earthquake of Feb. 9, 1971,", Report no. , no. EERC 73-2, University of California, Berkeley, Calif, USA, 1971.

[57] J. A. Knappett, R. F. Craig, Craig's Soil Mechanics , Spon Press, London, UK, 2012.

[58] R. N. Kostaschuk, J. M. T. Wilkinson, P. M. Byrne, "Predicting the deformation of tailings dams from earthquake liquefaction," in Proceedings of the 6th International Conference on Tailings and Mine Waste '99, pp. 237-248, AA Balkema, Fort Collins, Colo, USA, January 1999.

[59] R. B. J. Brinkgreve, H. L. Bakker, "Non-linear finite element analysis of safety factors," in Proceedings of the 7th International Conference on Computer Methods and Advances in Geomechanics, pp. 1117-1122, AA Balkema, Cairns, Australia, May 1991.

[60] ANCOLD, Guidelines on Tailings Dams: Planning, Design, Construction, Operation and Closure , 2012.