## Abstract

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 *K*_{s}/(*θ*_{s} − *θ*_{i}) under the constant-concentration condition (i.e., *K*_{s} – saturated hydraulic conductivity, *θ*_{s} – saturated volumetric moisture content, *θ*_{i} – initial volumetric moisture content).

## HIGHLIGHTS

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.

## INTRODUCTION

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.

## THEORY

### Governing equation and boundary conditions

#### Governing equation

*θ*

_{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):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 (mm

^{2}/h) at moisture content

*θ*, and is defined by: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.

*θ*(

*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: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;

*K*

_{s}represents the saturated hydraulic conductivity; and

*S*

_{e}represents effective saturation as given by

*S*

_{e}=(

*θ*−

*θ*

_{r})/(

*θ*

_{s}−

*θ*

_{r}).

#### Boundary conditions

It is generally assumed that when the rainfall intensity is less than the saturated hydraulic conductivity (*R* < *K*_{s}) 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 *t*_{p}) when the rainfall intensity is greater than the saturated hydraulic conductivity (*R* > *K*_{s}), 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 (*t*_{s},), where *t*_{s} 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, *t*_{p} and *t*_{s} 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

*i*

_{0}(

*t*) represents the flux density at the surface, i.e., infiltration rate;

*K*

_{i}denotes the unsaturated hydraulic conductivity (mm/h) at

*θ*=

*θ*

_{i}; and Θ represents relative moisture content, which can be expressed as follows: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).

*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:

#### Moisture migration model

*R* < *K*_{s}

*R*<

*K*

_{s}, the upper boundary condition is constant-flux (i.e.,

*i*

_{0}(

*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:

*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:where the

*δ*-limited range is shown in Section 4.3.

*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):where . If hydraulic conductivity at the initial soil moisture

*K*

_{i}is negligible, the wetting front depth

*z*

_{w}=

*z*(

*θ*

_{i},

*t*) can be obtained according to Equations (10) and (14), and expressed as follows:

*R* > *K*_{s}

*R*>

*K*

_{s}, the upper boundary condition is converted between constant-flux and constant-concentration. The boundary condition is the same as

*R*<

*K*

_{s}, and the moisture content distribution can be calculated by Equation (11) directly before ponding (0 <

*t*<

*t*

_{p}). The surface moisture content reaches saturation at the moment of ponding (

*t*=

*t*

_{p}), i.e., , the substitution of which into Equation (8) yields:The surface moisture content remains at the saturation state after ponding (

*t*>

*t*

_{p}), 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

*i*

_{p}. 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) :where

*N*and

*t*

_{0}can be expressed as follows:

*D*(

*θ*) is approximately a

*δ*-function, which is excessively large in a relatively small

*θ*region near

*θ*=

*θ*

_{s}(Philip 1973). This means that, as

*t*→ ∞,

*d*

*i*

_{p}/

*d*

*t*= 0,

*i*

_{p}(

*t*) =

*N*(

*θ*

_{s}) =

*K*

_{s}; the infiltrability after ponding decreases nonlinearly and tends to

*K*

_{s}with the rainfall duration. It is not possible to obtain the saturated time

*t*

_{s}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 (

*t*

_{s},), where

*t*

_{s}and can be expressed as Equation (21). denotes all values of infiltration rate in the right-neighborhood of the saturated hydraulic conductivity.where the

*δ*′-limited range is shown in Section 4.3.

*t*from

*t*

_{p}to

*t*and rearranging it yields:in which

*z*(

*θ*,

*t*

_{p}) and

*I*(

*t*) can be expressed as follows:

## VALIDATION

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
. | K_{s}/mm·h^{−1}
. | l
. |
---|---|---|---|---|---|---|

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^{−1}
. | N
. | K_{s}/mm·h^{−1}
. | l
. |
---|---|---|---|---|---|---|

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 |

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* < *K*_{s}. 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* > *K*_{s}, 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* < *K*_{s}, 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* = *t*_{p}), the change of *θ*_{0} is relatively small. Accordingly, the experimental determination of *t*_{p} 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 *K*_{s} 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 *t*_{p} under the conditions of *R* > *K*_{s}. Subsequently, the curves of *t**-*R* and *t*_{p}-*R* can be combined into a -line, and the curves of *t*_{s}-*R* are called the *t*_{s}-line. The -line as well as the *t*_{s}-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* > *K*_{s}), 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 *t*_{s}, 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 *z*_{w} 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 | 2K_{s1} | 8 | 2K_{s2} | 30 | 2K_{s3} |

Rainfall duration (h) | 30,100 | t_{s1},30 | 6,100 | t_{s2},6 | 3,100 | t_{s3},3 |

Soil texture . | Silt . | Loam . | Sandy Loam . | |||
---|---|---|---|---|---|---|

Infiltration rate (mm/h) | 2 | 2K_{s1} | 8 | 2K_{s2} | 30 | 2K_{s3} |

Rainfall duration (h) | 30,100 | t_{s1},30 | 6,100 | t_{s2},6 | 3,100 | t_{s3},3 |

*Note*: *K*_{s1}, *K*_{s2}, *K*_{s3} denote the unsaturated hydraulic conductivity for silt, loam, and sandy loam, respectively; *t*_{s1,}*t*_{s2,}*t*_{s3} 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., *t*_{s1,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.

*δ*′ is controlled from 0 to 0.2(i.e., the upper limit of

*δ*′ satisfies

*η*(

*R*/

*K*

_{s}− 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:

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

*x*′,

*z*′can be expressed as follows:

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*α*.

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

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.

## CONCLUSIONS

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* > *K*_{s} and *R* < *K*_{s}). 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

*K*_{s}/(*θ*_{s}−*θ*_{i}) in the stage of PIS.

## ACKNOWLEDGEMENTS

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.

## REFERENCES

*Infiltration Theory for Hydrologic Applications*. Water Resources Monograph 15, AGU, Washington, DC, USA