Numerical Analysis of Slope Stability Considering Grading and Seepage Prevention

The normal operation of Yulangpei tailings reservoir is affected by landslide stability. In this paper, taking the main and side slopes near the dam bank of the Yulangpei ditch as an example, water-soil coupling theory is applied to comprehensively evaluate the reliability of the side slopes of the tailings reservoir. Grading and seepage prevention (GSP) measures and the suction of the substrate are considered, as well as the infiltration of different rainfall and reservoir water levels. We numerically simulate the typical three forms of side slopes under the coupling conditions and conduct a reliable and comprehensive evaluation of tailings reservoir side slopes. The study shows that under six reservoir water level changes, the factor of safety (FS) of the bank slope shows a hysteresis effect. According to nine rainfall infiltration conditions and during rainfall, the greater the rainfall intensity, the greater the weakening effect. When rainfall stops, the FS rebounds. After GSP measures, the initial stability of the bank slope under different conditions is improved, but the main slope is more sensitive to changes in rainfall and water levels.


Introduction
Reservoir landslide disasters are one of the most common geological hazards in tailings reservoirs (Jiao et al. 2014). During tailings reservoirs operation, long-term water level changes and heavy rainfall will often cause groundwater periodicity on the banks. This fluctuation will further affect the stability of the bank slope (Cojean and Caï 2011). These changes will also cause hydrological changes inside the slope body, which in turn will lead to changes in the saturation state of the landslide, as well as changes in the geomechanical and physical properties of the groundwater table, which will affect the stability of the bank slope (Fourniadis et al. 2007). Due to changes in reservoir water levels and heavy rainfall, landslides may occur in younger stable slopes, while ancient landslides may once again become active. For these reasons, the influence of water level changes and heavy rainfall is considered to be an important hazard factor for the stability of tailings bank slopes. Increasing evidence suggests that approximately 90% of slope instability is related to some extent by pore water pressure (PWP). Due to the periodic changes in reservoir water levels and rainfall, the landslide soil undergoes periodic soil-water effects, namely: (1) changes in water content (Acharya et al. 2015); (2) alternating wet and dry effects (Reid and Parkinson, 1984); (3) groundwater seepage (Take et al. 2015;Pender et al. 2016); and (4) hydrochemistry (Wei et al. 2012;Wen and He 2012). At present, the stability evaluation of the bank slope mainly shows a steady state, i.e., a single event corresponds to a single safety factor. Therefore, it is beneficial to study the stability of the bank slope by introducing time as a factor, so as to discuss any dynamic changes occurring in the bank slope safety over time. Traditional analysis of slope events is often a comparative analysis of excavation (Chen et al. 2018;Ran et al. 2019;Singh and Lan 2017;Kang et al. 2017), earthquakes (Zhou et al. 2019;Javdanian and Pradhan 2018;Gomberg 2018), seismic activity , rainfall (Vanwoert et al. 2005;Zhuang et al. 2018;Zhao et al. 2018;Fu 2017), reservoir water level (Mandal et al. 2019;Sun et al. 2015), and freeze-thaw change effects (Fan et al. 2014;Cheng et al. 2017). Comparative analysis performed both before and after GSP is relatively rare. However, reservoir water level changes and heavy rainfall comparative stability analysis on bank slopes, before and after GSP, has been conducted more often. Previous studies have typically F. Yan • C. Xia () • G. Lu School of Geosciences and Info-Physics, Central South University, Changsha 410083, China e-mail: 2110409@tongji.edu.cn F. Yan e-mail: yanfuyu@csu.edu.cn G. Lu e-mail:luguangyin@csu.edu.cn yingagishi only considered steady-state saturated flow, and not the influence of matrix suction on landslide stability (Wei et al. 2008). The groundwater level is generally below the surface, while the PWP of the soil above the groundwater level is negative, i.e., suction occurs. This suction will increase the stability of the slope on unsaturated soils (Xiong et al. 2014). Kitamura and Sako (2010) reviewed studies of slope failure caused by rainfall over the past 50 years and emphasized the importance of unsaturated soil mechanics in elucidating the failure mechanism. In general, there are two physical mechanisms that cause soil suction stress, namely, the physical and chemical forces between particles, and capillary forces between particles. When the soil is near saturation, the suction stress is greatly reduced. This phenomenon may be the physical mechanism that instigates many shallow landslides under heavy rainfall conditions (Lu and Godt 2008). Existing research has not paid sufficient attention to the dynamic changes of soil pressure distribution , which is a common occurrence. Water-soil coupling is of great significance in landslide stability assessments.
In this study, the main and side slopes of the Yulangpei Gully, near the dam bank, are taken as an example. Based on the Morgenstern-Price method, using the theories of water-soil coupling and unsaturated seepage under the infiltration of different rainfalls and reservoir water level changes, a numerical simulation of the bank slope stability was performed under the coupling conditions. In addition, the factor of safety (FS), saturation distribution, and displacement of the slope with time, before and after grading and seepage prevention (GSP), of the bank slope were obtained. Finally, we performed a comprehensive evaluation on the stability of the bank slope.

Establishment of the control differential equation
This paper uses the SEEP/W module in Geo-studio to simulate seepage processes. In this study, the saturated-unsaturated seepage theory is used for numerical simulation, where rainfall infiltration is considered as a variable flow boundary in the unsaturated zone. According to the principle of mass conservation, the saturated-unsaturated seepage control differential equation is as follows (SEEP/W 2007): where h is the pressure head; ( ) r k h is the relative water permeability coefficient; ( ) 1 r k h  in the saturated region, and ( ) [0,1) r k h  in the unsaturated region; ij k is the saturated permeability tensor; s S is the unit water storage coefficient; ( ) C h is the water capacity; if in the positive pressure zone, ( ) 0 C h  ; if in the negative pressure zone, ( ) / C h h     ;  is the moisture content;  is the parameter for determining the saturated and unsaturated states; in the saturated region, 0   ; in the unsaturated region, 1   ; t is time; and S is the source exchange item.

Establishment of the unsaturated permeability coefficient
It can typically be expressed by the fitting equation in the Van Genuchten model (Genuchten 1980), which describes the water characteristic curve and unsaturated permeability coefficient: where  is the corrected volumetric water content; r  is the residual volumetric water content; s  is the saturated volumetric water content; s is the corrected saturation; r s is the residual saturation level; and a , n and m are the fitting parameters. According to the relationship between volume moisture content and saturation, Equation (3) can be obtained from Equation (2). According to the relationship of relative permeability and saturation, from Equation (4), Equation (5) can be obtained, where w k is the saturated permeability coefficient, r k is the relative permeability coefficient, and k is the adjusted coefficient of permeability.

Establishment of the finite element seepage equation
By applying the Galerkin method with a weighted margin to the governing equation, we can obtain the finite element format of the two-dimensional seepage equation:  is the interpolation function vector; q is the unit flow across the element boundary;  is the element thickness; t is time;  is the storage term; w w m  are equal to the transient downstream current; A is the summation sign on the element area; and L is the summation sign on the element boundary length.
In the axisymmetric analysis, the equivalent element thickness is the hoop distance at the radius R , which is relative to the axis of symmetry. The complete hoop distance is multiplied by 2 . Due to the fact that SEEP / W in Geo-studio is derived for a 1-radian surface, the equivalent thickness is R . The finite element equation in the case of axisymmetric is as follows: It is important to note that, unlike the thickness in the two-dimensional analysis, the distance of radial R in a unit is not a constant, and that the radial R changes within the integrand. The finite element seepage equation can be expressed in the following simplified form: is the element mass matrix; and { } Q is the flow vector applied on the element. The above equation is a general-purpose finite element equation used for transient seepage analysis in Geo-studio SEEP / W.
2.4 Establishment of the factor of safety for unsaturated soil and rock SLOPE/W uses the Morgenstern-Price method, based on limit equilibrium theory to calculate the FS. The modified method strictly satisfies the force balance and torque balance and is highly accurate. The expression is shown below: where i c is the cohesive strength for every soil slice; i is the soil slice number; i W is the weight of each soil slice; i P is the water pressure; i  is the angle of the bottom of the soil slice; i b is the length of each soil slice; i  is the frictional strength for each soil slice; i r is the radius of the sliding arc; and s F is the FS.

Establishment of the numerical calculation model
As shown in Fig. 1, the bank slope of the Pulang copper mine tailings bank is located to the northeast of Diqing Tibetan Autonomous Prefecture, in northwestern Yunnan Province, China. The Pulang copper mine is one of the largest underground copper mines in the country. The bank slope selected by the geological model is the bank slope near the rockfill dam of the Yulangpei tailings reservoir, as shown in Fig. 1c. The coordinate system of the numerical calculation model was selected as follows: The X (538 m) direction of the main slope is perpendicular to the surface runoff, and the Y (242 m) direction runs vertically upwards, as shown in Fig. 2a. The X (328 m) direction of the side slope is almost parallel to the surface runoff, while the Y (196 m) direction runs vertically upwards, as shown in Fig. 2c. The sizes of the slope and model after grading and seepage prevention (GSP) remain unchanged. The height of the main slope grading is about 165 m, and the length is about 290 m, as shown in Fig. 2b. The height and length of the side slope grading are roughly 135 and 240 m, respectively, as shown in Fig. 2d. In the numerical simulation, there are 59 monitoring points on the main slope and 66 monitoring points on the side slope. The thickness of the geomembrane is only 0.5-2 mm, and the actual permeability coefficient value is 10-13cm/s. To simulate the geomembrane, according to the equivalent principle of seepage, the total volume of seepage flow is constant, therefore, the equivalent permeability coefficient k  and the equivalent thickness tcan be expressed as: According to the research of Cen et al. (2012), when the equivalent thickness is 1000 mm (1m), the equivalent permeability coefficient is 10-13 cm/s, and the numerical simulation effect is optimized. For the HDPE geomembrane, the height and length of the impervious layer on the main slope are 125 and 220 m, respectively, as shown in Fig. 2c. The height and length of the impervious layer on the side slope are 75 and 130 m, respectively, as shown in Fig. 2d. The rock and soil constitutive model follows the Mohr-Coulomb model, and the FS is calculated using the Morgenstern-Price method. As shown in Table 1, the parameters for ascertaining the rock and soil layers of the bank slope are based on quality surveys and laboratory test data. Combined with the VG model parameter sensitivity analysis and literature (Chen et al. 2020), the parameters of the water retention curve of the rock mass can be attained by referring to Table 2. The relationship between the volumetric water content, permeability coefficient, and matrix suction in the soil-water characteristic curve is determined by geological survey data and Equations (2) -(5), as shown in Fig. 3.  Based on hydrogeological data, the initial water level at the far bank of the main slope is maintained at approximately 210 m, and the near-shore groundwater head is roughly 30 m. The initial water level at the side slope is fixed at roughly 150 m. The near-shore end groundwater head is approximately 40 m high, allowing the model to seep naturally for 600 days, and therefore the seepage-to-matrix suction is stable. The simulation time is 120 days. The numerical simulation considers three dynamic boundary conditions and their time-varying properties, which are defined as the three following operating conditions: Condition 1: Different reservoir water level height and assorted rising and falling rate combinations. It is assumed that the rising and falling water levels are linear, while the water level changes are generalized into three relative modes. Two water level heights are considered, and category 1 is fast rising and slow falling (water levels 1 and 4); category 2 is moderate rising and falling (water levels 2 and 5); category 3 is slow rising and fast falling (water level 3 and 6). There are six cases in each category, as shown in Fig. 4.
Condition 3: According to the real-time rainfall and monitored reservoir water level from 24:00 on March 1, 2019, to 24:00 on July 27, 2019, the influence of real rainfall and reservoir water level coupling on bank slope stability can be simulated, as shown in Fig. 5.

Analysis of condition 1
Figs. 6 and 7 show the calculation results for the factor of safety of the bank slope under working condition 1. For the different reservoir water level fluctuation modes, when the increase is small, whether the reservoir water level rises or falls, the valley of the safety factor (VFS) appears only once. Therefore, there is no factor of safety peak (PFS). When the increase is large, the VFS appears twice during the rise and fall of the reservoir water level, thus a PFS is present. The PFS of the three water level modes increases with the water level rising faster and higher.
The rise and fall of the bank slope FS is associated with the rising water level speed and height. When the water level rising height is larger, and the water level rises faster. the FS keeps increasing. This is because the pore water pressure growth rate in the slope is less than that of the reservoir water pressure on the main slope (Wu et al., 2017). This results in reverse pressure on the foot of the slope, and the difficulty of the slope body from saturation to saturation will gradually increase. The negative pressure zone decelerates, resulting in the bank slope FS increasing. The FS is also related to the decline rate of the reservoir water level. The faster the reservoir water level decreases, the faster the FS decreases. This is because the decline in the slope water level lags behind the decline in the reservoir water level. In addition, the pore water pressure decrease rate on the slope is less than the decrease rate of the water pressure on the main slope. Consequently, excess pore water pressure and hydrodynamic pressure are generated outside the inclined slope. The FS exhibits a lag effect with the water level change, that is, the FS corresponding to the highest water level is not the minimum value, yet the lag reaches a minimum after a certain period . This lag effect is further strengthened by the faster increase of water level, a slower decline, and a greater increase in water level. The water level 1 lag effect in mode 1 is the most apparent, and the water level 6 lag effect in mode 3 is the weakest. For the various reservoir water level fluctuation modes, during the water level rise of the reservoir main slope, the VFS appeared last in mode 3 (water levels 3 and 6), and first in mode 1 (water levels 1 and 4), as shown in Fig. 6.
Following GSP, the stability of the bank slope improved, and the FS change law remained basically the same as before GSP. Furthermore, the side slope lag effect was enhanced, but that of the main slope was weakened. The main slope VFS appeared later. The side slope FS was clearly improved, yet the main slope FS increase was not significant. The main slope PFS decreased, and that of the side slope increased significantly. show the progressive results of the slope body saturation and horizontal displacement at water level 1. Initially, the saturation above the water level surface is inversely proportional to the elevation, namely the further away from the water level surface the rock and soil body is, the less saturated it will be. At this time, the horizontal displacement of the slope is zero. After 20 days of rising water level, the infiltration line in the slope increases the slope internal saturation, and the rate of increase in the water level becomes greater than the slope infiltration rate. This causes the saturation curve to become "concave". The soil shear strength is low, and the horizontal displacement is large under the soil-water coupling. After 60 days, the water level drops. As the infiltration line in the slope decreases, the slope saturation decreases, and the slope seeps outwards, which in turn generates dynamic water pressure and excess pore water pressure. Consequently, the slope horizontal displacement further increases. After another 40 days, the water level drops again. Although the water level has returned to normal, the soil was rich in water and the saturation is still higher than its initial state. At this time, the displacement increases slowly. After grading, the initial slope unsaturated area decreases. Compared with pre-grading, the slope saturation and horizontal displacement are more sensitive to water level changes, the horizontal displacement area increases, and the maximum displacement area approaches the foot of the slope. Figs. 10-11 show the calculated result of the bank slope FS in condition 2. The stronger the rainfall intensity was, the more obvious the weakening effect is. After the rainfall ended, the FS rebounded. Following a heavy rainfall period, the slope FS, before GSP, decreased twice. The first drop was during rainfall. Under short-term and high-intensity rainfall, a strong outward tilting water pressure was generated inside the slope, which caused the slope stability to decrease . The second drop occurred after the rain stopped. The mine drainage lagged behind, and a large amount of water was stored in the slope. After the rainfall ceased, the drainage increased, the slope continued to seep and produced a strong outward water pressure, causing a decrease in the slope stability.
According to Article 5.3.1 of the Technical Specification for Construction Slope Engineering (GB50330-2013), the stability coefficient of the bank slope is 1.3, as shown in Fig. 10b. After GSP, the initial and final factors of safety of the main slope were enhanced. The sensitivity to rainfall was also improved, but during rainfall, with an intensity greater than 7 mm/h, the FS was less than 1.3, indicating that the main slope seepage protection was poor. After GSP, the initial and final side slopes safety coefficients were significantly improved. When rainfall intensity was no greater than 9 mm/h, the sensitivity to rainfall was reduced. When the rainfall intensity exceeded 13 mm/h, the FS was lower than previously, but remained in a stable state, as shown in Fig. 11b. Considering that the record for single-day maximum rainfall intensity in the Prang area is 8 mm/h, the side slope GSP can be considered safe and reliable. Figs. 12 and 13 show the progressive results of the saturation and horizontal displacement of the main slope under the rainfall intensity of 17 mm/h. After 5 days of heavy rainfall, the saturation of the slope surface increased sharply, and the horizontal displacement increased significantly. After 3 days without rain, the surface water infiltrated the slope, the slope surface saturation decreased, and the horizontal displacement also gradually decreased. After 5 days with no rain, the surface water continued to seep into the slope, and both the saturation zone and the horizontal displacement slowly increased. After grading and seepage, the main slope impervious layer slowed the water outflow and seepage of the slope, and the middle and lower parts of the slope drew close to saturation. The horizontal displacement area instigated by heavy rainfall increased. The smaller the monitoring point value, the higher the altitude. For the monitoring points above the highest point of the reservoir water level, the PWP is affected by changes in the reservoir water level. As the distance from the monitoring point to the highest water level point increases, the PWP gradually becomes smaller eventually reaching a fixed value. At monitoring points below the highest reservoir water level, the PWP is consistent with the changes in the reservoir water level. Furthermore, the closer the monitoring point is to the bottom of the slope, the greater the change in PWP will be. After GSP, a negative pressure zone appeared in the middle of the main slope, and the negative pressure zone increased in the lower and middle sides of the side slope, which is indicated by the red dotted area. Fig. 15 shows the calculated results of the main slope FS and displacement over time. It can be seen that two stages exist in the change of the FS from March 1 to July 27. The overall change trend is a small decline, followed by a sharp rise. This is correlated with the increased rainfall intensity in the rainy season, more concentrated rainfall, faster water level elevation, and greater water level increase. The displacement is the total displacement of each monitoring point, which increases steadily. When the rainfall intensity is greater and the water level rises faster, then the displacement growth rate increases. After GSP, the overall FS is increased by about 3%, and the total displacement is reduced by about 87%. Fig. 15a shows the FS and displacement of the side slope of working condition 3 with time. From March 1 to July 27, the safety factor was negatively correlated with the change in total displacement. The overall trend is an initial slight rise, followed by a sharp rise. This is related to the increased rainfall intensity during the rainy season, as well as more concentrated rainfall, faster water level elevation, and greater water level increases. After GSP, the overall safety factor of the side slope increased by approximately 41%, and the total displacement decreased by about 66%. When the rainfall intensity was greater than 4 mm/h the bank slope total displacement increased significantly.