Study of Gas Flow Characteristics in Tight Porous Media with a Microscale Lattice Boltzmann Model

To investigate the gas flow characteristics in tight porous media, a microscale lattice Boltzmann (LB) model with the regularization procedure is firstly adopted to simulate gas flow in three-dimensional (3D) digital rocks. A shale digital rock and a sandstone digital rock are reconstructed to study the effects of pressure, temperature and pore size on microscale gas flow. The simulation results show that because of the microscale effect in tight porous media, the apparent permeability is always higher than the intrinsic permeability, and with the decrease of pressure or pore size, or with the increase of temperature, the difference between apparent permeability and intrinsic permeability increases. In addition, the Knudsen numbers under different conditions are calculated and the results show that gas flow characteristics in the digital rocks under different Knudsen numbers are quite different. With the increase of Knudsen number, gas flow in the digital rocks becomes more uniform and the effect of heterogeneity of the porous media on gas flow decreases. Finally, two commonly used apparent permeability calculation models are evaluated by the simulation results and the Klinkenberg model shows better accuracy. In addition, a better proportionality factor in Klinkenberg model is proposed according to the simulation results.

recent years, LBM has also been adopted to investigate gas flow in microscale pores 10 . Compared with MD and DSMC methods, LBM is more efficient. So it can be used to simulate gas flow in relatively large porous media. However, although LBM can be obtained by discretizing the Boltzmann equation, it omits the high order terms. So it is the approximation of the Boltzmann equation and it is not so accurate as the MD and DSMC methods. With several years' development, LBM is considered to be accurate enough to do the microscale gas flow simulations now.
Since LBM is a N-S equation solver, it can be used to simulate gas flow in the continuum flow region by adopting the non-slip boundary condition, or in the slip flow region by adopting the slip boundary condition [11][12][13] . With the increase of Knudsen number (Kn > 0.1), the non-continuum effect becomes more pronounced and the N-S equation is no longer applicable. To extend LBM to simulate gas flow beyond the slip flow region, many efforts have been made, mainly including two ways. One is building the higher order LB models 14,15 , and the other is adopting the effective relaxation time modified by Knudsen number [16][17][18][19] . As the latter one is accurate without losing the advantages of traditional LB model, it is more popular. However, Suga et al. 20 pointed out that this kind of LB model can be adopted to simulate gas flow in microscale channels or tubes, but it cannot be used to simulate gas flow in microscale porous media because of the inadequate symmetry of the discrete model. To overcome this defect, the regularization procedure must be considered [20][21][22] .
In recent years, LBM has been adopted to simulate gas flow in shale gas reservoirs. Fathi et al. firstly introduced LBM into shale gas flow simulation in 2012. They adopted the bounce back boundary condition and Langmuir slip boundary condition to investigate the Klinkenberg slip phenomenon 23 and gas flow characteristics in inorganic and organic micro-channels 24 , respectively. However, their simulation results were questionable because the slip velocities were generated by the discretization error of their models 25 . Since then, more and more researchers adopted LBM to study gas flow characteristics in shale rocks. Zhang et al. 26 adopted the microscale LB model to simulate gas flow in micro-capillary and investigated the slip effects and permeability under different Knudsen numbers. Ning et al. 27 introduced the adsorption effect into the microscale LB model and adopted it to simulate shale gas flow in organic channels and two-dimensional porous media. Ren et al. 25 proposed a microscale LB model considering the effects of surface diffusion, gas slippage and adsorbed layer, and adopted it to investigate the effects of surface diffusion and gas slippage on gas flow in micro-channels. According to our literature review, most of the researchers adopted the microscale LB model to investigate shale gas flow characteristics in simple channels or tubes, because the slip boundary conditions are difficult to conduct on the complex solid boundaries in porous media, especially the 3D porous media. Although some other researchers also adopted LBM to simulate gas flow in microscale porous media [28][29][30] , the geometries of their physical models are simple compared with the real rocks. In addition, they didn't take the regularization procedure into account in their models. So their models are unsuitable for gas flow simulation in microscale porous media under relatively high Knudsen numbers. All the researches mentioned above adopted the common microscale LB model, in which the relaxation time is determined by Knudsen number and the slip boundary condition is adopted on the solid boundaries. In recent years, Chen et al. 31,32 proposed a different LB model to consider the microscale effect based on the dusty gas model. They simulated the viscous flow and Knudsen diffusion separately with different LB models and the total mass flux can be considered as a combined result of viscous flow and Knudsen diffusion. Although this model can be used to simulate microscale gas flow in real digital rocks under any Knudsen numbers, it is not based on the kinetic theory and it could not obtain the detailed velocity distributions in the porous media. As this work mainly focuses on the gas flow characteristics in microscale porous media, the common microscale LB model is adopted here.
In this work, the common microscale LB model is firstly combined with the digital rock technology to investigate gas flow characteristics in real tight rocks. The regularization procedure is introduced into the model to enable it to simulate gas flow in porous media under high Knudsen numbers, which could not be simulated in others' work [27][28][29][30] . Because of the random pore size distribution, the characteristic length is not a constant. The local characteristic length is introduced into the model which was always ignored in others' work [27][28][29][30] . In addition, the diffuse reflection boundary condition is adopted to deal with the random solid boundaries. The effects of pressure, temperature and pore size on microscale gas flow are investigated firstly at the pore scale and the influencing mechanisms are analyzed. Then the gas flow characteristics in digital rocks under different Knudsen numbers are studied and some new phenomena have been found. Finally, the simulation results are adopted to verify the accuracy of two commonly used apparent permeability calculation models and a modified Klinkenberg model is proposed according to the simulation results.

Results
The microscale LB model is adopted to simulate gas flow in digital rocks. The effects of pressure, temperature and pore size on microscale gas flow are investigated and the gas flow characteristics in tight porous media under different Knudsen numbers are analyzed.
Digital rocks. Two different digital rocks are used for the microscale gas flow simulations, as shown in Fig. 1(b,c). The shale digital rock is reconstructed by Markov Chain Monte Carlo (MCMC) method based on the SEM scanning image shown in Fig. 1(a) and it is adopted to study the effect of pressure and temperature on microscale gas flow. As this work mainly focuses on the microscale gas flow characteristics in tight porous media, the effect of kerogen is not taken into account here. To get the digital rock with the presence of kerogen, one can refer to the reference 33 . The sandstone digital rock is obtained by micro-CT scanning method. Although the real resolution of this digital rock is determined, by selecting different resolutions manually, we can investigate the effect of pore size on microscale gas flow. For more information about the MCMC method and micro-CT scanning method, one can refer to the references 33,34 .
Firstly, the traditional LB model not considering the microscale effect is adopted to simulate pressure driven gas flow in these two digital rocks. The top is inlet and the bottom is outlet. Pressure boundary condition is adopted in the z direction. For the x and y directions, periodic boundary condition is adopted. When the simulations reach the steady state, the intrinsic permeability and equivalent pore size of each digital rock are calculated based on the simulation results. The intrinsic permeability can be obtained based on Darcy's law and the equivalent pore size can be obtained by the following equation: where H equivalent is the equivalent pore size; k intrinsic is the intrinsic permeability; φ is the porosity. The equivalent pore sizes of the two digital rocks are 2.84 and 3.55 in lattice unit, respectively. They are adopted to calculate the equivalent Knudsen numbers in the following simulations. Then pressure driven gas flow in these two digital rocks is simulated by the microscale LB model to investigate the gas flow characteristics in tight porous media.
The effect of pressure on microscale gas flow. To investigate the effect of pressure on microscale gas flow, methane flow in the shale digital rock under different pressures is simulated by the microscale LB model. The temperature is 373 K and the inlet-outlet pressure differences are all 0.001 MPa. Although the pressures in shale gas reservoirs are always very high, the microscale effect is more obvious at low pressure. To investigate the effect of pressure on microscale gas flow more thoroughly, a wide range of pressure is adopted here. In the following simulations, the outlet pressures are 0. When the simulations reach the steady state, the volume flux across the digital rock can be obtained and the apparent permeability k a can be calculated based on Darcy's law. The effect of pressure can be reflected by k r which is defined as k r = k a /k intrinsic . The simulation results are shown in Fig. 2. Pressure has a significant effect on gas flow in microscale porous media. When the pressure is very high, the apparent permeability of the shale digital rock is similar to the intrinsic permeability. With the decrease of pressure, the apparent permeability increases gradually. Especially when the pressure is very low (smaller than 1.0 MPa), the apparent permeability increases dramatically with the decrease of pressure.  This can be explained by Fig. 3 which shows the velocity distributions in the shale digital rock under three different outlet pressures. As can be seen from Fig. 3, pressure has a significant effect on gas flow in the digital rock. When the pressure is high, gas velocity in the digital rock is very small; with the decrease of pressure, gas velocity in the digital rock increases dramatically. This is caused by the different gas flow mechanisms under different pressures. When the pressure is high, the gas molecules are close to each other. The molecular mean free path is small compared with the pore size. Therefore, the intermolecular collisions are much more frequent than the collisions between gas molecules and solid walls. There are almost no slip velocities on the solid walls. As a result, the velocity in the digital rock is small and the apparent permeability is similar to the intrinsic permeability. With the decrease of pressure, the molecular mean free path increases. The collisions between gas molecules and solid walls become increasingly obvious compared with the intermolecular collisions. Therefore, the slip velocities on the solid walls increase and the velocity in the whole digital rock increases. The apparent permeability increases accordingly.
The effect of temperature on microscale gas flow. The shale digital rock is also adopted to investigate the effect of temperature on microscale gas flow. The outlet pressure is 20.0 MPa and the inlet pressure is 20.001 MPa. The temperatures are 298 K, 323 K, 348 K, 373 K, 398 K, 423 K, 448 K and 473 K, respectively.
When the simulations reach the steady state, the results can be obtained, as shown in Fig. 4. The apparent permeability increases with the increase of temperature. As the intrinsic permeability remains unchanged, k r increases with the temperature. This is because the molecular mean free path is related to the temperature. With the increase of temperature, the kinetic energy of the gas molecules increase and they move faster. Then the collisions between gas molecules and solid walls become stronger. To keep the same pressure, the gas density must decrease. Thus the molecular mean free path increases and the collisions between gas molecules and solid walls become more obvious. Therefore, the microscale effect increases and the apparent permeability increases.
The effect of pore size on microscale gas flow. The sandstone digital rock is adopted to investigate the effect of pore size on microscale gas flow. Methane is adopted in the simulations and the temperature is 373 K. The physical size of the digital rock can be changed by selecting different resolutions. The following resolutions are adopted here: 3.9 nm/pixel, 5.5 nm/pixel, 7.8 nm/pixel, 11.0 nm/pixel, 15.6 nm/pixel, 22.0 nm/pixel, 39.0 nm/  pixel, 55.0 nm/pixel, 78.0 nm/pixel, 110.0 nm/pixel and 156.0 nm/pixel. As analyzed in above simulations, the microscale effect decreases with the increase of pressure. To magnify the effect of pore size, a relatively low pressure is adopted here. In the following simulations, the outlet pressures are all 1.0 MPa and the pressure gradients are all 6.25 × 10 −7 MPa/nm.
When the simulations reach the steady state, the results can be obtained, as shown in Fig. 5. The equivalent pore size is set as the x-axis. As can be seen from Fig. 5, the apparent permeability is always higher than the intrinsic permeability, and with the decrease of pore size, both the intrinsic permeability and apparent permeability decrease, but k r increases. When the pore size is large, the molecular mean free path is small compared with the characteristic length. So the microscale effect is not obvious and the apparent permeability is similar to the intrinsic permeability. As the pore size decreases, the collisions between gas molecules and solid walls become more important. So the microscale effect increases and the slip velocities on the solid walls increase. Then the difference between the apparent permeability and intrinsic permeability (k r ) increases.

Gas flow characteristics in digital rocks under different Knudsen numbers.
As stated in the introduction, Knudsen number is the characteristic parameter for gas flow in microscale porous media. So the Knudsen numbers of gas flow in the digital rocks under different conditions are calculated and the corresponding gas flow characteristics are analyzed. Figure 6 shows the proportions of volume flux in pores of different sizes under different Knudsen numbers. The pore volume distributions are also shown in this figure. As can be seen from Fig. 6, in small pores, the proportion of pore volume is higher than that of volume flux; while in large pores, the proportion of volume flux is higher. So the velocities in small pores are smaller than those in large pores, as shown in Fig. 7. With the increase of Knudsen number, the proportion distribution curve of volume flux in Fig. 6 moves to the lower left. So the relative volume fluxes in small pores increase while those in large pores decrease. As a result, the dimensionless average velocities in small pores increase while those in large pores decrease, as shown in Fig. 7. This means that  The following statements can explain the above phenomena. Gas flow in the porous media is motivated by the intermolecular collisions and collisions between gas molecules and solid walls. When Knudsen number is very small, gas flow is mainly dominated by the intermolecular collisions. There exist very small slip velocities on the solid walls in both small pores and large pores. In large pores there are more intermolecular collisions, while in small pores there are less intermolecular collisions. So the gas flow ability in large pores is much higher than that in small pores. The proportion of volume flux and dimensionless velocity in large pores are higher and the effect of heterogeneity of the porous media is obvious. With the increase of Knudsen number, the effect of the collisions between gas molecules and solid walls becomes increasingly obvious, so the slip velocities on the solid walls increase and the gas flow resistance decreases. In addition, the Knudsen numbers in small pores are higher than those in large pores because of their smaller sizes. So the microscale effect is more obvious and the gas flow resistance decreases more in small pores. Therefore, the difference of gas flow abilities in large pores and small pores decreases. Then gas flow in the whole digital rock becomes more uniform and the effect of heterogeneity of the porous media on gas flow decreases.
Evaluation of two apparent permeability calculation models. According to the simulation results, in tight porous media, such as shale gas reservoirs, because of the microscale effect caused by the ultra small sizes of pores, the permeability is no longer a constant. The apparent permeability is always used in macroscale shale gas reservoir simulations. Klinkenberg model 35 and Beskok-Karniadakis model (B-K model) 36 are two commonly used apparent permeability calculation models. For 3D fluid flow in porous media, they have the following formats.
Klinkenberg model: where c is the proportionality factor. Klinkenberg pointed out that the value of c seems to be slightly less than 1.0 35 . However for most researchers, they usually choose c = 1.0 37 . In this work, we also choose c = 0.8 for comparison. B-K model: − . LBM is adopted to evaluate the accuracy of these two models. Figure 8 shows the comparison of the calculated results of the two models with the simulation results of LBM. As shown in Fig. 8, when the Knudsen number is small, k r increases slowly with the increase of Knudsen number. Both models can well describe the microscale effect. But when Knudsen number gets higher, k r increases quickly as Knudsen number increases. The B-K model will obviously overestimate the microscale effect, while the Klinkenberg model gives a better prediction. In addition, our simulation results suggest that the proportionality factor c in Klinkenberg model should be modified to 0.8.

Kn
After getting the modified Klinkenberg model, our simulation results can be used in the following way. For any porous medium, it can be equivalent to the capillary bundles. Then the equivalent pore size can be obtained according to Eq. (1), and the corresponding Knudsen number can be calculated by Eq. (14). After obtaining the Knudsen number, the modified Klinkenberg model can be used to calculate the apparent permeability of the tight porous medium, which is a critical parameter in the macroscale numerical simulation of tight gas or shale gas reservoirs.

Discussion
In this work, a microscale LB model with the regularization procedure is adopted to simulate gas flow in 3D digital rocks. The diffuse reflection boundary condition is adopted to deal with the random solid boundaries and the local characteristic lengths of pores with different sizes are also introduced into the model. A shale digital rock and a sandstone digital rock are adopted to investigate the gas flow characteristics in microscale porous media. Gas flow simulations in the shale digital rock under different pressures and temperatures are conducted to investigate the effects of pressure and temperature on microscale gas flow. And the sandstone digital rock is used to study the effect of pore size on microscale gas flow by selecting different resolutions. The simulation results show that with the decrease of pressure or pore size, or with the increase of temperature, the difference between apparent permeability and intrinsic permeability increases. The decrease of pressure and the increase of temperature will increase the molecular mean free path while the decrease of pore size will decrease the characteristic length, all will make the collisions between gas molecules and solid walls become more obvious compared with the intermolecular collisions. As a result, the slip velocities on the solid walls increase and the apparent permeability increases. In addition, the results show that gas flow characteristics in the digital rocks under different Knudsen numbers are quite different. With the increase of Knudsen number, the dimensionless average velocities in small pores increase while those in large pores decrease. Gas flow in the digital rocks becomes more uniform and the effect of heterogeneity of the porous media on gas flow decreases. Finally, the Klinkenberg model and B-K model for calculating the apparent permeability of tight porous media are evaluated by the LBM simulation results. According to the simulation results, we suggest the Klinkenberg model in macroscale shale gas reservoir simulations because of its better accuracy. In addition, a better proportionality factor (c = 0.8) in Klinkenberg model is proposed according to the simulation results.
Kerogen is very common and important in shale gas reservoirs. However, this work mainly focuses on the gas flow characteristics in tight porous media. In our future work, we will investigate the adsorption effect of kerogen on microscale gas flow.

Methods
Lattice Boltzmann method. In this work, the 3D microscale LB model with the D3Q19 discrete velocity model 38 is adopted to simulate gas flow in digital rocks. The basic evolution equation with the Bhatnagar-Gross-Krook (BGK) collision approximation is shown as follows. Unless otherwise stated, all units in this paper are lattice units.
where f α is the density distribution function of α direction; α = 0, 1, 2, … , 18; r is the spatial location of the particles; e α is the velocity of α direction; t is time; δ t is time step; τ is the relaxation time; F α is the force term; α f eq is the local equilibrium distribution function of α direction: where τ e is the effective relaxation time considering the effect of Knudsen layer; ψ(Kn) is the modification function and we set ψ(Kn) = 1/(1 + 2Kn) 19 ; N = H/Δ x means the number of lattices occupied by the characteristic length, where H is the characteristic length, m, and Δ x is the size of one lattice, m.
As the solid boundaries in shale are quite rough, the diffuse reflection boundary condition is very appropriate for such boundaries. The discrete format of the diffuse reflection boundary condition in LBM is 11,13 :  where n is the inward unit normal vector; u w is the wall velocity; the subscript w means the solid walls; ′ α f is the distribution function after streaming. As long as the inward normal direction can be obtained, the diffuse reflection boundary condition can be adopted.
Regularization procedure. To apply the microscale LB model into gas flow simulation in porous media, the regularization procedure must be introduced 20,39 . With the regularization procedure, the evolution equation becomes: where α ∼ f noneq is the regularized distribution function and can be expressed as:  (2) is the second Hermite polynomial 21 . For more information about the regularization procedure one can refer to the references 20,39 . Local Knudsen number. According to the kinetic theory, for hard sphere molecules, the molecular mean free path is 16 :  where m is the molecular weight and m = 2.658 × 10 −26 kg for methane; σ is the molecular diameter and σ = 3.8 × 10 −10 m for methane; ρ is gas density, kg/m 3 . According to the definition of Knudsen number, it can be expressed as: 2 For gas flow in porous media, the characteristic length is the pore size. Because of the non-uniform pore size distribution, H is not a constant in porous media. The pore sizes in different positions can be obtained by the method in reference 40 . For a certain type of gas, m and σ are constants. As ρ is determined by the pressure, Knudsen number is related to P and H. Before the simulation, the reference characteristic length H ref  To study the effect of grid refinement, three slab models with different grid numbers are adopted here: slab model 1 with N x × N y × N z = 7 × 3 × 5, slab model 2 with N x × N y × N z = 12 × 3 × 10 and slab model 3 with N x × N y × N z = 22 × 3 × 20. The left and right nodes are solid nodes and the nodes in the middle area are pore nodes. Periodic boundary condition is adopted in y and z directions and a uniform pressure gradient of 5 × 10 −5 is applied in z direction. Gas flow in the slab models under four different Knudsen numbers are simulated by our model and the simulation results are compared with those of DSMC method 20 , as shown in Fig. 9. The results of the LB model match well with those of DSMC method, which verifies the accuracy of this model to simulate gas flow in different flow regions. In addition, the simulation results under different grid refinement almost fall on the same line, which demonstrates that the simulation results of this LB model is independent of grid refinement.
To further verify the accuracy of this model to simulate gas flow in microscale porous media, force driven gas flow in microscale porous media is simulated by this model and the simulation results are compared with those of MD method 41 . In addition, the microscale LB model not considering the regularization procedure is also adopted to do the same simulation. The physical model is a parallel slab model with a square cylinder in the middle 20,41 , as shown in Fig. 10. The resolution is 1.0 nm/pixel and methane is adopted in the simulation. The temperature is 373 K and the pressure is 2.413 MPa. Then the corresponding Knudsen number is 0.11. The driven acceleration is 5 × 10 −5 in x direction. The simulation results are shown in Fig. 11. The results of the LB model with regularization procedure match well with those of MD method, but the LB model without regularization procedure will generate unphysical velocity distributions. This demonstrates that for LBM, the regularization procedure is essential for gas flow simulation in microscale porous media. As the models adopted in others' research [27][28][29][30] did not take the regularization procedure into account, they could not be adopted to simulate gas flow in porous media under relatively high Knudsen numbers.