AllAboutShale.com

1Department of Energy and Resources Engineering, College of Engineering, Peking University, Beijing 100871, China.

2ERE & BIC-ESAT, College of Engineering, Peking University, Beijing 100871, China.

Given the complex nature of the interaction between gas and solid atoms, the development of nanoscale science and technology has engendered a need for further understanding of gas transport behavior through nanopores and more tractable models for large-scale simulations. In the present paper, we utilize molecular dynamic simulations to demonstrate the behavior of gas flow under the influence of adsorption in nano-channels consisting of illite and graphene, respectively. The results indicate that velocity oscillation exists along the cross-section of the nano-channel, and the total mass flow could be either enhanced or reduced depending on variations in adsorption under different conditions. The mechanisms can be explained by the extra average perturbation stress arising from density oscillation via the novel perturbation model for micro-scale simulation, and approximated via the novel dual-region model for macro-scale simulation, which leads to a more accurate permeability correction model for industrial applications than is currently available.

The transport mechanism of gas through nanopores is a fundamental concern in many fields, such as nanofluidics, fluid mechanics, material science, chemical engineering, and petroleum engineering, which has attracted considerable attention. For example, with the successful development of unconventional oil and gas1–3, obtaining a better understanding of the transport mechanisms in the nanopores of tight/shale matrix is important for predicting long-term gas production; on the other hand, the discovery of the rapid transport of gases in carbon nano- tubes and graphene membranes4–6 means that these materials have the potential to function as high-efficiency and low-cost materials for hydrogen or methane storage. All of these applications involve the movement of gas molecules through highly confined space where the gas molecules are strongly affected by gas–solid interaction in addition to purely gas–gas interaction7,8. Because the strong density oscillation of fluid atoms near the fluid–solid interface is a universal phenomenon9, the gas–solid interaction can be described with the help of the concept of physical adsorption10–12, which is a well-developed field in physical chemistry. The effects of adsorption cannot be neglected, since the surface-to-volume ratio is very high, and the critical dimension is comparable to the thickness of the adsorption layer in a nanoscale system9. Therefore, the mechanisms of gas transport under the influence of adsorption in nanopores need to be understood in detail.

Research on gas flow in confined spaces has a long history, dating back to the works of Maxwell13 and Knudsen14. They presented the phenomenon of gas slippage near a solid wall in the context of the gas–solid collision, and the models were subsequently refined by others15. Then, Klinkenberg addressed the enhanced gas flow in a porous medium and its effect on permeability16. As the adsorption phase contributes to transport of gases to a great extent in some cases, surface diffusion was introduced to explain the mass flow in excess to that predicted by Knudsen’s model, and on this bases more models were proposed and improved17–21. On the other hand, the adsorbing capacity were often modeled by various adsorption isotherms, such as the Langmuir model for monolayer adsorption22 and the Brunauer-Emmett-Teller (BET) model for multilayer adsorption23.

In general, the above models are expressed in the diffusional framework with a concentration dependent diffusivity. However, another way of investigating the transport mechanisms is based on the Navier-Stokes (N-S) equation framework with idealized models, such as the channel flow and the tube flow24. The flow under nanoscale exhibits non-continuum effects25, and its flow regimes can be classified as continuum flow, slip flow, transition flow, and free molecular flow26 by Knudsen number (Kn = λ/l, where λ is the molecular mean free path, and l is the characteristic length). As the characteristic length approaches the order of the molecular mean free path, the Kn of the gas flow in nanopores is in a similar range to that of the rarefied gas flow with large λ under low pressure. Then, the law of similarity can be applied, which is in the slip flow or transition flow regime24 with a low Reynolds number and a low Mach number. Further, the continuum model with corrections, such as the N-S equation, has been utilized to describe the gas flow in nanopores, which were proved through direct-simulation Monte Carlo (DSMC) and linearized Boltzmann equations24. After that, a variety of slip boundary conditions24,27–29, such as the extended N-S equation30 and empirical velocity profiles25,31,32 were developed and widely applied. In general, because of the weak contribution of adsorption in rarefied gas flow under ultra-low-pressure conditions, the adsorption effect is often neglected in these models. However, high pressure, corresponding to dense gas, is also very common in industrial applications, and the impact of adsorption may be significant on the gas transport process. In addition, transport behaviors are often described as a variety of mechanisms, such as Knudsen diffusion, slip flow, viscous flow, and/or surface diffusion, and the dominant mechanisms may vary under different conditions. Subsequently, some phenomenological models33–36 in a superposition form of the above mechanisms were proposed, and were utilized to interpret experimental results. However, the rationality of some phenomenological models is under debate, due to their insufficiency to address the coupled effect.

With the development of atomistic simulation methods and computing capacity, the Molecular Dynamics (MD) simulation has become a powerful tool for investigating the micro-scale mechanisms. Various MD simulation techniques, such as equilibrium molecular dynamics (EMD), non-equilibrium molecular dynamics (NEMD), and dual control volume grand canonical molecular dynamics (DCV-GCMD), were designed with different strengths and applied to verify theoretical models and discover novel phenomena37–39. For example, the large slip lengths in graphite slit-like pores and carbon nanotubes were presented, respectively, and an interfacial friction model was proposed, which indicated that the behavior was influenced by the surface corrugation and the fluid–wall interaction strength40,41; another friction-based model demonstrated that the combination of viscous flow and momentum exchange at the pore wall governs the transport42; the oscillator model, on the other hand, offered an exact theory for low density gas flow8; and sub-continuum mass transport of condensed hydrocarbon in disordered porous carbon, with pore sizes ranging from 0.3 to 1.3 nm, was investigated43. Overall, various phenomena have been discovered and modeled independently, however, the understanding of the mechanisms of gas transport coupled with adsorption in nanopores remains limited, and a tractable numerical model for large-scale simulation under an extensive range of pore size and pressure has yet to be proposed.

In the current paper, we present a theoretical study of the impact of adsorption on gas transport in nanopores, which aims at analyzing the coupled behavior in detail and proposing numerical models for the micro-scale phenomena description and macro-scale flow simulation. First, we perform MD simulations to show the special mass flow properties in the nano-channel made of illite and graphene, respectively. To account for the velocity behavior under the strongly heterogeneous gas density distribution along the cross-section of the channel, we propose a perturbation model, which decomposes the stress term into viscous stress and average perturbative stress. As the monolayer adsorption assumption remains valid for the majority of gas transport cases, we present a dual-region model in order to approximate the mass flow and demonstrate its application of permeability correction. Meanwhile, other classical models and frequently utilized models are compared within this paper.

We constructed the nano-channel (slit-like pore) structures with illite44 and graphene, respectively. Illite represents a common layered clay mineral with significant adsorption in the shale matrix. Graphene, on the other hand, represents the advanced material for various applications including hydrogen storage45 and methane storage46. In the case of illite, the silicon-oxygen layer was assigned as the pore surface, which can represent many frequently utilized solid structures. The center of the top layer’s atoms was set as the position of the solid wall. The effective pore width was smaller than the assigned value and may change according to the atoms’ properties. Two kinds of NEMD simulations were conducted for isothermal methane flow as a simple fluid at room temperature (T = 298 K). One was based on the reflecting particle method (RPM), which is easy to apply and the driving mechanism of which is close the actual pressure-driven flow38. As the pressure gradient in the highly confined space could not be well-defined or monitored, the other one, the external force method (EFM), was applied in the next section in order to accurately determine the pressure gradient. Supplementary Fig. S1 presents the comparison of RPM and EFM for this case. Additional details pertaining to the MD simulation methods and settings can be found in Supplementary Methods.

Figure 1 demonstrates the velocity profile and the accumulative mass flow across the 10-nm-wide channel under various pressure conditions. In the illite cases, the velocity profile can be divided into three regions. The first region is adjacent to the solid atoms where the gas atoms are strongly repulsed by the wall and no atoms or velocity appears. Therefore, this region can be eliminated by defining the effective width. The second region is around the adsorbed region where complex phenomena are present. The slip velocity is very small, and vanishes in the high-pressure case. The velocity gradient is smaller around the potential well, where a density peak exists, than that in the other regions. The third region is the inner region, the transport behavior of which mainly follows the classical parabolic type but with different curvatures caused by the rarefaction effect in the nanopores24. As the pressure decreases, the velocity profile transforms from the parabolic type to the plug type, due to the viscosity reduction as the Knudsen number increases. Additionally, a velocity inflection point always exists between the inner region and the adsorbed region. On the whole, due to the small velocity at the boundary, the adsorbed region contributes very little to the total mass flow, especially in the high-pressure case, which is consistent with the macroscopic behavior of shale gas production47. However, although the graphene cases also have three regions, a significant phenomenon arises. A large slip velocity, principally due to the ultra-smooth surface of the graphene, appears in the adsorbed region. Since the density peak is in the adsorbed region and the difference in velocity between the adsorbed region and the inner region is minute, the adsorbed region contributes substantial mass flux, even exceeding 40% in the relatively low-pressure case. This nice property has been explored in order to determine its potential in fast gas storage, the diffusivity of which was shown to be several orders of magnitude higher than that of normal nano-porous media4. The relative contribution of the adsorbed region decreases as the pressure increases, because the density ratio between the adsorbed region and the inner region decreases. Taken as a whole, the velocity profile fits a flat type, but with a similar trend under various pressures to that of the illite case (Fig. 1b inner plot).

Figure 1. Transport behavior under various pressure conditions in a 10-nm-wide channel from MD simulations. Velocity (dot) and accumulative mass flux (line) in the illitic channel (a) and the graphene channel (b) along the cross-section; inset: velocity profile after zooming in on the y-axis in graphene channel cases. The corresponding gas atoms number density distribution is shown in Supplementary Fig. S2. The pressure is the average value in the central free space with uniform density distribution. The velocity is normalized by the average velocity across the channel.

The molecular transport at the micro-scale can be easily decomposed into the streaming velocity and the thermal motion velocity or speed48,49, namely, u = u + u′. When we monitor the average streaming velocity in a specific bin, it is implicitly assumed that the velocity is a density-weighted average velocity, because the average density is also a statistical value and the heterogeneous density distribution exists in a nanopore. This definition can be expressed as u = ρu/ρ, where u is the velocity of the molecule, ρ is the density, and ρu and ρ denote the long-time average of momentum and of density in a specific bin, respectively. We use ux in the figures to represent the x-component of the average streaming velocity. If we take a long-time accumulation for the thermal velocity, the value should be zero. We can, therefore, treat the thermal velocity as a perturbative velocity. If we assume that the continuum theory is valid in this condition from a statistical perspective, it is reasonable to apply the Navier-Stokes equation framework30. In analogy to the turbulence theory, which often employs molecular transport as an analogy in return50, the two theoretical systems maintain the same decomposition scheme. Therefore, we propose a perturbation model to describe the flow behavior in nanopores:

where t represents time; p stands for the average pressure; τij signifies the stress term; u′ is the perturbation velocity; −ρu′i u′j is the extra average perturbation stress term; and i, j denotes the axis-direction. This model retains the same form as the compressible turbulence flow based on the mass-weighted average concept50. In order to determine a physical expression of the perturbation term, we can apply a simplified kinetic theory model for approximation. Then, if we define Jij= τij−ρui′uj′ for the stress term, the momentum transport across the unit area can be denoted as:

where <v> is the thermal motion speed. The derivation is given in Supplementary Methods. The first term can be expressed as the viscosity term in the relationship of η = 1/3ρλ<v> , where the viscosity η can be assumed to be independent of the density in gas flows49 and kept constant under certain temperature. For a Newtonian fluid, the stress due to viscosity can be expressed as:

Additionally, we can assume that:

Figure 2. Velocity fitting results from the perturbation model in the illitic channel. (a) Low-pressure case (5.06 MPa) in the 10-nm-wide channel. (b) High-pressure case (26.17 MPa) in the 10-nm-wide channel. (c) Low-pressure case (7.44 MPa) in the 2-nm-wide channel. (d) High-pressure case (27.74 MPa) in the 2-nm- wide channel. The pressure is the average value in the central free space with a uniform density distribution in (a,b). The pressure cannot be well defined in a 2-nm-wide channel, therefore, the pressures in (c,d) are representative values, just for comparison, from the average density via the Peng-Robinson equation of state. The velocity is normalized by the average velocity across the entire channel from the MD simulation. The post- processing methods for pressure and velocity are also applied in all other figures presented in this work.

Figure 3. Velocity fitting results from the perturbation model in a 10-nm-wide graphene channel. (a) Low- pressure case (5.16 MPa). (b) High-pressure case (26.11 MPa).

where c0 and c are scale coefficients to be fitted. Then, equation (4) indicates that the velocity oscillation may result from the extra average perturbation stress due to density oscillation.

According to the configuration of the model in the MD simulations, the equations can be reduced to:

Equation (5) is solved via the finite difference method with the no-slip boundary condition at the wall atoms’ center position and the symmetric boundary at the center of the pore. The density profile is the most important input information.

To verify the perturbation model, MD simulations were performed by imposing the pressure gradient through external force. Figure 2 demonstrates the velocity profile results from the illite cases with various channel widths. The perturbation model can thoroughly capture the features of velocity profiles, accurately describing the inflection point which is more significant in the 2nm cases. The scale coefficient c is 0.0370, 0.0105, 0.0175, and 0.0194 in the cases demonstrated in Fig. 2a–d, respectively. The viscosity was calculated with the LBC method and corrected with the rarefaction effect24. Because of the smooth surface of graphene, the gas density profile always maintains a very sharp peak adjacent to the wall (see Supplementary Fig. S2), which often leads to numerical instability at the boundary. However, if an average velocity in the adsorbed region near the wall is assigned to the boundary, the perturbation model can also capture the velocity properties in both the adsorbed region and the inner region. The scale coefficient c is 2.4×10−4 and 5.8 ×10−4 in the cases demonstrated in Fig. 3a,b, respectively. The 2nm graphene cases present flat velocity profiles (Fig. 4), because the slip velocity is very large and the effects of viscosity and density oscillation are covered by noise. In addition, the mass flux contribution corresponds to the density profile. Therefore, it is not necessary to fit the model in this case.

It should also be noted that the perturbation model offers an explanation of the slippage mechanism. The slip velocity or the small slip region near the wall can be described well by this model in the illite case (see Fig. 2c,d).

Figure 4. Transport behavior in a 2-nm-wide graphene channel under various pressure conditions. Velocity (line-dot) agrees with the flat type, and the oscillation is covered by noise; the accumulative mass flux curve (dashed line) is only affected by the density distribution.

In the graphene case, the large slip velocity may be related to the large density gradient near the wall, which is much larger than in the illite case (see Supplementary Fig. S2). Although no wall potential term appears in this model, the density profile is an effective indirect parameter containing not only the effect of the solid wall but also the interaction between the gas–gas molecules. In particular, this model avoids the problem encountered in the continuum method in which the solid wall’s effect is unmanageable. Furthermore, this model makes it possible to apply density profile information obtained from other methods, e.g., grand canonical molecular dynamics (GCMC) and density function theory (DFT), to gas flow simulation.

The perturbation model can describe the velocity profile very well at the micro-scale and offer a basis for gaining insights into and even explaining the mechanisms. However, applying the perturbation model in macro-scale simulations is challenging because it is difficult to obtain the density distribution within each of the pores for the macro-scale pore network model. From the MD results and the perturbation model, we found that the behavior of the velocity profile differs significantly between the adsorbed region and the inner region because of the density oscillation. We can assume two surfaces, corresponding to the positions of the large density valleys: one surface, which may have slip velocity, is between the solid wall and the adsorbed region, whereas the other surface, which is virtual in nature, is between the first adsorbed layer and the inner region. Although the second density peak appeared in the high-pressure case, it simply enhanced the velocity oscillation around the surface, and no obvious oscillation presented between the second density peak and the central region. In other words, only the monolayer adsorption atoms maintained a significant effect on gas flow. Therefore, it is reasonable to propose simplified models for macro-scale simulation to approximate this flow behavior.

Here, we propose a dual-region model (D-R model) for mass flow in nanopores, based on the partition of the adsorbed region and the inner region, the flow regimes of which differ. The density distribution is assumed to be uniform within each region. Additionally, the density of the adsorbed region is easy to determine from experiments, such as the adsorption isotherm test and the specific area test. To approximate the major velocity features in each region, we sacrificed the velocity precision around the virtual surface, which we treated it as a slip velocity. The slip boundary conditions were applied for both surfaces through the respective Knudsen number and the tangential momentum accommodation coefficient (TMAC), which may switch to the no-slip boundary under high pressure24. The viscosities were also corrected with Kn and density in each region, respectively. Therefore, the approximate velocity profile can be expressed as follows by taking the center of the pore as the zero point:

where H is the effective channel width; α and β are the ratios of the density and viscosity between the adsorbed region and the inner region, respectively; w1, w2 and w3 are coefficients in terms of α, β, Kn, TMAC, and the adsorbed region’s width (for expressions, see Supplementary Methods); the width of the adsorbed region is assigned as γσ, where σ corresponds to the distance at which the potential between the two atoms goes to zero, which can be approximately taken as the diameter of a gas atom, and γ is set to 0.8 for methane empirically; the subscripts in and ad denote the inner region and the adsorbed region, respectively; and Zad is the location of the virtual surface. Then, we can obtain the total mass flux as:

Figure 5. Comparisons of different models for velocity and accumulative mass flux in the illitic channel. Upper panel: velocity profile. Lower panel: accumulative mass flux. (a,b) 10-nm-wide, low-pressure (5.06 MPa); (c,d) 10-nm-width, high-pressure (26.17 MPa); (e,f) 2-nm-width, low-pressure (7.44 MPa); (g,h) 2-nm-width, high-pressure (27.74 MPa). The values of TMAC can be found in Supplementary Table S1.

where L • H is the cross area of the channel; ρg is the gas density in the inner region; and h* is the width ratio of the inner region and the entire channel. Other classical models, including the no-slip model and the Beskok-Karniadakis-Trimmer (BKT) model24,27 are introduced for comparison. The BKT model is widely applied for microflows and nanoflows, which are proposed based on rarefied gas flow and fail to account for adsorption. Notably, the dual-region model will reduce to the BKT model in the limiting cases, such as those in which no adsorption effect exists.

Figure 5 demonstrates the comparisons with the no-slip model and the BKT model in the illite case. In case 1 (Fig. 5a,b, width: 10 nm, low-pressure), both the dual-region model and the BKT model agree with the profile very well, whereas the no-slip model underestimates both the velocity and the mass flow. In case 2 (Fig. 5c,d, width: 10 nm, high-pressure), both the BKT model and the no-slip model overestimate the velocity and the mass flow rate. The BKT model always maintains slip velocity, because the bulk Kn is not small enough to obtain a no-slip boundary in such a narrow channel. The no-slip model cannot correct the viscosity in the adsorbed region like the dual-region model does. In contrast, the dual-region model maintains a very small Kn in the adsorbed region so that the slip velocity approaches zero, and can apply viscosity correction with each region’s Kn. However, in case 3 (Fig. 5e,f, width: 2 nm, low-pressure), the no-slip model agrees with the velocity profile very well but underestimates the flow rate, because this model cannot account for the density difference, which has little effect in cases 1 and 2 due to the small volume fraction of the adsorbed region. Although the BKT model could fit the mass flow profile by adjusting the TMAC, this model would have non-physical velocity. Case 4 (Fig. 5g,h, width: 2 nm, high-pressure) resembles case 2 but with a more sizable deviation, since the adsorbed region maintains a larger volume fraction.

In case 5 (Fig. 6a,b, width: 10 nm, low-pressure), because of the large slip velocity, the no-slip model significantly underestimates the velocity and the mass flow, whereas the BKT model and the dual-region model both agree very well with the MD results in the central major part. In addition, the dual-region model has an adjustment in the adsorbed region. However, the BKT model still underestimates the total mass flow because the enhanced density in the adsorbed region is not accounted for. Nevertheless, this deviation decreases, as the density ratio of the adsorbed region and inner region decreases with increasing pressure in case 6 (Fig. 6c,d, width: 10 nm, high-pressure).

Overall, the dual-region model can capture most of the features in nano-channels with the addition of three parameters (α, β, and h*), each of which has a definite physical meaning and can be estimated from theoretical calculations or measured from experiments. However, the classical no-slip model and the BKT model may either overestimate or underestimate the velocity and the mass flow. Although the BKT model can fit the total mass flow in some cases, this model may have a non-physical velocity. Therefore, the dual-region model offers a better approximate model for the macro-scale simulation than other models, since it accounts for the adsorption effect and the complex flow mechanisms.

Figure 6. Comparisons of different models for velocity and accumulative mass flux in a 10-nm-wide graphene channel. Upper panel: velocity profile. Lower panel: accumulative mass flux. (a,b) Low-pressure (5.16 MPa). (c,d) High-pressure (26.11 MPa). The values of TMAC can be found in Supplementary Table S2.

From the macroscopic perspective, the performance of gas transport through a porous medium is often denoted as the apparent permeability51, which can be expressed as the intrinsic permeability multiplied by the correction factor f. The intrinsic permeability is an inherent property of a porous medium, and does not depend on either fluid type or flow-regime. The BKT model can also convert to a correction factor function in terms of Kn, which can be expressed as24:

where a is a rarefaction coefficient for viscosity correction; b is a slip coefficient and often taken as −1 for channel flow and tube flow; ω = 6 in the channel flow and ω = 4 in the tube flow. For the dual-region model, the correction factor function is:

In the channel flow, the intrinsic permeability can be denoted as H2/12 based on the Stokes’ equation. Then, the MD simulation results can be interpreted in the form of the apparent permeability. As shown in Fig. 7a,b, the dual-region model maintains a more effective performance throughout all of the cases than the BKT model. The Langmuir-type adsorption isotherms are applied as input parameters for the apparent permeability calculation (see Supplementary Fig. S3). In addition, if we convert the apparent permeability to transport diffusion coefficient form, the diffusivity of the graphene channel maintains a similar value to the (10, 10) carbon nanotube4 (about 1×10−4m2s−1).

Figure 7. Apparent permeability of the nano-channel. (a) Illitic channel. (b) Graphene channel. Error bars of MD simulations are smaller than the symbol sizes.

To propose a correction factor model for industrial applications, the dual-region model is extended to the shale matrix permeability correction model, and a theoretical comparison with the BKT model is performed. The correction factor function based on the dual-region model is:

where r* is the radius ratio of the inner region and the entire pore and the width of the adsorbed region is set as γσ. The presented intrinsic permeability values are 10 nD, 100 nD, and 1,000 nD (1 nD = 10−9D = 9.869 × 10−22m2), and the corresponding representative pore radius is in the empirical relationship52 of:

where k∞ is the intrinsic permeability and φ is the porosity. With the adsorbed region’s effect taken into consideration, the Langmuir-type adsorption isotherm is also applied, and the effect of Langmuir adsorption volume VL is compared. TMAC is assigned as 1 for the cases. Other parameters are set to be equal with those of the BKT model, and the parameters are in a reasonable range for actual shale gas reservoirs47 (see Supplementary Table S3 and Supplementary Methods).

Figure 8a shows that the correction factor is larger than that predicted by the BKT model in the high-Kn regime (or transition regime) because of the extra contribution in the adsorbed region. The underestimation by the BKT model increases as the k∞ decreases, corresponding to the decrease of R, due to the increase of the adsorbed region’s volume proportion. On the other hand, a larger VL results in a larger correction factor with the same k∞, and the difference increases with the decrease of k∞ for the same reason. In general, the BKT model underestimates the apparent permeability under relatively low-pressure conditions in ultra-tight reservoirs by neglecting the effect of adsorption.

In contrast to the high-Kn regime, the BKT model, as shown in Fig. 8b, will overestimate the apparent permeability in the low-Kn regime (i.e., the slip or viscous regime). Further, the BKT model’s deviation also increases with the decrease of k∞ (Fig. 8b), which corresponds to the effects shown in Fig. 5c,d,g,h. However, VL does not have much influence in either the 100 nD cases or the 1,000 nD cases because the velocity in the adsorbed region is very small. Nevertheless, a small difference of the correction factor is presented in the 10 nD cases because the adsorbed region maintains a large enough proportion to show the minute differences in the results depending on different VL values. On the whole, adsorption has an obstructive effect on the transport under relatively high-pressure conditions, the magnitude of which is related to the volume proportion of the adsorbed region accordingly. Furthermore, the correction factor from the BKT model will never be less than 1 according to equation (8), because the BKT model only considers the enhancement of slippage and the rarefaction effect, which is not the case under the high-pressure and narrow-pore conditions. Therefore, the correction factor model based on the dual-region model establishes a more accurate correction method by accounting for adsorption, which can be utilized as an up-scaled model in macroscopic modeling and simulation.

Figure 8. Permeability correction factor for shale matrix. (a) The correction factor vs. the Knudsen number. The curve of the BKT model almost coincides with that of the 1,000 nD case. (b) Zoom in on the plot for the

low-Knudsen-number region with log-axis. Except for in the 10 nD case, the curves of the Langmuir adsorption volumes (1 scf ton−1=2.83×10−5m3kg−1) are almost coincident.

In the present study, we have demonstrated the complex mechanisms of gas transport under the influence of adsorption in nano-channels composed of different materials. As a novel alternative to the continuum method, the perturbation model accurately describes the microscopic phenomenon of velocity oscillation by taking into account the density distribution across the channel, which is the result of the combination of gas–solid and gas– gas interactions. Additionally, this model provides a link between the microscopic method and the macroscopic method, which can introduce the results from the non-continuum method to the continuum method, as well as constituting a potential method for large-scale high-precision simulations based on the continuum description. In reference to the mechanisms of velocity oscillation as expressed by the perturbation model, the dual-region model offers a simplified approach to larger-scale simulation in the context of contemporary computing capability. As a novel method to handle the adsorption effect, the dual-region model also offers a correction factor model for the quantitative prediction of total mass flow with performance that is superior to that of classical models. In particular, the shale apparent permeability correction case presents an up-scaled model that have potential to be employed in macroscopic modeling and simulation for industrial application. Nevertheless, this study is still limited to the gas transport of simple fluid, while generalized model involved complex fluid remains to be developed. For the up-scaled model, the systematical analysis of uncertainty, which usually comes from the heterogeneity, tortuosity, and determination of characteristic radius, is also worthy to be investigated in the future work for the industrial application. Overall, the present study offers an in-depth theoretical investigation of the mechanisms and the impact of adsorption on gas transport in nanopores from phenomena to models and their applications.

1. Cueto-Felgueroso, L. & Juanes, R. Forecasting long-term gas production from shale. Proc. Natl Acad. Sci. USA 110, 19660–19661 (2013).

2. Clarkson, C. R. et al. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel 103, 606–616 (2013).

3. Aguilera, R. Flow units: From conventional to tight-gas to shale-gas to tight-oil to shale-oil reservoirs. SPE Reserv. Eval. Eng. 17, 190–208 (2014).

4. Skoulidas, A. I., Ackerman, D. M., Johnson, J. K. & Sholl, D. S. Rapid transport of gases in carbon nanotubes. Phys. Rev. Lett. 89, 185901 (2002).

5. Sholl, D. S. & Johnson, J. K. Making high-flux membranes with carbon nanotubes. Science 312, 1003–1004 (2006).

6. Holt, J. K. et al. Fast mass transport through sub-2-nanometer carbon nanotubes. Science 312, 1034–1037 (2006).

7. Bhatia, S. K. & Nicholson, D. Hydrodynamic origin of diffusion in nanopores. Phys. Rev. Lett. 90, 016105 (2003).

8. Jepps, O. G., Bhatia, S. K. & Searles, D. J. Wall mediated transport in confined spaces: exact theory for low density. Phys. Rev. Lett. 91, 126102 (2003).

9. Karniadakis, G., Beskok, A. & Aluru, N. Microflows and Nanoflows: Fundamentals and Simulation. 365–387 (Springer Science & Business Media, 2006).

10. Ruthven, D. M. Principles of Adsorption and Adsorption Processes. 29–61 (Wiley, 1984).

11. Adamson, A. W. & Gast, A. P. Physical Chemistry of Surfaces. 599–676 (Wiley, 1997).

12. Roque-Malherbe, R. M. A. Adsorption and Diffusion in Nanoporous Materials. 39–56 (CRC Press, 2007).

13. Maxwell, J. C. On stresses in rarified gases arising from inequalities of temperature. Philos. Trans. R. Soc. London, 231–256 (1879).

14. Knudsen, M. Die Gesetze der Molekularströmung und der inneren Reibungsströmung der Gase durch Röhren. Annalen der Physik 333, 75–130 (1909).

15. Pollard, W. & Present, R. D. On gaseous self-diffusion in long capillary tubes. Phys. Rev. 73, 762 (1948).

16. Klinkenberg, L. In Drilling and Production Practice. 200–213 (American Petroleum Institute, New York, 1941).

17. Darken, L. S. Diffusion, mobility and their interrelation through free energy in binary metallic systems. Trans. Aime 175, 184–201 (1948).

18. Evans III, R., Watson, G. & Mason, E. Gaseous diffusion in porous media at uniform pressure. J. Chem. Phys. 35, 2076–2083 (1961).

19. Kapoor, A. & Yang, R. Contribution of concentration-dependent surface diffusion to rate of adsorption. Chem. Eng. Sci. 46, 1995–2002 (1991).

20. Hu, X., Do, D. & Rao, G. Experimental concentration dependence of surface diffusivity of hydrocarbons in activated carbon. Chem. Eng. Sci. 49, 2145–2152 (1994).

21. Medveď, I. & Černý, R. Surface diffusion in porous media: A critical review. Microporous Mesoporous Mater. 142, 405–422 (2011).

22. Langmuir, I. The adsorption of gases on plane surface of glass, mica and platinum. J. Am. Chem. Soc. 40, 1361–1403 (1918).

23. Brunauer, S., Emmett, P. H. & Teller, E. Adsorption of gases in multimolecular layers. J. Am. Chem. Soc. 60, 309–319 (1938).

24. Beskok, A. & Karniadakis, G. E. Report: a model for flows in channels, pipes, and ducts at micro and nano scales. Microscale Thermophys. Eng. 3, 43–77 (1999).

25. Lilley, C. R. & Sader, J. E. Velocity profile in the Knudsen layer according to the Boltzmann equation. Proc. R. Soc. London, Ser. A

464, 2015–2035 (2008).

26. Schaaf, S. A. & Chambré, P. L. Flow of Rarefied Gases. 1–55 (Princenton University Press, 1961).

27. Beskok, A., Karniadakis, G. E. & Trimmer, W. Rarefaction and compressibility effects in gas microflows. J. Fluids Eng. 118, 448–456 (1996).

28. Zhang, H., Zhang, Z., Zheng, Y. & Ye, H. Corrected second-order slip boundary condition for fluid flows in nanochannels. Phys. Rev.

E 81, 066303 (2010).

29. Lockerby, D. A., Reese, J. M. & Gallis, M. A. The usefulness of higher-order constitutive relations for describing the Knudsen layer.

Phys. Fluids 17, 100609 (2005).

30. Guo, Z., Shi, B. & Zheng, C. G. An extended Navier-Stokes formulation for gas flows in the Knudsen layer near a wall. Europhys. Lett. 80, 24001 (2007).

31. Lockerby, D. A. & Reese, J. M. On the modelling of isothermal gas flows at the microscale. J. Fluid Mech. 604, 235–261 (2008).

32. Myong, R. A generalized hydrodynamic computational model for rarefied and microscale diatomic gas flows. J. Comput. Phys. 195, 655–676 (2004).

33. Javadpour, F. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Pet. Technol. 48, 16–21 (2009).

34. Darabi, H., Ettehad, A., Javadpour, F. & Sepehrnoori, K. Gas flow in ultra-tight shale strata. J. Fluid Mech. 710, 641–658 (2012).

35. Akkutlu, I. Y. & Fathi, E. Multiscale gas transport in shales with local Kerogen heterogeneities. SPE J. 17, 1,002–001,011 (2012).

36. Fathi, E. & Akkutlu, I. Y. Lattice Boltzmann method for simulation of shale gas transport in kerogen. SPE J. 18, 27–37 (2012).

37. Cracknell, R. F., Nicholson, D. & Quirke, N. Direct molecular dynamics simulation of flow down a chemical potential gradient in a slit-shaped micropore. Phys. Rev. Lett. 74, 2463 (1995).

38. Li, J., Liao, D. & Yip, S. Coupling continuum to molecular-dynamics simulation: Reflecting particle method and the field estimator. Phys. Rev. E 57, 7259 (1998).

39. Arya, G., Chang, H.-C. & Maginn, E. J. A critical comparison of equilibrium, non-equilibrium and boundary-driven molecular dynamics techniques for studying transport in microporous materials. J. Chem. Phys. 115, 8112–8124 (2001).

40. Sokhan, V., Nicholson, D. & Quirke, N. Fluid flow in nanopores: an examination of hydrodynamic boundary conditions. J. Chem. Phys. 115, 3878–3887 (2001).

41. Sokhan, V. P., Nicholson, D. & Quirke, N. Fluid flow in nanopores: Accurate boundary conditions for carbon nanotubes. J. Chem. Phys. 117, 8531–8539 (2002).

42. Bhatia, S. K. & Nicholson, D. Friction based modeling of multicomponent transport at the nanoscale. J. Chem. Phys. 129, 164709 (2008).

43. Falk, K., Coasne, B., Pellenq, R., Ulm, F.-J. & Bocquet, L. Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nat. Commun. 6, 6949 (2015).

44. Drits, V. A., Zviagina, B. B., McCarty, D. K. & Salyn, A. L. Factors responsible for crystal-chemical variations in the solid solutions from illite to aluminoceladonite and from glauconite to celadonite. Am. Mineral. 95, 348–361 (2010).

45. Wang, L. et al. Graphene oxide as an ideal substrate for hydrogen storage. Acs Nano 3, 2995–3000 (2009).

46. Ning, G. et al. High capacity gas storage in corrugated porous graphene with a specific surface area-lossless tightly stacking manner. Chem. Commun. 48, 6815–6817 (2012).

47. Yang, T., Li, X. & Zhang, D. Quantitative dynamic analysis of gas desorption contribution to production in shale gas reservoirs. Journal of Unconventional Oil and Gas Resources 9, 18–30 (2015).

48. Bhatia, S. K. & Nicholson, D. Molecular transport in nanopores. J. Chem. Phys. 119, 1719–1730 (2003).

49. Blundell, S. & Blundell, K. M. Concepts in Thermal Physics. 76–83 (Oxford University Press, 2010).

50. Wilcox, D. C. Turbulence Modeling for CFD. 33–59 (DCW industries La Canada, CA, 1998).

51. Ziarani, A. S. & Aguilera, R. Knudsen’s permeability correction for tight porous media. Transport in Porous Media 91, 239–260 (2012).

52. Aguilera, R. Incorporating capillary pressure, pore throat aperture radii, height above free-water table, and Winland r35 values on Pickett plots. AAPG Bulletin 86, 605–624 (2002).

This work is partially funded by the National Natural Science Foundation of China (Grant No. U1262204, U1663208 and 51520105005) and the National Key Technolog y R&D Program of China (Grant No. 2012BAC24B02).

T.W. and D.Z. designed the study, analyzed the data, and wrote the paper together. T.W. performed the MD simulations and derived the models.

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Wu, T. and Zhang, D. Impact of Adsorption on Gas Transport in Nanopores. Sci. Rep. 6, 23629; doi: 10.1038/srep23629 (2016).

Correspondence and requests for materials should be addressed to Dongxiao Zhang (email: dxz@pku.edu.cn)

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/

AllAboutShale.com