Topographic Effects on Three-Dimensional Slope Stability for Fluctuating Water Conditions Using Numerical Analysis

With recent advances in calculation methods, the external factors that affect slope stability, such as water content fluctuations and self-configuration, can be more easily assessed. In this study, a three-dimensional finite element strength reduction method was used to analyze the stability of three-dimensional slopes under fluctuating water conditions. Based on soil parameter variations in engineering practice, the calculation models were established using heterogeneous layers, including a cover layer with inferior properties. An analysis of seepage, deformation and slope stability was carried out with 27 different models, including three different slope gradients and nine different corner angles under five different hydraulic conditions. The failure mechanism has been shown to be closely related to the change in matric suction of unsaturated soils and the geometric slope configuration. Finally, the effect of geometry (surface shape, turning corner and slope gradient) and water (fluctuations) on slope stability are discussed in detail. Emphasis is given to comparing safety factors obtained considering or ignoring matric suction.


Introduction
Frequent landslides have long motivated research into slope stability under various conditions [1][2][3][4][5][6][7][8][9]. Continuous research and field investigations have shown that geometric configuration is an important factor in the stability of both natural and constructed slopes [10][11][12][13][14][15]. However, the topographic effect on slope stability under fluctuating water conditions has rarely been investigated, even though slope deformation and instability occur frequently in reservoir areas [16][17][18][19][20][21][22]. As shown in Figure 1, with the Dagangshan Hydropower Station Reservoir area as an example, both the strike direction of the river and the slope angle around it change in a complex manner, often resulting in slope failures of different degrees. Moreover, the hydraulic conditions of the slopes in the reservoir area are very complex, as shown in Figure 2. On one hand, the poor properties of the slope mass cover layer reduce the overall slope stability, even for slope masses with reduced strength and complex seepage environments due to water level fluctuations [23]. On the other hand, most natural and constructed slopes exhibit a complex geometric configuration and three-dimensional (3D) state, which greatly increases the slope stability variability and calculation complexity. Coupled with the regular and sudden water level fluctuations, the topographic effect on slope stability considering saturated and unsaturated conditions urgently require research.  Slope stability analysis has been a research focus area for decades, using either the soil's slide/column method based on the limit equilibrium concept, or the shear strength reduction method based on the finite element method [6,9]. The feasibility of 3D slope stability analysis has been proved to be a commonly used method today [24][25][26][27][28]. Moreover, with the continuous deepening of research, more powerful numerical techniques, like stronger technology of information gathering and processing, more realistic technology of numerical calculations, and smarter technology of monitoring and early warning, are required. The advanced numerical techniques have been proved to be a powerful tool to study both the failure stage and run-out of landslides [29][30][31][32][33]. Based on this research backdrop, considering the complex configuration of natural slopes, several studies were conducted to investigate how topographic factors (e.g., turning corner, slope gradient, turning arc, and convex-and concave-shaped surface geometry) affect slope stability under different boundary conditions [11][12][13][14]34,35]. From the researches above, one of the most important factors identified in slope stability is the constraining effect of the turning corner compared to that of general flat slopes, which will also be discussed in this study.
Similarly, many scholars have analyzed the stability of slopes subjected to various changes in the external water environment, water fluctuations or rainfall [16,17,36,37]. One key issue was the effect of water on slope stability, which was early quantified by Fredlund and Morgernstern [38] by considering the slope under different degrees of submergence and drawdown. From then on, slope stability with water effects gradually became a research hotspot. Former researches on saturated and unsaturated soils have provided better calculating methods for closer analysis of actual slope stability in this study using the 3D finite element method (FEM) with strength reduction technique. Lane and Griffiths [7] noticed that traditional analysis methods cannot cover the full range of possible critical conditions or adequately represent inhomogeneous slopes and complex loading arrangements. Thus, they proposed a chart-based operating safety approach by utilizing the finite element method to identify potentially critical conditions from rapid drawdown and partial submergence. Wheeler et al. [39] considered the significant influence of the degree of saturation on suction and the stress-strain behavior of unsaturated soil, and proposed a new elasto-plastic framework for unsaturated soils, involving coupling of hydraulic hysteresis and mechanical behavior. Griffith and Lu [40] presented results of unsaturated slope stability analyses using elasto-plastic finite elements in conjunction with a novel analytical formulation for the suction stress above the water table. Cai and Ugai [41] used the finite element analysis of transient water flow through unsaturated-saturated soils to investigate effects of hydraulic characteristics, the initial relative degree of saturation and the method with shear strength reduction technique to evaluate the stability of slopes under rainfall. Huang and Jia [42] discussed the utilization of the finite element method with the shear strength reduction technique approach in the stability analysis of soil slopes subjected to transient saturated/unsaturated seepage and presented the comparison of safety factors obtained from limited equilibrium method (LEM) and strength reduction FEM. It has been proved that the finite element method with the shear strength reduction technique was an effective method for assessing the safety factor of soil slopes under transient unsaturated seepage.
On the basis of previous research, the slope stability considering both saturated and unsaturated condition have been deeply carried out. For instance, as one of the greatest hydraulic engineering projects, the Three Gorges reservoir in China has drawn much attention worldwide. The geological structure and evolution of landslides subjected to water fluctuations in the reservoir area have been fully discussed since the reservoir's impounding, using different types of methods, including limit equilibrium and finite element approaches [18][19][20][21]. Furthermore, Gao et al. [17] adopted a 3D rotational failure mechanism to investigate the influence of water drawdown on slope stability, and gave an example to demonstrate the difference in the safety factors obtained from 2D and 3D analyses. The results of this example indicated that the 3D effect can be neglected only when the slope has a large width (width to height ratio B/H ≥ 10.0).
These studies above showed that water fluctuations were an important factor in slope stability; in addition, the use of 3D analysis using the finite element method with strength reduction technique for different hydraulic conditions was validated. Meanwhile, field investigation and numerical simulation have become the most advantageous research methods. However, in engineering practice, the environment around the river or reservoir is complex. This requires additional attention, since slope deformation and instability in reservoir areas with varying geometric configurations and water levels occur frequently. In this study, the seepage, deformation and safety of slopes were analyzed using the finite element strength reduction method, considering both saturated and unsaturated slope mass conditions. Because soil parameters vary in engineering practice, the calculation models were set to be heterogeneous, including a cover layer with inferior properties. The topographic effects (surface shape, turning corner and slope gradient) on slope stability are then discussed in detail.

Fully Coupled Flow-Deformation Analysis
In previous studies, an effective stress-strain relationship was used to describe the mechanical behavior of saturated soil when the soil particles and water phase are assumed to be incompressible. However, soil behavior is highly nonlinear and irreversible, especially considering both saturated and unsaturated conditions with water level fluctuations [7,39,40,43]. In this study, a fully coupled flow-deformation analysis was conducted to analyze the simultaneous development of deformation and pore pressure in saturated and partially saturated soils as a result of time-dependent changes in the hydraulic boundary conditions. The fully coupled flow-deformation analysis considered the unsaturated soil behavior and the matric suction in the unsaturated zone above the phreatic level. The Soil Water Characteristic Curve (SWCC) proposed by Van Genuchten [44] was applied in this study to describe the hydraulic parameters of the groundwater flow in unsaturated zones, reflecting the soil's capacity to keep water at different stresses. The Van Genuchten function is a three-parameter equation that relates the saturation to the pressure head : where is the suction pore stress; is the unit weight of the pore fluid; is the residual saturation, which describes the part of the fluid that remains in the pores, even at high suction heads; is the saturation in a saturated situation, which defaults to 1.0; is a fitting parameter related to the soil's air entry value, which must be measured for a specific material; is a fitting parameter which is a function of the rate of water extraction from the soil once the air entry value has been exceeded; is a fitting parameter that is used in the general Van Genuchten equation.
The deformation model applied in this study was formulated within the framework of continuum mechanics. The continuum description was discretized according to the finite element method. Each element consists of a number of nodes. Each node has a number of degrees of freedom that correspond to discrete values of the unknowns in the boundary value problem to be solved. In the present case of deformation theory, the degrees of freedom correspond to the displacement components. Subsequently, the slope stability was investigated using finite element analysis with the shear strength reduction technique, with the inputs of pore-water pressure and degree of saturation from the seepage analysis.

Strength Reduction Method
In previous studies, the limited equilibrium method (LEM) is the most extensive means for evaluating slope stability [45][46][47][48][49]. However, its difficult application to the finite element method, together with the failure mechanism assumptions, limited its development. This type of difficulty can be overcome by introducing a FEM with the shear strength reduction technique [7,[50][51][52][53][54][55][56]. The strength reduction method (SRM) was first introduced in slope stability analysis by Zienkiewicz et al. [57], but was not widely applied due to the lack of finite element analysis procedures and strength criteria selection in the next few years [32,33,[50][51][52][53][54][55][56]. With the development of computer technology and continued scholarly research, this SRM has gradually been recognized in stability analysis. Furthermore, SRM has taken stability analysis to a new level.
The shear strength reduction factor is defined as the ratio of the maximum shear strength of the slope mass to the actual shear stress generated by the external load while the external load remains unchanged. When the slope reaches the failure state, the displacement on the sliding surface will change abruptly, resulting in large and unlimited plastic flow and non-convergence in the finite element calculation. Therefore, the non-convergence of force and displacement is used as the failure criterion for the slope. The sliding failure surface (the area with a sudden change in plastic strain and displacement) can be obtained automatically, based on the elasto-plastic calculation results. Simultaneously, the safety factor (FS) can be determined. The process of solving for the safety factor can be simplified to the following formulas. The total multiplier ΣMsf is used to define the value of the soil strength parameters at a given stage in the analysis: where the strength parameters with the subscript 'input' refer to the properties entered in the material sets, and parameters with the subscript 'reduced' refer to the reduced values used in the analysis. ΣMsf is set to 1.0 at the start of a calculation to set all material strengths to their input values. When slope failure occurs, the safety factor (FS) is given by:

Validation
Slope stability analysis has received scholarly attention for a long time. A typical slope model was used to verify the correctness of the calculations and settings. It was first proposed by Zhang [58] using an extended Spencer's method, and has been widely used by various investigators, as part of the validation of their particular 3D slope stability methods [11,13,27,[59][60][61][62][63]. In this study, a similar 3D model was established; the results are shown in Table 1. As shown, the result with the coarse mesh distribution model fits well with previous results; thus, the coarse form is used in the following analysis.

Numerical Model and Scheme
With the rapid development of computer hardware and numerical software, 3D slope analysis is gradually becoming mainstream for seepage, deformation and safety analysis [27,28]. In this paper, Plaxis 3D was used for numerical simulation. Furthermore, a fully coupled flow-deformation analysis was used to analyze the simultaneous development of deformation and pore pressure in saturated and partially saturated soil as a result of water fluctuations. The strength reduction method was used to determine global safety factors after a certain period of flow and deformation analysis. There were 27 different calculation conditions in this study, and the results are discussed in a later section.

Geometric Entities and Material Properties
First, a general 3D slope model was constructed by extending the 2D model longitudinally. The geometric and soil information is given in Figure 3. The height of each slope was set to 30 m, and the bank slope subjected to water level fluctuations was divided into two identical 50-m long slopes by the axis of symmetry at the turning point. Considering the heterogeneity of slope soil and the inferior cover layer properties, the slope and foundation consist of two parts with poor and good soil conditions, respectively. The model has four different types of soil mass: the relatively weak cover layer, the general slope mass, the relatively weak foundation and the general foundation. The cover layer is 5 m thick. The general soil mass outside the turning corner range was set to at least 85 m in thickness to avoid the effect of water level fluctuations on the seepage field. There is no flow in the furthest soil mass by default. The foundation is a total of 30 m high, with a 10-m high relatively weak foundation. Three different slope gradients (2:1, 1:1 and 1:2 from steep to gentle, respectively) were investigated to determine the effect of slope gradient on stability. Meanwhile, nine different corner angles (90°, 120°, 135° and 150° with convex shapes, and 210°, 225°, 240° and 270° with concave shapes, along with 180° for comparison) were investigated with the three slope gradients above. The 27 calculation conditions in this study are shown in Figure 4. Each model was established with the same structural layering as the 2D extending model (180°). Table 2 shows the mass properties of each layer. A simplified selection of the most common soil types (coarse, medium, medium-fine, fine and very fine nonorganic material and organic material) based on the Hypres topsoil classification series [64] was adopted for subsequent use in the Van Genuchten model.

Mesh Generation and Boundary Conditions
After fully defining the model geometry, it was divided into finite elements for finite element calculations. However, previous studies have proven that the results from the finite element method with the strength reduction technique are sensitive to the design of the mesh [65]. Referring to the verification model of Zhang 1988 [58], the coarse element distribution results in this study (shown in Table 1) fit better with those of Zhang. Therefore, the overall model geometry was set to a coarse element distribution in this study. Furthermore, since the model geometry, including the edges, corners and structural objects, was complicated, the relatively weak layer required a finer finite element mesh. To balance the accurate numerical results of a sufficiently fine mesh against the excessive calculation times of very fine meshes, local refinement was used during mesh generation. A local coarseness factor was used to indicate the relative element size; the coarseness factor value set to 1.0 by default for most geometry entities. The cover layer was refined four times and the relatively weak foundation was refined 2 times. An example of the slope element distribution for a turning corner of 90° and a slope gradient of 1:2 is shown in Figure 5. The element distribution on the turning surface was much finer than that in distant areas.
In addition, the boundary conditions assumed in the 3D numerical model were important; many scholars have investigated the boundary effects with different methods [27,28,66]. For a 3D slope with any complex shape, as shown in Figure 5, boundary surfaces generally fall into five categories (bottom surface, front surface, turning surface, back surface and side surface) [13]. In this study, all movements were restrained (in the normal, dip and strike directions) in the bottom plane of the model (i.e., ux, y, z = 0, where u is the displacement). The front surface and back surface were restrained in the normal direction, but free to move in the dip and strike directions. However, the boundary conditions for the side surfaces were flexible (generally classified into three types) as discussed in [11][12][13]27]. In this study, the side surface condition was set comparing the results and conclusions in previous studies. Combining the practices in the studied models, to reduce calculation time, the side surfaces were restrained only in the normal direction for the convex-shaped slopes (corner angles of 90°, 120°, 135° and 150°) to determine the actual slope failure mechanism. The side surfaces were completely restrained (in three directions) for concave-shaped slopes (with corner angles of 210°, 225°, 240° and 270°) to investigate the constraining effects of the turning corner and avoid a landslide near the side surfaces.

Hydraulic Conditions and Calculation Phases
After creating the model entity, a single borehole was established immediately to create a horizontal water surface (5 m as a regular condition) that extended to the model boundaries. The global water level then was specified as the established borehole water level to generate a simple hydrostatic pore pressure distribution (phreatic calculation type) for the full geometry. This global water level was used throughout the calculation process for the specific model. For the fluctuating water conditions, defining hydraulic conditions as groundwater flow boundary conditions was required. The specific hydraulic conditions have priority over global model conditions-if the global water level was used to define the hydraulic boundary conditions (groundwater head), the global water level would be replaced by the results of the groundwater flow calculation. In this study, the surface flow boundary conditions acting on the slope model's turning surface were selected as the hydraulic conditions throughout the calculation process to meet the fluctuating water level conditions. The specific surface flow in each calculation phase is listed by calculation phase in a later section.
To build the calculation phases, a contrasting initial phase was established under the condition of regular water level without fluctuation before comparing the topographic effects on slope stability. A steady-state pore pressure with a 5-m water head was generated as the input for the initial phase, in which the calculation type was set to 'Gravity loading'. The initial phase calculation determined the initial conditions of the pore water pressure and stress field, providing the basis for subsequent calculations. For the next condition (water rising), the time-dependent hydraulic boundary condition that water levels rise rapidly from 5 to 25 over 5 days was input before the transient groundwater flow analysis. The fully coupled flow-deformation analysis was selected for this phase to analyze the simultaneous development of deformation and pore water pressure in both saturated and unsaturated soils. Since the deformation due to gravity loading is physically meaningless, the 'Reset displacements to zero' option was chosen after 'Gravity loading' to remove these displacements. After the rapid water level rise, a five-day period of high water level was established to investigate the effect of pore water pressure due to infiltration. Similarly, the fully coupled flow-deformation analysis was selected for this and all later phases. The sharp drawdown of the water level was set to the same speed as that of the water rising (five days), after which a five-day period of low water level was established. The effects of the high water level infiltration and the low water level drainage are discussed later.

Results
Engineering practices have shown that geometric configurations have important effects on slope stability. To assess the topographic effect of the complicated terrain of different natural slopes on their stability under fluctuating water conditions, many numerical calculations were carried out. The three-dimensional geometric factors were simplified to different corner angles and slope gradients to control variables and build numerical calculation models. A total of 27 different conditions were analyzed using nine different corner angles (90°, 120°, 135° and 150° with convex shapes, 210°, 225°, 240° and 270° with concave shapes, along with 180° for comparison) and three different slope gradients (2:1, 1:1 and 1:2 for a very steep slope, a steep slope and a gentle slope, respectively). The above slope entities with different geometric conditions were simultaneously analyzed under different hydraulic conditions; the regular condition, water rising condition, water drawdown condition, and high and low water level conditions were evaluated.

Influence of Turning Corner
The turning corner of a natural 3D slope is one of the most representative geometric factors. Usually, convex-shaped and concave-shaped slopes will have different stability patterns, since they have different constraint effects.
For convex-shaped slopes, Figure 6 shows the safety factors of slopes with corner angles were very near to or smaller than those of the 2D extended model (with a turning corner of 180°). The safety factor in the initial phase decreased as the slope turning corner value decreased from 150° to 90°, which indicates that a prominently exposed soil mass exhibits greater stability risks. Figure 7 shows the potential failure mechanism of each slope under regular hydraulic conditions with steady seepage. As shown, a diminishing turning corner gradually exposes the originally flat slope, which increases the slope instability. The smaller the turning corner, the thinner the exposed slope mass is. The convex-shaped slopes with corner angles of 150° and 135° are sufficiently large that they showed similar security as the 2D extending model (180°). The near-flat corner angles increase the volume of the potential landslide mass, which requires a greater external disturbance. Furthermore, a similar action of favorable boundary constraints due to the symmetric slope mass distribution may result in a slight increase in slope stability as the turning corner decreases from 180° to 150°. This confirms the effect of the boundary constraints from Nian et al. [11] and Zhang et al. [13]. Unfortunately, as the turning corner continues to decrease, the slope stability declines and drops sharply at 90°, as shown in Figure 6, resulting in the instability of thinner and smaller landslide masses in Figure 7. For convex-shaped slopes, the slightly favorable boundary constraints due to a convex symmetric slope mass distribution were not sufficient to offset the negative effects of slope exposure.  As is commonly known, a concave slope can provide beneficial constraints, with each part divided by the axis of symmetry. This contributes to the higher stability of slopes with corner angles. To prevent the formation of a landslide body from a thinner slope mass on the sides and determine the degree of the constraint effect, in this study, the deformation boundary conditions on both side surfaces were set to completely restrained in all three (normal, dip, and strike) directions. For concave-shaped slopes, as shown in Figure 6, regardless of the degree of increase of the slope security, whenever there was a turning corner, the stability increased compared to the 2D extending model (180°). This indicated significantly favorable boundary constraints due to the concave symmetric slope mass distribution. Figure 8 shows the potential failure mechanism of each slope under the regular hydraulic condition with steady seepage. As shown in Figure 6, similar to the convex-shaped slopes, concave slopes with corner angles of 150° and 135° had the highest stability, in stark contrast to other slopes with corner angles of 90° and 120°. Comparing the potential failure mechanism of the two slope types with different stabilities, it was found that the failure mass was divided into two parts by the axis of symmetry. This reduced the beneficial constraint effect due to the symmetrical distribution and further explained why the stability did not increase, but decreased, as the turning corner continued to decrease further.

Influence of the Slope Gradient
In this study, three slope gradients (2:1, 1:1, and 1:2 from steep to gentle, respectively) were calculated and analyzed successfully considering the unsaturated zone and conditions above the phreatic line. As shown in Figure 9, one routine conclusion was confirmed that the steeper the slope was, the lower the level of safety was. This result verifies the disturbance sensitivity of steep slope soil. Additionally, the safety factor variations of all three slope gradients in each water fluctuation time step are shown in Figure 9. Slopes with gentler gradients always had the highest stability at any point in the water fluctuations. Furthermore, the slope gradient did not have a significant effect on the safety factor variation pattern due to water fluctuations. The most dangerous part of the potential failure mass is highlighted in Figure 10. The two areas that are most likely to slide in normal water level conditions were the most exposed mass on the axis of symmetry and the mass on the slope foot. Comparing the potential failure mechanism of the three slopes with 90° corner angles, the slip surface became deeper and wider as the slope gradient increased. This also confirms that stability drops sharply with increasing gradient. Figure 10. The difference in failure mechanism between slopes with three gradients: turning corner equal to 90°.

Effect of Water Level Fluctuations
In engineering practice, water level fluctuation is a common phenomenon acting on the bank slope of the reservoir area. Regular water storage and drainage in reservoirs may cause slope instability, such as the landslides that occurred at the Three Gorges Reservoir area in China [16][17][18][19][20][21][22][23]. It is then conceivable that rapid water level fluctuations may lead to higher risks to slope deformation and stability. In this study, rapid rise and sharp drawdown water level fluctuation conditions were considered, with the continuous infiltration at the high water level and the successive slope drainage at the low water level. The slope mass safety factor was obtained for calculated time steps; the results are shown in Figure 11. The pore water environment for different time steps is shown in Figure 12. In numerical calculation time, the zeroth day was set to calculate the seepage and stability of the bank slope for the normal water level condition considering steady-state groundwater flow. The slope calculation models remained stable initially. The rapid water level rise of 20 m was then set evenly from the first day to the fifth day under the transient (time-dependent) groundwater flow. Immediately after the rapid rise, a steady-state groundwater flow condition was assumed for the next five days (sixth to tenth days) considering continuous infiltration at a high water level. As shown in Figure 11, the stability increased as the water level rose, but decreased with continuous infiltration. The former result illustrates the beneficial effect of the water pressure acting as a back pressure. The latter result indicates the adverse effects of water infiltration caused by the increase in pore pressure. After a five-day high water level, from the eleventh day on, the water level dropped sharply for five days until the twentieth day, successively followed with a five-day drainage at the low water level. The slope stability decreased as the water level decreased, resulting in an even lower stability than the initial condition, due to the adverse effects of pore pressure. Thus, the slight increase in stability with the drainage of water over the next five days can be easily understood.  Figure 13 shows this study's model calculation nodes. The model surfaces are represented by many nodes, and the stress points can be seen in the model body. The displacement data for 17 monitoring nodes were obtained. Their positions in the model are highlighted in Figure 13, representing the condition of the turning surface. Divided into two parts by the axis of symmetry, 5 columns of monitoring nodes on the turning surface were symmetrically distributed. Four columns of the nodes (nodes 1-12) were positioned on the outmost and the middle of the two surface parts, at heights of 0, 10 and 20 m, respectively. The remaining columns of nodes (node [13][14][15][16][17] were distributed on the model's symmetry axis, at heights of 0, 10, 15, 20 and 25 m, respectively. The obtained node displacements at each phase are shown in Figures 14 and 15. The displacements of surface nodes are represented by different colors and icons for different heights (i.e., 0 m in white with ○, 10 m in pink with □, 15 m in green with ×, 20 m in blue with △, and 25 m in yellow with ◇), and different line types for different locations from the axis of symmetry. In general, the increasing trend of displacements with water level fluctuations was similar for all monitoring nodes. The general rule was that the node displacements increased rapidly during both the water level rise and fall processes, but increased very slowly as the water level remained unchanged.  As shown in Figure 14, the deformations of two nodes with a symmetric distribution were basically the same (i.e., deformations between nodes 1 and 2, nodes 3 and 4, nodes 5 and 6, nodes 7 and 8, nodes 9 and 10, and nodes 11 and 12). The deformations of nodes located different distances from the axis of symmetry and with different heights showed different changing characteristics for different turning corner values. One interesting phenomenon was the difference between the increasing trends of deformations for different turning corner values. Generally, when the water level was rising, the deformations increased slowly with a gradually increasing speed. However, the deformation increased rapidly with a gradually reducing speed for slopes with corner angles, while it increased rapidly with a gradually rising speed in 2D-extending slopes. Moreover, the location of the most dangerous part of the slope (i.e., the part with the largest deformation) was very different for slopes with different turning corner values. For a general shape slope (a 2D extending slope with a turning corner of 180°), the maximum deformation zone located lower, with a height of 10-15 m, while it located higher, at 15-20 m, on slopes with corner angles. This also confirms the slope failure mechanism shown in Figures 7 and 8. The reason was probably that the failure of general slopes without corner angles often requires a large sliding mass, resulting in higher stability. The slipping surfaces of these slopes are usually located deeper and wider, with the maximum deformation zone located at a lower elevation. As shown in Figure 15, the increase in deformation for slopes with different gradients shows a similar pattern. However, the very steep slope of 2:1 with a turning corner of 90° failed when the water level dropped below 17 m during the 13th to 14th day of the water level fluctuations, manifesting as a dramatic increase in displacement.

Engineering Significance
For convex-shaped slopes, a turning corner usually has negative effects on stability, especially when the corner is relatively small and the thinner exposed slope mass greatly reduced the slope stability. Only when the turning corners were sufficiently large (150° to 135°) that they showed similar security as the 2D extending model (180°). In engineering practice, the slopes with small turning corners may be more unstable. In the reservoir area, under the influence of common water level fluctuations, those slopes are prone to deformation and even failure. Thus, the very small turning corner should be avoided (e.g., <120°). When it is necessary to learn the impact of other turning corners in different valves on stability, engineers can determine the range of influence based on the interpolation method with Figure 6.
For concave-shaped slopes, whatever the valve of the turning corner was, it always had a beneficial constraint effect on slope stability contributing to the symmetrical distribution. In engineering practice, the concave-shaped slopes may be more stable than other cases. However, the extremely small concave corner should also be avoided since the failure mass may be divided into two parts by the plane of symmetry and the beneficial constraint effect is likely to be reduced.
One routine conclusion was confirmed in this study that the steeper the slope is, the lower the safety factor, which means, the slopes with gentler gradients always had the highest stability. In engineering practice, the high and steep slopes usually show bad stability itself and prone to slope failure because of their sensitive responses to external interference. Although the stability of the gentler slope is always better, the engineer still needs to maintain the balance between the safety and economy of the project to determine the location and gradient of the slope. If necessary, measures such as slope protection and excavation can be taken.

Discussions
As is commonly known, there are few high and steep soil slopes in engineering practice. The poor physical and mechanical properties of soil mass and sensitive responses to the external interference, together with the slope configuration contribute to the failure of them. For stability analysis of high and steep soil slopes, when soil mass below the phreatic level was considered to be fully saturated, and suction effects in the soil above the phreatic level were ignored, the slope may calculated to have failed. Thus, in numerical calculation, a more realistic value of the safety factor will be obtained for stability analysis when suction is considered, where the suction acted as a positive value (tension) in the pore water pressure. This value is generally higher than a conventional safety factor ignoring suction, which also explains the successful calculation of high slope models with very steep gradients in this study. In calculation, the slope mass is divided into saturated and unsaturated zones by the phreatic line. When considering suction, the soil saturation depended on the soil-water retention curve, as defined in the corresponding material data sets. A saturation level was defined in the calculation by setting the water conditions to partially saturated; the resulting pore stress (suction) is calculated accordingly.
In this study, a comparison of the two cases (considering or ignoring suction) was conducted under the regular steady-state groundwater flow condition. To improve the calculation efficiency, 2D-extending models in three gradients were applied in this discussion. The effect of suction on stability can be simplified to the increase in effective cohesion [67]. As shown in Table 3, the condition of ignoring suction was carried out, and the assumed compensation effect of suction on effective cohesion was determined by increasing the safety factor of the calculation model to an equal level. It was confirmed that the suction in the partially saturated zone did have beneficial effects on slope stability. For a very steep slope (2:1), the FS was 1.121 considering suction. Unfortunately, when ignoring suction, the slope failed. In other words, this high and steep slope may not exist if calculated without considering suction. For a gentle slope, the effect of suction was approximately equal to a 10 kPa increment in effective cohesion, and 15 kPa for a steep slope, 25 kPa for a very steep slope. With an increasing slope gradient, the beneficial effects of suction increased, while the slope stability decreased.
Taking suction into consideration, the topographic effects on slope stability were studied under transient unsaturated seepage. Parametric analysis was conducted with the results of combinations of multiple configurations and varied hydraulic conditions. With generalizing the results of this study, some constructive suggestions can be directly adapted and used by other engineers. However, more research is still required for an accurate assessment of suction's role in a partially saturated mass and the effects of more complex geometries on slope stability, which should be addressed in future studies.

Conclusions
In this study, a three-dimensional finite element strength reduction method was used to analyze the seepage, deformation and stability of slopes with different configurations. Combined with different soil parameters obtained from field investigation, the calculation models were established to be heterogeneous, with an inferior property cover layer. Considering both saturated and unsaturated zones of the slope mass, five hydraulic conditions varied with water level fluctuations (regular steady-state groundwater flow condition, water level rising and dropping conditions, high water level infiltration and low water level drainage) were conducted. Finally, the topographic effects (surface shape, turning corner and slope gradient) on slope stability were discussed from models with three different slope gradients (2:1, 1:1 and 1:2 from steep to gentle, respectively) and nine different corner angles (90°, 120°, 135° and 150° with convex shapes, and 210°, 225°, 240° and 270° with concave shapes, together with 180° for comparison). Based on the calculation results, the following conclusions are drawn.
(1) For convex-shaped slopes, the safety factor in the initial phase decreased as the slope turning corner decreased from 150° to 90°, indicating greater stability risks for the prominently exposed soil mass. The slope stability increases slightly from 180° to 150°, due to the favorable boundary constraints of the symmetric slope mass distribution, but these were not enough to offset the negative effects of the slope exposure.
(2) For concave-shaped slopes, whenever there was a turning corner, the stability significantly increased compared to the 2D extending model (180°). This indicated significantly favorable boundary constraints due to the concave symmetric slope mass distribution. However, the beneficial constraint effect of the symmetrical distribution was reduced for small corner angles (90-120°) because the failure mass was divided into two parts by the axis of symmetry.
(3) One routine conclusion was confirmed: the steeper the slope is, the lower the safety. Furthermore, the slope gradient did not significantly affect the safety factor variation pattern with water fluctuations; the slopes with gentler gradients always had the highest stability.
(4) The stability increased as the water level rose and decreased as the water level decreased, resulting in even lower stability than the initial condition due to the adverse effects of pore pressure. Conversely, stability decreased with continuous infiltration but increased slightly for a low water level with water drainage. The former result illustrates the beneficial effect of the water pressure acting as a back-pressure. The latter result indicates the adverse effects of water infiltration and the beneficial effects of water drainage, because of the variation in pore pressure.
(5) The node deformation increased rapidly during both water level rise and fall, but it increased very slowly as the water level remained unchanged. For a general shape slope (a 2D extending slope with a turning corner of 180°), the maximum deformation zone located lower, with a height of 10-15 m, while it located higher, at 15-20 m, for slopes with corner angles. Furthermore, some slopes failed when the water level dropped, manifesting as a dramatic displacement increase. (6) In engineering practice, concave-shaped slopes always show better stability than that of convex-shaped slopes. Convex slopes with small turning corners should be avoided in engineering since they are prone to deformations and even landslides under the influence of common water level fluctuations. Meanwhile, the engineer needs to maintain the balance between safety and economy to determine the location and gradient of the slopes in engineering projects.