Regular Article
A methodology to predict the gas permeability parameters of tight reservoirs from nitrogen sorption isotherms and mercury porosimetry curves
^{1}
Foundation for Research and Technology Hellas – Institute of Chemical Engineering Sciences, Stadiou Str., Platani, 26504 Patras, Greece
^{2}
Department of Petroleum Engineering, Curtin University, WA 6151, Australia
^{*} Corresponding author: ctsakir@iceht.forth.gr
Received:
25
March
2020
Accepted:
16
March
2021
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 poreandthroat 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 (N_{2}) 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 N_{2} 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.
© C.D. Tsakiroglou et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 socalled “unconventional” features, the traditional approaches are unable to determine the reservoir permeability, and both the microscale 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 3dimensional 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 3Dimensional (3D) images obtained from micro/nano Computed Tomography (CT) and FIBSEM 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).
Threedimensional nanoporous 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 LatticeBoltzman 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 FIBSEM technique was employed to observe and characterize the morphology of the nanopores in shales, Lattice Boltzmann Method (LBM) was used to simulate highKnudsen gas flows, and the nanoscale 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 nonDarcy effect is evident (Cao et al., 2019). The imagebased 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 μm), interconnected by smallscale 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 poreandthroat network reconstruction from 3D SEM images with resolution of 4 nm was developed by accounting for complex geometry and extremely small sizes of the porebodies 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.
2 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 μm–3.5 nm). The liquid permeability and all pertinent Klinkenberg parameters were measured with a modified steadystate 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.
Properties of shale samples.
2.2 Porescale 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 poreandthroat networks (Fig. 1b). In the present work, the porosity of shales is regarded as a stochastic network of pores (sites)andthroats (bonds), the sizes of which are sampled from bimodal site f _{ s }(r) and bond f _{ b }(r) radius distributions, both composed of lognormal component functions, namely:
(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.
Fig. 1 (a) SEM image of AC1. (b) Representation of the pore structure as a poreandthroat network. 
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 ^{2}. On the other hand, β _{ s } = 1.0, implies a network of slitshaped pores with their volume proportional to r ^{1}. 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 poreandthroat 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 crosssection of each pore is modeled by a highorder 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,
Fig. 2 (a) Polygonal pore; (b) adsorbed layer along with condensed N_{2} along a triangular pore; (c) capillary equilibrium at a pore corner. 
2.3 Models of N_{2} adsorption/desorption in poreandthroat 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 porewalls, 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_{2} 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)
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 soproduced equations to any β _{ s } value. The liquid N_{2}saturation in a pore includes the adsorbed and condensed phases and is given by,
(7)where F _{LG} is the fraction of the pore crosssection area occupied by liquid N_{2}, and depending on the liquid/gas/solid contact angle, θ _{LG} (Fig. 2c; we assume that θ _{LG} = 0°, for liquid N_{2}, unless otherwise stated). If equation (7) is generalized for any β _{ 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,
(12b)which in combination with equations (5) and (6) lead to,
(13b)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 desorptiondrainage) or to isolated vapor pockets (secondary desorptiondrainage). 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,
(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_{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),
2.4 Model of Hg intrusion in poreandthroat 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 nonwetting 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^{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 θ _{Hg}, and the produced datasets (V _{ p } vs. θ _{Hg}) are fitted to a nonlinear parametric function of the form,
(16)so that the values of A _{1}, A _{2}, θ _{0}, w are estimated.
The Hg intrusion in a poreandthroat 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,
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 nonwetting 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 _{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 and is obtained from the numerical solution of the following equation,
2.6 Critical Path Analysis (CPA)
The transport properties of porous materials, characterized in the terms of a parametric poreandthroat 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,
(24)where φ is the porosity, t is the universal scaling exponent of conductivity (t = 2.0 for 3D 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 pressuredependent 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 steadystate gas flow through a porous medium, the integrated form of Darcy equation, including the Klinkenberg effect, takes the form,
(27)where u _{ 2 }, P _{ 2 } 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,
2.7 Demonstration of methodology
For the nonlinear parameter estimation, numerical codes in the software platform of Athena Visual Studio (Stewart and Caracotsios, 2008) were developed, and a stepwise 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_{2} 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_{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.
Fig. 3 Multistep methodology of pore structure characterization. 
3 Results and discussion
3.1 Inverse modeling of Hg intrusion curve and N_{2} sorption isotherms
It is wellknown that for low porosity samples, surface irregularities (e.g. nooks, crannies, crevices) and intergranular pores that are filled with Hg over the low pressure range are interpreted as macropores. It is of key importance to exclude such characteristics from the deconvolution of Hg intrusion and N_{2} 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).
Fig. 4 Corrected Hg intrusion curves for (a) AC1, (b) AC5, (c) AC8. 
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 }(θ _{Hg}) relationship was incorporated into the set of equations modeling Hg intrusion and N_{2} sorption processes.
Fig. 5 Correlation of the total pore volume with Hg contact angle by matching Hg intrusion with N_{2} desorption datasets (Tab. 2). 
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.
Estimated parameter values for shale samples AC1, AC5, AC8.
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 μ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_{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.
Fig. 6 Measured vs. predicted (a) N_{2} sorption isotherms, (b) Hg intrusion curve for sample AC1. 
Fig. 7 Measured vs. predicted (a) N_{2} sorption isotherms, (b) Hg intrusion curve for sample AC8. 
Fig. 8 Measured vs. predicted (a) N_{2} sorption isotherms (simulated structure S1), (b) N_{2} sorption isotherms (simulated structure S2), (c) Hg intrusion curve (simulated structures S1 & S2) for sample AC5. 
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 percolationtype 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.
Fig. 9 The pore and throatradius distributions obtained from inverse modeling for (a) AC1, (b) AC8, (c) AC5. 
3.2 Prediction of transport properties
The parameter values of simulated throatandpore 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 sopredicted 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).
Comparison of predicted vs. measured flow properties of shales.
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,
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 diffusionbased formulations must be applied. Looking at the postcalculated 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_{2} adsorptiondesorption 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 poreandthroat 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.
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 cofinanced by Greece and the European Union (European Regional Development Fund).
Appendix
Fig. A1 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for the simulated pore structures of sample AC1. 
Fig. A2 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for the simulated pore structures of sample AC8. 
Fig. A3 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for simulated pore structures (S1, S2) of sample AC5. 
References
 Al Hinai A., Rezaee R. (2015) Pore geometry in gas shale reservoirs, in: Fundamentals of gas shale reservoirs, Rezaee R. (ed), John Wiley & Sons Inc., pp. 89–116. [Google Scholar]
 Al Hinai A., Rezaee R., Saeedi A., Lenormand R. (2013) Permeability prediction from mercury injection capillary pressure: An example from the Perth Basin, Western Australia, APPEA J. 53, 31–36. [Google Scholar]
 Cao G., Lin M., Ji L., Jiang W., Yang M. (2019) Characterization of pore structures and gas transport characteristics of Longmaxi shale, Fuel 258, 116146. [Google Scholar]
 Chen L., Zhang L., Kang Q., Viswanathan H.S., Yao J., Tao W. (2015) Nanoscale simulation of shale transport properties using the lattice Boltzmann method: Permeability and diffusivity, Sci. Rep. 5, 8089. [PubMed] [Google Scholar]
 Firincioglu T., Blunt M.J., Zhou D. (1999) Threephase and wettability effects in triangular capillaries, Colloids Surf. 155, 259–276. [Google Scholar]
 Gregg S.J., Sing K.S.W. (1982) Adsorption, surface area and porosity, Academic Press, London. [Google Scholar]
 Hill D.G., Nelson C.R. (2000) Gas productive fractured shales: An overview and update, Gas TIPS 6, 3, 4–13. [Google Scholar]
 Hunt A.G. (2001) Applications of percolation theory to porous media with distributed local conductances, Adv. Water Resour. 24, 279–307. [Google Scholar]
 Javadpour F., Fisher D., Unsworth M. (2007) Nanoscale gas flow in shale gas sediments, J. Can. Pet. Technol. 46, 10, 55–61. [Google Scholar]
 Jones K.L., Laudone G.M., Matthews G.P. (2018) A multitechnique experimental and modeling study of the porous structure of IG110 and IG430 nuclear graphite, Carbon 128, 1–11. [Google Scholar]
 Katz A.J., Thompson A.H. (1986) Quantitative prediction of permeability in porous rock, Phys. Rev. B 34, 8179–8181. [Google Scholar]
 Song W., Wang D., Yao J., Li Y., Sun H., Yang Y., Zhang L. (2019) Multiscale imagebased fractal characteristic of shale pore structure with implication to accurate prediction of gas permeability, Fuel 241, 522–532. [Google Scholar]
 Stewart W.E., Caracotsios M. (2008) Computeraided modeling of reactive systems, John Wiley & Sons, Hoboken, NJ. [Google Scholar]
 Sun H., Yao J., Cao Y.C., Fan D.Y., Zhang L. (2017) Characterization of gas transport behaviors in shale gas and tight gas reservoirs by digital rock analysis, Int. J. Heat Mass Transfer 104, 227–239. [CrossRef] [Google Scholar]
 Tsakiroglou C.D. (2014) Computation of the twophase flow properties of intermediatewet porous media: A pore network approach, Can. J. Chem. Eng. 92, 515–523. [Google Scholar]
 Tsakiroglou C.D., Payatakes A.C. (1993) Porewall roughness as a fractal surface and theoretical simulation of mercury intrusion/retraction in porous media, J. Colloid Interface Sci. 159, 287–301. [Google Scholar]
 Tsakiroglou C.D., Ioannidis Μ.Α. (2008) Dual porosity modeling of the pore structure and transport properties of a contaminated soil, Eur. J. Soil Sci. 59, 744–761. [Google Scholar]
 Tsakiroglou C.D., Burganos V.N., Jacobsen J. (2004) Pore structure analysis by using nitrogen sorption and mercury intrusion data, AIChE J. 50, 489–510. [Google Scholar]
 Tsakiroglou C.D., Ioannidis M.A., Amirtharaj E., Vizika O. (2009) A new approach for the characterization of the pore structure of dual porosity rocks, Chem. Eng. Sci. 64, 847–859. [Google Scholar]
 Valvante P.H., Piri M., Lopez X., Blunt M.J. (2005) Predictive porescale modeling of single and multiphase flow, Transp. Porous Media 58, 23–41. [Google Scholar]
 Zhang P., Hu L., Meegoda J.N. (2017) Porescale simulation and sensitivity analysis of apparent gas permeability in shale matrix, Materials 10, 104, 1–13. [CrossRef] [Google Scholar]
 Zhang P., Lu S., Li J. (2019) Characterization of pore size distributions of shale oil reservoirs: A case study from Dongying sag, Bohai Bay basin, China, Mar. Petrol. Geol. 100, 297–308. [Google Scholar]
 Zheng J., Wang Z., Gong W., Ju Y., Wang M. (2017) Characterization of nanopore morphology of shale and its effects on gas permeability, J. Nat. Gas Sci. Eng. 47, 83–90. [Google Scholar]
 Zheng D., Wang W., Reza Z. (2019) Porenetwork extraction algorithm for shale accounting for geometryeffect, J. Pet. Sci. Eng. 176, 74–84. [Google Scholar]
 Zhou D., Blunt M., Orr F.M. (1997) Hydrocarbon drainage along corners of noncircular capillaries, J. Colloids Interface Sci. 187, 11–21. [Google Scholar]
 Ziarani A.S., Aguilera R. (2012) Knudsen’s permeability correction for tight porous media, Transp. Porous Media 91, 239–260. [Google Scholar]
All Tables
All Figures
Fig. 1 (a) SEM image of AC1. (b) Representation of the pore structure as a poreandthroat network. 

In the text 
Fig. 2 (a) Polygonal pore; (b) adsorbed layer along with condensed N_{2} along a triangular pore; (c) capillary equilibrium at a pore corner. 

In the text 
Fig. 3 Multistep methodology of pore structure characterization. 

In the text 
Fig. 4 Corrected Hg intrusion curves for (a) AC1, (b) AC5, (c) AC8. 

In the text 
Fig. 5 Correlation of the total pore volume with Hg contact angle by matching Hg intrusion with N_{2} desorption datasets (Tab. 2). 

In the text 
Fig. 6 Measured vs. predicted (a) N_{2} sorption isotherms, (b) Hg intrusion curve for sample AC1. 

In the text 
Fig. 7 Measured vs. predicted (a) N_{2} sorption isotherms, (b) Hg intrusion curve for sample AC8. 

In the text 
Fig. 8 Measured vs. predicted (a) N_{2} sorption isotherms (simulated structure S1), (b) N_{2} sorption isotherms (simulated structure S2), (c) Hg intrusion curve (simulated structures S1 & S2) for sample AC5. 

In the text 
Fig. 9 The pore and throatradius distributions obtained from inverse modeling for (a) AC1, (b) AC8, (c) AC5. 

In the text 
Fig. A1 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for the simulated pore structures of sample AC1. 

In the text 
Fig. A2 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for the simulated pore structures of sample AC8. 

In the text 
Fig. A3 (a) Liquid permeability, and (b) Knudsen diffusion coefficient of the critical path as a function of throat radius for simulated pore structures (S1, S2) of sample AC5. 

In the text 