Moisture content distribution in soil is of great importance for understanding rainfall-induced slope failure and roadbed settlement. This study aims to develop a moisture migration model that improves the quantification of moisture content at an arbitrary point of the soil at any time for the whole process of infiltration under steady rainfall. The model was derived from the Richards equation using the flux-concentration relation, which was validated by numerical solutions calculated by Hydrus-1D software to evaluate the performance of the model. Results showed good accuracy and high adaptability for the moisture migration simulation of a wide range of soil types, and is applicable for short-duration and long-duration steady rainfall. Moreover, it can also reflect the stratification phenomenon for soil profile wetting by infiltration. Our analysis indicates that the flux and surface volumetric moisture content together can bound the boundary conditions of rainfall infiltration, which presents a shift from constant-flux to constant-concentration during long-duration steady rainfall. The migration rate of the wetting front in the later stage of infiltration positively correlates with rainfall intensity under the constant-flux condition, while it finally stabilizes at Ks/(θs − θi) under the constant-concentration condition (i.e., Ks – saturated hydraulic conductivity, θs – saturated volumetric moisture content, θi – initial volumetric moisture content).
Moisture migration model can improve the quantification of moisture content at an arbitrary point of the soil at any time for the whole process of infiltration under steady rainfall.
The model shows good accuracy and high adaptability for the moisture migration simulation of a wide range of soil types, and is applicable for short-duration and long-duration steady rainfall.
Rainfall is one of the most significant triggering factors of slope failure and roadbed settlement (Lee et al. 2009; Hedayati & Hossain 2015; Wu et al. 2015). In general, roadbeds or slopes consist of unsaturated soils with low initial moisture content, high matrix suction, and good stability. However, with the increase of wetting front depth, long-duration rainfall leads to matric suction decrease and results in the reduction of shear strength for unsaturated soils (Fang & Esaki 2012), and in addition it leads to the poor durability, lack of stability, and uneven settlement of the roadbed. The variation of moisture content under rainfall infiltration has been established, and implemented in different models to estimate the softened shear strength of soils (Jacobs et al. 2004; Li et al. 2005), resilient modulus of bound material (Khoury & Zaman 2004; Cary & Zapata 2011), and volumetric deformation of subgrade soil (Lytton et al. 2005). Hence, this study investigates the moisture migration law of the unsaturated soil profile under rainfall, i.e., moisture content distribution or wetting front depth, which are of crucial significance.
Nowadays, the investigation of the influence of rainfall infiltration on moisture content distribution has three approaches: field monitoring, numerical simulation, and analytical methods. Field monitoring is helpful for the study of rainfall infiltration effects on slope stability, however, it is a costly procedure (Zhan et al. 2013). The Richards equation (RE) (Richards 1931) is the governing equation to describe rainfall infiltration through unsaturated soil, which is the theoretical basis for numerical simulation and analytical analysis. Nonetheless, it is a highly nonlinear partial differential equation that is difficult to solve directly. In the numerical simulation approach, many scholars have proposed methods to improve the accuracy of numerical solutions of RE. These improvements have been integrated into software called Hydrus-1D. Although numerical solutions of RE are available for many conditions and have high accuracy, these are not universally robust, and thus are not expected to immediately replace the use of analytical solutions in hydrological models (Smith et al. 2002).
Among the analytical methods, the Philip infiltration model (Philip 1955; Philip 1957; Philip 1969) and Parlange infiltration model (Parlange 1971; Parlange 1972; Parlange et al. 1982; Parlange et al. 1985), both derived from RE, have been widely used, and apply the iterative method to approximate analytical solution of the RE under the initial and boundary conditions. Srivastava & Yeh (1991) developed an analytical solution to the linearized RE for vertical infiltration, which provides a reliable means to compare the accuracy of various numerical methods. In terms of flux-concentration relation (FCR) solution, Philip (1973) was the first to introduce the concept of FCR and proposed that, if the expression of FCR were exactly known for any given phenomenon, i.e., given boundary conditions and given soil characteristics, the relevant exact solution of RE could be found immediately. (Kutilek 1980) used the approximate expression of FCR to derive a model for the moisture content distribution in wetting profiles under vertical rainfall infiltration, although the model is only applicable for a constant flux boundary and short-term rainfall.
To date, analytical solutions of RE have been derived in many studies, however, few analytical solutions have been specifically derived to apply to the whole process of rainfall infiltration for a more accurate moisture distribution in wetting profiles. Rainfall infiltration is a dynamic process with the boundary conditions varying with rainfall duration. It can be observed that the surface moisture content gradually stabilizes with the increase of rainfall duration under a constant-flux boundary (i.e., as shown in Section 2.1.2), which has also been proved by experimental studies (Zeng et al. 2018). When the rainfall intensity is greater than the saturated hydraulic conductivity, ponding starts to occur and the upper boundary condition gradually shifts, i.e., the infiltration curve splits into pre-ponding and post-ponding stages (Zhu et al. 2006; Chiu et al. 2009; Chen et al. 2014). With an additional boundary condition, the initial condition gradually changes (i.e., the initial moisture content of the soil profile changes from uniform to non-uniform). In addition, the stratification phenomenon for the moisture distribution in wetting profiles comes into being gradually during rainfall infiltration (Colman & Bodman 1945). Therefore, the above analytical solutions of RE are not satisfactory to accurately describe the moisture content distribution in wetting profiles for the whole process of rainfall infiltration.
In this study, the boundary conditions satisfied by the whole process of infiltration under steady rainfall are analyzed by applying the conclusions of existing infiltration experiments, and the uncertainty of the demarcation point is confined on account of continuity in the transition from different infiltration boundary conditions (i.e., as shown in Section 4.3). In addition, the basic theory of FCR is introduced to derive semi-analytical solutions of RE for the whole infiltration process under rainfall conditions. The primary objectives of this study are thereby the following: (1) derive a moisture migration model (i.e., a time evolution of moisture content distribution and wetting front depth) which is applicable for short-duration and long-duration rainfall under a rainfall intensity that is greater than or less than the saturated hydraulic conductivity; (2) validate the accuracy of the moisture migration model by numerical simulation to provide robust solutions for the moisture migration simulation of different soil types under various rainfall intensities; (3) explore the relationship between the duration of different infiltration stages and the rainfall intensity during the whole process of rainfall infiltration.
Governing equation and boundary conditions
It is generally assumed that when the rainfall intensity is less than the saturated hydraulic conductivity (R < Ks) under constant rainfall intensity, no ponding is formed on the surface, and all of the rainwater is absorbed into the soil (i.e., the constant-flux boundary). The conclusions of existing infiltration experiments show that the surface volumetric moisture content gradually stabilizes (Zeng et al. 2018), and the different infiltration boundary then comes near the stage of stabilization. Furthermore, it is more reasonable if boundary conditions are synthetically considered between infiltration rate and surface volumetric moisture content. Next, the Type 1 critical point is defined (t*,), where t* is called Type 1 critical time and is called stable moisture content on the surface. The infiltration process can then advance into two stages of non-pressure infiltration transition (NPIT) and non-pressure infiltration stabilization (NPIS).
The ponding point is generated on the surface (where the ponding time is recorded as tp) when the rainfall intensity is greater than the saturated hydraulic conductivity (R > Ks), which is determined by the moisture content on the surface, i.e., θ = θs (Kutilek 1980). The surface is still the constant-flux boundary before ponding occurs, while it turns into the constant-concentration boundary after ponding. Based on the conclusion of existing infiltration experiments (Chen et al. 2014), the infiltration rate gradually tends to the saturated hydraulic conductivity under the constant-concentration boundary. The Type 2 critical point is then defined (ts,), where ts is called the Type 2 critical time or saturated time. The Type 2 critical point after ponding can simplify the infiltration process into two stages, i.e., transition and stabilization. Therefore, tp and ts reasonably divide the whole process of rainfall infiltration into three stages under rainfall intensity that is greater than the saturated hydraulic conductivity, i.e., non-pressure infiltration transition (NPIT), pressure infiltration transition (PIT), and pressure infiltration stabilization (PIS). The relevant boundary conditions are shown in Figure 2.
Moisture migration model of soil profile
The flux-concentration relation
Moisture migration model
R < Ks
R > Ks
In order to evaluate the performance of the proposed moisture migration model, it was compared with the solutions of RE under two types of rainfall conditions. For this purpose, three soils (Vogel et al. 2000) (i.e., sandy loam, loam, and silt) were selected for model evaluation, whose van Genuchten soil hydraulic parameters are shown in Table 1.
|Soil texture .||θr .||θs .||α/mm−1 .||N .||Ks/mm·h−1 .||l .|
|Soil texture .||θr .||θs .||α/mm−1 .||N .||Ks/mm·h−1 .||l .|
The Richards equation illustrated in Equation (1) was solved numerically using Hydrus-1D (Šimůnek et al. 2005). The basic principle of Hydrus-1D is to use the Galerkin finite element method to solve the Richards equation for soil moisture migration simulation under rainfall infiltration conditions. The Hydrus-1D simulation results are generally widely accepted, and can be used for evaluating the approach of this study.
For numerical simulation, rainfall intensities of 50 mm/h and 30 mm/h were selected for sandy loam, 15 mm/h and 8 mm/h for loam, and 5 mm/h and 2 mm/h for silt. In addition, the soil column depth was set to 1 m with dense mesh near the upper boundary and sparse mesh near the lower boundary. The upper boundary was considered as the atmospheric boundary with runoff, while the lower one was set as free drainage. This rainfall experiment simulates a duration of 48 hours, and the initial moisture content of the soil profile is 10%.
We add the relative errors to evaluate the performance of the model. The relative errors can be expressed as , where and denote calculation results and Hydrus-1D results when the moisture migrates to the same infiltration depth.
The moisture content distribution curves for different soil types under two types of rainfall intensity conditions were obtained according to Equations (14) and (30), where δ and δ′ are taken as 0.002. In this paper, Equations (15) and (31) were adopted to obtain the wetting front depth. Furthermore, Equations (8) and (17) were also applied to further discuss the upper boundary conditions.
RESULTS AND DISCUSSION
Moisture content distribution and wetting front depth
The moisture content distribution values calculated from Equations (14) and (30) were compared with those from Hydrus-1D for sandy loam, loam, and silt. It is observed that Equations (14) and (30) perform well compared with the Hydrus-1D results, with the calculation and simulation results shown in Figure 3(a)–3(c).
The calculation results of Equation (14) agree well with those of Hydrus-1D, with only 0.5% average relative errors when the rainfall intensity is less than the saturated hydraulic conductivity. Compared with Hydrus-1D, the difference between the calculation results of Equation (30) is the smallest for sandy loam in the whole process of rainfall infiltration (when the rainfall intensity is greater than the saturated hydraulic conductivity), the largest for silt, and loam is in between. For loam and silt, the relative errors mainly increase in the stages of PIT and PIS, which is mainly due to the function of F deviating slightly from the theoretical lower limit function in the late stage of rainfall for these two soil types. However, the relative error of soil moisture content calculation in the same infiltration depth is basically controlled at 2%, and it remains reasonable to use Equations (14) and (30) to calculate the soil moisture content distribution in each stage of rainfall infiltration. Thus, it is considered that Equations (14) and (30) apply to a wide range of soil types. In addition, these results further justify the stratification of the soil moisture content after the critical time or ponding time under two types of rainfall intensity conditions.
On the one hand, various rainfall intensities were selected to calculate the wetting front depth to verify the applicability of Equations (15) and (31) for sandy loam, loam, and silt. The calculation results of Equations (15) and (31) were referred to as the predicted wetting front depth, and the Hydrus-1D results were referred to as the simulated wetting front depth, with the results shown in Figure 3(d). The data points are all close to the reference line, indicating that Equations (15) and (31) show good consistency with Hydrus-1D results for the three soils, and the relative error is basically controlled at 3%. Therefore, Equations (15) and (31) can accurately predict the wetting front depth for various types of soils.
On the other hand, the wetting front migration rate was found to decrease nonlinearly and become gradually stable with infiltration time (Figure 3(d)). The influence of wetting front depth on rainfall intensity is obviously under the condition of R < Ks. The migration rate increases with the rainfall intensity at the stage of NPIS, indicating that the infiltration depth is greater and the migration rate is faster during the same rainfall duration. Furthermore, the wetting front depth is not sensitive to rainfall intensity under the condition of R > Ks, and the migration rate is fast, slow, and uniform during the stage of NPIT, PIT, and PIS, respectively.
Surface moisture content and infiltration rate
The results of Hydrus-1D and Equation (10) correspond very well (Figure 4(a)). It is apparent that Equation (10) can effectively characterize the relationship between the volumetric moisture content on the surface and rainfall duration under R < Ks, where the surface volumetric moisture content increases with rainfall duration in a nonlinear and rapid way and finally settles down to under various rainfall intensities. In addition, with increasing rainfall intensity, the surface moisture content tends toward earlier, thus it can be concluded that the asymptote of Equation (10) (i.e., K() = R) is true. In addition, it is established that the changes of boundary conditions in the entirety of the rainfall infiltration processes satisfy the shift of constant-flux → constant-concentration (constant-flux).
Compared with the Hydrus-1D results, both Equations (10) and (17) show good performance (Figure 4(b)), therefore Equation (10) is applicable to all rainfall intensities, which characterize the change of the volumetric moisture content on the surface, and the infiltration rate curve is expressed by Equation (17) accurately. In addition, the rapid change of moisture content close to rainfall duration t is equal to 0, while near to ponding (t = tp), the change of θ0 is relatively small. Accordingly, the experimental determination of tp at θ0 = θs will be very difficult. Furthermore, the surface infiltration rate is equal to the rainfall intensity before ponding, while it decreases with rainfall duration nonlinearly and slowly, and converges to Ks after ponding. It can be seen that the changes of boundary conditions in the entirety of the rainfall infiltration processes satisfy the shifts of constant-flux → constant-concentration → constant-flux (constant-concentration). With increasing rainfall intensity, the surface moisture content reaches saturation more quickly, i.e., surface runoff is generated. Furthermore, the rainfall intensity only affects ponding time and cumulative infiltration before ponding.
Combining Figure 4(a) and 4(b) yields the relationship between critical (or ponding) time and rainfall intensity for sandy loam, with the results shown in Figure 4(c). It is observed that the Type 1 critical time t* sharply decreases with rainfall intensity, and shows a smooth transition to ponding time tp under the conditions of R > Ks. Subsequently, the curves of t*-R and tp-R can be combined into a -line, and the curves of ts-R are called the ts-line. The -line as well as the ts-line divide the entirety of the infiltration processes under all rainfall intensity conditions into two or three periods, i.e., NPIT and NPIS; or NPIT, PIT and PIS. With increasing rainfall intensity, the period of NPIT gradually declines and PIT begins to appear (until R > Ks), where the period of PIT is extended and stabilized.
Limited scope of critical point
In effect, two types of critical point are not obvious in the curves of surface moisture content and infiltrability. Critical points, however, have limited scope if we consider a reasonable relative error range, which is related to δ and δ′ values, respectively. Accordingly, the wetting front depth was studied to explore the relationship between the relative error of wetting front depth and the value of δ or δ′ for all time points after critical time t*or ts, for the soil types of sandy loam, loam, and silt.
The relative error εr can be expressed as , where denotes the wetting front depth at any time t (i.e., after the critical time) when the critical point changes with δ or δ′, and zw denotes the wetting front depth at time t when the critical point tends to the stable point (i.e., the calculation step length is 0.001). Test runs on the rainfall intensity are not sensitive to the relative error of the wetting front in the same soil, thus two types of single steady rainfall intensity are chosen to participate in the discussion of the critical point. The relationship curves between the δ or δ′ values and the relative error εr of the wetting front depth are obtained according to the calculation scheme in Table 2. The calculation results are shown in Figure 5.
|Soil texture .||Silt .||Loam .||Sandy Loam .|
|Infiltration rate (mm/h)||2||2Ks1||8||2Ks2||30||2Ks3|
|Rainfall duration (h)||30,100||ts1,30||6,100||ts2,6||3,100||ts3,3|
|Soil texture .||Silt .||Loam .||Sandy Loam .|
|Infiltration rate (mm/h)||2||2Ks1||8||2Ks2||30||2Ks3|
|Rainfall duration (h)||30,100||ts1,30||6,100||ts2,6||3,100||ts3,3|
Note: Ks1, Ks2, Ks3 denote the unsaturated hydraulic conductivity for silt, loam, and sandy loam, respectively; ts1,ts2,ts3 denote the Type 2 critical time for silt, loam, and sandy loam, which all satisfy , where δ′ is equal to 0.001 to solve the ‘stable time’ problem, i.e., ts1,2,3 = 14.957 h, 3.784 h, 1.757 h .
If the range of δ can be reasonably controlled from 0 to 0.01, the relative error limit of the calculation result of the wetting front depth is controlled as 2.5% (Figure 5(a)). This applies to the calculation of wetting front depth for the wide range of soil types in the whole process of rainfall infiltration.
Moisture migration model application
Rainfall infiltration has potential triggering factors for slope failure as rainwater into the slope encounters moisture content increase in unsaturated soil, which is of practical significance for accurately conducting slope stability analysis under rainfall infiltration conditions. The model is only applicable to the horizontal surface and not to slopes.
Comparing Equation (37) with Equation (1), only the unsaturated hydraulic conductivity is multiplied by cosα in the governing equation for slope rainfall infiltration. We can obtain the moisture migration model applied for slopes, i.e., multiplying the rainfall intensity and unsaturated hydraulic conductivity simultaneously by cosα.
The proposed model is only applicable to the case where the initial moisture content is uniform in the soil column, and the rainfall intensity is constant over a long period. However, the model can reflect the different stages of rainfall infiltration, and improve the quantification of moisture content at an arbitrary point of the soil at any time for the whole process of infiltration under steady rainfall.
This paper developed a moisture migration model that calculates the volumetric moisture content at an arbitrary point of the soil at any time, and that estimates wetting front depth. The model was derived from RE using the FCR, and the whole process of rainfall infiltration is applicable for short-duration and long-duration rainfall under two types of rainfall intensity (i.e., R > Ks and R < Ks). Furthermore, the model was compared with numerical solutions calculated by Hydrus-1D, and the comparison results showed that the moisture content distribution obtained by Equations (14) and (30) agrees well with that from Hydrus-1D for moisture migration simulation during short-duration and long-duration steady rainfall. The main conclusions can thus be drawn as follows:
The proposed moisture migration model is reasonable and applicable to a wide range of soil types under both short-duration and long-duration steady rainfall.
The boundary conditions of rainfall infiltration can consider the shift from constant-flux boundary to constant-concentration boundary during long-duration steady rainfall. The infiltration rate and surface volumetric moisture content can divide rainfall infiltration into two stages of NPIT and NPIS under the constant-flux boundary, and into three stages of PIT and PIS under the constant-concentration boundary.
The model can reflect the stratification phenomenon for soil profile wetting by infiltration, which can be simplified into three zones (i.e., transition zone, wetting zone, and dry-soil zone) from top to bottom.
The moisture migration is more sensitive to rainfall intensity under the constant-flux boundary, while is not affected by the rainfall intensity under the constant-concentration boundary, i.e., the migration rate of the wetting front positively correlates with the rainfall intensity in the stage of NPIS, while it stabilizes at Ks/(θs−θi) in the stage of PIS.
This study was supported by the National Natural Science Foundation of China (NO. 51878560), and the National Science Fund for Distinguished Young Scholars (NO. 51408491, 52008341). The authors appreciate the anonymous reviewers for their valuable comments and suggestions that helped improve the manuscript.
CONFLICT OF INTEREST
We declare that there is no conflict of interests regarding the publication of this article.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.