Pore saturation model for capillary imbibition and drainage pressures

Leaching and transport of radionuclides from cementitious waste forms and from waste tanks is a concern at the Savannah River Site and other Department of Energy sites. Computer models are used to predict the rate and direction for migration of these through the surrounding soil. These models commonly utilize relative permeability and capillary pressure correlations to calculate migration rates in the vadose (unsaturated) zone between the surface and the water table. The most commonly used capillary pressure models utilize two parameters to relate the pressure to the relative saturation between the wetting (liquid) and nonwetting (gas) phases. The correlation typically takes the form of a power law relation or an exponential equation. A pore saturation model is used to derive the secondary drainage pressure and the bounding imbibition pressure as functions of a characteristic pore pressure and the liquid saturation. The model utilizes singularity analyses of the total energies of the liquid and gas to obtain residual saturations for the two phases. The model successfully correlates a selected set of laboratory imbibition and drainage data for sand. The capillary pressure model utilizes a single fitting parameter, a characteristic pore pressure, which is related to a characteristic pore diameter by the Laplace equation. This pore diameter approximately equals the diameters predicted by two different geometric pore models based on the particle diameter.


Background
A general parametric model for secondary drainage and bounding imbibition pressures for quasi-steady state flow is proposed. The model is derived by equating changes in the potential energy associated with capillary forces to the energy required for flow. The model does not include transient effects and therefore does not address scanning curves between the bounding pressures for imbibition and drainage. The proposed model differs from most previous capillary models in that it is based on a single parameter, a characteristic capillary pore pressure. The pore pressure can be related to the mean particle diameter through the Laplace equation and geometric pore models.
Previously, capillary pressures have been evaluated using either parametric correlations or pore-scale models. Parametric correlations typically express the capillary pressure, normalized with respect to a characteristic pore pressure, as a function of the liquid saturation, normalized with respect to the difference between the residual liquid saturation and full saturation. Most capillary pressure correlations utilize two or more parameters to relate the pressure to the relative saturation. As in the proposed model, the characteristic pore pressure serves as one of the parameters. Additional parameters are added as fitting constants.
The most commonly used parametric correlations are those of Brooks and Corey (1964) and van Genuchten (1980). Other correlations include the Brutsaert (1967), lognormal (Kosugi 1994(Kosugi , 1996, and Gardner-Russo (Russo 1988;Gardner 1958) correlations. The Brooks and Corey correlation is a two-parameter power law correlation. Both the van Genuchten and the Brutsaert correlations are comprised of asymptotic power law relationships; the Brutsaert correlation has two parameters and the van Genuchten has three parameters. The lognormal correlation takes a logarithmic form when solved for the normalized saturation as a function of the pressure and has two parameters. The Gardner-Russo correlation is a single-parameter correlation that is expressed as a mixed exponential function when solved for the saturation. Chen et al. (1999) present a comprehensive comparison of these parametric correlations.
Because these parametric correlations, as originally formulated, assume that the liquid can completely saturate the pores, in a strict sense they apply only to primary drainage from a completely saturated porous material. Modifications are required to address the hysteresis between drainage and imbibition and between primary drainage and secondary drainage from a state of residual gas saturation. Along these lines, Parker and Lenhard (1987) developed a model that includes effects of gas entrapment during imbibition. Their model defines limiting maximum capillary heads for drainage and minimum capillary heads for imbibition. They state that the actual pressure must lie at some value between these limits that depends on the flow history. Another way to incorporate hysteresis, applied by White and Oostrom (1996), is simply to define separate limiting capillary pressures for imbibition and drainage.
Pore-scale models encompass pore network models and lattice Boltzmann methods. Pore network models, introduced by Fatt (1956), use Monte Carlo methods to calculate pressures and flows in a network of capillary pores with a specified coordination number, loosely defined as the number of neighboring pores connected to each pore. In lattice Boltzmann methods, the fluid flow through the porous network is modeled as the flow of a stream of particles in a lattice of pores. In an early implementation of the lattice Boltzmann method, Gray (1990, 1993) included an interfacial area to enable the modeling of hysteresis (Porter et al. 2009). This concept was extended to a porous network model by Reeves and Celia (1996).
A complicating factor in modeling the capillary pressure is the existence of fractures through which flow preferentially occurs. To model fractured media, Klavetter and Peters (1986) and Nitao (1988) used dual porosity functions that assign different capillary pressure-saturation relations to the fractures and to the low permeability solid matrix. Their models equated the fracture and solid matrix pressures. Di Donato and Blunt (2004) incorporated dual porosities into a streamline flow model.
In this analysis, a probabilistic pore pressure model and a mechanical energy balance replace the semi-empirical approach of previous parametric correlations. The energy balances are solved for singularities in the driving forces for flow of the wetting and nonwetting phases, to obtain the residual wetting and nonwetting saturations, respectively. The pore pressure models are then used to calculate the capillary pressure as a function of saturation. An adjusted saturation is defined to account for preferential flow through fractures or fissures. Results are compared to selected capillary pressure measurements.

Derivation of energy balance for capillary pressure model
The singularity analyses for the residual saturations are based on comparisons of the potential energy associated with the difference between the pressures of the wetting and nonwetting phases in the porous material (hereafter referred to as the liquid and gas phases, respectively) and the work required for flow of each phase. The analysis is loosely analogous to the Gibbs' adsorption theorem, which relates work done by surface forces to the chemical potential of a system. The model replaces the chemical potential with a mechanical potential based on the stored interfacial energy and substitutes a flow work term for the work needed to extend an interface. The model is applied to an adiabatic control volume within a larger volume of the porous material, from which the only energy transfer occurs by flow of gas and liquid in or out. The model evaluates changes of the potential energy within the control volume and changes in the amount of flow work performed by the control volume on its surroundings.
The energy balance for the control volume equates changes in the enthalpy, ΔH, to changes in the potential energy associated with saturation of the porous material, ΔG s , and the work required for flow, ΔW The energy balance is formulated in terms of the enthalpy because it is an open system in which liquid and gas may cross the volume boundaries.
Because the model assumes that the soil is adiabatic, the enthalpy term is zero, and The mechanical potential and flow work terms can be expressed in terms of pressures much as the chemical potential in the Gibbs' adsorption equation. For a control volume initially comprised of one mole of gas, these expressions are and (1) where R g is the ideal gas constant, T is the temperature, P pot is the potential pressure associated with partial saturation of capillary pores, and P flow is a characteristic pressure arising due to capillary flow.

Derivation of pore pressure model
The model employs different flow pressures for imbibition and drainage. The flow pressure for imbibition is based on liquid phase flow, while the flow pressure for drainage is based on gas phase flow. To derive the pressures for the potential energy and the flow work, a model for the pore pressure as a function of saturation must be developed. The pore pressure model assumes that each pore is filled with either a liquid wetting or a gaseous nonwetting phase. The capillary pore pressure is distributed by gas-liquid interfaces between the pores. The distribution of gas-and liquid-filled pores and the connections between pores are assumed to be random. The difference between the average pore pressure and the liquid pressure is equal to the product of the capillary pore pressure, P c , and the probability that a liquid-filled pore is in contact with a gas-filled pore and not with liquid in another pore. This probability, in turn, is the gas saturation, s. Thus, the difference between the average pore pressure, P a , and the liquid pressure, P l , is given by Likewise, the difference between the gas pressure, P g , and the average pore pressure is equal to the product of the capillary pressure and the liquid saturation, or, The potential energy is calculated from the difference between P a and the reference pressure for flow through the soil, P o . For an influx of liquid into the porous material (imbibition), the reference pressure is the liquid phase pressure, and For flow of gas (drainage), the reference pressure is the gas phase pressure, so that These pressure differences apply only to the fractions of the pore volume occupied by liquid and gas, respectively. Consequently, for either liquid or gas flow, the potential energy is based on a pressure differential given by The pressure differential is defined to be positive for both imbibition and drainage. Substitution of this pore pressure in the equation for the potential energy yields

Analysis of residual saturations
At the residual wetting phase saturation, s wr , the resistance to flow of the liquid becomes infinite, and, in an overall energy balance, the resistance to flow of the gas can be ignored. It may be stated, then, that for a given change in the saturation, s, the change in the potential energy counteracts the change in the work needed to cause the liquid to flow. In other words, at s = s wr , The liquid work function is derived from the Darcy equation for flow of two immiscible phases. The Darcy equation for the liquid phase velocity, v w , is: or where ɛ is the porosity, k is Darcy permeability coefficient, k w is the relative liquid permeability, μ w is the liquid viscosity, and z is the direction of flow.
From this expression, the pressure acting on a given volume of liquid is seen to be P c ln (s). The change in the corresponding work function is given by As with the potential energy, this expression gives the change in the work function for small changes in the saturation.
Substitution of this work term and the change in the potential energy in the energy balance yields or This equation is satisfied for s wr = 0.236. A similar line of reasoning is used to calculate the residual nonwetting phase saturation, s nr . Here, the resistance to the gas flow is infinite, and the resistance to flow of the liquid can be ignored. Thus, The Darcy equation for the gas phase velocity, v n , is: or where k n is the relative gas permeability and μ n is the gas viscosity. From this expression, the pressure acting on a given volume of gas is seen to be − P c ln (1 − s). The capillary pressure is subtracted from this pressure. The capillary pressure is supplied to the gas to convert it from a stationary state in which it is in contact with a pore wall to a mobile state in which it is surrounded by liquid. The total pressure supplied to the gas phase is, therefore, − P c (1 + ln (1 − s)), and the corresponding change in the work function is: Substitution of this work term and the change in the potential energy in the energy balance gives or This equality is satisfied for s nr = 0.884. Measured residual saturations for a uniform sand (Demond and Roberts 1991;Morel-Seytoux and Khanji 1975;Morel-Seytoux et al. 1973) agree almost exactly with the derived values; the measured wetting saturation, estimated by graphical interpolation, was 0.230, and the measured nonwetting saturation was 0.884.

Calculation of capillary pressures for imbibition and drainage
The following section describes a capillary pressure model for imbibition and drainage, based on the residual nonwetting saturation. For imbibition, the liquid supplies the motive force for displacement of the gas.
Consequently, it may be argued that the effective capillary pressure within the porous material is the integral of the pressure gradient for the liquid. During imbibition, the force to displace the gas is applied over the entire surface of the porous material, but the porous material is aerated so that the air within the capillary pores remains at atmospheric pressure. A hydraulic advantage accompanies the force applied to the liquid in the pores, then, in inverse proportion to the liquid saturation. It follows that the measured liquid pressure is the integral of the liquid phase pressure gradient, divided by the saturation. For the measured capillary pressure for imbibition, P c,i , this gives: where P c,i,0 represents the minimum pressure required for imbibition of liquid into a saturated porous material containing non-displaceable gas, i.e., the entry head. During drainage, both liquid and gas displace the liquid that leaves the porous material. Hence it may be argued that the effective capillary pressure is the volume average of the integrals of the pressure gradients for the gas and the liquid. Again, because the air remains at atmospheric pressure during drainage, there is a hydraulic advantage applied to the measured pressure, so the effective pressure is this volume-average pressure integral, divided by the saturation.
At the maximum residual nonwetting phase saturation, the measured suction pressure drops by a step change from its value for drainage to its value for imbibition. This pressure change can be attributed to the change from a continuous gas phase at atmospheric pressure to a continuous liquid phase for which, at static equilibrium, the outside pressure equals the liquid phase pressure. Since the liquid phase pressure is less than the gas phase pressure by a pressure difference equal to the characteristic capillary pore pressure, the suction pressure for imbibition must be less than the pressure for drainage by just this pressure, divided by the saturation to account for the effective hydraulic advantage. These considerations, with adjustments to account for differences between the integration constants for the liquid and volume-average pressure gradients, give, for the capillary pressure for drainage, P c,d , (1 − s nr )(1 + ln (1 − s nr )) − s nr ln (s nr ) s nr The entry head most likely is a function of the capillary pore pressure. The most plausible explanation that fits the measured data is that an excess capillary force is required to displace air from an array of pores located on the surface in all three directions (one perpendicular to the surface and two in transverse directions). The required force is the three-dimensional vector sum of the forces required for unidirectional displacement, one component of which corresponds to the pressure. According to this interpretation, the entry head is the pore pressure multiplied by the square root of three: The entry head is applied at the outer surface of the porous material, so this pressure difference is not normalized with respect to the liquid saturation.
A combination of the capillary pressure equations for imbibition and drainage yields a critical saturation, s cr , where these two pressures coincide: Substitution of s nr in this expression gives s cr = 0.301. Below this saturation, the suction pressure that develops during drainage is less than the suction pressure that accompanies imbibition. The logical conclusion is that the porous material will not drain to a saturation below this critical value, unless the source of liquid is removed from contact with the porous material or there is an external gas flow.

Comparison of model predictions to capillary pressure data
Data from the UNsaturated SOil hydraulic DAtabase (UNSODA) maintained by the George E. Brown, Jr., Salinity Laboratory of the U. S. Department of Agriculture (Nemes et al. 1999(Nemes et al. , 2001 was selected for comparison with the model. The comparison is limited to data sets with paired imbibition and drainage pressures that exhibited a significant entry head P c,i,0 and that included particle size measurements. There were three such data sources, namely Shen and Jaynes (1988), Stauffer and Dracos (1986), and Poulovassilis (1970). An additional data source not included in UNSODA, that of Smiles et al. (1971) and Vachaud and Thony (1971), also is analyzed. All these data sources contained measurements made in sand. The fitting procedure entailed a calculation (1 − s cr )(1 + ln (s cr ) + ln (1 − s cr )) s cr of the capillary pore pressure P c , computed as the average of the measured capillary pressures divided by the capillary pressures predicted by the model at the relative saturations for each measurement. To eliminate measurements near saturation with either liquid or gas, the calculation was limited to relative liquid saturations between 0.4 and 0.75. Figure 1 compares these data with model predictions. As this figure shows, the Shen and Jaynes data and the Stauffer and Dracos data have minimum measured saturations about half the predicted critical minimum saturation s cr , while the Poulovassilis and Smiles et al. data have minimum measured saturations approximately equal to s cr .
The model fits the data closely for the Poulovassilis and Smiles et al. data. The most logical explanation for the apparent discrepancy between the minimum measured saturations for the Shen and Jaynes and the Stauffer and Dracos data and the model prediction is that the porous material develops fissures that fill with gas during drainage. According to this interpretation, the fissures are larger than the pores by a sufficient magnitude that they do not offer any flow resistance but rather simply conduct fluid among the pores.
The minimum saturation for a fissured material, s min , can be derived by applying a jump condition for the development of fissures. Following the logic used in the The minimum saturation is the saturation that will minimize the free energy associated with this jump, ΔG jump . This free energy, evaluated after the appearance of the fissures, is This condition is satisfied by A similar analysis, combined with the observation that the maximum saturation for a homogeneous porous material is the residual nonwetting saturation, gives for the maximum saturation for a fissured porous material, s max , The capillary pressure model describes the saturation for just the pores of a fissured material. Because the fissure volume merely conducts pressure among the pores, the overall saturation for the total volume of porous material is a linear function of the pore saturation. If the maximum saturation is defined to be equal to the residual nonwetting saturation to conform with the model, then a linear interpolation between this saturation and the minimum saturation yields the following expression for the overall saturation of a fissured material, s*, as a function of s, s nr , and s cr : Figure 2 compares the capillary pressure model with the Shen and Jaynes data and the Stauffer and Dracos data defined by s*. The model fits the adjusted imbibition data closely, but the Stauffer and Dracos drainage data deviate significantly from the model. The difference between the measured and predicted drainage pressures remains approximately constant as the saturation increases. This suggests the presence of a metastable equilibrium that originates from the upstream (dry) side of the capillary pressure gradient.
The most plausible explanation for the overshoot and apparent metastability phenomena is that at the minimum liquid saturation, the interstitial gas in the fissures begins to exert pressure on the residual liquid in the pores. As explained in the capillary pressure model derivation, this pressure would be the product of the capillary pore pressure P c and the fraction of the volume (27) �P jump = P c ((1 − s min ) − (1 − s cr )) = P c (s cr − s min ) (28) �G jump = −R g T ln (P c s min (s cr − s min )) (29) s min = 0.5s cr (30) s max = 0.5 + 0.5s nr (31) s * = 0.5s cr s nr − s cr s + s nr s s nr − 0.5s cr occupied by the interstitial gas, 0.5(1 − 0.5s cr ), divided by the residual pore saturation, 0.5s cr . Thus, the overshoot pressure, ΔP c,ov , is given by According to this interpretation, the apparent metastability results when the overshoot in the capillary pressure at minimum saturation propagates downstream in the nonwetting (gas) phase. The fraction of the overshoot pressure that propagates would be equal to the fraction of the pore volume occupied by gas, 1 − s cr . If this overshoot pressure propagates, then the metastable drainage pressure, P c,d,ms , would be given by As shown by Fig. 2, the predicted overshoot pressure approximately equals maximum measured pressures, and the Stauffer and Dracos data closely follow the predicted metastable drainage pressure curve.
It appears that the Stauffer and Dracos drainage data are not constrained to the maximum residual nonwetting phase saturation but instead extend to complete saturation. This suggests that the additional capillary (32) P c,ov P c = 1 − 0.5s cr s cr Shen and Jaynes, imbibition Shen Fig. 2 Variation of imbibition and drainage pressures with relative saturation, corrected for fissure effect pressure at the metastable equilibrium is sufficient to overcome the capillary pore pressure.

Evaluation of characteristic pore pressures
The final step in development of the capillary pressure model is to relate the characteristic pore pressure P c derived from the data fit to the particle size of the porous material. One may recall that, for each data set that was regressed, the pressure multiplier is 1 P c . This pore pressure is related to the capillary rise pressure or capillary head P c,i,0 by Eq. 25. The Laplace equation is used to relate the capillary head to the effective pore diameter, d pore , the interfacial tension, σ, and the liquid wetting angle, θ: Two geometric models have been developed to relate the effective pore diameter to the particle diameter and the porosity. From an analysis of the capillary rise in a powder, based on the specific surface area of the powder, White (1982) derived the following relationship between the particle diameter, d part , and the pore diameter: The Revil-Glover-Pezard-Zamora (RGPZ) model (Glover et al. 2006) was developed to correlate the permeability as a function of particle size based on the limiting flow rate through the necks connecting adjacent pores. According to the RGPZ model, for spherical particles, the effective pore diameter is related to the particle diameter by: The RGPZ model was developed to describe flows in low-porosity oil sands and therefore contains only the leading term for the porosity. The laminar term in the Ergun equation for flow through packed beds (Ergun 1952) indicates that, at low flow rates, the particle diameter scales as ε 1.5 1−ε with respect to pressure. This suggests that a modified RGPZ model for higher porosity materials should take the form: Both the White and RGPZ models were developed for a material comprised of uniform size spheres. For this study, it is assumed that the mean particle diameter calculated from sieve tray data is representative of the particle diameter in the models. The particle diameter data are binned by tray mesh size. Accordingly, mean particle diameters were calculated using an interval-censored method, based on an assumed log normal size distribution. The interval-censored averaging was performed using JMP ® statistical software 1 (see So et al. 2010). For each data set, the void fraction was either given, calculated from the bulk density of the sand, or estimated by fitting the residual nonwetting saturation to the capillary pressure data.
The Laplace equation for the pore diameter was evaluated using the properties of water in contact with air. The surface tension of water at room temperature is 0.0722 N/m (Lide 1994). The measured wetting angle is highly variable and depends on the measurement method. The calculated pore diameter is based on a dynamic wetting angle of 45°, measured for air-dried sand (Weisbrod et al. 2009). Table 1 compares pore diameters calculated using the White model, the RGPZ model, and the modified RGPZ model to the characteristic pore diameter for the capillary pressure model presented in this study. As may be noted, pore diameters calculated from the White and modified RGBZ models agree relatively closely with the effective pore diameter based on the capillary pressures, whereas the original RGBZ model underestimates the effective pore diameter.

Statistical comparison of pore diameter models
It is informative to compare the three geometric pore diameter models in terms of how closely they fit the calculated results from the capillary pressure model. Figure 3 plots the results from Table 1 in terms of the ratio of the pore diameter to the average particle diameter. The datum points in this figure are calculated by dividing the average pore diameter for each of the five analyzed tests by the pore diameter calculated from the model and the Laplace equation.
A statistical comparison of the three geometric models indicates that the White and modified RGPZ models fit the pore diameters calculated from the capillary pressure measurements significantly better than the original RGPZ model. To compare the models, the Bayesian Information Criterion (BIC) (Schwarz 1978) was applied to the prediction of the pore diameter to particle diameter ratio given the porosity and to the prediction of the porosity given the diameter ratio. Equal weights were given to each of these comparisons. The modified RGPZ provides the best fit; the BIC scores for the original RGPZ model and the White model relative to the modified RGPZ model were 8.9 and 0.9, respectively. According to the criteria proposed by Kass and Raftery (1995), the BIC score for the original RGPZ model provides a strong indication that this is not the correct model. The BIC scores for the modified RGPZ and White models do not differ sufficiently to distinguish which is the superior model.

Conclusions
A model has been developed that uses a characteristic pore pressure to predict residual wetting and nonwetting saturations and capillary imbibition and drainage pressures. This model has a closed form solution that does not contain any empirical or dimensional constants. The model is based on an energy balance that equates changes in potential surface energy to changes in pressure-volume work. A simple probabilistic distribution expresses capillary forces as a function of the relative saturation and the characteristic pore pressure. The model implicitly assumes that the porous material is homogeneous and isotropic.
The model predicts a residual wetting saturation, s wr , of 0.236 and a residual nonwetting saturation, s nr , of 0.884. These predicted residual saturations are in almost exact agreement with measurements. In addition, the model predicts that there is a critical saturation, s cr , of 0.301, below which the imbibition pressure exceeds the drainage pressure. A porous material will not drain to any lower saturation, provided that the porous material remains in contact with liquid and there is no external gas flow. The capillary pressure model successfully correlates limiting imbibition and drainage pressures for selected laboratory tests that used uniformly packed, sieved sand, although for certain tests an adjusted saturation is required to account for fissures. For the data that was correlated, the characteristic pore diameter from the capillary pressure model approximately equals pore diameters calculated from the average particle diameter using two different geometric models. This agreement suggests that, for uniform packed sands with a relatively narrow particle size distribution, a capillary pressure curve can be estimated from a measured particle size distribution. The applicability of the model to soils and rocks with widely varying pore size distributions has not been demonstrated.

Authors' contributions
The author read and approved the final manuscript.

Authors' information
James E. Laurinat PE is a Principal Engineer at the Savannah River National Laboratory with over thirty years combined experience in computer modeling, actinide chemical processing, including ion exchange and solvent extraction, nuclear reactor thermal hydraulics, and two-phase flow. Dr. Laurinat holds BS, MS, and PhD degrees in chemical engineering. He is author or coauthor of nine journal publications.  Fig. 3 Comparison of pore diameters from fit of capillary pressure data to pore diameters calculated from mean particle diameters