#### Abstract

Traditional slope stability analysis mostly adopts the limit equilibrium method, which predetermines the slope failure surface and assumes that failure occurs simultaneously at all points of the failure surface. The method is based on the balance of forces and torques. The slope stability is represented by the factor of safety. The lowest factor of safety obtained after repeated analysis indicates the most failure-prone slope surface. However, the factor of safety for only one slope failure surface is obtained when applying this method. The distribution and changes of factor of safety in the interior of the slope are not identified. In addition, the analysis of factor of safety is influenced by the uncertainty in soil mechanical parameters, whereas uncertainty is not quantified in the traditional deterministic analysis. Therefore, a probabilistic approach, which uses the probability distribution function to explain the randomness of parameters, is proposed for quantifying the uncertainty. Nonetheless, when the observation data are not sufficient for determining the probability distribution function, the fuzzy theory can be an alternative method for the analysis. The fuzzy theory is based on fuzzy sets. It expresses the ambiguity of incomplete sets of information using a membership function. Moreover, a correct judgment can be made without verbose iterations. Hence, the aim of this study is to examine the uncertainty in soil mechanical parameters. The membership functions between soil mechanical parameters, i.e., cohesion and angle of internal friction, were constructed based on the fuzzy theory. The fuzzy point estimation was used in combination with the hydrologic and mechanical coupling model on HYDRUS 2D and the Slope Cube Module. The local factor of safety at different depths of the slope was determined using the local factor of safety theory. The probability of failure at different depths was calculated through reliability analysis, which could serve as an early warning for subsequent slope failures.

#### 1. Introduction

Slope stability is affected by intrinsic and triggering factors. The intrinsic factors include soil, groundwater, vegetation, slope gradient, and lithology. The triggering factors include volcanic eruptions, earthquakes, and rainfall. A common trigger for natural slopes is rainfall [1–8]. Rainfall-induced slope failures are usually shallow, with a depth of failure not exceeding three meters, and they likely occur on slopes with a gradient of 30° to 40° [9]. Lu and Godt [10] suggested that the failure mechanism for rainfall-induced shallow failures is that, as the rainfall infiltrates into the soil, the soil matric suction declines and the pore pressure rises positively. As the soil matric suction decreases, there would be a nonlinear drop in soil shear strength. Hence, when the soil is nearly saturated, the matric suction approaches zero, resulting in slope instability and further inducing disasters such as landslides and debris flow.

Studies related to rainfall-induced slope failure can be divided into three types according to their theoretical basis: statistical-model-based [4, 11–19], contributing factor [20–24], and physical-model-based analyses [5, 25–31]. Among them, the physical-model-based analysis coupled with hydromechanical mechanism models has overcome the excessive dependence of statistical models on rainfall data. The method can describe the hydromechanical changes caused by transient rainfall in the interior of the slope, as well as the associated failure mechanism. With its higher predictive power and capability of quantifying the effect of each parameter on slope stability [32], the method is now widely used. Nevertheless, the analytical process is limited by uncertainty caused by measurement error, spatial variability, and limited information [33]. The result of slope stability analysis may deviate from reality owing to the uncertainty in model parameters [34].

Therefore, probabilistic analysis is used to quantify the uncertainty [7, 35–39]. Nawari and Liang [40] and Giasi et al. [41] suggested that an adequate number of reliable observation values are required for probabilistic analysis. Precise mean values and standard deviations are derived from the observation values to construct a reasonable probability density function [42]. In addition, Juang (in 1998) and Nawari and Liang [40] proposed that the uncertainty in parameters may be nonstochastic. Previous studies have shown that, when the data available are not sufficient for defining the probability density function, the uncertainty in rock mass parameters can be expressed effectively with the use of a fuzzy set [43, 44]. This method has been applied to some of the cases for slope stability analysis [30, 41, 45–48].

Traditional slope stability analysis adopts the limit equilibrium analysis, which discretizes the potential sliding soil mass into smaller vertical slices without considering soil deformation. It assumes that failure occurs simultaneously at all points of the failure surface. This method is based on the balance of forces and torques. The slope stability is represented by the factor of safety. Various analytical methods have been developed based on different assumptions on the balance of forces [49–52]. In recent years, the finite element method has been widely applied to slope stability analysis in order to calculate the factor of safety in slopes with high complexity (complex geometries, boundaries, and loading conditions) and to investigate the stress–strain relationship in soil [53–58]. Liu and Shao [59] introduced the finite element limit equilibrium analysis, which combines the limit equilibrium analysis and finite element analysis. It is used to examine the slope stability and evaluate the breaking load of a rigid foundation and retaining wall.

The above analytical methods based on the balance of forces or on the stress field usually seek a single general slope stability index. Hence, it is almost impossible to identify the changes in pore water pressure and effective stress owing to rainfall infiltration, or the actual slope failure surface and its geometry. Therefore, Lu et al. [60] proposed the theory of local factor of safety (LFS), which can calculate the factor of safety at discrete points in the soil mass and describe the geometry and position of the potential failure surface. Previous studies have revealed that the factor of safety (probability of failure) is highly dependent on the coefficient of correlation between cohesion and angle of internal friction [61–63]. It has been shown that the two parameters are not independent of each other and that the correlation between them is mostly negative [64–67]. Jiang et al. [63] noted that, when analyzing the probability of failure, a significant deviation may occur if we assume an independent relationship between cohesion and angle of internal friction (i.e., no correlation). Aladejare and Wang [68] also pointed out that neglecting the coefficient of correlation between cohesion and angle of internal friction may result in an order-of-magnitude difference in the result of the analysis. Moreover, the factor of safety does not necessarily reflect the actual safety level. With the use of reliability analysis, considering the variability of variables and calculating the probability of failure and reliability index will provide a more valid representation of the reliability of slope stability.

Hence, the aim of this study is to examine the uncertainty in soil mechanical parameters. The membership functions for the soil mechanical parameters, i.e., cohesion and angle of internal friction, were constructed based on the fuzzy theory. The fuzzy point estimation was used in combination with the hydromechanical coupling model on HYDRUS 2D and the Slope Cube Module. The local factor of safety at different depths of the slope was determined. The probability of failure at different depths was calculated through reliability analysis, which could serve as an early warning for subsequent slope failures.

#### 2. Materials and Methods

##### 2.1. Seepage Analysis

In this study, the analytic solution of transient seepage in an unsaturated layer developed by Šimůnek et al. [69] based on the Richards equation was used as the governing equation of the two-dimensional seepage as follows: where is the volumetric water content (-), is the time (), is the pore water pressure or hydraulic head (), is the total head (), is the source or sink (), is the hydraulic conductivity function (HCF) that varies with the pore water pressure (), and is the volumetric water content that varies with the pore pressure in the soil-water retention curve (SWRC) ().

The soil water content and HCF of an unsaturated zone vary with the hydraulic head and are highly nonlinear. In this study, the relationship between soil water content and matric suction was predicted using the closed-form analytic solution proposed by van Genuchten [70] (see equation (2)). It is also referred to as the SWRC. Based on the SWRC, Mualem [71] introduced the HCF for unsaturated layers (see equation (3)). where is the saturated soil water content (), is the residual soil water content (), is the matric suction (), is the reciprocal correlation of the air-entry value (), is related to the SWRC gradient (), is the hydraulic conductivity in saturated soil (), , is the coefficient of correlation of soil porosity (), and is the equivalent degree of saturation (), shown as

##### 2.2. Principle of Effective Stress in Unsaturated Soil

We adopted the principle of effective stress proposed by Lu and Likos [72], which unified the possible physical and chemical interparticle mechanisms in soil and proposed the concept of suction stress. The effective stress based on the concept of suction stress is shown as follows [73]: where is the suction stress (), is the Born repulsive force (), is the capillary force (), is the combined van der Waals attractive force and electric double-layer force (), is the degree of saturation in the soil (), and is also the matric suction (). The matric suction, capillary force, van der Waals attractive force, and electric double-layer force balance the Born repulsive force in the soil. However, as the grain size of the soil increases, the effect of the van der Waals attractive force and electric double-layer force becomes negligible.

As each of the stress components in soil can be expressed as a function of matric suction , degree of saturation , and water content , and as the suction stress in soil is mainly controlled by the soil water content, Lu et al. [74] derived the suction stress characteristic curve (SSCC) from the soil-water characteristic curve, based on the principle of thermodynamics and by considering suction stress as the energy stored in the pedon. The following analytical solution is shown: where is the equivalent degree of saturation (-), is the residual saturation (-), is the soil water content (-), is the saturated soil water content (-), and is the residual soil water content (-). Moreover, van Genuchten [70] calculated the equivalent degree of saturation using the following closed-form equation: where and are fitting parameters correlated to the air-entry value of SWRC and the gradient, respectively. Therefore, the suction stress can be expressed in the following forms. The change in soil suction stress with water content can be illustrated by estimating the SSCC:

##### 2.3. Theory of Local Factor of Safety

The local factor of safety is based on the Mohr–Coulomb failure criterion, and is defined by the ratio between the potential Coulomb stress and the current Coulomb stress as follows: where is the potential Coulomb stress and is the current Coulomb stress. The theory is illustrated in Figure 1, in which the current state of stress in the soil is represented by the realization of Mohr’s circle. The shear stress acting on the soil when a failure occurs is obtained by translating Mohr’s circle to the Mohr–Coulomb failure envelope. When the effective stress of the soil decreases owing to the increase in water content, Mohr’s circle is translated leftward, during which its size is almost unchanged. By extending the Coulomb stress, the potential Coulomb stress at the intersection point of Mohr’s circle and the Mohr–Coulomb failure envelope (point B) is determined. The local factor of safety is obtained by the calculation of similar triangles as follows: where is the effective cohesion of the soil, is the effective angle of friction of the soil, and and are the maximum and minimum effective stresses of the soil, respectively.

The following expression of LFS can be derived from equations (11) and (12):

Substituting equation (5) into equation (13) gives

Using modeling and finite element analysis, we can analyze the effect of changes in water content or suction stress on the stability of soil units at different locations or depths of the slope.

##### 2.4. Fuzzy Theory

The fuzzy theory is also called the fuzzy set theory. The fuzzy number is a special case in a fuzzy set. If no assumption is specified (when limited data are available), the fuzzy number is assumed to be triangular and comprises maximum (), minimum (), and modal or peak () values. The maximum and minimum values of a fuzzy number can be expressed as

The -value is determined by the actual engineering situation of the slope and ranges between 0.5 and 3 [45]. The larger the -value, the larger is the scope of distribution of the mechanical parameter, and hence, the lower is the reliability of the selected parameter, and vice versa. In this study, the -value was considered to be 2.

Fuzzy point estimation combines the vertex method and the point estimate method. The vertex method was proposed by Dong and Shah [75]. The method is based on -cut and interval analysis. It computes combinations of vertices of the variables, which replace the membership functions as an input variable. Therefore, given membership functions of the input variables, there would be combinations of vertices. The point estimate method developed by Rosenblueth [76] evaluates the uncertainty parameters of a performance function. The mean value and standard deviation of a performance function are assessed using a two-point estimate. The upper limits of variables obtained from the -cut sets are and . Four sets of vertex combinations are derived through modeling to yield four sets of output values . The contributions of each -cut set value to the result are compared. In this study, we adopted the concept of fuzzy weighted average. The mean value and standard deviation of the factor of safety are illustrated below: where is the number of -cut sets. Nine -cut sets ranging from 0.1 to 0.9 were considered in this study.

The probability of failure was calculated from the reliability index [77] assuming a normally distributed factor of safety. Therefore, the reliability index is normally distributed. The reliability index and the probability of failure () can be represented as follows:

#### 3. Results and Discussion

In this study, a two-dimensional numerical model was developed using HYDRUS 2D. We performed a transient seepage analysis based on the seepage theory proposed by Richards (1931). The Slope Cube Module was used to examine the stress change experienced by the soil. Slope stability analysis was performed using the local factor of safety theory. The probability of failure at different depths of the slope was calculated through reliability analysis. The slope is 18 m high, with a slope angle of 40°. Figure 2 illustrates the conceptual model of the slope. BCDE is the boundary of rainfall infiltration, AG and HF define the hydraulic head boundaries, and BG, EH, and AF are the zero-flow boundaries. Observation surfaces were set at the top, middle, and toe of the slope, whereas observation points were installed at the middle part of the slope. As shown in Figure 3, the grid consists of 5,661 nodes and 11,524 elements. We simulated the rainfall intensity with reference to the data from the Alishan Weather Station where the greatest rainfall was recorded during the 2009 Typhoon Morakot. The recorded 48 h cumulative rainfall was 2,361 mm. The simulation duration was set to 48 h and the rainfall intensity was set to 49.18 mm/h.

We have considered loam and silt as examples in this study. The soil hydraulic properties and mechanical parameters are listed in Tables 1 and 2. Using the empirical formula developed by van Genuchten [70], we estimated the SWRC and the HCF. The hydraulic characteristics of loam and silt are shown in Figure 4.

**(a)**

**(b)**

**(c)**

The variables in this study include cohesion and the angle of internal friction. The values from Table 2 were considered as the mean. The equations for triangular fuzzy numbers are as follows:

The degree of variation of parameters is described by the coefficient of variation (). The larger the cov, the greater is the degree of variation. The cov value of cohesion is approximately 25–30%, and that of the angle of internal friction is approximately 10–20% [33, 78, 79]. We have selected the maximum values for cov, which are 30% for cohesion and 20% for the angle of internal friction, to construct the triangular fuzzy numbers. The triangular fuzzy numbers for loam and silt are presented in Figures 5 and 6, respectively. The -cut values are listed in Table 3.

**(a)**

**(b)**

**(a)**

**(b)**

##### 3.1. Comparison of Probability of Failure for Different Types of Soil

In this study, the correlation between cohesion and angle of internal friction () was not considered. With the use of -cut sets, combinations with different degrees of membership were computed for modeling. We calculated the factor of safety and reliability index of loam and silt using data from the observation points at the middle of the slope under the same rainfall condition. Their relationship with the degree of membership is shown in Figure 7. The factor of safety of loam fluctuated within (1.2435, 1.3135) whereas that of silt fluctuated within (1.5588, 1.6225). For loam, the fuzzy reliability is determined to be 1.3357, and the probability of failure is 0.0908. For silt, the fuzzy reliability is determined to be 2.2299, and the probability of failure is 0.0129.

As the coefficient of permeability for loam was greater than that for silt in this study, the rainfall was likely to infiltrate into the interior of the slope, increasing the suction stress while decreasing the effective stress on the interior of the slope. Consequently, after 48 h of sustained rainfall, the factor of safety of loam was lower than that of silt at the observation points. The reliability index analysis reveals that, as the degree of membership increases, the reliability index increases. The results obtained from the observation points on the slope indicate that the probability of failure of a loam slope is 7.79% higher than that of a silt slope.

##### 3.2. Comparison of Probability of Failure at Different Times

We investigated the change in suction stress owing to the change in soil water content in the slope at different times, as well as the change in the probability of failure after a sustained infiltration of rainfall into the interior of the soil. Observations were obtained at the 12th, 24th, and 48th hours. The variations in water content, suction stress, and the probability of failure of loam and silt at the top, middle, and toe of the slope are presented as follows:

In the loam soil slope, under the effect of sustained rainfall infiltration, the rainfall intensity exceeded the coefficient of permeability. Consequently, at the 12th hour, the surface layer of loam approached saturation with a water content of 0.43. As shown in Figure 8, the thickness of the soil moisture belt at the 48th hour was 2.74 m at the top, 2.95 m at the middle, and 2.29 m at the toe of the slope. As the soil water content increased, the soil suction stress increased. At the 48th hour, the increase in suction stress was 18.40 kPa at the top, 14.08 kPa at the middle, and 10.88 kPa at the toe of the slope. Analysis of probability of failure reveals that there was a low probability of failure at the top of the slope. At the middle and toe of the slope, the greatest change in the probability of failure was observed to be on the soil surface. Such changes decreased with depth. The probability of failure of the loam slope varied with time. At the middle of the slope, it increased by 60.56% at the 12th hour, by 65.13% at the 24th hour, and by 67.34% at the 48th hour. At the toe of the slope, it increased by 20.97% at the 12th hour, by 21.93% at the 24th hour, and by 22.09% at the 48th hour.

**(a)**

**(b)**

**(c)**

In the silt soil slope, under the effect of sustained rainfall infiltration, the rainfall intensity exceeded the coefficient of permeability. Therefore, the surface layer of silt approached saturation, with a water content of 0.46 at the 12th hour. Figure 9 illustrates that the thickness of the soil moisture belt was 1.54 m at the top, 2.44 m at the middle, and 1.67 m at the toe of the slope. Changes in suction stress were determined to be 55.53 kPa at the top, 38.00 kPa at the middle, and 26.46 kPa at the toe of the slope. Analysis of probability of failure showed that there was a low probability of failure at the top of the slope. At the middle and toe of the slope, the greatest change in the probability of failure was observed to be on the soil surface. Such changes decreased with depth. The probability of failure of the silt slope varied with time. It increased by 32.86% at the 12th hour, by 42.07% at the 24th hour, and by 47.72% at the 48th hour. At the toe of the slope, it increased by 56.51% at the 12th hour, by 60.00% at the 24th hour, and by 61.63% at the 48th hour.

**(a)**

**(b)**

**(c)**

We observed that a greater change in probability of failure is associated with the infiltration depth and variation in suction stress. The variation of suction stress on the surface layer of silt was greater than that of loam. Nevertheless, the coefficient of permeability was lower for silt, limiting the rainfall infiltration depth. Consequently, under the same rainfall condition, the depth of the moisture band in silt was shallower than that of loam, as shown in Figures 10(a) and 10(c). Therefore, surface runoff owing to rainfall is likely to be formed in silt. As the suction stress of a slope is influenced by the depth of the moisture band, the effect of suction stress extended deeper in loam (Figure 10(b)) than in silt (Figure 10(d)). As the duration of rainfall increased, the probability of failure of the slope increased from the toe toward the middle of the slope. Figure 11 shows that the middle of the slope was affected by rainfall infiltration and recharge from the top of the slope simultaneously. Therefore, the probability of failure increased with time. Overall, after 48 h of rainfall, the area of the loam slope in which the probability of failure was greater than 50% was approximately twice as large as that of the silt slope.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 3.3. Effect of Parameter Correlation on the Probability of Failure

We assessed the effect of the coefficient of correlation of the parameters on the probability of failure. The top of the slope was not investigated because of its low probability of failure. Only the middle and toe of the slope were discussed. Previous studies have reported a correlation between cohesion and the angle of internal friction and that the correlation is mostly negative [64–66, 80]. In this study, the analysis of loam and silt slope stability reveals that (at the 48th hour) the stronger the negative correlation between the mechanical parameters the lower is the probability of failure. Observation points situated one meter below the surface at the middle and toe of the slope suggest a linear relationship between the probability of failure and correlation coefficient. Figure 12 illustrates the differences in the computed probability of failure when considering a negative coefficient of −0.8 (compared to that without considering the correlation). In the loam slope, the probability of failure decreases by 0.9% at the middle and by 0.7% at the toe of the slope. In the silt slope, the probability of failure decreases by 0.2% at the middle and by 0.5% at the toe of the slope. This is consistent with the observation of Aladejare and Wang [68] that when the correlation between mechanical parameters is ignored, the probability of failure obtained from the reliability analysis might be overestimated.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Conclusion

In this study, we have examined the uncertainty in parameters. Fuzzy transform was performed on the cohesion and the angle of internal friction. Fuzzy point estimation was used in combination with the hydromechanical coupling model on HYDRUS 2D and the Slope Cube Module to examine the slope stability. The result shows that the fuzzy theory can effectively evaluate the fluctuation interval, mean, and standard deviation of the factor of safety and the reliability index. The probability of failure in the interior of the slope was computed through reliability analysis. At our observation points on the loam slope, the fuzzy reliability of loam was determined to be 1.3357, and the probability of failure was 0.0908. For silt, the fuzzy reliability was observed to be 2.2299, and the probability of failure was 0.0129. The results of the slope failure mechanism investigation is that, after rainfall infiltrates into the soil, the change in water content causes an increase in suction stress (a decrease in its absolute value). The resulting decrease in soil effective stress leads to slope instability. It has been determined in this study that the change in the probability of failure is spatially related to the depth of the moisture band caused by the soil hydraulic conductivity and to the suction stress change controlled by the water content. After 48 hr of rainfall, the infiltration depth into the loam slope was deeper than that into the silt slope. The area of the loam slope in which the probability of failure exceeded 50% was approximately twice as large as that of the silt slope. It suggests that, as the rainfall infiltrates deeper, the area of instability in the slope increases. This study was also aimed at determining the effect of correlation between the parameters on the probability of failure. It was shown that a stronger negative correlation between the mechanical parameters yields a lower calculated probability of failure when performing slope stability analysis. When the correlation was considered, the computed probability of failure at observed points decreased by <1%. It suggests that the correlation between parameters may be ignored when a conservative estimate of slope stability is required.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors are grateful for the support of the Research Project of the Ministry of Science and Technology, Taiwan (MOST 106-2625-M-006-014).