Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity.

Porous structures of shales are reconstructed using the markov chain monte carlo (MCMC) method based on scanning electron microscopy (SEM) images of shale samples from Sichuan Basin, China. Characterization analysis of the reconstructed shales is performed, including porosity, pore size distribution, specific surface area and pore connectivity. The lattice Boltzmann method (LBM) is adopted to simulate fluid flow and Knudsen diffusion within the reconstructed shales. Simulation results reveal that the tortuosity of the shales is much higher than that commonly employed in the Bruggeman equation, and such high tortuosity leads to extremely low intrinsic permeability. Correction of the intrinsic permeability is performed based on the dusty gas model (DGM) by considering the contribution of Knudsen diffusion to the total flow flux, resulting in apparent permeability. The correction factor over a range of Knudsen number and pressure is estimated and compared with empirical correlations in the literature. For the wide pressure range investigated, the correction factor is always greater than 1, indicating Knudsen diffusion always plays a role on shale gas transport mechanisms in the reconstructed shales. Specifically, we found that most of the values of correction factor fall in the slip and transition regime, with no Darcy flow regime observed.

Gas-bearing mudstone (“shale”) reservoirs have become important sources of natural gas production in North America, and are expected to play increasingly important roles in Europe and Asia in the near future. Mudstone (“shale”) reservoirs refer to fine-grained sedimentary rocks that contain self-sourced hydrocarbons1. Experimental observations using advanced techniques such as focused ion beam-scanning electron microscopy (FIB-SEM) find that mudstones include pores, organic matter (kerogen) and nonorganic minerals234. Natural gas trapped in mudstone reservoirs is called shale gas (mainly methane), including free gas in void space as well as adsorbed gas on the solid surface. Shale gas is classified as an unconventional gas because nanosize pores are dominant in mudstones and the permeability of mudstones is extremely low5, leading to quite different transport mechanisms compared to that in conventional reservoirs6.

To extract shale gas from low permeability mudstone reservoirs economically, conductive pathways in the shale matrix need to be artificially generated to transport shale gas to the wellbore. Currently, hydraulic fracturing is a widely used technique to increase the permeability of a shale formation by extending and/or widening existing fractures and creating new ones through the injection of a pressurized fluid into shale reservoirs7. Due to hydraulic fracturing, both shale matrix and fracture system exist inside the shale reservoirs. Therefore, four levels of void space can be identified in shale reservoirs including mesoscale hydraulic fractures, meso/microscale natural fractures, micro/nanoscale inter and/or intra particle pores and nanoscale pores in organic matter5. The majority of pores in organic matter range from 2 nm to 50 nm24.

A fundamental understanding of the multiscale multiple physiochemical processes during shale extraction is crucial for improving the gas production and lowering production costs. While the fracture network greatly determines the productivity of shale reservoirs, the transport of shale gas within the matrix also plays an important role. Because of the nanoscale characteristics of the shale matrix, Darcy's law, which is widely adopted in the conventional formations, cannot realistically describe the variety of the relevant flow regimes other than the viscous flow regime6. Apparent permeability, which takes into account the effects of Knudsen diffusion and viscous-flow, is suggested to use in the Darcy equation68. Since the first study of gas slippage in porous media conducted by Klinkenberg9, several correlations between apparent permeability and intrinsic permeability have been proposed101112131415, and some are capable of describing the four flow regimes, namely, continuum flow regime, slip flow regime, transition regime and free molecular flow regime1415, classified by Knudsen number Kn14. Kn is defined as the ratio between gas mean free path to the pore size. When Kn < 0.001, flow is governed by continuum fluid equations, i.e., the Navier-Stokes equation, with the traditional no-slip velocity conditions on the walls. For larger Kn (0.001 < Kn < 0.1), the flow is in the slip regime where the non-continuum effect can be approximated by an empirical velocity slip near the walls with the bulk flow still governed by the Navier-Stokes equation. Further increase in Kn leads to the transition (0.1 < Kn < 10) and free molecular (Kn > 10) flow regimes. The flows in these two regimes can no longer be described by the Navier-Stokes equation, even with a modified boundary condition, and one has to rely on the Boltzmann equation of gas kinetic theory14.

To date, shale reservoirs and their transport properties are still poorly understood. There have been some studies experimentally observing the nanoscale structures of mudstone234. Due to the ultra-low permeability of mudstones, it is difficult to experimentally investigate the transport processes inside shales as well as to accurately measure the transport properties. Alternatively, advanced pore-scale numerical methods, such as the lattice Boltzmann method (LBM)1617, have been shown to be an efficient alternative181920212223. In the present study, the LBM is employed to simulate fluid flow and Knudsen diffusion processes in shales as well as to predict the transport properties. Three-dimensional (3D) nanoscale porous structures of shales are reconstructed using a reconstruction method called Markov chain Monte Carlo (MCMC) based on SEM images of shale samples. Structural characteristics of the reconstructed shales are analyzed in detail. Tortuosity, effective Knudsen diffusivity and intrinsic permeability of the reconstructed shales are numerically predicted. By considering the total flow flux as a combined result of viscous flow and Knudsen diffusion based on the dusty gas model (DGM)2425, the apparent permeability is determined. The correction factor between apparent and intrinsic permeability is calculated and compared with existing correlations in the literature. To the best of our knowledge, the present work is the first study of numerically investigating the effective Knudsen diffusivity and apparent permeability based on real porous structures of shales.

Results

Structural characterization of reconstructed shales

The first step towards performing the LB simulations is to obtain digitalized shale structures in which different components are labeled with different values. Four 3D nanoscale structures of shales are reconstructed using MCMC based on SEM images of shale samples from Sichuan Basin, China (see Methods). Fig. 1(a) shows the SEM image of Sample 1, and Fig. 1(b) is the corresponding reconstructed structure with size L × L × L of 1250 nm × 500 nm × 500 nm.

Figure 1

Reconstruction of the shale matrix.

(a) A SEM image of a shale rock; (b)3D reconstructed porous structures of shales; (c) pore size distribution of the reconstructed shale; (d) the “transport” (blue) and “dead” (read) portions of void space in the reconstructed shale.

Quantitative characterization of the shale structure provides the necessary information for the prediction of transport properties. Typically, the void space within a porous medium is an interconnected array of pores that can be characterized by geometrical quantities (porosity, pore size, specific surface area, etc) and topological descriptors (pore connectivity).

Porosities of the four reconstructed shales are 0.191, 0.226, 0.268 and 0.176 (dimensionless), respectively. Pore size distribution (PSD) is an important structure characteristic of porous media. The pore size of each local pore cell is also required for the calculation of the local Knudsen diffusivity (See Eq. (1)), which is determined by the 13 direction averaging method26. Fig. 1(c) shows the PSD of the reconstructed shale shown in Fig. 1(b), which presents the unimodal distributions. For all the reconstructed structures, the pore size is in the range 5~80 nm, indicating the nanoscale characteristics of void space in shales.

Specific surface area is defined as the ratio between the total surface area and the total volume of the solid phase. Although gas desorption is a relatively slow process, it can account up to 50% of the total gas production27. Desorption process is significantly affected by the specific surface area of shales. Total surface area is determined by counting the number of the interface cells between void and solid space. For the structures shown in Fig. 1(b), the total surface area and total volume of solid phase are 721438 (δx)2 and 2043351 (δx)3, respectively, leading to specific surface area as 0.353066(δx)−1. With resolution δx = 5 nm, the specific surface area is 7.06 × 107 m−1 (or 47.07 m2g−1 with kerogen density as 1.5 gcm−328). The specific surface area in the present study is higher than 14.0 m2 g−1 estimated by Chen et al.19, due to high resolution in the present study compared to 12 nm in their study19. The specific surface area for other three reconstructed structures is 46.72, 85.71 and 90.08 m2 g−1, respectively.

Connectivity of the void space is very important, as shale gas transports through the connected pores in the shale matrix. For the reconstructed shales, “transport” and “dead” pore cells are distinguished. A transport pore cell is a cell that belongs to a continuous percolation path throughout the entire domain. If a portion of void phase does not penetrate the entire domain, shale gas entering the shale from one end cannot arrive at the other end through this portion, and this portion is called “dead”. For identifying the “transport” and “dead” portions in the reconstructed shales, a connected phase labeling algorithm is adopted29. This algorithm scans the entire domain, checks the connectivity of a cell with its neighboring 18 points, labels the cell depending on the local connectivity, and finally assigns any portion with a distinct label if this portion is disconnected from other portions in the domain. Here, the neighboring 18 points are used, in coincidence with the D3Q19 lattice model in the LBM framework adopted in the present study. The connectivity of void space is defined as the ratio of the number of “transport” cells to the total numbers of pore cells. Fig. 1(d) shows the void phase cells after labeled. The blue portion represents the “transport” one while the red portions are “dead” which present as small and discrete blobs. For Sample 1, it is found that there are 358 portions of void space with one as “transport” and the remaining as “dead”. The “transport” portion occupies most of the void phase, with a connectivity of 98.0% (447460/456649, 456649 is the total number of pore cells). The connectivity for other reconstructed structures is 99.1%, 99.7% and 99.8%, respectively, indicating good connection of void space in shales of Sichuan Basin, China.

Effective Knudsen diffusivity

When the character length of a system is comparable to or smaller than the mean free path of the gas molecules, collisions between molecules and the solid wall are more frequent than that between gas molecules, leading to the Knudsen diffusion process. The Knudsen diffusivity in a local pore is30 where M (kg mol−1) is molar mass, R the universal gas constant and T (K) temperature. dp (m) is the effective pore diameter and is determined by the 13 direction averaging method26. The effective Knudsen diffusivity Dk,eff, taking into account the pore size distributions as well as structural characteristics of a porous medium, is predicted by the LB mass transport model (See Methods). Fig. 2 shows the simulated 3D distribution of methane concentration in the reconstructed structure Sample 1. It can be seen that the diffusion process is greatly affected by the local void space, showing very complicated concentration distributions. The effective Knudsen diffusivity can be determined based on the concentration fields predicted (See Methods). The ratio between Dk,eff and D0 is 0.0212, 0.0197, 0.0380 and 0.0130 for Samples 1–4, respectively. With D0 as 5.44 × 10−6 m2s−1 (the Knudsen diffusivity with dp = 25 nm and T = 383 K), Dk,eff for Samples 1–4 is 1.25 × 10−7 m2s−1, 1.17 × 10−7 m2s−1, 2.25 × 10−7 m2s−1 and 7.71 × 10−8 m2s−1, respectively.

Figure 2

Methane concentration distributions in the reconstructed shale.

Intrinsic permeability

Fluid flow inside the reconstructed shales is simulated by the LB fluid flow model (See Methods). Fig. 3 displays the 3D distribution of pressure and streamlines of velocity in Sample 1. The pressure difference between the inlet and outlet, which is expressed by the density difference in the simulation (in the LBM, pressure is related to density by p = ρcs2 where cs is the speed of sound), is 0.0002 in lattice units (in LB simulations, all the variables are in lattice units, and conversion between physical units and lattice units can be achieved by matching the dimensionless number, such as Reynolds number3132), leading to the maximum magnitude of velocity about 4 × 10−6 in lattice units, meeting the low Mach number (Ma, defined as the ratio between local velocity to the speed of sound) limit of LB fluid flow model33. As shown in Fig. 3, the streamlines are quite tortuous, and continuous pathways from the inlet to the outlet are quite few, resulting in quite low permeability of shales. The predicted permeability for Samples 1~4 is 240.5nD, 172.75nD, 367.5nD, 135.0nD, respectively, with 1nD ≈ 10−21 m2. According to KC equation34, for simple sphere-packed porous media, the higher the porosity is, the larger the permeability is; however, reverse trend of Dk,eff/D0 and kd for Sample 1 and 2 is obtained, indicating the complexity of the porous structures of shales. Javadpour experimentally measured 152 samples from nine reservoirs in United States and found that 90% of the measured permeability is less than 150 nD5. The predicted permeability of shales from Sichuan Basin of China is quite close to that from reservoirs in United States. Recently, Chen et al. adopted LBM to simulate fluid flow in a reconstructed shale, and the permeability predicted were about 1500~3920nD with the relatively high porosity (0.299)19.

Figure 3

Distribution of pressure and streamline in the reconstructed shale.

Discussion

Tortuosity

Continuum-scale models of transport processes in shale matrix highly depend on empirical relationships between macroscopic transport properties and statistical structural information of shales (permeability vs. porosity, effective diffusivity vs. porosity, etc). Bruggeman equation has been widely used in such models to calculate the effective diffusivity where τ is tortuosity of a porous medium, defined as the ratio between the actual flow path length and the thickness of a porous medium along the flow direction. Higher tortuosity indicates longer, more tortuous and more sinuous paths, thus resulting in greater transport resistance. In Bruggeman equation, tortuosity is commonly defined as an exponential function of porosity with α usually chosen as 0.5. However, Bruggeman equation was proposed empirically from simple sphere packed porous media, which therefore cannot reflect the complex and heterogeneous structures of shales. Based on the effective diffusivity predicted by LB simulations, Eq. (2) is adopted to calculate τ, and then Eq. (3) is then used to determine α. α determined based on the LB simulation results ranges 1.33 to 1.65 for the four samples, much higher than 0.5 in the Bruggeman equation. Pore-scale modeling, which is based on the realistic porous structures of shales, fundamentally reflects the actual transport processes. Therefore, the simulation results indicate that the pathway inside the shales is very tortuous. The Bruggeman equation overestimates the effective diffusivity, and a higher α needs to be adopted in the Bruggeman equation to more accurately calculate effective diffusivity of shales. In fact, for complex nanoscale porous medium, such as catalyst layer of proton exchange membrane fuel cells, the value of α is often found to be much higher than 0.5, where α in the range of 1.0~2.0 is recommended26.

Apparent permeability

Gas slippage in a porous medium leads to higher measured gas permeability (apparent permeability) compared to the liquid permeability (intrinsic permeability)9. Klinkenberg first conducted gas slippage in porous media and proposed a first-order correlation for apparent permeability ka9 where fc is the correction factor and kd is the intrinsic permeability. The parameter bk is the Klinkenberg's slippage factor, which depends on the molecular mean free path λ, character pore size of the porous media r and pressure p9. Based on Eq. (4), various expressions of bk have been proposed in the literature. Heid et al.10 and Jones and Owens11 proposed similar expressions relating bk to kd. Sampath and Keighin12 and Florence et al.13 developed a different form by relating bk to kd/ε. Klinkenberg's correlation is a first-order correlation. Beskok and Karniadakis14 developed a second-order correlation where b is the slip coefficient and equal -1 for slip flow. The rarefaction coefficient α(Kn) is given by the following expression proposed by Civan15 Eq. (5) combined with Eq. (6) is called Beskok and Karniadakis-Civan's correlation in the present study. This correlation has been shown to be capable of describing all four flow regimes including viscous flow, slip flow, transition flow and free molecular flow.

The apparent permeability can also be determined based on the Dusty gas model (DGM)2425. According to this model, the total flux of gas in a pore can be considered as a combined result of viscous flow and Knudsen diffusion2425 where J (kg m−2 s−1) is the mass flux per unit area, µ (Pa.s) is viscosity, ρ (kg m−3) is density and (Pa m−1) is pressure gradient. J is the viscous flow flux term and Jk is the Knudsen diffusion flux term where z is the gas compressibility factor, accounting for the effect of non-ideal gas. Therefore, the total flow flux can be expressed as Shale gas is a mixture of several components, but the methane dominates with mole fraction of about 90%5. Therefore, in this study shale gas is considered to be composed of pure methane, and binary gas diffusion is neglected. Based on Eq. (10), the apparent permeability can be determined Fig. 4 shows the correction factor in a cylinder predicted by Klinkenberg's correlation (Eq. (4)), Beskok and Karniadakis – Civan's correlation (Eqs. (5–6)) as well as DGM (Eq. (11)). For gas transport in a cylinder, the intrinsic permeability is With the mean free path λ given by15 the correction factor based on the DGM is as follows As shown in Fig. 4, when Kn is quite low (Kn<0.01), i.e., in the range of Darcy flow regime, the above equations predict quite similar values as the term with Kn can be neglected. However, as Kn increases, the discrepancy between Klinkenberg's correlation and other correlations becomes larger. The Beskok and Karniadakis-Civan's correlation is more accurate than Klinkenberg's correlation, especially at high Kn regime35. Values predicted by Eq. (11) agree well with Beskok and Karniadakis-Civan's correlation.

Figure 4

Correction factor between apparent permeability and intrinsic permeability predicted by LB simulations and empirical correlations for transport in a cylinder.

In the literature, for determining the apparent permeability of microchannels or micro porous media, in which the fluid flow falls in the slip or transition flow regime, Navier-stokes equations are solved with modified local viscosity and slip flow boundary conditions36, where delicate numerical schemes are very much required. Through such simulations, detailed distributions of velocity and pressure as well as the apparent permeability can be obtained. However, for reservoir simulations, the empirical relationships (permeability vs. porosity, effective diffusivity vs. porosity, etc.) rather than the detailed distributions are more required. The analysis related to Fig. 4 indicates that apparent permeability can be determined based on Eq. (11) by predicting the intrinsic permeability and Knudsen diffusivity of a porous medium. Standard governing equations are solved, and tedious numerical schemes in these models with modified viscosity and slip flow boundary conditions36 are not required.

Using the intrinsic permeability and effective Knudsen diffusivity predicted by the LB simulations, correction factor is calculated based on Eq. (11) and is compared with that predicted by Klinkenberg's and Beskok and Karniadakis-Civan's correlations. Pressure is chosen in the range of 10-4000 psi (1psi = 6894.75729 Pa) and temperature is fixed at 383K. Viscosity of methane is calculated using a online software called Peace Software. In Eqs. (4–6), the following equation is used to calculate Kn15 based on the intrinsic permeability, porosity and tortuosity of the shales predicted by LB simulations.

Fig. 5 shows the correction factor predicted by the LB simulations and empirical correlations under different pressures, and Fig. 6 displays correction factor under different Kn. For the wide pressure range investigated, fc is always greater than 1 as shown in Fig. 5, indicating Knudsen diffusion always plays a role on the shale gas transport mechanisms in the reconstructed shales. This is also confirmed in Fig. 6, where all data fall within the slip flow regime, transition flow regime and Knudsen diffusion regime, with no Darcy flow regime observed as Kn is always greater than 0.03. As the pressure decreases or the Kn increases, the Knudsen diffusion becomes increasingly dominant. For pressure lower than 100 psi, the correction factor can be as high as 10, implying that Knudsen diffusion dominates and the viscous flow can be ignored. Most of the values of correction factor are located in the transition and slip flow regimes. The simulation results are in better agreement with Beskok and Karniadakis-Civan's correlation, in consistence with that in Fig. 4. The values from our simulations are a little higher than that of Beskok and Karniadakis-Civan's correlation. By investigating the experimental data of tight gas sandstones from seven major tight gas basins in the United States, Ziarani and Aguilera35 also found Beskok and Karniadakis-Civan's correlation performed the best, especially for extremely tight porous media, and Klinkenberg's correlation underestimates the correction factor. Note that the Beskok and Karniadakis-Civan's correlation as well as Eq. (1) are derived based on the diffusion process in a cylinder1430. However, pores in the shales are complex and can hardly be described by a collection of cylinders. Using Eq. (1) will overestimate the Knudsen diffusivity in porous media. Therefore, in future simulations, complex structures of the pores as well as the nature of the redirecting collision between walls and gas molecules should be considered to improve Eq. (1)3738.

Figure 5

Correction factor predicted by numerical simulations and empirical correlations under different pressure.

Figure 6

Correction factor predicted by numerical simulations and empirical correlations under different Knudsen number.

Methods

Structure reconstruction

The shale samples came from a well located in the Sichuan Basin, China39. Hitachi S-4800 Scanning Electron Microscope (SEM) in China University of Petroleum was adopted to scan the shale samples, which can achieve a resolution as high as 1 nm. In the present study, the resolution of the SEM images is 5 nm. Totally 13 SEM images were obtained by scanning different samples. Using spatial information derived from these SEM images, a reconstruction method called Markov Chain Monte Carlo (MCMC) method, which can generate 3D pore structures based on 2D thin section images such as SEM images, was adopted to reconstruct the 3D nanoscale porous structures of shales40. The main idea of MCMC method is to construct a Markov Chain with stable distributions. 2D SEM images are analyzed and used to generate input data to determine the transition probabilities between voids and solids controlling the Markov Chain. MCMC takes into account spatial structural information that identifies all the transition probabilities for a given local training lattice stencil. A complicated multiple-voxel interaction scheme is developed to generate individual realizations which display structural characteristics matching the input data obtained by image analysis of 2D SEM images40. Such reconstruction approach and the structures it generates are referred to as “pore architecture models (PAM)”. For more details of MCMC method, one can refer to the work of Wu et al.40.

Lattice Boltzmann method

Owing to its excellent numerical stability and constitutive versatility, the LBM has developed into a powerful technique for simulating transport processes and is particularly successful in transport processes applications involving interfacial dynamics and complex geometries, such as multiphase flow and porous media flow1633. In the present study, intrinsic permeability and effective Knudsen diffusivity of the reconstructed shales are numerically predicted using the Multiple-relaxation-time LB fluid flow model and LB mass transport model, respectively.

A viscosity-dependent permeability is usually obtained when adopting SRT LBM for simulating fluid flow41. In order to overcome such defect, the multi-relaxation-time (MRT) model has been proposed to simulate fluid flow through porous media41. The MRT-LBM model transforms the distribution functions in the velocity space of the SRT-LBM model to the moment space by adopting a transformation matrix. MRT-LBM (See details in Supporting Information) is adopted in the present study to simulate fluid flow in the reconstructed shales.

There have been no numerical studies to predict effective Knudsen diffusivity of porous shales in the literature. In this work, the LBM is adopted to simulate the Knudsen diffusion in nanoscale porous structures of shales (See details in Supporting Information). Since SRT-LB mass transport model is sufficient to accurately predict the pure diffusion process4243, it is employed rather than MRT-LB one. According to Eq. (1), the Knudsen diffusivity is a function of local pore diameter. Based on the PSD of the reconstructed shales, a reference pore size dp,ref = 25 nm is chosen, and the corresponding relaxation time τg,ref in LBM (See details in Supporting Information) is set as 1.0. Hence, for any local pore with size dp, the relaxation time can be calculated as Through Eq. (16), the effects of PSD on the Knudsen diffusion process inside the nanoscale structures of porous shales are taken into account.

For MRT-LB fluid flow model, for fluid flow three kinds of boundary conditions are adopted: no-slip boundary condition at the fluid-solid interface inside the domain, pressure boundary condition at the inlet and outlet of the computational domain (x direction), and periodic boundary conditions for the other four boundaries (y and z directions). For methane diffusion, four types of boundary conditions are used: the no flux boundary condition at the fluid-solid interface inside the domain, concentration boundary condition at the inlet and outlet of the computational domain (x direction), and periodic boundary conditions for the other four boundaries (y and z directions). The corresponding boundary conditions in the LB framework can be found in Supporting Information.

When the simulation of fluid flow reaches steady state using the MRT-LBM fluid flow model, the intrinsic permeability is calculated using Eq. (8), where the flow flux Jd is expresses as with as the volume-averaged fluid velocity along flow direction.

The effective Knudsen diffusivity DK,eff along the x direction is calculated by where Cin and Cout is the inlet and outlet concentration, respectively. means local value at x = L.

Validation of our numerical models can be found in Supporting Information.

Author Contributions

L.C. wrote the main manuscript and prepared figures 1–6. J.Y. and L.Z. helped to reconstruct the structures of shales. H.S.V., Q.J.K. and W.Q.T. supervised the theoretical analysis and writing. All authors reviewed the manuscript.