A methodology to predict the gas permeability parameters of tight reservoirs from nitrogen sorption isotherms and mercury porosimetry curves

A methodology is suggested for the explicit computation of the absolute permeability and Knudsen diffusion coefficient of tight rocks (shales) from pore structure properties. The pore space is regarded as a pore-and-throat network quantified by the statistical moments of bimodal pore and throat size distributions, pore shape factors, and pore accessibility function. With the aid of percolation theory, analytic equations are developed to express the nitrogen (N2 ) adsorption/desorption isotherms and mercury (Hg) intrusion curve as functions of all pertinent pore structure parameters. A multistep procedure is adopted for the successive estimation of each set of parameters by the inverse modeling of N2 adsorption–desorption isotherms, and Hg intrusion curve. With the aid of critical path analysis of percolation theory, the absolute permeability and Knudsen diffusion coefficient are computed as functions of estimated pore network properties. Application of the methodology to the datasets of several shale samples enables us to evaluate the predictability of the approach.


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 2 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 . 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 2 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 lm), 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 . 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 2 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 2 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.

Materials and methodology 2.1 Shales and transport properties
The Perth Basin is a 100 000 km 2 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 2 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 lm-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.

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: f s r ð Þ ¼ c s f s1 r; hr s1 i; r s1 ð Þþ 1 À c s ð Þf s2 r; hr s2 i; r s2 ð Þ ; ð1aÞ where <r s1 >, r s1 , <r s2 >, r 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 >, r b1 , <r b2 >, r 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: where b s is a volume shape factor (b s ! 0). The parameter value b s = 2.0, implies a network of prismatic pores with their volume proportional to r 2 . On the other hand, b s = 1.0, implies a network of slit-shaped pores with their volume proportional to r 1 . So any value in the range 1 < b 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, and the pore radius, r, coincides to the radius of the inscribed circle ( Fig. 2a), given by, 2.3 Models of N 2 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),  where x is the relative vapor pressure, r is the diameter of N 2 molecule (r = 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) To get a visual view of the calculated values and fluid saturation, we first assume that b s = 2.0 (prismatic pores) and then we generalize the so-produced equations to any b s value. The liquid N 2 -saturation in a pore includes the adsorbed and condensed phases and is given by, where F LG is the fraction of the pore cross-section area occupied by liquid N 2 , and depending on the liquid/gas/solid contact angle, h LG ( Fig. 2c; we assume that h LG = 0°, for liquid N 2 , unless otherwise stated). If equation (7) is generalized for any b s value, then it is written as, With regard to the fluid saturation history, the N 2 capillary condensation resembles the imbibition of a wetting fluid (liquid N 2 ), and the critical radius of curvature for condensation, r c = r con , is given by, Equations (9a) and (9b), in conjunction with equations (5) and (6), yield the critical radius r = r s , given by, The saturation of a pore network with liquid N 2 can be expressed as a function of the relative vapor pressure by the relation, The capillary evaporation of liquid N 2 from a pore resembles the drainage of the wetting fluid (liquid N 2 ), and hence the critical radius of curvature for evaporation, r evp , is given by, which in combination with equations (5) and (6) lead to, where r b is the critical radius of throats which are "allowable" to N 2 evaporation (r ! r b ). N 2 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 2 desorption process, the fraction of bonds (throats) q b and sites (pores) q s which are "allowable" to the vapor phase are given by, 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(Tsakiroglou et al., , 2009). In the specific case of N 2 desorption, the accessibility function relates the fraction of pores which are accessible to vapor N 2 with the fraction of throats that are allowable to vapor N 2 . Finally, the liquid N 2 saturation during desorption is given by (Tsakiroglou et al., 2004), See Equation (15) top of the next page

Model of Hg intrusion in pore-and-throat network
In general, the value of pore volume calculated from N 2 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, h 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 3 /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 2 desorption curve to calculate the corresponding liquid N 2 volume. The sum of the two volumes yields the total pore volume. The procedure is repeated for various values of h Hg , and the produced datasets (V p vs. h Hg ) are fitted to a non-linear parametric function of the form, so that the values of A 1 , A 2 , h 0 , 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(Tsakiroglou et al., , 2009

by,
See Equation (17) bottom of the page where,

Accessibility functions
The accessibility functions can be described satisfactorily by the following functional form, where q 0 si is the initial fraction of sites occupied by the non-wetting fluid; q 0 si ¼ 0 for Hg intrusion (primary drainage, i = 0) and q 0 si ¼ q si for nitrogen desorption (secondary drainage initiating before the completion of adsorption, i = 1). Evidently, two different sets of parameters (a 0 , b 0 ) and (a 1 , b 1 ) are required to quantify the accessibility functions Y s0 (q b ) and Y s1 (q b ) that specify the Hg intrusion and N 2 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 d 2 Y si =dq 2 b À Á q b ¼q bci ¼ 0 and is obtained from the numerical solution of the following equation,

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, 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, where u 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, 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, 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, where u 2 , P 2 are the flow velocity and pressure at the outlet, L is the sample length, l g is the gas viscosity, and DP is the pressure drop across the sample. By combining equation (22) with equation (27), it is obtained,

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 , b s , and f s (r; <r s1 >, r s1 , <r s2 >, r s2 , c s ) were estimated with inverse modeling of N 2 adsorption isotherm. Then, the foregoing parameters were fixed, and f b (r; <r b1 >, r b1 , <r b2 >, r b2 , c b ) along with Y si (q b ) were estimated with inverse modeling of N 2 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 reestimated, with the simultaneous inverse modeling of N 2 adsorption/desorption isotherms and Hg intrusion curve.
3 Results and discussion

Inverse modeling of Hg intrusion curve and N 2 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 2 isotherms, and pay attention mostly on    (Fig. 4). For the three shale samples (Tab. 1), the Hg intrusion data were combined with N 2 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 (h Hg ) relationship was incorporated into the set of equations modeling Hg intrusion and N 2 sorption processes.
The application of the methodology of pore structure characterization from N 2 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 2 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 lm) 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 2 adsorption/desorption curves, and simultaneously the Hg intrusion curve. The initial guesses of TSD, PSD, and ACF are obtained from N 2 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 2 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 2 sorption to the estimation process (2 datasets against 1 dataset); (ii) the uncertainty of the "fixed values" of parameters included in N 2 adsorption/desorption equations (e.g. liquid N 2 contact angle); (iii) the lack of sufficient information concerning the large pore sizes that are not detected at all by N 2 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 h 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 , 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.

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) 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, k (nm) to the mean pore radius, r s (nm), and is given by,  (k L ) exp (nD) 3 Â 10 À4 61.0 61.0 6.6 (k L ) pred (nD) 0.048 36.8 (S1) 109.4 (S2) 0.39 (D kn ) exp = (bk L /l g ) exp (m 2 /s) 2.26 Â 10 À10 7.7 Â 10 À9 7.7 Â 10 À9 1.7 Â 10 À9 (D kn ) pred (m 2 /s) 1.3 Â 10 À10 6.7 Â 10 À9 (S1) 1.4 Â 10 À8 (S2) 9.5 Â 10 À10 k/<r s > 10.5 5.7 5.7 26.9 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.

Conclusion
Analytical models of Hg intrusion and N 2 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 2 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 2 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.