The Stability of Tailings Dams under Dry-Wet Cycles : A Case Study in Luonan , China

Xingang Wang 1,2,* ID , Hongbin Zhan 3,* ID , Jiading Wang 2 and Ping Li 2 1 Key Laboratory of Mine Geological Hazards Mechanism and Control, Shaanxi Institute of Geo-Environment Monitoring, Xi’an 710054, China 2 State Key Laboratory of Continental Dynamics, Department of Geology, Northwest University, Xi’an 710069, China; wangjd@nwu.edu.cn (J.W.); 20175080@nwu.edu.cn (P.L.) 3 Department of Geology and Geophysics, Texas A&M University, College Station, TX 77840, USA * Correspondence: xgwang@nwu.edu.cn (X.W.); zhan@geos.tamu.edu (H.Z.)


Introduction
Tailings dams are formed as a result of the accumulation of dump slag from mining activities [1,2].In China, there are more than 12,000 tailings dams and this number is still increasing [3,4].Instability of tailings dams will result in loss of life and property and serious environmental pollution [5][6][7].For example, 277 people died from the accident of the dam breaking at the tailings pond in Xiangfen County, China on 8 September 2008 [8].
The effect of water on the stability of tailings dams has been long recognized.For example, Blight (1997) [9] analyzed the effect of osmotic static force on the stability of tailings dams and reported that the wet ground surface of the tailings ponds could produce mudflows more easily than dry ground surfaces.More recently, Rico et al. [10] made a thorough investigation and evaluation of 147 known cases of tailings dam disasters around the world and concluded that water from sources such as extreme rainfall and inappropriate water discharge is considered as the most important factor that

Engineering Overview and Test Results
The tailings dam in this study is located in Luonan County, China (see Figure 1), at an altitude of 1392 m~1595 m.The maximum vertical length and horizontal length of the tailings dam is approximately 286.2 m and 166 m.Since a residential area, a highway, and a river are located downstream, instability of the tailings dam is a great hazard.A geological survey designated to evaluate the stability of the tailings dam was conducted, and samples were collected from the rear edge of the dam crest (see Figure 1b).After transporting the samples to the laboratory, the FDJ-20 unsaturated soil direct shear instrument (see Figure 2) was used to conduct shear tests under different DW cycles.Vertical pressures varying from 50 kPa to 300 kPa were applied in the present study.The test results are shown in Table 1.Table 1 shows that the cohesion and angle of shearing resistance of soil reduce non-linearly with the dry-wet cycles.Based on the results in Table 1, the determined cohesion as well as frictional angle is plotted in Figure 3 against the corresponding dry-wet cycles.By fitting the data in Figure 3, the relationships between the cohesion as well as frictional angle with the dry-wet cycles are expressed as follows: Table 1 shows that the cohesion and angle of shearing resistance of soil reduce non-linearly with the dry-wet cycles.Based on the results in Table 1, the determined cohesion as well as frictional angle is plotted in Figure 3 against the corresponding dry-wet cycles.By fitting the data in Figure 3, the relationships between the cohesion as well as frictional angle with the dry-wet cycles are expressed as follows: Table 1 shows that the cohesion and angle of shearing resistance of soil reduce non-linearly with the dry-wet cycles.Based on the results in Table 1, the determined cohesion as well as frictional angle is plotted in Figure 3 against the corresponding dry-wet cycles.By fitting the data in Figure 3, the relationships between the cohesion as well as frictional angle with the dry-wet cycles are expressed as follows: It is obvious from the fitting curves that cohesion and angle of shearing resistance reduce exponentially with the increse of dry-wet cycles, and the correlation coefficients of the fitting data for cohesion and angle of shearing resistance are 0.96 and 0.98, respectively.It is obvious from the fitting curves that cohesion and angle of shearing resistance reduce exponentially with the increse of dry-wet cycles, and the correlation coefficients of the fitting data for cohesion and angle of shearing resistance are 0.96 and 0.98, respectively.

Method for Calculating the Tailings Dam's Phreatic Line under DW Cycles
Numerical analysis has been widely used to solve complex issues about the slope stability of tailings dams [4,15,16].Affected by rainfall infiltration, tailings recycling water, dry beach face evaporation, and other factors, the position of the phreatic line in the tailings pond suffers in various ways in the DW cycles [1,17,18].Therefore, in the numerical analysis mode, the determination of the position of the phreatic line under DW cycles plays an important role [19][20][21][22][23][24][25][26][27]. Figure 4 shows the diagram for calculating the tailings dam's phreatic line.Hi is the starter dam elevation; Ht is the tailings dam elevation; Hn is the normal reservoir water level; Ha is the corrected amendatory reservoir water level; He is the water inlet elevation; Hl is the long-term eliminating seepage elevation; Hc is the calculated water head; L is the seepage length; La is the upstream projected length of the started dam; Lb is the distance between the starter dam and the tailings dam; L0 is the length of dry beach face; Lh is the corrected length of dry beach face; a0 is the water overflow height; a1 is the intercept of the phreatic line on the y-axis; A, B, and C is the intersection of the curve at the corresponding position in the graph.
As illustrated in the Chinese Standard Safety Technical Regulations for the Tailing Pond (AQ2006-2005) [28], the tailings dam's phreatic line is calculated as follows:  Numerical analysis has been widely used to solve complex issues about the slope stability of tailings dams [4,15,16].Affected by rainfall infiltration, tailings recycling water, dry beach face evaporation, and other factors, the position of the phreatic line in the tailings pond suffers in various ways in the DW cycles [1,17,18].Therefore, in the numerical analysis mode, the determination of the position of the phreatic line under DW cycles plays an important role [19][20][21][22][23][24][25][26][27]. Figure 4 shows the diagram for calculating the tailings dam's phreatic line.
the relationships between the cohesion as well as frictional angle with the dry-wet cycles are expressed as follows:  It is obvious from the fitting curves that cohesion and angle of shearing resistance reduce exponentially with the increse of dry-wet cycles, and the correlation coefficients of the fitting data for cohesion and angle of shearing resistance are 0.96 and 0.98, respectively.

Method for Calculating the Tailings Dam's Phreatic Line Under DW Cycles
Numerical analysis has been widely used to solve complex issues about the slope stability of tailings dams [4,15,16].Affected by rainfall infiltration, tailings recycling water, dry beach face evaporation, and other factors, the position of the phreatic line in the tailings pond suffers in various ways in the DW cycles [1,17,18].Therefore, in the numerical analysis mode, the determination of the position of the phreatic line under DW cycles plays an important role [19][20][21][22][23][24][25][26][27]. Figure 4 shows the diagram for calculating the tailings dam's phreatic line.H i is the starter dam elevation; H t is the tailings dam elevation; H n is the normal reservoir water level; H a is the corrected amendatory reservoir water level; H e is the water inlet elevation; H l is the long-term eliminating seepage elevation; H c is the calculated water head; L is the seepage length; L a is the upstream projected length of the started dam; L b is the distance between the starter dam and the tailings dam; L 0 is the length of dry beach face; L h is the corrected length of dry beach face; a 0 is the water overflow height; a 1 is the intercept of the phreatic line on the y-axis; A, B, and C is the intersection of the curve at the corresponding position in the graph.
As illustrated in the Chinese Standard Safety Technical Regulations for the Tailing Pond (AQ2006-2005) [28], the tailings dam's phreatic line is calculated as follows: Water 2018, 10, 1048 5 of 11 where L h is the value when the dry beach face is mostly covered by tailings water (i.e., when the phreatic line is in a "wet" state). H where θ is the slope of dry beach face; where m is the upstream slope ratio.
In Figure 4, the equation for segment AB (red curve) of the phreatic line is expressed as The equation for segment BC (blue curve from Figure 4) of the phreatic line of is expressed as Here, β value can be obtained from the two equations and the coordinates of the point of intersection between the two curves.The phreatic line can be calculated by combining Equations ( 2)-( 9) when the tailings dam is in a "wet" state.
After discharge of mine water (recycling water) and rainfall evaporation, the value in Equation ( 2) when a small area of dry beach face is covered (i.e., when the phreatic line is in a "dry" state) changes to [29]: The phreatic line equation varies with the change of L h value.Therefore, if we substitute Equation (10) into Equations ( 3)-( 9), we get the phreatic line when the tailings dam is in a "dry" state.

Numerical Model for the Tailings Dam under DW Cycles
Finite element numerical analysis method, of which the reliability has been verified, has been widely used in the stability analysis of tailings slope [3,4,14].According to a field survey of the tailings dam shown in Figure 1, a finite element numerical model was established for the A-A' profile in Figure 1b, and a mesh element was generated and is shown Figure 5.The numerical model is 75 m in height and 286.2 m in length.The mesh consists of 6456 triangles and 13,208 nodes.It is important to point out that two phreatic lines in the blue zone were calculated from the method presented in Section 3.1.The positions data of two phreatic lines were imported into numerical model, and then the two lines shaped a hydro-fluctuation belt with a certain area [23], which is shown in the water blue area in Figure 5.
The left boundary of the model is vertically constrained, while the bottom boundary is fully constrained.The strength reduction method developed by Dawson et al. [30] was applied to calculate the tailings dam's safety factor in this study.Previous research reported that the effect of the geotechnical parameters such as the modulus of elasticity and the Poisson's ratio of sandy soil on the dam stability is negligible when using the finite element strength reduction method to calculate the dam safety factor [4,31].Hence, in the simulations performed in this study, parameters for the materials in the numerical model were assigned based on field survey data and routine geotechnical tests of tailings sand, which are shown in Table 2.The angle of shearing resistance and the cohesion of tailings sand in the hydro-fluctuation belt under different DW cycles were assigned according to data in Table 1.
Finite element numerical analysis method, of which the reliability has been verified, has been widely used in the stability analysis of tailings slope [3,4,14].According to a field survey of the tailings dam shown in Figure 1, a finite element numerical model was established for the A-A' profile in Figure 1b, and a mesh element was generated and is shown in Figure 5.The numerical model is 75 m in height and 286.2 m in length.The mesh consists of 6,456 triangles and 13,208 nodes.It is important to point out that two phreatic lines in the blue zone were calculated from the method presented in 3.1.The positions data of two phreatic lines were imported into numerical model, and then the two lines shaped a hydro-fluctuation belt with a certain area [23], which is shown in the water blue area in Figure 5.The left boundary of the model is vertically constrained, while the bottom boundary is fully constrained.The strength reduction method developed by Dawson et al. [30] was applied to calculate the tailings dam's safety factor in this study.Previous research reported that the effect of the geotechnical parameters such as the modulus of elasticity and the Poisson's ratio of sandy soil on the dam stability is negligible when using the finite element strength reduction method to calculate the dam safety factor [4,31].Hence, in the simulations performed in this study, parameters for the materials in the numerical model were assigned based on field survey data and routine geotechnical tests of tailings sand, which are shown in Table 2.The angle of shearing resistance and the cohesion of tailings sand in the hydro-fluctuation belt under different DW cycles were assigned according to data in Table 1.

Results and Discussion
The Mohr-Coulomb criterion was chosen as the constitutive model for materials shown in Figure 5.Moreover, data in Tables 1 and 2 were substituted into the finite element program for calculation.After different dry-wet cycles, the maximum shear stress, maximum displacement, and safety factor are shown in Table 3.

Results and Discussion
The Mohr-Coulomb criterion was chosen as the constitutive model for materials shown in Figure 5.Moreover, data in Tables 1 and 2 were substituted into the finite element program for calculation.After different dry-wet cycles, the maximum shear stress, maximum displacement, and safety factor are shown in Table 3.In order to get insight into the results in Table 3, the calculated results including shear stress, displacement, and safety factor were analyzed and are discussed in the following sections.

Analysis of Shear Stress Contour
The shear stress contour is one of the criteria for judging the sliding surface caused by instability of the tailings dam, and the area where the sliding surface is located is often the area where the value of shear stress changes significantly [15].The shear stress of the numerical calculation model in Figure 5 after six DW cycles is shown in Figure 6.As shown in Figure 6a, after one DW cycle, local shear stress concentration was found at the front edge of the hydro-fluctuation belt and at the toe of the tailings dam, with the maximum shear stress of 0.56 MPa.After 2, 3, and 4 DW cycles, the area of shear stress concentration gradually expanded, with a continuous shear stress concentration zone in the hydro-fluctuation zone (see Figure 6b-d As shown in Figure 6a, after one DW cycle, local shear stress concentration was found at the front edge of the hydro-fluctuation belt and at the toe of the tailings dam, with the maximum shear stress of 0.56 MPa.After 2, 3, and 4 DW cycles, the area of shear stress concentration gradually expanded, with a continuous shear stress concentration zone in the hydro-fluctuation zone (see Figure 6b-d).At that moment, the corresponding values of maximum shear stress for 2, 3, and 4 DW cycles were 1.43, 1.75, and 1.96 MPa, respectively.After 5 and 6 DW cycles, the shear stress distribution gradually developed to the top of dry beach surface, accompanied by a potential sliding surface (see Figure 6e,f).At that moment, the values of maximum shear stress were 2.07 and 2.11 MPa, respectively for 5 and 6 DW cycles.It was observed that the maximum shear stress is mainly concentrated at the front edge of the hydro-fluctuation belt and at the toe of the tailings dam.Moreover, the maximum shear stress after 2 DW cycles was 2.55 times of that after 1 DW cycle.As the number of DW cycles increases, the difference in the maximum shear stress of adjacent cycles decreases.

Analysis of Displacement
The displacement of the numerical model in Figure 5 after 6 DW cycles is shown in Figure 7. and 1.96 MPa, respectively.After 5 and 6 DW cycles, the shear stress distribution gradually developed to the top of dry beach surface, accompanied by a potential sliding surface (see Figure 6e,f).At that moment, the values of maximum shear stress were 2.07 and 2.11 MPa, respectively for 5 and 6 DW cycles.It was observed that the maximum shear stress is mainly concentrated at the front edge of the hydro-fluctuation belt and at the toe of the tailings dam.Moreover, the maximum shear stress after 2 DW cycles was 2.55 times of that after 1 DW cycle.As the number of DW cycles increases, the difference in the maximum shear stress of adjacent cycles decreases.

Analysis of Displacement
The displacement of the numerical model in Figure 5 after 6 DW cycles is shown in Figure 7.As illustrated in Figure 7a, after 1 DW cycle, the maximum displacement of the model reached 0.5 mm at the toe of the tailings dam.After 2, 3, and 4 cycles, the maximum displacement increased to 12 mm, 16 mm, and 19 mm, respectively.It is clear that the maximum displacement occurred at the toe of the tailings dam.Furthermore, the area of maximum displacement gradually enlarged as the number of DW cycle increased (see Figure 7b-d).Additionally, after 5 and 6 DW cycles, the maximum displacement was 21 mm at the toe of tailings dam and 22 mm at the front edge of the hydro-fluctuation belt (see Figure 7e,f).After 6 DW cycles, the area of maximum displacement was large and the area of maximum displacement in the hydro-fluctuation belt zone was relatively continuous (see Figure 7f).

Analysis of Safety Factor
The proposed approach consists of five major steps to determine the safety factor of tailings dam: (i) As revealed by Lei et al. [11], the strength reduction factor F s is assumed, the cohesion c and the angle of shearing resistance ϕ on the sliding surface are reduced as follows: where c e and ϕ e are respectively the cohesion and the angle of shearing resistance after reduction.
(ii) Based on the elastoplastic theory, apply boundary condition and material parameters to determine the stress and deformation of the tailings dam.(iii) For elements on the sliding surface, the sliding resistance F r and the sliding force F t are calculated as follows [30]: where N e is the number of elements; N g is the number of Gaussian integration points in each element; σ nig is the normal stress at the gth Gaussian point in the ith element; τ nig is the shear stress at the gth Gaussian point in the ith element; v ig is the control volume at the gth Gaussian point in the ith element; t is the calculated thickness of elements on the sliding surface; c ei is the reduced cohesion; and ϕ ei is the reduced angle of shearing resistance.
(iv) Calculate the new safety factor F s * by Equation ( 11).
The calculation is considered finished when Equation ( 14) is satisfied and the safety factor is determined as F * s .Otherwise, the calculation should be repeated according to Equations ( 11) and (12) until the results meet the Equation (14).Based on the safety factor calculation equations (Equations ( 11)-( 14)) contained in the finite element software, calculation was performed for the model in Figure 5.Then, the tailings dam safety factors after different DW cycles were obtained and shown in Table 3.
As listed in Table 3, the tailings dam safety factor decreased continuously with the DW cycle.Moreover, it is obvious that the safety factor decreased significantly from 1 DW cycle to 2 DW cycles, while the difference in safety factor between adjacent DW cycles became small after 2 DW cycle.The results suggest the first and second DW cycles reduce the stability of the tailings dam remarkably.

Conclusions
In the present study, a direct shear test and finite element numerical calculations were conducted on sands collected from a tailings dam in Luonan, China.The following conclusions can be drawn: (i) A method that efficiently calculates the phreatic line of the tailings dam under DW cycles was developed.It is capable of obtaining the boundaries of the hydro-fluctuation belt of the tailings dam under DW cycles so that the results of laboratory tests under DW cycles can be applied in numerical simulation and calculation in practical engineering.
(ii) The maximum shear stress of the tailings dam after two DW cycles was 2.55 times that after one DW cycle.However, as the number of DW cycle increases, the difference in maximum shear stress decreases, accompanied by the weakened damage effect of DW cycles.After 5 and 6 DW cycles, the shear stress distribution gradually extends to the top of dry beach surface, with a potential sliding surface developed.(iii) With increasing DW cycles, the maximum displacement of the tailings dam increased from 0.5 mm to 22 mm, and the area of maximum displacement expanded mainly at the toe of the tailings dam and at the front edge of the hydro-fluctuation belt.Moreover, the corresponding safety factor of tailings dam after each cycle nonlinearly reduced, and DW cycles strongly affected the safety of the tailings dam, especially after 1 and 2 DW cycles.(iv) In the numerical simulation of this study, the tailings sand under DW cycles was identified as in fully saturated condition, without considering unsaturated tailings sand that is prevalent in real engineering.Therefore, our future study should take the boundaries between saturated and unsaturated tailings sand into account so as to calculate the tailings dam stability more realistically.

Figure 2 .
Figure 2. FDJ-20 soil direct shear instrument and the shear failure surface after shearing (a) Direct shear apparatus; (b) Shear failure surface after shearing.

Figure 3 .
Figure 3. Relationships between cohesion and angle of shearing resistance with dry-wet cycles.Note: c = cohesion, = angle of shearing resistance, n = numbers of dry-wet cycles.

Figure 4 .
Figure 4. Diagram for calculating the tailings dam's phreatic line.Hi is the starter dam elevation; Ht is the tailings dam elevation; Hn is the normal reservoir water level; Ha is the corrected amendatory reservoir water level; He is the water inlet elevation; Hl is the long-term eliminating seepage elevation; Hc is the calculated water head; L is the seepage length; La is the upstream projected length of the started dam; Lb is the distance between the starter dam and the tailings dam; L0 is the length of dry beach face; Lh is the corrected length of dry beach face; a0 is the water overflow height; a1 is the intercept of the phreatic line on the y-axis; A, B, and C is the intersection of the curve at the corresponding position in the graph.

Figure 3 .
Figure 3. Relationships between cohesion and angle of shearing resistance with dry-wet cycles.Note: c = cohesion, ϕ = angle of shearing resistance, n = numbers of dry-wet cycles.

3 .
Methods for Calculating and Modeling the Tailings Dam's Phreatic Line under DW Cycles 3.1.Method for Calculating the Tailings Dam's Phreatic Line under DW Cycles

Figure 3 .
Figure 3. Relationships between cohesion and angle of shearing resistance with dry-wet cycles.Note: c = cohesion, = angle of shearing resistance, n = numbers of dry-wet cycles.

Figure 4 .
Figure 4. Diagram for calculating the tailings dam's phreatic line.Hi is the starter dam elevation; Ht is the tailings dam elevation; Hn is the normal reservoir water level; Ha is the corrected amendatory reservoir water level; He is the water inlet elevation; Hl is the long-term eliminating seepage elevation; Hc is the calculated water head; L is the seepage length; La is the upstream projected length of the started dam; Lb is the distance between the starter dam and the tailings dam; L0 is the length of dry beach face; Lh is the corrected length of dry beach face; a0 is the water overflow height; a1 is the intercept of the phreatic line on the y-axis; A, B, and C is the intersection of the curve at the corresponding position in the graph.

Figure 4 .
Figure 4. Diagram for calculating the tailings dam's phreatic line.H i is the starter dam elevation; H t is the tailings dam elevation; H n is the normal reservoir water level; H a is the corrected amendatory reservoir water level; H e is the water inlet elevation; H l is the long-term eliminating seepage elevation; H c is the calculated water head; L is the seepage length; L a is the upstream projected length of the started dam; L b is the distance between the starter dam and the tailings dam; L 0 is the length of dry beach face; L h is the corrected length of dry beach face; a 0 is the water overflow height; a 1 is the intercept of the phreatic line on the y-axis; A, B, and C is the intersection of the curve at the corresponding position in the graph.

Table 1 .
Shear strength of tailings sand after dry-wet (DW) cycles.

Table 2 .
Calculation parameters for the numerical model.
Note:  is density; E is Young's modulus;  is Poisson's ratio;  is Tensile strength;  is angle of shearing resistance; c is Cohesion.

Table 3 .
Tailings dam calculated results after DW cycles.

Table 2 .
Calculation parameters for the numerical model.

Table 3 .
Tailings dam calculated results after DW cycles.