Numerical Simulation on the Evolution of Tailings Pond Dam Failure Based on GDEM Method

Because of the continuous exploitation of metal minerals and the limited availability of land resources, many tailings ponds have been expanded. It is of great importance to dynamically calculate the critical parameters of flooding and sand velocity after a collapse, and quickly determine the dangerous area downstream on reducing the hidden dangers for an accident. In view of the limitations of the single finite element and discrete element methods for the simulation of tailings pond instability, the authors use the finite element and discrete element coupling method (GDEM) for the first time to carry out a dynamical process for numerical simulation of the evolution of a tailings pond collapse. The sediment movement and submerged range were compared with the results of the Volume of Fluid (VOF) method. The comparison shows that the change in drainage flow at the dam foundation position and the depth change in the downstream sensitive point after the dam break are consistent with the VOF calculation results, which indicates that the GDEM can be used as an effective means for the analysis of tailings dam failure. The calculation results of this method can provide technical support for the determination of dangerous areas downstream of tailings ponds and the safety management of key areas.


INTRODUCTION
The majority of China's metal deposit resources are deeply buried low-grade ores, while high-grade ores are relatively scarce. As a result, a large number of tailings ponds are required to store mining byproducts. A total of 7,793 tailings ponds were present in the country by the end of 2017, with the highest tailings dam reaching 325 m (Yin et al., 2018). Because of the continuous exploitation of metal minerals and the limited availability of land resources, many tailings ponds have been expanded. The gradually increasing heights and slopes of tailings ponds have induced multiple failures of tailings dams throughout the country, leading to high loss of life, mineral resources, and state property (Li et al., 2015;Zhang et al., 2016). Thus, it is important to establish an effective numerical model that can simulate the unique flow of tailings pond slurry. Such a model would provide the means to quickly simulate tailings dam failure and immediately analyze the severity of failure incidents, and quantitatively divide at-risk areas downstream, thereby improving the safety management of downstream areas and reducing the severity of failure-induced damage.
The methods of numerically simulating tailings ponds can be divided into two types: The finite element methods and the discrete element methods. The finite element methods are based on continuum mechanics and are mainly used to simulate the continuous deformation and plastic failure of materials (Brzezinski, 2002;Zandarín et al., 2009). The failure analysis was carried out by establishing finite element model, and the damage of materials was evaluated numerically (Guo et al., 2020;Arafa et al., 2022;Tang et al., 2022). They can be further divided into the finite element method (Peter and Mahmood, 2003;Wang et al., 2011) and the finite difference method (Chakrabotry and Choudhury, 2009) and used to study the stability of tailings dam slopes and tailings ponds. The discrete element methods are based on discontinuous medium mechanics and are used to simulate the motion and collision characteristics of a bulk system. They can be further divided into two types, namely, block discrete element modelling and particle discrete element modelling. Presently, discrete element methods focus on the relationship between macro-and micromechanics (Fakhimi et al., 2002;Backstrom et al., 2008;Hsieh et al., 2008;Wang et al., 2014) and rock slope stability (Esaki et al., 1998;Xu et al., 2022), whereas limited research has been conducted on the stability of tailings ponds. The tailings pond discharge fluid is a mixture of tailings sand and water. The single finite element method does not easily realize the simulation of the movement process of microscopic particles in the process of tailings pond destruction. Although the single discrete element method can exhibit the dynamic process of damage, it does not easily reflect the whole process from the start of the dam failure to the end.
The continuum-based discrete element method (GDEM) software combines the continuum-based discrete element method (CDEM) and high-performance graphics processing unit (GPU) technology (Li et al., 2004;Zou et al., 2007;Zhang et al., 2019). CDEM couples the finite element method (FEM) with the discrete element method (DEM) to conduct finite element calculations inside block elements and discrete element calculations at block boundaries, which facilitates simulating the entire process from continuous deformation to fracture movement . The GDEM can be widely applied to progressive failure simulations and overall stability evaluations of materials in several engineering fields, including geotechnical, mining, tunnel, oil and gas, water conservancy, geological, and structural engineering.
To accurately analyze the tailings sediment movement, the inundation range, and the impact on sensitive targets downstream during the failure process of tailings dams, the present study applies GDEM numerical simulation for the first time to test the failure evolution of a tailings dam. The whole process of dam break of tailings pond was analyzed, thereby achieving a dynamic simulation of the entire failure process. In addition, the impact of tailings dam failure on sensitive targets downstream is determined by analyzing the movement of tailings slurry during the failure process, thereby providing theoretical and technical support for the safety management of tailings ponds, especially for downstream reservoir areas.
CDEM is a discrete element method based on continuum mechanics. By combining the finite element and discrete element methods, CDEM not only facilitates simulations of the continuous or discontinuous deformation and the movement of granular material under static and dynamic loading but also models the progressive failure of materials, from continuum to discontinuum (Nezami et al., 2004;Yeung et al., 2007).
The governing equation of the CDEM is the equation of particle motion (Li et al., 2007). In Eq. 1, m i is the mass of the ith node, and u i , _ u i , and € u i respectively represent the displacement, velocity, and acceleration vectors of the ith node. The corresponding calculation process is shown in Figure 1.
In CDEM, a discrete block can consist of one or more FEM elements (as shown in Figure 2). A continuum-based constitutive model and finite element approach are used inside the blocks, whereas a discontinuum-based constitutive model and discrete element approach are used for the block interfaces (Li et al., 2007;. The CDEM finite elements can be simple tetrahedrons, pentahedrons, and hexahedrons, or complex polyhedrons ( Figure 3). The discontinuous deformation is modeled by springs placed between the CDEM blocks, such that spring breakage simulates the cracking and slipping of material.
The CDEM uses the element stiffness matrix to obtain the deformation force generated by blocks. That is, the nodal force of the elements is obtained using the already known element stiffness matrix and the nodal displacement. Because the CDEM uses the dynamic relaxation method, the entire stiffness matrix need not be formed; rather, only obtaining the element stiffness matrix of each element is sufficient. At each iteration, the nodal force of the element itself is calculated using Eq. 2 and assigned to the corresponding node element. Among them, {F} e i is the nodal force vector of the ith element, [K] e i is the displacement vector of the ith element, and {u} e i is the element stiffness matrix of the ith element .
The normal and tangential springs of the CDEM contact surfaces are shown in Figure 4, where K n and K s represent the normal and tangential stiffness, respectively, C 1 and C 2 represent the normal and tangential damping coefficients, respectively.
The spring force can be calculated using Eq. 3, in which F j n and F j s represent the normal and tangential force of the jth spring, respectively, K j n and K j s represent the normal and tangential stiffness of the jth spring, respectively, and Δd j n and Δd j s represent the normal and tangential displacement of the jth spring, respectively (Li et al., 2007;.
When calculating element failure, the spring force is modified using the Mohr-Coulomb criterion and tensile criterion, as shown in Eq. 4.

FIGURE 2 | Blocks and interfaces in CDEM.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 901904 3 where, T is tensile strength, ϕ is internal friction angle, C is cohesion.

Principles of GDEM Method
The core of GDEM is continuous-discontinuous element method (CDEM) based on Lagrange equation, which organically combines finite element and discrete element. Block is used to characterize the continuum properties of materials, and interface is used to characterize the discontinuum properties of materials.
The gradual failure process of materials is simulated through the fracture of block boundary and interior. GDEM is used to simulate dam break, which can intuitively show the impact range of tailings dam break and effectively analyze the failure of tailings dam.

Descriptions on Geological Settings
The case tailings pond is surrounded by mountains on three sides. It is multiditched to form a claw-shaped valley. The surface of the mountain is a Quaternary slope, and some of the rocks are bare. The coverage by mountain vegetation is high, and the growth is general. There are Yingyun dioritic gneisses in the Shachang Formation of the Mishan Group in the Archean Period. Some parts of the red feldspar are more distinctive than others and have felsic veins.
The overall terrain of the site is relatively high and has little precipitation, the rock mass has no obvious water-collecting structure. The seismic fortification intensity is 7 degrees, and the basic acceleration of the designed earthquake is 0.15 g.

Working Conditions
The tailings pond has two initial dams, both of which are permeable rockfill dams. The elevation of the west dam is 149.3 m, and the elevation of the east dam is 143.5 m. The top width of the two dams is 5 m, the ratio of the inner and outer  The tailings accumulation dam is above the 163.5 m elevation of the dam crest at the beginning of the tailings pond. The dam slope has a total of 4 terraces at elevations of 173.5, 183.5, 193.5, and 203.5 m, and the mound width is 5 m, with a slope ratio of 1: 4.5, as shown in Figure 5.
Within 800 m downstream of the tailings pond, it is mainly farmland. There are railways and highways passing through approximately 1,000 m south of the tailings pond. There are three villages, machine repair workshops, mineral processing workshops, and comprehensive buildings downstream, with a total of 2,000 people.

Numerical Model Conditions
The tailings flow after dam failure is a free surface flow problem with a moving boundary and complex large deformation. The tailings flow after the failure is simulated using the smoothed particle hydrodynamics (SPH) method, which eliminates problems associated with entanglement and distortion; these are common to the traditional Eulerian grid method. As a result, the SPH method is uniquely advantageous in dealing with large-deformation free-surface and moving-boundary problems. It is highly suitable for simulating the tailings slurry flow after dam failure. To ensure the accuracy of the calculation results, the following assumptions were made about the model: 1) The tailings pond consists of countless discrete microparticles; 2) The tailings pond is made of a homogeneous material, and each part of the pond has the same physical properties; 3) The spatial model of the tailings pond is filled with continuous nonporous particles.
Assumed dam failure mode: the length of the breach is 350 m from the east-west dam junction, and the dam break depth to the bottom of the dam is 149.3 m at the bottom of the west dam and 143.5 m at the bottom of the east dam; the tailings pond water level is 207 m.
Assumed failure discharge: the total sediment discharge is taken as the total volume of water and tailings in the reservoir above 207 m; a total of 6 million cubic meters of tailings and water are discharged.
Flow resistance at the bottom of the total flow discharge: the change in the Manning coefficient is a function of the shear force; it simulates the flow characteristics of the tailings slurry, which are different from those of water; the relationship between the shear force and the Manning coefficient in the model is shown in Figure 6. Where, τ c is the critical shear force at which the tailings slurry begins to flow; τ cu is the critical shear force at which the tailings slurry changes from a yield-pseudoplastic fluid to a Bingham plastic fluid; M l is the minimum Manning coefficient, which characterizes the bottom friction when the tailings slurry begins to flow; M u is the maximum Manning coefficient, which characterizes the bottom friction when the tailings slurry changes from a yield-pseudoplastic fluid to a Bingham plastic fluid.

Modelling Procedure
During the analysis of the tailings dam failure, the area affected by the tailings pond is divided using unstructured triangular elements. The unstructured triangular mesh has the advantages of high adaptability for complex areas, flexible local encryption, and convenient self-adaptation. This method can effectively simulate natural boundaries and complex underwater topography and improves the accuracy of boundary simulations. Figure 7 shows the topographic map of the area and the designed stacked tailings dam generated in the CAD model.
The terrain surface and the dam model files are imported into the GDEM software. Next, the key surface points are stripped to form a NURBS surface. Finally, the surface model of the tailings dam is formed using the geometry-create-volume-by-contour command.
The terrain surface and the dam model are meshed. The terrain surface is divided into a 2D planar mesh, whereas the dam model is divided into a 3D particle mesh, ultimately generating 105,585 triangular terrain surface elements and 119,427 dam model particles.
The JavaScript script file is written from four aspects: as shown in Figure 8, setting the mechanical properties and material parameters of the particle model, calculating the critical parameters of the dam, and monitoring the data point, the solution step number, and the resulting output interval. The mechanical properties and material parameters of the case tailings library are shown in Table 1. The mechanical parameters of tailings were obtained by carrying out geotechnical tests on tailings, and were used as experimental parameters of this model.

Analysis of the Dam Failure Results
In Figure 8, Point one refers to the monitoring point at the foot of the dam, and Point two refers to the monitoring point of the comprehensive building. In the initial stage of the tailings pond failure, the mud gradually spreads from the dam site and then gradually diffuses into the low-lying area. Due to the high middle of the terrain and the low east and west sides, two mudflows gradually form downstream. After 1 min, the mud basically covers the entire calculation area. As the mud continues to flow to the lower area and out of the calculation area, the mud depth in the calculation area continuously decreases, and the submerged range shrinks, as shown in the flooding map at 70 s.
The failure process of a stacked dam can be divided into three stages: the initial stage of accelerated collapse, the intermediate stage of decelerated diffusion, and the final stage of stable accumulation. Steps 0 to 600 (0-6 s) are the initial stage of the accelerated collapse after the failure. During this stage, the tailings flow out of the tailings pond and diffuse rapidly westward, forming a fan-shaped collapse zone. Steps 600 to  During this stage, the tailings particles continue to pour out of the tailings pond; however, the particle velocity decreases gradually as their travel distance increases. At 5,500 steps, the velocity of the particles in the front of the flow tends to zero; however, the middle part of the tailings flow continues to  The calculation results of the submerged depth and flow rate of the case tailings pond reaching the downstream complex after the dam break are shown in Table 2. The flow velocity and quantity changes during the dam break process are shown in Figures 9, 10. The flow velocity and submerged depth during the tailings impact of the comprehensive building are shown in Figures 11, 12.
In Figures 9-12, the tailings instantaneous impact velocity after the dam breaks is 14.71 m/s with the maximum quantity of flow 51643 m 3 /s, and the maximum flow rate reaches 16.47 m/s in the 3rd second. As the volume of tailings in the reservoir continuously decrease, the velocity of the tailings flowing through the dam bottom gradually decrease. Finally, the tailings reach a steady-state in the 55th second, and the discharge flow rate is 0.97 m/s, which gradually approaches zero.
The tailings reach the downstream comprehensive building in the 16th second after the dam break, and the instantaneous flow rate is 24.68 m/s. Due to terrain effects, the impact tailing speed gradually decreases. After 80th second, the flow rate is gradually stabilized and finally reaches a steady state. The maximum overflow depth of the complex is 9.72 m, and the corresponding time is the 62 s after the occurrence of the collapse.

Comparative Analysis of the Calculation Results of Volume of Fluid
The Volume of Fluid (VOF) model tracks the volume of each fluid by solving a set of momentum equations and continuous   equations that simulate the motion of two or more fluids that are not blended. As shown in Figures 13, 14. The VOF adopts all the model breaks instantaneously to calculate and obtain the process line of the fracture flow. When the water level is 207 m and the upstream flow is not considered, the storage capacity is 6 million m 3 , and the time required is 150 s. The peak flow of the fracture is 270,000 m 3 /s, the maximum mud depth at the downstream comprehensive building is 1.45 m, and the corresponding time is 30 s after the occurrence of the collapse.
The comparison between the calculation results of the GDEM and the VOF method shows that the trend of the flow curve at the bottom is consistent, and there is a certain difference between the flow values. The change in the depth of the sediment that reaches the complex after the collapse is consistent, but the mud depth values are different.
By analyzing the calculation principle of the two methods, the two monitoring methods are shown to be different. The monitoring part of the GDEM adopts the arrangement monitoring ring method to monitor the discharge flow at the bottom of the dam body and in a certain area at the comprehensive building. The VOF calculates the discharge flow generated by the instantaneous discharge of all the broken sediments. The VOF method calculation of the discharge flow at the bottom of the tailings dam is higher than that of the GDEM. The VOF calculation does not consider the downstream topographic changes, and the downstream ground is set to the plane flow. The GDEM introduces the tailings pond and surrounding terrain into the original scale, which can truly reflect the change in the sediment flow after the tailings pond breaks, so the time of sediment flow in the VOF method is shorter than that in the GDEM. Comparison of the results of the two methods shows that the GDEM can effectively reflect the change in sediment flow after the tailings pond collapse and proves that the method can be used as an effective analysis method for tailings pond failure analysis.

CONCLUSION AND DISCUSSION
1) The process of tailings dam failure can be divided into three stages: accelerated collapse, decelerated diffusion and stable accumulation. During the initial stage of accelerated collapse, the tailings flow diffuses rapidly westward, forming a fan-   shaped collapse zone. During the intermediate stage of diffusion, the tailings particles continue to pour out from the tailings pond; however, the particle velocity gradually decreases and eventually stabilizes. During the final stage of stable accumulation, most of the tailings particles are discharged from the tailings pond, the tailings discharge flow is basically stable, and the distribution of the tailings flow does not have any significant changes; thus, the dam failure process ends. 2) The time of the dam-sludge mud-sand flow reaching the comprehensive building is 16 s, and the instantaneous impact velocity is 24.68 m 3 /s. The maximum submerged depth is 9.72 m at the 62 s. The final dam-sludge sand flow at the 83rd second tends to be stable at the comprehensive building. The experimental results show that once the tailings dam fails, it will have a serious impact downstream, seriously threatening the lives and property of the downstream residents. It is necessary to strengthen the safety management of the tailings pond to prevent the occurrence of dam failure. 3) Comparing the dam break calculation with the VOF method, the change in the discharge flow at the dam foundation position and the submerged depth at the downstream sensitive point after the dam break is consistent, which proves that the GDEM can be used for effective measurement in the tailings dam failure analysis. 4) A new large scale model is used to calculate the whole process of tailings dam failure. Using the GDEM to simulate dam failure can visually determine the impact range after tailings failure and then determine the dam break safety boundary. This method provides a quantitative approach for judging the stability of the tailings pond; in addition, it can effectively simulate the dam-break process and can be used as one of the methods to measure the stability of tailings pond failure. It also proves the feasibility of the GDEM in the evaluation of the stability of tailings ponds.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.