A Closed-Form Equation for Capillary Pressure in Porous Media for All Wettabilities

A saturation–capillary pressure relationship is proposed that is applicable for all wettabilities, including mixed-wet and oil-wet or hydrophobic media. This formulation is more flexible than existing correlations that only match water-wet data, while also allowing saturation to be written as a closed-form function of capillary pressure: we can determine capillary pressure explicitly from saturation, and vice versa. We propose Pc=A+Btanπ2-πSeCfor0≤Se≤1,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{{\text{c}}} = A + B\tan \left( {\frac{\pi }{2} - \pi S_{e}^{C} } \right)\,{\text{for}}\,0 \le S_{{\text{e}}} \le 1,$$\end{document}where Se\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{{\text{e}}}$$\end{document} is the normalized saturation. A indicates the wettability: A>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A>0$$\end{document} is a water-wet medium, A<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A<0$$\end{document} is hydrophobic while small A suggests mixed wettability. B represents the average curvature and pore-size distribution which can be much lower in mixed-wet compared to water-wet media with the same pore structure if the menisci are approximately minimal surfaces. C is an exponent that controls the inflection point in the capillary pressure and the asymptotic behaviour near end points. We match the model accurately to 29 datasets in the literature for water-wet, mixed-wet and hydrophobic media, including rocks, soils, bead and sand packs and fibrous materials with over four orders of magnitude difference in permeability and porosities from 20% to nearly 90%. We apply Leverett J-function scaling to make the expression for capillary pressure dimensionless and discuss the behaviour of analytical solutions for spontaneous imbibition.


Introduction
The capillary pressure is defined as the pressure difference between two fluid phases, normally the pressure of the less dense phase minus the pressure of the denser phase (which here we will assume is water). In local capillary equilibrium capillary pressure is related to the curvature of the fluid menisci through the Young-Laplace equation. In a porous medium, the capillary pressure is a function of the fluid arrangement in the pore space, the displacement history and wettability. Traditionally, however, the capillary pressure is normally fit to a closed-form function of saturation (conventionally the saturation of the denser phase).
A robust capillary pressure model should account for all wettabilities and, at the same time, provide a convenient analytical expression for capillary pressure as a function of saturation that can be inverted to find a closed-form expression for saturation as a function of capillary pressure. We propose the following expression: where the term tan( 2 − S C e ) can be replaced by cot( S C e ) , if preferred. There is no fundamental basis for the choice of this functional form, but, as we show later, it does match accurately a wide range of experimental data. The parameter A is an indicator of wettability. ( A > 0 represents a water-wet or hydrophilic, medium, A < 0 is hydrophobic, while A ≈ 0 represents mixed-wet media where locally both hydrophilic and hydrophobic regions of the pore space are present.) B is the curvature index, which quantifies the magnitude of the capillary pressure. (We show that it can be smaller in mixed-wet media compared to water-wet cases when the fluid menisci form approximately minimal surfaces.) C is the saturation exponent which controls the inflection point of the capillary pressure where the largest portion of the pore space is displaced. As we show later, all three parameters are influenced by the pore structure and wettability of the porous medium.
In many applications it is useful to be able to invert the expression to write saturation as a function of capillary pressure: where S e is the normalised saturation defined as, where S wi and S r are the irreducible water saturation and residual saturation of the other phase (here considered to be either a gas or oil) corresponding to the end points.
In hydrology, and other fields, capillary pressure is almost universally represented by the model of Van Genuchten (1980); where scales the capillary pressure, while m g and n g are coefficients which are a measure of the pore-size distribution and related to each other as m g = 1 − 1∕n g .
Another well-known correlation in the literature is the Brooks and Corey (1966) model. (1) where similar to the Van Genuchten (1980) P g is an entry or typical capillary pressure. is the so-called pore-size distribution index. Both the Brooks and Corey (1966) and Van Genuchten (1980) models are used principally for primary drainage, with S r = 0 ; imbibition processes can, however, be modelled using a finite residual saturation. The major restriction though is that these models are limited to water-wet systems and positive capillary pressure, while most porous media show different degrees of mixed wettability: this is true of soils and rocks exposed to crude oil or organic material, as well as manufactured porous media, such as gas diffusion layers where water-wet carbon fibres are partially coated with a hydrophobic plastic (Blunt and Lin 2021;Zhao et al. 2021;Chi et al. 2019;Niu et al. 2018;Ali et al. 2020;Eltoum et al. 2021;Yao et al. 2021;Deng et al. 2019;Li et al. 2019;Herring et al. 2021;Armstrong et al. 2021;Mahani et al. 2015;Shojaei et al. 2022;Qu et al. 2022.) There are correlations developed to study systems that are not water-wet such as the Huang et al. (1997) model, which uses different parameters to represent spontaneous imbibition ( P c > 0 ) and forced displacement ( P c < 0). A more general model was proposed by Skjaeveland et al. (2000) in the context of oil (o) and water (w) flow, where a w , c w , a o , and c o are constants matched to the data. This model is a general, flexible formulation that can also be used to match and predict arbitrary hysteresis loops representing repeated cycles of water and oil injection. However, it has two disadvantages over the model proposed here, Eq. (1): (1) As well as the end-point saturations, there are four adjustable parameters; as we show later our model with only three parameters can accurately fit a wide range of experimental data; and (2) for many applications, such as determining saturation-height functions and for pore-scale modelling studies where fluid configurations are determined from a known capillary pressure, it is convenient to have a closed-form expression for saturation as a function of capillary pressure -this is not possible, in general, using Eq. (6). For instance, consider multi-scale pore network models (Bultreys et al. 2015;Mehmani et al. 2013;Ruspini et al. 2021;Wang et al. 2022) where the saturation of micro-porous elements is found as a function of the prevailing capillary pressure. It is convenient to do this with a simple analytical expression, such as Eq. (2), to find saturation from capillary pressure.
In the next section, we will first perform a sensitivity analysis on the parameters of our model to examine what effect they have on capillary pressure and to emphasize their physical interpretation. We then fit Eq. (1) to a series of experimental data in different porous media, from rocks to gas diffusion layers in fuel cells, and for different wettabilities. Finally, we discuss the behaviour of analytical solutions for spontaneous imbibition using our model.

Impact of Parameters on the Model Behaviour
We can rewrite Eq. (1) in dimensionless form using Leverett J-function scaling by dividing each term, P c , A and B, by √ K∕ where is the interfacial tension, K is the permeability and is the porosity (Leverett 1941): Note that we deliberately do not include a contact angle in the scaling relationship: this is misleading in mixed-wet media in particular where instead of a single representative contact angle, there is locally a wide range of values (AlRatrout et al. 2018); wettability is implicitly encapsulated in the three parameters of our model.
Based on Eq. (7) we can see how the dimensionless parameters A , B and C affect the dimensionless capillary pressure. The variation of P c with A , B and C is shown in Figs. 1a, b and c, respectively. According to the profiles in Fig. 1a, it can be concluded that A represents the wettability state: A > 0 indicates water-wet media, A < 0 is hydrophobic, while small values of A represent a mixed-wet state. B is related to the pore size distribution and average curvature of the fluid menisci. As we show later the value of B is typically and C (c) on the dimensionless capillary pressure, P c much less than 1. Finally C controls the point of inflection in the capillary pressure, and the asymptotic behaviour near the end points: for C=1 the inflection is at S e =0.5. Note that inflection point is where the magnitude of the capillary pressure derivative with respect to saturation is lowest; this is quantifies the capillary pressure at which the largest portion of the pore space is filled. This model is designed primarily to study water invasion; however, as we show later, it can also accommodate drainage-type displacements where the water saturation decreases. If the model were used to study water both advancing and receding, hysteresis would mean that the parameters used to match the capillary pressure would be different for the two cases. We would expect, for instance, that A would be lower, representing less water-wet conditions for water invasion than for water receding due to contact angle hysteresis. B would also be lower for water injection, as the capillary pressure overall would be higher for the water receding cycle. Table 1 reports the details of various samples used in capillary pressure measurements including porosity ( ), permeability (K), and interfacial tension ( ) between the two phases from which the dimensionless scaling can be derived. For better comparison we fit the experimental data to the dimensionless form of the capillary pressure, Eq. (7). The fitting parameters A , B , and C are also reported in Table 1. In all cases the exponent C is close to 1. Also, the coefficient of determination, R 2 , is reported to show the goodness-offit. For most of these datasets the reported value for R 2 is higher than 0.80 which indicates that the model reliably fits the data. Figure 2 shows the model fit to ten groups of literature data, representing 29 experiments in total: the first group is measured capillary pressure for two water-wet Bentheimer samples (Lin et al. 2018;Raeesi et al. 2014): the first dataset found capillary pressure from measuring meniscus curvature in pore-scale images, while the second used the traditional porous plate technique. A has a positive value representing water-wet conditions. A single model fit is applied to match both datasets simultaneously. Figure 2b shows the fits for two mixed-wet samples for water flooding by Lin et al. (2019) and Alhammadi et al. (2020): here capillary pressure was also determined from directly measuring the interfacial curvature on pore-space images. In both cases A is small and negative indicating mixed-wet conditions. More remarkable is the value of B . The experiments of Lin et al. (2019) were also performed on Bentheimer sandstone, but the value of B is more than an order of magnitude lower that in the water-wet cases displayed in Fig. 2a. The pore structure is the same, but in the mixed-wet case the fluid menisci form approximately minimal surfaces of zero curvature with, in this case, oil bulging into water in one direction, while water bulges into oil in an orthogonal direction. This is not a single contact angle effect: there is a wide range of local contact angle in the system with values both above and below ∕2 . Note that B has a similar magnitude for the other dataset on a reservoir carbonate ) with a wide range of pore size: it is wettability and not the pore-size distribution that governs the capillary pressure behaviour in this case. For the data of Lin et al. (2019) our model fits through the centre of the points (mean saturation and capillary pressures for each fractional flow) with R 2 =0.91. Figure 2c shows the capillary pressure, also for water flooding, in two oil-wet samples by Scanziani et al. (2020) and Alhosani et al. (2020); again the capillary pressure is found from Table 1 The experimental data used to fit the capillary pressure model, and the dimensionless coefficients used  curvature measurements on pore-space images. A is negative consistent with the analysis of local contact angles (Foroughi et al. 2021). Again we see low values of B which is indicative once more of menisci with curvatures of opposite sign in orthogonal directions, as seen directly on the images . The fourth set of data, Fig. 2d, shows capillary pressure for drainage and water flooding in a carbon paper gas diffusion layer where some of the fibres have been coated in a hydrophobic plastic, polytetrafluoroethylene (PTFE) (Hao and Cheng 2012). The good fit to the data indicates that we can use our model for manufactured, fibrous, porous media, as well as porous rocks. Here the material appears water-wet during drainage (air injection) as A > 0 while the material appears slightly hydrophobic during water injection with A < 0 . In this case, the values of A and B have similar magnitude, and we hypothesize that in this case the fluid menisci do not form approximately minimal surfaces.

Match to Experimental Data
We also fitted our model to two other gas diffusion layer experiments (Nguyen et al. 2008), see Fig. 2e. These materials also exhibited both hydrophobic and hydrophilic properties. The behaviour is similar to that observed for the other gas diffusion layer dataset for water flooding, Fig. 2d.
In Fig. 2f the model is fit to two experiments where ethanol imbibed into a glass bead pack (Moseley and Dhir 1996). The model predicts the capillary pressure for these two samples well with R 2 approximately 0.9. Both reveal water-wet type behaviour as evident from the positive values of A. The next two groups are presented in Figs. 2g and h, which show the results of capillary pressure measurements in a series of carbonate samples taken from a producing oil reservoir (Masalmeh and Jing 2006) Fig. 2g shows primary drainage capillary pressures measured on cleaned samples by mercury injection, while Fig. 2h shows centrifuge water flood capillary pressures on samples after ageing -when the samples have been exposed to crude oil. During drainage A is positive, indicating water-wet conditions, while the samples are oil-wet after contact with crude oil during water flooding ( A is negative). B is a measure of the pore size distribution. For example, sample s85d has a uniform pore size distribution while s89a is a sample with wider range of pore size: note that the parameter B for sample s85d is much lower than s89a consistent with the analysis in the original paper (Masalmeh and Jing 2006). Figure 2i shows drainage and imbibition measurements performed on water-wet sandy soil (Zhuang et al. 2017). Finally, Fig. 2i corresponds to brine invasion into a quartz sand containing super-critical CO 2 at pressures of 8.5 and 12.0 MPa ( 45 • , C ), compared to air-brine measurements at 21 • C and 0.1 MPa (Tokunaga et al. 2013). The dimensionless capillary pressure for the air-brine system is higher than CO 2 -brine, indicating more strongly water-wet conditions in this case.

Behaviour of Solutions for Spontaneous Imbibition
The capillary pressure governs the behaviour of imbibition processes. It is possible to construct semi-analytical solutions for spontaneous imbibition for water-wet and mixedwet media (McWhorter and Sunada 1990;Schmid et al. 2011). In this section we study the asymptotic behaviour of these solutions following the approach in Blunt (2017) Chapter 9.
The situation considered is illustrated in Fig. 3: spontaneous counter-current imbibition in one dimension. At the inlet the capillary pressure is zero and the total velocity (the sum of the Darcy fluxes of the two phases) is also zero. The saturation can be found as a function of = x∕ √ t where x is distance and t is time. We will consider the shape of the profile S w ( ) at the inlet, x = 0 (where P c = 0 and we define S e = S * e ) and at the leading edge of the imbibition front, where S e = 0 . The slope of this profile, dS e ∕d is inversely proportional to the capillary dispersion; the capillary dispersion is in turn proportional to the product of the two-phase mobility (the relative permeability) times the saturation gradient of the capillary pressure, dP c ∕dS e . For all capillary pressure models, including the one proposed here, the advancing edge of the imbibition front has an infinite gradient ( dS e ∕d → ∞ ): this is not a shock, but a divergence in slope at S e = 0 since here the water mobility is zero. Even if the capillary pressure diverges, the water phase mobility, for most reasonable relative permeability models, falls to zero faster, giving a zero capillary dispersion at S e = 0 . Define S e = where ≪ 1 near the end point, and assume that the water relative permeability has a power-law type scaling k rw ∼ a . From Eq. (1) we find P c ∼ B∕ C and the capillary pressure derivative, which drives the imbibition process, scales as 1∕ 1+C . Therefore to have zero capillary dispersion at S e = 0 -governed by the product of the water mobility and the gradient of the capillary pressure -we require a > 1 + C . In our model we have C ≈ 1 in all cases, and so we require a Corey water exponent near the end point of 2 or more to reproduce this behaviour, which is usually the case.
The behaviour at the other end point is more interesting: this is the slope of the saturation profile at the inlet. In the van Genuchten model, the capillary pressure diverges at S e ≡ S * e = 1 : this provides an infinite capillary dispersion at the inlet, and the other phase, air, is assumed to be of infinite mobility for all saturations, giving a zero saturation gradient: the imbibing front moves into the porous medium with a flat profile, dS e ∕d = 0.
In contrast, for mixed-wet media, where S * e ≠ 1 , the capillary pressure and fluid mobility are all finite in our model, giving a finite capillary dispersion and a finite saturation gradient at the inlet. If the medium is strongly water-wet and S * e = 1 then the capillary pressure and its gradient diverge at the inlet. Defining now = 1 − S e the capillary pressure gradient scales as 1∕ 1+C for << 1 as before. Now the water mobility is finite, but the relative permeability of the non-wetting phase falls to zero. Assuming again a power-law scaling k rn ∼ b for the non-wetting phase, n, then the capillary dispersion is zero, giving, as at the leading edge of the profile, an infinite saturation gradient at the inlet for b > 1 + C . On the other hand for 1 + C > b the capillary dispersion is infinite and, as in the van Genuchten model, the saturation gradient is zero at the inlet. Both cases are plausible since, for waterwet media, b is typically in the range 1-2: experimental tests are necessary to determine the behaviour, but the advantage of this new model is that both types of spontaneous imbibition profile are possible.

Conclusions and Future Work
A new empirical capillary pressure model has been proposed. Its advantages over existing formulations are that it can accommodate different wettabilities, and there is a closed-form relationship between saturation and capillary pressure. The model can accurately match a wide range of data, measured on rocks, soils, bead and sand packs, as well as manufactured fibrous materials, for a range of wettability.
A dimensionless form for the capillary pressure was proposed with a physical interpretation of the parameters in the model. The scaling of capillary pressure near the end points was investigated and the consequences for predicted spontaneous imbibition profiles was discussed.
Future work could focus on exploring the capillary pressure behaviour for an even wider range of materials, and investigation of hysteresis loops.