Effect of Unimodal and Bimodal Soil Hydraulic Properties on Slope Stability Analysis

: Rainfall inﬁltration is the primary triggering factor of slope instability. The process of rainfall inﬁltration leads to changes in the water content and internal stress of the slope soil, thereby affecting slope stability. The soil water retention curve (SWRC) was used to describe the relationship between soil water content, matric suction, and the water retention characteristics of the soil. This characteristic is essential for estimating the properties of unsaturated soils, such as unsaturated hydraulic conductivity function and shear strength. Thus, SWRC is regarded as important information for depicting the properties of unsaturated soil. The SWRC is primarily affected by the soil pore size distribution (PSD) and has unimodal and bimodal features. The bimodal SWRC is suitable for soils with structural or dual-porous media. This model can describe the structure of micropores and macropores in the soil and allow the hydraulic behavior at different pore scales to be understood. Therefore, this model is more consistent with the properties of onsite soil. Few studies have explored the differences in the impact of unimodal and bimodal models on unsaturated slopes. This study aims to consider unimodal and bimodal SWRC to evaluate the impact of unsaturated slope stability under actual rainfall conditions. A conceptual model of the slope was built based on ﬁeld data to simulate changes in the hydraulic behavior of the slope. The results of seepage analysis show that the bimodal model has a better water retention capacity than the unimodal model, and therefore, its water storage performance is better. Under the same saturated hydraulic conductivity function, the wetting front of the bimodal model moves down faster. This results in changes in the pressure head, water content, and internal stress of the soil. The results show that the water content and suction stress changes of the bimodal model are higher than those of the unimodal model due to the difference in water retention capacity. Based on the stability of the slope, calculated using the seepage analysis, the results indicate that the potential failure depth of the bimodal model is deeper than that of the unimodal model. Furthermore, the SWRC parameters ﬁtted by ﬁeld data were used to simulate the conceptual slope model in this study. Before the simulation, the soil water content measured in the ﬁeld was calibrated to prove the rationality of the model. According to the observation data of the Babaoliao rainfall station (1 June 2019–1 August 2019), there was high-intensity rainfall from 13 June 2019 to 14 June 2019. Thus, this rainfall event was used to calibrate the soil water content. Based on the water content of ﬁeld observations and simulations, the value at a depth of 0.3 m is used as the basis for the calibration. The simulation results of the unimodal and bimodal SWRC parameters show that the change in water content is similar to the trend of the ﬁeld observations, as shown in Figure 5. After statistical calculation, the results showed that both models were consistent with the ﬁeld soil water content. The root mean square errors (RMSEs) were 0.0568 and 0.0057. The covariances were 0.000121 and 0.00015. The results demonstrate the rationality of the model used in this study. In the next section, seepage analysis is performed on the ﬁtted SWRC parameters. The difference between the two SWRC models is quantiﬁed and qualitatively analyzed in the subsequent section. in a decrease in the LFS. However, the LFS of the bimodal model decreased to a greater extent. The results show that under the same rainfall conditions, the slope of the bimodal model is considered to be affected by rainfall-induced instability to a greater degree. The slope of the unimodal model is considered to have a smaller range of inﬂuence and degree of change.


Introduction
In recent years, many studies have analyzed and explored rainfall-induced unsaturated slope failure and instability; thus, unsaturated soil mechanics are essential for relevant research [1][2][3]. Many studies have specifically focused on the hydraulic behavior of the soil water retention curve (SWRC), which contains basic unsaturated soil properties and can determine the relationship between soil water content and matric suction. The shape of the SWRC is primarily affected by the different soils and their pore size distribution (PSD). They are generally divided into unimodal and bimodal characteristics. Well-graded soil generally has a unimodal pore distribution, and its SWRC shows a unimodal characteristic, which means that the soil pore size is a single continuous pore domain. Residual or colluvial soils have a double pore structure, which can be divided into micropores in the aggregates and macropores between the aggregates. These SWRC usually have a bimodal characteristic and are related to the capillary theory [4,5]. In addition, related studies have pointed out that the soil pore size distribution has a significant impact on the water retention behavior of different pores under unsaturated conditions [6]. Scholars have proposed empirical and physical equations to describe SWRC to predict the hydraulic characteristics of unsaturated soils [7][8][9][10][11][12]. Furthermore, to explore the hydraulic properties of the dual-porous media in onsite soil, some scholars have proposed methods to quantify the interaction between the dual-porous structure (macropores and micropores) and the distribution of pores in the soil and preferential flow [5,[13][14][15]. In these methods, the bimodal SWRC was superimposed on two unimodal curves. The subcurves were related to the interaggregate pores and intra-aggregate pores of the soil. Many studies have confirmed that the SWRC can be used to calculate the soil hydraulic conductivity function (HCF) using the pore size distribution function [11,[16][17][18].
Previous studies have indicated that rainfall is the main triggering factor for slope failure. In unsaturated soil, rainfall infiltration increases the degree of saturation, which affects the shear strength characteristics, soil internal stress, and slope stability [3,[19][20][21][22]. The relevant literature has indicated that the accuracy of soil hydraulic parameters is essential for evaluating slope stability in the coupled hydromechanical model due to the limitations of field soil samples and measurement errors. The actual condition of the slope can be reflected by the field measurement, and the prediction model can be effectively calibrated using the estimated hydraulic parameters [23][24][25][26]. Therefore, the accuracy of SWRC fitting will impact the simulation and prediction of the coupled hydromechanical model. If a difference is observed between the experimental data and the best fitting curve (unimodal), the SWRC may have a bimodal characteristic. At this time, the bimodal SWRC can be fitted to a more accurate curve and can be used to accurately estimate the hydraulic conductivity function and shear strength of unsaturated soil [27,28].
Many studies have used experiments to measure SWRC to obtain the hydraulic properties of unsaturated soils. It is also used to estimate the shear strength and hydraulic conductivity function of unsaturated soil to solve the soil hydraulic properties required to evaluate slope instability caused by rainfall [29][30][31][32]. Lu et al. [33] proposed a relationship between suction stress and the effective degree of saturation. The SWRC was used to describe the contribution of matric suction to the evaluation of unsaturated soil shear strength. The stress state of the unsaturated soil is expressed using the suction stress characteristic curve (SSCC). The suction stress is regarded as the energy stored in the soil unit [34,35]. When rainfall infiltrates the soil, it can easily cause changes in the internal hydraulic behavior of the slope, which is the key to the instability of shallow slopes. However, few studies have considered the impact of rainfall infiltration on the internal stress of unsaturated slopes in dual-porous media. The stability of the slope is estimated by using the bimodal SWRC and SSCC to predict the change between the water content and the soil stress.
The limit equilibrium method (LEM) is mostly used in traditional slope stability analysis and research because of its validity and reliability [36][37][38][39][40][41]41]. In the LEM, the potential failure surface of the slope is preset, and the potential sliding soil is discretized into smaller vertical slices. The factor of safety (FS) is used as an indicator of slope stability based on the principle of force balance and moment balance. However, this method only seeks a single stability index for the slope as a whole. It is not possible to understand the changes in pore water pressure and effective stress caused by rainfall infiltration, changes in slope stability with temporal and spatial distribution, and location of initial failure. The local factor of safety (LFS) theory is based on the Mohr-Coulomb failure criterion that evaluates the stability of each point in the slope [42]. It is an extension of the method described by Iverson and Reid (1992) for viscous, variable saturation porous materials [43] and provides a physical value to evaluate the stability of a partially saturated slope. Additionally, the distribution of soil suction stress calculated by considering dualporous media under rainfall conditions can be used as the LFS and the degree of stability of each point in the hillslope. Therefore, the slope stability can be effectively analyzed using the LFS at different depths, and the potential failure locations can be evaluated. This method can effectively solve the challenges encountered in the traditional slope stability analysis. The slope stability and the soil hydraulic behavior are essential for the sustainable hillslope management and geotechnical engineering. The above approaches can be used to evaluate and predict the risk of landslides which may lead to a disaster for humans. The slope stability modeling can assess the risk and prevent hazard exposure for human well-being.
This study explores the influence of considering dual-porous media on the stability of unsaturated slopes caused by rainfall. First, the unimodal and bimodal SWRC were fitted based on field soil data and used to estimate the unsaturated hydraulic conductivity function. Second, a conceptual model of the slope was built based on the field data of the landslide event. We selected actual rainfall events and the measurement of field water content to calibrate the accuracy of the model. Then, unimodal and bimodal SWRC were used to analyze the slope seepage process and the differences between the two models were compared. Finally, according to the results of the seepage analysis, the slope stability was further evaluated to quantify the changes in the two models.

Conceptual Model of Hillslope
This study explores the influence of soil dual-porous media on the seepage and stability of unsaturated slopes caused by rainfall. A conceptual model of the slope profile based on the field hydrogeological data of the Babaoliao landslide area was established, as shown in Figure 1. The Babaoliao landslide area is located in Dongxing Village Zhongpu Township, Chiayi County, Taiwan. In the field, it is found that the hillslope is affected by creep and has many features, such as tension cracks and subsidence, and these cracks may form an excellent soil channel for rainfall infiltration when rainfall occurs. Therefore, this study considers the soil dual-porous media for evaluation. The conceptual model of the slope includes the upper colluvial soil and the lower sandstone, and there is a thin layer of surface soil covering the upper colluvial soil. The surface soil is formed by rainfall alluvial from the upper area and its particle is small. The colluvial is composed of the debris and deposits with larger particles. The surface of the slope is the boundary of rainfall infiltration, and the left and right sides are the constant head boundaries, respectively, and the others are the zero flow boundaries. The grid of the conceptual slope was divided into 17,261 nodes and 33,913 elements. Two observation wells (BH-05CI and BH-02) and the observation profiles are set in the slope. We mainly intend to explore the hydraulic behavior of the shallow soil layer; thus, it will be especially discussed in the hydraulic behavior of the observation profile as shown in Figure 1 to avoid the influence of boundary conditions. When the slope is infiltrated by rainfall, the observation profile is used to analyze the soil seepage situation and the influence of internal stress and stability. In this study, the rainfall event was based on the observation data of the Babaoliao rainfall station (1 June 2019-1 August 2019), which show that there was high-intensity rainfall from 13 June 2019 to 14 June 2019. Thus, the 48-h rainfall data were used as the rainfall input value of the model.
In this study, the slope soil materials and pressure cooker test data of soil sampling points refer to the report from the Soil and Water Conservation Bureau, Council of Agriculture, Executive Yuan (2017). The mechanical parameters of the three soils were determined from the field soil test results, as shown in Table 1. Field surveys showed that this area was dominated by silt-bearing sandstone, and the neighboring sites were affected by soil creep and had many features, such as tension cracks and subsidence. These cracks may become a channel for rainfall to infiltrate the soil and cause discontinuities in the distribution of the pressure head and water content inside the slope. This leads to changes and unstable development of slope soil heterogeneity and internal hydraulic behavior. Therefore, the bimodal SWRC considering the soil dual-porous media was adopted for subsequent analysis and discussion. In this study, the slope soil materials and pressure cooker test data of soil sampling points refer to the report from the Soil and Water Conservation Bureau, Council of Agriculture, Executive Yuan (2017). The mechanical parameters of the three soils were determined from the field soil test results, as shown in Table 1. Field surveys showed that this area was dominated by silt-bearing sandstone, and the neighboring sites were affected by soil creep and had many features, such as tension cracks and subsidence. These cracks may become a channel for rainfall to infiltrate the soil and cause discontinuities in the distribution of the pressure head and water content inside the slope. This leads to changes and unstable development of slope soil heterogeneity and internal hydraulic behavior. Therefore, the bimodal SWRC considering the soil dual-porous media was adopted for subsequent analysis and discussion.

Seepage Model of Unsaturated Soil
In this study, the HYDRUS 2D module [44] was used to evaluate and analyze the seepage process in unsaturated soil. This module can be used to analyze unsaturated soil slopes effectively. Based on the Richards equations, transient and two-dimensional seepage governing equations were developed to estimate the unsaturated soil layer and soil hydraulic properties. The seepage governing equation is as follows: where H [L] is the total pressure head, t [T] is the time, h [L] is the pressure head or suction head, and W [T −1 ] is the source or sink, which refers to seepage caused by pumping or infiltration recharge.
) h (  [-] is the volumetric water content with the change in the pore water pressure and is dimensionless. is the unsaturated hydraulic conductivity function. In this study, the SWRC proposed by van Genuchten (1980) was used to describe the relationship between soil water content and matric suction. This method is in

Seepage Model of Unsaturated Soil
In this study, the HYDRUS 2D module [44] was used to evaluate and analyze the seepage process in unsaturated soil. This module can be used to analyze unsaturated soil slopes effectively. Based on the Richards equations, transient and two-dimensional seepage governing equations were developed to estimate the unsaturated soil layer and soil hydraulic properties. The seepage governing equation is as follows: is the pressure head or suction head, and W [T −1 ] is the source or sink, which refers to seepage caused by pumping or infiltration recharge. θ(h) [-] is the volumetric water content with the change in the pore water pressure and is dimensionless. K(h) is the unsaturated hydraulic conductivity function. In this study, the SWRC proposed by van Genuchten (1980) was used to describe the relationship between soil water content and matric suction. This method is in good agreement with the experimental results for fitting the SWRC parameters [45]. Additionally, the hydraulic conductivity function was proposed by Mualem (1976) based on the SWRC, and both can be described by the following formulas: where θ s [[-] is the saturated soil water content and is dimensionless. θ r [-] is the residual soil water content and is dimensionless. h [L] is the matric suction, α [L −1 ] is related to the soil air entry value, n [-] is related to the soil pore size distribution and is dimensionless, and l [-] is the pore connectivity parameter and is dimensionless, generally expressed as is the effective saturation and is dimensionless. K s [LT −1 ] is the saturated hydraulic conductivity function. α, n, and l are the empirical parameters that define the shape of the hydraulic conductivity function.

Bimodal Soil Water Retention Curve
Some soils have a macrostructure and microstructure, which may exhibit bimodal characteristics in the SWRC and a bimodal pore size distribution. However, soil with a bimodal pore size distribution cannot use unimodal SWRC to fully describe the bimodal characteristics. Therefore, many studies have proposed a bimodal SWRC equation [11,16,17,46] to describe the hydraulic characteristics of macropores and micropores in structured soils. In this study, the bimodal SWRC proposed by Durner (1994) was used to fit onsite soil data. To correctly represent the water retention characteristics of soils with heterogeneous pore structures, a multipeak water retention function is overlapped by the subcurves (associated with macropores and micropores) in the van Genuchten model as illustrated in Equation (4): , and m i [-] are the empirical parameters of the two pore domains, respectively, and m = 1 − 1/n. w i [-] is the weighting factor of the single pore domain and is dimensionless ( The bimodal unsaturated hydraulic conductivity function is combined with the Mualem (1976) and van Genuchten (1980) models, and the bimodal SWRC can be written as follows [18]: where k is the mode of the model (k = 2, which is a double peak) and w i is the weighting factor for each pore domain. K(h i ) is the unsaturated hydraulic conductivity function defined by the steady-state flow method at pressure head of i = −1 hPa.

Mechanics Model of Unsaturated Soil
In this study, the finite element analysis model HYDRUS 2D combined with the slope cube module is used to calculate the unidirectional coupling variable saturation flow and stress problems. Based on the seepage analysis in the previous section, the slope cube module was used to analyze the changes in soil stress to describe the impact of the unsaturated layer on soil hydraulic behavior after rainfall infiltration. In the slope cube module, the finite element method (FEM2D) was adopted to solve the stress distribution of various points inside the slope based on the momentum balance [47]. This method is based on plane stress linear elasticity to simulate the stress change caused by the transient unit weight.
The linear elastic material total stress governing equation is given by Equation (6): where σ is the independent stress variable in two-dimensional space, γ is the unit weight of the slope soil material, and b is the unit vector of the body force.
In this study, the effective stress principle is consistent with Terzaghi's effective stress [48]. By modifying the contribution of the saturation to the effective stress, the effective stress is regarded as Bishop's extension [49,50] and Terzaghi's extension for all saturations [35]. The effective stress is given by Equation (7): , and S e [-] is the effective saturation and is dimensionless. σ s [ML −1 T −2 ] is the suction stress, which represents all the physical and chemical mechanisms that may occur between soil particles as illustrated in Equation (8) [35,51]: is the capillary force and σ pc [ML −1 T −2 ] is the combination of the van der Waals attractive force and electric double layer force. S [-] is the degree of saturation and is dimensionless, is the Born repulsive force, which is the interatomic resistance that prevents two contacting particles from penetrating each other. The Born repulsive force balances the capillary force, van der Waals attractive force, electric double layer force, and matric suction in the soil. Among them, the influence of the van der Waals force and electric double layer force is ignored as the soil particle size increases. The soil suction stress is primarily controlled by the water content because each stress component of the soil can be described as a function of matric suction and water content. Lu et al. [33], based on the theory of thermodynamics, regarded suction stress as the energy stored in a unit of soil. The suction stress characteristic curve (SSCC) was derived using the van Genuchten (1980) SWRC model. The derivation process is as follows:

Local Factor of Safety Theory
In this study, the local factor of safety (LFS) proposed by Lu et al. [42] was used to evaluate slope stability. This method is based on the Mohr-Coulomb failure criterion and the process of changing the stress state toward the failure condition caused by rainfall. It is used to evaluate the current stress and stable state of each location point inside the slope. This method differs from the traditional limit equilibrium method (LEM). In the LEM, each slice is given a safety factor (FS) value as a possible potential failure surface in the slope. However, the LFS can calculate the FS of each point to find the potential failure location, rather than being limited to a single slice. The soil strength of the current stress state can be estimated by the intercept of the Mohr-Coulomb failure envelope to define the LFS of each point in the slope. The LFS is represented by the ratio of the potential stress to the current stress, as shown in Equation (11): where τ * is the shear strength which represents the value of the potential Coulomb stress. τ is the shear stress which represents the current Coulomb stress. ϕ is the effective friction angle, c is the effective cohesion, σ 1 is the maximum principal stress, σ 3 is the minimum principal stress, and σ s is the suction stress, as shown in Figure 2. When the LFS is less than 1, this means that the point may fail. When the LFS is greater than 1, the point is relatively stable. This explains the locations of potential failure areas and relatively stable areas. This method cannot only overcome the challenges of traditional slope stability analysis effectively but also allow the influence of the stability change of the slope at each position after rainfall infiltration to be analyzed.
where *  is the shear strength which represents the value of the potential Coulomb stress.  is the shear stress which represents the current Coulomb stress. '  is the effective friction angle, ' c is the effective cohesion, 1  is the maximum principal stress, 3  is the minimum principal stress, and s  is the suction stress, as shown in Figure 2. When the LFS is less than 1, this means that the point may fail. When the LFS is greater than 1, the point is relatively stable. This explains the locations of potential failure areas and relatively stable areas. This method cannot only overcome the challenges of traditional slope stability analysis effectively but also allow the influence of the stability change of the slope at each position after rainfall infiltration to be analyzed.

Results of Soil Hydraulic Fitting Parameters and Field Data Calibration
In this study, three types of soil data measured in the Babaoliao landslide area were used: soil, colluvial soil, and sandstone. The unimodal and bimodal SWRC were fitted based on their water content and matric suction, and the differences between the two models were compared. The SWRC parameters are obtained through the unimodal Equation (2) and the bimodal Equation (4) based on the field measurement information, as shown in Table 2. The fitting results of the three soils show that the bimodal SWRC has high accuracy, can completely fit the field measurement data, and present a bimodal shape, as shown in Figure 3. Fittings to experimental data are good for both unimodal and bimodal at the range of matric suction less than 1000 kPa. A difference appears at a high matric suction range (larger than 1000 kPa) in which there is no experimental data. It causes the phenomenon of uncertainty in high matric suction range, which affected the hydraulic behavior of the micropores. The values of saturated water content and residual water content are slightly different for unimodal and bimodal functions due to insufficient experimental data at the low and high matric suction range. In this study, the saturated hydraulic conductivity functions of the two models were the same. The SWRC parameters were fitted into Equations (3) and (5) to estimate the unsaturated hydraulic conductivity functions of the three soils in the unimodal and bimodal models, as shown in Figure 4. The bimodal SWRC is affected by different pore size distributions; therefore, the permeability will be different in the micropores and macropores regions, resulting in a difference between the unimodal and bimodal hydraulic conductivity functions. The results show that under the same saturated hydraulic conductivity function, the unsaturated soil in the bimodal model has a higher hydraulic conductivity function. Seepage and stability analyses will be affected by this result.     Furthermore, the SWRC parameters fitted by field data were used to sim conceptual slope model in this study. Before the simulation, the soil water cont ured in the field was calibrated to prove the rationality of the model. Accord observation data of the Babaoliao rainfall station (1 June 2019-1 August 2019), Furthermore, the SWRC parameters fitted by field data were used to simulate the conceptual slope model in this study. Before the simulation, the soil water content measured in the field was calibrated to prove the rationality of the model. According to the observation data of the Babaoliao rainfall station (1 June 2019-1 August 2019), there was high-intensity rainfall from 13 June 2019 to 14 June 2019. Thus, this rainfall event was used to calibrate the soil water content. Based on the water content of field observations and simulations, the value at a depth of 0.3 m is used as the basis for the calibration. The simulation results of the unimodal and bimodal SWRC parameters show that the change in water content is similar to the trend of the field observations, as shown in Figure 5. After statistical calculation, the results showed that both models were consistent with the field soil water content. The root mean square errors (RMSEs) were 0.0568 and 0.0057. The covariances were 0.000121 and 0.00015. The results demonstrate the rationality of the model used in this study. In the next section, seepage analysis is performed on the fitted SWRC parameters. The difference between the two SWRC models is quantified and qualitatively analyzed in the subsequent section. Furthermore, the SWRC parameters fitted by field data were used to simulate the conceptual slope model in this study. Before the simulation, the soil water content measured in the field was calibrated to prove the rationality of the model. According to the observation data of the Babaoliao rainfall station (1 June 2019-1 August 2019), there was high-intensity rainfall from 13 June 2019 to 14 June 2019. Thus, this rainfall event was used to calibrate the soil water content. Based on the water content of field observations and simulations, the value at a depth of 0.3 m is used as the basis for the calibration. The simulation results of the unimodal and bimodal SWRC parameters show that the change in water content is similar to the trend of the field observations, as shown in Figure 5. After statistical calculation, the results showed that both models were consistent with the field soil water content. The root mean square errors (RMSEs) were 0.0568 and 0.0057. The covariances were 0.000121 and 0.00015. The results demonstrate the rationality of the model used in this study. In the next section, seepage analysis is performed on the fitted SWRC parameters. The difference between the two SWRC models is quantified and qualitatively analyzed in the subsequent section.

Comparison of the Influence of Unimodal and Bimodal Models on Seepage Analysis
This section discusses the effects of unimodal and bimodal SWRC on unsaturated soils after rainfall infiltration. The two models were analyzed in terms of the amount of rainfall infiltration and the degree of change in the pressure head. In this study, unimodal and bimodal SWRC parameters were used for seepage analysis. The change in rainfall infiltration and pressure head were simulated in the slope over a time of 48 h. In the actual rainfall event used in this study, the average rainfall intensity is 2.2 mm/h, the total rainfall volume is 106 mm within 48 h, and the maximum rainfall intensity is 29.5 mm/h at 26 h. First, this study focuses on rainfall infiltration at the atmospheric boundary (surface). The change in rainfall infiltration over time and the cumulative rainfall infiltration are shown in Figure 6. In Figure 6, it can be observed that under the same rainfall conditions, the rainfall infiltration of the two models is different. The rainfall infiltration in the bimodal model was relatively larger than that in the unimodal model, and it showed a significant difference with the maximum rainfall intensity at 26 h. This phenomenon is mainly affected by unsaturated hydraulic conductivity. The bimodal unsaturated hydraulic conductivity function is larger than that of the unimodal. Hence, the water flow rate was faster, and the soil water retention capacity was better in the bimodal model after rainfall infiltration. This result is attributed to the hydraulic behavior of the macropores and micropores being considered in the bimodal model. The better permeability in the macropores leads to the difference of the unsaturated hydraulic conductivity function in the two models. This result affects the subsequent seepage and stability analyses. In addition, the change in the pressure head was explored in the slope. In the initial state (t = 0 h), the pressure head is linearly distributed with the depth of the slope. According to rainfall events, the simulation results were divided into 6, 12, 24, 36, and 48 h. The pressure head changes in the observation profile of the unimodal and bimodal SWRC parameters are observed, as shown in Figure 7. Figure 7a,b shows the pressure head changes in the unimodal and bimodal models, respectively. The rainfall event selected in this study had heavy rainfall from 23 to 36 h; as a result, the pressure head changed significantly after 24 h. The slope is infiltrated by rainfall, which causes the wetting front to move downward, resulting in a drop in the pressure head. The results show that the pressure head drop on the unimodal and bimodal models primarily occurs in shallow areas under the same rainfall conditions. This is attributed to the poor permeability of the soil, and the deep part of the soil was not affected. As shown in Figure 7a, after 48 h of rainfall, the wetting front of the unimodal model moved down to a depth of approximately 1 m. The wetting front of the bimodal model moves down to a depth of approximately 1.2 m, and no longer affects the deeper soil, as shown in Figure 7b. The two models have different speeds and depths of the wetting front during the wetting process, which means that the rainfall behavior for slope seepage analysis should be related to the SWRC and hydraulic conductivity function. From Figure 8, it can be observed that the pressure head of the bimodal model changes in deeper areas. Although the saturated hydraulic conductivity functions of the two models are the same, the unsaturated hydraulic conductivity functions calculated by Equations (3) and (5) are different, which leads to this situation. Therefore, under the conditions of the same rainfall intensity, time and saturated hydraulic conductivity function, the wetting front of the bimodal model moves downward faster, and the pressure head dissipates faster in the soil. Evidently, the soil heterogeneity considered by the bimodal model affects the change in seepage in the slope.

Influence of Rainfall Conditions on Unimodal and Bimodal Models
Under the same rainfall conditions, the speed and depth of the wetting front estimated by the unimodal and bimodal SWRC parameters will have an impact. Figure 8 shows the pressure head curve of the slope observation profile in the unimodal and bimodal models at 48 h. The soil water retention capacity of the bimodal model was larger than that of the unimodal model. The wetting front of the bimodal model moves down faster, which is approximately 0.2 m deeper than that of the unimodal model. This result makes it easier for the pressure head to dissipate. This means that the difference in seepage behavior of different models on slopes dominates the stability of the slope.
This study explores the change in the pressure head over time at a single depth in the observation profile, as shown in Figure 9. It shows the pressure head curve at a depth of 0.85 m in the slope observation profile of unimodal and bimodal models. The value of the pressure head decreased with increasing rainfall time. The pressure head curve indicates

Influence of Rainfall Conditions on Unimodal and Bimodal Models
Under the same rainfall conditions, the speed and depth of the wetting front estimated by the unimodal and bimodal SWRC parameters will have an impact. Figure 8 shows the pressure head curve of the slope observation profile in the unimodal and bimodal models at 48 h. The soil water retention capacity of the bimodal model was larger than that of the unimodal model. The wetting front of the bimodal model moves down faster, which is approximately 0.2 m deeper than that of the unimodal model. This result makes it easier for the pressure head to dissipate. This means that the difference in seepage behavior of different models on slopes dominates the stability of the slope.

Comparison of the Influence of Unimodal and Bimodal Models on Slope Stability
The pressure head and soil water content distribution values obtained from the seepage analysis were used to calculate and evaluate the slope stability in this section. In this study, the LFS was adopted to calculate the stress state and the degree of stability at different points in the slope. It can effectively overcome the challenges encountered in the traditional slope stability analysis. The stability of the slope is mainly affected by rainfall infiltration, which results in a change in the pressure head and an increase in water content. The amount and distribution of water in the soil affect the change in stress and cause the stability of the slope to change and decline. Therefore, this section mainly discusses the relationship between soil water content, suction stress, and the LFS and analyzes the changes in the above three factors in the unimodal and bimodal models with the rainfall conditions. This study explores the change in the pressure head over time at a single depth in the observation profile, as shown in Figure 9. It shows the pressure head curve at a depth of 0.85 m in the slope observation profile of unimodal and bimodal models. The value of the pressure head decreased with increasing rainfall time. The pressure head curve indicates that the pressure head of the unimodal model drops from −22.8 m to −8.56 m over time, representing a total decrease of 14.24 m, and the rate of change is 62.4%. The pressure head of the bimodal model drops from −22.8 m to −1.85 m over time, representing a total decrease of 20.95 m, and the rate of change is 91.8%. From the above results, it can be seen that the pressure head simulated by the bimodal model has a greater degree of change than the simulation of the unimodal model. The drop in the pressure head between the unimodal and bimodal models can be observed in Figure 8. As the rainfall time increased, there was a difference in the drop in the pressure head between the two models. The above results are mainly affected by soil hydraulic properties because the unsaturated hydraulic conductivity function estimated by the bimodal SWRC parameters has a higher value. This situation allows the rainfall to flow faster in the soil and results in a greater drop in the pressure head. Thus, the hydraulic conductivity function in the unsaturated state has a significant influence on slope seepage analysis. If the hydraulic conductivity function that is consistent with the field situation can be accurately obtained, the accuracy of the simulation can be improved.

Comparison of the Influence of Unimodal and Bimodal Models on Slope Stability
The pressure head and soil water content distribution values obtained from the seepage analysis were used to calculate and evaluate the slope stability in this section. In this study, the LFS was adopted to calculate the stress state and the degree of stability at different points in the slope. It can effectively overcome the challenges encountered in the traditional slope stability analysis. The stability of the slope is mainly affected by rainfall infiltration, which results in a change in the pressure head and an increase in water content. The amount and distribution of water in the soil affect the change in stress and cause the stability of the slope to change and decline. Therefore, this section mainly discusses the relationship between soil water content, suction stress, and the LFS and analyzes the changes in the above three factors in the unimodal and bimodal models with the rainfall conditions.

Comparison of the Influence of Unimodal and Bimodal Models on Slope Stability
The pressure head and soil water content distribution values obtained from the seepage analysis were used to calculate and evaluate the slope stability in this section. In this study, the LFS was adopted to calculate the stress state and the degree of stability at different points in the slope. It can effectively overcome the challenges encountered in the traditional slope stability analysis. The stability of the slope is mainly affected by rainfall infiltration, which results in a change in the pressure head and an increase in water content. The amount and distribution of water in the soil affect the change in stress and cause the stability of the slope to change and decline. Therefore, this section mainly discusses the relationship between soil water content, suction stress, and the LFS and analyzes the changes in the above three factors in the unimodal and bimodal models with the rainfall conditions. First, we can observe the change in soil water content in the unimodal and bimodal models over time, as shown in Figure 10a. The water content of the unimodal model increased from 0.257 to 0.293, and the rate of change was 14%. The water content of the bimodal model increased from 0.258 to 0.325, and the rate of change was 25.9%. The results show that the increase in the water content in the bimodal model is greater than that in the unimodal model. This result will affect the change in the internal stress of the soil, which represents the change in the soil suction stress in the unimodal and bimodal models with rainfall, as shown in Figure 10b. The result shows that the suction stress of the unimodal model drops from −155.2 kPa to −64.7 kPa, and the rate of change is 58.3%, whereas the suction stress of the bimodal model drops from −151.7 kPa to −15.5 kPa, and the rate of change is 89.7%. The soil of the bimodal model has better water capacity and can store more rainfall; therefore, the water content of the slope increases, and the soil suction stress decreases to a greater extent.
In this study, LFS was calculated and analyzed for slope stability. Figure 11 shows the variation in the LFS of the unimodal and bimodal models with depth at 48 h. From the results, it can be observed that both unimodal and bimodal slopes are unstable in the shallow layer (the LFS is less than 1), indicating the possibility of failure. The potential failure depths of the slopes in the unimodal and bimodal models are approximately 0.6 and 0.8 m, respectively. Additionally, we analyzed the stability of a single depth, as shown in Figure 12. It shows the variation of the LFS with the rainfall time at a depth of 0.85 m. The LFS of the unimodal model decreased from 12.16 to 3.38, and the change rate was 72.1%. The LFS of the bimodal model decreased from 11.32 to 1.40, and the change rate was 87.6%. After the slope was infiltrated by rainfall, the soil water content increased. The downward movement of the wetting front caused the pressure head and suction stress to decrease, resulting in a decrease in the LFS. However, the LFS of the bimodal model decreased to a greater extent. The results show that under the same rainfall conditions, the slope of the bimodal model is considered to be affected by rainfall-induced instability to a greater degree. The slope of the unimodal model is considered to have a smaller range of influence and degree of change.
with rainfall, as shown in Figure 10b. The result shows that the suction stress of the unimodal model drops from −155.2 kPa to −64.7 kPa, and the rate of change is 58.3%, whereas the suction stress of the bimodal model drops from −151.7 kPa to −15.5 kPa, and the rate of change is 89.7%. The soil of the bimodal model has better water capacity and can store more rainfall; therefore, the water content of the slope increases, and the soil suction stress decreases to a greater extent. In this study, LFS was calculated and analyzed for slope stability. Figure 11 shows the variation in the LFS of the unimodal and bimodal models with depth at 48 h. From the results, it can be observed that both unimodal and bimodal slopes are unstable in the shallow layer (the LFS is less than 1), indicating the possibility of failure. The potential failure depths of the slopes in the unimodal and bimodal models are approximately 0.6 and 0.8 m, respectively. Additionally, we analyzed the stability of a single depth, as shown in Figure 12. It shows the variation of the LFS with the rainfall time at a depth of 0.85 m. The LFS of the unimodal model decreased from 12.16 to 3.38, and the change rate was 72.1%. The LFS of the bimodal model decreased from 11.32 to 1.40, and the change rate was 87.6%. After the slope was infiltrated by rainfall, the soil water content increased. The downward movement of the wetting front caused the pressure head and suction stress to decrease, resulting in a decrease in the LFS. However, the LFS of the bimodal model decreased to a greater extent. The results show that under the same rainfall conditions, the slope of the bimodal model is considered to be affected by rainfall-induced instability to a greater degree. The slope of the unimodal model is considered to have a smaller range of influence and degree of change.
These results are consistent with the expected outcomes. The bimodal unsaturated soil hydraulic conductivity function was relatively large, which resulted in a faster rainfall infiltration rate. This affects the downward movement rate of the wetting front. Considering infiltration, the surface infiltration volume of the bimodal model is greater than that of the unimodal model. This causes the soil water content to increase further and leads to instability in the bimodal model. Thus, the bimodal SWRC will affect the stability of the slope.

Conclusions
In this study, the soil dual-porous media was used to evaluate the SWRC parameters. The bimodal SWRC was adopted to compare and discuss the unimodal SWRC. The results were used to quantify the influence of soil dual-porous media on the slope stability caused These results are consistent with the expected outcomes. The bimodal unsaturated soil hydraulic conductivity function was relatively large, which resulted in a faster rainfall infiltration rate. This affects the downward movement rate of the wetting front. Considering infiltration, the surface infiltration volume of the bimodal model is greater than that of the unimodal model. This causes the soil water content to increase further and leads to instability in the bimodal model. Thus, the bimodal SWRC will affect the stability of the slope.

Conclusions
In this study, the soil dual-porous media was used to evaluate the SWRC parameters. The bimodal SWRC was adopted to compare and discuss the unimodal SWRC. The results were used to quantify the influence of soil dual-porous media on the slope stability caused by rainfall and to compare them with homogeneous soils. The conclusions of this study can be summarized as follows: 1.
Considering the dual-porous soil media, the fitting SWRC parameters of the field soil yielded better results. This result describes the changes in the hydraulic behavior of dual-porous soil media caused by rainfall. It can also be used to quantify and predict the seepage process and stability status of unsaturated slopes.

2.
Under the same rainfall conditions, the infiltration of the surface soil of the bimodal model was larger than that of the unimodal model. This means that the soil water retention capacity of bimodal soil is greater than that of the unimodal model, which is related to the hydraulic conductivity function of unsaturated soil. It is attributed to the fact that the hydraulic behavior of the macropores and micropores are considered in the bimodal model, leading to the difference of the unsaturated hydraulic conductivity function in the two models. Under the same saturated hydraulic conductivity function, the unsaturated hydraulic conductivity function of the bimodal model was larger than that of the unimodal model, implying that the water flow in the soil of the bimodal model was faster than the water flow in the unimodal model.

3.
After the bimodal model is infiltrated by rainfall, the wetting front moves down faster, which affects the deeper soil layer. Rainfall infiltration increases the soil water content, leading to changes in the internal stress and stability of the slope.

4.
From the results of this study, it can be found that if the field soil has structural or preferential flow characteristics, the bimodal model can be used for analysis. The soil dual-porous media can be considered in this model, and the results are consistent with the actual situation. Therefore, this model can be used to analyze engineering evaluations based on the characteristics of field soil in the future.

Funding:
The authors are grateful for support from the Headquarters of University Advancement at the National Cheng Kung University, sponsored by the Ministry of Education, Taiwan. This research received no external funding.