Stability Assessment of the Dam of a Tailings Pond Using Computer Modeling—Case Study: Coroiești, Romania

Abstract

Anthropogenic activities related to mining generate both progress and a vast amount of waste that is responsible for environmental degradation. The Jiu Valley is one of the areas of Romania where mining has affected large areas of land, used to build mines and tailings ponds. The former Coroiesti coal processing plant (CCPP) is such a location with a total area of 25 ha containing approximately 5.5 million tons of tailings. The assessment of the stability of tailings dams is extremely important from safety and environmental aspects. This study proposes a solution based on numerical methods for determining the stability of a section of the dam of a tailings pond. The model of tailings pond no. 1, compartment B, from the Coroieşti Coal Preparation was built using COMSOL Multiphysics. Two scenarios of stability analysis were conducted on a section of the tailings dam: the FOS was determined using the shear strength reduction (SSR) method for both the initial and the current state of this TP. This method is a modern alternative to the limit equilibrium method, and its implementation by COMSOL is new to our country, thus aligning this methodology with current worldwide trends and developments in the field. The results obtained proved to be in line with those calculated in the past with traditional analytical methods, proving that the safety criteria of the studied TP/TD are being met.

1. Introduction

The modern world could not function without the raw materials produced by the mining industry. As the demand for mineral products is increasing, mining produces vast amounts [1]—of the order of billion tons yearly—of various types of mining waste [2], including tailings. These are mixtures of crushed rock and processing fluids from mills, washeries or concentrators that remain after the extraction of the economic metals, minerals, mineral fuels or coal [3]. The tailings are usually stored under water in order to prevent the formation of surface dusts and acid. There are several approaches towards the handling and storage of tailings (i.e., submarine or riverine disposal, wetland retention, backfilling), but the main method currently used is their storage behind dammed impoundments, with the resulting structures termed “tailings ponds” or “tailings dams” (TP/TD).

The creation and management of these structures is a process that begins in the early stages of a mining operation but is continued throughout the mine’s life until its closure and beyond (post-closure), with the main objective being to safely store tailings and reduce the risks over the long term [4,5,6]. There are a large number of recorded failures of tailing ponds present in the literature [7,8,9,10], and most of them produced severe environmental damage and/or loss of human life.

Several scholars, notably McColl [11], Clarkson et al. [12,13] reviewed the causes and mechanisms leading to TP/TD failures, and categorized them as physical/structural, functional or environmental, each with several subcategories [14,15]. These failures are not usually produced by a single factor, but rather are the result of multiple factors acting together [16]. Most of the time, the influence of the external environment (earthquakes, rainfall, flooding, dam foundation subsidence, etc.,) cause changes in the stress and seepage fields, leading to the destabilization of the structure.

The most frequent causes of failure were determined to be linked with groundwater level and its influence on dam foundation-soil permeability [17,18], the erosion produced by heavy rainfalls [19,20,21,22], seismic activity [23,24,25]; all of these causes leading to failure of the dam’s slope due to weakened and liquefied underlying layers, overtopping of the liquefied material, and an increase in lateral pressure on the dams.

In order to prevent these situations, engineers developed several procedures and tools that can be used to assess the stability of TP/TD. Stability analysis methods play a crucial role in ensuring the safety and reliability of these structures [26]. By utilizing advanced techniques and software, engineers are able to assess factors such as the critical slip surface of slopes, plastic zone distribution, deformation characteristics of slopes, and factor of safety (FOS). Various methods are used for determining the FOS of TP/TD using different approaches and tools. First, there are conventional analytical approaches based on the limit equilibrium methods (LEM) [27], with multiple techniques belonging to LEM: Swedish slip circle, ordinary method of slices, i.e., Bishop, Morgenstern–Price, Spencer, Janbu, Fredlund, etc., [28,29,30,31,32] which usually establish very conservative and safe estimations of FOS due to their simplification assumptions. These methods have been used for several decades and proven to be effective in many practical engineering problems. Next, there are more modern approaches for assessing the safety of TP/TD including cloud model analysis [33], set pair analysis [34,35], or fuzzy logic [36,37].

Finally, there are the numerical methods of analysis [38,39,40]. As computer performance and processing power have increased in recent years, these have allowed the simulation of problems which otherwise were too complex for conventional methods, because of complex geometry, material anisotropy, nonlinear behavior, in situ stresses, etc. These types of methods can model both continuous and discontinuous slopes using the techniques of the finite element method (FEM) for continuum modelling [31,32,33,34,35,36,37,38,39,40,41,42,43], discrete elements method (DEM) for discontinuum modelling [44,45,46] and the hybrid/coupled method which combines FEM and DEM [47,48,49,50,51]. Stead et al. [52] reviewed the wide range of computer programs available to engineers for the analysis of slopes, emphasizing the emerging powerful numerical modelling codes that enable the realistic simulation of slope failure. However, it was found that each method and every software has its own advantages and limitations, so there is no accepted standard imposed, to date, requiring the choice of a particular software product or analysis method.

The COMSOL Multiphysics software and its geomechanics module [53] (part of the structural mechanics module) proved to be a reliable tool for the scope of the paper, given the type of simulation involved. This is a versatile multiphysics piece of software [54,55] that differs from niche programs in user interface, model building capabilities, solvers available, mesh generation complexity and also the possibility to solve coupled models (fluid flow and structural mechanics in this case). COMSOL provides multiple ways to define geometries using strong user or physics-controlled mesh generation tools; the soil plasticity can be simulated either by internal available models or by user-defined yield functions; it also offers either direct or iterative solvers and can generate results as charts, graphs, tables, animated outputs, or data that can be exported and handled in postprocessing tools. There are several simulations in geoengineering carried out with COMSOL. Zhou et al. [56] calculated several safety factors of a high southwest rock slope of the Dagushan iron mine under various rainfall conditions using the nonlinear strength reduction method. A study for the assessment of slope stability was developed by Wu et al. [57] where the COMSOL numerical model considered the coupled thermal (temperature variation), hydraulic (pore water pressure), and mechanical (stress and displacement) processes, obtaining results that coincided well with those given by the limit equilibrium method. Again, COMSOL Multiphysics software was used by Shao et al. [58] to compute a coupled dual-permeability model with a soil mechanics model for landslide stability evaluation on a hillslope and to compare the model with two existing solutions. Memon [59] conducted a study that focused on comparisons between safety factors calculated utilizing computer modeling (COMSOL) for simple slopes using the extensively used limit equilibrium methods and the less common finite element approach. The comparative analysis is based on a multitude of slope geometry and material property combinations. In their paper, Mukhlisin et al. [60] proposed a numerical approach using the same software of seepage analysis where soil moisture distribution, water movement phenomena and slope stability characteristics in an unsaturated soil slope were analyzed based on the strength reduction method using a two-step approach: pressure-head analysis, and static analysis. A multifield coupling theory was proposed by Wu et al. [61] to solve the coupled partial differential equations for soil slopes under rainfall infiltration using COMSOL. Sysala et al. [62] combined Plaxis, Matlab and COMSOL Multiphysics for an optimized approach to the modified shear strength reduction method (SSRM) used for slope stability analysis and determination of the factor of safety and related failure zones.

The purpose of this work is to assess the stability of a TP/TD for a tailings pond situated in Coroiesti, Romania. In order to do so, a stability analysis was conducted on a section of the tailings dam, with the determination of the FOS in various scenarios using the commercial software COMSOL Multiphysics.

2. Site Description

Several TP/TD resulted from the processing of the underground mining of coal in the Jiu Valley coal basin of Romania (Figure 1) and still exist today, whether in operation or under conservation. Due to their age, they all pose a risk for structural failure and thus of accidental environmental contamination or even human loss of life.

The coal processing plant from Coroiești was shut down several years ago, when three of the four mines that delivered the run-of-mine coal to it were closed. However, during the years of operation from the processing of coal after the flotation of the raw slurry, waste resulted in the form of tailings that was decanted and stored in two tailings ponds located on the right bank of the Western Jiul River, at a distance of approx. 1 km downstream from the processing plant.

Both these ponds were made by enclosing surfaces with perimetral dams made of mining and processing tailings. There are two tailing ponds situated east of the plant, named Tailings Pond no. 1 and Tailings Pond no. 2.

Tailings Pond 1 consists of two compartments (A) and (B) as shown in Figure 2, both of sedimentation type, which have functioned alternately over time, with an area of 140,000 m² and a volume of 3,000,000 m³. The elevation of the natural land at the base of TP no. 1 is between 579 and 589 m above sea level, and the top of the dam was at 600.50 m to 608.90 m above the sea level, when last surveyed in 2014. Since then, the dike of compartment B has been raised by about 2 m. Due to the exhaustion of the storage capacity of TP no. 1, it was decided to reactivate and store the tailings in TP no. 2.

Tailings Pond no. 2 is presented in Figure 3; it is located in continuation of the two compartments of TP no. 1 and was commissioned in 1968. It has a single compartment, with an area of 110,000 m² and a volume of 2,000,000 m³. The elevation of the natural land at the base varies between 576 m and 584 m above sea level, and the last measured elevation of the dam top was between 594.50 m and 599.30 m above the sea level.

The last assessment of the safety status of Tailings Ponds Nos. 1 and 2 from Coroieşti dates back to 2001 [63] and included the determination of the FOS, which was determined using the classic analytical methods of slices, with the calculations carried out with a simple computer program, for static loads as well as for pseudostatic loads. This resulted in several FOS, depending on the height and angle of the dam slopes, with values of FOS ranging from 1.48 at the dike height of 20 m and at the slope angle of 30 degrees to 1.33 at the slope angle of 35 degrees. The numerical modeling within this study based on finite elements method (FEM) was carried out for a two-dimensional section of the dam, with a total height of 30 m and a median slope of 30 degrees. The medium was considered as isotropic in terms of permeability.

3. Description of the Simulation Model

3.1. Theoretical Aspects of the Method Used

As stated in the Introduction section, slope stability analysis is important for the prediction of soil slippage, slides and deformation caused by various mechanical loads, and it is also imposed by safety concerns. For a section of the dam of the TP no. 1 Compartment B described in Section 2, such an analysis was conducted using the geomechanics module of COMSOL Multiphysics version 5 software. This software was chosen for the analysis, in consideration of its complexity, its extensive use in the literature, as presented in Section 1, and also based on the expertise of the authors’ team [64,65] in solving engineering problems using the same software.

Within the model studied, the software used Darcy’s law for the computation of the pore pressure in the soil and the fluid flow through the porous medium, while the Mohr–Coulomb criterion was used to represent the soil elastic–plastic behavior during the analysis. The Mohr–Coulomb criterion is a mathematical model which describes the response of a certain material to normal and shear stresses, and it is widely used in geoengineering.

This technique of determining the factor of safety is called the shear strength reduction (SSR) method. According to [66,67], the SSR technique for slope stability analysis involves systematic finite element analysis to determine a stress reduction factor (SRF) that brings the slope to the verge of failure, thus obtaining the FOS value. By definition [68], the factor of safety of a slope is the “ratio of actual soil shear strength to the minimum shear strength required to prevent failure,” or the factor by which soil shear strength must be reduced to bring a slope to the verge of failure.

The shear strengths of all materials in the model were progressively reduced by the SRF until collapse occurred during a multistep approach carried out automatically by the software, based on the algorithm: (1) Develop the model of the slope. Add the material properties established (strength and deformation). Compute and record the maximum deformation; (2) Increase the value of the SRF, add the reduced Mohr–Coulomb material parameters to the model. Re-compute and record the maximum deformation; (3) Repeat step 2 until the model does not converge. In other words, the strength parameters of the material are reduced until the slope fails. This critical SRF that is obtained just before the failure, is the actual FOS for the slope.

In this method, the Mohr–Coulomb material parameters are defined as functions of the SRF. The Mohr–Coulomb yield function defining the associated plastic potential can be expressed as [55,69]:

$$F=Q=m\cdot \sqrt{J_2}+\alpha \cdot I_1-k$$

(1)

where

$$\alpha =\frac{\sin \mathsf{\varphi }}{3}; k=C\cdot \cos \mathsf{\varphi }$$

(2)

With this technique, the FOS affects the cohesion as well as the angle of internal friction. Cohesion describes how strongly a material will stick together. The angle of internal friction describes the frictional shear resistance of the soil. The parametrized angle of internal friction $\mathsf{\varphi }$ and cohesion C are expressed [67,70] as:

$$\begin{array}{l}\mathsf{\varphi }=\arctan {\left(\frac{\tan {\mathsf{\varphi }}_u}{\mathrm{SRF}}\right)}_{(\mathrm{for}p<0)}+\arctan {\left(\frac{\tan {\mathsf{\varphi }}_s}{\mathrm{SRF}}\right)}_{(\mathrm{for}p\ge 0)}\\ C=\frac{c}{\mathrm{SRF}}\end{array}$$

(3)

where c is the cohesion of the material, ${\mathsf{\varphi }}_u$ and ${\mathsf{\varphi }}_s$ are the angles of internal friction for unsaturated and saturated soils, and p is the pore pressure as defined by Darcy’s law. The dam section is modelled in 2D in order to reduce the computation time, using the plane strain approximation method and including the hydrostatic pressure and gravity effects.

3.2. Geometry of the Model

Two scenarios were simulated. The first determined the FOS for the initial geometric characteristics of the TP, when it was commissioned in use, while the second determined the FOS for the current geometric characteristics of the TP, taking into consideration the slurry deposition on the bottom of the pond and thus the increased level of water. Figure 4 and Figure 5 present the geometry of the dams for the two scenario models.

4. Scenario 1: Determination of the Stability of the Tailings Dam in Its Initial State

4.1. Simulation Model Construction and Methodology

According to Darcy’s law, hydraulic conductivity has a random character being defined by means of a probability function with normal distribution, with two parameters and a standard deviation of 0.001 as in Figure 6.

The simulation parameters are defined in COMSOL Multiphysics using the model builder window under global definitions, as presented in Figure 7. They comprise the geometric dimensions of the model as well as the characteristics of the dam material.

Next, the geometry of the model was constructed using the geometry toolbar, as can be seen in Figure 8. In addition to the polygon representing the outline of the dam model, two points were defined. The first defines the water level with the coordinates ((Hw)/tan(alfa1*pi/180), Hw), while the second defines the possible seepage level (L1 + L2 + L3−Hs/tan(alfa2*pi/180), Hs).

In the model builder window, several variables were defined:

-

two of Boolean type: saturated and unsaturated, for the saturated and unsaturated dam materials;

-

K representing the hydraulic conductivity;

-

C for the parametrized cohesion;

-

PHI for the parametrized internal friction angle.

This step is highlighted in Figure 9.

The hydrostatic pressure was defined according to the y-coordinate of the coordinate system as:

$$H_{p0}=H_w-y$$

(4)

where

-

H_w is the height of the water layer (defined in the simulation parameters);

-

y is the current height relative to the origin of the coordinate system.

In the following steps, the water pressure is added as a boundary load that acts on the submerged face of the dam and the possible seepage, using the Darcy’s law interface of the COMSOL model builder (Figure 10).

Also, in the solid mechanics module, constraints were defined for the boundary surfaces as either (fixed) for the base of the geometry or (roller) for the faces that could move vertically. These operations are shown in Figure 11 and Figure 12.

The finite element mesh is of the user-controlled mesh type, for which the size of the finite elements depends on their position in relation to the adopted model. This was achieved by using the size expression feature with the spatial distribution function:

$$\mathrm{if} \left(Y>-L4/5 \& \& X>L1\ast 0.9 \& \& X<\left(L1+L2+L3\ast 1.3\right),1,10\right)$$

(5)

Figure 13 shows the resulting mesh.

4.2. Studies, Results of Simulation, and Discussion

Obtaining the results involved performing calculations by running separate studies:

(1)

A first study using the Darcy’s law interface;

(2)

A second study using the solid mechanics interface with the in situ stress initialization;

(3)

A third study using the solid mechanics interface taking into account the FOS.

The computing of the first study enabled the pressure head (hydrostatic pressure) graphic representation (Figure 14) as it acts on the dam.

A positive pressure head indicates a positive pore pressure and thus a saturated soil, while a zero pressure head means unsaturated soil. The phreatic surface that divides the unsaturated and saturate soil is indicated by the zero-pressure head line. The hydrostatic pressure can be expressed as:

$$\psi =\frac{p}{\delta \cdot g}$$

(6)

where p is the pressure [N/m²]; $\rho$ is the density [kg/m³]; g is the gravitational acceleration [m/s²].

Figure 15 presents the slip area just before the slope failure, with the arrows highlighting the direction of displacement. For a more suggestive image of the soil slipping phenomenon, a combined contour and arrow representation diagram was used.

It must be noted that Figure 15 presents the slipping for a FOS value of 1.45. This value was imposed in the third aforementioned study where the FOS was defined in the [1–1.45] interval with an increment step of 0.001.

For higher values, the simulation of the adopted model did not converge, which means that the 1.45 value of FOS is the upper limit where the slope fails due to increased plastic stresses and a subsequent reduction in shear strength. The arrows show the direction of the soil particles’ displacement. In the bottom-right corner, the soil did not slip because of the base boundary fixed constraints.

The pattern in Figure 16 shows the actual plastic strain just before the collapse, indicating the failure mechanism. The displacement is also shown in a tridimensional representation in Figure 17, obtained by using the extrusion datasets in the results menu.

This result means that the soil slippage prediction using the plane strain approximation is a precise and effective way, with the tridimensional visualization being available after postprocessing thus avoiding the processing-intensive task of solving complex 3D problems.

Figure 18 presents the displacement amplitude in a 2D view. Finally, the graphic of the maximum displacement function of the FOS is plotted in Figure 19, which indicated after further analysis that the displacement significantly increases when the FOS exceeds 1.4, indicating the start of the slope failure.

As mentioned, three separate stationary studies were created and run in order to obtain relevant results for the pore pressure effect on the slope stability. In the first study, the hydrostatic pressure (pressure head) shape was obtained by computing Darcy’s Law, considering the pore pressure. The obtained hydrostatic pressure combined with gravity was then added to simulate the slope stability in the second study. Finally, in the third study, the stress and strain features were added, along with the pore pressure from the first study and the initial stresses generated by gravity from the second study.

By including the Mohr–Coulomb criterion, we could examine the elastic–plastic failure of the soil as a result of the gravity and pore pressure variation combined, and thus the FOS was determined for the slope.

The internal friction angle for unsaturated and saturated soils is different, resulting in the use of pore pressure instead. Thus, in regions with unsaturated soil, due to the fact that pores are considered interconnected and subject to constant atmospheric pressure, no external stresses were applied.

5. Scenario 2: Determination of the Stability of the Tailings Dam in Current State, Taking into Consideration the Slurry Deposition on the Bottom of the Pond

This scenario was simulated using the model constructions and methodology similar to those described in Section 4.1.

The difference was in the simulation parameters definition in COMSOL Multiphysics, when a new parameter is introduced (Figure 20), along with the other geometric dimensions of the model and characteristics of the dam material; the height of the slurry layer denoted by H_sl.

Studies, Results of Simulation, and Discussion

The studies were constructed, defined and run similarly to the first scenario. By computing, the simulation allowed for the representation of the same data: the pressure head (hydrostatic pressure) graphic representation (Figure 21) as it acts on the dam; the slip area just before the slope failure with the arrows highlighting the direction of displacement (Figure 22); the actual plastic strain just before the collapse, indicating the failure mechanism (Figure 23); the displacement of the dam in a tridimensional representation (Figure 24); the displacement amplitude in a 2D view (Figure 25); and finally, the plot of the maximum displacement function of the FOS (Figure 26), which indicated in the case of this scenario that the displacement significantly increases when the FOS exceeds 1.15.

6. General Conclusions

By comparing the results obtained in the two scenarios that were simulated, a series of aspects were observed.

The FOS is higher in the initial state of the TP, leading to the conclusion that although the height of the water is higher in its initial state, the deposition of slurry on the bottom of the pond determines over time a decrease in the stability factor. The comparison of the pressure heads (hydrostatic pressure) diagram, shows a pressure distribution in almost the entire dam in the current state, as compared to the initial state, where the pressure is concentrated towards the bottom of the dam. In both situations analyzed, there is the same tendency of the outer slope of the dam to fail. The maximum plastic strain just before the collapse is four times higher in the initial state as compared to the current state of the tailings dam. Also, the displacement of the dam is five times greater in the original state as compared to the current state. In both cases analyzed, the plots of the variation in the maximum displacement function of the FOS had a rapid increase after FOS values of 1.40 for the initial state and 1.15 for the current state.

Comparatively, the results obtained by simulation for the FOS are consistent with those obtained in the last known study conducted [63] in which the safety factor varied in the range of 1.33–1.48. For the initial state of the TP/TD, the difference between the value of the FOS of 1.48 obtained using the classic analytical method and 1.40 obtained with COMSOL was 5.71%. For the current state of the TP/TD with slurry deposition on the bottom of the pond, the FOS value of 1.15 is still within the safe limits.