# Influence of wettability on water retention curves in unconsolidated porous media

Gerhard Schäfer, Kevin Hernandez Perez, Panav Hulsurkar, Magda Ibrahim

Youssif, François Lehmann, Mohammad Piri

► To cite this version:

Gerhard Schäfer, Kevin Hernandez Perez, Panav Hulsurkar, Magda Ibrahim Youssif, François Lehmann, et al.. Influence of wettability on water retention curves in unconsolidated porous media. Journal of Contaminant Hydrology, 2025, 269, 10.1016/j.jconhyd.2025.104502. insu-05082255

HAL Id: insu-05082255

https://insu.hal.science/insu-05082255v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Influence of wettability on water retention curves in unconsolidated porous media

Gerhard Schäfer<sup>a,*</sup>, Kevin Hernandez Perez<sup>a</sup>, Panav Hulsurkar<sup>b</sup>, Magda Ibrahim Youssif<sup>b</sup>, François Lehmann<sup>a</sup>, Mohammad Piri<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Hanging column

Aging

Fractionally wet

Contact angles

Leverett J-function

Bartell-Osterhof equation

## ABSTRACT

Controlled laboratory experiments were carried out using the hanging column method. Prior to the experiments, three uniform silica sands, which were originally water-wet, were aged in contact with crude oil until they were moderately oil-wet. Five fractionally wet sands were obtained by mixing the water-wet sands with oil-wet sands containing 25, 50 and 75 vol% oil-wet sands. A total of 11 tests formed the basis for the present study. The measured water retention curves showed that the capillary pressure heads were greatly reduced in sands that were oil-wet or fractionally wet. Changes in the wettability of the sands also affected their irreducible water saturation: The higher the proportion of oil-wet sand in the sand mixtures, the lower the irreducible water saturation. To quantify the characteristics of the measured water retention curves, the Van Genuchten model was used. The two optimized parameters seem to indicate a general trend: The higher the volume fraction of oil-wet sand, the higher is $\alpha$ and the lower is $n$. For the three unaged sands and the aged medium-sized sand, it was found that each of the two branches of the measured retention curves can be suitably scaled to a unique curve if, in addition to the petrophysical parameters (intrinsic permeability, porosity, surface tension, gas-water contact angle), the irreducible residual water saturation and the residual air saturation are taken into account. To quantify the observed deviations of the other two aged sands from the unified Leverett J-function, a theoretical fit function was used to match the experimental data of the three unaged sands. The experimental data sets for P2040ag and P100ag were found to be overestimated overall by the fit function. However, when the petro-physical parameters of the unaged sands were used instead of the actual measured parameters, the individual experimental Leverett J-functions came closer to the uniform J-curve. Based on this, it could be concluded that the apparent differences in pore structure between aged and unaged sands in addition to wettability, expressed by the cosine of the contact angle, may have contributed to a further reduction in the capillary pressure plateaus of the aged sands, which was particularly visible and significant in the P2040 and P100 sand. Using the measured static contact angles for two-phase gas-oil and oil-water systems and the measured interfacial tensions when the porous medium is either water-wet or moderately oil-wet, it was shown that the Bartell-Osterhof equation overestimates the measured gas-water contact angles. Reasonable agreement was achieved when a calibrated roughness factor of the solid surface was considered in the Young's contact angle.

## 1. Introduction

Numerical simulation of multiphase flow has become a crucial method for solving water-resource management and other environmental problems such as CO₂ sequestration in saline aquifers and coupled geomechanical processes in the subsurface (Helmig et al., 2013; Celia et al., 2015; Benisch et al., 2020). While multiphase flow models

for water, non aqueous phase liquids (NAPLs) such as petroleum hydrocarbons, and gas phases were initially adopted in petroleum research for reservoir simulation, they were later successfully exploited for hydrosystem monitoring or remediation. Pressure and saturation of fluids are the key variables in the analysis of multiphase flow in porous media (Schäfer et al., 2020). Furthermore, secondary variables such as relative permeabilities and capillary pressures, which are functions of

<footnote><sup>a</sup> Institut Terre et Environnement de Strasbourg (ITES) UMR 7063, Université de Strasbourg, CNRS, Strasbourg, France</footnote>

<footnote><sup>b</sup> Center of Innovation for Flow through Porous Media (COIFPM), Department of Energy and Petroleum Engineering, University of Wyoming, Laramie, WY, USA</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: schafer@unistra.fr (G. Schäfer).</footnote>
fluid saturations, are required for the simulation of displacement processes in porous media.

Conventional macroscopic equations for two-phase flows in porous media require a constitutive relationship for capillary pressure ($p_c$). In hydrology, these constitutive relationships for a air-water fluid system are commonly referred to as water retention curves. The capillary pressure between two immiscible fluids can be defined as (Morrow, 1976)

$$p_c = p_{nr} - p_r \quad (1)$$

where $p_{nr}$ is the pressure in the non-reference (or non wetting) fluid (air) and $p_r$ is the pressure in the reference (or wetting) fluid (water). Drainage is defined as a decrease in the saturation of the reference phase and imbibition as an increase in the saturation of the reference phase. Under near-surface pressure conditions, air is generally close to atmospheric pressure. Therefore, in unconfined porous aquifers, where water is retained in the unsaturated soils by capillary forces above the water table, the water pressure is lower than the atmospheric pressure, resulting in capillary pressure values between air and water higher than zero.

In a typical laboratory experiment, a functional macroscopic capillary pressure-saturation relationship ($p_c-S$) is measured as the difference between the pressure of the non-wetting fluid at the inflow edge of the porous medium and the pressure of the wetting fluid downstream as function of the saturation of the reference fluid in the porous medium. It is known that macroscopic $p_c-S$ relationships accurately describe physics only under equilibrium conditions and also when both fluids are in hydraulic contact. Based on the parameterization procedure proposed by Parker et al. (Parker et al., 1987), the known water retention curves of a porous sample are often used to obtain capillary pressure-saturation curves for the different fluid pairs present in a multiphase system (e.g. Helmig, 1997; Bohy et al., 2006). As found in the early 1960s, water retention curves may also depend on the direction and on the 'displacement history' of the water-air front, further on referred as hysteresis effects (Poulavissilis, 1962; Mualem, 1984; Parker and Lenhard, 1987; Stauffer and Kinzelbach, 2001). Later studies have demonstrated that hysteresis effects are not the only additional influence on this relationship and they show that there are considerable dynamic effects on the measured capillary pressure curves (Hassanizadeh and Gray, 1988; Bourgeat and Panfilov, 1998; Hassanizadeh et al., 2002; Dahle et al., 2005; Løvoll et al., 2011).

Capillary pressure has a significant effect on the behavior and prediction of flow in porous media. When a porous medium has a clear wetting preference for one of the fluid phases, the wetting fluid will fill the smallest available pores. At any given wetting phase saturation, the capillary pressure as well as the relative permeability of the two fluid phases will be controlled by the radius of the smallest pore invaded by the non-wetting fluid. This is reason why capillary pressure measurements are commonly used to estimate two-phase relative permeabilities (e.g. Van Genuchten, 1980). Capillary pressure is a critical parameter in the study of multiphase flow in porous media, especially in the context of reservoir engineering, hydrogeology, and soil science. Several factors influence capillary pressure-saturation relationships. Main key parameters are the pore size distribution, pore geometry and heterogeneity of the porous medium, grain surface roughness, fluid properties, wettability of the porous medium (e.g. Morrow, 1976; Lenormand et al., 1983; Blunt, 2017). The size and distribution of pores in the porous medium including their shape and connectivity play a significant role in capillary pressure-saturation relationships. Tortuosity, throat size, and pore size distribution and the ratio between the pore and throat sizes are important aspects of pore geometry. Smaller pores tend to have higher capillary pressures, which affect the capillary pressure-saturation curve. The degree of heterogeneity in the porous medium can significantly affect $p_c-S$ relationships. Variations in pore size and wettability at different locations within the medium can lead to complex capillary

pressure curves. The interfacial tension between the wetting and non-wetting phases of the fluids involved impact the capillary pressure through the Young-Laplace equation. Furthermore, in the case of brine and oil systems, the salinity and ionic strength of the brine can influence the capillary pressure-saturation relationship due to changes in interfacial tension. Variations in temperature can also affect the interfacial tension and thus have an effect on the capillary pressure curves (Sinnokrot et al., 1971). Roughness of grain surfaces can influence the wetting properties and contact angles defined as fluid-fluid contact angles on the grain surface within a porous medium (Morrow, 1976).

The wettability of a porous medium affects how fluids interact with the solid surfaces. It is generally defined as the tendency of one fluid to spread on or adhere to a solid surface in the presence of other immiscible fluids (Craig, 1971), and is commonly expressed in terms of the contact angle at the NAPL/water/mineral interface (Cassie and Baxter, 1944). In the oilfield terminology (Schön, 2015), a porous medium is strongly water-wet if the oil-water contact angle ($\theta_{ow}$) is smaller than 70 degrees, has an intermediate wettability if the $\theta_{ow}$ is between 70 and 110 degrees, and is strongly oil-wet if the contact angle is larger than 110 degrees. The intermediate state between water-wet and oil-wet can be caused by a mixed-wet system, in which some surfaces or grains are water-wet, and others are oil-wet, or a neutral-wet system, in which the surfaces are not strongly wet by either water or oil. Wettability influences the fluid distribution in the porous medium because it affects the quantity displaced as well as how the displacement proceeds. The wetting fluid will fill the smallest pores and will be present in the largest pores as films on the solid surface. The existence of such films enhances the continuity of the wetting phase. In oil reservoirs, water is believed to be the original wetting fluid, but oil can become the wetting fluid through the deposition of hydrocarbon residues on pore walls or the presence of charged molecules in oil or water. In the case of hydrocarbon spills in aquifers, there is relatively little time for the pore surfaces to be modified and the hydrocarbons should be non-wetting, especially if the mineral is quartz. The presence of calcite can sometimes lead to non-wetting conditions because the calcite surface is positively charged, whereas the quartz (and other silicate) surface is negatively charged. Although most minerals are naturally water-wet, the adsorption of surface-active constituents in oil may alter their wettability to neutral- or oil-wet (Buckley, 1998).

The traditional method of determining water retention curves in the laboratory is the porous plate method or also called hanging water column. A summary of methods used in soil science to determine water retention curves can be found in (Klute, 1986). The hanging column method is based on a series of equilibria established between water in the porous medium and a body of water at known water suctions. It is a simple and effective technique to measure both drainage and imbibition capillary pressure curves in sandy soils at low suctions of up to around -15 kPa. Initially water saturated soil samples are placed on a semi-permeable porous plate, where they are in hydraulic contact with bulk water. By the use of hanging water columns, a defined negative pressure is applied to the water phase of the samples. The capillary pressure is controlled by creating suction under the plate (Haines, 1930) or adjusting the air phase pressure surrounding the sample in the chamber (e.g., Bruce and Welge, 1947; Raeesi et al., 2014). During drainage, the actual water saturation in the sample is generally determined from the outflow liquid volume. Water retention curves obtained during drainage differ from those measured during imbibition. During drainage, water saturation is larger than during imbibition for the same capillary pressure; that is, the soil water retention curve is not unique and is affected by the saturation history. Hysteresis can be attributed to the ink-bottle effect, entrapped air, and the difference in contact angle between an advancing and a receding water-air front over a solid surface (Kutílek and Nielsen, 1994). Instabilities in interface configurations are a primary cause of hysteresis, which are referred to as Haines' jumps (Haines, 1930). These interfaces can be seen to move suddenly from one stable position to another. The spontaneous movement often occurs in
cascades involving large clusters of individual pores. The time of movement is estimated to be on the order of milliseconds (Berg et al., 2013). No matter how small the pressure steps and how much time is allowed for equilibration after each step, the irreversible jumps in the positions of interfaces cannot be avoided. They are inherent to drainage and imbibition.

Leverett (1941) was the first to publish results showing how rocks with the same lithology but of different porosities and permeabilities and for different fluid couples characterized by interfacial tensions and contact angles of the fluid couples could be described by a single function that normalizes their capillary pressure-saturation curves. The so-called Leverett J-function can then be related to water saturation of a porous medium. More recently, Desouky (2003) applied the J-function to capillary pressure measurements in a shale-sandstone reservoir by using normalized “effective” water saturations, taking into account an irreducible water saturation, to summarize the capillary pressure data into a single curve. Validation for various porous media characterized by different intrinsic permeabilities and porosities yielded excellent agreement between measured and calculated data with average relative errors ranging from 3 to 5.5 %. Dury et al. (1998) showed that the air-water *p*<sub>c</sub>-S curve for a sandy soil can be predicted from a single set of capillary pressure-saturation curves if the corresponding surface tensions are known. How the water retention curve of an unconsolidated porous medium under various wettability conditions could be scaled to its known petrophysical properties including interfacial tension and air-water contact angle within a porous medium remains an open question.

While great efforts have been made in the past to specifically investigate the influence of wettability on capillary pressure-saturation relationships in the oil-water system or oil-gas system in solid rocks (e.g. Leverett, 1941; Morrow, 1976; Anderson, 1987), to our knowledge, only a few laboratory studies have been conducted to date to investigate the influence of wettability on water retention curves in unconsolidated porous media (e.g. Liu et al., 2012). The question to what extent the capillary pressure saturation relations in the air-water system in unconsolidated porous media can be described mathematically by changing the wettability is still topical. Fatt and Klikoff (1959) introduced first the concept of fractional water wettability of rock which is defined as the fraction of the solid surface area that is in contact with water. In their studies they showed that oil-water capillary pressure and relative permeability of unconsolidated sands are functions of fractional wettability. More recently, Geistlinger et al. (2021) used mixed hydrophilic and hydrophobic 1-mm sand and 1-mm glass beads as model systems for fractional wet porous media. Using a µ-CT study, they showed the trapping behavior of air in artificial unconsolidated porous media. While the effect of interfacial tension has been studied in the work of Dury et al. (1998) and Nouri et al. (2014) for sandy soils, the influence of changing air-water contact angle on the water retention has not yet been analyzed in detail, especially when the wettability of the unconsolidated porous medium changes from being water-wet to intermediate or oil-wet. Addressing these scientific gaps will lead to a better understanding of the relationship between wettability and water retention curves, which is essential for managing water resources and optimizing environmental conservation.

The aim of the present study was to overcome this limitation by investigating the influence of wettability on the water retention curves in unconsolidated porous media. We therefore performed controlled laboratory experiments using the hanging column method on three uniform quartz sands with a mean grain size diameter between 0.2 and 0.65 mm. The originally water-wet sands were aged in contact with crude oil to become moderately oil-wet. Prior to the main drainage and main imbibition studies of the air-water *p*<sub>c</sub>-S relationships, we quantified their petrophysical properties such as porosity and intrinsic permeability and used in situ µ-CT to measure the static contact angles between gas-water, gas-oil, and oil-water for each of the six sand samples. Five fractionally wet sand samples were obtained from the water-wet and oil-wet samples to complete the series of experiments. A total of

11 experiments to measure the water retention curve formed the basis for the present study.

Our study aimed to: (1) quantify the characteristics of the measured water retention curves using the Van Genuchten model (Van Genuchten, 1980) to assess the value of irreducible water saturation, residual air saturation, parameters α and n for each of the 11 sand samples; (2) verify the scaling of the measured *p*<sub>c</sub>-S relationships based on a modified Leverett J-function, focusing on the effect of both the measured air-water contact angle and irreducible water saturation; and (3) test the ability of the Bartell-Osterhof equation (Bartell and Osterhof, 1927) to predict the measured gas-water contact angles based on the measured contact angles for two-phase gas-oil and oil-water systems when the porous medium is either water-wet or moderately oil-wet.

## 2. Materials and methods

### 2.1. *Sands and fluids used*

The porous medium consisted of three quartz sands, named in this study P100, P4070, and P2040, with a mean grain diameter ($d_{50}$) of 0.21, 0.34, and 0.65 mm, and a uniformity coefficient (defined by $d_{60}/d_{10}$) of approximately 1.59, 1.59, and 1.45, respectively. Fig. 1 shows the grain size distribution of the three soil samples. According to the soil classification of the United States Department of Agriculture, P100 is a fine sand, P4070 a medium sand and P2040 a coarse sand. The grain size distributions were obtained using the laser diffraction method, wherein particle sizes are reported as a volume-equivalent sphere diameter.

The sands originally water-wet were treated to become oil-wet. Results of the contact angle measurements are presented in section 2.2. To study the influence of wettability of the porous medium on the water retention curves, we performed controlled experiments on 11 sand samples: three water-wet sands (P2040, P4070, P100), three oil-wet sands by aging of the initially water-wet sands (P2040ag, P4070ag, P100ag) and five fractionally wet sands obtained by mixture of the water-wet sands with oil-wet sands containing 25, 50 and 75 % volume fractions of oil-wet sands (P2040ag50, P4070ag25, P4070ag50, P4070ag75, P100ag50). The fractionally wet sands are characterized by a heterogeneous wettability due to the locally varying contact angles.

Prior to the aging of the sands in crude oil, the impact of water films initially adsorbed on the quartz surfaces of our sands has been taken into consideration by establishing an initial water saturation (Alotaibi et al., 2010). The dried sands were mixed with brine solution made of distilled water, 2 wt% CaCl₂, 15 wt% NaI and 0.01 wt% NaN₃ for three days; then dewatering the wet sands by spreading them on filter paper and covering them with filter paper for four hours at room temperature so that only the bound water remained on the grains. Then, the sands were put in a stainless-steel high-pressure-high-temperature (HPHT) chamber filled with crude oil at 1000 psi and 120 °C for three weeks. The oil-impregnated sand was then placed on filter paper for three days in a fume hood to desorb and evaporate the excess oil from the grains (Le et al., 2023; Youssif et al., 2024).

The total porosity and particle density of the three water-wet and oil-wet sands were determined using the Archimedes' liquid displacement method (e.g., Flint and Flint, 2002) (Table 1). It is worth noting that the porosity of P2040ag, P4070ag and P100ag is reduced by 10, 23 and 13 %, respectively, compared to their unaged counterparts. The water-repellent behavior of the aged sands in these static experiments was accompanied by greatly reduced soil particle densities, which undoubtedly indicate an increased volume of the sand grains due to the aging process, which was also shown to be a controlling factor in the permeability measurements.

The hydraulic conductivity of the six sands were measured under steady state flow conditions using Darcy's Law in a vertical Plexiglass column (ID = 3 cm, L = 20 cm). To achieve complete initial water saturation of the sands, the sands were poured into the column in small portions under water and, in the case of the aged sands, stirred manually
<figure><img src="image_4.png" /><figcaption>Fig. 1. Cumulative particle size distribution of the experimental porous media.</figcaption></figure>

<table><caption>**Table 1**

Physical properties of aged and unaged sands used in the study (ag = aged).</caption><thead><tr><th>Porous media</th><th>P2040</th><th>P2040ag</th><th>P4070</th><th>P4070ag</th><th>P100</th><th>P100ag</th></tr></thead><tbody><tr><td>Porosity φ (−)</td><td>0.366</td><td>0.333</td><td>0.398</td><td>0.307</td><td>0.368</td><td>0.319</td></tr><tr><td>Soil particle density ρ<sub>s</sub> (gcm<sup>-3</sup>)</td><td>2.59</td><td>2.33</td><td>2.60</td><td>2.25</td><td>2.55</td><td>2.34</td></tr><tr><td>Hydraulic conductivity K (ms<sup>-1</sup>)</td><td>2.26<br/>×<br/>10<sup>-3</sup></td><td>1.32 ×<br/>10<sup>-3</sup></td><td>4.66<br/>×<br/>10<sup>-4</sup></td><td>3.12 ×<br/>10<sup>-4</sup></td><td>1.52<br/>×<br/>10<sup>-4</sup></td><td>1.21 ×<br/>10<sup>-4</sup></td></tr><tr><td>Permeability k (m<sup>2</sup>)</td><td>2.20<br/>×<br/>10<sup>-10</sup></td><td>1.29 ×<br/>10<sup>-10</sup></td><td>4.54<br/>×<br/>10<sup>-11</sup></td><td>3.04 ×<br/>10<sup>-11</sup></td><td>1.48<br/>×<br/>10<sup>-11</sup></td><td>1.18 ×<br/>10<sup>-11</sup></td></tr></tbody></table>

to break up agglomerated sand particles. A thin layer of water was left over the sand pack and a manual tamper was used to compact and achieve a uniform spatial distribution of the sand particles within each sample. The pressure measurements were performed using two highly precise pressure transmitters PR41X (Keller pressure sensors, Switzerland) placed on the column at a vertical distance of 10 cm and recorded on the PC via a serial interface cable. During the experiments, a constant flow rate of water was applied to the top of the water-saturated sand column using a peristatic pump. The water flow rate was monitored in-line by weighing the cumulative outflow differences within a given time period via a digital Sartorius CPA 34001S balance (Sartorius Weighing Technology GmbH, Germany). Four hydraulic gradients were applied in each flow experiment, varying from 0.15 to 0.75, 0.5 to 3.3 and 1.5 to 5.5 for the sands P2040 and P2040ag, P4070 and P4070ag, and P100 and P100ag respectively. To quantify the hydraulic conductivity of the individual sands, an automatic linear fit of the relationship between Darcy velocity and hydraulic gradient was performed. The regression line determined using the least mean squares method yielded the slope of the straight line corresponding to the hydraulic conductivity with coefficients of determination (R squared) close to one. As expected, the fine sand P100 has a finer pore structure than P2040 and P4070, which contributes to a lower permeability compared to the coarse and medium sand (Table 1). Aging of the originally water-wet sands significantly reduced their hydraulic conductivity by 42, 33 and 21 % for P2040, P4070 and P100, respectively. It is worth noting that replicate measurements of the permeability and porosity of the aged P2040 sand were carried out and gave exactly the same values as shown in Table 1. The physical properties of the fluids used in the experiments are provided in Table 2.

<table><caption>**Table 2**

Physical properties of fluids used in the study (at 21 °C).</caption><thead><tr><th>Two phase fluid system</th><th>Air-Brine</th><th>Oil-Brine</th><th>Air-Oil</th></tr></thead><tbody><tr><td>Interfacial tension [mNm<sup>-1</sup>]</td><td>62.45<br/>61.56*</td><td>11.61</td><td>24.44</td></tr><tr><td>Fluid phase</td><td>Air</td><td>Oil</td><td>Brine</td></tr><tr><td>Density [kgm<sup>-3</sup>]</td><td>1.2</td><td>813</td><td>1124</td></tr></tbody></table>

## 2.2. Contact angle measurements

In two-phase fluid systems, at equilibrium state, the contact angle $\theta$ is defined by the Young equation (Young, 1805) based on the corresponding interfacial tension of the fluid couple and the interfacial tensions of both fluids at the solid surface (see Fig. 2). $\theta_{gw}$, $\theta_{go}$ and $\theta_{ow}$ are fluid-fluid contact angles within a porous medium and measured through water in the presence of gas, through oil in the presence of gas and through water in presence of oil, respectively. It should be noted that the indices used for the contact angles follow the commonly used

<figure><img src="image_5.png" /><figcaption>Fig. 2. Definition of contact angles in two-phase fluid systems (adapted from Blunt et al., 2021). $\sigma_{ws}$, $\sigma_{gs}$ and $\sigma_{os}$ are the interfacial tension between water and solid, gas and solid, and between oil and solid.</figcaption></figure>

<footnote>* Under thermodynamic equilibrium conditions of brine with oil.</footnote>
nomenclature, but are the reverse of those chosen in the study of Blunt et al. (2021). In three-phase fluid systems, by combining Young's equation written for each of the three fluid pairs, once two contact angles are known, the third can then be predicted using the Bartell-Osterhof equation (Bartell and Osterhof, 1927):

$$ \sigma_{gw} \cos\theta_{gw} = \sigma_{go} \cos\theta_{go} + \sigma_{ow} \cos\theta_{ow} \quad (2) $$

where $\sigma_{gw}(mNm^{-1})$, $\sigma_{go}(mNm^{-1})$ and $\sigma_{ow}(mNm^{-1})$ are the interfacial tensions of fluid couples gas/water, gas/oil and oil/water.

Aged and unaged sand samples were prepared to assess static contact angles in the case of 2-phase gas-brine, gas-oil and oil-brine fluid systems. A 10-cm-long fine glass tube (ID = 0.8 cm) was used as a sample holder. Once the sand sample was prepared, we used the Thermo Scientific™ HeliScan™ µ-CT system for imaging the porous media containing the distribution of two fluids. The specimen was placed at a distance of 12 mm to the X-ray; a voxel size of about 2 µm was chosen in all scans resulting in a total number of 2920 × 2920 × 2160 voxels. The reconstructed X-ray images of the different sand samples were then used to measure manually the contact angles $\theta_{gw}$, $\theta_{go}$ and $\theta_{ow}$ (Table 3) using the ImageJ software (Schneider et al., 2012). The detailed results of the contact angle measurements are shown in the histograms presented in Figs. A1, A2 and A3 (Appendix A). The oil-water contact angles measured for the three aged sand samples P2040ag, P4070ag and P100ag were approximately 99, 95 and 93 degrees, respectively. The increased oil-water contact angles determined for the aged sands clearly show the change in wettability from initially water-wet to moderately oil-wet. Note that in the oilfield terminology, porous media with oil-water contact angles between 70 and 110° are commonly called intermediate wet. However, in our study, we consistently refer to the three aged sand samples with contact angles greater than 90 degrees as oil-wet. As a result of the aging process, the gas-oil contact angles for P2040ag, P4070ag and P100ag decrease considerably, from around 47° to 26°, 46° to 26° and 45° to 27° respectively. The gas-water contact angles are however significantly increased, by a factor of nearly 2, which emphasizes that water now has a lower affinity to the solid grains than in the water-wet sand samples.

## 2.3. Interfacial tension measurements

Characterizing interfacial tension (IFT) is essential in numerous scientific and engineering fields, as it provides insights into fluid interface behavior and phase interactions. In this study, a custom-built apparatus was utilized to measure the IFT of gas-water-oil systems using the pendant/sessile drop technique. IFT measurements were performed between gas-oil ($\sigma_{go}$), gas-brine ($\sigma_{gw}$), and oil-brine ($\sigma_{ow}$) (Table 2). The drop shape analysis was carried out using the Laplace transform method integrated into the Axisymmetric Drop Shape Analysis (ADSA) software (Saraji et al., 2013).

To prepare the fluids, crude oil, brine, and nitrogen were placed in a 1 L beaker and vigorously shaken at 180 rpm (RPM) for 24 h using a VWR Advanced Digital Shaker. After shaking, the solutions were allowed to stabilize and equilibrate before being used in the main interfacial tension experiments. It should be noted that, for gas-oil and gas-brine surface tension (ST) measurements, the cell was filled with the

<table><caption>**Table 3**
In situ static contact angles quantified for the six aged and unaged sands (ag = aged).</caption><thead><tr><th></th><th>P2040</th><th>P2040ag</th><th>P4070</th><th>P4070ag</th><th>P100</th><th>P100ag</th></tr></thead><tbody><tr><td>oil-water</td><td>39.0 ± 8.8</td><td>98.5 ± 12.9</td><td>37.7 ± 5.9</td><td>94.5 ± 12.1</td><td>35.8 ± 5.2</td><td>93.3 ± 9.4</td></tr><tr><td>θ<sub>ow</sub> (°)</td><td>47.0 ± 6.0</td><td>26.4 ± 4.2</td><td>45.8 ± 5.8</td><td>25.8 ± 4.2</td><td>45.3 ± 6.2</td><td>27.5 ± 4.7</td></tr><tr><td>gas-oil</td><td>24.1 ± 5.7</td><td>46.7 ± 6.4</td><td>29.2 ± 4.8</td><td>48.8 ± 7.6</td><td>28.2 ± 6.2</td><td>48.8 ± 5.5</td></tr><tr><td>θ<sub>go</sub> (°)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>gas-water</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>θ<sub>gw</sub> (°)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

equilibrated gas, and the equilibrated oil or brine droplet was injected through a needle of the selected size from the top of the cell. In contrast, for oil-brine interfacial tension measurements, the cell was filled with the equilibrated brine, and an equilibrated oil droplet was suspended from the bottom of the cell through a needle.

Before commencing the tests, the cell was vacuumed for 30 min to remove any contaminants and then filled with the designated fluids, depending on the two phases being measured, under a pressure of 10 psi. To clarify further, during the brine-gas ST measurements, the cell was inverted and thoroughly vacuumed to eliminate any potential contaminants, ensuring a clean environment. Following this, the cell was filled with nitrogen gas ($N_2$) to establish the gas phase. A droplet of equilibrated brine was then carefully injected through a 2.1 mm diameter needle, allowing for precise control over its dimensions. The external surface of the brine droplet, which formed an interface with the surrounding nitrogen gas, was analyzed to measure surface tension. This interface was evaluated using the Young-Laplace equation to accurately determine the ST values.

## 2.4. Experimental setup

Controlled capillary pressure measurements were made in the experimental device shown in Fig. 3. In the hanging water column experiments, a 1-cm thick sand sample was placed on the semi-permeable membrane of a 600 mL sintered glass Buchner funnel Kontes® (inner diameter (ID) = 10 cm, pore sizes 10–15 µm, threshold pressure for gas breakthrough 1 bar) corresponding to a total volume of approximately 78 mL occupied by the porous medium. All experiments were performed at a constant room temperature of 20 ± 1 °C. It is worth noting that the thickness of the sand sample was kept small to minimize the influence of the vertical non-uniform saturation distribution in the sand sample on the quantification of its water content (Moura et al., 2019) and to reduce the effects of gravity on the water retention curve (Ayaz et al., 2020). The diameter of the sand pack was chosen to be large in relation to its thickness in order to increase the pore volume of the sands and thus increase the accuracy of the amount of water that can be drained.

The brine chamber was connected to a brine reservoir via a flexible hose, 2-way and 3-way valves, and a stainless-steel needle (ID = 1 mm). The tip of the needle was dipped into the brine through a hole in the Teflon cap, the diameter of which was just large enough to avoid contact with the needle. In addition, the upper part of the Buchner funnel was covered with a plastic film to reduce the evaporation of water from the sand samples, while maintaining atmospheric pressure. With this arrangement, the evaporation rate from the reservoir was extremely low

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Experimental device of the hanging water column</figcaption></figure>
at less than 0.001 mLh⁻¹ and thus could be neglected when calculating the water saturation of the 1-cm thick sand sample. The brine reservoir was a glass flask set on a digital Mettler Toledo balance PL 1502-S with a precision of 0.01 g. During the primary drainage of the water saturated sand samples, the brine chamber was initially raised in cm-increments relative to the initial level of brine solution in the reservoir to establish a defined capillary pressure, and the mass of drained brine was measured to equilibrium saturation (which was reached in the range of 5 to 10 min during the experiments) using the digital balance. The sequence of capillary pressures was imposed monotonically to produce the primary drainage branch of the capillary-saturation relation. The increase of the liquid level in the brine reservoir during drainage or, in the case of imbibition, the decrease of the brine level was measured with a digital Eberbach® cathometometer telescope. The capillary pressure acting on the sand sample was determined on the basis of the measured height difference Δz between the center of the sand sample and the liquid level in the brine reservoir using a ruler and the cathometometer telescope (see Fig. 3).

Once the maximum capillary pressure level of approximately 110 cm was reached, secondary imbibition of the sand samples was performed by gradually lowering the position of the brine chamber to return to its original starting position. Subsequently, the same experimental procedure was used to determine the capillary pressure-saturation curves during the main drainage and the main imbibition. The total duration of the experiment, which included the four branches, was approximately 12 h.

Based on the measured porosities and densities of the sands quantified in preliminary tests (see Table 1), the individual sand volumes for each of the experiments were fixed. The sand pack was superimposed by a load of 534 g placed on a perforated hard PVC plate to ensure a constant pressure on the sand pack. Prior to the filling procedure, the sintered glass funnel was completely filled with brine solution to avoid any air bubbles below the semi-permeable membrane. The pore volume of each sand pack was individually measured by the weight of brine solution used for its saturation with brine solution, varying between 28 and 37 g corresponding to 24.9 and 32.9 mL, respectively. The packing process of the sands on the semi-permeable membrane was conditioned by the fact that the aged sands were highly water repellent. In order to achieve complete initial water saturation of the sands and a uniform distribution of sand particles where all pores have more or less the same size, the air-dried aged sands were poured onto the membrane in small portions under water and stirred manually to break up the agglomerated sand particles. A thin layer of water was left over the sand pack and a manual tamper was used to compact and achieve a uniform spatial distribution of the sand particles within each sample. Once a final sand layer thickness of approximately 1.1 cm was reached, the perforated rigid PVC sheet was placed and loaded with the ballast to compress the sand pack and maintain a constant pressure on the sand layer throughout the experiment. At this stage, excess water that appeared on the perforated PVC plate was then extracted from the sand pack and measured in order to quantify the exact pore volume of the sand pack.

## 2.5. Analysis of $p_c-S$ curves in a two-phase water-gas system

In this study, we used the model introduced by (Van Genuchten, 1980). It is based on a simple analytical expression to describe the capillary pressure curves. In a two-phase gas-water system, the capillary pressure $p_{cgw}$ (Pa) as function of the effective saturation of the wetting fluid phase (water) is described by

$$p_{cgw} = \frac{1}{\alpha} \left[ S_{we}^{\frac{n}{n-1}} - 1 \right]^{\frac{1}{n}} \quad (3)$$

where $\alpha$ (Pa⁻¹) and $n$ (–) are empirical parameters. Parameter $n$ controls the “curvature” of the $p_c-S$ function; $\alpha$ quantifies the “level” of the

plateau of the capillary pressure. The effective saturation $S_{we}$ (–) of the wetting fluid phase is defined by:

$$S_{we} = \frac{S_w - S_{wi}}{S_{max} - S_{wi}} \quad (4)$$

where $S_w$, $S_{wi}$ and $S_{max}$ are the actual, the irreducible, and the maximum saturation of the wetting fluid phase (here water). During primary drainage of the porous medium, water is displaced by gas; the maximum saturation $S_{max}$ is equal to one.

On the other hand, when gas is displaced by water, in the case of secondary imbibition, $S_{max}$ is limited to values lower than one because the non-wetting fluid phase will remain in the porous medium at residual saturation $S_{gr}$. The maximum water saturation is therefore complementary to the residual saturation of the nonwetting fluid phase: $S_{max} = 1-S_{gr}$.

# 3. Results and discussion

## 3.1. Capillary pressure-saturation relationships

Fig. 4 shows the measured water retention curves for the untreated water-wet sand P100 and oil-wet sand P100ag as examples. The $p_c-S$ curves were determined during the four successive phases of the experiment: Primary drainage, secondary imbibition, main drainage and main imbibition. The major differences between the results of the measured water retention curves become visible when comparing the water-wet sand P100 and oil-wet sand P100ag: The so-called plateau of the capillary pressure is significantly lower with P100ag compared to P100 as in the irreducible water saturation. This will be discussed in more detail later on the basis of the results obtained during the main drainage and the main imbibition. It can also be seen that the $p_c-S$ curves obtained during the secondary imbibition and main imbibition almost overlap.

Figs. 5, 6 and 7 show the main drainage and imbibition branches of the capillary pressure-saturation relationships obtained for the three unaged and aged sands. The effect of hysteresis between main drainage and main imbibition was pronounced in both water-wet, fractionally wet sands and oil-wet sands; the hysteresis loop was closed. Sands becoming oil-wet or fractionally wet strongly lowered the capillary pressure heads at which a given water saturation was reached. Apparently, the absolute hysteresis effect was not significantly reduced.

A “shoulder” effect is noticeable in the water retention curves, which is particularly evident in the drainage curves for water saturations greater than 75 %. It is likely that the load resting on the top of the sand samples caused this premature decrease in water saturation as soon as the pore water pressure decreased. When the capillary pressure is zero, the load exerted on the sample is balanced by the normal stress in the solid matrix and the pore water pressure. The so-called Terzaghi principle states that a decrease in water pressure increases the normal stress in the solid matrix, which is accompanied by a consolidation of the porous medium and thus contributes to an increased desaturation of the sand. The effect of loading on the dewatering curves of aged and unaged sands was clearly noted at high water saturations in this study. For unaged sands, loading was not required to perform controlled $p_c-S$ measurements. However, due to the water-repelling nature of the sand particles and the preferential formation of agglomerates on contact with water, it was necessary to apply a well-defined upper boundary condition to the aged sand packs by means of loading. It is worth noting that in the reference case of unaged sand P4070, where replications were performed, the main drainage plateau was not affected by the application of a load to the sand pack compared to the no load configuration.

The wettability of the sands significantly affected the irreducible water saturation the higher the fraction of oil-wet sand in the sand mixtures, the lower the irreducible water saturation (Fig. 6). This confirms the experimental results of Ustohal et al. (1998) on water retention curves with artificial mixtures of beads composed of pure quartz sand
<figure><img src="image_7.png" /><figcaption>Fig. 4. Water retention curves measured during primary drainage, secondary imbibition, main drainage and main imbibition for (a) water-wet sand P100 and (b) oil-wet sand P100ag. The capillary pressure is expressed as pressure head of the brine solution.</figcaption></figure>

with hydrophilic behavior and pure PFA grains with hydrophobic behavior where the irreducible water saturation correlated roughly with the weighted mean of irreducible water saturations for the single packings. However, our findings do not confirm results of earlier studies by Morrow (1970) and Anderson (1987) where the irreducible wetting phase saturation in glass beads and sand packs was primarily controlled by the pore structure and did not significantly depend on wettability of the porous medium. On the other hand, the wettability did not significantly affect the residual gas saturation (approximately 8 % for sand P2040 compared to 7 % for P2040ag).

The water retention curves measured during main drainage and main imbibition were analyzed using the Van Genuchten (VG) model (Eq. 3). To quantify its input parameters, we used the experimentally determined end-point irreducible water saturations ($S_{wi}$) and maximum water saturation ($S_{max}$) as known. The only remaining model parameters, VG parameters $n$ and $\alpha$, were then obtained by a nonlinear fit of the theoretical $p_c-S$ function against the experimentally observed retention curve using the NLFiT tool of Origin 2019 with the Levenberg-Marquardt iteration algorithm. The optimized values are listed in Tables 4, 5 and

6. All water retention curves could be described reasonably well with the VG model when $\alpha$ and $n$ were optimized. The coefficient of determination ($R^2$) of the sand samples P2040, P4070 and P100 was between 0.82 and 0.89, 0.85 and 0.90 and between 0.83 and 0.94, respectively (see Tables 4, 5 and 6).

In all cases, an increase in the VG parameter $\alpha$ can be observed when the volume fraction of oil-wet sand in the sample increases. The increase in $\alpha$ corresponds to a reduction in the height of the capillary fringe in an unconfined porous aquifer. Analyzing the fitted water retention curves for all water-wet, fractionally wet and oil-wet sands, the two optimized VG parameters seem to indicate a general trend: The higher the volume fraction of oil-wet sand, the higher is $\alpha$ and the lower is $n$. The findings on the influence of volume fraction of oil-wet sand on the reduction of the capillary pressure plateau confirm the experimental results of Motealleh et al. (2013) conducted in oil-wet sand packs. They reported a systematic variation of the capillary pressure curves during the oil-water displacement when the proportion of oil-wet grains varied, with much larger differences than the ones observed in this study. It should be noted that in their study, there was a significantly larger difference
<figure><img src="image_8.png" /><figcaption>Fig. 5. Water retention curves measured and fitted for water-wet coarse sand P2040, oil-wet sand P2040ag and a mixture of both sands during (a) main drainage and (b) main imbibition. The capillary pressure is expressed as pressure head of the brine solution. The solid lines represent the drainage and imbibition curves according to the Van Genuchten model. The fitted Van Genuchten parameters are given in Table 4.</figcaption></figure>

between aged and unaged contact angles, which may have caused the apparent large variation in capillary pressure curves.

While the measured irreducible water saturations of the mixture correlate quite well with the weighted average of the irreducible water saturations of the individual sands, this is not fully achieved for the VG parameter α. As can be seen from Tables 4 and 6, the quantified VG parameters α correlate quite well with the weighted mean values of the unaged and aged sands in the case of a 50 % mixture of P2040 and P100. However, this is not the case for the mixtures of P4070. For example, the water retention curves of P4070ag50 are strongly dominated by the behavior of the oil-wet sand in the mixture (see Fig. 6). The predomi- nance of the aged sand fraction is also reflected in the VG parameter α, which is 1.8 % and 14 % lower than the aged P4070 sand counterparts for the main drainage and main imbibition curves, respectively. At this stage, we have not found a satisfactory explanation for the high influ- ence of the oil-wet grains on the water retention curve.

3.2. Normalization of capillary pressure-saturation relationships

Attempts were made to correlate capillary pressure-saturation curves with the physicochemical properties of the porous medium. Leverett (1941) was the first to introduce a dimensionless capillary pressure correlation function. To summarize capillary pressure data into a single curve for porous media with various irreducible water saturations, Desouky (2003) extended Leverett’s approach by using effective water saturations S<sub>we</sub>. In the case of a gas-water fluid system, the so-called Leverett J-function can then be written as:

$$
J = \frac{p_c}{\sigma_{gw} \cos \theta_{gw}} \sqrt{\frac{k}{\phi}} \qquad (5)
$$

where p<sub>c</sub> (Pa<sup>-1</sup>) denotes capillary pressure, σ<sub>gw</sub> (Nm<sup>-1</sup>) is the interfacial tension, θ<sub>gw</sub> (°) is the contact angle, k is the intrinsic permeability (m<sup>2</sup>), and φ is the porosity (—). The effective water saturation is expressed by Eq. 4, which normalizes the water saturation with respect to the indi- vidually observed irreducible saturation S<sub>wi</sub> and maximum water
<figure><img src="image_9.png" /><figcaption>Fig. 6. Water retention curves measured and fitted for water-wet medium sand P4070, oil-wet sand P4070 ag and three mixtures of both sands containing 25, 50 and 75 volume fractions of P4070ag during (a) main drainage and (b) main imbibition. The capillary pressure is expressed as pressure head of the brine solution. The solid lines represent the drainage and imbibition curves according to the Van Genuchten model. The fitted Van Genuchten parameters are given in Table 5.</figcaption></figure>

saturation ($S_{max}$), resulting in a value between zero and one. Eq. (5) accounts for change of permeability, porosity and wettability of a porous medium as long as the general pore geometry remains constant. This fact was confirmed by Rose and Bruce (1949) and Brown (1951).

Eq. 5 was applied to the measured main drainage and imbibition branches of the retention curves obtained for the six aged and unaged sands (P2040, P2040ag, P4070, P4070ag, P100 and P100ag), taking into account the experimentally quantified values for $S_{wi}$ and $S_{max}$ (see Tables 4–6), intrinsic permeability $k$ and porosity $\phi$ (see Table 1), gas-water interfacial tension $\sigma_{gw}$ (see Table 2) and gas-water contact angle $\theta_{gw}$ (see Table 3). Fig. 8 illustrates the obtained (dimensionless) Leverett J-function as function of the effective water saturation. Each of these individual curves should theoretically coalesce in one unified J-curve that collectively considers all the sands. However, for sands P2040ag and P100ag the collapse is not complete. To quantify the observed deviations, we employed a detailed analysis using the normalized capillary pressure curves of the three unaged sands as reference curves.

To establish a relationship that better matches the properties of the individual experimental Leverett J-functions, we used a curve fitting

technique to model the relationship between the observed data of the three unaged sands and a theoretical function. Several attempts were made to find a theoretical function that most accurately reflects the key characteristics of the data sets. For the drainage case, the best results were achieved with a fit function that combines the properties of an exponential function and a cubic polynomial function and is based on six optimization parameters:

$$J_{\text{fit}}(S_{\text{we}}) = a \exp^{b} S_{\text{we}} + cS_{\text{we}}^{3} + dS_{\text{we}}^{2} + eS_{\text{we}} + f \quad (6)$$

In contrast, for the imbibition case, the same general expression was applied, but the cubic term ($cS_{we}^3$) was omitted, as it was found to not contribute significantly to the data fit. Thus, the best-fit function for the imbibition case simplified to an exponential and quadratic function with five optimization parameters.

To minimize the difference between the experimental data and the fitted values obtained with Eq. (6), we employed the curve_fit function from the scipy.optimize library in Python. The six parameters of $J_{fit}$ were iteratively adjusted to obtain the most accurate representation of the experimental data (see Fig. 8).
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Water retention curves measured and fitted for water-wet fine sand P100, oil-wet sand P100 ag and a mixture of both sands during (a) main drainage and (b) main imbibition. The capillary pressure is expressed as pressure head of the brine solution. The solid lines represent the drainage and imbibition curves according to the Van Genuchten model. The fitted Van Genuchten parameters are given in Table 6.</figcaption></figure>

**Table 4**

Measured irreducible water saturation $S_{wi}$ and maximum water saturation $S_{max}$ of the main drainage and imbibition branches of the retention curves obtained for sand P2040, sand P2040ag and the 50 % mixture, as well as optimized numerical values of the Van Genuchten model.

<table><thead><tr><th rowspan="2"></th><th colspan="2">P2040</th><th colspan="2">50% mixture</th><th colspan="2">P2040ag</th></tr><tr><th>drainage</th><th>imbition</th><th>Drainage</th><th>imbition</th><th>drainage</th><th>imbition</th></tr></thead><tbody><tr><td>S<sub>wi</sub><sup>a</sup></td><td>0.081</td><td>0.081</td><td>0.046</td><td>0.046</td><td>0.029</td><td>0.029</td></tr><tr><td>S<sub>max</sub><sup>a</sup></td><td>0.921</td><td>0.915</td><td>0.933</td><td>0.930</td><td>0.933</td><td>0.931</td></tr><tr><td>α (Pa<sup>-1</sup>)<sup>b</sup></td><td>6.2 × 10<sup>-4</sup> ± 7.3 × 10<sup>-5</sup></td><td>8.9 × 10<sup>-4</sup> ± 1.5 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-3</sup> ± 2.6 × 10<sup>-4</sup></td><td>1.9 × 10<sup>-3</sup> ± 6.5 × 10<sup>-4</sup></td><td>1.6 × 10<sup>-3</sup> ± 5.3 × 10<sup>-4</sup></td><td>3.7 × 10<sup>-3</sup> ± 1.5 × 10<sup>-3</sup></td></tr><tr><td>n<sup>b</sup></td><td>5.02 ± 0.32</td><td>4.41 ± 0.32</td><td>3.89 ± 0.34</td><td>3.21 ± 0.29</td><td>3.35 ± 0.31</td><td>2.63 ± 0.20</td></tr><tr><td>R<sup>2</sup></td><td>0.89</td><td>0.85</td><td>0.86</td><td>0.82</td><td>0.84</td><td>0.84</td></tr></tbody></table>

The values obtained for the optimized fit parameters with the associated standard deviations and their correlation coefficients are summarized in Tables B1 and B2 (Appendix B). The correlation matrix shows

the relationships between the different pairs of parameters. For example, in the case of the main drainage branch, parameters *c* and *e* are highly correlated (correlation coefficient of 0.948), which means that if one

<footnote>a Measured.</footnote>

<footnote>b Optimized Van Genuchten parameters (see Eq. 3).</footnote>
<caption>**Table 5**</caption>

Measured irreducible water saturation $S_{wi}$ and maximum water saturation $S_{max}$ of the main drainage and imbibition branches of the retention curves obtained for sand P4070, sand P4070ag and three mixtures of both sands containing 25, 50 and 75 volume fractions of P4070ag, as well as optimized numerical values of the Van Genuchten model.

<table><thead><tr><th rowspan="2"></th><th colspan="2">P4070</th><th colspan="2">75 % unaged - 25 % aged mixture</th><th colspan="2">50 % mixture</th></tr><tr><th>drainage</th><th>imbition</th><th>drainage</th><th>imbition</th><th>drainage</th><th>imbition</th></tr></thead><tbody><tr><td>S<sub>wi</sub><sup>a</sup></td><td>0.081</td><td>0.081</td><td>0.069</td><td>0.069</td><td>0.054</td><td>0.054</td></tr><tr><td>S<sub>max</sub><sup>a</sup></td><td>0.898</td><td>0.875</td><td>0.826</td><td>0.809</td><td>0.902</td><td>0.855</td></tr><tr><td>α (Pa<sup>-1</sup>)<sup>b</sup></td><td>3.6 × 10<sup>-4</sup> ± 4.0 × 10<sup>-5</sup></td><td>4.9 × 10<sup>-4</sup> ± 9.4 × 10<sup>-5</sup></td><td>4.2 × 10<sup>-4</sup> ± 5.4 × 10<sup>-5</sup></td><td>7.8 × 10<sup>-4</sup> ± 2.0 × 10<sup>-4</sup></td><td>5.6 × 10<sup>-4</sup> ± 8.3 × 10<sup>-5</sup></td><td>9.4 × 10<sup>-4</sup> ± 2.2 × 10<sup>-4</sup></td></tr><tr><td>n<sup>b</sup></td><td>5.18 ± 0.36</td><td>4.47 ± 0.42</td><td>4.67 ± 0.34</td><td>4.19 ± 0.41</td><td>3.82 ± 0.26</td><td>3.09 ± 0.23</td></tr><tr><td>R<sup>2</sup></td><td>0.89</td><td>0.85</td><td>0.88</td><td>0.87</td><td>0.90</td><td>0.90</td></tr><tr><td colspan="7">25 % unaged - 75 % aged mixture<sup>a</sup></td></tr><tr><th></th><th colspan="2">P4070ag</th><th colspan="2"></th><th colspan="2"></th></tr><tr><th></th><th>drainage</th><th>imbition</th><th>drainage</th><th>Imbibition</th><th colspan="2"></th></tr><tr><td>S<sub>wi</sub><sup>a</sup></td><td>0.045</td><td>0.045</td><td>0.034</td><td>0.034</td><td colspan="2"></td></tr><tr><td>S<sub>max</sub><sup>a</sup></td><td>0.878</td><td>0.859</td><td>0.903</td><td>0.894</td><td colspan="2"></td></tr><tr><td>α (Pa<sup>-1</sup>)<sup>b</sup></td><td>5.6 × 10<sup>-4</sup> ± 8.4 × 10<sup>-5</sup></td><td>9.1 × 10<sup>-4</sup> ± 2.3 × 10<sup>-4</sup></td><td>5.7 × 10<sup>-4</sup> ± 8.6 × 10<sup>-5</sup></td><td>1.1 × 10<sup>-3</sup> ± 3.0 × 10<sup>-4</sup></td><td colspan="2"></td></tr><tr><td>n<sup>b</sup></td><td>3.74 ± 0.25</td><td>3.09 ± 0.25</td><td>3.69 ± 0.25</td><td>2.78 ± 0.21</td><td colspan="2"></td></tr><tr><td>R<sup>2</sup></td><td>0.89</td><td>0.88</td><td>0.89</td><td>0.89</td><td colspan="2"></td></tr></tbody></table>

<caption>**Table 6**</caption>

Measured irreducible water saturation $S_{wi}$ and maximum water saturation $S_{max}$ of the main drainage and imbibition branches of the retention curves obtained for sand P100, sand P100ag and the 50 % mixture, as well as optimized numerical values of the Van Genuchten model.

<table><thead><tr><th rowspan="2"></th><th colspan="2">P100</th><th colspan="2">50 % mixture</th><th colspan="2">P100ag</th></tr><tr><th>drainage</th><th>imbition</th><th>drainage</th><th>imbition</th><th>drainage</th><th>imbition</th></tr></thead><tbody><tr><td>S<sub>wi</sub><sup>a</sup></td><td>0.064</td><td>0.064</td><td>0.054</td><td>0.054</td><td>0.035</td><td>0.035</td></tr><tr><td>S<sub>max</sub><sup>a</sup></td><td>0.861</td><td>0.863</td><td>0.860</td><td>0.828</td><td>0.874</td><td>0.873</td></tr><tr><td>α (Pa<sup>-1</sup>)<sup>b</sup></td><td>2.6 × 10<sup>-4</sup> ± 1.7 × 10<sup>-5</sup></td><td>3.8 × 10<sup>-4</sup> ± 2.6 × 10<sup>-5</sup></td><td>3.8 × 10<sup>-4</sup> ± 3.8 × 10<sup>-5</sup></td><td>5.6 × 10<sup>-4</sup> ± 7.7 × 10<sup>-5</sup></td><td>4.7 × 10<sup>-4</sup> ± 4.0 × 10<sup>-5</sup></td><td>8.7 × 10<sup>-4</sup> ± 1.6 × 10<sup>-4</sup></td></tr><tr><td>n<sup>b</sup></td><td>7.00 ± 0.53</td><td>5.58 ± 0.30</td><td>6.69 ± 0.53</td><td>5.21 ± 0.40</td><td>5.11 ± 0.28</td><td>3.82 ± 0.28</td></tr><tr><td>R<sup>2</sup></td><td>0.87</td><td>0.94</td><td>0.83</td><td>0.87</td><td>0.92</td><td>0.88</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Water retention curves normalized using Eq. 5 for the six selected aged and unaged sands and compared to the fit function (Eq. 6): (a) main drainage branch, (b) main imbibition branch.</figcaption></figure>

parameter increases (or decreases), the other is likely to increase (decrease) as well. On the other hand, parameters c and d are negatively correlated (correlation coefficient of -0.990), which means that if one parameter increases, the other is likely to decrease.

To quantify the effectiveness of the prediction of the measured J-values, we quantified two coefficients: the coefficient of residual mass (CRM) and the coefficient of determination R-squared. CRM is defined as

$$CRM = \frac{\sum_{i=1}^{n} M_i - \sum_{i=1}^{n} P_i}{\sum_{i=1}^{n} M_i} \qquad (7)$$

where $M_i$ and $P_i$ represent the measured and predicted J-values, and n is

<footnote>a Measured.</footnote>

<footnote>b Optimized Van Genuchten parameters (see Eq. 3).</footnote>

<footnote>a Measured.</footnote>

<footnote>b Optimized Van Genuchten parameters (see Eq. 3).</footnote>
the number of experimental data. The CRM illustrates the tendency of a model to overestimate or underestimate the measurements (positive values for CRM indicate that the model underestimates the measurements and negative values for CRM indicate a tendency to overestimate).

In the present study, the quantified coefficient of determination (R-squared) is expressed by:

$$R^2 = 1 - \frac{\sum_{i=1}^{n} [M_i - P_i]^2}{\sum_{i=1}^{n} [M_i - \bar{M}]^2}, \text{ with } \bar{M} = \frac{1}{n} \sum_{i=1}^{n} M_i \quad (8)$$

R-squared provides a measure of how well the measurements are likely to be predicted by the fit function. A positive value of R-squared indicates a reasonable fit, and the closer it is to one, the better the fit function captures the variation in the data. The best possible score is one. It is worth noting that R-squared can be negative if the fit function is arbitrarily bad.

Figs. 9 and 10 show the detailed comparison of each data set with the quantified fit function for the main drainage branch and main imbibition branch. The quality of the prediction is illustrated by the two coefficients, CRM and R-squared. The R-squared values generally indicate a good fit for the three unaged sands and P4070ag. However, negative values of R-squared are recorded for sands P2040ag and P100ag. The Leverett J-function of the aged medium-sized sand P4070 is close to the fit function, while the normalized $p_c$-S curves of the other two aged sands are largely below the fit function (expressed by high negative CRM values). Using the measured data presented in Tables 1 and 3 and calculating the Leverett scaling factors $\sqrt{\frac{k}{\partial}} \frac{1}{\cos n\theta_{gw}}$ for each pair of a ratio between aged and unaged sands, the normalized Leverett J-functions should differ by a factor of 1.07 for the P2040 sands, 1.23 for the P4070 sands, and 1.28 for the P100 sands. However, if we take the $p_c$ values measured at a water saturation of 0.5 as the reference value for the capillary pressure plateaus, we obtain $p_c$ ratios of unaged sand to aged sand of approximately 2.5 for P2040 sands, 1.4 for P4070 sands and 2.0 for P100 sands in the case of the main drainage. Only for P4070, the Leverett scaling factor is quite close to the estimated $p_c$ ratios, which explains why the normalized capillary pressure curves of P2040ag and P100ag are lower than the fit function (expressed by negative CRM values) and do not coincide with the other curves. One possible reason could be the assumption that the relationship between permeability and porosity was taken into account when normalizing the capillary pressure curves for the aged sands. Using the permeability to represent the pore size and porosity data in the Leverett scaling factor assumes an isomorphic pore structure of the aged sands, characterized by uniform pore size, uniform pore distribution, and the same geometric shape of the pores. In the present study, the pore structure of the aged sands appears to be significantly different from that of the unaged sands, as both permeability and porosity have decreased due to the aging process.

The fact that the soil particle densities of the aged sands are lower than those of the unaged sands indicates that the particle size distribution of the aged sands is not exactly the same, as the particle size of the individual silica grains has increased due to the aging process. Possible changes in the pore structure and texture, which are probably more pronounced in the P2040 and P100 sands, have certainly superimposed the influence of wettability, expressed by the ratio of the cosine values of the measured gas-water contact angles and quantified with a factor of about 1.3. Since the permeability of the aged sands in this study appears to be a poor approximation for the characteristic pore size, we adopted the same ratio of permeability to porosity for the aged sands P2040 and P100 as for the unaged sands in a modified normalized Leverett J-function. As shown in Figs. 9 and 10, the normalized capillary curves of P2040ag, labeled “Test 1”, are now much closer to the other normalized curves. In the case of P100ag, the CRM values improved insignificantly from -0.77 to -0.70 for the drainage branch and from -1.28 to -1.19

for the imbibition branch.

A second test was carried out to reduce the observed differences. Using test 1, the gas-water contact angles were increased by the measured standard deviation. In the case of P2040ag, the contact angle was increased by 6.4°, while the contact angle of P100ag was increased by 6.2°. The normalized capillary curves of P2040 obtained in Test 2 still underestimate the values of the fit function, but the deviations are further reduced, and the quantified R-squared value of 0.69 and 0.72 for the drainage branch and the imbibition branch respectively, reflects a reasonable calibration quality. For P100ag, however, the R-squared values remain negative. The reduction in the magnitude of the negative CRM values in Test 1 and Test 2 supports the hypothesis that these adjustments mitigate the problem of global overestimation.

### 3.3. Gas-water contact angles

Table 3 summarizes the measured contact angles obtained for the three unaged and three aged sands. The contact angles determined for oil-water fluid system indicate that the three aged sands are now intermediate-wet as the $\theta_{ow}$ contact angle is slightly higher than 90°. The corresponding $\theta_{gw}$ contact angles are increased indicating that water becomes less wetting with regard to air.

Based on Eq. (2) and using the measured gas-oil and oil-water contact angles for each of the six aged and unaged sands and the measured IFT values (see Table 2), gas-water contact angles can be predicted

$$\theta_{gw} = \cos^{-1} \left[ \frac{\sigma_{go} \cos\theta_{go} + \sigma_{ow} \cos\theta_{ow}}{\sigma_{gw}} \right] \quad (9)$$

Uncertainties of the predicted gas-water contact angles were quantified from Eq. (9) based on the total derivative expansions for the correlated variables $\theta_{go}$ and $\theta_{ow}$ as follows:

$$\Delta\theta_{gw} = \frac{1}{\sqrt{1 - \cos^2\theta_{gw}}} \left[ \frac{\sigma_{go} \sin\theta_{go} \Delta\theta_{go} + \sigma_{ow} \sin\theta_{ow} \Delta\theta_{ow}}{\sigma_{gw}} \right] \quad (10)$$

where $\Delta\theta_{go}$ and $\Delta\theta_{ow}$ are the uncertainties of the measured contact angles between gas-oil and oil-water, respectively.

Table 7 summarizes the static gas-water contact angles predicted by eqs. (9) and (10) for each of the six aged and unaged sands and compared with the measured gas-water contact angles. It can be seen that the predicted gas-water contact angles are significantly overestimated by the theoretical relationship (Eq. 9). However, the increased contact angles of the aged sands compared to the originally water-wet sands are well predicted. A possible reason for the overestimation could be the fact that the Bartell-Osterhof equation is only valid if all three fluid phases are present and in thermodynamic equilibrium. When three fluid phases are present simultaneously, the interfacial tensions differ from those in two-fluid phase systems. For example, the gas-water interfacial tension can be significantly lower than its two-phase value because oil can cover the interface with a thin film of molecular thickness (Zhou and Blunt, 1997). The other two-phase interfacial tensions can also vary when the third phase is present. When measuring the gas-brine interface in the presence of oil, the value decreased only insignificantly from the original 62.45 mNm⁻¹ to 61.56 mNm⁻¹.

Another possible reason for the theoretical values being higher than the observed contact angles could be the ideal conditions of the Young equation used in the Bartell-Osterhof equation: Young’s equation is valid only for ideal solids whose surfaces are homogeneous, isotropic, smooth, rigid, and are inert to surrounding fluids (Drelich, 2019). To take into account the roughness of the solid surface, Wenzel (1949) proposed the following relationship to express the measured apparent contact angle ($\theta^a$) as function of the intrinsic contact angle ($\theta^i$), which is defined by Young’s law (Young, 1805) and measured on an ideal, smooth and inert surface of the same material

$$\cos\theta^a = r \cos\theta^i, \quad (11)$$
<figure><img src="image_12.png" /><figcaption>Fig. 9. Analysis of the main drainage branch: comparison of each data set with the optimized fit function. The quality of the prediction of the measured experimental data with the fitted function is expressed by the coefficients CRM and R-squared.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Analysis of the main imbibition branch: comparison of each data set with the optimized fit function. The quality of the prediction of the measured experimental data with the fitted function is expressed by the coefficients CRM and R-squared.</figcaption></figure>
<table><caption>**Table 7**
Predicted and measured gas-water contact angles for each of the six aged and unaged sand samples.</caption><thead><tr><th></th><th>P2040</th><th>P2040ag</th><th>P4070</th><th>P4070ag</th><th>P100</th><th>P100ag</th></tr></thead><tbody><tr><td>Predicted gas-water<br>θ<sub>gw</sub> [°]</td><td>65.7 ± 3.0</td><td>71.2 ± 3.3</td><td>65.2 ± 2.5</td><td>70.3 ± 3.1</td><td>64.8 ± 2.5</td><td>70.4 ± 2.8</td></tr><tr><td>Predicted gas-water<br>θ<sub>gw</sub> [°]<sup>*</sup></td><td>65.3 ± 3.1<sup>*</sup></td><td>70.9 ± 3.3<sup>*</sup></td><td>64.8 ± 2.6<sup>*</sup></td><td>70.0 ± 3.2<sup>*</sup></td><td>64.5 ± 2.6<sup>*</sup></td><td>70.1 ± 2.8<sup>*</sup></td></tr><tr><td>Predicted<br>gas-water<br>θ<sub>gw</sub> [°]<sup>**</sup></td><td>31.2 ± 5.4<sup>**</sup></td><td>47.8 ± 4.2<sup>**</sup></td><td>29.2 ± 4.8<sup>**</sup></td><td>45.4 ± 4.2<sup>**</sup></td><td>27.8 ± 5.0<sup>**</sup></td><td>45.7 ± 3.7<sup>**</sup></td></tr><tr><td>Measured gas-water<br>θ<sub>gw</sub> [°]</td><td>24.1 ± 5.7</td><td>46.7 ± 6.4</td><td>29.2 ± 4.8</td><td>48.8 ± 7.6</td><td>28.2 ± 6.2</td><td>48.8 ± 5.5</td></tr></tbody></table>

where $r$ is the roughness factor of the solid surface, which is generally defined as the true coarse surface area divided by the projected area. The numerical value of the roughness factor is larger than unity for rough surfaces, and is equal to one for smooth surfaces.

In the present experimental study, the contact angle measurements along the three-phase contact line were performed on a curved surface of the solid grains, which was assumed to be smooth. Therefore, to verify Eq. 11, the gas-water contact angles previously predicted by the Bartell-Osterhof equation were considered as intrinsic gas-water contact angles ($\theta_{gw}^i$) and the gas-water contact angles measured in the experiments were considered as measured apparent contact angles ($\theta_{gw}^a$). Assuming that the roughness of the different solid grains is uniform, the cosine of $\theta_{gw}^a$ was then plotted against the cosine of $\theta_{gw}^i$ for each of the six sand samples (Fig. 11). The data bloc with three high values of cosines corresponds to the contact angles measured for the three unaged sand samples. As can be seen from the cross plot (Fig. 11), the data set as a whole is quite well described by a linear regression giving a fit result with an R-squared of 0.91; the slope of the regression line corresponds to the roughness factor of Eq. 11. The obtained calibrated roughness factor of 2.05 was then used, based on Eq. 11, to estimate corrected values of the predicted gas-water angle for each of the six sand samples (Table 7). It is certainly not surprising that the calibrated contact angles are now quite close to their measured counterparts.

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Cross-plot of the cosine of $\theta_{gw}^a$ against the cosine of $\theta_{gw}^i$ for each of the six aged and unaged sand samples. The data bloc with three high values of cosines corresponds to the contact angles measured for the three unaged sand samples.</figcaption></figure>

## 4. Summary and conclusions

The influence of wettability on capillary pressure-saturation relationships in two- or three-phase fluid systems, especially in solid rocks such as sandstones or carbonates, has been the subject of numerous experimental investigations. Only few experimental studies in the past have focused on the effects of wettability on water retention curves in unconsolidated porous media. This was the motivating starting point for the present work, especially since wettability is a key factor controlling multiphase flow behavior in porous media, particularly in applications in environmental and engineering sciences, such as soil remediation, CO₂ sequestration and groundwater management.

Our experimental investigations in quartz sands show that wettability significantly influences capillary pressure and fluid distribution. Sands that became oil-wet or fractionally wet greatly reduced the capillary pressure heights at which a certain water saturation was reached; the absolute hysteresis effect was not significantly reduced. In the originally water-wet sands, the pressure for displacement of water by air was high, as water adheres to the pore walls as a wetting phase. In our oil-wet sands, characterized by a higher gas-water contact angle, the capillary pressure-saturation curve was more gradual. The air as a nonwetting phase spread over the sand surfaces so that the water could be displaced more easily, resulting in lower capillary pressures at similar saturations. This is consistent with previous studies on solid rock (Anderson, 1987; Morrow, 1990) and with recent observations of two-phase water-oil systems at the pore level using X-ray µ-CT (e.g. Arshadi et al., 2020), which state that the shape of the $p_c-S$ curve changes as the wettability of the solid changes from water-wet to oil-wet, affecting the residual oil and water saturations. Changes in the wettability of our sands also affected their irreducible water saturation. The higher the proportion of oil-wet sand in the sand mixtures, the lower the irreducible water saturation. This is due to the fact that in oil-wet sands the gas-water contact angle is increased and water tends to remain in the larger pores and thus be more easily displaced by air, resulting in a lower irreducible water saturation of the sand sample.

Quantitative analysis of the measured water retention curves for water-wet, fractionally wet and oil-wet sands was performed using the Van Genuchten model. Based on the irreducible water saturation and residual air saturation observed for each main drainage and main imbibition branch of the water retention curve, the model parameters $n$ and $\alpha$, were then determined by a nonlinear fit of the theoretical $p_c-S$ function to the experimentally observed water retention curve. In all cases, an increase in the VG parameter $\alpha$ can be observed when the volume fraction of oil-wet sand in the sample increases. The increase in $\alpha$ corresponds to a reduction in the height of the capillary fringe in an unconfined porous aquifer. When analyzing the fitted water retention curves for all water-wet, fractionally wet and oil-wet sands, the two optimized VG parameters seem to show a general trend: The higher the volume fraction of oil-wet sand, the higher is $\alpha$ and the lower is $n$.

<footnote>* $\sigma_{gw} = 61.56 \text{ mNm}^{-1}$.</footnote>

<footnote>** $\theta_{gw}$ calibrated using roughness factor $r = 2.05$ (Eq. 11).</footnote>
We tested the scalability of the measured main drainage and imbibition branches of the retention curves obtained for the six aged and unaged sands (P2040, P2040ag, P4070, P4070ag, P100 and P100ag). The extended Leverett J-function (Eq. 5) using effective water saturations was applied to the two measured branches of the retention curves. For the three unaged sands and the aged P4070 sand, it was found that each of the two branches of the measured retention curves can be scaled to a single curve if, in addition to the petrophysical parameters (intrinsic permeability, porosity) and surface tension, the gas-water contact angle, irreducible residual water saturation and residual air saturation are also taken into account. Theoretically, each of the six individual curves should have collided into a uniform J-curve that takes all sands together into account. However, the collapse of the retention curves into a unified J-curve was not fully achieved for the aged P2040 and P100 sands. To quantify the observed deviations, we calculated as a reference, a theoretical fit function that best matches the properties of the individual experimental Leverett J-functions of the three unaged sands. The experimental data sets for P2040ag and P100ag were found to be overestimated overall by the fit function. However, when the petrophysical parameters of the unaged sands were used instead of the actual measured parameters, the individual experimental Leverett J-functions came closer to the uniform J-curve. From this it could be concluded that the apparent differences in pore structure between aged and unaged sands in addition to wettability, expressed by the cosine of the contact angle, may have contributed to a further reduction in the capillary pressure plateaus of the aged sands, which was particularly visible and significant in the P2040 and P100 sand.

Our experimental results also shed light on the ability to predict measured gas-water contact angles from the contact angles we measured in situ for the other two fluid pairs. To our knowledge, there are no comparable studies in the literature in which all three pairs of interfacial tensions and all three pairs of static contact angles have been measured in detail. One of the rare publications comparing the predicted contact angles from the other two is the review paper by Siddiqui et al. (2018) on the measurement of oil-water contact angle in shale reservoirs. In their approach to the Bartell-Osterhof equation, they used data from the literature on the corresponding interfacial tension for the fluid pairs when no measured data were available, measured air-water contact angles and assumed that the gas-oil contact angle is 0° (i.e. complete spreading of the oil on the solid surface in the presence of air). The oil-water contact angles they predicted for two shale formations differ by about three to five degrees from their measured counterparts. However, when applying the Bartell-Osterhof equation (Eq. 9) to our data, the predicted values for the gas-water contact angle were much higher, about 20 to 25° higher than the measured values.

What could be the reason for the significant overestimation of the measured gas-water contact angles? We found that the 3-phase equilibrium condition for interfacial tensions was not the reason. We recognized as a physical explanation that the roughness of the soil could be the actual cause, since Young's law, which underlies the Bartell-Osterhof equation, assumes ideal smooth and inert solid surfaces. We tested the approach proposed by Wenzel (1949), which considers a roughness factor that converts the ideal contact angle (used in Young's law) into an apparent measured contact angle. In our study, we assumed that the roughness factor was the same for all six aged and unaged sands and corrected the original predicted values for the gas-water contact angles of the six sands with a fitted roughness factor of 2.05. As expected, this allowed a near perfect overlap between the predicted contact angles and the measured values. The quantified roughness factor, defined as ratio of the real rough surface of the sand grain and the projected grain surface appears to be physically plausible. Indeed, recent studies have been conducted to measure the roughness of quartzite grains with relatively spherical local shapes using the Fogale Nanotech

optical microscope, where the surface topography was described by an interferogram (Yang and Baudet, 2016; Yang et al., 2016). They investigated three uniform particle size distributions varying between 0.6 and 1.18, 1.18–2 and 2–5 mm; the surface measurement were always performed on a sub-area of the grain surface with dimensions 106.6 µm × 106.6 µm. Yang and Baudet (2016) reported average values of estimated surface areas for a grain size distribution of 1.18–2 mm of approximately 26,000 µm², corresponding to an equivalent roughness factor of 2.29. In further studies, Yang et al. (2016) presented detailed results for the grain size distribution of 0.6–1.1 mm, where the estimated surface area of the grain was even higher, approaching 300,000 µm². However, the extent to which a constant roughness factor can be assumed for all six sands in our study remains to be seen. In particular, as shown by Yang et al. (2016), the grain size can have an influence on the surface roughness, and furthermore, there are still no findings on how aging of quartzite sand grains changes the roughness of the irregularities distributed on the grain surface. Further detailed investigations under the optical microscope could provide information on this.

The new insights gained in this study could be used in advanced numerical codes to improve their predictive capability by implementing the possibility of wettability-dependent two-phase capillary pressure curves that can change in time at a given location. This can be of great importance in field applications in environmental science, where the wettability of the porous medium influences the dispersion of contaminants (e.g. hydrocarbons) in sandy soils and aquifers. Oil-wet conditions can be achieved in old, long-standing oil spills in soil, which can hinder the removal of oil spills as the oil adheres to the sand grains and resists flushing with water.

### CRediT authorship contribution statement

**Gerhard Schäfer:** Writing – original draft, Visualization, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Kevin Hernandez Perez:** Writing – review & editing, Visualization, Software, Investigation. **Panav Hulsurkar:** Writing – review & editing, Methodology, Investigation, Data curation. **Magda Ibrahim Youssif:** Writing – review & editing, Investigation. **François Lehmann:** Writing – review & editing, Formal analysis, Data curation. **Mohammad Piri:** Resources, Project administration, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgment

This research was granted from the French National Centre for Scientific Research (CNRS) within the International Research Project « CONTINUUM » established between ITES and COIFPM. We are also grateful for the financial support by the University of Strasbourg and the administrative and technical support by the University of Wyoming. We extend our gratitude to Wattana Chaisoontornyotin, Henry Plancher and Van Si Le for their technical and scientific support in planning and conducting the experiments. We also acknowledge Dr. Soheil Saraji, who enriched the discussion of our IFT measurements under thermodynamic equilibrium of three fluid phases. We would also like to thank Martine Trautmann of the EOST soil analysis laboratory, University of Strasbourg for measuring the particle size curves of the sands using the laser diffraction method.
<figure><img src="image_15.png" /><figcaption>Fig. A.1. Histogram of in-situ contact angles measured in sands P2040 and P2040ag. The red dash dot line shown in each histogram represents the calculated normal distribution whose mean value and standard deviation are summarized in Table 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. A.2.** Histogram of in-situ contact angles measured in sands P4070 and P4070ag. The red dash dot line shown in each histogram represents the calculated normal distribution whose mean value and standard deviation are summarized in Table 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. A.3. Histogram of in-situ contact angles measured in sands P100 and P100ag. The red dash dot line shown in each histogram represents the calculated normal distribution whose mean value and standard deviation are summarized in Table 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
## Appendix B. Quantified parameters of the fitted J-function

Table B1 and Table B2 summarize the optimized parameters of function $J_{fit}$ (Eq. 6) quantified for both the main drainage branch and main imbibition branch of the experimental data and the associated correlation matrix. The tables specify uncertainty intervals (standard errors) associated with the parameter estimates.

### Table B1

Optimized parameters of $J_{fit}$ and associated correlation matrix quantified for the main drainage branch.

<table><thead><tr><th colspan="7">Fit parameter</th></tr><tr><th>a<br>amplitude</th><th>b<br>exponential coefficient</th><th>c<br>cubic coefficient</th><th>d<br>quadratic coefficient</th><th>e<br>linear coefficient</th><th>f<br>constant offset</th></tr></thead><tbody><tr><td>2.23 ± 0.01</td><td>-34.3 ± 0.3</td><td>-0.86 ± 0.1</td><td>0.31 ± 0.2</td><td>-0.19 ± 0.09</td><td>0.74 ± 0.01</td></tr><tr><th colspan="7">Correlation matrix</th></tr><tr><td></td><td>a</td><td>b</td><td>c</td><td>d</td><td>e</td><td>f</td></tr><tr><td>a</td><td>1.</td><td>0.702</td><td>0.778</td><td>-0.828</td><td>0.880</td><td>-0.915</td></tr><tr><td>b</td><td>0.702</td><td>1.</td><td>0.671</td><td>-0.730</td><td>0.802</td><td>-0.888</td></tr><tr><td>c</td><td>0.778</td><td>0.671</td><td>1.</td><td>-0.990</td><td>0.948</td><td>-0.836</td></tr><tr><td>d</td><td>-0.828</td><td>-0.730</td><td>-0.990</td><td>1.</td><td>-0.983</td><td>0.893</td></tr><tr><td>e</td><td>0.880</td><td>0.802</td><td>0.948</td><td>-0.983</td><td>1.</td><td>-0.953</td></tr><tr><td>f</td><td>-0.915</td><td>-0.888</td><td>-0.836</td><td>-0.893</td><td>-0.953</td><td>1.</td></tr></tbody></table>

### Table B2

Optimized parameters of $J_{fit}$ and associated correlation matrix quantified for the main imbibition branch.

<table><thead><tr><th colspan="7">Fit parameter</th></tr><tr><th>a<br>amplitude</th><th>b<br>exponential coefficient</th><th>d<br>quadratic coefficient</th><th>e<br>linear coefficient</th><th>f<br>constant offset</th></tr></thead><tbody><tr><td>2.44 ± 0.01</td><td>-50.6 ± 0.3</td><td>-0.75 ± 0.03</td><td>0.24 ± 0.03</td><td>0.48 ± 0.01</td></tr><tr><th colspan="7">Correlation matrix</th></tr><tr><td></td><td>a</td><td>b</td><td>d</td><td>e</td><td>f</td></tr><tr><td>a</td><td>1.</td><td>0.390</td><td>-0.622</td><td>0.678</td><td>-0.746</td></tr><tr><td>b</td><td>0.390</td><td>1.</td><td>-0.650</td><td>0.725</td><td>-0.839</td></tr><tr><td>d</td><td>-0.622</td><td>-0.650</td><td>1.</td><td>-0.984</td><td>0.822</td></tr><tr><td>e</td><td>0.678</td><td>0.725</td><td>-0.984</td><td>1.</td><td>-0.900</td></tr><tr><td>f</td><td>-0.746</td><td>-0.839</td><td>0.822</td><td>-0.900</td><td>1.</td></tr></tbody></table>

## Data availability

Data will be made available on request.

## References

* Alotaibi, M.B., Azmy, R., Nasr-El-Din, H.A., 2010. Wettability challenges in carbonate reservoirs. In: SPE Improved Oil Recovery Symposium. https://doi.org/10.2118/129972-ms.

* Anderson, W.G., 1987. Wettability literature survey—part 4: the effects of wettability on capillary pressure. J. Pet. Technol. 39, 1283–1300. https://doi.org/10.2118/15271-PA.

* Arshadi, M., Gesho, M., Qin, T., Goual, L., Piri, M., 2020. Impact of mineralogy and wettability on pore-scale displacement of NAPLs in heterogeneous porous media. J. Contam. Hydrol. 230. https://doi.org/10.1016/j.jconhyd.2020.103599.

* Ayaz, M., Toussaint, R., Schäfer, G., Måløy, K.J., 2020. Gravitational and finite-size effects on pressure saturation curves during drainage. Water Resour. Res. 56 (10). https://doi.org/10.1029/2019WR026279.

* Bartell, F.E., Osterhof, H.J., 1927. Determination of the wettability of a solid by a liquid. Ind. Eng. Chem. 19 (11), 1277–1280.

* Benisch, K., Wang, W., Delfs, J.-O.D., Bauer, S., 2020. The OGS-eclipse code for simulation of coupled multiphase flow and geomechanical processes in the subsurface. Comput. Geosci. 24, 1315–1331. https://doi.org/10.1007/s10596-020-09951-8.

* Berg, S., Ott, H., Klapp, S.A., Schwing, A., Neiteler, R., Brussee, N., Makurat, A., Leu, L., Enzmann, F., Schwarz, J.-O., Kersten, M., Irvine, S., Stampanoni, M., 2013. Real-time 3D imaging of Haines jumps in porous media flow. Proc. Natl. Acad. Sci. 110 (10), 3755–3759. https://doi.org/10.1073/pnas.1221373110.

* Blunt, M.J., 2017. Multiphase Flow in Permeable Media: A Pore-Scale Perspective. Cambridge University Press, Cambridge. https://doi.org/10.1017/9781316145098.

* Blunt, M.J., Alhosani, A., Lin, Q., Scanziani, A., Bijeljic, B., 2021. Determination of contact angles for three-phase flow in porous media using an energy balance. J. Colloid Interface Sci. 582, 283–290. https://doi.org/10.1016/j.jcis.2020.07.152.

* Bohy, M., Dridi, L., Schäfer, G., Razakarisoa, O., 2006. Transport of a mixture of chlorinated solvent vapors in the vadose zone of a sandy aquifer. Vadose Zone J. 5, 539–553.

* Bourgeat, A., Panfilov, M., 1998. Effective two-phase flow through highly heterogeneous porous media: capillary nonequilibrium effects. Comput. Geosci. 2 (3), 191–215.

* Brown, H.W., 1951. Capillary pressure investigation. J. Pet. Technol. 3, 67–74. https://doi.org/10.2118/951067-G.

* Bruce, W.A., Welge, H.J., 1947. Restored state method of determination of oil in place and connate water. Oil Gas J. 46 (12), 223–228.

* Buckley, J.S., 1998. Wetting alteration of solid surfaces by crude oils and their asphaltenes. Rev. Inst. Fr. Pétrol. 53 (3), 303–312.

* Cassie, A.B., Baxter, S., 1944. Wettability of porous surfaces. Trans. Faraday Soc. 40, 546–551.

* Celia, M.A., Bachu, S., Nordbotten, J.M., Bandilla, K.W., 2015. Status of CO₂ storage in deep saline aquifers with emphasis on modeling approaches and practical simulations. Water Resour. Res. 51, 6846–6892. https://doi.org/10.1002/2015WR017609.

* Craig, F.F., 1971. The reservoir engineering aspects of waterflooding. In: Monograph Series SPE, vol. 3. H.L. Doherty Memorial Fund of AIME, New York.

* Dahle, H.K., Celia, M.A., Hassanizadeh, S.M., 2005. Bundle-of-tubes model for calculating dynamic effects in the capillary-pressure-saturation relationship. Transp. Porous Media 58 (1–2), 5–22.

* Desouky, S.E., 2003. A new method for normalization of capillary pressure curves. Oil Gas Sci. Technol. Rev. Inst. Fr. Pétrol. 58 (5), 551–556. https://doi.org/10.2516/ogst:2003038.

* Drellich, J.W., 2019. Contact angles: from past mistakes to new developments through liquid-solid adhesion measurements. Adv. Colloid Interf. Sci. 267, 1–14. https://doi.org/10.1016/j.cis.2019.02.002.

* Dury, O., Fischer, U., Schulin, R., 1998. Dependence of hydraulic and pneumatic characteristics of soils on a dissolved organic compound. J. Contam. Hydrol. 33 (1–2), 39–57. https://doi.org/10.1016/S0169-7722(98)00064-3.

* Fatt, I., Klikoff, W.A., 1959. Effect of fractional wettability on multiphase flow through porous media. J. Pet. Technol. 11, 71–76. https://doi.org/10.2118/1275-G.
Flint, A.L., Flint, L.E., 2002. 2.2 Particle density. In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis, Part 4, Physical Methods. Soil Science Society of America, Madison, WI, USA, pp. 229–240.

Geistlinger, H., Zulfiqar, B., Schlueter, S., Amro, M., 2021. New structural percolation transition in fractional wet 3D-porous media: a comparative µ-CT study. Water Resour. Res. 57. https://doi.org/10.1029/2021WR030037.

Haines, W.B., 1930. Studies in the physical properties of soil. V. The hysteresis effect in capillary properties and the modes of moisture distribution associated therewith. J. Agric. Sci. 20, 97–116. https://doi.org/10.1017/S002185960008864X.

Hassanizadeh, S.M., Gray, W., 1988. High-velocity flow in porous-media by Hassanizadeh and Gray—reply. Transp. Porous Media 3 (3), 319–321.

Hassanizadeh, S.M., Celia, M.A., Dahle, H.K., 2002. Dynamic effect in the capillary pressure—saturation relationship and its impacts on unsaturated flow. Vadose Zone J. 1, 38–57.

Helmig, R., 1997. Multiphase Flow and Transport Processes in the Subsurface, A Contribution to the Modeling of Hydrosystems. Springer (Environmental Engineering), Berlin.

Helmig, R., Flemisch, B., Wolff, M., Ebigbo, A., Class, H., 2013. Model coupling for multiphase flow in porous media 51. Adv. Water Resour. 51, 52–66.

Klute, A., 1986. Water retention: Laboratory methods. In: Klute, A. (Ed.), Methods of Soil Analysis. SSSA Book Series, pp. 635–662. https://doi.org/10.2136/ sssabookser5.1.2ed.c26.

Kutflek, M., Nielsen, D.R., 1994. Soil Hydrology. Catena Verlag, Cremlingen-Destedt, Germany.

Le, V.S., Sharma, K.V., Hanamertani, A.S., Youssif, M.I., Elkhatib, O., Rane, K., Piri, M., Katiyar, A., Nagarajan, N., 2023. Methane foam performance in oil-wet unconsolidated porous media: a systematic experimental investigation at reservoir conditions. Fuel 344. https://doi.org/10.1016/j.fuel.2023.128002.

Lenormand, R., Zarcone, C., Sarr, A., 1983. Mechanisms of the displacement of one fluid by another in a network of capillary ducts. J. Fluid Mech. 135, 337–353. https://doi. org/10.1017/S00221120830003110.

Leverett, M.C., 1941. Capillary behavior in porous solids. Trans. AIME 142, 152–169. https://doi.org/10.2118/941152-G.

Liu, H., Ju, Z., Bachmann, J., Horton, R., Ren, T., 2012. Moisture-dependent wettability of artificial hydrophobic soils and its relevance for soil water desorption curves. Soil Sci. Soc. Am. J. 76, 342–349. https://doi.org/10.2136/sssaaj11.0081.

Løvoll, G., Jankov, M., Måløy, K.J., Toussaint, R., Schmittbuhl, J., Schäfer, G., Méheust, Y., 2011. Influence of viscous fingering on dynamic saturation-pressure curves in porous media. Transp. Porous Media 86, 305–324. https://doi.org/ 10.1007/s11242-010-9622-8.

Morrow, N.R., 1970. Irreducible wetting-phase saturations in porous media. Chem. Eng. Sci. 25 (11), 1799–1815. https://doi.org/10.1016/0009-2509(70)80070-7.

Morrow, N.R., 1976. Capillary pressure correlations for uniformly wetted porous media. J. Can. Pet. Technol. 15. https://doi.org/10.2118/76-04-05.

Morrow, N.R., 1990. Wettability and its effect on oil recovery. J. Pet. Technol. 42 (12), 1476–1484.

Motealleh, S., Ashouripashaki, M., DiCarlo, D., Bryant, S., 2013. Unified model of drainage and imbibition in 3D fractionally wet porous media. Transp. Porous Media 99, 581–611. https://doi.org/10.1007/s11242-013-0201-7.

Moura, M., Flekkoy, E.G., Måløy, K.J., Schäfer, G., Toussaint, R., 2019. Connectivity enhancement due to film flow in porous media. Phys. Rev. Fluids 4, 094102. https:// doi.org/10.1103/PhysRevFluids.4.094102.

Mualem, Y., 1984. A modified depend-domain theory of hysteresis. Soil Sci. 137, 283–291.

Nouri, M., Homaee, M., Bybordi, M., 2014. Quantitative assessment of LNAPL retention in soil porous media. Soil Sediment Contam. Int. J. 23 (8), 801–819. https://doi.org/ 10.1080/15320383.2014.887650.

Parker, J.C., Lenhard, R.J., 1987. Model for hysteretic constitutive relations governing multiphase flow 1. Saturation- pressure relations. Water Resour. Res. 23, 2187–2196. https://doi.org/10.1029/WR023i012p02187.

Parker, J.C., Lenhard, R.J., Kuppusami, T., 1987. A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res. 23, 618–624. https://doi.org/10.1029/WR023I004P00618.

Poulavissilis, A., 1962. Hysteresis of pore water, an application of concept of independent domains. Soil Sci. 93, 405–412.

Raeesi, B., Morrow, N.R., Mason, G., 2014. Capillary pressure hysteresis behavior of three sandstones measured with a multistep outflow-inflow apparatus. Vadose Zone J. 13 (3), 1–12. https://doi.org/10.2136/vzj2013.06.0097.

Rose, W., Bruce, W.A., 1949. Evaluation of capillary characters in petroleum reservoir rocks. J. Pet. Technol. 1 (5), 127–142. https://doi.org/10.2118/949127-G.

Saraji, S., Goual, L., Piri, M., Plancher, H., 2013. Wettability of supercritical carbon dioxide/water/quartz systems: simultaneous measurement of contact angle and interfacial tension at reservoir conditions. Langmuir 29, 6856–6866. https://doi. org/10.1021/la3050863.

Schäfer, G., di Chiara Roupert, R., Alizadeh, A.H., Piri, M., 2020. On the prediction of three-phase relative permeabilities using two-phase constitutive relationships. Adv. Water Resour. 145. https://doi.org/10.1016/j.advwatres.2020.103731.

Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9 (7), 671–675. https://doi.org/10.1038/nmeth.2089.

Schön, J.H., 2015. Physical properties of rocks. In: Fundamentals and Principles of Petrophysics. Elsevier. ISBN: 978-0-08-100404-3.

Siddiqui, M.A.-Q., Ali, S., Fei, H., Roshan, H., 2018. Current understanding of shale wettability: a review on contact angle measurements. Earth Sci. Rev. 181, 1–11. https://doi.org/10.1016/j.earscirev.2018.04.002.

Sinnokrot, A.A., Ramey Jr., H.J., Marsden Jr., S., 1971. Effect of temperature level upon capillary pressure curves. Soc. Petrol. Eng. J. 11 (1), 13–22. https://doi.org/ 10.2118/2517-PA.

Stauffer, F., Kinzelbach, W., 2001. Cyclic hysteretic flow in porous medium column: model, experiment, and simulations. J. Hydrol. 240, 264–275.

Ustohal, P., Stauffer, F., Dracos, T., 1998. Measurement and modeling of hydraulic characteristics of unsaturated porous media with mixed wettability. J. Contam. Hydrol. 33 (1-2), 5–37. https://doi.org/10.1016/S0169-7722(98)00063-1.

Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.

Wenzel, R.N., 1949. Surface roughness and contact angle. J. Phys. Chem. 53, 1466–1467. https://doi.org/10.1021/j150474a015.

Yang, H., Baudet, B.A., 2016. Characterisation of the roughness of sand particles. Procedia Eng. 158, 98–103. https://doi.org/10.1016/j.proeng.2016.08.412.

Yang, H., Baudet, B.A., Yao, T., 2016. Characterization of the surface roughness of sand particles using an advanced fractal approach. Proc. R. Soc. A 472. https://doi.org/ 10.1098/rspa.2016.0524.

Young, T., 1805. An essay on the cohesion of fluids. Philos. Trans. R. Soc. Lond. 95, 65–87. https://doi.org/10.1098/rstl.1805.0005.

Youssif, M.I., Sharma, K.V., Goual, L., Piri, M., 2024. Experimental evaluation of foam-assisted gas injection in proppant-packed fractured oil-wet carbonate. Energy Fuel 38 (4). https://doi.org/10.1021/acs.energyfuels.3c04322.

Zhou, D., Blunt, M., 1997. Effect of spreading coefficient on the distribution of light nonaqueous phase liquid in the subsurface. J. Contam. Hydrol. 25 (1–2), 1–19. https://doi.org/10.1016/S0169-7722(96)00025-3.