Seepage, Deformation, and Stability Analysis of Sandy and Clay Slopes with Di ﬀ erent Permeability Anisotropy Characteristics A ﬀ ected by Reservoir Water Level Fluctuations

: Evaluation of slope stability under water level ﬂuctuations is an important topic in the Three Gorges Reservoir (TGR) in China. However, most of the previous studies regarded slope soil as isotropic material, or only considered the inﬂuence of anisotropy ratio ( k r = k x / k y ) but ignored the anisotropy direction ( α ). Meanwhile, the pore pressure–stress coupling was rarely considered in the previous numerical simulations. In the present study, the SIGMA / W and SLOPE / W modules in Geo-studio are utilized to carry out the numerical simulation of Caipo slope under the drawdown of the reservoir water level, and the anisotropy ratio ( k r ) as well as the anisotropy direction ( α ) of two kinds of soils (clay and sand) are included. Results show that the anisotropy ratio k r and anisotropy direction α decrease the inﬁltration capacity of the soil, which increases the inﬁltration line hysteretic elevation (ILHE) as well as maximum horizontal displacement (MHD), and reduces the minimum safety factor (MSF). The slope toe ﬁrstly fails with the drawdown of water level. The inﬂuence of reservoir water level drop on seepage, deformation, and stability of the sand slope is less than that of the clay slope. For the sandy soil slope, it is not only necessary to consider the inﬂuence of k r , but also the inﬂuence of α . For the soil slope, we can only consider α in order to simplify calculation.


Introduction
With the rapid development of China's economy, large amounts of hydropower stations have been built, which provides a great guarantee for China's clean energy production, and the Three Gorges Project is the typical one. Since the completion of the Three Gorges Reservoir (TGR) in 2003, landslide disasters in the reservoir area have occurred continuously under the cycle of reservoir water level fluctuations, as is shown in Figure 1a [1,2]. The TGR began to impound water in June 2003, and the elevation of water level was 135 m, which was maintained until October 2006, and the number of landslides was the largest in this period. During the period from October 2006 to October 2008, the October 2008, the reservoir water level fluctuated between 145 m and 165 m, and the number of landslides in this period was more concentrated. From October 2008 to December 2015, the reservoir water level stabilized between 145 and 175 m, and the number of landslides decreased significantly compared with the completion period but still happened occasionally. Therefore, it can be found that the fluctuations of reservoir water level is an important inducement of landslides in the TGR. The main reasons for slope instability caused by reservoir water level fluctuations are reflected in these three aspects: (1) The increase of reservoir water level leads to the increase of pore water pressure in slope soil, which contributes to the reduction of strength and effective stress of slope soil [3][4][5]. (2) The decrease of reservoir water level leads to the change of pore pressure in slope soil, which lags behind the reservoir water level and causes the dynamic seepage force, which drives the slope soil to move away [6][7][8]. (3) The decrease of reservoir water level also leads to the reduction of the confining pressure over the slope, which reduces the slope stability. The main consequences of slope instability caused by reservoir water level fluctuations are reflected in these two aspects: (1) Landslides cause damage to buildings in the area where the slope is located. (2) Landslides threaten the lives and properties of the residents around the disaster area. Typical landslide events caused by the change of reservoir water level such as the Bazimen landslide, which is located in Xiangxi Village, Guizhou Town, Zigui County, Hubei Province, China, and happened in June 2003 [9], as shown in Figure 1b, the Taping landslide, which is located in Quzhi Township, Wushan County, Chongqing City, China, and happened in June 2012 [10,11], as shown in Figure 1 (c), all caused great losses of lives and properties. Therefore, it is important to grasp the laws and contributing factors of slope instability caused by reservoir water level fluctuations in order to correctly understand the mechanisms of reservoir water seepage and to prevent and control landslide disasters. Bazimen landslide [10]. (b) Taping landslide [11]. (c) Taping landslide [11]. (a) Bazimen landslide [10]. (b) Taping landslide [11]. (c) Taping landslide [11].
Scholars have conducted large amounts of research studies on the mechanical and seepage responses of slope soil to the reservoir water level fluctuations; the research results mainly focused on field investigation, model tests, and numerical simulation. Field investigations are the most direct and accurate way to study the influence of reservoir water level fluctuations on landslides, but the cost and time are the largest. Wang et al. [12] carried out a geological survey and analyzed the monitoring data of the Wushaxi landslide; the deformation characteristics under the combination of reservoir water level and rainfall were obtained and a new analysis method was put forward. Javed et al. [13] carried out a drilling survey and analyzed the monitoring data to find the position of the sliding surface and studied the relationship between landslides and surrounding environmental factors. Yan et al. [14] conducted a geological survey on the Donglingxing landslide and analyzed the seepage characteristics under rainfall and reservoir water level fluctuations; the drainage measures were also studied according to the specific characteristics. The above geological investigation studied the mechanical mechanism and seepage characteristics of the landslide, but the investigating period was too long and the whole landslide could not be completely evaluated. Model tests have solved this problem. Luo et al. [15] conducted a model test of the failure process of a soil slope with different angles under the drawdown of the water level and found that the position of the slide surface became deeper with the increase of the slope angles. Hu et al. [16] carried out laboratory model tests of the deformation mechanism of landslides under water level fluctuations and concluded that the increase of pore water pressure inside the slope soil was the key factor triggering the slope failure. Yang et al. [17] studied the failure process of a silt slope under the sudden drop of reservoir water level and found that failure mode and crack propagation were closely related to the mechanical properties of unsaturated soil. However, the similarity between the model test and the real landslide is still a problem to be solved. Numerical simulations are an important means to verify the laboratory test results and predict landslide instability, which has the advantage of low cost and is less time consuming. Previous studies mostly utilized the unsaturated limit equilibrium method to calculate the slope safety factors; for example, Hu et al. [18] used the finite element method to simulate the stability of the Binjiang landslide based on the unsaturated seepage theory. Jian et al. [19] calculated the seepage and stability of the Qianjiangping landslide under a combination of rainfall and reservoir water level fluctuations. Song et al. [20] carried out a numerical analysis of the seepage characteristics and slope stability under different reservoir water drawdown rates. Huang et al. [21] analyzed and discussed the seepage and stability characteristics of a slope with different soil permeability. However, the seepage variation under reservoir water drawdown is a typical saturated-unsaturated process; the reservoir water not only has static water pressure and buoyancy effect on slope soil [22][23][24][25] but also causes matrix suction change, dynamic water pressure, and excess pore water pressure [26][27][28][29][30][31]. The altering of pore water pressure caused by the flow of fluid in unsaturated soil and external water load changes the effective stress and seepage force, which means the stress distribution inside the soil changes, which causes the deformation and meanwhile causes a change of seepage. Therefore, the analysis of slope under the reservoir water fluctuations includes three aspects, namely seepage, deformation, and stability, but previous numerical simulation rarely used a unified finite element method for seepage calculation, deformation analysis, and stability prediction.
Meanwhile, most of the previous studies regarded the slope soils as isotropy materials. According to the SEM microcosmic study of Song et al. [32], the anisotropy of the permeability coefficient is caused by the flocculation microstructure of clay and other soils; meanwhile, soil permeability coefficient anisotropy is greatly affected by dry density and freeze-thaw cycles [33]. Generally speaking, the hydraulic conductivity anisotropy ratio (k x /k y ) can reach 2~10, and can possibly reach 100 for clay soil [34]. There also exist joint cracks in rock slopes [35,36], which lead to strong anisotropy of seepage characteristics. The anisotropy coefficient not only has great influence on the transient seepage, but it also has an impact on the safety factor of the slope. According to Mahmood et al. [37], the difference of slope safety factors between considering and not considering soil anisotropy is about 40% [38]. However, most of the previous studies ignored the seepage anisotropy of slope soil, and the research results of soil anisotropy of slopes were few and incomplete. Yeh et al. [39] took into account the hydraulic conductivity anisotropy ratio of soil slope, and simulated the seepage characteristics and local safety factor, but ignored the hydraulic conductivity anisotropy direction. In fact, the horizontal permeability coefficient k x and vertical permeability coefficient k y coincide with the natural coordinate axis only in some special cases such as in layered crushed earth dams or naturally deposited layered soil. In nature, there are more cases where the anisotropy principal direction does not coincide with the coordinate axis. The actual requirements of the project cannot be accurately reflected only by considering the hydraulic conductivity anisotropy ratio.
In view of the shortcomings of previous studies, the seepage force is considered in the traditional slice method to analyze the slope stability, which is also combined with stress-strain analysis to establish a new numerical method. The relationship between suction and permeability coefficient is obtained according to the soil-water characteristic curves (SWCCs) to calculate the seepage field, and the calculation equation is discretized in time and space. The unified finite element method can quantitatively reflect the interaction between stress and seepage in soil and the change of pore water pressure. The Capo slope in the TGR is taken as a case study in this paper [40]. The numerical analysis is carried out under the condition of reservoir water drawdown, and two kinds of soil (clay and sand) and the hydraulic conductivity anisotropy ratio and direction are included. Although the influences of vegetation on slope stability have been declared in previous studies [41][42][43], we have neglected this for ease of reading and understanding. The research results provide some reference for the understanding of the seepage anisotropy law and prevention of landslides.

Methods and Theory
The SIGMA/W and SLOPE/W modules in Geo-studio were utilized to calculate the pore pressure-stress coupling and slope stability. The SIGMA/W module creates three equations for each node, two of which are the equilibrium equations (displacement) and the third one is the continuous flow equation (pore water pressure). Then the variations of displacement and pore water pressure can be given by solving three equations at the same time. The slope safety factors can be obtained by SLOPE/W based on the pore water pressure calculated by SIGMA/W.

Theory of Unsaturated Seepage
The seepage control equation is derived from the saturated and unsaturated Darcy law [44], which can be expressed as: where x and y are the coordinates in the direction of x and y, respectively; k x is the hydraulic conductivity in the x-direction; k y is the hydraulic conductivity in the y-direction; H is the total head; Q is the applied boundary flux; t is the time; m w is the slope of the storage curve; and γ w is the unit weight of water.
Applying the Galerkin method of weighed residual to the governing differential equation, the finite element for two-dimensional seepage equation can be derived as: where [B] is the gradient matrix; [C] is the element hydraulic conductivity matrix; [H] is the vector of nodal heads; <N> is the vector of interpolating function; q is the unit flux across the edge of an element; τ is the thickness of an element; λ is the storage term for a transient seepage equal to m w γ w ; A is a designation for summation over the area of an element; and L is a designation for summation over the edge of an element. In an abbreviated form, the finite element seepage equation can be expressed as: where [K] is the element characteristic matrix; [M] is the element mass matrix; and {Q} is the element applied flux vector.

Theory of Pore Pressure-Stress Coupling
For a two-dimensional space, the incremental stress-strain relationship can be written as: where ε is the normal strain; γ is the engineering shear strain; σ is the normal stress; τ is the shear stress; u a is the pore-air pressure; u w is the pore-water pressure; [D] is the drained constitutive matrix; and {m H } can be expressed as where H is the unsaturated soil modulus for soil structure with respect to matrix suction (u a − u w ).

Safety Factor for Unsaturated Soil
SLOPE/W adopts the Morgenstern-Price method based on limit equilibrium theory to calculate the safety factor; the modified method strictly satisfies the force balance and torque balance, and the calculation accuracy is high. The expressions are as follows: where c i is the cohesive strength for every soil slice; i is the soil slice number; W i is the weight of every soil slice; P i is the water pressure; β i is the angle of the bottom of the soil slice; b i is the length of every soil slice; ϕ i is the frictional strength for every soil slice; r i is the radius of the sliding arc; and F s is the safety factor.

Numerical Model and Boundary Conditions
The Caipo slope in the TGR was taken as a case study [40], which is located in the Yanzi Town, Hefeng County. The distance of the slope from the dam was about 1.32 km, the shape of which was similar to that of the rectangle. The slope angle was 15-40 • . The elevation of the front edge was 160 m, and the elevation of the trailing edge was 220 m. The thickness of the sliding mass was 2-4 m, the total area was about 1.2 × 10 4 m 2 , and the total volume was about 3.6 × 10 4 m 2 . In order to reduce the influence of boundary conditions on the calculation results, the ranges of the slope foot and top were extended. The whole model was divided into 1084 nodes and 1001 elements, which is shown in Figure 2.
For the purpose of investigating the seepage and deformation characteristics at different positions, three sections were set to reflect the effect of hydraulic conductivity anisotropy, whose positions were x = 100 m (top of the slope), x = 200 m (middle of the slope), and x = 300 m (toe of the slope). The monitoring sections ran through the landslide. The initial conditions were the steady seepage field calculated by the fixed water level boundaries of AB 9 m and DF 175 m. The boundary conditions are as follows: AB is the fixed water level boundary of 190 m, DEF is the reservoir water level drawdown boundary of 175 m to 145 m, and AF and CD are the impervious boundaries.

Unsaturated Soil Properties
The soil-water characteristic curves (SWCCs) adopted the Fredlund and Xing model, which can be written as [45] where kw is the calculated conductivity for a specified water content or negative pore-water pressure; ks is the measured saturated conductivity; Θs is the volumetric water content; e is the natural number 2.71828; y is a dummy variable of integration representing the logarithm of negative pore-water pressure; i is the interval between the range of j to N; j is the least negative pore-water pressure to be described by the final function; N is the maximum negative pore-water pressure to be described by the final function; Ψ is the suction corresponding to the jth interval; Θ' is the first derivative of the equation. Θ can be described as where a is the air-entry value of the soil; n is a parameter that controls the slope at the inflection point in the volumetric water content function; and m is a parameter that is related to the residual water content. C(Ψ) is a correcting function defined as where Cr is a constant related to the matric suction corresponding to the residual water content.

Unsaturated Soil Properties
The soil-water characteristic curves (SWCCs) adopted the Fredlund and Xing model, which can be written as [45] where k w is the calculated conductivity for a specified water content or negative pore-water pressure; k s is the measured saturated conductivity; Θ s is the volumetric water content; e is the natural number 2.71828; y is a dummy variable of integration representing the logarithm of negative pore-water pressure; i is the interval between the range of j to N; j is the least negative pore-water pressure to be described by the final function; N is the maximum negative pore-water pressure to be described by the final function; Ψ is the suction corresponding to the jth interval; Θ is the first derivative of the equation. Θ can be described as where a is the air-entry value of the soil; n is a parameter that controls the slope at the inflection point in the volumetric water content function; and m is a parameter that is related to the residual water content. C(Ψ ) is a correcting function defined as where C r is a constant related to the matric suction corresponding to the residual water content. The sandy soil and clay soil were selected for analysis in this study, which represents the high and low permeability [46], as shown in Figure 3. The unsaturated parameter values are shown in Table 1 [47], and the SWCCs are shown in Figure 4. The sandy soil and clay soil were selected for analysis in this study, which represents the high and low permeability [46], as shown in Figure 3. The unsaturated parameter values are shown in Table 1 [47], and the SWCCs are shown in Figure 4.

Physical and Mechanics Parameters
The mechanical parameters of sand and clay were determined according to the laboratory test [47]. The Mohr-Coulomb elastoplastic criterion was applied, and the parameters are shown in Table  2.  The sandy soil and clay soil were selected for analysis in this study, which represents the high and low permeability [46], as shown in Figure 3. The unsaturated parameter values are shown in Table 1 [47], and the SWCCs are shown in Figure 4.

Physical and Mechanics Parameters
The mechanical parameters of sand and clay were determined according to the laboratory test [47]. The Mohr-Coulomb elastoplastic criterion was applied, and the parameters are shown in Table  2.

Physical and Mechanics Parameters
The mechanical parameters of sand and clay were determined according to the laboratory test [47]. The Mohr-Coulomb elastoplastic criterion was applied, and the parameters are shown in Table 2.

Definition of Anisotropy and Calculation Conditions
Previous studies mostly ignored the anisotropy ratio and direction. In fact, anisotropy widely exists in the soil, and for the hydraulic conductivity matrix [C] in Equation (2), it can be expressed as where C 11 = k x cos 2 α + k y sin 2 α, C 22 = k x sin 2 α + k y cos 2 α, and C 12 = C 21 = k x sinαcosα + k y sinαcosα. The k x , k y anisotropy direction α can be defined according to Figure 2; k x is the horizontal permeability coefficient, k y is the vertical permeability coefficient, and α is the direction between the k y and y axes.
Equation (11) is adopted in reference [39], only considering the anisotropy ratio k r = k x /k y ; however, the definition of anisotropy not only includes the ratio k r but also the direction α. However, previous investigations ignored the anisotropy and especially the anisotropy direction α. The anisotropy ratio and direction can be obtained from the consolidation test of soil [32], which can be applied in the real situation.
To completely discuss the anisotropy of sandy soil and clay soil, including the anisotropy ratio k r and the anisotropy direction α, different calculation conditions were used, as shown in Table 3, which includes the anisotropy ratios k r = 1, 10, 50, and 100 and the anisotropy directions α = 0 • , 15 • , 30 • , 45 • , 60 • , 75 • , and 90 • . The reservoir water elevation changed from 175 m to 145 m, the drawdown rate of which was 1.2 m/d. Meanwhile, the 55 d after the reservoir water level drop was also considered, and the total duration was 80 d.

Results
Based on the calculation conditions in Table 3, we carried out a total of 56 numerical simulations. For ease of reading, this section carries out the discussion based on the classification of sandy soil and clay soil slopes. The infiltration line, x-displacement, and safety factors of the 25th day with k r = 10 are displayed to illustrate the influence of anisotropy direction α on seepage, deformation, and stability of the Caipo slope, and those of the 25th day with α = 45 • are also displayed to show the impact of the anisotropy ratio k r .

Sandy Soil
The variation of infiltration lines under different permeability anisotropy degrees k r and directions α of sandy slope is shown in Figure 5. It can be seen that there existed a delay phenomenon in the infiltration line close to the reservoir area, which means that the variation of the infiltration line lagged behind the change of reservoir water level. The anisotropy direction α had a great impact on the seepage characteristics; the infiltration line elevation increased with the increase of α, but when α was more than 45 • , the elevations of different conditions were almost the same.
Based on the calculation conditions in Table 3, we carried out a total of 56 numerical simulations. For ease of reading, this section carries out the discussion based on the classification of sandy soil and clay soil slopes. The infiltration line, x-displacement, and safety factors of the 25th day with kr = 10 are displayed to illustrate the influence of anisotropy direction α on seepage, deformation, and stability of the Caipo slope, and those of the 25th day with α = 45° are also displayed to show the impact of the anisotropy ratio kr.

Sandy Soil
The variation of infiltration lines under different permeability anisotropy degrees kr and directions α of sandy slope is shown in Figure 5. It can be seen that there existed a delay phenomenon in the infiltration line close to the reservoir area, which means that the variation of the infiltration line lagged behind the change of reservoir water level. The anisotropy direction α had a great impact on the seepage characteristics; the infiltration line elevation increased with the increase of α, but when α was more than 45°, the elevations of different conditions were almost the same. The anisotropy ratio kr equally also had a significant influence on the variations of infiltration lines. With the increase of kr, the elevation of the infiltration line increased. When kr was larger than 50, the infiltration line of different conditions was almost the same.

Clay Soil
The variation of infiltration lines under different permeability anisotropy degrees kr and directions α of clay slope is shown in Figure 6.
As can be seen from Figure 6, the infiltration line elevation of the clay slope was higher than that of the sandy slope under the same conditions. What is more, the differences under different kr and α values were smaller than those of the sandy slope. The elevation increased with the increase of kr and α, which was the same as the sandy slope. The anisotropy ratio k r equally also had a significant influence on the variations of infiltration lines. With the increase of k r , the elevation of the infiltration line increased. When k r was larger than 50, the infiltration line of different conditions was almost the same.

Clay Soil
The variation of infiltration lines under different permeability anisotropy degrees k r and directions α of clay slope is shown in Figure 6.
As can be seen from Figure 6, the infiltration line elevation of the clay slope was higher than that of the sandy slope under the same conditions. What is more, the differences under different k r and α values were smaller than those of the sandy slope. The elevation increased with the increase of k r and α, which was the same as the sandy slope.

Analysis of Infiltration Line Hysteretic Elevation (ILHE)
It can be inferred from Section 4.1.1 and 4.1.2 that not only the sandy slope but also the clay slope all experienced the delay phenomenon, which led to a steep increase in the hydraulic gradient in the slope close to the reservoir area and the seepage force's direction out of the slope, thus aggravating the instability. In order to quantitatively describe this phenomenon, the infiltration line hysteretic elevation (ILHE) was defined, as shown in Figure 7a

Analysis of Infiltration Line Hysteretic Elevation (ILHE)
It can be inferred from Sections 4.1.1 and 4.1.2 that not only the sandy slope but also the clay slope all experienced the delay phenomenon, which led to a steep increase in the hydraulic gradient in the slope close to the reservoir area and the seepage force's direction out of the slope, thus aggravating the instability. In order to quantitatively describe this phenomenon, the infiltration line hysteretic elevation (ILHE) was defined, as shown in Figure 7a, as the height of the intersection point between the infiltration line and the toe section of the slope to the 145 m elevation. Figure 7c show the ILHE of sandy and clay slopes, respectively, under different k r and α values.

Analysis of Infiltration Line Hysteretic Elevation (ILHE)
It can be inferred from Section 4.1.1 and 4.1.2 that not only the sandy slope but also the clay slope all experienced the delay phenomenon, which led to a steep increase in the hydraulic gradient in the slope close to the reservoir area and the seepage force's direction out of the slope, thus aggravating the instability. In order to quantitatively describe this phenomenon, the infiltration line hysteretic elevation (ILHE) was defined, as shown in Figure 7a Generally speaking, with the increase of anisotropy ratio k r and direction α, the ILHE of sandy and clay slopes increased to the maximum. This is because as α increased, the permeability coefficient parallel to the x axis was getting smaller and smaller for k r > 1, and the interior water in the slope soil was more difficult to infiltrate into the reservoir, thus leading to the larger ILHE. Meanwhile, the permeability coefficient parallel to the y axis decreased with the increase of k r , and the infiltration water line inside the slope was more difficult to draw down, which also contributed to the increase of ILHE. Compared with the clay slope, the permeability anisotropy characteristics had a greater impact on ILHE of the sandy slope. For the sandy slope, the difference between the maximum and minimum ILHE, only considering k r , was 6.51%, which was 35.53% when considering both k r and α. For the clay slope, the difference only considering k r was 4.7%, which was 9.94% considering both k r and α.

Sandy Soil
The variations of x-displacement of the sandy slope top, slope middle, and slope toe were obtained based on the pore pressure-stress coupling method, which are shown in Figure 8.
had a greater impact on ILHE of the sandy slope. For the sandy slope, the difference between the maximum and minimum ILHE, only considering kr, was 6.51%, which was 35.53% when considering both kr and α. For the clay slope, the difference only considering kr was 4.7%, which was 9.94% considering both kr and α.

Sandy Soil
The variations of x-displacement of the sandy slope top, slope middle, and slope toe were obtained based on the pore pressure-stress coupling method, which are shown in Figure 8.
It can be seen from the figure that the variations of x-displacement at different positions were different. For slope top, the x-displacement increased first, then decreased with the elevation. The maximum x-displacement was about 5 m from the surface. For slope middle and slope toe, the x-displacement decreased with the elevation, and the maximum x-displacement was on the slope surface. The overall x-displacement increased when the distance of the slope to the reservoir decreased, which was because the decrease of reservoir water level led to the delay phenomenon in slope soil, the most dramatic of which was the slope toe, and it caused the dynamic seepage force, which drove the slope soil to move away [6][7][8]. Therefore, we can infer that slope toe first failed with the drawdown of the reservoir water level. The permeability anisotropy characteristics also had a great impact on the variation of x-displacement. The x-displacement value of the slope top was very small compared to the slope middle and slope toe, such that it could be ignored. The x-displacement of slope middle decreased with the increase of kr, and that of the slope toe increased with the increase of kr. This may due to the fact that the increase of kr and α increased the hydraulic gradient in the slope toe but decreased the hydraulic gradient in the slope middle.  It can be seen from the figure that the variations of x-displacement at different positions were different. For slope top, the x-displacement increased first, then decreased with the elevation. The maximum x-displacement was about 5 m from the surface. For slope middle and slope toe, the x-displacement decreased with the elevation, and the maximum x-displacement was on the slope surface. The overall x-displacement increased when the distance of the slope to the reservoir decreased, which was because the decrease of reservoir water level led to the delay phenomenon in slope soil, the most dramatic of which was the slope toe, and it caused the dynamic seepage force, which drove the slope soil to move away [6][7][8]. Therefore, we can infer that slope toe first failed with the drawdown of the reservoir water level. The permeability anisotropy characteristics also had a great impact on the variation of x-displacement. The x-displacement value of the slope top was very small compared to the slope middle and slope toe, such that it could be ignored. The x-displacement of slope middle decreased with the increase of k r , and that of the slope toe increased with the increase of k r . This may due to the fact that the increase of k r and α increased the hydraulic gradient in the slope toe but decreased the hydraulic gradient in the slope middle.

Clay Soil
The variations of x-displacement of clay slope top, slope middle and slope toe are shown in Figure 9.

Clay Soil
The variations of x-displacement of clay slope top, slope middle and slope toe are shown in Figure 9.  As can be seen from Figure 9, the variations of x-displacement under different kr and α values for clay slope were similar to those for the sandy slope. What is to be noticed is that the influence of permeability anisotropy on the clay slope was much smaller than that of the sandy slope. The x-displacement for the clay slope in the slope top and slope middle was smaller than that of the sandy slope, which was larger in the slope toe. This is easy to understand that for the clay slope, the dramatic change of infiltration line concentrated in the slope toe, but rarely changed in the slope middle and slope top. For the sandy slope, the reservoir water drawdown not only influenced the slope toe but also had some impacts on the slope middle.

Analysis of Maximum Horizontal Displacement (MHD)
For the purpose of quantitative study, the maximum horizontal displacement (MHD) for sandy and clay slopes under different permeability anisotropy characteristics and the variations of MHD under different conditions are shown in Figure 10.
The impact of permeability anisotropy on MHD changed with the positions of slope sections (slope top, slope middle, and slope toe). The MHD was at the maximum in the slope toe, was smaller in the slope middle, and was the smallest in the slope top. For the sandy slope, the difference between the maximum and minimum MHD, only considering kr for the slope top, slope middle, and slope toe was 69.5%, 20.7%, and 6.9%, respectively, while it became 75.05%, 77.2%, and 21.4% for the clay slope. Taking both kr and α into consideration, the difference was 80.0%, 94.6%, and 45.0% for the sandy slope, while it became 75.1%, 87.1%, and 32.0% for the clay slope. As can be seen from Figure 9, the variations of x-displacement under different k r and α values for clay slope were similar to those for the sandy slope. What is to be noticed is that the influence of permeability anisotropy on the clay slope was much smaller than that of the sandy slope. The x-displacement for the clay slope in the slope top and slope middle was smaller than that of the sandy slope, which was larger in the slope toe. This is easy to understand that for the clay slope, the dramatic change of infiltration line concentrated in the slope toe, but rarely changed in the slope middle and slope top. For the sandy slope, the reservoir water drawdown not only influenced the slope toe but also had some impacts on the slope middle.

Analysis of Maximum Horizontal Displacement (MHD)
For the purpose of quantitative study, the maximum horizontal displacement (MHD) for sandy and clay slopes under different permeability anisotropy characteristics and the variations of MHD under different conditions are shown in Figure 10.

Sandy Soil
The variations of safety factors (SFs) for the sandy slope under different permeability anisotropy degrees kr and directions α are shown in Figure 11. The SF first decreased then increased with the reservoir water draw down; what should be stressed is that the permeability anisotropy not only influenced the value of the SF but also had some impact on the critical time when it reached the minimum safety factor (MSF). The SF decreased with the increase of kr and α, which was closely related to the delay phenomenon discussed in Section 4.1. Meanwhile, the increase of kr and α made the critical time of the MSF increase, due mainly to three reasons. One was the discharge of the water pressure acting on the slope surface, which reduced the anti-sliding force and was regarded as an unfavorable factor for slope stability. The second reason was the delay phenomenon of the infiltration line inside the slope soil, which contributed to the sharp increase in the hydraulic gradient in the slope toe and was also regarded as an unfavorable factor. The last reason was the increase of effective stress and soil strength with the decrease of infiltration line, which was regarded as a positive factor. Combined with the analysis in Section 4.1.1, the smaller the α and kr, the more easily the interior water in the slope soil infiltrated into the reservoir, and the influence of water pressure unloading and the hydraulic gradient increase was less than the soil strength increase, which led to the earlier critical time. This was the reverse when α and kr increased and the critical time was delayed. The impact of permeability anisotropy on MHD changed with the positions of slope sections (slope top, slope middle, and slope toe). The MHD was at the maximum in the slope toe, was smaller in the slope middle, and was the smallest in the slope top. For the sandy slope, the difference between the maximum and minimum MHD, only considering k r for the slope top, slope middle, and slope toe was 69.5%, 20.7%, and 6.9%, respectively, while it became 75.05%, 77.2%, and 21.4% for the clay slope. Taking both k r and α into consideration, the difference was 80.0%, 94.6%, and 45.0% for the sandy slope, while it became 75.1%, 87.1%, and 32.0% for the clay slope.

Sandy Soil
The variations of safety factors (SFs) for the sandy slope under different permeability anisotropy degrees k r and directions α are shown in Figure 11. The SF first decreased then increased with the reservoir water draw down; what should be stressed is that the permeability anisotropy not only influenced the value of the SF but also had some impact on the critical time when it reached the minimum safety factor (MSF). The SF decreased with the increase of k r and α, which was closely related to the delay phenomenon discussed in Section 4.1. Meanwhile, the increase of k r and α made the critical time of the MSF increase, due mainly to three reasons. One was the discharge of the water pressure acting on the slope surface, which reduced the anti-sliding force and was regarded as an unfavorable factor for slope stability. The second reason was the delay phenomenon of the infiltration line inside the slope soil, which contributed to the sharp increase in the hydraulic gradient in the slope toe and was also regarded as an unfavorable factor. The last reason was the increase of effective stress and soil strength with the decrease of infiltration line, which was regarded as a positive factor. Combined with the analysis in Section 4.1.1, the smaller the α and k r , the more easily the interior water in the slope soil infiltrated into the reservoir, and the influence of water pressure unloading and the hydraulic gradient increase was less than the soil strength increase, which led to the earlier critical time. This was the reverse when α and k r increased and the critical time was delayed.

Clay Soil
The variations of safety factors (SFs) for the clay slope under different permeability anisotropy degrees kr and directions α are shown in Figure 12. The variations of SFs under different kr and α values for the clay slope were different in that the SF first decreased then remained unchanged. What is more, the influence of permeability anisotropy on the clay slope was much smaller than that of the sandy slope. What is to be noticed is that the SF of the clay slope was smaller than that of the sandy slope under the same conditions, which was because of the more obvious delay phenomenon in the clay slope.   Figure 13. The MSF decreased to the minimum with the increase of kr and α. What is to be stressed is that the overall MSF of the clay slope was smaller than that of the sandy slope. The difference between the maximum and minimum MSF, only considering kr, was 2.72% for the sandy slope, which was 9.29% for the clay slope. Taking both kr and α into consideration, the difference was 11.58% and 12.31% for the sandy and clay slopes, respectively.

Clay Soil
The variations of safety factors (SFs) for the clay slope under different permeability anisotropy degrees k r and directions α are shown in Figure 12. The variations of SFs under different k r and α values for the clay slope were different in that the SF first decreased then remained unchanged. What is more, the influence of permeability anisotropy on the clay slope was much smaller than that of the sandy slope. What is to be noticed is that the SF of the clay slope was smaller than that of the sandy slope under the same conditions, which was because of the more obvious delay phenomenon in the clay slope.

Clay Soil
The variations of safety factors (SFs) for the clay slope under different permeability anisotropy degrees kr and directions α are shown in Figure 12. The variations of SFs under different kr and α values for the clay slope were different in that the SF first decreased then remained unchanged. What is more, the influence of permeability anisotropy on the clay slope was much smaller than that of the sandy slope. What is to be noticed is that the SF of the clay slope was smaller than that of the sandy slope under the same conditions, which was because of the more obvious delay phenomenon in the clay slope.   Figure 13. The MSF decreased to the minimum with the increase of kr and α. What is to be stressed is that the overall MSF of the clay slope was smaller than that of the sandy slope. The difference between the maximum and minimum MSF, only considering kr, was 2.72% for the sandy slope, which was 9.29% for the clay slope. Taking both kr and α into consideration, the difference was 11.58% and 12.31% for the sandy and clay slopes, respectively.  Figure 13. The MSF decreased to the minimum with the increase of k r and α. What is to be stressed is that the overall MSF of the clay slope was smaller than that of the sandy slope. The difference between the maximum and minimum MSF, only considering k r , was 2.72% for the sandy slope, which was 9.29% for the clay slope. Taking both k r and α into consideration, the difference was 11.58% and 12.31% for the sandy and clay slopes, respectively.  Figure 14 shows the dynamic variation process of the infiltration line of the sandy slope with kr = 1 and α = 0°, which is regarded as isotropic, as adopted by most of the previous studies. Under the action of reservoir water drawdown, the drawdown of the infiltration line lags behind that of the reservoir water level, which is due to the holding water capacity of soils. Therefore, the impact of reservoir water drawdown on the slope mainly reflects the following three aspects: (1) The delay phenomenon leads to higher pore water pressure (excess pore water pressure), which cannot be dissipated in time and forms an unsteady seepage field. The excess pore water pressure contributes to the sharp increase in the hydraulic gradient in the slope toe and is regarded as an unfavorable factor. (2) The pore pressure acting on the slope surface prevents the slope from sliding when the reservoir water level is high, but the reservoir water drawdown contributes to the unloading of the water pressure, which is also regarded as an unfavorable factor. (3) The drawdown of the reservoir water level results in the decrease of the pore water pressure interior of the slope, which increases the soil strength and the effective stress and is a positive factor. Therefore, the slope stability is affected by the above three factors, and the slope soil as well as the permeability anisotropy are important contributing factors. Meanwhile, it is found that the x-displacement of the slope toe is the largest during the drawdown of the reservoir water level, which is also due to the delay phenomenon. In this paper, the pore pressure-stress coupling method is utilized to obtain the variation of x-displacement contour, as is shown in Figure 15a. The x-displacement in the slope toe is the largest on the 5 th day. With the further drawdown of reservoir water level, the x-displacement in the slope middle and slope top gradually gets bigger, which illustrates that the slope toe firstly fails during the reservoir water drawdown. Jiang et al. [48] conducted a series of model tests of landslides and found that the  Figure 14 shows the dynamic variation process of the infiltration line of the sandy slope with k r = 1 and α = 0 • , which is regarded as isotropic, as adopted by most of the previous studies. Under the action of reservoir water drawdown, the drawdown of the infiltration line lags behind that of the reservoir water level, which is due to the holding water capacity of soils. Therefore, the impact of reservoir water drawdown on the slope mainly reflects the following three aspects: (1) The delay phenomenon leads to higher pore water pressure (excess pore water pressure), which cannot be dissipated in time and forms an unsteady seepage field. The excess pore water pressure contributes to the sharp increase in the hydraulic gradient in the slope toe and is regarded as an unfavorable factor. (2) The pore pressure acting on the slope surface prevents the slope from sliding when the reservoir water level is high, but the reservoir water drawdown contributes to the unloading of the water pressure, which is also regarded as an unfavorable factor. (3) The drawdown of the reservoir water level results in the decrease of the pore water pressure interior of the slope, which increases the soil strength and the effective stress and is a positive factor. Therefore, the slope stability is affected by the above three factors, and the slope soil as well as the permeability anisotropy are important contributing factors.  Figure 14 shows the dynamic variation process of the infiltration line of the sandy slope with kr = 1 and α = 0°, which is regarded as isotropic, as adopted by most of the previous studies. Under the action of reservoir water drawdown, the drawdown of the infiltration line lags behind that of the reservoir water level, which is due to the holding water capacity of soils. Therefore, the impact of reservoir water drawdown on the slope mainly reflects the following three aspects: (1) The delay phenomenon leads to higher pore water pressure (excess pore water pressure), which cannot be dissipated in time and forms an unsteady seepage field. The excess pore water pressure contributes to the sharp increase in the hydraulic gradient in the slope toe and is regarded as an unfavorable factor. (2) The pore pressure acting on the slope surface prevents the slope from sliding when the reservoir water level is high, but the reservoir water drawdown contributes to the unloading of the water pressure, which is also regarded as an unfavorable factor. (3) The drawdown of the reservoir water level results in the decrease of the pore water pressure interior of the slope, which increases the soil strength and the effective stress and is a positive factor. Therefore, the slope stability is affected by the above three factors, and the slope soil as well as the permeability anisotropy are important contributing factors. Meanwhile, it is found that the x-displacement of the slope toe is the largest during the drawdown of the reservoir water level, which is also due to the delay phenomenon. In this paper, the pore pressure-stress coupling method is utilized to obtain the variation of x-displacement contour, as is shown in Figure 15a. The x-displacement in the slope toe is the largest on the 5 th day. With the further drawdown of reservoir water level, the x-displacement in the slope middle and slope top gradually gets bigger, which illustrates that the slope toe firstly fails during the reservoir water drawdown. Jiang et al. [48] conducted a series of model tests of landslides and found that the Meanwhile, it is found that the x-displacement of the slope toe is the largest during the drawdown of the reservoir water level, which is also due to the delay phenomenon. In this paper, the pore pressure-stress coupling method is utilized to obtain the variation of x-displacement contour, as is shown in Figure 15a. The x-displacement in the slope toe is the largest on the 5 th day. With the further drawdown of reservoir water level, the x-displacement in the slope middle and slope top gradually gets bigger, which illustrates that the slope toe firstly fails during the reservoir water drawdown. Jiang et al. [48] conducted a series of model tests of landslides and found that the model slope failed from local slip failure to leading-edge failure under declining water levels, showing a progressive retrogressive mode, which is consist with our numerical results and verifies the correctness and rationality of the pore pressure-stress coupling method. Therefore, in the actual reservoir area landslide treatment project, we should strengthen the reinforcement measures at the slope toe of the landslide.

Effect of Reservoir Water Level Drawdown on Seepage Deformation and Stability of Slope
Water 2020, 11, x FOR PEER REVIEW 16 of 20 model slope failed from local slip failure to leading-edge failure under declining water levels, showing a progressive retrogressive mode, which is consist with our numerical results and verifies the correctness and rationality of the pore pressure-stress coupling method. Therefore, in the actual reservoir area landslide treatment project, we should strengthen the reinforcement measures at the slope toe of the landslide.

Effect of Permeability Anisotropy on Seepage Deformation and Stability of Slope
It is also found that the permeability anisotropy has a great influence on the seepage, deformation, and stability of the slope. Previous studies only took kr into consideration but ignored the effect of α. Table 4 shows the statistics of the differences between only considering kr and considering both kr and α compared with the isotropy condition. As can be seen from Table 4, the values of ILHE, MHD, and MSF for the sandy slope are quite different under the conditions with only considering kr and with considering both kr and α. Therefore,

Effect of Permeability Anisotropy on Seepage Deformation and Stability of Slope
It is also found that the permeability anisotropy has a great influence on the seepage, deformation, and stability of the slope. Previous studies only took k r into consideration but ignored the effect of α. Table 4 shows the statistics of the differences between only considering k r and considering both k r and α compared with the isotropy condition. As can be seen from Table 4, the values of ILHE, MHD, and MSF for the sandy slope are quite different under the conditions with only considering k r and with considering both k r and α. Therefore, for the sandy slope, which represents the kind of soil slope with larger permeability, k r and α must be considered in the seepage, deformation, and stability analysis. However, for the clay slope, which represents the kind of soil slope with lower permeability, the differences are much smaller. Thus, we can infer that for soil with smaller permeability, we can consider only the effect of k r but ignore α to simplify the calculation.

Advice for Future Studies
In our analysis, the 2D model was applied to carry out the sensitivity analysis of the effect of the permeability anisotropy characteristics on the seepage, deformation, and stability of a slope, which has the advantage of simply and straightforwardly illustrating the rules. The results show that the effect of k r and α is huge and can determine the overall slope behavior, which provides a reference for the analysis of a 3D model in real cases.

Conclusions
In this paper, the numerical simulation on the seepage, deformation, and stability of a slope is carried out based on the pore pressure-stress coupling model, and the Capo slope in the TGR is taken as a case study with consideration of the effect of permeability anisotropy ratio k r and direction α, which includes the analysis of two kinds of soil type, namely a sandy slope and clay slope. The following conclusions can be drawn: 1.
The increase of k r and α decreases the infiltration capacity of slope soil, which increases the pore water pressure, as well as the x-displacement, and decreases the SF. The delay phenomenon occurs during the drawdown of the reservoir water level. 2.
The infiltration line hysteretic elevation (ILHE), the maximum horizontal displacement (MHD), and the minimum safety factor (MSF) are defined to characterize the delay phenomenon, the deformation, and the stability. The values of ILHE, MHD, and MSF increase with the increase of k r and α. The ILHE of the sandy slope is smaller than that of the clay slope, but the MHD and MSF of the sandy slope is larger than that of the clay slope.

3.
There are mainly three aspects of the impact of reservoir water drawdown on the slope, which include two unfavorable factors containing the delay phenomenon and the water pressure unload and one positive factor containing the increase of soil strength and effective stress. The x-displacement reaches a maximum in the slope toe, which indicates the slope toe firstly fails during the reservoir water drawdown, and it is consistent with the test results in the existing literature. 4.
The differences of ILHE, MHD, and MSF between only considering k r as well as both considering k r and α are statistically analyzed. For the sandy slope, the values of ILHE, MHD, and MSF are quite different under the conditions with only considering k r and with both considering k r and α, while they become relatively small for the clay slope. Therefore, we must consider k r and α for the sandy slope; however, for the clay slope, we can consider only k r to simplify computation. 5.
In this paper, the permeability anisotropy ratio k r and direction α are only considered in our simulation; future studies should combine the permeability anisotropy with the soil strength anisotropy in numerical analysis and develop new corresponding anisotropic material and carry out the model tests to verify the numerical simulation results.
Author Contributions: S.Y. completed the statistical analysis and wrote the paper; X.R., J.Z., and H.W. provided the writing ideas and supervised the study; J.W. and W.Z. did the numerical analysis. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Fund, grant number U1765204.

Conflicts of Interest:
The authors declare no conflict of interest.