Integration of DSM and SPH to Model Tailings Dam Failure Run-Out Slurry Routing Across 3 D Real Terrain

Tailings dam failure accidents occur frequently, causing substantial damage and loss of human and animal life. The prediction of run-out tailings slurry routing following dam failures is of great significance for disaster prevention and mitigation. Using satellite remote sensing digital surface model (DSM) data, tailings pond parameters and the advanced meshless smoothed particle hydrodynamics (SPH) method, a 3D real-scale numerical modelling method was adopted to study the run-out tailings slurry routing across real downstream terrains that have and have not been affected by dam failures. Three case studies, including a physical modelling experiment, the 2015 Brazil Fundão tailings dam failure accident and an operating high-risk tailings pond in China, were carried out. The physical modelling experiment and the known consequences were successfully modeled and validated using the SPH method. This and the other experiments showed that the run-out tailings slurry would be tremendously destructive in the early stages of dam failure, and emergency response time would be extremely short if the dam collapses at its full designed capacity. The results could provide evidence for disaster prevention and mitigation engineering, emergency management plan optimization, and the development of more responsible site plans and sustainable site designs. However, improvements such as rheological model selection, terrain data quality, computing efficiency and land surface roughness need to be made for future studies. SPH numerical modelling is a powerful and advanced technique that is recommended for hazard assessment and the sustainable design of tailings dam facilities globally.


Introduction
A tailings pond is a structure made up of one or more dams, built for the purpose of storing unwanted and currently uneconomical mine waste tailings.Unlike water-retaining dams, where the construction materials usually consist of concrete, rock or soil, tailings dams are mostly designed to be built using tailings themselves, in order to minimize costs [1].However, the proportion of tailings dam failure accidents between 1910 and 2010 has been estimated to be 1.2%, more than two orders of magnitude higher than water-retaining dams (0.01%) [2,3].Tailings dam failure accidents, with their limited emergency response times, large volumes of released materials, high velocity slurry flows and long run-out distances, can often lead to heavy casualties and irrecoverable losses [4].For example, on 4 June 2018, the Cieneguita mine tailings dam in Chihuahua, Mexico, collapsed, releasing 439,000 m 3 of tailings and dam material that travelled 29 km downstream and killed at least three workers [5].On 12 March 2017, the Tonglvshan mine tailings dam failure in Hubei province, China, released approximately 200,000 m 3 of tailings, flooding approximately 27 hectares of fish ponds and killing at least two people [5].On 4 August 2014, the Mount Polley mine tailing dam failure in British Columbia, Canada, released 24.4 Mm 3 of tailings and water that flowed into adjacent Polley Lake, through Hazeltine Creek, and eventually entered Quesnel Lake approximately 10 km away [5].
As contemporary mining is focusing on the extraction of large volumes of low grade ore, larger amounts of mine waste are being produced than previously [6].For example in China, one of the world's largest developing economies with vast mineral resource exploitation and consumption, more than 1.5 billion tons of tailings were produced in 2016, and this figure is still increasing [7].Despite significant efforts by the Chinese government to improve the reuse and recycling of tailings, as many as 8869 tailings ponds were documented in 2016.Among these, 1425 tailings ponds were designated 'overhead tailings ponds' (OTPs), which are high-risk tailing ponds located upstream and within only 1 km of residential areas, industrial plants, schools, hospitals and other important facilities [8].Furthermore, these dams were mostly constructed by the economical and practical upstream method, which has been statistically proven to be less reliable than other methods [2].Therefore, the government and stakeholders are under tremendous pressure to implement measures to avoid or minimize casualties, environmental impacts, financial losses and social problems that would be caused by unpredictable tailing dam failures.Most previous studies have focused on the environmental impacts and dam geotechnical properties [4,[8][9][10][11].Hence, the aim of this study is to fill this knowledge gap by using numerical modelling to determine the impact area of the tailings slurry run-out and provide evidence to initiate more responsible and sustainable site plans and designs, as well as improving tailings dam failure emergency response plans.
Numerical modelling is playing an increasingly important role in solving complex and devastating geo-disaster problems, such as landslides, debris flows, dam failures and soil liquefaction [12][13][14], due to its high efficiency, flexible setting and low cost.However, tailings run-out routing is considered to be difficult to analyze and predict due to the large-scale deformation caused and the complex downstream terrains [15,16].Grid methods such as the finite difference method (FDM) and the finite element method (FEM) with mesh can lead to inevitable numerical difficulties including mesh tangling, twisting and distortion in solving problems with large deformation and free surfaces [17].Smoothed particle hydrodynamics (SPH) is a fully mesh-free, particle-based, pure Lagrangian method that is considered to be an approach circumventing mesh tangling problem [18].It has been decades since it was first proposed to solve astrophysical problems [19,20].It is a relatively new and advanced computational fluid dynamics (CFD) method with industrial applications, such as in the aerospace industries, car industries, energy production industries, marine industries and mining industries [18].Recently, SPH has been successfully applied to the geo-disaster field [17,18].For example, Huang et al. [21] analyzed run-out, flow-like landsides triggered by an earthquake in China using the SPH method and the results showed good agreement with the observed field characteristics.Vacondio et al. [22] adopted SPH to simulate falling, slide-generated flow in a water reservoir, and the flow maximum run-up distances and surface elevations were reproduced satisfactorily.McDougall and Hungr [23] developed a SPH-adapted numerical model to analyze rapid landside motion across 3D terrain and carried out a series of laboratory flume tests for verification.
The open-source DualSPHysics code is implemented in C++ and the compute unified device architecture (CUDA) language to allow either central processing unit (CPU) or graphics processing unit (GPU) modelling, so that the SPH method can be applied to industrial cases with large domains and long computing runtimes [24].GPU modelling has been shown to achieve speeds of up to two Water 2018, 10, 1087 3 of 15 orders of magnitude higher than those of single-core CPUs, thus achieving the accuracy of the CPU modelling and the efficiency of the GPU computing [25].Therefore, by using the DualSPHysics code, the large-scale 3D modelling of the run-out slurry routing caused by tailings dam failure over real terrain becomes possible with reasonable costs and high efficiency.

DSM Data
The digital surface model (DSM) data used in this study were taken from the freely available "Advanced Land Observing Satellite (ALOS) World 3D 30 m mesh" (AW3D30) dataset provided by Japan Aerospace Exploration Agency (JAXA), with 1 arc second (approximately 30 m) resolution [26].The dataset was generated using the data acquired from 2006 to 2011 by the ALOS.The tailings pond 3D geometry was built based on the dam height, tailings volume and nearby terrain, and was integrated with the AW3D30 DSM data.A 3D-geometry model including the tailings pond and nearby real terrain was then generated for SPH modelling.

SPH Method
In the SPH method, grids are completely abandoned and the continuum is represented as a set of discrete particles characterized by their own physical properties such as mass, density, pressure and energy.A particle is determined by its neighboring particles in the supporting domain, instead of by a grid relationship, as in traditional mesh methods.This step is usually called the kernel approximation.The interpolation function is designated the smoothing kernel function, which can be written in a continuous form as: where h is the smoothing length that 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, the origin of the SPH approach.Equation (1) can be expressed in a non-continuous form after a discrete approximation step as: where N is the amount of neighboring particles, m is the mass and ρ is the density.The smoothing function W(r a − r b , h) is correlated with the smoothing length h and (r a − r b ) i.e., the distance between particles a and b.This study used the fifth-order quantic kernel by Wendland [27], where the weighting function vanishes for inter-particle distances greater than 2h.The Wendland kernel is defined in 3D as: where q = (r a − r b )/h and α D is the normalization constant that equals 7/(4πh 2 ) in 2D and 21/(16πh 3 ) in 3D.

State Equation
The weakly compressible state equation proposed by Monaghan [28] was adopted for this study.The equation that relates the hydrostatic pressure to the local density is as follows: where the constant B sets a limit for the range of density values.B equals 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: where P represents the pressure, g = (0, 0, −9.81) m/s 2 is the gravity acceleration, and Ψ is the dissipative term.This study was carried out under the assumption of maximum capacity, and that dam failure would happen instantaneously and release a huge volume of slurry.The slurry is usually reported to have a concentration of less than 50% and is described as a fluid with a hydraulic behavior that is little different from that of water [29,30].Therefore, the artificial viscosity scheme proposed by Monaghan [31] was adopted in this study, due primarily to its simplicity and applicability in free surface water simulation [32][33][34][35].The scheme is described as: where v is the velocity vector and Π ab is the viscous term according to the artificial viscosity.
The artificial viscosity can be calculated as: where r ab = r a − r b , v ab = v a − v b , ρ ab = (ρ a + ρ b )/2, µ ab = hv ab • r ab /(r 2 ab + η 2 ), η 2 = 0.01h 2 , α is the coefficient with the main purpose of preventing instability and spurious oscillations in the numerical scheme, c is the sound speed and c ab = (c a + c b )/2.
The fluid density fluctuation is computed by solving the mass conservation, as the particle mass is kept constant through weakly compressible SPH (WCSPH) computing.The continuity equation in SPH is as follows: The particle position is updated using the following XSPH momentum equation [36]: being constant ε value ranges 0~1 and ρ ab = (ρ a + ρ b )/2.
A delta-SPH formulation [37] is also available from DualSPHysics, which introduces a diffusive term in the continuity that is then applied to reduce density fluctuations.This gives more reliable results for the WCSPH.The equation can be written as: where δ is the delta-SPH coefficient.

DualSPHysics GPU Execution
According to the DualSPHysics code [24], the input parameters of the system configuration and its execution are defined by a series of .xml(EXtensible Markup Language), .stl(STereoLithography) and .sh(SHell script) files.The code was executed using the NVIDIA Tesla K80 GPU card (4992 processor cores with a dual-GPU design, 560 MHz, NVIDIA, Santa Clara, CA, USA).

Research Background
A laboratory physical modelling experiment was carried out by Jing et al. [30] to study the consequences of the potential failure of the Yangtianqing tailings pond located in Yuxi, China (Figure 1).This was an operating tailings pond that was constructed using the upstream method in the nearby valley, with a designed maximum capacity of 108.9 Mm 3 .As Figure 1b shows, Mimao village is located only about 0.8 km downstream of this typical 'overhead tailings pond' (OTP), occupying both the northern and southern banks of the valley's sharp bend.A laboratory study [30] was adopted here as a reference to validate the SPH method of computing tailings slurry routing.
Water 2018, 10, x FOR PEER REVIEW 5 of 15

DualSPHysics GPU Execution
According to the DualSPHysics code [24], the input parameters of the system configuration and its execution are defined by a series of .xml(EXtensible Markup Language), .stl(STereoLithography) and .sh(SHell script) files.The code was executed using the NVIDIA Tesla K80 GPU card (4992 processor cores with a dual-GPU design, 560 MHz, NVIDIA, Santa Clara, CA, USA).

Research Background
A laboratory physical modelling experiment was carried out by Jing et al. [30] to study the consequences of the potential failure of the Yangtianqing tailings pond located in Yuxi, Yunnan Province, China (Figure 1).This was an operating tailings pond that was constructed using the upstream method in the nearby valley, with a designed maximum capacity of 108.9 Mm 3 .As Figure 1b shows, Mimao village is located only about 0.8 km downstream of this typical 'overhead tailings pond' (OTP), occupying both the northern and southern banks of the valley's sharp bend.A laboratory study [30] was adopted here as a reference to validate the SPH method of computing tailings slurry routing.

Model Setup
A simplified model of the tailings pond and downstream terrain was constructed at 1:400 scale of the experimental laboratory flume.The volumetric concentration of tailings was set to 40%, according to in situ tailings disposal parameters and debris mudflow accidents [30].As shown in Figure 2, tailings slurry was mixed and pumped into a 3000 mm × 3048 mm storage tank that represented the tailings pond.The run-out slurry was designed to flow towards the downstream flume shown in Figure 2b through a sharp bend after the tank gate was opened.A piezometer and data sampling system were set up at the sharp bend to capture the impact pressure at various time steps.Slurry wave shape and velocity data were collected by a series of digital cameras and rulers.

Model Setup
A simplified model of the tailings pond and downstream terrain was constructed at 1:400 scale of the experimental laboratory flume.The volumetric concentration of tailings was set to 40%, according to in situ tailings disposal parameters and debris mudflow accidents [30].As shown in Figure 2, tailings slurry was mixed and pumped into a 3000 mm × 3048 mm storage tank that represented the tailings pond.The run-out slurry was designed to flow towards the downstream flume shown in Figure 2b through a sharp bend after the tank gate was opened.A piezometer and data sampling system were set up at the sharp bend to capture the impact pressure at various time steps.Slurry wave shape and velocity data were collected by a series of digital cameras and rulers.
The rebuilt 3D geometry and fluid properties for SPH modelling were set to correspond to those of the flume experiment.The smoothing length was set to 0.015 m according to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA).Eventually, 315,983 boundary particles and 894,768 fluid particles were generated.The rebuilt 3D geometry and fluid properties for SPH modelling were set to correspond to those of the flume experiment.The smoothing length was set to 0.015 m according to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA).Eventually, 315,983 boundary particles and 894,768 fluid particles were generated.

Results
According to the SPH numerical results shown in Figure 3, the slurry flow reached x = 2 m, x = 6.5 m, x = 9 m and x = 15.35 m at t = 1.6 s, t = 3.8 s, t = 4.8 s and t = 7.4 s, respectively.The slurry flow velocity distribution law was in good agreement with the phenomenon described in the literature [30].For example, the front flow velocity was higher and a decrease was observed along the flow's long tail.The maximum velocity of the slurry flow at x = 6.5 m was 2.41 m/s, with only 3.6% difference to the experimental result of 2.5 m/s.Flow pattern variations were observed at the sharp bend, and these are ascribed mainly to the centrifugal force and slurry reflection wave.The velocity was much lower at the outer side of the sharp bend, likely because the centrifugal force caused high-velocity slurry to run up to the outer side, slowing down as its kinetic energy was transformed into gravitational potential energy.When the slurry reached a steady flow stage at t = 18 s, a clear vortex was observed in both the experiment snapshot shown in Figure 4a and the numerical particle trajectories shown in Figure 4b This vortex was likely caused by the interaction between the northeastern high velocity reflection waves and the southwestern low velocity flow.

Results
According to the SPH numerical results shown in Figure 3, the slurry flow reached x = 2 m, x = 6.5 m, x = 9 m and x = 15.35 m at t = 1.6 s, t = 3.8 s, t = 4.8 s and t = 7.4 s, respectively.The slurry flow velocity distribution law was in good agreement with the phenomenon described in the literature [30].For example, the front flow velocity was higher and a decrease was observed along the flow's long tail.The maximum velocity of the slurry flow at x = 6.5 m was 2.41 m/s, with only 3.6% difference to the experimental result of 2.5 m/s.Flow pattern variations were observed at the sharp bend, and these are ascribed mainly to the centrifugal force and slurry reflection wave.The velocity was much lower at the outer side of the sharp bend, likely because the centrifugal force caused high-velocity slurry to run up to the outer side, slowing down as its kinetic energy was transformed into gravitational potential energy.When the slurry reached a steady flow stage at t = 18 s, a clear vortex was observed in both the experiment snapshot shown in Figure 4a and the numerical particle trajectories shown in Figure 4b This vortex was likely caused by the interaction between the northeastern high velocity reflection waves and the southwestern low velocity flow.The rebuilt 3D geometry and fluid properties for SPH modelling were set to correspond to those of the flume experiment.The smoothing length was set to 0.015 m according to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA).Eventually, 315,983 boundary particles and 894,768 fluid particles were generated.

Results
According to the SPH numerical results shown in Figure 3, the slurry flow reached x = 2 m, x = 6.5 m, x = 9 m and x = 15.35 m at t = 1.6 s, t = 3.8 s, t = 4.8 s and t = 7.4 s, respectively.The slurry flow velocity distribution law was in good agreement with the phenomenon described in the literature [30].For example, the front flow velocity was higher and a decrease was observed along the flow's long tail.The maximum velocity of the slurry flow at x = 6.5 m was 2.41 m/s, with only 3.6% difference to the experimental result of 2.5 m/s.Flow pattern variations were observed at the sharp bend, and these are ascribed mainly to the centrifugal force and slurry reflection wave.The velocity was much lower at the outer side of the sharp bend, likely because the centrifugal force caused high-velocity slurry to run up to the outer side, slowing down as its kinetic energy was transformed into gravitational potential energy.When the slurry reached a steady flow stage at t = 18 s, a clear vortex was observed in both the experiment snapshot shown in Figure 4a and the numerical particle trajectories shown in Figure 4b This vortex was likely caused by the interaction between the northeastern high velocity reflection waves and the southwestern low velocity flow.Figures 5 and 6 show the comparisons between the numerical and experimental results of the slurry submerged depth (at x = 6.5 m and x = 9 m) and impact pressure (at the sharp bend) curves versus time.The depth versus time curves generally agree with the experimental results, but the calculated depth peak time is slightly delayed.Nevertheless, the calculated depth peak values of 10.80 cm at x = 6.5 m and 8.19 cm at x = 9 m generally agree with the experimental results of 10.88 cm at x = 6.5 m and 7.80 cm at x = 9 m, respectively.The experimental impact pressure was measured by preset piezometer.Due to limitations in the measurement resolution and sensitivity, only pressure data from t = 0 s to t = 21 s were collected.The calculated curve peak between 0 s and 5 s was much sharper than that of the experimental result, but the curve variation trends were generally similar.The numerical pressure peak values of 21.67 kPa at t = 1.8 s and 22.96 kPa at t = 2 s were slightly higher than, and delayed compared to, the experimental results (19.05 kPa at t = 1.32 s and 18.92 kPa at t = 1.76 s).Both the numerical and experimental pressure peak values were approximately 3.5 times greater than the pressure at the stable flow stage after t = 7.5 s.Figures 5 and 6 show the comparisons between the numerical and experimental results of the slurry submerged depth (at x = 6.5 m and x = 9 m) and impact pressure (at the sharp bend) curves versus time.The depth versus time curves generally agree with the experimental results, but the calculated depth peak time is slightly delayed.Nevertheless, the calculated depth peak values of 10.80 cm at x = 6.5 m and 8.19 cm at x = 9 m generally agree with the experimental results of 10.88 cm at x = 6.5 m and 7.80 cm at x = 9 m, respectively.The experimental impact pressure was measured by preset piezometer.Due to limitations in the measurement resolution and sensitivity, only pressure data from t = 0 s to t = 21 s were collected.The calculated curve peak between 0 s and 5 s was much sharper than that of the experimental result, but the curve variation trends were generally similar.The numerical pressure peak values of 21.67 kPa at t = 1.8 s and 22.96 kPa at t = 2 s were slightly higher than, and delayed compared to, the experimental results (19.05 kPa at t = 1.32 s and 18.92 kPa at t = 1.76 s).Both the numerical and experimental pressure peak values were approximately 3.5 times greater than the pressure at the stable flow stage after t = 7.5 s.Figures 5 and 6 show the comparisons between the numerical and experimental results of the slurry submerged depth (at x = 6.5 m and x = 9 m) and impact pressure (at the sharp bend) curves versus time.The depth versus time curves generally agree with the experimental results, but the calculated depth peak time is slightly delayed.Nevertheless, the calculated depth peak values of 10.80 cm at x = 6.5 m and 8.19 cm at x = 9 m generally agree with the experimental results of 10.88 cm at x = 6.5 m and 7.80 cm at x = 9 m, respectively.The experimental impact pressure was measured by preset piezometer.Due to limitations in the measurement resolution and sensitivity, only pressure data from t = 0 s to t = 21 s were collected.The calculated curve peak between 0 s and 5 s was much sharper than that of the experimental result, but the curve variation trends were generally similar.The numerical pressure peak values of 21.67 kPa at t = 1.8 s and 22.96 kPa at t = 2 s were slightly higher than, and delayed compared to, the experimental results (19.05 kPa at t = 1.32 s and 18.92 kPa at t = 1.76 s).Both the numerical and experimental pressure peak values were approximately 3.5 times greater than the pressure at the stable flow stage after t = 7.5 s.Figures 5 and 6 show the comparisons between the numerical and experimental results of the slurry submerged depth (at x = 6.5 m and x = 9 m) and impact pressure (at the sharp bend) curves versus time.The depth versus time curves generally agree with the experimental results, but the calculated depth peak time is slightly delayed.Nevertheless, the calculated depth peak values of 10.80 cm at x = 6.5 m and 8.19 cm at x = 9 m generally agree with the experimental results of 10.88 cm at x = 6.5 m and 7.80 cm at x = 9 m, respectively.The experimental impact pressure was measured by preset piezometer.Due to limitations in the measurement resolution and sensitivity, only pressure data from t = 0 s to t = 21 s were collected.The calculated curve peak between 0 s and 5 s was much sharper than that of the experimental result, but the curve variation trends were generally similar.The numerical pressure peak values of 21.67 kPa at t = 1.8 s and 22.96 kPa at t = 2 s were slightly higher than, and delayed compared to, the experimental results (19.05 kPa at t = 1.32 s and 18.92 kPa at t = 1.76 s).Both the numerical and experimental pressure peak values were approximately 3.5 times greater than the pressure at the stable flow stage after t = 7.5 s.

Modelling the 2015 Brazilian Fundão Tailings Dam Failure Accident
On 5 November 2015, the Fundão tailings dam of Samarco Mine, Minas Gerais, Brazil, collapsed and released approximately 32 Mm 3 of tailings.The downstream community of Bento Rodriguez, located 5 km away was submerged, resulting in at least 17 people losing their lives and more than 650 km river being contaminated [5].Based on the JAXA AW3D30 DSM dataset, the released tailings volume of 32 Mm 3 , the dam crest elevation before the collapse of 900 m and the dam shape [29], a 3D geometrical model of the research area (Figure 7) was built.This case study was carried out in order to determine the applicability of the SPH method to a real, very large-scale tailings dam failure.According to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA), the smoothing length was set to 3 m, which eventually generated 18,132,290 boundary particles and 2,987,759 fluid particles.Some typical numerical result snapshots are shown in Figure 8.At t = 300 s, the peak value of the run-out slurry velocity was more than 20 m/s and appeared at about 1 km downstream of the tailings dam.This is due to the large terrain undulation at this position and the conversion of the gravitational potential energy of the slurry into kinetic energy, causing the velocity to increase sharply.At t = 300 s, the slurry front flow reached the downstream Santarem dam 3 km away.According to the panel investigation report [29], the slurry accumulated and overtopped, but did not breach the Santarem dam, then entered Bento Rodriguez before being routed to its ultimate destination of the Atlantic Ocean.At t = 600 s, the slurry overtopped Santarem dam and its velocity reduced considerably before the dam.Soon after, the slurry velocity increased to about 16 m/s due to elevation differences.At t = 800 s, the slurry front flow approached Bento Rodriguez, 5 km downstream of the Fundão dam.At this time, the flow become stable, moving slowly downstream

Modelling the 2015 Brazilian Fundão Tailings Dam Failure Accident
On 5 November 2015, the Fundão tailings dam of Samarco Mine, Minas Gerais, Brazil, collapsed and released approximately Mm 3 of tailings.The downstream community of Bento Rodriguez, located 5 km away was submerged, resulting in at least 17 people losing their lives and more than 650 km river being contaminated [5].Based on the JAXA AW3D30 DSM dataset, the released tailings volume of 32 Mm 3 , the dam crest elevation before the collapse of 900 m and the dam shape [29], a 3D geometrical model of the research area (Figure 7) was built.This case study was carried out in order to determine the applicability of the SPH method to a real, very large-scale tailings dam failure.According to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA), the smoothing length was set to 3 m, which eventually generated 18,132,290 boundary particles and 2,987,759 fluid particles.

Modelling the 2015 Brazilian Fundão Tailings Dam Failure Accident
On 5 November 2015, the Fundão tailings dam of Samarco Mine, Minas Gerais, Brazil, collapsed and released approximately 32 Mm 3 of tailings.The downstream community of Bento Rodriguez, located 5 km away was submerged, resulting in at least 17 people losing their lives and more than 650 km river being contaminated [5].Based on the JAXA AW3D30 DSM dataset, the released tailings volume of 32 Mm 3 , the dam crest elevation before the collapse of 900 m and the dam shape [29], a 3D geometrical model of the research area (Figure 7) was built.This case study was carried out in order to determine the applicability of the SPH method to a real, very large-scale tailings dam failure.According to the memory size and computing efficiency of the NVIDIA Tesla K80 GPU card (NVIDIA, CA, USA), the smoothing length was set to 3 m, which eventually generated 18,132,290 boundary particles and 2,987,759 fluid particles.Some typical numerical result snapshots are shown in Figure 8.At t = 300 s, the peak value of the run-out slurry velocity was more than 20 m/s and appeared at about 1 km downstream of the tailings dam.This is due to the large terrain undulation at this position and the conversion of the gravitational potential energy of the slurry into kinetic energy, causing the velocity to increase sharply.At t = 300 s, the slurry front flow reached the downstream Santarem dam 3 km away.According to the panel investigation report [29], the slurry accumulated and overtopped, but did not breach the Santarem dam, then entered Bento Rodriguez before being routed to its ultimate destination of the Atlantic Ocean.At t = 600 s, the slurry overtopped Santarem dam and its velocity reduced considerably before the dam.Soon after, the slurry velocity increased to about 16 m/s due to elevation differences.At t = 800 s, the slurry front flow approached Bento Rodriguez, 5 km downstream of the Fundão dam.At this time, the flow become stable, moving slowly downstream Some typical numerical result snapshots are shown in Figure 8.At t = 300 s, the peak value of the run-out slurry velocity was more than 20 m/s and appeared at about 1 km downstream of the tailings dam.This is due to the large terrain undulation at this position and the conversion of the gravitational potential energy of the slurry into kinetic energy, causing the velocity to increase sharply.At t = 300 s, the slurry front flow reached the downstream Santarem dam 3 km away.According to the panel investigation report [29], the slurry accumulated and overtopped, but did not breach the Santarem dam, then entered Bento Rodriguez before being routed to its ultimate destination of the Atlantic Ocean.At t = 600 s, the slurry overtopped Santarem dam and its velocity reduced considerably Water 2018, 10, 1087 9 of 15 before the dam.Soon after, the slurry velocity increased to about 16 m/s due to elevation differences.At t = 800 s, the slurry front flow approached Bento Rodriguez, 5 km downstream of the Fundão dam.At this time, the flow become stable, moving slowly downstream with an average velocity of about 8 m/s.At t = 1800 s, Bento Rodrigues was mostly submerged by the run-out slurry.At this time, the slurry flow was divided into two branches in the 'T-shaped' valley on the south-western side of the community.One branch flowed south into a valley, and the other flowed in a south-easterly direction over the village and into the Gualaxo do Norte River.As the terrain became flatter and the flow more stable, the slurry velocity decreased to less than 5 m/s.In general, the SPH numerical results of the slurry flow directions and impact area were in good agreement with the actual results, shown in Figure 7, which was drawn based on the satellite image.
Water 2018, 10, x FOR PEER REVIEW 9 of 15 with an average velocity of about 8 m/s.At t = 1800 s, Bento Rodrigues was mostly submerged by the run-out slurry.At this time, the slurry flow was divided into two branches in the 'T-shaped' valley on the south-western side of the community.One branch flowed south into a valley, and the other flowed in a south-easterly direction over the village and into the Gualaxo do Norte River.As the terrain became flatter and the flow more stable, the slurry velocity decreased to less than 5 m/s.In general, the SPH numerical results of the slurry flow directions and impact area were in good agreement with the actual results, shown in Figure 7, which was drawn based on the satellite image.The curves of the slurry submerged depth, velocity and impact pressure at the lowest altitude point in Bento Rodriguez are illustrated in Figure 9.These show that the slurry reached this location at as early as t = 825 s.The velocity and impact pressure decreased greatly as the village is located 5 km away downstream on the north-east side of the 'T-shaped' valley with a low slope.The velocity peak value was only 4.5 m/s and the impact pressure fluctuated between 13.7 kPa and 19 kPa.The submerged depth increased slowly with time, reaching 20.4 m at t = 1800 s.It can be inferred that, although the velocity and pressure had decreased greatly when the slurry reached the village, catastrophic damage was still unavoidable because of the huge volume of run-out tailings.The curves of the slurry submerged depth, velocity and impact pressure at the lowest altitude point in Bento Rodriguez are illustrated in Figure 9.These show that the slurry reached this location at as early as t = 825 s.The velocity and impact pressure decreased greatly as the village is located 5 km away downstream on the north-east side of the 'T-shaped' valley with a low slope.The velocity peak value was only 4.5 m/s and the impact pressure fluctuated between 13.7 kPa and 19 kPa.The submerged depth increased slowly with time, reaching 20.4 m at t = 1800 s.It can be inferred that, although the velocity and pressure had decreased greatly when the slurry reached the village, catastrophic damage was still unavoidable because of the huge volume of run-out tailings.
at as early as t = 825 s.The velocity and impact pressure decreased greatly as the village is located 5 km away downstream on the north-east side of the 'T-shaped' valley with a low slope.The velocity peak value was only 4.5 m/s and the impact pressure fluctuated between 13.7 kPa and 19 kPa.The submerged depth increased slowly with time, reaching 20.4 m at t = 1800 s.It can be inferred that, although the velocity and pressure had decreased greatly when the slurry reached the village, catastrophic damage was still unavoidable because of the huge volume of run-out tailings.

Modelling of an Operating Tailings Pond in China
Using the method described above, the modelling and prediction of the run-out slurry routing from the potential failure of an operating 'overhead tailings pond' (OTP) in Northern China was carried out.As Figure 10 shows, communities are located within only 0.5 km and industrial plants are located within 1 km downstream of this typical OTP.Based on the designed maximum capacity of 33 Mm 3 , dam height, dam shape and DSM data, the real-scale 3D geometry of the OTP and its surrounding area were rebuilt.A total of 3,900,144 boundary particles and 4,462,575 fluid particles were generated using the smoothing length of 3 m.The dam failure mode was set up to represent an instantaneous full collapse under the worst conditions, in order to predict the possible consequences of a major tailings dam failure.
Water 2018, 10, x FOR PEER REVIEW 10 of 15

Modelling of an Operating Tailings Pond in China
Using the method described above, the modelling and prediction of the run-out slurry routing from the potential failure of an operating 'overhead tailings pond' (OTP) in Northern China was carried out.As Figure 10 shows, communities are located within only 0.5 km and industrial plants are located within 1 km downstream of this typical OTP.Based on the designed maximum capacity of 33 Mm 3 , dam height, dam shape and DSM data, the real-scale 3D geometry of the OTP and its surrounding area were rebuilt.A total of 3,900,144 boundary particles and 4,462,575 fluid particles were generated using the smoothing length of 3 m.The dam failure mode was set up to represent an instantaneous full collapse under the worst conditions, in order to predict the possible consequences of a major tailings dam failure.The SPH numerical results in Figure 11 show that the slurry flow velocity increased sharply from the dam toe downstream in the first 100 s.Its peak value exceeded 30 m/s due to the steep slope under the dam toe.At t = 40 s, the run-out slurry submerged nearby buildings and pools, and continued to flow downstream at a high velocity of around 30 m/s.A small portion of the slurry flowed towards the south-east at a relatively low speed of around 15 m/s, but began to submerge the main communities 0.5 km away from the dam toe.At t = 100 s, most of the communities were influenced by the run-out slurry.At the same time, the slurry front flow divided into two branches at about 1 km downstream, flowing towards northern and south-eastern low-lying areas separately at a velocity of more than 20 m/s.At t = 200 s, the slurry flow velocity decreased greatly and the low-lying area on the north side was completely submerged.At this time all the buildings within 1 km downstream of the dam toe were also affected.Finally, at t = 400 s the slurry velocity at 1.5 km downstream became stable at about 10 m/s.The frontal portion of the flow beyond 1.5 km downstream was divided into two branches by a ground bulge, and flowed into the low-lying open pits on both the south and north sides.The SPH numerical results in Figure 11 show that the slurry flow velocity increased sharply from the dam toe downstream in the first 100 s.Its peak value exceeded 30 m/s due to the steep slope under the dam toe.At t = 40 s, the run-out slurry submerged nearby buildings and pools, and continued to flow downstream at a high velocity of around 30 m/s.A small portion of the slurry flowed towards the south-east at a relatively low speed of around 15 m/s, but began to submerge the main communities 0.5 km away from the dam toe.At t = 100 s, most of the communities were influenced by the run-out slurry.At the same time, the slurry front flow divided into two branches at about 1 km downstream, flowing towards northern and south-eastern low-lying areas separately at a velocity of more than 20 m/s.At t = 200 s, the slurry flow velocity decreased greatly and the low-lying area on the north side was completely submerged.At this time all the buildings within 1 km downstream of the dam toe were also affected.Finally, at t = 400 s the slurry velocity at 1.5 km downstream became stable at about flow downstream at a high velocity of around 30 m/s.A small portion of the slurry flowed towards the south-east at a relatively low speed of around 15 m/s, but began to submerge the main communities 0.5 km away from the dam toe.At t = 100 s, most of the communities were influenced by the run-out slurry.At the same time, the slurry front flow divided into two branches at about 1 km downstream, flowing towards northern and south-eastern low-lying areas separately at a velocity of more than 20 m/s.At t = 200 s, the slurry flow velocity decreased greatly and the low-lying area on the north side was completely submerged.At this time all the buildings within 1 km downstream of the dam toe were also affected.Finally, at t = 400 s the slurry velocity at 1.5 km downstream became stable at about 10 m/s.The frontal portion of the flow beyond 1.5 km downstream was divided into two branches by a ground bulge, and flowed into the low-lying open pits on both the south and north sides.To gain further insight into the properties of the slurry flow motion, the submerged depth, velocity and impact pressure versus time curves at 0.5 km downstream, where the communities are located, were plotted (Figure 12).These show that the slurry reached this location as early as t = 40 s.The submerged depth increased sharply due to the dam break, reaching a peak value of 37 m at t = 98 s.Subsequently, this depth decreased, reaching 18.9 m at t = 600 s.The slurry velocity peak value appeared at t = 64 s and the impact pressure peak value at t = 88 s.Thus, all the peak values appeared at a very early stage, before t = 100 s, which indicates that the downstream industrial plants and communities would suffer serious damage in an extremely short time if the tailings dam collapsed.This could be due to the huge designed capacity, short distance, large height difference and the narrow downstream topography.After t = 100 s, with the rapid reduction of the pond capacity, the impact pressure and velocity both show a significant downward trend.The submerged depth at this location also decreased.To gain further insight into the properties of the slurry flow motion, the submerged depth, velocity and impact pressure versus time curves at 0.5 km downstream, where the communities are located, were plotted (Figure 12).These show that the slurry reached this location as early as t = 40 s.The submerged depth increased sharply due to the dam break, reaching a peak value of 37 m at t = 98 s.Subsequently, this depth decreased, reaching 18.9 m at t = 600 s.The slurry velocity peak value appeared at t = 64 s and the impact pressure peak value at t = 88 s.Thus, all the peak values appeared at a very early stage, before t = 100 s, which indicates that the downstream industrial plants and communities would suffer serious damage in an extremely short time if the tailings dam collapsed.This could be due to the huge designed capacity, short distance, large height difference and the narrow downstream topography.After t = 100 s, with the rapid reduction of the pond capacity, the impact pressure and velocity both show a significant downward trend.The submerged depth at this location also decreased.To gain further insight into the properties of the slurry flow motion, the submerged depth, velocity and impact pressure versus time curves at 0.5 km downstream, where the communities are located, were plotted (Figure 12).These show that the slurry reached this location as early as t = 40 s.The submerged depth increased sharply due to the dam break, reaching a peak value of 37 m at t = 98 s.Subsequently, this depth decreased, reaching 18.9 m at t = 600 s.The slurry velocity peak value appeared at t = 64 s and the impact pressure peak value at t = 88 s.Thus, all the peak values appeared at a very early stage, before t = 100 s, which indicates that the downstream industrial plants and communities would suffer serious damage in an extremely short time if the tailings dam collapsed.This could be due to the huge designed capacity, short distance, large height difference and the narrow downstream topography.After t = 100 s, with the rapid reduction of the pond capacity, the impact pressure and velocity both show a significant downward trend.The submerged depth at this location also decreased.

Significance of This Study
Serious tailings dam failure accidents have occurred globally in recent years.The disaster prevention and mitigation topic is attracting more and more attention from both government and stakeholders, especially in China, with 8869 tailings ponds all over the country, including 1429 potentially dangerous OTPs.Routing of the run-out slurry flow from tailings dam failure is known to be a complex, large deformation, destructive and uncontrollable process.Numerical simulations, especially the advanced, meshless SPH method, which has the advantages of low costs, flexible settings and high efficiency, is playing an important role in disaster prevention and mitigation.The results could provide evidence for disaster mitigation and protection engineering layout and emergency management plan optimization.Globally, the sustainability and viability of long-term mining communities are considered to be increasingly critical [38,39].The OTP problems in China now are mainly caused by the short-term considerations of tailings disposal plans when the mines were first built in the 20th century.Therefore, in contemporary mining, the mine site location, tailings dam construction method and mine waste disposal plan must be more responsible and sustainable, to cover the mine's whole service life.Knowledge of the potential consequences of dam failure at designed capacity through modelling, as is presented in this study, could be helpful for better and more sustainable site design.

Proposed Method Applicability
The proposed method was adopted to simulate the laboratory physical modelling experiment based on the operating Yangtianqing tailings pond in China and the known Brazil Fundão tailings dam failure accident for validation.The comparison confirmed the validity of the SPH method.
In the Yangtianqing tailings pond case study, the numerical calculated results of the submersion depths, velocity, impact pressure and vortex phenomenon at different time steps were all in line with the data measured in the laboratory experiment.Furthermore, it can be inferred that the Mimao community, located 0.8 km downstream of the tailings dam, would suffer serious damage at the early stage of a potential tailings dam failure, due to the large impact pressure and submerged depth modelled.Therefore, more analysis and assessment of the current design need to be carried out before further operation.However, due to the limitations of the measurement resolution and sensitivity, as well as the complex dam shape and failure triggering modes, more advanced and sophisticated experiments are needed for further validation.
In the Brazil Fundão tailings dam failure case study, the calculated results of the direction of the tailings run-out flow and the impact area generally agree with the panel report and the actual consequences shown in the satellite images.More details and data would be needed for further validation.The significant agreements above indicate the applicability of the SPH method when modelling tailings slurry flow motions.

Serious Risks of the OTPs
Modelling the routing of the run-out of tailings from a typical OTP in China under the worst conditions and with the maximum designed capacity was carried out.The results showed the high risks of this typical OTP.Downstream communities and other important facilities would be submerged in an extremely short period of time, leaving very limited time to evacuate residents and conduct further emergency response.The potentially destructive power of the tailings slurry was emphasized by the high impact pressure and flow velocity of the tailings slurry at 0.5 km downstream, where the communities are located.
The impacts of OTPs on the local ecosystem and resident health should also be considered in future studies, for example through agricultural irrigation, potable water, air pollution and dust emission.The good news is that both the government and stakeholders in China are paying increasing attention to OTP safety and environmental risks.A series of effective campaigns have been launched since 2016, such as OTP closure, tailings reuse, land reclamation, strict surveillance standards, emergency response plan exercises and community migrations.

Inadequacies of This Study
There are still some inadequacies in this study that need to be considered in future studies.Firstly, the run-out slurry from tailings dam failure is known to be a multi-phase complex fluid.The rheological model and parameters are determined by the dam failure mode, slurry concentration, tailings particle size, spatial distribution and other physical parameters.The parameters of the tailings slurry run-out lacked field measured data and were mostly simplified in previous physical and numerical studies.
Secondly, the resolution of the satellite remote sensing DSM data was relatively low and the data sampling process usually lasts for several years.For example, the JAXA AW3D30 DSM data used in this study was collected between 2006 and 2011, with a resolution of approximately 30 m [26].With this poor resolution, infrastructure is not well-represented.Higher resolution and more recent terrain data such as those sampled by the unmanned aerial vehicle (UAV) photogrammetry [40] would be helpful to produce more reliable and detailed numerical results.Thirdly, as tailings ponds and their potential impact regions cover very wide areas, millions of boundary and fluid particles would be generated when these are modelled.Computational efficiency is thus becoming one of the bottlenecks restricting broader future industrial applications of this study.
Finally, downstream land surface roughness, vegetation and sediment may also have great influence on the flow pattern and run-out slurry volume.These interacting mechanisms need to be considered in future studies.

Conclusions
In this study, a 3D real-scale numerical modelling method based on remote sensing digital surface model (DSM) data and the advanced meshless smoothed particle hydrodynamics (SPH) method was used to study the routing of tailings slurry run-out across downstream terrains, to determine its impact area during potential tailings dam failures.Case studies using a physical experimental test, and based on the 2015 Brazil Fundão tailings dam failure accident and one operating tailings pond in China were carried out.Key conclusions are:

•
The physical modelling flume experiment and the 2015 Brazil Fundão accident consequences were restored and validated using the SPH method.This work showed that the tailings dam failure run-out slurry routing across the real terrain could be solved using the proposed method.Parameters such as impact area, submersion depths, impact pressure and slurry velocity should be obtained for further analysis.

•
Modelling of high-risk 'overhead tailings ponds' (OTPs) showed that the tailings slurry run-out would be tremendously destructive and the emergency response time would be extremely short if tailings dam failure occurred with the pond at full designed capacity.Further analysis and assessment is recommended before further operation of these OTPs.

•
The results could provide evidence for disaster prevention and mitigation engineering layout, emergency management plan optimization and the development of more responsible site plans and sustainable site designs.However, improvements such as the rheological model selection, terrain data quality, computing efficiency and land surface roughness need to be considered in future studies.

•
SPH numerical modelling is a powerful and advanced technique that is recommended for hazard assessments and the sustainable design of tailings dam facilities globally.

Figure 1 .
Figure 1.(a) Locations of Yunnan Province and the Yangtianqing tailings pond; (b) Locations of the Yangtianqing tailings pond and downstream villages (Google Earth).

Figure 1 .
Figure 1.(a) Locations of Yunnan Province and the Yangtianqing tailings pond; (b) Locations of the Yangtianqing tailings pond and downstream villages (Google Earth).

Figure 2 .
Figure 2. Experiment flume setup and dimensions (mm): (a) the slurry mixing device and storage tank; (b) the downstream flume; (c) the piezometer and data sampling system.

Figure 2 .
Figure 2. Experiment flume setup and dimensions (mm): (a) the slurry mixing device and storage tank; (b) the downstream flume; (c) the piezometer and data sampling system.

Water 2018 , 15 Figure 2 .
Figure 2. Experiment flume setup and dimensions (mm): (a) the slurry mixing device and storage tank; (b) the downstream flume; (c) the piezometer and data sampling system.

Figure 3 .
Figure 3. Top view snapshots of the SPH modelled tailings slurry flow motion in the flume.

Figure 4 .
Figure 4. Comparison between the numerical and experimental snapshots of slurry flow at the sharp bend: (a) experimental result; (b) numerical particle trajectories.

Figure 5 .
Figure 5.Comparison of slurry depth versus time at x = 6.5 m (from t = 1.6 s) and x = 9 m (from t = 3.8 s).

Figure 3 . 15 Figure 3 .
Figure 3. Top view snapshots of the SPH modelled tailings slurry flow motion in the flume.

Figure 4 .
Figure 4. Comparison between the numerical and experimental snapshots of slurry flow at the sharp bend: (a) experimental result; (b) numerical particle trajectories.

Figure 5 .
Figure 5.Comparison of slurry depth versus time at x = 6.5 m (from t = 1.6 s) and x = 9 m (from t = 3.8 s).

Figure 4 .
Figure 4. Comparison between the numerical and experimental snapshots of slurry flow at the sharp bend: (a) experimental result; (b) numerical particle trajectories.

Water 2018 , 15 Figure 3 .
Figure 3. Top view snapshots of the SPH modelled tailings slurry flow motion in the flume.

Figure 4 .
Figure 4. Comparison between the numerical and experimental snapshots of slurry flow at the sharp bend: (a) experimental result; (b) numerical particle trajectories.

Figure 5 .
Figure 5.Comparison of slurry depth versus time at x = 6.5 m (from t = 1.6 s) and x = 9 m (from t = 3.8 s).

Figure 5 .
Figure 5.Comparison of slurry depth versus time at x = 6.5 m (from t = 1.6 s) and x = 9 m (from t = 3.8 s).

Figure 6 .
Figure 6.Comparison of slurry impact pressure vs time at the start of the sharp bend.

Figure 7 .
Figure 7. Locations of the Samarco collapsed tailings pond, submerged area and community.

Figure 6 .
Figure 6.Comparison of slurry impact pressure vs. time at the start of the sharp bend.

Water 2018 , 15 Figure 6 .
Figure 6.Comparison of slurry impact pressure vs time at the start of the sharp bend.

Figure 7 .
Figure 7. Locations of the Samarco collapsed tailings pond, submerged area and community.

Figure 7 .
Figure 7. Locations of the Samarco collapsed tailings pond, submerged area and community.

Figure 8 .
Figure 8. SPH numerical results of the Samarco tailings dam failure accident.

Figure 9 .Figure 8 .
Figure 9. Numerical results of submerged depth, velocity and impact pressure at Bento Rodriguez.

Figure 9 .
Figure 9. Numerical results of submerged depth, velocity and impact pressure at Bento Rodriguez.

Figure 9 .
Figure 9. Numerical results of submerged depth, velocity and impact pressure at Bento Rodriguez.

Figure 10 .
Figure 10.Locations of the 'overhead tailings pond' and nearby important facilities.

Figure 10 .
Figure 10.Locations of the 'overhead tailings pond' and nearby important facilities.
s.The frontal portion of the flow beyond 1.5 km downstream was divided into two branches by a ground bulge, and flowed into the low-lying open pits on both the south and north sides.

Figure 12 .
Figure 12.Numerical results of submerged depth, velocity and impact pressure at 0.5 km downstream of the dam toe.

Figure 12 .
Figure 12.Numerical results of submerged depth, velocity and impact pressure at 0.5 km downstream of the dam toe.

1 .Figure 12 .
Figure 12.Numerical results of submerged depth, velocity and impact pressure at 0.5 km downstream of the dam toe.