Soil moisture migration model based on Flux-concentration relation under steady rainfall

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 quanti ﬁ cation of moisture content at an arbitrary point of the soil at any time for the whole process of in ﬁ ltration 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 re ﬂ ect the strati ﬁ cation phenomenon for soil pro ﬁ le wetting by in ﬁ ltration. Our analysis indicates that the ﬂ ux and surface volumetric moisture content together can bound the boundary conditions of rainfall in ﬁ ltration, which presents a shift from constant- ﬂ ux to constant-concentration during a long-duration steady rainfall. The migration rate of the wetting front in the later stage of in ﬁ ltration positively correlates with rainfall intensity under the constant- ﬂ ux condition, while it ﬁ nally stabilizes at K s /( θ s (cid:1) θ i ) under the constant-concentration condition (i.e., K s -saturated hydraulic conductivity, θ s -saturated volumetric moisture content, θ i -initial volumetric moisture content). Moisture migration model was derived to improve the quanti ﬁ cation of moisture content at an arbitrary point of the soil at any time for the whole process of in ﬁ ltration under steady rainfall, which 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.

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) (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) ;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) (Kutilek 1980) used the approximate expression of FCR to derive a model for the moisture content distribution in wetting profiles under vertical rainfall infiltration, while the model is only applicable for the constant flux boundary and the short-term rainfall.
To date, the analytical solutions of RE have been derived in many studies, however, few analytical solutions were 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 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 becomes 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 the 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 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
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., Richards' equation) (Richards 1931), which applies to homogenous 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) generally difficult to solve directly, while it can be numerically solved by solvers such as Hydrus 1D with given initial and boundary conditions (Simunek 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: 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 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 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 constantflux 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, Type 1 critical point is defined (t*, u 0 ), where t* is called Type 1 critical time and u 0 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. Type 2 critical point is then defined (t s , K s ), where t s is called 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.

The Flux-concentration relation
The Flux-concentration relation, introduced by Philip (1973), satisfies the following form (Philip 1973): where 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).
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 limit 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 deltadiffusivity (Sivapalan & Milly 1989;Wang et al. 2018). Furthermore, Sivapalan et al. (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 timedependence 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: 2.2.2. Moisture migration model 2.2.2.1. R , K s . When 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: Substituting F(Θ,t) ¼ Θ 0 into Equations (8) and (9) yields The derivative of Equation (10) yields According to Equation (12), The surface moisture content tends to stabilize with the duration of rainfall, and the numerical solution may not converge when θ 0 ¼ u 0 . For convenience, the Type 1 critical point is improved as (t*, u Ã 0 ), where u Ã 0 is the critical moisture content, which denotes all moisture content in the left neighborhood of u 0 , and can be expressed in the following form: where the δ-limited range is shown in Section 4.3. Equation (10) is classified as the segmented function using 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 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: 2.2.2.2. R . K s . when R . K s , the upper boundary condition is converted between constant-flux and constantconcentration. 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., u 0 (t) ¼ u s , 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. 2002b;Brutsaert 2005;Wang et al. 2017), the relationship of infiltrability after ponding can be expressed by Equation (17) (Wang et al. 2018 where N and t 0 can be expressed as follows: The derivative of Equation (17) with respect to rainfall duration is obtained in the same manner as Equation (20), thus it can be introduced that di p =dt , 0. 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 → ∞, di p /dt ¼ 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 , K Ã s ), where t s and K Ã s can be expressed as Equation (21). K Ã s denotes all values of infiltration rate in the right-neighborhood of the saturated hydraulic conductivity.
where the δ 0 -limited range is shown in Section 4.3. F(Q, t) ¼ Q s 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 t p , t t s , thus Equation (22) can be reduced as Equation (24).
Substituting the second equation into the first equation in Equation (24) and subsequent integration yields Integrating Equation (25) with respect to t from t p to t and rearranging it yields in which z(θ, t p ) and I(t) can be expressed as follows: Rainfall infiltration is in the stage of PIS(i.e., t . t s ), where i p (t s ) ≈ K s . Equation (26) can be reduced to Equation (29).
Therefore, when the rainfall intensity is greater than saturated hydraulic conductivity, the moisture content distribution in the whole process of rainfall infiltration can be expressed by Equation (30).
The wetting front depth z w ¼ z(θ i ,t) can be obtained directly according to Equation (30). Considering K i ≈ 0, the wetting front depth can be reduced to Equation (31).
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).

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. 2001) (i.e., sandy loam, loam, and silt) were selected for model evaluation, whose van Genuchten soil hydraulic parameters are shown in Table 1. (1) was solved numerically using Hydrus 1D (Simunek 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.

The Richards' equation illustrated in Equation
For numerical simulation, the 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 1 0 r can be expressed as 1 0 r ¼ (u a À u s )=u s , where u a and u s 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 δ 0 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).
The calculation results of Equation (14) agree well with that 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 for 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 in nonlinear 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, indicting that the infiltration depth is Uncorrected Proof 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 u 0 under various rainfall intensities. In addition, with increasing rainfall intensity, the surface moisture content tends toward u 0 earlier, thus it can be concluded that the asymptote of Equation (10) (i.e., K( u 0 ) ¼ R) is true. In addition, it is established that the changes of boundary conditions in the entire rainfall infiltration processes satisfy the shift of constantflux → 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 converge to K s after ponding. it can be seen that the changes of boundary conditions in the entire 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 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 t Ã p -line, and the curves of t s -R are called as t s -line. The t Ã p -line as well as the t s -line divide the entire 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. 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 δ 0 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 δ 0 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 1 r ¼ (z Ã w À z w )=z w , where z Ã w denotes the wetting front depth at any time t (i.e., after the critical time) when the critical point changes with δ or δ 0 , 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 δ 0 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.
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 δ 0 is controlled from 0 to 0.2(i.e., the upper limit of δ 0 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 δ 0 is related to the rainfall intensity, where η can be expressed as follows: Figure 5 | The curves of ε r -δ 0 (a) and ε r -δ(b) 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%.
based on the derived model, which can be involved in the analysis of the stability for 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 the short-duration and long-duration rainfall under two types of rainfall intensities (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 the 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.