**1. 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 [4–28]. 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 *<sup>r</sup>*(*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.

**Figure 1.** 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*).

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–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–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 Monte Carlo (MC) simulation techniques imposes an enormous computational burden for routine applications. MC simulations were used, for instance, by References [55–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–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–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].

Most of the above-mentioned models assumed zero or small surface slope that does not affect the infiltration process. However, in most real situations, infiltration occurs over surfaces characterized by different gradients [76,77] and the role of surface slope on infiltration is not clear. In fact, the results obtained by some theoretical and experimental investigations [78–97] lead to contrasting conclusions, suggesting that an improved understanding and modeling of infiltration on sloping surfaces is required.

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.

#### **2. 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:

$$q = -K\left(\frac{\partial\Psi}{\partial z} - 1\right),\tag{1}$$

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:

$$
\frac{\partial \theta}{\partial t} = \frac{\partial q}{\partial z} \; \; \; \tag{2}
$$

where *θ* is volumetric water content and *t* is time. The substitution of Equation (1) into Equation (2) leads to the well-known Richards equation:

$$\mathbf{C}\_{1}\frac{\partial\psi}{\partial t} = \frac{\partial}{\partial z}\left[\mathbf{K}\frac{\partial\psi}{\partial z}\right] - \mathbf{C}\_{2}\frac{\partial\psi}{\partial z'}\tag{3}$$

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 *ψ* = *ψ<sup>i</sup>*, and the upper boundary conditions at the soil surface, where *z* = 0, are:

$$q\_0 = r, \; 0 < t \le t\_P; \; \theta\_0 = \theta\_{5\prime}, t\_P < t \le t\_{\bar{r}}; \; q\_0 = 0, \; t\_{\bar{r}} < t\_{\prime} \tag{4}$$

where *r* is the rainfall rate, *tp* is the time to ponding, and *tr* 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 *zb* which is not reached by the wetting front is *ψ*(*zb*) = *ψi* for *t* > 0. The soil water hydraulic properties can be represented by the following parameterized forms [22]:

$$
\psi = \psi\_b \left[ \left( \frac{\theta - \theta\_r}{\theta\_s^\* - \theta\_r} \right)^{-c/\lambda} - 1 \right]^{1/c} + d \,, \tag{5a}
$$

$$K = K\_s^\bullet \left[ 1 + \left( \frac{\Psi - d}{\Psi\_b} \right)^c \right]^{-(b\lambda + a)/c},\tag{5b}$$

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:

$$
\psi\_1(Z\_\mathfrak{c}) = \psi\_2(Z\_\mathfrak{c}) = \psi\_\mathfrak{c} \, , \tag{6a}
$$

$$K\_1 \left[ \left( \frac{\partial \psi\_1}{\partial z} \right)\_{Z\_c} - 1 \right] = K\_2 \left[ \left( \frac{\partial \psi\_2}{\partial z} \right)\_{Z\_c} - 1 \right],\tag{6b}$$

where hereafter the subscripts 1, 2, and *c* denote variables in the upper layer, in the lower layer, and at the interface, respectively, and *Zc* is the interface depth.

In particular cases, it can be useful to express Equation (3) in terms of *θ* obtaining:

$$\frac{\partial\theta}{\partial t} = \frac{\partial}{\partial z}\left[D(\theta)\frac{\partial\theta}{\partial z} + K(\theta)\right],\tag{7}$$

where *D*(*θ*) is the soil water diffusivity equal to *K*(*θ*)(*dψ*(*θ*)/*dθ*).

#### **3. 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–11,17,18, 20,22,23,26,28–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.

#### *3.1. Horton Empirical Equation*

In the empirical equation proposed by References [29,30], the infiltration capacity, *fc*, exponentially decreases as follows (see also Figure 2):

$$f\_{\mathbb{C}} = f\_f + \left(f\_0 - f\_f\right) \exp(-\text{at})\,,\tag{8}$$

where *f* 0 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).

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