Three-dimensional numerical simulation of mud flow from a tailing dam failure across complex terrain

A tailing dam accident can cause serious ecological disaster and property loss. Simulation of a tailing dam accident in advance is useful for understanding the tailing flow characteristics and assessing the possible extension of the impact area. In this paper, a three-dimensional (3-D) computational fluid dynamics (CFD) approach was proposed for reasonably and quickly predicting the flow routing and impact area of mud flow from a dam failure across 3-D terrain. The Navier–Stokes equations and the Bingham–Papanastasiou rheology model were employed as the governing equations and the constitutive model, respectively, and solved numerically in the finite volume method (FVM) scheme. The volume-of-fluid (VOF) method was used to track the interface between the tailings and air. The accuracy of the CFD model and the chosen numerical algorithm were validated using an analytical solution of the channel flow problem and a laboratory experiment on the dam-break problem reported in the literature. In each issue, the obtained results were very close to the analytical solutions or experimental values. The proposed approach was then applied to simulate two scenarios of tailing dam failures, one of which was the Feijão tailing dam that failed on 25 January 2019, and the simulated routing coincided well with the in situ investigation. Therefore, the proposed approach does well in simulating the flow phenomenon of tailings after a dam break, and the numerical results can be used for early warning of disasters and emergency response.

1 Introduction

A tailing dam is one of the three basic projects at a mine, and it is a dangerous source of artificial mud flow with high potential energy. The stability of tailing dams is highly susceptible to adverse natural and human-induced factors. According to the International Commission on Large Dams (ICOLD) (Dams, 2001), the accident rate of tailing dams in the past 100 years (1900–2000) was approximately 1.2 %, which is 2 orders of magnitude higher than the accident rate of reservoirs (0.01 %). Since the beginning of the 20th century, there have been hundreds of tailing dam accidents in the world (Rico et al., 2008b). In 2008, the Xiang Fen tailing dam collapsed in China, causing 381 deaths. In 2015, the Fundão mine tailing dam failure in Brazil released 32 million m³ of toxic tailings, polluted 650 km of rivers and flowed into the Atlantic Ocean, causing serious environmental consequences (Burritt and Christ, 2018). About 4 years later, more than 232 people died in the failure of the Feijão tailing dam, which belongs to the same company as the Fundão mine (Santamarina et al., 2019).

To date, the number of tailing dams in the world exceeds 10⁵, and ensuring the stability of such reservoirs is an extremely challenging task for the mining industry. By documenting 147 historical cases of tailing dam failure worldwide, Rico et al. (2008b) found that the most common causes of failure are related to unusual rain and seismic liquefaction. Furthermore, dozens of factors, such as management operations, slope instability, fluvial undermining and foundation failure, are also likely to cause dam failures. As a consequence, even if all the tailing dams were constructed correctly, the occurrence of pond failure cannot be avoided.

Dam-break simulation and run-out analyses can provide useful information for assessing the risk associated with tailing dam collapses, such as run-out path and quantification of potential losses, which have become more important in the geotechnical assessment of tailing dams. Nevertheless, little work has been conducted to analyze the dam-break process and investigate the flow of tailings released following dam failure. These related dam-break analyses can be divided into empirical models and numerical models. The empirical model predicts the volume of released tailings (V_F) and the run-out distance (D_max) by establishing a regression equation based on historical cases of collapse, dam height and the impounded volume of tailings (Rico et al., 2008a; Larrauri and Lall, 2018). However, owing to the difference in the actual state of the tailing reservoir and the surrounding environment, and without considering the rheological properties of the tailings such as the viscosity coefficient and yield stress (Martin et al., 2015), the predicted results may deviate greatly from the actual dam-break situation. Additionally, the flow pattern and movement regularity of tailing fluid cannot be accurately described.

Since numerical modeling based on computational fluid dynamics (CFD) has the advantages of sophisticated theory and low cost, it is considered to be a more appropriate solution for analyzing geographic hazards such as water reservoir failures, landslides and debris flows (Dai et al., 2014; Issakhov and Imanberdiyeva, 2019; Han et al., 2019). Using CFD to study and analyze a tailing run-out path is rare. In many cases, the depth-averaged shallow water model has been widely used to compute the fluxes on rectangular grids. It can be numerically solved by different schemes such as smoothed particle hydrodynamics (SPH) (Cascini et al., 2014), the finite volume method (FVM) (Medina et al., 2008; Li et al., 2013; Wang et al., 2020) and the finite difference method (FDM) (Wu et al., 2013). However, the shallow water equations are two-dimensional simplifications of Navier–Stokes equations, neglecting the observed nonlinear behavior through the mud depth (Han et al., 2019). For more reliable results, the 3-D Navier–Stokes equations are recommended to control the flow dynamic (Wang et al., 2016; Lee et al., 2010).

One key issue limiting the application of the CFD method in tailing flow simulation is the rheological properties. Unlike water, which is a Newtonian fluid, the tailing fluid formed by mixing tailings and water is in nature a non-Newtonian fluid (Henriquez and Simms, 2009) with complex rheological properties. In the early research, the Voellmy frictional model, the frictional model and the general two-phase model were used to describe the rheological behavior of a mixture of water- and clay-like tailing fluids and were successfully applied in 2-D shallow water equations, but these constitutive models may present difficulty in discretizing in the 3-D FVM scheme (Han et al., 2019). Currently, the Bingham model is considered a better way to describe the dynamic behavior of tailing fluids and easy to apply in a 3-D case (Mizani et al., 2013; Henriquez and Simms, 2009; Liao and Zhou, 2015; Gao and Fourie, 2015). Another key factor that makes it difficult to analyze and predict tailing run-out routing by numerical modeling is the difficulty of building the boundary of large-scale complex downstream terrains. Recently, some initial work, such as laboratory experiments and relatively simple deposition simulations of tailing fluid, has been attempted (Mizani et al., 2013; Zhang et al., 2015; Babaoglu and Simms, 2017; Gao and Fourie, 2019; Luo et al., 2017). These studies discussed the flow and rheological properties of tailing fluids in detail but did not take into account the influence of real terrain on the movement distance of tailing fluids. Wang et al. (2018) integrated the SPH and the 30 m resolution digital elevation model (DEM) of the ALOS satellite to simulate the dam-break process of a tailing reservoir, while the complex rheological properties of the tailing fluids were not mentioned at all.

In this paper, a three-dimensional (3-D) model based on the FVM method for simulating mud flow from a tailing dam break across 3-D terrain was presented to provide quantitative information for risk assessment and disaster prevention, and the accuracy of the model was tested using small-scale laboratory experiments of a dam-break-type flow. Aiming at the difficulty of accurately reproducing large-scale downstream terrain, unmanned aerial vehicle (UAV) aerial photogrammetry technology was used to obtain fine DEM data. Furthermore, the calculation of the model was efficiently executed on the supercomputer cluster of the Fujian high-performance computing (HPC) center. Finally, we simulated two scenarios of tailing dam failure.

2 Methodology

2.1 Mathematical and numerical models

In this study, the mixture of tailings and water is regarded as an incompressible fluid, and the tailing-flow dynamic behavior can be described by the 3-D continuity and Navier–Stokes (N–S) governing equations, which follow the physics laws of conservation of mass and momentum. The governing equations (Versteeg and Malalasekera, 2007) can be written in the tensor form as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }{u}_i\right)=\mathrm{0}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }{u}_i\right)& +{\displaystyle \frac{\partial }{\partial x_j}}\left(\mathit{\rho }{u}_i{u}_j\right)=\\ & -{\displaystyle \frac{\partial p}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left({\mathit{\mu }}_{\mathrm{e}}\left({\displaystyle \frac{{\partial }_{\stackrel{\mathrm{‾}}{u_i}}}{\partial x_j}}+{\displaystyle \frac{{\partial }_{\stackrel{\mathrm{‾}}{u_j}}}{\partial x_i}}\right)\right),\end{array}\end{array}$$

where ρ is the density, t is the time, u_i is the velocity, p is the pressure and μ_e is the effective viscosity.

Most natural or engineering fluids are in a turbulent state that causes fluctuations in velocity, pressure, temperature and other transport scalars in space and time. Tailing reservoirs are mostly located in mountainous areas, with a large terrain gradient. Under the influence of gravity, tailing fluids flow out of the break at high speed and will be in a turbulent state (Blight, 1997). The standard k-ε Reynolds-averaged Navier–Stokes (RANS) equations (Jones and Launder, 1972) are used to address the turbulence effect. This method introduces the kinetic energy k and the turbulent dissipation rate ε with good convergence, which is the most widely used method in the engineering field. The RANS equations can be written as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }\stackrel{\mathrm{‾}}{u_i}\right)& +{\displaystyle \frac{\partial }{\partial x_j}}\left(\mathit{\rho }\stackrel{\mathrm{‾}}{u_i{u}_j}\right)=\\ & -{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left({\mathit{\mu }}_{\mathrm{e}}\left({\displaystyle \frac{{\partial }_{\stackrel{\mathrm{‾}}{u_i}}}{\partial x_j}}+{\displaystyle \frac{{\partial }_{\stackrel{\mathrm{‾}}{u_j}}}{\partial x_i}}\right)-\mathit{\rho }\stackrel{\mathrm{‾}}{u_i^{\prime }}\stackrel{\mathrm{‾}}{u_j^{\prime }}\right),\end{array}\end{array}$$

where – denotes the average value and ^′ denotes the pulsation value.

In the present model, the tailing fluid and air are assumed to be two incompressible phases, and the surface of the tailing phase movement is captured using the volume-of-fluid (VOF) method (Hirt and Nichols, 1981) that uses a step function F whose value is unity at any control volume occupied by tailing fluid and zero otherwise. The time dependence of F is governed by Eq. (4):

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathit{F}}{\partial t}}+\mathrm{\nabla }\cdot \left(u_i\mathit{F}\right)=\mathrm{0}.\end{array}$$

The above governing equations are solved based on the FVM, which is the most commonly used method in the field of fluid engineering (Schraml et al., 2015; Yoon et al., 2014). We split the computational region into series of control volumes, which is represented by a mesh of hexahedral cells in the Eulerian framework according to the FVM. The governing equations are discretized into a set of algebraic expressions involving the unknown values of the physical quantities at the center of each control volume by integrating the N–S equation over the control volume. The following high-order discretization schemes (Table 1) are used to interpolate the physical quantities of the centers of the control volumes onto the faces, and this setup is second-order accurate and fully bounded.

Table 1 Discretization schemes of the model.

Here ∇ is the gradient operator, ∇⋅ is the nabla operator and ∇² is the Laplacian operator.

Preconditioning conjugate gradients (PCGs) and smooth linear algebraic solvers are selected to calculate symmetric and asymmetric matrices separately in the work. Additionally, in terms of pressure–velocity coupling, we use the pressure implicit with the splitting of operators (PISO) (Issa et al., 1991) algorithm to achieve it. The PISO algorithm uses prediction–correction–recorrection steps, speeding up the convergence speed in the iterative process and improving computational efficiency. To guarantee model stability, the Courant–Friedrichs–Lewy (CFL) condition must be satisfied. For the n-dimensional case, the CFL number is defined in the following general form:

$$\begin{array}{}\text{(5)}& \mathrm{CFL}=\mathrm{\Delta }t\sum _{i=\mathrm{1}}^n{\displaystyle \frac{u_i}{\mathrm{\Delta }{x}_i}}.\end{array}$$

The CFL number is a measure of how much information (u_i) traverses a control volume (Δx_i) in a given time step (Δt). The CFL condition can be restricted by a maximum CFL coefficient between 0 and 1, which is chosen to be 0.5 in the present work.

2.2 Rheological behavior of tailing fluid

Tailing fluid is similar to debris flow and is generally considered to be an incompressible and viscous fluid. Currently, most existing studies on debris flows and landslides assume that the mud–water mixture is a homogeneous non-Newtonian fluid with a constant density (Dai et al., 2014; Wang et al., 2016; Han et al., 2019). Several authors concluded that a fluid with solid concentrations above 10 % is a non-Newtonian fluid (Komatina and Jovanovic, 1997; Sloff, 1993), and their rheological properties can be described by the Bingham model (Pastor et al., 2014; Komatina and Jovanovic, 1997; Han et al., 2019) that can be represented by the following constitutive equation:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{ij}=\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{\dot{\mathit{\gamma }}}}+{\mathit{\mu }}_{\mathrm{b}}\right){\dot{\mathit{\gamma }}}_{ij},\end{array}$$

where τ_b is the yield stress, μ_b is the viscosity coefficient, τ_ij is the shear stress tensor and ${\dot{\mathit{\gamma }}}_{ij}$ is the rate-of-deformation tensor:

$$\begin{array}{}\text{(7)}& {\dot{\mathit{\gamma }}}_{ij}={\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}.\end{array}$$

The symbol $\dot{\mathit{\gamma }}$ is the second invariant of the ${\dot{\mathit{\gamma }}}_{ij}$, which is given by

$$\begin{array}{}\text{(8)}& \dot{\mathit{\gamma }}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\dot{\mathit{\gamma }}}_{ij}{\dot{\mathit{\gamma }}}_{ij}}.\end{array}$$

Bingham fluid is a kind of plastic fluid; as long as the shear stress exceeds the threshold τ_b specific to the material, the material flows like a Newtonian fluid. Where the threshold is not exceeded, the plastic fluid behaves like a solid. The average diameter of the tailings is generally less than 0.1 mm. At present, several laboratory rheological tests have shown that the rheological properties of a water–tailing mixture conform to those of a Bingham fluid (Mizani et al., 2013; Henriquez and Simms, 2009; Liao and Zhou, 2015; Gao and Fourie, 2015); hence, the Bingham model is adopted as the constitutive model to describe the fluidization characteristics of flow-like tailings in this work. In the Bingham model, the unknown parameters are τ_b and μ_b, which can be determined by laboratory measurements. The apparent viscosity in a Bingham model can be represented as

$$\begin{array}{}\text{(9)}& {\mathit{\mu }}_{\mathrm{e}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{\dot{\mathit{\gamma }}}}+{\mathit{\mu }}_{\mathrm{b}}.\end{array}$$

Understandably, the effective viscosity will become unbound when $\dot{\mathit{\gamma }}$ tends to be infinitely small, such as in the core region of a Bingham flow, which may cause numerical divergence or even crash in the solution procedure (Shao and Lo, 2003). To overcome this inherent discontinuity issue, several regularized approaches have been proposed, where the infinite viscosity presented in the rigid zone can be approximated by a highly viscous fluid (Papanastasiou, 1987; Hannani et al., 2007; Frigaard and Nouar, 2005). The most popular regularization is the Bingham–Papanastasiou model (Papanastasiou, 1987), which is chosen as the constitutive equation of tailing fluid in the present work. Hence, the apparent viscosity is expressed as

$$\begin{array}{}\text{(10)}& {\mathit{\mu }}_{\mathrm{e}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{\dot{\mathit{\gamma }}}}\left[\mathrm{1}-e^{(-q\dot{\mathit{\gamma }})}\right]+{\mathit{\mu }}_{\mathrm{b}},\end{array}$$

where q denotes the stress growth parameter, such that exceeding the yield stress τ_b finite shear stress is allowed for small shear rates, while the stress grows linearly with shear stress beyond the yield stress. Moreover, the limit of m→∞ is fully equivalent to an ideal Bingham fluid.

In this study, mathematical and rheological models are numerically built and implemented using the release of version 6 (https://github.com/OpenFOAM/OpenFOAM-6, last access: 3 March 2020) of the open-source Field Operation and Manipulation (OpenFOAM) C++ libraries under a Linux environment. OpenFOAM is a framework for developing application executables that use packaged functionality contained within a collection of approximately 100 C++ libraries (Greenshields, 2018). In the original OpenFOAM code, it is impossible to cope with the Bingham flow. To make it possible for the tailing-flow simulation, the regularized Bingham–Papanastasiou model without numerical divergence was introduced into the original code.

2.3 Experimental verification

To investigate the performance of the model in predicting the behavior of a Bingham fluid, we compared the simulation results with the analytical solutions and experimental results for the flow of a Bingham fluid between two parallel plates and dam-break-type flow. The purpose of the flow test of a Bingham fluid between two parallel plates is to check whether the Bingham–Papanastasiou regularization method can correctly describe the flow phenomena of a Bingham fluid, and the 3-D dam-break experiment is used to test the ability of the above scheme to capture the free surface of a fluid and calculate pressure and velocity.

2.3.1 Flow between two parallel plates

The test of flow between two parallel plates was used for the validation of the implementation of the rheological model in OpenFOAM since the analytical solution can be easily acquired from

$$\begin{array}{}\text{(11)}& u\left(y\right)={\displaystyle \frac{\mathrm{\Delta }PD}{\mathrm{2}{\mathit{\mu }}_{\mathrm{b}}L}}\left(y-{\displaystyle \frac{y^{\mathrm{2}}}{D}}\right)-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}y}{{\mathit{\mu }}_{\mathrm{b}}}}\mathrm{for}\mathrm{0}\le y<S.\end{array}$$

Equation (11) gives the expression for velocity for τ>τ_b (Gopala et al., 2011; Waarde, 2007), and the velocity at the plug is given by Eq. (12):

$$\begin{array}{}\text{(12)}& u_{\mathrm{plug}}=u\left(S\right)\mathrm{for}S\le y<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}D.\end{array}$$

For the simulation of flow between parallel plates, a 2-D computational domain was constructed, which consisted of a rectangular channel with length L=2 m and width D=0.1 m. The test configuration is schematically shown in Fig. 1. The domain was discretized with 400×40 cells limited by three types of boundary faces: inlet, outlet and no-slip walls, and the Bingham fluid flowed from the inlet to the outlet at a constant velocity of 0.069 m s⁻¹. Two types of Bingham fluids were simulated, namely, a self-compacting cement (SCC) mixture and grout, and their rheological parameters are shown in Table 2 with reference to Gopala et al. (2011). In Fig. 2, the velocity profiles were compared between the numerical simulation and the analytical solution for the SCC and the grout. The simulated results and analytical solutions appear to be in very close agreement, with the implemented results of the Bingham–Papanastasiou model fitting the analytical velocity curves very well. In terms of velocity, the average error of simulated values compared with analytical solutions is 0.000077 m s⁻¹, and the maximum error is 0.0016 m s⁻¹. For the plug flow region, the yield stress of the SCC is 2 orders of magnitude higher than that of the grout, so there is an obvious plug region in the flow of the SCC, as shown in Fig. 2.

Figure 1 Schematic of the flow between parallel plates.

Figure 2 Comparison of velocity variation between analytical solutions (Gopala et al., 2011) and numerical results simulated in two parallel plates based on different rheological parameters, i.e., (a) SCC and (b) grout.

2.3.2 3-D dam-break experiment

The second experiment presented here is one of the dam-break-type problems to validate the model. The dam-break problem is a classical validation case for the assessment of free surface modeling methods. As shown in Fig. 3, the experiment is carried out in a large tank of 3.22 m × 1 m × 1 m with an open roof, and the right part of the tank is a rectangular box of water between a fixed wall and a temporary door. At time t=0 s, the temporary door is removed instantaneously, allowing the fluid to collapse under the influence of gravity g=9.81 m s⁻². During the collapse, the fluid impacts an obstacle at the bottom of the tank and creates a complicated flow structure. The geometry and the probe positions are briefly described in Fig. 3; points P1 to P8 are pressure sensors on the obstacle box, while probes H1 to H4 are used to observe the fluid heights. In this simulation, grids are used in a structural pattern with an initial spacing of 0.015 m, so a set of 971 498 cells is employed, and the boundary consists of the no-slip walls and the top boundary that is free to the atmosphere so that both outflow and inflow are permitted according to the internal flow. For the rheology, the typical values used in the simulations are μ=0.001 Pa s, τ_b=0 Pa and ρ=1000 kg m⁻³ for water-based suspensions at 20 ^∘C.

Table 2 Rheological parameters of a simulated Bingham fluid.

Figure 3 Measurement positions for water heights and pressures in the dam-break experiment (from Kleefsman et al., 2005). P1 to P8 are pressure probes, and H1 to H4 are height probes.

Figure 4 Screenshots of the simulation compared with the experiment (Kleefsman et al., 2005) at t=0.4 s (a) and t=0.56 s (b). The color map of the simulation results is mapped by velocity characteristics.

The experiment was performed by the Maritime Research Institute Netherlands (MARIN), and all detailed experimental data can be downloaded free from this website (https://app.spheric-sph.org/sites/spheric/tests/test-2, last access: 3 March 2020). In this case, we compared the simulation results against the experimental data obtained by MARIN. The free surface and velocity field of the simulation against different physical times t=0.4 s and t=0.56 s are shown in Fig. 4, which is similar to the results in Kleefsman et al. (2005) and Lee et al. (2010). The small pictures on the top right of the screenshots represent the fluid behind the gate. As this figure suggests, the simulation results are in good agreement with the experimental data. In Fig. 5, characteristic curves of pressure and vertical water heights varying with time at P1 and H4, respectively, are depicted. Similarly, the global behaviors of the simulation and experiment are fairly consistent but slightly shifted about 0.2 s after 2.55 s at probe H4 and after 4.6 s at probe P1, as shown in Fig. 5. Kleefsman et al. (2005) pointed out that an underpredicted difference could be reduced by increasing the grid resolution in the comments of their VOF simulations. At probe H4, the water height experienced a period of decline from the initial 0.55 m, and the experimental and simulated values reached a minimum of about 0.1 and 0.075 m in approximately 0.25 and 0.27 s, respectively. At the pressure probe P1, the difference between the simulated pressure value and the experimental value reaches a maximum of about 1617 Pa when the pressure reaches a maximum at about 0.41 s, and the average error in other times is 115 Pa. Consequently, the ability of the above scheme to capture the free surface of a liquid and calculate the pressure and flow rate is indicated.

Figure 5 Comparison of vertical water heights in the reservoir at probe H4 (a) and the characteristic curve of pressure changing with time at probe P1 (b) between experimental data and simulation results.

3 Simulation of tailing dam failures – Feijão and A'xi tailing dams

3.1 Feijão tailing dam

The Córrego do Feijão iron mine complex, with two tailing dams, a cargo terminal, an administrative office and a small railway network for the transport of iron ore (Fig. 6a), is located in Brumadinho, Minas Gerais state, Brazil, in the upper valley of the Paraopeba River. The tailing dam I, measuring 86 m high and 720 m long, was used for disposal of approximately 12.7 million m³ of iron tailings. On 25 January 2019, dam I suffered structural damage and eventually collapsed catastrophically, causing at least 232 deaths (Santamarina et al., 2019). According to a field survey (Porsani et al., 2019) that was conducted after the collapse, approximately 11.7 million m³ of tailings were spilled from the dam, destroying office buildings, a railway bridge and a small community. Under the action of gravity, the tailing fluids traveled 9 km and eventually flowed into the Paraopeba River, the region's main river. The downstream area submerged by tailing fluids was more than 2.54 million m² based on measurements from remote sensing images from Sentinel Satellite S2 (see Fig. 6b). According to surveillance cameras around the accident area, the velocity of the tailing fluid when the tailing dam failed was more than 10 m s⁻¹. Figure 6 shows the area impacted by the dam I collapse event.

Figure 6 Location of the Feijão iron mine complex with reconstructions of the area impacted by the 2019 dam I failure event: (a) location and downstream area of the Feijão tailing dam. A – location of dam I. B – location of dam VI which was destabilized by the event. C & E – locations of the rail network and bridge. D – location of the office center and a small community. F – location of the Paraopeba River. (b) G – the destroyed downstream area. Map data: © OpenStreetMap contributors 2019. Distributed under a Creative Commons BY-SA License and Sentinel Satellite S2 imagery.

3.2 A'xi gold tailing dam

The A'xi gold mine tailing reservoir (Fig. 7) is 58 m high and located 46 km from the Ili River, Xinjiang, China, in the upper valley of the A'xi River, a tributary of the Ili River. Designed for disposal of approximately 3.6 million m³ of gold tailings, this tailing dam measures 210 m long and occupies 170 000 m² and is the largest tailing dam in Xinjiang. The tailing storage facility of the A'xi gold mine was constructed using the “paddock” system of construction, where the subsequent dams raised were built on the tailings. In addition, it is located on the south slope of the main ridge of Kogurqin Mountain in the Borokunu Mountains of the western Tian Shan with an elevation of 1513 m and a large downstream topographic gradient. The mining area belongs to the continental climate of the north temperate zone, where the overall rainfall is low but unusual rainfall events occur frequently. At present, the tailing dam has entered the tail period of its service life with a storage capacity of approximately 2 million m³ of tailings. Additionally, the area is in an earthquake-prone zone and is prone to abnormal rainfall. Therefore, the possibility of tailing dam failure cannot be ruled out. If a tailing accident occurs, it will cause environmental disasters and large property losses and even pollute the large reservoirs downstream. Considering the serious consequences if the A'xi gold mine tailing dam experiences an accident, we chose it as a disaster assessment case study to simulate and analyze its dam-break process and run-out path.

Figure 7 Location of A'xi gold mine tailing dam: the red area is the A'xi gold mine complex; the blue line represents the valley of the A'xi River; the green area represents a large water reservoir. Image data: Sentinel Satellite S2 imagery and UAV imagery.

3.3 Key factors and model setup

The presented and validated CFD model was used to simulate the run-out path and range of tailing fluid from the Feijão tailing dam failure and to simulate the breaking of the A'xi tailing dam. The features of run-out are greatly affected by three significant factors: the fineness of the downstream topography, the rheological characteristics of the tailing fluid and the initial state of the tailing fluid.

Figure 8 The 3-D terrain model around the A'xi Gold Mine. (a) The fine DEM generated by UAV (right) and ALOS satellite (left), (b) reconstructed 3-D terrain model, and (c) the computational meshes.

First, the terrain has a great influence on the tailing fluid's velocity and run-out path due to the gradient variations in the terrain. The fine terrain boundary, which is closer to the real landform, has more features and obvious topographic fluctuations than the coarse boundary, making the numerical simulation results more in line with the realistic phenomenon. Conversely, the coarse terrain boundary has a greater loss of detailed topographic characteristics. If the resolution of the terrain boundary is too low, it may cause a large difference between the simulation results and the real event. The topographical boundaries of a large area possibly affected by the tailing dam failure were constructed using a DEM. For the A'xi tailing dam, we went to the field to collect a group of aerial photographs with a certain overlapping range using the UAV. Based on these overlapped images, a DEM was reconstructed with a resolution of approximately 0.5 m × 0.5 m, which can be further generated into a 3-D terrain model. However, due to the large altitude drop and the limited endurance of the drone used in this area, images of some areas that may be affected by tailing fluids have not been collected. For these areas, we chose the ALOS PALSAR RTC (radiometrically terrain corrected) DEM with a 12.5 m × 12.5 m resolution released by the Alaska Satellite Facility (https://www.asf.alaska.edu/sar-data/palsar/terrain-corrected-rtc/, last access: 3 March 2020) to replace them. Then, the UAV DEM is mosaicked into the ALOS DEM and forms a broader DEM with of cell size of 0.5 m × 0.5 m. After that, the inlaid DEM was georeferenced, cropped and reconstructed to form the final 3-D terrain geometry around the A'xi gold mine. Because the Feijão tailing dam has collapsed, the DEM data collected by the drone before the accident could not be obtained; therefore, we also used the ALOS PALSAR RTC DEM to make the 3-D terrain geometry. Figure 8 shows the DEM of the A'xi tailing reservoir and the reconstructed three-dimensional model computational grid. The 3-D model of the terrain was built to make the meshes of the computational domain, which were the data that must be used in numerical calculations based on the FVM. For mesh generation, an OpenFOAM utility, snappyHexMesh, was employed to generate 3-D meshes containing hexahedra and split-hexahedra automatically from a 3-D file (STL or OBJ format).

Furthermore, the travel distance and extent of a tailing outflow are affected by the rheological equation. Several studies (Pastor et al., 2004; Henriquez and Simms, 2009; Gao and Fourie, 2015, 2019; Babaoglu and Simms, 2017) have established that Bingham constitutive law provides a reasonably good compromise between accuracy and simplicity, and it has been frequently adopted to describe the rheology of tailing fluid. With respect to the determination of the rheological parameters, Henriquez and Simms (2009) determined the yield stress and viscosity of tailing flow using rheometer and slump tests; the results showed that yield stress ranges roughly from 40 to 300 Pa, and dynamic viscosity ranges from approximately 0.1 to 1.12 Pa s when the solid concentration is between 60 % and 70 %. Gao and Fourie (2015) concluded that the range of yield stress values is from approximately 18 to 60 Pa, which covers the majority of cases of actual tailing disposal operations, with a viscosity ranging from around 0.3 to 0.8 Pa s. The yield stress of the tailing flow generally has more influence on the final profile than does velocity during the deposition process in flume tests, while the importance of viscosity for the formation of the final profile of tailings increases with a rise of inertia effects (Gao and Fourie, 2019). In this study, the Bingham viscosity and yield stress were determined based on a previous investigation report conducted in testing the rheological properties of tailings (Liao and Zhou, 2015). Liao and Zhou (2015) collected tailing samples from three different tailing ponds and mixed these tailings with water to form different concentrations of tailing fluid for rheological tests. A total of 15 rheological tests were carried out using a LAMY RM100 rotational viscometer, with a range of dams (no. 1–no. 3) and volume concentrations of mixtures of 60 %–80 %. The no. 1 tailing sample and no. 2 tailing sample were mainly fine sand with a max particle size of less than 2 mm, but the powder content (0.075 mm < diameter <0.25 mm) of no. 1 tailing sand (68.68 %) was more than that of no. 2 tailing sand (42.95 %). The no. 3 tailing sample was mainly composed of coarse sand with a max particle size of less than 5 mm. Because Feijão iron tailings and A'xi tailings are similar to no. 1 tailings, they are both silt and have a small particle size, and the main mineral composition of no. 1 tailings is iron ore. In this work, the 80 % solid concentration of the no. 1 tailing dam is selected, corresponding to a density of 1830 kg m⁻³, a viscosity of 0.741 Pa s and yield stress of 59.82 Pa. The values of this set of parameters are in accordance with the recommended viscosity and yield stress range of Henriquez and Simms (2009) and Gao and Fourie (2015).

Figure 9 The time sequence of the numerically computed event up to the tailings reaching the Paraopeba River. The base map is the satellite image (Sentinel S2) after the dam break. The color map of the simulation results is mapped by velocity characteristics.

Undoubtedly, the initial conditions of the released tailing fluid can have an important effect on the velocity and out-flow distance of the tailings, especially at the time of the dam break, which is mainly caused by the initial geometry and volume of the tailings. For the Feijão tailing dam, according to the monitoring video of the accident area, the structure of the dam body completely collapsed instantaneously, and the same hypothesis was also used for the A'xi tailing pond. The volume of tailings released by the Feijão tailing accident was set to 11.7 million m³ based on subsequent reports. While the Feijão tailing dam almost discharged all the tailings it had stored, tailing dam accidents will generally only release part of the tailings stored, according to many historical accidents. The volume of tailings released after the dam break at the A'xi gold mine tailing reservoir is estimated to be approximately 756 000 m³ based on the regression equation proposed by Larrauri and Lall (2018):

$$\begin{array}{}\text{(13)}& \log \left(V_{\mathrm{F}}\right)=-\mathrm{0.477}+\mathrm{0.954}\times \log \left(V_{\mathrm{T}}\right),\end{array}$$

where V_T is the total impounded volume and V_F is the volume of tailings that could potentially be released. For the initial 3-D model of the tailing dam, its volume is equal to the total volume of tailings released, which is constructed based on the dam height and geometric shape of the tailing pond. Formerly, the 3-D topographic model was modified to suit the initial geometry of the tailings.

Generally, a smaller cell spacing provides a better simulation result, but it will dramatically increase the time for calculation. In our simulation, the cell spacing of the Feijão tailing dam and the A'xi tailing dam are fixed to 10 m × 10 m × 3 m and 3 m × 3 m × 3 m, generating 3 241 998 cells and 6 657 418 cells, respectively. The terrain elevation within each cell is evenly distributed, so it is necessary to supplement relatively accurate surface information such as terrain roughness. Owing to the difference in the surface cover, the roughness length of the downstream terrain of the Feijão tailing dam and the A'xi tailing dam is set to 1 and 0.1 m, respectively, based on a criterion given by Henderson-Sellers et al. (1993). Finally, the geometry and associated fields are decomposed into pieces using the decomposePar utility, and each separate part of the decomposed domain was run on the cluster of the Fujian HPC by using the IntelMPI (an implementation of the standard message passing interface) and Slurm workload manager.

3.4 Results and discussion

The proposed 3-D CFD approach was used to simulate the dynamic behavior of tailing flow over the 3-D topography based on the above settings, and the discussion of the propagation features of the tailing outflow is presented as follows.

Figure 10 Displacement time history for the front of Feijão failure event in the simulation.

Figure 11 Average velocity time history for the free surface simulated of the Feijão failure event.

For the Feijão tailing reservoir, we simulated the travel of tailings for 2500 s at the time of the disaster. The typical simulation snapshots of the fluidization process of the 2019 Feijão dam-break event are presented in Fig. 9. In addition, to analyze the movement characteristics of the tailings after failure, the progress of the front displacement within time sequence is shown in Fig. 10. After dam failure, the gravitational potential energy of the tailing fluid was converted into kinetic energy, causing the velocity to increase suddenly, even more than 20 m s⁻¹. Subsequently, it took 150 s for the tailing fluid to travel approximately 1300 m and reach the office center and small community, with a velocity of 15 m s⁻¹. After 150 s, the front flow smashed the railway bridge at approximately 3150 m downstream, with a decreasing average velocity of less than 5 m s⁻¹. Finally, it flowed into the Paraopeba River at approximately 2400 s, with a small average velocity of less than 2 m s⁻¹, and the total travel distance was approximately 9 km. Figure 11 shows the mean velocity time history of the free surface. It can be obtained that the average velocity of the tailing flow increases continuously to more than 6 m s⁻¹ in the first 100 s, and then it gradually slows down, and the average speed stabilizes at around 1 m s⁻¹ from 1000 s. Furthermore, to check the quality of the results simulated on the Feijão event, a comparison of the CFD simulated and measured submerging area is shown in Fig. 12, from which we can see that the numerical simulation results of the tailing flow directions and impact area are also modeled satisfactorily. There is a significant difference in the location of the railway network (see Figs. 6a and 12), which is mainly due to the acquisition time of the DEM used in 2011 when the railway network had not been built. The simulated inundation area is about 2.572 km², which is only 0.036 km² different from the actual measured inundation area of 2.536 km². It can be observed that the simulation results are fairly consistent with the on-site historical runoff area, with only minor differences in some locations, and we believe that more accurate DEM data and smaller grid size could be helpful to obtain a better result. However, the collapse of the tailing pond is a very complicated geographical phenomenon. It is impossible to fully satisfy the historical situation in the simulation because uncertainty factors are not considered, so it might not be closely matched in terms of time.

Figure 12 Comparison of simulation results and measurement results in the area destroyed by the Feijão accident.

Figure 13 The time sequence of the numerically computed event up to the tailings reaching the left boundary of the computational domain. The base map is the satellite image (Sentinel S2). The color map of the simulation results is mapped by velocity characteristics.

According to the relevant settings in Sect. 3.3, if the A'xi gold mine tailing pond experiences a dam failure, the flow path and submerged range of the tailings were predicted within 800 s. The simulated propagation of the tailing fluid released from the A'xi tailing dam failure is shown in Fig. 13. In addition, to analyze the movement characteristics of the tailings after failure, the progress of the front displacement within time is shown in Fig. 14. After the tailing fluid is released from the breach, it flows down along the terrain under the action of gravity. At 20 s, it first submerges the environmental protection depot and power distribution room (see Fig. 13) approximately 250 m downstream of the tailing pond, and the speed exceeds 15 m s⁻¹. Immediately afterward, the tailing fluid enters the A'xi Valley, approximately 320 m from the dam site, and flows downstream along the terrain in the valley. Although the elevation difference in the valley is more than 155 m (in the calculation domain), the valley is so tortuous that the flow rate of the tailing fluid gradually decreases. Eventually, at 760 s, the front reaches the simulated boundary with an average velocity of less than 2 m s⁻¹ and a flow distance of approximately 3.5 km. Figure 15 shows the mean velocity time history of the free surface. It can be obtained that the average velocity of the tailing flow increases continuously to more than 11 m s⁻¹ along the steep channel in the first 10 s, and then it gradually slows down due to the continuous curved valley, and the average speed stabilizes at around 2 m s⁻¹ from 200 to 550 s. After that, there is a slight increase in speed due to the large topographic elevation changes. It can be seen from the numerical results that almost 75 600 m³ of tailings released after the break of the A'xi gold mine tailing dam is almost poured into the A'xi Valley except that some of the tailings are blocked by the environmental protection depot. The normal runoff volume of the A'xi River is 0.5–2 m³ s⁻¹, and the runoff volume during the dry season is 0.002 m³ s⁻¹. Overall, the runoff volume is relatively small, so the influence of water flow is not considered in the simulation. However, the region is prone to abnormal rainfall events, which lead to an increase in river runoff. Under this situation, the toxic tailings that flow into the river valley would be transported by the river to the large reservoir downstream, or even eventually to the Ili River, causing serious environmental pollution.

Figure 14 Displacement time history for the front flow of A'xi failure in the simulation.

Figure 15 Average velocity time history for the free surface simulation of the A'xi dam failure.

From the simulated results shown above, it can reasonably be concluded that the proposed method can accurately represent the whole dynamic processes of tailing movement after dam failure. In addition, parallel technology was used to speed up computation. We tested the parallel speedup for different numbers of processors on the supercomputing cluster of the Fujian HPC, with a computational load of more than 3 million cells, and the speedup was nearly a factor of 50. When the number of processors was more than 50, the calculation time was no longer reduced because the acceleration ratio also depends on the amount of computation. If the number of cells increases, the speedup will be greater. Figure 16 shows the parallel speedup of different numbers of processors.

Figure 16 Parallel speedup with different numbers of processors. The green line is the ideal speedup, and the purple line is the actual speedup.

4 Conclusions

The damage caused by tailing dam failures occurring frequently around the globe can be costly and catastrophic. The topic of disaster prevention and mitigation is being discussed more in government departments and the mining industry. The debris flow from a tailing dam failure is a very complex phenomenon involving physical mechanisms and the evolution of fluid properties. In this paper, a 3-D CFD model based on the FVM was proposed to predict the routing and impact area of tailings across complex topography during tailing dam failure. The Navier–Stokes equations and the Bingham rheology model are introduced into the FVM framework as the governing equations and the constitutive model, respectively. To construct a wide range of detailed downstream terrain, UAV photogrammetry and satellite imagery are used to sample surface information. Moreover, parallel technology was used to simulate the complex problem of the supercomputing cluster to achieve rapid and reasonable predictions for disaster prevention and mitigation.

To verify the performance of the reported scheme, the numerical model is first tested against an analytical solution in the literature for the flow of a Bingham fluid between two parallel plates, demonstrating the correct implementation of the Bingham–Papanastasiou constitutive model in OpenFOAM. Next, the numerical model is validated against the experimental value for the 3-D dam-break experiment, thereby ensuring a good prediction of the free surface profile. To further illustrate the performance of the presented method, the 2019 Feijão dam-break event in Brazil was selected as a case study, and the numerical simulation results coincide well with the in situ investigation. Finally, considering the serious consequences of the A'xi gold mine tailing dam if it experiences an accident, we chose it as a disaster assessment case study to simulate and analyze its dam-break process and run-out path. Therefore, CFD numerical modeling can provide an effective means for mapping hazardous areas, estimating hazard intensity and designing appropriate protective measures.