SPH modelling of dam breach run out flow for a site planning tailings storage facility

Tailings storage facilities (TSFs) are being built globally for containing the chief solid waste stream from mining industry. Catastrophic TSF breach accidents have occurred frequently since the beginning of the 21 century, causing severe impacts on the environment, economy and community safety. The recent example is the 2019 Brumadinho accident in Brazil that released 12 million m of tailings and killed more than 249 people. The foreknowledge of the TSF breach run out overland flow can be crucial to prevent or minimize possible losses. Using the Digital Surface Model (DSM) terrain data and the smoothed particle hydrodynamics (SPH) numerical method, this study proposed a procedure to predict the routings of hypothetical TSF breach run out flow over downstream complex terrain. A case study of a planning TSF site in Guizhou Province of China was carried out to evaluate its applicability. The results suggested that the maximum routing distance of the TSF breach run out flow was 1.45 km. At 240 s, the run out flow began to impact the downstream viaduct piers with the maximum submerged depth of 3.3 m and the maximum impact force of 21.8 kPa. Essential protective measures were recommended before the TSF site construction. The proposed procedure is then recommended for the safety management of the TSFs globally.


Introduction
A tailings storage facility (TSF) is an essential structure in mining industry built for the purpose of storing unwanted and currently uneconomical solid waste tailings behind one or more tailings dams. Unlike water-retaining dams, where construction materials usually consist of concrete, rock or soil, tailings dams are mostly being built using tailings themselves in order to minimize costs [1]. Therefore, the safety management are often more complicated, in terms of higher risks of dam breach, debris flow or overtopping accidents. The literatures [2,3] present that the rate of tailings dam breach (TDB) accident between 1910 and 2010 was estimated to be 1.2%, more than two orders of magnitude higher than that of normal water-retaining dams (0.01%). During the recent two decades, catastrophic TDB accidents have occurred frequently [4,5]. For example, on 25 January 2019, the Brumadinho dam accident at Córrego de Feijão iron ore mine in Minas Gerais, Brazil, released almost its complete holdings of 12 million m 3 of tailings. The mining loading station, buildings and railways were destroyed by the devastating run out flow, more than 249 people were killed [4]. On 5 November 2015, another catastrophic TDB accident occurred in Minas Gerais, Brazil [4]. The Fundão tailings dam of Samarco iron mine collapsed and released 32 million m 3 of tailings. The downstream community was flooded, at least 17 people lost their lives and more than 650 km of rivers were contaminated [4]. On 4 June 2018, the Cieneguita mine tailing dam in Chihuahua, Mexico, collapsed, released 439,000 m 3 of tailings and dam material that travelled 29 km downstream and at least 3 workers weredead [4].
The TDB accidents often initiate in very short time, release large volume of materials, transport in high velocities and travel long distances. Therefore, the foreknowledge of the run out overland flow routings is important to prevent or minimize possible losses. Especially for the site planning TSFs, the prediction of TDB run out flow routings can be helpful for the safety assessment and decision making. As on-site experiment is nearly impossible due to its uncontrollable and devastating process. Numerical method is considered to be an irreplaceable tool in modelling TDB problems, in terms of its high efficiency, flexible settings and reasonable costs.
However, traditional grid methods such as the Finite Difference Method (FDM) and the Finite Element Method (FEM) can lead to inevitable numerical difficulties including mesh winding, twisting and distortion in modelling problems with large deformation and free surfaces [6]. The fully mesh-free, particle-based, Lagrangian method of Smoothed Particle Hydrodynamics (SPH) proposed by Lucy [7], Gingold and Monaghan [8] is becoming more commonly-used in the latest geo-disaster studies [6,9]. For example, Huang, et al. [10] analyzed the flow-like landsides triggered by anearthquake in China with the SPH method and the results indicated good agreement with observed field characteristics. Vacondio, et al. [11] adopted SPH to simulate falling slide-generated flow in a water reservoir. The flow maximum run-up distances and surface elevations were reproduced satisfactorily. McDougall and Hungr [12] developed a SPH-adapted numerical model to analyze rapid landside motion across 3D terrain and designed a series of laboratory flume tests for verification. On the other hand, the run out flow properties such as velocities, submerged depths and directions can be directly affected by the complex terrains. However, to the best of the authors' knowledge, the downstream terrain was mostly simplified or ignored in the previous studies. Therefore, the aim of this study is to determine the potential for using the satellite DSM terrain data and the SPH numerical method to predict a hypothetical TDBrun out flow over downstream real terrain. And a case study of a site planning TSF in southwestern China was carried out to evaluate its applicability.

Study area
The study site of the Dalong TSF is situated inTongren city, Guizhou Province of southwestern China, as Figure 1 shows. The climate of study area is the warm and humid monsoon climate in the middle subtropics, with the annual precipitation of 1199.5 mm. According the records, the daily maximum precipitation was 226.2 mm. The site location is to chosen be surrounded by the mountain valley in order to minimizecosts. It was designed to store the manganese slag and mining waste from the local mines. And its capacity will augment by upstream raising the main dam with the slope ratio of 1/4. The capacity of the TSF is designed to be 2.15 million m 3 , with the main dam height of 30 m. The annual disposal of mining waste is estimated to be 340 thousand tons.

Fig. 1. The Dalong TSF site landform (left) and its location in China (right)
For the safety concerns, the downstream commuties within 600 m of the TSF will be relocated before the Dalong TSF is built. And approximately 600 m away of the TSFin the downstream area, several viaduct piers are being built, as Figure 2 shows. Therefore, the main concerns of this case study will be paid on the safety of the viaductpiers. The study area was defined to be within 3 km of the Dalong TSF.

Smoothed Particle Hydrodynamics (SPH) method
The applicability of using SPH method to model tailings dam breach routings wasproved by Wang, et al. [13] . In SPH, the continuum is represented by defining a set of discrete particles, rather than by generating mesh in traditional grid methods. Each particle is initially characterized by its specific physical properties, such as mass, density, pressure and energy. And then a particle is influenced by its neighboring particles by using the kernel approximation. The interpolation function designated as the smoothing kernel function can be written in a continuous form as: (1) where h is the smoothing length which indicates the interaction distance between two neighboring particles, r represents the particle location vector, Ω is the total supporting domain determined by h, and W(r-r', h) is the Dirac delta weighting function which is E3S Web of Conferences 303, Clean Coal Technologies: Mining, Processing, Safety, and Ecology 2021 fundamental to the SPH approach.
Equation (1) can be expressed in a non-continuous form after a discrete approximation step as: (2) where N is the amount of neighboring particles, m is the mass and ρ is the density.The smoothing function W(ra-rb, h) is correlated with smoothing length h and (ra-rb) i.e. the distance between particles a and b. In this study, the fifth-order quantic kernel by Wendland [14] where the weighting function vanishes for inter-particle distances greater than 2h is used. The Wendland kernel is defined in 3D as: where q=(ra-rb)/h and αD is the normalization constant that equals to 7/(4πh 2 ) in 2D and 21/(16πh 3 ) in 3D.

State Equation
The weakly compressible state equation proposed by Monaghan [15] is adopted in this study. The equation which relates the hydrostatic pressure to local densities is as follows: where constant B sets a limit for the range of density values. B equals to 200(ρ0)gH/γ when the liquid level is H. Constant γ is usually taken as 7, ρ0 is the reference density.

Governing Equations
The momentum equation in the Lagrange coordinate system can be written as: (5) where P represents the pressure, g= (0, 0, -9.81) m/s 2 is the gravity acceleration, and ψ is the dissipative term.
In most reported tailings dam breach cases, the slurry concentration is normally less than 50% and described as a fluid of hydraulic behavior little different from water [16,17]. Therefore the artificial viscosity scheme proposed by Monaghan [18] is adopted in this study due primarily to its simplicity and applicability in free surface water simulation [19][20][21][22]. The scheme is described as: E3S Web of Conferences 303, Clean Coal Technologies: Mining, Processing, Safety, and Ecology 2021 (6) v being the velocity vector and the viscous term according to the artificial viscosity. The artificial viscosity can be calculated as: (7) where α is the coefficient with the main purpose of preventing instability and spurious oscillations in the numerical scheme, c is the sound speed and / The fluid density fluctuation is computed by solving the mass conservation as theparticle mass keeps constant in the Weakly Compressible SPH (WCSPH) computing. The continuity equation in SPH is as follows: (8) Particle position is updated using the following XSPH momentum equation [23]: (9) ε being a constant value in the range of 0~1 and To account for the WCSPH, the Delta-SPH equation proposed by Molteni and Colagrossi [24] was implemented. A diffusive term is introduced to reduce the density fluctuations. The Delta-SPH equation can be written in the following form: (10) where is the diffusive coefficient of Delta-SPH.

Numerical modelling execution
The open-source SPH solver DualSPHsysics v4.0 [25] (www.dual.sphysics.org) was used in the present study. A series of .xml (EXtensible Markup Language), .stl (STereoLithography) and .sh (SHell script) files were compiled to define the material parameters, system configuration and modelling execution. Based on the ©JAXA ALOS global Digital Surface Model (DSM) data, the real-scale 3D geometry of the research area was constructed after a series of transitions. Firstly, the DSM raster data were imported into the open-source QGIS software (https://qgis.org) and converted to the .stl (STereoLithography) format 3D geometry. Thereafter, according to the previous studies, the dam breach widths in most earthen dam breach cases are 3 times of the breach depths [26], and the ratio of the top breach width to the bottom breach width ranges from 1.06 to 1.74, with an average value of 1.29 [27]. The TDB mode was simplified to be E3S Web of Conferences 303, Clean Coal Technologies: Mining, Processing, Safety, and Ecology 2021 under the extreme conditions of one-in-two-centuries flood and major earthquake, that the breach initiates instantaneously at the first step of the simulation. On the basis of available computational capacities, the smoothing length was set to 2 m which eventually generated 821,009boundary particlesand 48,565fluid particles. Thereafter, with the rapid decrease the TSF capacity, the velocities of run out flow significantly reduced at 240 s, 600 s and 1200 s. The flow velocity peak value had reduced to less than 10 m/s at 240 s and began to impact the downstream commuties and viaduct piers. Thereafter at 600 s and 1200 s, the flow tended to be stable and wasspreading at less than 6 m/s through the downstream valley. The maximum routing distance of the TSF breach run out flow was calculated to be 1.45 km. According to the measured data at the viaduct piers that 600 m away from the planning TSF site, the maximum submerged depth was 3.3 m and the maximum impact force was 21.8 kPa. The safety of the viaduct piers should be further considered. Essential protective measures such as diversion or cut-off ditches were recommended.
In addition, the maximum submergence depth of the run out flow at 1 km downstream of the TSF was 3.5 m, and the maximum impact force was 15.2 kPa. Lastly, at the farthest submergence range that about 1.45km downstream of the TSF, the maximum submerged depth was 2.6m and the maximum impact force was 9.8 kPa.

Conclusions
Over the past decade, catastrophic tailings dam breach accidents highlighted the importance of TSF risk assessment and emergency management. As contemporary mining industry is exploiting more low-grade ores and will discharge more tailings [28],the capacity and amount of TSF will continue to grow in the foreseeable future.
This study presents a novel procedure of using satellite DSM terrain data and SPH numerical method to determine the run out overland flow of a TSF breach. A case studyof a hypothetical TSF breach in Dalong TSF, Guizhou Province, China was then carried out. Under the worst conditions of one-in-two-centuries flood and major earthquake, the run out flow peak velocity of 18 m/s was identified at its dam toe. At 600 s the run out flow began to impact the viaduct piers that 600 m away from the planning TSF site,the maximum submerged depth was measured to be 3.3 m and the maximum impact force was 21.8 kPa. The safety of the viaduct piers and essential protective measures were recommended to be further considered. In addition, the maximum submerged depth and impact force at different interest area can be directly obtained from the numerical results.
The proposed procedure can be helpful for the safety assessment and decision making of the TSFs construction and management, especially for the site planning TSFs. The foreknowledge of the key features such as the peak velocity, submerged area, depths and impact force can be crucial for the stakeholders to mitigate possible consequences and ensure that catastrophic accidents will not occur.
Funding: This research was funded by the National Natural Science Foundation of China (grant number 51674150 and 51974051).