Analyzing the Effect of Soil Hydraulic Conductivity Anisotropy on Slope Stability Using a Coupled Hydromechanical Framework

Abstract

In studies on the effect of rainfall on slope stability, soil hydraulic conductivity is usually assumed to be isotropic to simplify the analysis. In the present study, a coupled hydromechanical framework based on transient seepage analysis and slope stability analysis is used to investigate the effects of hydraulic conductivity anisotropy on rainfall infiltration and slope safety at various slope locations (the top of the slope, the slope itself and the toe of the slope). The results show that when the vertical hydraulic conductivity (K_y) is constant, the horizontal hydraulic conductivity (K_x) increases (i.e., anisotropy increases). This occurs because rainfall tends to infiltrate into the interior of the slope, resulting in the soil on top of the slope and on the slope itself being easily influenced by rainfall, leading to soil instability. The change of rainfall infiltration at the slope itself is the most significant. When the anisotropic ratio K_r (=K_x/K_y) increased from 1 to 100, the depth of the wetting zones for loam, silt and clay slopes increased by 23.3%, 33.3% and 50%, respectively. However, increased K_r led to a slower infiltration rate in the vertical direction at the toe of the slope. Compared to the results for K_r = 1 and for K_r = 100, the thickness of the wetting zones at the toe of loam and silt slopes decreased by 23.3% and 30.0%, respectively. For the clay slope, K_r changes did not significantly affect the wetting zones because of poor permeability. The results of this study suggest that the effect of soil hydraulic conductivity anisotropy should be considered when estimating slope stability to better understand the effect of rainfall on slopes.

1. Introduction

Studies on landslides, which occur in many parts of the world, have found that slope stability is influenced by not only internal factors, such as topography, geological structure of the slope and groundwater conditions but also external factors, such as precipitation, earthquakes and human activity [1,2,3,4,5]. Rainwater infiltration is the most common cause of landslides [6,7]. Mountainous areas all over the world are subjected to heavy rainfall in the wet season [8,9]. During rainfall events, rainwater infiltrates into the pores of the soil, increasing the water content of the soil in the dry zone (above the groundwater table and near the surface), which becomes a wetting zone or transient saturated zone [10,11,12,13]. Water that infiltrates into the soil increases the shear stress and the pore water pressure of the soil. When the pore water pressure of the soil increases, the effective stress decreases, resulting in a decrease in the shear strength of the soil and a decrease in slope stability. Therefore, rainfall seepage is considered to be an important factor in slope failure [14,15,16,17,18,19,20].

Numerous studies have focused on the mechanism of rainfall-triggered landslides [21,22,23,24,25,26,27,28,29,30]. Results show that when rainfall occurs on slopes, the hydraulic properties of soil and rainfall patterns control the seepage behavior. The spatial variability of soil’s properties is also an important factor that can influence rainfall infiltration and slope stability [31,32]. The saturated hydraulic conductivity of soil, K_s, limits the rainfall infiltration rate and rainfall characteristics, such as intensity and duration, control the rainfall amount that can infiltrate into the soil [33,34]. However, most research has assumed that the hydraulic conductivity of soil is isotropic. In a natural environment, the hydraulic conductivity of soil is usually anisotropic; however, few studies have considered the effect of soil hydraulic conductivity anisotropy on slope stability. Studies that discuss hydraulic anisotropy usually define the anisotropic ratio as K_x/K_y to describe soil hydraulic conductivity anisotropy, where K_x and K_y are the hydraulic conductivity in the horizontal and vertical directions, respectively. The values of the anisotropic ratio for alluvium range from 2 to 10 and that for clay soil can be higher than 100 [35]. Dong and Hsu [36] used a numerical model to estimate the effect of hydraulic conductivity anisotropy on the stability of stratified, poorly cemented rock slopes. They showed that neglecting hydraulic conductivity anisotropy resulted in overestimation of the safety factor of the slope by about 40%. Mahmood et al. [37] discussed the effect of anisotropic conductivity on the pore water pressure and reliability index of an unsaturated embankment. When the horizontal conductivity is fixed and the vertical conductivity is reduced, the rainfall infiltration rate in the vertical direction decreases and the slope reliability index remains unchanged. These studies show the significance of hydraulic anisotropy on the mechanism of rainfall-induced landslides. However, previous studies have not clarified the effect of hydraulic conductivity anisotropy on rainfall infiltration and slope stability at various locations on a slope (the top of the slope, the slope itself and the toe of the slope). Therefore, a coupled hydromechanical framework based on transient seepage analysis and slope stability analysis is established in this study to investigate the effects of hydraulic conductivity anisotropy on rainfall infiltration and slope safety at various locations on a slope.

Reliable limit equilibrium methods such as the ordinary method [38], Bishops’ simplified method [39], Janbu’s simplified method [40] and the Morgenstern and Price method [41] are commonly used to investigate slope stability. Although these methods can be used to easily calculate the safety factor of a slope by presuming the slope’s failure surface, it is difficult to identify where the slope failure begin and the real failure surface. Thus, the local factor of safety method [42], which is a new stability analysis method based on the Mohr-Coulomb failure criterion, is used in the present study. The advantage of using the local factor of safety (LFS) method is that this method can compute the factor of safety of each soil element, this is helpful to compare the effect of rainfall infiltration on soil stability at various locations or depths of a slope. In this study, a coupled hydromechanical framework based on transient seepage analysis and slope stability analysis is utilized for analyzing the effect of soil hydraulic conductivity anisotropy on slope stability.

2. Methodology

2.1. Seepage Analysis

Transient seepage analysis was conducted using the HYDRUS 2D [43] numerical model. The solution for transient unsaturated flow is based on Richards equation [44]:

$$\frac{\partial \theta \left(h\right)}{\partial t}=\nabla \cdot K\left(h\right)\nabla H+W$$

(1)

where $\theta$ [-] is the volumetric water content, t [T] is the time, h [L] is the pore water pressure or suction head, H [L] is the total head, W [L³T⁻¹] is the boundary flux, K(h) is the hydraulic conductivity function (HCF) and $\theta \left(h\right)$ is the soil-water retention curve (SWRC). In this study, the van Genuchten [45] SWRC model—shown in Equation (2)—is used to describe the relationship between soil volumetric water content and soil suction and the Mualem model [46]—shown in Equation (3)—is used for the HCF:

$$\theta \left(h\right)={\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left[1+{\left(\alpha \left|h\right|\right)}^n\right]}^{1-1/n}}$$

(2)

$$K=K_s\cdot {\left(S_e\right)}^l{\left\{1-{\left[1-{\left(S_e\right)}^{1/m}\right]}^m\right\}}^2$$

(3)

where ${\theta }_s$ [-] represents the saturated water content, ${\theta }_r$ [-] represents the residual water content, $\alpha$ [L⁻¹], n [-] and l [-] are the fitting parameters related to the inverse of the air entry value (suction), pore size distribution of soil and pore connectivity of soil, respectively and m = 1 − 1/n. $K_s$ [LT⁻¹] represents the saturated hydraulic conductivity and $S_e$ [-] is the equivalent degree of saturation, defined as:

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(4)

2.2. Theory of Effective Stress and Local Factor of Safety

In the stability analysis, the concept of suction stress [47] was used to calculate the soil effective stress. Suction stress includes all the inter-particle stresses such as capillary stress, the electric double-layer force, the van der Waals attractive force and the matric suction of soil. Lu et al. [48] postulated that suction stress represents the energy stored by the soil water and provided the closed-form equations for the effective stress and suction stress in soils, respectively, expressed as:

$${\sigma }^{\prime }=\sigma -u_a-{\sigma }^s$$

(5)

$${\sigma }^s=-\left(u_a-u_w\right)$$

(6)

$${\sigma }^s=-\frac{\left(u_a-u_w\right)}{{\left\{1+{\left[\alpha \left(u_a-u_w\right)\right]}^n\right\}}^{\left(n-1\right)/n}}$$

(7)

where ${\sigma }^s$ is the suction stress [ML⁻¹T⁻²], $u_a$ is the pore air pressure [ML⁻¹T⁻²], $u_w$ is the pore water pressure [ML⁻¹T⁻²], $\left(u_a-u_w\right)$ is the matric suction of soil [ML⁻¹T⁻²] and $\alpha$ [L⁻¹] and n [-] are fitting parameters identical to those in the van Genuchten [45] SWRC model.

The local factor of safety method [42] was used in this study to estimate slope stability. This method is based on the Mohr-Coulomb failure criterion (Figure 1). It defines the local factor of safety (LFS) as the ratio of the potential Coulomb stress to the current Coulomb stress, shown in Equation (7). In this theory, the solid Mohr circle is assumed to be the current state of stress of a soil element. When the water content of the soil is increased, the solid Mohr circle shifts leftward because of the suction stress increasing. The size of the Mohr circle remains nearly unchanged because it is related to the difference in the principal total stresses $\left({\sigma }_1-{\sigma }_3\right)$, which are determined by the slope geometry and the self-weight of the soil. If the Mohr circle shifting leftward, the Mohr circle is closer to the Mohr-Coulomb failure envelope, then, the state of stress of the soil is closer to the failure condition.

$$LFS=\frac{\tau *}{\tau }=\frac{BE}{CD}=\frac{\cos \varphi }{{\sigma }_1^{\prime }-{\sigma }_3^{\prime }}\left(2c+\left({\sigma }_1^{\prime }+{\sigma }_3^{\prime }\right)\tan \varphi \right)$$

(8)

In Equation (8), $\tau *$ is the shear strength or potential Coulomb stress [ML⁻¹T⁻²], $\tau$ is the shear stress or current Coulomb stress [ML⁻¹T⁻²], c is the effective cohesion of the slope material [ML⁻¹T⁻²], ${\sigma }_1^{\prime }$ and ${\sigma }_3^{\prime }$ are the maximum effective stress and minimum effective stress [ML⁻¹T⁻²], respectively and $\varphi$ is the angle of shearing resistance of the slope material [LL⁻¹].

In addition, referring to the effective stress theory of Lu and Likos [47], substituting Equation (5) into Equation (8) yields:

$$LFS=\frac{\tau *}{\tau }=\frac{\cos \varphi }{{\sigma }_1-{\sigma }_3}\left(2c+\left({\sigma }_1+{\sigma }_3-2{\sigma }_s\right)\tan \varphi \right)$$

(9)

Then, the solution of LFS can be utilized for calculating the stability of the slope and is combined with the finite element method to estimate the safety factor variation with space and time. If the suction stress of the soil is increased, indicating that the LFS of the soil is reduced, which means the stability of the soil is decreased.

2.3. Coupled Hydromechanical Framework

The finite element software HYDRUS 2D [43] and the Slope Cube Module [49] were used to establish a coupled hydromechanical framework to investigate the effect of rainfall infiltration on slope stability. The geometry and boundary conditions of the two-dimensional (2D) slope model are shown in Figure 2. The boundaries of the slope surface (BC, CD and DE boundaries) were set as the rainfall infiltration boundaries and the sides of the slope (AG and FH boundaries) were set as the boundaries of the fixed water level to maintain the depth of the original groundwater level in the slope. BG, EH and AF were the no-flow boundaries; the flux at these three boundaries is zero.

In this study, three types of soil were considered, namely loam, silt and clay, representing high to low hydraulic conductivity residual soils. The hydraulic properties of the materials include saturated water content (${\theta }_s$), saturated hydraulic conductivity ($K_s$) and the fitting parameters $\alpha$ and n for SWRC [44], as shown in

Table 1. Hydraulic properties of slope material [50].

Material	${\mathit{\theta }}_{\mathit{s}}$	$\mathit{\alpha }$ (kPa−1)	n	Ks (m/s)
Loam	0.43	0.36	1.56	$1.0\times {10}^{-5}$
Silt	0.46	0.16	1.37	$5.0\times {10}^{-6}$
Clay	0.38	0.08	1.09	$1.0\times {10}^{-7}$

. From the theories of SWRC, HCF and suction stress mentioned above, the corresponding pore water pressure, soil water permeability and soil suction stress of the slope materials can be obtained for various water content conditions (Figure 3).

Table 2. Mechanical properties of slope material [51].

Material	Gs (-)	c (kPa)	$\mathit{\varphi }$ (°)	E (kPa)	ν
Loam	2.65	10	35	15,000	0.3
Silt	2.70	15	30	10,000	0.35
Clay	2.80	20	25	7000	0.45

shows the compiled mechanical characteristic parameters of the soils, where G_s represents the specific weight of the soil, E represents the Young’s modulus of the soil and v is Poisson’s ratio.

In the coupled hydromechanical framework, the solution for transient unsaturated flow [43] was used to conduct seepage analysis of 2D slopes to examine the effects of rainfall infiltration on the water content, soil suction and soil unit weight in unsaturated slope layers. Then, the solutions of the seepage analysis were input into FEM2D [52] which is one part of code wrote in the Slope Cube Module and is used to compute the soil stresses and displacements by using the field of moist unit weight as the gravity term [49]. A linear momentum balance equation is used to describe the total stress changes in slope materials:

$$\nabla \cdot \left(\sigma \right)+\gamma \left(\theta \right)b=0$$

(10)

where $\sigma$ represents the stress tensor in 2D space, $\gamma \left(\theta \right)$ represents the unit weight of the slope material related to water content and b is the unit vector of body forces with two components. The field of soil suction stress and effective stress can be estimated using Equations (5) and (7), respectively, with the solutions of the field of soil saturation, soil suction and total stress. The LFS method can then be used to calculate the soil stability field.

3. Results and Discussion

In a nature environment, soil always sediment by gravity and shows layered distribution, resulting in the hydraulic conductivity of soil in all directions are different (i.e., hydraulic conductivity anisotropy). According to previous studies, the conductivity in the horizontal direction (K_x) is usually higher than that in the vertical direction (K_y) and the anisotropic ratio (K_r = K_x/K_y) ranges from 1 to 100 [35]. The K_y value of slope materials (loam, silt and clay) was defined as constant in this research and the K_x value was varied to design different anisotropic cases for the analysis. The settings of each case are shown in

Table 3. Case design for hydraulic conductivity anisotropy.

Rainfall Intensity (mm/h)	Ky (m/s)	Kr
14.5	$\mathrm{loam}:1.0\times {10}^{-5}$	$\left[\begin{array}{c}1\\ 10\\ 50\\ 100\end{array}\right]$
	$\mathrm{silt}:5.0\times {10}^{-6}$	$\left[\begin{array}{c}1\\ 10\\ 50\\ 100\end{array}\right]$
	$\mathrm{clay}:1.0\times {10}^{-7}$	$\left[\begin{array}{c}1\\ 10\\ 50\\ 100\end{array}\right]$

. The rainfall intensity was set at 14.5 mm/h, which is the threshold value for long-duration extreme events, as defined by the National Science and Technology Center for Disaster Reduction. The rainfall duration is longer than 1 day and the rainfall amount in 1 day is about 350 mm [53]. In this study, the analysis is performed for 48 h (rainfall duration is 48 h) and was average divided into 48 time steps.

The analysis results show a significant difference in the rainfall effect on the slope saturation and stability at different locations (Figure 4). When K_r is equal to 100 (high anisotropy case), the thickness of the wetting zone caused by rainfall increases for the slope. However, at the toe of the slope, with increasing anisotropic ratio, the thickness of the wetting zone decreases for a given rainfall intensity. Thus, the results of the X-X,’ Y-Y’ and Z-Z’ sections (Figure 2), representing the top of the slope, the slope itself and the toe of the slope, respectively, are discussed here to clarify the effect of hydraulic conductivity anisotropy.

3.1. Top of Slope (X-X’ Section)

The results of water content, soil suction stress and LFS at the top of the slope (X-X’ section) after a rainfall event in different cases are first discussed. The results for the loam slope show that the thickness of the wetting zone caused by rainfall infiltration increases by about 12.5% when K_r is increased from 1 to 100 (Figure 5). When K_r is increased and K_y remains constant, the soil permeability in the horizontal direction increases. The rainfall in the soil tends to infiltrate in the horizontal direction and thus the rainwater permeates deep into the slope. The soil at the top of the slope is also affected by this phenomenon. The soil water content at depths in the range of 2.8 to 3.2 m is increased in the high anisotropy case and the change of water content also affect soil suction stress and LFS, as shown in Figure 5b,c.

For silt and clay slopes, although high K_r increases soil permeability in the horizontal direction, the original hydraulic conductivity of these two slope materials is so small that the results of modeling do not show the effect of anisotropy changes on the saturation and stability at the top of the slope (X-X’ section).

3.2. Slope (Y-Y’ Section)

The results of the Y-Y’ section are used to discuss the effect of anisotropy changes on the saturation and stability of the slope. The results of water content field, suction stress field and LFS field are shown in Figure 6. In each case, the results show that the thickness of the wetting zone caused by precipitation increases with increasing K_r. This is because K_x increases when K_r is increased, rainfall infiltration flow tends to be in the horizontal direction and seepage can easily affect soil saturation and stability. Compared with the K_r = 1 results, the depth of the wetting zone increased by about 23.3%, 33.3% and 50% for loam, silt and clay slopes, respectively, for K_r = 100. It also shows that the LFS tends to easily drop down when K_r is high.

3.3. Toe of Slope (Z-Z’ Section)

For K_r = 100, the soil below depths of 2.8, 3.0 and 5.0 m for the loam, silt and clay slopes, respectively, is affected by groundwater table rises, as shown in Figure 7. The groundwater table was originally inclined by 5°; however, when the horizontal hydraulic conductivity increases, the initial groundwater table is not under a steady state, it tends to become level with time, resulting in a groundwater table rise at the toe of the slope. This study also conducts these analyses without rainfall, the results are the same, show that the groundwater table rising is because of the horizontal hydraulic conductivity increases.

As mentioned above, when K_r increases, rainfall tends to infiltrate in the horizontal direction. Thus, the rainwater can flow into the slope body quickly in the horizontal direction in high K_r cases, which decreases the vertical infiltration rate. As shown in the results for the loam and silt slopes, when K_r = 100, the thickness of wetting zone at the toe of the slope (Z-Z’ section) decreases by about 23.3% and 30.0%, respectively, compared with that for K_r = 1. And the changes of the wetting front also affect soil stability, the results show that the safety factor of the soil at 2 m below the surface is higher when K_r increases. The hydraulic conductivity of clay is so small that the anisotropy changes do not affect the soil saturation and stability of the clay slope.

The analysis results show that soil hydraulic conductivity anisotropy can influence the behavior of rainfall infiltration at different locations of the slope, especially when K_r is higher than 10. However, as mentioned above, the values of the anisotropic ratio for alluvium usually range from 2 to 10 and that for clay soil can be much higher. Therefore, slope stability analyses in the further research could consider hydraulic conductivity anisotropy to better understand the impact of rainfall on the slope when the anisotropic ratio is higher than 10. In addition, this study shows the advantages of the LFS method, which allows the effect of rainfall on soil stability to be evaluated at each point within a slope. However, the study considering a slope made up of homogeneous material is a strong limitation that limits the application to real cases of study. In a real slope, the hydrology of a hillslope is usually influenced by the internal discontinuities (layering, bedding, pedology, etc.). Thus, in the future research, conducting geological surveys to determine slope geometry, structure and using probabilistic method to determine soil properties in the numerical model could make the results of the analyses more reliable.

4. Conclusions

A coupled hydromechanical framework was used in this study to investigate the effects of hydraulic conductivity anisotropy on rainfall infiltration and slope safety. Changes in the water content field, suction stress field and LFS field show that the anisotropy of soil permeability has an effect on rainfall infiltration. When the vertical hydraulic conductivity is fixed, an increase in K_r represents an increase in soil seepage mobility in the horizontal direction, making the soil on the slope itself more susceptible to rainfall infiltration, decreasing the soil safety factor. When K_r was increased from 1 to 100, the depth of the wetting zone on the slope itself increased by 23.3%, 33.3% and 50.0% for loam, silt and clay slopes, respectively and resulting in the soil LFS more easily to decrease. Although the results of loam slope also show the effect of anisotropy changes on the saturation and stability at the top of the slope (X-X’ section), the original hydraulic conductivity of silt and clay slope is so small that the results of modeling do not show the effect of anisotropy changes. In addition, when the horizontal hydraulic conductivity increased, the rainfall at the toe of the slope tended to flow into the slope body in the horizontal direction, decreasing the vertical infiltration rate. Compared to the results for K_r = 1, for K_r = 100, the thickness of the wetting zone at the toe of the loam and silt slopes decreased by 23.3% and 30.0%, respectively. Thus, it can be seen that rainfall infiltration behavior was most affected by hydraulic conductivity anisotropy at the slope (Y-Y’ section). The rainfall infiltration behavior at the top of the slope and the toe of the slope was less sensitive to the anisotropy of hydraulic conductivity for slope materials with poor permeability, such as clay. Therefore, hydraulic conductivity anisotropy could be considered in analyses of slope stability to better predict the effect of rainfall infiltration on slope safety, especially when K_r is higher than 10.

However, the study considering a slope made up of homogeneous material is a strong limitation that limits the application to real cases of study. Thus, future research could conduct geological surveys to determine the slope geometry, structure and soil properties in the numerical model to improve the reliability of the study results.