Sensitivity Analyses of the Seepage and Stability of Layered Rock Slope Based on the Anisotropy of Hydraulic Conductivity: A Case Study in the Pulang Region of Southwestern China

: In the study of the seepage characteristics of layered rock slope under rainfall conditions, the majority of previous research has considered the hydraulic conduction to be isotropic, or only considered the anisotropy ratio of the hydraulic conductivity, ignoring the anisotropy angle. In the current study, a layered rock slope in the Pulang region was selected as an example. Then, based on the ﬁtting parameters of the Van Genuchten model, pore water pressure sensitivity analyses of the layered rock slope were carried out. The anisotropy ratio and anisotropy angle were used to analyze the sensitivity of the seepage and stability of the layered rock slopes. The results show that as the anisotropy angle of hydraulic conductivity of layered rock slope decreased, the maximum volume water content of surface (MWCS) of layered rock slope gradually increased. Additionally, as the anisotropy ratio decreased and the anisotropy angle increased, the rising heights of the groundwater (RHG) of layered rock slope gradually increased. When the hydraulic conduction of layered rock slope was considered isotropic, the factor of safety (FS) tended to be overestimated. As the anisotropy ratio decreased and the anisotropy angle increased, the factor of safety (FS) of layered rock slope decreased. Prevention should be the objective for rock slopes with larger dip angles in the bedding plane in the Pulang region. This study provides feasible schemes for the evaluation of the seepage and stability of layered rock slopes in Pulang region of southwestern China.


Introduction
According to available statistical data, more than 90% of rock slope failures are related to fluid hydraulic conduction [1]. Collins and Znidarcic [2] pointed out that the depths of slope failures are not only controlled by strength characteristics, but also by hydraulic characteristics, and both should be used to assess the stability of slopes. Dong et al. [3] found out that the anisotropic hydraulic conduction of geotechnical material affects slope pore water pressure distributions and factors of safety.
Water 2020, 12 Dong et al. [3] also concluded that the use of isotropic and hydrostatic assumptions to estimate pore water pressure levels will produce erroneous results. The majority of previous anisotropy studies of layered rock slopes have focused on the anisotropy of rock strength. For example, Zbek et al. [4] conducted anisotropy ratio and anisotropy angle studies on the strength of four types of metamorphic rock. In a study by Lee et al. [5], the anisotropy ratio and anisotropy angle of slate strength were taken into consideration, and excavation deformations of slate tunnels were successfully simulated. Current anisotropy studies regarding the hydraulic conduction of rock masses are generally focused on rock masses with joint fissures [6][7][8]. The results have shown that the deformations and hydraulic conduction of layered rock masses display anisotropic characteristics [9]. Yeh et al. [10] found that layered rock masses exhibit heterogeneity in hydraulic conduction and strength. It was considered that the effects of anisotropy of hydraulic conductivity and strength on groundwater flow and slope stability were mainly reflected in the direction of the bedding plane. Dong et al. [11] also investigated the anisotropy ratio of hydraulic conductivity and strength; the authors studied the stability of layered and weakly cemented sandstone slopes, applying two-dimensional seepage models and steady-state seepage analysis methods to determine that the anisotropy ratio of the bedding plane had significant effects on the pore pressure distributions and slope stability. The ratio of the parallel to the vertical hydraulic conductivity of the bedding plane in the western foothills of Taiwan was 10 to 100 [10]. However, the aforementioned studies did not examine the anisotropy angles of the hydraulic conductivity. The anisotropy angles of the rock layers are known to be closely related to the dip angles of the bedding plane. Yu et al. [12] pointed out that the horizontal hydraulic conductivity and the vertical conductivity only coincide with the natural coordinate axis in some special cases. It has been observed that in nature, there are many examples where the anisotropic angle does not coincide with the coordinate axis. Therefore, research on the seepage characteristics of layered rock slope needs to consider the dip angles of the rock layers.
The hydraulic conductivity in rock masses is not a constant value and can be affected by saturation degree [13]. In particular, for highly weathered rock masses, due to the loose porous media, saturated-unsaturated states will significantly affect the hydraulic conductivity. The unsaturated seepage of the rock masses has been reflected in many previous studies. However, the majority of those studies have focused on the unsaturated flow of rock fractures [6,14,15], and the unsaturated flow of the rock layers was rarely involved. Since the unsaturated parameters of rock layers are difficult to measure under laboratory conditions, classic unsaturated porous media water retention curves are generally substituted [16][17][18], along with relative hydraulic conduction models [19]. Additionally, determinations of the water retention curves of unsaturated rock and porous media tend to be divided into parameter sensitivity analysis and inversion analysis processes. Yu et al. [20] studied the effects of different unsaturated porous media parameters on the seepage characteristics and stability levels of slopes under different types of rain events using the Fredlund and Xing models. Chen et al. [6] utilized orthogonal designs, finite element forward analysis, artificial neural networks, and genetic algorithms to perform inverse analyses on unsaturated rock masses; it was found that the unsaturated parameters of the Van Genuchten model, with the calculated flow and water levels, were in good agreement with the field observations.
Because of the shortcomings of previous studies, firstly, a sensitivity analysis of pore water pressure based on the parameters of the Van Genuchten model, was performed. Compared with the field monitoring data, a reasonable water retention curve of the layered rock mass was produced. Secondly, a reasonable initial matric suction was determined, taking a strong weathered layered rock slope in the Pulang region, China as an example. Finally, this study analyzed the sensitivity of the seepage and stability of the layered rock slope under different anisotropy ratios and anisotropy angles of hydraulic conductivity. This was specifically ascertained from the analysis results of the maximum volume water content of surface (MWCS), rising heights of the groundwater (RHG) and the factors of safety (FS). Field investigations verified the feasibility of numerical simulation. The research results Water 2020, 12, 2314 3 of 20 provide some references for the understanding of the seepage anisotropy law and the prevention of rock landslides.

Control Differential Equation
To simulate the seepage processes of rock slopes, rainfall infiltration was regarded as the variable flow boundary of the unsaturated zone, and a saturated-unsaturated seepage theory was applied to perform the simulations. Under the principle of mass conservation, the saturated-unsaturated seepage control differential equation was as follows [21]: where H is the total head, the unit is m; x and y are the coordinates in the direction of x and y; k x is the hydraulic conductivity in the x-direction, the unit is m/s; k y is the hydraulic conductivity in the y-direction, the unit is m/s; Q is the applied boundary flux, m w is the slope of the storage curve, t is the time, the unit is s; γ w is the unit weight of water, the unit is N/m 3 .

Factors of Safety for Unsaturated Layered Rock Slope
This study adopted a Morgenstern-Price Method based on a limit equilibrium theory to calculate the factors of safety. Since Morgenstern and Price [22] developed a method for the analysis of rock slope stability which offered various user-specified inter-slice force functions and also considered both shear and normal inter-slice forces, these properties make it possible to use for the rock slope stability problem with the actual situation of complex engineering characteristics. Therefore, Morgenstern-Price method is considered, for many rock slope stability analyses, a more reasonable and comprehensive analysis method [13,23,24]. It was found that the modified method strictly satisfied the force balance and torque balance, and the calculation accuracy was found to be high. The expressions were as follows [22,25]: where u a is the pore-air pressure and the unit is kPa; u w is the pore-water pressure and the unit is kPa; c i is the effective cohesion for every rock slice and the unit is kPa; i indicates the rock slice number; W i represents the weight of every rock slice; the unit is N/m 3 ; P i indicates the water pressure and the unit is kPa; ϕ b is an angle defining the increase in shear strength for an increase in suction, β i is the angle of the bottom of the rock slice; b i denotes the length of every rock slice and the unit is m; ϕ i represents the shear strength angle for every rock slice; r i is the radius of the sliding arc; F s denotes the factors of safety.

Finite Element Model and Boundary Conditions
The case study is a layered rock slope in the Pulang Region of Southwestern China, which is shown in Figure 1. In Figure 2, I and II represent strongly weathered carbonaceous slate and moderately weathered carbonaceous slate, respectively, which were characterized by bedding structures. This study's finite element model of the slope body was 240 m high and 280 m long. The left and right borders extended horizontally for 50 m, and the lower border extended downward for 50 m, to eliminate any border effects. To investigate the strongly weathered layered slate, the grid size of the strongly weathered area was established as 2 m, and the grid size of the moderately weathered area was set as 10 m. The grid model was divided into 6333 nodes and 6208 elements, as shown in Figure 2. To study the characteristics of the seepage at different positions, three monitoring lines and three monitoring points were set at x = 150 m (top of the slope), x = 250 m (middle of the slope), and x = 330 m (bottom of the slope), and the monitoring points were located 10 m below the slope. In Figure 2, the boundary conditions are as follows: AB and FD are the fixed water level boundaries of 200 m and 90 m, respectively; CD indicates the boundary of the rainfall infiltration; BC and DE are the boundaries of small liquidity; AF is an impervious boundary. To restore the rainfall boundary, the average daily rainfall from 1 July to 1 August 2019 in the Pulang Region was adopted as the rainfall conditions, as detailed in Figure 3.
Water 2020, 12, x FOR PEER REVIEW 4 of 20 slope), and x = 330 m (bottom of the slope), and the monitoring points were located 10 m below the slope. In Figure 2, the boundary conditions are as follows: AB and FD are the fixed water level boundaries of 200 m and 90 m, respectively; CD indicates the boundary of the rainfall infiltration; BC and DE are the boundaries of small liquidity; AF is an impervious boundary. To restore the rainfall boundary, the average daily rainfall from 1 July to 1 August 2019 in the Pulang Region was adopted as the rainfall conditions, as detailed in Figure 3.

Determination of the Maximum Initial Matric Suction
Since the studied rock slope has a strong degree of weathering and the rock layers are so intensely fractured that the properties of soil or a porous medium can be assumed, it has a certain degree of matric suction. The extent of the initial maximum matric suction may vary, and the seepage states and stability levels of slopes may also be quite different. Therefore, to make the initial matric suction more in line with the actual situation (refer to previous studies [25,26]), this study set different maximum initial matric suctions as −50 kPa, −75 kPa, and −100 kPa, respectively, and the annual average rainfall of the Pulang area as 0.06 mm/h in the steady-state comparative analysis process. As

Determination of the Maximum Initial Matric Suction
Since the studied rock slope has a strong degree of weathering and the rock layers are so intensely fractured that the properties of soil or a porous medium can be assumed, it has a certain degree of matric suction. The extent of the initial maximum matric suction may vary, and the seepage states and stability levels of slopes may also be quite different. Therefore, to make the initial matric suction more in line with the actual situation (refer to previous studies [25,26]), this study set different maximum initial matric suctions as −50 kPa, −75 kPa, and −100 kPa, respectively, and the annual average rainfall of the Pulang area as 0.06 mm/h in the steady-state comparative analysis process. As

Determination of the Maximum Initial Matric Suction
Since the studied rock slope has a strong degree of weathering and the rock layers are so intensely fractured that the properties of soil or a porous medium can be assumed, it has a certain degree of matric suction. The extent of the initial maximum matric suction may vary, and the seepage states and stability levels of slopes may also be quite different. Therefore, to make the initial matric suction more in line with the actual situation (refer to previous studies [25,26]), this study set different maximum initial matric suctions as −50 kPa, −75 kPa, and −100 kPa, respectively, and the annual average rainfall of the Pulang area as 0.06 mm/h in the steady-state comparative analysis process. As detailed in Figure 4, the changes in the matric suction of the average annual rainfall were similar to the setting of the maximum initial matric suction of −75 kPa. Therefore, the initial maximum matric suction of all this study's simulation experiments was set as −75 kPa, which was determined to be in line with the real situation in the study area.
Water 2020, 12, x FOR PEER REVIEW 6 of 20 detailed in Figure 4, the changes in the matric suction of the average annual rainfall were similar to the setting of the maximum initial matric suction of −75 kPa. Therefore, the initial maximum matric suction of all this study's simulation experiments was set as −75 kPa, which was determined to be in line with the real situation in the study area.

Van Genuchten Model Parameter Sensitivity Analysis
Due to the difficulties encountered when measuring the unsaturated characteristic curves in rock masses, the water retention curves of classic unsaturated porous media [16][17][18]27] and relative hydraulic conduction models [19] are often used to describe the unsaturated flow in fissures and weakly permeable rock masses. Therefore, this study adopted the Van Genuchten Model to describe unsaturated features. The Van Genuchten Model is as follows [18]: where p is the suction, unit is kPa; θ indicates the adjusted volumetric water content, unit is m −1 ; θr is the residual volumetric water content, unit is m −1 ; θs indicates the saturated volumetric water content, unit is m −1 ;.a is the fitting parameter closely related to the air-entry value of the unsaturated rock mass, and the unit is kPa; n and m (m = 1 − 1/n and n > 1)are fitting parameters that control the slope at the inflection point in the volumetric water content function [18]. kw is the saturated hydraulic conductivity, unit is m/s; k is the adjusted hydraulic conductivity, unit is m/s. Parameters such as a, n, m, and θr are known to have major influences on the unsaturated flow of rock slope. For example, based on the analysis results of the seepage and stability of different slopes under the Fredlund and Xing parameters, Yu et al. [20] determined that the fitting parameters a, n, and m have major influences on the seepage characteristics of the slope. Additionally, Chen et al. [6] found that the unsaturated fitting parameters of the VG model were in the range of a = 0.2 to 100, m = 0.31 to 0.69, and n = 1.45 to 3.23 for rock masses. In a study by Jiang et al. [28], for general soil, the fitting parameters a = 100 and n = 2 were adopted in the Geostudio software, and for rock mass, the fitting parameter a = 10.
Therefore, this study referred to the recommendations of previous, related studies [6,18,20,28] and divided the parameters into three groups, as detailed in Table 1. Since m = (1 − 1/n), this study chose to omit discussion of m.

Van Genuchten Model Parameter Sensitivity Analysis
Due to the difficulties encountered when measuring the unsaturated characteristic curves in rock masses, the water retention curves of classic unsaturated porous media [16][17][18]27] and relative hydraulic conduction models [19] are often used to describe the unsaturated flow in fissures and weakly permeable rock masses. Therefore, this study adopted the Van Genuchten Model to describe unsaturated features. The Van Genuchten Model is as follows [18]: where p is the suction, unit is kPa; θ indicates the adjusted volumetric water content, unit is m −1 ; θ r is the residual volumetric water content, unit is m −1 ; θ s indicates the saturated volumetric water content, unit is m −1 ; a is the fitting parameter closely related to the air-entry value of the unsaturated rock mass, and the unit is kPa; n and m (m = 1 − 1/n and n > 1)are fitting parameters that control the slope at the inflection point in the volumetric water content function [18]. k w is the saturated hydraulic conductivity, unit is m/s; k is the adjusted hydraulic conductivity, unit is m/s. Parameters such as a, n, m, and θ r are known to have major influences on the unsaturated flow of rock slope. For example, based on the analysis results of the seepage and stability of different slopes under the Fredlund and Xing parameters, Yu et al. [20] determined that the fitting parameters a, n, and m have major influences on the seepage characteristics of the slope. Additionally, Chen et al. [6] found that the unsaturated fitting parameters of the VG model were in the range of a = 0.2 to 100, m = 0.31 to 0.69, and n = 1.45 to 3.23 for rock masses. In a study by Jiang et al. [28], for general soil, the fitting parameters a = 100 and n = 2 were adopted in the Geostudio software, and for rock mass, the fitting parameter a = 10. Therefore, this study referred to the recommendations of previous, related studies [6,18,20,28] and divided the parameters into three groups, as detailed in Table 1. Since m = (1 − 1/n), this study chose to omit discussion of m. In Figure 5a, parameter a represents the quantity closely related to the air-entry value of the unsaturated rock mass; the smaller the value of a, the more easily the water in the pores of the rock mass will be discharged. At the same time, due to the entry of air into the pores of the rock body, the water in the rock body will have difficulty flowing. Therefore, with decreases in the fitting parameter a, the hydraulic conductivity will drop earlier, which is shown in Figure 5d. As detailed in Figure 5b, the fitting parameter n controls the slope at the inflection point in the volumetric water content function. Therefore, with increases in parameter n, the slope part of the hydraulic conductivity of the rock will also increase, as shown in Figure 5e. It can be seen in Figure 5f that the residual volume of the water content has only a small effect on the hydraulic conductivity of rock. Therefore, θ r can take any value from 0.001 to 0.05, so θ r is set as 0.001 in this study.  In Figure 5a, parameter a represents the quantity closely related to the air-entry value of the unsaturated rock mass; the smaller the value of a, the more easily the water in the pores of the rock mass will be discharged. At the same time, due to the entry of air into the pores of the rock body, the water in the rock body will have difficulty flowing. Therefore, with decreases in the fitting parameter a, the hydraulic conductivity will drop earlier, which is shown in Figure 5d. As detailed in Figure 5b, the fitting parameter n controls the slope at the inflection point in the volumetric water content function. Therefore, with increases in parameter n, the slope part of the hydraulic conductivity of the rock will also increase, as shown in Figure 5e. It can be seen in Figure 5f that the residual volume of the water content has only a small effect on the hydraulic conductivity of rock. Therefore, θr can take any value from 0.001 to 0.05, so θr is set as 0.001 in this study. Field monitoring data of pore water pressure were measured by the vibrating string type pore water pressure gauge (BGK10-BGK-4500S). The time and positions of field monitoring measurement are consistent with time and position the numerical simulation. It can be seen from Figure 6a,c that during the period of rainfall, the closer the fitting parameter a is to 10, the closer to field monitoring data the pore water pressure from the numerical simulation is. From Figure 6d,f, it can be seen that the closer the fitting parameter n is to 1.5, the closer to field monitoring data the pore water pressure from the numerical simulation is. So the fitting parameters of the water retention curve of the rock masses were determined, as detailed in Table 2. Accordingly, the water retention curve of the rock mass is shown in Figure 7. Field monitoring data of pore water pressure were measured by the vibrating string type pore water pressure gauge (BGK10-BGK-4500S). The time and positions of field monitoring measurement are consistent with time and position the numerical simulation. It can be seen from Figure 6a,c that during the period of rainfall, the closer the fitting parameter a is to 10, the closer to field monitoring data the pore water pressure from the numerical simulation is. From Figure 6d,f, it can be seen that the closer the fitting parameter n is to 1.5, the closer to field monitoring data the pore water pressure from the numerical simulation is. So the fitting parameters of the water retention curve of the rock masses were determined, as detailed in Table 2. Accordingly, the water retention curve of the rock mass is shown in Figure 7.

Definition of Anisotropy and the Calculation Conditions
It was found that the previous related studies had generally ignored the anisotropy ratio and angle. However, anisotropy widely exists in rock masses. In this study, the hydraulic conductivity matrix [C] was expressed as follows: In Equation (5), 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 and k y , as well as the anisotropy angle α, were defined according to Figure 2. In the present study, k x represents the horizontal hydraulic conductivity; k y is the vertical hydraulic conductivity; α is the direction between k x and x-axis. Therefore, if α = 0 • , then [C] will be reduced to: Equation (6) was adopted in this study [29] with only the anisotropy ratio k r = k y /k x considered. However, for the layered rock masses, due to the dip angles of the bedding plane, the conditions under which the anisotropy angle α is not equal to 0 should also be considered. Therefore, in this study, to restore a reasonable rock formation state, a dip angle range of between 50 • and 70 • in the layered rock slope in the Pulang area was taken into consideration. As shown in Figure 2, the anisotropic angle of hydraulic conductivity is equivalent to the dip angle of the layered rock slope. Combining the recommendations of previous related research [6,12,18,[30][31][32], it was determined that the anisotropy ratio k r was 0.01, 0.02, 0.1, and 1, and the anisotropy angle α was 50 • , 55 • , 60 • , 65 • , and 70 • , respectively, neglecting the anisotropy of the moderately weathered carbonaceous slate, as detailed in Table 3. The failure criterion of the rock layered slope simulation adopted the Mohr-Coulomb Criterion, and the rock mass strength parameters were obtained from the geotechnical test results, as shown in Table 4.

Analysis of the Volumetric Water Content
According to the calculation conditions in Table 3, we carried out a total number of 20 numerical simulations, and obtained, in all, 60 sections of volumetric water content variation. The variations of the volumetric water content of different α values with k r = 0.01 and k r = 0.1 are shown in Figure 8 to illustrate the influence of anisotropy direction α on seepage characteristics, and different k r values with α = 50 • and α = 70 • are also displayed in Figure 9 to show the impacts of anisotropy ratio k r .
Water 2020, 12, x FOR PEER REVIEW 10 of 20 According to the calculation conditions in Table 3, we carried out a total number of 20 numerical simulations, and obtained, in all, 60 sections of volumetric water content variation. The variations of the volumetric water content of different α values with kr = 0.01 and kr = 0.1 are shown in Figure 8 to illustrate the influence of anisotropy direction α on seepage characteristics, and different kr values with α = 50° and α = 70°are also displayed in Figure 9 to show the impacts of anisotropy ratio kr.   0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16  In Figure 8, the initial state was the distribution of the volumetric water content of the rock slope before rainfall, and the remainder was the distribution of the volumetric water content after the rain had ceased. Rainfall increased the volumetric water content of the surface of the rock slope, and the volumetric water content of the surface of the rock slope was affected by the anisotropy angle. Regarding the rock slope body, as the anisotropy angle decreased, the volumetric water content of the surface of the slope gradually increased. This was because a smaller the inclination in the plane led to a smaller hydraulic conductivity in the vertical direction, and a larger hydraulic conductivity in the horizontal direction, which prevented the rainwater from infiltrating into the deep part of the slope body and increased the surface volumetric water content of the rock slope. Figure 9 shows that rainfall increased the volumetric water content of the slope, and that the volumetric water content of the slope was affected by anisotropy ratio k r . It is worth mentioning that k r had a great influence on the variation of water content in the shallow part of the bottom of the slope, as shown in Figure 9c,f. It can be seen in Figures 8 and 9 that the turning point of the volumetric water content appeared within the shallow part of the rock slope. The moisture condition of the region above this point began to change from unsaturated to saturated, which indicated that the moist front had reached that turning point after the rain had ceased. The turning point of volumetric water content also appeared within the deep part of the rock slope. The moisture condition of the region above this point began to change from saturated to unsaturated, indicating that the groundwater line rose to this turning point after the rain had ceased.

Analysis of the Changes in the Groundwater Levels
As detailed in Figure 10a,c, the longer the duration of the rainfall event, the higher the groundwater levels rose. The increases in the groundwater level near the free face of the rock slope were found to be greater than those far from the free face. It is worth noting that when the rock masses were isotropic, the groundwater at the bottom of the slope decreased, since the hydraulic conductivity of the horizontal direction reached the maximum, and the horizontal drainage capacity the bottom of the slope was stronger than the vertical infiltration capacity of the rainwater, and the groundwater levels decreased. It can be seen in Figure 10d,e that when k r was constant, it was found that with the increases of α, the differences in the heights of the groundwater levels near the free face of the rock were smaller. Then, when α was constant, with the anisotropy being increased, the differences in the heights of groundwater levels near the free face of the rock slope were larger. (e) different kr values with α = 70° after rain had ceased.

Analysis of the Maximum Water Content of the Surface and the Rising Height of the Groundwater
To more prominently examine the seepage characteristics of the rock slope under different kr and α conditions, the maximum water content on the surface (MWCS) and the rising heights of the groundwater (RHG) were determined [25]. Figure 11 shows the variations of the volumetric water content on the slope over time. It can be seen in the figure that as the rainfall continued, the surface volumetric water content in the slope changed and gradually become saturated. Additionally, the maximum was reached at the end of the rainfall event. Therefore, the MWCS was used to describe the maximum water content on the slope surface during the rainfall event. As the rainfall infiltrated the slope, the groundwater level increased and the saturated areas expanded. Therefore, the RHG was used to describe the maximum height of the groundwater level increase when rainfall had subsided. To more prominently examine the seepage characteristics of the rock slope under different k r and α conditions, the maximum water content on the surface (MWCS) and the rising heights of the groundwater (RHG) were determined [25]. Figure 11 shows the variations of the volumetric water content on the slope over time. It can be seen in the figure that as the rainfall continued, the surface volumetric water content in the slope changed and gradually become saturated. Additionally, the maximum was reached at the end of the rainfall event. Therefore, the MWCS was used to describe the maximum water content on the slope surface during the rainfall event. As the rainfall infiltrated the slope, the groundwater level increased and the saturated areas expanded. Therefore, the RHG was used to describe the maximum height of the groundwater level increase when rainfall had subsided. volumetric water content in the slope changed and gradually become saturated. Additionally, the maximum was reached at the end of the rainfall event. Therefore, the MWCS was used to describe the maximum water content on the slope surface during the rainfall event. As the rainfall infiltrated the slope, the groundwater level increased and the saturated areas expanded. Therefore, the RHG was used to describe the maximum height of the groundwater level increase when rainfall had subsided.  It can be seen in Figure 12a,c that MWCS increases gradually with the decrease of α. The reason for this is mentioned in Section 3.1.1. It is worth noting that with the change of k r , MWCS has different change rules at different positions of slope body. At the top of the slope, MWCS gradually decreases with the increase of α. At the bottom of the slope, MWCS gradually increases with the increase of α. In the middle of the slope, the variations in α have little connection with MWCS under different k r and α conditions. The variations of the MWCS at the top of the slope were found to range between 0.154 and 0.163 and the change rate was 5.8%. The variations of the MWCS at the middle of the slope ranged between 0.155 and 0.160, and the change rate was 3.2%. The variations of the MWCS at the middle of the slope ranged between 0.156-0.160, and the change rate was 2.6%, The maximum and minimum values of the MWCS both appeared at the top of the slope, and the change rate of the MWCS at the top of the slope was the largest, which indicated that the MWCS at the top was most greatly affected by the k r and α. The change rate of the MWCS at the bottom of the slope was the smallest, indicating that the MWCS at the bottom was the least affected by the k r and α.
It can be seen from Figure 12d,f that the RHG gradually increased with decreases in k r and increases in α, since vertical penetration was enhanced and horizontal discharge was weakened, resulting in a large increase in the height of the water level. It was observed that under different k r and α conditions, the variations of the RHG at the top of the slope ranged between 11.32 and 12.47 m and the change rate was 10.16%. The variations of the RHG in the middle of the slope ranged between 12.40 and 15.05 m and the change rate was 21.37%. Also, the variations of the RHG at the bottom of the slope ranged from 3.88 to 12.22 m and the change rate was 214.95%.
Therefore, it was found that the maximum value of the RHG appeared in the middle of the slope, and the minimum appeared at the bottom of the slope. The change rate of the RHG at the bottom of the slope was the largest, indicating that the RHG at the bottom was the most affected by the k r and α. The RHG change rate at the top of the slope was observed to be the smallest, indicating that the RHG at the top was the least affected by the k r and α. Water 2020, 12, x FOR PEER REVIEW 14 of 20

Effects of the Hydraulic Conduction Anisotropy on the Stability of the Rock Slope
In the present study, the stability of rock slopes under different sliding surface conditions was comprehensively analyzed. Four arc sliding surfaces were assumed, referred to in this study as sliding surfaces 284, 293, 313, and 112, as shown in Figure 13. These types of sliding surfaces mainly penetrated through the strongly weathered layered slate. For the four types of surfaces, it was observed that as the rainfall continued, FS gradually decreased, as shown in Figure 14. It was found that under various anisotropy conditions, the differences in FS increased over time. When kr = 1 (for example, when the rock layers were isotropic to seepage), the FS tended to be overestimated. When analyzing the stability of the layered rock slope, only the anisotropy of the shear strength was incomplete, and the anisotropy of seepage could not be ignored for the layered rock slope.

Effects of the Hydraulic Conduction Anisotropy on the Stability of the Rock Slope
In the present study, the stability of rock slopes under different sliding surface conditions was comprehensively analyzed. Four arc sliding surfaces were assumed, referred to in this study as sliding surfaces 284, 293, 313, and 112, as shown in Figure 13. These types of sliding surfaces mainly penetrated through the strongly weathered layered slate. For the four types of surfaces, it was observed that as the rainfall continued, FS gradually decreased, as shown in Figure 14. It was found that under various anisotropy conditions, the differences in FS increased over time. When k r = 1 (for example, when the rock layers were isotropic to seepage), the FS tended to be overestimated. When analyzing the stability of the layered rock slope, only the anisotropy of the shear strength was incomplete, and the anisotropy of seepage could not be ignored for the layered rock slope.

Effects of the Hydraulic Conduction Anisotropy on the Stability of the Rock Slope
In the present study, the stability of rock slopes under different sliding surface conditions was comprehensively analyzed. Four arc sliding surfaces were assumed, referred to in this study as sliding surfaces 284, 293, 313, and 112, as shown in Figure 13. These types of sliding surfaces mainly penetrated through the strongly weathered layered slate. For the four types of surfaces, it was observed that as the rainfall continued, FS gradually decreased, as shown in Figure 14. It was found that under various anisotropy conditions, the differences in FS increased over time. When kr = 1 (for example, when the rock layers were isotropic to seepage), the FS tended to be overestimated. When analyzing the stability of the layered rock slope, only the anisotropy of the shear strength was incomplete, and the anisotropy of seepage could not be ignored for the layered rock slope.   To complete a more in-depth examination of the stability of the different sliding surfaces of layered slopes under different kr and α conditions, the FS after the rainfall had ceased were discussed. As detailed in Figure 15, it was observed that as the kr decreased and the α increased, the FS of each sliding surface decreased. When kr = 0.01 and α = 70°, the FS was found to be minimum. The FS of sliding surface No. 284 was determined to range between 1.41 and 1.47, and the rate of change was 4.26%. The FS of sliding surface No. 112 ranged between 1.49 and 1.60, and the rate of change was 7.38%. The FS of sliding surface No. 293 was observed to range from 1.57 to 1.67, and the rate of change was 6.37%. Additionally, the No. 293 sliding surface displayed a factor of safety of 1.66 to 1.70, with a rate of change of 2.41%. It was found that the smallest factor of safety occurred in the No. 284 sliding surface.  To complete a more in-depth examination of the stability of the different sliding surfaces of layered slopes under different k r and α conditions, the FS after the rainfall had ceased were discussed. As detailed in Figure 15, it was observed that as the k r decreased and the α increased, the FS of each sliding surface decreased. When k r = 0.01 and α = 70 • , the FS was found to be minimum. The FS of sliding surface No. 284 was determined to range between 1.41 and 1.47, and the rate of change was 4.26%. The FS of sliding surface No. 112 ranged between 1.49 and 1.60, and the rate of change was 7.38%. The FS of sliding surface No. 293 was observed to range from 1.57 to 1.67, and the rate of change was 6.37%. Additionally, the No. 293 sliding surface displayed a factor of safety of 1.66 to 1.70, with a rate of change of 2.41%. It was found that the smallest factor of safety occurred in the No. 284 sliding surface.  Figure 16 shows some rock landslides in this area and their stereographic projections. It can be seen that most of the sliding surfaces of rock landslides are circular arcs in the strongly weathered layer, which verifies the assumptions regarding the arc sliding surface in Section 3.3. Although some stereographic projections show that the rocky slopes are stable, they still slip. It can be seen that they all present large dip angles. Figure 17 shows the results of the field survey of landslides in the Pulang region. The number of field landslides increases as the anisotropy angle increases, which verifies that as the dip angle of the bedding plane increases, the FS of rock slope decreases.  Figure 16 shows some rock landslides in this area and their stereographic projections. It can be seen that most of the sliding surfaces of rock landslides are circular arcs in the strongly weathered layer, which verifies the assumptions regarding the arc sliding surface in Section 3.3. Although some stereographic projections show that the rocky slopes are stable, they still slip. It can be seen that they all present large dip angles. Figure 17 shows the results of the field survey of landslides in the Pulang region. The number of field landslides increases as the anisotropy angle increases, which verifies that as the dip angle of the bedding plane increases, the FS of rock slope decreases.

Conclusions
A layered rock slope in the Pulang Region, China was numerical simulated, considering different anisotropy ratios and anisotropy angles of hydraulic conductivity. Anisotropy ratios are the ratios of vertical hydraulic conductivity to horizontal hydraulic conductivity; anisotropy angles are the dip angles of the bedding plane. The research results were as follows: (1) The residual volumetric water content of the VG model had little effect on the hydraulic conductivity. The closer the fitting parameter a of the VG model was to 10, and the closer the fitting parameter n of the VG model was to 1.5, the closer the pore water pressure of the rock slope was to the field monitoring data.
(2) The maximum initial matric suction was determined to be of major significance to the unsaturated seepage of the rock slope and the subsequent calculations. This study set the maximum initial matric suction to −75 kPa, which was determined to be consistent with the actual situation.
(3) The different anisotropy ratios and dip angles of the bedding plane were found to have major impacts on the seepage in the layered rock slope.
(4) The MWCS and RHG were determined to characterize the seepage characteristics of the rock slope. As the dip angles of the bedding plane decreased, the MWCS gradually increased. As the anisotropy ratios decreased and the dip angles increased, the RHG gradually increased.
(5) When the seepage of layered rock slope was considered isotropic, the FS tended to be overestimated. As the anisotropy ratios decreased and the dip angles of the bedding plane increased, the FS of each sliding surface was reduced. When the dip angles of the bedding plane of the rock slope are larger, more essential protections are needed.

Conclusions
A layered rock slope in the Pulang Region, China was numerical simulated, considering different anisotropy ratios and anisotropy angles of hydraulic conductivity. Anisotropy ratios are the ratios of vertical hydraulic conductivity to horizontal hydraulic conductivity; anisotropy angles are the dip angles of the bedding plane. The research results were as follows: (1) The residual volumetric water content of the VG model had little effect on the hydraulic conductivity. The closer the fitting parameter a of the VG model was to 10, and the closer the fitting parameter n of the VG model was to 1.5, the closer the pore water pressure of the rock slope was to the field monitoring data. (2) The maximum initial matric suction was determined to be of major significance to the unsaturated seepage of the rock slope and the subsequent calculations. This study set the maximum initial matric suction to −75 kPa, which was determined to be consistent with the actual situation. (3) The different anisotropy ratios and dip angles of the bedding plane were found to have major impacts on the seepage in the layered rock slope. (4) The MWCS and RHG were determined to characterize the seepage characteristics of the rock slope. As the dip angles of the bedding plane decreased, the MWCS gradually increased. As the anisotropy ratios decreased and the dip angles increased, the RHG gradually increased. (5) When the seepage of layered rock slope was considered isotropic, the FS tended to be overestimated.
As the anisotropy ratios decreased and the dip angles of the bedding plane increased, the FS of each sliding surface was reduced. When the dip angles of the bedding plane of the rock slope are larger, more essential protections are needed.