HYPROP measurements of the unsaturated hydraulic properties of a carbonate rock sample

The unsaturated hydraulic properties provide important theoretical and practical information about fluid flow in soils and rocks for a range of soil, environmental and engineering applications. In this study we used the evaporation (HYPROP) and chilled-mirror dew point (WP4C) methods to estimate the water retention and unsaturated hydraulic conductivity curves of an Indiana Limestone carbonate rock sample. The obtained data were analyzed in terms of unimodal and bimodal VG (van Genuchten, 1980) and PDI (Peters et al., 2015) type functions. Bimodal functions were found to produce excellent descriptions for the unsaturated hydraulic data we measured, slightly better than the standard unimodal formulations. For our particular test using an Indiana Limestone rock sample, we did not find much improvement when accounting for film and corner flow using the PDI formulation, relative to the VG model. As far as we know, this is the first time the HYPROP methodology was applied to a rock sample.


Introduction
The hydraulic properties of unsaturated porous media provide important information about single-and multi-phase fluid flow in the subsurface. They are critical to addressing many theoretical and practical studies in the soil, hydrogeologic, agricultural, civil and petroleum engineering disciplines. A large number of experimental methodologies has been developed and tested over the years to estimate the water retention, θ(h), and unsaturated hydraulic conductivity, K(h) or K(θ), relationships, where θ is the volumetric water content, h is pressure head, and K the hydraulic conductivity. Most of the standard techniques, especially for the hydraulic conductivity, are suitable for intermediate or relatively wet conditions (Dane and Hopmans, 2002;Durner and Lipsius, 2005;Looney and Falta, 2000), although some can also be applied to the dry range, such as hot-air, centrifugation, or dew-point methodologies (Arya, 2002;Nimmo et al., 2002;Scanlon et al., 2002;Durner and Lipsius, 2005).
One popular experimental approach has long been the evaporation method, starting with early studies by Gardner and Miklich (1962) and Wind (1968), and many others later (e.g., Wendroth et al. (1993) and Iden and Durner (2008), and references therein). Measurements of the evaporation rate and pressure head at multiple depths in the sample permit then the simultaneous estimation of the water retention and hydraulic conductivity curves, either using direct measurements or inverse procedures (e.g., Simunek et al., 1998). A popular recent version of the evaporation method involving semi-automated direct measurements is the HYPROP system Peters and Durner, 2008;Schindler et al., 2010;among others), commercialized by the METER group (München, Germany). The HYPROP system involves pressure head measurements versus time at two depths within a short 5-cm sample as water evaporates from its surface, with the evaporation rate determined by continuous weighing of the column setup.
Several studies successfully tested the HYPROP evaporation approach against multistep outflow and other methods (e.g., Schelle et al., 2010;Zhuang et al., 2017), including tests against synthetic data generated with the Richards unsaturated flow equation (e.g., Bezerra-Coelho et al., 2018;Iden et al., 2019). Most of these and other studies used van Genuchten-Mualem (VG) type expressions (van Genuchten, 1980;van Genuchten and Nielsen, 1985) for the unsaturated hydraulic functions, while Peters et al. (2015) used the equations of Kosugi (1996) as well as the PDI model of Peters and colleagues (Iden and Durner, 2014;Peters, 2014) to account for the effects of film and corner flow in very dry soils.
Most of the above techniques, including the HYPROP approach, typically work well for unconsolidated media, generally soils. While some of the techniques can be applied also to consolidated rocks, different approaches are more suitable for such media, such as mercury intrusion porosimetry, water adsorption or desorption measurements, centrifuge methods, or core flooding (Purcell, 1949;Dullien, 1979;Liu et al., 2016;Haghi et al., 2020). In this study we extend the HYPROP approach to a carbonate rock sample. The HYPROP results, augmented with WP4C (chilled-mirror dew point) water retention data in the dry range, are analyzed in terms of van Genuchten type hydraulic functions, assuming the presence of single-porosity (unimodal) or double-porosity (bimodal) pore systems. The data are for an Indiana Limestone carbonate rock sample extracted from a quarry in the USA, similar to those found in Brazilian and other oil-bearing pre-salt reservoirs.

Experimental procedures
Here we briefly review the HYPROP measurement system, as well as summarize WP4C and related experimental protocols used for this study, including how the data were analyzed in terms of functional descriptions of the hydraulic properties. Our studies were carried out using two samples of Indiana Limestone, a calcite-cemented grainstone carbonate rock. One standard rock sample was used for the HYPROP measurements, and a much smaller separate sample for the WP4C experiments (taken from the same larger Indiana Limestone core we had). In addition, another standard sample was used to drill out plugs needed to cap the holes drilled for the two tensiometers (see Fig. 1). These plugs also originated from the same larger Indiana Limestone core used for the permeability and routine core analysis. The rock consisted of fossil fragments and concentric lamellar calcium carbonate particles (oolites) (Indiana Limestone Handbook, 1975). Petrographic and X-ray diffraction (XRD) analyses indicated that Indiana Limestone is made up nearly exclusively by calcite (99%), with very small amounts (about 1%) of quartz (Churcher et al., 1991;Drexler et al., 2019).

The HYPROP measurement system
To evaluate the relationship between volumetric water retention and the pressure head (negative capillary pressure, given in cm in this study) of the rock samples, we used the HYPROP evaporation method for the wet range (between saturation and pressure heads of about − 850 cm) and the WP4C psychrometer method for the dry range. Instrumentation for both approaches were obtained from the METER Group AG (München, Germany). Fig. 1A provides a schematic of the HYPROP setup involving pressure head measurements at two depths within a 5-cm long fluid-saturated sample. Water contents and fluid fluxes versus time are determined by automated weighing of the sample as water evaporates from the surface. Observed pressure heads, water contents, and evaporation fluxes are subsequently used to estimate the water retention and unsaturated hydraulic conductivity functions (Pertassek et al., 2015;Schindler et al., 2010). The measurement range of the hydraulic conductivity at the wet side is restricted by limitations of the tensiometers to register very small pressure head differences, while limitations at the dry side are due to water outgassing, bubble formation and subsequent bubble expansion in the tensiometers, usually at about − 850 cm.
Pressure heads versus time within the sample are recorded at two locations (1.25 cm and 3.75 cm from the evaporating surface), while the evaporation rate is obtained by automated weighing. Water retention and hydraulic conductivity data points are subsequently estimated from the average water content, the pressure heads provided by the two tensiometers, and the mean total head gradient between the two tensiometers. The approach assumes that linear vertical distributions are present of the pressure head and the water content within the sample. These assumptions have been shown to be very much acceptable (Bezerra-Coelho et al., 2018;Peters et al., 2015).
The standard HYPROP setup is suitable for cylindrical samples measuring 5 cm in height, and either 5 or 8 cm in diameter. For our experiments we used an Indiana Limestone sample having a diameter of only 3.85 cm. Since the sample's diameter was smaller than the inner diameter of 5 cm, we used a 3D printer to construct an impermeable outside ring. A small rotary saw was used to drill circular holes with a diameter of 0.5 cm vertically through the rock sample to receive the tensiometers and a holder pin at the bottom (Fig. 1B). The tensiometers shafts were placed within the holes such that the middle of the tensiometers cups were located at 1.25 and 3.75 cm from the bottom. Very small, fine silt-sized, crushed rock particles that remained from drilling, were used to ensure tight lateral contact between the tensiometer cups and the sample. The remaining open parts of the two holes above the tensiometers were capped with drilled-out plugs (Fig. 1C) from a neighboring sample of the same overall Indiana Limestone core we used. Wetted ground sample material was first placed around the tensiometer cups using a small syringe to again guarantee close contact between the plugs and the tensiometer cups, as well as laterally with the rock sample. The ground rock particles were compressed as much as possible to avoid creating artificial macropores. Additional details about the standard HYPROP system and the evaporation methodology are provided by Pertassek et al. (2015) and METER (2015). For our experiments we used the HYPROP setup for hydraulic measurements down to pressure heads of about − 850 cm. For measurements in the dry range we used the WP4C chilled-mirror dew point method.

WP4C and other measurements
Psychrometer (WP4C) measurements were performed on two 0.5 cm long Indiana Limestone samples having the same lateral dimensions as used for the HYPROP measurements. Using drippers, the sample were wetted to a water content at or slightly below full saturation to avoid hysteresis effects in the drying curves by starting at relatively wet conditions. Drying was done stepwise using a desiccator. All WP4C measurements were preceded by equilibrating the temperatures of the cup and the samples using a thermal equilibration plate as recommended by METER (München, Germany). WP4C evaluations started at approximately pF = 3.5 (about − 3000 cm pressure head), where pF = − log|h|, with h expressed in cm. The WP4C measurements continued until pF = 5.5 (about − 27,000 cm) when the weight variations could no longer be detected accurately. After the last WP4C measurements the samples were dried in an oven at 105 • C to check the dry weight of the samples.
A DV-4000 Poropermeameter, an automatic steady-state gas permeameter-porosimeter system from Weatherford Laboratories (Houston, USA), was further used to measure the absolute permeability and porosity of the sample. The porosity was measured using Helium (He) gas and the absolute permeability, k, using Nitrogen gas (N 2 ), with corrections for the Klinkenberg effect to account for the slippage of air in the sample (Klinkenberg, 1941). The corrected air permeability was automatically obtained internally by the DV-4000 equipment by measuring the permeability at several pressures and extrapolating to infinite pressure (WL, 2017). Measurements had a reliability range between 0.001 and 40,000 mD (millidarcies) for the absolute permeability. The saturated hydraulic conductivity for water was subsequently obtained using the relation K s = ρ w gk/μ, where k is the permeability, ρ w is the fluid density, g is gravitational acceleration and μ the fluid viscosity.
These basic petrophysical measurements of porosity and the hydraulic conductivity at full saturation were fixed during analysis of the unsaturated data in terms of the hydraulic functions.

Unsaturated hydraulic functions
The measured data pairs of water content (θ) versus pressure head (h) as obtained with the HYPROP and WP4C experiments were fitted jointly using the HYPROP-FIT (Pertassek et al., 2015). Soil hydraulic properties for most scenarios were initially described using the standard van Genuchten-Mualem (VG) formulation given by (van Genuchten,

1980)
where S e is effective saturation, θ s and θ r are the saturated and residual water contents, respectively, K s is the saturated hydraulic conductivity, α and n are semi-empirical shape parameters, m = 1-1/n, and L is a poreconnectivity parameter. We additionally explored the performance of a dual-porosity (bimodal) extension of the hydraulic functions to account for the presence of distinct but interacting macropore and micropore regions. The functions we used are given by Durner (1994) and Priesack and Durner (2006): where S i , α i , n i , and m i (i = 1, 2) are the same as in Eqs.
(1) and (2) for the macropore and micropore regions, respectively (or for the fracture and matrix regions if interpreted for unsaturated fractured rock), while w defines the division of the porous medium in macropore and micropore regions. We further tested the van Genuchten formulation assuming variable m and n values in Eq. (1) as described by van Genuchten and Nielsen (1985), as well as the PDI formulation by Peters et al. (2015), the latter accounting for capillary flow, film and corner flow and vapor flow. Because of their complex nature, we decided not to restate here the complete mathematical formulations of the extensions involving the variable m and n cases and the PDI formulation. Once the pressure heads and actual evaporation rates (estimated from the monitored sample weights) were obtained, the HYPROP-FIT analysis was applied to the measured data. HYPROP-FIT uses the root mean square error (RMSE) to quantify differences between the measured (y i ) and calculated (y i c ) water retention and hydraulic conductivity data: where n p is the number of data points, and y i refers to either water content, θ i , or the logarithm of the unsaturated hydraulic conductivity, i. e., log(K i ). The HYPROP-FIT analysis also provides values of the corrected Akaike Information Criterion (AICc), usually a negative value. The larger the absolute AICc number, the more appropriate the model. Fig. 2 shows the observed HYPROP and WP4C water retention data, and the HYPROP derived hydraulic conductivity data, the latter plotted versus volumetric water content as well as the pF (=− log|h|, with h expressed in cm). The Indiana Limestone retention data reflect a somewhat bimodal pore size distribution, presumably in part due to fine calcite crystals lining the pores and creating microporosity, and in part perhaps due to intra-particle microporosity in some of the fossil fragments and oolites. A study by Churcher et al. (1991) using thin sections and scanning electron microscopy showed that the porosity and permeability of these rocks are controlled mostly by the distribution of coarse pore-filling calcite cement. As compared to the water retention data, the hydraulic conductivity data did not show a similar clear bimodal behavior, mostly because they are outside of the micropore range identified with the WP4C data. The dual-porosity nature of our Indiana Limestone samples is consistent with several previous studies showing the bimodal nature of Indiana Limestone rocks. For example, Mercury Injection Capillary Pressure (MICP) tests by Churcher et al. (1991) on Indiana Limestone samples showed very clear bimodal poresize distributions. An NMR study by Dunn et al. (1994) similarly showed bimodality in the relaxation time and associated pore-size distributions, with the peak of the smaller pores being associated with irreducible saturation.

Results and discussion
We next analyzed the observed water retention and conductivity data in terms of the functional hydraulic descriptions in Section 2.2. Fig. 3 shows results for the water retention curve. Clearly, the traditional unimodal van Genuchten and PDI models (including the van Genuchten model with variable m and n parameters) did not match the data well due to the bimodal nature of the curves, while the bimodal equivalents provided excellent fits to the data. Fitted parameter values of the various models are listed in Table 1, while Table 2 compares the statistical results of the different optimizations. We also included in Table 1 the PDI parameters. While the unimodal PDI model improved slightly upon the unimodal van Genuchten equation, far better results were obtained using the bimodal functions. For these reasons we focus below only on the standard unimodal and bimodal van Genuchten functions assuming m = 1-1/n. The VG and PDI scenarios assuming variable m,n parameters are hence not further shown. Fig. 4 shows results for the standard unimodal and bimodal van Genuchten hydraulic conductivity functions, assuming m = 1-1/n, plotted again versus volumetric water content as well as the pressure head (the latter in terms of pF values). Results for the variable m,n functions were only marginally better, which was the case also for the PDI functions compared to the bimodal VG functions. The negligible differences between the classical VG functions and the PDI modifications within the range of available conductivity data is not surprising, since the expected differences between the VG and PDI formulations should occur primarily in the dry range of the conductivity function (>pF 3), which is not covered by our measurements.
We note here that the various optimizations were carried out assuming fixed values of the saturated water content (θ s ) and the saturated hydraulic conductivity (K s ) of the Indiana Limestone sample. These values were obtained by means of poropermeameter measurements on samples IH2 and IH3. The values obtained for IH2 were 311 mD (or 25.9 cm/d) for air, and 287 mD (or 23.9 cm/d) for water. IH3 presented values of 400 mD (33.3 cm/d) for air, and 372 mD (30.9 cm/ d) for water. IH2 and IH3 further showed helium porosities of 0.185 and 0.201 respectively. For our studies we used the IH3 data.
The relative accuracy of the different formulations is best demonstrated by comparing RMSE values of the fitted water content (RMSE θ ) and the hydraulic conductivity (RMSE logK ) data, as well as AICc values. Results for all scenarios are shown in Table 2. The data indicate indeed very little improvement when adopting variable m,n van Genuchten and PDI models, at least for our Indiana Limestone sample. This was the main reason for us not to further pursue those two approaches in this paper, which are far more complicated numerically since the formulations lead to incomplete beta functions or hypergeometric functions (Dourado Neto et al., 2011;van Genuchten and Nielsen, 1985).   Fig. 4. Observed HYPROP hydraulic conductivity curves fitted with the unimodal and bimodal van Genuchten equations assuming m = 1-1/n.

Conclusions
To the best of our knowledge, this is the first time HYPROP and WP4C were used in conjunction to obtain estimates of the water retention and hydraulic conductivity functions of a rock sample. The tested rock, Indiana Limestone, is a carbonate formation sample showing bimodal porosity behavior as evidenced from both our measured retention data, and confirmed by several previous studies. Goodness of fit values showed better results when the van Genuchten bimodal models were fitted to the data, as reflected visually from the fitted data and also seen from RMSE and AICc statistical analyses. Bimodal models indicate the presence of interacting macropore and micropore parts of the medium, in our case each occupying approximately 50% of the pore space. The HYPROP and WP4C data relating volumetric water content and pressure head could be described well with the bimodal van Genuchten hydraulic functions.

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.