Numerical Analysis of an Upstream Tailings Dam Subjected to Pond Filling Rates

Abstract

One of the challenges in upstream tailings dam projects is to ensure the allowable rate of deposition of tailings in the pond (i.e., pond filling rate) while maintaining the stability of the dam. This is due to the fact that an upstream tailings dam is constructed by placing dikes on top of previously deposited soft tailings, which could lead to a decrease in dam stability because of the build-up of excess pore water pressure. The main purpose of this work is to investigate the effects of pond filling rates on excess pore water pressure and the stability of an upstream tailings dam by a numerical study. A finite element software was used to simulate the time-dependent pond filling process and staged dam construction under various pond filling rates. As a result, excess pore water pressure increased in each raising phase and decreased in the subsequent consolidation phase. However, some of the excess pore water pressure remained after every consolidation phase (i.e., the build-up of excess pore water pressure), which could lead to a potentially critical situation in the stability of the dam. In addition, the remaining excess pore water pressure varied depending on the pond filling rates, being larger for high filling rates and smaller for low filling rates. It is believed that the approach used in this study could be a guide for dam owners to keep a sufficiently high pond filling rate but still ensure the desirable stability of an upstream tailings dam.

1. Introduction

This study focuses on the effects of pond filling rates on excess pore water pressure and the stability of an upstream tailings dam. Owing to a low initial investment (as a small amount of building material is required), the upstream method is the most popular design for a raised tailings embankment, especially in low-risk seismic areas [1]. It should be noted that there are over 3500 tailings dams worldwide, of which 50% are of the upstream design [2]. One of challenges in an upstream tailings dam is to ensure the allowable rates of deposition in the pond (i.e., pond filling rates) and still maintain the stability of a tailings dam [3]. This is due to the fact that an upstream tailings dam is constructed by placing dikes on top of previously deposited soft tailings, which could lead to a decrease in dam stability [4] because of the build-up of excess pore water pressure. During the operation of an upstream tailings dam, a small starter dike is first built at the extreme downstream toe. The dam wall is then progressively raised on the upstream side, mainly founded on the tailings beach. The pond filling operation process is then gradually accomplished as the height of the dam increases during the stage of construction [5,6]. In general, an adequate pond filling rate allows excess pore water pressure to dissipate gradually during the consolidation process. However, if the consolidation process is not completed when a new dike is built, i.e., in case of high pond filling rates, the excess pore water pressure will not have enough time to dissipate. In this sense, excess pore water pressure could build up under the dikes, which could lead to a critical stability situation of the dam [7]. A close relationship can thereby be seen between the consolidation process and the stability of the dam. In the present study, an approach which can simulate both consolidation processes and stability for a tailings dam is needed. Somogyi [8] investigated a slow deposition process in impoundments of slimes and tailings based on a nonlinear one-dimensional numerical approach. Gassner and Fourie [9] used a simple one-dimensional numerical approach for optimizing the allowable rate of deposition on tailings dams.

However, due to the complexities of the consolidation process coupled with the stability analysis as well as the complicated behavior of tailings materials, an advanced numerical approach would be appropriate. Note that the actual consolidation in the facility is two-dimensional under the assumption of plane strain or axisymmetric conditions, at least at the straight portions of the tailings impoundments and dams [10]. With regard to the use of the advanced numerical approach for a tailings dam, there have been a number of studies addressing two-dimensional numerical simulations [7,11,12,13,14,15] and three-dimensional numerical simulations [16,17,18]. In the study by Psarropoulos and Tsompanakis [12], an investigation on the mechanical behavior and the stability of a tailings dam under static and dynamic loading was performed. Two-dimensional numerical simulations were used in their study for three typical types of tailings dams. Zandarín et al. [13] developed a numerical model to study the stability of tailings dams subjected to the role of capillarity. Their numerical model considered the consolidation processes under self-weight loads, continuous addition of water with the tailings discharge, infiltration of rainwater, and evaporation. Ormann et al. [7,11] and Zardari et al. [15] boused the two-dimensional finite element method for the simulations of both consolidation processes and the stability of an upstream tailings dam. Ormann et al. [7,11] focused on static aspects, i.e., stability of a curved embankment and strengthening by rockfill embankments, whereas Zardari et al. [15] focused on earthquake-induced liquefaction of an upstream tailings dam. However, none of the abovementioned references have closely evaluated the effects of the pond filling rates on excess pore water pressure and the stability of an upstream tailings dam. Vick [3] addressed the need to manage the pond filling rates of an upstream tailings dam to prevent the build-up of excess pore water pressure that can reduce the shear strength of the fill material. Excessive rates of the pond filling would cause a trigger for static liquefaction that has been the underlying cause for many upstream tailings impoundment failures [19].

In this study, a numerical approach (two-dimensional finite element method) was used to investigate the effects of pond filling rates on excess pore water pressure and the stability of an upstream tailings dam. The time-dependent pond filling process and the stage of dam construction under various pond filling rates were considered in the simulations. The UBCSAND model was used for the tailings materials, whereas the Mohr–Coulomb (M-C) model was applied for moraine and rockfill materials. The results show that the finite element method can be a useful tool for studying how an upstream tailings dam should be built in order to be stable enough for different pond filling rates. With the finite element method, a finite model of a tailings dam can be built easily, e.g., a gradual raising of a dam in its geometry.

2. Numerical Modelling and Methods

2.1. General Description

An idealized upstream tailings dam was simulated using the finite element software Plaxis 2D [20]. Note that the plane strain condition is commonly adopted in simulations of tailings dams as it is a long straight section (except for corners) [7,12,21]. A full model was 360 m wide and 55 m high after the last construction stage, as shown in Figure 1a. In general, the model size should be large enough to avoid boundary effects on simulation outputs. In this study, the left vertical boundaries were 205 m away from the dam in all models. The slope inclination of 2H:1V (horizontal to vertical) was adopted as in the previous work of a simplified upstream tailings dam simulation by Psarropoulos and Tsompanakis [12]. The model comprised six soil types: moraine (foundation), rockfill (downstream support), initial dike (starter dike), layered tailings, compacted sand tailings (dikes), and rockfill (Figure 1b). Figure 1b shows the initial pond filling step, which was located at 10 m above the initial ground level. From this level, five construction stages comprising a raising phase (R) and a consolidation phase (C) were simulated in the model. In each raising phase, a new dike was constructed on top of previously layered tailings, assumed to be completed in 10 days and followed by a 10-day raising period for both new layered tailings and rockfill. The consolidation analysis was carried out to investigate the transient system response allowing excess pore water pressure to dissipate gradually during the consolidation process. The time in the consolidation phase (consolidation time) depends on the pond filling rate. In this current work, four case studies corresponding to four pond filling rates of 10.0 m/year, 5.0 m/year, 3.3 m/year, and 2.5 m/year were used to investigate the effects of the pond filling rates on excess pore water pressure and the stability of the upstream tailings dam. Therefore, the consolidation times in each stage were set over 163 days, 345 days, 528 days, and 710 days corresponding to the pond filling rates of 10.0 m/year, 5.0 m/year, 3.3 m/year, and 2.5 m/year, respectively. The phreatic line was assumed to be located at the surface of the tailings impoundment (worst-case scenario). The locations of the phreatic lines were adapted depending on the pond filling steps. In each analysis step, a new set of hydraulic boundary conditions was imposed, i.e., closed boundaries at the bottom (Y_min) and the left vertical edge (X_min), and open boundaries for others. As for the mechanical boundary conditions, the model was assumed to be fully fixed at its bottom. The horizontal displacements were assumed to be zero along the lateral edges (i.e., both left and right vertical boundaries). These hydraulic and mechanical boundary conditions can be found in many previous works dealing with two-dimensional numerical simulations of upstream tailings dams [7,11,12,13,14,15]. The finite element mesh of the upstream tailings dam after the last construction stage is shown in Figure 2. The finite element mesh in each cluster was composed of 15-node triangular elements. These elements give a fourth-order (quartic) interpolation for displacements [20]. A massive number of elements were generated in the areas of interest (structural zone), providing the finer mesh near the slope with embankments and rockfill. This is due to the fact that these areas would be affected by large strains during the stage of construction. The coarser mesh was then generated at the far-field areas to minimize computation time.

Apart from the fully coupled analysis on deformation and consolidation, the global factor of safety (FoS) for slope failure was computed for every stage of construction. The global factor of safety is computed in the finite element software Plaxis 2D using the shear strength reduction method. The idea of this method is that the soil strength is gradually reduced, and when a failure occurs, the corresponding strength reduction factor can be considered as the factor of safety of soil strength [20,22,23,24,25]. In this study, the factor of safety should have a value of at least 1.5 under normal operation conditions according to the Swedish dam safety guidelines [26]. In order to maintain the factor of safety of 1.5 at every stage of the construction, a strengthening by rockfill berms was added on the downstream side of the slope. The procedure for the construction of the rockfill berms was simulated in the finite element model. The rockfill berms were placed immediately after each dike construction in the finite element model. The width of the rockfill was then optimized in this study by a parametric sensitivity analysis where the rockfill width was varied until no FoS became smaller than approximately 1.5 at every stage of construction. Four rockfill widths of 5.0 m, 10.0 m, 15.0 m, and 20.0 m were applied to reach the target (discussed in Section 3.1, Figure 6).

2.2. Constitutive Models

In this study, the UBCSAND constitutive model was used to simulate the tailings materials, whereas the Mohr–Coulomb (M-C) model was applied for moraine and rockfill materials. The M-C model is a linear elastic–perfectly plastic model. The UBCSAND model is an effective stress plasticity model used in advanced stress deformation analyses of geotechnical materials. A fully coupled analysis (mechanical and groundwater flow calculations) can be performed simultaneously in this model. It can predict the shear stress–strain behavior of soil using an assumed hyperbolic relationship [27]. Therefore, this model can overcome the drawback of the elastic–perfectly plastic M-C model with a possibility to capture the nonlinearity in the elastic part of the soil. In addition, the UBCSAND model is appropriate for tailings materials as excess pore water pressure can build up in tailings during pond filling and dam raising. Figure 3 presents the behavior of tailings materials during drained simple shear tests and simulations. The laboratory tests were performed by Wiklund [28] at the Luleå University of Technology. Several simulation models were performed in this study to evaluate parameters for the UBCSAND model. As can be observed, while the stress–strain curves followed the nonlinear behavior of the experimental data well in the UBCSAND model, they overestimated the shear strength of the tailings materials significantly in the M-C model, regardless of the applied normal stresses (50 kPa, 100 kPa, 150 kPa) or the tailings types (layered tailings, compacted tailings). The UBCSAND model has also been used for materials in the development of excess pore water pressure-induced static liquefaction [29] or seismic liquefaction [15,30].

The M-C input parameter values of moraine and rockfill materials were adopted from previous works [31,32], as summarized in

Table 1. Parameters of Mohr–Coulomb (M-C) model.

Parameter	Symbols	Moraine (Foundation)	Moraine (Initial Dike)	Rockfill	Unit
Soil unit weight above phreatic level	γunsat	20	20	18	kN/m3
Soil unit weight below phreatic level	γsat	22	22	20	kN/m3
Young’s modulus	E	20,000	20,000	40,000	kN/m2
Cohesion	c’	1	1	1	kN/m2
Friction angle	φ’	37	35	42	degree
Horizontal hydraulic conductivity	kx	4.98 × 10−8	1.0 × 10−7	1 × 10−1	m/s
Vertical hydraulic conductivity	ky	1.0 × 10−8	4.98 × 10−8	1 × 10−1	m/s

, except for the Poisson’s ratio and the dilatancy angle. In this study, the Poisson´s ratio (υ) and the dilatancy angle (ψ) were assumed to be 0.33 and 0, respectively. These assumed values have been used in many simulation works of tailings dams in Sweden by Ormann et al. [7,11], Knutsson et al. [21], and Zardari et al. [15]. The laboratory tests were performed on tailings materials taken from an iron mine in northern Sweden. In addition, the parametric sensitivity analyses (curve fitting method) for obtaining elastic shear modulus number (K_G^e), plastic shear modulus number (K_G^p), and elastic bulk modulus number (K_B^e) were performed. In the curve fitting method, the friction angles and cohesion, which were evaluated based on the direct shear test results, were fixed. Elastic shear modulus number (K_G^e) and plastic shear modulus number (K_G^p) were varied to capture the stress–strain curves from the experiments. Finally, elastic bulk modulus number (K_B^e) was estimated as 70% of elastic shear modulus number (K_G^e), as recommended in [20] (material models part). As can be seen in Figure 3, there is a good agreement between the simulations (UBCSAND model) and experimental results, regardless of the tailings types used. The index parameters m_e, n_e, and n_p were assigned the values of 0.5, 0.5, and 0.4, respectively, as suggested in the original UBCSAND report by Beaty and Byrne [27]. The hydraulic conductivity of the tailings was adopted from previous studies [28,31]. In the present study, the hydraulic conductivity values in the horizontal direction were assumed to be 10 times higher than those in the vertical direction due to the layered nature of the tailings [26,31]. All input parameter values of the UBCSAND model used in the numerical analyses are tabulated in

Table 2. Parameters of UBCSAND model.

Parameter	Symbols	Layered Tailings	Compacted Tailings (Dikes)	Unit
Soil unit weight above phreatic level	γunsat	17	18	kN/m3
Soil unit weight below phreatic level	γsat	20	20	kN/m3
Elastic shear modulus number	KGe	130	130	-
Plastic shear modulus number	KGp	50	75	-
Elastic bulk modulus number	KBe	91	91	-
Elastic bulk modulus index	me	0.5	0.5	-
Elastic shear modulus index	ne	0.5	0.5	-
Plastic shear modulus index	np	0.4	0.4	-
Constant volume friction angle	ɸv	25	28	degree
Peak friction angle	ɸp	28	36	degree
Cohesion	c’	0	0	kN/m2
Horizontal hydraulic conductivity	kx	1.0 × 10−7	1.0 × 10−6	m/s
Vertical hydraulic conductivity	ky	1 × 10−8	1.0 × 10−7	m/s

.

3. Simulation Results

3.1. Stability Analyses

Stability analysis is an important component in the design of any earth structure, including tailings dams. Stability is usually expressed in terms of the factor of safety (FoS), which is defined as the ratio between the available shear strength of the soil and the minimum shear strength required against failure [12]. Figure 4 illustrates the factor of safety of the dam (pond filling rate of 5.0 m/year) during the pond filling, comprising raising phases (R) and consolidation phases (C). As shown, the FoS decreased gradually and dropped to less than 1.5 after the construction of three dikes of a height of 20 m from the initial pond filling level, i.e., FoS was 1.37 for R3 and 1.48 for C3. These values kept falling and even dropped to around 1.0 after the last construction stage. In addition, it was found that the FoS values decreased after every raising phase and then increased again due to the subsequent consolidation phase. This can be explained by the increase in effective stresses at the same rate as the excess pore pressure dissipates. The shear strength increases gradually during the consolidation phase. In this sense, the lowest stability of the dam is expected immediately after a new dike construction (a raising phase). This observation is similar to that of many previous studies on the stability of an upstream tailings dam [7,12,21,33]. The possible failure mechanism of the dam is shown in Figure 5. The slip surface most likely associated with the lowest FoS occurred according to Figure 5. In this case, the FoS did not meet the adopted code requirements of the Swedish dam safety guidelines. Therefore, a slope strengthening was introduced in this study to maintain the factor of safety of 1.5 at every stage of construction.

Strengthening an upstream tailings dam with rockfill berms at the downstream slope is a common technique that has been introduced and used in many previous works [7,15,21,25,33]. In this study, the width of the rockfill was optimized. The rockfill width was increased in the models (with the rockfill widths of 5 m, 10 m, 15 m, and 20 m) until no FoS became smaller than approximately 1.5 at every stage. The factors of safety for all the cases can be observed in Figure 6. These observed factors of safety were higher than those from the upstream tailings dam without rockfill berms. The weight of the rockfill berms placed next to the dikes could provide a resisting moment which increases the slope stability. As expected in all the cases, the FoS values decreased after every raising phase and then increased again due to the subsequent consolidation phase, regardless of rockfill widths. This trend is found to be similar to that in Figure 4. In conclusion, the rockfill width of 20 m should be used to keep the dam stable enough (FoS ≥ 1.5) for the pond filling rate of 5.0 m/year.

3.2. Excess Pore Water Pressure Analyses

Figure 7 shows the excess pore water pressure at the last construction stage (filling rate of 2.5 m/year and rockfill width of 20 m). There was a relatively high excess pore pressure in some zones beneath the dikes after the raising phase R5 (Figure 7a). However, this pressure decreased significantly after the consolidation phase C5 (Figure 7b). The remaining high excess pore pressure after the consolidation phase C5 occurred in the lower part of the impoundment but further away from the dam. This is due to the presence of the impermeable base in the finite element model. The trend of the excess pore water pressure distribution observed in this study is similar to that reported by Ormann et al. [7] and Saad and Mitri [14], i.e., larger remaining excess pore pressure after a consolidation phase occurred in the lower part of the impoundment but farther away from the dam.

3.3. Effects of the Pond Filling Rates on the Stability and Excess Pore Water Pressure

Figure 8 shows the effect of the pond filling rates (PFR) on the stability of the upstream tailings dam with the rockfill width of 20 m during the pond filling, comprising raising phases (R) and consolidation phases (C). As shown, the FoS decreased very slightly when the PFR increased from 2.5 m/year to 5.0 m/year, regardless of construction phases. At the pond filling rate of 5.0 m/year, the factor of safety of the dam at every stage of construction fulfilled the adopted code requirement of the Swedish dam safety guidelines, except for the raising phase R5 (FoS at R5 was 1.472). However, when the PFR increased from 5.0 m/year to 10 m/year, the FoS suddenly dropped to 1.40 at the raising phase R5. Thus, the FoS did not satisfy the abovementioned requirement, which might be a sign of potential instability issues for the dam.

Figure 9 shows excess pore water pressure after the last consolidation phase C5 with respect to the pond filling rates. For the ease of interpretation, the shading views of excess pore water pressure in Figure 9 are fixed at the same color intervals, except for the model with the pond filling rate of 10 m/year. As shown, an effect of the pond filling rate on excess pore water pressure was discovered: an increase in pond filling rate can lead to an increase in excess pore water pressure in the model. In particular, the maximum excess pore water pressure in the model with the pond filling rate of 2.5 m/year was 29 kPa. This value rose up to 43 kPa, 68 kPa, and 138 kPa in the models with the pond filling rates of 3.3 m/year, 5.0 m/year, and 10 m/year, respectively. In cases of high pond filling rates, there was less time for the excess pore pressure to dissipate during consolidation.

In addition, excess pore water pressure which developed along cross sections A–A’ and B–B’ (defined in Figure 1b) are plotted in Figure 10. Excess pore water pressure increased with an increase in the PFR, regardless of the cross sections analyzed. Take the elevation of 5 m from the initial ground level as an example: at the section A–A’, excess pore water pressure in the model with the pond filling rate of 2.5 m/year had the lowest value of 4.65 kPa, followed by 7.50 kPa (PFR = 3.3 m/year), 12.25 kPa (PFR = 5.0 m/year), and 27.25 kPa (PFR = 10 m/year). Similarly, at the cross section B–B’ at 5 m height, excess pore water pressure in the models with the pond filling rates of 2.5 m/year, 3.3 m/year, 5.0 m/year, and 10 m/year were 7.5 kPa, 12.1 kPa, 19.6 kPa, and 43.1 kPa, respectively. Interestingly, the same trend was observed at both cross sections, even though excess pore water pressure at the cross section B–B’ was higher than the corresponding excess pore water pressure at the cross section A–A’. It should be noted that cross section A–A’ had a draining distance shorter than cross section B–B’. It was also found from Figure 10 that excess pore water pressure decreased with increased heights from the initial ground level, regardless of the PFR and sections. This is because of the presence of the impermeable base in the finite element model.

Finally, the build-up of excess pore pressure (EPP) was also considered in this study. Figure 11 shows the basic concept of excess pore water pressure build-up during the stage of construction at point B and PFR = 10 m/year. As shown, the excess pore water pressure increased after every raising phase and then decreased after every consolidation phase which followed. Again, this finding supports the abovementioned stability analyses, that is, the FoS values decrease after raising phases and then increase again after consolidation phases. More importantly, it is worth observing that EPP remained after every pond filling step (i.e., after consolidation phases C1, C2, C3, C4, and C5). The remaining EPP were built up until the last phase of consolidation C5. Based on this concept, the EPP build-up during the stage of construction with respect to pond filling rates (PFR) is presented in Figure 12. Typical points A and B (in Figure 1b) were selected for the analyses. As shown, the build-up of excess pore water pressure was dependent on the PFR. Take point B in Figure 12b as an example: when the PFR was small (i.e., PFR = 2.5 m/year), the remaining EPP were EPP_C1 = 3.9 kPa, EPP_C2 = 4.0 kPa, EPP_C3 = 5.1 kPa, EPP_C4 = 6.3 kPa, and EPP_C5 = 7.6 kPa after the consolidation phases C1, C2, C3, C4, and C5, respectively. The remaining EPP, in this case, was built up, reaching only a small value of 7.6 kPa (EPP_C5) after the last consolidation phase C5. However, higher values of EPP_C5 were found when the PFR was higher (i.e., 12.0 kPa, 19.7 kPa, and 43.2 kPa for the cases with PFR = 3.3 m/year, PFR = 5.0 m/year, and PFR = 10.0 m/year). The same effect of the pond filling rates on the EPP build-up during the stage of construction can be observed in the results of point A (Figure 12a).

4. Discussion

It is concluded by Ormann et al. [7] that the stability of an upstream tailings dam could be improved by utilizing rockfill berms as supports on the downstream side. The stability analyses obtained in this study agree with this conclusion; refer, for example, to Figure 4 and Figure 6. In addition, it was found from this study that, in order to keep the stability of an upstream tailings dam, the exceptionally high pond filling rate is not recommended. However, the recommendation is only applicable in the case of 20 m rockfill width and the tailings type used in this study. In practice, larger rockfill width could be added next to dikes or on the downstream side to enhance the slope stability and hold the production speed. In this case, key questions remain to be answered: (i) which design of rockfill berms is needed for different production speeds to maintain the stability of a tailings dam, and (ii) how can the volume of rockfill berms be minimized (i.e., for minimum cost) while still maintaining the stability and production speed?

In practice, an upstream tailings dam is typically constructed by placing dikes on top of previously deposited soft saturated tailings, which could lead to the excess pore water pressure gradually increasing during pond filling and dam raising. Therefore, an adequate pond filling rate is required to allow excess pore water pressure to dissipate gradually during the consolidation process. The numerical results from this study showed that excess pore water pressure increased in each raising phase and decreased in the subsequent consolidation phase. However, some of the excess pore water pressure remained after every consolidation phase (i.e., the build-up of excess pore water pressure) (Figure 11). The build-up excess pore water pressure varied depending on the pond filling rates. The dam with a high pond filling rate had a shorter time for the excess pore pressure to dissipate during every consolidation phase than the dam with low ones. Therefore, the highest value was found in the dam with the pond filling rate of 10 m/year in this study (Figure 12).

Even though the numerical results obtained are encouraging, some additional work will need to be performed in the future to simulate other conditions involved in the responses of an upstream tailings dam under cyclic and seismic loading.

5. Conclusions

In this study, a numerical analysis of an upstream tailings dam subjected to pond filling rates was conducted. A finite element software was used to simulate the time-dependent pond filling process and the staged dam construction under various pond filling rates. An advanced constitutive model was adopted in the finite element models to capture the nonlinearity in the elastic deformations of the tailings materials. Based on the results of this study, the following conclusions can be drawn:

• Excess pore water pressure increased after every raising phase and then decreased after every subsequent consolidation phase in the numerical models. However, excess pore water pressure remained after every consolidation phase, leading to the build-up of excess pore water pressure and, eventually, a potentially critical situation in the stability of the dam (i.e., especially for the dam without strengthening).
• The stability of the dam decreased during the raising phase but increased during the consolidation phase. This is due to the changes of excess pore water pressure, i.e., increased values with the raising phase and decreased values with the consolidation phase. In addition, it is also vital to mention that the strengthening of an upstream tailings dam with rockfill berms at the downstream slope could improve the stability of the dam.
  The basic concept of excess pore water pressure build-up during the stage of construction was presented based on the remaining excess pore water pressure after every consolidation phase. The remaining excess pore water pressure varied depending on the pond filling rates, being larger in high filling rates and smaller in low filling rates. A very high pond filling rate is, therefore, not recommended for an upstream tailings dam as the dissipation of excess pore water pressure usually takes a long time.
• This study provided an in-depth analysis of the build-up of excess pore water pressure and slope stability during the stage of construction (pond filling and dam raising) of an upstream tailings dam. It is believed that the finite element method could be used to study how an upstream tailings dam should be built in order to be stable for different pond filling rates.