Skip to main content

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.


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, 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, Hassanizadeh and 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, ΔGs, and the work required for flow, ΔW

$$\Delta H = \Delta G_{s} + \Delta 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

$$\Delta G_{s} = - \Delta W$$

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

$$\Delta G_{s} = - R_{g} T\Delta \rm{{ln}} \,\left( {P_{pot} } \right)$$


$$\Delta W = - R_{g} T\Delta \left( { - {\rm{ln}}\, \left( {P_{flow} } \right)} \right)$$

where Rg is the ideal gas constant, T is the temperature, Ppot is the potential pressure associated with partial saturation of capillary pores, and Pflow 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, Pc, 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, Pa, and the liquid pressure, Pl, is given by

$$P_{a} - P_{l} = P_{c} \left( {1 - s} \right)$$

Likewise, the difference between the gas pressure, Pg, and the average pore pressure is equal to the product of the capillary pressure and the liquid saturation, or,

$$P_{g} - P_{a} = P_{c} s$$

The potential energy is calculated from the difference between Pa and the reference pressure for flow through the soil, Po. For an influx of liquid into the porous material (imbibition), the reference pressure is the liquid phase pressure, and

$$P_{a} - P_{o} = P_{c} \left( {1 - s} \right)$$

For flow of gas (drainage), the reference pressure is the gas phase pressure, so that

$$P_{a} - P_{o} = - P_{c} s$$

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

$$P_{pot} = P_{c} s\left( {1 - s} \right)$$

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

$$\Delta G_{s} = - R_{g} T\Delta \it{ln} \left( {P_{c} s\left( {{\text 1} - s} \right)} \right)$$

Analysis of residual saturations

At the residual wetting phase saturation, swr, 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 = swr,

$$\frac{{{\rm{d}}G_{s} }}{{\rm{d}}s} = - \frac{{{\rm{d}}W_{l} }}{{\rm{d}}s}$$

The liquid work function is derived from the Darcy equation for flow of two immiscible phases. The Darcy equation for the liquid phase velocity, vw, is:

$$\varepsilon sv_{w} = - \frac{{kk_{w} }}{{\mu_{w} }}\frac{{{\rm{d}}P_{l} }}{{\rm{d}}z} = - \frac{{kk_{w} }}{{\mu_{w} }}P_{c} \frac{{\rm{d}}s}{{\rm{d}}z}$$


$$v_{w} = - \frac{{kk_{w} }}{{\mu_{w} \varepsilon }}P_{c} \frac{{\rm{d}}{\rm{ln}} \,\left( s \right)}{{\rm{d}}z}$$

where ɛ is the porosity, k is Darcy permeability coefficient, kw 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 Pc ln (s). The change in the corresponding work function is given by

$$\Delta W_{l} = - R_{g} T\Delta {\rm{ln}}\, \left( { - P_{c} {\rm{ln}}\, \left( s \right)} \right)$$

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

$$- R{}_{g}T\frac{{{\rm{d}}{\rm ln}\, \left( {P_{c} s\left( {{\text 1} - s} \right)} \right)}}{{\rm{d}}s} = - R_{g} T\frac{{{\rm{d}}{\rm{ln}} \,\left( { - P_{c} {\rm{ln}} \,\left( s \right)} \right)}}{{\rm{d}}s}$$


$${\rm{ln}}\, \left( s \right) = \frac{s - {\text 1}}{{\text 1} - 2s}$$

This equation is satisfied for swr = 0.236.

A similar line of reasoning is used to calculate the residual nonwetting phase saturation, snr. Here, the resistance to the gas flow is infinite, and the resistance to flow of the liquid can be ignored. Thus,

$$\frac{{{\rm{d}}G_{s} }}{{\rm{d}}s} = - \frac{{{\rm{d}}W_{g} }}{{\rm{d}}s}$$

The Darcy equation for the gas phase velocity, vn, is:

$$\varepsilon \left( {{\text 1} - s} \right)v_{n} = - \frac{{kk_{n} }}{{\mu_{n} }}\frac{{{\rm{d}}P_{g} }}{{\rm{d}}z} = - \frac{{kk_{n} }}{{\mu_{n} }}P_{c} \frac{{\rm{d}}s}{{\rm{d}}z}$$


$$v_{n} = - \frac{{kk_{n} }}{{\mu_{n} \varepsilon }}\left( { - P_{c} \frac{{{\rm{d}}{\rm{ln}} \,\left( {{\text 1} - s} \right)}}{{\rm{d}}z}} \right)$$

where kn 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 − Pc 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, − Pc(1 + ln (1 − s)), and the corresponding change in the work function is:

$$\Delta W_{g} = - R_{g} T\Delta {\rm{ln}} \,\left( { - P_{c} \left( {{\text 1} + {\rm{ln}} \,\left( {{\text 1} - s} \right)} \right)} \right)$$

Substitution of this work term and the change in the potential energy in the energy balance gives



$${\rm{ln}}\,\left( {{\text 1} - s} \right) = \frac{3s - {\text 1}}{{\text 1} - 2s}$$

This equality is satisfied for snr = 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, Pc,i, this gives:

$$\frac{{P_{c,i} }}{{P_{c} }} = \frac{{P_{c,i,0} }}{{P_{c} }} - \frac{{\rm{ln}}\, \left( s \right)}{s}$$

where Pc,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, Pc,d,

$$\begin{aligned} \frac{{P_{c,d} }}{{P_{c} }} & = \frac{{P_{c,i,0} }}{{P_{c} }} + \frac{1}{{s_{nr} }} - \frac{{{\rm{ln}}\, \left( {s_{nr} } \right)}}{{s_{nr} }} + \frac{{\left( {{\text 1} - s} \right)\left( {{\text 1} + {\rm{ln}}\, \left( {{\text 1} - s} \right)} \right) - s{\rm{ln}}\, \left( s \right)}}{s} \\ & \quad - \frac{{\left( {{\text 1} - s_{nr} } \right)\left( {1 + {\rm{ln}}\, \left( {{\text 1} - s_{nr} } \right)} \right) - s_{nr} {\rm{ln}}\, \left( {s_{nr} } \right)}}{{s_{nr} }} \\ \end{aligned}$$

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:

$$\frac{{P_{c,i,0} }}{{P_{c} }} = \sqrt 3$$

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, scr, where these two pressures coincide:

$$\frac{{1 - \left( {1 - s_{nr} } \right)\left( {1 + {\rm{ln}}\, \left( {s_{nr} } \right) + {\rm{ln}}\, \left( {{\text 1} - s_{nr} } \right)} \right)}}{{s_{nr} }} = - \frac{{\left( {1 - s_{cr} } \right)\left( {1 + {\rm{ln}}\, \left( {s_{cr} } \right) + {\rm{ln}}\, \left( {{\text 1} - s_{cr} } \right)} \right)}}{{s_{cr} }}$$

Substitution of snr in this expression gives scr = 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, 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 Pc,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 of the capillary pore pressure Pc, 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 scr, while the Poulovassilis and Smiles et al. data have minimum measured saturations approximately equal to scr.

Fig. 1

Variation of imbibition and drainage pressures with relative saturation

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, smin, can be derived by applying a jump condition for the development of fissures. Following the logic used in the analysis of residual wetting saturations, the pressure applied to the liquid phase across such a jump, ΔPjump, is

$$\Delta P_{jump} = P_{c} \left( {\left( {1 - s_{\text{min} } } \right) - \left( {1 - s_{cr} } \right)} \right) = P_{c} \left( {s_{cr} - s_{\text{min} } } \right)$$

The minimum saturation is the saturation that will minimize the free energy associated with this jump, ΔGjump. This free energy, evaluated after the appearance of the fissures, is

$$\Delta G_{jump} = - R_{g} T\ln \left( {P_{c} s_{\text{min} } \left( {s_{cr} - s_{\text{min} } } \right)} \right)$$

This condition is satisfied by

$$s_{\text{min} } = 0.5s_{cr}$$

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, smax,

$$s_{\text{max} } = 0.5 + 0.5s_{nr}$$

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, snr, and scr:

$$s* = \frac{{0.5s_{cr} s_{nr} - s_{cr} s + s_{nr} s}}{{s_{nr} - 0.5s_{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.

Fig. 2

Variation of imbibition and drainage pressures with relative saturation, corrected for fissure effect

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 Pc and the fraction of the volume occupied by the interstitial gas, 0.5(1 − 0.5scr), divided by the residual pore saturation, 0.5scr. Thus, the overshoot pressure, ΔPc,ov, is given by

$$\frac{{\Delta P_{c,ov} }}{{P_{c} }} = \frac{{1 - 0.5s_{cr} }}{{s_{cr} }}$$

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 − scr. If this overshoot pressure propagates, then the metastable drainage pressure, Pc,d,ms, would be given by

$$\frac{{P_{c,d,ms} }}{{P_{c} }} = \frac{{P_{c,d} }}{{P_{c} }} + \frac{{\left( {1 - s_{cr} } \right)\left( {1 - 0.5s_{cr} } \right)}}{{s_{cr} }}$$

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 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 Pc 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 \(\frac{1}{{P_{c} }}\). This pore pressure is related to the capillary rise pressure or capillary head Pc,i,0 by Eq. 25. The Laplace equation is used to relate the capillary head to the effective pore diameter, dpore, the interfacial tension, σ, and the liquid wetting angle, θ:

$$d_{pore} = \frac{4\sigma \cos \left( \theta \right)}{{P_{c,i,0} }}$$

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, dpart, and the pore diameter:

$$d_{pore} = \frac{2\varepsilon }{{3\left( {1 - \varepsilon } \right)}}d_{part}$$

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:

$$d_{pore} = \frac{{2\sqrt 3 \varepsilon^{1.5} }}{3}d_{part}$$

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 \(\frac{{\varepsilon^{1.5} }}{1 - \varepsilon }\) with respect to pressure. This suggests that a modified RGPZ model for higher porosity materials should take the form:

$$d_{pore} = \frac{{2\sqrt 3 \varepsilon^{1.5} }}{{3\left( {1 - \varepsilon } \right)}}d_{part}$$

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 Footnote 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.

Table 1 Comparison of pore diameters based on particle diameters with pore diameters from capillary pressure measurements

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.

Fig. 3

Comparison of pore diameters from fit of capillary pressure data to pore diameters calculated from mean particle diameters

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.


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, swr, of 0.236 and a residual nonwetting saturation, snr, 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, scr, 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.


  1. 1.

    JMP is a registered trademark of The SAS Institute of Cary, North Carolina, USA.



Bayes information criterion


Revil-Glover-Pezard-Zamora (pore diameter model)


UNsaturated SOil hydraulic DAtabase


  1. Brooks RH, Corey AT (1964) Hydraulic properties of porous media. Hydrol Pap 3. Civ Eng Dept, Colorado State University, Fort Collins

    Google Scholar 

  2. Brutsaert W (1967) Some methods of calculating unsaturated permeability. T ASEA 10(3):400–404

    Article  Google Scholar 

  3. Chen J, Hopmans JW, Grismer ME (1999) Parameter estimation of two-fluid capillary pressure—saturation and permeability functions. Adv Water Resour 22(5):479–493

    Article  CAS  Google Scholar 

  4. Demond AH, Roberts PV (1991) Effect of interfacial forces on two-phase capillary pressure-saturation relationships. Water Resour Res 27(3):423–437

    Article  Google Scholar 

  5. Di Donato G, Blunt MJ (2004) Streamline-based dual-porosity simulation of reactive transport and flow in fractured reservoirs. Water Resour Res 40(4):W04203

    Article  Google Scholar 

  6. Ergun S (1952) Fluid flow through packed columns. Chem Eng Prog 48(2):89–94

    CAS  Google Scholar 

  7. Fatt I (1956) The network model of porous media: I. capillary pressure characteristics. Petrol Trans AIME 207:144–159

    Google Scholar 

  8. Gardner WR (1958) Some steady state solutions of unsaturated moisture flow equations with application to evaporation from water table. Soil Sci 85(4):228–232

    Article  Google Scholar 

  9. Glover PWJ, Zadjali II, Frew KA (2006) Permeability prediction from MICP and NMR data using an electrokinetic approach. Geophysics 71(4):F49–F60

    Article  Google Scholar 

  10. Hassanizadeh SM, Gray WG (1990) Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv Water Resour 13(4):169–186

    Article  Google Scholar 

  11. Hassanizadeh SM, Gray WG (1993) Thermodynamic basis of capillary pressure in porous media. Water Resour Res 29(10):3389–3404

    Article  Google Scholar 

  12. Kass RE, Raftery AE (1995) Bayes factors. J Am Stat Assoc 90(430):773–795

    Article  Google Scholar 

  13. Klavetter EA, Peters RR (1986) Estimation of hydrological properties of unsaturated fractured rock mass. U.S. DOE Sandia National Laboratory Report SAND84-2642

  14. Kosugi K (1994) Three-parameter lognormal distribution model for soil water retention. Water Resour Res 30(4):891–901

    Article  Google Scholar 


    Article  Google Scholar 

  16. Lide DR (ed) (1994) CRC handbook of chemistry and physics, 75th edn. CRC Press, Boca Raton

    Google Scholar 

  17. Morel-Seytoux HJ, Khanji J (1975) Prediction of imbibition in a horizontal column. Soil Sci Soc Am J 39(4):613–617

    Article  Google Scholar 

  18. Morel-Seytoux HJ, Khanji J, Vachaud G (1973) Prediction errors due to uncertainties in the measurement and extrapolation of the diffusivity function. Civil Engineering Report, Colorado State University, Fort Collins, Colorado 80523, CEP72–73, HJM48

  19. Nemes A, Schaap MG, Leij FJ (1999) The UNSODA unsaturated soil hydraulic database, version 2.0. U. S. Salinity Laboratory, U. S. Department of Agriculture, Agricultural Research Service, Riverside

    Google Scholar 

  20. Nemes A, Schaap MG, Leij FJ, Wosten JHM (2001) Description of the unsaturated soil hydraulic database UNSODA version 2.0. J Hydrol 251(3–4):151–162

    Article  Google Scholar 

  21. Nitao JJ (1988) Numerical modeling of the thermal and hydrological environment around a nuclear waste package using the equivalent continuum approximation: horizontal emplacement. U.S. DOE Lawrence Livermore National Laboratory Report UCID-2144

  22. Parker JC, Lenhard RJ (1987) A model for hysteric constitutive relations governing multiphase flow: 1. Saturation-pressure relations. Water Resour Res 23(12):2187–2196

    Article  CAS  Google Scholar 

  23. Porter ML, Schaap MG, Wildenschild D (2009) Lattice-Boltzmann simulations of the capillary pressure—saturation—interfacial area relationship for porous media. Adv Water Resour 32(11):1632–1640

    Article  Google Scholar 

  24. Poulovassilis A (1970) Hysteresis of pore water in granular porous bodies. Soil Sci 109(1):5–12

    Article  Google Scholar 

  25. Reeves PC, Celia MA (1996) A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour Res 32(8):2345–2358

    Article  Google Scholar 

  26. Russo D (1988) Determining soil hydraulic properties by parameter estimation: on the selection of a model for the hydraulic properties. Water Resour Res 24(3):453–459

    Article  Google Scholar 

  27. Schwarz GE (1978) Estimating the dimension of a model. Ann Stat 6(2):461–464

    Article  Google Scholar 

  28. Shen R, Jaynes DB (1988) Effect of soil water hysteresis on simulation, infiltration, and distribution. Shuili Xuebao J Hydr Eng (China) 10:11–20

    Google Scholar 

  29. Smiles DE, Vachaud G, Vauclin M (1971) A test of the uniqueness of the soil moisture characteristic during transient, nonhysteric, flow of water in a rigid soil. Soil Sci Soc Am J 35(4):534–539

    Article  Google Scholar 

  30. So Y, Johnston G, Kim SH (2010) In: Proceedings, SAS global forum 2010, Analyzing interval-censored data with SAS® Software. paper 257-2010, Seattle, Washington, April 11–14

  31. Stauffer F, Dracos T (1986) Experimental and numerical study of water and solute infiltration in layered porous media. J Hydrol 84(1–2):9–34

    Article  Google Scholar 

  32. Vachaud G, Thony J-L (1971) Hysteresis during infiltration and redistribution in a soil column at different initial water contents. Water Resour Res 7(1):111–127

    Article  Google Scholar 


    Article  Google Scholar 

  34. Weisbrod N, McGinnis T, Rockhold ML, Niemet MR, Selker JS (2009) Effective Darcy-scale contact angles in porous media imbibing solutions of various surface tensions. Water Resour Res 45(4):W00D39

    Article  CAS  Google Scholar 

  35. White LR (1982) Capillary rise in powders. J Colloid Interf Sci 90(2):536–538

    Article  CAS  Google Scholar 

  36. White MD, Oostrom M (1996) STOMP, subsurface transport over multiple phases: theory guide. US DOE Pacific Northwest National Laboratory Report PNNL-11217, UC-2010

Download references

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.

Competing interests

The author declares no competing interests.

Availability of data and materials

Version 2.0 of the UNSODA database has been made publicly available by the US Department of Agriculture National Agricultural Library in Riverside, California. The website for the database can be accessed by entering the following address into a search engine: (Accessed June 2018).

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.


This work was funded by the US Department of Energy Office of Environmental Management under contract number DE-AC09-08SR22470.

Publisher’s Note


Author information



Corresponding author

Correspondence to James E. Laurinat.

Rights and permissions


Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Laurinat, J.E. Pore saturation model for capillary imbibition and drainage pressures. Environ Syst Res 7, 18 (2018).

Download citation


  • Capillary pressure
  • Imbibition
  • Drainage
  • Hysteresis
  • Model