Freely available online through the author-supported open access option.

Abstract

Optimality principles have been used to investigate physical processes in different areas. This work applied an optimal principle (that water flow resistance is minimized for the entire flow domain) to steady-state unsaturated flow processes. Based on the calculus of variations, under optimal conditions, hydraulic conductivity for steady-state, gravity-dominated unsaturated flow is proportional to a power function of the magnitude of water flux. This relationship is consistent with an intuitive expectation that for an optimal water flow system, locations where relatively large water fluxes occur should correspond to relatively small resistance (or large conductance). This theoretical result was also consistent with observed fingering-flow behavior in unsaturated soils and an existing model.

This work presents a new theory that hydraulic conductivity for steady-state unsaturated flow is proportional to a power function of water flux when water flow resistance is minimizedfor the entire flow domain. This result was also demonstrated to be consistent with observed fingering-flow behavior.

Optimality principles refer to that state of a physical process that is controlled by an optimal condition that is subject to physical or resource constraints. For example, Eagleson (2002) demonstrated that under natural conditions and in water-limited areas, vegetation tends to grow under maximum-productivity and unstressed conditions. He called the function and forms of vegetation, following the optimality principle, the results of “Darwinian expression.” After studying a variety of natural phenomena characterized by tree-like structures, Bejan (2000) proposed a “constructal theory,” which states that “for a finite-size open system to persist in time (to survive) it must evolve in such a way that it proves easier and easier access to the currents that flow through it.” While the definition of “easy access” is not always clear, Bejan (2000) demonstrated that tree-like structures are the direct results of the minimization of flow resistance across whole flow systems under consideration. During the past 30 yr, the maximum entropy production (MEP) principle has been successfully applied, in a heuristic sense, to the prediction of steady states of a wide range of systems (Niven, 2010; Kleidon, 2009). The MEP principle states that a flow system subject to various flows or gradients will tend toward a steady-state position of maximum thermodynamic entropy production (Niven, 2010). The theoretical connections between these optimality principles and the currently existing fundamental laws, however, are not well established. The alternative point of view is that these principles are actually self-standing and do not follow from other known laws (Bejan, 2000).

The role of optimality principles in forming complex natural patterns has been recognized for many years in the surface hydrology community (Leopold and Langbein, 1962; Howard, 1990; Rodriguez-Iturbe et al., 1992; Rinaldo et al., 1992; Liu, 2010). For example, Leopold and Langbein (1962) proposed a maximum entropy principle for studying the formation of landscapes. Rodriguez-Iturbe et al. (1992) postulated principles of optimality in energy expenditure at both local and global scales for channel networks. While previous studies mainly use spatially “discrete” approaches as a result of considering energy dissipation through channel networks only, Liu (2010) developed a group of (partial differential) governing equations for steady-state optimal landscapes (including both channel networks and associated hillslopes) using the calculus of variations.

The importance of optimality principles has also been intuitively recognized in the vadose zone hydrology community for a long time. For example, it seems to be well known that fingering flow is due to the fact that unsaturated water tends to form flow paths corresponding to the minimized flow resistances. (Note that for a given water flux, fingering flow always gives lower flow resistance [or higher conductance] compared with uniform flow because fingering flow paths generally have higher local water saturations that correspond to larger unsaturated conductivities.) Rigorous applications of this optimality principle have not been fully explored, however. Because of fingering flow, water propagates quickly to significant depths while bypassing large portions of the vadose zone, and solute travel times from a contamination source (located on the soil surface or in the vadose zone) to the groundwater are shorter than a priori expected. As a result of the important effects of this flow process on groundwater contamination (an important issue for water resources management), preferential flow has been a major research area in the vadose zone hydrology community for a number of years and considered probably the most frustrating processes in terms of hampering accurate predictions of contaminant transport in the vadose zone (e.g., Glass et al., 1988; Flury and Flühler, 1995; Liu et al., 2003; Šimůnek et al., 2003; Nimmo, 2010).

This study developed a conductivity relationship for gravity-dominated unsaturated flow derived from a principle that the energy dissipation rate (or flow resistance) is minimized for the entire flow system. Preliminary evaluation of this relationship was conducted by comparing it with relevant experimental observations and the currently existing models. The potential limitations and further improvements of this work are also briefly discussed.

Theory

As the first step, we consider a relatively simple, steady-state unsaturated flow system associated with a homogeneous and isotropic porous medium. From water mass (volume) conservation, the steady-state water flow equation is given by 
\[\frac{{\partial}q_{x}}{{\partial}x}{+}\frac{{\partial}q_{y}}{{\partial}y}{+}\frac{{\partial}q_{z}}{{\partial}z}{=}0\]
[1]
where x and y are two horizontal coordinate axes, z is the vertical axis, and qx, qy, and qz are volumetric fluxes of water along the x, y, and z directions, respectively. (Water flux is used here to mean the volumetric flux of water.)
The parameter E (a function of x, y, and z) represents the total energy, including both potential (corresponding to elevation z) and (capillary) pressure energy: 
\[E{=}z{+}\frac{P}{{\rho}g}{=}z{+}h\]
[2]
where g is gravitational acceleration, P is capillary pressure, ρ is water density, and h is capillary pressure head. Accordingly, the energy expenditure rate for a unit control volume, ΔE, can be expressed as 
\[\mathrm{{\Delta}}E{=}\frac{{\partial}\left(q_{x}E\right)}{{\partial}x}{+}\frac{{\partial}\left(q_{y}E\right)}{{\partial}y}{+}\frac{{\partial}\left(q_{z}E\right)}{{\partial}z}\]
[3]
This equation simply states that for a given unit volume, the energy expenditure rate at that location is equal to the energy carried by the water flowing into the volume minus the energy carried by the water flowing out of the volume.
A combination of Eq. [1] and [3] yields 
\[\mathrm{{\Delta}}E{=}q_{x}\frac{{\partial}E}{{\partial}x}{+}q_{y}\frac{{\partial}E}{{\partial}y}{+}q_{z}\frac{{\partial}E}{{\partial}z}\]
[4]
Throughout this development, Darcy's law is assumed to applied to unsaturated flow: 
\[q_{x}{=}{-}K\frac{{\partial}E}{{\partial}x}\]
[5a]
 
\[q_{y}{=}{-}K\frac{{\partial}E}{{\partial}y}\]
[5b]
 
\[q_{z}{=}{-}K\frac{{\partial}E}{{\partial}z}\]
[5c]
where K is hydraulic conductivity and is given by 
\[K{=}K\left(h\mathrm{,}S\right)\]
[5d]
 
\[S{=}\left(\frac{{\partial}E}{{\partial}x}\right)^{2}{+}\left(\frac{{\partial}E}{{\partial}y}\right)^{2}{+}\left(\frac{{\partial}E}{{\partial}z}\right)^{2}\]
[5e]
In Eq. [5d], hydraulic conductivity is assumed to be a function of both the capillary pressure head (h) and the square of the energy gradient (S). A previous study for an optimal landscape indicated that water-flow conductance is a function of the water flux in the system (Liu, 2010). That result is mathematically identical to the empirical relations between water depth and surface slope observed from many river basins (e.g., Leopold and Langbein, 1962). Assuming K to be a function of water flux is equivalent to assuming it to be a function of the energy gradient because water flux, energy gradient, and K are related through Darcy's law. The function form of K (Eq. [5d]) was a subject of study in this work. Note that our theory was developed for a macroscopic scale that may include a number of fingering or preferential flow paths. Local scale refers to the continuum scale within each finger. The unsaturated flow process at the local scale is mainly controlled by pore-scale physics. We also need to emphasize that assuming K to be a function of S at the macroscopic scale does not exclude the possibility that K may have nothing to do with S in our final results. In this case, the mathematically derived function form of K would not include S as an independent variable.
When we combine Eq. [4] and [5], the global energy expenditure rate through domain Ω is given by 
\[{\int}{{\int}_{\mathrm{{\Omega}}}}\mathrm{{\Delta}}Edxdydz{=}{\int}{{\int}_{\mathrm{{\Omega}}}}\left({-}KS\right)dxdydz\]
[6]
The optimality principle in our problem is to minimize the absolute value of the above integral. To do so, we use the calculus of variations that seeks optimal (stationary) solutions to a functional (a function of functions) by identifying unknown functions (Weinstock, 1974).
Based on Eq. [5] and [6], the Lagrangian for the given problem is given by 
\[L{=}{-}KS{+}{\lambda}_{1}\left[S{-}\left(\frac{{\partial}E}{{\partial}x}\right)^{2}{-}\left(\frac{{\partial}E}{{\partial}y}\right)^{2}{-}\left(\frac{{\partial}E}{{\partial}z}\right)^{2}\right]\]
[7]
Note that the first term is from Eq. [6] and the second term is a constraint from Eq. [5]. The use of the constraint term allows consideration of related functions to be independent when the optimal solution to Eq. [6] is determined. The l function is the Lagrangian multiplier. A mathematically equivalent way to define L to avoid the use of some (or all) constraints is to directly insert Eq. [5] into the first term of Eq. [7] (Pike, 2001). In this case, the number of independent functions will be reduced; however, the use of Eq. [7] is more straightforward and easier to handle for the given problem. Also, note that the other constraint should be the continuity equation and is dealt with below.
The following Euler–Lagrangian equation is used to determine an unknown function w associated with L to minimize the integral defined in Eq. [6] (Weinstock, 1974): 
\[\frac{{\partial}L}{{\partial}w}{-}\frac{{\partial}}{{\partial}x}\left(\frac{{\partial}L}{{\partial}w_{x}}\right){-}\frac{{\partial}}{{\partial}y}\left(\frac{{\partial}L}{{\partial}w_{y}}\right){-}\frac{{\partial}}{{\partial}z}\left(\frac{{\partial}L}{{\partial}w_{z}}\right){=}0\]
[8]
where wx, wy, and wz are partial derivatives with respect to x, y, and z, respectively. In this study, w corresponds to S and h (or E). (Also note that application of the Euler–Lagrangian equation to Lagrangian multipliers will recover Eq. [5e]).
Replacing w with S in Eq. [8] yields 
\[{\lambda}_{1}{=}\frac{{\partial}\left(KS\right)}{{\partial}S}\]
[9]
Replacing w with h (or E) in Eq. [8] and using Eq. [9] and the continuity equation, we have 
\[\frac{{\partial}\left[\left({\partial}K\mathrm{/}{\partial}\mathrm{log}S\right)\left({\partial}E\mathrm{/}{\partial}x\right)\right]}{{\partial}x}{+}\frac{{\partial}\left[\left({\partial}K\mathrm{/}{\partial}\mathrm{log}S\right)\left({\partial}E\mathrm{/}{\partial}y\right)\right]}{{\partial}y}{+}\frac{{\partial}\left[\left({\partial}K\mathrm{/}{\partial}\mathrm{log}S\right)\left({\partial}E\mathrm{/}{\partial}z\right)\right]}{{\partial}z}{=}\frac{S}{2}\frac{{\partial}K}{{\partial}h}\]
[10]
In general, it is difficult to obtain an analytical solution to the above equation. For some special case in which the term on the right-hand side is small compared with the other terms, however, a closed-form solution can be obtained. This may be true for a gravity-dominant flow based on a rough order-of-magnitude analysis [S(∂K/∂h)]/(∂K/∂logS) = ∂S/∂h. Note that for gravity-dominant flow, the energy gradient is close to one in the vertical direction and not a strong function of local capillary pressure at locations occupied by flow patterns. In other words, for two different gravity-dominated flow situations for a given flow system, changes in the energy gradient are much smaller than changes in the capillary pressure. The argument is also consistent with an observation that ∂K/∂h → 0 for large water saturation (corresponding to the gravity-dominated conditions) (van Genuchten, 1980). A comparison between Eq. [10] (without the term on the right-hand side) with the continuity equation (Eq. [1] and [5]) yields 
\[\frac{{\partial}K}{{\partial}\mathrm{log}S}{=}AK\]
[11]
where A is a constant.
To get practically useful results, we consider K(h,S) to be further expressed by 
\[K\left(h\mathrm{,}S\right){=}f\left(h\right)g\left(S\right)\]
[12]
Substituting Eq. [12] into [11] results in 
\[g\left(S\right){\propto}S^{A}\]
[13]
Based on Darcy's law, Eq. [13] can be rewritten as 
\[g\left(S\right){\propto}\left(\frac{\left|q\right|}{K}\right)^{A\mathrm{/}2}\]
[14]
where |q| is the magnitude of water flux given by 
\[\left|q\right|{=}\left[q_{x}^{2}{+}q_{y}^{2}{+}q_{z}^{2}\right]^{1\mathrm{/}2}\]
[15]
Combining Eq. [14] and [12] gives our final conductivity relationship: 
\[K{=}F\left(h\right)\left(\frac{\left|q\right|}{K_{\mathrm{sat}}}\right)^{a}\]
[16]
where a = A/(2 + A) and Ksat is the saturated hydraulic conductivity. It is very interesting to note that although the mathematical derivation processes are considerably complex and involve solving a group of partial differential equations, the final result (Eq. [16]) is amazingly simple for gravity-dominant flow problems.

Discussion

Under optimal flow conditions corresponding to the minimum energy dissipation rate (or flow resistance), the derived conductivity is a power function of water flux (Eq. [16]) for gravity-dominated unsaturated flow. This result physically makes sense. For the positive power values, the smallest flow resistance occurs within flow paths with the largest water flux. Intuitively, it is easy to understand that this conductivity distribution will result in minimized total flow resistance globally. This finding is also consistent with our daily life experiences. For example, to maximize traffic transportation efficiency, our highways always have more lanes (or higher “conductance”) in locations with high traffic fluxes. (Highway networks may be considered to be analogous to fingering flow paths.)

There may be different interpretations of Eq. [16]. One interpretation is that F(h) is the local-scale hydraulic conductivity within the fingering flow zone and that the power function of flux in the equation represents the fraction of the fingering flow zone in an area normal to the water flux direction. This is justified because h is a local-scale variable. In this case, our result is supported by the analysis results of Wang et al. (1998). On the basis of a number of laboratory experimental observations of vertical fingering flow in homogeneous soils, Wang et al. (1998) presented a relation between flow conditions and a parameter, Fa, defined as the ratio of the horizontal cross-sectional area occupied by gravity fingers to the total cross-sectional areas:  
\[F_{\mathrm{a}}{=}\left[\frac{\left|q\right|}{K_{\mathrm{sat}}}\right]^{0.5}\]
[17]
Obviously, Eq. [17] is identical to our theoretical result with a = 0.5. In other words, our theoretical result agrees with the laboratory observations cited by Wang et al. (1998).
Our theoretical result is also consistent with the active region model (ARM) proposed by Liu et al. (2005) for dealing with unsaturated flow in soils. The ARM is an extension of the active fracture model developed for modeling unsaturated water in fractured rock (Liu et al., 1998). Both the active fracture model and the ARM have been evaluated with a variety of experimental data and remarkable agreements between the models and the data have been observed (Liu et al., 1998, 2003, 2005; Sheng et al., 2009); however, the ARM has been tested only for soils that are relatively homogeneous (e.g., Sheng et al., 2009). The ARM assumes a flow domain to be divided into an active region (fingering flow zone) and an inactive region. Flow occurs only in the active region. The volumetric portion of the active region is given as 
\[f{=}\mathrm{{\theta}}^{{\gamma}}\]
[18]
where θ is the average effective water saturation across the whole flow domain (including both active and inactive regions), and γ is a constant factor between zero and one. Note that f is equivalent to Fa in Eq. [17].
By definition, the average water saturation is related to the effective water saturation (θa) within the active region by 
\[\mathrm{{\theta}}{=}f\mathrm{{\theta}}_{\mathrm{a}}\]
[19]
For gravity-dominated flow, the energy gradient approximately equals one and the vertical water flux is about the same as the hydraulic conductivity. Using the well-known Brooks–Corey relationship (Brooks and Corey, 1964) to describe the hydraulic conductivity within the active region, we can write the total vertical water flux as 
\[\frac{q}{K_{\mathrm{sat}}}{=}\mathrm{{\theta}}_{\mathrm{a}}^{\mathrm{{\beta}}}f\]
[20]
where β is the Brooks–Corey exponent. Combining Eq. [18–20] yields 
\[f{=}\left(\frac{\left|q\right|}{K_{\mathrm{sat}}}\right)^{\frac{1}{1{+}\left[\mathrm{{\beta}}\left(1{-}{\gamma}\right)\mathrm{/}{\gamma}\right]}}\]
[21]
Thus, Eq. [21] derived from the ARM is equivalent to our conductivity relationship with 
\[a{=}\frac{1}{1{+}\left[\mathrm{{\beta}}\left(1{-}{\gamma}\right)\mathrm{/}{\gamma}\right]}\]
[22]
In the other words, we demonstrate the equivalence between our Eq. [16] and the ARM for gravity-dominated unsaturated flow under the condition that the power function in Eq. [16] is interpreted as a volumetric fraction of the fingering flow zones within the soil.

It is of interest to note that for typical values of β = 4 and γ = 0.7 (Brooks and Corey, 1964; Sheng et al., 2009), a = 0.4 is close to the value of 0.50 given in Eq. [17]. Whether or not a single value for parameter a is valid for different soils needs further research based on experimental observations.

Finally, this work is the first step to incorporate the optimality principle into unsaturated flow. Consequently, some limitation of the current work still exists. For example, detailed pore-scale unsaturated-flow physics is not adequately incorporated yet. This physics requires that the upper limit of K in Eq. [16] should be F(h), which, however, is not reflected in our theory. This can be approximately accounted for in practice by limiting the K value calculated from Eq. [16] to the corresponding F(h) value. Nevertheless, the major focus of this study was to highlight the potential for developing new unsaturated water flow theories based on the optimal principle. This principle may hold the key to resolving a number of problems associated with emerging patterns in unsaturated soils.

Conclusions

Based on the calculus of variations, this work showed that under optimal conditions, hydraulic conductivity for steady-state, gravity-dominated unsaturated flow is proportional to a power function of the magnitude of water flux. It is consistent with an intuitive expectation that for an optimal water flow system, locations where relatively large water fluxes occur should correspond to relatively small resistance (or large conductance). Consistence between this theoretical result with observed fingering flow behavior in unsaturated soils and the ARM was also demonstrated. Finally, it is important to note that the classic unsaturated-flow theory is applicable to capillarity-dominated cases while the current work focused on unsaturated flow under gravity-dominated conditions. Whether the optimality principle can be used to develop a general theory (that includes both cases as two special ones) deserves further research.

The initial version of the paper was carefully reviewed by Dr. Jim Houseworth and Dr. Dan Hawkes. The work was performed under DOE Contract DE-AC03–76F00098.