Nonlinear Seepage Model of Gas Transport in Multiscale Shale Gas Reservoirs and Productivity Analysis of Fractured Well

Shale is abundant in nanoscale pores, so gas flow in shales cannot be simply represented by Darcy formula anymore. It is crucial to figure out the influence of gas flow in nano/micro pores on actual productivity, which can provide basic theories for optimizing parameters and improving the gas production from engineering perspective. This paper considers the effects of slippage and diffusion in nanoscale based on Beskok-Karniadakis (BK) equation, which can be applicable for different flow regimes including continuum flow, slip flow, transition flow, and free-molecule flow. A new non-Darcy equation was developed based on the analysis of effects of high order terms of BK equation on permeability correction factor. By using the conformal transformation principle and pressure coupling method, we established the productivity formula of fractured well (infinite and limited conductivity) satisfying mass variable seepage flowing in fractures. The simulation results have been compared with field data and influencing parameters are analyzed thoroughly. It is concluded that slippage effect affects gas production of fractured well when wellbore pressure is less than 15MPa, and the effects of slippage and diffusion have greater influence on gas production of fractured well for reservoir with smaller permeability, especially when permeability is at nano-Darcy scale.


Introduction
Shale has been found to contain numerous well-developed nanoscale pores in recent researches: pore diameter of Haynesville basin in North America mainly distributes from 2 nm to 20 nm [1]; pore diameter of Mississippian basin is between 5 nm and 750 nm [2]; pore diameter of Sichuan basin in China is around 100 nm [3].Core analysis shows that the permeability of shale is in a range of 1 × 10 −9 -1 × 10 −3 m 2 [4], and flow in extremely low permeability shales undergoes a transition from Darcy regime to other regimes owing to the significant effect of collisions between molecules and pore walls on gas transport.Darcy Formula cannot describe all the flow regimes in shale gas reservoirs, so it is necessary to establish a transport equation valid for all flow regimes in multiscale shale gas reservoirs.
Gas flow in shales cannot be simply represented by Darcy equation anymore owing to pore throat diameters on the order of nanometers [4][5][6][7][8][9][10][11].Javadpour et al. [5] employed a group of Gaussian distribution curves to model nonlinear flow in a combination of nano/micropore networks and was validated with the results of canister desorption tests.Wang and Reed [6] proposed that gas flow in shale matrix can be non-Darcy as a result of the slippage effect but a Darcy type in fractures.Huang et al. [7,8] used apparent permeability which combines Knudsen diffusion coefficient and slippage factor to describe gas flow in nanoscale pores within shale matrix.Michel et al. [9] established a model to correct the permeability to account for all flow regimes including continuum, such as slip, transition, and molecular flow in shales based on Beskok and Karniadakis equation [10,11], but the slip coefficient  is assumed to be constant.Deng et al. [4] developed a seepage model in regard to different slip coefficients on the basis of Michel's model, and diffusion is considered as well.
However, the slippage factor is only assigned to some different constants in previous articles which cannot really simulate the entire flow regimes for gas flow in shale gas reservoirs, but actually slippage factor is changing with pressure, pore radius, and other parameters.This paper considers the effects of slippage and diffusion based on Beskok-Karniadakis (BK) equation, and slippage factor and diffusion coefficient are varying with pressure and shale permeability in the process of simulation.A new non-Darcy equation of motion is obtained based on the analysis of the effects of high order terms of BK equation on permeability correction factor.By employing conformal transformation principle and pressure coupling method, a second order differential equation satisfying variable mass flowing in fractures is derived, and productivity formula of fractured well (infinite and limited conductivity) for shale gas reservoir is established.The simulation results are compared with field data, and effects of slippage and diffusion on productivity of fractured well have been analyzed thoroughly as well as other parameters.

Flow Regimes in Shale Gas Reservoirs
Knudsen number Kn is defined as the ratio of the fluid mean free path  and pore throat diameter , which is a widely recognized dimensionless parameter to determine the degree of appropriateness of continuum model: Gas flow regimes can be classified into four categories [12] depending on the Knudsen number (Table 1): (1) continuum flow, (2) slip flow, (3) transition flow, (4) and free-molecule flow.In continuum flow regime, no-slip boundary condition is valid and gas flow is linear.As Kn increases, the rarefaction effects become more pronounced and the continuum assumption breaks down eventually.So in flow regimes other than continuum flow, Darcy law is not applicable.
To analyze the flow regimes in shale gas reservoirs, Figure 1 presents the Knudsen number under different pore sizes between 1 nm and 100 m and different pressure ranging from atmosphere to 100 MPa where gas flow from formation to surface undergoes.Knudsen number increases when pressure drops and pore throat diameter decreases.The horizontal plane in Figure 1 represents that Kn equals 0.001.Gas flow regime in shales is continuum flow when Knudsen number is beneath the horizontal plane, but it changes into other flow regimes when Knudsen number is above the horizontal plane.
For example, when formation pressure is 100 MPa, gas flow in shale gas reservoirs is slip flow when pore diameter is between 1 nm and 63 nm and is transition flow when pore diameter is smaller than 1 nm, but when formation pressure drops to 10 MPa, gas flow in shale gas reservoirs is slip flow when pore diameter is between 10 nm and 631 nm and is transition flow when pore diameter is smaller than 10 nm.Gas  flow in shale gas reservoirs is free-molecule flow when pore diameter is smaller than 6.3 nm and pressure is below 8 MPa.So we can conclude that gas flow in shale gas reservoirs is a multiscale flow process which includes continuum flow, slip flow, transition flow, and free-molecule flow.

Multiscale Non-Darcy Seepage Model in Shale Matrix
Volume flow rate of gas flow in reservoir can be represented by Darcy Law: According to the Beskok-Karniadakis equation which accounts for different flow regimes of the continuum flow, slip flow, transition flow, and free-molecule flow, the relationship between flow rate and pressure gradient can be described as Combining ( 2) and (3), permeability correction factor can be written as We can see from (4) that the apparent permeability is close to the value of absolute permeability when Kn approaches 0. When Knudsen number gets bigger, which  means that the flow in microtube is no longer Darcy flow, the apparent permeability should be corrected.The rarefaction coefficient is defined by Beskok and Karniadakis (1999) as To simplify Beskok-Karniadakis equation [4,9,13], permeability correction factor can be rewritten as The first two terms of (6) represent the first order term to no-slip flow in Beskok-Karniadakis equation.When Kn < 0.1, the following approximations are valid [4,9]: When Kn > 0.1, the value of rarefaction coefficient calculated by ( 7) is negative, which is very different from rarefaction coefficient calculated by ( 5), so we use another approximation instead: By combining ( 8) and ( 9), the permeability correction factor can be rewritten as Here, the influence of high order terms of (10) on permeability correction factor will be discussed.In Figures 2(a) and 2(b), the term of Kn 2.4 has little influence on permeability correction factor whether the value of slippage factor  is big or small; the effect of term Kn 2 on permeability correction factor can be neglected when the value of slippage factor is small, but this effect becomes greater when slippage factor gets bigger.For slip flow and continuum flow, the high order (>2) correction term of (10) can be neglected, and (3) can be derived as The definition of mean path is given as [14] Knudsen diffusion coefficient is mainly changing with pore diameter, which is defined as [5,12,15]: Slippage factor  is changing with pressure and pore diameter, which is given by [16][17][18]  = ( 8  ) Knudsen number can be represented by Knudsen diffusion coefficient according to (12) and ( 13): Taking ( 15) into (11), we can get According to Hagen-Poiseuille model, relationship between permeability and pore diameter can be written as Taking ( 17) into ( 16), we can derive Considering both sides of ( 18) multiplied by seepage area and divided by volume factor, we can get the volume flow rate at standard condition: where the gas volume factor   can be derived according to the equation of state:

Productivity Model of Fractured Well
The natural productivity of shale gas reservoirs is very low, so formation has to be artificially stimulated (hydraulically fractured) to render commercial production.So this paper establishes steady seepage model of fractured well for shale gas reservoirs by conformal transformation.According to the principle of conformal transformation, the gas production and pressure of the boundary remains the same, only the length of line and flow forms change [19,20].

Physical Assumptions.
To make this mathematical model more tractable, the following simplified assumptions are made.
(1) The reservoir is isotropic and gas flow in shale gas reservoir is planar (steady state).
(2) Flow in nanopores of shale matrix is assumed to be non-Darcy, and fluid in shale gas reservoirs is single phase.
(3) Fractured vertical well is produced at a constant wellbore pressure and is penetrated fully in shale gas reservoir.The height of hydraulic fractures is equal to the thickness of formation.
(4) The fractures are vertical and symmetrically located at two sides of wellbore.Fracture half-length is   ; fracture width is negligible when the vertical fractures are assumed to have infinite conductivity and is   when vertical fractures are assumed to have finite conductivity.

Fractured Well with Infinite Conductivity. Conformal transformation function is defined as
Taking  =  + ,  =  + V into (21), we can get We can see from Table 2 that the upper half plane formation in  plane (Figure 3(a)) is transformed to a half infinite formation with width of  at the right side of fractured well in  plane (Figure 3(b)); elliptical constant pressure boundary with radius of   is transformed to a linear one with length of ; fractured well with length of 2  is transformed to a line sink in  plane.
Equipotential line equation of gas flow in  plane can be obtained through (22) The boundary can be viewed as circular when the boundary is far from fractured well, so Then the equipotential line equation of boundary can be rewritten as Solving (25), we can get the relationship between the length of drainage area  0 in  plane and the radius of boundary   in  plane: Combined with (19), production rate of drainage area in  plane can be derived: Gas production of fractured well in  plane can be obtained according to (26) (28)

Fractured Well with Finite Conductivity.
Compared with hydraulic fracture with infinite conductivity, fractured vertical well with finite conductivity is more close to real condition, and the effects of parameters of hydraulic fractures on well performance can be analyzed as well.
Considering an element of hydraulic fracture in  plane in Figure 3(b), we assume that gas only flows in V direction (flows to wellbore) in hydraulic fracture (Figure 4).So flow rate in fracture V V is a constant value in  direction and a variable value in V direction; flow rate from reservoir to the element V  equals a constant value.
The mass balance equation of gas flow in hydraulic fracture can be obtained by principle of mass conservation: Considering both sides of (29) divided by ΔV, the partial difference equation can be derived when ΔV takes an infinitesimal: Pseudopressure function of gas flow in hydraulic fracture is defined as Equation (30) can be The boundary conditions of hydraulic fracture are given as follows: (33) Combining with (33), we can get the solution of (32): where the parameters  1 ,  2 , and  are defined as follows: The gas production of fractured vertical well can be written as Combined with (31), we can get the productivity formula of fractured vertical well with finite conductivity: where the first term in the bracket is the gas production calculated by Darcy formula, which is defined as   .

Results and Discussion
According to the productivity formula of fractured well with finite conductivity, the effects of permeability correction factor, slippage factor, Knudsen diffusion coefficient, matrix permeability, fracture half-length, and fracture conductivity on productivity of fractured well are analyzed by combining the data from a single well of shale gas reservoirs in China.
The data for production simulation are presented in Table 3.

Effect of Permeability Correction Factor on Productivity.
Figure 5 presents the effect of permeability correction factor on productivity of fractured well, and we can see from this figure that permeability correction factor has a great  influence on productivity.Production of fractured well is very small when our model only considers Darcy flow in nano/microscale pores ( = 1), but gas production increases a lot when matrix apparent permeability is corrected by permeability correction factor . Field data shows that gas production of a fractured vertical well is 1.2 × 10 4 m 3 when production pressure drop keeps 7 MPa and fracture halflength is 200 m [4].The calculation results of multiscale non-Darcy seepage model are much more close to the field data compared with the production rate calculated by Darcy formula and production rate calculated by productivity formula only considering first order term of permeability correction factor.

Effects of Formation and Fracture Parameters on Productivity.
The effect of matrix permeability  on productivity of fractured well is illustrated in Figure 6.Gas production of fractured well increases with the matrix permeability, and gas production increases to 1.71 times, 2.32 times, and 2.87 times when matrix permeability increases to 2 times, 3 times, and 4 times, respectively.Figure 7 reflects productivity of fractured well affected by fracture half-length.Gas production increases to 1.54 times, 2.06 times, and 2.65 times when fracture half-length increases to 3 times, 4 times, and 5 times, respectively.We can find that the increase of gas production rate affected by matrix permeability is still more than production rate affected by fracture half-length even if the increase ratio of fracture half-length is bigger than matrix permeability, which means that longer fracture half-length only can improve local seepage ability of formation adjacent to hydraulic fractures, but improvement of flow capacity in the whole formation can increase the gas production substantially.
Effect of fracture conductivity on productivity of fractured well has been analyzed under different fracture penetration ratio (  /  ).We can see from Figure 8 that gas production increases with the fracture conductivity and fracture penetration ratio.But gas production increase becomes slow when the fracture conductivity increases to a certain value, and this value becomes bigger when fracture penetration ratio increases.For examples, this value will be 2 m 2 ⋅cm when fracture penetration ratio is 0.1, but this value will increase to 4 m 2 ⋅cm when fracture penetration ratio is 0.3.This chart can help to optimize fracture conductivity under different fracture penetration ratio.

Effects of Slippage and Knudsen Diffusion on Productivity.
The relationship between wellbore pressure of fractured well with finite conductivity and gas production under different slippage factor  has been shown in Figure 9. Slippage effect has negligible influence on productivity of fractured well when wellbore pressure is high, but slippage effect begins to affect gas production when wellbore pressure is lower than 15 MPa and productivity of the fractured well increases with the slippage factor.So decreasing the wellbore pressure of fractured well properly can improve gas production when exploiting shale gas reservoirs.We can see from Figure 10 that gas production rate increases with Knudsen diffusion coefficient, which has a greater impact on gas production than slippage effect.As we can see in Figure 11, the productivity of fractured well derived in this paper (  ) has been compared with production calculated by Darcy formula (  ) under different matrix permeability and wellbore pressure.Figure 11(a) shows that bigger matrix permeability and a lower wellbore  pressure lead to a higher productivity of fractured well, and gas production   is much bigger than Darcy production rate   .The effects of Knudsen diffusion and slip flow on productivity of fractured well can be observed by taking the ratio of   and   in Figure 11(b).The effects on gas production increase when matrix permeability becomes smaller and wellbore pressure declines, which become most significant under the smallest matrix permeability (1 × 10 −5 mD) and the lowest wellbore pressure (1 MPa) in Figure 11(b).Figure 12 represents the relationship of gas production rate   and matrix permeability  at different wellbore pressure.For conventional reservoirs, production rate decreases with formation permeability; for shale gas reservoirs, gas production rate decreases with matrix permeability first but then increases with permeability owing to the effects of Knudsen diffusion and slip flow after a certain permeability value, which is called turning point in this paper.The value of turning point decreases when wellbore pressure rises, which is because Knudsen diffusion and slippage effect are more likely to occur in formation with smaller permeability and lower pressure, so a smaller turning point would be required when pressure rises.
Table 4 shows the ratio of gas production rate at minimum permeability   min in Figure 12 and production at turning point under different wellbore pressures  TP .The ratios decrease when wellbore pressure increases, which means that Knudsen diffusion and slippage effects have a greater influence on productivity of fractured well under a lower wellbore pressure.But the gas production rate   ( = 1 × 10 −5 mD,   = 1 MPa) is only 1.6 times the production   at turning point (Table 4) even though the production   is almost 150 times the Darcy production rate   in Figure 11(b).Knudsen diffusion and slippage effects can improve the productivity at extremely low permeability significantly compared with gas production calculated by Darcy formula (without considering effects of diffusion and slip flow), but the improvement of productivity is still small compared with higher permeability when Knudsen diffusion and slippage effects are considered in both situations.

Conclusions
In this paper, a multiscale comprehensive mathematical model had been established, which can simulate different flow regimes including continuum flow, slip flow, transition flow, and free-molecule flow.We deduced the productivity formula for fractured well with limited conductivity for shale gas reservoirs, and the influencing parameters were analyzed thoroughly.The following conclusions can be drawn.
(1) Flow regimes in shale gas reservoir were analyzed by calculating Knudsen number under different pressure and different pore throat diameter.Gas flow in shale gas reservoirs with nano/microscale pores is a combination of continuum flow, slip flow, transition flow, and free-molecule flow, which cannot simply be represented by Darcy formula anymore.
(2) A new non-Darcy equation was developed based on Beskok-Karniadakis equation which is applicable for all flow regimes, and effects of high order terms of BK equation on permeability correction factor were analyzed under different Knudsen number.Permeability correction factor is increasingly deviating from a Darcy type as Knudsen number increases.
(3) Productivity formula of fractured well satisfying variable mass flowing in fractures was obtained, and the simulated results showed that production rate calculated by non-Darcy seepage model with higher order terms of permeability correction factor is more close to the field production data, compared with the production rate calculated by Darcy formula and production rate calculated by productivity formula only considering first order term of permeability correction factor.
(4) Effects of slippage and diffusion which are changing with pressure and pore sizes in nanoscale pores were considered in simulation.Simulation results showed that slippage effect affects gas production of fractured well only when wellbore pressure less than 15 MPa,  and the effects of slippage and diffusion have a greater influence on gas production of reservoirs with smaller matrix permeability, especially when permeability is at nano-Darcy scale.Matrix permeability, fracture half length, fracture penetration ratio, and fracture conductivity can improve gas production at different degrees as well.

P
or e th ro at d ia m et er (n m ) Pr es su re (k Pa )

Figure 1 :
Figure 1: The relationship between Knudsen number, pressure, and pore radius.

Figure 2 :
Figure 2: Relationship between permeability correction factor and Knudsen number.

Figure 3 :
Figure 3: Sketch map of conformal transformation to fractured well.

Figure 4 :
Figure 4: Sketch map of hydraulic fracture element.

Figure 5 :
Figure 5: Effect of permeability correction factor  on productivity.

Figure 10 :
Figure 10: Effect of Knudsen diffusion coefficient  K on productivity.

2 M
at ri x p er m ea bi lit y (m D ) W e ll b o r e p r e s s u r e ( M P a ) Effects of slip flow and Knudsen diffusion Comparison between   and   (at different permeability and wellbore pressure) ix p e r m e a b il it y ( m D ) W e ll b o re p re ss u re (M P a The value of   /  at different matrix permeability and wellbore pressure

Figure 11 :Figure 12 :
Figure 11: Effects of Knudsen diffusion coefficient  K and slip factor  on productivity.

Latin𝑏 :
Slippage coefficient   : Gas volume factor  K : Knudsen diffusion coefficient, m 2 /s : Fraction of molecules striking pore wall which are diffusely reflected ℎ: Thickness of shale gas reservoir, m : Apparent permeability Kn: Knudsen number   : Permeability of hydraulic fracture  0 : Absolute permeability   : Fracture half length : Molecular mass, kg/mol : Pseudopressure   : Initial pseudopressure   : Pseudopressure of fractured well : Pressure, Pa   : Formation pressure, Pa   : Wellbore pressure, Pa  avg : Average pressure, Pa  sc : Pressure at standard condition, Pa   : Volume flow rate calculated by Darcy formula, m 3 /s

Table 2 :
Conformal transformations of the fractured well.

Table 3 :
Data of shale gas reservoir for production simulation.

Table 4 :
Comparison of production at turning point and production at minimum permeability under different wellbore pressure (the gas production rate unit: ×10 4 m 3 /d).