Effect of Surface Type on the Flow Characteristics in Shale Nanopores

Key Laboratory of Unconventional Oil & Gas Development, Ministry of Education, Qingdao 266580, China School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, Shandong 266580, China Petroleum Engineering Technology Research Institute of Shengli Oilfield, SINOPEC, Dongying 257067, China Postdoctoral Scientific Research Working Station of Shengli Oilfield, SINOPEC, Dongying 257067, China Shengli Oilfield Exploration and Development Research Institute, SINOPEC, Dongying 257015, China Research Institute of Exploration and Development, Tarim Oilfield Company, PetroChina, Korla 84100, China Exploration and Development Research Institute of North China Oilfield Company, PetroChina, Renqiu 062550, China


Introduction
With the depletion of conventional energy resources and the expanding energy demands, the unconventional energy resources such as the shale oil and gas have attracted more and more attention [1][2][3][4][5][6][7]. Unlike the conventional gas reservoir, numerous pores in shale are in nanometer scale [7][8][9][10][11]. According to the value of pore diameter, d, the pores could be divided as micropore (d < 2 nm), mesopore (2 nm ≤ d ≤ 50 nm), and macropore (d > 50 nm) suggested by the International Union of Pure and Applied Chemistry (IUPAC) [1,2,12]. In one aspect, such amounts of nanopores provide large internal surface areas and result in great adsorption of shale gas; the adsorbed gas occupies a significant percentage of the gas-in-place in shale [13][14][15]. The bulk density and the viscosity cannot represent the physical properties actually since the distribution of the gas confined in nanopore may be inhomogeneous [16,17]. In another aspect, as the mean free path of methane molecules changes to be comparable to the pore size in such tiny pores [18][19][20], the traditional Darcy Law may be invalidated and make the prediction of shale gas flow be challenging [6,21,22].
Basically, as shown in Figure 1, the gas flows are generally classified into four different regimes referred to the key parameter, Knudsen (Kn) number, which is defined as the ratio of mean free path of gas molecule (λ) to the characteristic length of the pore (such as the pore width, W), Kn = λ/W [12,21,23,24]. The mean free path of the molecule λ could be derived by λ = R g T/ ffiffi ffi 2 p πσ 2 N A P, where R g is the gas constant, T and P mean the temperature and the pressure, respectively, σ is the diameter of the molecule, and N A is the Avogadro number [23]. Thus, for a given reservoir pressure and temperature condition and a certain type of fluid molecule, the Kn is only related to the pore size. When the pore size is enough large, Kn is sufficiently small and lower than 0.001; the fluid flow is assumed as continuum flow and could be described by the no-slip condition combined with the Navier-Stokes equation [25,26]. The Poiseuille flow under slit-like pores could be predicted as the following [22,23]: in which η is viscosity and ∂P/∂x means the pressure gradient in the flow direction, x-direction. When the pore size decreases and 0:001 < Kn < 0:1, the flow is regarded as in a slip flow regime, and the slip boundary condition is often applied with the Navier-Stokes equation. Such slip boundary is often characterized by the slip length L s , and the flux could be corrected by incorporating the slip length into the Navier-Stokes equation [21,27], At the same time, to measure the flow enhancement owing to the velocity slip on the pore surface, the enhancement factor ε is raised and defined as the ratio of the real flux to the value predicted by Equation (1), Then, the relationship between the enhancement factor ε and the slip length L s is given as [22,23] If the pore size is small enough and 0:1 < Kn < 10, the fluid flow regime is assumed to be transition flow. In such cases, the gas interaction with the pore wall as well as the gas slippage plays more important roles to the fluid flow. Generally, in shale pore systems under the reservoir conditions, the flow regimes are mainly in the range of slip flow and the transition flow [28].
To have a better understanding on the flow mechanisms and a much better prediction of such abnormal fluid flow in small nanopores, numerous research has been finished. Some researchers believe that the transport mechanisms of the shale gas are composed of the viscous flow, slip flow, Knudsen diffusion, and surface diffusion in the current [29][30][31][32][33], while some other scholars agree that the slip flow and Knudsen diffusion describe the same phenomenon, which is mainly the migration mode caused by the collision between the gas particles and the pore surface wall [24,34,35]. Thus, the calculation by coupling only the viscous flow, slip flow (or Knudsen diffusion), and surface diffusion is sufficiently accurate in the predict of the gas flow in shale [34]. Sheng et al. established one mathematical theoretical model in which considering these triple shale gas migration modes and the physical model are given in Figure 2 [24]. The flow region in nanopores could be divided into three types, adsorption layer, Knudsen layer, and bulk gas region. In the adsorption layer, the gas molecules are assumed as mainly driven by the chemical potential gradient. Whereas in the Knudsen layer, the collisions between the free gas molecules and the adsorption layer play the dominant role. In the bulk gas region, the free gas molecules are mainly driven by the pressure gradient and the traditional viscous flow forms. They also showed that the surface diffusion of the gas molecules in the adsorption layers plays a key role in the gas slippage, especially for the fluids in sub-5 nm nanopores [24]. Such previous models shed light on the gas flow mechanisms and provide various methods to predict the gas flux, while it is still confusing in the methods selection from so many models or methods [9,12,24,28,31,[36][37][38].
Furthermore, to have an accurate understanding of gas transport characteristics in nanoporous/microporous media, the molecular dynamic (MD) simulation is also widely applied to understand the gas flow behaviors, as summarized and depicted in one recent review by Yu et al. [39,40]. The nanochannels are generally used to mimic the fluids flow confined in shale nanoporous media in these MD simulations [22,41,42]. Since the shale organic matter is mainly composed of carbon atoms and the pore wall is hydrophobic; thus, the graphene sheets, carbon nanotubes, and arrayed carbons are abstracted as organic pore wall materials in numerous previous work [39]. Employing parallel carbon walls to model the complex organic pore structure, Chen 2 Geofluids et al. investigated the pressure-driven shale gas flow behavior; they found that the velocity profiles translate from plug-like to parabolic when the pore width increases from 2 nm to 10 nm [41]. The similar carbon structure was also used to study the slip behavior of shale gas by Nan et al. [22]. Recently, the slit-like pore confined in amorphous kerogen structure was constructed and used to reveal the transport behavior of geo-fluids. Unlike the obvious slippage flow confined in graphene nanochannel, the velocity profiles in such slit-like kerogen pore display no-slip parabolic shape [40]. The roughness in atomic-scale roughness plays a nonnegligible role in interrupting the molecules moving along the pore wall. Furthermore, Yu et al. also found that the roughness has great effect on the mass-transport of shale gas using a series of rough kerogen structures [43]. While in the MD simulation, it is not easy to control the roughness of the kerogen structure [43]. Actually, the pore wall in shale includes various minerals, such as inorganic minerals (i.e., quartz, feldspar, calcite, and clays), organic minerals (i.e., kerogen) [2,11,44]; it is still hang on doubt that whether the surface diffusion of gas molecules is great of importance in all these types of shale nanopores or not. To reveal the effect of the type of pore surface to the gas migration and provide the reference to the application and selection of the various theoretical models, we use the molecular dynamic (MD) simulation to study the gas flow behavior in organic and inorganic nanopores under the typical reservoir conditions. In addition, the flux contributions of each flow region in various pore sizes are also compared. Here, the pure methane molecules are used to represent the shale gas since the methane is the most one component of shale gas usually.
The paper is organized as follows. In Section 2, the molecular models, force fields, and the simulation procedure are introduced briefly. In Section 3, the static properties (i.e., density profiles) and the hydrodynamic performance are presented firstly; then, the effect of pore surface to the flow enhancement and flux contributions of each region are introduced in detail. At last, the key conclusions are summarized, and the potential instructions are discussed in Section 4.

Molecular Models and Simulation Methods
2.1. Organic and Inorganic Nanopore Models. In this work, to weed out the effect of random roughness of the wall structure to the gas transport, a full atomic smooth structure of three graphite layers is used to represent the organic nanopores in shale reservoir (as shown in Figure 3). Such a simplified method is widely used in the MD simulations which are focused on the research of geofluids adsorption and flow performance in shale organic pores [15,[21][22][23]45]. The dimension of the graphene layers in x-and y-direction are L x = 20:27 nm and L y = 10:07 nm, respectively, and such a box size is enough to make the transport behavior of the gas molecules reach a steady state in the NEMD simulation. The thickness of the three-layer carbon sheets, h, is 0.67 nm, which is in agreement with the previous works [21][22][23]. In addition, the distance in z-direction between the inner plane of the top carbon sheets and the bottom sheet are defined as the pore width, H. Here, the 5 nm, 10 nm, and 15 nm are selected since that the mesopore are abundant in shale according to the pore size distribution equipment [44,46].
Meanwhile, the fully hydroxylated quartz sheet is adopted to indicate the pore wall of inorganic nanopore [47]. As a brittle mineral, the quartz makes the shale formation easy to be fractured and obtain the commercial industry production [7]. While in the process of diagenetic evolution during the shale sedimentary, the quartz surfaces could capture the -OH groups and be fully hydroxylated. Thus, the (100) crystallographic orientation of a-quartz is prepared firstly, and then, the flexible -OH groups are added on such cleaved surface fully. The 3D snapshot of the inorganic quartz nanopore is presented in Figure 4.
In this work, the pressure and the temperature of the molecular dynamic simulation are selected as 350 K and 50 MPa, which is in the range of typical shale reservoir condition [22]. Then, referring to the bulk methane density at the reservoir condition from NIST database [48], ρ bulk = 0:2336 g/cm 3 , the number of methane molecules confined in the nanopore is evaluated firstly. While due to the interaction between the pore wall and the methane, the excess adsorption may make the bulk methane density observed by the

Geofluids
MD simulation, ρ MD bulk deviated from the standard value. Then, to make ρ MD bulk equals to ρ bulk , according to the measurement discrepancy between the ρ MD bulk in the MD simulation and standard ρ bulk , the corresponding number of the methane molecules could be calculated combined with the pore volume. And then, the renewable number of methane molecules is inserted in the nanopore and compare ρ MD bulk with ρ bulk . By this loop, the suitable number of the methane molecules could be obtained until ρ MD bulk = ρ bulk , as stated in the previous work [49]. The detailed information of each MD simulation system is illustrated in Table 1. Periodic boundary conditions are all applied in the x-, y-, and z-directions.

Force Fields.
The inorganic quartz surface is described by CLAYFF force field [50,51], and the potential of organic graphene is listed in the previous work [52]. At the same time, the TraPPE model is used to describe the methane [53]. All the interactions between methane and pore surface are represented by the Lennard-Jones (LJ) 12-6 potentials and Coulomb electrostatic interactions [53,54], where q i and q j are the partial charges of site i and j, ε 0 is the dielectric constant of vacuum, r ij indicates the separation distance between particle i and j, and ε ij and σ ij are the LJ energy and size parameters, which are deducted from the standard Lorentz-Berthelot mixing rules, in which σ i and σ j are the LJ size parameters of atom i and j, and ε i and ε j are the LJ energy parameters of atom i and j. The Coulomb interactions are computed using standard three-dimensional particle-mesh Ewald (PME) method. All nonbonded simulation parameters used in this work are listed in Table 2. In addition, the cut-off radius of the LJ interaction is set as 1.5 nm, and the shift corrections are applied in all the simulation.

Simulation
Procedure. The pore surfaces are set as rigid, and a certain amount of methane molecules are inserted into the corresponding nanopore according to the Table 1. 3 ns equilibrium molecular dynamic (EMD) simulation in NVT

Geofluids
ensemble is conducted to make the simulation system reach an equilibrium state [22]. In this EMD simulation process, the Nosé-Hoover thermostat with a relaxation time of 0.1 ps is utilized to control the temperature of the fluid system by coupling the CH 4 velocities in all x-, y-and z-directions. The following 1 ns simulation is sampled to obtain the methane density distribution confined in various confined nanopore.
To make the methane molecules reach a full-developed flow, the external force nonequilibrium molecular dynamic (EF-NEMD) simulation is adopted and one external force [8,16], and f is applied to each CH 4 particle. Such external force f has the functional formula with the equivalent pressure gradient along the flow direction (x-direction), as follows [22,55]: where ∇PðxÞ is the pressure gradient along x-direction and N CH 4 is the number of CH 4 molecules. In addition to the EF-NEMD simulation, DCV-GCMD and boundary-driven NEMD are also widely used to make the fluids reach the steady-state flow state [21,23,56,57]. Compared to the other methods, the EF-NEMD could be not only time and computer resources saving but also predicts the accurate geofluids flux under similar pressure gradient gradients very well [22,58]. In order to make sure the temperature is demonstrated by the thermal motion of gas molecules rather than the movement of center-of-mass when coupling the thermostat in EF-NEMD simulation [41,59], only the velocity in y-direction, which is perpendicular to the driving force, is considered to control the temperature of the con-fined geofluids [22,45,55]. The steady-state flow in all the simulation systems could be achieved after 5 ns usually; the following 15 ns is used to collect the methane flow velocity and density distributions for 3 times. The bins of thickness 0.02 nm are used to sampling the data (i.e., the density distributions and velocity profiles). It was already found that owing to lower friction of the perfect smooth atomic graphene surface, the alkane flow enhancement in the organic graphene nanopore could be 60 times faster than that in inorganic quartz nanopore [45]. Owing to such unique properties of the organic graphene surface and the inorganic quartz surface, the flow characteristic of the methane appears much more different, which will be discussed later. To make the flow velocity profile much more comparable, the equivalent pressure gradient ∇PðxÞ is set as the 0.04 MPa/nm in the quartz nanopore and 0.01 MPa/nm in the graphene nanopore. All MD simulation jobs are finished by the large-scale atomic/molecular massively parallel simulator (LAMMPS) [60], and the snapshots of the molecular system are viewed by the visual molecular dynamics (VMD) [61].

Results and Discussion
3.1. The Density and Velocity Distributions of Methane. In Figure 5, we show the snapshots of the methane confined in 5 nm, 10 nm, and 15 nm graphene nanopores. The density and velocity profiles of methane are also provided in Figure 5(d). The methane molecules form obvious layering structures near the graphene surfaces. The first density peak is 1.411 g/cm 3 in our simulation system under the reservoir condition 350 K and 50 MPa, which is line with the results at the same temperature from Nan et al., 1.334 g/cm 3 under the condition of 40 MPa and 1.452 g/cm 3 under 60 MPa [22]. The second density peak is 0.393 g/cm 3 . Such density values of the first peak and the second peak could be 5.52 and 1.54 times than that of bulk, ρ bulk . And such layer structures near the graphene surface are independent on the pore size as long as the reservoir pressure and temperature conditions are the same, which is in agreement with the methane confined in various graphene nanopore [22] and alkane confined in the nanopores [12,45,47]. In addition, the methane flow with the external force almost has a negligible influence on the methane density distributions; thus, the density profiles of CH 4 molecules at static and steady-state flow conditions coincide with each other exactly.
Although the CH 4 layering structure near the graphene surface are the same in the various pore size, while the methane velocity at the wall is various when the pore size is different, as shown in Figure 5(d), even though the pressure gradient is consistent with each other. It shows that the methane molecules could have a great slip velocity on the graphene surface, and such slip velocity could be up to 8.9 m/s, 16.3 m/s, and 28.6 m/s when the pore size is 5 nm, 10 nm, and 15 nm, respectively. Actually, in the slit-like graphene nanopore, such obvious surface slip velocity appears not only for CH 4 , but also in the flow of water, alkane, and some other geofluids [45,62,63]. The result also indicates that the pore size has nonnegligible effect on the slip velocity

Geofluids
when CH 4 is confined in organic graphene nanopore. With the pore size increases, the slip velocity of CH 4 increases. Similar tendency has been also reported for the alkane flow in graphene nanopore [22].
When it comes to the cases of methane confined in inorganic quartz nanopores, the configurations of the simulation system, the density, and the velocity profiles of methane are shown in Figure 6. The obvious methane layer structures are also formed near the quartz surface, while the density peak is much lower than that near the graphene surface. The first density peak is only 0.722 g/cm 3 , which is almost half of the value in the graphene nanopore. The second density peak is 0.331 g/cm 3 , which is lower than that in the graphene nanopore slightly. As to the velocity profile, the methane slip velocity on the quartz surface is almost zero, which indicates the slip flow mode plays negligible role for the methane flow in quartz nano-channels. Furthermore, unlike the slug-like velocity profiles in the graphene nanopore, the velocity profiles in quartz nanopore are parabolic, which are in agreement with the descriptions of the traditional Navier-Stokes equation [23,46]. Such different slip behavior of methane molecules on the quartz and graphene surface is interesting; the different liquid-solid interactions are generally regarded as the main reason, which is mainly dominated by the atom's parameters σ and ε [64,65] in Table 2 for the CH 4 molecules confined in various nanopores. In addition, the breakdown of gas molecules slippage on the quartz surface is especially owing to the inherent atomic roughness of the surface by Yu et al. [43,66], especially, the potential energy roughness of the quartz surface is higher than that on graphene surface up to two magnitude of orders. Thus, the gas molecules collide with the quartz surface and just go round and round and not move along the walls as what the molecules act on the graphene surface [40].
It is emphasized that although the surface types are different, such as the graphene and the quartz, the Kn numbers are the same when the methane molecules flow under the same reservoir pressure (50 MPa) and temperature (350 K) conditions. The Kn numbers are around 0.1 in these nanopores of this work, and the slip flow mode is obvious only in the organic graphene nanopore; meanwhile, the viscous continuum flow still dominant the flow behavior in the quartz nanopores. Thus, the surface type is also needed to be considered carefully in the prediction of the flow performance of the gas in nanoscale. From the view of shale gas in which J MD is the flux obtained by the MD simulation. And then, referring to Equation (3), the enhancement factor could be given by Here, when calculating the J NS , the effective pore width W needs to be provided. Suggested by Botan et al. and Bourg et al. [68,69], the position where the Gibbs dividing surface, z surf is set as the wall surface exactly, according to Then, the distance between both Gibbs dividing surface is defined as the effective pore width, W. Therefore, the slip length L s could be obtained using Equation (4). Figure 7 presents the enhancement factor and the slip length of methane flow confined in graphene and quartz nanopores with various pore size. The enhancement of methane flow in organic nanopore is very obvious, and ε could up to be 14.4, 7.2, and 5.1 when the pore size are 5 nm, 10 nm, and 15 nm, respectively. With the increase of the pore size, such enhancement becomes more and more weaker. At the same time, when the pore size is ranging from 5 nm to 15 nm, the slip length L s is around 10 nm, which is comparable to the pore size. Also, with the pore size increases from 5 nm to 10 nm, L s decreases slightly firstly, then L s tends to be a constant value, which is around 10 nm. Similar tendency The equivalent pressure gradient from the inlet to the outlet of the quartz nanopore is 0.04 MPa/nm. 7 Geofluids of slip length is also found for the methane flow confined in the graphene nanopore under high-pressure conditions (i.e., 40 MPa and 60 MPa) [22]. Such obvious enhancement is mainly attributed to the high adsorption and slip flow of the gas molecules on the graphene surface [24]. Whereas for the methane molecules confined in the inorganic quartz nanopores, such enhancement flow disappears. The ε is closely to 1, and L s is almost 0, which indicates that the traditional Navier-Stokes equation combined with the no-slip boundary condition could capture the methane flow behavior in the quartz nanopores fully. The attribution of the surface diffusion or the slippage on surface could be negligible when the methane molecules transport in the inorganic quartz nanopores.

Flux Contribution of Each Region.
Combined with the density distributions and the flow region dividing criteria [24], the methane flow region is divided into three types, adsorption layer, Knudsen layer, and the bulk gas region. According to the migration characteristic of the gas molecules confined in the nanopores [24,70], the first adsorption layer is defined as the adsorption layer. The region far away the wall is classified as the bulk gas region where the gas shows bulk density continuously. And the region between the adsorption layer and the bulk region is defined as the Knudsen layer. Given the pore type (i.e., graphene or quartz) and the reservoir condition (i.e., pressure and temperature), it is found that the thickness of the adsorption layer and Knudsen layer, which is noted as h Ad and h Kn , respectively, is a constant, and only the thickness of the bulk gas region h bulk changes with the pore size. Thus, only the clarification of the flow region in 5 nm pores is presented here. As depicted in Figure 8, the h Ad and h Kn in graphene nanopores are always 0.26 nm and 0.76 nm, respectively. Similarly, the h Ad and h Kn in the quartz nanopores are 0.28 nm and 0.78 nm, respectively. Such a tiny deviation may owe to the different interaction between the surface type and the methane molecules; in addition, the precision of the bin scale (0.02 nm) also has a limit to the characterize the region boundary. Overall, it is reasonable to assume the thickness of the adsorption layer as the diameter of the methane molecules, and the width of the Knudsen layer is around 0.77 nm.
In Figure 9, we present the flux contribution ratio of each flow region in the organic graphene and inorganic quartz nanopores. Firstly, the flux of each region J i is evaluated by the calculation [55],  where z i and z j are the boundaries of each region. And then, the flux ratios of each region could be obtained by J MD i /J MD . In the 5 nm graphene nanopore, the contributions of the adsorption layer and Knudsen layer could be up to 11.5% and 32.72%, and the total amount of both regions is almost half of the total flux, which could not be neglected. While in the 5 nm quartz nanopore, the total contribution is only around 30%. With the increase of the pore size, the flux contribution ratio of the adsorption layer and the Knudsen layer decreases, and the percent of the bulk gas increases. When the pore size is 15 nm, the total contributions of the adsorbed layer and the Knudsen layer are 12.52% in the graphene nanopore and 3.51% in the quartz nanopore. Thus, the flux of the adsorption layer and Knudsen layer could be neglected in the inorganic nanopore when the pore size is larger than 15 nm. While if we do not consider the flux contribution of the adsorption layer and Knudsen layer, the flux would be underestimated by 12.5% in the 15 nm graphene nanopore.

Conclusions
In this work, based on molecular dynamic simulation, we presented the density distributions and the flow performance of methane in organic graphene nanopores and inorganic quartz nanopores. It was observed that the layering structure forms on all the graphene and quartz surface, and the peak density of the first adsorption layer on the former is much higher than that on the latter. The slippage of the methane molecules is obvious on the graphene surface while disappears on the quartz surface. Owing to such great slip on the graphene surface, the enhancement to the gas flux is nonnegligible in the organic graphene nanopores, and the enhancement factor could be up to 14.43, 7.16, and 5.12 when the pore sizes are 5 nm, 10 nm, and 15 nm, respectively. Correspondingly, the slip length is always larger than 10 nm when the methane flow in graphene nanochannels, which is comparable to the pore size. While in the quartz nanopores, the enhancement to the gas flow is almost none, and the slip length is around zero. In addition, the flux contributions of the adsorbed layer, the Knudsen layer, and the bulk gas region were also presented. The flux contribution ratios of the adsorbed layer and the Knudsen layer in the quartz nanopore are much lower than those in the graphene nanopores when the pore size is the same. The total flux contribution of the adsorbed layer and the Knudsen layer is lower than 10% in the quartz nanopore if the pore size is larger than 10 nm, which could be neglected in the prediction of the gas flux confined in quartz nanopores. It was also noted that in addition to the Knudsen number, the surface type also has effect on the characteristics of the gas migration in shale nanopores. Our work could provide a molecular view and important reference to the prediction of the gas flow in various types of shale nanopores.

Data Availability
All the data has been presented in the manuscript.

Conflicts of Interest
The authors declare that they have no conflicts of interest.