A novel numerical model of gas transport in multiscale shale gas reservoirs with considering surface diffusion and Langmuir slip conditions

Multiflow mechanisms coexist in shale gas reservoirs (SGRs) due to the abundant nanopores and the organic matter as a medium of gas souring and storage. The gas transport mechanisms in nanopores including bulk gas transfer and adsorption‐gas surface diffusion were already investigated in pore‐scale models, but their effects on actual gas production of multistage fractured horizontal wells in SGRs are not clearly understood, which are crucial for the economic development of unconventional resources. Therefore, a comprehensive apparent permeability (AP) model which couples the surface diffusion of adsorbed gas, slippage flow considering the additional flux generated by surface diffusion based on Langmuir's theory, and Knudsen diffusion is established. The presented model is validated with the experimental data and lattice Boltzmann method (LBM) simulation results. Then, we propose a numerical model which combines multiflow mechanisms in microscale pores and a multistage fractured horizontal well (MSFHW) in macroscale shale gas reservoirs together. The effects of different transport mechanisms on both AP of nanopores and gas production are analyzed thoroughly. The results show that the effect of surface diffusion on the apparent permeability of nanopores is much greater than that on the actual gas production of MSFHW, and the influence of high‐pressure condition must be considered when calculating the surface diffusion coefficient. The presented numerical model has important implications for accurate numerical simulation and efficient development of shale gas reservoirs.


| INTRODUCTION
The gas transfer mechanisms in shale gas reservoirs (SGRs) are significantly different from that of conventional gas reservoirs, which is mainly caused by the abundant nanopores in shales and the organic matter as a medium of gas sourcing and storage. So multiflow mechanisms coexist when gas transport in shales, thus a comprehensive model for gas transfer in nanopores is essential for accurate numerical simulation and effective development of shale gas reservoirs. Nanopores of organic matter are shown to be abundant in shales, which have a large surface area and contain a significant amount of adsorbed gas. The surface diffusion of adsorbed molecules occurs under the gradient of chemical potential, which may significantly contribute to mass transfer in porous media. 1,2 The mobility of adsorbed benzophenone molecules was first demonstrated by Volmer and Adhikari through experiments, and then, various analytical and simulation models of surface diffusion have been studied by researchers. The flux equation of surface diffusion was proposed by Do et al 2 for hydrocarbons in activated carbon. Sheng et al 3 investigated the effect of surface diffusion on shale gas flow in Kerogen pores based on Do's model, but they applied a constant value to represent the surface diffusion coefficient. Hwang et al 4 and Guo et al 5 both proposed the equation of surface diffusion coefficient, which is influenced by temperature, adsorbent, and adsorbate. However, the influence of pressure was not considered and cannot be applied for surface diffusion calculation of adsorbed gas for a high-pressure condition. Wu et al 6 derived the model of surface diffusion under a low-pressure condition and considered the effect of adsorbed gas coverage under high pressure based on Hwang's model, and they found that the contribution of surface diffusion to the gas mass transfer is dominant in micropores. These researchers mainly focused on porescale modeling or rate-transient analysis. 7 Nevertheless, little work has been done to investigate the effect of surface diffusion on gas production of MSFHWs for shale gas reservoirs.
In addition, the mechanisms of free gas flow in nanopores are the key issue in shale gas reservoir research, which are viscous flow and Knudsen diffusion. 8,9 And the viscous flow must be modified to satisfy the slip boundaries due to the breakdown of the continuum assumption in nanopores. An empirical formula or a theoretical coefficient F 10-12 is used to account for the slippage effect. Another type of slip model is based on slip boundary conditions, including the first-and second-order modifications ( Table 1). The Maxwell slip model 13 is a first-order approximation from kinetic theory, and second-order [14][15][16][17][18] slip models have been successively proposed to improve calculation precision. In summary, slip velocity models can uniformly be expressed as 19 : Table 1 shows the different slip velocity models and the disadvantages when these models are applied in shale system. To avoid these disadvantages and their impact on the permeability models developed based on Maxwell slip model and other second-order slip models, a Langmuir slip boundary condition more accurately describes the real slip phenomenon from the perspective of a physical meaning. Langmuir slip condition considering the gas-surface molecular interaction based on the theory of surface chemistry was first proposed by Myong 20 and applied to the permeability model for shale gas reservoir by Singh and Javadpour. 21 But many scholars 21,24-29 did not fully consider Knudsen diffusion, surface diffusion and slip flow considering Langmuir slip condition in their apparent permeability (AP) model for nanopores of shales.
Although most studies proposed the pore-scale models for gas transfer in nanopores of shales which couple bulk gas transfer and adsorption-gas surface diffusion, little work has been done to investigate these gas transfer mechanisms on actual gas production of MSFHWs in macroscale shale gas reservoir. And the effects of the gas transfer mechanisms on the permeability of a single nanopore are different from that on gas production, and the latter one will be more complicated. There are two methods to study the rate-transient analysis of MSFHWs, which include analytical methods (the analytical trilinear flow model 30,31 and semi-analytical models 32-34 based on Green's function and the source/sink method) and numerical simulation method. But analytical methods require many simplification assumptions which are not completely based on the actual physical model compared

T A B L E 1 Different slip velocity models
Note: C 1 and C 2 are two slip coefficients, and σ v is the tangential momentum accommodation coefficient.
with the numerical models, so researchers [35][36][37][38] usually use numerical simulation methods or commercial software to study the performance of MSFHWs. In this work, a comprehensive AP model which couples the surface diffusion of adsorbed gas and its effect on slippage flow, Knudsen diffusion with considering solid surface roughness is proposed. The AP model is validated with the experimental data and lattice Boltzmann method (LBM) simulation results. Then, a multiscale numerical model is established by combining the AP model that considers gas transfer mechanisms in microscale pores with a multistage fractured horizontal well (MSFHW) in a macroscale shale gas reservoir. The effects of different gas transport mechanisms both on apparent permeability of nanopores and gas production of a MSFHW were both analyzed thoroughly. The presented numerical model has important implications for evaluating and predicting the productivity of fractured wells, which is of great significance for effective and economic development of shale gas reservoirs.

SHALE NANOPORES
Shale gas reservoir is self-generating and self-storing, so free gas and adsorbed gas coexist in the reservoir. Due to the submicron effects of the nanopores, bulk phase free gas is mainly transferred by slippage viscous flow under pressure gradient and Knudsen diffusion under the impact of collisions between molecules and solid walls ( Figure 1). Besides, adsorbed gas molecules can hop from one adsorption site to another by the means of surface diffusion, which was validated with Lattice Boltzmann simulation. 39

| Surface diffusion
Surface diffusion is mainly caused by the interaction between adsorbed gas molecules and solid surface, which is a continuous process in which gas molecules randomly hop from one adsorption site to another under the action of activation energy, as shown in Figure 2.
According to the Maxwell-Stefan method, the driving force of surface diffusion is the chemical potential gradient of the adsorbed phase 40 : The chemical potential can be written as 40,41 : Combined with Equation (3), Equation (2) can be written as: where D 0 s = L m RT. Based on the assumption that Langmuir adsorption is monolayer adsorption, the coverage can be expressed as: The number of adsorbed gas molecules per unit surface area can be written as: The mass of adsorbed gas per unit surface area can be written as: The concentration of adsorbed gas can be obtained due to the single adsorbed layer: Combined with Equations (6), (7) and (8), the concentration of adsorbed gas can be written as: The surface diffusion coefficient can be also expressed as 5 : (2) J s = −L m s C s u x .
(3) u = u o + RT ln p. The surface diffusion coefficient in Equation (10) is obtained through experiments under low pressure, which is only applicable to describe the surface diffusion under a low-pressure condition. According to the kinetic method proposed by Chen and Yang, 42 the surface diffusion under a high-pressure condition is described by considering the influence of gas coverage on surface diffusion: where H(1−κ) is a Heaviside function, dimensionless.
Therefore, the surface diffusion mass flux in shale nanopores can be expressed as:

Langmuir's theory
Since gas molecules are adsorbed on pore walls, the adsorbed layer becomes the flow boundary of bulk phase gas transport, and the moving molecules of adsorbed layer accelerate the gas slippage and further increase the viscous flow.
Considering laminar flow in a capillary tube with radius r t , the velocity profile can be derived from Navier-Stokes Equation 4 : where C 1 and C 2 are integral constants obtained from boundary conditions.
The boundary conditions of gas flow in the tube can be expressed as: By combining Equations (15) and (16), we can get: So the velocity profile can be obtained as: The slippage velocity was derived by Myong 20 on the basis of Langmuir's theory: The velocity of gas adsorbed on the pore surface can be derived from Equation (14): The local velocity a mean free path away from the pore surface can be obtained according to Equation (18): So the viscous flux with considering slippage model based on Langmuir's theory in a single pore can be derived as: We can get the equivalent permeability according to the viscous flux with considering slippage velocity model: where Kn is the Knudsen number, which can be expressed as Kn = λ/r t .
So the permeability modification coefficient F(p, Kn) can be obtained through Equation (23): where the first term in Equation (24) represents viscous flow, the second and the third terms represent slip flow near the solid surface, and the last term is caused by the effect of moving boundary of adsorbed layer.
The viscous flux with considering slippage velocity model in shale nanopores can be written as:

| Knudsen diffusion with considering solid surface roughness
As the pore diameter decreases or molecular mean free path increases (at low pressures), gas molecules tend to collide more with the pore wall compared with the collisions between gas molecules. The Knudsen diffusion flux can be written as: According to the gas state equation, the gas concentration can be expressed as: The influence of solid surface roughness on Knudsen diffusion cannot be negligible, because the Knudsen diffusion decreased due to the long gas retention in the vicinity of a pore wall under severe wall roughness. Combined with Equation (27), the Knudsen diffusion flux with considering solid surface roughness 19,25 in shale nanopores can be written as: The Knudsen diffusion coefficient is defined as 43 :

| Pore confinement and real gas effect
Knudsen number is a function of the mean free path of gas molecules and equivalent hydraulic radius, which is defined as: The gas compressibility factor can be calculated by 44 : The viscosity of real gas can be calculated by 44 : So the Knudsen number of real gas is obtained as 45 : The parameters used in Equations (31) and (32) are listed in Table 2.

| Adsorbed layer on pore radius
The effective radius can be calculated based on the assumption of monolayer adsorption on solid surfaces: So the effective permeability considering adsorption layer thickness can be written as:

| Apparent permeability modeling
The contributions of slippage viscous flow and Knudsen diffusion to total mass flux are not a linear assembly, which is influenced by the ratio of intermolecular collision frequency to the total collision frequency. 46 And surface diffusion contributes to total mass flux independently because the surface diffusion is not caused by molecular collisions at the pore wall. Thus, the total mass flux can be written as: ω is the ratio of intermolecular collision frequency to total collision frequency, which is the function of Knudsen number 3,46 : For surface diffusion, the equivalent permeability can be written as: For viscous flow with considering slippage velocity model based on Langmuir's theory, the equivalent permeability can be written as: For Knudsen diffusion with considering solid surface roughness, the equivalent permeability can be written as: Therefore, the apparent permeability of shale nanopores can be expressed as the sum of bulk gas transport and adsorption-gas surface diffusion according to Equation (36):

| Validation of the AP model
To examine the accuracy of our AP model, the permeability modification coefficient F(p, Kn) is validated with the experimental data and LBM simulation results proposed by Fathi et al, 47 and the experimental results are obtained from gas uptake measurements using crushed nanoporous material. The parameters required to validate the slip model are listed in Table 3. We can see from Figure 3 that the calculated F(p, Kn) is proved to match well with the experimental data and the LBM results of Fathi et al.
Then, our AP both considering bulk phase transporting and surface diffusion and AP only considering bulk phase transport are validated with the generalized lattice Boltzmann model (LBM) simulation results of Wang et al, 26 and the parameters required to validate the AP model are listed in Table 4. As can be seen from Figure 4, the results of our AP model match the LBM results very well under different pore size regardless of whether the surface diffusion is taken into account. We also find that surface diffusion begins to affect AP when the pore radius is <20 nm.  (41) is only applicable to calculate the apparent permeability of a single nanopore. The AP of shale reservoir should be modified as follows.

| SIMULATION OF SHALE GAS
The surface diffusion mass flux of shale reservoir can be expressed as: where γ s is the correction factor of surface diffusion of adsorbed gas in shale reservoir, which is expressed as 6 : The viscous flux with considering slippage velocity model in shale reservoir can be written as: where γ b is the correction factor of AP of bulk gas in shale reservoir, γ b = ϕ/τ. The Knudsen diffusion flux with considering solid surface roughness in shale reservoir can be written as: Therefore, the apparent permeability of shale reservoir can be expressed as the sum of bulk gas transport and adsorbed gas diffusion according to Equation (36):  The volume of adsorbed gas in shale reservoir can be calculated by: where k app in Equation (47) can be calculated by Equation (46), the apparent permeability model in shale reservoir which considers surface diffusion of adsorbed gas and bulk gas transport.

| Gas flow in hydraulic fractures
The effective conductivity of both propped and unpropped fractures declines with the decrease of effective stress, 48 and the relationship between the effective permeability and pressure can be calculated by the following equation based on the experimental results of Zhang et al 49 : By applying the theory of equivalent fracture conductivity, the effective permeability of hydraulic fracture which considers the stress-dependence effect can be calculated by 48 : Gas flow in hydraulic fractures can be described by Darcy' law, so the following equation is achieved:

| RESULTS AND DISCUSSION
The numerical simulator for shale gas reservoirs was validated in our previous work, 50 and the parameters needed for simulation are listed in Table 5. In this part, the effects of different gas transport mechanisms such as bulk gas transport and surface diffusion on apparent permeability of nanopores have been studied, and the effects of surface diffusion under different pore radius, surface diffusion coefficient and stimulated reservoir volume on the production rate, and cumulative production of multifractured horizontal well are analyzed thoroughly. Figure 5 shows the ratios of apparent permeabilities of different gas transport mechanisms (surface diffusion, Knudsen diffusion, and viscous flow with considering slippage effect) to total apparent permeability k app . From Figure 5A, we can find that the ratio of the AP of surface diffusion to total AP (k sur /k app ) increases with the decrease of pore radius and increases rapidly when pore radius is <10 nm. Surface diffusion becomes the main gas transport type when pore radius gets smaller than 5 nm, and the value of k sur /k app almost reaches 1 when the pore radius decreases to 1 nm. Although pore radius is <5 nm, the value of k sur /k app decreases when pressure increases. This is because that the pore radius approaches the diameter of methane molecule (0.38 nm) when pore radius gets smaller than 5 nm, so surface diffusion plays a more important role than other gas transport mechanisms.

| Effects of different gas transport mechanisms on apparent permeability of nanopores
But when pressure increases, more methane molecules occupied in pores, so bulk diffusion will exist as well as surface diffusion.
The ratios of AP of Knudsen diffusion to total AP (k Kn / k app ) under different pore radius and pressure are illustrated in Figure 5B. It can be seen that the value of k Kn / k app increases with the decrease of pore radius and pressure when the pore radius is >10 nm, while the value of k Kn / k app decreases with the decrease of pore radius when the pore radius is <10 nm. This is because surface diffusion is the main gas transport type when pore radius is <10 nm. Figure 5C presents the ratios of AP of viscous flow to total AP (k v-s /k app ) under different pore radius and pressure, and the value of k v-s /k app increases with the increase of pore radius and pressure. Figure 6 presents the equivalent permeability of different gas transport mechanisms under different pore radius and pressure. The equivalent permeability of surface diffusion increases with the decrease of pore radius and pressure. The equivalent permeability of Knudsen diffusion increases with the decrease of pressure in Figure 6B, which is in accord with Equation (40), and k Kn decreases with the decrease of pore radius according to Equations (29) and (40), but the decrease of k Kn is slowly with the decrease of pore radius because that the ratio of collisions between molecules and pore wall to total collisions (1−ω) increases with the decrease of pore radius. The equivalent permeability of viscous flow increases with the increase of pore radius and pressure, which is shown in Figure 6C. Figure 7 shows the ratio of AP of different gas transport mechanisms to total AP under different pore radius and pressure. When the pore radius is smaller than 100 nm, surface diffusion and Knudsen diffusion both play important roles, but the effect of surface diffusion on apparent permeability is much greater than that of Knudsen diffusion. When pore radius is bigger than 10 nm, viscous flow with considering slippage effect is the main gas transport type when pressure is higher than 15 MPa. It can also be seen from Figure 7 that the pressure between 1 and 15 MPa has a greater influence on apparent permeability than the pressure between 15 MPa and 30 MPa. Figure 8 shows the k sur /k app value of each grid block in shale gas reservoir after 600 days of production when the pore radius is 5, 15, and 50 nm, respectively. The k sur /k app value of each grid block in shale gas reservoir increases with the decrease of pore radius. When the pore radius of shale reservoir is smaller, the pressure of grid blocks in shale gas reservoir declines faster, which leads to a higher value of k sur /k app . Figure 8 shows the effects of different pore radius on the value of k sur /k app of each grid block in shale gas reservoir, and then, Figure 9 will further investigate the effect of surface diffusion under different pore radius on gas production of MSFHWs. Figure 9 presents the effect of surface diffusion under different pore radius on production rate and cumulative production of MSFHW. When pore radius is 5 nm, the effect of surface diffusion on production rate and cumulative production cannot be ignored. But the effect of surface diffusion on production rate and cumulative production becomes smaller when pore radius is bigger than 15 nm. And the effects of Knudsen diffusion and slip flow on production rate and cumulative production increases with the decrease of pore radius. Although the proposed study shows the importance of matrix as a significant source of gas in the total gas production, there are studies 24,51-54 that do not endorse these results in that the importance of matrix relative to fractures in total production is insignificant.

| Effect of surface diffusion coefficient on gas production of MSFHWs
The coverage of adsorbed gas on pore walls varies with Langmuir pressure according to Equation (5), and the coverage increase with the decrease of Langmuir pressure. Surface diffusion coefficient in Equation (10) is obtained under low pressure, and surface diffusion coefficient calculated by Equation (11) considers the adsorbed gas coverage under a high-pressure condition. Figure 10 presents the effect of different surface diffusion coefficient under different Langmuir pressure on production rate and cumulative production of MSFHW. When D s is calculated by Equation (10) the effect of Langmuir pressure on gas production can be ignored. But the effect of Langmuir pressure on gas production becomes significant when D s is calculated by Equation (11), and the difference between surface diffusion coefficients calculated by Equation (10) and Equation (11) becomes bigger when the Langmuir pressure gets lower, because Equation (10) does not consider the influence of pressure, which cannot describe surface diffusion under high pressure. The adsorbed gas coverage is bigger under high-pressure condition compared with low-pressure condition, which increases the surface diffusion flux according to Equations (4) and (9). And when Langmuir pressure is bigger, the adsorbed Knudsen diffusion viscous flow with considering slippage model gas coverage becomes smaller according to Equation (5), which leads to a smaller surface diffusion flux.

| Effect of stimulated reservoir volume (SRV) on gas production of MSFHWs
The effects of different stimulated reservoir volumes on gas production rate and cumulative production of MSFHW are investigated. Six different SRV schemes were proposed, and the case description of SRV schemes was given in Table 6. Figure 11 shows the schematic diagram of six different SRV schemes. Figure 12 shows the effects of different stimulated reservoir volumes on production rate and cumulative production. From Figure 12, we can find that the production rate and cumulative production increase with the increase of the length of SRV region. But the growth rates of production rate and cumulative production of MSFHW decline when the length of SRV region increases to a certain degree.
It can also be seen from Figure 13 that the performance of MSFHW is greatly enhanced once a fracture network is formed around the main hydraulic fracture. The production can be increased by 28.1% when a 30 m wide fracture network is formed around the main hydraulic fracture. The production can be increased by 44.4% when a 50 m wide fracture network is formed. The production can be increased by 55.7% when a 70 m wide fracture network is formed, which is more than half of the production of MSFHW without volume fracturing. But the growth rate of the production of MSFHW increases slowly when the width of the fracture network is greater than half of the hydraulic fracture spacing. For example, the production only increased by 14.9% when the width of the fracture network increases from 70 to 130 m. So there is no need to maximize the volume of the fracture network when volume fracturing is carried out, as long as the width of fracture network reaches half of the hydraulic fracture spacing.

| CONCLUSIONS
This paper presents a comprehensive apparent permeability model which couples the surface diffusion of adsorbed gas, slip flow considering additional flux produced by surface diffusion, and Knudsen diffusion. The AP model was combined with the multistage fractured horizontal well in a macroscale shale gas reservoir simulator, and the influencing parameters were analyzed thoroughly. The following conclusions can be drawn: 1. The AP model is validated with the experimental data and lattice Boltzmann method (LBM) simulation results, and the effects of different transport mechanisms on AP of nanopores are analyzed. When pore radius is smaller than 100 nm, surface diffusion and Knudsen diffusion both play important roles, but the effect of surface diffusion on apparent permeability is much bigger than Knudsen diffusion. When pore radius is bigger than 10 nm, viscous flow with considering slippage effect is the main gas transport type when pressure is higher than 15 MPa. 2. The effect of surface diffusion on gas production of MSFHW is investigated. When pore radius is 5 nm, the effect of surface diffusion on production rate and cumulative production cannot be ignored. But the effect of surface diffusion on production rate and cumulative production becomes smaller when pore radius is bigger than 15 nm. 3. When surface diffusion coefficient is calculated without considering the effect of high-pressure condition, the effect of Langmuir pressure on gas production can be ignored. If  F I G U R E 9 Effects of surface diffusion under different pore radii on production rate and cumulative production of MSFHW (decline curves represent gas production rate curves, increasing curves represent cumulative production curves). A, pore radius is 5 nm; B, pore radius is 15 nm; and C, pore radius is 50 nm the high-pressure condition is considered, the gas production increases with the decrease of Langmuir pressure. 4. The effect of stimulated reservoir volume on gas production of MSFHW is also discussed. The performance of MSFHW is greatly enhanced once fracture network is formed around the main hydraulic fracture, and there is no need to maximize the volume of fracture network when volume fracturing is carried out, as long as the width of fracture network reaches half of the hydraulic fracture spacing.