## Micro/Nano-pore Network Analysis of Gas Flow in Shale Matrix

#### Pengwei Zhang^{1}, Liming Hu^{1}, Jay N. Meegoda^{2} & Shengyan Gao^{2}

## Abstract

The gas flow in shale matrix is of great research interests for optimized shale gas extraction. The gas flow in the nano-scale pore may fall in flow regimes such as viscous flow, slip flow and Knudsen diffusion. A 3-dimensional nano-scale pore network model was developed to simulate dynamic gas flow, and to describe the transient properties of flow regimes. The proposed pore network model accounts for the various size distributions and low connectivity of shale pores. The pore size, pore throat size and coordination number obey normal distribution, and the average values can be obtained from shale reservoir data. The gas flow regimes were simulated using an extracted pore network backbone. The numerical results show that apparent permeability is strongly dependent on pore pressure in the reservoir and pore throat size, which is overestimated by low-pressure laboratory tests. With the decrease of reservoir pressure, viscous flow is weakening, then slip flow and Knudsen diffusion are gradually becoming dominant flow regimes. The fingering phenomenon can be predicted by micro/nano-pore network for gas flow, which provides an effective way to capture heterogeneity of shale gas reservoir.

## Introduction

According to the U.S. Energy Information Administration (EIA) data, as of 2012 shale gas is more than 40% of U.S. domestic total natural gas production and it is expected to reach 50% in 2039. The Chinese government is also determined to exploit shale gas resources and is planning to reach an annual production rate of 60 to 100 billion cubic meters by 2020^{1}. Hence it is essential to understand shale reservoir properties and gas flow mechanisms in order to achieve above shale gas production goals.

Shale is an ultra-tight rock with relatively low pore connectivity^{2,3}, and extremely low permeability of 10^{−18 } to 10^{−21}m^{2 4,5}. Hydraulic fracturing breaks up shale matrix and connects natural fractures to create flow paths to improve the reservoir connectivity. Thereafter, free gas can quickly release along hydraulic fractures, and gas stored in the matrix slowly transport to connected fractures^{6,7}. Finally, the adsorbed gas in pore surfaces begins to desorb due to reduction in pressure gradient. Therefore, gas storage in the shale matrix and matrix gas flow pattern are important for prediction of gas production.

Pore size of shale matrix ranges from several nanometers^{7-12} to several hundred nanometers, and pore throat size is even much smaller^{9-13}, which is closer to the mean free path of gas molecule. Therefore, gas flow in such nano-scale channels is not governed by viscous flow, where the drag force along rough pore surface and gas molecular collision with pore wall become non-ignorable. Knudsen number (Kn) is defined to separate different dominant flow regimes, and expressed as Kn=λ/r, which is the ratio between gas molecular mean free path (λ) and gas flow characteristic size (r) in porous media^{14-18}. The Knudsen number accounts for the frequency of molecule-molecule and molecule-wall collisions based on gas molecular dynamics theory^{19,20}, which also can be used to separate different flow regimes. For very small values of Knudsen number (Kn<0.01), the characteristic size of gas flow channels is relatively larger than the mean free path of gas molecule, which enables gas flow in continuous state and viscous flow type is available. As the Knudsen number becomes larger (0.01<Kn< 0.1), no-flow assumption near pore surface becomes invalid, and slipping along pore surface occurs^{21,22}. Thereafter, in the transition flow regime (0.1<Kn<10) gas slippage along pore surface and collision with pore wall simultaneously occur^{23}. When Kn>10, gas flow regime becomes free molecular flow, in which gas molecule-molecule interaction turns weak, and molecule-wall interaction is intensified. In real shale reservoir gas flow, flow regime changes with pore gas pressure, and hence a single flow regime cannot accurately predict mass flux and pressure distributions.

Therefore, a theoretical model to capture dynamic gas flow under various flow regimes would be of valuable, and has attracted interests of many researchers. Beskok and Karniadakis (1999) developed a unified gas flow model by defining a rarefaction factor, which explained flow regime changes with variation in gas molecular interaction frequency. Florence et al. (2007) based on large volume of industrial data, defined a pseudo Knudsen number and modified the model proposed by Beskok and Karniadakis (1990), which provided refinements to previous models, especially for ultra-tight porous media24. Additionally, Civan (2010) proposed a simple inverse power-law relationship of rarefaction factor and correlated with experiment data. These unified models were obtained from straight tube and the validity for nano-scale shale matrix is uncertain. The flow regimes in typical shale reservoir are mainly slip and transition flows^{25–28}, hence the gas flow model for shale matrix should combine slip flow and Knudsen diffusion^{23,26,29,30}.

A true molecular simulation using Molecular Dynamics (MD) or Lattice Boltzmann Method (LBM) would be a preferred approach. Chen et al. (2015) simulated the three-dimensional (3D) nanoscale porous structures of shale using a reconstruction method called Markov chain Monte Carlo (MCMC) based on SEM images of shale samples. Then they used LBM to simulate nanoscale shale gas flow. The method proposed by Chen et al. (2015) is quite fundamental but computationally very intensive, hence a fairly accurate but simpler method using typical shale reservoir data is proposed in this research.

The pore-scale simulation is an efficient way to understand micro/nano scale fluid flow^{31-42}. The shale matrix is porous media of low connectivity and it is time consuming to perform permeability tests. Hence pore network models may provide an effective way to obtain basic hydraulic parameters for shale reservoir, to understand dynamic migration of shale gas and to predict reservoir gas production. The pore size distribution in shale matrix has bimodal characteristics^{10,43}. Mehmani et al. (2013) established a multi-scale pore-network with a constant coordination number of 4 by extracting pore-network from a dense random pack of spheres by Delaunay tessellation method, and the results show that gas flow in shale matrix is mainly determined by nano-pore fraction.

In this study, a mathematical micro/nano-pore network model is proposed, where pore sizes, pore throat sizes and coordination numbers satisfy typical shale reservoir data. Isolated pores and clusters also exist in this pore network, while the connected backbone is the research focus of this study. Furthermore, the pore-network model has to satisfy the law of statistics, which is verified in a numerical simulation. Based on the proposed pore network, numerical experiments of apparent permeability tests were conducted for two types of pore throat distributions, i.e., a constant pore throat size and that obeys normal distribution. Then the dynamic gas migration after considering Klinkenberg slip effect and Knudsen diffusion is analyzed. The slip flow and diffusion theories coexist in the shale matrix pore network as the flow regimes range from slip to transition flow^{23,25}. In this work, the two flow regimes for gas mass flux contribution at different stages are considered and dynamic gas flow is analyzed. Finally, the fingering effect in the pore-scale network model is incorporated, hence the scaled up pore network model can effectively capture the heterogeneity of the porous media.

## Results

### Structure of shale matrix pore-network model

Figure 1(a) presents a generated fully-connected micro/nano pore network. In order to match the measured data, the interconnected coordination bonds were eliminated at inlets and outlets to avoid planar flow. A dilution procedure was developed to reflect the low connectivity of shale matrix, and the backbone of the pore-network after the isolated pores and isolated pore clusters were eliminated as shown in Fig. 1(b). The gas flow is along X direction as shown in Fig. 1(c). Pressure difference was applied to the inlet and outlet boundary surface. A constant pressure (1MPa) was maintained at the outlet boundary during shale gas exploitation.

The other boundary surfaces are set no flow conditions. In this research, the pore size distribution can be divided into two parts: large pores and small pore throats.

Both large pores and small pore throats obey normal distribution, where the average pore size is 300 nm and average pore throat size is 6 nm^{13}. The pore diameter ranges from 100 nm to 500 nm, and pore throat size ranges from 1 nm to 10 nm. As reported in many literatures^{44–46}, the average coordination number for sandstones is approximately 4. The coordination number decreases with decreases in porosity, thus average coordination number is assumed as equal to 3 for shale matrix in the proposed pore network model is reasonable, and it also obeys normal distribution. The upper bound of coordination number is 26 as each pore can be connected to 26 adjoining pores for gas transport^{47,48}.

## Theoretical gas flow model of shale matrix

During the dynamic gas flow, the variation in Knudsen number results in the change of flow regimes. Figure 2 illustrates the relationship between the Knudsen number and the pore throat size as well as reservoir gas pressure. The Pore throats were used in this work instead of larger pores for gas flow computation in shale matrix. When the pore throat sizes are used for the flow computation the Knudsen number ranges from 0.01 to 1, and the Knudsen number may larger than 1 when the reservoir pressure decreases. Therefore, the typical gas flow in shale matrix is mainly slip and transition flow. Hence it is necessary to combine the slippage effect and Knudsen diffusion in the theoretical model^{23}.

The single fluid phase flow rate in a pipe can be described by Hagen-Poiseuille equation, if it is applied to gas phase considering gas compressibility, the mass flux equation is written as:

where ρ is the gas density equals pM/RT, p is the reservoir gas pressure in Pa, M is gas molar mass in kg/mol, R is universal gas constant in J/mol/K, T is the thermodynamic temperature in K, μ is the gas dynamic viscosity in Pa s. For slippage effect, the gas slip along pore wall was considered by proposing a correction factor (F)^{22},

where c is a collision proportionality factor, normally given as equal to 1.0, λ is the gas molecular mean free path, which can be obtained by the gas kinetic theory:

Then the modified gas mass flux is

The diffusion part can be described by the Knudsen diffusion model^{23}:

Combining Eq. (4) and (5) to obtain:

where J is the mass flux in kg/m^{2}/s. Eq. (6) was used to calculate the discharge through each coordination bond connecting two adjacent pores. It is in essence similar to Javadpour (2009) flow rate equation, except a minor difference in the slip part as shown in Fig.3.

Figure 3 shows several typical apparent model simulation results. As it indicates, apparent permeability is higher than intrinsic permeability, especially for higher Knudsen numbers. The apparent permeability models of Florence et al. (2007) and Civan (2010) are multi-flow regimes, which are obviously higher than results of Klinkenberg (1941). The model proposed by Brown et al. (1946) is similar to Klinkenberg (1941), while the former considered the wall roughness by bringing in the tangential momentum accommodation coefficient. When the tangential momentum accommodation coefficient (α ) increases the apparent permeability decreases. Specifically, when α approaches 1, the two apparent permeability models are the same. Additionally, when α approaches 0.8 the Brown (1946) model is very close to those of Florence et al. (2007) and Civan (2010) even in transition flow region.

This indicates that Brown (1946) model can account for the gas molecular interactions with pore walls. Therefore, in present mass flux model Klinkenberg (1941) slip flow combined with Knudsen diffusion (Eq. (6)) was applied in the transient flow region, which may accurately account the contribution from the two regimes separately. Javadpour (2009) studied the gas flow regimes in nanoscale shale matrix, and the Brown’s slippage effect and Knudsen diffusion were considered in the apparent permeability model. The apparent permeability model proposed by Chen et al. (2015) was based on the Dusty gas model (DGM), which considered the viscous flow and Knudsen diffusion.

## Variation of Dominant flow regimes

In shale reservoir, free gas release quickly from micro fracture and connected shale matrix after hydraulic fracturing. With the decrease of reservoir gas pressure, he dominant flow regimes for gas exploitation varied at each stage. Based on gas kinetic theory, the variation of macro flow regimes is due to the interaction of gas molecules with themselves and with the wall. At each stage, the flow regimes include viscous flow, slip flow and Knudsen diffusion, and they are mixed with different proportions as shown in Fig. 4.

During the viscous flow range and slip flow range, Knudsen diffusion is almost negligible and the pressure driven flow is the main gas production source. However, the mass flux ratio contributed by typical Hagen-Poiseuille viscous flow model gradually decreases with the increase of Knudsen number, especially in the transition flow range and the viscous pipe flow almost extinct. In contrast, Klinkenberg slippage effect and Knudsen diffusion tend to play a major role in transition flow range. Specifically, total flow in Fig. 4 includes Hagen-Poiseuille flow and Klinkenberg slip flow, and the latter takes up almost 40% of total flux when Knudsen number is larger than 1. The Knudsen diffusion starts to enhance sgnificantly when the Knudsen number is larger than 0.1, and it contributes to more than 60% of total mass flux at the later stage of transition flow range. Therefore, at the later stage of gas exploitation, gas production is mainly due to Knudsen diffusion and slip flow.

## Dynamic gas flow in pore-network model

Gas flow through the pore network is a dynamic process and each pore satisfies conservation of mass equation. Take pore i for example, it connects with several adjacent pores. During each time step, pore throat pressure is equal to the average pressure of connected two pores. From the time step k to k+1, the mass changes in pore i can be balanced by the summation of all the mass flux related to pore i:

For ideal gas, the gas mass store in pore i at k time step is:

Substitute Eq. (8) into Eq. (7) yields:

Eq. (9) shows the iteration process for pressure. At the start, pressure difference only applied to the boundary pores and pores connected to it, subsequently it will dynamically cover all pores in the pore network. During this computation, time step is crucial for convergence. The Courant number which reflects the relation between time step and space step was used to verify the validity of numerical simulation. Courant number should smaller than 1 to ensure the computational stability, thus the range of time step can be determined by Eq. (10).

where the lcb is the coordination bond, which connects, pore i and pore j in the pore network.

The advantage of pore network simulation is that the dynamic gas transfer can be captured. Figure 5(a) presents the variation of gas flow rate of whole pore network. At the time of hydraulic fracture, gas flow rate near the outlet or dowmstream fracture is higher because of the higher initial pressure gradient. With the passage of time after hyraulic fracture, the gas flow rate in the inner layer shows an ascending tendency. After relatively balanced gas flow rate is reached, like the time node of 1 ms, gas rate of all layers of the pore network gradually decrease. In order to explore the relation of gas flow rate changes of three typical pore network layers (downstream, middle and upstream layers) are analyzed in Fig. 5(b).

The pressure gradient for downstream layer is continously decreasing, which leads to steadily decreasing gas flow rate during gas exploration. While the gas flow rate for middle and upstream layers experience an ascending rates initially followed by a declining rates. This observation can be explained using the initial condition of the reservoir, which is constant pressure for all pores, therefore, the initial pressure gradient for inner layer is zero, and it experiences increase initially and then decrease later. Additionally, at the end of gas exploration the larger pressure gradient isoline moves into inner layers, as shown in Fig. 5(b).

## Discussion

### Apparent permeability in shale matrix.

Figure 6 shows the apparent gas permeability of shale matrix using the proposed pore network, where two types of pore throat size distribution at different reservoir pressure are compared. Scenario 1 is for constant pore throat size, and Scenario 2 is for pore throat size obeying normal distribution. The results indicate that apparent permeability is sensitive to reservoir gas pressure and pore throat size. The gas density increases if reservoir gas pressure increases, which leads to decrease in kinematic viscosity (dynamic viscosity divides gas density) and decrease in apparent permeability. As it shows, the apparent permeability of constant pore throat size model is larger than that for the normally distributed one. Although the apparent permeability for Scenario 1 is approximately 2 times larger than Scenario 2, it still indicates that apparent permeability value of shale depends on the proportion small pore throat. Additionally, the results also show that the apparent permeability obtained from laboratory tests under lower pressures is higher than the permeability values of the real reservoir.

### Gas flow rate using different theoretical models

Figure 7 compares the gas flow rates from three theoretical models (Hagen-Poiseuille flow, Klinkenberg slip flow, and flow by present model) using the proposed pore network model, the results show that the gas flow rate of proposed model considering slip flow and Knudsen diffusion is significantly higher than viscous flow model and Klinkenberg slip flow at the early time. However, according to the initial and boundary conditions for the pore network, the stored mass of gas in the shale matrix is assumed and hence the dynamic gas flow during gas extraction can be divided into three stages. At the initial stage, the reservoir gas pressure is relatively high and the pressure gradient decreases rapidly (larger slope) for the proposed model than that for other two models. Thus the gas flux ratio, defines as the flux calculated from the proposed model to that of viscous model, is gradually decreases from a starting value of slightly over 3.2. In the middle stage, with decreases of pore gas pressure in the proposed micro/nano network model, there will be a transition region where the gas flux calculated from the proposed model is less than that from other models (gas flux ratio is less than 1). In the final stage, the slope of gas flux ratio at the transition region is gradually decreasing, and at the last period of gas exploitation the low gas pressure leads to a higher Knudsen number, therefore the slip flow and Knudsen diffusion dominates and hence the gas flow ratio begins to increase (the ratio tends to increase slightly).

### Fingering effect

Heterogeneity is quite common in rock and soil. Macro simulation of heterogeneity in hydrogeology researches often assumes the hydraulic conductivity satisfy certain distribution^{49,50}. Unfortunately, the hydraulic conductivity of in-situ shale reservoir is too difficult to measure due to deep burial conditions and extremely low permeability values. The Micro/Nano pore network analysis of gas migration in shale matrix may provide insights to the study the heterogeneity.

Figure 8 shows a pressure distribution of the same pore network layer where the pore size distribution is according to the real shale matrix pore size range and the gas pressure ratio (gas pressure of pore i divides maximum pore pressure of this layer) is a curved surface demonstrating the fingering effect. It means that gas flow occurs through least resistant paths corresponding to larger pore throats in this pore network model of the shale matrix. In the proposed pore-network model, the micro fractures are not considered and the model size is only several microns. The up-scaled pore network model considering micro fracture should be able to accurately account for the fingering effect of preferential flow in real shale reservoir gas flow. However, it is out of the scope of this work.

The pore-scale simulation in shale matrix is an attractive way to study gas flow and predict long-term gas production. The impact of isolated pores in the pore network is not considered in this paper. Furthermore, with the gas release, the effective stress of shale reservoir increases, and this may lead to compression of reservoir and change of the permeability of shale matrix, which is also not considered in this paper. Hence it is quite promising to further develop the pore network model to perform more realistic simulations.

### Methods

#### Porosity of the pore-network model

The methodology to generate pore network for typical porous media has been clearly presented^{46-48}, however the low connectivity is the key characteristic of shale matrix, which should be adequately considered during the construction of pore structure. As mentioned before, the pore size, porosity and pore connectivity are very small for gas shale matrix. These parameters are essential for micro/nano-pore network analysis.

According to data from Barnett, Marcellus, Haynesville, and Eagle Ford shale, the shale porosity ranges from 1% to 10%^{7,8,10,51}. To generate the regular micro/nano-pore network, the initial porosity was assumed to be 7% for the fully connected network, which is consistent with data for the Marcellus shale7. For a regular pore network, the pore spacing (pore center to center distance) is constant and determined by porosity. First, the number of pores in each direction is assumed in this model (n_{x}=10, n_{y}=15, n_{z}=10), which should satisfy representative property of gas flow^{32,40,48}. Then, the volume of all pores is calculated, and with the porosity the volume of the model V is obtained. Finally, the boundary length in each direction can be obtained using model volume V (l_{x}*l_{y}*l_{z}= V), and pore spacing by dividing length by number of pore numbers in each direction. Although pore center distance is constant for the whole model, the coordination bond length is varied by connecting to different pores where the coordination bond length (lcb) is equal to pore center distance minus the radius of the adjacent two pores.

where l_{ij} is the pore center distance between pores i and j in m, r_{i } and r_{j} are the radius of pores i and j respectively in m. When two adjacent pores are interconnected, which means there is no coordination bond, then the two pores are merged into one larger pore^{47}.

#### Pore connectivity

Narrow pore throat is disordered in shale matrix. Therefore, a multi-directional pore throats may capture this property. Raoof and Hassanizadeh (2009) accomplished this by assigning 26 coordination number for each nodes. The 26 coordination bonds were not suitable for low connective shale and dilution coordination bond is necessary. Gao et al. (2012b, a) proposed rigorous dilution procedure for dynamic two-phase flow in porous media. However, this dilution procedure is not suitable for shale matrix as shale matrix is a relatively low interconnected porous media. The dilution procedure in shale matrix is similar to that in percolation theory where each coordination bond has a probability threshold which determines if the bond is opened or blocked^{48,52}. The probability threshold for bond connect pore i and pore j is as follows:

where N_{ri}, N_{rj} are the random distribution coordination number for pores i and j respectively, N_{ai}, N_{aj} are the assigned coordination number for pores i and j in the fully connected pore network. Similar to the term “elimination number” defined by Raoof and Hassanizadeh (2009), if the elimination number for a bond is larger than the corresponding probability threshold, then that bond will be eliminated. As average coordination number for shale matrix is relatively small, the probability threshold for each bond will be much smaller. Therefore, elimination number obeying random distribution between 0 and η, where η is the reduction factor ranging from 0 to 1 was used. The reduction factor depends on the connectivity of pore network and it is a function of average coordination number and porosity. If the connectivity is high, η can be set smaller, otherwise, there will be more isolated pores. Besides, η can be different in each direction, which makes our model quite flexible to simulate material anisotropy. The pore-network model in this work can also be simplified to a constant coordination number model^{47}. For shale matrix pore network model, η value has to satisfy statistics law, which is verified in a numerical simulation. Figure 4(b) shows the interconnected pore network after dilution. For this average coordination number and porosity in shale matrix pore network, 0.45 is an upper bound value for reduction factor. The statistic value from 25 groups tests were conducted, and the coefficients of variation (total pores and pore coordination numbers) are almost within 5%, which shows a good stability.

#### Permeability tests

Steady state flow is usually assumed in permeability tests of pore network models^{3,47}, as it can establish mass balance equation

for each pore. Unfortunately, for shale matrix pore network in this work, with low connectivity directly solving the equation set will lead to a large sparse matrix which is difficult to invert. Furthermore, as it is nonlinear equation and the convergence for iteration method like Newton-Raphson is time consuming. Therefore, a transient permeability test method was applied, which set constant pressure difference between inlet and outlet, the gas flow will reach steady state and the permeability can be obtained at that state. The convergence criteria is defined as that searching for the maximum pressure difference during time step k and k+1 for all pores, and let it be less than the allowable error: Max [p_{ki}-p_{kj}]≤ err. The allowable error to reach steady state was set as: err= poulet/104, it is 0.1 kPa and much less than the MPa range in reservoir gas pressure.

After the steady state was reached, the permeability for the whole shale matrix pore network can be calculated as follows:

## References

1. Shi, M. et al. Bromide: A pressing issue to address in China’s shale gas extraction. Environ Sci Technol 48, 9971–9972, doi: 10.1021/es502848p10.5942/jawwa.2013.105.0093 (2014).

2. Hu, Q., Ewing, R. P. & Dultz, S. Low pore connectivity in natural rock. J Contam Hydrol 133, 76–83, doi: 10.1016/j. jconhyd.2012.03.006 (2012).

3. Mehmani, A., Prodanović, M. & Javadpour, F. Multiscale, Multiphysics network modeling of shale matrix gas flows. Transport Porous Med 99, 377–390, doi: 10.1007/s11242-013-0191-5 (2013).

4. Amann-Hildenbrand, A., Ghanizadeh, A. & Krooss, B. M. Transport properties of unconventional gas systems. Mar Petrol Geol31, 90–99, doi: 10.1016/j.marpetgeo.2011.11.009 (2012).

5. Sakhaee-Pour, A. & Bryant, S. L. Gas permeability of shale. SPE Reserv Eval Eng 15, 401–409, doi: http://dx.doi.org/10.2118/146944PA (2012).

6. Lee, D. S. et al. A critical evaluation of unconventional gas recovery from the marcellus shale, northeastern United States. KSCE J Civil Eng 15, 679–687, doi: 10.1007/s12205-011-0008-4 (2011).

7. Soeder, D. J. Porosity and permeability of Eastern Devonian gas shale. Transport Porous Med 3, 116–124, doi: http://dx.doi. org/10.2118/15213-PA (1988).

8. Curtis, M. E., Sondergeld, C. H., Ambrose, R. J. & Rai, C. S. Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging. AAPG Bull 96, 665–677, doi: 10.1306/08151110188 (2012).

9. Loucks, R. G., Reed, R. M. & C.Ruppel, S. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett shale. J Sediment Res 79, 848–861, doi: 10.2110/jsr.2009.092 (2009).

10. Chalmers, G. R., Bustin, R. M. & Power, I. M. Characterization of gas shale pore systems by porosimetry, pycnometry, surface area, and field emission scanning electron microscopy/transmission electron microscopy image analyses: Examples from the Barnett, Woodford, Haynesville, Marcellus, and Doig units. AAPG Bull 96, 1099–1119, doi: 10.1306/10171111052 (2012).

11. Firouzi, M., Rupp, E. C., Liu, C. W. & Wilcox, J. Molecular simulation and experimental characterization of the nanoporous structures of coal and gas shale. Int J Coal Geol 121, 123–128, doi: 10.1016/j.coal.2013.11.003 (2014).

12. Zhang, X., Xiao, L., Shan, X. & Guo, L. Lattice Boltzmann simulation of shale gas transport in organic nano-pores. Sci Rep 4, 4843, doi: 10.1038/srep04843 (2014).

13. Gao, Z. & Hu, Q. Estimating permeability using median pore-throat radius obtained from mercury intrusion porosimetry. J Geophys Eng 10, 025014, doi: 10.1088/1742-2132/10/2/025014 (2013).

14. Freeman, C. M., Moridis, G. J. & Blasingame, T. A. A Numerical Study of Microscale Flow Behavior in Tight Gas and Shale Gas Reservoir Systems. Transport Porous Med 90, 253–268, doi: 10.1007/s11242-011-9761-6 (2011).

15. Civan, F. Effective correlation of apparent gas permeability in tight porous media. Transport Porous Med 82, 375–384, doi: 10.1007/s11242-009-9432-z (2009).

16. Ziarani, A. S. & Aguilera, R. Knudsen’s permeability correction for tight porous media. Transport Porous Med 91, 239–260, doi:

10.1007/s11242-011-9842-6 (2011).

17. Roy, S. et al. Modeling gas flow through microchannels and nanopores. J Appl Phys 93, 4870, doi: 10.1063/1.1559936 (2003).

18. Gilron, J. & Soffer, A. Knudsen diffusion in microporous carbon membranes with molecular sieving character. J Membrane Sci 209, 339–352, doi: 10.1016/S0376-7388(02)00074-1 (2002).

19.Thomson, S. L. & Owens, W. R. A survey of flow at low pressures. Vac 25, 151–156, doi: 10.1016/0042-207X(75)91404-9 (1975).

20. Beskok, A. & Karniadakis, G. E. Report: A Model for flows in channels, pipes, and ducts at micro and nano scales. Nanosc Microsc Therm 3, 43–77, doi: 10.1080/108939599199864 (1999).

21. Brown, G. P., DiNardo, A., Cheng, G. K. & Sherwood, T. K. The flow of gases in pipes at low pressures. J Appl Phys 17, 802, doi: 10.1063/1.1707647 (1946).

22. Klinkenberg, L. J. The permeability of porous media to liquids and gases. Paper presented at Drilling and Production Practice New York, New York, USA (1941).

23. Javadpour, F. Nanopores and apparent permeability of gas flow in Mudrocks (Shales and Siltstone). J Can petrol Technol 46, 55–61, doi: http://dx.doi.org/10.2118/09-08-16-DA (2009).

24. Florence, F. A., Rushing, J. A., Newsham, K. E. & Blasingame, T. A. Improved permeability prediction relations for low permeability sands. SPE 107954, 1–18, doi: http://dx.doi.org/10.2118/107954-MS (2007).

25. Shi, J., Zhang, L., Li, Y. & Yu, W. Diffusion and flow mechanisms of shale gas through matrix pores and gas production forecasting.SPE Unconventional Resources Conference Canada, Calgary, Alberta, Canada, doi: http://dx.doi.org/10.2118/167226-MS (2013).

26. Song, H. Q. et al. Dynamic characteristics of gas transport in nanoporous media. Chinese Phys Lett 30, 014701, doi: 10.1088/0256307x/30/1/014701 (2013).

27. Javadpour, F. Nanoscale gas flow in shale gas sediments. J Can petrol Technol 48, 16–21, doi: http://dx.doi.org/10.2118/07-10-06 (2007).

28. Chen, L. et al. Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity. Sci Rep 5, 8089, doi: 10.1038/srep08089 (2015).

29. Shabro, V., Torres, V. & Javadpour, F. Numerical simulation of Shale gas production: from pore scale modeling of slip flow, Knudsen diffusion, and Langmuir desorption to reservoir modeling of compressible fluid. SPE 144355, 1–11, doi: http://dx.doi. org/10.2118/144355-MS (2011).

30. Darabi, H., Ettehad, A., Javadpour, F. & Sepehrnoori, K. Gas flow in ultra-tight shale strata. J Fluid Mech 710, 641–658, doi:10.1017/jfm.2012.424 (2012).

31. Gao, S., Meegoda, J. N. & Hu, L. Simulation of dynamic two-phase flow during multistep air sparging. Transport Porous Med 96, 173–192, doi: 10.1007/s11242-012-0081-2 (2012).

32. Gao, S., Meegoda, J. N. & Hu, L. A dynamic two-phase flow model for air sparging. Int J Numer Anal Met 37, 1801–1821, doi:10.1002/nag.2109 (2013).

33. Blunt, M. J., Jackson, M. D., Piri, M. & Valvatne, P. H. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv Water Resour 25, 1069–1089, doi: 10.1016/S0309-1708(02)00049-0 (2002).

34. Ryazanov, A. V., Dijke, M. I. J. & Sorbie, K. S. Two-phase pore-network modelling: existence of oil layers during water invasion. Transport Porous Med 80, 79–99, doi: 10.1007/s11242-009-9345-x (2009).

35. Acharya, R. C., Sjoerd, E. A. T. M. & Leijnse, A. Porosity–permeability properties generated with a new 2-parameter 3D hydraulic pore-network model for consolidated and unconsolidated porous media. Adv Water Resour 27, 707–723, doi: 10.1016/j. advwatres.2004.05.002 (2004).

36. Valvatne, P. H., Piri, M., Lopez, X. & Blunt, M. J. Predictive pore-scale modeling of single and multiphase flow. Transport Porous Med 58, 23–41, doi: 10.1007/s11242-004-5468-2 (2005).

37. Sheng, Q. & Thompson, K. Dynamic coupling of pore-scale and reservoir-scale models for multiphase flow. Water Resour Res 49, 5973–5988, doi: 10.1002/wrcr.20430 (2013).

38. Gao, S., Meegoda, J. N. & Hu, L. Microscopic modeling of air migration during air sparging. J Hazard Toxic Radioact 15, doi: 10.1061/(asce) (2011).

39. Reeves, P. C. & Celia, M. A. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour Res 32, 2345–2358, doi: 10.1029/96wr01105 (1996).

40. Nogues, J. P., Fitts, J. P., Celia, M. A. & Peters, C. A. Permeability evolution due to dissolution and precipitation of carbonates using reactive transport modeling in pore networks. Water Resour Res 49, 6006–6021, doi: 10.1002/wrcr.20486 (2013).

41. Celia, M. A., Reeves, P. C. & Ferrand, L. A. Recent advances in pore scale models for multiphase flow in porous media. Rev Geophys 33, 1049–1057, doi: 10.1029/95RG00248 (1995).

42. Nordhaug, H. F., Celia, M. & Dahle, H. K. A pore network model for calculation of interfacial velocities. Adv Water Resour 26,1061–1074, doi: http://dx.doi.org/10.1016/S0309-1708(03)00100-3 (2003).

43. Ross, D. J. K. & Marc Bustin, R. The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs. Mar Petrol Geol 26, 916–927, doi: 10.1016/j.marpetgeo.2008.06.004 (2009).

44. Øren, P. E. & Bakke, S. Process based reconstruction of sandstones and prediction of transport properties. Transport Porous Med 46, 311–343, (2002).

45. Lindquist, W. B., Venkatarangan, A., Dunsmuir, J. & Wong, T. F. Pore and throat size distributions measured from synchrotron X-ray tomographic images of Fontainebleau sandstones. J Geophys Res 105, 21509, doi: 10.1029/2000jb900208 (2000).

46. Arns, J. Y. et al. Effect of network topology on relative permeability. Transport Porous Med 55, 21–46, (2004).

47. Gao, S., Meegoda, J. N. & Hu, L. Two methods for pore network of porous media. Int J Numer Anal Met 36, 1954–1970, doi: 10.1002/nag.1134 (2012).

48. Raoof, A. & Hassanizadeh, S. M. A new method for generating pore-network models of porous media. Transport Porous Med 81, 391–407, doi: 10.1007/s11242-009-9412-3 (2009).

49. Feehley, C. E., Zheng, C. & Molz, F. J. A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: application to the macrodispersion experiment (MADE) site. Water Resour Res 36, 2501, doi: 10.1029/2000wr900148 (2000).

50. Bianchi, M., Zheng, C., Tick, G. R. & Gorelick, S. M. Investigation of small-scale preferential flow with a forced-gradient tracer

test. Ground water 49, 503–514, doi: 10.1111/j.1745-6584.2010.00746.x (2011).

51. Vega, B., Dutta, A. & Kovscek, A. R. CT imaging of low-permeability, dual-porosity systems using high X-ray contrast gas.

Transport Porous Med 101, 81–97, doi: 10.1007/s11242-013-0232-0 (2014).

52. Loucks, R. G., Reed, R. M., Ruppel, S. C. & Hammes, U. Spectrum of pore types and networks in mudrocks and a descriptive

classification for matrix-related mudrock pores. AAPG Bull 96, 1071–1098, doi: 10.1306/08171111061 (2012).

#### Acknowledgements

The authors would like to express their sincere gratitude to the financial supports from National Basic Research Program of China (Grant No. 2012CB719804), and National Natural Science Foundation of China (Project No. 51323014, 41372352).

#### Author Contributions

P.W.Z. wrote the main manuscript and prepared the figures. L.M.H. and J.N.M. supervised the theoretical analysis and writing. P.W.Z. and S.Y.G. developed the program of pore-network model. All authors reviewed the manuscript.

#### Additional Information

Competing financial interests: The authors declare no competing financial interests. How to cite this article: Zhang, P. et al. Micro/Nano-pore Network Analysis of Gas Flow in Shale Matrix. Sci. Rep. 5, 13501; doi: 10.1038/srep13501 (2015).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

1 State Key Laboratory of Hydro-Science and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, P. R. China. 2 Department of Civil and Environmental Engineering, New Jersey Institute of Technology, Newark, NJ 07102, USA. Correspondence and requests for materials should be addressed to L.M.H. (email: gehu@tsinghua.edu.cn)