Averaged water potentials in soil water and groundwater , and their connection to menisci in soil pores , field-scale flow phenomena , and simple groundwater flows

Abstract. The movement of subsurface water is mostly studied at the pore scale and the Darcian scale, but the field and regional scales are of much larger societal interest. Volume-averaging has provided equations at these larger scales, but the required restrictions rendered them of little practical interest. Others hypothesized a direct connection at hydrostatic equilibrium between the average matric potential of a subsurface body of water and the average pressure drop over the menisci in the soil pores. The link between the volume-averaged potential energy of subsurface water bodies and large-scale fluxes remains largely unexplored. This paper treats the effect of menisci on the potential energy of the water behind them in some detail, and discusses some field-scale effects of pore-scale processes. Then, various published expressions for volume-averaged subsurface water potentials are compared. The intrinsic phase average is deemed the best choice. The hypothesized relationship between average matric potential and average meniscus curvature is found to be valid for unit gradient flow instead of hydrostatic equilibrium. Still, this restriction makes the relationship hold only for a specific depth range in the unsaturated zone under specific conditions, and certainly not for entire fields or catchments. In the groundwater, volume-averaged potential energy is of more use: for linearized, steady flows with flow lines that are parallel, radially diverging, and radially converging, proofs are derived for proportionality between averaged hydraulic potentials and fluxes towards open water at a fixed potential. For parallel flow, a simplified but relevant transient flow case also exhibits this proportionality.


Introduction
In the saturated zone of the subsurface, groundwater fills the entire pore space.If the solid phase is assumed rigid and unchanging, the liquid-solid interface is static.The single fluid flow can be simplified to Darcian flow for most practical situations, and the continuum approach (Bear and Bachmat, 1991, 14-42 pp.) facilitates the quantitative analysis without having to resort to pore-scale details.The potential energy of water has a gravitational and a pressure component which values and gradients are well defined at the scale of the representative elementary volume (Darcian scale).
In the unsaturated zone, the introduction of the gas phase creates solid-liquid, solid-gas, and liquid-gas interfaces, the latter of which are mobile (if we assume the solid phase to be rigid and its architecture time invariant), while the total area of the other two can vary.Natural conditions often lead to situations where the liquid phase in an extended area (a field of several hectares or a catchment of several square kilometers) is confined behind its solid-liquid and liquid-gas interfaces, possibly with enclosed pockets of the gas phase trapped within it.The mobility of the gas-liquid interface nevertheless allows ample movement of and exchange between subsurface and atmospheric water through infiltration and evapotranspiration.Still, the bulk of the gas phase is above the bulk of the liquid phase, and most of the gas phase is generally well-connected to the atmosphere.The continuum approach can be applied in the unsaturated zone, and the potential energies of water (gravitational and matric potential) again are well-defined at the Darcian scale, but phenomena such as hysteresis and sharp wetting fronts in the unsaturated zone require an understanding of the behavior of the liquid-gas interfaces.
G. H. de Rooij: Averaged water potentials in soil water and groundwater Hassanizadeh and Gray (1990) provided a comprehensive thermodynamic treatment of the interaction between the various phases, their interfaces, and the contact lines between the interfaces.Their analysis produces 28 equations and unknowns, and moves from the pore scale to the Darcian scale (in the order of several hundreds of pore diameters).It demonstrates that the generally accepted proportionality between the curvature of a meniscus and the pressure jump across it is only valid at equilibrium.Hassanizadeh and Gray (1990) also offered a modified version of Darcy's Law that accounts for the masses and energies of the interfaces in situations where these move and deform rapidly.Most issues of societal interest arise at the field and regional scales, but Hassanizadeh and Gray's (1990) full analysis is difficult to carry out beyond the Darcian scale.
A formal volume-averaging over a wider range of scales was attempted by various authors (e.g., Crapiste et al., 1986;Quintard andWhitaker, 1988, 1990a, b;Gray, 2002).The mathematical rigor of this approach required considerable limitations in the type of problems that could be analyzed in order to solve the closure problem (see Quintard and Whitaker, 1988), leading to such requirements as clearly separated hierarchical spatial scales of heterogeneities and absence of gravity, which severely hampered the application to realistic subsurface environments.
The influence of the atmosphere and the vegetation on the unsaturated zone creates large fluctuations of the matric potential and in the curvature of the menisci.Through this curvature, the configuration of liquid and gas in the pore space is affected.Thus, large-scale processes such as rainfall and evapotranspiration from vegetated areas interact with porescale water distribution and meniscus curvatures.Zehe et al. (2006) were among the first to attempt to quantify and interpret the potential energy of bodies of subsurface water at the catchment scale by volume-averaging local Darcianscale values.Zehe et al. (2006) only looked at the scale of the entire integration volume without relating the variations of the averaged variable within the integration volume to its volume-average, thereby avoiding the closure problem.They postulated that the volume-average of the pore water matric potential in the unsaturated zone should reflect the pressure drop over the gas-liquid interface averaged over the unsaturated zone of an entire catchment, and thus proposed a very direct link between two very disparate scales.
Both Hassanizadeh and Gray (1990) and Zehe et al. (2006) focused more on the status of the liquid phase than on its movement.De Rooij (2009) used area-averaged flux densities and hydraulic potentials to investigate scale limitations to the validity of Darcy's Law.It is worthwhile investigating if simple averaging techniques like Zehe et al.'s (2006) and de Rooij's (2009) can be expanded to flows beyond the Darcian scale.In view of Roth's (2008) finding that processes operating on small time scales are difficult to upscale, the dynamic behavior of the liquid phase in the unsaturated zone makes the saturated zone better suited to explore the link between averaged potential energies and large-scale fluxes.For the unsaturated zone, the direct connection between the mean curvature of the liquid-gas interfaces and the volume-averaged matric potential as hypothesized by Zehe et al. (2006) has not been tested yet.
The objectives of this paper address the potential of straight-forward volume-averaging as used by Zehe et al. (2006) andde Rooij (2009) for application to realistic problems in subsurface hydrology.The first objective is to examine the relationship between the volume-averaged potential energy of the soil solution and the curvatures of the gas-liquid interfaces, thereby advancing and refining the analysis by Zehe et al. (2006).In pursuing the first objective, the response of the pressure potential to addition or removal of water will be treated in some detail and some practical cases of field scale effects of the pore-scale behavior of the interfaces will be discussed.The second objective is to establish a proof of principle for a link between large scale fluxes and average potential energy of the groundwater by investigating simple saturated flow problems.

Theory
The three phases and their interactions are simplified according to Sposito (1981, 187-190 pp.), which implies that all phases are uniform in composition and do not mix or react with one another.Furthermore, we consider an isothermal system with a rigid solid phase, an inelastic liquid phase, and a continuous gas phase (no entrapped soil air unless stated otherwise).The three surface tensions are assumed constant, and flows are slow enough to make the kinematic energy negligible.The potential energy of the liquid at any location is assumed to be fully quantified by its gravitational potential and its pressure potential (below the phreatic level) or matric potential (between the soil surface and the phreatic level).The conventional definitions of these quantities are used (e.g., Jury et al., 1991, 48-52 pp.).The formal unit for potential energy is J kg −1 (potential energy per unit mass), but for clarity the potential energy per unit volume is used here (J m −3 = Pa).As Sposito (1981, 193-194 pp.) points out, these definitions are interchangeable, and the potential should not be considered a true pressure, even if it has the same dimensions.In the following, the terms "potential" and "potential energy" refer to the volume-based quantities unless stated otherwise.For actual liquid or gas pressures the term "pressure" is used; potentials expressed in units of pressure carry the qualifier "potential" (e.g., pressure potential).An important difference is the absolute nature of pressure (cannot be negative) and the relative nature of the pressure potential and matric potential, both of which equal zero when the liquid is at atmospheric pressure.
In view of the simplifying assumptions, the use of the Laplace-Young Law (e.g., Brutsaert, 2005, 255-257 pp.;Jury et al., 1991, p. 41;Or and Wraith, 2000) is Hydrol.Earth Syst.Sci., 15, 1601Sci., 15, -1614Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1601/2011/ permitted to determine the pressure jumps across the liquidgas interfaces: with P lg (Pa) the positive pressure jump from the wetting to the non-wetting phase, γ lg (N m −1 ) the liquid-gas surface tension, and r 1 and r 2 (m) the principle radii of curvature of the interface.Similarly, flow is described by the classical form of Darcy's Law.The use of these classical equations restricts the analysis at this time to problems in which the energy changes and movements of the interfaces are negligible, i.e., saturated flows or unsaturated flows in which the liquid-gas interfaces move only slowly.Note that the use of the conventional form of Darcy's Law requires this restriction to be met; the restriction thus applies to the vast majority of modeling efforts and applications published to date.
Processes such as infiltration in dry soils are not amenable to the analysis presented here.This restriction is in line with Roth's (2008) finding that phenomena with small characteristic times (e.g., infiltration fronts) lend themselves poorly to upscaling.Furthermore, we will not consider the forces exerted upon the soil solution by the porous matrix, such as Van der Waals forces and hydrogen bonds (see Or and Wraith, 2000, for more detail), other than through their combined effect as expressed in the matric potential.

Hydrostatic pressure potential below a liquid-gas interface
In soils, the capillary fringe above the phreatic level is defined here as the region where all pores other than macropores like cracks and biopores are liquid-saturated (except for possible pockets of enclosed air).Thus, the capillary fringe extends to the depth where the largest capillary pore empties and a continuous path of soil gas extends from that pore to the atmosphere.The elevation of the liquid-gas interface of that pore, and thus of the top of the capillary fringe is denoted x lg (m).For a continuous liquid phase below the top of the capillary fringe in a satiated porous medium, the pressure at hydrostatic equilibrium is given by: with x 3 (m) the vertical coordinate (positive upwards), P (Pa) the pressure in the liquid, P atm (Pa) the gas pressure above x lg , ρ (kg m −3 ) the density of the liquid phase, and g (m s −2 ) the gravitational acceleration.The sign before the last term is positive for non-wetting liquids (in which case the capillary fringe does not exist) and negative for wetting liquids.If a small quantity of water V (m 3 ) is added or removed and x lg changes as a result, we have: where A lg (m 2 ) is the liquid-occupied horizontal area at x lg .Thus, the pressure change in the liquid phase is equal to the change in gravitational potential of the water at the liquidgas interface plus its pressure change caused by changes in the curvature of the interface.The difference-form of Eq. ( 3) is generally applicable, even when the adding/removal of V involves the sudden emptying or filling of pores, e.g. by Haynes jumps.Only when the pores are tubular or their radii change gradually and monotonically, can the differential form of Eq. ( 3) can be used.The approximate equality in Eq. (3) then holds exactly.
From Eqs.
(2) and (3) it follows that the pressure at an arbitrary depth below the top of the capillary fringe can be manipulated by changing x lg (requiring adding or removing significant volumes of water), or by changing the curvature of the interface (which requires only minute volumes of water, depending on the pore architecture).An obvious third method would be to change the pressure in the gas phase above the column.
The sometimes very substantial changes in the potential energy of large bodies of liquid caused by adding or removing small quantities of liquid seem to suggest that energy is added and dissipated in quantities that are only determined by the volume of liquid experiencing these changes in matric/pressure potential, not by the possibly very much smaller volume of liquid that causes the perturbation.It should be noted that the perturbation changes the potential energy of the liquid at the interface, either by a change in the position of the interface in the gravitational field (affecting the gravitational potential), a change in the interface curvature (affecting the matric/pressure potential), or both.If the liquid below the capillary fringe is enclosed by no-flow boundaries, hydrostatic equilibrium after the perturbation can only be restored if the potential energy of the rest of the liquid follows suit.Since the location of the liquid is fixed, only the matric/pressure potential is available to accommodate the required change.If the soil solution can flow freely, the change in potential energy will generate a flow that changes the position or curvature of the meniscus.This flow will usually be in the direction that tends to reverse the effect of the initial perturbation.In most cases, the liquid displacement necessary to restore equilibrium is comparable to the liquid volume that caused the perturbation, for instance by allowing just enough soil solution to drain from a saturated soil profile to allow the liquid-gas interface to recede from the soil surface (where the interface is flat) into the pores just below the soil surface (see Sect. 3).
G. H. de Rooij: Averaged water potentials in soil water and groundwater As a practical example, relatively rapid pressure changes as a result of water extractions can occur in a confined aquifer with a small phreatic intake area, where the aquifer's only airwater interface exists.Such an aquifer can only sustain small transfers of water to a deeper aquifer before the hydraulic potential drops to equilibrium with the lower neighbor and flow ceases.If the aquifer over some extent is artesian, too many artesian wells may cause the hydraulic head to drop below the soil surface and the wells will then run dry.

Menisci in the unsaturated zone
The liquid-gas interface of a body of subsurface water consists of a multitude of menisci with pressure jumps from the atmosphere to the liquid phase determined by the direction and magnitude of the curvatures of the interfaces, which in turn depend on the pore architecture, the various surface tensions, which determine the wetting angle (Koorevaar et al., 1983, 64-65 pp.), and the detailed liquid-gas configuration (the degree of wetting of the pore walls and the associated hydraulic connectivity).If the gas phase is well connected to the atmosphere, the reservoir providing the gas pressure effectively consists of the atmosphere of the entire planet.Thus, the atmospheric pressure is independent of the configuration of the soil solution in the unsaturated zone.The liquid pressure, on the other hand, has a component that is determined by the atmosphere: the atmospheric pressure is passed on to the hydrosphere, but not the other way around.At the Darcian scale this matters little -changes in the atmospheric pressure are simply passed on but mostly do not affect the configuration of the soil solution and the liquid-gas interfaces.At the scale of basins that are large enough to have spatial variations of the air pressure within their area, the passage of weather fronts and areas of high and low atmospheric air pressure can affect flows in larger connected systems, such as aquifers that are partially phreatic and partially confined.
In the unsaturated zone, the larger pores are filled with air (with the soil solution possibly present as a film on the pore walls), while smaller pores can be entirely filled with soil solution.In very dry soils, the soil solution is predominantly present in pendular rings around the contact points between soil particles and possibly in liquid films on the solid particle surfaces.Soil solution films in dry soil are more strongly affected by der Waalsforces, H-bonds, and, in diffuse double layers, electrostatic forces than the solution at larger distance from the solid phase behind a meniscus in filled capillaries.Therefore the degree of curvature of their liquid-gas interfaces not only depends on the liquid-gas surface tension, but also on the various forces acting upon the water molecules.The relationship between pressure and gas-liquid interface curvatures becomes ambiguous (Iwata et al., 1995, p. 7 and 13).
In saturated small pores and pendular rings, there is little if any flow (often driven by condensation and evaporation at the interface).The matric potential in these pockets of soil solution is entirely governed by the air-water interface.Its curvature determines the matric potential directly behind the interface and -through the hydrostatic equilibrium within the pore or the pendular ring -anywhere else in the pore or the pendular ring at a sufficient distance from the solid phase to be unaffected by any forces it exerts on the water molecules.For a hypothetical pore with a single exit in which the meniscus is located (Fig. 1), the intrinsic phase average (Whitaker, 1986, see below) of the matric/pressure potential when no flow occurs is: where ψ m l (Pa) is the intrinsic phase average matric potential in the pore, x * 3 (m) is the vertical position of the meniscus, A(x 3 ) (m 2 ) is the horizontal liquid-filled cross-section of the pore at x 3 , P g (Pa) is the pressure in the gas phase immediately in front of the interface (often equal to P atm ), and P lg is given by Eq. (1).For wetting liquids P lg acquires a minus sign, for non-wetting liquids a plus sign.It is interesting to note that the liquid pressure behind the interface affects and even determines the matric/pressure potential everywhere in the pore, but yet the average matric potential is determined by both the liquid pressure behind the interface and the geometry of the pore.
For a pocket of liquid enclosed by multiple liquid-gas interfaces, the pressure in the surrounding continuous gas phase will differ much less between the menisci than the liquid pressure because of the density contrast.In the pore system the menisci will be located in places where the curvatures allowed by the pore geometry and contact angles produce pressure jumps that vary with position of the menisci according to: known pressure jump can serve as the reference interface that defines x * 3 .The pocket of liquid for which Eq. ( 5) describes the pressure jumps at its interfaces is in principle unlimited in size, as long as the gas phase is continuous.The body of liquid can therefore be conceived to be the entire connected body of subsurface water in a catchment (groundwater and the soil solution), with the menisci obviously only present between the soil surface and the top of the capillary fringe.Equation (5) does not apply to menisci enclosing entrapped air bubbles because the gas pressure inside these bubbles is not determined independently from the menisci but the result of the interplay between the amount of gas, liquid pressure, pore space architecture, and interfacial tensions.
Particularly in the unsaturated zone (for which it is most relevant), the assumption of hydrostatic equilibrium will limit the direct application of Eq. ( 5) to such large scales.If the matric potential field is known, Eq. ( 5) can be generalized for non-equilibrium conditions: where x (m) is the location vector, ψ m (Pa) is the matric potential, and x * (m) denotes the location of the reference meniscus.(For clarity, only the version for wetting liquids is given.)

Averaging potentials
The values of the matric or pressure potential and the gravitational potential at a point are well defined, but if one is interested in the average potential energy of non-zero liquid volumes, these point values need to be volume-integrated and averaged.Several methods have been applied.Zehe et al. (2006) used direct volume-averaging: where ψ (Pa) can represent the gravitational, matric/pressure or hydraulic potential, and V (m 3 ) is the averaging volume within the porous domain (in the case of Zehe et al., 2006, V is the subsurface volume of an entire catchment).Whitaker (1986) only considered the volume of the phase of interest within the averaging volume and defined the phase average of the potential as: where the integration is only carried out over the volume occupied by the phase of interest within V .Since we are interested here in the liquid phase this volume is denoted by V l (m 3 ).Whitaker (1986) also defined the intrinsic phase average, in which the averaging is performed over V l only: At the pore-scale, where the distributions of the solid, liquid, and gas phases are known, integration of ψ over V l in Eqs. ( 8) and ( 9) can be achieved by integrating over V while including in the integrand an indicator function that assumes the value of one in locations within V l and zero elsewhere.The volume V l in Eq. ( 9) equals the volume integral of that indicator function over V .At the field and catchment scales, the continuum approach (Bear and Bachmat, 1991, 14-42 pp.) is invoked.In that case, to arrive at the phase average, the local value of ψ needs to be multiplied by the local volume fraction occupied by its phase, which is the volumetric water content θ for the liquid phase.The resulting expression for the phase average is: The intrinsic phase average then becomes: Equation (7a) suffers from undefined values of ψ outside V l when applied to the pore scale (for which Zehe et al., 2006, did not intend it to be used).Only when ψ is set to zero outside V l does its value there not influence the value of the volume integral.Equation (7a) then becomes equivalent to the phase average (Eq.7b), because the indicator function has the same effect.When applied to the field and catchment scale (its intended application), Eq. (7a) weighs every local value of the potential in the solid-liquid-gas continuum equally, irrespective of the local water content.As a consequence, Eq. (7a) does not satisfy Gray's (2002) additivity property that ensures that the potential energy is conserved during the volume averaging, and is therefore not recommended.The inability of Eq. (7a) to conserve potential energy can be illustrated by volume-averaging the gravitational potential energy of the water in a vertical soil column of area A extending between heights − x 3 to 0. Irrespective of the water distribution in the column, its average gravitational potential will be calculated as −ρg x 3 /2 J m −3 , which obviously is only true if the vertical soil water distribution is symmetrical around − x 3 /2.
The expressions for the phase average (Eq.7b and d) and intrinsic phase average (Eq.7c and e) both satisfy the energy conservation property, even though the numerical values of the spatial averages will differ.The phase average expresses the potential energy as Joule per volume for the entire averaging volume V , while the intrinsic phase average gives the potential energy per volume over the smaller volume V l , with its value correspondingly higher.When the respective averages are multiplied by the volumes for which they hold, they yield the same numerical value of potential energy in Joules.When ψ is uniform over V l , the intrinsic phase average ψ l will be equal to the local value ψ, which is an attractive property that makes the intrinsic phase average the most representative of the three averages of the conditions in the phase of interest (Whitaker, 1986).
The above volume integrations are all carried out for potentials expressed as potential energy per volume (Pa), and therefore weighting by volume is appropriate.For potentials expressed by mass (J kg −1 ) or weight (potential head, m), the weighting needs to be adjusted accordingly.The corresponding expressions for the intrinsic phase averages are given below (compare e.g.Gray, 2002, Eq. 5e).Expressions for the phase average are analogous.Note that the various expressions are equivalent if the phase density and the gravitational acceleration are uniform within V .
Here, µ (J kg −1 ) denotes the potential expressed per unit mass.When the potential head is used instead, the correct expression is: where H denotes the hydraulic head (m) and can be replaced by the matric head h or the gravitational head x 3 .
3 Illustrative cases and large-scale consequences

Rapid rise of groundwater levels during rainfall
In fine-textured soils, a fully saturated capillary fringe can extend over several decimeters above x p .Above the capillary fringe, the largest pores are filled with air, but most of the pore space is still filled with soil solution.A small amount of rainfall suffices to saturate the pore space over a much larger vertical extent than the equivalent water layer of the rainfall, and x p as well as the capillary fringe shoot upward (e.g., Sklash and Farvolden, 1979;Rosenberry and Winter, 1997;Seibert et al., 2003).For soils with x p only a few decimeters deeper than the height of the capillary fringe, most of the pore space below the soil surface is saturated.As a consequence, the rising capillary fringe may extend to the soil surface, causing saturation-excess overland flow (see observations by van der Velde et al., 2010, among others).
The menisci at the soil surface can no longer go up and can only flatten (Fig. 2).The pressure increase by adding a water layer then combines with the reduction of the curvature of the menisci to make x p rapidly rise to the soil surface, even though the phreatic level prior to rainfall may have been decimeters below the soil surface.
Table 1 illustrates the magnitude of this effect for hypothetical hydrophilic soils with capillaries of uniform radius.The table gives the amount of water needed to let the matric potential at the soil surface go from the air entry value to zero, thereby bringing the phreatic level to the soil surface.The air-entry value is calculated as − P lg from Eq. (1) for a zero contact angle between solid and liquid.Under hydrostatic equilibrium, with the soil surface being at the air-entry value, the phreatic level x p would be at P lg /ρg m below the soil surface.For fine-textured soils, the minute amount of water needed to bring the phreatic level to the soil surface is even smaller than for coarse soils (3rd column of Table 1), yet the matric head jump P lg /ρg that it causes (2nd column) is much larger and occurs much faster (4th column).Since fine-textured soils have the most extensive capillary fringes (2nd column), their ratio of the rise of the phreatic level and the thickness of the added water layer that causes the rise is enormous (last column).

Shallow infiltration of small amounts of rainfall/irrigation in dry soil
Small amounts of rainfall or irrigation water on well-sorted, dry soils not always percolate to the subsoil but instead only wet the top soil (Youngs, 1958).Raats (1973) explained this by the difference between the water-entry pressure at the wetting front and the air-entry pressure at the soil surface.His explanation can be cast within the framework outlined above.Incomplete wetting at the wetting front leads to nonzero contact angles and meniscus radii larger than the pore radii.The pressure jump P lg across the liquid-gas interface Fig. 2. A fine-textured soil at hydrostatic equilibrium with different phreatic levels.In the left pane the phreatic level is at such a depth that the capillary fringe exactly reaches the soil surface.The menisci at the soil surface are curved, but all pores are saturated.In the right pane, just enough water was added to flatten these menisci.
As a result the phreatic level (where the water is at atmospheric pressure) moved up to the soil surface.
at wetting front depth F (m) below the soil surface is therefore limited (Eq.1), and the water pressure behind the front only slightly subatmospheric.For the soil surface to dry, water-filled pores must empty, and these therefore will be completely wetted (meniscus radii equal to the pore radii).
The pressure jump across the interface is quite large, and the water behind the drying front will be at a considerably lower pressure than the water at the wetting front.For slow-moving fronts, the water pressure profile within the wetted layer will be nearly hydrostatic, and Eq. ( 5) can be applied to the problem.The pressure difference from the soil water to the atmosphere at the soil surface will be P lg + ρgF .According to Eq. (1) this pressure difference must be larger than the air entry pressure difference 2 γ lg /r m required for air to enter the soil, with r m the radius of the largest pore throats giving access to a continuous pore network that can be invaded by air.This will only occur if F > (2γ lg /r m − P lg )/ρg: the thickness of the wetted layer F must be larger than the difference between the matric potential heads at water entry and air entry.

Relation between large-scale average matric potential and the average curvature of the menisci
The total volume-averaged hydraulic potential of a body of water at equilibrium locked in the porous medium behind gas-liquid interfaces can be calculated from Eq. ( 4) and the expression for the pressure-equivalent of its intrinsic phase average gravitational potential ψ g l (Pa) (compare de Rooij, 2009): G. H. de Rooij: Averaged water potentials in soil water and groundwater This equation is more general if A(x 3 ) is interpreted as the water-filled portion of the horizontal cross-section at x 3 : it then applies to any body of subsurface water, stagnant or mobile.But since Eq. ( 4) is valid for a single connected body of stagnant water, the total intrinsic phase averaged hydraulic potential ψ H l (Pa) is subject to the more limiting constraints of Eq. ( 4).Its expression is: with + P lg for hydrophobic soils and − P lg for hydrophilic soils.
The inclusion of the water-filled cross-section A in the equations above limits the volume integration to the liquid phase only.This is a refinement of the volume integration of Zehe et al. (2006) who integrated over the entire subsurface volume, i.e., over the gas, solid, and liquid phases.As discussed above and in line with Gray (2002), integrating only over the liquid phase ensures that the upscaling operation conserves potential energy.
Since hydrostatic equilibrium requires the hydraulic potential to be the same everywhere, the simple expression for the volume-averaged hydraulic head expressed in the final result of Eq. ( 10) is not surprising.But it does point out that under equilibrium conditions, the gas pressure and interface curvature at any gas-liquid interface, and its vertical position, suffice to describe the total volume-averaged hydraulic potential in a body of subsurface water at equilibrium, but not its volume-averaged matric/pressure potential ψ m l , not even for a water pocket behind a single pore.This is contrary to Zehe et al.'s (2006) hypothesis that the volumeaveraged capillary pressure should reflect the average pressure drop across the liquid-gas interfaces in the pores at equilibrium.For multiple-pore systems, the situation becomes even more troublesome, since the average pressure drop over the menisci will not only be governed by the matric potential field (see Eqs. 5 and 6), but also by the distribution of menisci.Below the top of the capillary fringe, for instance, the number of menisci drops to zero while a significant portion of the soil water is likely to reside in the capillary fringe; the other extreme is the dry top layer of coarse soils, where a meniscus is present at each pendular ring of water that retracted around a contact point of solid grains.A large population of menisci thus represents only a small quantity of water.The relation between average matric potential and average meniscus curvature is analyzed in detail for a bundle of capillary tubes reflecting a soil water characteristic in the Appendix.Zehe et al.'s (2006) hypothesis gains credibility under unit gradient flow conditions.Then, the flow is driven by gravity only, and the matric potential is uniform.(Such conditions are unlikely in saturated conditions, and we therefore disregard the pressure potential here.)In this case, the difference between the matric potential terms in Eq. ( 6) is zero.The matric potential determined behind any individual interface is representative for the entire body of water, and the average meniscus curvature is obviously equal to each of the individual curvatures.This warrants the conclusion that Zehe et al.'s (2006) hypothesis is correct for unit gradient flow conditions, but not for the equilibrium conditions as stipulated by the authors.Zehe et al. (2006) formulated their assertion for a population of interfaces and large bodies of water, and for the case of unit gradient flow that is appropriate: at the pore scale, and even at the Darcy-scale, unit gradient flow cannot exist in heterogeneous media because it precludes lateral flows by assuming horizontal matric potential gradients to be zero.Thus, the redistribution of water needed to keep the local vertical flux densities equal to the local hydraulic conductivities (a necessity if the vertical matric potential gradient is to vanish) is not allowed.Nevertheless, unit gradient conditions have been observed in the field at depths sufficient to dampen out fluctuations in the boundary conditions at the soil surface.At scales slightly larger than the Darcy scale (1 m and beyond) unit gradient flow can safely be assumed to occur if conditions are favorable (e.g., Davidson et al., 1969, Reichardt et al., 1998): the lateral matric gradients inevitably caused by soil heterogeneity often turn out to be relatively small owing to the dissipative nature of matric potential gradients.
Neither hydrostatic equilibrium nor unit gradient flow is likely to occur over an entire catchment.But even for smaller systems (e.g., a field of a few hectares discharging into a ditch or stream) the required conditions are problematic: hydrostatic equilibrium implies there is no precipitation, no evapotranspiration (and hence to plant growth), no groundwater recharge or lateral groundwater flow, etc.Under circumstances with a small downward water flux in the unsaturated zone and a conductive, sufficiently thick phreatic aquifer, the lateral groundwater flow necessary to discharge the groundwater recharge form the unsaturated zone may induce only minor lateral gradients and hardly affect the vertical gradients within the groundwater and the capillary fringe.Under such circumstances, near-hydrostatic equilibrium conditions may prevail below the top of the capillary fringe, but there is no obvious relationship between the average potential in this area and the average curvature of the liquid-gas interfaces in the depth range above it (the only region where such interfaces occur).
Hydrol.Earth Syst.Sci., 15, 1601Sci., 15, -1614Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1601/2011/ Unit gradient flow conditions at least allow water to move, but only downward.Unlike hydrostatic equilibrium, unit gradient conditions can be found for prolonged periods of time in the field: arid areas with some rainfall and deep groundwater tend to develop unit gradient conditions at some depth below the soil surface.Still, these conditions apply only to the unsaturated zone below the depth at which the signals of precipitation and evapotranspiration are fully damped, and not to the shallow soil or the groundwater.The vertical flux is typically small.Horizontal groundwater flow is necessary to carry off the recharge from the unsaturated zone.Below the top of the capillary fringe, lateral gradients will therefore be present.For conductive and sufficiently thick aquifers, this lateral gradient is very small.In that case, within the groundwater and the capillary fringe, as well as in a poorly defined region above it, hydrostatic equilibrium might be a better approximation of the potential field below the depth range where unit gradient conditions prevail.If the aquifer is less conductive and/or thin, appreciable lateral gradients may develop that complicate the potential field.
In either case, the presumed conditions can only exist in a part of the subsurface flow domain.Only in the case of unit gradient flow does the desired condition appear in the region which also contains the liquid-gas interfaces.Still, the hypothesis that the average pressure drop over liquid-gas interfaces reflects the average potential of all the subsurface water in a catchment does not seem tenable.Only in the case of unit gradient flow does the (uniform) curvature of these interfaces reflect the (uniform) matrix potential, but only within the depth range in which the unit gradient conditions apply.

Relation between large-scale average hydraulic potential and saturated flow across groundwater system boundaries
In the unsaturated zone, the concept of average potential energy appears to be of limited value: the matric potential can change rapidly in response to minor changes in the configuration of the liquid-gas interfaces without major or lasting effects on the flow.But in groundwater flows, average potential energies may still provide a useful tool to describe flows at a more integrative scale if two criteria are met: the flow must be confined within a finite volume to be able to compute the intrinsic phase average of the groundwater body, and the groundwater body must have a well-defined boundary across which a flux can be determined.Many classical groundwater problems for which analytical solutions are available meet these criteria.For problems that are described by linear differential equations, the superposition principle applies.Such systems can therefore be expected to exhibit a linear relationship between the intrinsic phase average of the hydraulic potential and any flux across a system boundary.
To establish a proof of principle, three classical flow problems with parallel, radially diverging, and radially converging flow are recast in terms of the intrinsic average hydraulic  potential and the flow across a system boundary -in line with the desirability to quantify fluxes at larger-than-Darcy scales the focus thus is on the flow across system boundaries, rather than the details of the flow within the system.In the analysis, the hydraulic head H instead of the hydraulic potential is used to adhere to established convention.The first system to be considered is a confined, non-sloping aquifer of finite thickness D (m) and width L (m), and infinite length.The aquifer has a uniform saturated hydraulic conductivity K (m d −1 ) and is recharged from above at constant and uniform recharge rate R (m d −1 ).The horizontal coordinate x 1 (m) runs from 0 to L. At x 1 = 0, the aquifer is bounded by a no-flow boundary; at x 1 = L, the aquifer is in contact with open water with a constant level H 1 above the bottom of the aquifer (a no-flow boundary), which is also the reference height where the vertical coordinate x 3 is defined to be zero.The Dupuit assumptions apply (zero vertical gradient in H , recharge instantly divided over the full vertical cross-section).Note that this flow pattern also approximates that in a field drained by ditches that are spaced 2L apart and reach down to the impermeable layer (Fig. 3a).The aquifer then is phreatic and the vertical saturated extent D(x 1 ) is determined by the local groundwater level: D(x 1 ) = H (x 1 ).In that case, the difference between the maximum groundwater level located halfway between the ditches (where symmetry considerations impose a vertical no-flow boundary) and H 1 must be negligible compared to This system leads to a parabolic H (x 1 ) relationship: For uniform porosity θ s , and D assumed constant, the intrinsic phase average of H is: During steady state flow, the flux Q(0) per unit width (m 2 d −1 ) into the open water is equal to the total recharge RL.Solving Eq. ( 12) for RL then gives: This expression is the system-scale equivalent of Darcy's Law: it expresses a flux across the system boundary as a function of the difference in average potential energy of the water in the systems on either side of the boundary and the system properties.
In case of a phreatic aquifer, flow towards the ditches will continue after recharge ceased, until the groundwater level equals H 1 everywhere.We can hypothetically extend the steady-state proportionality between the flux across the boundary and the difference between averaged hydraulic heads to the non-steady receding phase that commences after recharge stops at a given time t 0 (d).As the groundwater table lowers, some water will remain behind in the previously saturated portion of the soil, and the unsaturated zone above the receding phreatic level will respond to the changing conditions.The net effect is that the delivery of water to the ditch will be less than the water stored in the soil volume that became unsaturated.Thus, the drainable storage at any time is smaller than the water in the saturated zone above H 1 , and is simplified as θ d L( H l − H 1 ).Here, θ d can be viewed as the volume fraction of the initially saturated soil in the field that will release its water when the groundwater level drops to H 1 .For deep groundwater tables and very dry conditions the theoretical maximum of θ d is θ s − θ r (the latter being the residual water content), but normally it will be smaller.Thus, maintaining the proportionality of Eq. ( 13) for zero recharge translates into proportionality between the flux and drainable storage, which leads to an exponentially decreasing flux: in which the right hand side without the proportionality constant 3KD/L describes the temporal evolution of H l − H 1 .
The value of H l (t 0 ) is given by Eq. ( 12).This linear reservoir-type behaviour of a drained field was already reported by de Zeeuw and Hellinga (1958).A more rigourous treatment of non-stationary discharge by Dumm and Glover can be approximated well by an exponential decay, particularly for prolonged times (Dumm, 1954).So, for this case of transient flow, the field-scale analysis relying on averaged hydraulic heads generates results consistent with classical Darcian approaches.
The second case (Fig. 3b) involves radially divergent flow in a circular aquifer surrounded by open water at H 1 : an aquifer in a circular island in a freshwater lake.The lake bottom and the impermeable layer underneath the aquifer are located at x 3 = 0 m.Invoking the Dupuit assumptions and considering only steady flows leaves only the radial coordinate r (m), which runs from 0 to L. Again, we assume the difference between H (0) at the centre of the island and H (L) = H 1 to be negligible (D(r) ≈ H 1 ).The expression of the phreatic level becomes: The intrinsic phase average of the hydraulic head is: The flux Q(0) discharging into the lake is π L 2 R. It can be expressed as: Not surprisingly, the radially diverging flow has a considerably larger proportionality constant than the parallel flow system for realistic values of L. Interestingly, the radius of the island does not appear in the constant (compare with Eq. 13).
Both cases above have a parabolic phreatic level.The third case is that of a radially converging flow in a circular aquifer of thickness D enclosed on the outer rim by a no flow boundary at radius L. The constant recharge R flows through the aquifer towards a fully penetrating well in the center, with radius r w .The hydraulic head inside the well is maintained at H 1 , which is such that the aquifer remains saturated over its full thickness, even for small r.In this system, the expression for the hydraulic head becomes: which, for uniform porosity, leads to the intrinsic phase average: which, for r w L, simplifies to (with the right-hand-side included to facilitate comparison with the expression for diverging flow in Eq. 16 In comparison to the proportionality constant for diverging flow in Eq. ( 17), this constant is considerably smaller, reflecting the need for larger gradients to achieve the same flux (not flux density) during converging flow.This also demonstrates that the system-scale proportionality constant comprises properties of the system as well as the flow.In general, a proportionality constant for a groundwater system will only be valid for flow patterns similar to that for which the constant was determined.
This flow-specific nature of the proportionality constants may also apply to aquifer transmissivities determined from groundwater and contour lines of the hydraulic head.The transmissivities could change when the direction of the flow changes, or when the flow pattern changes, e.g. from parallel flow to radial flow after installation of a drinking water pumping well.

Summary and conclusions
To study the energy status and flow of subsurface water at scales larger than the Darcian scale, the relationship between changes in location and curvature of gas-liquid interfaces in soils and the matric/pressure potential of the water behind those interfaces was examined.The relevance of porescale processes at larger scales was illustrated by the role of menisci in highly dynamic matric/pressure potential responses of fine-textured soils to rainfall: one well-known part of this response is the result of limited available storage in the pore space, another part is caused by extremely rapid water pressure changes when curved interfaces at the soil surface are flattened.A detailed analysis of the relationship between the curvature of a single meniscus or a population of menisci and the volume-averaged matric potential of a body of subsurface water did not offer theoretical support for the direct relationship between average matric potential and average meniscus curvature at equilibrium that was hypothesized in the literature.However, for strictly gravity-driven flow, a direct proportionality between average meniscus curvature and average matric potential is plausible.Still, this is of limited value for practical situations.
Various methods for volume averaging water potentials were reviewed and the intrinsic phase average deemed the most appropriate.For saturated conditions, a proof of principle was established that demonstrates that for flows that are governed by linear differential equations, the flux across a boundary of a hydrological system is proportional to the difference between the intrinsic phase averaged hydraulic potentials on either side of the boundary.This can greatly simplify the description of the interaction of groundwater bodies with their surrounding elements of the hydrological cycle: average groundwater levels may suffice to generate discharge estimates.
G. H. de Rooij: Averaged water potentials in soil water and groundwater tubes with radii within the range defined by the two values.Letting this range go to zero gives the derivative: If N (m −2 ) denotes the number of water-filled tubes per square meter and n (m −3 ) its derivative dN/dr, Eq. (A4) provides the number of capillary tubes per square meter with radius r by noting that dθ/dr is the fraction of the total tube area occupied by the combined inner area of these tubes, i.e., n • πr 2 : In a bundle of capillary tubes with its lower end submerged in a water reservoir (mimicking the groundwater), all tubes will contain water, with the level rising to a/r above the water level (Eq.A1).Equation (A5) then gives the number of menisci per square meter with curvature 2/r (m −1 ).The average curvature of the menisci in the entire bundle of tubes can then be found by deriving the distribution function of 2/r from Eq. (A5).By changing variables to make y = 2/r (thereby letting y (m −1 ) denote the meniscus curvature), Eq. (A5) can be written as: Note that the minimum curvature is −ρgh ae /γ and belongs to the largest capillary, emptying at h ae .The curvatures in tubes that are water filled over their entire length (where h = −x p if we set x 3 equal to zero at the soil surface) are ignored, since they represent water pockets in fine pores behind larger pores with menisci represented by the larger capillaries.The arithmetic mean of the curvature y a (m −1 ) then is: which reduces to: to give The radius of a meniscus in a fully wetted cylindrical capillary corresponding to this curvature would be 2 y a .The matric head of water just behind this meniscus is: We can also calculate the average matric head of all the water in the population of tubes.For a tube of radius r, the water rises to 2γ /ρgr m above x p , and the tube therefore contains 2π rγ /ρg m 3 of water.The total volume of water per square meter that is stored in tubes in which the water rises to 2γ /ρgr m above x p is n(r) • 2π rγ /ρg, and that volume of water has an average matric head of −γ /ρgr m.The overall average matric head equals the weighted average of these average matric heads for all r, with the weighting factor being the amount of water residing in tubes of radius r.The tubes that are so narrow that they are filled with water all the way to the soil surface all have an average matric head of −x p /2.The largest of these tubes has radius −2γ /ρgx p .
The average matric head of the population of tubes smaller than this is: To arrive at the overall average matric head h a (m), the average matric head for the partially filled tubes h 1 a and that for the full tubes (−x p /2) have to be averaged, weighted by their respective portions of the total water they represent: Clearly, Eqs.(A10) and (A16) are unequal, and for the case of a soil simplified to a bundle of capillary tubes, the matric potential derived from the mean curvature of the menisci does not represent the volume-averaged matric potential of the water above the groundwater table if the water is at hydrostatic equilibrium with a fixed groundwater table.
One may wonder if the analysis above would not give completely different results if the vertical capillaries were connected and could exchange water.This can be tentatively assessed by Fig. A1.Only when a horizontal throat connects a filled capillary to an empty one does it create an extra interface.The right-hand side figure represents a throat that is wider than one of the capillaries it connects, which would create an interface that is less curved than it could have been in that capillary.This situation does not seem probable in much more complex natural pore networks.The extra interface in the connecting throat in the left hand side of the figure would not occur in a bundle of isolated capillaries.Still, the population of throats in the natural pore space would hold water at low matric potentials and thus be visible in the water retention curve.Since the radii of the bundles are chosen such that they collectively represent the entire pore space, the interfaces present in these narrow throats are represented by tubes with small radii in the bundle model.
The capillary bundle representation would probably not change very dramatically if the vertical capillaries were connected by throats, since only those throats that happen to connect two capillaries at an elevation where precisely one of them is water-filled would add an interface.Nevertheless, it remains doubtful that the bundle model accurately represents the total number of interfaces in a soil, even if it does accurately reproduce the soil water characteristic.

Fig. 1 .
Fig. 1.A hypothetical liquid-filled pore with a single exit in which the liquid-gas interface is located at vertical position x * 3 .The horizontal liquid-filled cross-section A at arbitrary elevation x 3 is also indicated.

Fig. 3 .
Fig. 3. Two of the three groundwater flow cases governed by linear differential equations, as given in the figure.Figure (A) reflects parallel flow in cartesian coordinates, Figure (B) depicts radial flow in cylindrical coordinates.The dashed vertical line at the left of each figure is the axis of symmetry.At L, the groundwater body is in contact with open water.The open water and the aquifers are bounded below by an impermeable base.Further details are in the text.
Fig. 3. Two of the three groundwater flow cases governed by linear differential equations, as given in the figure.Figure (A) reflects parallel flow in cartesian coordinates, Figure (B) depicts radial flow in cylindrical coordinates.The dashed vertical line at the left of each figure is the axis of symmetry.At L, the groundwater body is in contact with open water.The open water and the aquifers are bounded below by an impermeable base.Further details are in the text.

−
average matric head relates to a volume of water present in those tubes in a 1 m 2 area that are not entirely filled over their full length −x p .This water volume is given Hydrol.Earth Syst.Sci., 15, 1601-1614, 2011 www.hydrol-earth-syst-sci.net/15/1601/2011/ by the denominator of Eq. (A11).The remaining portion of the water resides in the capillaries that are filled with water up to the soil surface.The volume of water per square meter of soil stored in these capillaries with radii between 0 and −2γ /ρgx p is: −2γ /ρgx p 0 −x p π r 2 n(r)dr = ρgλh ae (θ s − θ r )x p 2γ

Fig. A1 .
Fig. A1.Hypothetical distributions of liquid, gas, and liquid-gas interfaces for a simple system of vertical capillaries connected by horizontal throats.The lower end of the vertical capillaries is in contact with open water at a constant level (not drawn).

Table 1 .
Amount of water needed to bring the phreatic level up to the soil surface when initially the capillary fringe extends to the soil surface, and the time needed to deliver this water during rainfall of moderate intensity.The soils are assumed to have cylindrical pores of uniform radius.The soil water is 15 • C (with corresponding density and air-water surface tension), and the gravitational acceleration is 9.812 m s −2 .For the third column, a porosity of 0.60 is assumed for pores ≤0.1 mm, and of 0.40 for larger pores.Note that the last row gives the limiting radius where the capillary rise equals the pore radius.