Next Article in Journal
Expert-Based Assessment of the Potential of Agroforestry Systems in Plain Regions across Bihor County, Western Romania
Next Article in Special Issue
Spatial and Temporal Patterns of Groundwater Levels: A Case Study of Alluvial Aquifers in the Murray–Darling Basin, Australia
Previous Article in Journal
Regenerative Agriculture as a Sustainable System of Food Production: Concepts, Conditions, Perceptions and Initial Implementations in Poland, Czechia and Slovakia
Previous Article in Special Issue
Assessing Impacts of Mining-Induced Land Use Changes on Groundwater and Surface Water Quality Using Isotopic and Hydrogeochemical Signatures
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments

Dipartimento di Scienze Della Terra “Ardito Desio”, Università Degli Studi di Milano, Via Botticelli 23, 20133 Milano, Italy
Sustainability 2023, 15(22), 15723; https://doi.org/10.3390/su152215723
Submission received: 8 September 2023 / Revised: 4 November 2023 / Accepted: 6 November 2023 / Published: 8 November 2023
(This article belongs to the Special Issue Groundwater, Soil and Sustainability)

Abstract

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

1. 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).
This paper is not a comprehensive review of modeling unsaturated flow in soils and sediments, which can be found elsewhere (e.g., [7,8,9,10,11,12,13]). Moreover, this paper does not consider numerical models (for this aspect, see, e.g., [14,15]).
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, physically-based 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 K- based 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.

2. Equations of Water Flow in Unsaturated Soils and Sediments

2.1. 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
θ t = div q ,
where:
  • θ (dimensionless) is the volumetric water content,
  • q ( [ length · time 1 ] ) 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
q = K grad H ,
where:
  • K ( [ length · time 1 ] ) is the ( θ dependent) hydraulic conductivity,
  • H = h + z ( [ length ] ) is the hydraulic head or soil water potential per unit weight (z axis is taken positive upwards),
  • h < 0 ( [ length ] ) 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:
θ t = div K grad H .
It must be stressed that Equation (3) is based on the following assumptions:
  • isothermal conditions;
  • 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;
  • the soil or sediment is assumed to be non-deformable;
  • the influence of the air phase on water flow is negligible;
  • 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;
  • 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;
  • 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:
  • the Fokker–Planck equation, obtained when (3) is expressed as an equation for θ ;
  • the Richards equation, when (3) is transformed in an equation for h;
  • Φ -based equation, when the matric flux potential, Φ , is introduced and used as a dependent variable;
  • 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 = h ^ ( θ ) indicate the functional dependence of K and h on θ , through the functions K ^ and h ^ , respectively. Finally, let the hydraulic diffusivity, D ( [ length 2 · time 1 ] ), be defined as
D = D ^ ( θ ) K ^ ( θ ) d h ^ d θ ( θ ) ,

2.2. The Fokker–Planck Equation

If H = h ^ ( θ ) + z is substituted in (3), then
θ t = div K ^ ( θ ) d h ^ d θ ( θ ) grad θ + K ^ ( θ ) grad z ,
i.e.,
θ t = div D ^ ( θ ) grad θ + d K ^ d θ ( θ ) θ z .
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.

2.3. The Richards Equation

If
d h ^ d θ ( θ ) 0 ,
the functional relationship h = h ^ ( θ ) can be inverted to find θ = θ ^ ( h ) such that h ^ ( θ ^ ( h ) ) = h and θ ^ ( h ^ ( θ ) ) = θ . K can also be expressed as a function of h in the following way: K = K ˜ ( h ) K ^ ( θ ^ ( h ) ) . Let the moisture capacity function, C ( [ length 1 ] ), be defined as
C ( h ) d θ ^ d h ( h ) .
Then, (3) can be written in the following form:
C ( h ) h t = div K ˜ ( h ) grad h + d K ˜ d h ( h ) h z .
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.

2.4. Generalized Matric Flux Potential Equation

Several analytical solutions can be obtained for the steady-state case when
K ˜ ( h ) = K 0 exp ( α h ) ,
where K 0 and α are constants, by introducing the matric flux potential [24,25,26], Φ ( [ length 2 · time 1 ] ), defined by:
Φ θ 0 θ D ^ ( θ ) d θ = h 0 = h ^ ( θ 0 ) h = h ^ ( θ ) K ˜ ( h ) d 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:
θ t = d ω d Θ Θ t , grad h = d h ^ d θ d ω d Θ grad Θ , grad K ˜ = d K ˜ d h d h ^ d θ d ω d Θ grad Θ
and, if Δ denotes the Laplacian operator,
Δ h = div grad h = grad d h ^ d θ d ω d Θ · grad Θ + d h ^ d θ d ω d Θ Δ Θ = d 2 h ^ d θ 2 d ω d Θ 2 grad Θ 2 + d h ^ d θ d 2 ω d Θ 2 grad Θ 2 + d h ^ d θ d ω d Θ Δ Θ .
Then
div K ˜ grad h + K ˜ z = grad K ˜ · grad h + K ˜ Δ h + K ˜ z = d K ˜ d h d h ^ d θ d ω d Θ 2 grad Θ 2 + K ˜ d 2 h ^ d θ 2 d ω d Θ 2 + d h ^ d θ d 2 ω d Θ 2 grad Θ 2 + K ˜ d h ^ d θ d ω d Θ Δ Θ + d K ˜ d h d h ^ d θ d ω d Θ Θ z
and finally
d ω d Θ Θ t = d d Θ K ˜ d h ^ d θ d ω d Θ grad Θ 2 + K ˜ d h ^ d θ d ω d Θ Δ Θ + d K ˜ d h d h ^ d θ d ω d Θ Θ z .
Equation (11) can be simplified, if ω is chosen such that the coefficient of grad Θ 2 disappears. This means that ω satisfies
d d Θ K ˜ d h ^ d θ d ω d Θ = 0 ,
i.e.,
K ˜ ( h ˜ ( ω ( Θ ) ) ) d h ^ d θ ( ω ( Θ ) ) d ω d Θ ( Θ ) = c ,
with c = constant . From K ˜ ( h ˜ ( ω ( Θ ) ) ) = K ^ ( ω ( Θ ) ) and by the definition of hydraulic diffusivity, one has
D ^ ( ω ( Θ ) ) d ω d Θ ( Θ ) = c ,
and, as a consequence,
θ 0 = ω ( Θ 0 ) θ = ω ( Θ ) D ^ ( θ ) d θ = Θ 0 Θ D ^ ( ω ( Θ ) ) d ω d Θ ( Θ ) d Θ = c Θ Θ 0 .
The following property can now be stated
h 0 = h ^ ( θ 0 ) h = h ^ ( θ ) K ˜ ( h ) d h = θ 0 θ K ˜ ( h ^ ( θ ) ) d h ^ d θ ( θ ) d θ = θ 0 θ D ^ ( θ ) d θ = c Θ Θ 0 ,
and finally one obtains
Θ = Θ 0 + 1 c θ 0 θ D ^ ( θ ) d θ = Θ 0 + 1 c h 0 h K ˜ ( h ) d h .
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), if
c = 1 , Θ 0 = 0 .
Equation (11) can now be simplified and rewritten as follows:
Θ t = D ^ ( ω ( Θ ) ) Δ Θ + d K ˜ d h ( h ^ ( ω ( Θ ) ) ) d h ^ d θ d ω d Θ ( ω ( Θ ) ) Θ z .
If (9) holds, then
d K ˜ d h = α K ˜
and (15) reads
Θ t = K ˜ d h ^ d θ Δ Θ + α Θ z ,
which, under steady-state conditions, reduces to the following simple linear equation:
Δ Θ + α Θ z = 0 .
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:
q = c grad Θ K k ,
where k is the upward-directed vertical unit vector. Once again, if (9) holds, then
K = α c Θ Θ 0 + K ˜ ( h 0 ) .
If (14) is applied and h 0 is chosen to be such that K ˜ ( h 0 ) 0 , then q = grad Θ α Θ k .

2.5. K-Based Equation

Since hysteresis effects are disregarded and experimental data show that K is strictly increasing with θ and h, both h ^ and K ˜ functions can be inverted, so that one can write θ = θ ^ ( h ) and h = h ˜ ( K ) , with h ^ ( θ ^ ( h ) ) = h and K ˜ ( h ˜ ( K ) ) = K .
By expressing (3) in terms of the K parameter, one obtains
d θ ^ d h d h ˜ d K K t = div K d h ˜ d K grad K + K z .
If the following notation is introduced
A ( K ) d h ˜ d K ( K ) = d K ˜ d h ( h ˜ ( K ) ) 1 and B ( K ) d θ ^ d h ( h ˜ ( K ) ) d h ˜ d K ( K ) ,
then
B ( K ) K t = div K A ( K ) grad K + K z .
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
B ( K ) K t = 1 α Δ K + K z .
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).

3. 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.

3.1. 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 h ^ ( θ ) , 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
h ^ ( θ ) h e < 0 , θ [ θ irr , θ sat ] ,
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 h ^ ( θ ) 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 h ^ ( θ ) = .
This means that the function h ^ ( θ ) shows a vertical asymptote for θ irr and, as a consequence,
lim θ θ irr d h ^ d θ ( θ ) = + , lim θ θ irr d 2 h ^ d θ 2 ( θ ) = .
Both experimental evidence and studies on capillary bundle models show that h ^ is a monotonously non-decreasing function (see, e.g., [40]), so that
d h ^ d θ ( θ ) 0 , θ [ θ irr , θ sat ) ,
and
lim θ θ sat d h ^ d θ ( θ ) = + , lim θ θ sat d 2 h ^ d θ 2 ( θ ) = +
i.e., the function h ^ ( θ ) has a vertical tangent while approaching θ sat , a value of θ for which the function h ^ is not uniquely defined, as mentioned in the remark following (24) above.
Equations (26) and (28) imply that
d 2 h ^ d θ 2 ( θ i . p . ) = 0 , for some θ i . p . ( θ irr , θ sat ) ,
i.e., h ^ ( θ ) 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 :
K ^ ( θ ) > 0 , θ ( θ irr , θ sat ] ; K ^ ( θ irr ) = 0 ; K ^ ( θ sat ) = K sat .
It is experimentally verified that K increases with soil water content [42,43], so that
d K ^ d θ ( θ ) > 0 , θ [ θ irr , θ sat ] .
Since hydraulic diffusivity is always positive,
D ^ ( θ ) > 0 , θ ( θ irr , θ sat ] ,
Equations (30) and (4) imply
d h ^ d θ ( θ ) > 0 , θ ( θ irr , θ sat ) .
Equation (33) implies condition (6), so that h ^ ( θ ) is invertible in the interval ( θ irr , θ sat ) , and therefore a single-valued function θ ^ ( h ) such that θ ^ ( h ^ ( θ ) ) = θ and h ^ ( θ ^ ( h ) ) = h can be found, as it is required to introduce the Richards and the K-based equations (see Section 2.3 and Section 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
d D ^ d θ ( θ ) 0 , θ [ θ irr , θ sat ] .
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 θ θ sat D ^ ( θ ) = + .
The formula
d D ^ d θ = d K ^ d θ d h ^ d θ + K ^ ( θ ) d 2 h ^ d θ 2 0 , θ [ θ irr , θ sat ]
implies
d 2 h ^ d θ 2 K ^ 1 d K ^ d θ d h ^ d θ , θ [ θ irr , θ sat ] .
In particular,
lim θ θ irr K ^ 1 d K ^ d θ d h ^ d θ = + .
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.
Gilding [47,48] showed that two conditions are sufficient to create a wetting front in non-linear systems:
D ^ ( θ ) θ θ irr be integrable for θ [ θ irr , θ sat ] ;
K ^ ( θ ) m ( θ θ irr ) , for some m > 0 and θ θ irr .
Condition (36) implies
D ^ ( θ irr ) = 0 ,
whereas (37) implies
K ^ ( θ irr ) = 0 ,
in perfect agreement with (38).
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.

3.2. 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]:
h ^ ( θ ) = h e S 1 / λ ,
where
S ( θ ) = θ θ irr θ sat θ irr
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]:
h ^ ( θ ) = h 0 S 1 / m 1 1 / n ,
where h 0 < 0 and m and n are two positive quantities.
If n < 1 , (42) violates (28).
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
K ^ ( θ ) = K sat S α ,
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.

3.3. 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.

3.3.1. 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:
  • 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;
  • 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 0 h e K ˜ ( h ) d h 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.

3.3.2. 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
q · n = K ^ ( θ ) grad ( h ^ ( θ ) + z ) · n = D ^ ( θ ) grad θ · n K ^ ( θ ) k · n
for the Fokker-Planck Equation (5),
q · n = K ˜ ( h ) grad ( h + z ) · n
for the Richards Equation (8),
q · n = grad Θ + K ^ ( ω ( Θ ) ) k · n
for the Θ -based Equation (15), as obtained from (18) and (14),
q · n = K A ( K ) grad K + k · n
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).

3.3.3. 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
K wb ( H wb h ^ ( θ ) z ) = D ^ ( θ ) grad θ · n K ^ ( θ ) k · n
for the Fokker Planck Equation (5),
K wb ( H wb h z ) = K ˜ ( h ) grad ( h + z ) · n
for the Richards Equation (8),
K wb ( H wb h ˜ ( ω ( Θ ) ) z ) = grad Θ + K ^ ( ω ( Θ ) ) k · n
for the Θ -based Equation (15), as obtained from (18) and (14),
K wb ( H wb h ˜ ( K ) z ) = K A ( K ) grad K + k · n
for the K-based Equation (22). From the mathematical point of view, these equations appear in the form of non-linear Robin boundary conditions.

3.3.4. 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.

4. Conclusions

The comparison among the four Equations (5), (8), (15) and (22) shows several interesting properties which are briefly summarized here.
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 saturated-unsaturated 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]).

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No data were created or analyzed in this study.

Acknowledgments

The guest editors of the Special Issue on “Groundwater, soil and sustainability” are warmly acknowledged for the invitation to submit this contribution. Three anonymous reviewers are sincerely acknowledged, because their comments and questions were very useful to improve the quality of the manuscript during the revision process. This research has been partially conducted in the framework of the project “Geosciences for society: resources and their evolution” supported by the Italian Ministry of University and Research (MUR) through the funds ‘Dipartimenti di Eccellenza 2023/2027’. Multidisciplinary Digital Publishing Institute is warmly acknowledged for having waived APC and the editorial office is kindly acknowledged for the high professional level with which this submission has been managed.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. UNESCO World Water Assessment Programme. The United Nations World Water Development Report 2022: Groundwater: Making the Invisible Visible. 2022. Available online: https://unesdoc.unesco.org/ark:/48223/pf0000380721.locale=en (accessed on 5 November 2023).
  2. Bierkens, M.F.P.; Wada, Y. Non-renewable groundwater use and groundwater depletion: A review. Environ. Res. Lett. 2019, 14, 063002. [Google Scholar] [CrossRef]
  3. Bredehoeft, J. Safe Yield and the Water Budget Myth. Groundwater 1997, 35, 929. [Google Scholar] [CrossRef]
  4. Paniconi, C.; Putti, M. Physically based modeling in catchment hydrology at 50: Survey and outlook. Water Resour. Res. 2015, 51, 7090–7129. [Google Scholar] [CrossRef]
  5. Vassena, C.; Rienzner, M.; Ponzini, G.; Giudici, M.; Gandolfi, C.; Durante, C.; Agostani, D. Modeling water resources of a highly irrigated alluvial plain (Italy): Calibrating soil and groundwater models. Hydrogeol. J. 2012, 20, 449–467. [Google Scholar] [CrossRef]
  6. Martin-Ortega, J.; Ferrier, R.C.; Gordon, I.J.; Khan, S. Water Ecosystem Services: A Global Perspective; UNESCO: Paris, France; Cambridge University Press: Cambridge, UK, 2015; Available online: https://unesdoc.unesco.org/ark:/48223/pf0000244743.locale=en (accessed on 5 November 2023).
  7. Nielsen, D.R.; Van Genuchten, M.Th.; Biggar, J.W. Water flow and solute transport processes in the unsaturated zone. Water Resour. Res. 1986, 22, 89S–108S. [Google Scholar] [CrossRef]
  8. van Genuchten, M.T.; Jury, W.A. Progress in unsaturated flow and transport modeling. Rev. Geophys. 1987, 25, 135–140. [Google Scholar] [CrossRef]
  9. Feddes, R.; Kabat, P.; Van Bakel, P.; Bronswijk, J.; Halbertsma, J. Modelling soil water dynamics in the unsaturated zone—State of the art. J. Hydrol. 1988, 100, 69–111. [Google Scholar] [CrossRef]
  10. Raats, P.A. Developments in soil–water physics since the mid 1960s. Geoderma 2001, 100, 355–387. [Google Scholar] [CrossRef]
  11. Diamantopoulos, E.; Durner, W. Dynamic Nonequilibrium of Water Flow in Porous Media: A Review. Vadose Zone J. 2012, 11, vzj2011.0197. [Google Scholar] [CrossRef]
  12. Subbaiah, R. A review of models for predicting soil water dynamics during trickle irrigation. Irrig. Sci. 2013, 31, 225–258. [Google Scholar] [CrossRef]
  13. Vereecken, H.; Schnepf, A.; Hopmans, J.; Javaux, M.; Or, D.; Roose, T.; Vanderborght, J.; Young, M.; Amelung, W.; Aitkenhead, M.; et al. Modeling Soil Processes: Review, Key Challenges, and New Perspectives. Vadose Zone J. 2016, 15, vzj2015.09.0131. [Google Scholar] [CrossRef]
  14. Reeves, M.; Baker, N.A.; Duguid, J.O. Review and Selection of Unsaturated Flow Models; TRW Environmental Safety Systems Inc.: Las Vegas, NV, USA; INTERA Inc.: Las Vegas, NV, USA, 1994. [Google Scholar] [CrossRef]
  15. Zha, Y.; Yang, J.; Zeng, J.; Tso, C.H.M.; Zeng, W.; Shi, L. Review of numerical solution of Richardson–Richards equation for variably saturated flow in soils. WIREs Water 2019, 6, e1364. [Google Scholar] [CrossRef]
  16. Bear, J. Dynamics of Fluids in Porous Media; American Elsevier Publishing Company: Amsterdam, The Netherlands, 1972. [Google Scholar]
  17. Carrillo, F.J.; Bourg, I.C. A Darcy-Brinkman-Biot Approach to Modeling the Hydrology and Mechanics of Porous Media Containing Macropores and Deformable Microporous Regions. Water Resour. Res. 2019, 55, 8096–8121. [Google Scholar] [CrossRef]
  18. Hassanizadeh, M.; Gray, W.G. General conservation equations for multi-phase systems: 1. Averaging procedure. Adv. Water Resour. 1979, 2, 131–144. [Google Scholar] [CrossRef]
  19. Nimmo, J.R. The processes of preferential flow in the unsaturated zone. Soil Sci. Soc. Am. J. 2021, 85, 1–27. [Google Scholar] [CrossRef]
  20. Fokker, A.D. Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld. Ann. Der Phys. 1914, 348, 810–820. [Google Scholar] [CrossRef]
  21. Planck, M. Über einen Satz der statistischen Dynamik und seine Erweiterung in der Quantentheorie. Sitzungsberichte Der K. Preuss. Akad. Der Wiss. Berl. 1917, 24, 324–341. [Google Scholar]
  22. Philip, J.R. The theory of infiltration: 1. The infiltration equation and its solutions. Soil Sci. 1957, 83, 345–358. [Google Scholar] [CrossRef]
  23. Richards, L.A. Capillary conduction of liquids through porous mediums. Physics 1931, 1, 318–333. [Google Scholar] [CrossRef]
  24. Raats, P.A.C. Steady Infiltration from Line Sources and Furrows. Soil Sci. Soc. Am. J. 1970, 34, 709–714. [Google Scholar] [CrossRef]
  25. Braester, C. Moisture variation at the soil surface and the advance of the wetting front during infiltration at constant flux. Water Resour. Res. 1973, 9, 687–694. [Google Scholar] [CrossRef]
  26. Warrick, A.W.; Amoozegar-Fard, A. Soil water regimes near porous cup water samplers. Water Resour. Res. 1977, 13, 203–207. [Google Scholar] [CrossRef]
  27. Pullan, A.J. The quasilinear approximation for unsaturated porous media flow. Water Resour. Res. 1990, 26, 1219–1234. [Google Scholar] [CrossRef]
  28. van Lier, Q.d.J.; Metselaar, K.; van Dam, J.C. Root Water Extraction and Limiting Soil Hydraulic Conditions Estimated by Numerical Simulation. Vadose Zone J. 2006, 5, 1264–1277. [Google Scholar] [CrossRef]
  29. Kirchhoff, G. Vorlesungen über Mathematische Physik: Bd. Vorlesungen über die Theorie der Wärme; B.G. Teubner: Berlin, Germany, 1894; Volume 4. [Google Scholar]
  30. Philip, J.R. Steady Infiltration From Buried Point Sources and Spherical Cavities. Water Resour. Res. 1968, 4, 1039–1047. [Google Scholar] [CrossRef]
  31. Wooding, R.A. Steady Infiltration from a Shallow Circular Pond. Water Resour. Res. 1968, 4, 1259–1273. [Google Scholar] [CrossRef]
  32. Davidson, M.R. Numerical calculation of saturated-unsaturated infiltration in a cracked soil. Water Resour. Res. 1985, 21, 709–714. [Google Scholar] [CrossRef]
  33. Batu, V. Time-dependent linearized two-dimensional infiltration and evaporation from nonuniform and periodic strip sources. Water Resour. Res. 1983, 19, 1523–1529. [Google Scholar] [CrossRef]
  34. Srivastava, R.; Yeh, T.C.J. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resour. Res. 1991, 27, 753–762. [Google Scholar] [CrossRef]
  35. Shan, C.; Stephens, D.B. Steady Infiltration Into a Two-Layered Soil from a Circular Source. Water Resour. Res. 1995, 31, 1945–1952. [Google Scholar] [CrossRef]
  36. Ursino, N. Linear Stability Analysis of Infiltration, Analytical and Numerical Solution. Transp. Porous Media 2000, 38, 261–271. [Google Scholar] [CrossRef]
  37. Romano, N.; Brunone, B.; Santini, A. Numerical analysis of one-dimensional unsaturated flow in layered soils. Adv. Water Resour. 1998, 21, 315–324. [Google Scholar] [CrossRef]
  38. Freeze, R.; Cherry, J. Groundwater; Prentice-Hall: Englewood Cliffs, NJ, USA, 1979; ISBN 0-13-365312-9. [Google Scholar]
  39. Hillel, D. Environmental Soil Physics: Fundamentals, Applications, and Environmental Considerations; Elsevier Science: Amsterdam, The Netherlands, 1998. [Google Scholar]
  40. Hohenbrink, T.L.; Jackisch, C.; Durner, W.; Germer, K.; Iden, S.C.; Kreiselmeier, J.; Leuther, F.; Metzger, J.C.; Naseri, M.; Peters, A. Soil water retention and hydraulic conductivity measured in a wide saturation range. Earth Syst. Sci. Data 2023, 15, 4417–4432. [Google Scholar] [CrossRef]
  41. Nemes, A.; Schaap, M.; Leij, F.; Wösten, J. Description of the unsaturated soil hydraulic database UNSODA version 2.0. J. Hydrol. 2001, 251, 151–162. [Google Scholar] [CrossRef]
  42. Peters, A. Simple consistent models for water retention and hydraulic conductivity in the complete moisture range. Water Resour. Res. 2013, 49, 6765–6780. [Google Scholar] [CrossRef]
  43. Luo, Z.; Kong, J.; Shen, C.; Lu, C.; Hua, G.; Zhao, Z.; Zhao, H.; Li, L. Evaluation and application of the modified van Genuchten function for unsaturated porous media. J. Hydrol. 2019, 571, 279–287. [Google Scholar] [CrossRef]
  44. Mantovani, F.; Menziani, M.; Pugnaghi, S.; Rivasi, M.R.; Santangelo, R. Approximating functions in hydrology of unsaturated soils—I. Nuovo C. C 1988, 11, 527–535. [Google Scholar] [CrossRef]
  45. Giudici, M.; Menziani, M.; Pugnaghi, S.; Rivasi, M.R.; Santangelo, R. Approximating functions in hydrology of unsaturated soils—II. Nuovo C. C 1989, 12, 781–786. [Google Scholar] [CrossRef]
  46. Giudici, M.; Menziani, M.; Rivasi, M.R. Determination of soil water content by mathematical models. In Proceedings of the Monitorare l’Ambiente Agrario e Forestale, Porto Conte, Italy, 4–6 June 1991; Accademia dei Georgofili: Florence, Italy, 1991; pp. 769–782. [Google Scholar]
  47. Gilding, B.H. The occurrence of interfaces in nonlinear diffusion-advection processes. Arch. Ration. Mech. Anal. 1988, 100, 243–263. [Google Scholar] [CrossRef]
  48. Gilding, B.H. Qualitative mathematical analysis of the Richards equation. Transp. Porous Media 1991, 6, 651–666. [Google Scholar] [CrossRef]
  49. Witelsky, T.P. Perturbation Analysis for Wetting Fronts in Richard’s Equation. Transp. Porous Media 1997, 27, 121–134. [Google Scholar] [CrossRef]
  50. Morelli, S.; Santangelo, R.; Vincenzi, S. Confined diffusion. Nuovo C. C 1991, 14, 377–390. [Google Scholar] [CrossRef]
  51. Kool, J.; Parker, J.; van Genuchten, M. Parameter estimation for unsaturated flow and transport models—A review. J. Hydrol. 1987, 91, 255–293. [Google Scholar] [CrossRef]
  52. Hopmans, J.W.; Šimůnek, J.; Romano, N.; Durner, W. 3.6.2. Inverse Methods. In Methods of Soil Analysis; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2002; pp. 963–1008. [Google Scholar] [CrossRef]
  53. Durner, W.; Lipsius, K. Determining Soil Hydraulic Properties. In Encyclopedia of Hydrological Sciences; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2006; tg 4 r 75. [Google Scholar] [CrossRef]
  54. Giudici, M.; Baratelli, F.; Cattaneo, L.; Comunian, A.; Filippis, G.D.; Durante, C.; Giacobbo, F.; Inzoli, S.; Mele, M.; Vassena, C. A conceptual framework for discrete inverse problems in geophysics. arXiv 2021, arXiv:1901.07937. [Google Scholar] [CrossRef]
  55. Dai, Z.; Samper, J. Inverse problem of multicomponent reactive chemical transport in porous media: Formulation and applications. Water Resour. Res. 2004, 40, W07407. [Google Scholar] [CrossRef]
  56. Zheng, L.; Samper, J. Formulation of the inverse problem of non-isothermal multiphase flow and reactive transport in porous media. In Computational Methods in Water Resources: Volume 2; Developments in Water Science; Miller, C.T., Pinder, G.F., Eds.; Elsevier: Amsterdam, The Netherlands, 2004; Volume 55, pp. 1317–1327. [Google Scholar] [CrossRef]
  57. Brooks, R.; Corey, A. Hydraulic Properties of Porous Media; Hydrology Paper No. 3; Civil Engineering Department, Colorado State University: Fort Collins, CO, USA, 1964. [Google Scholar]
  58. van Genuchten, M.T. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef]
  59. Kozeny, J. Über kapillare Leitung des Wassers im Boden. Sitzungsberichte Der Kais. Akad. Der Wiss. Wien 1927, 136, 271–306. [Google Scholar]
  60. Carman, P.C. Fluid flow through a granular bed. Trans. Inst. Chem. Eng. 1937, 15, 150–156. [Google Scholar] [CrossRef]
  61. Burdine, N. Relative Permeability Calculations From Pore Size Distribution Data. J. Pet. Technol. 1953, 5, 71–78. [Google Scholar] [CrossRef]
  62. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 1976, 12, 513–522. [Google Scholar] [CrossRef]
  63. Leij, F.J.; Russell, W.B.; Lesch, S.M. Closed-Form Expressions for Water Retention and Conductivity Data. Groundwater 1997, 35, 848–858. [Google Scholar] [CrossRef]
  64. Facchi, A.; Ortuani, B.; Maggi, D.; Gandolfi, C. Coupled SVAT—Groundwater model for water resources simulation in irrigated alluvial plains. Environ. Model. Softw. 2004, 19, 1053–1063. [Google Scholar] [CrossRef]
  65. Romano, E.; Giudici, M. Experimental and modeling study of the soil-atmosphere interaction and unsaturated water flow to estimate the recharge of a phreatic aquifer. ASCE J. Hydrol. Eng. 2007, 12, 573–584. [Google Scholar] [CrossRef]
  66. Romano, E.; Giudici, M. On the use of meteorological data to assess the evaporation from a bare soil. J. Hydrol. 2009, 372, 30–40. [Google Scholar] [CrossRef]
  67. Comunian, A.; Giudici, M.; Landoni, L.; Pugnaghi, S. Improving Bowen-ratio estimates of evaporation using a rejection criterion and multiple-point statistics. J. Hydrol. 2018, 563, 43–50. [Google Scholar] [CrossRef]
  68. Liu, L.; Dai, S.; Ning, F.; Cai, J.; Liu, C.; Wu, N. Fractal characteristics of unsaturated sands—Implications to relative permeability in hydrate-bearing sediments. J. Nat. Gas Sci. Eng. 2019, 66, 11–17. [Google Scholar] [CrossRef]
  69. Chen, H.; Yang, M.; Chen, K.; Zhang, C. Relative Permeability of Porous Media with Nonuniform Pores. Geofluids 2020, 2020, 5705424. [Google Scholar] [CrossRef]
  70. Chen, H.; Chen, K.; Yang, M.; Xu, P. A fractal capillary model for multiphase flow in porous media with hysteresis effect. Int. J. Multiph. Flow 2020, 125, 103208. [Google Scholar] [CrossRef]
  71. Solazzi, S.G.; Rubino, J.G.; Jougnot, D.; Holliger, K. Dynamic permeability functions for partially saturated porous media. Geophys. J. Int. 2020, 221, 1182–1189. [Google Scholar] [CrossRef]
  72. Nguyen, V.N.A.; Jougnot, D.; Luong, D.T.; Phan, V.D.; Tran, T.C.T.; Dang, T.M.H.; Nguyen, M.H. Predicting water flow in fully and partially saturated porous media: A new fractal-based permeability model. Hydrogeol. J. 2021, 29, 2017–2031. [Google Scholar] [CrossRef]
  73. Carrillo, F.J.; Bourg, I.C. Capillary and viscous fracturing during drainage in porous media. Phys. Rev. E 2021, 103, 063106. [Google Scholar] [CrossRef] [PubMed]
  74. Stults, J.; Illangasekare, T.; Higgins, C.P. The Mass Transfer Index (MTI): A semi-empirical approach for quantifying transport of solutes in variably saturated porous media. J. Contam. Hydrol. 2021, 242, 103842. [Google Scholar] [CrossRef] [PubMed]
  75. Ladd, A.J.; Szymczak, P. Reactive Flows in Porous Media: Challenges in Theoretical and Numerical Methods. Annu. Rev. Chem. Biomol. Eng. 2021, 12, 543–571. [Google Scholar] [CrossRef] [PubMed]
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.
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.
FunctionProperties
h ^ ( θ ) h ^ ( θ ) h e < 0 , for θ [ θ irr , θ sat ] ; h ^ ( θ ) , as θ θ irr .
h ^ ( θ ) h ^ ( θ ) > 0 , for θ [ θ irr , θ sat ] ; h ^ ( θ ) + , as θ θ irr ; h ^ ( θ ) , as θ θ sat .
h ^ ( θ ) h ^ ( θ ) , as θ θ irr ; h ^ ( θ ) + , as θ θ sat ; h ^ ( θ i . p . ) = 0 , for some θ i . p . ( θ irr , θ sat ) .
K ^ ( θ ) K ^ ( θ ) > 0 , for θ ( θ irr , θ sat ] ; K ^ ( θ irr ) = 0 ; K ^ ( θ sat ) = K sat ; K ^ ( θ ) m ( θ θ irr ) for some m > 0 and θ θ irr .
K ^ ( θ ) K ^ ( θ ) > 0 , for θ [ θ irr , θ sat ] .
D ^ ( θ ) D ^ ( θ ) > 0 , for θ ( θ irr , θ sat ] ; ( θ θ irr ) 1 D ^ ( θ ) be integrable for θ [ θ irr , θ sat ] .
D ^ ( θ ) D ^ ( θ ) 0 , for θ ( θ irr , θ sat ] ; D ^ ( θ irr ) = 0 ; D ^ ( θ ) + , as θ θ sat .
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Giudici, M. Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments. Sustainability 2023, 15, 15723. https://doi.org/10.3390/su152215723

AMA Style

Giudici M. Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments. Sustainability. 2023; 15(22):15723. https://doi.org/10.3390/su152215723

Chicago/Turabian Style

Giudici, Mauro. 2023. "Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments" Sustainability 15, no. 22: 15723. https://doi.org/10.3390/su152215723

APA Style

Giudici, M. (2023). Modeling Water Flow in Variably Saturated Porous Soils and Alluvial Sediments. Sustainability, 15(22), 15723. https://doi.org/10.3390/su152215723

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop