Statistical mechanics of unsaturated porous media

We explore a mean-field theory of fluid imbibition and drainage through permeable porous solids. In the limit of vanishing inertial and viscous forces, the theory predicts the hysteretic “retention curves” relating the capillary pressure applied across a connected domain to its degree of saturation in wetting fluid in terms of known surface energies and void space geometry. To avoid complicated calculations, we adopt the simplest statistical mechanics, in which a pore interacts with its neighbors through narrow openings called “necks,” while being either full or empty of wetting fluid. We show how the main retention curves can be calculated from the statistical distribution of two dimensionless parameters λ and α measuring the specific areas of, respectively, neck cross section and wettable pore surface relative to pore volume. The theory attributes hysteresis of these curves to collective first-order phase transitions. We illustrate predictions with a porous domain consisting of a random packing of spheres, show that hysteresis strength grows with λ and weakens as the distribution of α broadens, and reproduce the behavior of Haines jumps observed in recent experiments on an ordered pore network.


I. INTRODUCTION
Unsaturated porous media are ubiquitous in geophysical and industrial processes.They include, for example, soils partially filled with water, fuel cells, or oil reservoirs holding several gas and liquid phases.The principal challenge is to predict imbibition and drainage, viz.how fluid penetrates into-or emerges from-a porous solid matrix upon applying a macroscopic gradient in the pressure p of the wetting liquid relative to the ambient pressure p g .In general, porous solids delimit a complex void network consisting of a large number of small individual pores interconnected through several narrow openings (Fig. 1).When a porous sample traps one kind of liquid at a volume fraction θ , a simple measure of its partial filling state is the degree of saturation S ≡ θ/(1 − ν), where ν is the solid volume fraction of the dry sample.
Because the capillary energy of a pore depends upon the saturation state of its connected neighbors, the establishment of a local equilibrium derives from many-body interactions similar to those handled by statistical mechanics [1], such as lattice gases [2], neural networks [3,4], bird flocks [5], and spin glasses [6,7].In this context, we propose a mean-field theory of fluid retention in porous media partially filled with a connected wetting fluid.
In this regime, the porous network exhibits a nonlinear hysteretic behavior, whereby an applied macroscopic capillary pressure ψ ≡ p g − p delivers a higher saturation when ψ rises to expel the wetting liquid than when ψ is reduced to imbibe the sample.The relation S = f (ψ) then includes two limiting "main fluid retention curves" that describe, respectively, the imbibition of an initially dry porous solid S = f w (ψ) and the draining of a completely saturated sample S = f d (ψ).
Analyses capturing this phenomenon typically adopt one of two approaches.The objective of the first is to describe the porous medium on a scale large enough for practical applications.To that end, this traditional approach uses partial differential equations (PDEs) that incorporate retention curves as a constitutive law.The second approach exploits recent progress in direct numerical simulations and threedimensional x-ray tomography to observe geometry and behavior of a smaller sample.
To establish the constitutive retention behavior of the medium, the first approach models unsaturated porous media as a collection of independent pores described, for example, as a bundle of capillary tubes [8], or by analogy with other physical processes like magnetism [9][10][11] or neural networks [12].Notably, Poulovassilis [13] and Mualem [14] invoked the "ink-bottle" effect, whereby the total liquid content of a single bulging pore at a given ψ is greater if capillarity thwarts expulsion of the liquid, than it is when the liquid is drawn into the bulge.This observation recognized the role played by connected pores in setting the local saturation, but mainly attributed hysteresis to a distribution in pore size, rather than to the collective behavior that we now explore.
In the independent pore model, parameters are ultimately fitted to experimental data.For example, to characterize the water retention of soils, Van Genuchten [15] put forth a convenient expression with five parameters.The latter include the liquid volume fraction θ s at saturation, its residual value θ r after the liquid ceases to percolate between the two boundaries across which ψ is applied, a characteristic capillary pressure ψ a , and two exponents that are adjusted to reproduce the shape of the fluid retention curve.If this simple model does not capture data satisfactorily, other mathematical forms of the retention curve are chosen to minimize the number of additional empirical parameters needed to describe the entire range of S from saturation to complete dryness [16].In the terminology of wet granulation [17], this range begins with the "funicular" regime near saturation, where percolation allows ψ to be felt through most of the liquid.It then transitions to the "pendular" regime, where liquid congregates near grain contacts [18].Finally, liquid films condense or evaporate from surfaces FIG. 1. (Color online) Section from the x-ray computed tomography (CT) of an unsaturated porous soil aggregate, with superimposed qualitative identification of pore volumes v p , pore wettable surface area a p , and neck area a n .Continuous and dotted blue (white) lines show, respectively, air-water interfaces and boundaries between two connected water-filled pores.Continuous and dotted orange (light gray) lines mark, respectively, pore surface area in contact with water or air.Photo courtesy of Valérie Pot.within the porous solid [19].Hysteresis is typically handled by fitting different parameters for imbibition and drainage.
Traditionally, once established experimentally across a macroscopic sample, retention curves are used to predict the effective permeability K of the porous medium to the wetting fluid in terms of its saturated value K 0 , for example through the integral models of Burdine [20], Brooks and Corey [21], or Mualem [22,23].Retention curve and permeability are incorporated into PDEs like that of Richards [24], which combines Darcy's law and mass conservation to predict the evolution of S over space and time.Other treatments also include effects such as "non-Darcy" inertial behavior [25].The PDEs are then solved subject to boundary conditions for drainage or imbibition [26], which account for the relative roles of viscous, inertial, and capillary forces at and near the free surface of the porous medium.
Unfortunately, experiments show that, if a large sample is subject to a pressure drop imposed across distant boundaries, the resulting hysteresis is more complicated than based upon the two unique "main" retention curves S = f d (ψ) and S = f w (ψ).For example, if capillary pressure decreases before the sample drains completely, the overall degree of saturation follows a path in the S vs ψ diagram that invades the region between f d and f w [13,[27][28][29].Although the function f (ψ) appears continuous upon reversal, its derivative is not, and a new path S = f (ψ) is opened each time the rate of change of ψ switches sign.In general, this "return-point memory" [30] predicates the current state of a large sample on details of its past history, which complicates an already challenging numerical integration [31,32].
To address this difficulty, other treatments subsume hysteresis by introducing constitutive expressions that include state variables other than ψ and S, such as microscopic interface curvature [33], or by involving dynamic relaxation [34].Averaging procedures based on moments of the Boltzmann equation [35] or methods consistent with thermodynamics [36] are then employed to build macroscopic governing equations for the evolution of these variables in space and time from physics at the microscopic scale [37][38][39].However, to solve practical problems, the averaging still requires a measurable closure, such as a retention curve, to capture the constitutive behavior of the system.
The second approach begins with a detailed geometric description of the porous medium.Like direct numerical simulations, this method aims at reproducing the retention behavior by integrating on the pore scale PDEs that incorporate capillary, inertial, and viscous forces [26,[40][41][42][43].Other treatments operate on an intermediate scale to handle uneven fluid distribution, for example, in the pendular regime [44].To interpret the complex ramification of pore networks, a fertile literature invokes percolation theories that handle collective interactions naturally [45][46][47].However, the percolation framework, while fruitful in describing geometrical domain size and shape, does not directly account for surface energies and displacement work, unlike the approach outlined here.
Recent experimental techniques, such as x-ray microtomography, have the potential to inform numerical simulations by revealing the internal liquid distribution among pores within a solid matrix in detail [48][49][50][51][52][53].Here a challenge is to relate the complex microscopic geometry of triple contact lines where gas, solid, and liquid meet to the macroscopic behavior of the porous medium [54].Because the location and structure of these contact lines are affected by surface roughness and impurities that are difficult to discern or control, data interpretation and reproducibility are challenging [55].
As with other complex multiphase flows, the choice among the two approaches is mandated by overall system size.Because direct numerical simulations or three-dimensional experiments are rarely large enough to handle realistic applications, a formulation based on PDEs remains of practical interest.However, the complicated history-dependent hysteresis observed on macroscopic samples has called into question the meaningfulness of the retention curve and the judiciousness of describing the filling state of an unsaturated porous media by the degree of saturation alone [33,56].
Several fundamental questions therefore persist.Are the traditional retention curves S = f (ψ) meaningful?If so, at what scale?Can they be derived from porous geometry and surface energies?What causes sudden "Haines jumps" in capillary pressure and fluid speed that are observed as a pressure gradient or a flow rate are imposed on macroscopic boundaries of the system [57]?
To address these questions, we explore a framework that underscores collective interactions of pores in unsaturated permeable solids.Without parametric fitting, the theory predicts the behavior of the main fluid retention curves for a porous domain of known geometry and surface energies.To present the analysis as clearly as possible, we deliberately adopt the simplest mean-field statistical mechanics and illustrate it with generic examples.Specifically, we adopt an "Ising model" that characterizes the state of a single pore as being either full or empty of wetting fluid.Adjacent pores interact through the narrow openings or "necks" that connects them (Fig. 2), either by possessing a gas-liquid interface if their filling states are different or by not having one if they contain the same fluid.
We begin with a derivation of the theory.We then show how hysteresis of the main imbibition and drainage curves is associated with collective first-order phase transitions in the void network.To illustrate how the geometry of pores and necks affects hysteresis, we consider generic statistical distributions of these quantities, and we examine the collective behavior of a porous medium created by a dense packings of spheres.Finally, by identifying viscous forces as the mechanism that absorbs the "latent energy" released in the phase transitions, we predict recent observations of Haines jump by Armstrong and Berg [58] with an ordered pore network.

II. MEAN-FIELD THEORY
We consider a void space of mean volume fraction (1 − ν) delimited by internal surfaces of a porous solid, filled with two immiscible wetting and nonwetting fluids.For convenience, we refer to these fluids as a liquid and a gas, although the nonwetting fluid may also be a liquid.Although the theory could be extended to a hydrophobic situation such as mercury porosimetry [59,60], thereby reversing the sign of capillary pressure, we illustrate it with hydrophilic solids relevant to fuel cells and soils.
The void space consists of many pores of volume v p and wettable surface area a p interconnected to N adjacent pores of index i by narrow "necks" of interfacial area a n,i (Fig. 2).For simplicity, we adopt the Ising model [1] and assume that pores contain only gas or liquid, with known surface energies γ s , γ g , and γ gs of the solid-liquid, liquid-gas, and gas-solid interfaces, respectively.Consistent with this framework, we assign a binary filling state variable σ to each pore, whereby σ = +1 denotes a pore full of gas, while σ = −1 denotes one with liquid only.
This Ising framework ignores liquid films.Therefore, it does not capture their role in the "hydraulically connected transition zone" found in evaporative drying [61,62].However, it is relevant to "capillary pumping" [63].
We begin by evaluating the energy E(−σ → +σ ) that must be supplied to a pore to change its state from −σ to +σ .Without loss of generality, we calculate the amount E(−1 → +1) that is needed to empty out a pore initially full of liquid.This energy input has three contributions.The first two are independent of the state of neighboring pores.They include the volume work that integrates infinitesimal energies needed to counteract the net pressure (p − p g ) = −ψ resisting the substitution of a gas for a liquid of volume shrinking from v p to zero.The second is the energy required to lift the liquid from the wetted pore surface of area a p , p = (γ gs − γ s )a p = γ g a p cos θ c , where θ c is the static contact angle at the triple line where solid, liquid, and gas meet.The third and most notable contribution to E depends on the filling of all N connected pores of index i.If the adjacent pore i has no liquid (σ i = +1), then emptying the pore of interest involves the "exothermic" destruction of a gas-liquid interface.If instead pore i is full (σ i = −1), then the process entails the "endothermic" creation of a new such interface.Overall, the energy required is a sum over all adjacent pores, A similar argument shows that, for the opposite transformation (σ = +1 → σ = −1), all signs are flipped in Eqs. ( 1)-(3).Therefore, in general, With the convenient (and inconsequential) choice of a ground state at σ = 0 halfway between a completely filled or a completely empty porous domain, the energy of a single pore is then In the mean-field theory [1], σ i is approximated by the domain-average filling state σ , recognizing that the relative fluctuation in σ among adjacent pores is small if their number is large.Unlike problems of granular mechanics, what matters here is pore space contained outside the solid matrix.Therefore, in this article, an overbar denotes the volume average of any state quantity over the domain of total void volume V = v p .With the self-consistent mean-field assumption σ i σ , Eq. ( 4) then becomes where a n ≡ a n,i is the total cross-sectional area of all necks connected to the pore of interest.In a disordered porous material, v p , a p , and a n are random.For a crystal of identical pores, they instead possess discrete values.
In the framework of statistical mechanics, we consider a large ensemble of identical, yet distinguishable, copies of a porous sample.Each copy has a distinct distribution of filling states.Copies with similar total energy in Eq. ( 5) form a microstate with Bose-Einstein statistics (i.e., there is no limitation on the number of pores holding a given filling state).Then, to maximize the entropy of the system at equilibrium, the probability to find a pore at the filling state σ , subject to the constraints of a given ensemble-average energy and a normalized probability, conforms to the Maxwell-Boltzmann distribution [1], where β is a Lagrange multiplier and Z is the partition function The expected filling state in the pore of interest is then where ϕ denotes the ensemble-average of any state quantity ϕ over all copies of the sample.This expression suggests how the problem should be made dimensionless.First, we define a characteristic length ¯ 0 based on the mean dry volume of N p pores in the domain, in which the subscript 0 represents the dry void geometry.
Then the dimensionless capillary pressure can be written Similarly, the dimensionless specific neck cross section relative to pore volume is λ ≡ a n ¯ 0 /v p (11) and the specific wettable pore surface area is such that Eq. ( 8) becomes where β ≡ βγ g v p / ¯ 0 .[Because we only use ¯ 0 to make capillary pressure dimensionless, other convenient definitions are possible.For media with a wider distribution in v p , it would be prudent to adopt a ¯ 0 that is not as skewed toward larger populations of small pores as Eq. ( 9)].
In Appendix A we show that β is a very large number in the limit of small inertial and viscous forces.In other words, the porous medium exhibits "frozen disorder" [64].This lets us simplify the hyperbolic tangent in Eq. ( 13), which effectively becomes the Heaviside function H(ψ − ψ c ) jumping from −1 to +1 at a dimensionless threshold pressure, Then, in a domain where pores and necks are not uniform, the domain-average expected filling state is where F (λ,α) is the normalized joint distribution function by volume of λ and α in Eqs.(11) and (12), such that dv = V F dλdα (16) is the elementary pore volume distributed within the domain of total void volume V that has λ ∈ [λ,λ + dλ] and α ∈ [α,α + dα].For an ergodic system at equilibrium, σ = σ , (17) so that Eq. ( 15) can be written where I is the integral function and + and − are complementary nonintersecting regions filling the parameter space (λ,α) and satisfying Meanwhile, because σ = ∓1 for pores filled with liquid and gas, respectively, the domain-average σ is related to the liquid volume fraction through or, equivalently, S = (1 − σ )/2.A stable solution of Eq. ( 18) yields a point on the retention curve of θ vs ψ .

III. PHASE TRANSITION AND HYSTERESIS
We now illustrate how the theory predicts hysteretic retention curves and first-order phase transitions.The latter arise as a saturated sample abruptly drains most of its liquid or, conversely, as a dry one suddenly jumps toward saturation.To that end, we examine the simplest case of an ordered hypothetical "foam" with single-valued v p , a p , and a n , consisting of identical hollow spherical pores with circular necks carved from a solid on a regular cubic lattice.(A hexagonal close-packed crystal of solid spheres is analyzed in Appendix H. ) For this foam, because F is a normalized δ function located at (λ,α) (square symbol in Fig. 3), the two integrals on the right of Eq. ( 19) are unity if the δ function belongs to the domain of integration and zero otherwise.Therefore, their combination can only take one of three values (−1,0,+1).
To derive the resulting retention curves [Fig.3(a)], consider first an initially saturated sample with σ = −1 and ψ = 0 [dashed line in Fig. 3(b)], for which the entire (λ,α) space is occupied by − [Eq.(20)].At small capillary pressure, the δ function at (λ,α) remains within − until ψ reaches Then, because I = −1, ∀ ψ < ψ − , Eq. ( 18) is satisfied and, consequently, σ = −1 remains its stable solution throughout the range ψ < ψ − .However, for ψ > ψ − , σ = −1 ceases to be a solution, forcing the system to jump to the other state σ = +1 marked by an open triangle in Fig. 3.In short, the porous domain remains saturated until ψ reaches ψ − .A similar argument [Fig.3(c)] implies that a sample having initially σ = +1 and ψ → ∞ stays dry as ψ is progressively decreased within the range ψ + < ψ < +∞, where Because ψ + < ψ − , an initially saturated sample with σ = −1 transitions abruptly to the other solution σ = +1 of Eq. ( 18) as ψ increases beyond ψ − .Conversely, an initially dry sample jumps to saturation as ψ falls below ψ + .In short, the porous domain with single-valued λ and α undergoes a hysteresis loop marked by two abrupt phase transitions with for drainage and imbibition, respectively [Fig.3(a)].Note that the singularity of ∂S/∂ψ at ψ ± is a consequence of the meanfield and Ising assumptions.In a physical system, the derivative should not change sharply, but continuously.If one constructed a more elaborate Potts statistical mechanics [65] on groups of m pores, as suggested in Appendix B, then the greater number of available energy states would produce transitions with rounded edges.
Because these transitions arise as an external field (viz.the capillary pressure ψ ) passes through a critical point at which the total system energy per unit void volume H ≡ (1/V ) E(σ i ) is discontinuous, they are classified as "first-order" phase transitions [30].
In this ordered geometry, the "air-entry potential," i.e., the capillary pressure at which a saturated sample begins to drain, is ψ − = γ g (α cos θ c + λ)/ ¯ 0 .Meanwhile, the dimensionless separation ψ ≡ (ψ − − ψ + ) between the two transition pressures is a measure of hysteresis strength, Without geometrical disorder, Eqs. ( 11) and ( 25) thus imply that hysteresis strength grows as the neck area relative to the (2/3) power of the pore volume.

IV. COLLECTIVE BEHAVIOR IN A DISORDERED MEDIUM
Traditional treatments regard unsaturated porous media as a collection of independent pores, each acting as a separate domain, and each having an individual critical capillary pressure for imbibition and another for drainage.In that view, a pore contributes to decreasing the average retention of the whole medium by emptying once its own critical pressure is reached; conversely, upon wetting the medium, each pore fills up after ψ decreases below another, lower threshold [14].
We explore an alternative approach underscoring the collective nature of pore interactions seen in experiments  20) when the mesoscopic domain undergoes its draining phase transition at ψ − 4.98 and σ − −0.925 (S 0.96).The white dashed line is the wetting phase transition with ψ + 0.35 and σ + +0.820 (S 0.09).The solid and dashed lines (respectively L − and L + in Fig. 3) mark the corresponding transitions for a hypothetical porous medium with single-valued pore and neck sizes having the same mean λ and ᾱ. [66][67][68].Specifically, the self-consistent integral formulation of Eqs. ( 18)- (20) implies that the entire unsaturated porous domain contributes collectively to its equilibrium solution, rather than as a superposition of individual transitions.In our mean-field analysis, pores are interconnected to the bulk filling index σ through the necks they each possess.Therefore, the hysteresis of our retention curve is determined by the distribution function F , rather than by averaging retention curves of individual pores.
To show how such collective behavior arises in the statistical mechanics, we now consider a porous domain with frozen disorder embodied in a broader distribution F (λ,α), which we illustrate with a random dense packing of spheres obtained in numerical simulations (Fig. 4).Appendix B outlines how F is calculated from the diameter of spheres and the position of their centers.
Unlike the crystal example in Sec.III, the integral function I no longer takes on discrete values, but adopts instead a sigmoidal shape.Figure 5 illustrates the search for solutions to Eq. ( 18) by superposing I vs σ and the diagonal representing the ergodic condition σ = σ .In general, there can be one, two, or three solutions at intersections of the sigmoidal and diagonal lines in Figs.5(b) and 5(c).
To establish whether any solution is stable, we must first determine causality among state variables.On the one hand, the mean filling state σ determines the regions ± of integration of I or, in short, σ ⇒ I. Conversely, ergodicity implies that knowledge of the expected filling state σ = I leads to knowledge of σ , i.e., σ ⇒ σ .By determining whether the domain returns to a solution upon small excursions away from it, these causal relations indicate that, out of three solutions, the middle one is unstable, while the others at low and high filling are both stable.(Instability of the middle solution explains why we ignored σ = 0 in the example of Sec.III.) Consider an initially saturated domain ( σ = −1) without capillary pressure (ψ = 0) [Fig.5(b)].Here λ σ + ψ − α cos θ c < 0,∀ (λ,α), so that − represents the entire (λ,α) space shown in Fig. 4, and I = −1.This saturated state is a single solution of Eq. ( 18) represented by the lower left corner in Fig. 5(b).As capillary pressure is increased, the I curve shifts leftward, until one, then two, new intersections arise beside the lower left solution near σ ∼ −1.Yet, because this solution is stable, any other solution is ignored and the domain remains near saturation.However, as ψ increases, the sigmoidal curve eventually moves too far leftward to intersect the diagonal near saturation.At the pressure ψ − (solid triangle), the whole domain undergoes a first-order phase transition whereby σ changes sign and jumps to the dryer solution on the upper right (open triangle).Figure 5 sigmoidal curve at ψ = ψ + intersects the diagonal at only one point (open circle).At lower capillary pressures, the domain must then jump to its other solution closer to saturation (solid circle).
As with the crystal in Sec.III, the theory suggests that the transition pressure upon wetting is greater than that upon draining, ψ − > ψ + , so hysteresis arises once again.However, as Appendix C shows, the difference ψ ≡ ψ − − ψ + is always smaller with a broadly distributed F than with singlevalued λ and α, Meanwhile, the form of Eq. ( 20), which marks the border between the domains − and + in the (λ,α) parameter space, implies that the air entry potential ψ − shifts toward higher values as ᾱ increases, consistent with the ordered example in Eqs. ( 22) and ( 23). Figure 6 confirms these trends with retention curves calculated from artificial distributions F (λ,α) having a simple analytical form.In short, increasing the mean λ strengthens the hysteresis; i.e., it widens the gap between imbibition and drainage curves without shifting their midposition along the pressure axis [in Fig. 6(b), compare curves for larger λ with the base case].In contrast, increasing the mean ᾱ translates the midposition toward higher ψ without affecting hysteresis strength [Fig.6(a), as ᾱ is raised].Last, spreading F attenuates hysteresis (see curves for "wider F ").However, doing so along the α axis mitigates hysteresis more than along λ, to the point that hysteresis and both abrupt transitions disappear altogether with a sufficiently wide F α [Fig.6(a), retention curve marked "wider F α "].
Even in the example of Figs. 4 and 5 with monodisperse spheres, the retention curves exhibit a more gradual draining and wetting transitions than in Sec.III, as the F (α,λ) is no longer a δ function.More generally, a wider particle-size distribution (PSD) should induce greater spread in F and therefore less hysteresis and a more gradual transition roll-off.Then, for example, this theory suggests that a sand with relatively narrow PSD should have significant hysteresis.However, if the same sand also included fine particles, the strength of its hysteresis should be diminished.
Because the theory is built upon actual areas and volumes through λ and α, it tacitly accounts for interface deformations caused by a contact angle θ c = π/2, in addition to the explicit cos θ c appearing in Eq. ( 20).However, it is difficult to determine the distribution F (λ,α) with complicated gasliquid interfaces.Instead, it is more straightforward to find its counterpart F 0 (λ 0 ,α 0 ) on a dry sample (which is roughly equivalent to θ c = π/2).Then one can estimate how a contact angle = π/2 affects interface area and pore volume.As Appendix E shows, unless λ 0 is relatively large, values of θ c < π/2 mainly affect λ by raising it uniformly, as estimated in Eq. (E4), thereby exacerbating hysteresis strength.

V. CRITICAL DISORDER
Ji and Robbins [69] and Sethna et al. [30] emphasized the role of randomness and scale in disordered systems of interconnected sites.In this section, we suggest that the (b) From lowest to highest insets: "base case" for F α = F λ , with s λ = s α = 0.32 and λ 1.2 (solid transition lines); case of a "wider F λ ," with s λ = 0.42 (dashed transition lines); case of a "larger λ," where λ alone is increased to 1.7 from the base case (dotted lines).magnitude of geometrical disorder matters to the retention characteristics of unsaturated porous media.
Although these authors considered random fluctuations in applied field, rather than our frozen disorder in pore and neck geometry, their Hamiltonians resembled our Eq.( 4).Specifically, their H played the role of our ψ, their J ij was our γ g a n,i , and their random field f i (or h i ) was our γ g a p cos θ c .They then defined the parameter R as the ratio of the full width at half maximum of their Gaussian distribution of f i to the fixed magnitude of their J .Translated to our nomenclature, this relative randomness R represents approximately R ∼ = 2 √ 2 ln(2)s α cos θ c / λ, where s α is the standard deviation of α.
Notably, those authors found a critical relative disorder R = R c in the applied field below which a draining system undergoes a single abrupt "avalanche" similar to our collective phase transitions in Fig. 5.As they increased R toward its critical value, their magnetic equivalent of our S vs ψ curve became progressively less abrupt, eventually leading to a smooth curve beyond R c and to percolationlike behavior of the system [70][71][72][73][74][75].For R > R c , their simulations broke down into a myriad fractal domains similar to capillary fingering [76] and undergoing multiple avalanches with sizes exhibiting a power-law distribution [64,77].
Experiments conducted on artificial porous media with weak geometrical disorder [58,78,79] confirm the sharp transition behavior observed at small R [30,69].In our analysis, the critical disorder marks the disappearance of abrupt transitions, which arises as F α widens in Fig. 6(a).This occurs roughly when s α 0.72 and λ 1.2 or, equivalently, when R c cos θ c .Whereas the relatively ordered media shown in Figs.4-6 undergo relatively abrupt transitions as their relative disorder lies below R c , a more complicated fractal behavior of liquid penetration, albeit with a smoother retention curve, is expected as F widens along α.In short, experiments on weakly disordered systems may not behave like porous media with high disorder.

VI. MESOSCOPIC DOMAIN
Although Ji and Robbins [69] and Sethna et al. [30] noted avalanches of wide spatial extent below the critical disorder, experiments [58,80] and simulations [81] at relatively weak disorder suggest that invasion events only involve mesoscopic domains with limited number of successive nearest neighbors.To estimate the size of this domain, we calculate in Appendix D the probability Pr c (m) that the mth nearest neighbor has the same filling state as the original pore of index m = 0. Then a measure of the mesoscopic domain size is the value of m where Pr c (m) approaches zero.If tortuosity of the porous medium is known [82,83], then m can be converted to a correlation length.We find As expected, Pr c decays with m ever more slowly as the system becomes either dry or saturated.Equation ( 27) also indicates that the mesoscopic domain is smallest for S = 1/2.Meanwhile, the integral scale is a more objective measure of the correlation neighbor index m c at which Pr c has decayed substantially.For the dense spherical packing with retention behavior in Figs. 4 and 5, m c − 25 and m c + 10 at the draining and wetting transitions S − and S + , respectively.Here imbibition occurs in a smaller mesoscopic region than drainage.Because both m c − and m c + just fit within the simulation shown in Fig. 4, the inset illustrates how small a mesoscopic domain of identical spheres can be.Such small mesoscopic size challenges volume averaging.
In practice, the latter also requires a limited magnitude of capillary pressure gradients, particularly if these are directed along the flow, u • ∇ψ > 0, an unstable configuration that triggers the onset of fingering [84], like unfavorable pressure gradients in a boundary layer.Because our framework, like equilibrium thermodynamics, does not directly involve space or time, it needs insight on ramification of the mesoscopic domains upon which it is built.For example, our estimate of the correlation index in Eq. ( 28) could be refined using percolation theory.Although the latter also invokes statistical mechanics, its objectives differs from ours.By focusing on pore geometry, percolation theory predicts how size and shape of connected liquid clusters obey power laws in the difference between the degree of saturation S and a critical value beyond which the first infinite cluster spans the entire system [85].Here our aim is instead to predict the hysteretic relation between S and applied capillary pressure ψ, which also exhibits critical phase transitions that are no longer purely geometrical.Nonetheless, for a gas-liquid mixture, percolation theory complements our analysis by predicting the minimum S at which individual pores feel the gradient ∇ψ imposed on the liquid across distant boundaries [86].For mixtures of two liquids [58], ψ propagates instead throughout the system, as both components can sustain negative pressures.In that case, percolation theory could lend insight into the geometrical structure of liquid clusters within the mesoscopic domain, thus refining Eq. ( 28).[However, unlike percolation models tracking only the minority phase, our analysis must involve interpenetrating clusters of both phases, for example, yielding a mesoscopic domain size in Eq. ( 28) that is symmetric in S and 1 − S.]

VII. HAINES JUMPS
Because our equilibrium theory ignores time and gradients, it considers infinitesimal increments in capillary pressure and fluid saturation up to irreversible first-order phase transitions of drainage or imbibition, thereby only addressing the limit of negligible inertial and viscous forces.Nonetheless, we take here a first step toward nonequilibrium by considering rapid fluid rearrangements called "Haines jumps" [77], which we associate with dissipation of the latent energy released in phase transitions.
For inertial forces to be significant [87], the Weber number We ≡ ρ u 2 ¯ 0 /γ g typically exceeds unity [88], where u is a characteristic interstitial flow speed of a wetting fluid of density ρ .For example, inertial forces become significant only if water seepage during the wetting or draining transitions produces a speed u 1 m/s in a sand bed with ¯ 0 ∼ 50 μm.Neglecting viscous forces mandates lower capillary numbers Ca ≡ μ u /γ g < 10 −5 [89] based on the dynamic viscosity μ of the wetting fluid and, therefore, much lower speeds.(An alternative definition Ca * lets rapid fluid mobilization occur when Ca * ∼ O(1) [90].)Then we expect viscosity to play a role during Haines jumps seen in experiments [52,80,91] and numerical simulations [55].Conversely, we do not expect viscosity to matter during the slower, reversible redistribution of fluid [92], which Berg et al. [52] also observed for a substantial fraction of the displaced volume, and which we interpret as the reversible approach to a phase transition.
Their experiments had a viscosity ratio μ ≡ μ g /μ of the nonwetting and wetting fluids that renders viscous fingering unimportant [76].Because their artificial pore network had no geometrical disorder (R R c ), phase transitions should be sharp.As expected, sudden Haines jumps arose as water was drained by injecting immiscible nonwetting decane at controlled flow rates.Unlike the air-water system in which the gas cannot hold tensile stresses, negative pressures percolated through both liquids at any S. Therefore, our theory should hold even at low water saturation.Because the porous solid medium was etched into glass on a hexagonal pattern with uniform pores and necks, its retention curves possessed the Heaviside shape in Fig. 3. From its known geometry and contact angle (quoted in the caption of Fig. 7), we calculate λ 0.627 using Eq.(E4).Such relatively small λ should also produce modest hysteresis.
To predict the behavior in Haines jumps, we first calculate the total energy H in a unit volume of the mesoscopic domain by summing Eq. ( 5) over all pores.For this monodisperse network, it is, in dimensionless form, Then the draining first-order phase transition produces a dimensionless "latent energy" per unit volume Because the overall energy in Eq. ( 29) is proportional to σ to leading order, the average rate of latent energy produced is approximately −(∂S/∂t)|L |γ g / 0 during a draining phase transition in which S remains uniform in the mesoscopic domain.It balances the combined energy dissipation rate μ u 2 (1 − ν)S/K + μ g u 2 g (1 − ν)(1 − S)/K g from viscous forces on both fluids averaged over the whole mesoscopic domain.Because u ∝ x, the domain average of u 2 is (4/3) ū2 .In short, the energy balance is Substituting the continuity Eq. ( 31) and introducing the ref- 2 and time τ 0 ≡ 0 /u 0 , Eqs. (31) and (32) are, in dimensionless form where t ≡ t/τ 0 and u ≡ ū /u 0 , ∂ ln S/∂t = −2u /m c (33) and where μ ≡ μ g /μ is the ratio of viscosities of the nonwetting and wetting fluids.Equation (34) predicts that the mean dimensionless velocity peaks at a value u max given by ln u max −0.036 20 ln μ 2 − 0.4206 ln μ − 2.745 in the range 10 −3 < μ < 10 3 , at a degree of saturation given by ln S max −0.004 72 ln μ 2 − 0.0849 ln μ − 0.7340.In the experiments of Armstrong and Berg [58], the maximum Weber number is ∼1.2 × 10 −3 , which is too low for inertial forces to matter.However, if Armstrong and Berg [58] had used air (μ 0.02) rather than decane (μ 0.96) as a nonwetting fluid, the peak speed should have been about 2.7 times faster, yet not large enough for inertia to become important.
To solve this problem, we first eliminate u from Eqs. ( 33) using (34) and obtain an ODE for S in terms of t .We then impose S max as an initial condition and solve the ODE both forward and backward from the time of peak mean speed.As Fig. 7 shows, results compare well with experiments without resorting to parametric fitting.In particular, they reproduce the different behaviors that Armstrong and Berg [58] reported before and after the peak, namely a rapid rise in speed, followed by more gradual deceleration.This drawn-out approach to a new drained equilibrium, which is revealed by model and experiment, gives the impression that the system reaches a residual degree of saturation ∼0.1 on short time scales.
From the peak speed u , we estimate the magnitude of the largest expected drop δψ max in capillary pressure by integrating its gradient across the mesoscopic domain, δψ max ∼ m c 0 μ ū (1 − ν)/K .Using Eq. ( 34), this is, in dimensionless form, max ].Then we expect capillary pressures to vary in the range ψ − ± δψ max during a Haines jump.From γ g of the water-decane interface [94], we find ψ − ± δψ max 15 100 ± 1700 Pa, which is consistent with the pressure jump of 15 790 ± 2820 Pa that Armstrong and Berg [58] reported.

VIII. DISCUSSION
Our analysis suggests that hysteresis is the natural collective behavior of a mesoscopic domain exhibiting phase transitions [30], as it is in other instances where statistical mechanics is useful, such as magnetism [95] or shape-memory alloys [96].In this view, it is not possible to subsume hysteresis by introducing another state variable.
The principal question is to what mesoscopic scale this collective hysteretic behavior applies.Recent experiments on a single plane layer of sintered glass beads [80] or in more complicated three-dimensional Berea sandstone [52] indicate that sudden collective rearrangements in liquid distribution only involve a few near -neighbors, thereby implying that the mesoscopic domain size should be relatively small, as Eq. ( 28) implies.
As derived, our mean-field theory predicts that, at the mesoscopic level, the region bound by the two curves in (ψ ,S) state space cannot be invaded upon a reversal of capillary pressure.Nonetheless, practical applications often stage porous media that are larger than the scale in Eq. ( 28).For such macroscopic systems, a counterexample based on the Preisach [97] model (Appendix F) shows how a porous medium with inhomogeneous saturation can produce such an invasion in (ψ ,S) state space.Because, in general, there is not a unique combination of degrees of saturations S i that produces an overall S = S i χ i in a medium composed of distinct mesoscopic domains with volume fraction χ i , any measurement (ψ , S) on a macroscopic sample that falls within the two main retention curves cannot describe the state of the system unambiguously.In other words, the existence of points within the main curves in (ψ ,S) state space implies a sensitivity to past conditions, which jeopardizes experimental reproducibility and solution unicity.Therefore, to sidestep these difficulties, one should avoid measuring retention curves with a pressure ψ imposed across distant boundaries, and one should integrate PDEs like Richards' equation [24] at the mesoscopic scale [98].
However, it is yet unclear whether this mesoscopic approach is sufficient to proscribe an invasion of the region within the two main curves, thereby guaranteeing experimental reproducibility and restoring unicity of a numerical solution.At present, our mean-field theory attributes randomness to geometry alone.However, because an inhomogeneous macroscopic sample may be composed of individual mesoscopic domains with different transition pressures, local fluctuations in the applied field ψ could arise even if a steady capillary pressure is imposed across a large sample.Such fluctuations would effectively raise the relative randomness R arising from Eq. ( 4).Therefore, although our analysis predicts retention curves that are reversible until a phase transition occurs, there may be irreversible precursors to the main phase transition.If such precursors existed for R < R c , they could lead to an invasion of the main retention curves upon a reversal in applied pressure.Future numerical simulations of Eq. ( 4) on porous media with a distribution F (λ,α) should establish the relative roles of geometrical randomness, fluctuations in applied pressure, and domain size.
Unlike the problem of vapor sorption on solid surfaces [99], in which ψ is written in terms of vapor pressure by eliminating interface curvature between the Kelvin and Young-Laplace Eqs.[100,101], capillary pressure appears explicitly through its volume work in Eq. ( 4).Therefore, vapor pressure is not as crucial here as it is in analyses of hysteretic isotherms of sorption on nanometric porous structures [102][103][104][105].
However, for our description to hold, negative pressure must be transmissible through a connected liquid.While this assumption holds when wetting and nonwetting fluids are both able to sustain a tensile pressure (e.g., the experiments of Armstrong and Berg [58]), a gas phase restricts applications to a large enough S that allows ψ to percolate across the mesoscopic domain.In the pendular regime where narrow liquid bridges congregate across necks surrounded by gas on both sides, saturation can be too low for ψ to connect, thus producing the residual volume fraction θ r observed in experiments [106].Nonetheless, morphological observations on spherical packings at the pore scale [107] indicate that liquid arranges in pore-filling clusters for S as low as 0.15, followed by a percolation threshold for long-range connectivity around S 0.2.At lower liquid volume fractions, pendular bridges do not disappear until desorption takes place.
In short, because pore energy in Eqs. ( 1)-( 4) ignores pendular bridges forming around necks at low saturation, we do not expect our mean-field theory to capture the transition from dry to wet as quantitatively as its converse from wet to dry, unless the predicted degree of saturation S + at ψ + is high enough to uphold the form of Eq. ( 4), or unless both fluids can sustain tensile pressures.Nonetheless, it may be possible to extend our analysis to lower saturation of a gas-liquid system by combining it with statistical mechanics of liquid sorption on nanoscopic surfaces [108][109][110][111]. Similarly, pores of the "F type" [61] coated with a residual liquid film could be handled in this framework by considering other energy states and modifying Eq. ( 4) accordingly.
Finally, our analysis has restricted attention to immutable geometries in which mechanical forces do not contribute to the system's energy.However, pressure forces exerted on the porous medium could produce microscopic rearrangements, for example, in unsaturated soils [112,113].In this case, the statistical mechanics should be refined to incorporate mechanical energy.Irreversible changes to the geometry would then lead to another mechanism for invasion of the main retention curves that is not yet included in our analysis.

IX. CONCLUSIONS
We presented a mean-field statistical mechanics of porous media filled with two immiscible fluids, one wetting and the other nonwetting.The model assumed that interstitial voids can transmit tensile capillary pressure applied on domain boundaries.Whereas such connectivity is readily achieved if both fluids are liquid, this assumption requires saturation above a percolation threshold if the nonwetting fluid is a gas.Although the theory requires no fitting parameter, it ignores phase trapping and liquid films in this case.
To avoid complicated derivations, we used a mean-field analysis based on the simplest Ising assumption that the void space is made up of pores connected to their neighbors through narrow necks, where fluid interfaces typically reside, and that pores are either full or empty of wetting fluid.We derived the energy E of an individual pore in terms of its filling state σ and average σ in its surroundings.The form of E prescribed how the applied capillary pressure ψ should be made dimensionless with mean pore size and interface energy of the two fluids.
With a simple example, we showed that a saturated mesoscopic sample having a disordered void space drains progressively as ψ rises, until it experiences a collective first-order phase transition that empties its wetting fluid abruptly.We then attributed hysteresis of the retention curve to another phase transition that the sample experiences as ψ is subsequently returned to low values.
We calculated the retention curve from statistical moments of the porous geometry, represented by the specific pore wettable surface area α, and the specific neck interfacial area λ, both made dimensionless relative to individual pore volume and average pore size.Having showed how these two parameters could be evaluated in spherical packings, we found that the retention curve is shifted to higher capillary pressures as ᾱ grows and that its hysteresis loop gets wider as the volume-averaged λ rises or as the hydrophilic contact angle decreases.We also found that, like other mesoscopic systems undergoing phase transitions, hysteresis is inevitable, but could disappear with a sufficiently broad distribution of α.
Finally, we showed that a velocity time history recorded in Haines jumps could be attributed to viscous dissipation of the latent energy released in the collective first-order transition predicted by the theory.

APPENDIX A: LAGRANGE MULTIPLIER
We show that the dimensionless Lagrange multiplier β in Eq. ( 13) is typically very large.As a scale for pore energy fluctuations, 1/β depends on fluid speed u .At vanishing u , energy fluctuations are only due to thermal agitation.They then grow with u as they do in agitated granular gases [114], thus possibly reducing β.
To evaluate β as u ∼ 0, we consider a reversible approach to the draining phase transition.(A similar calculation applies to its wetting counterpart.)Because at least one liquid phase is involved, the speed of sound is high enough for a typical mesoscopic domain to reequilibrate quickly and thoroughly to any small reduction dS < 0 in the degree of saturation.The resulting change in dimensionless total entropy ℵ t , which involves all active molecular degrees of freedom, has three contributions, dℵ t = β dH + dℵ g − dℵ .The first is from the change dH in dimensionless total energy [1], the second is the entropy input of nonwetting fluid dℵ g = −[ρ g 3 0 ŝg /(kM g )]dS, and the third is the corresponding entropy output of wetting fluid, dℵ = −[ρ 3 0 ŝ /(kM )]dS.In these expressions, k is Boltzmann's constant, and ρ g , ŝg , and M g are, respectively, the density, absolute molar entropy, and molecular weight of nonwetting fluid, while ρ , ŝ , and M are the corresponding quantities for the wetting fluid.Because the process is reversible until the phase transition, dℵ t = 0. Combining these relations, In typical situations, the term in square brackets, which we call A, is very large.For the air-water system at standard temperature with ŝ 70 J/mol K and ŝg 194 J/mol K, A 3 × 10 11 with 0 = 1 μm.If decane is the nonwetting fluid (ρ g = 730 kg/m 3 , ŝg 425 J/mol K), then A 10 11 .In the spirit of the mean-field theory, we estimate the total energy from Eq. ( 29) as H Toward the draining phase transition where |∂ψ /∂S| 1 and S 1, β ∼ A/2 λ is therefore very large.A consequence is that the retention curve should be insensitive to thermal temperature, as typically observed.

APPENDIX B: SPHERE PACKINGS
We consider the porous void space formed by random, dense packings of spheres.Studies of such granular assemblies often focus on compaction of the solid [115], rather than on the statistical distribution of voids and openings [116].For example, Edwards and Oakeshott [117] applied statistical mechanics to powder compaction.However, the mean solid volume fraction ν matters less to the retention curve than it does to the mechanics of granular materials.Whereas a small increment in ν can induce jamming in random packings [118], such increment does not affect the void space as much and, consequently, the retention curve is not changed substantially.
To our knowledge, the joint distribution F (λ,α) has not yet been studied.However, unless the PSD is wide, Luchnikov et al. [119] identified necks and pores with a Voronoi-Delaunay partitioning of void space in random packings of spheres.As a first attempt to characterize F , we exploit their Delaunay triangulation, which subdivides space into irregular tetrahedra having vertices on the centers of spheres of diameter d j .This triangulation is unique, as long as centers lie in "general position"; i.e., no three centers are aligned, no four centers lie on a plane or a circle, and no five centers are arrayed on a sphere [120].
Here necks have planar cross sections obtained by removing three circular sectors from each tetrahedral face.Similarly, the pore of index i is what remains of a tetrahedron of volume v t i after excision of N = 4 spherical sectors centered on each apex j with solid angles j subtending the opposite face.Then pore volume and area are, respectively v p i = v t i − j d 3 j /24 and a p i = j d 2 j /4 for j = 1,N .Although triangulation is a robust method for identifying unambiguously pores and necks in a packing of spheres, it only allows N = 4 near neighbors around an individual pore.Because accuracy of the mean-field theory grows with N [1], one could contemplate a Potts model [65] grouping m adjacent tetrahedral pores into a new unit cell with N = 3m + 1 external necks and 2 m possible filling states, each with a different unit cell energy E. The partition function in Eq. ( 7) would then include 2 m terms.Unfortunately, there would now be three kinds of dimensionless ratios arising from the energy of a unit cell, namely m parameters similar to α, another m similar to λ, and m − 1 independent ratios of pore and unit cell volumes.For example, a unit cell with two pores of index 1 and 2 would have seven necks, four possible filling states (σ 1 ,σ 2 ) = (+1,+1), (+1,−1), (−1,+1), or (−1,−1), and a five-dimensional joint statistical distribution of independent geometrical parameters.While being more accurate, this higher-order mean-field treatment would trade off the relative simplicity that we have exploited with m = 1.However, it would make it straightforward to handle rare pores that a triangulation subdivides artificially, such as a five-sphere pyramid with a square base.For soils with more complicated pore geometry (e.g., Fig. 1), techniques inspired from Monga et al. [121] could be developed to calculate F .

APPENDIX C: HYSTERESIS STRENGTH
Consider a R 2+ parameter space (λ,α), like that shown in Fig. 4. From Eq. ( 20), − and + are divided by the straight line of slope σ / cos θ c and intercept ψ / cos θ c , and are therefore complementary (i.e., + ∪ − = R 2+ ; + ∩ − = ∅).Meanwhile, integrals in Eq. ( 19) can be interpreted as "masses" contained within each domain with "surface density" F .In that view, the "center of mass" of the entire plane has coordinates Consider now a porous domain with broader surface density F (Fig. 4).Without loss of generality, assume that it begins saturated ( σ = −1).Because the magnitude | σ− | of the transition filling state is necessarily <1, the slope of the dividing line between − and + at phase transition is not as steep as L − .In addition, this dividing line can no longer pass through the center of mass, since doing so would imply an equal "mass" on both sides, which, from Eq. ( 19), would lead to I = 0.However, because σ− = 0, I = 0 could not be a solution to Eq. (18).Instead, the dividing line must be further displaced downward, so the − domains captures more mass, bringing σ− to the negative sign that we expect.In short, there are two reasons why the new transition ψ − is lower with a broader F .First, the slope of the dividing line is not as steep, thus reducing the intercept ψ / cos θ c .Second, the latter moves farther down to satisfy Eq. (18).
A similar argument shows that ψ + increases with a broader F .Overall, if hysteresis strength is measured as the difference between the two capillary transition pressures, it is always weaker with a broad distribution of pore and neck areas than the single-valued case in Eq. ( 25), hence justifying Eq. ( 26).

APPENDIX D: CORRELATION PROBABILITY
We calculate the probability that the mth nearest neighbor of a pore holds the same filling state, and we extract at which nearneighbor index m this probability has decayed substantially.We restrict attention to the case where all pores and necks have a single size.(For complicated geometries, Monte Carlo simulations constrained to uphold the overall voidage could refine predictions, albeit without the benefit of simplicity.) In a system with uniform pore volume, the probability for a pore to hold a filling state σ of the same sign as the volume average σ is (D1) (For example, if σ = 0, pores have equal probability to hold one or the other state; for | σ | = 1, all pores have the same sign.)Equivalently, the probability Pr(σ ) for the original pore to hold the state σ is Pr 0 if σ > 0 or 1 − Pr 0 if σ < 0. First, consider a pore filled with liquid (σ = −1).If σ > 0, then the probability for its first-nearest neighbor to have the same filling state is Pr 1 = 1 − Pr 0 .If instead σ < 0, then Pr 1 = Pr 0 .A similar argument applies to a gas-filled pore.In general, Pr 1 (σ ) = (1 + σ σ )/2.In the spirit of the mean-field theory, whereby the energy of a pore is only affected by the state of its first neighbors, Pr 1 is generalized recursively.Then, the probability for the mth neighbor to have the same σ is where σ can be either −1 or +1.Therefore, the probability to have a mth neighbor of the same filling state as the original pore is Pr c = Pr m (+1)Pr(+1) + Pr m (−1)Pr(−1).Overall, substituting S for σ with Eq. ( 21), )

APPENDIX E: ROLE OF CONTACT ANGLE
Because a neck involves a gas-liquid interface if σ changes sign across it, its area can deform away from the dry opening cross section any time the contact angle θ c differs from π/2.In principle, capillary pressure could also modify the interface shape.However, because effects discussed in this article occur in a relatively narrow range of ψ around ∼γ g / ¯ 0 , and because necks are small, ψ should contribute negligibly to interface distortion.Similarly, pore size typically makes the Bond number Bo ≡ g ¯ 2 0 (ρ − ρ g )/γ g too small for the gravitational acceleration g to matter.
Ignoring details of interface geometry [54], we estimate the deformation for θ c = π/2 by analogy with a cylindrical capillary, in which the meniscus is small enough to be nearly spherical.In this analogy, the ratio of neck interface area a n and its dry cross section a n 0 is similar to the ratio of meniscus area to capillary cross section, where the subscript 0 represents "dry" geometrical quantities measured in the absence of liquid.The interface is equally distorted if gas resides on one side or the other.If instead a neck connects two pores with the same phase (gas or liquid), then its area does not matter to pore energy, and may be taken to satisfy Eq. (E1).Then, in this estimate, all necks have greater cross section than the dry value a n 0 by a common factor 1 < a n /a n 0 < 2. Similarly, the deformation of a neck matters to the volume v p in Eq. ( 1), only if it connects two pores of a different phase.Therefore, the contribution of a single neck to the difference (v p − v p 0 ) is either zero or a volume increment of magnitude v p .Exploiting a similar capillary approximation than toward Eq.(E1), and assuming for simplicity that all N necks around a given pore have the same area, interface distortion amounts to Consider a pore filled with liquid (σ = −1).Here distortion of gas-liquid interfaces shrinks its dry volume v p 0 by an amount σ N = v p contributed by its N = necks connected to gas-filled neighbors.A similar argument applies to a gas-filled pore, where N = necks to liquid-filled neighbors swell its dry volume, again by the amount σ N = v p .In all cases, using results in Appendix D, we have N = /N = 1 − Pr 1 , so that (v p − v p 0 )/v p 0 = σ N v p (1 − σ σ )/2 = N v p (σ − σ )/2.On average, (v p − v p 0 )/v p 0 δ(σ − σ )/2. (E3) Therefore, by the ergodic condition σ = σ , this correction in average pore volume vanishes whenever δ is constant, for example with single-valued dry neck area and pore volume.It also vanishes in the general case if, as expected, the domain statistics of a 3/2 n 0 /v p 0 are uncorrelated with those of σ .However, the mean absolute excursion |v p − v p 0 |/v p 0 δ(1 − σ σ )/2 does not vanish, but it is usually small, since λ 0 is typically ∼O(1) and the terms in parentheses in Eq. (E2) is 2/3 √ Nπ,∀ θ c .For example, tetrahedral pores in a hexagonal close packing have a (E4)

APPENDIX F: MACROSCOPIC INHOMOGENEITIES
To illustrate the role of inhomogeneities, we consider a Preisach model [97] for a hypothetical medium consisting of two adjacent mesoscopic subdomains a and b occupying different parts V a = χ a V and V b = χ b V of the total pore volume V = V a + V b with distinct mean degrees of saturation S a and S b , but subject to the same dimensionless applied pressure ψ .For clarity, the entire medium has single-valued pore and neck sizes.Therefore, a and b share the transition pressures ψ + and ψ − (Sec.III).
Figure 8 shows the response to generic initial conditions.If the two subdomains a and b have an initial capillary pressure ψ 0 < ψ + (saturated), or ψ 0 > ψ − (dry), they both experience the same saturation history and effectively behave as a single domain.If instead they start at an intermediate pressure ψ + < ψ 0 < ψ − , then overall saturation depends on initial conditions S a 0 and S b 0 : If an initially saturated a (S a 0 = 1) is juxtaposed with a dry b (S b 0 = 0), then a stays saturated as ψ rises from ψ 0 until the phase transition at ψ − .This produces a horizontal segment on the joint retention retention curve at an apparent mean S = θ/(1 − ν) = χ a .Conversely, if a starts dry while b is saturated, then the retention curve features a similar segment at S = χ b = 1 − χ a .In either case, these segments invade the region in (ψ ,S) state space between the main curves f d and f w .Similar invasions arise with a more general distribution F (λ,α), with positions in (ψ ,S) that are sensitive to the initial inhomogeneous distribution of liquid.

FIG. 2 .
FIG. 2. (Color online) Pore geometry characterization for a porous assembly of identical spheres.(a) Single pore wedged between an irregular tetrahedron of four spheres.The contribution of the top sphere to the wettable pore area is orange (light gray).Three spheres in the foreground delimit the dry neck cross section of index i = 1.Cartesian coordinates of the four spheres relative to their diameter and measured from their barycenter are (x,y,z) = (0.05,−0.50,0.26);(−0.50,0.05,−0.36);(−0.14,0.46,0.51);(0.58,−0.01,−0.42).For this pore, λ 0 5.23 and α 0 3.20.(b) Section through centers of the three darker spheres showing traces of three of the four dry neck cross-sectional area a n,i , wettable pore surface area a p 0 and pore volume v p 0 .

FIG. 3 .
FIG. 3. (Color online) (a) Retention curves of degree of saturation S = (1 − σ )/2 vs dimensionless capillary pressure ψ for a hypothetical foam of identical spherical pores on a square lattice of solid volume fraction ν = 0.3, having single-valued λ 0 = 1.530 . . .and α 0 = 4.735 . . . .For the contact angle θ c = 50 • of this example, λ 1.73 and α 4.74 (Appendix E).Arrows show directions of changes in ψ .(b) Draining transition in (λ, α) space.For this foam, F is a δ function centered on the square symbol.The dashed line marks the border in Eq. (20) between domains + below the line and − above, with ψ = 0.The red (light gray) line L − is the corresponding border for the draining phase transition, which occurs when ψ , rising in the direction of the red arrow, reaches ψ − .(c) As ψ decreases in the direction of the blue (dark gray) arrow, it eventually reaches the wetting transition at ψ + 1.31.

FIG. 4 .
FIG. 4. (Color online) Contour plot of F (λ,α) for the random packing shown in the inset with 20 000 spheres forming 123 914 pores at ν 0.604 (detail in Fig. 2).The Delaunay triangulation outlined in Appendix B yields ¯ 0 /d 0.389, λ0 3.81, and ᾱ0 3.55.For a contact angle θ c = 50 • , λ 4.32.The open square marks the "center of mass" of F defined in Eqs.(C1) and (C2).The white solid line of slope σ − / cos θ c and intercept ψ − / cos θ c is the boundary between the regions + and − in Eq. (20) when the mesoscopic domain undergoes its draining phase transition at ψ − 4.98 and σ − −0.925 (S 0.96).The white dashed line is the wetting phase transition with ψ + 0.35 and σ + +0.820 (S 0.09).The solid and dashed lines (respectively L − and L + in Fig.3) mark the corresponding transitions for a hypothetical porous medium with single-valued pore and neck sizes having the same mean λ and ᾱ.
FIG. 5. (Color online) (a) Imbibition and drainage retention curves charting degree of saturation vs dimensionless capillary pressure for the distribution F (λ,α) in Fig. 4 and θ c = 50 • .(b) Integral I and domain-averaged expected filling state σ vs σ as ψ is increased from zero (I, black line) to its value ψ − 4.98 at the draining phase transition, beyond which the filling state jumps from σ− −0.925 (or S − 0.963, solid triangle) to σ +1 (open triangle), as I (left sigmoidal line) no longer intersects the diagonal to satisfy Eq. (18).(c) The corresponding graphs of I and σ vs σ as ψ is decreased from +∞ to ψ + 1.31 at the wetting transition (open to solid circles) with σ+ 0.820 (or S + 0.090).

FIG. 6 .
FIG. 6. (Color online) Dependence of retention curves S = S(ψ ) on F (λ,α), illustrated with artificial normalized γ distributions of the form F = F λ × F α , where F λ ∝ (λ − d λ ) a λ c λ −1 exp{−[(λ − d λ )/b λ ] c λ } for λ > d λ and zero otherwise, and F α ∝ (α − d α ) aα cα −1 exp{−[(α − d α )/b α ] cα } for α > d α and zero otherwise.Insets show the corresponding contours of F (λ,α).Vertical arrows and lines mark, respectively, wetting and draining phase transitions.(a) White arrows show how hysteretic retention curves are affected by an increase in the mean value of α from ᾱ 1.7 to 4.7 (respective phase transitions represented by dotted and dashed lines), while keeping λ 1.7 and the standard deviations s λ and s α of F λ and F α at 0.32.The retention curve without hysteresis (black line without a phase transition) is obtained from F in the "no hysteresis" inset having ᾱ 4.7, a wider s α = 0.72, and the same F λ as above.(b)From lowest to highest insets: "base case" for F α = F λ , with s λ = s α = 0.32 and λ 1.2 (solid transition lines); case of a "wider F λ ," with s λ = 0.42 (dashed transition lines); case of a "larger λ," where λ alone is increased to 1.7 from the base case (dotted lines).
Similarly, the wetting transition has L + = +ψ + − α cos θ c .Substituting transition pressures in Eqs.(22) and(23), both transitions have the same volumetric latent energy,L = −λ.(30)Therefore,they are both "exothermic" and should occur spontaneously, but irreversibly.Such latent energy is absorbed by viscous dissipation in the two fluids.To model this mechanism, we consider an open rectilinear mesoscopic domain of uniform saturation with coordinate −m c 0 < x < m c 0 , involving m c neighbors from Eq. (28), as Armstrong and Berg[58] observed.In a transition, volume conservation ∂S/∂t = −S∂u /∂x relates the gradient in water velocity u to temporal variations in S, and it binds decane and water velocities through u g = −u S/(1 − S).Taking u = 0 midway through the domain, the ODE integrates to u = −x∂ ln S/∂t or u = 2 ū x/(m c 0 ), where ū = −(m c 0 /2)∂ ln S/∂t (31) is the domain-average speed of water.Meanwhile, viscous forces exerted by water and decane in a unit pore volume are μ u (1 − ν)S/K and μ g u g (1 − ν)(1 − S)/K g , where K and K g are Carman-Kozeny permeabilities corrected for the respective incomplete degrees of saturation S and (1 − S) of the wetting and nonwetting fluids (here water and decane) and μ and μ g are their respective viscosities.In the integral model of Burdine [20], a Heaviside-shaped retention curve yields K /K 0 = S n b and K g /K 0 = (1 − S) n b , where K 0 /(1 − ν) (d 2 /180)[(1 − ν)/ν] 2 and n b 3. From the periodic unit cell of volume V cell , we calculate the Carman-Kozeny equivalent "particle diameter" d = (6νV cell /π ) 1/3 .

3 / 2 n
0 /v c 0 2.49 (see Appendix H), so|v p − v p 0 |/v p 0 < 0.23(1 − σ 2 ),∀ θ c < π/2 and <0.1(1 − σ 2 )for typical θ c > 35 • .Because phase transitions occur near | σ | 1, fluctuations in pore volume due to θ c can therefore be ignored in most cases.Finally, unless necks are unusually long, the position of the triple contact line on their periphery should be relatively insensitive to the contact angle, and therefore the pore surface area a p should not significantly depend on θ c either.Overall, hydrophilic contact angles should mostly affect the magnitude of λ.From Eq. (E1), λ/λ 0 2/(1 + sin θ c ).

FIG. 8 .
FIG. 8. Apparent macroscopic retention curve for two adjacent mesoscopic domains subject to the same capillary pressure but different initial mean degree of saturation, and drawn for conditions of Fig. 3. (Right) The open square represents initial degrees of saturation S a 0 = 0 and S b 0 = 1, and an initial capillary pressure ψ 0 between ψ + and ψ − .(Left) The solid square is the same initial pressure but with S a 0 = 1 and S b 0 = 0.Both apparent retention curves invade the region between the main curves (dashed lines) shown in Fig. 3(a).