Rainfall Infiltration Modeling : A Review

Infiltration of water into soil is a key process in various fields, including hydrology, hydraulic works, agriculture, and transport of pollutants. Depending upon rainfall and soil characteristics as well as from initial and very complex boundary conditions, an exhaustive understanding of infiltration and its mathematical representation can be challenging. During the last decades, significant research effort has been expended to enhance the seminal contributions of Green, Ampt, Horton, Philip, Brutsaert, Parlange and many other scientists. This review paper retraces some important milestones that led to the definition of basic mathematical models, both at the local and field scales. Some open problems, especially those involving the vertical and horizontal inhomogeneity of the soils, are explored. Finally, rainfall infiltration modeling over surfaces with significant slopes is also discussed.


Introduction
Brutsaert [1] offers a concise definition of infiltration as "the entry of water into the soil surface and its subsequent vertical motion through the soil profile".Infiltration plays an important role in the partitioning of applied surface water into surface runoff and subsurface water-both of these components govern water supply for agriculture, transport of pollutants through the vadose zone, and recharge of aquifers [2].The infiltration of rain and surface water is influenced by many factors, including soil depth and geomorphology, soil hydraulic properties, and rainfall or climatic properties [3].Researchers have understood for ages that rainfall wets the soil and may produce runoff.Our understanding of the physics of the process and the dynamics of porous media hydraulics has come rather recently.Our ability to mathematically describe the response of a soil to rainfall and to understand the parameters that affect infiltration, has developed only in the last few decades.
The spatio-temporal evolution of infiltration rate under natural conditions cannot be currently deduced by direct measurements at all scales of interest in applied hydrology, and infiltration modeling with the aid of measurable quantities is of fundamental importance.
Even though the representation of the main natural processes in applied hydrology requires areal infiltration modeling for both flat and sloping surfaces, research activity has been limited to the development of local or point infiltration models for many years.A variety of local infiltration models for vertically homogeneous soils with constant initial soil water content and over horizontal surfaces has been proposed .A milestone paper describing the process of infiltration from a ponded surface condition was published by Reference [4].Successively, Reference [29] published his assessment of the role of infiltration in flood generation, defining "infiltration capacity" (IC) as a hyetograph separation rate that was generally applicable as a threshold for application to a rainfall intensity graph.A few years later, Reference [30] refined this concept by referring to it as an infiltration rate that declines exponentially during a storm, and then published a conceptual derivation of the exponential decay infiltration equation.On this basis, if the rainfall overcomes the IC, only a portion of it may infiltrate while the remaining quantity ponds over the surface or moves depending of the local slope.Therefore, the IC can be regarded as an important soil characteristic.As an example, Figure 1 shows several infiltration regimes when rainfall occurs over a soil surface.The line represents the IC(t) curve for a silty loam soil.The dark columns represent the rainfall rate, r(t), observed in the Umbria Region (Central Italy) during a generic frontal system.The dashed columns represent the simulated behavior of infiltration rate, f (t), during the rainfall event.With low values of rainfall intensities, the total amount of r(t) infiltrates through the soil surface and r(t) is equal to f (t).With high r(t) values only a part of the r(t) can infiltrate, while the difference r(t)−f (t) becomes runoff.As it can be observed, there exists no straightforward relationship between IC and runoff as the effective f (t) is a function of the specific rainfall intensity.
Water 2018, 10, x FOR PEER REVIEW 2 of 20 hyetograph separation rate that was generally applicable as a threshold for application to a rainfall intensity graph.A few years later, Reference [30] refined this concept by referring to it as an infiltration rate that declines exponentially during a storm, and then published a conceptual derivation of the exponential decay infiltration equation.On this basis, if the rainfall overcomes the IC, only a portion of it may infiltrate while the remaining quantity ponds over the surface or moves depending of the local slope.Therefore, the IC can be regarded as an important soil characteristic.As an example, Figure 1 shows several infiltration regimes when rainfall occurs over a soil surface.The line represents the IC(t) curve for a silty loam soil.The dark columns represent the rainfall rate, r(t), observed in the Umbria Region (Central Italy) during a generic frontal system.Natural soils are rarely vertically homogeneous.In hydrological simulations, the estimate of effective rainfall can be reasonably schematized by a two-layered vertical profile [31,32].Soils representable by a sealing layer over the parent soil or a vertical profile with a more permeable upper layer are found frequently in nature.Some models for infiltration into stable crusted soils were developed under significant approximations by adapting the Green-Ampt model [33][34][35][36][37][38] or the twostage infiltration equations proposed by Mein and Larson [39] for homogeneous soils.An efficient approach which represents transient infiltration into crusted soils was proposed by Reference [40], while Reference [41] formulated a model which describes upper-layer dynamics but in the limits of a ponded upper boundary.A more general semi-analytical/conceptual model for crusted soils was later formulated by Reference [42], and was extended by Reference [43] to represent infiltration and re-infiltration after a redistribution period under any rainfall pattern and for any two-layered soil where either layer may be more or less permeable than the other.For a much more permeable upper layer, and under more restrictive rainfall patterns, a simpler semi-empirical/conceptual model was presented by Reference [44].Under conditions of surface saturation, a simple Green-Ampt-based model was proposed [45].
In applied hydrology, upscaling of point infiltration models to the field scale is required to estimate areal-average infiltration.This is a challenging task because of the natural spatial heterogeneity of hydraulic soil properties [46][47][48][49][50][51][52], and particularly of the soil saturated hydraulic conductivity [53,54] that may be assumed as a random field with a lognormal univariate probability distribution.The mathematical problem is not analytically tractable, whereas the use of accurate Natural soils are rarely vertically homogeneous.In hydrological simulations, the estimate of effective rainfall can be reasonably schematized by a two-layered vertical profile [31,32].Soils representable by a sealing layer over the parent soil or a vertical profile with a more permeable upper layer are found frequently in nature.Some models for infiltration into stable crusted soils were developed under significant approximations by adapting the Green-Ampt model [33][34][35][36][37][38] or the two-stage infiltration equations proposed by Mein and Larson [39] for homogeneous soils.An efficient approach which represents transient infiltration into crusted soils was proposed by Reference [40], while Reference [41] formulated a model which describes upper-layer dynamics but in the limits of a ponded upper boundary.A more general semi-analytical/conceptual model for crusted soils was later formulated by Reference [42], and was extended by Reference [43] to represent infiltration and re-infiltration after a redistribution period under any rainfall pattern and for any two-layered soil where either layer may be more or less permeable than the other.For a much more permeable upper layer, and under more restrictive rainfall patterns, a simpler semi-empirical/conceptual model was presented by Reference [44].Under conditions of surface saturation, a simple Green-Ampt-based model was proposed [45].
In applied hydrology, upscaling of point infiltration models to the field scale is required to estimate areal-average infiltration.This is a challenging task because of the natural spatial heterogeneity of hydraulic soil properties [46][47][48][49][50][51][52], and particularly of the soil saturated hydraulic conductivity [53,54] Water 2018, 10, 1873 3 of 20 that may be assumed as a random field with a lognormal univariate probability distribution.The mathematical problem is not analytically tractable, whereas the use of accurate Monte Carlo (MC) simulation techniques imposes an enormous computational burden for routine applications.MC simulations were used, for instance, by References [55][56][57] to describe field-scale infiltration.MC simulations were also used by Reference [58] to corroborate a relation between areal average and variance of infiltration rate under a time-invariant rainfall rate; however, the averaging procedure was applied in space over a single realization and the behavior of the resulting errors was not specified.A variant to MC sampling for the representation of the random variability of a soil property is Latin Hypercube sampling [59] adopted by Reference [60] to develop a simple parameterized approach for areal-average infiltration.
Even though MC simulations, performed by many realizations of the random variable, are rather expensive in terms of computational effort for practical applications, they are useful as a tool to parameterize simple semi-empirical approaches or to serve as a benchmark for validating semi-analytical models.Along these lines, Reference [61] developed three versions of a semi-analytical/conceptual model for estimating the expected areal-average infiltration into vertically homogeneous soils under spatially uniform rainfall, but with random horizontal values of the saturated hydraulic conductivity (see also [62]).A shortcoming of this model is that the process of infiltration of overland flow running over pervious downstream areas (run-on process) is neglected.On the other hand, the importance of run-on was shown in a few investigations concerning the effects of horizontal variability of saturated hydraulic conductivity on Hortonian overland flow [57,[63][64][65][66][67], but this process has generally been disregarded in hydrological models.
In addition to the heterogeneity of saturated hydraulic conductivity, rainfall is also characterized by spatial variability [68,69].Some studies combining the random variability of these two quantities were performed by References [70][71][72].The latter paper presented a semi-analytical model developed under less restrictive conditions, even though run-on was not incorporated, and was based upon the use of cumulative infiltration as the independent variable that was linked with an expected time.Subsequently, Reference [73] formulated a more complete mathematical model for the expected areal-average infiltration, which considers both the saturated hydraulic conductivity and rainfall rate as random variables, and then combines the aforementioned semi-analytical approach with a semi-empirical/conceptual component to represent the run-on process.
A model of the expected areal-average infiltration into a much more permeable upper layer of a two-layered soil was also proposed by Reference [74] considering only a spatial horizontal random field of the saturated hydraulic conductivity.It involves the solution of a set of algebraic equations obtained by upscaling simple local infiltration equations to the field scale.The areal-average infiltration for a vertically non-uniform soil characterized by a saturated hydraulic conductivity decreasing with depth according to a power law was derived by Reference [75] and upscaled to the field scale by the same semi-analytical technique used in previous papers [61,72].
Finally, when macropore flow plays a significant role in determining infiltration, amendments of the Darcy-Richards approach may be used, especially at the local scale [98,99].
The main intent of this review paper is to critically assess the complexities of the infiltration process and to provide guidance for developments related to open problems.This paper will provide classical approaches developed for rainfall events with continuous saturation at the soil surface, and a more general formulation suitable for any type of rainfall pattern for applications at the local (point) scale in homogeneous soils.Simple models for point infiltration into a two-layered soil with a more permeable upper layer and a more complex model for any two-layered soil type are presented.Semi-empirical and semi-analytical field-scale infiltration models are analyzed.Finally, problems linked to rainfall infiltration into surfaces with significant slopes are also discussed.

Basic Physical Models for Infiltration
By considering a horizontally homogeneous soil, water movement in the vertical direction is governed by one-dimensional soil water flow and continuity equations.The flow rate, q, per unit cross-sectional area is described by Darcy's law, actually proposed by Reference [100], as: where K is the hydraulic conductivity, ψ is the soil water matric capillary head, and z is the vertical soil depth assumed positive downward.The infiltration rate, q 0 , is given by Equation ( 1) applied at the soil surface.
In the absence of changes in the water density and soil porosity as well as of sinks and sources, the continuity equation is: where θ is volumetric water content and t is time.The substitution of Equation ( 1) into Equation ( 2) leads to the well-known Richards equation: where C 1 = dθ/dψ and C 2 = dK/dψ with a typical assumption that θ and K are unique functions of ψ thus neglecting hysteresis in these functions.The initial condition at time t = 0 for z ≥ 0 is ψ = ψ i , and the upper boundary conditions at the soil surface, where z = 0, are: where r is the rainfall rate, t p is the time to ponding, and t r is the duration of rainfall.Hereafter the subscripts i and s denote initial and saturation quantities, respectively, while 0 stands for quantities at the soil surface.The lower boundary conditions at a depth z b which is not reached by the wetting front is ψ(z b ) = ψ i for t > 0. The soil water hydraulic properties can be represented by the following parameterized forms [22]: where θ * s and K * s are used as scaling quantities, ψ b is the air entry head, θ r is the residual volumetric water content, c, λ, and d are empirical coefficients, and b = 3 and a = 2 according to Burdine's method [101].For particular values of the parameters, Equations (5a) and (5b) reduce to the well-known equations proposed by References [101,102].For two-layered soils, two additional conditions are required at the interface between the two layers: Water 2018, 10, 1873 5 of 20 where hereafter the subscripts 1, 2, and c denote variables in the upper layer, in the lower layer, and at the interface, respectively, and Z c is the interface depth.
In particular cases, it can be useful to express Equation (3) in terms of θ obtaining: where D(θ) is the soil water diffusivity equal to K(θ)(dψ(θ)/dθ).

Point Infiltration Modeling for Homogeneous Soils
Many local infiltration models for vertically homogeneous soils with constant initial soil water content and over horizontal surfaces are widely recognized in the scientific literature [4,5,[8][9][10][11]17,18,20,22,23,26,[28][29][30]103].Furthermore, for isolated storms and when ponding is not achieved instantly, extended forms of the Philip model [45], the Green-Ampt model [104,105], and the Smith and Parlange model [106] have been widely used, whereas for arbitrary rainfall patterns, the model presented in [26] serves as a useful method.These four models, together with a widely adopted equation suggested by References [29,30], have been extensively used in applied hydrology and as building blocks in the development of infiltration approaches at the field scale.

Horton Empirical Equation
In the empirical equation proposed by References [29,30], the infiltration capacity, f c , exponentially decreases as follows (see also Figure 2): where f 0 and f f represent the initial and final values of f c , respectively, and α is the decay constant.
When t→∞, f f can be considered equal to the saturated hydraulic conductivity of the soil.If K and D are independent of θ, References [107,108] demonstrated that Equation ( 8) can be obtained from Equation (7).
Water 2018, 10, x FOR PEER REVIEW 5 of 20

Point Infiltration Modeling for Homogeneous Soils
Many local infiltration models for vertically homogeneous soils with constant initial soil water content and over horizontal surfaces are widely recognized in the scientific literature [4,5,[8][9][10][11]17,18,20,22,23,26,[28][29][30]103].Furthermore, for isolated storms and when ponding is not achieved instantly, extended forms of the Philip model [45], the Green-Ampt model [104,105], and the Smith and Parlange model [106] have been widely used, whereas for arbitrary rainfall patterns, the model presented in [26] serves as a useful method.These four models, together with a widely adopted equation suggested by References [29,30], have been extensively used in applied hydrology and as building blocks in the development of infiltration approaches at the field scale.

Horton Empirical Equation
In the empirical equation proposed by References [29,30], the infiltration capacity, fc, exponentially decreases as follows (see also Figure 2): where f0 and ff represent the initial and final values of fc, respectively, and α is the decay constant.When t→∞, ff can be considered equal to the saturated hydraulic conductivity of the soil.If K and D are independent of θ, References [107,108] demonstrated that Equation ( 8) can be obtained from Equation (7).

Philip Equation
A widely adopted analytical solution of the Richards equation was proposed by References [9][10][11]109] under the conditions of vertically homogeneous soil, constant initial moisture content, and saturated soil surface with immediate ponding.For early to intermediate times the semi-analytical solution for the infiltration capacity, based on a series expansion and truncation after the first two terms, is expressed as: where S is the sorptivity, depending on soil properties and initial moisture content, and A is a quantity ranging from 0.38 Ks to 0.66 Ks.For t→∞, Equation ( 9) is replaced by fc = Ks.The integration of Equation ( 9) yields the cumulative infiltration:

Philip Equation
A widely adopted analytical solution of the Richards equation was proposed by References [9][10][11]109] under the conditions of vertically homogeneous soil, constant initial moisture content, and saturated soil surface with immediate ponding.For early to intermediate times the semi-analytical solution for the infiltration capacity, based on a series expansion and truncation after the first two terms, is expressed as: where S is the sorptivity, depending on soil properties and initial moisture content, and A is a quantity ranging from 0.38 K s to 0.66 K s .For t→∞, Equation ( 9) is replaced by f c = K s .The integration of Equation ( 9) yields the cumulative infiltration: Philip's model was extended for applications to less restrictive conditions.For a constant rainfall rate r > K s , surface saturation occurs at a time t p > 0 and, following [45], infiltration can be described through an equivalent time origin, t 0 , for potential infiltration after ponding as: For unsteady rainfall under the condition of a continuously saturated surface for t > t p , the infiltration process can be represented adopting a similar procedure.Generally, S and A are derived from the calibration of hydrological models; however, S can also be approximated as [110]: where φ is the soil porosity and ψ av is the soil water matric capillary head at the wetting front.

Green-Ampt Model
This model represents infiltration into homogeneous soils under the conditions of continuously saturated soil surface and uniform initial soil moisture as: where F the cumulative depth of infiltrated water.To express the infiltration as a function of time, this equation can be solved after the substitution f c = dF/dt.The resulting equation [45] is: Equation ( 16) assumes immediate ponding and an infinite supply of water at the surface.Under more general conditions, with a constant rainfall rate r > K s , that begins at the time t = 0, surface saturation is reached at a time t p > 0. For t ≤ t p, the infiltration rate q 0 is equal to r and later to the infiltration capacity.Mein and Larson [104] formulated this process through Equation ( 15) as: Water 2018, 10, 1873 7 of 20 which leads to the determination of t p as: Furthermore, for t > t p, Equation ( 16) becomes: Equation ( 19) can be solved at each time, for example, by successive substitutions of F which is then substituted into Equation ( 15) to obtain the corresponding value of f c .

Parlange-Lisle-Braddock-Smith Model
A 3-parameter model obtained through an analytical integration of the Richards equation has been formulated by Reference [106] and can be expressed as: where F = F − K i t is the cumulative dynamic infiltration rate, α is a parameter linked to the behavior of hydraulic conductivity and soil water diffusivity as functions of θ, and G is the integral capillary drive defined by: Equation ( 20) includes as limiting forms [60] the Reference [17] equation (α = 1) and the Green-Ampt equation (α→0).It can be applied to determine t p and f c for any rainfall pattern, and for t > t p can be rewritten under the condition of surface saturation in a time dependent form as: where K d = K s − K i and F p = F (t p ).The quantities F p and t p are the values of F and t at ponding, respectively, at which Equation (20) with f c = r(t p ) is first satisfied.The value of α usually ranges from 0.8 to 0.85 [3].

Corradini-Melone-Smith Semi-Analytical/Conceptual Model
When the hypothesis of uniform initial soil moisture cannot be guaranteed, the models presented in the previous subsections cannot be applied.For complex rainfall patterns involving rainfall hiatus periods or a rainfall rate after time to ponding less than soil infiltration capacity, an approach for the application of the aforementioned classical models was developed [111][112][113] starting from the time compression approximation proposed by Reference [114] for post-hiatus rainfall producing immediate ponding (see also [1]).However, by comparing results with the Richards equation, Reference [22] showed that this approach was not sufficiently accurate because it neglects the soil water redistribution process (Figure 3) which is particularly important when long periods with a light rainfall or a rainfall hiatus occur.
approach for the application of the aforementioned classical models was developed [111][112][113] starting from the time compression approximation proposed by Reference [114] for post-hiatus rainfall producing immediate ponding (see also [1]).However, by comparing results with the Richards equation, Reference [22] showed that this approach was not sufficiently accurate because it neglects the soil water redistribution process (Figure 3) which is particularly important when long periods with a light rainfall or a rainfall hiatus occur.Models combining infiltration and redistribution are therefore the best solution when soils are subjected to complex rainfall patterns, which are prevalent under natural conditions (see also [115,116]).Dagan and Bresler [20] developed an analytical model along this line starting from depth integrated forms of Darcy's law and the continuity equation and using simplifications in the initial and surface boundary conditions that make easier areal investigations but reduce practical applications at the local scale.In any case, the model is not applicable to local studies with erratic rainfalls producing successive infiltration-redistribution cycles.A more general model was formulated by Reference [26] starting from the same integrated equations, and then combined with a conceptual representation of the wetting soil moisture profile.
To demonstrate the structure of the model presented in [26], a specific rainfall pattern to describe all the involved components is presented here.The application to erratic rainfall is then a straightforward extension.Let us consider a stepwise rainfall pattern involving successive periods of rainfall with constant r > Ks, separated by periods with r = 0 (see Figure 4).We denote by t1 the duration of the first pulse, t2 and t3, the beginning and end of the second pulse, respectively, and t4 the beginning of the third pulse.The model was derived considering a soil with a constant value of z θ i θ 0 θ Models combining infiltration and redistribution are therefore the best solution when soils are subjected to complex rainfall patterns, which are prevalent under natural conditions (see also [115,116]).Dagan and Bresler [20] developed an analytical model along this line starting from depth integrated forms of Darcy's law and the continuity equation and using simplifications in the initial and surface boundary conditions that make easier areal investigations but reduce practical applications at the local scale.In any case, the model is not applicable to local studies with erratic rainfalls producing successive infiltration-redistribution cycles.A more general model was formulated by Reference [26] starting from the same integrated equations, and then combined with a conceptual representation of the wetting soil moisture profile.
To demonstrate the structure of the model presented in [26], a specific rainfall pattern to describe all the involved components is presented here.The application to erratic rainfall is then a straightforward extension.Let us consider a stepwise rainfall pattern involving successive periods of rainfall with constant r > K s , separated by periods with r = 0 (see Figure 4).We denote by t 1 the duration of the first pulse, t 2 and t 3 , the beginning and end of the second pulse, respectively, and t 4 the beginning of the third pulse.The model was derived considering a soil with a constant value of θ i and combining the depth-integrated forms of Darcy's law and the continuity equation.In addition, as the event progresses in time, a dynamic wetting profile, of lowest depth Z and represented by a distorted rectangle through a shape factor β(θ 0 ) ≤ 1, was assumed.The resulting ordinary differential equation is: where p is a parameter linked with the profile shape of θ and G(θ i ,θ 0 ) is expressed by Equation ( 21) modified by the substitution of θ s with θ 0 and K s with K 0 .Equation ( 23) can be applied for 0 < t < t 2 , and the profile shape of θ(z) is approximated [23] by: Functional forms for β and p were obtained by calibration using results provided by the Richards equation applied to a silty loam soil, specifically: Equation ( 23) can be solved numerically.For q0 = r, with F′ = (r-Ki)t, it gives θ0(t) until time to ponding, t′p, corresponding to θ0 = θs and dθ0/dt = 0.Then, for t′p < t ≤ t1, with θ0 = θs and dθ0/dt = 0, it provides the infiltration capacity (q0 = fc) and for t1 < t < t2, with q0 = 0, it gives dθ0/dt<0 thus describing the redistribution process.Equation ( 23) can be solved numerically.For q 0 = r, with F = (r-K i )t, it gives θ 0 (t) until time to ponding, t p , corresponding to θ 0 = θ s and dθ 0 /dt = 0.Then, for t p < t ≤ t 1 , with θ 0 = θ s and dθ 0 /dt = 0, it provides the infiltration capacity (q 0 = f c ) and for t 1 < t < t 2 , with q 0 = 0, it gives dθ 0 /dt < 0 thus describing the redistribution process.
The second rainfall pulse leads to a new time to ponding, t" p , but reinfiltration occurs according to two alternative approaches determined by a comparison of r and the downward redistribution rate, D F (t = t 2 ), expressed by: More specifically, for r ≤ D F , the reinfiltrated water is distributed to the whole dynamic profile and θ 0 (t) can be still computed by Equation ( 23), whereas for r > D F , the profile of θ(z, t = t 2 ) is assumed temporarily invariant and starts a superimposed secondary wetting profile which advances alongside the pre-existing profile according to Equation (23) modified by substituting θ i with θ 0 (t 2 ) and F with F 2t accumulated for t ≥ t 2 .If the secondary profile reaches the depth of the first steady one, the compound profile reduces to a single profile and then Equation ( 23) can be again applied (see Figure 3 for a schematic representation of θ(z,t)).On the other hand, if at t = t 3 the secondary profile has not caught up with the first one, redistribution is first applied to the secondary profile and then re-established to the single profile in the successive rainfall hiatus.Finally, in the case at t = t 4 , the θ(z) profile is still compound and r is larger than D F (t 4 ), and a procedure of consolidation that merges the composite profile is applied early to avoid the formation of a further additional profile.
This model incorporates all the components required for application to any natural rainfall pattern.It was calibrated by Reference [26] for a silty loam soil, and then tested using different soils from clay loam to sandy loam soil types.Weighted implicit finite difference solutions of the Richards equation were used as a benchmark.For each soil, the model accuracy was found to be acceptable in terms of both infiltration rate and soil water content, even though better results were obtained for fine-textured soils.In addition, the model was found to simulate fairly well the θ(z,t) profiles observed in laboratory (Figure 5) and field experiments [117,118].
the θ(z) profile is still compound and r is larger than DF(t4), and a procedure of consolidation that merges the composite profile is applied early to avoid the formation of a further additional profile.
This model incorporates all the components required for application to any natural rainfall pattern.It was calibrated by Reference [26] for a silty loam soil, and then tested using different soils from clay loam to sandy loam soil types.Weighted implicit finite difference solutions of the Richards equation were used as a benchmark.For each soil, the model accuracy was found to be acceptable in terms of both infiltration rate and soil water content, even though better results were obtained for fine-textured soils.In addition, the model was found to simulate fairly well the θ(z,t) profiles observed in laboratory (Figure 5) and field experiments [117,118].

Point Infiltration Modeling for Vertically Non-Uniform Soils
A two-layer approximation, with each layer being schematized as homogeneous, is frequently used to set up models of infiltration for natural soils.The process of formation of a sealing layer was accurately examined by References [31,119], and that of disruption was considered by References [120][121][122].Evidence of the role of crusted soils in semi-arid regions has been recently provided by Reference [123].Green-Ampt-based models for infiltration into stable crusted soils were proposed by References [33,34,36,39].An efficient approach that represents transient infiltration into crusted soils was proposed by Reference [40], whereas Reference [41] formulated a model that describes upper-layer dynamics for a ponded upper boundary.On the other hand, vertical profiles with a more permeable upper layer are observed in hydrological practice and can be also used, for example, as a first approximation in the representation of infiltration into homogeneous soils with grassy vegetation [124].In the latter layering type, the simple model presented by Reference [45] can be usefully applied for infiltration into a saturated soil surface.Under more general conditions, the model by Reference [43] appears to be accurate with modest computational effort.

Green-Ampt-Based Model for a Layered Soil
A model for infiltration into a two-layered soil with a more permeable upper layer under the condition of continuously saturated soil surface has been formulated by Reference [45].The classical Green-Ampt equation is applied until the wetting front is in the upper layer, then the following equations are used: Water 2018, 10, 1873 where L 2 is the depth of the wetting front below the interface.Equation ( 27) can be solved at each time for L 2 by successive substitutions; L 2 is then used in Equations ( 28) and ( 29) to determine F and f c , respectively.
4.2.Corradini-Melone-Smith Semi-Analytical/Conceptual Model for a Two-Layered Soil A semi-analytical/conceptual model applicable to any horizontal two-layered soil where either layer may be more permeable has been formulated by Reference [43].It relies on the same elements previously used by Reference [26] but adopted here in each layer and integrated at the interface between the two layers by the boundary conditions expressing continuity of flow rate and capillary head (Equations (6a) and (6b)).In addition, at the lower boundary for t > 0, we have ψ 2 = ψ i .The initial condition is ψ 1 = ψ 2 = ψ i constant at t = 0, and at the interface q(Z c ) is approximated through the downward flux in the upper layer as: where G(ψ c ,ψ 10 ) is expressed by Equation ( 21) modified by the substitutions of D(θ)dθ with K(ψ)dψ, θ s with ψ 10 , and θ i with ψ c .As long as water does not infiltrate in the lower layer, the model presented in [26] is used, then starting from the time t c when the wetting front enters the lower layer, the following system of two ordinary differential equations may be applied: with P L (ψ c , t) and F 2 defined as: and C 1 (ψ 10 ) = dθ 10 /dψ 10 , C 1 (ψ c ) = dθ 1c /dψ 1c .The quantity γ represents a conceptual proportion of the upper layer where θ is increasing due to rainfall and is assumed equal to 0.85 [42], β 2 and p 2 are determined by Equations (25a-c) but applied using θ 20 , θ 2s , θ 2i , and θ 2r and substituting r with q(Z c ).
On the basis of the same stepwise rainfall pattern earlier adopted to explain the model presented in [26], Equations ( 31) and ( 32) may be used for t c < t < t 2 .Then, the two-layer model has to be applied in each layer by analogy with the procedure described for homogeneous soils; in particular, compound and consolidated profiles develop in each layer.In the underlying soil, the generation of additional profiles occurs through q(Z c ) in substitution of r.On the basis of the described steps, model application to arbitrary rainfall patterns is straightforward.The solution of the above system, Equations ( 31) and (32), may be obtained by a library routine for the Runge-Kutta-Verner fifth-order method with a variable time step.Calibration and testing of the model were performed through a comparison with numerical solutions of the Richards equation.Three soils (clay loam, silty loam, and sandy loam) with a variety of thicknesses were combined to realize two-layered soils, where either layer was more permeable, that were selected as test cases.In all instances, the simulations involved the cycle of infiltration-redistribution-reinfiltration.The infiltration rate as well as the water content at the surface and interface were found to be very accurately estimated by the semi-analytical/conceptual model.

Areal Infiltration Models over Soil with Variable Hydraulic Properties
At the field scale, considering the significant spatial variability of the main soil hydraulic properties, rainfall infiltration modeling is not analytically tractable, whereas the use of accurate Monte Carlo (M C) simulation techniques imposes an enormous computational burden for routine applications.
In the Introduction section, the evolution of studies for field-scale infiltration models for spatial variability in soil hydraulic properties and rainfall was presented.Models for infiltration at the field scale have been recently developed and represent a useful support for practical hydrological purposes.Two models characterized by significant differences in complexity and application area are presented here.

Smith and Goodrich Approach
A semi-empirical model to determine the areal-average infiltration rate into areas with random spatial variability of K s has been proposed by Reference [60].The authors assumed a lognormal probability density function (PDF) of K s with a mean value <K s > and a coefficient of variation CV(K s ), and considered one realization of the random variable.Then, adopting the model presented in [106] (see also [22]) and the Latin Hypercube sampling method, and through a large number of simulations performed for many values of CV(K s ) and rainfall rates, they developed the following effective areal relation for the scaled areal-average infiltration rate, I * e , linked with the corresponding scaled cumulative depth, F * e : with c a 1 + 0.8 [CV(K s )] 1. 3 1 − e (0.85(r * b −1)) , (36) where I * e = I e /K e , F * e = F e /[G(θ s − θ i )], r * e = r/K e and, r * b = r/ < K s > .The quantity K e denotes the areal effective value of K s given by: where f K s (K) is the PDF of K s .Finally, Equation ( 35) may be also applied for r variable with time.

Govindaraju-Corradini-Morbidelli Semi-Analytical/Conceptual Model
Govindaraju et al. [72] formulated a semi-analytical model to estimate the expected field-scale infiltration rate <I n (F)> under the condition of negligible effects of the run-on process.The model incorporates heterogeneity of both K s and r assumed as random variables with a lognormal and a uniform PDF, respectively.The quantity <I n (F)> is estimated through the averaging procedure over the ensemble of two-dimensional realizations of K s and r.For the sake of simplicity, we first examine the model under a steady-rainfall condition, and then provide the guidelines for applications to a time-varying rainfall rate.
Starting from the extended Green-Ampt model and choosing F as the independent variable, <I n (F)> at a given F can be written as: where f r (r) and f K s (K) are the PDFs of r and K s , respectively, with K c which denotes the maximum value of K s leading to surface saturation in the i-th cell, determined by: Equation ( 38) may be expressed as: with r min and r min + R extreme values of the PDF of r and M K s given by: where K a and ω stand for the first and the second argument, respectively, of the M K s function.To relate time to F, an implicit relation between the expected value of t, <t>, and F is provided as: To extend the model by incorporating the run-on effect, an additional empirical term is used in the form presented in [73]: where t pa is the time to ponding (see Equation ( 18)) associated with <r> and <K s >.The parameters a, b, and c are expressed by: a where Equation ( 44) holds for θ i θ s and for θ i →θ s we have a→0.Furthermore, Equation ( 46) is undefined for CV(r) and/or CV(K s ) equal to 0. In Equations ( 43)-( 46), length units are in mm and time is expressed in hours.
When the spatial heterogeneity of r can be neglected, with spatially uniform and steady rainfall and negligible run-on, Equations ( 40) and ( 42) can be replaced by [61]: The last formulation was also extended for applications involving r variable in time.The same method can be used to adapt Equations ( 40) and ( 42) for unsteady rainfall patterns.Furthermore, the additional empirical term for run-on could be adapted for unsteady rainfalls following the guidelines indicated by Reference [73].
The solution of the model even in the conditions of coupled spatial variability of r and K s is fairly simple and requires limited computational effort.
The model for coupled heterogeneity of r and K s (Equations ( 40), (42), and ( 43)) was validated by comparison with the results derived starting from MC sampling and using a combination of the extended Green-Ampt formulation at the local scale with the kinematic wave approximation [125] that is required to represent run-on.Through a wide variety of simulations it was shown that: (1) the model produced very accurate estimates of <I> over a clay loam soil and a sandy loam soil; (2) the spatial heterogeneity of both r and K s can be neglected only when <r> <K> or for storm durations much greater than t pa ; (3) the effects on <I> produced by significant values of CV(K s ) and CV(r) are similar; (4) run-on plays a significant role for moderate storms and high values of CV(K s ) and CV(r); and (5) the model can be simplified using Equations ( 47) and (48) for CV(r) substantially less than CV(K s ) and steady rainfalls.

Conclusions and Open Problems
In spite of the continuous developments in infiltration modeling, challenges regarding infiltration exist for many scales of hydrologic interest.The conflicting results from different studies on the effect of slope on infiltration over homogeneous surfaces show that a solution continues to elude researchers.Careful experimental and theoretical investigations in this regard are needed to fully comprehend this fundamental infiltration behavior over sloping surfaces.
The estimate of infiltration at different spatial scales (i.e., from the local to watershed scales) is a complex problem as further challenges are imposed by the natural spatial variability of soil hydraulic characteristics and that of rainfall.An important issue to be addressed when areal estimates are involved is that concerning the determination of <K s >, CV(K s ), <r>, and CV(r) together with the corresponding quantities for soil moisture content [126].
The models presented here apply to infiltration into a soil matrix.When macropore flow is significant, the problem becomes much more complicated even though simplified approaches have been proposed.For example, two practical approximations to describe infiltration into soils with macropores rely upon the use of modified values of the saturated hydraulic conductivity of the soil matrix [127] or upon the representation of the two processes of infiltration controlled by the matrix potential and the macropore volume [128].However, there exist many other models to represent the macropore flow [129], based on the single continuum approach (e.g., [130]), the dual continuum approach, the dual permeability approach (e.g., [131]), and the dual porosity approach (e.g., [132]).Notwithstanding the high number of simplified models, and given the difficulty to investigate the Navier-Stokes equations over finite soil portions, macropore flow still needs a convincing physical theory for the scales of practical interest.
Finally, all these models are formulated for horizontal land surfaces.Extensions of the classical infiltration theory to inclined surfaces were proposed by References [81,87]; however, these theories do not explain the results of laboratory experiments, for example those performed on bare soils in the absence of erosion and a sealing layer [88,92].The modeling of the slope effects has to be therefore considered as an open problem [97], in particular when surfaces with vegetation are involved [95].
Further challenges exist in our ability to independently measure soil properties at the local scale to identify the true nature of field-scale variability.The use of common measuring instruments for soil hydraulic properties does not yield consistent estimates of variability (e.g., [133]).Both the inter-and intra-instrument errors contribute to a level of uncertainty that has not been understood.Even though measurement techniques are not the topic of this review, measurement errors and uncertainty nevertheless influence modeling efforts that rely on these data.Efforts to appropriately combine disparate measurement results are also needed to realize the full worth of the data that are generated from expensive and time-consuming experimental campaigns.
Field-scale experiments measuring runoff and deep flow are influenced by spatial variability as well, but as yet no theory exists for elucidating the underlying spatial variability from these experiments.Current approaches (e.g., [134]) rely on brute force calibration techniques; however, such methods are often plagued by identifiability and non-uniqueness problems.Moreover, calibration is known to compensate for various model deficiencies, and authors deriving parameter estimates from these approaches must be cognizant of the role played by the underlying model structure.Better conceptual and theoretical underpinnings are needed to move the science of infiltration forward.

Figure 1 .
Figure1.Infiltration capacity, IC(t), and infiltration rate, f (t), of a silty loam soil under a natural rainfall event (observed in Central Italy) characterized by a variable intensity, r(t).

Figure 2 .
Figure 2. Graphical representation of the Horton empirical equation.

Figure 2 .
Figure 2. Graphical representation of the Horton empirical equation.

Figure 3 .
Figure 3. Graphical representation of the soil water redistribution process.

Figure 3 .
Figure 3. Graphical representation of the soil water redistribution process.

)Figure 4 .
Figure 4. (a) Rainfall pattern selected to describe the Corradini-Melone-Smith semi-analytical model for point infiltration.(b,c) Profiles of soil water content at various times indicated in (a) and associated with different infiltration-redistribution stages.For symbols, see the text. rainfall

Figure 4 .
Figure 4. (a) Rainfall pattern selected to describe the Corradini-Melone-Smith semi-analytical model for point infiltration.(b,c) Profiles of soil water content at various times indicated in (a) and associated with different infiltration-redistribution stages.For symbols, see the text.

Figure 5 .
Figure 5. Laboratory experimental system adopted by References [117,118] to verify the Corradini-Melone-Smith model.

Figure 5 .
Figure 5. Laboratory experimental system adopted by References [117,118] to verify the Corradini-Melone-Smith model.
The dashed columns represent the simulated behavior of infiltration rate, f(t), during the rainfall event.With low values of rainfall intensities, the total amount of r(t) infiltrates through the soil surface and r(t) is equal to f(t).With high r(t) values only a part of the r(t) can infiltrate, while the difference r(t)−f(t) becomes runoff.As it can be observed, there exists no straightforward relationship between IC and runoff as the effective f(t) is a function of the specific rainfall intensity.
Figure1.Infiltration capacity, IC(t), and infiltration rate, f(t), of a silty loam soil under a natural rainfall event (observed in Central Italy) characterized by a variable intensity, r(t).