© C.D. Tsakiroglou et al., published by IFP Energies nouvelles, 2021

1 Introduction

Tight and shale reservoir systems pose a tremendous potential resource for future development, and the determination of their microscopic and macroscopic properties has long attracted the attention of scientific society. In particular, shale gas reservoirs possess many so-called “unconventional” features, the traditional approaches are unable to determine the reservoir permeability, and both the micro-scale analysis of flow regimes and macroscopic estimation of extremely low permeability pose emerging challenges (Hill and Nelson, 2000; Javadpour et al., 2007).

Depending on the pore size range, and prevailing conditions, the flow of natural gas through conventional porous rocks, tight formations and shale systems, is controlled by a variety of flow regimes, including Knudsen, transition, slip and viscous Darcy flow (Ziarani and Aguilera, 2012). Having known and understood the rock properties, including pore geometry, flow/transport parameters, the determination of hydrocarbon storage and recovery from shales becomes tractable. The microstructural characterization of shale samples could be achieved by combining Mercury Intrusion Porosimetry (MIP), with low‐field Nuclear Magnetic Resonance (NMR), N₂ adsorption/desorption, and high resolution Focused Ion Beam–Scanning Electron Microscopy (FIB–SEM) images (Al Hinai and Rezaee, 2015).

The gas permeability of shales can be correlated with the capillary pressure data from mercury injection measurements, by using theoretical and empirical equations (Al Hinai and Rezaee, 2015). On the other hand, the pore network approach is a practical way to explain the macro flow behavior of porous media from a microscopic point of view, so that the shale matrix is simulated by 3-dimensional pore network models that include typical bimodal pore size distribution, anisotropy and low connectivity (Zhang et al., 2017). Moreover, the digital shale rock could be developed by 3-Dimensional (3D) images obtained from micro/nano Computed Tomography (CT) and FIB-SEM images reconstructed from 2D SEM images of the pore structure, and using the Lattice Boltzmann method to compute the intrinsic permeability, porosity and tortuosity; these parameters are used to calculate the apparent permeability by accounting for different gas transport mechanisms (Sun et al., 2017).

Three-dimensional nano-porous structures of shales reconstructed from SEM images of shale samples enabled the computational analysis of structural characteristics and calculation of tortuosity, effective Knudsen diffusivity and intrinsic permeability by Lattice-Boltzman method. It was found that Knudsen diffusion is always of key importance on shale gas transport mechanisms in the reconstructed shales, while no Darcy (viscous) flow regime was observed (Chen et al., 2015). The FIB-SEM technique was employed to observe and characterize the morphology of the nanopores in shales, Lattice Boltzmann Method (LBM) was used to simulate high-Knudsen gas flows, and the nano-scale pore morphology was found to play an important role on the gas transport properties of shale matrix (Zheng et al., 2017). A method was developed to identify organic and inorganic pores from SEM images of shale samples collected from the wells of a productive formation, and it was found that: (i) most of inorganic pores are mesopores; (ii) the sizes for organic pores tend to be lognormal distributions with peak in the range 10–30 nm; (iii) the gas storage space obtained by SEM is some orders of magnitude larger than that obtained by gas adsorption method; (iv) the permeability is strongly anisotropic, and a significant non-Darcy effect is evident (Cao et al., 2019). The image-based fractal characteristics of shale pore structure were examined from 100 SEM images, at multiscale resolutions ranging from 15 to 420 nm, their effect on the accurate prediction of gas permeability was analyzed, and it was found that the box counting method enables the accurate computation of fractal dimension, compared to Sierpinski carpets analytical solution, and the gas permeabilities approach to the exact values when the resolution is better than 50 nm (Song et al., 2019). Four major techniques, including N₂ adsorption, SEM, MIP, and NMR, were employed to characterize the Pore Size Distributions (PSDs) of a shale oil reservoir, and it was found that there is a wide range of pore sizes ranging from micropores to macropores (>1 μm), interconnected by small-scale throats, while the comprehensive characterization of shale oil reservoir pore structures requires the combination of more than one techniques, like NMR and MIP (Zhang et al., 2019). A sophisticated algorithm of pore-and-throat network reconstruction from 3D SEM images with resolution of 4 nm was developed by accounting for complex geometry and extremely small sizes of the pore-bodies and throats, and validated for a Marcellus shale core sample (Zheng et al., 2019).

In the present work, with the aid of percolation theory, analytic models of the N₂ adsorption–desorption isotherms and Hg intrusion curve are formulated as functions of the pore and throat radius distributions and pore network connectivity (Tsakiroglou et al., 2004). Then, a stepwise methodology is adopted to estimate the pore network parameters with inverse modeling of N₂ sorption and Hg intrusion curves and applied to datasets of several shale samples. The critical path analysis of percolation theory is used to calculate the liquid permeability and Knudsen diffusion coefficient of porous media explicitly from the aforementioned pore structure parameters and compare the predictions with experimental measurements.

2 Materials and methodology

2.1 Shales and transport properties

The Perth Basin is a 100 000 km² area covering the Western Australian margin between Augusta and Geraldton. Samples were collected at various depths from the same well of Carynginia shale gas formation (Al Hinai and Rezaee, 2015). The N₂ adsorption/desorption isotherms were measured with Micromeritics TriStar II 3020 (pore sizes probed: 1.3–200 nm), whereas the Hg intrusion curve was measured with Micrometrics Autopore IV Porosimeter (pressure range: 1–60 000 psi; pore sizes probed: 400 μm–3.5 nm). The liquid permeability and all pertinent Klinkenberg parameters were measured with a modified steady-state method, where the sample was embedded in a resin disc, the gas flow rate was measured at the outlet, the average gas pressure was in the range of 1–7 bar, and the net confining pressure was 70 bar (Al Hinai et al., 2013). The properties of samples, measured by the various techniques, are shown in Table 1.

2.2 Pore-scale models of drainage and imbibition

The actual pore space of sedimentary rocks, including shales (Fig. 1a) is complex, and it is commonly characterized in terms of pore-and-throat networks (Fig. 1b). In the present work, the porosity of shales is regarded as a stochastic network of pores (sites)-and-throats (bonds), the sizes of which are sampled from bimodal site f _s (r) and bond f _b (r) radius distributions, both composed of log-normal component functions, namely:

(1a)

(1b)where <r _s1>, σ _s1, <r _s2>, σ _s2 are the mean radii and standard deviations of the two component site distributions with contribution fractions to the overall one c _s and 1 − c _s , respectively, and <r _b1>, σ _b1, <r _b2>, σ _b2 are the mean radii and standard deviations of the two component bond distributions with contribution fraction to the overall one c _b and 1 − c _b , respectively.

The throats (bonds) are considered as narrow constrictions without volume, while the pores are assigned the entire pore volume according to the relation:

(2)where β _s is a volume shape factor (β _s  ≥ 0). The parameter value β _s  = 2.0, implies a network of prismatic pores with their volume proportional to r ². On the other hand, β _s  = 1.0, implies a network of slit-shaped pores with their volume proportional to r ¹. So any value in the range 1 < β _s  < 2 corresponds to a mixture of pore shapes between the foregoing limits.

The representation of the pore shapes in reservoir rocks with angular geometries contributes to their realistic characterization in terms of pore-and-throat networks (Tsakiroglou and Payatakes, 1993). On the other hand, angular geometries allow the rational computation of the oil and water flow rates, and wettability effects on simulated relative permeability curves (Firincioglu et al., 1999; Tsakiroglou, 2014; Valvante et al., 2005; Zhou et al., 1997). The cross-section of each pore is modeled by a high-order regular polygon with the number of sides (angles), area, and perimeter equal to n _s , A, and P, respectively (Fig. 2a). The pore shape factor, G, is given (Tsakiroglou et al., 2004) by,

(3)and the pore radius, r, coincides to the radius of the inscribed circle (Fig. 2a), given by,

(4)

2.3 Models of N₂ adsorption/desorption in pore-and-throat networks

The graual saturation of pores with liquid nitrogen takes place through two parallel mechanisms (Fig. 2b): (1) gas adsorption in multiple layers across the pore-walls, and (2) capillary condensation of liquid at the corners of adsorbed layer. The thickness of the adsorbed layer, t _c , is governed by the generalized Frenkel–Halsey–Hill (FHH) equation (Gregg and Sing, 1982),

(5)where x is the relative vapor pressure, σ is the diameter of N₂ molecule (σ = 3.54 A), the parameter b depends on the properties of adsorbent and adsorptive, and s is an index with values ranging from 2 to 3. The radius of curvature of the condensate/vapor interfaces, r _c , is given by Kelvin equation (Gregg and Sing, 1982)

(6)

To get a visual view of the calculated values and fluid saturation, we first assume that β _s  = 2.0 (prismatic pores) and then we generalize the so-produced equations to any β _s value. The liquid N₂-saturation in a pore includes the adsorbed and condensed phases and is given by,

(7)where F _LG is the fraction of the pore cross-section area occupied by liquid N₂, and depending on the liquid/gas/solid contact angle, θ _LG (Fig. 2c; we assume that θ _LG = 0°, for liquid N₂, unless otherwise stated). If equation (7) is generalized for any β _s value, then it is written as,

(8)

With regard to the fluid saturation history, the N₂ capillary condensation resembles the imbibition of a wetting fluid (liquid N₂), and the critical radius of curvature for condensation, r _c  = r _con, is given by,

(9a)

(9b)

Equations (9a) and (9b), in conjunction with equations (5) and (6), yield the critical radius r = r _s , given by,

(10a)

(10b)

The saturation of a pore network with liquid N₂ can be expressed as a function of the relative vapor pressure by the relation,

(11)

The capillary evaporation of liquid N₂ from a pore resembles the drainage of the wetting fluid (liquid N₂), and hence the critical radius of curvature for evaporation, r _evp , is given by,

(12a)

(12b)which in combination with equations (5) and (6) lead to,

(13a)

(13b)where r _b is the critical radius of throats which are “allowable” to N₂ evaporation (r ≥ r _b ). N₂ empties a pore if its radius is less than the critical size (r _b or r _s ) and is accessible either to the bulk vapor phase (primary desorption-drainage) or to isolated vapor pockets (secondary desorption-drainage). During the N₂ desorption process, the fraction of bonds (throats) q _b and sites (pores) q _s which are “allowable” to the vapor phase are given by,

(14)respectively. In percolation theory, the accessibility function depends on pore network, topology, size, and statistical spatial pore size correlations, and is defined as the fraction of allowable to a fluid pores that are accessible to this fluid (Tsakiroglou et al., 2004, 2009). In the specific case of N₂ desorption, the accessibility function relates the fraction of pores which are accessible to vapor N₂ with the fraction of throats that are allowable to vapor N₂. Finally, the liquid N₂ saturation during desorption is given by (Tsakiroglou et al., 2004),

(15)

2.4 Model of Hg intrusion in pore-and-throat network

In general, the value of pore volume calculated from N₂ sorption does not coincide to the corresponding value calculated from Hg intrusion curve (Tab. 1). Although Hg is a non-wetting fluid for the majority of materials, by “convention” the Hg contact angle, θ _Hg, is measured with respect to the wetting fluid (low pressure air) and is considered always lower than 90°. In order to estimate a unique pore volume value, V _p (cm³/g), the following method is adopted. For a given Hg intrusion contact angle, the maximum intrusion capillary pressure is converted to an equivalent relative vapor pressure. This value is interpolated in N₂ desorption curve to calculate the corresponding liquid N₂ volume. The sum of the two volumes yields the total pore volume. The procedure is repeated for various values of θ _Hg, and the produced datasets (V _p vs. θ _Hg) are fitted to a non-linear parametric function of the form,

(16)so that the values of A ₁, A ₂, θ ₀, w are estimated.

The Hg intrusion in a pore-and-throat network is analogous to capillary evaporation for q _si  = 0 (primary desorption) and is controlled by throat sizes. The mercury saturation in a pore network is given (Tsakiroglou et al., 2004, 2009) by,

(17)where,

(18)

(19a)

(19b)

2.5 Accessibility functions

The accessibility functions can be described satisfactorily by the following functional form,

(20)where is the initial fraction of sites occupied by the non-wetting fluid; for Hg intrusion (primary drainage, i = 0) and for nitrogen desorption (secondary drainage initiating before the completion of adsorption, i = 1). Evidently, two different sets of parameters (a ₀, b ₀) and (a ₁, b ₁) are required to quantify the accessibility functions Y _s0(q _b ) and Y _s1(q _b ) that specify the Hg intrusion and N₂ desorption in pore networks, respectively. A “physical” parameter characterizing the network accessibility function is the “percolation threshold”, q _b  = q _bci , which corresponds to the inflection point of the curve and is obtained from the numerical solution of the following equation,

(21)

2.6 Critical Path Analysis (CPA)

The transport properties of porous materials, characterized in the terms of a parametric pore-and-throat network, with pore sizes spanning some orders of magnitude, can be computed either explicitly with simulation of the corresponding process (e.g. single phase flow) in a reconstructed network (Jones et al., 2018), or implicitly with the use of approximate methods, such as the CPA of percolation theory (Tsakiroglou and Ioannidis, 2008; Tsakiroglou et al., 2009). The main idea hidden behind the CPA is that for a pore system with a broad pore size distribution most of the flow is transferred by a network spanning cluster consisting of pores with conductance greater than a critical value, g _c (Katz and Thompson, 1986). This critical conductance g _c depends on pore space characteristics (e.g. pore shape, pore size distribution, connectivity, dimensionality, etc), and CPA enables the calculation of the corresponding critical pore dimensions. The transport properties of the critical path may be governed by universal scaling laws allowing the approximate computation of accessibility function at the vicinity of percolation threshold (Tsakiroglou et al., 2009). In mathematical terms, the viscous flow or Knudsen flow problem reduces to the calculation of the critical pore dimensions that maximize the total conductance of the critical path (Hunt, 2001; Tsakiroglou and Ioannidis, 2008; Tsakiroglou et al., 2009). In general, the dependence of the local viscous flow conductance on pore size may be different from that of the local Knudsen diffusion conductance, so that the critical pore dimension may depend on the transport mechanism (Tsakiroglou et al., 2009). The gas molar flux in a porous medium for the transition region between Knudsen diffusion and viscous flow is given by the sum of the individual contributions, namely,

(22)where k _L is the liquid permeability, and D _kn is the Knudsen diffusion coefficient of the critical path. According to the CPA (Tsakiroglou et al., 2009), these properties are approximated by the relationships,

(23)

(24)where φ is the porosity, t is the universal scaling exponent of conductivity (t = 2.0 for 3-D systems), and mean pore radius, <r _sq > can be expressed as function of throat radius, r, by,

(25)

The values of k _L and D _kn are maximized for r = r _max,f and r = r _max,d , respectively, and are considered approximately equal to the permeability and Knudsen diffusion coefficient of the investigated porous medium.

For low permeable porous media, the measured gas permeability is pressure-dependent and is commonly expressed by the relationship,

(26)where b is the Klinkenberg factor, and P _m is the mean gas pressure (Al Hinai et al., 2013). For the steady-state gas flow through a porous medium, the integrated form of Darcy equation, including the Klinkenberg effect, takes the form,

(27)where u ₂ , P ₂ are the flow velocity and pressure at the outlet, L is the sample length, μ _g is the gas viscosity, and ΔP is the pressure drop across the sample. By combining equation (22) with equation (27), it is obtained,

(28)

2.7 Demonstration of methodology

For the non-linear parameter estimation, numerical codes in the software platform of Athena Visual Studio (Stewart and Caracotsios, 2008) were developed, and a step-wise procedure was adopted (Fig. 3). First, the parameters n _s , β _s , and f _s (r; <r _s1>, σ _s1, <r _s2>, σ _s2, c _s ) were estimated with inverse modeling of N₂ adsorption isotherm. Then, the foregoing parameters were fixed, and f _b (r; <r _b1>, σ _b1, <r _b2>, σ _b2, c _b ) along with Y _si (q _b ) were estimated with inverse modeling of N₂ desorption isotherm. Then, Y _s0 (q _b ) was estimated from inverse modeling of the Hg intrusion curve. Finally, the throat size distribution, secondary and primary pore accessibility functions, and Hg contact angle were re-estimated, with the simultaneous inverse modeling of N₂ adsorption/desorption isotherms and Hg intrusion curve.

3 Results and discussion

3.1 Inverse modeling of Hg intrusion curve and N₂ sorption isotherms

It is well-known that for low porosity samples, surface irregularities (e.g. nooks, crannies, crevices) and inter-granular pores that are filled with Hg over the low pressure range are interpreted as macro-pores. It is of key importance to exclude such characteristics from the deconvolution of Hg intrusion and N₂ isotherms, and pay attention mostly on mesopores and micropores. In this context, the Hg intrusion curves of all samples were corrected by removing the part of the curve extending over low pressures (Fig. 4).

For the three shale samples (Tab. 1), the Hg intrusion data were combined with N₂ desorption isotherms to calculate the total pore volume as a function of the Hg contact angle, and the results were fitted with equation (16) to estimate all pertinent parameters (Tab. 2, Fig. 5). The V _p (θ _Hg) relationship was incorporated into the set of equations modeling Hg intrusion and N₂ sorption processes.

The application of the methodology of pore structure characterization from N₂ adsorption/desorption isotherms and Hg intrusion curves was applied to the datasets of the three shale samples (Tab. 1), and the results are shown in Table 3.

The measured N₂ isotherms are predicted satisfactorily by simulated ones for all cases (Figs. 6a, 7a, 8a and 8b), indicating that the estimated statistics of pore and throat sizes over the low size range (<0.5 μm) are consistent with the estimated accessibility function. On the other hand, regarding the Hg intrusion curves, there is a discrepancy between the measured and predicted curve over the low and very high pressure ranges (Figs. 6b, 7b and 8c). It seems that, in spite of the large number of “unknown” parameters and the iterative method of inverse modeling, we are unable to track precisely the N₂ adsorption/desorption curves, and simultaneously the Hg intrusion curve. The initial guesses of TSD, PSD, and ACF are obtained from N₂ adsorption (PSD) and desorption (TSD and ACF) curves, and unavoidably the finally estimated PSD and TSD might have been biased toward the low radii sizes. On the other hand, large throat radii that are reflected in the low pressure range of Hg intrusion curve (Figs. 6b, 7b and 8c) are overlooked since they are not “shown” at all in the N₂ desorption curve. Other factors that prevent the simultaneous prediction of all datasets and favor the shift of PSD and TSD toward small sizes are: (i) the higher contribution fraction of N₂ sorption to the estimation process (2 datasets against 1 dataset); (ii) the uncertainty of the “fixed values” of parameters included in N₂ adsorption/desorption equations (e.g. liquid N₂ contact angle); (iii) the lack of sufficient information concerning the large pore sizes that are not detected at all by N₂ adsorption isotherm, and are shadowed in Hg intrusion curve; such large pore sizes might be correlated with the Hg retraction curve (Tsakiroglou et al. 2009) which was overlooked in the present approach.

In the case of AC5, the significant discrepancy observed between observed and simulated Hg intrusion curve (Fig. 8c) led us to repeat the last step of parameter estimation [f _b (r), Y _si (q _b ), Y _s0(q _b ), θ _Hg] by modifying the initial guesses. Finally, the estimated TSD widens toward large radii, and ACFs shift to higher q_b values, characterized by higher percolation thresholds. Not having imposed any limitation on the ratio of pore to throat sizes, the overlapping between PSD and TSD may be discernible (Figs. 9a–9c), with the mean value of estimated TSD exceeding that of estimated PSD (Tab. 3, AC5, S1/S2). The significant overlapping of TSD and PSD (Fig. 9) leads to comparable values of q _bc0 and q _bc1 (for AC1 and AC8), while in the case of AC5, with the throat sizes being greater than the pore sizes, q _bc0 is much less than q _bc1, which is opposite to conventional predictions of percolation-type approaches. The physical interpretation of such a quantitative description is that the critical capillary sizes for Hg penetrating in throats are comparable to those of Hg penetration in interconnecting pores. Probably, a simpler model like a capillary channel network, described by only one pore size distribution and one accessibility function, might be more convenient than the present model for the representation of the pore space.

3.2 Prediction of transport properties

The parameter values of simulated throat-and-pore networks (Tab. 3, Fig. 9) were used to compute k _L and D _kn as functions of throat radius, r, through the approximate equations (23) and (24), respectively (Figs. A1–A3). The maximum k _L,max or D _kn,max occur at the minimum throat radius, over which the interconnected critical path of throats and pores transfers most of the viscous (Figs. A1a, A2a and A3a) or Knudsen (Figs. A1b, A2b and A3b) flow, respectively. The so-predicted permeability and Knudsen diffusion coefficient are compared with corresponding experimental values (Tab. 4). It seems that D _kn is predicted satisfactorily for almost all cases (Tab. 4), while the permeability is almost predictable only for sample AC5, while it fails for AC1 and AC8 (Tab. 4).

The Knudsen number is commonly used to classify flow regimes in fine pores where deviation from continuum flow is important (Ziarani and Aguilera, 2012). It is defined as the ratio of molecular mean free path, λ (nm) to the mean pore radius, r _s (nm), and is given by,

(29)

If K _n  < 0.01, then viscous flow dominates, and the conventional fluid dynamics equations and Darcy equation are applicable. If 0.01 < K _n  < 0.1, then slip flow regime occurs and Darcy equation with Klinkenberg or Knudsen modification is applicable. If 0.1 < K _n  < 10, then the transition flow regime prevails, where both slip (continuum) and diffusion flows can occur, and the traditional fluid dynamics equations start to fail. The higher the Knudsen number, the higher the chance of failure. Even though conventional equations could be applied (i.e., Darcy with Knudsen correction), the validity of such a formulation is questionable. It is safer to use Knudsen diffusion equations, especially for flow at higher Knudsen numbers (close to 10). If K _n  > 10, then the Knudsen (free molecular) flow regime prevails, where the continuum fluid flow equations fully break down, and diffusion-based formulations must be applied. Looking at the post-calculated Knudsen number, K _n (Tab. 4), the transition (AC5) and Knudsen (AC1, AC8) flow regimes (AC1, AC8) prevail. Therefore, Knudsen flow is dominant over all cases (AC1, AC5, AC8), but only in the one case (AC5), viscous flow becomes evident. Based on the aforementioned interpretation, failing to predict the value of k_L for AC1 and AC8, and succeeding to predict satisfactorily D _kn for all cases seem reasonable.

4 Conclusion

Analytical models of Hg intrusion and N₂ adsorption-desorption in porous materials are used for the inverse modeling of experimental datasets of shale samples and estimation a full set of parameters for the pore space, regarded as a stochastic pore-and-throat network. Using the critical path analysis of percolation theory, analytical relationships are developed for the explicit computation of the Klinkenberg parameters (liquid permability and Knudsen diffusion coefficient) of the gas permeability from pore space parameters. The methodology is evaluated with demonstration to the datasets of three shale samples:

• The estimated N₂ adsorption/desorption curves are predicted satisfactorily, while the predicted Hg intrusion curve deviates from the experimental one, over low and high pressure range.

• The initial guesses of parameters are based on N₂ sorption isotherms, so that the estimated throat and pore size distributions are biased over the low sizes.

• The overlapping between the pore and throat size distributions is significant for the investigated shales.

• Application of the methodology to datasets of shales indicates that, based on the estimated pore space parameters, the liquid permeability is predictable at moderate K _n values, over which the transition regime is dominant, while the Knudsen diffusion coefficient is also predictable at high K _n values, over which Knudsen flow is dominant.

Acknowledgments

Dr. C.D. Tsakiroglou acknowledges support of this work by the project POLITEIA II (MIS 5002478) which is implemented under the “Action for the Strategic Development on the Research and Technological Sector”, funded by the Operational Programme “Competitiveness, Entrepreneurship and Innovation” (NSRF 2014–2020) and co-financed by Greece and the European Union (European Regional Development Fund).

Appendix

References

All Tables

All Figures