References

In order to establish a thermodynamic justification of the theoretical relationship between the hydraulic, pressure/matric, and gravitational head in subsurface hydrology, the thermodynamic literature pertaining to subsurface flow processes is reviewed. The incompressibility of liquids negates a thermodynamic definition of pressure, which gives rise to several inconsistencies in pore scale theories. At larger scales, the gravitational potential and fluid pressure are treated as additive potentials. This superposition principle is replicated in the well-established relationship between the various heads according to subsurface hydrological theory. The necessary requirement that the superposition be maintained across scales is combined with conservation of energy during volume integration to establish consistent upscaling equations for the various heads. The power of these upscaling equations is demonstrated by the derivation of an upscaled water content-matric head relationship and the resolution of an apparent paradox reported in the literature that is shown to have arisen from a violation of the superposition principle.


Introduction
In saturated and unsaturated water flow in porous media, the continuum approach is generally invoked (e.g., Bear and Bachmat, 1991, p. 14-31), and the resulting laws of mass conservation and movement, as well as the variables appearing in there, are defined for a mathematical point centered within a Representative Elementary Volume 20 (REV). The resulting equation for water movement is Darcy's Law. Its general acceptance is illustrated by the fact that the REV-scale has frequently been referred to as the Darcy scale (e.g., Bhattacharya and Gupta, 1983;Jacquin and Adler, 1987;Kavvas, 2001, among many others).
In order to develop theories (and their associated measurable variables) at scales Introduction selves, there is a drive to upscale these Darcy-scale relationships and variables to larger spatial scales. In doing so, pertinent questions arise as to the physical soundness of upscaling variables such as fluid pressure, matric and hydraulic potential, and water flux. It must be stressed that this problem differs from upscaling from the pore scale to the Darcy scale. Nonetheless, the rigour with which the fundamental ques-5 tions surrounding the nature of the upscaling operation and the physical relevance of the resulting upscaled properties and variables have been addressed for scale transfers starting at the pore scale may provide valuable insights for the next-larger scale transfer addressed here. At the pore scale, the solid phase, the soil solution, any other fluids present, and 10 the soil gas share interfaces where two phases meet, contact lines where interfaces meet, and contact points were contact lines meet. The various surface tensions and pressures together with gravity govern the distribution of the various fluids and soil gas over the pore space. A sizeable body of literature exists on the thermodynamic description of such systems at the pore scale, of which Gray and Miller (2007) provided 15 an insightful overview. Gray (2002) presented a systematic method to upscale the pore-scale thermodynamic expressions for internal, kinetic, and potential energy of the various phases and their interfaces and contact lines to larger scales (Darcy scale and beyond) by integrating over a soil volume. The most recent culmination of this on-going work appeared in 20 Gray andMiller (2005, 2006) and Miller and Gray (2005), where a systematic thermodynamically constrained averaging theory was laid out. The resulting equations bear some similarity with the early work of Groenevelt and Bolt (1969) on the non-equilibrium thermodynamics of subsurface water flow, although the former preferred to work in a Langrangian framework while the latter favored an Eulerian approach. 25 Nordbotten et al. (2007) recently proposed a volume-average of the Navier-Stokes equations for saturated flow over a REV, without the term accounting for the kinetic energy. While integration over a soil volume appears acceptable for thermodynamic quantities, a volume-averaged flux has no physical meaning, and Nordbotten et al. (2007) Gray and Miller (2005) argued against the inconsistencies in and the hysteretic nature of many Darcy-and macroscale models and propose instead a rigorous development from pore scale physics (with the assumptions required to achieve the desired upscaling identified explicitly instead of being tacitly incorporated in the model). The general upscaling methodology of Miller and Gray (2005) and Gray and Miller (2006) 10 does not have the limitations of Nordbotten et al. (2007). But its complexity is such that conventional modeling practices will likely persevere for some time, particularly for the unsaturated zone. Fortunately, a full thermodynamic treatment may not always be necessary; at the Darcy-and the macroscale, the various forces and pressures acting in the pores are efficiently unified in the core-scale matric and gravitational potential 15 while the kinetic energy can generally be ignored. The configuration of the fluids in the pore space as described by their respective indicator functions coalesce into the core-scale volumetric fractions of the bulk medium occupied by the various phases (readily calculated from the volume integrals of their indicator functions; Nordbotten et al., 2007), at the cost of these fractions becoming hysteretic. The practical applicability 20 of existing models, their conceptual and thermodynamic shortcomings notwithstanding, is supported by the single-fluid case elaborated by Gray and Miller (2006), which yielded a result very similar to the classical groundwater flow equation.
Whatever the method of upscaling may be, gravity gains importance as the vertical dimension of the system under study expands. In natural systems, gravity manifests 25 itself in such phenomena as the finite height of capillary rise in the unsaturated zone, and in fluid pressure in the saturated zone. The connection between this 'gravitational' pressure and the thermodynamic definition of pressure as the negative of the partial derivative of internal energy with respect to volume has been treated in various ways HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential in the literature and will be elaborated on in Sect. 2.1. The role of gravity and pressure is apparent in the work of Gray and Miller (2004), who studied the upscaling of Darcy's Law in an idealized system involving a configuration with horizontal flow in a volume containing two porous media with different porosities. Gray and Miller (2004) averaged over a cross-section perpendicular to the main 5 direction of flow, thereby bypassing the difficulties associated with averaging fluxes over volumes (Nordbotten et al., 2007). They reported inconsistencies in Darcy's Law if the layering of the media was tilted. Their analysis also discussed the averaging of fluid pressure. In the unsaturated zone, the matric potential is the equivalent to fluid pressure in the saturated zone. Therefore, the work of Gray and Miller (2004) may   10 also have ramifications for averaging the matric potential of a soil volume, and thus for the determination of soil water retention curves at various scales. The effect of gravity on the nature of the soil water retention curve was experimentally addressed in another context by Liu and Dane (1995), and their findings will be made to bear on the upscaling formalism. 15 This paper first discusses the treatment of pressure in various thermodynamic models for the distribution and movement of fluids in porous media, and tries to reconcile that with established soil physical and groundwater hydrological theory. It then connects some of the results of the thermodynamic analysis to more traditional soil physics by developing consistent volume averaging equations for various soil physical 20 variables. Problems associated with upscaling Darcy's Law itself will be clarified. The paper then examines the inconsistencies reported by Gray and Miller (2004) between formulations of Darcy's Law and its upscaled equivalent, which was illustrated for a special case. In the thermodynamic analysis of Gray (2002), the total energy of a phase that occupies a portion of the pore space is the sum of internal, kinetic, and potential energy. In 5 classical thermodynamics, the differential of the internal energy of a single-component system is (Callen, 1985, p. 35): The definition of pressure is intimately connected to the definition of work W [ML 2 T −2 ] done by the system by expanding its volume against an external pressure p: with the minus-sign arising from the convention that work is positive if it is performed by the system, thereby reducing its internal energy (Callen, 1985, p. 19-35). Gray's (2002) definition of the internal energy of a unit volume of fluid or gas at the pore scale is (excluding the energy of its interfaces with other phases): , ρ the density [ML −3 ] of the phase (assumed by Gray to be a single component of perfect purity, but an extension to solutions and mixtures is straightforward), and µ M the chemical potential of the phase per mass [L 2 T −2 ]. This is the Euler relation (Callen, 1985, p. 59-60), normalized by V . While the change in the units of 5 mass and the resulting dimensions of the chemical potential compared to Eq. (1) is trivial, the normalization by volume is key. This implies incompressibility of the fluid, just as the more common normalization by mass implies constant mass (vz. Eqs. 2.23 and 2.24 of Callen, 1985, p. 39). The incompressibility assumption is quite accurate for the pore solution but renders the thermodynamic definition of pressure meaningless: a volume that cannot expand or contract cannot perform work according to Eq.
(2), leaving p undefined. Since, for such a fluid, dV ≡0, the pressure term drops out of Eq. (1), and a definition of p is not required. Gray (2002) consequently left out the term involving pressure in his Eq. (13a) (essentially equivalent to Eq. 1 above). This is consistent with the equation that Groenevelt and Bolt (1969) obtained by applying the normaliza-15 tion by volume directly to Eq. (1). Obviously, the presence of molecules at non-zero absolute temperature ensures the existence of a pressure in the fluid volume comprising the system under study, but its value cannot be established from the principles of classical thermodynamics. This leaves the pressures in Gray's (2002) thermodynamic derivations somewhat ambiguous, particularly since a flow velocity field is imposed. 20 The only condition under which the pressure is not affected by flow is unit gradient or gravity flow, when the positive downward flux density is equal to the hydraulic conductivity (e.g., Wagenet, 1984;Jury et al., 1985, p. 100) and loss in gravitational potential exactly matches the energy dissipation by viscous flow.
In contrast with Gray's (2002) assumption of fluid incompressibility, Gray and Miller 25 (2007) presented the classical thermodynamic definition of pressure as in Eq. (2) above. Since they normalized the Euler relation by mass instead of volume, they did not implicitly assume incompressibility. However, they explicitly addressed the thermodynamics of interfaces between a wetting and a non-wetting phase, implying that the HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential Interactive Discussion two did not mix well, which could only be the case if at least one of the two was a liquid. The probable incompressibility of at least one phase complicates the interpretation of the derivations with respect to the density that occur in the equations for fluid pressures, interfacial tensions, and the differentials of internal energy. Both Gray (2002) and Gray and Miller (2007) strived to upscale pore scale variables 5 and relationships by rigorous volume averaging. Consequently, they could not describe the effect of the macroscale pressure field on the flow since these pressures do not emanate from local conditions in the pores. Instead, they are imposed by large-scale factors like climate, land use, aquifer structure, location within the flow domain, etc., and therefore cannot and should not be inferred from pore-scale information alone. 10 While Gray and Miller (2007) focused on equilibrium conditions and did not consider flow, Gray (2002) accommodated the large-scale forcing by imposing a macroscopic velocity field, which he explicitly incorporated in the equations. Groenevelt and Bolt (1969) and Gray and Miller (2006) aimed to develop macroscale equations and therefore were forced to include the macroscale pressure field and cal- 15 culate, rather then impose, the resulting macroscale flow velocity field. Both did so by treating the pressure as a potential and adding it to the gravitational potential. Groenevelt and Bolt (1996) simply declared the pressure a potential that was additive to the gravitational potential. Gray and Miller (2006) allowed the volume fraction of a phase to change, which in turn allowed its differential to be non-zero. They equated the macro-20 scopically observable pressure to the thermodynamic pressure according to Eq. (2). By doing so they could replace µ M by the pressure potential by applying the Gibbs-Duhem relation (Callen, 1985, p. 60-62) for isothermal conditions. Thus, from Eq. (100) in Gray and Miller (2006): where ψ is the total potential per volume [ML −1 T −2 ], g is the gravitational acceleration [LT −2 ], x 3 is the vertical coordinate [L], and x 3,ref [L] is an arbitrary reference height.
Note that the inclusion of x 3,ref makes the potential relative. Since the gradient of ψ is HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential Interactive Discussion all that is needed, this is inconsequential. The superposition of the potentials created by the pressure field and by the gravitational field suggests that methods developed to average one potential should be applicable to the other as well. Furthermore, the superposition principle implies that, in any upscaling operation through systematic averaging, the averaging manipulations be 5 identical for both potentials to maintain the superposition property across scales. Note that, as long as x 3,ref is kept constant, it does not impair the superposition property.
In subsurface hydrology, the hydraulic and matric head are related as (e.g., Brutsaert, 2005, p. 268;Jury et al., 1985, p. 77): It is immediately clear that this textbook equation can be obtained by dividing the thermodynamically established Eq. (5) by ρg and replacing p by p minus the atmospheric pressure to allow h to be zero at the phreatic level. In this way the pressure is made relative in the same fashion as the elevation in Eq. (5). Without loss of generality we can set the arbitrary reference height x 3,ref [L] to zero: For h≥0 in a multiphase system with air at instantaneous equilibrium with the atmosphere and a fluid, the pressure in the fluid exceeds the atmospheric pressure, and h is routinely termed the pressure head. For wetting fluids this implies that the medium is saturated with this fluid. For matric potentials <0, the fluid pressure is smaller than 20 atmospheric pressure, and h represents the matric and interfacial forces that retain the fluid in the pores. It is therefore termed matric head. Note that the magnitude of the matric forces can cause values of h that are so low that the pressure equivalent would be negative. Since negative pressures are physically unacceptable, Hassanizadeh and Gray (1990) proposed a wettability potential in addition to the fluid pressure to represent the effect of matric and interfacial forces on a phase. This wettability potential quantifies the energy change of the phase associated with a change in saturation. The 5,2008 Thermodymanics and upscaling of subsurface water potential sum of the fluid pressure potential (absolute or relative to the atmospheric pressure potential) and the wettability potential gives h. The wettability potential satisfies the superposition property of the pressure and the gravitational potential. The excessive pressure of the air for h<0 allows it to invade the pores. The extent to which it will do so is described by the fluid retention curve, which indicates how much 5 fluid resides in the pore space at a given matric head. Thus, this generally hysteretic curve, which is specific for a combination of medium and fluid, facilitates a consistent characterisation of the system in terms of fluid content and energy status of the fluid (Jury et al., 1985, p. 45-68;Hillel, 1998, p. 129-168).

HESSD
2.2 Upscaling by spatial averaging 10 The discussion here pertains to upscaling within a soil volume that is large enough to be heterogeneous in a non-stationary sense (Cressie, 1993, p. 52-53), i.e., trends in soil hydraulic properties may be present. At the scale of the analysis the continuum approach is assumed to be valid. The pore architecture is such that the flow is laminar anywhere, e.g., at the REV-scale, the validity of Darcy's Law is not disputed: water content θ(x 1 ,x 2 ,x 3 ,t), h, H, q (all functions of x 1 , x 2 , x 3 , and t), and K (θ). For saturated flow, K (θ) simplifies to K (x 1 ,x 2 ,x 3 ). By letting n vary in space but not in time we exclude from the analysis porous media with varying pore space (i.e., swelling and shrinking soils), and implicitly limit the discussion to time scales much smaller than those of interest for most geological and soil morphological processes. 5,2008 Thermodymanics and upscaling of subsurface water potential Large-scale averages of porosity and volumetric water content are easily found and their physical meaning is immediately clear:

HESSD
where V [L 3 ] denotes an arbitrary volume occupied by a porous medium, and the sub-5 script V denotes the volume-averaged value of the subscripted variable. Note that Eqs. (8) and (9) are themselves upscaled values of volume averages of indicator functions that take on the value of 1 whenever they are located in a pore or in the water phase, respectively (see also Nordbotten et al., 2007). The main difference is observational: measuring n and θ is feasible in realistic porous media, while the indicator 10 functions can only be determined in small, simplified systems. Averages of n and θ over an arbitrary area or line along the principle directions within the porous medium can be found by a corresponding reduction of the dimensions over which the integrals in Eqs. (8) and (9) are performed. As shown above, the hydraulic head represents the total energy of the water at a 15 given location. Its spatial average should therefore reflect the total energy of the water present in the porous medium volume for which the average is determined. Therefore, the local values should be weighted by the local water amount. With H being the energy per unit weight (Hillel, 1998, p. 153), weighting by the local weight (θρg) would be consistent. But since ρg is assumed constant, weighting by θ is acceptable: HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential with the subscript V denoting a volume average as above. The weighting by θ ensures that Eq. (10) satisfies the additivity condition advocated by Gray (2002) in that it conserves the total energy irrespective of the size of V . Because of the identical averaging manipulations required by the superposition property of the components h and x 3 of H and, we also have: and which identifies the horizontal plane around the center of gravity of the water in V . Note that the numerators of Eqs. (10)(11)(12) have dimensions L 4 and are measures of the total, 10 matric, and gravitational energy, respectively, stored in a body of subsurface water at a given time. If desired, the relative nature of H, h, and x 3 can be removed, and the resulting absolute values can be multiplied by ρg to obtain the respective energies in Joules.
As with n and θ, averages of the various potentials over an arbitrary area or line 15 along the principle directions within the porous medium can be found by modifying the dimensions over which the integrals in Eqs. (10)(11)(12) are performed. In many cases, a flow has a well-defined macroscopic direction, for instance because the flow domain is enclosed by impermeable barriers on all but two opposite sides. Integrations over cross-sections perpendicular to the macroscopic flow direction can then provide aver-20 age potentials that are local with respect to the coordinate in the flow direction. This is particularly useful for calculating upscaled head gradients.
HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential Any of the tensorial components of K can be volume-averaged, but the resulting average does not bear a relationship to any flux at the scale of V . For K this direct approach seems to be of limited value. Averaging the components of q is more fruitful. Here, volume averages are less useful (as Nordbotten et al., 2007, also noted) than averages over some area (e.g., the boundary of V , or a cross-section through V ). For where the off-diagonal elements of K were assumed to be equal to zero. With averaged flux densities thus derived, the question arises whether a hydraulic conductivity can be found that relates the flux densities to some average gradient in H, effectively enforcing a scale-invariant formulation of Darcy's Law: 15 where K j,A [LT −1 ] is the areally averaged hydraulic conductivity in the direction of x j , and f (x 1 ,x 2 ,x 3 , t) is an as yet undefined weighting function. An intuitively appealing weighting function is θ (x 1 ,x 2 ,x 3 ,t), in analogy with Eqs. (10)(11)(12). However, Darcy's Law (Eq. 7) does not have a term with θ. It describes flux densities, for which K (related to water movement) is a more relevant quantity than θ, which is related to the HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential If we tentatively set f to 1 everywhere (uniform weighting), Eqs. (13) and (14) combine to: The definition of K j,A follows immediately: Equation (13) arrives at an areally averaged flux density by arithmetically averaging the local flux densities over A. According to Eq. (16), the corresponding hydraulic conductivity of A is found by weighting the local values of K j by the local hydraulic head gradient ∂H/∂x j . The derivation involves an arithmetic areal average of ∂H/∂x j in Eq. (15).
In the above we saw that the total energy in a volume of subsurface water (or any cross-section thereof) can be found by a θ-weighted average of the H, but that the areal average of the flux density involves arithmetic averaging of the gradient of H. The local energy gradient is closely related to the local flux density, as described by Darcy's Law. The connection between H and q at the non-point scale is not similarly defined. The 15 first conclusion is that upscaling Darcy's Law through consistent averaging of H and K while keeping the differential equation itself identical to its point-scale form (Eq. 5) is impossible. Furthermore, for both saturated and unsaturated flow K becomes a function of H, creating the non-linearity that at the Darcy scale only manifests itself in the unsaturated flow equation (Richards' equation; Jury et al., 1985, p. 105-109;Hillel, 20 1998, p. 212-214). A third conclusion is that the total energy of a body of water is of little use in describing its tendency to generate subsurface flow.
From a practical point of view, the limited use of H V (or the areal average H A ) for large-scale flow calculations is not problematic. In many aquifers, H is fairly uniform 1150 Introduction mation about the distribution of n and θ, and about the boundaries of the subsurface water body of interest, all of which is difficult to obtain.

Upscaling the water content-matric head relationship
For the upscaled version of the soil water characteristic h(θ), the relation h V (θ V ) is 10 an obvious candidate, but the experimental conditions of h(θ) need to be taken into account. At the Darcy scale, an undisturbed sample is typically exposed to a fixed H long enough to establish hydrostatic equilibrium (e.g., Dane and Hopmans, 2002a, b, c, d;Romano et al., 2002). From the elevation of the sample with respect to the chosen x 3,ref the value of h at the center of the sample is derived. The total water content of 15 the sample is determined (usually by weighing) to calculate a point (h, θ) on the h(θ)curve. The height of the sample (usually ∼5 cm) is assumed to be too small to let the value of θ be affected by the vertical gradient of h that counters the gradient in the gravitational potential.
Even at this small scale, assuming the vertical variation in h to be negligible may not 20 be permitted. Liu and Dane (1995) developed a criterion to test the validity of the assumption for the general two-phase case, involving the densities of the wetting and the non-wetting fluid and the geometry of the experimental set-up. In an elegant analysis, Liu and Dane (1995)  in the curve's derivative, resulted in a curve more similar to van Genuchten's (1980) continuously differentiable shape. They calculated the smooth average curve by numerically treating the sample as a stack of thin slices of soil with identical properties, but with different values of h corresponding to their respective elevations. By summing the water contents of the slices for their respective values of h, the sample water 5 content could be found. For larger and possibly irregularly shaped averaging volumes V , the vertical extent of V needs to be accounted for in conjunction with porous medium heterogeneity within V . The procedure proposed by Liu and Dane (1995) can be adapted by assuming hydrostatic equilibrium within V , characterized by uniform H or by h V . Applying the 10 hydrostatic equilibrium condition H=h+x 3 to Eq. (11) gives: where A(x 3 ) is the horizontal area [L 2 ] of V as function of x 3 . Note that either H or h V is sufficient to fully characterize the distribution of h over V if the boundary of V is known. The corresponding average water content is (see Eq. 8): An upscaled h V (θ V ) relationship according to Eqs. (17) and (18) incorporates spatial heterogeneity and allows h to vary with elevation under hydrostatic equilibrium conditions. The relationship for macroscopic V (e.g., ranging from a soil layer within a field plot to an entire field) will be of little use to calculate actual flow, but by comparing HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential (e.g., the groundwater, or a stream) can be established with a more or less quantitative measure.

Macroscopically horizontal flow in a container filled with a heterogeneous porous medium
This is a generalization of the hypothetical case discussed by Gray and Miller (2004), 5 where a horizontal, closed rectangular container of length L [L], height B [L], and crosssectional area Ω [L 2 ] was filled with a porous medium. Two water-filled reservoirs in contact with the porous medium over the full cross-section at opposite ends at x 1 =0 and x 1 =L established fixed hydraulic head boundary conditions. A uniform layer of vertical extent b [L] and porosity n 2 was sandwiched between uniform layers with porosity 10 n 1 . The stripe with porosity n 2 tilted from the container bottom at x 1 =0 to the top at x 1 =L to make the distance c 2 [L] between the vertical midpoint of the layer and the container bottom: Gray and Miller (2004) derived an area-averaged pressure (equivalent to an averaged 15 h) over a vertical cross-section using the θ-weighted averaging underlying Eq. (11) above. If the integration in Eq. (11) were carried out over Ω instead of a volume (as explained below Eq. 12) and the head were converted to a pressure, the resulting equation would be identical to Gray and Miller's. According to the superposition principle (Sect. 2.1) that requires identical averaging manipulations for all components of H, 20 the gravitational head x 3 should have been averaged in the same way to arrive at the correct gradient in the averaged H Ω (with the subscript indicating the cross-sectional average). For n 2 >n 1 , the depth interval b centered on c 2 would overcontribute to the areally averaged θ-weighted h. Initially, Gray and Miller (2004) did not carry out the averaging of the gravitational head arguing that gravitation plays no part in horizontal flow. However, even in horizontal flow, the gravitational head contributes to H. Consequently, the calculated ∂H Ω /∂x 1 was erroneous, resulting in an obviously incorrect non-zero water flux while the water levels in the reservoirs at both ends of the column were the same.
In the next section of their paper, Gray and Miller (2004) developed an averaging equation for the gravitational potential that is the areal averaging equivalent of Eq. (12) above. The consistency between the averaging manipulations of h and x 3 was thereby restored, and the correct gradient of H Ω was obtained.
The corrected equations by Gray and Miller (2004) are only applicable to the spe-10 cific configuration of their hypothetical set-up, and their methodology rapidly becomes intractable for more complicated set-ups. In contrast, the superposition principle outlined above naturally leads to the generally applicable set of Eqs. (10)(11)(12), that reduces to Gray and Miller's (2004) equations for their special case. Thus, consistent application of the superposition principle resolves the paradox discussed by Gray and Miller 15 (2004), and the validity of Darcy's Law for the scale for which Darcy himself formulated it remains unchallenged.

Conclusions
The relationship between the hydraulic, pressure/matric, and gravitational heads is supported by the thermodynamic interpretation of fluid pressure and the gravitational po-20 tential. The superposition property of the heads constituting the hydraulic head translates into a consistency requirement for the upscaling manipulations of all heads. With the added constraint that the amount of energy must be conserved during volume integrations, the superposition property produces a set of consistent upscaling equations. Application of the equations to two cases demonstrates the usefulness of the equation set. In one case, an apparent paradox reported in the literature that threatened the validity of Darcy's Law was elucidated by demonstrating that it emerged from a violation HESSD 5,2008 Thermodymanics and upscaling of subsurface water potential