Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments

: The sustainable exploitation of groundwater resources is a multifaceted and complex problem, which is controlled, among many other factors and processes, by water flow in porous soils and sediments. Modeling water flow in unsaturated, non-deformable porous media is commonly based on a partial differential equation, which translates the mass conservation principle into mathematical terms. Such an equation assumes that the variation of the volumetric water content ( θ ) in the medium is balanced by the net flux of water flow, i.e., the divergence of specific discharge, if source/sink terms are negligible. Specific discharge is in turn related to the matric potential ( h ), through the non-linear Darcy–Buckingham law. The resulting equation can be rewritten in different ways, in order to express it as a partial differential equation where a single physical quantity is considered to be a dependent variable. Namely, the most common instances are the Fokker–Planck Equation (for θ ), and the Richards Equation (for h ). The other two forms can be given for generalized matric flux potential ( Φ ) and for hydraulic conductivity ( K ). The latter two cases are shown to limit the non-linearity to multiplicative terms for an exponential K -to- h relationship. Different types of boundary conditions are examined for the four different formalisms. Moreover, remarks given on the physico-mathematical properties of the relationships between K , h , and θ could be useful for further theoretical and practical studies.


Introduction
The most recent report on groundwater from the UNESCO World Water Assessment Programme [1] recognizes the huge uncertainty about the estimates of the volume of fresh water on the Earth, but nevertheless considers that the vast majority (possibly up to 99%) of freshwater is stored in aquifers.Furthermore, the same source of information reports that groundwater is estimated to provide almost one-half of the volume of water withdrawn for domestic use by the global population, and about one-fourth of water withdrawn for irrigation, serving more than one-third of the world's irrigated land.
These numbers, although affected by large uncertainty, are nevertheless sufficient to demonstrate the relevance of groundwater to life and the development of the planet Earth and its population.
These remarks prompted the scientific community to consider the protection and the proper management of groundwater resources, also within the realm of the Sustainable Development Goals (SDGs) identified by the UN General Assembly, with the adoption of Resolution A/RES/70/1 on "Transforming Our World: The 2030 Agenda for Sustainable Development" in 2015.One of the SDGs (6-Clean water and sanitation) is explicitly devoted to the availability and sustainable management of freshwater for all the world's population, but the other SDGs are strictly related to the availability of good quality water.
For a general review of the discussion of the sustainability of groundwater resources, the reader is referred to a recent paper by Bierkens and Wada [2] and to the references therein.In that paper, "physically sustainable groundwater use" is defined as the "prolonged (multi-annual) withdrawal of groundwater from an aquifer in quantities not exceeding average annual replenishment, resulting in dynamically stable water tables or hydraulic heads".However, the authors rightly criticize the simplistic statement that the sustainable average annual withdrawal is equal to the annual aquifer recharge.This approach was already mentioned as the "safe yield myth" by Bredehoeft [3].Instead, the interaction between surface-and ground-water is very important in many cases and a good review on this topic is given by Paniconi and Putti [4].Contaminant transport in soil and groundwater should also be considered, as it can limit the availability of good quality water, especially for drinkable purposes.In regions characterized by a well-developed and widespread irrigation network, watering methods and practices interfere in a rather complex way with groundwater, as shown, e.g., by Vassena et al. [5].Finally, water, obviously including groundwater, impacts several ecosystem services [6].
The above remarks show that sustainable management of groundwater resources is very challenging, as it is a multifaceted issue.However, a crucial point remains the behavior of water in soils and in shallow alluvial sediments, which can be considered to be partially saturated porous media.Water flow in variably saturated porous soils and alluvial sediments affects crop production, recharge of phreatic aquifers, transport of fertilizer and contaminant transport, and many other processes.Therefore, its study is of utmost importance in soil physics, hydrology, water resources management, agronomy, and so on.
Within this framework, this paper provides a theoretical view of groundwater flow in partially saturated porous media.In particular, the basic equation, founded on the mass conservation principle and the Darcy-Buckingham law, yields a non-linear partial differential equation, which involves the volumetric water content (θ) and the matric potential (h).Such an equation can be expressed in several different forms, according to the physical quantity considered to be the dependent variable, with respect to which the equation is solved.The fundamental goal of this work is to review four of the different forms of the flow equations: the Fokker-Planck equation, which is written in terms of θ; the Richards equation, which is written in terms of h; the equation for the generalized matric flux potential (Φ); the equation for hydraulic conductivity (K).
Instead, this paper is devoted to giving some physical or mathematical insights into the different formulations of the flow equation for unsaturated porous media, with specific attention, for each formulation, to the physical assumptions at its basis and to the form of the boundary conditions, including the conditions at the sharp interface between two media characterized by different physical properties.
Moreover, physical and mathematical conditions, which should be satisfied by the phenomenological relations between the physical variables used to determine the state and the flow of water in porous media, are discussed.In the author's opinion, these conditions should be carefully considered by researchers working on different methods (e.g., inverse modeling, capillary tube, and fractal-based models, etc.) to predict retention curves, conductivity curves, and other semi-empirical relationships.In fact, they are fundamental for a proper, physicallybased simulation of water flow in variably saturated porous media.
The most innovative content of this work is twofold and it is mainly related to Φand Kbased equations.First, Φ is shown to be a specific case within a family of variables, in particular, it permits slightly simplifying the non-linearity of the flow equations.Second, the K-based equation is proposed for a general case, whereas its introduction in the literature is usually limited to the application of an exponential dependence of conductivity and matric potential on water content.

Basic Equations and Notation
The fundamental, basic equation, that describes moisture content distribution and flow in soils and unsaturated sediments, is derived from the continuity equation for water phase [16], which, in the absence of source or sink terms, reads where: ) is the water flow through a unit surface of a porous medium.
The second ingredient for modeling groundwater flow in partially saturated porous media is the Darcy-Buckingham law, which states the proportionality between water flux and hydraulic head gradient and reads where: ) is the hydraulic head or soil water potential per unit weight (z axis is taken positive upwards), ) is the matric potential, often referred to as suction; recall that for saturated porous media, it is substituted by the pressure head (p − p a ) • (ρ w g) −1 , where p is the water pressure, p a is the atmospheric pressure, ρ w is the water density and g is the gravity acceleration.
From Equations (1) and (2), one obtains the following fundamental equation: It must be stressed that Equation (3) is based on the following assumptions: 1. isothermal conditions; 2.
water and, even more acceptable, solid grains are considered to be incompressible, i.e., their density is assumed to be constant in space and time; 3.
the soil or sediment is assumed to be non-deformable; 4.
the influence of the air phase on water flow is negligible; 5.
hysteresis effects, which would give a non-unique dependence between h and θ, are disregarded; such an approximation is acceptable either if the medium does not show any difference in the h-to-θ relationship during the wetting and draining phases or if only one of the two phases (wetting or draining) is considered; 6.
the medium is supposed to be homogeneous, i.e., the functions relating K and h to θ are the same at every point of the domain; 7.
the medium is assumed to be porous.
With respect to the last point, notice that (1) holds also for a fractured medium, which possibly behaves as anisotropic, so that (2) is often used, but K is a second-order symmetric tensor.In the case of a fractured medium, the physical meaning of the quantities appearing in (1)-(3), should be modified and adapted to this kind of medium.
With respect to the point 3 above, recent advances in modeling flow in deformable media, by taking into account the presence of macropores also, are given by Carrillo 128 and Bourg [17], who apply a volume-averaging approach to upscale physics from the microscopic pore-and-grain scale to the Representative Elementary Volume (REV) scale [18].Moreover, the study of preferential flow, due to funneling, fingering, and macropores, has been recently reviewed by Nimmo [19].
The main difficulties in solving (3) are due to its non-linearity.
Formulations of (3) found in the scientific literature differ mainly for the dependent variable, which the equation is solved for.In the following subsections, four instances are considered: 1.
the Fokker-Planck equation, obtained when (3) is expressed as an equation for θ; 2.
the Richards equation, when (3) is transformed in an equation for h; 3.
Φ-based equation, when the matric flux potential, Φ, is introduced and used as a dependent variable; 4.
K-based equation, when the dependent variable is the K function.
Before discussing these equations in detail, some notations and definitions have to be introduced.Let K = K(θ) and h = ĥ(θ) indicate the functional dependence of K and h on θ, through the functions K and ĥ, respectively.Finally, let the hydraulic diffusivity, D ([length 2 • time −1 ]), be defined as

The Fokker-Planck Equation
If H = ĥ(θ) + z is substituted in (3), then This is a well-known equation of mathematical physics, often referred to as the Fokker-Planck Equation [20,21].An interesting and fundamental example of a solution of this equation was found by [22].
It is apparent that non-linearity involves the right-hand side of (5), where spatial derivatives are found.

If
d ĥ the functional relationship h = ĥ(θ) can be inverted to find θ = θ(h) such that ĥ( θ(h)) = h and θ( ĥ(θ)) = θ.K can also be expressed as a function of h in the following way: ), be defined as Then, (3) can be written in the following form: The first presentation of this equation goes back to Richards' work in 1931 [23].Its advantages against the Fokker-Planck equation are that h is a continuous function across boundaries between different soils or sediments, whereas θ is discontinuous at the interface separating two media characterized by different physical properties (see Section 3.3.4).In addition, this equation directly controls the variable h and, as a consequence, H, whose gradient is the "driving force" for water to flow in porous media.Anyway, whereas ( 5) is non-linear because of the non-linearity of physical parameters appearing under the spatial derivatives, in (8) the coefficient multiplying the temporal derivative of the solution depends on the solution itself, thus increasing the non-linearity of the problem.

Generalized Matric Flux Potential Equation
Several analytical solutions can be obtained for the steady-state case when where K 0 and α are constants, by introducing the matric flux potential [24][25][26], Φ ([length 2 • time −1 ]), defined by: h= ĥ(θ) This approach has been called "quasilinear approximation" by Pullan [27], who proposed an interesting discussion also about the K-to-h relationship.This quantity has been used also to determine the onset of limiting hydraulic conditions for root water uptake [28].
Here, it is shown that the definition of Φ given by ( 10) is the most natural choice among a family of new variables, which satisfies a condition for simplifying (3).
In particular, let a new variable Θ be searched for, such that θ = ω(Θ) with ω an invertible and regular, but otherwise arbitrary, function.
From theorems about differentiation of composite functions, one obtains: and by the definition of hydraulic diffusivity, one has and, as a consequence, The following property can now be stated h= ĥ(θ) and finally one obtains It is trivial to recognize that Θ coincides with the transformation originally proposed by Kirchhoff i in 1894 [29] and reduces to Φ, as defined by (10) Equation ( 11) can now be simplified and rewritten as follows: If (9) which, under steady-state conditions, reduces to the following simple linear equation: A long list of papers presenting analytical solutions for this linear steady-state equation, with different geometries of the domain and boundary conditions, could be cited.Among the others, it is worth recalling some seminal works [24,30,31].Studies about numerical solutions go back at least to Davidson [32] and time-dependent solutions were obtained by Braester and Batu [25,33].
Two final remarks.From the mathematical point of view, among all regular invertible transformations, ω(Θ), the family of functions implicitly defined by ( 12) collects all the transformations that cancel the coefficient of |gradΘ| 2 appearing in (11).
On the other hand, water flow in the porous medium can be expressed through a simple equation, in which Θ appears explicitly: where k is the upward-directed vertical unit vector.Once again, if (9) holds, then If ( 14) is applied and h 0 is chosen to be such that K(h 0 ) 0, then q = −gradΘ − αΘk.

K-Based Equation
Since hysteresis effects are disregarded and experimental data show that K is strictly increasing with θ and h, both ĥ and K functions can be inverted, so that one can write θ = θ(h) and h = h(K), with ĥ( θ(h)) = h and K( h(K)) = K.
By expressing (3) in terms of the K parameter, one obtains If the following notation is introduced then In analogy with Section 2.4, ( 22) can be simplified when ( 9) is valid.In fact, from (21) it follows A(K) = (αK) −1 and Equation ( 22) reads Non-linear terms now appear only in front of the time derivative of the solution, a remark analogous to that related to (16); this property significantly simplifies the problem, above all under steady-state conditions.In fact, the steady-state equation is easily seen to be equivalent to (16), thanks to (19), valid when the exponential K-to-h relation holds.
Notice that a similar approach was applied by Srivastava and Yeh [34], who introduced an even more simple version of (23) for an exponential dependence of θ on h, and then followed also by other researchers for theoretical studies or to validate numerical models (see, e.g., [35][36][37], among the others).

Discussion about Initial and Boundary Conditions and about Phenomenological
Relationships among θ, h and K Equations ( 5), ( 8), ( 15) and ( 22) must be supported by phenomenological laws relating K, h, and θ and must be accompanied by proper initial and boundary conditions.Therefore, this section is devoted to a discussion of the physical and mathematical properties behind some of the most common phenomenological laws and behind the different types of boundary conditions.

Properties of the Relationship of K and h to θ
The physical, chemical, and biological processes which occur in unsaturated soils and sediments are quite complex and they obviously have non-trivial impacts on the relationship among h, θ, and K. Indeed, several simplifications have already been, implicitly or explicitly, included in the derivation of Equations ( 5), ( 8), ( 15) and (22).Experimental data about the relationships among h, θ, and K have been reviewed in several textbooks (e.g., [16,38,39]) and in a large set of papers.Here it is worth mentioning the data set compiled very recently by Hohenbrink et al. [40], the UNSODA database [41], and the papers by Peters [42], who tested empirical models, which distinguish between capillary, adsorptive and film components, with data from ten different soils, and by Luo et al. [43], who showed results for five representative soils, out of 27 different soils based on a compilation of data from eight papers.
The first part of this sub-section summarizes some of the results presented in papers, which have not found wide diffusion in the international scientific literature [44][45][46].In order to do this, it is assumed that the functions ĥ(θ), K(θ) are sufficiently regular, so that they can be differentiated as many times as required.
First of all, θ varies in the interval between θ irr , the irreducible water content, and θ sat , the water content at saturation.Both quantities can be defined for a REV of the porous medium: θ irr is the ratio between the volume of water that adheres on the surface of soil grains inside the REV, so that it can be extracted only if the energy given to the medium is much higher than that involved in normal, purely mechanical processes, and the volume of the REV itself; θ sat is the ratio between the volume of water at saturation in the REV and the volume of the REV itself.Notice that if total porosity is denoted by , then 0 ≤ θ irr ≤ θ ≤ θ sat ≤ .
The matric potential is a negative quantity, so that ĥ where h e is called the air entry value and is a threshold, such that air can enter the porous medium only when suction is lower than h e .The function ĥ(θ) can assume any value between 0 and h e for θ = θ sat .When θ is small, i.e., if θ → θ irr , it is necessary to provide a lot of energy to the porous medium in order to dry it.In fact, the "free" water inside the pores can flow out from a sample of the porous medium by gravity only, but heating with an oven is necessary to extract the water involved in capillary processes.In other words, the energy required to extract small quantities of water from an almost dry porous medium is much greater than that involved in the mechanical processes occurring in most of the field conditions.Therefore, this could be expressed in mathematical form as lim θ→θ irr ĥ(θ) = −∞. ( This means that the function ĥ(θ) shows a vertical asymptote for θ irr and, as a consequence, lim Both experimental evidence and studies on capillary bundle models show that ĥ is a monotonously non-decreasing function (see, e.g., [40]), so that d ĥ and lim i.e., the function ĥ(θ) has a vertical tangent while approaching θ sat , a value of θ for which the function ĥ is not uniquely defined, as mentioned in the remark following (24) above.Equations ( 26) and ( 28) imply that d 2 ĥ dθ 2 (θ i.p. ) = 0, for some θ i.p. ∈ (θ irr , θ sat ), i.e., ĥ(θ) has an inflection point at θ = θ i.p. .The first trivial condition for K is that it should be positive for the whole range of θ values, with the exception of θ = θ irr .Moreover, it should tend to K sat , the value of hydraulic conductivity for a saturated medium, when θ θ sat : It is experimentally verified that K increases with soil water content [42,43], so that Since hydraulic diffusivity is always positive, Equations ( 30) and ( 4) imply d ĥ Equation ( 33) implies condition (6), so that ĥ(θ) is invertible in the interval (θ irr , θ sat ), and therefore a single-valued function θ(h) such that θ( ĥ(θ)) = θ and ĥ( θ(h)) = h can be found, as it is required to introduce the Richards and the K-based equations (see Sections 2.3 and 2.5).
Moreover, it is worth recalling that (32) implies that also ( 13) is invertible and it makes sense to consider θ = ω(Θ).
Experiments show that D is a monotonically non-decreasing function, so that Possible discrepancies from this behavior are disregarded here, but could occur in small intervals of θ values, where vapor movement could make processes more complex [39].
Equation (28) implies lim The formula In particular, lim The second part of this subsection is devoted to some properties, which are based on the results of studies [47,48] which show that water infiltration in initially dry soil creates a wetting front, separating the wet region from the dry region.These properties were used by other authors to study the propagation of wetting fronts (e.g., [49]).For linear diffusion equations, it is well known that an infinite velocity of propagation is obtained, which is not physically acceptable.A modification of the diffusion equation, which accounts for the possible presence of a front, has been proposed for linear systems [50], but has not yet been given great attention in the scientific literature.
Table 1 provides a complete list of the properties of different functions.These properties can be used to derive basic properties of the other functions that have been defined in Section 2: θ(h), K(h), C(h), A(K), and B(K).
The conditions discussed in this subsection should be verified by any model used to fit field data, because they are necessary to guarantee that a given model for the relationships among soil hydraulic properties correctly reproduces the physical behavior of natural systems.Recall that the functional relationships among soil hydraulic properties proposed in the scientific literature and/or the phenomenological parameters appearing in the corresponding formulas are usually determined with experimental procedures (field tests; laboratory experiments on samples) or with inverse modeling [51][52][53].With specific reference to field data inversion, general reviews of the discrete inverse problem (see, e.g., [54], and references therein) or specific applications to flow and transport in unsaturated media (e.g., [55,56]) can be found in the scientific literature.
Table 1.Properties of different functions.Notice that in this table f and f denote, respectively, the first-and second-order derivatives of a function f , with respect to its argument.

Critical Analysis of the Most Common Phenomenological Laws
Starting from the relationship between suction and volumetric water content, namely, the h-θ relationship, one of the oldest fitting expressions is that proposed by Brooks and Corey in 1964 [57]: ĥ where is the saturation function and λ is a parameter often called the pore-size distribution index.The limitations of ( 40) are due to the fact that it does not obey (28), so it can give problems close to saturation, where S = 1.Moreover, it is basically a power function, which does not show any inflection point, contrary to what is expected from (29).
The most common fitting relationship is that proposed by van Genuchten in 1980 [58]: where h 0 < 0 and m and n are two positive quantities.
Moving to the K function, i.e., to the K-to-θ relationship, the classical approach by Kozeny in 1927 and revised by Carman in 1937 [59,60] yields where α is a fitting parameter.If α > 1, the conditions introduced in Section 3.1 are satisfied.This is in agreement with the typical values obtained from fitting experimental data.Burdine's and Mualem's approaches [61,62] are often coupled with the van Genuchten h-θ relationship (42).Leij et al. [63] provided a long list of closed-form expressions that could be considered and analyzed under the framework described in this work.

Initial and Boundary Conditions
Initial conditions necessary to solve evolution equations should describe the state of the system at the starting time, conventionally chosen as t = 0, of the modeled time interval.
In the case of partially saturated porous media, under the conditions mentioned in Section 3.1, initial conditions can be expressed by fixing the values of θ, h, Φ, or K at time t = 0, at any point of the domain under study.
The discussion is more complex for boundary conditions and requires further remarks.From the physical point of view, three cases can be considered.

Prescribed Value of Water Content at Saturation
In groundwater hydrology and soil physics, one can prescribe the water content at the border of the domain if the evolution of the physical system does not affect the state of the external water bodies, with which the system under study is in contact.The assignment of this boundary condition basically requires perfect hydraulic contact between the porous medium and the external water bodies.In principle, this condition can be prescribed at the water table, where θ = θ sat or, equivalently, S = 1.Notice, however, that two requisites should be met for rigorous use of this kind of boundary condition: 1.
the position of the water table should be independent of the water flow in the vadose zone, i.e., the unsaturated portion of the subsurface; 2.
the capillary fringe should have a negligible thickness, otherwise, more complex approaches are necessary for a rigorous treatment [16].
For this case, the boundary condition can be written as θ = θ sat for the Fokker-Planck equation, h = h e for the Richards equation, Φ = h e h 0 K(h ) dh for the Φ-based equation, and K = K sat for the K-based equation.From the mathematical point of view, these are Dirichlet boundary conditions for all the four types of equations considered in this work.

Prescribed Value of q • n
Here n denotes the inward-pointing unit vector perpendicular to the border of the domain.In this case, one assumes that the value of the incoming flux per unit surface can be prescribed.The most common situations, where this condition can be applied, are the contact with impermeable materials and at the ground surface.For the latter instance, the downward water flux (or upward moisture flux in case of prevailing evapotranspiration) can be evaluated by assessing the energy and mass balance at the ground, possibly with Soil-Vegetation-Atmosphere-Transfer (SVAT) models (see [64][65][66][67], for examples showing the use of different methods, some of which are physically-based, and of meteorological data).Thus, it is possible to estimate q • n at the ground surface.This condition, when coupled with the Darcy-Buckingham law, means that one prescribes the value of for the Fokker-Planck Equation ( 5), for the Richards Equation ( 8), for the Θ-based Equation ( 15), as obtained from ( 18) and ( 14), for the K-based Equation (22).These equations cannot be assimilated in a straightforward way to Neumann or Robin boundary conditions.They involve both the normal component of the derivative of the solution, but with non-linear coefficients for ( 44), (45), and (47), and the solution itself, through non-linear dependencies for ( 44)-( 46).

Interaction between Groundwater and Surface Water
The third kind of boundary condition, that is considered in this work, corresponds to the contact of soils or sediments with water bodies, but with the fundamental difference, with respect to the case of Section 3.3.1,that the hydraulic contact is not perfect.Therefore, the water saturation, and the other quantities for which the equations are solved for, do not perfectly and instantaneously follow the saturation of the water body with which they are in contact.The most common example, for which this boundary condition can be used, is when soil or sediments are in contact with a river, but the water exchange is controlled by low-permeability sediments on the stream-bed.If H wb is the hydraulic level of the external water body, which is assumed to be independent of the level inside the soils and sediments, and if K wb ([time −1 ]) is the conductance of the poorly permeable sediments, which connects the water body with the vadose zone, then q • n = K wb (H wb − H).As a consequence, these boundary conditions transform into for the Fokker Planck Equation ( 5), for the Richards Equation ( 8), for the Θ-based Equation ( 15), as obtained from ( 18) and ( 14), for the K-based Equation (22).From the mathematical point of view, these equations appear in the form of non-linear Robin boundary conditions.

Conditions at the Interface between Media with Different Physical Properties
In the presence of a sharp interface between two media with different physical properties, i.e., for which the relationship between soil hydraulic quantities is different, two basic physical conditions can be introduced: continuity of h and continuity of the normal component of q across the interface.
These conditions are quite complex, once again, due to non-linearity, with the exception of the continuity of h for the Richards equation, which is directly related to the dependent variable.
Notwithstanding the fact that the Fokker-Planck and Richards equations are perhaps the most natural formulations of (3), they are affected by strong non-linearity.The Richards equation is possibly the best choice for problems involving heterogeneous porous media, e.g., for layered media (see Section 3.3.4),or simultaneous saturatedunsaturated flow.
The generalized matric flux potential equation is a simplified version of (3) and generalizes the matric flux equation introduced by other Authors.In particular, matrix flux potential is the most natural choice among the transformations which allow canceling one of the terms of the spatial derivatives out of the fundamental equation.In particular, when (13) holds, then one obtains the simple form of the Θ-based equation, given by (15), where non-linear terms appear only as multiplicative factors of the laplacian and of the vertical derivative of Θ.
When an exponential K-to-h relation is valid, both Θ-based and K-based equations present non-linear coefficients only as multiplicative terms in front of time or spatial derivatives, thus giving linear steady-state equations for which several analytical solutions are available.Unfortunately, this dependence usually does not fit the trend of K over the whole range of h values, but it can be applied only for problems that consider a small range of h and θ values.
The non-linearity of the equations makes it quite complex the definition of boundary conditions.Dirichlet boundary conditions can be prescribed if the water table is independent of the flow in the vadose zone and if the capillary fringe can be neglected; if this is the situation, a saturation condition can be imposed at the water table.Non-linear Robin conditions correspond to the cases when water flux through the border of the study domain can be prescribed or when the interaction between the porous medium and superficial water bodies is examined.
Finally, it is important to recall that some remarks introduced in the previous sections could be relevant both in theoretical and practical studies.For instance, they can drive the research on the mathematical properties of the inverse problem for unsaturated flow in porous media, or on innovative methods for the determination of retention curves (see, e.g., [68][69][70][71][72], for fractal-based models) but they can also have an impact on the dynamics of porous media [73] and on the modeling of solute transport (see, e.g., [74,75]).