J. Hydrol. Hydromech., 72, 2024, X, X–X
https://doi.org/10.2478/johh-2024-0012© 2024 Luwen Zhuang et al., published by Sciendo. This work is licensed under theCreative Commons Attribution-NonCommercial-NoDerivatives 4.0 International Licence.

Analysis of the unsaturated hydraulic properties of rocks using multiple
laboratory methods
Luwen Zhuang1, Fernanda O. Hoerlle2, Hao Chen1, Elizabeth M. Pontedeiro2,3,Martinus Th. van Genuchten3,4, Paulo Couto2, Chao-Zhong Qin5, Kairong Lin1*
1 Center for Water Resources and Environment, and Guangdong Key Laboratory of Marine Civil Engineering, School of Civil Engineering,

Sun Yat-sen University, Guangzhou 510275, China.
2 Department of Civil Engineering, LRAP, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazil.
3 Department of Earth Sciences, Utrecht University, Utrecht, Netherlands.
4 Department of Nuclear Engineering, LASME, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazil.
5 School of Resources and Safety Engineering, Chongqing University, Chongqing, China.
* Corresponding author. E-mail: linkr@mail.sysu.edu.cn

Abstract: Proper characterization of the unsaturated hydraulic properties in rocks is significant for predicting fluid flow
in soil, hydrogeologic, and petroleum science and engineering problems. In this study, we contributed rigorous analysis of
the unsaturated hydraulic properties of three reservoir rock samples (Berea Sandstone, Guelph Dolomite, and Indiana
Limestone). An improved version of the standard evaporation method (HYPROP) was developed to cater specifically torock samples. The improved HYPROP setup enables measurements of local water pressures within rock samples without
disturbing the upper portion of the samples. The obtained results were compared with those obtained using the conventionalpressure plate method and a state-of-the-art nuclear magnetic resonance (NMR) method. Observed data were analyzed interms of four different unimodal and bimodal hydraulic functions. The HYPROP data were found to be relatively close to
the pressure plate data of two carbonate rocks. The NMR-based data were reasonably consistent with the HYPROP data,with differences likely due in part to the fact that they were obtained using two different 5-cm long plugs taken from thesame core. Heterogeneity along the rock cores from which the samples were taken could be a major reason for the observed
differences, and hence should be considered in reservoir analyses.
Keywords: Hydraulic properties; Pressure plate; HYPROP; WP4C; NMR.
INTRODUCTION
The unsaturated hydraulic properties of rock formations are
crucial to understanding fluid flow and contaminant transport
processes in many vadose zone hydrogeologic, geotechnical, andpetroleum science and engineering problems, among other applications. The capillary pressure-saturation and relative permeability functions are the two most important properties of rocks.
Several laboratory techniques (such as pressure plate, centrifuge,
and mercury intrusion methods) are commonly used to quantify
capillary pressure as a function of fluid saturation of rock samples (Shikhov and Arns, 2015). Pressure plate methods are often
assumed to give accurate data without a saturation gradient along
the samples. Unfortunately, pressure plate methods are generally
very time-consuming, often requiring up to several months. Mercury intrusion and centrifuge methods on the other hand are relatively fast (Turturro et al., 2022). With mercury injection capillary pressure (MICP) methods it is possible to measure a broad
range of pore throat sizes, including the distribution and connectivity of pores, but samples are contaminated by mercury andgenerally cannot be used anymore (Liu et al., 2019; Yao and Liu,2012). Centrifuge methods allow one to estimate the drainage
and imbibition capillary pressure curves by applying a centrifugal acceleration on saturated samples with stepwise increasing
rotation speeds (Albuquerque et al., 2018; Chen and Balcom,2005). A drawback of the methods is inaccurate average saturation estimation at each rotation speed (Caputo and Nimmo,2005), although more dynamic methods have now become alsoavailable (Bui and Mori, 2021; van den Berg et al., 2009).

Various methods are also available to measure the relative
permeability of reservoirs rocks. Generally, they can be divided
into steady-state and unsteady-state flow based methods (Favvas
et al., 2009). Regarding steady-state methods, fluid saturation,
flow rate, and pressure data are obtained stepwise at equilibrium
conditions after fluid injection at each step. Relative permeability values can then be calculated using the Darcy-Buckingham
law. These methods can be time-consuming because of the needto reach steady state for each mean pressure. By comparison, unsteady-state laboratory methods are relatively fast, leading
potentially to relative permeability and saturation data that canbe obtained in a few hours depending upon sample size andpermeability. However, many uncertainties still influence theaccuracy of the unsteady-state methods, such as non-equilibrium
capillary effects and viscous fingering (Kumahor et al., 2015;
Zhuang et al., 2017a, 2017b). Almost all of these methods exclusively measure pressures within the reservoir without capturing
internal pressures within the rock samples themselves.
We note that reliable hydraulic property measurements areequally important for applications involving unconsolidated
materials, including soils. Many approaches have beendeveloped over the years to measure capillary pressure andrelative permeability data of soils (e.g., Dane and Hopmans,
2002; Schelle et al., 2013; Zhuang et al., 2017a). Among these,
the evaporation method has become especially popular(Schindler, 1980; Šimůnek et al., 1998; Wendroth et al., 1993),
including the HYPROP semi-automated evaporation systemdesigned by the METER group (METER, 2015). HYPROP isdesigned to continuously collect tensiometer readings at two

1

Luwen Zhuang et al.

depths from the evaporation surface as well as corresponding
sample weights during the full evaporation process (Schindler etal., 2010), from which the hydraulic properties can be estimated.
Lipovetsky et al. (2020) first applied the HYPROP methodology
to reservoir rocks; as far as we know no subsequent studies havebeen reported.
For our studies we improved the HYPROP system used by
Lipovetsky et al. (2020) to measure the unsaturated hydraulicproperties of reservoir rock samples without disturbing the upper
portion of the samples. Different from traditional methods, ourapproach can measure local fluid pressure inside samples. To
validate our approach, we furthermore compared the measured
data obtained using the traditional pressure plate method and an
updated methodology using nuclear magnetic resonance (NMR)
data.
MATERIALS AND METHODSRock samples
Three different rock outcrop samples, obtained from Kocurek
Industries, were used for our studies. One was a sandstone

sample (Berea Sandstone, BS) and two were carbonate rocksamples (Guelph Dolomite, GD, and Indiana Limestone, IL).
Prior to the HYPROP experiments, the plugs were cleaned using
methanol and toluene to remove any salt or impurities from thepore system. After cleaning, the samples were dried for theroutine core analyses. The samples for the HYPROP
experiments were 2 inches (5.08 cm) in diameter andapproximately 5 cm in length. The NMR sample had a diameterof 1.5 inches (4.3 cm) and a length of approximately 3 cm.
Petrophysical properties of the samples are shown in Table 1.
Permeability values of the different rock samples wereobtained each from the same but longer rock core. The porositiesof the two Indiana Limestone plugs were nearly identical.
However, it is worth noting that the measured porosity of one
Indiana Limestone plug was greater than that of the other one
taken from the same core. This variation in porosity may suggest
potential heterogeneity within a longer rock core (Manjunath,2022; Manjunath et al., 2023). Microscope images (Vega 3,
TESCAN, Czech Republic) of the selected samples (Guelph
Dolomite and Indiana Limestone) are shown in Fig. 1. Mineralcompositions of the GD and IL samples are listed in Table 2.

Table 1. Porosities and permeabilities of the rock samples used in this study.
Sample
BS1
GD1
GD2
IL1
IL2

Porosity
(%)
19.12
8.01
11.20
16.24
16.48

Permeability
Klinkenberg (mD)
171.6
1.18
4.37
38.6
68.1

Length
(cm)
5.0
4.9
3.02
4.9
3.03

Diameter
(cm)
5.1
5.1
4.3
5.1
4.3

Rock type
Berea Sandstone
Guelph Dolomite
Guelph Dolomite
Indiana Limestone
Indiana Limestone

Experimental
procedure
HYPROP
HYPROPNMR
HYPROP
NMR

Table 2. Mineral composition of the selected rock samples.
SampleCalciteMagnesian Calcite
Guelph Dolomite–
0.5
Indiana Limestone
99.3–
Total mineralogical composition for rock sample in mass %Dolomite
99.3–

Fig. 1. Microscopic images of the Guelph Dolomite (GD) and Indiana Limestone (IL) samples.

2

K-feldspar
0.1–

Quartz0.1
0.7

Analysis of the unsaturated hydraulic properties of rocks using multiple laboratory methods

Evaporation method
The standard HYPROP system (Meter Group AG, München,
Germany) was previously modified by Lipovetsky et al. (2020)
to enable evaporation experiments of rock samples. The standard
HYPROP system is often employed to measure water retention
and hydraulic conductivity of soils. The soil samples are thensaturated and exposed at the top for evaporation. During theprocess, soil water pressure head readings from the tensiometers
at two depths (1.25 and 3.75 cm) and weight losses of thesamples are collected. Geometric means of the pressure heads
and saturations are used to obtain water retention data, while
unsaturated hydraulic conductivities are estimated from theweight loss and pressure head differences between the twotensiometers (Peters and Durner, 2008). Due to the air-entry
value of tensiometer cups, the HYPROP evaporation method isable to obtain retention data in the relatively wet andintermediate range to about 1,000 cm capillary pressure head
(Schindler et al., 2010). Details about the standard HYPROP
system can be found in METER (2015, 2021).
The standard sample sizes of the HYPROP system are 5 or 8
cm in diameter and 5 cm in height, which is not always suitable
for rock samples. In our work, sample rings for the rock samples
were made using a 3-D printer, with the exact inner size being
the same as the diameter of the rock samples. Silicon grease wasapplied to the inner wall of the sample rings to avoid preferential
flow along the wall. In the system used by Lipovetsky et al.
(2020), openings for the tensiometers were drilled completely
through the samples, with two plugs subsequently used to seal
the holes above the tensiometers after they were installed. Their
setup potentially creates preferential flow paths between theplugs and the sample. In our approach we drilled two holes withthe same dimensions as the two tensiometers from the bottom of
the samples only to the correct positions of the two tensiometers(Fig. 2), thereby minimizing the risk of preferential flow throughany gaps. Ground material from drilled-out rock parts was used
to fill gaps between the tensiometers and the rock samples. The
samples were saturated with deaired water from the bottom
overnight, with air allowing to escape from the top of the sample.
We emphasize that our approach hence differed from the oneused by Lipovetsky et al. (2020), who drilled the vertical holes
through the entire rock sample and backfilled the upper parts of
the holes above the tensiometers with drilled-out 0.5-cm
diameter plugs from a neighboring sample of the same rock core.
Our approach prevents uneven or preferential flow of moisture
between the tensiometers and the top boundary of the samples.
Pressure plate method
We conducted pressure plate experiments using the SWC-150
Fredlund SWCC device, developed by GCTS Testing Systems(Arizona, USA). Saturated rock samples were placed in a sealedcontainer for this purpose. The air pressure around the samples
was increased step by step for the primary drying process. The
bottom of the samples was in contact with a saturated ceramicporous plate. At each step, the air pressure was increased to acertain value to allow water to move out from the sample throughthe ceramic plate at the bottom. A burette was connected to thereservoir below the ceramic plate to collect outflowing water andmeasure changes in the water volume. For each pressure step, thewater volume was monitored to ensure that equilibrium wasreached. The initial water content was calculated based on the
difference between the weights of saturated rock samples andoven-dried rock samples. The rock samples and the ceramic
plates were saturated prior to the experiments. The maximum
capillary pressure head h could be extended to about 1.5 × 104
cm (pF = log(h) = 4.2) before desaturating the ceramic plate. All
experiments were conducted in a constant-temperature room at
21±0.5 °C. Distilled and degassed water was used in theexperiments.
Chilled-mirror dewpoint method
The chilled-mirror dewpoint method using WP4C
psychrometers (Meter Group AG, München, Germany) was used
to obtain water retention data in the relatively dry moisture range(Gee et al., 1992). The method measures the sum of the matricand osmotic (negligible in our study) potentials during theequilibration process between liquid-phase water and vapourphase water in a sealed container (Bittelli and Flury, 2009;
Solone et al., 2012). Neighbouring rock samples were cut into
pieces with dimensions of 1 cm in height and 4 cm in diameter,
which is the exact inner size of the WP4C sample container. The
pieces were saturated by evacuation prior to the measurements.
The saturated samples were subsequently put into the WP4C
container with the top surface exposed to evaporation. The
samples were under primary drying conditions, which isconsistent with the measurement conditions of the HYPROP
system. During evaporation, water pressure heads, of thesamples were measured every 15 min, and water contents
gravimetrically, until there were almost no changes in the sample
weights. The WP4C measurements were reliable from aroundpF = 4.0 (about 10,000 cm capillary pressure head) to pF = 6.5(about 3,100,000 cm capillary pressure head), where pF = log(h).
NMR method
Nuclear magnetic resonance (NMR) is now widely used toimage core body structures to predict capillary pressure andrelative permeability curves of reservoir rocks in petroleum
engineering (see e.g., Dias et al., 2023; Hosseinzadeh et al.,2020; Lima et al., 2022), but increasingly also in the soils

Fig. 2. Schematic of the rock sample with drilled holes for the two tensiometers.

3Luwen Zhuang et al.

literature (e.g., Rahmati et al., 2019). The method uses the
response to hydrogen protons in cores. Petrophysical properties
are evaluated based on different decay rates of the NMR signal.
Capillary pressure and relative permeability curves can then be
obtained from the images using an advanced interpretation
method of NMR. NMR measurements were carried out only on
the Guelph Dolomite and Indiana Limestone samples. Thesesamples, taken from the same core as the HYPROP plugs, weresaturated with NaCl (sodium chloride) brine having a salinity of
30,000 ppm. For this, the plugs were submitted to a vacuum
chamber, which was then invaded with brine and subsequently
pressurized to 2,000 psi for 24 hours. The samples were weighted
regularly during saturation to verify the effectiveness of theprocedure. A saturation of 100% is usually expected, but an errorof ±5 % is generally thought to be acceptable (API RP 40, 1998).
The Guelph Dolomite sample had a saturation of 100% and theIndiana Limestone sample 97%.
The samples were wrapped by Teflon tape to keep themsaturated and avoid fluid losses. The NMR measurements weredone in a Geospec 2 MHz device (Oxford Instruments, UK)
using a Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence(Carr and Purcell, 1954). The CPMG consists of a 90°
radiofrequency pulse succeeded by a 180° pulse train lagged at
π/2 of the first pulse (Lima et al., 2020). To acquire the time
distribution curves of the transverse relaxation (T2), the inverse
Laplace transform was used with the LithoMetrix™ softwarefrom Green Instruments Technology (Frederiction, NB,
Canada). The T2 signal is directly proportional to the poresurface-volume ratio, which is a measure of pore size when
the formation is fully saturated with a single and well-behaved
fluid, such as brine or light oil. In this way T2 can be associatedto the pore-size distribution, porosity, and fluid type at thepore space (Coates et al., 1998). Smaller pores are related withshorter relaxation times and larger pores with longer relaxationtime.
In our studies we used NMR techniques to obtain transverse
relaxation times (T2) of the Guelph Dolomite and Indiana
Limestone samples. From the acquired data, we estimated thepartitioned porosity by applying the superficial relativity for eachtype of rock (Chen, 2014; Souza et al., 2013) and converted theT2 values to equivalent pore radii (Hoerlle et al., 2018). The
approach allows one to obtain the pore-size distribution ofsamples with visual heterogeneities in their pore system. Fig. 3

shows NMR derived T2 relaxation time, which were converted topore radii using surface relaxations of 8 for the IL sample and 13
for the GD sample. The pore radii distributions in turn wereconverted to equivalent capillary pressure – saturation curves
using the approach of Lima et al. (2020) based on the YoungLaplace equation. The T2 curves we obtained for Guelph Dolomite
and Indiana Limestone (Fig. 3) showed two distinct pore size
distribution regions (meso- and macro-pores), which in terms of
the unsaturated hydraulic properties can be represented by dual
porosity functions. We used the RETC optimization program (vanGenuchten et al., 1991) to estimate the hydraulic parameters ofvarious uni- and bimodal hydraulic models for the two samples.
Hydraulic functions
The water retention and unsaturated hydraulic conductivity
data were analyzed in terms of several hydraulic models. We
mostly used van Genuchten (1980) type models in our analyses.
For unimodal data we used the standard van Genuchten-Mualem(VG) equations given by:
S e (h) =S (h) − S r
1
=
m
1− Sr
1 + (α h ) n 



(

2

(2)

where S and Sr denote water saturation and residual water
saturation, respectively, Se is effective water saturation, h is thecapillary pressure head, Ks is the saturated hydraulicconductivity, and α, n, m (m = 1–1/n), and l are mostly semiempirical fitting parameters. Here the values of Ks were fixed at
the corresponding values of ρμ/g times the permeability shownin Table 1, where ρ denotes water density, μ is viscosity, and g
represents gravitational acceleration. Bimodal functions were
further employed to account for the dual-porosity nature ofseveral samples (Durner, 1994; Priesack and Durner, 2006):
S e (h) =S (h) − S r
w1
w2
=
+
m1
m
n
1− Sr
1 + (α h ) 1 
1 + (α h ) n 2  2
1
2





Fig. 3. T2 relaxation times for the Guelph dolomite (GD) and Indiana limestone (IL) rock samples.

4)m
K ( S e ) = K s S e l 1 − 1 − S e1/ m 



(1)

(3)

Analysis of the unsaturated hydraulic properties of rocks using multiple laboratory methods( w1S1 + w2S2 ) l w1α1 1− (1− Se1/m1 )  + w2α2 1− (1− Se1/m2 )m1

K ( Se ) = K s







( w1α1 + w2α2 ) 2

where αi, ni, and mi (i = 1, 2) are the same parameters as those in
Eqs. (1) and (2) but for the macropore and micropore regions,
respectively, while wi (i = 1, 2) denotes the weight factors of thetwo regions, with w1+w2 = 1. The van Genuchten (VG) model
was additionally extended with the Peters-Durner-Iden
formulation (the VG-PDI model). In the VG-PDI formulation,water inside pores is divided into capillary water and absorbed
water, with capillary water being described using the VG modeland adsorbed water using the PDI approach. The absorbed water
mainly refers to film flow, angular flow, and isothermal vapour
flow in relatively dry regions (Peters et al., 2015). Details about
the VG-PDI formulation are given by Peters et al. (2015) and the
HYPROP user manual (METER, 2015).
The HYPROP-FIT software (METER, 2021) was used toanalyse the measured data. The software used the root meansquare error (RMSE) to evaluate deviations between themeasured and fitted water retention and unsaturated hydraulicconductivity values. The software also provides values of thecorrected Akaike information criterion value (AICc), which isalways a negative value. Large absolute AICc values reflect
relatively close agreement of the model with the data. The RMSE
and AICc values are respectively given by:
RMSE =b1

 b ( y ic − y i )

2

(5)

i =1

and
AICc =

2k − 2 L  2k (k + 1) +b
 b − k −1 

(6)in which b is the number of measured data, k is the numberof parameters in the fitted model, L is logarithmic likelihood,and yi and yic are the measured and calculated values,
respectively.
RESULTS AND DISCUSSION

The measured water retention data obtained using the threeexperimental methods are shown in Fig. 4. The HYPROP
measurements for each sample lasted about 10 days when the



2
m2 



(4)

capillary pressure head reached a maximum of approximately
1000 cm. To also cover the relatively low moisture range, we
complemented the HYPROP data with WP4C measurements
above capillary pressure heads of about 5000 cm. The retention
curve obtained from the HYPROP and WP4C measurements of
the BS sample showed a remarkably steep structure with weak
bimodal features. The steep structure suggests a relatively
uniform pore size distribution within the sample. The HYPROP
and WP4C retention data of Guelph Dolomite distributed
relatively evenly over the entire moisture range (Fig. 4b), while
Indiana Limestone showed more bimodal distribution (Fig. 4c).
The retention data obtained using the pressure plate (PP)method are shown as rectangles in Fig. 4. They took around one
month or more for the drying branches of each sample.
Extremely long equilibrated times were especially needed for themeasurements at low moisture. The pressure plate data of Berea
Sandstone (BS) were shifted to the right compared to the
HYPROP data. One reason could be uneven water distribution
within the rock sample caused by a non-uniform pore size distribution, but imperfect equilibration of the moisture regime
appears more likely (which is often the case with PP measurements). Far fewer deviations were observed for Guelph Dolomite. The deviations between the pressure plate and the WP4C
data around 10000 cm of the Indiana Limestone sample could
also be due to the measured pressure plate data not reaching equilibrium. It took more than one week for each Guelph Dolomitepressure plate measurement at relatively low moisture to reachequilibrium. Because of its apparent shortcomings, we will not further consider the pressure plate data.
The NMR retention data for the GD and IL samples (shownas crosses in Fig. 4) were slightly steeper than the HYPROP and
WP4C data. The NMR and HYPROP measurements wereconducted using different plugs taken from the same rock core.
Slightly different porosity values of the two rock plugs could
cause differences in the retention data (see Table 1). Overall, of
the three rocks, only Indiana Limestone showed a clear bimodal
form of its retention curve.
The combined HYPROP and WP4C retention andunsaturated conductivity data were fitted using four hydraulicfunctions, i.e., the original unimodal van Genuchten-Mualem(VG) model, the bimodal van Genuchten (BiVG) model,

Fig. 4. Water retention data of the three rock samples (BS, GD, IL) obtained using the evaporation method (HYPROP, circles), the chilledmirror dewpoint method (WP4C, triangles), the pressure plate method (PP, rectangles). Also shown as cross marks are retention data basedon NMR measurements of the two carbonate rocks (GD and IL).

5

Luwen Zhuang et al.

the van Genuchten-PDI (VGPDI) model, and the bimodal van
Genuchten-PDI (BiVGPDI) model. The saturated hydraulicconductivities were fixed at their calculated values derived fromthe measured permeabilities shown in Table 1.
The measured and fitted retention curves for Berea Sandstone
are shown in Fig. 5a, with the optimized hydraulic parametervalues and statistical results listed in Table 3. The fitted
unimodal curve VG function deviated from the measurements
over the full moisture range. However, the BiVG, VGPDI, and
BiVGPDI functions closely matched the measurements, with
only minor differences in the relatively dry range. Fig. 5b further
shows the measured and fitted hydraulic conductivity curves for
the Berea Sandstone. Results obtained using the unimodal VG
and VGPDI functions deviated slightly from the available
unsaturated conductivity data. The optimized bimodal BiVG and
BiVGPDI curves almost exactly overlapped the measurements.
As expected, the BiVGPDI function, with more fitting
parameters, described the observations best among the fourmodels, which is seen also from the statistical analyses (Table4). However, the BiVG results were also acceptable with fewer
fitting parameters. The calculations reflected the only slightly

bimodal pore structure of the Berea Sandstone rock sample.
Please note that pFdry, ω and A in Table 4 are VGPDI
parameters (METER, 2021), with pFdry representing the pF
value where relative saturation in the very dry region becomeszero (oven dryness).
Observed and calculated retention and unsaturated hydraulicconductivity data of Guelph Dolomite are shown in Fig. 6. The
water retention results obtained with the four hydraulic functions
were almost identical over the full moisture range. The optimized
hydraulic conductivity models almost overlapped with eachother over the available data range, with excellent descriptionsof the HYPROP data. The resulting fitted parameter values arelisted in Table 5, with Table 6 listing the corresponding accuracy
parameters. The VGPDI model described the measurements
minimally better based on the accuracy analyses, although theother functions also matched the experiments over the available
data range. The deviating line of the hydraulic conductivity nearsaturation obtained with VGPDI is an irregularity and could be
avoided by fixing pFdry at 4.8 as is generally assumed following
Schneider and Goss (2012) among others. In our analyses we
optimized the value of pFdry.

Fig. 5. Measured and fitted water retention (a) and unsaturated hydraulic conductivity (b) curves of Berea Sandstone (BS) using the VG
model (black solid lines), the bimodal VG model (black dashed lines), the VG-PDI model (blue solid lines), and the bimodal VG-PDI model(blue dashed lines).
Table 3. Fitted hydraulic parameter values of the unimodal and bimodal models for Berea Sandstone.
Parameter (unit)α (1/cm)n (–)
Sr (–)
Ks (cm/d)l (–)
w2 (–)α2 (1/cm)
n2 (–)
pFdry
ω (–)
A (–)

VG model0.0077
1.93
0.05
14.63
4.38
–
–
–
–
–
–

BiVG model
0.0077
15.00
0.00
14.63
3.04
0.51
0.0054
1.33–––

VGPDI model
0.0077
12.66
0.42
14.63
1.28
–
–
–
5.17
4.7×10–4–1.70

BiVGPDI model
0.0045
6.22
0.37
14.63
2.70
0.75
0.007915.00
5.84
1.2×10–3–2.90

VGPDI model
0.03630.4581–1230

BiVGPDI model
0.03050.3318–1286

Table 4. Statistical results of the unimodal and bimodal models for Berea Sandstone.
Statistical criterion
RMSES
RMSElogKAICc

6

VG model0.10470.6126–1027

BiVG model
0.03680.3140–1269

Analysis of the unsaturated hydraulic properties of rocks using multiple laboratory methods

Fig. 6. Measured and fitted water retention (a) and unsaturated hydraulic conductivity (b) curves of Guelph Dolomite (GD) using the VG
model (black solid lines), the bimodal VG model (black dashed lines), the VG-PDI model (blue solid lines), and the bimodal VG-PDI model(blue dashed lines).
Table 5. Fitted parameter values for different unimodal and bimodal models of rock sample GD.
Parameter(unit)α (1/cm)n (–)
Sr (–)
Ks (cm/d)l (–)
w2 (–)α2 (1/cm)
n2 (–)
pFdry
ω (–)
A (–)

VG
model
0.0124
1.45
0.00
0.10
–0.69
–
–
–
–
–
–

HYPROP+WP4C
BiVG
VGPDI
model
model
0.01240.0121
7.60
1.55
0.00
0.25
0.10
0.10
–0.92
0.35
0.89
–
0.0112
–
1.34
–
–
6.24
–
8.94×10–4––0.21

BiVGPDImodel
0.0448
1.01
0.50
0.10
1.49
0.89
0.0096
1.97
5.47
1.00×10–4
–1.51

VG
model
0.0071
1.93
0.00
–
–
–
–
–
–
–
–

BiVG
model
0.0014
1.56
0.00
–
–
0.73
0.0067
3.20
–
–
–

NMR
VGPDI
model
0.0067
2.78
0.03
–
–
–
–
–
6.80
–
–

BiVGPDImodel
0.0067
3.21
0.00
–
–
0.27
0.0015
1.56
6.82
–
–

BiVG
model
0.0014
–
–1733

NMR
VGPDI
model
0.0016––1544

BiVGPDImodel
0.0007––1732

Table 6. Statistical analysis for different unimodal and bimodal models of rock sample GD.
Statistical
criterion
RMSES
RMSElogK
AICc

VG
model
0.17130.1765–1869

HYPROP+WP4C
BiVG
VGPDI
model
model
0.17130.17130.14990.1515–1873–1874

BiVGPDImodel
0.17130.1456–1871

The results for Indiana Limestone obtained with the HYPROP
and WP4C techniques are shown in Fig. 7, with the corresponding
parameter values listed in Table 7. For the retention curves, theVG and BiVG models described the entry part reasonably well
but not the lower moisture part. In contrast, the PDI functions
gave better agreement in the low moisture range. Regarding theunsaturated hydraulic conductivity, the results of all hydraulicfunctions were almost identical over the available data range. But
they deviated significantly from each other between saturationand the first data points. Based on the statistical analyses data
(Table 8), the VGPDI and BiVGPDI functions described themeasurements better in the low moisture range.
Figs. 4b and 4c show that the NMR-derived water retention
curves of Guelph Dolomite and Indiana Limestone deviated fromthe HYPROP and WP4C data, with the NMR curves being

VG
model
0.0042––1309

steeper than those shown by the full HYPROP and WP4C data.
Still, the NMR-based retention curves are in the general areaof the measured data and hence provide a useful initial indicationof the hydraulic properties of the rock samples. For completeness,
we also analysed the NMR-based curves in terms of thevarious hydraulic models. The NMR derived and optimized
retention curves of Guelph Dolomite and Indiana Limestone areshown in Fig. 8. The VG model deviated from the NMR curves
at relatively high and low moisture contents, while the othermodels gave almost identical descriptions with generally
excellent agreement with the NMR based curve. Based on theaccuracy analyses provided in Tables 6 and 8, the BiVG and
BiVGPDI models described the NMR curves minimally better
than the other models, with only weak hints of having bimodal
distributions.

7Luwen Zhuang et al.

Fig. 7. Measured and fitted HYPROP and WP4C water retention (a) and unsaturated hydraulic conductivity (b) curves of Indiana Limestone(IL) using the VG model (black solid lines), the bimodal VG model (black dashed lines), the VG-PDI model (blue solid lines), and thebimodal VG-PDI model (blue dashed lines).
Table 7. Fitted parameter values for different unimodal and bimodal models of rock sample IL.
Parameter(unit)α (1/cm)n (–)
Sr (–)
Ks (cm/d)l (–)
w2 (–)α2 (1/cm)
n2 (–)
pFdry
ω (–)
A (–)

VG
model
0.0193
1.30
0.00
3.29
–2.69
–
–
–
–
–
–

HYPROP+WP4C
BiVG
VGPDI
model
model
0.50000.0137
1.64
3.23
0.00
0.69
3.29
3.29
–3.39
5.002
0.98
–
0.0136
–
1.36
–
–
5.01
–
1.62×10–2––1.44

BiVGPDImodel
0.0282
1.05
0.00
3.29
4.23
0.06
0.0093
8.33
5.03
5.43×10–2
–1.50

VG
model
0.0130
1.59
0.00
–
–
–
–
–
–
–
–

BiVG
model
0.0173
3.25
0.00
–
–
0.53
0.0019
2.08
–
–
–

NMR
VGPDI
model
0.0135
1.55
0.03
–
–
–
–
–
5.04
–
–

BiVGPDImodel
0.0169
3.07
0.01
–
–
0.48
0.0016
2.38
6.09
–
–

Fig. 8. Measured and fitted NMR-based water retention curves of Guelph Dolomite (GD) and Indiana Limestone (IL) samples using the VG
model (black solid lines), the bimodal VG model (black dashed lines), the VG-PDI model (blue solid lines), and the bimodal VG-PDI model(blue dashed lines).
Table 8. Statistical analysis for different unimodal and bimodal models of rock sample IL.
Statistical
criterion
RMSES
RMSElogK
AICc

8

VG
model
0.03690.1713–1538

HYPROP+WP4C
BiVG
VGPDI
model
model
0.03450.03060.14200.1386–1559–1592

BiVGPDImodel
0.02190.1417
–1651

VG
model
0.0044
–
–1298

BiVG
model
0.0012
–
–1609

NMR
VGPDI
model
0.0042––1303

BiVGPDImodel
0.0008––1683

Analysis of the unsaturated hydraulic properties of rocks using multiple laboratory methods

CONCLUSIONSIn this work, we measured the unsaturated hydraulicproperties of three reservoir rock samples (Berea Sandstone,
Guelph Dolomite, and Indiana Limestone) using an improved
HYPROP experimental setup. The improved HYPROP wasdeveloped to cater specifically to rock samples. This improved
HYPROP setup enabled measurements of the local waterpressures within the rock samples without disturbing the upper
portion of the samples. The WP4C (chilled-mirror dewpoint)
method was employed to extend the capillary pressure data to thedry region. We compared our results with conventional pressure
plate data and nuclear magnetic resonance (NMR) based data.
We have found the pressure plate data shifted slightly to the rightof the HYPROP data, while no such clear deviation wasobserved for Guelph Dolomite and Indiana Limestone. By
comparison, the NMR retention data were steeper than the
HYPROP and WP4C data. Only a slight bimodal nature wasobserved for Indiana Limestone.
Four hydraulic functions (i.e., the original van Genuchten(VG), the bimodal van Genuchten (BiVG), the van GenuchtenPDI (VGPDI), and the bimodal van Genuchten-PDI (BiVGPDI)
functions), were employed to describe the HYPROP and WP4C
retention and unsaturated hydraulic conductivity data, as well as
the NMR-based retention curve. For the HYPROP and WP4C
measurements, the bimodal models were found to reproduce thedata of Berea Sandstone relatively well. The PDI functions
matched the data of Guelph Dolomite and Indiana Limestone.
All models closely described the NMR measurements. Some ofthe deviations between the pressure plate, HYPROP, WP4C andNMR-based hydraulic properties may have been due toheterogeneity along the core from which the respective plugs
were taken. This possible heterogeneity within the cores mayneed to be recognized in reservoir analyses in the future.
Acknowledgements. The authors would like to acknowledge thefinancial support received from Projects of National Key R&D
Program of China (2022YFC3002903). The first author received
funding from the National Natural Science Foundation of China(Grant Nos. 42007165 and 42377057). This research was carriedout in association with the ongoing R&D project registered as
ANP nº 23020-1, “Caracterização Experimental, Modelagem e
Otimização de Processos de Injeção de Água Alternada Com Gás– WAG-EX Fase II" (UFRJ/Shell Brasil/ANP), sponsored by
Shell Brasil Petróleo Ltda under the ANP R&D levy as “Compromisso de Investimentos com Pesquisa e Desenvolvimento”. This study was financed in part by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior- Brasil (CAPES).
Paulo Couto also acknowledges support by the Brazilian Research Council (CNPq) through Grant 311143/2022-9.
REFERENCES

Albuquerque, M., Eler, F., Camargo, H., Compan, A., Cruz, D.,
Pedreira, C., 2018. Estimation of capillary pressure curves
from centrifuge measurements using inverse methods. In:
Proc. Oil & Gas Expo and Conference, Rio de Janeiro, Brazil.
API RP 40, 1998. American Petroleum Institute (API) RP-40
recommended practices for core analysis. Washington, DC,
USA, 220 p.
Bittelli, M., Flury, M., 2009. Errors in water retention curves
determined with pressure plates. Soil Sci. Soc. Am. J., 73,
1453–1460. https://doi.org/10.2136/sssaj2008.0082
Bui, L.T., Mori, Y., 2021. Pinhole multistep centrifuge outflow
method for estimating unsaturated hydraulic properties withsmall volume soil samples. Water, 2021, 13, 1169.
https://doi.org/10.3390/w13091169
Coates, G.G., Borrello, J.A., McFarland, E.G., Mirowitz, S.A.,
Brown, J.J., 1998. Hepatic T2‐weighted MRI: A prospective
comparison of sequences, including breath‐hold, half‐Fourier
turbo spin echo (HASTE). J. Magn. Reson. Imaging, 8, 3,
642–649. https://doi.org/10.1002/jmri.1880080319
Caputo, M.C., Nimmo, J.R., 2005. Quasi‐steady centrifuge
method for unsaturated hydraulic properties. Water Resour.
Res., 41, 11, W11504. https://doi.org/10.1029/2005WR003957
Carr, H.Y., Purcell, E.M., 1954. Effects of diffusion on free
precession in nuclear magnetic resonance experiments. Phys.
Rev., 94, 3, 630. https://doi.org/10.1103/PhysRev.94.630
Chen, Q., Balcom, B.J., 2005. Measurement of rock-corecapillary pressure curves using a single-speed centrifuge andone-dimensional magnetic-resonance imaging. J. Chem.
Phys., 122, 21, 214720. https://doi.org/10.1063/1.1924547
Chen, P., 2014. Enhanced oil recovery in fractured vuggy
carbonates. PhD Thesis. The University of Texas at Austin,USA.
Dane, J.H., Hopmans, J.W., 2002. Water retention and storage:
Laboratory. In: Dane, J.H., Topp, G. (Eds.): Methods of Soil
Analysis. Part 4. Physical Methods. SSSA Book Ser. 5. Soil
Science Society of America, Madison, WI, USA, pp. 675–
720. https://doi.org/10.2136/sssabookser5.4.c25
Dias, C.H., Eler, F.M., Cordeiro, C., Ramirez, M.G., Soares,
J.A., Nunes, D., Lima, M.C., Couto, P., 2023.Effects of poresize and pore connectivity on trapped gas saturation. J. Hydrol. Hydromech., 71, 1, 11–21. https://doi.org/10.2478/johh2022-0042
Durner, W., 1994. Hydraulic conductivity estimation for soils
with heterogeneous pore structure. Water Resour. Res. 30,
211–223. https://doi.org/10.1029/93WR02676
Favvas, E.P., Sapalidis, A.A., Stefanopoulos, K.L., Romanos,
G.E., Kanellopoulos, N.K., Kargiotis, E.K., Mitropoulos,
A.C., 2009. Characterization of carbonate rocks bycombination of scattering, porosimetry and permeability
techniques. Microporous Mesoporous Mater., 120, 109–114.
https://doi.org/10.1016/j.micromeso.2008.09.015
Gee, G.W., Campbell, M.D., Campbell, G.S., Campbell, J.H.,
1992. Rapid measurement of low soil water potentials using awater activity meter. Soil Sci. Soc. Am. J., 56, 1068–1070.
https://doi.org/10.2136/sssaj1992.03615995005600040010x
Hoerlle, F.O., Rios, E.H., Silva, W.G.A.L., Pontedeiro, E.M.B.D.,
Lima, M.C.O., Corbett, P.W.M., Alves, J.L.D., Couto, P., 2018.
Nuclear magnetic resonance to characterize the pore system ofcoquinas from Morro do Chaves Formation, Sergipe-Alagoas
basin, Brazil. Revista Brasileira de Geofísica, 36, 3, 317–324.
http://dx.doi.org/10.22564/rbgf.v36i3.1960
Hosseinzadeh, S., Kadkhodaie, A., Yarmohammadi, S., 2020. NMR
derived capillary pressure and relative permeability curves as an
aid in rock typing of carbonate reservoirs. J. Pet. Sci. Eng., 184,106593. https://doi.org/10.1016/j.petrol.2019.106593
Kumahor, S.K., de Rooij, G.H., Schlüter, S., Vogel, H.-J., 2015.
Water flow and solute transport in unsaturated sand: A
comprehensive experimental approach. Vadose Zone J., 14,
2. https://doi.org/10.2136/vzj2014.08.0105
Lima, M.C.O., Pontedeiro, E.M., Ramirez, M.G., Boyd, A., vanGenuchten, M.Th., Borghi, L., Couto, P., Raoof, A., 2020.
Petrophysical correlations for permeability of coquinas (carbonate rocks). Transp. Porous Media 135, 287–308.
https://doi.org/10.1007/s11242-020-01474-1
Lima, M.C.O., Pontedeiro, E.M., Ramirez, M.G., Favoreto, J.,
Santos, H.N., van Genuchten, M.Th., Raoof, A., 2022.
Impacts of mineralogy on petrophysical properties. Transp.

9

Luwen Zhuang et al.

Porous Media, 145, 103–125. https://doi.org/10.1007/s11242022-01829-w
Lipovetsky, T., Zhuang, L., Teixeira, W.G., Boyd, A., May
Pontedeiro, E., Moriconi, L., Alves, J.L.D., Couto, P., vanGenuchten, M.Th, 2020. HYPROP measurements of the unsaturated hydraulic properties of a carbonate rock sample. J.
Hydrol.,591,
125706.
https://doi.org/10.1016/j.jhydrol.2020.125706
Liu, Y., Liu, Y., Zhang, Q., Li, C., Feng, Y., Wang, Y., Xue, Y,
Ma, H., 2019. Petrophysical static rock typing for carbonate
reservoirs based on mercury injection capillary pressure
curves using principal component analysis. J. Pet. Sci. Eng.,
181, 106175. https://doi.org/10.1016/j.petrol.2019.06.039
Manjunath, G.L., 2022. Estimating stress percolation patterns inhydraulic fractured Gondwana coal using Raman
spectroscopy. J. Nat. Gas Sci. Eng., 100, 104469.
https://doi.org/10.1016/j.petrol.2019.06.039
Manjunath, G.L., Akono, A. T., Haljasmaa, I., Jha, B., 2023. Roleof CO2 in geomechanical alteration of Morrow Sandstone
across micro–meso scales. Int. J. Rock Mech. Min. Sci., 163,
105311. https://doi.org/10.1016/j.ijrmms.2022.105311
METER, 2015. Operation Manual HYPROP 2. METER GroupAG, Munich, Germany.
METER, 2021. WP4C Manual. METER Group AG, Munich,Germany.
Peters, A., Durner, W., 2008. Simplified evaporation method for
determining soil hydraulic properties. J. Hydrol., 356, 1–2,
147–162. https://doi.org/10.1016/j.jhydrol.2008.04.016
Peters, A., Iden, S.C., Durner, W., 2015. Revisiting thesimplified evaporation method: Identification of hydraulicfunctions considering vapor, film and corner flow. J. Hydrol.,527, 531–542. https://doi.org/10.1016/j.jhydrol.2015.05.020
Priesack, E., Durner, W., 2006. Closed-form expression for themulti-modal unsaturated conductivity function. Vadose Zone
J., 5, 1, 121–124. https://doi.org/10.2136/vzj2005.0066
Rahmati, M., Pohlmeier, A., Abasiyan, S.M.A., Weihermüller,
L., Vereecken, H., 2019. Water retention and pore size distribution of a biopolymeric-amended loam soil. Vadose Zone J.,
18:180205. https://doi.org/10.2136/vzj2018.11.0205
Schelle, H., Heise, L., Janicke, K., Durner, W., 2013. Water
retention characteristics of soils over the whole moisture
range: a comparison of laboratory methods. Eur. J. Soil Sci.,
64, 814–821. https://doi.org/10.1111/ejss.12108
Schindler, U., 1980. Ein Schnellverfahren zur Messung derWasserleitfähigkeit
imteilgesättigten
Boden
an
Stechzylinderproben. Arch. Acker-u. Pflanzenbau u.
Bodenkd. Berlin, 24, 1–7.
Schindler, U., Durner, W., Von Unold, G., Mueller, L., Wieland,R., 2010. The evaporation method: Extending themeasurement range of soil hydraulic properties using the air‐entry pressure of the ceramic cup. J. Plant Nutr. Soil Sci., 173,
563–572. https://doi.org/10.1002/jpln.200900201
Schneider, M., Goss, K.U., 2012. Prediction of the water
sorption isotherm in air dry soils. Geoderma, 170, 64–69.
https://doi.org/10.1016/j.geoderma.2011.10.008

10Shikhov, I., Arns, C.H., 2015. Evaluation of capillary pressure
methods via digital rock simulations. Transp. Porous Media,
107, 623–640. https://doi.org/10.1007/s11242-015-0459-z
Šimůnek, J., van Genuchten, M.Th., Wendroth, O., 1998. Parameter
estimation analysis of the evaporation method for determining
soil hydraulic properties. Soil Sci. Soc. Am. J., 62, 4, 894–905.
https://doi.org/10.2136/sssaj1998.03615995006200040007x
Solone, R., Bittelli, M., Tomei, F., Morari, F., 2012. Errors inwater retention curves determined with pressure plates:
Effects on the soil water balance. J. Hydrol., 470, 65–74.
https://doi.org/10.1016/j.jhydrol.2012.08.017
Souza, A., Carneiro, G., Zielinski, L., Polinski, R., Schwartz, L.,
Hürlimann, M.D., Boyd, A., Rios, E.H., dos Santos, B.C.C.,
Trevizan, W.A., Machado, V.F., de Vasconcellos Azeredo,
R.B., 2013. Permeability prediction improvement using 2D
NWR diffusion-T2 maps. In: Proc. PWLA 54th Annual
Logging Symposium, New Orleans, Louisiana, USA.
Turturro, A.C., Caputo, M.C., Gerke, H.H., 2022. Mercury
intrusion porosimetry and centrifuge methods for extendedrange retention curves of soil and porous rock samples.
Vadose Zone J., 21, 1–11. https://doi.org/10.1002/vzj2.20176
van den Berg, E.H., Perfect, E., Tu, C., Knappe, P.S.K., Leao,
T.P., Donat, R.W., 2009. Unsaturated hydraulic conductivity
measurements with centrifuges: A review. Vadose Zone J., 8,
3, 531–547. https://doi.org/10.2136/vzj2008.0119
van Genuchten, M.Th., 1980. A closed-form equation forpredicting the hydraulic conductivity of unsaturated soils.
Soil
Sci.
Soc.
Am.
J.,
44,
892–898.
https://doi.org/10.2136/sssaj1980.03615995004400050002x
van Genuchten, M.Th., Leij F.J., Yates, S.R., 1991. The RETC
code for quantifying the hydraulic functions of unsaturated
soils. Report No. EPA/600/2-91/065. R.S. Kerr
Environmental Research Laboratory, U.S. Environmental
Protection Agency, Ada, OK. 85. https://www.pcprogress.com/Documents/programs/retc.pdf
Wendroth, O., Ehlers, W., Hopmans, J.W., Kage, H.,
Halbertsma, J., Wösten, J.H.M., 1993. Reevaluation of theevaporation method for determining hydraulic functions in
unsaturated soils. Soil Sci. Soc. Am. J., 57, 6, 1436–1443.
https://doi.org/10.2136/sssaj1993.03615995005700060007x
Yao, Y., Liu, D., 2012. Comparison of low-field NMR andmercury intrusion porosimetry in characterizing pore size
distributionsofcoals.
Fuel,95,
152–158.
https://doi.org/10.1016/j.fuel.2011.12.039
Zhuang, L., Bezerra Coelho, C.R., Hassanizadeh, S.M., vanGenuchten, M.Th., 2017a. Analysis of the hysteretic
hydraulic properties of unsaturated soil. Vadose Zone J., 16,
5. https://doi.org/10.2136/vzj2016.11.0115
Zhuang, L., Hassanizadeh, S.M., Qin, C.-Z., de Waal, A., 2017b.
Experimental investigation of hysteretic dynamic capillarity
effect in unsaturated flow. Water Resour. Res., 53.
https://doi.org/10.1002/2017WR020895
Received 5 February 2024
Accepted 15 May 2024