Numerical Simulation of Non-Darcy Flow in Naturally Fractured Tight Gas Reservoirs for Enhanced Gas Recovery

1. Introduction

Fluid flow in porous media is present in several areas of scientific and technological knowledge [1]. In this context, due to technological advances, so-called unconventional reservoirs have aroused great interest in the oil industry. However, even with the most recent knowledge and technologies, the economically viable production of these reservoirs still presents challenges when aiming to increase gas recovery rates [2]. On the other hand, they have great potential regarding the geological storage of carbon dioxide. Therefore, the numerical simulation of flow in tight and shale gas reservoirs is a topic of great importance for the energy sector on a global scale.

Specifically, it is relevant to mention unconventional reservoirs of low permeability, tight gas, or those with very low permeability, that is, shale gas. Shale is a fine-grained sedimentary rock typically composed of clay, silt, and organic matter. Unlike conventional reservoirs, shale gas reservoirs have low porosity and permeability, which makes commercial exploitation of the gas more difficult. For example, sandstone-type formations have low permeability.

Usually, to extract gas from shales, a process called hydraulic fracturing (fracking) is used. It involves injecting a mixture of water, sand, and chemicals into the formation to create fractures, which allow gas to flow more easily. This extraction method differs from traditional natural gas production methods, such as drilling a well and allowing the gas to flow to the surface, which may not be effective for shale-type gas reservoirs [3].

1.1. Brief Review: Shale and Tight Gas Reservoirs

In recent years, investigations into low and very-low permeability reservoirs containing natural gas have intensified. Among the topics explored, it is possible to highlight the studies: flow mechanisms in shale gas reservoirs [4,5]; correlations for absolute permeability as a function of pressure [6,7,8]; fractured horizontal wells in shale gas reservoirs [9]; about naturally fractured reservoirs [10]; two-phase flow including gas slippage and production through fractured horizontal wells [11,12,13]; and the effects of heterogeneity and non-Darcy in tight gas reservoirs [14].

Guria [15] addressed the effects of non-ideal gas behavior and slippage when measuring permeability in porous media. The author used cubic equations of state as the van der Waals, Soave-Redlich-Kwong, and Peng-Robinson to reproduce the non-ideal behavior of gases, including or not hydrocarbons. He also proposed a comprehensive mathematical model for evaluating the apparent permeability, dependent on pressure and temperature, in porous media when real gases flow.

Rubin et al. [16] performed experiments on shales of the Marcellus formation. They showed that the effects of gas slippage and matrix compaction are significant when evaluating gas production due to the substantial depletion of reservoir pressure, especially at the end of production. However, their impact on gas recovery quantification and hydraulic fracturing design has not yet been understood and systematically investigated. According to the authors, the results showed that ignoring both effects when simulating flow in the reservoir leads to an incorrect estimate of gas production.

In the work by Ding et al. [17], they carried out experiments to understand the nonlinear behavior of flow in tight gas reservoirs, including the effects of slippage and inertial/turbulent effects. The authors also attested that experimental methods applicable to conventional reservoirs lead to significant errors when applied to these unconventional reservoirs.

Friedel and Voigt [18] applied numerical reservoir simulation to study the inertial effect in a fractured reservoir. They conducted simulations considering a fully implicit formulation and investigated synthetic cases for gas production. The results indicated that non-Darcy flow influenced productivity, even for relatively low gas production rates. The authors also mention that they evaluated a scenario in which they verified a reduction between 21% and 40% in gas production.

In 2012, Clarkson et al. [19] discussed the rapid evolution of the technology used to evaluate the properties of unconventional gas reservoirs and hydraulic fracturing. The authors proposed a model for optimizing the development of unconventional gas reservoir fields, raising critical issues related to the analysis of reservoir samples, production data analysis, and transient flow responses. The authors aimed, for example, to review existing methods for reservoir evaluation and to introduce new methods for evaluating heterogeneities, focusing on tight and shale reservoirs and discussing the influence of non-Darcy inertial effects.

On the other hand, Ye et al. [20] considered non-Darcy inertial effects using Forchheimer’s model and unconventional coal-bed methane reservoirs. According to experimental results, the so-called non-Darcy flow factor, $\beta$, associated with the inertial effects, can be expressed as a power law dependent on absolute permeability. The authors point out that for certain conventional gas reservoirs, $\beta$ can be considered constant, but this may not be suitable for coal-bed methane with a notable change in permeability. They also report that when the paper was published, few studies considered $\beta$ as a function of coal permeability in reservoir simulations. The results of the simulations showed that considering a constant $\beta$ factor can significantly underestimate or overestimate the gas production rate for coal-bed methane wells.

Other authors have also pointed to the possibility of errors in calculating the productivity of gas wells if the inertial effects are not considered, including cases with and without fractures and the use of horizontal wells. [21] dealt with this issue, proposing the concept of the apparent equivalent radius of the well to consider the non-Darcy inertial flow in a fracture. Firstly, the authors presented an iterative process to calculate the productivity of vertical fractured gas wells using the proposed concept, followed by an analysis of the productivity of horizontal multi-fractured gas-producing wells. The results showed that non-Darcy inertial effects can significantly decrease production even in low-permeability reservoirs (less than 1 mD) and reduce the potential interference between hydraulic fractures.

According to Farahani et al. [22], two main effects for accurately modeling the productivity of gas reservoirs are slippage and inertial flow in porous media. These effects lead to an introduction of an apparent permeability for the gas phase, different from the absolute permeability measured when we have a single liquid phase flow. Therefore, we must consider these elements when performing pressure tests to characterize gas reservoirs. The authors studied gas flow conditions in laboratory measurements using Iranian Kangan and Dalan carbonates, developing a correlation to estimate the inertia factor of the samples. In addition, they introduced an equation to estimate the maximum allowable pressure variation per length to maintain the laminar flow regime (Darcy flow) in the samples.

Furthermore, Fu et al. [12] proposed a new productivity model for multi-fractured horizontal wells and considered the beta inertial factor to incorporate the effects of high velocity. In addition, the authors considered the consequences of stress reduction (negative on effective permeability) and slip (positive on absolute permeability). The authors studied tight gas reservoirs, including the impacts of the presence of the water phase in the porous medium.

Assessing the economic viability of unconventional gas reservoirs has been challenging due to the inherent difficulties in accurately quantifying absolute permeabilities on the micro/nano-darcy scale [23]. It greatly influences the behavior of gas mass transport in porous media with low permeability. However, correctly measuring the permeability of porous media of this type is difficult. Gas flow differs from what occurs with liquids, for example, due to its high compressibility, the effect of slippage, and, sometimes, adsorption [24].

The study by Luo et al. [25] aimed to analyze temperature variation and quantitatively diagnose water production in multi-fractured horizontal wells for two-phase flows in tight gas reservoirs. The work led to a model for determining temperature considering effects such as, for example, Joule-Thomson and thermal expansion. They simulated synthetic cases to illustrate temperature behavior. The sensitivity analysis indicated that the gas phase slippage was responsible for reducing the well temperature. They also concluded that it was significantly impacted depending on the reservoir permeability values and the length of the fractures. Subsequently, they applied the proposed model to a field case, and the comparison showed a good agreement between the results.

Liu et al. [26] experimentally examined the effect of the pore-throat structure on the behavior of gas-water flow in a formation composed of low-permeability sandstone bearing natural gas. More specifically, they carried out mercury injection tests under controlled pressure to measure capillary pressure and identify the characteristics of a pore-throat structure and its connectivity and distribution of pore diameters. Subsequently, the specificities of the flow, including the variation in the saturation of the mobile fluid, the relative permeability, and the degree of water blockage, were quantified using displacement experiments and the nuclear magnetic resonance technique. The pore-throat structure, especially its connectivity, was estimated to be the dominant factor in the distribution of mobile fluid saturation and the infiltration characteristics of the two phases in a tight gas reservoir composed of sandstone. This study allowed us to understand how the structure of the pore throat influences the behavior of gas-water flow in a low-permeability gas reservoir, emphasizing the importance of identifying and classifying flow patterns.

The gas reservoir in the Kelasu field, in the Tarim basin (China), is of the tight gas type, ultra-deep, with existing fractures at various scales within the matrix, including faults. The flow inside it was studied by Sun et al. [27] via numerical simulations of the pressure test type in vertical wells, considering the porous matrix, fractures, and faults, combining the random generation of natural fracture networks with the use of discrete unstructured fractures. The numerical methodology used the mixed finite element method, and they obtained typical curves using different random fracture networks. Based on the observed results, they classified the distribution of the fracture network of fractured low-permeability sandstone reservoirs. Furthermore, they also discussed the influence of the random generation of fracture networks on well tests.

1.2. Carbon Sequestration

Therefore, as we can see, tight gas and shale gas reservoirs have become increasingly important within unconventional reserves. Among the applications involving them, we can cite carbon dioxide injection into the reservoir, aiming at carbon sequestration and natural gas recovery [28]. For example, due to the competitive adsorption between carbon dioxide and methane, the capture and carbon storage in tight and shale gas reservoirs also constitutes an opportunity to improve gas recovery. Furthermore, we already use the injection of CO₂ in oil recovery.

To address the impacts of greenhouse gas emissions into the atmosphere, one of the strategies to mitigate the problem of carbon dioxide release is its capture and sequestration in underground rock formations, which can have other purposes besides storage [29]. One example would be the consideration of low-permeability reservoirs, as the controlled injection of CO₂ into these formations improves natural gas recovery and proves to be an efficient solution for storing carbon dioxide [30].

According to Liu et al. [31], CO₂ storage to enhance natural gas recovery (CO₂ Storage with Enhanced Gas Recovery, CSEGR) is a promising option. It can reduce the concentration of carbon dioxide in the environment by sequestering it in gas reservoirs and simultaneously increase the production of natural gas, which can be used, for example, in sustainable hydrogen production through methane reforming [32].

Aminu and Manovic [33] expect that carbon capture and storage will have a role in reducing greenhouse gas emission. For example, in the southern North Sea, the Bunter sandstone formation has been identified as a potential reservoir that can store significant amounts of carbon dioxide. The formation has favorable rock and geometry properties, making CO₂ storage feasible on an industrial scale. However, they noted that captured CO₂ typically contains impurities, which can alter the boundaries of the carbon dioxide phase diagram, resulting in higher costs for storage operations. The authors also modeled the effect of CO₂ and impurities (NO₂, SO₂, and H₂S) on reservoir production performance.

From the point of view of the presence of heterogeneities, Wang et al. [34] report that natural or induced fractures are typically present in sub-surface geological formations. They, therefore, need to be studied for a reliable estimate of long-term carbon dioxide storage. Fractures can undermine storage security as they increase the risk of carbon dioxide escaping. In addition, fractures can act as flow barriers, causing significant pressure gradients in relatively small regions near fractures. However, despite their high sensitivities, the impact of fractures on the complete storage process is still not understood. In the cited study, the authors developed and applied a numerical model to analyze the role of discrete fractures in the flow mechanism and transport of CO₂ plumes in simple and complex fracture geometries. The results indicate that the fracture has effects that depend on the different trapping mechanisms.

Therefore, carbon capture and storage can contribute to climate change mitigation if implemented on a large scale. Achieving high injection rates in deep saline aquifers, for example, requires a detailed assessment of injectivity performance and an evaluation of the processes that alter rock properties in the region of nearby wells, as discussed by Parvin et al. [35]. One of the most common forms of injectivity loss in the context of carbon dioxide storage in saline aquifers is salt precipitation driven by brine evaporation into the relatively dry flow of injected CO₂, according to the authors.

2. Non-Darcy Flow

For flows inside a porous medium corresponding to high Reynolds numbers, the relationship between the flow rate and the pressure gradient [36] is no longer linear. The nonlinearity, in this case, comes from inertial effects. Although it may originate from other causes, in petroleum engineering it is essentially related to flow values. According to Amao [37], different behavior from that predicted by Darcy’s law due to inertial effects appears mainly in the following cases: near the wellbore; hydraulically fractured wells; gas reservoirs; condensate reservoirs; high flow potential wells; naturally fractured reservoirs; and gravel packs. Therefore, in these cases it would be recommended that specific flow models be used, taking into account inertial effects.

Barree and Conway [38] proposed a model, which introduces the notion of apparent permeability, as a way of including inertial effects [39,40],

$$\mathbf{v}=-{\displaystyle \frac{{\mathbf{k}}_{app}}{\mu }}\left(\nabla p-\rho g\nabla z\right),$$

(1)

where $\mathbf{v}$ is the apparent velocity, $\mu$ is the fluid viscosity, p is the pressure, $\rho$ is the density, g is the magnitude of the acceleration due to gravity, z is the depth, and the apparent permeability tensor is given by

$${\mathbf{k}}_{app}={\mathbf{k}}_{min}+{\displaystyle \frac{\mathbf{k}-{\mathbf{k}}_{min}}{{\left[1+{\left({\displaystyle \frac{\rho \left|\mathbf{v}\right|}{\mu {\tau }_c}}\right)}^F\right]}^E}},$$

(2)

and ${\mathbf{k}}_{min}$ is a minimum permeability value at high flow rates, $\mathbf{k}$ is the permeability tensor, F and E are dimensionless constants and ${\tau }_c$ is the inverse of characteristic length. It is an advantageous alternative to the Forchheimer model [36].

Usually we take $F=E=1$ and introduce ${\mathbf{k}}_{mr}={\mathbf{k}}_{min}{\mathbf{k}}^{-1}$ [38], to obtain the form used here

$${\mathbf{k}}_{app}={\mathbf{k}}_{BC}\mathbf{k}=\left({\mathbf{k}}_{mr}+{\displaystyle \frac{\mathbf{I}-{\mathbf{k}}_{mr}}{1+Re}}\right)\mathbf{k},$$

(3)

where

$$Re={\displaystyle \frac{\rho \left|\mathbf{v}\right|}{\mu {\tau }_c}}$$

(4)

is the Reynolds number, and $\mathbf{I}$ is the identity tensor.

As already mentioned, this model has some advantages, such as providing a limit value for apparent permeability at high flows [39]. Furthermore, it would allow a better representation of transitions between flow regimes in the porous medium, according to Lai et al. [39].

In this work, we also consider a simplified model to account for the effects of effective stress on permeability variation [41]

$$\mathbf{k}=exp[-\gamma (p_0-p)]{\mathbf{k}}_0,$$

(5)

where $\gamma$ is the permeability modulus, and $p_0$ and ${\mathbf{k}}_0$ are the initial pressure and absolute permeability tensor.

Furthermore, we know that the slippage effects can be introduced in terms of the Knudsen number [42,43]. Therefore, the combined effects of gas slippage and effective stress on permeability variation are determined in the form

$$\mathbf{k}=exp[-\gamma (p_0-p)]\left(1+{\displaystyle \frac{4{\overline{K}}_n}{1+{\overline{K}}_n}}\right){\mathbf{k}}_0,$$

(6)

where ${\overline{K}}_n$ is the Knudsen number of the mixture and, according to Wang et al. [44], for a two-component flow

$${\overline{K}}_n={\displaystyle \frac{1}{r_h}}\left(x_1{\lambda }_1+x_2{\lambda }_2\right),$$

(7)

where $r_h$ is the hydraulic radius and, for each component, ${\lambda }_1$ and ${\lambda }_2$ are the mean free path of gas molecules [45], and $x_1$ (CO₂) and $x_2$ (CH₄) represent the molar fraction, respectively.

3. Two-Component Flow

In this article, we are interested in simulating the enhanced recovery process when injecting carbon dioxide into a naturally fractured tight gas reservoir (in general, with matrix permeability less than 0.1 mD and porosity less than 10%) containing methane. When developing our numerical code, we considered the following hypotheses: three-dimensional one-phase isothermal flow; newtonian fluids; compressible fluids; small compressibility rock [40]; no adsorption effects; no chemical reactions; and fractured media modeled using Discrete Fracture Networks (DFN) [46].

The governing equations for two-component single-phase flow can be obtained from the continuity equation [47] and modified Darcy’s law (1) [39,40], since we assume that the flow is isothermal. Therefore, after combining these two equations we obtain [48,49]

$${\displaystyle \frac{\partial }{\partial t}}\left(\varphi \xi \right)=\nabla ·\left[{\displaystyle \frac{\xi }{\mu }}{\mathbf{k}}_{app}\left(\nabla p-\rho g\nabla z\right)\right]+q,$$

(8)

where $\varphi$ represents the porosity, $\xi$ is the molar density, given by the real gas law

$$\xi ={\displaystyle \frac{p}{RTZ}},$$

(9)

q is a source term, R is the universal gas constant, T the temperature and Z is the compressibility factor.

On the other hand, when it comes to molar fractions, the governing equation is given by [48,49],

$${\displaystyle \frac{\partial }{\partial t}}\left(\varphi x_m\xi \right)=\nabla ·\left[{\displaystyle \frac{x_m\xi }{\mu }}{\mathbf{k}}_{app}\left(\nabla p-\rho g\nabla z\right)+\varphi \xi \mathbf{D}\nabla x_m\right]+q_m$$

(10)

where $m=1,2,\cdots ,N_c$, $x_m$ is the molar fraction, $q_m$ is a source term, $N_c$ is the number of components and $\mathbf{D}$ is the effective dispersion tensor, which can be introduced as [49,50]:

$$\mathbf{D}\left(\mathbf{v}\right)=\tau D_m\mathbf{I}+D_{\perp }\left|\mathbf{v}\right|\mathbf{I}+\left(D_{\Vert }-D_{\perp }\right){\displaystyle \frac{\mathbf{v}\otimes \mathbf{v}}{\left|\mathbf{v}\right|}},$$

(11)

where $\tau$ is the tortuosity, $D_m$ is the molecular diffusion coefficient, $D_{\perp }$ and $D_{\Vert }$ are the transversal and longitudinal hydrodynamic dispersion tensor coefficients, and $\mathbf{v}$ is given by Equation (1). In our simulations we will consider that [50]

$$D_l=\zeta {\left({log}_{10}{L}_c\right)}^{\eta },$$

(12)

$$D_t=0.1D_l,$$

(13)

where $\zeta$ and $\eta$ are parameters that must be provided, and $L_c$ represents a characteristic length.

The dependent variables, the pressure and mole fraction of one of two components, are obtained by solving the partial differential Equations (8) and (10) numerically.

The Peng-Robison cubic equation of state is employed to calculate the compressibility factor Z [49,50]. For this purpose, we start from its original form

$$p={\displaystyle \frac{RT}{v-b}}-{\displaystyle \frac{a}{v^2+2bv-b^2}}$$

(14)

where v is the molar volume and for each component

$$a_m=0.45724\left({\displaystyle \frac{R^2{T}_{cm}^2}{p_{cm}}}\right)\left\{1+{\kappa }_m{\left[1-{\left({\displaystyle \frac{T}{T_{cm}}}\right)}^{1/2}\right]}^2\right\},$$

(15)

$$b_m=0.07780\left({\displaystyle \frac{R{T}_{cm}}{p_{cm}}}\right),$$

(16)

$${\kappa }_m=0.37464+1.54226{\omega }_m-0.26992{\omega }_m^2,$$

(17)

where ${\omega }_m$ is the acentric factor of component, and $p_{cm}$ and $T_{cm}$ are the respective critical pressure and temperature.

For a mixture of $N_c$ components [51]

$$a=\sum _{m=1}^{N_c}\sum _{n=1}^{N_c}{x}_m{x}_n(1-k_{mn})\sqrt{a_m{a}_n}$$

(18)

$$b=\sum _{m=1}^{N_c}{x}_m{b}_m$$

(19)

where $k_{mn}$ is a binary interaction parameter, considered to be null in our simulations.

Thus, once known the values of the universal gas constant (R), pressure (p) and temperature (T), we can determine the compressibility factor (Z) from the Peng-Robison equation (14) rewritten as

$$Z^3-\left(1-B\right)Z^2+\left(A-2B-3B^2\right)Z-\left(AB-B^2-B^3\right)=0,$$

(20)

where

$$A={\displaystyle \frac{ap}{R^2{T}^2}},B={\displaystyle \frac{bp}{RT}}$$

(21)

where a and b are calculated using Equations (18) and (19).

Finally, the molar density value can be calculated using Equation (9) and the viscosity using the correlation proposed by Lohrenz et al. [52].

4. Finite Volume Formulation

We used the Finite Volume Method (FVM) [53] in the numerical resolution of the governing equations, since their respective theoretical solutions are not known.

As is usual in the case of this method, a finite number of control volumes are superimposed on the domain and the governing equations are discretized by integrating them in time (from $t^n$ to $t^{n+1}$) and space, over a finite volume [53].

To obtain the discrete equations, we use time-implicit formulations, approximations such as three points centered differences for the spatial derivatives, with the exception of advective terms, which are approximated through a first-order upwind type scheme [49].

Furthermore, we employ a structured three-dimensional mesh when partitioning the domain and the planar fractures are discretized and inserted into it.

By way of illustration, Figure 1 contains the representation of a three-dimensional structured mesh, where the nodes are positioned in the centers of finite volumes and we highlight the one belonging to the central finite volume as well as those of its six neighbors.

We also adopted the established nomenclature, in which the integer subscripts i, j and k represent the indices of the center of one of the cells (or finite volumes), according to their position along the directions of the x-, y- and z-axes respectively. Besides that, the subscripts $i\pm 1/2$ ,$j\pm 1/2$, and $k\pm 1/2$ indicate the faces of the finite volumes in the direction of the Cartesian x-, y- and z-axes.

The steps relating to the integration of the governing equations will be omitted, since the usual procedures inherent to the method were employed. More details can be found in the book authored by Versteeg and Malalasekera [53]. Therefore, considering that the apparent permeability (${\mathbf{k}}_{app}$) and effective dispersion ($\mathbf{D}$) tensors are diagonal [47], the discrete final form of Equation (8) is [48,49]

$$\begin{array}{cc}\hfill & {\mathbb{T}}_{x,i+1/2,j,k}^{n+1}\left(p_{i+1,j,k}^{n+1}-p_{i,j,k}^{n+1}\right)-{\mathbb{T}}_{x,i-1/2,j,k}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i-1,j,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\mathbb{T}}_{y,i,j+1/2,k}^{n+1}\left(p_{i,j+1,k}^{n+1}-p_{i,j,k}^{n+1}\right)-{\mathbb{T}}_{y,i,j-1/2,k}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i,j-1,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\mathbb{T}}_{z,i,j,k+1/2}^{n+1}\left(p_{i,j,k+1}^{n+1}-p_{i,j,k}^{n+1}\right)-{\mathbb{T}}_{z,i,j,k-1/2}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i,j,k-1}^{n+1}\right)\hfill \\ \\ \hfill -& {\left(\rho g\mathbb{T}\right)}_{x,i+1/2,j,k}^{n+1}\left(z_{i+1,j,k}-z_{i,j,k}\right)+{\left(\rho g\mathbb{T}\right)}_{x,i-1/2,j,k}^{n+1}\left(z_{i,j,k}-z_{i-1,j,k}\right)\hfill \\ \\ \hfill -& {\left(\rho g\mathbb{T}\right)}_{y,i,j+1/2,k}^{n+1}\left(z_{i,j+1,k}-z_{i,j,k}\right)+{\left(\rho g\mathbb{T}\right)}_{y,i,j-1/2,k}^{n+1}\left(z_{i,j,k}-z_{i,j-1,k}\right)\hfill \\ \\ \hfill -& {\left(\rho g\mathbb{T}\right)}_{z,i,j,k+1/2}^{n+1}\left(z_{i,j,k+1}-z_{i,j,k}\right)+{\left(\rho g\mathbb{T}\right)}_{z,i,j,k-1/2}^{n+1}\left(z_{i,j,k}-z_{i,j,k-1}\right)+Q_{i,j,k}^{n+1}\hfill \\ \\ \hfill =& {\left[V_{b}c\left(p^{n+1}\right){\displaystyle \frac{p^{n+1}-p^n}{\Delta t}}\right]}_{i,j,k},\hfill \end{array}$$

(22)

where $p^{n+1}$ must be calculated numerically, knowing that

$$c\left(p^{n+1}\right)={\varphi }^n{\displaystyle \frac{\partial \xi }{\partial p}}+{\xi }^{n+1}{\varphi }^0{c}_{\varphi },$$

(23)

$${\displaystyle \frac{\partial \xi }{\partial p}}={\displaystyle \frac{1}{RTZ}}-{\displaystyle \frac{p}{RT{Z}^2}}{\displaystyle \frac{\partial Z}{\partial p}},$$

(24)

where Equation (23) was obtained from a conservative expansion [54], the last partial derivative is calculated from Equation (20), and we consider a constant rock compressibility.

On the other hand, to determine the molar fraction, from Equation (10) we obtain

$$\begin{array}{cc}\hfill & {\mathbb{D}}_{x,i+1/2,j,k}^{n+1}\left(x_{m,i+1,j,k}^{n+1}-x_{m,i,j,k}^{n+1}\right)-{\mathbb{D}}_{x,i-1/2,j,k}^{n+1}\left(x_{m,i,j,k}^{n+1}-x_{m,i-1,j,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\mathbb{D}}_{y,i,j+1/2,k}^{n+1}\left(x_{m,i,j+1,k}^{n+1}-x_{m,i,j,k}^{n+1}\right)-{\mathbb{D}}_{y,i,j-1/2,k}^{n+1}\left(x_{m,i,j,k}^{n+1}-x_{m,i,j-1,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\mathbb{D}}_{z,i,j,k+1/2}^{n+1}\left(x_{m,i,j,k+1}^{n+1}-x_{m,i,j,k}^{n+1}\right)-{\mathbb{D}}_{z,i,j,k-1/2}^{n+1}\left(x_{m,i,j,k}^{n+1}-x_{m,i,j,k-1}^{n+1}\right)\hfill \\ \\ \hfill +& {\left({\mathbb{T}}_m\right)}_{x,i+1/2,j,k}^{n+1}\left(p_{i+1,j,k}^{n+1}-p_{i,j,k}^{n+1}\right)-{\left({\mathbb{T}}_m\right)}_{x,i-1/2,j,k}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i-1,j,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\left({\mathbb{T}}_m\right)}_{y,i,j+1/2,k}^{n+1}\left(p_{i,j+1,k}^{n+1}-p_{i,j,k}^{n+1}\right)-{\left({\mathbb{T}}_m\right)}_{y,i,j-1/2,k}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i,j-1,k}^{n+1}\right)\hfill \\ \\ \hfill +& {\left({\mathbb{T}}_m\right)}_{z,i,j,k+1/2}^{n+1}\left(p_{i,j,k+1}^{n+1}-p_{i,j,k}^{n+1}\right)-{\left({\mathbb{T}}_m\right)}_{z,i,j-1/2,k}^{n+1}\left(p_{i,j,k}^{n+1}-p_{i,j,k-1}^{n+1}\right)\hfill \\ \\ \hfill -& {\left(\rho g{\mathbb{T}}_m\right)}_{x,i+1/2,j,k}^{n+1}\left(z_{i+1,j,k}-z_{i,j,k}\right)+{\left(\rho g{\mathbb{T}}_m\right)}_{x,i-1/2,j,k}^{n+1}\left(z_{i,j,k}-z_{i-1,j,k}\right)\hfill \\ \\ \hfill -& {\left(\rho g{\mathbb{T}}_m\right)}_{y,i,j+1/2,k}^{n+1}\left(z_{i,j+1,k}-z_{i,j,k}\right)+{\left(\rho g{\mathbb{T}}_m\right)}_{y,i,j-1/2,k}^{n+1}\left(z_{i,j,k}-z_{i,j-1,k}\right)\hfill \\ \\ \hfill -& {\left(\rho g{\mathbb{T}}_m\right)}_{z,i,j,k+1/2}^{n+1}\left(z_{i,j,k+1}-z_{i,j,k}\right)+{\left(\rho g{\mathbb{T}}_m\right)}_{z,i,j,k-1/2}^{n+1}\left(z_{i,j,k}-z_{i,j,k-1}\right)+Q_{m,i,j,k}^{n+1}\hfill \\ \\ \hfill =& {\left[V_b{\displaystyle \frac{{\left(\varphi x_m\xi \right)}^{n+1}-{\left(\varphi x_m\xi \right)}^n}{\Delta t}}\right]}_{i,j,k},\hfill \end{array}$$

(25)

where $m=1,2,\cdots ,N_c-1,$. The new molar fraction $x_m^{n+1}$ must be determined. Here, $Q_{m,i,j,k}=V_{b_{i,j,k}}{q}_{m,i,j,k}$, $Q_{i,j,k}=V_{b_{i,j,k}}{q}_{i,j,k}$, $V_{b_{i,j,k}}={\left(\Delta x\Delta y\Delta z\right)}_{i,j,k}$, where $\Delta x_{i,j,k}$, $\Delta y_{i,j,k}$ and $\Delta z_{i,j,k}$ are the grid sizes in the x-, y- and z-directions.

In these equations, we define the transmissibilities $\mathbb{T}$ by

$${\mathbb{T}}_{x,i\pm \frac{1}{2},j,k}^{n+1}\equiv {\left({\displaystyle \frac{\xi k_{ap{p}_{xx}}{A}_x}{\mu \Delta x}}\right)}_{i\pm \frac{1}{2},j,k}^{n+1},$$

(26)

$${\mathbb{T}}_{y,i,j\pm \frac{1}{2},k}^{n+1}\equiv {\left({\displaystyle \frac{\xi k_{ap{p}_{yy}}{A}_y}{\mu \Delta y}}\right)}_{i,j\pm \frac{1}{2},k}^{n+1},$$

(27)

$${\mathbb{T}}_{z,i,j,k\pm \frac{1}{2}}^{n+1}\equiv {\left({\displaystyle \frac{\xi k_{ap{p}_{zz}}{A}_z}{\mu \Delta z}}\right)}_{i,j,k\pm \frac{1}{2}}^{n+1},$$

(28)

and ${\mathbb{T}}_m=x_m\mathbb{T}$.

We also introduce the discrete dispersion coefficients $\mathbb{D}$ as

$${\mathbb{D}}_{x,i\pm 1/2,j,k}\equiv {\left({\displaystyle \frac{A_x\varphi \xi D_{xx}}{\Delta x}}\right)}_{i\pm 1/2,j,k},$$

(29)

$${\mathbb{D}}_{y,i,j\pm 1/2,k}\equiv {\left({\displaystyle \frac{A_y\varphi \xi D_{yy}}{\Delta y}}\right)}_{i,j\pm 1/2,k},$$

(30)

$${\mathbb{D}}_{z,i,j,k\pm 1/2}\equiv {\left({\displaystyle \frac{A_z\varphi \xi D_{zz}}{\Delta z}}\right)}_{i,j,k\pm 1/2},$$

(31)

for the three spatial directions, respectively.

4.1. Auxiliary Conditions

In terms of boundary conditions, the reservoir is sealed (zero flow across its boundaries) and gas is injected, at a prescribed flow rate ($Q_{sc}$), in the lower left corner and produced in the upper right at the same flow rate ($-Q_{sc}$).

The reference pressure ($p_0$) and reference temperature ($T_0$) values are provided at the initial time and the reservoir contains only methane gas at the beginning of the simulations.

4.2. Mesh Refinement

When simulating discrete fractures, depending on their dimensions and properties, we must prevent the appearance of discontinuities through the use of mesh refinement, introducing a transition region in their neighborhood [40,55]. For this purpose, the spatial dimensions of the cells neighboring the fractures will vary according to a logarithmic function and they will increase in size as they distance themselves from the fractures [49].

4.3. Solving Systems of Algebraic Equations

Because the discretized governing equations are nonlinear, we cannot directly employ algebraic methods formulated specifically for solving systems of linear algebraic equations. For us to use them, we must first linearize these equations. As a result, we will have a coupled system of equations for determining pressure and mole fraction. Aiming to decouple them [49], we employ an operator-splitting method [56].

Once the system of equations was decoupled, we chose two well-known and efficient methods to obtain their numerical solutions. As the coefficient matrix is symmetric in the case of the algebraic system whose dependent variable is pressure, we decided to use the Conjugate Gradient Method (CG) [57]. As this is not the case when determining the mole fraction, we chose to use the Biconjugate Gradient Stabilized Method (BiCGSTAB) [57].

4.4. Numerical Validation

Concluding this section, we must highlight that the numerical results obtained with this simulator were previously quantitatively validated in the case of flow in a square domain containing two orthogonal fractures that intersect each other, in the region to the right of the reservoir, forming an inverted “L”. In this test case, we considered the tracer transport in a single-phase flow of an incompressible fluid in a two-dimensional domain in the presence of these two fractures. Additional information about the domain geometry, dimension and positioning of fractures, in addition to physical properties, can be found in Debossam et al. [49].

5. Results and Discussion

The main objective of this work is to quantify the influence of inertial, gas slippage and effective stress effects on production in naturally fractured tight gas reservoirs. In advanced gas recovery techniques, we know that a fluid can be injected into the reservoir to extract a greater amount of gas. As the environmental issue is currently of great importance, we decided to inject carbon dioxide (CO₂) into a reservoir containing methane (CH₄). Therefore, in addition to increasing the volume of gas produced, we can use the reservoir to store this fluid, which is harmful to the environment.

As we are considering two-component single-phase flow in a naturally fractured reservoir, we need a model to represent the fracture network within the porous medium. Among the existing possibilities, for structured meshes, we opted for the arrangement known as Sugar Cube [58]. It is based on Barenblatt et al. [59]’s double porosity model, which predicts a system of fractures with high permeability, a porous matrix with low permeability and the existence of flow from the porous matrix to the fractures.

As previously stated, fractures are represented using the Discrete Fracture Networks (DFN) [46] technique. In it, fractures are inserted into the computational mesh and represented by a set of finite volumes. They are differentiated from the porous matrix, in addition to their dimensions, by their different physical properties, that is, porosity and absolute permeability. Furthermore, due to the discrepancy between the thickness of the fracture and the dimensions of the finite volumes, a mesh refinement is introduced in the regions around the fractures [40].

With regard to the dimensions of the reservoir, as well as other properties and parameters common to all simulations, they are presented in Table 1. We remind that the subscripts r refer to the porous matrix while f refers to the fractures.

Next, the specific properties of the injected and resident fluids, respectively carbon dioxide and methane, can be seen in Table 2.

Finally, we present in Table 3 the length of the fractures ($L_f$), their thickness (w), permeability ($k_f$) and porosity (${\varphi }_f$). In the Sugar Cube configuration, the fractures are positioned in the vertical and horizontal directions, equidistant and orthogonal to each other, and have the same length and thickness. In each $xy$-plane, we positioned nine fractures, parallel and equidistant from each other, along the longitudinal and transverse directions.

5.1. Mesh Refinement

From the results of previous simulations for a mesh refinement study [49], we found that numerical convergence is achieved and numerical diffusion is minimized to acceptable levels when we employ a mesh containing 320 × 320 ×3 finite volumes. Therefore, all simulations are performed considering this mesh.

5.2. Numerical Results

In this work, we chose to use a sealed reservoir and flow rate conditions imposed for the gases that are being injected and produced. Therefore, production curves will not be presented, but rather the advance front of CO₂ inside the reservoir, in pre-established plans. As a result of the flow rate being constant, the inertial effects, gas slippage and effective stress will be felt in the variation of the molar fraction field inside the reservoir. Subsequently, in a future work, we intend to introduce a well-reservoir coupling model [60] and work with the pressure prescribed in the producing well.

As mentioned, our interest is to study the influence of inertial effects, gas slippage, and effective stress on the displacement of the advance front of the gas injected into the reservoir. We did this by monitoring the variation in the molar fraction of CO₂ inside the reservoir. Comparisons are made based on results obtained considering classical Darcy’s law without incorporating the mentioned effects.

For all simulations, the maximum production time is equal to 6,000 days, and we presented the molar fraction values in the $xy$-plane for a value of z corresponding to half of $L_z$ and for four different selected time instants.

Therefore, let’s start with the advance of carbon dioxide in the case of flow governed by classical Darcy’s law, see Figure 2. These are the reference results against which we will analyze the impacts caused by including the effects already referenced.

When we observe the distribution of the molar fraction, we can see that when we attained the maximum time, a portion of the CO₂, thanks to the presence of fractures, has already reached the production region, with its values being approximately equal to 0.2.

In the injection region, we found that the first square region delimited by the crossing of the fractures has not yet been filled by the injected gas, and the maximum molar fraction value is higher than 0.8.

Next, we move on to the case of Darcy’s law, modified to account for inertial effects. With the introduction of apparent permeability, in the model proposed by [38] its values are limited between $\mathbf{k}$ (for Reynolds number approaching zero) and ${\mathbf{k}}_a=k_{mr}\mathbf{k}$ in the limit of $Re$ tending to infinity, with $k_{mr}$ less than the unit. In short, in regions where the flow has higher Reynolds number values, the tendency is to have apparent permeability values lower than those considered in classical Darcy’s law.

We can see this trend in the fields presented in Figure 3. As the apparent permeability tends to be lower in fractures, the gas that advances closer to the injection will take longer to reach the production region. We also verify that the mole fraction is less than 0.2, contrary to the case in the previous example. So, we can say that inertial effects are slowing down the displacement of the injected gas.

Continuing, we focus on the effects arising from the stress field. In it, we know that permeability values can change depending on pressure variations. Its variation will grow exponentially and will be higher when the reservoir pressure is higher than the initial one.

In regions of the reservoir not yet disturbed by gas injection, the permeability value will tend to be equal to that of Darcy’s classical law. On the other hand, in those where the pressure is higher than the initial, we will have a higher value. We must remember that this effect is not taken into account in fractures.

In practical terms, we note (Figure 4) that, for the evaluated parameters, there were only little variations concerning the distribution of CO₂ when we compare the values to those in the first example. We can distinguish a smaller filling of the region delimited by the fractures in the lower left corner and a smaller molar fraction of the injected gas than in the case where the permeability is constant, upper right corner.

The last effect incorporated was that arising from gas slippage. In this model, if the Knudsen number of the mixture of components tends to zero, the permeability value would be the same as the initial. On the other hand, when it tends to infinity, its value would become five times greater than the initial value.

However, we must remember that the slip regime is characterized by $0.001\le Kn\le 0.1$. Thus, the highest viable variation would be for an apparent permeability equal to 1.4${\mathbf{k}}_0$.

Unfortunately, depending on the parameters and properties chosen, the variation range of the mixture’s Knudsen number does not provide changes that differentiate this case from that of flow using Darcy’s classical law, see Figure 5.

We can not visually notice the differences in the mole fraction values when we compare these two cases. We believe that cumulative effects appear in the long term in the amount of gas produced. However, we cannot verify this in the present work.

Finally, we would like to show the results of the simulations encompassing all the effects added to our model. Although separately, they may have a greater or lesser influence on the flow of injected carbon dioxide, their combination will lead to a different methane gas production scenario in the long term.

From what we have discerned, inertial effects most affected the flow, followed by those arising from effective stress. Differently, the gas slip did not introduce changes that could modify the flow from that obtained with Darcy’s law without any modification incorporated.

The reader can see for themselves, in Figure 6, how these effects changed the mole fraction field due to the combination of them all.

To begin with, we found that the first region between the fractures in the injection region, bottom left corner, is practically filled by CO₂, unlike when we did not consider any effect. In the methane production region, as we are injecting and producing the same volume of gas, CO₂ may arrive earlier in the production zone, and we recover the injected CO₂ along with the CH₄.

Therefore, we can say that gas production will not be the same when we disregard the combined effects and should be higher than in reality.

6. Conclusions

We address the numerical simulation of single-phase two-component gas flow in a naturally fractured tight reservoir having low permeability and porosity values.

Our main objective was to show the implications for the production of using classical Darcy’s law in simulating flow in these types of reservoirs, especially in the presence of fractures. We know that they are responsible for introducing preferential paths for flow, presenting Reynolds numbers higher than those of regions without fractures.

Therefore, we incorporated inertial, effective stress, and gas slippage effects. In the first case, the apparent permeability values vary as a function of the Reynolds number and have an upper limit value. When it comes to the effects of effective stress, permeability values grow exponentially depending on the difference between the current and initial pressure values prevailing in the reservoir. Finally, the correction due to gas slip depends on the value of the flow’s Knudsen number and has an upper limit.

We chose a reference case to represent a naturally fractured tight gas reservoir (sugar cube configuration) and to allow us to perform simulations in a domain with zero flow conditions at the boundaries and injection flow rates (CO₂) and production (CH₄) prescribed.

Compared to the case of simulation using unmodified Darcy’s law, we verified that the distribution of the molar fraction of carbon dioxide presented a different pattern when we considered all effects simultaneously. This fact indicates that production is affected and that future predictions based on simulations made without taking their influence into account may be incorrect.

From the specific case considered, we found that the inertial and effective stress effects were dominant, with the first being the one that most influenced the flow. On the other hand, we did not observe variations when considering gas slippage.

Finally, in future work, we intend to do a quantitative study on the importance of these effects, besides including the possibility of gas adsorption and a well-reservoir coupling technique.