3D Numerical Modelling of Tailings Dam Breach Run Out Flow over Complex Terrain: A Multidisciplinary Procedure

Tailings dams, as essential mining structures, are being built globally for containing the chief waste stream of the mining industry. Catastrophic tailings dam breaches have occurred frequently over the past decade, causing severe impacts on the environment, economy, and human health. The foreknowledge of the tailings dam breach overland flow is crucial for the risk assessment and emergency response planning in order to prevent or minimize possible losses. Using unmanned aerial vehicles (UAV) photogrammetry and smoothed particle hydrodynamics (SPH) numerical method, this study proposed a multidisciplinary procedure for modelling a hypothetical tailings dam breach run out flow over the downstream complex terrain. A case study on a 97-m-height tailings dam in Shandong Province of China was carried out. The proposed procedure was proven applicable to determine the overland tailings flow. The submerged area and flow velocities suggested that the downstream G2 highway would hardly be threatened and more concerns should be paid on the factory plants and workers deployed between the dam toe and the highway. Additionally, the application of UAV photogrammetry in the mining industry as a supplementary surveying method can be further expanded, especially for the numerous small-scale mining sites. The proposed procedure is then recommended for the safety management of the tailings’ storage facilities globally.


Introduction
A tailings storage facility (TSF) is an essential structure in the mining industry built for the purpose of storing unwanted and currently uneconomical mine 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 operation and emergency 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 over the 100 years between 1910 and 2010 was estimated to be 1.2%, which is more than two orders of magnitude higher than that of normal water-retaining dams (0.01%). In the last two decades, with the great improvements in damming practices, regulation, and modern construction technology, the TDB accident rate has sufficiently

Study Area and Proposed Procedure
The study site of the Yujiaquan TSF is situated in Laiwu city, Shandong Province of Eastern China, as Figure 2 shows. The TSF has been operating since 1985 in order to store mining waste from the Luzhong iron mine. It was designed to deposit tailings behind a 29-m-high compacted rockfill starter dam (El. + 256 m ~El. + 285 m) and then augment capacity by upstream raising the main dam with the slope ratio of 1/4. The annual disposal of tailings is estimated to be 2.6 million tons, which rises the TSF reservoir surface by 2.4 m every year. To ensure a safe and sustainable tailings disposal operation, the mining waste rocks, the relatively coarse tailings from the hydrocyclone underflow, and the permeable geotextile are utilized for the dam construction, as Figure 3 shows. The finegrained hydrocyclone overflow slurry is then discharged into the TSF. The TSF designed maximum capacity was 52.55 Mm 3 with the dam height of 114 m (El. 370 m). It can also be seen in Figure 2 that a densely populated community with more than 1000 residents is located only 1.2 km downstream of the TSF. In addition, the G2 highway, which links the capital city Beijing in Northern China and the central major city Shanghai, is crossing through the study area from the northwest to the southeast, which is approximately 1 km from the dam toe. Several factory plants and a water recovery pond are deployed on the flat area between the dam toe and the G2 highway.

Study Area and Proposed Procedure
The study site of the Yujiaquan TSF is situated in Laiwu city, Shandong Province of Eastern China, as Figure 2 shows. The TSF has been operating since 1985 in order to store mining waste from the Luzhong iron mine. It was designed to deposit tailings behind a 29-m-high compacted rockfill starter dam (El. + 256 m~El. + 285 m) and then augment capacity by upstream raising the main dam with the slope ratio of 1/4. The annual disposal of tailings is estimated to be 2.6 million tons, which rises the TSF reservoir surface by 2.4 m every year. To ensure a safe and sustainable tailings disposal operation, the mining waste rocks, the relatively coarse tailings from the hydrocyclone underflow, and the permeable geotextile are utilized for the dam construction, as Figure 3 shows. The fine-grained hydrocyclone overflow slurry is then discharged into the TSF. The TSF designed maximum capacity was 52.55 Mm 3 with the dam height of 114 m (El. 370 m). It can also be seen in Figure 2 that a densely populated community with more than 1000 residents is located only 1.2 km downstream of the TSF. In addition, the G2 highway, which links the capital city Beijing in Northern China and the central major city Shanghai, is crossing through the study area from the northwest to the southeast, which is approximately 1 km from the dam toe. Several factory plants and a water recovery pond are deployed on the flat area between the dam toe and the G2 highway.

Study Area and Proposed Procedure
The study site of the Yujiaquan TSF is situated in Laiwu city, Shandong Province of Eastern China, as Figure 2 shows. The TSF has been operating since 1985 in order to store mining waste from the Luzhong iron mine. It was designed to deposit tailings behind a 29-m-high compacted rockfill starter dam (El. + 256 m ~El. + 285 m) and then augment capacity by upstream raising the main dam with the slope ratio of 1/4. The annual disposal of tailings is estimated to be 2.6 million tons, which rises the TSF reservoir surface by 2.4 m every year. To ensure a safe and sustainable tailings disposal operation, the mining waste rocks, the relatively coarse tailings from the hydrocyclone underflow, and the permeable geotextile are utilized for the dam construction, as Figure 3 shows. The finegrained hydrocyclone overflow slurry is then discharged into the TSF. The TSF designed maximum capacity was 52.55 Mm 3 with the dam height of 114 m (El. 370 m). It can also be seen in Figure 2 that a densely populated community with more than 1000 residents is located only 1.2 km downstream of the TSF. In addition, the G2 highway, which links the capital city Beijing in Northern China and the central major city Shanghai, is crossing through the study area from the northwest to the southeast, which is approximately 1 km from the dam toe. Several factory plants and a water recovery pond are deployed on the flat area between the dam toe and the G2 highway.

UAVs SfM-Photogrammetry and Field Surveys
UAVs have recently been used for the surveillance of heritage, vegetation, landslides, and infrastructure engineering [17][18][19][20][21]. In addition to UAV technology, structure from motion (SfM) photogrammetry is a low-cost topographical survey technique for estimating three-dimensional structures from a two-dimensional image sequence. The SfM technique is described in detail by Westoby, et al. [22]. Recent advances in photogrammetry software, high-performance computing, and consumer-grade UAVs have made SfM photogrammetry a reliable topographic survey technique in terms of reasonable costs, great flexibility, and efficiency. UAVs photogrammetry has become a matter of growing interest among both researchers and industries. Rauhala et al. [23] had successfully adopted the technique for the surveillance of a closed TSF in a subarctic area in Finland. The obtained results indicated that UAVs remote sensing can produce sufficiently accurate data for supporting tailings management operations and tracking annual surface displacements in the decimeter range. Barreiro, et al. [24] applied UAV photogrammetry and the SPH method to compute the trajectory of the water flow across a close road during extreme weather events under different protection scenarios.
In the present study, the photogrammetric field surveys were performed in four stages: (1) mission planning, (2) acquisition of ground control points (GCPs), (3) flights implementation and images acquisition, and (4) photogrammetry processing. The aerial inspection devices used in this work were two four-propeller UAV drones from DJI (DJI, Shenzhen, China) (Shown in Table 1). One was the Phantom 4 Advanced model equipped with a digital FC6310 camera (1-inch 20-megapixel CMOS sensor, 24-mm f/2.8 wide-angle lens). The other one was the Phantom 3 Professional model equipped with a digital FC300X camera (1/2.3-inch 12-megapixel sensor, 20-mm f/2.8 wide-angle lens). The field images acquisition was performed on 1 August 2018. The sky condition was generally clear. As shown in Figure 4a,b, the UAVs were deployed on site. During the flights, the UAVs were controlled autonomously using the Pix4Dcapture mobile application (Pix4D, Prilly, Switzerland). The flying height was 100 m with the frontal overlap of 75% and side overlap of 60%. The GCPs distributed evenly on the main dam crest and the dam body were recorded with a Hi-Target V30 GNSS RTK receiver (Hi-Target Surveying Instrument, Guangzhou, China) and several self-made GCP mark boards, as Figure 4c,d shows. The static positioning accuracy of 2.5 mm and the elevation accuracy of 5 mm can be achieved. Eventually, 3002 images and 33 GCPs data covering the 2.28 km 2 study area were collected. The photogrammetry processing was then performed on the highperformance cluster with 69 GB RAM and Intel ® Xeon ® Platinum 8124M CPU (Intel Corporation, Santa Clara, CA, USA) at 3.00 GHz.

UAVs SfM-Photogrammetry and Field Surveys
UAVs have recently been used for the surveillance of heritage, vegetation, landslides, and infrastructure engineering [17][18][19][20][21]. In addition to UAV technology, structure from motion (SfM) photogrammetry is a low-cost topographical survey technique for estimating three-dimensional structures from a two-dimensional image sequence. The SfM technique is described in detail by Westoby et al. [22]. Recent advances in photogrammetry software, high-performance computing, and consumer-grade UAVs have made SfM photogrammetry a reliable topographic survey technique in terms of reasonable costs, great flexibility, and efficiency. UAVs photogrammetry has become a matter of growing interest among both researchers and industries. Rauhala et al. [23] had successfully adopted the technique for the surveillance of a closed TSF in a subarctic area in Finland. The obtained results indicated that UAVs remote sensing can produce sufficiently accurate data for supporting tailings management operations and tracking annual surface displacements in the decimeter range. Barreiro et al. [24] applied UAV photogrammetry and the SPH method to compute the trajectory of the water flow across a close road during extreme weather events under different protection scenarios.
In the present study, the photogrammetric field surveys were performed in four stages: (1) mission planning, (2) acquisition of ground control points (GCPs), (3) flights implementation and images acquisition, and (4) photogrammetry processing. The aerial inspection devices used in this work were two four-propeller UAV drones from DJI (DJI, Shenzhen, China). One was the Phantom 4 Advanced model equipped with a digital FC6310 camera (1-inch 20-megapixel CMOS sensor, 24-mm f/2.8 wide-angle lens). The other one was the Phantom 3 Professional model equipped with a digital FC300X camera (1/2.3-inch 12-megapixel sensor, 20-mm f/2.8 wide-angle lens). The field images acquisition was performed on 1 August 2018. The sky condition was generally clear. As shown in Figure 4a,b, the UAVs were deployed on site. During the flights, the UAVs were controlled autonomously using the Pix4Dcapture mobile application (Pix4D, Prilly, Switzerland). The flying height was 100 m with the frontal overlap of 75% and side overlap of 60%. The GCPs distributed evenly on the main dam crest and the dam body were recorded with a Hi-Target V30 GNSS RTK receiver (Hi-Target Surveying Instrument, Guangzhou, China) and several self-made GCP mark boards, as Figure 4c,d shows. The static positioning accuracy of 2.5 mm and the elevation accuracy of 5 mm can be achieved. Eventually, 3002 images and 33 GCPs data covering the 2.28 km 2 study area were collected. The key features of UAV photogrammetry data acquisition are shown in Table 1. The photogrammetry processing was then performed on the high-performance cluster with 69 GB RAM and Intel ® Xeon ® Platinum 8124M CPU (Intel Corporation, Santa Clara, CA, USA) at 3.00 GHz.

Smoothed Particle Hydrodynamics (SPH) Method
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. 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: 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 fundamental to 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(ra − rb, h) is correlated with smoothing length h and (ra − rb), i.e., the distance between

Smoothed Particle Hydrodynamics (SPH) Method
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. 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: 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 fundamental to 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 smoothing length h and (r a − r b ), i.e., the distance between particles a and b. In this study, the fifth-order quantic kernel by Wendland [25] where the weighting function vanishes for inter-particle distances greater than 2h is used. The Wendland kernel is defined in 3D as: where q = (r a − r b )/h and α D is the normalization constant that equals to 7/(4πh 2 ) in 2D and 21/(16πh 3 ) in 3D.
To simplify the calculation of the pressure term, the weakly compressible state equation proposed by Monaghan [26] 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 and ρ 0 is the reference density. 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.
In most reported tailings dam breach cases and studies [27][28][29], the run out tailings flow is normally considered to be the homogeneous mud-water mixture similarly to debris flow at a low concentration of less than 40%. In addition, its rheological behavior is often described as homogeneous non-Newtonian fluid [29][30][31]. The generalized Herschel-Bulkley-Papanastasiou (HBP) model proposed by Papanastasiou [32] is adopted in this study. The effective viscosity (µ e f f ) can be defined as: where K and m are the constant coefficients, .
y is the shear rate and τ y is the yield stress. When n = 1, the HBP model reduces to a simple Bingham model. The HBP model is described in detail in the research [31][32][33]. In the present study, the tailings' run out flow was regarded as 40% solid concentration homogeneous slurry with the yield stress of 2.5 Pa and the viscosity of 0.15 Pa·s.
The fluid density fluctuation is computed by solving the mass conservation as the particle mass keeps constant in the weakly compressible SPH (WCSPH) computing. The continuity equation in SPH is as follows: To account for the WCSPH, the Delta-SPH equation proposed by Molteni and Colagrossi [34] was implemented. A diffusive term is introduced to reduce the density fluctuations. The Delta-SPH equation can be written in the following form.
where δ is the diffusive coefficient of Delta-SPH.
The open-source SPH solver DualSPHsysics v4.0 [35] (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. The basic SPH constants used were based on the suggested values of the 3D simulation cases in Duals physics. The code was then executed on the Graphics Processing Unit (GPU) node of a high-performance cluster (HPC) with the NVIDIA Tesla K80 GPU (4992 processor cores with a dual-GPU design, 560 MHz, NVIDIA, Santa Clara, CA, USA).
Based on the produced high-resolution UAV SfM DSM result, the real-scale 3D geometry of the research area was constructed after a series of transitions. First, the SfM 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 three times of the breach depths [36], 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 [37]. Based on the fact that the elevation of the 29-m-high starter dam was from +256 m to +285 m and the measured average elevation of the reservoir was +349.27 m. The TDB width was then set to be 279.81 m (top breach width) and 216.91 m (bottom breach width). The TDB mode was simplified to be under extreme conditions where the breach initiates instantaneously at the first step of the simulation. On the basis of available computational capacities, the particle distance was set to 2 m, which eventually generated 936,381 boundary particles and 3,494,863 fluid particles.

UAV SfM Photogrammetry
The orthophoto and Digital Surface Model (DSM) of the research area are shown in Figures 5 and 6b. Several factory plants, a water recovery pond, and the G2 highway can be clearly identified on the southwest downstream direction. The TSF waste water delivered by the underdrains is discharged into the 2.2 × 10 4 m 2 water recovery pond from where it will be pumped back to the processing plant for reuse. The perimeter and area of the TSF are measured to be 2934 m and 5.11 × 10 5 m 2 , respectively. The average elevations of the dam crest and the reservoir are calculated to be +353.75 m and +349.27 m. The beach width, which indicates the proximity of the supernatant water to the tailings dam and is often regarded as an important surveillance factor can also be clearly observed, as well as the ongoing greening work on the rock-soil embankment in order to control ambient environmental pollution. Meanwhile, the elevation, width, and length of the concrete auxiliary dam crest on the north side are measured to be +355 m, 6 m, and 225 m, respectively. The result ties well with the field survey records. The extracted data can provide evidence for the tailings' discharge planning and operation.
Water 2020, 12, x FOR PEER REVIEW 7 of 14 simulation cases in Duals physics. The code was then executed on the Graphics Processing Unit (GPU) node of a high-performance cluster (HPC) with the NVIDIA Tesla K80 GPU (4992 processor cores with a dual-GPU design, 560 MHz, NVIDIA, Santa Clara, CA, USA). Based on the produced high-resolution UAV SfM DSM result, the real-scale 3D geometry of the research area was constructed after a series of transitions. First, the SfM 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 three times of the breach depths [36], 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 [37]. Based on the fact that the elevation of the 29-m-high starter dam was from +256 m to +285 m and the measured average elevation of the reservoir was +349.27 m. The TDB width was then set to be 279.81 m (top breach width) and 216.91 m (bottom breach width). The TDB mode was simplified to be under extreme conditions where the breach initiates instantaneously at the first step of the simulation. On the basis of available computational capacities, the particle distance was set to 2 m, which eventually generated 936,381 boundary particles and 3,494,863 fluid particles.

UAV SfM Photogrammetry
The orthophoto and Digital Surface Model (DSM) of the research area are shown in Figure 5 and Figure 6b. Several factory plants, a water recovery pond, and the G2 highway can be clearly identified on the southwest downstream direction. The TSF waste water delivered by the underdrains is discharged into the 2.2 × 10 4 m 2 water recovery pond from where it will be pumped back to the processing plant for reuse. The perimeter and area of the TSF are measured to be 2934 m and 5.11 × 10 5 m 2 , respectively. The average elevations of the dam crest and the reservoir are calculated to be +353.75 m and +349.27 m. The beach width, which indicates the proximity of the supernatant water to the tailings dam and is often regarded as an important surveillance factor can also be clearly observed, as well as the ongoing greening work on the rock-soil embankment in order to control ambient environmental pollution. Meanwhile, the elevation, width, and length of the concrete auxiliary dam crest on the north side are measured to be +355 m, 6 m, and 225 m, respectively. The result ties well with the field survey records. The extracted data can provide evidence for the tailings' discharge planning and operation.

Tailings Dam Breach Overland Flow
The snapshots of simulated hypothetical dam breach ROOF velocities at 50 s, 150 s, 300 s, and 600 s are presented in Figure 7. The flow velocity increased greatly at 50 s, as the gravitational potential energy of the tailings flow is swiftly transferred into kinetic energy. In addition, the velocity peak value of 22.6 m/s was identified at the dam toe. Thereafter, with the rapid decrease, the TSF capacity and the velocities of ROOF significantly reduced at 150 s, 300 s, and 600 s. The velocity peak value had reduced to less than 15 m/s at 150 s and began to impact the downstream water recovery pond. The front flow velocity was lower than 5 m/s. Thereafter, at 300 s and 600 s, the flow was spreading at less than 4 m/s over the flat area where the water recovery pond and factory plants are located in. At 600 s, a fan-shaped submerged area had formed and the front flow was approaching the G2 highway at less than 2.5 m/s. The maximum run out distance from the dam crest was measured to be 1224 m with only 78 m from the G2 highway. At 600 s, the maximum submerged depth at the low-lying water recovery pond (655 m downstream from the dam crest) was 24 m. Furthermore, the obtained submerged depth of the front flow was 4~10 m. In addition, the front flow was accumulating at several low-lying spots and several earthworks were built along the highway. Therefore, it can be inferred that the concerned G2 highway was at a safety distance from the Yujiaquan TSF, as the volume and velocity of the approaching front flow had reduced greatly, the road bed altitude was relatively higher, and the emergency response time was sufficient. However, the safety of the factory plants and workers deployed on the downstream flat area should be further considered. In addition, the freely available "Advanced Land Observing Satellite (ALOS) World 3D 30 m mesh" (AW3D30) dataset provided by Japan Aerospace Agency (JAXA) are extracted and illustrated in Figure 6a. The dataset was produced based on the data acquired from 2006 to 2011 by the ALOS satellite and the resolution was 1 arc second (approximately 30 m) [38]. A comparative discussion between the AW3D30 and UAV SfM DSMs is presented in Section 4.2.

Tailings Dam Breach Overland Flow
The snapshots of simulated hypothetical dam breach ROOF velocities at 50 s, 150 s, 300 s, and 600 s are presented in Figure 7. The flow velocity increased greatly at 50 s, as the gravitational potential energy of the tailings flow is swiftly transferred into kinetic energy. In addition, the velocity peak value of 22.6 m/s was identified at the dam toe. Thereafter, with the rapid decrease, the TSF capacity and the velocities of ROOF significantly reduced at 150 s, 300 s, and 600 s. The velocity peak value had reduced to less than 15 m/s at 150 s and began to impact the downstream water recovery pond. The front flow velocity was lower than 5 m/s. Thereafter, at 300 s and 600 s, the flow was spreading at less than 4 m/s over the flat area where the water recovery pond and factory plants are located in. At 600 s, a fan-shaped submerged area had formed and the front flow was approaching the G2 highway at less than 2.5 m/s. The maximum run out distance from the dam crest was measured to be 1224 m with only 78 m from the G2 highway. At 600 s, the maximum submerged depth at the low-lying water recovery pond (655 m downstream from the dam crest) was 24 m. Furthermore, the obtained submerged depth of the front flow was 4~10 m. In addition, the front flow was accumulating at several low-lying spots and several earthworks were built along the highway. Therefore, it can be inferred that the concerned G2 highway was at a safety distance from the Yujiaquan TSF, as the volume and velocity of the approaching front flow had reduced greatly, the road bed altitude was relatively higher, and the emergency response time was sufficient. However, the safety of the factory plants and workers deployed on the downstream flat area should be further considered. Water 2020, 12, x FOR PEER REVIEW 9 of 14

Quality of the SfM Photogrammetry Results
The accuracy and quality of the photogrammetry results are determined by the Ground Sampling Distance (GSD), which means the distance between the center of two consecutive pixels on the ground surface. A lower GSD value indicates a better and more detailed spatial resolution of the results. The GSD is decided by the real focal length FR (mm), the flight height H (m), and the covered width of one single image DW (m). FR (mm) can be calculated by: where F35 (mm) is the 35-mm equivalent focal length and SW (mm) is the real sensor width. Furthermore, constant k equals to 34.6 in the 4:3 aspect ratio and 36 in the 3:2 aspect ratio. The covered distance on the ground surface by one single image DW (m) is given by: being Wim is the image width (pixel) and GSD is the expected Ground Sampling Distance (cm/pixel). Given the fact that: = Therefore, combining Equations (9)-(11), the GSD can be calculated by:

Quality of the SfM Photogrammetry Results
The accuracy and quality of the photogrammetry results are determined by the Ground Sampling Distance (GSD), which means the distance between the center of two consecutive pixels on the ground surface. A lower GSD value indicates a better and more detailed spatial resolution of the results. The GSD is decided by the real focal length F R (mm), the flight height H (m), and the covered width of one single image D W (m). F R (mm) can be calculated by: where F 35 (mm) is the 35-mm equivalent focal length and S W (mm) is the real sensor width. Furthermore, constant k equals to 34.6 in the 4:3 aspect ratio and 36 in the 3:2 aspect ratio. The covered distance on the ground surface by one single image D W (m) is given by: being W im is the image width (pixel) and GSD is the expected Ground Sampling Distance (cm/pixel). Given the fact that: Therefore, combining Equations (9)-(11), the GSD can be calculated by: According to the manufacturer's specifications, the two UAV drones (DJI Phantom 4 Advanced and DJI Phantom 3 Professional) used in the present study are equipped with a 24-mm FC6310 camera and a 20-mm FC300X camera (35 mm equivalent focal length), respectively. The aspect ratio of the FC6310 camera is 3:2 (5472 × 3648) and the aspect ratio of FC300X camera is 4:3 (4000 × 3000). Therefore, the GSD from the DJI Phantom 4 Advanced is calculated to be 2.74 cm/pixel. In the same manner, the GSD from the DJI Phantom 3 Professional is 4.33 cm/pixel. According to the processing report, the average GSD is 3.8 cm/pixel.
As Equation (12) suggests, the result quality is highly impacted by the image acquisition parameters and aerial inspection devices. Generally, lower flight heights, longer focal lengths, and higher image widths could decrease the GSD value and produce better results. Focal lengths of 18-24 mm are frequently used for aerial inspection campaigns because lower flight heights require a wider-angle lens for aerial coverage and image overlap. Moreover, longer focal length lenses are generally heavier, which could shorten the UAV endurance. In addition, proper camera settings, e.g., a low ISO value, a broad aperture, a high shutter speed, RAW image file format, and good weather conditions are also critical to produce better results.

Comparison with Spaceborne DSM Product
The freely available AW3D30 dataset shown in Figure 6a was used for the comparison study. The resolution of the AW3D30 is 1 arc second (approximately 30 m), which is much lower than the UAV SfM DSM. Moreover, significant elevation changes of the study area could be observed, as the AW3D30 DSM was produced based on the dataset sampled from 2006 to 2011, 7~12 years earlier than the present UAV field survey. Subsequently, DSMs differences between the AW3D30 and SfM datasets and the contour lines ( Figure 8a) were generated by calculating the pixel-by-pixel elevation differences. The elevation differences distribution is plotted and shown in Figure 8b. The mean and median of the histogram are 15.4 m and 16 m, respectively. The middle 80% of the values range between 1 m and 29 m. It clearly shows that most of the study area, specifically the TSF reservoir surface, the main dam, the G2 highway, and its adjacent area had significantly risen by 10 m~40 m due to the tailings' disposal activities and highway roadbed constructions occurring during the 7~12 years period. In addition, a negative anomaly could be observed on the north-western and northern edges mainly due to the ambient safety maintenance activities during the TSF expansion.
Water 2020, 12, x FOR PEER REVIEW 10 of 14 = 100 × × × (12) According to the manufacturer's specifications, the two UAV drones (DJI Phantom 4 Advanced and DJI Phantom 3 Professional) used in the present study are equipped with a 24-mm FC6310 camera and a 20-mm FC300X camera (35 mm equivalent focal length), respectively. The aspect ratio of the FC6310 camera is 3:2 (5472 × 3648) and the aspect ratio of FC300X camera is 4:3 (4000 × 3000). Therefore, the GSD from the DJI Phantom 4 Advanced is calculated to be 2.74 cm/pixel. In the same manner, the GSD from the DJI Phantom 3 Professional is 4.33 cm/pixel. According to the processing report, the average GSD is 3.8 cm/pixel.
As Equation (12) suggests, the result quality is highly impacted by the image acquisition parameters and aerial inspection devices. Generally, lower flight heights, longer focal lengths, and higher image widths could decrease the GSD value and produce better results. Focal lengths of 18-24 mm are frequently used for aerial inspection campaigns because lower flight heights require a widerangle lens for aerial coverage and image overlap. Moreover, longer focal length lenses are generally heavier, which could shorten the UAV endurance. In addition, proper camera settings, e.g., a low ISO value, a broad aperture, a high shutter speed, RAW image file format, and good weather conditions are also critical to produce better results.

Comparison with Spaceborne DSM Product
The freely available AW3D30 dataset shown in Figure 6a was used for the comparison study. The resolution of the AW3D30 is 1 arc second (approximately 30 m), which is much lower than the UAV SfM DSM. Moreover, significant elevation changes of the study area could be observed, as the AW3D30 DSM was produced based on the dataset sampled from 2006 to 2011, 7~12 years earlier than the present UAV field survey. Subsequently, DSMs differences between the AW3D30 and SfM datasets and the contour lines ( Figure 8a) were generated by calculating the pixel-by-pixel elevation differences. The elevation differences distribution is plotted and shown in Figure 8b. The mean and median of the histogram are 15.4 m and 16 m, respectively. The middle 80% of the values range between 1 m and 29 m. It clearly shows that most of the study area, specifically the TSF reservoir surface, the main dam, the G2 highway, and its adjacent area had significantly risen by 10 m~40 m due to the tailings' disposal activities and highway roadbed constructions occurring during the 7~12 years period. In addition, a negative anomaly could be observed on the north-western and northern edges mainly due to the ambient safety maintenance activities during the TSF expansion.

Discussion of the SPH Numerical Results
As presented in Section 3.2, the hypothetical tailings dam breach ROOF on the downstream complex terrain was simulated using UAV photogrammetry and an SPH numerical method. The concerned G2 highway was considered to be at a safety distance from the TSF due to the accumulation of the front flow at several low-lying spots and the blocking of both the surrounding earthworks and the water recovery ponds near the dam toe. At 600 s, both the front flow velocity and the ROOF submerged depth approaching the G2 highway significantly reduced. In the engineering practice aspect of the TSF disaster prevention and mitigation, the retaining dam or diversion gully built on the path of ROOF is important to minimize the consequences. Therefore, the simulation of hypothetical tailings dam breach ROOF can provide evidence for designing such prevention and mitigation structures. On the other hand, it can be seen that the ROOF routings and submerged area were greatly influenced by the downstream terrain as its gradient variations. In our previous study [16], the procedure of integrating spaceborne remote sensing DSM (approximately 30 m*30 m resolution) and the SPH method was introduced to predict tailings dam breach routings over downstream terrain. In the present study, the precise and fine timeliness terrain represented by the UAV photogrammetry generated high-resolution DSM is then closer to the real terrain and expected to produce more reliable ROOF simulation results.
Since the amount of disposal tailings will continue to grow for the foreseeable future, and major tailings dam breach accidents are frequently occurring since the 21st century, TSF safety still ranks as the major concern in contemporary mining [16]. Taking China as an example, despite the ongoing efforts of underground tailings backfill mining and utilization actions, the annual tailings' disposal amount is still more than 1.5 billion tons. The documented TSF amount is up to 8869 and among them 1425 were designated as the high-risk TSF, which locate within only 1 km upstream of communities, industrial plants, or other important facilities. The prediction of tailings dam breach ROOF is then of great significance for the stakeholders and administrators to implement better emergency response plans (such as efficient evacuation routes, rational deployments of emergency shelters and rescue equipment, important facilities protection measures, etc.).

Significance of the Study
Numerical simulation with the advantages of flexible settings, reasonable costs, and high efficiency can be the ideal method for the risk assessment of the devastating tailings dam breach accidents. In order to solve the numerical difficulties when modelling large deformation problems using traditional grid methods, the advanced mesh-free SPH method was used in this study. The proposed multidisciplinary procedure was proven effective by a hypothetical tailings dam breach case study in Shandong Province, China. The proposed procedure is then recommended for the risk assessment of the high-risk TSFs.
In addition, UAV SfM photogrammetry is deemed to be a reliable topographic survey technique and has been widely used in the surveillance of heritage, landslides, infrastructure engineering, and so on. In the present study, the reproduced orthophoto and DSM results of Yujiaquan TSF showed significant accuracy and resolution. Moreover, compared with the spaceborne DSM data, the UAV photogrammetry results have the advantages of great timeliness, a repeatable process, flexible schedule, and reasonable costs. It can be inferred that the multi-temporal UAV photogrammetry dataset could provide strong support for the TSF operation during its whole life cycle of site planning, tailings depositing, capacity expansions, site closure, greening, and reclamation.

Conclusions and Future Work
Over the past decade, catastrophic tailings dam breach accidents highlighted the importance of TSF risk assessment and emergency management. This paper presents a novel multidisciplinary procedure of using UAV photogrammetry and SPH numerical method to determine the run out overland flow of a hypothetical TDB. A case study of a hypothetical tailings dam breach in Yujiaquan TSF, Shandong Province, China was then carried out. Some key conclusions are:

1.
UAV photogrammetry is an advanced and powerful surveying technique in terms of great timeliness, repeatable process, high resolution, flexible schedule, and reasonable costs, when compared with spaceborne remote sensing. It has broad application scenarios in the mining industry, especially for the numerous small-scale mining sites worldwide, where professional surveillance equipment and topographical surveys are not affordable.

2.
The proposed procedure is presented and successfully applied in a case study in Shandong Province, China. The flow peak velocity of 22.6 m/s was identified at its dam toe at 50 s and a fan-shaped submerged area was formed at 600 s. The safety of the factory plants and workers deployed between the dam toe and the G2 highway was suggested to be primarily concerned.

3.
The proposed procedure is recommended for the risk assessment of the TSFs. The foreknowledge of the ROOF key features such as the peak velocity and submerged area is critical for the stakeholders to mitigate possible consequences and ensure that catastrophic accidents will not occur.
On the other hand, despite the applicability of the proposed method, dam breach is regarded as very complex phenomena involving a multitude of physical mechanisms, such as sedimentation, rheological property evolution, and more. The dam breach parameters and flow properties are decided by multiple factors, such as the dam breach causes (overtopping, earthquake, seepage corrosion, etc.), tailings saturation rate, damming method, and materials. In the present study, the implemented rheological properties of the ROOF were simplified as low concentration and homogeneous fluid. In addition, the dam breach width and mode were set according to the empirical formula in previous studies. Moreover, the flow interacting mechanism with downstream attachments also need to be considered in future studies.