Gas Transport in Shale Nanopores with Miscible Zone

Based on the results of molecular dynamics simulation, in a gas-water miscible zone, the velocity profiles of the flowing water film do not increase monotonously but increase first and then decrease, which is due to the interaction between water and gas molecules. This exhibits a new physical mechanism. In this paper, we firstly propose a gas-water flow model that takes into account the new physical phenomena and describes the distribution of gas-water velocity in the whole pore more accurately. In this model, a decreasing factor for water film in the gas-water miscible zone is used to describe the decrease of water velocity in the gas-water miscible zone, which leads to the gas velocity decrease correspondingly. The new flow model considers the interaction among gas and water molecules in the miscible zone and can provide more accurate velocity profiles compared with the flow models not considering the miscible region. Comparison calculation shows that the previous model overestimates the flow velocity, and the overestimation increases with the decrease of the pore radius. Based on the new gas-water flow model, a new permeability correction factor is deduced to consider the interaction among gas and water molecules.


Introduction
Shale formation is a great energy source. As a result of rising energy demand and prices, research into unconventional oil and gas is deepening, and successful development of shale formation in the United States has led to a boom in research and development [1]. Nowadays, there are about 213 tcm of recoverable shale gas reserves in the world [2]. Shale gas is a very important unconventional energy, and its exploitation capacity is mainly determined by gas storage and flow capacity [3,4]. Shale gas has a big advantage over other energy sources in terms of carbon dioxide emissions and prices and is becoming a bigger and bigger part of total energy production [5,6]. However, because shale is mainly stored in natural fractures, micropores, nanopores, and other complex and relatively impermeable shale formations, its exploitation has many technical and flow mechanism problems [1,[7][8][9]. Therefore, it is important to study the flow mechanism in shale formation in order to provide more theoretical support for the improvement of mining technology and efficiency.
Shale formation pores are mainly nanoscale. The pore radius varies widely and is mainly on the nanometer scale, which is quite different from the transport mechanism in conventional formation and has a significant effect on the gas flow performance [1,7]. When considering the transport of single-phase gas in the nanometer pore channel, there is slip effect on the wall surface due to the interaction between solid and gas molecules [10]. Surface diffusion of adsorbed gas cannot be neglected in smaller pores and increases significantly with the increase of the maximum adsorption capacity [11]. By defining the ratio of the average free path of gas to the characteristic dimension as the Knudsen number K n , the gas transport mechanism can be divided into viscous flow, slip flow, transition flow, and molecular free flow [12][13][14]. Therefore, when establishing the flow equation, the flow boundary cannot be regarded as a nonslip boundary, and it needs to be corrected by adding gas slip [15]. Based on the Knudsen number K n , two kinds of models are established: continuity equation obtained by converting the boundary condition to nonslip boundary [16], several transport mechanisms combined by weight coefficients with a unified gas transport equation [17,18]. Based on a unified diffusion coefficient, Cai et al. [19] proposed an apparent permeability model of shale which investigated the gas transport mechanism in a shale nanopore by considering convective flow, gas diffusion, and surface diffusion. However, the above results only apply to single-phase gas transport. In a real shale reservoir, the injection of fracturing fluid and initial water saturation make the real pore flow more than a single-phase gas flow [20]. Therefore, when studying nanopores, the gas-water two-phase velocity model should be established, which makes the model more consistent with the actual situation.
In order to establish a more accurate gas flow model in nanometer pores, a lot of literatures have discussed gas and water flow models. Li et al. [21] proposed an analytical method based on the Hagen-Poiseuille formula and capillary pressure curve to model and analyze the gas-water relative permeability of nanopores with interface effect. Singh and Cai [8] proposed a new method which predicts permeability of fractured shale by discretizing the medium into matrix and can estimate permeability at any scale of interest and is used to predict relative permeability estimates of two-phase flow in fractured shale samples. Li et al. [22] proposed an analytical model to consider the flowing water film by changing the boundary conditions of the gas-liquid interface and the solid-liquid interface. The model has two different nanoscale characteristic sizes of crack hole and circular hole and is in good agreement with the experimental data, which proves that the high-viscosity flowing water film can improve the flow capacity of gas.
However, the gas-water two-phase flow in nanopores is more complex than the above literatures. Molecular dynamics simulation plays an important role in studying the mechanism of gas-water flow in pores. Xu et al. [23] performed a comprehensive study on the two-phase transport characteristics of shale gas and water through hydrophilic and hydrophobic nanopores combined with molecular dynamics (MD) simulation and analysis model. Hao et al. [24] used nonequilibrium molecular dynamics to simulate the mixed flow behavior of water and methane gas in shale pores. Based on Hao et al.'s research, we found that the velocity profiles of the gas-water in the nanometer pore are different from the previous studies, as shown in Figure 1. At the gas-water interface, due to the interaction between gas and water molecules, the water film velocity profile has an obvious downward trend near the gas-water interface, which is partly enlarged as shown in Figure 1(a), but the water velocity profile in a previous paper is increasing monotonically as shown in Figure 1(b) [25]. Based on the molecular dynamics results, there are gas and water molecules in a zone near the gaswater interface which is called the miscible zone in the next section. The effect of the miscible zone on flow behavior is not considered in the above-mentioned models. This shows that there are some physical mechanisms that are not considered in previous papers.
Based on the research gap, in this paper, a gas-water twophase model was proposed to study the influence of the miscible zone at the gas-water interface. The flow model takes a decreasing factor into account to describe the effect of the miscible zone on flow behaviors.

Model Establishment
2.1. Miscible Zone. Figure 2 gives the density distribution of gas and water molecules in pores, the density 2D distribution of gas and water molecules, and an elevation of water and gas at equilibrium in a pore in Hao et al.'s molecular dynamics simulation results [20]. As shown in Figure 2(a), at the gaswater interface, there is a region where gas-water molecules coexist. Figures 2(b) and 2(c) show the density 2D distribution of gas and water molecules. Near the pore wall are all water molecules, and in the middle of the pore are all gas molecules. At the gas-water interface, there is a region where the density of gas and water molecules is not zero, and in this region, the gas and water molecules coexist. In order to show the region where gas-water molecules coexist clearly, the 1D density distribution of gas and water is shown in Figures 2(d) and 2(e), which show that there are three zones: zone 1, zone 2, and zone 3. Only water molecules are in zone 1, which represents the water film. In zone 2, there are gas and water molecules at the same time. The number of water molecules decreases while gas molecules increase (indicated by the shaded part). Combined with the analysis of 1D and 2D density distribution of gas-water molecules, water and gas molecules form interfacial regions with the 1-99 thickness (1%-99% of bulk gas density) around 0.3 nm [23,26], and we call the interfacial regions as the miscible zone (zone 2). In zone 3, the water molecule number of zone 3 becomes less and is about five or six times less than that in zone 2, but the number of gas molecules reaches its maximum. Therefore, in zone 3, gas molecules are those that dominated and it is considered as pure gas phase flow. Figure 2(b) shows that the density distribution of water and gas molecules is not strictly symmetrical. In order to model the flow behavior, Figure 2 is simplified to Figure 3, in which there are three flow zones in the circular pore.
They are a high-viscosity water film zone near the wall, a gas zone in the middle, and a miscible zone between the gas zone and the water film zone. Hao et al.'s results show that due to the gas-water molecule interaction, the farther away from the wall of the pore, the slower the velocity of water, as shown in Figure 2. In this paper, we use a decreasing factor β to describe the effect of molecule interaction in the miscible region on the water velocity as follows: where h m and h w are the thickness of the miscible zone and water film zone, respectively, and r 0 is the pore radius.

Mobile
High-Viscosity Water Film. The walls of the nanopores are mainly composed of hydrophilic or hydrophobic materials [27]. However, because of the diversity of wall materials, hydrophobic substances will also be doped. Due to the solid attraction of the hydrophilic channel wall, the water molecules are trapped on the surface of the hydrophilic solid and arranged in an orderly manner within a few molecular diameters near the channel wall [27,28]. A large number of molecular dynamics simulation results and experimental data illustrate that the thickness of water film is about 0.7 nm [29][30][31]. When the pressure gradient in the pores reaches a certain height, the water film will flow. Such a water film will show the characteristics of high viscosity and slow flow rate, which is of great significance in establishing the gas-liquid two-phase velocity model in nanometer pores.
Because of the interaction of gas and water molecule, the real slip of confined water can be calculated as [22] where θ is the contact angle, l s is the slip length, and C is assigned to be 0.41, dimensionless.

Flow Equation.
This model is based on the Hagen-Poiseuille equation for steady-state laminar flow through circular pores. As shown in Figure 3, the high-viscosity water film is distributed on the pore surface evenly, and the miscible zone exists at the gas-water interface. Based on the model of Mattia and Calabrò [25], gas and water velocity are, respectively, where μ w and μ g are viscosity of water and gas, respectively; L is the length of the pore; and ΔP is the pressure difference between entrance and exit. The boundary conditions for velocity continuity are as follows: where v s describes the gas-water momentum transport and the interaction between gas molecules, which can be defined as [32] v where v s is the slip velocity between water and gas, σ v is the tangential momentum accommodation coefficient, b 4 Geofluids is a slip coefficient, and λ is the mean free path, which can be defined as where R is the gas constant, M is the molecular weight, T is the temperature, and Z is the gas compressibility factor, which is calculated [22]: where T c is the critical temperature and P c is the critical pressure. (e) z direction Figure 2: The elevation of water and gas at equilibrium in a pore, the density 2D distribution of gas and water molecules, and an average density distribution curve of gas and water in the yz plane under a static equilibrium state of equilibrium molecular dynamics [24]. (a) The elevation of water and gas at equilibrium in a pore. The gray ball represents the carbon atoms in the methane molecule, which can represent gas molecules. The green ball represents the oxygen atoms in the water molecule, which can represent water molecules. (b) The density 2D distribution of gas. (c) The density 2D distribution of water. (d) The density distribution curve of water and gas molecules along the y direction, and (e) is the density distribution curve of water and gas along the z direction. The shaded areas are estimated to be highly interactive between the gas and water molecules.

Geofluids
By the combination of Equations (3) and (4), the velocity profiles are obtained: By integrating Equation (8) along the r direction, the gas flow equation can be further deduced as where Q g is the flow rate of gas in nanopores. In Darcy's equation, the flow rate Q d of gas is as follows: The flow rate of porous media is modified by the ratio of tortuosity to porosity [22]: Substituting Equations (9) and (10) into Equation (11) to get the formula for calculating the gas apparent permeability, where K Ag is the gas apparent permeability when considering the flow of water film in the nanopores. When the miscible zone and flowing water film with high viscosity are not considered, it reduces to a single gas flow model, and its boundary condition is given as follows: Using the same derivation process in the formula of velocity profile and apparent permeability of gas without considering the miscible zone and flowing high-viscosity water film can be deduced as follows: where V Agi is the gas velocity in the nanometer pores without considering the miscible zone and the mobile water film and K Agi is the apparent permeability without considering the miscible zone and the high-viscosity flowing water film.

Model Validation
In order to verify the model in this paper, a comparison was made between the velocity profile of gas-water calculated by the proposed model and the result obtained by molecular dynamics simulation. Figure 4 shows the velocity profile of gas and water in a pore that contains a local magnification of the velocity of the water film on the left and right sides. The velocity profile of the gas in the middle of the pore is parabolic, while the velocity of the water film near the pore wall first rises and then falls at the gas-water interface. Comparing Figure 1 with Figure  and the proposed model in this paper with the pore radius of 2 nm, 3 nm, and 5 nm, respectively. It shows that gas flow velocity of the single gas model is smaller than that of the proposed gas-water flow model. We can see that when the pore radius increases, the difference of gas velocity profile between the proposed model and single gas flow becomes larger. When the pore radius is 2 nm, their gas velocity is approximately equal. For the smaller nanopores, although the water film thickness occupies a relatively large pore radius and enhances gas velocity at the interface of water film and gas, the miscible zone reduces the gas velocity, so that the combined effect is not obvious.
However, when the pore radius is 5 nm, the difference of gas velocity profile between the proposed model and single gas flow model becomes more obvious. According to Equation (1), we know that the value of the decreasing factor at the interface of water film and gas increases with the increase of the pore radius. The mobile water film increases the gas flow capacity bigger, which leads to larger gas velocity profile.

Gas Transport Capacity.
In this paper, the ratio of K Ag to K Agi which is defined as the apparent permeability enhancement factor K c , is used to evaluate the enhanced gas flow capacity by considering the miscible zone and high-viscosity flowing water film compared with the single gas model. Figure 6(a) shows that the value of K c is always larger than 1 and that K c increases with the increase of the pressure,   8 Geofluids which means that the gas transport capacity is underestimated if the miscible zone and mobile water film are not considered. Figure 6(b) shows that K c is smaller for larger pore radius. When flowing water film and miscible zones are con-sidered coherently, enhancement factor K c is very small, for example, K c = 1:013 for r 0 = 50 nm. The mobile water film leads to the increase of the apparent permeability enhancement factor. The miscible zone lowers the apparent permeability enhancement factor. The two factors are canceled   [22]. The reason is that the miscible zone lowers the gas velocity and thus lowers the apparent permeability enhancement factor.
We can see that when the pore radius increases, the difference of gas velocity profile between the proposed model and Li et al.'s model becomes smaller. This means that for the smaller nanopores, the miscible zone has larger effect on the gas velocity. Figure 5(d) gives the relative differences of gas velocity between Li et al.'s model and the proposed model in this paper. It clearly shows that error would increase for the smaller nanopores when the miscible zone is neglected. At the centerline, the differences are 2.85%, 1.5%, and 0.82%

Apparent Permeability Enhancement
Factor. In order to exhibit the effect of the miscible zone on permeability correction under the mobile water film, Figure 7 gives the comparison of the flow enhancement factor between with and without consideration of the miscible zone. When the miscible zone is neglected, the gas velocity will be larger at the interface of mobile water film and gas. This leads to a larger flow enhancement factor as shown in Figure 7.

Gas Velocity Profile.
In order to exhibit the effect of the miscible zone on gas velocity profile under the flowing water film, Figure 5 gives the water velocity profile comparison between with and without consideration of the miscible zone.
The gas velocity at the center line is the highest, and the overall gas velocity is parabolic. The gas velocity profile comparison in Figure 5 shows that gas velocity is smaller when the miscible zone is considered. The reason is that when the miscible zone is considered, the water velocity profile near the boundary increases first and then decreases (Figure 1(b)) and is like a parabola. According to boundary conditions in Equation (13), the gas velocity at the boundary will be smaller compared with the velocity of that not considering the miscible zone, which will lead to a smaller gas velocity profile.

Water Flow Rate.
In order to exhibit the effect of the miscible zone on the water flow rate under the flowing water film, Figure 9 gives the gas velocity comparison between with and without consideration of the miscible zone. As shown in Figure 9, the water flow rate in the proposed model is smaller than that in Li et al.'s paper, but all of them increase with the increase of pressure. In this paper, we use the proposed model to describe the velocity distribution of water film, which is different from Li et al.'s model. Due to the existence of the miscible zone, there is an obvious downward trend around the gas-water interface instead of a monotonic rise. Therefore, the average velocity of water film in this paper is less than the average velocity of water film not considering the miscible zone in Li et al.'s paper. This indicates that the water flow rate of the proposed model is smaller than that of Li et al.'s model. This means that the water flowback ratio is smaller due to the existence of the miscible zone.

Conclusions
In shale nanopores, the miscible zone and mobile highviscosity water film are potential influencing factors for gas transport. Inspired by the results of molecular dynamics simulation, the miscible zone is considered into the gas flow model in this paper. The calculation results of Li et al.'s model are compared to verify that the model in this paper is more comprehensive. According to the results of the study and discussion, the following conclusions are drawn: (1) Simulating the gas flow pattern in actual shale formations, miscible zones should be taken into account. Ignoring the influence of miscible regions leads to an overestimation of the velocity of the gas in the pores. If the flow of water film and miscible zone are ignored, the velocity of gas in pores will be underestimated (2) The flow enhancement factor reflects the transport capacity of shale gas under different conditions, which can be concluded as follows: in the case of larger pressure and smaller pores, the enhancement factor is larger and the flow enhancement is more significant (3) By comparing with Li et al.'s paper and single gas model in the aspects of velocity, gas flow rate, water flow rate, and enhancement factor, it is verified that ignoring the miscible zone and only considering the mobile high-viscosity water film will overestimate the flow rate of water and gas and the flow capacity of gas will also be overestimated in the calculation process The above results show that proposing the miscible zone has a significance in modifying the gas flow model. In the study of the miscible zone, the paper only discusses the indirect influence on gas velocity through the influence on water film velocity and the value of the thickness cannot be used to every situation. Therefore, the specific mechanism of action of the miscible zone needs to be explained, and the calculation of the thickness of the miscible zone needs to be more accurate which will be the focus of the next work.

Nomenclature b:
Slip coefficient, dimensionless C: Constant, dimensionless h m : Thickness of the miscible zone (m) h w : Thickness of the water film zone (m) K Ag : Gas apparent permeability with the miscible zone and mobile high-viscosity water film ðm 2 Þ K Agi : Gas apparent permeability with the static highviscosity water film (m 2 ) K c : Enhancement factor, dimensionless L: Length of the pore model (m) l s : True slip length (m) M: Molecular weight (kg/mol) P c : Critical pressure (Pa) ΔP: Pressure difference between the entrance and exit (Pa) Q d : Darcy rate of gas flow (m 3 /s) Q g : Flow rate of gas (m 3 /s) r 0 : Pore radius (m) R: Gas constant (8.314 J/(K · mol)) T: Temperature (K) T c : Critical temperature (K) V Agi : Gas velocity with the static high-viscosity water film (m/s) v s : Gas slip velocity (m/s) Z: Gas compressibility factor, dimensionless λ: Mean free path (m)