© F. Bignonnet et al., published by IFP Energies nouvelles, 2016

Introduction

Low permeability caprock characterisation is important for both subsurface storage activities (Carbon Capture & Sequestration, gas, etc.) and unconventional resource production. In most published papers, emphasis is put on the threshold capillary pressure characterisation. However, the available literature data about the associated relative permeabilities is scarce although it also contributes to either the confinement of storage facilities or the production rates of wells. This is mainly due to the complexity of deriving these data in a representative manner from laboratory experiments in a reasonable time frame. This work is a contribution to provide such data from a tight carbonate rock with clay content and to propose an approach to interpret their shapes.

Relative gas permeability measurements can be performed in a reasonable time frame on standard permeable media, such as conventional sandstones [1–3]. However, relative gas permeability measurements on weakly permeable media remain challenging. Recent developments have allowed the study of low permeability materials, such as concrete [4] and tight sandstones [5].

The present study focuses on a tight carbonate rock with clay content, extracted from a caprock facies interval. First, the relative gas permeability is assessed experimentally with respect to partial water saturation. To get a deeper understanding of the porous network properties, two porosity measurement methods are used and a capillary pressure curve is determined (Sect. 1 and 2). Second, an analytical permeability upscaling model is proposed and compared to experimental data in Section 3.

1 Material and methods

1.1 Material

1.1.1 Sample Origin and Macroscopic Description

The studied material comes from a single core extracted at a depth of 450 m. In the study, a confining pressure of 9 MPa has been assumed to be representative of the in-situ hydrostatic pressure. At the macroscopic scale, the material is of the type “mudrock”, very fine grained, mostly carbonated and fairly homogeneous (Fig. 1). Geochemical analysis indicates that the carbonate content is high with 276,000 parts per million in weight (ppm) of Ca, located in the light-gray regions. It also features a minor component of pyrite (1.1%) as well as silicates and clay as deduced from the measured values of 7,030 ppm of Mg, 7,980 ppm of K, 11,250 ppm of Fe and 6,130 ppm of S. EDX analysis indicates that the clay content is located in the dark-grey regions and is most likely illite. The total organic carbon is 0.57%.

1.1.2 Microscopic Description

A petrographic analysis using Scanning Electron Microscopy (SEM) shows at a low magnitude of amplification (Fig. 2a) a microstructure of the type “grain-supported”, that is with carbonated grains (light gray) in contact. The size of the carbonated grains ranges from 1 to 10 μm approximately. Pyrite spheres (in white) and silicate silts (dark grey) smaller than 10 μm are embedded in this matrix. Observation at a finer scale (Fig. 2b) shows that the carbonates (light grey) tend to form clusters or aggregates with low porosity. The porosity is rather located in micro-metric regions with a high concentration of clay particles (dark grey) and inter-particle porosity (black). These high porosity regions are seemingly connected through thin “canals” crossing the carbonate aggregates.

1.1.3 Sample Preparation

Eight vertical cylindrical samples with diameter 37 mm and height 10 mm have been prepared from the main vertical core, so that the axis of the cylinders is assumed perpendicular to the bedding. First, longer cylinders are obtained by coring. Second, the long cylinders are cut into several 10 mm height cylinders with a wire saw. The upper and lower faces of the samples are then polished with sandpaper. The height h and diameter d of each sample are then measured with a calliper with a precision of 0.02 mm. The total volume V of each sample is simply assessed by V = πhd²/4.

1.2 Methods

1.2.1 Partial Water Saturation by Relative Humidity Control

All samples are first oven-dried at 60°C until mass stabilisation is achieved, which occurs after about two weeks. This state is considered as the reference dry state, and the water saturation S_w is assumed null. The mass in the dry state m_dry is measured with a precision of 0.01 g.

To obtain a partial water saturation, the samples are stored in a chamber comprising air and water vapour in equilibrium with a saturated saline solution at 20°C [4]. The relative humidity of the atmosphere in the chamber is imposed by thermo-dynamical equilibrium and depends on the type of saline solution. The relative humidities achieved by this method are reported in Table 1.

Since the chemical activities of liquid water and vapour at equilibrium are equal, the capillary pressure P_cap = P_g − P_l between the gas (here at atmospheric pressure) and the liquid water (which condenses in the pore network of the sample) is related to the relative humidity rh of the atmosphere by Kelvin’s law:$$P_{\mathrm{cap}}=-\frac{{\rho }_l\mathrm{RT}}{M}\ln \left(rh\right)$$(1)

where ρ_l is the density of water, M the molar mass of the water, R the perfect gas constant and T the temperature. For water at 20°C, the capillary pressures corresponding to the imposed relative humidities are reported in Table 1.

Moreover, the capillary pressure is related to the curvature C of the liquid-gas meniscii by Young-Laplace’s law:$$P_{\mathrm{cap}}=\gamma C$$(2)

where γ is the coefficient of surface tension, equal to 0.073 N.m⁻¹ under ambient conditions. If the contact angle between the water and the solid phase is null, the curvature of the liquid-gas interface in a cylinder of radius r is C = 2/r. In the case of a meniscus in-between two planes spaced by an aperture e (the case of a “flat pore”), the curvature is C = 2/e. For water at 20°C, the cylinder radii and flat pore apertures corresponding to the imposed relative humidities through the Kelvin-Laplace law are reported in Table 1.

A sample successively equilibrated in chambers with increasing relative humidity undergoes a progressive imbibition, from the smaller to the larger pores. At each step, equilibrium is assumed once mass stabilisation of the sample is achieved, which occurs after about two weeks. The mass m(rh) in the partially water saturated state due to the relative humidity rh is measured with a precision of 0.01 g. Note that the water saturation is not directly known at this stage.

At the end of the progressive imbibition process by successive equilibrium in chambers with increasing relative humidities up to 100%, the samples are immersed in water (partially immersed for the 6 first hours to allow the air to drain, then fully immersed) in a vacuumed chamber for 24 hours. The mass in this saturated state m_sat is measured.

Assuming the final immersion step guarantees a total water saturation of the sample, a way to determine the porosity ϕ is:$$\varphi =\frac{m_{\mathrm{sat}}-m_{\mathrm{dry}}}{{\mathrm{V\rho }}_l}$$(3)

The relation (3) provides a preliminary way of measuring the porosity. Porosity measured using this method will be referred to as “water porosity”. This method will be compared to another method later in the paper.

Additionally, the partial water saturation S_w (rh) achieved at any of the previous steps by equilibrium with relative humidity rh may now be assessed using:$$S_w\left(rh\right)=\frac{m\left(rh\right)-m_{\mathrm{dry}}}{m_{\mathrm{sat}}-m_{\mathrm{dry}}}$$(4)

The S_w − rh curve obtained from Equation (4) can be used to derive the capillary pressure curve by means of Kelvin’s law (1). The latter may be used with caution to estimate the pore size distribution through Young-Laplace’s law (2), provided strong assumptions on the shape (cylindrical or flat pores for example) and connectivity of the pores are made.

1.2.2 Gas Permeability Measurements

Steady State Apparatus

Gas permeability K is measured using a uniaxial steady state gas flow apparatus (Fig. 3). The apparatus consists of a confining cell which is allowed to reach confining pressures P_c as high as 60 MPa, together with a gas injection device [1, 4–9]. The gas used in the current study is 99%-pure Argon with a viscosity of μ = 2.2 × 10⁻⁵ Pa.s.

Let x = 0 be the coordinate of the injection (or upstream) side of the sample and x = h the drainage (or downstream) side. The injection pressure P_i at x = 0 is set to values ranging from 0.5 MPa to 1.5 MPa, whereas the back-pressure P_d at x = h is equal to the atmospheric pressure. The injection pressure is measured with a pressure sensor with a precision of 0.002 MPa. The drainage flow rate q_d at x = h is measured with a set of flow-meters, the finer of them allowing measurement of flow rates from 0.1 mL/min to 1 mL/min with a precision of 0.01 mL/min.

No-Slip Flow

In the case where gas-slippage effects can be neglected, combining Darcy’s law and the perfect gas law to mass conservation yields the uniaxial diffusivity equation describing the pressure p(x, t) in the sample:$$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\varphi p(x,t)\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(p(x,t)\frac{K}{\mu }\frac{\mathrm{\partial }p(x,t)}{\mathrm{\partial }x}\right)$$(5)

In a steady state test, boundary conditions constrain the pressure profile inside the sample to be:$$\begin{array}{cc}p\left(x\right)=\sqrt{P_i^2-(P_i^2-P_d^2)\frac{x}{h}}& (0\le x\le h)\end{array}$$(6)

so that the pressure gradient is . Combining Darcy’s law and flow continuity at x = h yields the expression for the permeability K as a function of the measured drainage flow rate q_d:$$K=\frac{2\mu h{P}_d{q}_d}{A(P_i^2-P_d^2)}$$(7)

where A = πd²/4 is the area of the sample section.

Slip Flow

In the case where the local permeability is governed by Klinkenberg’s slippage law [10], the apparent local permeability depends on the gas pressure as:$$K^{\mathrm{app}}\left(p\left(x\right)\right)=K^{\mathrm{int}}\left(1+\frac{\beta }{p\left(x\right)}\right)$$(8)

where K^int is the intrinsic gas permeability and β is Klinkenberg’s coefficient. Carefully introducing (8) into (5) and considerations similar to the no-slip case rigorously leads to:$$K_m^{\mathrm{app}}=K^{\mathrm{int}}\left(1+\frac{\beta }{P_m}\right)=\frac{2\mu h{P}_d{q}_d}{A(P_i^2-P_d^2)}\hspace{1em}\mathrm{with}\hspace{1em}{P}_m=\frac{P_i+P_d}{2}$$(9)

where P_m is a mean pressure and is a mean apparent permeability, in an averaged way along the sample. The similarity of the relations (9) with (8) and (7) is striking and convenient, but was non-obvious a priori and is merely due to the particular form of the pressure dependence of the permeability in Klinkenberg’s law. Hence, the plot of with respect to 1/P_m allows determination of the intrinsic gas permeability and the Klinkenberg coefficient, irrespective of the fact that P_i and P_d may be very different. In practice, three different gas injection pressures are successively applied: 0.5, 1 and 1.5 MPa.

Effective Gas Permeability

In order to assess the relative gas permeability as a function of the water saturation S_w, effective gas permeability measurements are performed at each step of saturation obtained by the procedure described in Section 1.2.1. As the Argon gas injected in the partially water saturated sample is dry, the mass of the sample is checked before and after the measurement to monitor potential desaturation of the sample. In practice for the current study, changes in water saturation due to dry gas injection have been observed to be negligible.

1.2.3 Gas Porosity Measurements

The porosity has been measured at three steps of loading/unloading of the dry gas permeability measurements using a dedicated method [6], similar in principle to gas pycnometry.

The measurement is performed in the same confining cell as for the permeability test (Fig. 4). The room temperature is set to 20°C. Three volumes have to be distinguished:

• a reference tank with a calibrated volume V_t of 70 mL,

• the pore volume ϕ × V to be measured,

• the volume V_c of the tubes connecting the reference tank to both sides of the sample.

The reference tank is initially at a pressure P_t of 1 MPa which is monitored using a pressure sensor with a precision of 100 Pa, whereas the porosity and the connecting tubes are at atmospheric pressure P_atm. The valves between the reference tank and the sample are then opened while the pressure in the reference tank is monitored. Once gas pressure and temperature are equilibrated, application of the law of perfect gases in isothermal conditions links the final pressure P_f to the sample porosity through:$$P_t{V}_t+P_{\mathrm{atm}}(V_c+\varphi V)=P_f(V_t+V_c+\varphi V)$$(10)

The measure of porosity deduced from Equation (10) is referred to as “gas porosity”. The volume V_c (≈ 5 mL) is measured based on the same principle using a non-porous steel sample.

2 Experimental Results

2.1 Porosity

The measured “water porosities” (Eq. 3) and “gas porosities” (Eq. 10) are reported in Table 2. Water porosities are measured without confining pressure, whereas gas porosities are measured at three steps of a confining pressure cycle (2 MPa, 9 MPa (≈ in-situ pressure) and 2 MPa after unloading).

Gas porosities exhibit a low dispersion and range from 11.7 to 12.7% at 2 MPa of confining pressure. Water porosities are more dispersed, mostly because of the presence of the two specimens MIL-20 and MIL-21 which have been sampled closer to the upper part of the received core than all other samples.

The water porosities appear to systematically underestimate the gas porosities by 1 to 2%. As water molecule and Argon atom diameters are respectively 0.28 nm and 0.38 nm, their size cannot explain this discrepancy. We rather believe that the water saturation method by immersion did not guarantee a total water saturation because of trapped gas. Another source of difference may be the sorption of Argon on pore walls during the porosity measurement experiment, which has not been quantitatively assessed.

The confining pressure increase from 2 MPa to 9 MPa leads to a decrease of 1% in porosity. After unloading, the initial porosity values are retrieved up to the measurement errors. This indicates that the evolution of the samples along the imposed stress path is mostly elastic.

2.2 Capillary Curve

In order to determine the length scales associated with the porosity, the water sorption and desorption curves are measured based on the method described in Section 1.2.1. The S_w − rh curves obtained from relation (4) have been used to derive by means of Kelvin’s law (1) the capillary pressure curves presented in Figure 5.

The desorption curve, which has been measured on a single sample, is close to the sorption curves, tending slightly towards higher water saturations. As for the set of eight sorption curves, the dispersion increases with decreasing capillary pressure (S_w = 0.68 ± 0.10 at P_cap = 2.7 MPa; S_w = 0.115 ± 0.015 at P_cap = 114 MPa). Again, the two specimens MIL-20 and MIL-21, which have the highest total porosities, have also the weakest water saturation at P_cap = 2.7 MPa. This indicates that the extra porosity of these two specimens is mostly a macro-porosity.

The capillary pressure curves presented in Figure 5 are re-interpreted as cumulative pore size distributions in Figure 6 using Young-Laplace’s law (2). To do so, the pores have been assumed flat, i.e. described as the space between two parallel planes with an aperture e, which is of course a rough approximation.

Under these assumptions, around 70% ± 5% of the porosity corresponds to pores with an aperture below 50 nm. This fine porosity is attributed to the clay phase, which is located either in the canaliculi spanning carbonate aggregates or in the high-porosity regions (Fig. 2). The pore size distribution above 50 nm cannot be evaluated by the present method, but mainly corresponds to the inter-particular pores observed in the high-porosity regions.

Finally, in the studied range of pore sizes, a log-uniform repartition law provides an appropriate fit of the desorption curve. The minimum and maximum flat pore apertures of the fitted log-uniform law are respectively 0.6 nm and 370 nm. As these extremal values are out of the range of validity of the present method (which is ≈ 2 nm to 50 nm), they should be considered with care.

2.3 Dry Gas Permeability

The gas permeability has been initially measured on the dried samples using the steady state method presented in Section 1.2.2. Each sample undergoes the following confining pressure cycle: 2 MPa, 9 MPa and 2 MPa. At each confining pressure step, the apparent gas permeability is measured for three different gas injection pressures: 0.5 MPa, 1 MPa and 1.5 MPa. Except for the two specimens MIL-20 and MIL-21, the apparent gas permeability for an injection pressure of 1 MPa is 1.25 ± 0.15 × 10⁻¹⁸ m² and Klinkenberg corrected intrinsic permeability is 2.4 ± 1 × 10⁻¹⁹ m² (Tab. 3). The specimens MIL-20 and MIL-21, which also have more macro-porosity, have permeabilities 2 to 8 times greater than all other samples.

Reported Klinkenberg coefficients (except for specimen MIL-20) correspond to flat pore apertures in the order of 15 to 50 nm [11], according to a slip-flow model based on the kinetic theory of gases. These values are in excellent agreement with the estimated cumulative pore size distributions presented in Figure 6.

The moderate increase in confining pressure from 2 to 9 MPa results in a decrease in permeability of only 10%. After unloading to 2 MPa, the initial permeability is almost retrieved with less than 5% relative difference. As it is well known that small changes of the pore network can lead to significant changes in permeability, these relatively low fluctuations indicate that the deformation of the samples is mostly elastic during the present load cycle.

2.4 Relative Gas Permeability

In order to assess the relative gas permeability, effective gas permeability measurements are carried out on the partially water saturated samples. At each imposed relative humidity step, as in the dry case, the samples undergo a confining pressure cycle from 2 MPa to 9 MPa and back to 2 MPa. At each confining pressure step, the apparent effective gas permeability is measured for three different gas injection pressures: 0.5 MPa, 1 MPa and 1.5 MPa. For water saturations above 0.5, Klinkenberg’s law appeared inappropriate to model the effect of gas pressure. As a result, we choose to present relative gas permeabilities which correspond to apparent permeabilities (computed using Eq. 7) for an injection pressure of 1 MPa. The measured relative gas permeabilities are reported as a function of water saturation in Figure 7, for confining pressures of 2 and 9 MPa.

For water saturations below 0.5, effective gas permeabilities exhibited limited hysteresis to a cycle in confining pressure (2-9-2 MPa). However at the highest water saturation tested (corresponding to rh = 98%), the effective gas permeability after the confining pressure cycle ranged only from 1/2 to 1/3 of the value before the cycle, for the same reference value (2 MPa) of confining pressure. Two hypotheses may be proposed to explain this hysteresis:

• pore compression due to hydrostatic loading induces a spatial reorganisation of the water or,

• the saturation cycle has induced micro-cracking of the samples.

Subsequent permeability measurements after a second drying of the samples suggest that limited micro-cracking might indeed have occurred.

Excepted for the two outlier specimens MIL-20 and MIL-21 which suffered from severe micro-cracking induced by the saturation process, all the measured relative permeabilities in Figure 7 lie on a master curve, for both confining pressures (2 and 9 MPa). The main feature of the relative permeability curve is the presence of a low critical water saturation S_w,c around 0.71 ± 0.04 for which gas permeability nearly vanishes. Such so-called hydraulic cut-off has also been observed on tight sandstones [5].

The influence of the saturation history on the relative permeability curve has not been assessed and requires further investigation. This effect has been studied on two industrial concretes in [4], comparing desorption and re-saturation paths. It has been observed for these concretes that the saturation history influences relative gas permeability mainly in the case of a strong hysteresis in desorption/re-saturation isotherms (i.e. S_w − rh relations). The higher relative gas permeability for the re-saturation path has been attributed to bottle neck effects within the pore network or capillary pressure effects. However, as the imbibition and drainage capillary curves (Fig. 5) do not significantly differ, one may expect little influence of the saturation history on the relative gas permeability for the material presently under consideration.

3 Multi-Scale Modeling

In this section, a method to estimate the relative gas permeability of the present tight carbonate rock with clay content is proposed. We assume the existence of a Representative Volume Element (RVE) for the permeability of this porous material. Two distinct scales are hence considered:

• the macroscopic scale, which is supposed to correspond to the scale at which the previous experiments have been carried out;

• the microscopic scale at which one may distinguish elementary particles such as carbonate aggregates, clay particles and pores.

The upscaling is performed within the framework of random media or so-called continuum micromechanics [12], widely used to derive effective mechanical properties. This strategy has been applied to determine the effective permeability of specific porous media such as fractured media [13–17], argillite [18] or tight sandstones with inter-granular hydraulic joints [19]. Due to the disordered or polycrystalline-like arrangement of the elementary particles at the pore scale, the self-consistent homogenisation scheme initially proposed by [20] in linear elasticity has been extensively used to estimate the macroscopic permeability tensor K^hom. The self-consistent scheme naturally exhibits a so-called percolation threshold in the case of a strong contrast in permeability of the elementary particles. We wish to take advantage of this feature to model the hydraulic cut-off which has been observed experimentally at the critical water saturation.

3.1 Description of the Microscopic Scale

A direct way to upscale the permeability would be to resort to numerical homogenisation on explicit representations of the microstructure, either obtained by 3D imaging techniques or computer synthesised [21–28]. We instead turn to analytical modeling with an implicit description of the microstructure.

Our first modeling step lies in the extraction of the most salient properties of the organisation of the microstructure. We purposely wish to avoid a too exhaustive description as it would lead to a model with too many parameters, thus of limited practical use. From a hydraulic point of view, the pore structure may be simplified as follows (Fig. 2): most of the pores are concentrated in high-porosity regions together with clay particles. These regions of micro-metric extent (≈2 μm) lie in between larger carbonate aggregates or silts (10 to 20 μm), which are assumed too dense to be significantly permeable. The key point is that the connection between the high-porosity regions is only ensured by canaliculi (resp. joints) at the outskirt of the large carbonate aggregates (resp. silts). In other words, the flow is assumed to be mostly controlled by these canaliculi and joints, while the high-porosity regions play a minor role.

A RVE Ω of volume |Ω| of this material is now considered. The high-porosity regions occupy the domain Ω_hp with volume fraction f = |Ω_hp|/|Ω|. To simplify the discussion, carbonate aggregates surrounded by canaliculi and silts with hydraulic joints will be generically called grains with (hydraulic) interfaces. Each interface, which may contain several clay particles, is assumed to have an equivalent aperture e for the fluid flow. To describe the grain-interface phases (Fig. 8), the composite elementary particle of interest is a grain surrounded by an interface (more precisely, just a half-interface as will be explained later). Let G_i(e) be the domain occupied by a grain i, such that its surrounding hydraulic interface has an actual aperture e. The boundary S_i(e) of each grain i corresponds to a half interface. The total domain occupied by all the grains is denoted Ω_g = ∪_iG_i. The volume fraction of all the grains is then 1 − f. As the interfaces are very thin as compared to the radius of the grains, their contribution to the pore volume is neglected.

The apertures e of the interfaces are not all equal; instead, they are assumed to follow a distribution α(e) in the range [e_min, e_max] where e_min and e_max are respectively the smallest and the largest interface apertures. The distribution is normalised such that . The description of the pore space based on size distribution has been classically used in statistical modes based on bundles of cylindrical pores [29–32].

In the case of a partial water saturation S_w resulting from an imposed capillary pressure P_cap, we denote e^⋆ the interface aperture which corresponds to P_capvia Young-Laplace’s law. All the interfaces with apertures e < e^⋆ are assumed fully water saturated, while all interfaces with apertures e > e^⋆ are assumed fully gas saturated. It is then convenient to introduce the fraction s ∈ [0, 1] of water saturated interfaces, which is related to the interface aperture distribution by:$$s={\int }_{e_{\min }}^{e^{\star }}\mathrm{}\alpha \left(e\right)\mathrm{d}e$$(11)

3.2 Hydraulic Properties of the Phases

In the framework of periodic homogenisation, [33] has shown that the steady state, laminar, single phase flow of an incompressible Newtonian viscous fluid in a porous medium obeys Darcy’s law at the macroscopic scale. However, the direct resolution of the Stokes equations on a true RVE of the microstructure is a difficult task. First, the very disperse pore sizes involved and the complexity of their connections make it almost impossible to obtain a well defined, truly representative element of the pore network of an actual rock such as considered in this paper – despite all the advanced imaging techniques currently available. Further, in the case where such RVE were available, it would be far too computationally intensive to solve Stokes equations with a reasonable resolution.

This has motivated authors to propose simplifications of the problem by replacing the different constituents of the micro-scale such as micro-cracks [13–17], inter-granular joints [19] or flat pores [18] by fictitious Darcy media. The simplified problem then amounts to the homogenisation of a heterogeneous Darcy medium. Based on these considerations, let us propose the following fictitious equivalent permeabilities of the elementary particles of each phase.

3.2.1 Dry Case

To start with, let us consider the case where all pores are gas saturated. In what follows, we systematically consider a gas viscosity μ = 1, with no loss of generality.

High-Porosity Regions

At the micro-metric scale, an isotropic apparent permeability K_hp = K_hp1 (where 1 is the second order identity tensor) is attributed to the high-porosity regions which occupy the domain Ω_hp. It is assumed that this permeability could be obtained from the upscaling of a Stokes flow on a representative volume of the high-porosity region from a lower scale. Dimensional analysis indicates that K_hp scales as the square of a characteristic pore size. As will be seen later, further dimensional analysis will indicate that this scaling for K_hp is much greater than the effective permeability K^hom of the RVE, so that exact value of K_hp is not needed in a first approximation. Physically speaking, this corresponds to the case where the pressure drop in the high-porosity regions can be neglected as compared to the pressure drop along the interfaces.

Grains with Interfaces

As proposed for tight sandstones in [19], a two dimensional mathematical representation of the interfaces will be adopted. This is justified by the very thin openings e of the canaliculi or of the joints (e smaller than a few hundreds of nanometres, generally much smaller) as compared to the grain radii R which are around 10 μm (i.e. e ≪ R). At the micro-metric scale, the core of a grain G_i is assumed impermeable (K → 0). In turn, based on Poiseuille’s flow between two planes with aperture e, the flow in the half interface S_i(e) surrounding the grain G_i(e) is described by a surface Darcy’s law:$$\mathit{q}\left(\mathit{z}\right)=-\eta \left(e\right)\mathbf{gra}{\mathbf{d}}_{s}p\left(\mathit{z}\right)\hspace{1em}(\mathit{z}\in S_i(e\left)\right)$$(12)

where p(z) is the pressure at a point z of the interface, η is a surface permeability and q(z) is the surface velocity vector defined as the cumulative velocity across the half interface$$\mathit{q}\left(\mathit{z}\right)={\int }_0^{e/2}\mathrm{}\mathit{v}(\mathit{z}+Z\mathit{n})\mathrm{d}Z$$(13)

where n is the unit external normal to the boundary of the grain at point z. The reason to consider at this step only one half of the interface surrounding a grain G_i(e) is to avoid double counting in the estimation of the total flow through all the interfaces, which will later be retrieved by summation over all grains.

The classical solution to Poiseuille’s flow between two planes with aperture e and the definition of the surface velocity (Eq. 13) lead to the definition of the surface permeability in the law in Equation (12) as:$$\eta \left(e\right)=\frac{e^3}{24}$$(14)

for a half interface [19].

3.2.2 Partially Saturated Case

In the case of a partial water saturation S_w of the pore network, we assume that the water does not move due to gas flow. This decoupling in the extension of Darcy’s law to two-phase flow is frequently assumed in the capillary dominated flow regime (i.e. when viscous forces are negligible in comparison to capillary forces) [2]. It is here justified on a macroscopic level by the above presented experiments in which gas flow did not significantly modify the overall water saturation of the samples.

The interfaces which are water saturated are assumed to have a vanishing gas permeability. Consequently, the grains surrounded by water saturated interfaces behave as fully gas impermeable particles.

Further, it is assumed for simplification that a subdomain of volume fraction fS_w of high-porosity regions is totally water saturated and impermeable to gas (Fig. 8). The complementary subdomain of volume fraction f(1 − S_w) is assumed fully gas saturated, with unchanged gas permeability. These are strong assumptions which may have to be reconsidered in future works.

3.3 Determination of the Macroscopic Permeability by the Self Consistent Estimate

3.3.1 Darcy Medium Homogenisation Problem

Due to the previous simplifications consisting of replacing elementary particles by fictitious Darcy media, we are left with the homogenisation of a heterogeneous Darcy media. The homogenisation problem is defined on the RVE Ω by:$$\begin{array}{llll}& \mathrm{div}\mathit{v}=0& & (\mathit{z}\in \mathrm{\Omega })\\ & \mathit{v}=-{\mathit{K}}_{hp}\cdot \mathbf{grad}p& & (\mathit{z}\in {\mathrm{\Omega }}_{hp}^g)\\ & \mathit{v}=0& & (\mathit{z}\in {\mathrm{\Omega }}_{hp}^w\cup {\mathrm{\Omega }}_g)\\ & \mathit{q}=-\eta \left(e\right)\mathbf{gra}{\mathbf{d}}_{s}p& & (\mathit{z}\in S_i(e),e>e^{\star })\\ & \mathit{q}=0& & (\mathit{z}\in S_i(e),e<e^{\star })\\ & p=\nabla P\cdot \mathit{z}& & (\mathit{z}\in \mathrm{\partial }\mathrm{\Omega })\end{array}$$(15)

where z is the position vector, ∇P is the macroscopic pressure gradient which is applied uniformly at the boundary ∂Ω of the RVE and q is the cumulative velocity across the half interface defined by Equation (13).

A direct consequence of the uniform pressure gradient boundary condition is the so-called pressure gradient averaging rule, which relates the local pressure gradients to the macroscopic pressure gradient by:$$\nabla P=\frac{1}{\left|\mathrm{\Omega }\right|}{\int }_{\mathrm{\Omega }}\mathrm{}\mathbf{grad}p\mathrm{d}V\equiv \overline{\mathbf{grad}p}$$(16)

Decomposition of the volume average over the elementary particles yields:$$\nabla P=f(1-S_w){\overline{\mathbf{grad}p}}^{hp,g}+f{S}_w{\overline{\mathbf{grad}p}}^{hp,w}+(1-f){\int }_{e_{\min }}^{e^{\star }}\mathrm{}{\overline{\mathbf{grad}p}}^{G_i\left(e\right)}\alpha \left(e\right)\mathrm{d}e+(1-f){\int }_{e^{\star }}^{e_{\max }}\mathrm{}{\overline{\mathbf{grad}p}}^{G_i\left(e\right)}\alpha \left(e\right)\mathrm{d}e$$(17)

where the notation:$${\overline{\mathbf{grad}p}}^{hp,w}=\frac{1}{\left|{\mathrm{\Omega }}_{hp}^w\right|}{\int }_{{\mathrm{\Omega }}_{hp}^w}\mathrm{}\mathbf{grad}p\mathrm{d}V$$(18)

indicates the intrinsic volume average over the phase , and similarly for the phase or the grains G_i(e).

In turn, the macroscopic velocity V corresponds to the volume average of the local velocity. By linearity of the homogenisation problem (15), V is linearly related to the macroscopic pressure gradient by the macroscopic permeability K^hom:$$\mathit{V}=\bar{\mathit{v}}=-{\mathit{K}}^{\mathrm{hom}}\cdot \nabla P$$(19)

The overall volume average of the velocity can be decomposed on phases and grains similarly to Equation (17). Hence, the macroscopic gas permeability K^hom can be deduced from the knowledge of the phase averages (or elementary particles averages) of the pressure gradient and velocity.

3.3.2 Generalised Eshelby Problems

The self-consistent scheme relies on the estimate of the velocity and pressure gradient phase averages involved (e.g. in Eq. 17 and Eq. 19) based on the solution to auxiliary problems: the so-called Eshelby problems.

The auxiliary problems are devised as follows. For each phase, an inclusion I (or elementary particle representative of the phase) is embedded in an infinite medium (Fig. 8). The infinite medium has a uniform permeability K₀. An auxiliary uniform pressure gradient ∇P₀ is applied at infinity (remote boundary conditions). The resolution of this problem yields the average of the pressure gradient and of the velocity in the inclusion linearly as a function of the auxiliary pressure gradient ∇P₀. If the inclusion I is an ellipsoid with homogeneous permeability K_I, the pressure gradient and the velocity are actually constant inside the inclusion. In this case, the pressure gradient inside the inclusion is given by:$$\mathbf{grad}p\left(\mathit{z}\right)={\left[\mathbf{1}+{\mathit{P}}_0^I\cdot \left({\mathit{K}}_I-{\mathit{K}}_0\right)\right]}^{-1}\cdot \nabla P_0\hspace{1em}\left(\mathit{z}\in I\right)$$(20)

where is the so-called Hill tensor of the inclusion I in the infinite medium of permeability K₀, which is related to the Eshelby tensor by . In the case where K₀ is isotropic (i.e. K₀ = K₀1), the Eshelby tensor is independent of K₀. For a spherical or spheroidal (oblate or prolate) inclusion I, the Eshelby tensor decomposes as [1], with:$$\begin{array}{ll}{\mathit{S}}^{\mathrm{sphere}}& =\frac{\mathbf{1}}{3}\\ {\mathit{S}}^{\mathrm{prolate}}& =\frac{-\sqrt{{\omega }^2-1}+\omega \mathrm{arccosh}\left(\omega \right)}{({\omega }^2-1{)}^{3/2}}\mathit{n}\otimes \mathit{n}+\frac{\omega (\omega \sqrt{{\omega }^2-1}-\mathrm{arccosh}(\omega \left)\right)}{2({\omega }^2-1{)}^{3/2}}(\mathbf{1}-\mathit{n}\otimes \mathit{n})\\ {\mathit{S}}^{\mathrm{oblate}}& =-\frac{-\sqrt{1-{\omega }^2}+\omega \arccos \left(\omega \right)}{(1-{\omega }^2{)}^{3/2}}\mathit{n}\otimes \mathit{n}-\frac{\omega (\omega \sqrt{1-{\omega }^2}-\arccos (\omega \left)\right)}{2(1-{\omega }^2{)}^{3/2}}(\mathbf{1}-\mathit{n}\otimes \mathit{n})\end{array}$$(21)

where ω is the aspect ratio of the spheroid (ω > 1, prolate; ω < 1, oblate) and n is the unit vector collinear with the axis of revolution of the spheroid. The velocity inside the inclusion is then simply deduced from the combination of Darcy’s law with the Eshelby localisation rule (20) as v(z) = −K_I⋅grad p(z).

The case of a generalised Eshelby problem with a composite inclusion made of an impermeable core surrounded by a permeable interface (in the sense of Eq. 12) as been dealt with by [19]. A remarkable result is that the intrinsic average of the pressure gradient and of the velocity in the composite inclusion (i.e. including the interface) can be retrieved from Eq. 20 in the uniform spherical inclusion case for an appropriate choice of equivalent fictitious uniform permeability K^eq = K^eq(1) of the inclusion. If η denotes the surface permeability (Eq. 12) and R the radius of the composite inclusion, this equivalent fictitious permeability is:$$K^{\mathrm{eq}}=\frac{2\eta }{R}$$(22)

Another remarkable result is that K^eq does not depend on the permeability K₀ of the surrounding medium. This last feature is specific to a Darcy medium as it is untrue for linear elasticity (see e.g. [34], Eq. 37). Moreover, this equivalent permeability corresponds to the limiting case of a composite sphere made of a rigid core and a fluid shell of vanishing thickness in which an incompressible Newtonian fluid flows, with uniform normal velocity/no tangential stress boundary conditions on the outside of the composite sphere [11, 35, 36].

Alternatively, the result in Equation (22) can be retrieved from a generalised Eshelby problem with a composite spherical inclusion made of an impermeable core surrounded by a shell with thickness h and isotropic permeability η/h in the limit as h → 0. Furthermore, the resolution of a generalised Eshelby problem with a composite inclusion made of an impermeable spheroidal core surrounded by a confocal spheroidal shell with isotropic permeability η/h and vanishing average thickness h indicates that an equivalent transverse isotropic permeability tensor independent of K₀ can be similarly defined.

3.3.3 Self-Consistent Estimate

The principle of the self-consistent estimate is to choose the permeability K₀ of the infinite medium in the Eshelby problems equal to the homogenised permeability K^hom. In turn, the auxiliary pressure gradient ∇P₀ in the Eshelby problem is determined as a linear function of the actual macroscopic pressure gradient ∇P by inserting the Eshelby localisation rule of the pressure gradient in Equation (20) in the pressure gradient averaging rule (17). Combination with Equation (19) classically leads to an implicit equation on the self-consistent estimate of the macroscopic permeability K^hom:$${\mathit{K}}^{\mathrm{hom}}=\stackrel{̅}{\mathit{K}\cdot {\left[\mathbf{1}+{\mathit{P}}_{\mathrm{hom}}^I\cdot \delta \mathit{K}\right]}^{-1}}\cdot {\left(\stackrel{̅}{{\left[\mathbf{1}+{\mathit{P}}_{\mathrm{hom}}^I\cdot \delta \mathit{K}\right]}^{-1}}\right)}^{-1}$$(23)

where δK = K − K^hom. Equation (23) can be simplified into:$$\overline{\delta \mathit{K}\cdot [\mathbf{1}+{\mathit{P}}_{\mathrm{hom}}^I\cdot \delta \mathit{K}{]}^{-1}}=\mathbf{0}$$(24)

In the following developments, all the grains are assumed spheroidal with the same aspect ratio ω. The distribution of the orientations of the spheroidal grains is assumed isotropic. Based on the results recalled in the previous subsection, an equivalent gas permeability tensor K^eq (e) is associated with grains surrounded by an interface of thickness e > e^⋆ whereas the grains with water saturated interfaces are gas-impermeable. Denoting the unit vector collinear with the axis of revolution of the spheroid by n, the transverse isotropic equivalent permeability tensor decomposes as . The Eshelby and Hill tensors of the grains are respectively denoted and . In turn, the high porosity regions are described by spherical inclusions with gas permeability K^hp when they are gas saturated or 0 when they are water saturated. The Eshelby and Hill tensors of the spherical high porosity regions are respectively denoted and . Consequently, the volume average in Equation (24) can be split as in Equation (17) into:$$-f{S}_w{\mathit{K}}^{\mathrm{hom}}\cdot {\left[\mathbf{1}-{\mathit{S}}_{\mathrm{hom}}^{hp}\right]}^{-1}+f(1-S_w)({\mathit{K}}^{\mathrm{hp}}-{\mathit{K}}^{\mathrm{hom}})\cdot {\left[\mathbf{1}+{\mathit{P}}_{\mathrm{hom}}^{hp}\cdot ({\mathit{K}}^{\mathrm{hp}}-{\mathit{K}}^{\mathrm{hom}})\right]}^{-1}-(1-f){\int }_{e^{\min }}^{e^{\star }}\mathrm{}<{\mathit{K}}^{\mathrm{hom}}\cdot {\left[\mathbf{1}-{\mathit{S}}_{\mathrm{hom}}^g\left(\omega \right)\right]}^{-1}>\alpha \left(e\right)\mathrm{d}e+(1-f){\int }_{e^{\star }}^{e^{\max }}\mathrm{}<\left({\mathit{K}}^{\mathrm{eq}}\right(e)-{\mathit{K}}^{\mathrm{hom}}){\left[\mathbf{1}+{\mathit{P}}_{\mathrm{hom}}^g\cdot \left({\mathit{K}}^{\mathrm{eq}}\right(e)-{\mathit{K}}^{\mathrm{hom}})\right]}^{-1}>\alpha \left(e\right)\mathrm{d}e=\mathbf{0}$$(25)

where < … > indicates the angular average over all directions of the orientations of the spheroids.

Under the hypothesis of isotropy, the macroscopic permeability is K^hom = K^hom1 and Equation (24) reduces to a single scalar equation. Making use of:

• the orthogonality properties of the tensors n ⊗ n and 1 − n ⊗ n ,

• the angular average of n ⊗ n for an isotropic distribution of orientations < n ⊗ n > = 1/3,

• the definition of the fraction s of water saturated interfaces (11) allows simplification of Equation (25) to:

$$-f{S}_w\frac{3}{2}+f(1-S_w)3\frac{K^{\mathrm{hp}}/K^{\mathrm{hom}}-1}{2+K^{\mathrm{hp}}/K^{\mathrm{hom}}}-(1-f)s\frac{1}{3}\left[\frac{1}{1-S_n\left(\omega \right)}+\frac{2}{1-S_t\left(\omega \right)}\right]+(1-f)\frac{1}{3}{\int }_{e^{\star }}^{e^{\max }}\mathrm{}[\frac{K_n^{\mathrm{eq}}\left(e\right)/K^{\mathrm{hom}}-1}{1+S_n\left(\omega \right)\left(K_n^{\mathrm{eq}}\right(e)/K^{\mathrm{hom}}-1)}+2\frac{K_t^{\mathrm{eq}}\left(e\right)/K^{\mathrm{hom}}-1}{1+S_t\left(\omega \right)\left(K_t^{\mathrm{eq}}\right(e)/K^{\mathrm{hom}}-1)}]\alpha \left(e\right)\mathrm{de}=0$$(26)

The percolation threshold of the self-consistent scheme is reached for the critical water saturation and interface water saturation s^c such that K^hom → 0. Inserting this limiting case in Equation (26) indicates that the critical water saturations for the present model verify:$$3f\left[1-\frac{3}{2}{S}_w^c\right]+\frac{1-f}{3}[(1-s^c)\left(\frac{1}{S_n\left(\omega \right)}+\frac{2}{S_t\left(\omega \right)}\right)-s^c\left(\frac{1}{1-S_n\left(\omega \right)}+\frac{2}{1-S_t\left(\omega \right)}\right)]=0$$(27)

irrespective of the values of K^hp and K^eq (e), as long as they are not null. The influence of the grains aspect ratio ω and of the high porosity regions volume fraction f is illustrated in Figure 9 under the assumption that the interface water saturation s is equal to the total water saturation S_w.

In order to provide an explicit expression of the relative gas permeability, we choose to use a spherical shape of the inclusions in the Eshelby problems. The grains are assumed to all have the same radius R. Further, the experimental results presented in Figure 6 hint that a log-uniform distribution of interface apertures is a satisfactory approximation, that is:$$\alpha \left(e\right)=\frac{1}{\ln \left(e_{\max }/e_{\min }\right)}\times \frac{1}{e}$$(28)

For the log-uniform distribution, the median interface aperture e_med, the maximum to minimum aperture ratio ρ and the critical aperture e^⋆ which separates water saturated from gas saturated interfaces are:$$\begin{array}{ll}{e}_{\mathrm{med}}& =\sqrt{e_{\min }{e}_{\max }}\\ \rho & =e_{\max }/e_{\min }\\ e^{\star }& =e_{\min }^{1-s}{e}_{\max }^s=e_{\max }{\rho }^{s-1}\end{array}$$(29)

In the case where the high porosity region permeability K^hp ≫ K^hom, the solution to equation (26) is then explicit:$$K^{\mathrm{hom}}=\frac{e_{\mathrm{med}}^3}{24R}\frac{{\rho }^{\frac{3}{2}}-{\rho }^{\frac{1-3f(1-S_w)}{1-f}+\frac{3}{2}(2s-1)}}{{\rho }^{\frac{1-3f(1-S_w)}{1-f}}-1}$$(30)

The approximation K^hp ≫ K^hom can be a posteriori justified by the fact that the high porosity region permeability scales as the square of a characteristic pore size () whereas the homogenised permeability K^hom (Eq. 30), like the equivalent grain permeability (Eq. 22), scales as with e_med ≪ R.

3.4 Comparison to Experimental Results

The explicit estimate Equation (30) of the homogenised gas permeability is compared to experimental results in Figure 10, under the hypothesis that the interface water saturation s is equal to the total water saturation S_w and that the high porosity region volume fraction is equal to the total porosity ϕ. Based on the experimental results, we used ϕ = 0.11 (Tab. 2), e_min = 0.6 nm and e_max = 360 nm (Fig. 6) such that no parameter has to be fitted. The radius R of the grains has no impact on the relative permeability curve as results are normalised by the permeability corresponding to the dry case. It is emphasised that the threshold saturation is an output of the model. In the present case, where spheres have been simply chosen to represent the grains, .

The agreement is good for water saturations below 0.3, but beyond that the model tends to over-estimate the relative gas permeability due to an inflexion point in the curve of the model. The percolation threshold given by the self-consistent scheme for spherical grains and S_w = s is slightly below the actual critical water saturation (≈ 0.7). Choosing prolate (elongated) spheroids with aspect ratio ω = 2 to describe the grains would have led to a better agreement for the critical water saturation. Note however that all these remarks are very sensitive to the choice of the s − S_w relation, which here has a priori been set to s = S_w for the sake of simplicity.

Conclusion

An experimental characterisation of the porosity, capillary pressure curve, dry and relative gas permeabilities has been carried out on a set of samples of a low permeability caprock. Two methods have been used to measure the porosity: with water or gas. Water porosities, around 9%, are systematically below gas porosities, around 11%. This difference could be attributed to trapped gas in the “water porosity” measurements or to gas sorption on pore walls when performing “gas porosity” measurements, and has been observed on other materials such as argillite. The interpretation of the capillary pressure curves as approximate pore size distributions indicates that 2/3 of the pore volume is comprised of very small pores, below 100 nm opening. Based on SEM imaging, these small sizes are attributed to the clay content in high porosity regions, to canaliculi crossing carbonate aggregates and to joints around silts.

The Klinkenberg-corrected intrinsic gas permeabilities on dry samples are low, around 2 × 10⁻¹⁹ m² for most samples. The relative gas permeability measurements demonstrate a strong repeatability as all samples lie on the same curve. The low confining pressures which have been used (9 MPa ≈ as in-situ and 2 MPa) did not have a significant impact on the porosity or relative gas permeability measurements. The important feature of the relative permeability curve is a gas permeability cut-off for critical water saturations around 0.71 ± 0.04. This has strong implications for gas storage applications, since it implies that the caprock can prevent gas flow even in the case of low desaturation. It means that a caprock exhibits two interesting features in term of confinement properties. The first is the standard one which consists of a static barrier related to the threshold capillary pressure value. The second one evidenced by this work consists of a dynamic barrier, related to the gas relative permeability curve shape, which first prevents and then delays the gas flow at low saturations.

The proposed micro-mechanical model relies on a strong simplification of the microstructure organisation. Based on several SEM images, gas flow has been assumed to be controlled by hydraulic interfaces corresponding to canaliculi surrounding impermeable carbonate aggregates or joints along impermeable silts. These interfaces have been assumed to connect high porosity regions which are partially clay-filled. The contribution of the hydraulic interfaces has been taken into account thanks to a composite grain model recently proposed in the literature. Based on a simplification of the Newtonian flow problem into an equivalent heterogeneous Darcy problem, the self-consistent scheme has been used to assess the macroscopic effective gas permeability. The model includes an interface aperture size distribution and allows study of the effect of progressive water saturation. The percolation threshold naturally arising in the self-consistent scheme allowed modeling of the critical water saturation corresponding to the gas permeability cut-off. A closed-form model has been proposed, whose parameters may all be deduced from porosimetry measurements. The agreement with experimental data is fair for the relative gas permeability at low water saturations (S_w < 0.3) and for the value of the critical water saturation. However, the relative gas permeability is still poorly modeled at intermediate water saturations (0.3 < S_w < 0.6). Further improvements could be made based on a more accurate description of the progressive and perhaps differential localisation of the water in the interfaces or in the high porosity regions.

Acknowledgments

The authors acknowledge financial support from ENGIE E&P International and Storengy. The authors would like to thank Bertrand van de Moortèle and Gilles Dromart from ENS Lyon for the petrographic analysis of the tested rock.

References

All Tables

All Figures