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

$$f\_c = K\_s \left[ 1 + \frac{a}{\exp\left(\frac{aF'}{G(\theta\_s - \theta\_i)}\right) - 1} \right],\tag{20}$$

where *F* = *F* − *Kit* 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:

$$G = \frac{1}{\mathcal{K}\_s} \int\_{\theta\_i}^{\theta\_s} D(\theta) d\theta \tag{21}$$

Equation (20) includes as limiting forms [60] the Reference [17] equation (*α* = 1) and the Green–Ampt equation (*α*<sup>→</sup>0). It can be applied to determine *tp* and *fc* for any rainfall pattern, and for *t* > *tp* can be rewritten under the condition of surface saturation in a time dependent form as:

$$\ln\left[(1-a)K\_{\rm s} - K\_{\rm i}\right](t - t\_{\rm P}) = F^{\prime} - F\_{\rm p}^{\prime} - \frac{(\theta\_{\rm s} - \theta\_{\rm i})K\_{\rm s}G}{K\_{\rm d}} \ln\left\{\frac{\exp\left[\frac{aF^{\prime}}{(\theta\_{\rm s} - \theta\_{\rm i})G}\right] - 1 + \frac{aK\_{\rm e}}{K\_{\rm d}}}{\exp\left[\frac{aF\_{\rm p}^{\prime}}{(\theta\_{\rm s} - \theta\_{\rm i})G}\right] - 1 + \frac{aK\_{\rm e}}{K\_{\rm d}}}\right\},\tag{22}$$

where *Kd* = *Ks* − *Ki* and *<sup>F</sup>p* = *F*(*tp*). The quantities *<sup>F</sup>p* and *tp* are the values of *F* and *t* at ponding, respectively, at which Equation (20) with *fc* = *r*(*tp*) is first satisfied. The value of *α* usually ranges from 0.8 to 0.85 [3].

#### *3.5. 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–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.

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

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

$$\frac{d\theta\_0}{dt} = \frac{(\theta\_0 - \theta\_i)\beta(\theta\_0)}{F'\left[ (\theta\_0 - \theta\_i)\frac{d\beta(\theta\_0)}{d\theta\_0} + \beta(\theta\_0) \right]} \left[ q\_0 - K\_0 - \frac{(\theta\_0 - \theta\_i)G(\theta\_i, \theta\_0)\beta(\theta\_0)pK\_0}{F'} \right],\tag{23}$$

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 *Ks* with *K*0. Equation (23) can be applied for 0 < *t* < *t*2, and the profile shape of *θ*(*z*) is approximated [23] by:

$$\frac{\theta(z) - \theta\_i}{\theta\_0 - \theta\_i} = 1 - \exp\left[\frac{\beta z(\theta\_0 - \theta\_i) - F'}{(\beta - \beta^2) - F'}\right].\tag{24}$$

Functional forms for *β* and *p* were obtained by calibration using results provided by the Richards equation applied to a silty loam soil, specifically:

$$\beta(\theta\_0) = 0.6 \frac{\theta\_s - \theta\_i}{\theta\_s - \theta\_r} + 0.4,\tag{25a}$$

$$
\beta p = 0.98 - 0.87 \exp\left(-\frac{r}{K\_s}\right) \qquad \frac{d\theta\_0}{dt} \ge 0,\tag{25b}
$$

$$
\beta p = 1.7 \qquad \frac{d\theta\_0}{dt} < 0 \,. \tag{25c}
$$

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

Equation (23) can be solved numerically. For *q*0 = *r*, with *F* = (*r-Ki*)*<sup>t</sup>*, it gives *θ*0(*t*) until time to ponding, *<sup>t</sup>p*, corresponding to *θ*0 = *θs* and *dθ*0/*dt* = 0. Then, for *tp* < *t* ≤ *t*1, with *θ*0 = *θs* and *dθ*0/*dt* = 0, it provides the infiltration capacity (*q*0 = *fc*) 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, *<sup>t</sup>"p*, but reinfiltration occurs according to two alternative approaches determined by a comparison of *r* and the downward redistribution rate, *DF*(*t* = *t*2), expressed by:

$$D\_F = \left(-\frac{1}{\beta} \frac{d\beta}{d\theta\_0} - \frac{1}{\theta\_0 - \theta\_i}\right) F' \frac{d\theta\_0}{dt} - K\_i \quad \text{t} = \ t\_2. \tag{26}$$

More specifically, for *r* ≤ *DF*, the reinfiltrated water is distributed to the whole dynamic profile and *θ*0(*t*) can be still computed by Equation (23), whereas for *r* > *DF*, the profile of *θ*(*<sup>z</sup>*, *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 *<sup>F</sup>*2*t* 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 *θ*(*<sup>z</sup>*,*<sup>t</sup>*)). 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 *DF*(*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 *θ*(*<sup>z</sup>*,*<sup>t</sup>*) profiles observed in laboratory (Figure 5) and field experiments [117,118].

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

#### **4. 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–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.

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

$$L\_{2}\frac{(\theta\_{2s}-\theta\_{2i})}{K\_{2s}} + \frac{\left[ (\theta\_{2s}-\theta\_{2i})Z\_{c}K\_{2s} - (\theta\_{2s}-\theta\_{2i})K\_{1s}(\psi\_{av\ 2}+Z\_{c}) \right]}{K\_{1s}K\_{2s}}\ln\left[1+\frac{L\_{2}}{\psi\_{av\ 2}+Z\_{c}}\right] = t\,,\tag{27}$$

$$F = Z\_c(\theta\_{1s} - \theta\_{1i}) + L\_2(\theta\_{2s} - \theta\_{2i}) \, , \tag{28}$$

$$f\_c = \frac{K\_{1s}K\_{2s}}{Z\_4K\_{2s} + L\_2K\_{1s}}(\psi\_{av\ 2} + Z\_c + L\_2)\tag{29}$$

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 *fc*, 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 = *ψ<sup>i</sup>*. The initial condition is *ψ*1 = *ψ*2 = *ψi* constant at *t =* 0, and at the interface *q*(*Zc*) is approximated through the downward flux in the upper layer as:

$$q(Z\_{\varepsilon}) = K\_{1s} \frac{G\_1(\psi\_{\varepsilon}, \psi\_{10})}{Z\_{\varepsilon}} + K\_{1\varepsilon} \tag{30}$$

where *G*(*ψ<sup>c</sup>*,*ψ*10) is expressed by Equation (21) modified by the substitutions of *D*(*θ*)*dθ* with *<sup>K</sup>*(*ψ*)*dψ*, *θs* with *ψ*10, and *θi* with *ψ<sup>c</sup>*. As long as water does not infiltrate in the lower layer, the model presented in [26] is used, then starting from the time *tc* when the wetting front enters the lower layer, the following system of two ordinary differential equations may be applied:

$$\frac{d\psi\_{10}}{dt} = \frac{1}{\gamma Z\_{\text{c}} \mathbb{C}\_{1}(\psi\_{10})} \left[ q\_{10} - K\_{1\varepsilon} - \frac{K\_{1\varepsilon} G\_{1}(\psi\_{\varepsilon}, \psi\_{10})}{Z\_{\text{c}}} \right] - \frac{(1 - \gamma) \mathbb{C}\_{1}(\psi\_{\varepsilon})}{\gamma \mathbb{C}\_{1}(\psi\_{10})} \frac{d\psi\_{\varepsilon}}{dt} t \ge t\_{\varepsilon} \,, \tag{31}$$

$$\frac{d\psi\_{\rm c}}{dt} = \frac{1}{P\_{\rm L}(\psi\_{\rm c}, t)} \left[ K\_{\rm lc} + \frac{K\_{\rm lc} G\_{\rm l}(\psi\_{\rm c}, \psi\_{\rm l0})}{Z\_{\rm c}} - K\_{\rm 2c} - \frac{\beta\_{2}(\theta\_{2c}) p\_{2}(\theta\_{2c} - \theta\_{2i}) K\_{2s} G\_{2}(\psi\_{\rm i}, \psi\_{\rm c})}{F\_{2}} \right] t \geq t\_{\rm c}, \tag{32}$$

with *PL*(*ψ<sup>c</sup>*, *t*) and *F*2 defined as:

$$P\_L(\psi\_c, t) = \left[\beta\_2(\theta\_{2c}) + \frac{d\beta\_2}{d\theta\_{2c}}(\theta\_{2c} - \theta\_{2i})\right] \frac{F\_2}{(\theta\_{2c} - \theta\_{2i})\beta\_2(\theta\_{2c})}C\_2(\psi\_c) \,, \tag{33}$$

$$F\_2 = F - Z\_{\varepsilon}[\gamma(\theta\_{1s} - \theta\_{1i}) + (1 - \gamma)(\theta\_{1c} - \theta\_{1i})] - K\_{2i}t\_{\prime} \tag{34}$$

and *C*1(*ψ*10) = *dθ*10/*dψ*10, *C*1(*ψc*) = *dθ*1*c*/*dψ*<sup>1</sup>*c*. 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, *θ*2*s*, *θ*2*i*, and *θ*2*r* and substituting *r* with *q*(*Zc*). 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 *tc* < *t<t2*. 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*(*Zc*) 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.

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

#### *5.1. Smith and Goodrich Approach*

A semi-empirical model to determine the areal-average infiltration rate into areas with random spatial variability of *Ks* has been proposed by Reference [60]. The authors assumed a lognormal probability density function (PDF) of *Ks* with a mean value < *Ks*> and a coefficient of variation *CV*(*Ks*), 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*(*Ks*) 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* :

$$I\_{\varepsilon}^{\*} = 1 + (r\_{\varepsilon}^{\*} - 1) \left\{ 1 + \left[ \frac{r\_{\varepsilon}^{\*} - 1}{a} \left( \varepsilon^{a F\_{\varepsilon}^{\*}} - 1 \right) \right]^{\varepsilon\_{d}} \right\}^{-1/\varepsilon\_{d}} \qquad r\_{\varepsilon}^{\*} > 1,\tag{35}$$

with

$$\mathcal{L}\_a \cong 1 + \frac{0.8}{\left[\mathcal{C}V(\mathcal{K}\_s)\right]^{1.3}} \left[1 - e^{\left(0.85(r\_b^\*-1)\right)}\right],\tag{36}$$

where *I*∗ *e* = *Ie*/*Ke*, *F*∗ *e* = *Fe*/[*G*(*<sup>θ</sup>s* − *<sup>θ</sup>i*)], *r*∗ *e* = *r*/*Ke* and, *r*∗ *b* = *r*/ < *Ks* > . The quantity *Ke* denotes the areal effective value of *Ks* given by:

$$K\_{\mathfrak{k}} = \int\_0^r K f \chi\_{\mathfrak{k}}(K) dK + [1 - \int\_0^r f \chi\_{\mathfrak{k}}(K) dK] r \,\,\,\,\,\tag{37}$$

where *fKs*(*K*) is the PDF of *Ks*. Finally, Equation (35) may be also applied for *r* variable with time.

#### *5.2. Govindaraju–Corradini–Morbidelli Semi-Analytical/Conceptual Model*

Govindaraju et al. [72] formulated a semi-analytical model to estimate the expected field-scale infiltration rate <*In*(*F*)> under the condition of negligible effects of the run-on process. The model incorporates heterogeneity of both *Ks* and *r* assumed as random variables with a lognormal and a uniform PDF, respectively. The quantity <*In*(*F*)> is estimated through the averaging procedure over the ensemble of two-dimensional realizations of *Ks* 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:

$$\langle I(F)\rangle \,\, = \int\_0^\infty \int\_{K\_c}^\infty r f\_r(r) f\_{K\_s}(K) dr \, dK + \int\_0^\infty \int\_0^{K\_c} \left(1 + \frac{\psi(\theta\_s - \theta\_i)}{F}\right) K f\_r(r) f\_{K\_s}(K) dr \, dK \,\, , \tag{38}$$

where *fr*(*r*) and *fKs* (*K*) are the PDFs of *r* and *Ks*, respectively, with *Kc* which denotes the maximum value of *Ks* leading to surface saturation in the *i*-th cell, determined by:

$$K\_{\mathbb{C}} = \frac{Fr\_i}{\psi(\theta\_{\mathbb{S}} - \theta\_i) + F} = F\_{\mathbb{C}}r\_i \,. \tag{39}$$

Equation (38) may be expressed as:

$$\begin{aligned} \langle I\_n(F) \rangle &= \frac{1}{2K\_c^2} \{ M\_{K\_s}[(r\_{\min} + R)F\_{c'}, 2] - M\_{K\_s}[(r\_{\min})F\_{c'}, 2] \} \\ &- \frac{r\_{\min}^2}{2\mathbb{R}} \{ M\_{K\_s}[(r\_{\min} + R)F\_{c'}, 0] - M\_{K\_s}[(r\_{\min})F\_{c'}, 0] \} \\ &+ \left( r\_{\min} + \frac{R}{2} \right) \{ 1 - M\_{K\_s}[(r\_{\min} + R)F\_{c'}, 0] \} \\ &+ \frac{1}{F\_c} \left( \frac{r\_{\min} + R}{R} \right) \{ M\_{K\_s}[(r\_{\min} + R)F\_{c'}, 1] - M\_{K\_s}[(r\_{\min})F\_{c'}, 1] \} \\ &- \frac{1}{F\_c} \left( \frac{1}{RT\_c} \right) \{ M\_{K\_s} \left[ (r\_{\min} + R)F\_{c'}, 2 \right] - M\_{K\_s}[(r\_{\min})F\_{c'}, 2] \} \\ &+ \frac{1}{F\_c} \{ M\_{K\_s} \left[ (r\_{\min})F\_{c'}, 1 \right] \} \end{aligned} \tag{40}$$

with *rmin* and *rmin* + *R* extreme values of the PDF of *r* and *MKs*given by:

$$M\_{\mathcal{K}\_\bullet}(\mathcal{K}\_a, \omega) = \int\_0^{\mathcal{K}\_a} \mathcal{K}^\omega f\_\mathcal{K}(\mathcal{K}) d\mathcal{K} \,\,\,\,\tag{41}$$

where *Ka* and *ω* stand for the first and the second argument, respectively, of the *MKs* function. To relate time to *F*, an implicit relation between the expected value of *t*, <*t*>, and *F* is provided as:

$$\begin{split} \left< t(F) \right> &= \frac{F}{\left< r \right>} \left\{ 1 - M\_{\mathcal{K}\_{\boldsymbol{\theta}}} \left[ \left< r \right> F\_{\boldsymbol{\varepsilon}}, 0 \right] \right\} + \left[ F + \psi(\theta\_{\boldsymbol{s}} - \theta\_{\boldsymbol{i}}) \ln \left( \frac{\psi(\theta\_{\boldsymbol{s}} - \theta\_{\boldsymbol{i}})}{\overline{\psi(\theta\_{\boldsymbol{s}} - \theta\_{\boldsymbol{i}}) + F}} \right) \right] \\ & \left\{ M\_{\mathcal{K}\_{\boldsymbol{\varepsilon}}} \left[ \left< r \right> F\_{\boldsymbol{\varepsilon}}, -1 \right] \right\} + \psi(\theta\_{\boldsymbol{s}} - \theta\_{\boldsymbol{i}}) \sum\_{j=1}^{\infty} \frac{1}{(j+1) \left< r \right>^{j+1}} \left\{ M\_{\mathcal{K}\_{\boldsymbol{\varepsilon}}} \left[ \left< r \right> F\_{\boldsymbol{\varepsilon}}, j \right] \right\} \ . \end{split} \tag{42}$$

To extend the model by incorporating the run-on effect, an additional empirical term is used in the form presented in [73]:

$$<\langle I(t)\rangle \stackrel{\approx}{\equiv} \langle I\_n(t)\rangle \, + \,\,\langle r\rangle \, a\left(\frac{t}{t\_{pa}}\right)^b \exp\left(-c\frac{t}{t\_{pa}}\right) \,. \tag{43}$$

where *tpa* is the time to ponding (see Equation (18)) associated with <*r*> and <*Ks*>. The parameters *a*, *b*, and *c* are expressed by:

$$a = 2.8[\mathbb{C}V(r) + \mathbb{C}V(K\_{\mathbb{S}})]^{0.36},\tag{44}$$

$$b = 5.35 - 6.32[\mathcal{CV}(r)\mathcal{CV}(\mathbb{K}\_s)]\_\prime \tag{45}$$

$$\varepsilon = 2.7 + 0.3 \left[ \frac{\langle r \rangle / \langle \mathbf{K\_s} \rangle}{\langle \mathbf{C} V(r) \mathbf{C} V(\mathbf{K\_s}) \rangle} \right]^{0.3},\tag{46}$$

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*(*Ks*) 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]:

$$\left< I(F) \right> = r[1 - M\_{K\_\ast}(K\_{\mathbb{C}}, 0)] + \frac{\psi(\theta\_\ast - \theta\_i) + F}{F} M\_{K\_\ast}(K\_{\mathbb{C}}, 1) \,, \tag{47}$$

*Water* **2018**, *10*, 1873

$$
\begin{split}
\langle t \rangle &= \frac{F}{r} \{ 1 - M\_{\mathcal{K}\_{\mathbf{s}}}[\mathcal{K}\_{\mathbf{c}}, 0] \} + \left[ F + \psi(\theta\_{\mathbf{s}} - \theta\_{\mathbf{i}}) \ln \left( \frac{\psi(\theta\_{\mathbf{s}} - \theta\_{\mathbf{i}})}{\overline{\psi(\theta\_{\mathbf{s}} - \theta\_{\mathbf{i}}) + F}} \right) \right] \\ & \quad \left\{ M\_{\mathcal{K}\_{\mathbf{s}}}[\mathcal{K}\_{\mathbf{c}\prime} - 1] \right\} + \psi(\theta\_{\mathbf{s}} - \theta\_{\mathbf{i}}) \sum\_{i=1}^{\infty} \frac{M\_{\mathcal{K}\_{\mathbf{s}}}[\mathcal{K}\_{\mathbf{c}}, i]}{(i+1)r^{i+1}}.
\end{split} \tag{48}
$$

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 *Ks* is fairly simple and requires limited computational effort.

The model for coupled heterogeneity of *r* and *Ks* (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 *Ks* can be neglected only when <*r*> <*K*> or for storm durations much greater than *tpa*; (3) the effects on <*I*> produced by significant values of *CV*(*Ks*) and *CV*(*r*) are similar; (4) run-on plays a significant role for moderate storms and high values of *CV*(*Ks*) and *CV*(*r*); and (5) the model can be simplified using Equations (47) and (48) for *CV*(*r*) substantially less than *CV*(*Ks*) and steady rainfalls.

#### **6. 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*(*Ks*), <*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 ye<sup>t</sup> 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.

**Author Contributions:** Investigation, Writing—original draft preparation and Writing—review and editing, R.M., C.C., C.S., A.F., J.D., R.S.G.

**Funding:** This research was partly financed by the Italian Ministry of Education, University and Research (PRIN 2015).

**Conflicts of Interest:** The authors declare no conflict of interest.
