Numerical Analysis of Slope Stability under Reservoir Water Level Fluctuations Using a FEM-LEM-Combined Method

Large-scale slopes at the banks of reservoirs pose a serious threat to the safety of hydropower stations. The fluctuation of the reservoir water level is a key factor in the slope stability. However, the parameters to describe the relationship among water content, matric suction, and soil strength are difficult to measure using unsaturated soil strength theory. To solve this problem, a simple FEM-LEM-combined scheme considering pore pressure, seepage force, and strength weakening is presented to calculate the safety factor. A numerical study on the impact of reservoir water level fluctuations on stability of a glaciofluvial deposit slope is implemented. Two typical profiles are used to estimate the stability of the glaciofluvial deposit slope in response to rising and lowering water levels. The results indicate that this method proposed a simple and efficient tool for water level-induced slope stability analysis.


Introduction
Compared with general landslides, landslides occurring at the banks of reservoirs are unique, in that their occurrence is correlated to fluctuations in reservoir water levels, wave erosion effects, and water immersion [1][2][3]. Data indicates nearly 50% of reservoir-induced slope failures occur during the first impoundment, while slope failures mainly occurred within 3-5 years of a dam's construction [4,5]. Many studies have demonstrated that fluctuating water levels are one of the essential causes of reservoir landslides [6][7][8][9][10].
Thus, from the perspective of hydropower engineering, the potential of reservoir landslides plays a determining role in dam site selection, most especially in the wake of the catastrophic Vaiont slope failure [11]. During the past several decades, this potentiality has routinely been consid-ered in the design and construction of hydropower stations throughout southwest China, a region abundant in hydropower resources, owing to the presence of ideal geological phenomena: great rivers, high mountains, deep and narrow river valleys, and ample water drops in river beds [12,13]. However, previous studies have indicated this region is particularly prone to landslides [14][15][16][17][18][19][20][21]. Some case studies have already analyzed the impact of unstable reservoir slopes and landslides on the construction of hydropower stations [2,17,[22][23][24][25]. Therefore, the slope stability analysis under reservoir water level fluctuations is of high significance for the safety of hydropower engineering.
In the analysis of slope stability with water level fluctuation, unsaturated soil strength theory is often used. However, the parameters to describe the relationship among water content, matric suction, and soil strength are difficult to measure [5,[26][27][28][29][30][31][32][33]. Mostly, the project design and construction unit just give the mechanical parameters for the geomaterials in the natural and saturated states. For most widely used codes, the unsaturated theory is used, and the slope stability in this situation cannot be calculated directly. To solve this problem, a simple FEM-LEM coupling method for slope stability under reservoir water level fluctuations is proposed. Two profiles of a glaciofluvial deposits are analyzed in the impoundment and drawdown condition to indicate the ability of this method.

FEM-LEM-Combined Analysis Method
2.1. Reservoir Water Impact Mechanism. Before the numerical modeling, a mechanism analysis on the instabilities of the glaciofluvial deposit slope is conducted in the conditions of the reservoir water level fluctuations. Generally speaking, in this problem, there are largely three significant factors affecting the slope's stability, namely, interspace pore pressure, seepage force, and strength weakening of materials.
2.1.1. Influences of Pore Pressure. In order to analyze the mechanism of the influence of the water level on the slope stress states, we assume two extreme cases for a homogeneous slope, which are completely above and below the reservoir water level. For simplicity and clarity, the strength parameters of the slope geomaterials in natural states and wet states are assumed to be the same. To begin with, the slope is assumed to be in a natural state. The stress state of a random point in the slope is schematically shown in Figure 1(a).
The balance partial differential equations are as follows: where γ is the natural volume weight of the slope's geomaterial and α is the dip angle of the slope.
In fact, this is a classical problem and the solution has been presented by Chen [34]. Therefore, for simplicity, this paper gives the results directly: where σ 1 and σ 3 are the major and the minor principal stress at the point, respectively. λ is a variable for representing the lateral soil pressure coefficient, which relates σ x and σ y , in the equation σ x = λσ y . The variable h means a vertical distance from the slope surface to the point. Thus, the principal stresses at the point can be simplified as where k 1 = ðγ/2Þððλ + 1Þ cos 2 α + ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðλ − 1Þ 2 cos 4 α + sin 2 2α q Þ and k 3 = ðγ/2Þððλ + 1Þ cos 2 α − ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðλ − 1Þ 2 cos 4 α + sin 2 2α q Þ with the relationship k 1 > k 3 . If the point is located on a slice in the state of critical instability, according to the Mohr-Coulomb criterion, the strength curve will lie tangent to the Mohr stress circle and the cohesion can written as where c and φ are the geomaterial cohesion and internal friction angle, respectively. Substituting Equations (3) and (4) into Equation (5), the critical value of cohesion c can be expressed as follows: As the whole slope is submerged into water, the volume weight will be changed to the floating one γ ′ . Correspondingly, the principal stress σ 1 and σ 3 will decrease to be σ 1 ′  2 Geofluids and σ 3 ′ , and the reductions are k 1 hðγ − γ ′ Þ and k 3 hðγ − γ ′ Þ, respectively, which can be seen in Figure 1(b). Just like the analysis for natural unsaturated conditions, a similar relationship is proposed for saturated states.
where γ ′ and c ′ are the floating volume weight and cohesion of the geomaterial, respectively. Here, it is assumed that the geomaterial has common internal friction angles in both natural and saturated conditions. Comparing the expression of c with that of c′, it can be reasonably inferred that a slope will be safe in the impoundment situation and be dangerous in the drawdown circumstance.
2.1.2. Influence of Seepage Force. A concept slope for glaciofluvial deposit is presented in Figure 2(a). According to the flow net, the average hydraulic gradient of each grid is set to be j i and the related seepage force can be expressed as follows.
Here, γ w means volume weight of water, α i represents the area of the i-th grid, and j i represents the hydraulic gradient.
Considering the moment arm for J i is d i , the corresponding moment can be written as γ w j i α i d i . Assuming there are several grids in a slice (Figures 2(b)-2(c)) and the sliding force caused by the seepage field is distributed on slip surface uniformly, thus, the value of this force is where L s is the arc length of the slice at the sliding surface.
When the slope is in the condition of impoundment, the influence of the seepage force on the slope stability is in favor of increasing the safety factor. Otherwise, it will increase the sliding fore and reduce the slope safety.

Influence of Strength
Weakening. The glaciofluvial deposit slice can be divided into two categories on the sliding surface of CD and DA, and the number of each kind of slices is m and n. Neglecting the force in between the slices, the formula of the safety factor can be written as follows: where c and φ are, respectively, the cohesion and friction angle of the glaciofluvial deposit in the natural condition and c ′ and φ ′ are the effective parameters in the saturated condition.
With swift decreases in water levels, the strength parameters cannot recover to their natural condition, and the effects of seepage force as well as pore pressure work to increase slipping force. Thus, according to the analysis of Equation (10),     4 Geofluids the stability of glaciofluvial deposit slope tends to drop markedly. However, if the water level is low for a long time, the mechanical parameter may recover which increases the stability of the slope. In contrast, during impoundment, the wetting area will be changed soon if the permeability is high, and it will be bad for slope stability.

FEM-LEM-Combined
Scheme. The analysis above indicates that the glaciofluvial deposit slope subject to the water level fluctuations is a complex system and its stability is controlled by three main factors: interspace pore pressure, seepage force, and weakening degree of strength. In traditional limit equilibrium method (LEM) software, like Slide and GEOSLOPE, these factors cannot be taken into consideration simultaneously.
To solve this problem, we proposed a novel FEM-LEMcombined numerical scheme considering these three factors, as shown in Figure 3, for estimating the stability of the glaciofluvial deposit slope. The main idea of this method can be composed of three parts: first, calculating the seepage of the glaciofluvial deposit slope; second, changing the material partition of the model based on the result of seepage field and assigning weakened parameters for the wetting region; and finally, estimating the safety factor of the slope with material repartition. In this way, FEM is employed for the analysis of seepage field and LEM is employed for safety factor analysis based on the result of FEM.

Geofluids
For the calculation of the safety factor of the glaciofluvial deposit slope, the problem is divided into two conditions. As the reservoir water level rises, first, based on the hydraulic boundary condition, the seepage field in the slope is determined. Then, according to the seepage field, the strength parameters of the glaciofluvial deposits under infiltration lines are changed. On the basis of the seepage field calculated with a finite element program, geomaterial zone repartitions are reassigned new mechanical values, then the Morgenstern-Price method is applied to compute the safety factor of the glaciofluvial deposit slope. In the rapid drawdown condition, considering that the strength parameters of the glaciofluvial deposit cannot recover suddenly, their mechanical properties are assumed to be constant in the subsequent calculation.
As for the implementation of the FEM-LEM-combined method, GeoStudio software is employed in this work. Two modules, namely, Slope/W and Seep/W, are used alternately. Seepage analysis is firstly carried out to get the groundwater saturation line. Then, the mechanical property is reassigned based on this line, and the model is transferred to the Slope/W module. The limit equilibrium analysis is implemented easily. In this way, making the calculation of the dynamic seepage field and combined with the limit equilibrium analysis, the stabilities of glaciofluvial deposit slope under various reservoir water levels can be evaluated.

Geological Condition.
A glaciofluvial deposit slope located on the left bank of the Gushui Reservoir is employed in this study. This slope is the largest and requires engineering treatment, insofar as it could potentially threaten the construction and future operations of the proposed hydropower plant. The said slope is close to the dam, within a distance of about 300 m from its downstream boundary to the left abutment.
Presented in Figure 4 is a contour map showing the boundaries, topography, and geological features of the studied glaciofluvial deposit slope. From this, one can see that the Meilishi-Cigu fault passes through the lower part of the glaciofluvial deposit and the Lancang river, trending northeast to southwest. Geological surveys confirmed this fault was the largest running beneath the slope. In the lower and upper parts of the slope, several exploratory audits, named PD15 and PD17, were drilled into the mountain. It was ascertained that the glaciofluvial deposit is approximately 50 meters thick with a volume of 700 × 10 4 m 3 .
For the sake of detailed analysis, two engineering geological profiles I and II, extracted from the contour map (Figure 4), are shown in Figure 5. There is a wide distribution of the deposit in these two profiles, and the slope stability can represent the whole model. It can be seen that the glaciofluvial deposit is composed of layers of quaternary diluvium and glaciofluvial soil-stone-blended medium. Under these two layers lies a layer of toppling rock mass, which stretches continuously from the upper reaches of the reservoir to the front of the dam and trends northeast to southwest, toppling towards the river. The green imaginary curves, seen in Figure 4, show the three obvious beddings, which correspond to the three thick layers of the toppling failure rock masses shown in Figure 5. In geological surveys of the dam, a large area of glaciofluvial deposits covering a thick layer of toppling rock masses can be found, isolated not just to the study slope. To some extent, these signs explain some of the region's evolutionary features in the formations of the Lancang river valleys.
It is interesting to find that the toppling deformation is rather serious in the contact zones between the bottom of the glaciofluvial deposits and the top of the toppling rock masses. At greater depths, the toppling deformational phenomenon of the underlying rock masses becomes less obvious. Although toppling deformation in the contact zones is evident, geological surveys indicate that the toppling rock masses remain stable. This paper mainly focuses on the stability analysis of the glaciofluvial deposit slope that lies above them.
The study area is a region with a mountain monsoon climate, characterized by low air temperature, low rainfall, low humidity, evaporation, and a long frost period. The annual mean temperature was about -2.6 to 12.0°C, and annual average precipitation is approximately 494.7 to 522.4 mm. The lack of vegetation on both sides of the banks of the Lancang   6 Geofluids river, as indicated by Figure 6, highlights the fragility of the study area.

Numerical Models and Parameters.
For the purposes of analyzing and evaluating the stabilities of the studied glaciofluvial deposit slope, a numerical modeling study of the slope's dynamic seepage process in response to water level fluctuations has been performed. The study's aim is to provide appropriate methods for increasing the stability of the slope. The study models a dynamic cyclic process, involving water seepage into the glaciofluvial deposits during increases in the reservoir's water levels, as well as water seeping out from the glaciofluvial deposits during periods where the water levels decrease as the reservoir drains. This is a complex seepage-mechanical coupling problem and very common with reservoirs associated with hydropower stations. In the following study, a lot of emphasis has been placed on the numerical simulation of the seepage-mechanical coupling process and the failure mechanism, as well as the stability of the glaciofluvial deposit slope. Both of the two typical profiles are chosen to be the geological models for analyzing the stabilities of the slope, which are shown in Figure 7. In this study, the analysis emphasizes the stabilities of the glaciofluvial deposits, so the underlying toppling rock masses and bedrocks are simply considered as one homogeneous material.
For seepage analysis, it is assumed that if the elevation of the point is above the reservoir water level, the flux on the point is zero. Other far-field boundaries are assumed to be nonpermeable. To analyze the influences of the water level fluctuations on the safety factor, the water level elevations are assumed to range from 2230 to 2265 m.
The periods of reservoir water drawdown and impoundment are assumed to last 14 days, implying the corresponding velocity is 2.5 m a day. It should be noted, however, that this is a rather extreme condition and rarely encountered in reality. A part of the engineering geology report [35] indicates that the altitude of the far-field underground water level is assumed to be 2235 m.
Compared with simple soils, special structures of the glaciofluvial deposits have rather different mechanical and hydraulic properties [25,[36][37][38]. Because of the large number of stones and gravels contained in the glaciofluvial deposits, they are more permeable, and thus their permeability parameters range more. Due to the loosened structural states [10,39,40], the glaciofluvial deposits have greater saturation and less evident unsaturation above the underground water table.
Thus, soil-water characteristic curves for the glaciofluvial deposit slope cannot be as easily obtained as those for soil. For simplicity's sake, the glaciofluvial deposits are assumed to be a kind of saturated-only material. All the materials in the two models are simplified into two kinds, namely, rock mass and glaciofluvial deposit. Their conductivity coefficients are 3 × 10 −3 cm/s and 5 × 10 −5 cm/s, respectively, according to estimates obtained from geological surveys conducted by HydroChina Kunming Engineering Corporation [41]. Correspondingly, the parameters for water tolerance are 0.76 and 0.53, which are given in Table 1.
In terms of stability analysis of the reservoir bank slopes, water erosion and degradation effects on the geomaterial mechanical properties should be taken into account. Some seepage-mechanical tests [26,42] on the soil-stone-blended medium or glaciofluvial deposits indicate that their mechanical properties when saturated with water deviate from their natural conditions. These effects unavoidably impair the stability of the glaciofluvial deposit slopes.
At the same time, water seeping into the pores lead to saturation of the glaciofluvial deposits, reducing their gravity with buoyancy. This effect, to a degree, can greatly decrease the slide moment along the slip surfaces and are conducive  to slopes' stabilities. Based on this qualitative analysis, one can see that this study on the slope's stability subject to fluctuating water levels is complicated and dependent on the model's boundaries and the associated mechanical parameters, as well as the solution strategy. According to several field and laboratory test results, the physical and mechanical parameters of the glaciofluvial deposits and underlying toppling rock masses are listed in Table 2. These are used in the following numerical simulations.

Numerical Simulation on Deposit Slope
Glaciofluvial deposit slopes have two kinds of potential failure modes, namely, local and overall failures, respectively, shown in Figures 8(a) and 8(b). In the case of local failures, a part of the slope submerged underwater creates a local slide into the reservoir. In the case of overall failure, a local slide occurs underwater, triggering a landslide above water. Generally speaking, overall failure involves considerably more slide masses into the reservoir, proportionally causing greater damage. In the realization, the local failure surface is that with a minimum safety factor obtained using a search algorithm. The overall failure mode is computed using a specified slide surface. Hence, two safety factors for the glaciofluvial deposit slope corresponding to the different failure modes must be calculated.

Slope Stabilities of the Impoundment Condition.
During the period of impoundment, there are three factors affecting the slope's stabilities: strength of the glaciofluvial deposits softened by water, variations of the pore pressure, and seepage force through the slope masses. Therefore, the schemes proposed in Section 2 are employed for numerical calculation.
Firstly, the dynamic seepage field varying with water level fluctuations is calculated based on the finite element method, with the analysis results in the form of the seepage lines at different times presented in Figure 9. Then, the material partitions of the glaciofluvial deposits change and are reassigned with new values at every step varying with the dynamic seepage field. Finally, after the repartition of the slope materials, the safety factors are computed by the Morgenstern-Price method.
The safety factors for profiles I-I and II-II in the fast impoundment condition are computed and listed in Table 3. For profile I-I, when the water level changes from 2230 m to 2265 m, the local safety factor increases from 1.106 to 1.163, while the overall safety factors (Fs) decrease from 1.204 of 1.191. They show opposite changing trends. In terms of profile II-II, the overall safety factor has a drop from 1.196 to 1.184, and the local safety factor falls from 1.169 to 1.135.
The changing trends of the two profiles for the glaciofluvial deposit slope are presented in Figure 10 clearly. In general, the safety factor of the glaciofluvial deposit slope is always larger than 1.0, which means that the slope can keep stable during the reservoir storage. In spite of the fact that the slope can keep a bit stable during the fast impoundment stage, increases in the risk of overall failures cannot be neglected. A necessary support measure is rather needed to stabilize the slope in this condition.

Slope Stabilities of the Drawdown Condition.
In terms of the rapid drawdown condition, it is assumed that the mechanical parameters of the glaciofluvial deposits submerged underwater are still constant and cannot recover as quickly as the descent of the water level. Thus, in the rapid drawdown conditions of the reservoir water level, the   Figure 11. After calculation of the changing stabilities of the glaciofluvial deposit slope with dynamic seepage field, the safety factors with time are summarized in Table 4. In the rapid drawdown condition, both of the local and overall safety factors have a downward trend. When the reservoir water level is at the altitude of 2265 m, the safety factors of the overall failure mode for profiles I-I and II-II are 1.193 and 1.186, respectively. After a 7-day rapid drawdown, the altitude of the water level drops to 2230 m in elevation. Then, the overall safety factors at profiles I-I and II-II change to 1.185 and 1.175, respectively.
As for another failure mode, the safety factor of profile I-I drops evidently from 1.169 to 1.073 while that of profile II-II changes from 1.135 to 1.121 (in Figure 12). Although most of the safety factor is larger than 1.1, the safety redundancy is little. Compared with Table 3, the Fs in the condition of the reservoir water level rapid drawdown are a bit less than that in the fast impoundment condition. Therefore, due to the rapid descent of the reservoir water level during discharge, the glaciofluvial deposit slope is more prone to slide.

Discussions
In the analysis of overall failure mode, the slide surface of the glaciofluvial deposit slope is fixed in numerical calculation. The safety factor decreases in both fast impoundment and rapid drawdown conditions. In the rapid drawdown condition, all the factors affecting the slope stability are unfavorable. However, the seepage force and interspace pore pressure are good for slope stability in the fast impoundment conditions. Decrease of the safety factor is mainly ascribed to the degradation of the strength parameters of the glaciofluvial deposits, which dominates over other factors. This may be the reason for the soil slope often failing as the reservoir impounds.
As for the local failure mode, the slide surface is dynamic. In the fast impoundment stage, the safety factor of profile I-I decreases while that of profile II-II increases. In order to 0 1.08  9 Geofluids illustrate the difference of this phenomenon, the dynamic slide surfaces are illustrated in Figure 13. In profile I-I, the reservoir water submerges a large amount of the glaciofluvial deposit slope, and the slide surface is under the reservoir level. As for profile II-II, a major part of the slide surface is above the water. Hence, the difference changing law in safety factors is mainly caused by the shape and elevation levels of the glaciofluvial deposit.
In both rapid drawdown and fast impoundment conditions, the smallest safety factor is 1.073, which is slightly larger than 1.0. Besides, from the perspective of risk-based design, the upper and lower bounds of the stability interval are the overall and local safety factors, respectively. In the fast impoundment stage, though the local safety factor increases, the failure risk still does not lessen due to the decreasing overall safety factor. Hence, monitoring equipment and more detail in situ tests should be implemented before and while the reservoir is in service.

Conclusions
This manuscript proposed a FEM-LEM-combined method to compute slope stability under reservoir water level fluctuations. The glaciofluvial deposit slope on the left bank of the  10 Geofluids Gushui Reservoir is employed to analyze the variation of the safety factor.
The results indicate that the overall failure safety factor falls in both the fast impoundment and rapid drawdown conditions. The local increasing safety factor of profile I-I and the decreasing one of profile II-II imply that the shape and elevation levels of the glaciofluvial deposit affect the stability of the glaciofluvial deposit. The minimum safety factor occurs in the rapid drawdown condition. When the reservoir water level fluctuates in the altitude range from 2230 m to 2265 m, the most hazardous slide surface has a safety factor of 1.073, which means a small redundancy. The suggestion that monitoring equipment and more detailed in situ tests should be implemented before the projects are in service is given in the end.
The proposed FEM-LEM-combined method can reflect the impact mechanism of the water level on slope stability considering the interspace pore pressure, seepage force, and change of strength parameters. This article proposed a simple tool for the slope stability analysis under reservoir water level fluctuation and be helpful for hydraulic engineering projects in Southwest China.

Data Availability
Data is available on request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.