Pressure-driven supercritical CO2 transport through a silica nanochannel

A thorough understanding of supercritical CO2 (scCO2) transport through nanochannels is of prime significance for the effective exploitation of shale resources and the mitigation of greenhouse gas emission. In this work, we employed the non-equilibrium molecular dynamics simulations method to investigate the pressure-driven scCO2 transport behavior through silica nanochannels with different external forces and pore sizes. The simulations reveal that the capability of scCO2 diffusion enhances both in the bulk region and the surface adsorbed layer with the increasing of pressure gradient or nanochannel size, in addition, the slip length increases nonlinearly with the external acceleration or nanochannel width increases and finally reaches a maximum value. The negative slippage occurs at lower pressure gradient or within the narrower nanochannel. Overall, it is the combined effect of strong adsorption, surface diffusion and slippage that causes the nonlinear relation between flow rate and pressure gradient or nanochannel size. The present work would provide theoretical guidance for the scCO2 enhanced shale oil/gas recovery, CO2 storage, and mass transport in nanoporous materials.


Introduction
Recently, the behaviors of supercritical CO 2 (scCO 2 ) in narrow pores and conned spaces have attracted extensive attention in the exploration and development of hydrocarbons stored in shale, 1,2 such as scCO 2 fracturing 3,4 and injecting. 5,6 On the other hand, CO 2 injection into tight shales is one of the promising solutions for mitigating carbon emissions 7,8 and global warming. 9,10 The nanopores are widely present in shale reservoirs. [11][12][13][14] The basis for scCO 2 geosequestration 15, 16 and enhanced hydrocarbon recovery depends on its unique transport properties induced by the connement of shale nanopores. The evaluation of scCO 2 enhancing hydrocarbon recovery, enhancement of fracturing efficiency, and estimation of storage capacity require understanding of microscopic behaviors of scCO 2 through shale nanopores, such as adsorption, diffusion and ow.
Conned uids in nanopores demonstrate microstructural, dynamical and thermophysical behavior that differ signicantly from their bulk counterparts due to the molecular asymmetry between uid-uid and solid-uid interactions. 17 This asymmetry leads to inhomogeneous local uid density distributions, where the superimposition of the attractive potentials from two walls induces conned environments. 18 Extensive molecular simulation studies were reported on CO 2 sorption in the graphite slit-pores [18][19][20][21][22][23][24] usually treated as the organic matrix of shale reservoirs. In addition, several molecular simulations were addressed on the adsorption structure of CO 2 conned in clay-like 25,26 and silica 17,27 slit pores. The results revealed that the nanopore size dramatically inuences the adsorption structure of CO 2 , which would has potentially signicant effects on the transport properties of CO 2 inside the nanopore. 28 The diffusion of conned CO 2 obviously weakens as compared with that in bulk uid 25,27 and CO 2 diffusion coefficient reduces with the decrease in the pore size. 25 The diffusion and migration of CO 2 in silica nanopore occur predominantly along the surface due to hydrogen bonds formed between CO 2 and the surface -OH group. 27,29 Besides, the widely developed nanopores and ultralow permeability in shale reservoirs raise concerns on the applicability of Darcy's law to account for mass transport in shale. However, it breaks down for ow through very narrow pores because the continuum ow hypothesis is invalid. 23,28 Efforts have been made to moderate for the breakdown of Darcy's law via gas slippage theory, 23 for example, the Klinkenberg effect. 30,31 Such empirical corrections can not capture the complex adsorption and transport behavior of the uid in ultra-conning porous materials. 32 Falk et al. demonstrated that the non-Darcy behavior results from strong adsorption in kerogen and the breakdown of hydrodynamics at the nanoscale. 32 Monteiro et al. suggested a hydrodynamic model of gas ow in nanoporous media by introducing a pressure gradientdependent permeability of kerogen, which obeys power function relationship. 33 Wang et al. reported the pressure-driven ow behavior of scCO 2 conned in shale organic nanopores. 24 They observed that the streaming velocity prole tends to a plug ow and the magnitude far exceeds the prediction of the no-slip Poiseuille equation. Firouzi and Wilcox investigated the gas slippage and Klinkenberg effects of CO 2 conned in carbon slit pores. 23 They reported a plug ow occurs at pore size less than 2 nm, a parabolic velocity prole at pore size larger than 10 nm, and the dramatic underestimation of the CO 2 permeability using the bulk phase viscosity. In summary, the transport of scCO 2 in inorganic nanochannel remains unclear, and there is a strong need for a reliable theoretical framework of scCO 2 transport inside nanoporous matrix due to the constraint of the lowest permeability in the uid path to the overall permeability of the shale.
Herein, non-equilibrium molecular dynamics (NEMD) simulations were conducted to investigate the pressure-driven scCO 2 transport through the silica nanochannel. The effects of ow driving pressure gradient and nanochannel size on scCO 2 transport were determined by exploring the adsorption, diffusion and ow of scCO 2 molecules. The density proles and residence autocorrelation functions were employed to describe adsorption. The velocity proles, the slippage effect, and the ow rate as the function of pressure gradients or nanochannel sizes were applied to characterize the ow behavior.

Models and methodology
The nanochannel with the size of 71.60 Â 50.12 Â HÅ 3 was formed by two 10.5Å cristobalite plates. H denoted the interlayer space of two plates providing the nanochannel for scCO 2 transport. The silica nanochannels with hydroxylated surfaces 34 were xed in all simulations. Periodic boundary conditions were applied in three dimensions. To obtain saturated adsorption of scCO 2 in nanochannels under 373.15 K and 35 MPa, equilibrium molecular dynamics (EMD) simulations were carried out for 2 ns to get the isothermal adsorption curve 35 (ESI Fig. S1 †). The initial model was shown in Fig. 1(a).
As applying an acceleration or a force on all atoms directly 24,36-38 was improper to simulate desorption from surfaces, the whole nanochannel was divided into three segments along X dimension, named G-region, Buffer and S-region, respectively, as shown in Fig. 1(b). A driving pressure gradient was exerted on the system by applying an acceleration along X direction on scCO 2 molecules in G-region in order to generate the Poiseuille ow through the nanochannel. A bigger acceleration corresponds to a greater pressure gradient. The effect of pressure gradient on transport behavior was examined by altering the acceleration ranged from 0.001 to 0.05 kcal (Å À1 g À1 ) as H is 28.50Å. The effect of nanochannel sizes on transport behavior was studied by altering H from 7.25 to 50.00Å as the acceleration is 0.01 kcal (Å À1 g À1 ). All calculations were executed in S-region. The NEMD simulations were performed to simulate mass transport for 6 ns which is enough to obtain a steady scCO 2 ow.
Both the EMD and NEMD simulations were performed using large-scale atomic/molecular massively parallel simulator (LAMMPS) 39 soware in the NVT ensemble with a time step of 1 fs. Visualization of the dynamics process was carried out using the program VMD. 40 The Nosé-Hoover 41,42 thermostat was employed to control the temperature. The silica nanochannel was modeled with the CLAYFF 43 force eld. CO 2 molecules were described by the three-site EPM2 (ref. 44) model which was successfully demonstrated to reproduce critical point of CO 2 in experiment. All the VDW interactions were described by the 12-6 Lennard-Jones potential. The cutoff of non-bonded interaction was set to be 9.5Å. Long-range electrostatic interaction was calculated by PPPM 45 summation method with the accuracy of 0.0001e. The interaction between different types of atoms was determined by the Lorentz-Berthelot rule:3 ij ¼ ffiffiffiffiffiffi ffi 3 i 3 j p and s ij ¼ (s i + s j )/2. All parameters were given in the ESI Table S1. †

Results and discussions
3.1 Effect of pressure gradient 3.1.1 Effect of pressure gradient on adsorption of scCO 2 . In Fig. 2(a), we presented the density proles of scCO 2 under different pressure gradients to explore the inuence of pressure gradient on adsorption. As shown, the density proles along Z direction are symmetrical about the central axis of the nanochannel. There are two primary peaks and two secondary peaks appear near the surfaces, indicating that the scCO 2 form four structuring adsorption layers. This shows that the scCO 2 molecules pack more closely near the surface than those far away from the surfaces. The monolayer peak positions are almost independent of pressure gradient, while the second layer peaks shi slightly as the gradient pressure increases and almost disappears at the highest gradient pressure. The density in the center of the channel increases slightly under relatively low pressure gradients while increases obviously at higher pressure gradients. The results indicate that a high pressure gradient dramatically affects the distribution of scCO 2 inside the nanochannel. Under the high pressure gradients, some of adsorbed scCO 2 desorb from the surfaces and enter into the channel central areas. Such behavior could improve the effective permeability 46 of scCO 2 in channel, which contributes to the scCO 2 ow. However, the shape of the density proles remains almost unchanged under the lower pressure gradients indicating a little inuence on the distribution of scCO 2 , which is in agreement with the results obtained by Kasiteropoulou. 47 In order to quantify the residence time of scCO 2 molecules remaining near the nanochannel surface, residence autocorrelation function, C R (t), was calculated by: 48-50 where O w (t) describes whether one scCO 2 molecule is in the rst adsorption layer at time t. O w (t) equals one indicating a scCO 2 molecule within the rst adsorption layer. O w (t) equals zero meaning that a scCO 2 molecule leaves the rst adsorption layer. The angular bracket denotes the ensemble average. Carbon atom in scCO 2 is considered to identify the position of one scCO 2 molecule. C R (t) decays from 1 to 0. The faster decay of C R (t) implies the shorter time for scCO 2 residing in the rst adsorption layer. The C R (t) for scCO 2 at the different pressure gradients is shown in Fig. 2(b). C R (t) decreases faster as pressure gradient increases, implying a shorter residence time of scCO 2 in the adsorption layer at bigger pressure gradients. The high pressure gradients could facilitate the uid mobility and enhance the occurrence of molecular collisions, and then driving the adsorbed scCO 2 molecules away from the surface. The C R (t) plateaus occurs aer about 120 ps, indicating that some of the adsorbed scCO 2 molecules are not desorbed within the simulation time because of the electrostatic interaction and the hydrogen bond between scCO 2 and the surface. 51 These scCO 2 molecules that can not be desorbed from the surfaces will reduce the scCO 2 ow in shale reservoirs but benet the CO 2 sequestration.
3.1.2 Effect of pressure gradient on diffusion coefficient of scCO 2 . The diffusion coefficient is a signicant parameter revealing the mechanism of mass transfer. The self-diffusion coefficient can be obtained by Einstein's relation where r j is the position vector of the j scCO 2 molecule, d is the dimension of the system, N is the total scCO 2 molecules and the angular bracket represents the ensemble average. In this work, D was calculated for motions in Y and Z dimensions, and hence d ¼ 2. Fig. 3 illustrates the diffusion coefficients of scCO 2 at different pressure gradients and areas in the nanochannel. As  This journal is © The Royal Society of Chemistry 2018 pressure gradient increases, the diffusion coefficient of scCO 2 in the adsorption layers increases while it decreases rstly and then rises for bulk scCO 2 . The variation of diffusion coefficients of scCO 2 under lower pressure gradients was mainly affected by the variation in the density distribution. As shown in Fig. 2(a), the density near the surface reduces with the increasing pressure gradient as the acceleration is less than 0.007 kcal (Å À1 g À1 ). This indicates that the partial pressure in adsorbed layers decreases, leading to the increase of scCO 2 diffusion coefficient. The increasing density of bulk scCO 2 indicates the enhancement of the partial pressure in the central area of nanochannel, resulting in the decrease of scCO 2 diffusion coefficient. In this case, the density distribution plays the dominant role in the diffusibility of scCO 2 within the nanochannel. As the acceleration is greater than 0.007 kcal (Å À1 g À1 ), the largely undermined adsorption strength with increasing pressure gradient indicates the increase of the scCO 2 velocity within the nanochannel. The raising velocity may lead to the decrease of scCO 2 pressure, which induces the increase of the diffusion coefficient. Herein, the increasing velocity plays the pivotal role in scCO 2 diffusibility inside the nanochannel. As indicated from the Stokes-Einstein relation, the increase in the diffusion coefficient implies the decrease in the viscosity. Therefore, we conclude that such changes could affect the ow of scCO 2 within the nanochannel.
3.1.3 Effect of pressure gradient on ow behavior of scCO 2 . In Fig. 4(a), we depicted the velocity proles at different pressure gradients in order to characterize the ow behavior of scCO 2 within the nanochannel. The velocity proles were obtained by tting data using the quadratic function: where a, b, and c are the tting parameters describing the ow process, a characterizes the curvature. The lager the absolute value of a is, the steeper the curve and the greater the velocity gradient. c represents the maximum velocity. The tting parameters were provided in the ESI Table S4. † In Fig. 4(a), the velocity proles with the parabolic shape 24,36,38 are symmetrical about the central axis of the nanochannel. The larger pressure gradients induce the greater axial velocities and sharper parabolic velocity curves. Moreover, we notice that there are non-zero velocities near the surfaces, which leads to the slippage phenomenon. The slip velocity increases as the pressure gradient increases. Here, we employed the slip length L s to describe the slippage phenomenon, as dened by 52 where v(z surf ) is the slip velocity, v is the velocity prole in the ow direction, and (dv/dz) z surf is the velocity gradient with respect to z on the surface. Slip length is the distance extrapolated into the surface where the velocity vanishes as supposed by no-slip boundary condition. 53 Fig. 4(b) demonstrates the slip lengths for scCO 2 ow in silica nanochannel under different pressure gradients. The inset shows the derivatives of the velocity with respect to z on the surface and the slip velocities. Clearly, a negative slip length 54-57 occurs as the acceleration is less than 0.004 kcal (Å À1 g À1 ). The main reason can be explained as follows. Under small pressure gradients, the shear force, as reected by velocity gradient in the inset in Fig. 4(b), is too small to overcome the attractive interactions between scCO 2 and the surface. This leads to the zero velocity for scCO 2 molecules near the surfaces, indicating the presence of scCO 2 boundary layer whose thickness equals slip length. In this case, adsorbed scCO 2 are trapped on nanochannel surfaces, which is so-called sticking phenomenon, 58 conning the ow and diffusion of scCO 2 . In addition, the slip length shows a non-monotonous increase with pressure gradients. The increasing pressure gradient promotes scCO 2 ow and detaches scCO 2 molecules from the surfaces due to the enhancement of shear rate, which stimulates surface diffusion 59 of scCO 2 . Moreover, the rate of increase in slip velocity is greater than that of velocity gradient. Consequently, slip length increases rapidly with increasing pressure gradient. As pressure gradient further increases, the ratio of slip velocity to velocity gradient decreases gradually which cause the slip length increases slowly and reaches a maximum value nally. Aer that, the velocity of bulk scCO 2 increases much faster than slip velocity with the increase in pressure gradient due to the strong adsorbed potential of the surfaces. This leads to a large increase in velocity gradient greater than the rate of increase in slip velocity. Therefore, slip length will decrease gradually as pressure gradient increases. Fig. 5 depicts the scCO 2 ow rate as a function of pressure gradient, revealing a nonlinear relationship different from Darcy's law where permeability is dened as a material property of the rock. This indicates that Darcy's law dramatically fails to describe scCO 2 transport in nanochannel. Form the analysis of Fig. 2 and 3 mentioned above, we conclude that the permeability increases and the viscosity decreases as pressure gradient increases. Hence, the ratio of the permeability to the viscosity, i.e., mobility, increases with the increasing of pressure gradient, characterizing the seepage capacity of scCO 2 in nanochannel. It is observed that the mobility is linearly related to pressure gradient shown by the inset in Fig. 5. This corresponds to the hydrodynamic model of gas ow in nanoporous media where permeability is power function of pressure gradient. 33 3.2 Effect of nanochannel width 3.2.1 Effect of nanochannel width on adsorption of scCO 2 . The density proles of scCO 2 in different nanochannels were presented in Fig. 6(a), showing the formation of several structuring layers for scCO 2 molecules, which is consistent with previous reports. 18,24,37,60 For H ¼ 7.25Å, the two symmetrical density peaks indicate that scCO 2 form two primary adsorption layers due to partial superimposition of the attractive potentials from two surfaces. 37 Accordingly, the peaks of the density proles are the highest. As H ¼ 14.5Å, the appearing of two primary peaks and two secondary peaks indicates that there are four adsorption layers form within the nanochannel. This arises from the further separation of the attractive potentials from two surfaces. Accordingly, the primary peaks of the density proles decrease. For H ¼ 28.5Å, bulk scCO 2 occurs in the central area of the nanochannel apart from the four adsorption layers. The density peaks further decrease because the attractive potentials from two surfaces have been effectively separated into two single potential system. A similar phenomenon can be observed as H further increases from 35.5Å to 50Å. The results show that nanochannel size has a signicant inuence on the adsorption mechanism of scCO 2 especially in very narrow channels. For narrower nanochannels, the superposition of attractive potentials from two surfaces results in strong adsorption of scCO 2 on the surfaces. This can reduce the permeability of scCO 2 in shale porous medium, which contributes to CO 2 storage rather than its ow in shale reservoirs. Fig. 6(b) shows the residence autocorrelation functions, C R (t), as a function of the residence time of scCO 2 molecules remaining near the surface. As H increases, C R (t) decays fast, implying the shorter time that scCO 2 remains in the absorbed layer. As H increases, the overlapping attractive potentials from two surfaces gradually separate from each other and nally form two single potential systems near the surfaces. This weakens the interaction between scCO 2 and the surfaces resulting in the shorter residence time of scCO 2 on the surface, and also implies a more frequent exchange of scCO 2 appeared between the adsorption layer and bulk phase in the center of nanochannel inside the wider nanochannel. For H ¼ 7.25Å, C R (t) always equals to one indicating that all scCO 2 molecules reside in the adsorption layer. As H further increases from 14.5 A to 50.0Å, C R (t) plateaus appears aer about 120 ps indicating that a part of scCO 2 molecules always reside near the surface within the simulation time. Value of C R (t) plateaus decrease as H increases, also showing the stronger adsorption for narrower nanochannels. As H ranges from 35.5Å to 50.0Å, C R (t) curves are basically coincident, suggesting that scCO 2 adsorption is independent of H larger than 35.5Å. This is mainly due to that attractive potentials from two surfaces have been effectively divided into two single potential systems near the channel surface. The results also demonstrate that the narrower nanochannel storages CO 2 more securely while the wider nanochannel favors scCO 2 transport.
3.2.2 Effect of nanochannel width on diffusion coefficient of scCO 2 . Fig. 7 demonstrated the dependence of the diffusion coefficient on H. The diffusion coefficient increase monotonically with H, indicating that scCO 2 molecules migrate more easily in wider channels. Compared to the less constraint of bulk phase with a larger diffusion coefficient, the adsorption phase experience a stronger constraint from the surfaces which leads to a smaller diffusion coefficient. The order of diffusibility for scCO 2 at different areas in the nanochannel is bulk phase > second adsorption layer > rst adsorption layer, as described in Fig. 7. The diffusion coefficient of scCO 2 molecules in the rst adsorption layer is nonzero, which may indicate that the surface diffusion occurs by hopping of scCO 2 molecules from one site to another. 61 Therefore, surface diffusion is one factor that scCO 2 molecules move more quickly in the nanochannel. This mainly results from two aspects. One is the decrease of scCO 2 adsorption density (Fig. 6(a)), leading to the viscosity reduction. The other is the increase of scCO 2 velocity gradient (the inset in Fig. 8(b)), indicating the enhancement of shear effect and then results in the shear thinning, which reduces scCO 2 viscosity.
3.2.3 Effect of nanochannel width on ow behavior of scCO 2 . To understand the effect of nanochannel width on scCO 2 ow behavior, the velocity proles and the slip lengths were plotted in Fig. 8(a) and (b), respectively. The related tting  Table S5. † In Fig. 8(a), the scCO 2 velocity proles exhibit a parabolic shape and are symmetrical about the central axis of the nanochannel. Moreover, as H increases, the velocity curve becomes sharper and the axial velocity of scCO 2 increased. This is caused by the intensive constraint on the scCO 2 molecules near the surface in the smaller nanochannel, as shown in Fig. 6. However, the bulk scCO 2 in the central area of the wider nanochannel are rare constrained.
In Fig. 8(b), the slip length exhibits non-monotonous increase with a maximum value as H increases. For 7.25Å < H < 14.5Å, the slip velocity increases fast while the velocity gradient increment is small as H expands, as indicated by the inset in Fig. 8(b). Then the ratio of the slip velocity to the velocity gradient, i.e., slip length, increases rapidly. For 14.5Å < H < 28.5Å, both the slip velocity and the velocity gradient increase rapidly. Ratio of slip velocity to velocity gradient enhances gradually and reaches the maximum value as H ¼ 28.5Å. Therefore, the slip length increases slowly and its maximum value occurs. For 35.5Å < H < 50.0Å, the increase of the slip velocity becomes slower while the velocity gradient increase rapidly. The ratio of them weakens gradually. Consequently, the slip length decreases slowly. Fig. 8(b) also shows a negative slip length 54-57 occurring at H less than about 10.0Å. In this case, the overlapping between attractive potentials from the surfaces is so strong that a part of scCO 2 molecules stick on the surfaces. The external pressure gradient cannot completely overcome the attractive interactions between scCO 2 and the surface. This leads to the zero velocity near the surfaces, indicated the presence of scCO 2 in boundary layer.     9 describes the dependence of ow rate on H. As H is less than 30.0Å, ow rate increases slowly due to the strong constraint on scCO 2 molecules from the superimposition of attractive potentials of two surfaces. The attractive potential from the surface has a signicant effect on scCO 2 transport through nanochannel. As H is larger than 30.0Å, ow rate increases rapidly because of the less constraint on scCO 2 molecules where attractive potentials from two surfaces separate completely. Here, the attractive potential from the surface has little effect on the transport of scCO 2 in nanochannel. Therefore, nanochannels with H less than 30.0Å contribute to CO 2 storage in shale. To quantitatively describe the dependence of the ow rate on H, the relationship between them was obtained by tting data applying the power function polynomial The tting parameters were shown by the inset in Fig. 9. Compared to the continuum ow with slip effect, 52 the fourth power of H appears in the tted formula. This could be mainly due to the fact that both the viscosity and the slip length of scCO 2 in nanochannel relate to H, as indicated in Fig. 7 and 8(b).

Conclusions
NEMD simulations were employed to study the scCO 2 transport behavior within the silica nanochannel in shale. The inuences of ow driving pressure gradient and nanochannel size were systematically investigated. These effects were understood by exploring density distribution of scCO 2 molecules, residence time of scCO 2 in the rst adsorption layer, the diffusion coef-cient, the velocity prole, the slippage effect, and the ow rate.
Our simulations reveal that the structuring adsorption layers of scCO 2 molecules occur within the silica nanochannel yet under the driving of pressure gradient. As either pressure gradient or nanochannel width increases, the peaks of scCO 2 density proles decrease and the density of the bulk phase increases. The higher pressure gradients result in the larger shear force, which leads to desorbing scCO 2 molecules near the nanochannel surfaces. The wider nanochannel leads to the weaker overlapping of the attractive potentials from the two surfaces, attenuating the adsorption of scCO 2 molecules near the surfaces. The higher pressure gradients and the wider nanochannel can promote scCO 2 diffusibility for both bulk region and the adsorption layers near the surfaces. The improvement of scCO 2 diffusion ability may be mainly determined by the decrease of partial pressure near the surfaces, the increase of velocity of the central area in the nanochannel, and the surface diffusion of scCO 2 hopping from one site to another. The increased diffusion ability also indicates the reduction in scCO 2 viscosity. The negative slippage occurs at lower pressure gradients or within the narrower nanochannel is mainly due to the stronger adsorption of scCO 2 molecules near the surfaces. The slip length exhibits non-monotonic increase and existence of a maximum value with the increase in the pressure gradient or the nanochannel size. The ow rate of scCO 2 through the nanochannel shows a nonlinear characteristic with pressure gradient or nanochannel size. It is also found that the mobility is proportional to the pressure gradient. This implies the breakdown of Darcy's law at nanoscale because of the strong adsorption effect, surface diffusion and slippage effect.

Associated content
Additional simulation results to determine the initial state of scCO 2 in channel, parameters of force eld and a part of tting parameters are described in the text.

Conflicts of interest
There are no conicts to declare. Fig. 9 The dependence of flow rate on channel width for scCO 2 molecular.
This journal is © The Royal Society of Chemistry 2018