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.

flux-concentration relation, moisture content distribution, rainfall infiltration, wetting front depth

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

Governing equation

This study assumes a uniform distribution of initial soil moisture content θi, constant and continuous rainfall intensity R, and disregards the effects of groundwater. The governing equation for one-dimensional vertical infiltration in unsaturated soils is subsequently given by the following equation (i.e., the Richards equation) (Richards 1931), which applies to homogeneous and isotropic soil (Figure 1):
formula
(1)
where z represents the vertical coordinate, positive upward; θ represents the volumetric moisture content of a soil; t represents time (h); K(θ) represents unsaturated hydraulic conductivity (mm/h), which is a function of θ; and D(θ) represents soil water diffusivity (mm2/h) at moisture content θ, and is defined by:
formula
(2)
In general, Equation (1) is difficult to solve directly, while it can be numerically solved by solvers such as Hydrus-1D with given initial and boundary conditions (Šimůnek et al. 2005). The soil water retention curve θ(h) and unsaturated hydraulic conductivity function K(θ) are required for the solution.
The van Genuchten equation is the most widely used function for θ(h), while the modified Mualem–van Genuchten formulation for K(θ) is adopted in this study (Mualem 1976; van Genuchten 1980). This can be expressed in the following form:
formula
(3)
formula
(4)
where θr, θs represent residual moisture content and saturated moisture content, respectively; α, m, n are empirical parameters, and the parameters m, n satisfy the relationship given by m = 1 − 1/n, n > 1; l represents the pore-connectivity parameter and is evaluated to be about 0.5 as an average for many soils; Ks represents the saturated hydraulic conductivity; and Se represents effective saturation as given by Se=(θ−θr)/(θs−θr).

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.

Figure 1
Schematic representation of rainfall infiltration and definition of the coordinate system: x, z are Cartesian rectangular space coordinates (mm) with x parallel to the surface direction, and z normal to the surface; L denotes the soil thickness (m).
View largeDownload slide

Schematic representation of rainfall infiltration and definition of the coordinate system: x, z are Cartesian rectangular space coordinates (mm) with x parallel to the surface direction, and z normal to the surface; L denotes the soil thickness (m).

Figure 1
Schematic representation of rainfall infiltration and definition of the coordinate system: x, z are Cartesian rectangular space coordinates (mm) with x parallel to the surface direction, and z normal to the surface; L denotes the soil thickness (m).
View largeDownload slide

Schematic representation of rainfall infiltration and definition of the coordinate system: x, z are Cartesian rectangular space coordinates (mm) with x parallel to the surface direction, and z normal to the surface; L denotes the soil thickness (m).

Close modal
Figure 2
Boundary conditions of infiltration behavior under rainfall, where represents the constant-flux boundary, and represents the constant-concentration boundary.
View largeDownload slide

Boundary conditions of infiltration behavior under rainfall, where represents the constant-flux boundary, and represents the constant-concentration boundary.

Figure 2
Boundary conditions of infiltration behavior under rainfall, where represents the constant-flux boundary, and represents the constant-concentration boundary.
View largeDownload slide

Boundary conditions of infiltration behavior under rainfall, where represents the constant-flux boundary, and represents the constant-concentration boundary.

Close modal

Moisture migration model of soil profile

The flux-concentration relation

The flux-concentration relation, introduced by Philip (1973), satisfies the following form:
formula
(5)
where i0(t) represents the flux density at the surface, i.e., infiltration rate; Ki denotes the unsaturated hydraulic conductivity (mm/h) at θ = θi; and Θ represents relative moisture content, which can be expressed as follows:
formula
(6)
where θ0 is the soil surface moisture content, which is time-dependent; after surface saturation, infiltration near the surface will be equal to the infiltrability (mm/h).
Philip has proposed that if the exact expression of F(Θ,t) is known for any given process, the relevant exact solution of Equation (1) can be found immediately. In general, F(Θ,t) is not normally exactly known, but is obtained by using a guessed F(Θ,t) function. In addition, the theoretical upper and lower limits of F(Θ,t) proposed by Philip are widely accepted, which correspond to two types of soils, namely, ‘delta-function’ soil and ‘linear’ soil. For natural soils, however, soil water diffusivity D(θ) is usually close to a delta function, i.e., the ‘delta-function’ soil, which means that F(Θ,t) = Θ is exact for soils with delta-diffusivity (Sivapalan & Milly 1989; Wang et al. 2018). Furthermore, Sivapalan & Milly (1989) stated that, for a variety of flow processes, the FCR is nearly the same regardless of the actual boundary conditions. For constant-flux absorption, the time-dependence of FCR appears to be important at short times only (White et al. 1979). In this paper, F(Θ,t) = Θ is used for the whole process of rainfall, and for the convenience of expression, F(Θ,t) is expressed as follows:
formula
(7)

Moisture migration model

R < Ks
When R < Ks, the upper boundary condition is constant-flux (i.e., i0(t) = R). The expressions for the soil surface moisture content θ0 and the moisture content profiles z(θ,t) are derived based on FCR (Kutilek 1980), and can be formulated as follows:
formula
(8)
formula
(9)
Substituting F(Θ,t) = Θ0 into Equations (8) and (9) yields:
formula
(10)
formula
(11)
The derivative of Equation (10) yields:
formula
(12)
According to Equation (12), dθ0/dt > 0. In the case of dθ0/dt = 0, K(⁠⁠) = R. The surface moisture content tends to stabilize with the duration of rainfall, and the numerical solution may not converge when θ0 = ⁠. For convenience, the Type 1 critical point is improved as (t*,⁠), where is the critical moisture content, which denotes all moisture content in the left-neighborhood of ⁠, and can be expressed in the following form:
formula
(13)
where the δ-limited range is shown in Section 4.3.
Equation (10) is classified as a segmented function using the Type 1 critical point, including a nonlinear function and a constant function. Equation (11) is independent of the rainfall duration after the surface moisture content reaches the critical moisture content, which is not applicable to calculate the soil moisture content distribution after the critical duration. The established rainfall infiltration tests show that the wetting front depth increases linearly with rainfall duration in the late period of infiltration (Qin et al. 2017). Accordingly, the moisture distribution of the soil profile occurs in two zones after the critical time, i.e., wetting zone and transition zone, where the thickness of the wetting zone tends to be stable and the moisture content of the transition zone is equal to the critical moisture content. Therefore, Equation (9) can be advanced into two stages of NPIT and NPIS, as shown in Equation (14):
formula
(14)
where ⁠. If hydraulic conductivity at the initial soil moisture Ki is negligible, the wetting front depth zw = z(θi,t) can be obtained according to Equations (10) and (14), and expressed as follows:
formula
(15)
R > Ks
When R > Ks, the upper boundary condition is converted between constant-flux and constant-concentration. The boundary condition is the same as R < Ks, and the moisture content distribution can be calculated by Equation (11) directly before ponding (0 < t < tp). The surface moisture content reaches saturation at the moment of ponding (t = tp), i.e., ⁠, the substitution of which into Equation (8) yields:
formula
(16)
The surface moisture content remains at the saturation state after ponding (t > tp), i.e., past the constant-concentration boundary, and runoff is produced. Furthermore, the infiltration rate gradually moves away from the rainfall intensity, which is called the soil infiltrability ip. For vertical rainfall infiltration conditions, based on the FCR and the time compression approximation (TCA), which uses shift time and ponding time to calculate the post-ponding infiltration (Smith et al. 2002; Brutsaert 2005; Wang et al. 2017), the relationship of infiltrability after ponding can be expressed by Equation (17) (Wang et al. 2018) (Equation (19), in which γ = 0) :
formula
(17)
where N and t0 can be expressed as follows:
formula
(18)
formula
(19)
The derivative of Equation (17) with respect to rainfall duration is obtained in the same manner as Equation (20), and thus it can be introduced that ⁠. In addition, D(θ) is approximately a δ-function, which is excessively large in a relatively small θ region near θ = θs (Philip 1973). This means that, as t → ∞, dip/dt = 0, ip(t) = N(θs) = Ks; the infiltrability after ponding decreases nonlinearly and tends to Ks with the rainfall duration. It is not possible to obtain the saturated time ts using the inverse solution of Equation (17), where the infiltration rate is equal to the saturated hydraulic conductivity. Therefore, the Type 2 critical point is improved as (ts,⁠), where ts and can be expressed as Equation (21). denotes all values of infiltration rate in the right-neighborhood of the saturated hydraulic conductivity.
formula
(20)
formula
(21)
where the δ′-limited range is shown in Section 4.3.
and the continuity equation are combined to obtain the moisture content distribution after ponding, i.e., Equation (22). Rainfall infiltration is in the stage of PIT, and the initial conditions and upper boundary conditions are shown in Equation (23) in the period of ⁠, thus Equation (22) can be reduced as Equation (24):
formula
(22)
formula
(23)
formula
(24)
Substituting the second equation into the first equation in Equation (24) and subsequent integration yields:
formula
(25)
Integrating Equation (25) with respect to t from tp to t and rearranging it yields:
formula
(26)
in which z(θ, tp) and I(t) can be expressed as follows:
formula
(27)
formula
(28)
Rainfall infiltration is in the stage of PIS (i.e., ⁠), where ip(ts) ≈ Ks. Equation (26) can be reduced to Equation (29):
formula
(29)
Therefore, when the rainfall intensity is greater than the saturated hydraulic conductivity, the moisture content distribution in the whole process of rainfall infiltration can be expressed by Equation (30):
formula
(30)
The wetting front depth zw = z(θi,t) can be obtained directly according to Equation (30). Considering Ki ≈ 0, the wetting front depth can be reduced to Equation (31):
formula
(31)
formula
(32)
where ⁠, and and indicate that the results are obtained by replacing R in Equations (16) and (19) with ⁠, i.e., ⁠.

The migration depth of the wetting front is related to the soil type, initial moisture content, and rainfall intensity according to Equations (15) and (31).

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.

Table 1

Average hydraulic parameters of the three main soil texture groups

Soil textureθrθsα/mm−1NKs/mm·h−1l
Sandy loam 0.065 0.41 0.0075 1.89 44.21 0.5 
Loam 0.078 0.43 0.0036 1.56 10.40 0.5 
Silt 0.034 0.46 0.0016 1.37 2.50 0.5 
Soil textureθrθsα/mm−1NKs/mm·h−1l
Sandy loam 0.065 0.41 0.0075 1.89 44.21 0.5 
Loam 0.078 0.43 0.0036 1.56 10.40 0.5 
Silt 0.034 0.46 0.0016 1.37 2.50 0.5 
View Large

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.

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).

Figure 3
Comparisons of moisture content distribution (a–c) and wetting front depth (d) simulated by the proposed moisture migration model for three typical soils (sandy loam, loam, and silt).
View largeDownload slide

Comparisons of moisture content distribution (a–c) and wetting front depth (d) simulated by the proposed moisture migration model for three typical soils (sandy loam, loam, and silt).

Figure 3
Comparisons of moisture content distribution (a–c) and wetting front depth (d) simulated by the proposed moisture migration model for three typical soils (sandy loam, loam, and silt).
View largeDownload slide

Comparisons of moisture content distribution (a–c) and wetting front depth (d) simulated by the proposed moisture migration model for three typical soils (sandy loam, loam, and silt).

Close modal

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).

Figure 4
The surface moisture content and infiltration rate curve predicted by Hydrus-1D (solid lines), Equations (10) and (17) (scatter-dot) for three types of soil under different rainfall intensities, where (a) R < Ks, (b) R > Ks. (c) Relationship between critical (ponding) time and rainfall intensity.
View largeDownload slide

The surface moisture content and infiltration rate curve predicted by Hydrus-1D (solid lines), Equations (10) and (17) (scatter-dot) for three types of soil under different rainfall intensities, where (a) R < Ks, (b) R > Ks. (c) Relationship between critical (ponding) time and rainfall intensity.

Figure 4
The surface moisture content and infiltration rate curve predicted by Hydrus-1D (solid lines), Equations (10) and (17) (scatter-dot) for three types of soil under different rainfall intensities, where (a) R < Ks, (b) R > Ks. (c) Relationship between critical (ponding) time and rainfall intensity.
View largeDownload slide

The surface moisture content and infiltration rate curve predicted by Hydrus-1D (solid lines), Equations (10) and (17) (scatter-dot) for three types of soil under different rainfall intensities, where (a) R < Ks, (b) R > Ks. (c) Relationship between critical (ponding) time and rainfall intensity.

Close modal

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.

Table 2

Calculation scheme design

Soil textureSiltLoamSandy 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 textureSiltLoamSandy 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 ⁠.

View Large
Figure 5
The curves of (a) εr-δ′ and (b) εr-δ for three types of soil (silt, loam, and sandy loam), where the relative error limit of the wetting front depth is controlled as 2.5%.
View largeDownload slide

The curves of (a) εr-δ′ and (b) εr-δ for three types of soil (silt, loam, and sandy loam), where the relative error limit of the wetting front depth is controlled as 2.5%.

Figure 5
The curves of (a) εr-δ′ and (b) εr-δ for three types of soil (silt, loam, and sandy loam), where the relative error limit of the wetting front depth is controlled as 2.5%.
View largeDownload slide

The curves of (a) εr-δ′ and (b) εr-δ for three types of soil (silt, loam, and sandy loam), where the relative error limit of the wetting front depth is controlled as 2.5%.

Close modal

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.

When the range of δ′ is controlled from 0 to 0.2(i.e., the upper limit of δ′ satisfies η(R/Ks − 1) = 0.2, where 0 < η ≤ 1), the relative error limit of wetting front depth is similarly controlled as 2.5% (Figure 5(b)). The upper limit of δ′ is related to the rainfall intensity, where η can be expressed as follows:
formula
(33)

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.

We add the simplification conditions to extend the application scope of the model, i.e., assume that the slope is infinitely long and the moisture content contour in the soil is parallel to the slope. The slope profile is shown in Figure 6, where x′, z′can be expressed as follows:
formula
(34)
Figure 6
The schematic representation of slope rainfall infiltration and the definition of the coordinate system, where α is the slope inclination angle; x, z are Cartesian rectangular space coordinates (mm); and x′, z′ are slope surface space (rotated) (mm) with x′ taken parallel to the slope surface direction and z′ taken normal to the slope surface.
View largeDownload slide

The schematic representation of slope rainfall infiltration and the definition of the coordinate system, where α is the slope inclination angle; x, z are Cartesian rectangular space coordinates (mm); and x′, z′ are slope surface space (rotated) (mm) with x′ taken parallel to the slope surface direction and z′ taken normal to the slope surface.

Figure 6
The schematic representation of slope rainfall infiltration and the definition of the coordinate system, where α is the slope inclination angle; x, z are Cartesian rectangular space coordinates (mm); and x′, z′ are slope surface space (rotated) (mm) with x′ taken parallel to the slope surface direction and z′ taken normal to the slope surface.
View largeDownload slide

The schematic representation of slope rainfall infiltration and the definition of the coordinate system, where α is the slope inclination angle; x, z are Cartesian rectangular space coordinates (mm); and x′, z′ are slope surface space (rotated) (mm) with x′ taken parallel to the slope surface direction and z′ taken normal to the slope surface.

Close modal
The governing equation is the two-dimensional Richards equation for a homogeneous soil as given by the following equation:
formula
(35)
Substituting Equation (34) into Equation (35) yields:
formula
(36)
It can be seen that the moisture content in the soil is only related to the rainfall duration t and z′ according to the assumption, and then Equation (36) can be reduced to Equation (37):
formula
(37)

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 moisture migration model for slopes is designed to reflect the moisture content distribution in the soil by rainwater into the soil. Thus the matrix suction at an arbitrary point of the soil can be derived from the van Genuchten equation, i.e., Equation (38). Generally, the matrix suction has a great influence on the shear strength of unsaturated soils. Fredlund et al. (1978) proposed that matrix suction in unsaturated soils can increase cohesion, and its influence on shear strength is reflected by Equation (39). It can be seen that the shear strength at an arbitrary point in the slope can be further obtained based on the derived model, which can be involved in the analysis of the stability of the slope under the rainfall infiltration.
formula
(38)
formula
(39)

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:

  • (1)

    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.

  • (2)

    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.

  • (3)

    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.

  • (4)

    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.

We declare that there is no conflict of interests regarding the publication of this article.

All relevant data are included in the paper or its Supplementary Information.

Brutsaert
W.
2005
Hydrology: An Introduction
.
Cambridge University Press
,
Cambridge, UK
.
Google Scholar
Crossref
Search ADS
 
Cary
C. E.
 &
Zapata
C. E.
2011
Resilient modulus for unsaturated unbound materials
.
Road Materials and Pavement Design
12
(
3
),
615
–
638
.
Google Scholar
Crossref
Search ADS
 
Chen
H. K.
,
Tang
H. M.
,
Tang
Y. H.
,
Tan
L.
,
Wei
L.
,
He
X. Y.
 &
Gao
Y. H.
2014
Experimental study of infiltration characteristics of strongly weathered mudstone under action of heavy rainfall
.
Rock and Soil Mechanics
35
(
10
),
2755
–
2760
.
Google Scholar
 
Chiu
A. C. F.
,
Zhu
W.
 &
Chen
X.
2009
Rainfall infiltration pattern in an unsaturated silty sand
.
Journal of Hydrologic Engineering
14
(
8
),
882
–
886
.
Google Scholar
Crossref
Search ADS
 
Colman
E. A.
 &
Bodman
G. B.
1945
Moisture and energy conditions during downward entry of water into moist and layered soils
.
Soil Science Society of America Journal
9
(
C
),
3
–
11
.
Google Scholar
Crossref
Search ADS
 
Fang
W.
 &
Esaki
T.
2012
Rapid assessment of regional superficial landslide under heavy rainfall
.
Journal of Central South University
19
(
9
),
2663
–
2673
.
Google Scholar
Crossref
Search ADS
 
Fredlund
D. G.
,
Morgenstern
N. R.
 &
Widger
R. A.
1978
The shear strength of unsaturated soils
.
Canadian Geotechnical Journal
15
(
3
),
313
–
321
.
Google Scholar
Crossref
Search ADS
 
Jacobs
J. M.
,
Mohanty
B. P.
,
Hsu
E. C.
 &
Miller
D.
2004
SMEX02: field scale variability, time stability and similarity of soil moisture
.
Remote Sensing of Environment
92
(
4
),
436
–
446
.
Google Scholar
Crossref
Search ADS
 
Kutilek
M.
1980
Constant-rainfall infiltration
.
Journal of Hydrology
45
(
3–4
),
289
–
303
.
Google Scholar
 
Lee
M. L.
,
Gofar
N.
 &
Rahardjo
H.
2009
A simple model for preliminary evaluation of rainfall-induced slope instability
.
Engineering Geology
108
,
272
–
285
.
Google Scholar
Crossref
Search ADS
 
Li
A. G.
,
Yue
Z. Q.
,
Tham
L. G.
,
Lee
C. F.
 &
Law
K. T.
2005
Field-monitored variations of soil moisture and matric suction in a saprolite slope
.
Canadian Geotechnical Journal
42
(
1
),
13
–
26
.
Google Scholar
Crossref
Search ADS
 
Lytton
R.
,
Aubeny
C.
 &
Bulut
R.
2005
Design Procedure for Pavements on Expansive Soils
.
Texas Transportation Institute, Texas A&M University System
,
College Station, TX, USA
.
Google Scholar
 
Parlange
J.-Y.
,
Lisle
I.
,
Braddock
R. D.
 &
Smith
R. E.
1982
The three-parameter infiltration equation
.
Soil Science
133
(
6
),
337
–
341
.
Google Scholar
Crossref
Search ADS
 
Philip
J. R.
1969
Theory of infiltration
.
Advances in Hydroscience
5
,
215
–
296
.
Google Scholar
Crossref
Search ADS
 
Qin
X. H.
,
Liu
D. S.
,
Song
Q. H.
,
Du
C. L.
 &
Wang
X.
2017
Experimental study on one-dimensional vertical infiltration in soil column under rainfall and the derivation of permeability coefficient
.
Chinese Journal of Rock Mechanics and Engineering
36
(
2
),
475
–
484
.
Google Scholar
 
Šimůnek
J.
,
Van Genuchten
M. T.
 &
Šejna
M.
2005
The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media
,
Version 3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA
.
Sivapalan
M.
 &
Milly
P. C. D.
1989
On the relationship between the time condensation approximation and the flux concentration relation
.
Journal of Hydrology
105
(
3–4
),
357
–
367
.
Google Scholar
 
Smith
R. E.
,
Smettem
K. R. J.
,
Broadbridge
P.
 &
Woolhiser
D. A.
2002
Infiltration Theory for Hydrologic Applications. Water Resources Monograph 15, AGU, Washington, DC, USA
.
van Genuchten
M. T.
1980
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils
.
Soil Science Society of America Journal
44
(
5
),
892
–
898
.
Google Scholar
Crossref
Search ADS
 
Vogel
T.
,
van Genuchten
M. T.
 &
Cislerova
M.
2000
Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions
.
Advances in Water Resources
24
(
2
),
133
–
144
.
Google Scholar
Crossref
Search ADS
 
Wang
K.
,
Yang
X.
,
Liu
X.
 &
Liu
C.
2017
A simple analytical infiltration model for short-duration rainfall
.
Journal of Hydrology
555
,
141
–
154
.
Google Scholar
Crossref
Search ADS
 
White
I.
,
Smiles
D. E.
 &
Perroux
K. M.
1979
Absorption of water by soil: the constant flux boundary condition
.
Soil Science Society of America Journal
43
(
4
),
659
–
664
.
Google Scholar
Crossref
Search ADS
 
Wu
L. Z.
,
Huang
R. Q.
,
Xu
Q.
,
Zhang
L. M.
 &
Li
H. L.
2015
Analysis of physical testing of rainfall-induced soil slope failures
.
Environmental Earth Sciences
73
,
8519
–
8531
.
Google Scholar
Crossref
Search ADS
 
Zeng
L.
,
Li
G. Y.
,
Shi
Z. N.
,
Liu
D. S.
,
Liu
J.
 &
Li
D. K.
2018
Experiment on seepage characteristics of unsaturated soil under rainfall infiltration
.
China Journal of Highway and Transport
31
(
2
),
191
–
199
.
Google Scholar
 
Zhan
T. L. T.
,
Jia
G. W.
,
Chen
Y.-M.
,
Fredlund
D. G.
 &
Li
H.
2013
An analytical solution for rainfall infiltration into an unsaturated infinite slope and its application to slope stability analysis
.
International Journal for Numerical and Analytical Methods in Geomechanics
37
(
12
),
1737
–
1760
.
Google Scholar
Crossref
Search ADS
 
Zhu
W.
,
Chen
X. D.
 &
Zhong
X. C.
2006
Observation and analysis of rainfall infiltration
.
Rock and Soil Mechanics
27
(
11
),
1873
–
1879
.
Google Scholar
 
© 2022 The Authors
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).