A Triple Pore Network Model (T-PNM) for Gas Flow Simulation in Fractured, Micro-porous and Meso-porous Media

In this study, a novel triple pore network model (T-PNM) is introduced which is composed of a single pore network model (PNM) coupled to fractures and micro-porosities. We use two stages of the watershed segmentation algorithm to extract the required data from semi-real micro-tomography images of porous material and build a structural network composed of three conductive elements: meso-pores, micro-pores, and fractures. Gas and liquid flow are simulated on the extracted networks and the calculated permeabilities are compared with dual pore network models (D-PNM) as well as the analytical solutions. It is found that the processes which are more sensitive to the surface features of material, should be simulated using a T-PNM that considers the effect of micro-porosities on overall process of flow in tight pores. We found that, for gas flow in tight pores where the close contact of gas with the surface of solid walls makes Knudsen diffusion and gas slippage significant, T-PNM provides more accurate solution compared to D-PNM. Within the tested range of operational conditions, we recorded between 10 and 50% relative error in gas permeabilities of carbonate porous rocks if micro-porosities are dismissed in the presence of fractures.


Single Pore Network Model
A pore network model (PNM) is a simplified proxy model to simulate various transport phenomena in porous materials. PNMs can be statistically generated (Raoof and Hassanizadeh 2010;Karadimitriou et al. 2012) or extracted from realistic tomography images of porous media (Dong and Blunt 2009;Rabbani et al. 2014;Gostick 2017). Despite PNM simplifications, it is still interesting due to the low computational cost, especially with the growing size and resolution of the tomography images of porous materials . Single PNM has been initially introduced by Fatt (1956) (Fig. 1a). He discussed several approaches to tackle the porous media transport problems from single-size capillary tube to dynamic displacement of two-phase fluids. Since then many developments  (Blunt 2001;Fatt 1956), b D-PNM with meso-pores coupled to a lattice micro-pores (Bekri et al. 2005), c dual with meso-pores and parallel micro-pores and throats (Bauer et al. 2012), d a fully twoscale PNM with micro-and meso-pores (Jiang et al. 2013), e dual with extra micro-throats to connect the meso-pores (Bultreys et al. 2015), f dual with meso-pores and fractures with variable aperture (Hughes and Blunt 2001), g dual with meso-pores and fractures with fixed aperture (Erzeybek and Akin 2008), h D-PNM with extra bypassing fracture links (Liu et al. 2017), i dual with unstructured network and variable geometry of fracture (Jiang et al. 2017), j the presented T-PNM which includes meso-and micro-pores, meso-and micro-throats, fracture nodes, and fracture links have been accomplished by researchers to enhance the capability and reliability of PNMs (Baychev et al. 2019) to realistically simulate and predict the properties of porous materials (Garboczi 1990;Joekar-Niasar and Hassanizadeh 2012;Xiong et al. 2016), from hydraulic permeability (Dong and Blunt 2009) and mass diffusivity (Burganos and Sotirchos 1987;Erfani et al. 2019) to electrical (Friedman and Seaton 1998) and thermal conductivities (Plourde and Prat 2003).
In the present study, the main focus of PNM is on simulation of gas transport in porous media. One of the first applications of PNM for gas flow simulation was presented by Millington and Quirk (1961). They used the physical concept of capillaries network demonstrated by Fatt (1956) and studied the effects of Knudsen diffusion (Wakao and Smith 1962) on gas permeability of the porous solids. From 1960s to 1990s, most of the gas PNM studies were dedicated to chemical processes and catalysis applications (Beeckman and Froment 1979;Webb et al. 1991;Coppens and Froment 1995), while gradually by rising the role of the natural gas in world energy market (Victor et al. 2006), many petroleumoriented studies emerged in the literature (Wang and Reed 2009;Javadpour et al. 2007;Zhang et al. 2015). Specifically, in the past decade, the prosperity of unconventional natural gas resources such as shale gas, tight gas, and coal-bed methane triggered a significant research trend to model and simulate gas transport in these irregular materials (Javadpour 2009;Ghanizadeh et al. 2015;Zhang et al. 2015;Naraghi and Javadpour 2015;Chen 2016;Wang et al. 2018;Tahmasebi 2018;Song et al. 2019). As a prominent example, shale gas is pioneering among the unconventional gas research subjects and many numerical studies have been recently presented to include a wide range of physics into a gas transport model such as adsorption (Chai et al. 2019;Yang et al. 2019), multiple gas components (Bai et al. 2019), multi-level of porosity (Chen et al. 2018;Ma et al. 2014), super critical condition (Chen et al. 2018), and arbitrary pore shapes (Wu et al. 2015). Furthermore, in parallel to the natural gas PNM studies, many fuel cell researchers started to use PNMs as a tool to model gas diffusion in porous membrane of proton exchange cells (Gostick et al. 2007;Sinha and Wang 2007) even with capability of including micro-porosities in the network (Sadeghi et al. 2017) which is one of the main focus of the present study.

Meso-pore and Micro-pore Coupling
A single continuum pore network model defines a void-solid structure in which the granular portion of the material does not contribute to the fluid flow (Blunt 2001). In the more recent pore network studies, many researchers have started to consider an implicit contribution for the solid phase assumed in the transport properties of porous media (Tahmasebi and Kamrava 2018;Bekri et al. 2005;Bauer et al. 2012Bauer et al. , 2011Khan et al. 2019;Bultreys et al. 2015;Jiang et al. 2013;de Vries et al. 2017). The main logic behind this assumption is the fact that in the heterogeneous material, micro-pores residing in the solid section can affect the macroscopic properties of porous material in certain physical conditions. In some physical processes such as gas transport and two-phase imbibition, micro-pores can make a statistically significant change in the simulation results as well as in the experiments (Javadpour 2009;Bultreys et al. 2015).
Many geometrically different structures for dual pore network of meso-and micro-pores have been suggested by researchers. Bekri et al. (2005) presented a D-PNM with latticelike networks of micro-pores embedded between the neighbouring meso-pores (Fig. 1b).
Then, they predicted relative permeabilities with a quasi-static assumption and compared them with experimental values. In another effort, Bauer et al. (2012) considered that some portion of the meso-pores are connected to the adjacent pores with both meso-and microthroats in a paralleled way (Fig. 1c). Using an approach similar to Bekri et al. (2005), they compared their relative permeability results obtained from modelling and experiment. Jiang et al. (2013) presented the next comprehensive work on D-PNMs and developed a fully connected PNM with two scales of pores and throats to consider the effect of unresolved porosities on PNMs constructed based on micro-tomography images. They studied the effects of changing the density of micro-throats on two-phase properties of porous media such relative permeability. Inspired by their work, Bultreys et al. (2015) used micro-CT images of several tight and heterogeneous porous carbonate rocks to extract a pore network that in addition to the meso-throats, utilizes some micro-throats to connect adjacent meso-pores (Fig. 1e). Some of their non-fractured raw images of porous carbonated rocks with multi-levels of porosity are used in this paper as case studies.

Meso-pore and Fracture Coupling
Similar to the smaller-scale features such as micro-porosities that can be plugged into pore networks, some larger-scale features like fractures can be coupled to regular PNMs. Researchers must answer two important questions when coupling these components: how should they specify a void section as a fracture or simply a pore body, and what network properties, such as hydraulic permeability, should those sections exhibit? To answer the second question, we usually assume that a pore throat's hydraulic permeability obeys Hagen-Poiseuille law (Schmid and Henningson 1994;Mortensen et al. 2005), and equals r 2 ∕8 , where r is the throat hydraulic radius. Similarly, we can calculate the absolute permeability of an ideal fracture segment using h 2 ∕12 where h is the distance between two planes of the fracture (Witherspoon et al. 1980;Zimmerman and Bodvarsson 1996). Hughes and Blunt (2001) presented one of the first studies to simulate fluid flow through the fractures using PNM. They arranged several box-shaped conductors in a row to simulate the aperture variation of a fracture which is in contact with a regular porous space (Fig. 1f). Erzeybek and Akin (2008) studied different arrangements of such fracture nodes in different pore and fracture configurations to evaluate the contribution of each element in the fluid flow simulations. Jiang et al. (2017) presented a more versatile approach to calculate the absolute permeability of a pore-fracture system. They utilized the medial axis method to extract the regular PNM then devised a shrinking algorithm to locate the fractures, answering the important question of how to specify a void section as a fracture or a pore body. In this shrinking algorithm, they removed non-planar structures in a stepwise manner to obtain the remaining fracture body and simulated it with a virtual network of locally variable throats (Fig. 1i). As mentioned, this virtual network for simulating fracture behaviour has been introduced by Hughes and Blunt (2001), and it is used in the present study with modifications. Details of fracture modelling and method validation will be discussed later in this paper.
In another effort to couple fractures and porous space, Weishaupt et al. (2019) used a Navier-Stokes model for free flow and a PNM for the porous domain to simulate pressure drop and velocity map of a porous-fractured system in two dimensions. They studied the mass transport of a component in the porous space due to the free flow in fractures for different Reynolds numbers. Their developed method is dynamic, robust and flexible for modelling of unstructured networks coupled with fractures; however, it could become computationally intensive for complex 3D geometries.

Triple Pore Network Model
To the best of our knowledge, there is no previous study to couple micro-pores, mesopores, and fractures at the same time in a triple pore network model (T-PNM). Here, we present an approach to extract a triple pore network model based on the semi-real microtomography images of heterogeneous porous material. Then by simulating liquid and gas flow through the porous-fractured geometries of the selected samples, we demonstrate the capability of the proposed method and justify the conditions that a T-PNM would give more accurate results than single PNMs or D-PNMs. The diagram of the previous D-PNM structures in the literature and the presented triple structure are illustrated in Fig. 1. Some of the illustrations in this figure are inspired by Tahmasebi and Kamrava (2018).

Methodology
In this research, we extract a T-PNM from micro-tomography images of heterogeneous porous media and simulate steady-state liquid and gas flow through the network to investigate the viability/applicability/advantages of the proposed method. For this purpose, we describe the porous material preparation and image processing techniques that lead to extracting a PNM consisting of (1) meso-pores, (2) micro-pores which are unresolved in the images, and (3) fractures which have been generated synthetically on the volumetric images of porous carbonated rocks. We define micro-pores as pores so small that they are not explicitly visible in the tomography images. However, a population of the micro-pores causes local density reduction that locally decreases the intensity of the X-ray tomography image (Taud et al. 2005). On the other hand, based on our definition, meso-pores are clearly visible in the images, and pore boundaries can be distinct from the solid or partially solid background. After extracting the regular PNM, we need to include fractures and micro-pores. Finally, we clarify the equations that we use to simulate steady-state liquid and gas flow through the proposed T-PNM by assuming no surface reaction or physical adsorption.

Materials
The International Union of Pure and Applied Chemistry (IUPAC) defines meso-pores as between 2 and 50 nm in size (Everett 1972), and states that size of micro-pores should not exceed 2 nm. However, in the present study, we assume a dynamic threshold level for distinguishing between the micro-pores and meso-pores. This dynamic threshold is equal to the imaging spatial resolution (see Table 1), in which pores with sizes less than spatial resolution are considered as micro-pores and larger ones are assumed to be meso-pores.
Three carbonate porous rocks are selected as case studies. These samples are heterogeneous with two scales of porosity imaged and published by Bultreys et al. (2015) in the public domain. Several artificial fractures are created though these samples to mimic the conditions of a triple-pore system. These fractures are carved into the samples with a costminimizing random walk approach inspired by Mhiri et al. (2015) and coded for our specific application. The method dictates the fracture pathway to grow in the direction with highest porosity and consequently with less mechanical strength. This approach is physically plausible, considering the fact that breaking a material at loose connections requires less amount of energy (Alaboodi and Sivasankaran 2018). Details of generating fracture realizations are beyond the scope of this study; however, more information is provided in "Appendix 1". Estaillades, Savonnieres and Massangis are the conventional names of three types of carbonated porous rocks that we use as the raw material ( Fig. 2a-c, respectively). Dimensions and spatial resolution of the porous samples studied are provided in Table 1. Figure 2 illustrates the 3D map of CT numbers for the porous-fractured samples studied and the corresponding porosity map of each sample. We assume that mineralogy and consequently mineral density remain constant throughout various sections of the samples. Based on this assumption, bulk density of a section is the major factor that controls the energy of the photons received during the tomography process (Bryant et al. 2012;Lagravère et al. 2006). Bulk density of the pure materials is inversely related to porosity (Athy 1930). Consequently, we can roughly assume that any increase in porosity, linearly decreases the local CT number of the material (Vega et al. 2014). Therefore, in order to obtain a porosity map for the unresolved part of the image we simply consider a linear relationship with a unit slope between the normalized CT number and micro-porosity. Then, the darkest portion of the CT number map is considered to be completely void and the lightest parts of the geometry are assumed to be completely solid.
The segmentation between the void space and porous space of the samples is conducted using Otsu algorithm (Vala and Baxi 2013) which is based on the minimization of standard deviation in each segment of the data. So, in Fig. 2d-f, the pure yellow segments of the samples are assumed to be completely porous and the porosity of other parts of the samples varies between 0 and 1. In addition to tomography images, an experimental curve for the pore size distribution of each porous material has been obtained by mercury intrusion porosimetry. These curves are adopted from Bultreys et al. (2015), Gibeaux et al. (2018) and Neveux et al. (2014), respectively, for Estaillades, Savonnieres and Massangis samples. Dashed part of the curves in Fig. 2g represent the size distribution of the sub-resolution pores.

Extraction of the Single Pore Network Model
Extraction of a single-continuum pore network model is a routine procedure reported in the literature of modeling fluid flow in porous materials (Xiong et al. 2016). In order to extract the regular PNM, we use the watershed segmentation algorithm (Mangan and Whitaker 1999) which has been extensively employed to segment objects with morphological connections or overlapping (Wildenschild and Sheppard 2013;Gostick 2017;Rabbani and Ayatollahi 2015). It has been shown in the literature that watershed segmentation is reasonably accurate and computationally efficient for extraction of the pore bodies and pore  (Gostick 2017;Wildenschild and Sheppard 2013). This algorithm takes the distance transform of the pore space to locate the narrowest pathways of pore network and record them as pore throats (Rabbani et al. 2014;Sheppard et al. 2005;Arns et al. 2007). Distance transform gives the minimum distance of a zero voxel to the nearest non-zero voxel (Breu et al. 1995) (Fig. 3b). This transform is a function available in most of the image processing libraries and packages. We used MATLAB 2018a Image Processing Toolbox to analyse 2D and 3D images. In order to avoid over-segmentation of the structure, it is recommended to apply a smoothing transform such as Gaussian filtering or image opening on the distance map (Gostick 2017;Baychev et al. 2019;. Figure 3 illustrates a simple workflow for extracting single PNMs (a-d), D-PNMs (e-h), and T-PNMs (i-l). Initially, we consider a simple bed of spherical solids (Fig. 3a). Applying the distance transform on this geometry generates a map visualized in Fig. 3b. After smoothing, watershed algorithm uses this map to generate the distinct map of the pore bodies labelled with random shades of grey in Fig. 3c. This algorithm simulates a hypothetical flooding process starts from the local maximum values of the distance map. In a stepwise  Mangan and Whitaker 1999;Angulo and Jeulin 2007). Watershed segmentation inherently calculates the labelled map of pore bodies. In Fig. 3c, we have assigned random shade of grey to each pore body to demonstrate the segmentation. Now, using a 3 × 3 sliding window we scan the whole area of the image and detect connections between different pore labels. Consequently, a pore network can be constructed based on the interpore connections (Fig. 3d).
To simulate fluid flow within the extracted single PNM, we should determine the magnitude of pressure drop occurring when fluid is moving between a pair of connected pore bodies. If we ignore the pore body pressure drop, the major element that controls the fluid transport will be the pore throat that has minimal opening area compared to its connected pores. Hydraulic permeability in a throat can be simply calculated by Hagen-Poiseuille law (Schmid and Henningson 1994), and it is equal to r 2 ∕8 that r is the throat hydraulic radius. A more realistic approach to calculate the absolute permeability of a throat with an arbitrary cross-sectional shape is presented by Rabbani and Babaei (2019) and it is based on the Lattice-Boltzmann method (LBM) to simulate a steady-state single-phase flow through the throat. They showed that numerical results obtained for permeability of arbitraryshaped throats are highly correlated with averaged distance map of the throat image and permeability can be obtained as the following quadratic form: where K p is the absolute permeability of the porous space links such as a meso-throat with pixel 2 unit, and D is the mean distance of the throat cross-sectional image. This approach for obtaining absolute permeability of a single throat is referred to as image-based throat permeability model (ITPM) (Rabbani and Babaei 2019) and its code is available in the public domain. 1 As mentioned previously, a distance map shows the minimum distance between a zero voxel to the nearest non-zero voxel (Fig. 3b) and to obtain D , we perform an averaging on the non-zero values of a distance map. More details of statistical and physical justifications for this empirical equation are provided in Rabbani and Babaei (2019). Also, the basics and modelling assumptions of LBM simulation used in Rabbani and Babaei (2019) and extended to the current study are described in "Appendix 2". By describing the hydraulic behaviour of the links in a single continuum PNM, we are able to model fluid flow through the extracted networks that consist of meso-pores and meso-throats.

Fracture Inclusion
The common approach to model a D-PNM containing meso-pores and fractures is to separate the single PNM from fracture elements and then use different flow equations for each type of the elements (Jiang et al. 2017). Various morphological operations can be used to specify the location of fracture elements (Jiang et al. 2017); however, this distinction could become difficult in cases with complicated morphologies. For instance, a meso-throat can have a fracture-like elongated cross-sectional shape that casts shadow on the distinct definition of the elements. In order to avoid this issue, we present a unified permeability model for both meso-scale and fracture networks. A fracture network is assumed to be composed of an interconnected structure of nodes and links. As a similar definition to the meso-scale network, fracture links represent geometrical bottlenecks between fracture nodes so that these bottlenecks create a pressure drop when fluid is moving between two connected nodes.

Shape Issues
In this subsection, the capability of ITPM (Eq. 1) to calculate the absolute permeability of meso-throats as well as fracture links will be investigated. In this regard, we generate a set of cross-sectional images with a wide range of elongations from 1 to 10 and different roundnesses from 0.1 to 1 (Fig. 4). Elongation is the ratio between the major axis and the minor axis of the shape, and roundness is a measure that shows the deviation of the shape radii from a perfect circle, thus a shape with pointy corners will have less roundness (Diepenbroek et al. 1992). By analysing the sets of geometries, we observed that the hydraulic permeability values obtained by LBM simulation and ITPM are consistent and averaged relative error is around 3.3%. Relative error percentage does not necessarily increase or decrease when objects are more elongated. This observation is validated by performing Student's t-test between the relative errors and elongation for 20 different cross sections illustrated in Fig. 4. The obtained Student's t-test p-value is 1.5 × 10 −4 that verifies the independence of these two variables.

Lateral Segments
Another issue when we are extracting a fracture-included porous material is unfavourable lateral segments. Consider a channel cross section as visualized in Fig. 5a with flow direction perpendicular to the image surface. If our pore segmentation method mistakenly divides this elongated channel into 3 sections, its permeability is likely to be underestimated due to the added unrealistic zero-velocity boundaries as shown in Fig. 5b. Considering that we are using averaged distance values of throat images ( D ) for calculating their permeability (Eq. 1), Fig. 5b exhibits a lower permeability compared to Fig. 5a. The simple solution to avoid the aforementioned issue is to calculate the distance map prior to the pore network segmentation. Thus, the distance map will not be affected in the case that lateral segmentation occurs (Fig. 5c).

Coupling to the Meso-pore Network
We showed that ITPM can be used for permeability estimation of both meso-throats and fracture links. So these elements are effortlessly coupled to each other with no extra calculation to distinguish each of them. Figure 3e-h illustrates the steps of network extraction process for our defined D-PNMs and as it can be seen, the workflow is similar to the single PNMs.

Micro-porosity Inclusion
As described, the solid elements of the porous material could contain micro-porosities. As mentioned in Sect. 2.1, we are able to generate an approximated porosity map for micropores based on the mineralogy and normalized CT numbers. In this section, we describe the hydraulic properties of partially solid elements by knowing their approximated porosity and pore size distribution. Our proposed configuration of T-PNM is illustrated in Fig. 6a in which two micro-pores are connected to a pair of meso-pores and one fracture node. We define two types of micro-throats: (1) micro-throats that connect two micro-pores (Fig. 6c), and (2) micro-throats that connect a micro-pore into a meso-pore or fracture node ( Fig. 6b and d). Each micro-throat is consisted of a bundle of micro-tubes with different range of radii. These micro-tubes intersect at the centre of the solid element that is called main micro-pore ( Fig. 6b-2). Also, some isolated cavities could exist within the solid element and they affect the image-based porosity we obtain from CT numbers while not contributing to fluid flow ( Fig. 6b-1). If we assume that the radial density of micro-tubes is constant, based on a probabilistic permeability model of a bundle of tubes developed by Juang and Holtz (1986), the absolute radial permeability of the solid elements ( K s ) are: where m is micro-porosity of the solid element, r m is the radius of the micro-tubes with a distribution described by f (r m ) , f (r m ) is a function that gives the void fraction of the solid element occupied by micro-tubes with the radius of r m . Parameter is porosity correction factor between 0 and 1 that subtracts the portion of the total micro-porosity occupied by isolated pores. In other words, when is 0 it means that we have no interconnected microporosity in the solid element and when it is 1, it means there is no isolated micro-porosity within the solid element. This parameter can be obtained from nano-tomography, highresolution FIB-SEM images (Ma et al. 2018), or by density analysis of the crushed particles of the porous sample (Luffel and Guidry 1992). Here, we have simply assumed that is 0.95. The term Meso-Pore Meso-Throat Micro-Pore Micro-Throat Frac-Node Frac-Link Solid Fig. 6 Different configurations of micro-porosities in connection with other elements in a T-PNM, a a system of two micro-pores in connection with two meso-pores and one fracture node, b the structure of a micro-throat type 2 in connection with a fracture node, c the structure of a micro-throat type 1, d structure of a micro-throat type 2 in connection with a meso-pore (details of figure elements are described in Table 2) et al. 2005). The integral used in Eq. 2 can be simply solved by numerical discretization.
Considering that the overall pore size distribution of samples is experimentally obtained by mercury intrusion, we only need to extract and normalize the sub-resolution part of the pore size curve ( f (r m ) ) (Fig. 2g). Then by descretizing f (r m ) over different values of r m , the integral will be converted to a set of summations. Then, to calculate the absolute permeability of the micro-throats type 1 (Fig. 6c), we use Eq. 2 but instead of a single value for micro-porosity ( m ), we put the averaged value of micro-porosities from two adjacent solid elements. This averaging needs to be weighted relative to the radii of the elements, which means that the larger size of the solid elements makes them more significant in permeability averaging for micro-throats type 1. Additionally, in order to calculate the absolute permeability of the micro-throats type 2, we ignore the pressure drop happening in the void section of the micro-throat (Fig. 6b-4 and Fig. 6d-4). Considering that the absolute permeability of a series of elements is mainly controlled by the smallest permeability value (Ahmed 2018), this assumption is reasonable.

Coupling to the Dual Pore Network Model
In the previous sections, we developed a D-PNM consisted of meso-pores, and fractures. Now, in order to couple micro-pores into that D-PNM, we need to revise the pore network extraction method. In order to clarify the workflow of T-PNM extraction, Fig. 3i-l is presented. Fig. 3i illustrates the porosity map of a system with three types of porosities: micro-pores, meso-pores and fractures. In the next step, we apply Euclidean distance transform on both void and solid spaces separately. Then, the obtained distance maps are summed as follow: That DM is the overall Euclidean distance map (Fig. 3j), DT is the distance transform function that takes a binary map and gives back its distance map, and PM is porosity map of the media as appears in Fig. 3i. Consequently, PM < 1 implies a binary map in which all the pixels with porosities less than 1 are set to one and other pixels are set to zero. Similarly PM = 1 represents a binary map in which all the pixels that are completely porous, are set to 1 and the rest of the pixels are set to zero. On the basis of the obtained overall distance map, we apply watershed segmentation and generate the node map that includes all three types of aforementioned nodes (Fig. 3k). Finally, based on the generated node map, we detect the neighbour nodes and interfaces between them. Then by determining the types of micro-throats, we assign proper hydraulic properties to them. Finally, the T-PNM could be visualized as shown in Fig. 3l.

Results and Discussion
In this section, we first verify the developed method using simplified geometries. Then, statistics of the T-PNMs will be compared to D-PNMs in order to provide a more quantitative measure of differences. Next, steady-state gas flow (formulation presented in "Appendix 3") through the extracted PNMs will be simulated to investigate the effects of micro-porosity presence on the governing flow mechanisms. Finally, the results of gas flow simulation will be extended to a wide range of pressure, temperature, and molecular mass  Interface area between a micro-pore and a fracture node in which a bundle of micro-tubes with variable sizes represents a type 2 micro-throat 2 Interface area between two micro-pores in which a bundle of micro-tubes with variable sizes represent a type 1 micro-throat 3 Interface area between a micro-pore and a meso-pore in which a bundle of micro-tubes with variable sizes represents a type 2 micro-throat   The length of the micro-tubes that form the micro-throat type 1 Figure 6d 1-6 Same as Fig. 6b-1-6 conditions and the observed trends in the apparent permeability of the samples will be discussed.

Method Verification
In order to measure accuracy of the proposed methodology, we have constructed six idealized geometries that contain a combination of the three porous elements: meso-pores, micro-pores and fractures (Fig. 7a). Then, absolute permeability of each geometry is calculated using the developed methodology and analytical equations. The analytical absolute permeability of the geometries is calculated based on expressions r 2 ∕8 and h 2 ∕12 , respectively, for permeabilities of capillary tubes and ideal fractures as discussed before. Also, analytical permeability of the partially solid section is obtained by assuming ideal microtubes with an identical radius. The size of all cubic geometries is 200 3 μm and spatial resolution is 1 μm per voxel. Four of the geometries (Fig. 7a-2, a-4, a-5 and a-6) contain six capillary tubes with the diameter of 20 μm , equally distributed and aligned to the direction of flow. Similarly, in four of the geometries ( Fig. 7a-1, a-3, a-5 and a-6), an ideal fracture with no roughness is present in the middle section with the aperture of 10 μm . Additionally, in geometries a-3, a-4, and a-6, we have considered 30% of micro-porosity in the shape of parallel micro-tubes with diameter of 2 μm (Fig. 7d).
Watershed segmentation cannot be used to extract the pore network of these simplified geometries due to their idealized shapes. The reason is that, most of the pore network extraction methods including watershed is useful to detect the pathways with minimal cross-sectional area, while in an idealized uniform geometry, there is no minimal crosssectional area in the capillary tubes nor in the fracture. Consequently, to generate the node map we have divided the void-space and solid-space into equal cubic segments with the size of 10 μm (Fig. 7c) and then we have superimposed these maps to generate a complete node map which includes all flow domains. As can be seen in the magnified part of Fig. 7c, each node is labelled with a random colour.
Considering the porous elements mentioned in this paper, we investigate six types of geometries that include a combination of the different elements, and then using the developed pore network model, we calculate the absolute or liquid permeability of each geometry as presented in Fig. 7a. We assumed 1 Pa pressure gradient between two opposite faces of the PNM to calculate the absolute permeability of the system. The average relative error between the estimated permeability using the present method and the analytical solution is less than 2%. Figure 7b illustrates the pore pressure obtained by simulating single-phase steady-state flow on the geometry with triple porosity.
As can be seen in the magnified section of Fig. 7c, capillary tubes contain lateral segments. These segments did not affect the permeability calculations (Fig. 7a), because we have calculated the distance transform of the geometry prior to the segmentation. To calculate the analytical permeability of the whole geometry, we assume a system of parallel conductors with the same length. Thus, the equivalent analytical permeability of the whole geometry ( K e ) is Ahmed (2018): where s is the volume fraction of the geometry occupied by the partially or fully solid elements, K f is the absolute permeability of the fracture, f is the volume fraction of the geometry occupied by the fracture, and K p is the overall permeability of the six capillary tubes.

Network Statistics
In this subsection, we discuss the statistics of the extracted triple pore networks and compare them with D-PNMs that contain only meso-pores and fractures. In other words, the effects of including micro-porosities on the statistical features of the PNMs are investigated. Verification of estimated absolute permeabilities on a simplified geometry, a comparison of absolute permeabilities obtained by the present model and analytical solution for six combinations of porosity types with an average relative error less than 2%, b pore pressure distribution in a single-phase steady-state flow modelling with 1 Pa pressure difference between two opposite faces of the T-PNM, c extracted node map with a random colour for each node, d original geometry of the medium with a fracture, six capillary tubes, and 0.3 porosity at the solid sections Some statistics of these networks, including network porosity, average link permeability, specific surface, number of links, and number of nodes are presented in Table 3. As it can be seen, T-PNM shows higher porosity up to two times higher than the D-PNM which does not contain micro-porosity. On the contrary, the average link permeability in a T-PNM is commonly four to five times lower than the D-PNM due to the effect of multitudes of low permeable micro-throats which are present in a T-PNM. This ratio is virtually the same as the number of links of the T-PNMs compared to the D-PNMs. It is notable that microthroats can contain up to thousands of micro-tubes with minuscule radii. Therefore, the actual number of links can be two or three orders of magnitude larger than the reported values. Specific surface of the pore network is defined as the ratio between total internal surface area and the bulk volume of the whole geometry. As it can be seen in Table 3, specific surface of the T-PNM can be up to 20-40 times greater than a D-PNM. This property is significant in the simulation of the surface processes such as reactions or diffusion layers (Peigney et al. 2001;Bacsa et al. 2000).
Additionally, distributions of two pore network parameters are illustrated in Fig. 8 for the triple and dual approaches. Figure 8a shows the distribution of the node connectivity in both networks of Estaillades sample. T-PNM has considerably higher coordination number due to the presence of partially solid elements which are in touch with many adjacent nodes and increase the maximum number of connections per node up to 50. The other effect of micro-porosities is shown in Fig. 8b that presents the distribution of node surface area of both PNMs. The presence of partially solid elements exhibit nodes with surface area with two to three orders of magnitudes larger than the regular nodes. Consequently, simulation of the surface dependant processes can give significantly different results by using a T-PNM instead of a dual or single PNM.
In order to provide more insight into the physical features of the extracted PNMs, Fig. 9 is presented for Estaillades sample. Figure 9a illustrates the top view of the original volumetric image in which the level of darkness intuitively indicates the porosity. Figure 9b and c show two different sides of the PNM, respectively, with dual and triple architecture, stacked to each other for better comparison. The T-PNM half is denser and contains many large size nodes that indicate the presence of the large partially solid elements, while the D-PNM half is less dense with relatively smaller size of nodes. Large nodes of the triple side are highly connected to the neighbouring nodes and this is the source of high coordination number as shown in Fig. 8a. In terms of the surface area, we expect to see large values for partially solid elements ( Fig. 9d and e) due to the presence of hundreds to thousands of micro-tubes in each of the links. Finally, if we calculate the absolute permeability of each link based on Eqs. 1 and 2, the distribution of this parameter is visualized in Fig. 9f and g. As it can be seen, despite the large sizes of the partially solid nodes, their connected links have low permeability. However, due to the large number of connections, their role can become significant in some cases. Roughly, it can be stated that the highest absolute permeabilities of both PNMs belong to fractures by comparing the original geometry (Fig. 9a) and permeability network of Fig. 9f and g.

Gas Flow Simulation
In this subsection, we present the gas flow simulation results (formulation presented in "Appendix 3"). Initially, effective flow mechanisms are investigated and then flow pattern is discussed among D-PNMs and T-PNMs. Finally, we present a sensitivity analysis for conditions with different pressure and temperature, followed by a discussion on the effect of gas molecular mass.

Flow Mechanisms
As discussed in "Appendix 3", three main types of gas transport mechanisms happen in porous media: Knudsen, slippage and viscous flow. Each of these mechanisms could become significant in a specific range of Knudsen number. This sensitivity of gas flux ratio to Knudsen number is plotted for each type of the mechanisms and their combinations in Fig. 10a. In order to calculate the gas flux fractions, we have assumed an ideal tube geometry and then by gradual changing of the tube radius, flux fractions are obtained in a wide range of Knudsen number. In Fig. 10a and Knudsen fluxes, and finally, non-slip flux is the summation of the Knudsen and viscous fluxes. Assuming the definition of the Knudsen number which is the ratio of the gas mean free path to the tube radius, it is reasonable that the viscous flux fraction decreases in higher Knudsen numbers. Under conditions where gas molecules frequently collide with tube walls rather than other gas molecules. When Knudsen number is greater than 1, Knudsen flux and slippage are dominant flow mechanisms. On the other hand, when Knudsen number is smaller than 0.01, non-viscous fluxes are minimal. It is noteworthy that the fraction of the slip-flux is always less than the Knudsen flux at any Knudsen number range. In order to investigate the effective flow mechanisms in a more complex geometry, we have considered four sets of gas flow mechanisms to simulate the full pore network model of three porous materials. Figure 10b shows the ratio between the triple and dual gas permeabilities calculated with different mechanistical assumptions for three rock samples at the standard temperature-pressure condition. Simulated gas is methane and we have considered 1 Pa pressure difference between two opposite faces of the pore network model in the x direction. As readers can see in Fig. 10b, the ratio between triple-and dual-porosity permeabilities increases as we integrate more flow mechanisms into our model. The final columns indicates the condition that all three mechanisms have been applied, and it is found that triple permeability can be 1.3-1.55 times larger than the dual permeability which ignores the micro-porosities. Also, for the simulated conditions, it can be concluded that Knudsen flux, is more significant compared to slippage, since the major jump in permeability ratio happens when Knudsen flow is added to the system (Fig. 10b).

Flow Pattern
Considering a steady-state condition we have simulated gas flow in the direction of Y axis through D-PNM and T-PNM of Estaillades sample. The pressure map and flow rate distribution are visualized in Fig. 11. Pressure difference between two opposite faces of the PNM is 1 Pa and we have assumed that in this short range of pressure, the gas density and viscosity remains constant. The average pressure for calculating gas properties is 101 kPa and the temperature is 25 • C . As it can be inferred from Fig. 11a and b, iso-pressure lines does not show sharp changes when changing the geometry from a T-PNM to a D-PNM. Figure 11c and d illustrates the flow rates of links within D-PNM and T-PNM partly imposed on each other. Despite the pressure distribution, flow pattern is significantly different between two PNMs and many micro-throats are contributing in fluid transport. Due to the numerous connections, micro-throats make a difference in the overall permeability results, while each single of them does not pass a significant amount of flow compared to the meso-throats or fractures.

Sensitivity Analysis
In order to provide an insight into the application of a T-PNM and clarify the conditions that ignoring the partially porous part of the media could create a significant error, we analyse the sensitivity of gas permeability by changing three independent parameters: average pressure, temperature, and gas molecular mass. We simulate gas flow in D-PNMs and T-PNMs of all three porous samples for a wide range of each parameter by keeping the other two as constant. The results are presented in Fig. 12. Figure 12a shows the gas permeability changes in different average pressures for methane gas and constant temperature of 25 • C . Low pressures could dramatically increase the gas permeability of T-PNMs, while D-PNMs are relatively less sensitive to the pressure change. For the tested range of average pressure, lowering the pressure two orders of magnitude increases the gas permeability of T-PNMs up to two times. This is mainly due to the increment of gas mean free path in low pressures and consequently higher Knudsen flow rates. Gas permeability is relatively less sensitive to the temperature variations compared to the pressure (Fig. 12b).
By increasing the average temperature of the gas and keeping constant pressure of 101 kPa, methane permeability linearly increases. This change is more significant for T-PNMs compared to D-PNM. However, the change hardly exceeds 20% for 500 • of temperature variation. It should be noted that in the case of subsurface flow modelling in heterogeneous porous material such as the case in tight carbonate reservoirs or gas shales, due to the higher range of environmental pressure, the effect of micro-pores in permeability will become less significant. In contrary due to a higher subsurface temperature, molecules will have more interaction with solid walls and non-viscous flow mechanisms will be still important.
Finally, we have examined the effect of changing the gas molecular mass which is hypothetically equivalent to changing the gas type (Fig. 12c). It has been found that for lighter gases in the standard temperature and pressure condition, T-PNM permeability is around 20% higher than the heavy gases with molecular mass above 100 g/mol. This change is less significant in D-PNMs and its almost less than 10%. Gas permeability gradually approaches to a constant value when molecular mass increases and this means that mean free path of the molecules becomes so small that they less frequently interact with the solid walls of the throats compared to interacting with themselves.

Conclusions
In this study, we introduced the concept of triple pore network models that incorporated a meso-pore network coupled with fractures and micro-porosities. The proposed triple pore network presents a more detailed model of a porous material in which two levels of porosity exist in connection with fracture or channel-like elements. After describing the method for developing these networks, we simulated steady-state gas flow through three porous samples in different range of pressure, temperature and gas molecular masses. The a b c Fig. 12 Sensitivity analysis of gas permeability due to changes in average pressure, temperature and molecular mass, a sensitivity to pressure for methane in 25 • C , b sensitivity to temperature for methane in 101 kPa pressure, c sensitivity to molecular mass in 101 kPa pressure and 25 • C developed method was validated by comparing the gas permeability results obtained in simple geometries and analytical solutions. We showed that using a T-PNM to model a fractured porous media with two scales of porosity made substantive differences compared to a D-PNM. Concisely, the main conclusions and findings of this research can be drawn as follow: • Triple pore network models are introduced and constructed based on the semi-real tomography images in order to simulate gas and liquid permeabilities. • A unified approach is invented to simulate fluid flow within a network with both mesopores and fractures. In this approach we used distance map of the channels with arbitrary cross sections and estimated the LBM permeabilities using an empirical correlation. The benefits of using this approach are that we bypass the lateral segmentation problem in pore network extraction as well as no requirement to locate the fracture elements explicitly. • A hypothetical micro-network structure is presented for taking into the account the hydrodynamical effects of micro-porosities in liquid and gas flow through porous samples with dual scales of porosity. • Effects of micro-porosities on gas permeability are investigated and operational conditions in which ignoring the presence of micro-porosities make a significant error are discussed. • Within the tested range of operational conditions, we have recorded between 10 and 50% relative error in gas permeabilities if micro-porosities are dismissed in the presence of the fractures.

Nomenclature
Here is the list of different variables or parameters used in the main text and in the appendices: Variable Definition Function that gives back the probability of micro-tubes with the radius equal to r m K s Absolute permeability of a solid element with micro-porosity K a Absolute permeability of the network link K n Knudsen permeability K e Total equivalent analytical permeability of the geometry K f The absolute permeability of the fracture F r Gas slippage correction factor Variable Definition K p Overall permeability of the meso-PNM Mean free path of gas

A i
The effective surface area of a network link in the direction of flow L i The length of a network link and its equal to the distance between two connected nodes M Molecular mass of gas Average gas density T Average temperature g Gas viscosity k b Gas Boltzmann constant Additionally, considering the differences between the definition of technical terms in the literature and in order to clarify the vocabulary used/defined in this paper, the following terminology list is provided:

Term Definition
Pore Short way to imply pore body Throat Short way to imply pore throat Micro-pore A hypothetical volume-less joint that connects the micro-throats, unresolved in the tomography images Micro-throat A bundle of micro-tubes that connects micro-pores to other elements in the network Micro-tubes Capillary tubes with sub-resolution radii that form a micro-throat Meso-pore Adequately resolved pore body based on the images from tomography Meso-throat The interface between two meso-pore bodies that is the tightest pathway between them Fracture node A similar element to meso-pores but with an elongated structure Fracture-link The interface between two fracture node that is the tightest pathway between them Throat length The Euclidean distance between the centre of masses of two connecting pore bodies Throat area Cross-sectional area of a throat image which is projected on a 2D surface perpendicular to the throat axis Throat axis The straight line that connects the centre of masses of two adjacent pore bodies Micro-porosity Existence of a micro-pore and micro-throat system in a partially solid element within a tomography image of porous material Solid element A partially solid part of the image that hosts micro-porosities

Appendix 1: Fracture Realization
In each sample we have created 12 fractures which half of them are deviated by 30 • and the rest by 60 • from y axis and the angle is measured in the xy plane. These angles only indicate the general trend of the fractures and fracture local curvature changes due to porous media texture. These fractures are generated via a random walk approach which is inspired by the idea used in Mhiri et al. (2015). In this approach, we consider some arbitrary walkers that can move to their neighbouring cells/pixel in each time step in a random manner. In 2D, each walker, has 8 neighbouring pixel to travel. In order to control the general trend of the pathway, a bias weight is put on the pixels which are aligned to the desired fracture direction. So, the walkers wobble randomly but in a larger picture they are gradually moving towards the biased direction. Then, we run this process for several times and calculate a cost function for each trial. If a walker path cuts hard grains/solid sections of the porous media in which the CT number is higher, the path cost would increase. Finally, in a random search approach, we select the least expensive path to establish a fracture in that area. Figure 13 illustrates the cost minimization process which is used to generate a horizontal fracture in the Estaillades sample. As can be seen, at the initial generations of the random search, fracture path has a high cost due to hitting many high density zones in the image (Fig. 13a) Fig. 13 Cost minimization process used to generate realizations of fractures within porous material, a fracture path at the beginning of optimization with high cost and several solid hits, b fracture path at the middle of the optimization with partially passed solids, c final fracture path that bypasses most of the solid regions, and d best function cost during the optimization process to find the least expensive fracture pathway bypass some expensive zones (Fig. 13b). And finally, after 1000 generations, a plausible trend for an artificial fracture is obtained in which the sample is cut at the points with the lowest cost while maintaining the optimization limits such as the maximum number of steps (Fig. 13c).
It is noteworthy that the cost function is defined as the averaged values of the voxels within the fracture path. Also, the voxel cost of the whole map is normalized to be between 0 and 1 (Fig. 13d). A simple version of this fracture approximation code is published in the public domain. 2

Appendix 2: Lattice-Boltzmann Simulation
We used LBM to numerically calculate absolute permeability of pores and fractures, and compared it with ITPM results. In order to simulate a steady-state, single-phase, and fully developed flow in a channel-like link in the network (pore throat or fracture), we assume it has arbitrary cross section and a uniform shape along its length. A 2D cross section of the fluid conductor is taken 6 times and they are stacked upon each other to create the simulation geometry. Image voxels will be simulation grids so if the conductor image is 50 by 50 voxels, the simulation grid is 50 by 50 by 6 grids. The scheme of LBM is D3Q19 which means there will be 18 possible directions for fluid flow within each block. The code that has been modified and used to calculate LBM permeability is originally published by Haslam et al. (2008). In this simulation, fluid is assumed to be Newtonian with BGK collision model (Flekkoy 1993). Additionally, bounce-back is only in the direction normal to the geometry boundary. The same as the regular formulation of LBM, we have assumed that the particle distribution function f i evolves in the directions of the distribution vectors e i at each time step ( t + Δt ) and location (x) as Haslam et al. (2008): where f i eq is a truncated Maxwell-Boltzmann equilibrium distribution and it can be written as a function of local velocity (v) in 18 directions and is equilibrium state time. Also i is the distribution function index for different neighbouring cells. Simulation starts by imposing a fixed flow velocity at the inlet of geometry (first layer) and periodic boundary condition is considered in the direction of flow. We have defined the convergence criteria when calculated permeability value is stabilized and its relative error compared to the previous time step is less than 10 −6 . More details on LBM permeability calculation is provided in Rabbani and Babaei (2019). Freeman et al. 2011). Absolute permeability of the network links is commonly addressed as the liquid permeability. However, gas permeability is partially different than Newtonian liquids due to the molecular interactions with solid walls (Wu and Pruess 1998). In this paper, we calculate the absolute permeability of gas in different temperature-pressure conditions for the whole pore network with a steady-state approach. Due to the assumed small pressure changes from inlet to the outlet of the PNM, gas properties such as viscosity and density are taken constant. Based on Javadpour et al. (2007), Javadpour (2009) formulation for gas flow in nano to micro size channels, and by assuming three flow mechanisms of Knudsen, slippage and viscous, the gas permeability of a channel can be calculated as Javadpour et al. (2007), Javadpour (2009), Naraghi and Javadpour (2015): where M is molecular mass of the gas, k b is the gas Boltzmann constant, T is average temperature, ̄ is the average density of the gas, g is gas dynamic viscosity, r is tube equivalent radius, F(r) is a correction factor to include slip flow in gas permeability as a function of r, and K a is the absolute or liquid permeability of the channel which is independent to the fluid properties. K n (r) describes the channel permeability increment due to the Knudsen diffusion. This term is a function of gas properties and channel radius (r) and its presence increases the apparent gas permeability of any type of network links including mesothroats, micro-throats and fracture links.
We have previously calculated the absolute permeability of meso-throats and fracture links ( K p ) in Eq. 1, so this value can be replaced in Eq. 6 as K a . We cannot do the same replacement for the gas permeability of micro-throats since they contain a range of microtubes with different sizes and we need to calculate slippage ( F(r m ) ) and Knudsen terms ( K n (r m ) ) for each size of the tubes, and then integrate the calculated permeabilities over different micro-tube radii ( r m ). Thus, the gas permeability of all different links within the T-PNM can be calculated as: where R s is the spatial resolution of the image in μm∕voxel . Additionally, gas slippage term for all types of channels can be calculated as Javadpour et al. (2007), Javadpour (2009), Zhang et al. (2015): where r can be replaced by r m in the case of micro-throats, c is a collision proportionality factor, and is commonly set equal to 1, and ̄ is average mean free path of the gas and can be calculated as Bird (1983): K g = K n (r) + F(r)K p (D)R 2 s meso-throat and fracture-link ∫ f (r m )K n (r m ) + 1 8 f (r m )F(r m ) m r 2 m dr m micro-throat (9) F(r) = 1 + 4cr where P is gas average pressure. Equations 9 and 10 can be used for all types of network links including micro-tubes with the radius of r m and meso-throats or fracture links with the equivalent radius of r. By assuming that gas properties remain constant in a small pressure range of the flow, we can write a steady-state gas flow balance for each pore based on Darcy's law (Whitaker 1986;Huang et al. 2016): where A i is the effective surface area of a network link in the direction of flow, ΔP i is the pressure difference between two sides of the link, L i is the length of the i th link, and n is the total number of network links connected to the specified node and g is gas viscosity obtained by empirical correlation developed by Lee et al. (1966) at the average temperature-pressure condition: where X, Y and Z are correlation constants used to estimate the gas viscosity. By writing Eq. 11 for each of the nodes in the network, a set of equations is given in which the only unknown parameter is the node pressure ( P i ). In order to solve this set of equations, we apply a small pressure difference of 1 Pa between two opposite faces of the PNM as a boundary condition. Other sides of the geometry are assumed as no-flow boundaries. Then, using the biconjugate gradients method (Fletcher 1976) coupled with a lower-upper (LU) factorization (Saad 1994), we solve the system of equations obtained in an iterative manner to find the pore pressure at the centres of each node. Now, using Darcy's law, we find the volumetric flow rate of each network link as well as the overall flow rate. Finally, writing Darcy's law for the whole network, the equivalent gas permeability of the T-PNM is obtained.