A Comprehensive Model Coupling Embedded ... - ACS Publications

Article pubs.acs.org/EF

A Comprehensive Model Coupling Embedded Discrete Fractures, Multiple Interacting Continua, and Geomechanics in Shale Gas Reservoirs with Multiscale Fractures Kun Wang,* Hui Liu, Jia Luo, Keliu Wu, and Zhangxin Chen

Downloaded via UNIV OF LOUISIANA AT LAFAYETTE on September 25, 2018 at 20:59:44 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Department of Chemical and Petroleum Engineering, University of Calgary, Calgary, Alberta T2N 1N4, Canada ABSTRACT: Shale gas has become one of the primary energy resources during the past few years, and its impact has been profound in many countries. Hydraulic fracturing treatments are required for the development of shale gas reservoirs, and the consequent hydraulic fractures usually connect with the original small-scale natural fractures forming complex fracture networks in these reservoirs. Therefore, a model for numerical simulation, which is capable of accurately modeling naturally and hydraulically fractured reservoirs, is essential in optimization and management of such reservoirs. In this paper, we develop a comprehensive model that couples embedded discrete fractures, multiple interacting continua, and geomechanics to accurately simulate the fluid flow in shale gas reservoirs with multiscale fractures. Large-scale hydraulic fractures are described by an embedded discrete fracture method, while middle-scale and small-scale natural fractures are modeled by a multiple interacting continua method. Usually, the connection of matrix−natural fractures−hydraulic fractures−wells is considered as the main pathway for the gas flow from a reservoir to a production well. However, geomechanics effects are significant in fractured reservoirs, which may lead to the closure of fractures and a dramatic decrease in gas conductivity in the fractures during the depletion of pressure. When the permeability of fractures is close to that of matrix, the gas production directly from the pathway of matrix−hydraulic fractures−wells as well as the fluid flow in shale matrix cannot be ignored. To accurately predict the fluid flow and well performance, the geomechanics effects and all the possible connections between different regimes must be taken into account. We implement this comprehensive model in our in-house reservoir simulator and study its behavior by numerical experiments. According to the numerical simulation results, accurate and comprehensive production prediction is performed, and reasonable physical phenomena is captured. Sensitivity studies are also performed to show the impacts of different parameters on the prediction of well performance. phenomena and multiphase flow behavior cannot be accurately handled by them. The multiple continua models and the discrete fracture models are two of the most commonly used numerical methods to simulate the fluid flow in fractured reservoirs. The first multiple continua model, a dual-porosity (DP) model, was developed by Warren and Root,4 which separates each grid block into two sub-blocks (Figure 1a) and employs them to represent fractures and matrix, respectively. A DP model assumes that there is no connection between neighboring matrix sub-blocks, and these sub-blocks act as a sink/source term to a fracture system. A dual-porosity dual-permeability (DPDP) model is an enhancement of a DP model, which allows a fluid transfer between matrix sub-blocks. In refs 5 and 6, a multiple interacting continua (MINC) model was developed, and DP and DPDP models can be seen as two special situations of a MINC model. For a MINC model, each grid block is divided into a sequence of nested sub-blocks (Figure 1b), and the nested sub-blocks have their own variables, such as pressure and saturation. Multiple continua models are the conventional models for simulating naturally fractured reservoirs and have been implemented by many commercial simulators. All these models assume that fractures are regular

1. INTRODUCTION Shale gas has become the focus of an important energy resource around the world, and its impact has been profound in the United States. A shale gas reservoir consists of complicated gas flow regimes, including gas adsorption/desorption, surface diffusion, and other fluid flow regimes. Since the permeability of shale matrix is extremely small, hydraulic fracturing techniques are widely used for the development of shale gas reservoirs. Complex fracture networks with a wide range of spacial scales are often associated with hydraulic fracturing. The conductivity of gas in fractures is sensitive to a stress change. These features make the mathematical modeling and numerical simulation of shale gas reservoirs a big challenge. Researchers have made significant efforts in this area in the past few years. Analytical and semianalytical approaches have been developed to predict the fluid flow and well performance in shale gas reservoirs with complex fracture networks. Zhou et al.1 proposed a semianalytical model that combines a numerical solution and an analytical solution on discretized fracture panels. Ye et al.2 developed a semianalytical model with complex nonplanar fractures, and Yang et al.3 performed studies on various gas flow mechanisms. These analytical and semianalytical methods are capable of predicting well production faster than numerical simulations. However, significantly simplified assumptions are applied in these methods, and the coupling of highly nonlinear physical © 2017 American Chemical Society

Received: February 9, 2017 Revised: April 29, 2017 Published: July 11, 2017 7758

DOI: 10.1021/acs.energyfuels.7b00394 Energy Fuels 2017, 31, 7758−7776

Article

Energy & Fuels

matrix, natural/induced fractures. Ding et al.20 have a similar idea and develop an EDFM-MINC approach. Their models assume that the large-scale fractures represented by an EDFM first interact with the small-scale fractures represented by the fracture medium of a MINC model, and the direct mass transfer to the large-scale fractures from the matrix medium is neglected. If geomechanics effects are taken into account, the permeability of fractures decreases dramatically and approaches the permeability of matrix. Therefore, the mass transfer from the matrix medium to large-scale fractures must be considered. In this paper, we develop a comprehensive model coupling a MINC model and an EDFM. The geomechanics effects and the mass transfer between large-scale fractures and matrix as well as between neighboring matrix blocks are included in this model. We present simulation results to demonstrate the applicability and accuracy of the developed model. Through the simulation results, we also compare our model with other models in the literature and perform sensitivity studies to show the factors that may affect the differences in the predicted production between different coupled models.

2. GAS STORAGE AND TRANSPORT MECHANISMS IN SHALE RESERVOIRS It is generally agreed that a shale matrix is comprised of organic matter and inorganic matter. Organic and inorganic matter have different gas storage and transport mechanisms,21 such as gas adsorption/desorption, non-Darcy’s flow, surface diffusion of adsorbed gas,22 and stress dependence of shale rock and fractures.23,24 In this paper, we consider the gas adsorption/ desorption, which happens on the surface of organic matter, the gas flow in nanoscale pores, which is included by an apparent permeability correction in a general form, and the stressdependent permeability of shale matrix and fractures. 2.1. Gas Adsorption/Desorption. Gas adsorption/desorption is one of the features of shale gas reservoirs and makes significant contributions to original gas storage in place. The quantity of the adsorbed gas on pore walls of organic matter in shale is

Figure 1. (a, b) Multiple continua models.

and well connected in a reservoir. In ref 7, a MINC model was developed for the simulation of single-phase and multiphase flow in naturally fractured vuggy reservoirs. A dual-porosity model proposed in ref 8 was developed for the fluid flow in microscale porous media. In ref 9, Wu et al. adopted a MINC model, and a generalized mathematical framework model and numerical approach was discussed. However, this assumption is often not satisfied, and the approach is not applicable to modeling large-scale disconnected fractures. Compared with multiple continua models, discrete fracture models employ unstructured grids to describe fractures and treat the fractures in a lower dimension.10−12 To avoid difficulties in the generation of unstructured grids, embedded discrete fracture models (EDFM) are proposed in refs 13−15. An EDFM discretizes fractures into segments, treats the mass transfer between fracture segments and matrix through transfer terms, and accurately predicts well production. A simulation of fractured reservoirs with an EDFM explicitly representing all the small-scale and large-scale fractures, however, is not practical due to its expensive computational cost. To efficiently model reservoirs with multiscale fractures, Jiang et al.16 proposed two hybrid models: The first one couples an EDFM with a MINC model and the second one uses a DFM to describe large-scale fractures with an unstructured grid and employs the MINC method to simulate small-scale fractures. Then Jiang et al.17,18 adopted their models for simulating multiphase flow in shale gas reservoirs with a complex fracture system. Zhang et al.19 employ a discrete fracture model with unstructured grids employed to describe hydraulic fractures, and a dual-porosity dual-permeability model was adopted to represent the stimulated region, including

mg = (1 − ϕ)ρρ V r gs E

(1)

where ϕ is the rock porosity, ρr is the rock bulk density, ρgs is the gas density at standard condition, and VE is the adsorption isotherm function (m3/kg). The adsorption isotherm function can be expressed by pg VE = VL pg + pL (2) where VL is the Langmuir volume (the maximum adsorption capacity at a given temperature, m3/kg), pg is the gas pressure, and pL is the Langmuir pressure (the pressure corresponding to one-half of the Langmuir volume, Pa). If real gas is considered, eq 2 can be written as pg

VE = VL

Z pg Z

+ pL

(3)

where Z is the compressibility factor. 2.2. Apparent Permeability. Shale gas formations are characterized by pore sizes from one to a hundred nanometers. Darcy’s law is not adequate to describe the various nonviscous gas flow regimes. In this paper, we used apparent permeability25 7759


Article

Energy & Fuels to present the non-Darcy gas flow mechanism in nanopores, which can be represented by K am = K im·f (K n)

employ an EDFM to handle the large-scale hydraulic fractures. On the basis of different assumptions, three coupled models are introduced. 3.1. Multiple Interacting Continua Model in the Shale Matrix and Natural Fractures. Since the permeability is ultralow in a shale formation, the transient flow may occur for a long period of time inside the matrix, and it is improper to treat the flow between fractures and matrix as a pseudosteady state. Thus, it is necessary to divide each matrix block into several sub-blocks to capture its large potential gradients. A MINC model can resolve the gradients in detail and give more accurate production prediction. In addition, the MINC model can employ different sub-blocks to represent organic and inorganic matter separately, and treat the different flow mechanisms in organic and inorganic matter in a realistic way. Figure 2 shows a two-dimensional (2D) FIB/SEM shale

(4)

Km i

where is the intrinsic permeability of matrix at the initial condition, Kn is the Knudsen number, and f(Kn) is a correction factor. The correction factor is defined by ⎛ 4K n ⎞ f (K n) = (1 + αK n)⎜1 + ⎟ 1 + Kn ⎠ ⎝

(5)

where α is the rarefaction parameter and can be expressed as α=

128 tan−1(4K n0.4) 15π 2

(6)

The Knudsen number Kn is the ratio of the mean free path of gas at reservoir conditions and a hydraulic radius:

Kn =

λ Rh

(7)

The mean free path λ of real gas is λ=

μ pg

πZRT 2M

(8)

and the hydraulic radius Rh is R h = 2 2τ

K im ϕ

(9)

where R is the gas constant (8314 J/(kmol K)), T is the temperature (K), τ is the tortuosity, and M is the molecular weight (kg/kmol). 2.3. Geomechanics Effects. The effect of geomechanics is more significant in fractured reservoirs than that in conventional reservoirs because fractures are much more stresssensitive than the rock matrix. During production, the depletion of pressure leads to a stress change and fracture closure, which may significantly reduce the fluid transmissibility in the fractures. In spite of field and experimental research that has demonstrated the stress-dependent fracture properties, rock porosity and permeability are still often treated as static parameters in fractured reservoir simulations. In ref 24, it was assumed that a natural fracture permeability change is related to the mean normal stress and can be expressed as

Figure 2. Two-dimensional FIB/SEM image: The black depicts pores; the dark gray is organic matter; the light gray is inorganic matter.26

image, from which we can see that the volume of organic matter is smaller than that of inorganic matter, and the organic matter is usually surrounded by the inorganic matter. Hence it is reasonable to use the inner sub-blocks to represent the organic matter and the outer sub-blocks to describe the inorganic matter. Since a large number of natural fractures are distributed in a shale reservoir, it is often unlikely to locate the pre-existing natural fractures, and a huge computational cost makes it unpractical to explicitly represent the natural fractures. Therefore, it is more reasonable to use a most outer sub-block in the MINC model to represent natural fractures. Figure 3 illustrates the MINC model for shale matrix with natural fractures. The sub-blocks in the MINC model connect with each other by the matrix−matrix and matrix−fracture transfer terms. Let M be the number of matrix sub-blocks in each grid block; the matrix−matrix transfer terms between the kth matrix sub-block and the (k + 1)th matrix sub-block as well as between the matrix sub-block and fracture sub-block for the gas phase can be defined by

0

K f = K if ecf {(1 + vf )/(1 − vf )}αf (pg,f − pg,f )

(10)

where Kfi is the fracture permeability at an initial condition, cf is the fracture compressibility (1/Pa), vf is Poisson’s ratio, αf is a Biot coefficient, and p0g,f is the initial reservoir pressure. We employ a similar formula to calculate the permeability change in shale matrix and hydraulic fractures: 0 ⎧ m m c {(1 + vm)/(1 − vm)}αm(pg,m − pg,m ) ⎪ K = Ki e m ⎨ 0 ) ⎪ K F = K FecF{(1 + vF)/(1 − vF)}αF(pg ,F − pg,F ⎩ i

⎛ K rgρ ⎞ g ⎟ (p − p ) qg, l l = σ ⎜⎜ ⎟ l2 l1 12 μ ⎝ g ⎠t

(11)

3. A COMPREHENSIVE MODEL IN FRACTURED SHALE GAS RESERVOIRS In order to accurately predict the production of shale gas reservoirs with multiscale fractures, we adopt a MINC model to describe the shale matrix and small-scale natural fractures and

(l1 , l 2) = (f , m1), (m1 , m2), ..., (mM − 1 , mM )

(12)

where t in (·)t is equal to l1 or l2 depending on which sub-block is upstream, f is the fracture sub-block, mi, i = 1, 2, ..., M are the 7760


Article

Energy & Fuels

from the fracture. For a 3D problem, see Figure 4, and ⟨d⟩ is calculated by

Figure 3. MINC model for the shale matrix. Figure 4. A grid block that contains a fracture.

matrix sub-blocks, Krg is the gas relative permeability, ρg is the gas density, μg is the gas viscosity, and σ is the shape factor. In this paper, the shape factor developed by Gilman and Kazemi27 is employed. 3.2. Embedded Discrete Fracture Model in Hydraulic Fractures. The assumption in a MINC model that fractures are distributed uniformly and are well connected does not conform to the fracture distribution resulted from hydraulic fracturing. Consequently, a multiple continua model is not adequate to model large-scale and disconnected hydraulic fractures. A discrete fracture model (DFM) is a new type of method in fractured reservoir simulation, which employs unstructured grids to provide more realistic representations of fractures. Fractures are represented by the boundaries of grid blocks. The flow through fractures is simulated in a lower dimension, i.e., 2D in a three-dimensional (3D) reservoir simulation case and one-dimensional (1D) in a 2D reservoir simulation case. The flow in matrix is coupled with that of fractures by a reasonable boundary condition. However, the generation of grids that conform to each fracture network can be a big challenge. An embedded DFM (EDFM) proposed in ref 13 uses a structured grid for the matrix representation, introduces additional grids for fractures, and models the flow flux between fractures and matrix through transfer terms. With an EDFM, it is not necessary to generate complex unstructured grids, and it is easy for this approach to incorporate with the commonly used finite difference (volume) methods. To couple an EDFM with a MINC model, four types of connections are taken into consideration: (1) Type I: connections between matrix and hydraulic fractures; (2) Type II: connections between two intersecting hydraulic fractures; (3) Type III: connections between two fracture cells of an individual hydraulic fracture; and (4) Type IV: connections between natural fractures and hydraulic fractures. The transmissibility in these connections can be calculated as follows: • Transmissibility in connection Type I: The EDFM assumes that the distribution of pressure around a fracture is linear. Then the transmissibility between fractures and matrix in each grid block can be calculated by

∫V n ⃗ ·x ⃗ dv

⟨d⟩ =

(14) V where n⃗ is the unit normal vector, x⃗ is the normal distance of the element from the fracture, dv is the volume element, and V is the volume of a grid block. • Transmissibility in connection Type II: For the fracture cells interacting in a grid block, see Figure 5, and the transmissibility is approximated by

TF1F2 =

TF1TF2 TF1 + TF2

(15)

and TF1 = TF2 =

K F1ωF1L F1 d F1

(16)

K F2ωF2L F2 d F2

(17)

where ωF and KF are the fracture aperture and permeability, respectively; LF is the length of the fracture bounded in a grid block; and dF is the average distance from the center of the fracture bounded in a grid block to the center of the intersection line. • Transmissibility in connection Type III: For two fracture cells of an individual fracture in two separate grid blocks, see Figure 6, and the transmissibility is calculated by TF1F2 =

KFA F12 DF12

(18)

where DF12 is the distance between two fracture segments, which is DF12 = DF1 + DF2

(19)

and AF12 is the fracture aperture times the length of the intersection line, which is A F12 = L FωF

KA ̅ |⟨d⟩| (13) where K̅ is the harmonic average for the permeability of the fracture and matrix, A is the fracture surface area in the grid block, and ⟨d⟩ is the average normal distance TmF =

(20)

• Transmissibility in connection Type IV: On the basis of the assumption in the MINC model that natural fractures are distributed uniformly, by extending eq 13, the transmissibility in connection Type IV is calculated by 7761


Article

Energy & Fuels

Figure 5. (a, b) Transmissibility of interacting fractures.

Figure 6. (a, b) Transmissibility of one fracture in two grid blocks.

• The coupled model I. In this model, it is assumed that the hydraulic fractures only interact with the natural fracture blocks in the MINC model. Since the permeability of the matrix is much lower than that of natural fractures, a direct mass transfer between the matrix medium and the hydraulic fractures is much less than that between natural fractures and hydraulic fractures and thus can be neglected. The connection between the neighboring matrix sub-blocks is also neglected, and the shale matrix only behaves as a sink/ source term of the natural fracture sub-blocks. The connection list of the grid blocks in this model is shown in Figure 8a. Actually, the connection list from the models proposed in the literature16−20 can be included by the coupled model I. • The coupled model II. Because of an irregular distribution of natural fractures, some grid blocks may not contain natural fractures but have hydraulic fractures passing through them. In this situation, the mass transfer to the hydraulic fractures is only possible from the matrix medium. More importantly, during the gas production, if the geomechanics effects are taken into consideration, the depletion of pressure leads to the closure of fractures and a dramatic decrease in the fracture permeability. Then the fracture permeability is close to the matrix permeability, and the mass transfer directly from shale matrix to hydraulic fractures should not be neglected. Hence the connection between shale matrix and hydraulic fractures should be taken into consideration. The grid connection list in the second coupled model is shown in Figure 8b.

K̅ fFA |⟨d⟩| (21) where K̅ fF is the harmonic average for the permeability of a natural fracture and a hydraulic fracture. 3.3. Coupled Models. Taking the grid shown in Figure 7 as an example, there are four grid blocks Ci, i = 1, 2, 3, 4 and one TfF =

Figure 7. Example of grid blocks for the coupled models.

hydraulic fracture passing through three blocks, C1, C2, and C4. Each grid block is divided into four sub-blocks, in which the most outer sub-blocks, f i, i = 1, 2, 3, 4, represent the natural fractures, the most inner sub-blocks, mi,3, i = 1, 2, 3, 4, are for organic matter, and inorganic matter is represented by the subblocks mi,j, i = 1, 2, 3, 4, j = 1, 2. Three different block connecting relations are listed in Figure 8a−c. On the basis of these connection relations, we have the following three coupled models: 7762


Article

Energy & Fuels Oi Oi ∂(ϕOisgOiρgOi ) + (1 − ϕOi)ρρ V r gs E

∂t

= qgOi − 1 − Oi − qgOi − Oi + 1

+ qgOi − F , i = 1, ..., MO

(22)

O

where M is the number of the organic sub-blocks in each grid block, qOi−Oi+1 is the mass transfer term between the ith organic sub-block and (i + 1)-th organic sub-block, and qOi−F is the mass transfer term between the ith organic sub-block and an adjacent hydraulic fracture. qOMO−OMO+1 represents the mass transfer term between the most outer organic sub-block and the most inner inorganic sub-block, and qO0−O1 is equal to zero. The gas equation for inorganic matter is ∂(ϕ IisgIiρgIi ) ∂t

⎛ K Ii ρ Ii ⎞ rg g = ∇·⎜⎜KaIi I ∇ΦIgi⎟⎟ + qgIi − 1− Ii − qgIi − Ii + 1 + qgIi − F , μg i ⎝ ⎠

i = 1, ..., MI

(23)

where MI is the number of the inorganic sub-blocks in each grid block, qIi−Ii+1 is the mass transfer term between the ith inorganic sub-block and (i + 1)-th inorganic sub-block, and qOi−F is the mass transfer term between the ith inorganic sub-block and an adjacent hydraulic fracture. qIMI−IMI+1 represents the mass transfer term between the most outer inorganic sub-block and an adjacent natural fracture sub-block, and qI0−I1 represents the mass transfer term between the most inner inorganic sub-block and the most outer organic sub-block, which has qI0−I1 = −qOMO−OMO+1. The gas equation in natural fractures is

∂(ϕf sgf ρgf ) ∂t

⎛ K f ρf ⎞ rg g f = ∇·⎜⎜K f f ∇Φgf ⎟⎟ + qgf − IM + qgf − F − qg,W μg ⎝ ⎠

(24)

f−IMI

is the mass transfer term between the most outer inorganic where q sub-block and an adjacent natural fracture, qf−F is the mass transfer term between adjacent natural and hydraulic fractures, and qfg,W is the well production from natural fractures. Similarly, we have the water equations for organic matter: Figure 8. Example of the connection list for the coupled models.

∂(ϕOiswOiρwOi )

• The coupled model III. In the above two models, the matrix sub-blocks disconnect with their neighboring matrix sub-blocks. However, the shale matrix may have equal conductivity for fluid flow as natural fractures, then the shale matrix becomes one of the main media for fluid flow, and the connection of the matrix sub-blocks with their neighboring grid blocks is also needed to consider, which results in the coupled model III. Since the shale matrix has been divided into a series of nested sub-blocks in the MINC model, it is reasonable to assume that the fluid can flow between the most outer neighboring matrix sub-blocks. Actually, the coupled model III has the same consideration as a dual permeability model. The connection list of the coupled model III can be seen in Figure 8c.

∂t

= qwOi − 1 − Oi − qwOi − Oi + 1 + qwOi − F i = 1, ..., MO

(25)

for inorganic matter: ∂(ϕ IiswIi ρwIi ) ∂t

⎛ K Ii ρ Ii ⎞ rw = ∇·⎜⎜K Ii I w ∇ΦIwi ⎟⎟ + qwIi − 1− Ii − qwIi − Ii + 1 + qwIi − F , i μw ⎝ ⎠

i = 1, ..., MI

(26)

and in natural fractures:

∂(ϕf swf ρwf ) ∂t

⎛ K f ρf ⎞ rw w f = ∇· ⎜⎜K f ∇Φfw ⎟⎟ + qwf− IM + qwf− F − qw,W f μ ⎝ ⎠ w (27)

The governing equations for gas and water in hydraulic fractures are

We will perform numerical experiments and study the behavior of the three coupled models in 5, from which we can conclude the coupled model III is the most comprehensive one to capture reasonable physical phenomena and accurately predict production.

∂(ϕFsαFραF ) ∂t −

⎛ K F ρF ⎞ rα = ∇·⎜⎜K F F α ∇ΦαF ⎟⎟ − μα ⎝ ⎠

∑

∑

qαOi − F

i = 1, ..., MO

qαIi − F − qαf − F − qαF,W , α = g , w

i = 1,..., MI

(28)

There are also other constraint equations:

4. GOVERNING EQUATIONS AND NUMERICAL METHODS

Φα = pα + ρα , Nz , α = w , g sw + sg = 1

In this paper, a two-phase gas−water system is considered. The gas

pw = pg − pcgw

equation for organic matter is 7763


Article

Energy & Fuels ⎧A Ii − 1 Ii Ii + 1 F ⎪ IiIi − 1δ V + A IiIi δ V + A IiIi + 1δ V + A IiFδ V ⎪ = b Ii , i = 1, ..., MI ⎪ ⎪ A fI M δ V I MI + A ff δ V f + A fFδ V F = bf ⎨ I ⎪ ⎪ ∑ AFOi δ V Oi + ∑ AFIi δ V Ii + AFf δ V f ⎪ i = 1,..., MO i = 1, ..., MI ⎪ F F ⎩ + AFFδ V = b

In eqs 22−28, O, I, f, and F are for organic matter, inorganic matter, natural fractures, and hydraulic fractures, respectively, ϕ and K are porosity and permeability, for phase α (α = w, g) Φα is the phase potential, and sα, μα, pα, Bα, ρα, Krα, and qα are the saturation, viscosity, pressure, format volume factor, density, relative permeability and production rate, respectively, pcgw is the capillary pressure between the water and gas phases, N is the gravitational constant, and z is the reservoir depth. With proper boundary and initial conditions, a close system is obtained. To solve the nonlinear eqs 22−28, a fully implicit Newton− Raphson iterative method is employed with the finite difference (volume) method. We choose the pressure (pOi, i = 1, ···, MO, pIi, i = 1, ···, MI, pf, and pF) and gas saturation (sOg i, i = 1, ···, MO, sIgi, i = 1, ···, MI, ∂ sfg, and sFg ) as the primary variables. For a cumulative term ∂t (C), it can be expressed as

(33)

Therefore, a linear system is established AV = b

(34)

where

1 n+1 ∂ (C) = [C − C n] ∂t Δt where Cn+1 and Cn are the variables at the (n + 1)-th and nth time steps. For a flow term ∇(K T∇Φ), it can be discretized as 1 n+1 n+1 ∇(KT ∇Φ) = K i + 1/2, j , kTin++1/2, j , k(Φi + 1, j , k − Φi , j , k ) − K i − 1/2, j , k 1 n+1 n+1 n+1 Tin−+1/2, j , k(Φi − 1, j , k − Φi , j , k ) + Ti , j + 1/2, kTi , j + 1/2, k

and

(Φin,+j +11, k − Φin,+j , k1) − Ti , j − 1/2, kTin, j+−11/2, k(Φin,+j −11, k − Φin,+j , k1)

⎛ δ V O1 ⎞ ⎜ ⎟ ⎜ ⋯ ⎟ ⎜ δ V O MO ⎟ ⎜ ⎟ ⎜ δ V I1 ⎟ V=⎜ ⋯ ⎟ ⎜ ⎟ ⎜ δ V I MI ⎟ ⎜ f ⎟ ⎜ δV ⎟ ⎜ ⎟ ⎝ δ VF ⎠

+ Ti , j , k + 1/2Tin, j+, k1+ 1/2(Φin,+j , k1+ 1 − Φin,+j , k1) − Ti , j , k − 1/2Tin, j+, k1− 1/2 (Φin,+j , k1− 1 − Φin,+j , k1), where (ζ)i, j, k is the variable ζ at the grid block (i, j, k) and the transmissibility K i ± 1/2, j , k , K i , j ± 1/2, k , and K i , j , k ± 1/2 can be calculated as shown in 3.2. An upwind method is used to calculate Ti + 1/2, j , k as

⎧ Ti , j , k , if Φi , j , k > Φi − 1, j , k ⎪ Ti + 1/2, j , k = ⎨ ⎪ ⎩Ti + 1, j , k , if Φi , j , k ≤ Φi − 1, j , k

(29)

and the same method for Ti − 1/2, j , k , Ti , j ± 1/2, k , and Ti , j , k ± 1/2 . At each nonlinear iteration, the increment of the primary variables is to be solved. Taking the pressure as an example, we represent the pressure at the (n + 1)-th time step as

pn + 1, l + 1 = pn + 1, l + δp

The linear system will be solved at each nonlinear iteration and the solutions will be used to update the variables until the nonlinear residual reaches the stop criterion.

5. NUMERICAL EXPERIMENTS AND DISCUSSIONS We implement the three coupled models in our in-house simulator.28−30 Our simulator is built on a parallel platform PRSI.31 PRSI is written by C language and Message Passing Interface (MPI), which is designed for large-scale reservoir simulation and provides a grid management module, a parallel input and output module, a linear solver module, a matrix/ vector module, and a well management module. The constrained pressure residual (CPR)-type preconditioners,28,32−36 which are considered as the most efficient and parallel scalable preconditioners in reservoir simulation, are also implemented in our simulator. The simulator is capable of simulating large-scale problems with hundreds of millions of grid blocks. The simulator has also been coupled with a parallel geomechanics simulator37,38 for large-scale fractured unconventional reservoir simulations. 5.1. Behavior of the Coupled Models. The first case describes a reservoir with natural fractures and considers one horizontal well with four stages of hydraulic fracturing. The detailed parameters of the reservoir are summarized in Table 1, and the position of the horizontal well and hydraulic fractures is shown in Figure 9a−b. The operation constraints of the horizontal well are listed in Table 2, in which “STG” and

(30)

where pn+1,l is pressure from the lth nonlinear iteration. Other variables ζ that depend on the pressure, such as porosity and density, are represented as

⎛ ∂ζ ⎞n + 1, l ζ n + 1 = ζ n + 1, l + ⎜ ⎟ δp ⎝ ∂p ⎠

(31)

Consequently, from eqs 22 and 25, we can have the following equations for organic matter A OiOi − 1δ V Oi − 1 + A OiOi δ V Oi + A OiOi + 1δ V Oi + 1 + A OiFδ V F = bOi , i = 1, ..., MO

(36)

(32)

where δVm = (δpm, δsmg ) are the increment of the primary variables for different media m; AOiOi−1, AOiOi+1, and AOi F are the matrices resulted from the discretization of the transfer terms qOα i−1−Oi, qOα i−Oi+1, and qOα i−F, α = g, w; AOiOi is the matrix resulted from the discretization of the cumulative term and transfer terms; and bOi is the right-hand side. Similarly, we can also draw the equations for inorganic material, natural fractures and hydraulic fractures from eqs 23, 26, 24, 27, and 28 as 7764


Article

Energy & Fuels Table 1. Parameters of the Base Case parameter grid blocks block dimension (x, y, z) organic material porosity inorganic material porosity organic material permeability inorganic material permeability organic/inorganic volume ratio V ( O)

Table 2. Well Constraints

value 60 × 51 × 1 (12.5, 12.5, 12.5) 0.4 0.2 (1 × 10−4, 1 × 10−4, 1 × 10−4) (1 × 10−3, 1 × 10−3, 5 × 10−4) 1:1

unit

time

constraint

0

Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP:

ft The 50th day md

The 100th day

md

The 150th day The 200th day

VI

hydraulic fracture porosity hydraulic fracture aperture hydraulic fracture permeability natural fracture spacing natural fracture aperture natural fracture permeability natural fracture shape factor Langmuir volume Langmuir pressure organic material density temperature Biot coefficient Poisson’s ratio Simulation time

0.9 0.01 (5000, 5000, 2000) (10, 10, 10) 0.001 (125, 125, 125) GK 0.08 10000 400 157.73 1.0 0.33 365

The 250th day ft md ft ft md

The 300th day

2000 MMScf 6000 psia 1800 MMScf 4000 psia 1600 MMScf 3000 psia 1400 MMScf 2000 psia 1200 MMScf 1000 psia 1000 MMScf 500 psia 800 MMScf 200 psia

Table 3. Geomechanics Parameters I

ft3/lbm psi lbm/ft3 F

Parameters I

material

value

unit

matrix natural fracture hydraulic fracture

1e-7 3e-4 3e-5

1/psi 1/psi 1/psi

we expect, almost the same cumulative gas production is predicted by the coupled model II and III, which is slightly larger than that by the coupled model I; see Figure 10. Figure

day

“BHP” are the gas production rate and well bottom-hole pressure, respectively. On the basis of this case, we compare the simulation results of the three coupled models discussed in 3. We first compare the behavior of the coupled models for the situations without and with geomechanics effects. In this test, we use one sub-block to represent the natural fracture and two sub-blocks to represent the matrix, in which the inner one is for the organic matter and the outer one is for the inorganic matter. When the geomechanics effects are taken into consideration, it is assumed that the shale matrix has much smaller compressibility than the fractures and that the compressibility of natural fractures is larger than that of hydraulic fractures due to the existence of proppants in hydraulic fractures. The geomechanics parameters shown in Table 3 are used. When the geomechanics effects are not considered, the permeability of fractures is constant and much larger than that of matrix, and most of the produced gas is from the pathway of matrix-NF (natural fractures)-HF (hydraulic fractures)-well. As a result, the connection between hydraulic fractures and matrix contributes little to the gas production and can be neglected. As

Figure 10. Cumulative gas production without geomechanics effects.

11a−c shows the profile of gas saturation in natural fractures, from which we can see a similar phenomenon presented by the three different coupled models: The grid blocks containing hydraulic fractures have lower gas saturation than their neighboring blocks. The minimum gas saturation in the natural fractures from the coupled model I is approximately 0.164, while the lowest values from the coupled models II and III are 0.386 and 0.348, respectively. That is because the connection

Figure 9. (a, b) Reservoir description of the first case. 7765


Article

Energy & Fuels

Figure 11. (a−c) Gas saturation profiles in the cases without geomechanics effects.

Figure 12. (a, b) Permeability change in the cases with geomechanics parameters I.

of natural fractures is even smaller than that of organic and inorganic matter. The simulation results of the coupled models II and III show that the gas in the hydraulic fractures mainly comes from the direct connection of NF-HF at the beginning of production, while the direct connection of matrix-HF contributes an equal amount of gas to the hydraulic fractures as the connection of NF-HF does after pressure drops. We calculate the gas mass transfer between different media and summarize it in Figure 14a−b. We can see that (1) if no geomechanics effects are taken into consideration, only a small proportion of gas comes from the connection between the hydraulic fractures and matrix, Figure 14a; and (2) for the

between matrix and hydraulic fractures is considered by the coupled model II and III, and thus the shale matrix contributes more gas to the hydraulic fractures. Therefore, the coupled models II and III have more gas production and higher gas saturation in natural fractures. In general, the geomechanics effects are significant for fractured unconventional reservoirs and should be modeled in numerical simulation; otherwise, the closure of fractures will not be captured by simulation, the gas conductivity in fractures will be constant during a pressure drop, and thus the cumulative gas production is overestimated. Figure 18a shows the cumulative gas production when the geomechanics effects are taken into consideration. By comparing Figures 10 and 18a, we can see a big difference on production. During the depletion of reservoir pressure, when the geomechanics parameters I are used, the permeability of fractures decreases dramatically. The intrinsic permeability of natural fractures drops from 100 md to 1 md during production; see Figure 12a. Since the apparent permeability correction is used to simulate nonviscous gas flow in nanoscale pores in this paper, the effective permeability of matrix increases as shown in Figure 12b, and the effective permeability

Figure 13. (a, b) Diagram of gas flow in fractured reservoirs: The solid line means a major mass transfer while the dash line means a minor mass transfer. 7766


Article

Energy & Fuels

Figure 14. (a, b) Gas transfer between hydraulic fractures and natural fractures, hydraulic fractures and organic matter, and hydraulic fractures and inorganic matter.

gas in the hydraulic fractures comes from the connection of matrix-HF if there is no geomechanics effect, while the number increases to about 35% for the coupled model II and 45% for the coupled model III. The diagram of gas flow is shown in Figure 13a−b, in which the solid lines mean a major mass transfer, while the dash lines mean a minor mass transfer. Hence the gas saturation profile from the coupled model I is totally different from that from the coupled models II and III; see Figure 16a−c. The gas saturation change history in three selected grid blocks (the locations of the three blocks are shown in Figure 9a) is plotted in Figure 17. As we can see, compared with the coupled model I, the coupled models II and III give a significantly different prediction on the gas saturation change. At point I, natural fractures play the role as the source of well production. The gas flow rate coming from the shale matrix to the natural fractures is smaller than that flowing from the natural fractures to the hydraulic fractures, which results in the decreasing of gas saturation in natural fractures. During production, the reservoir pressure decreases, and consequently

geomechanics parameters I, a huge amount of gas in the hydraulic fractures directly comes from the matrix according to the simulation results of the coupled models II and III, Figure 14b. It can be seen more clearly in Figure 15: Only 15% of the

Figure 15. Percentage of gas transfer between hydraulic fractures and natural fractures, hydraulic fractures and organic matter, and hydraulic fractures and inorganic matter.

Figure 16. (a−c) Gas saturation profile in the cases with geomechanics parameters I. 7767


Article

Energy & Fuels

Figure 17. (a−c) Gas saturation history in the cases with geomechanics parameters I at selected points.

Figure 18. (a, b) Geomechanics effects on cumulative gas and water production when geomechanics effects are considered.

Figure 19. (a, b) Reservoir description of the second case.

source of production, and the gas saturation should decrease. When the coupled model II is employed, due to the connection of shale matrix and hydraulic fractures and a larger transmissibility between them (compared with the transmissibility at points I and II), the gas flowing into natural fractures is more than that flowing out from natural fractures and then the gas saturation increases. For the cumulative gas production, shown in Figure 18a, there is 10% percent more cumulative gas production predicted by the coupled model III. As seen from Figure 18b, the cumulative water production resulting from the coupled model I is totally different from that from the other

the fracture permeability dramatically decreases. In this situation, the shale matrix becomes the main source of the well production. The gas flowing from the shale matrix to the natural fractures is more than that from the natural fractures to the hydraulic fractures, which leads to the increasing of gas saturation in natural fractures. The pressure drops faster at point III than at point I, and hence the gas saturation decreases in a shorter period of time and then starts to increase when the coupled models II and III are used. At point II, if the coupled model I is used, the gas can only flow through the pathway of matrix-NF-HF, the natural fractures always play the role as the 7768


Article

Energy & Fuels

predicted by the coupled model III is the highest, and much more water is produced by the coupled model I; (3) the gas saturation change history resulting from different models at a selected point is significantly different; and (4) when the coupled model III is used, the gas flow which is directly from the shale matrix to the hydraulic fractures is close to that from the natural fractures to the hydraulic fractures, and the gas that comes from the connection of matrix-HF is about one-third of the total produced gas. On the basis of the above numerical experiments and discussion, some conclusions are made. During the pressure drops, the effective permeability of natural fractures becomes close to or even smaller than that of matrix. As a result, the pathway of matrix-HF-well makes a nonignorable contribution to the well production. Thus, it is essential to consider the connection of matrix-HF in the coupled models; otherwise, a coupled model may underestimate the production and it leads to unreasonable saturation distribution. Compared with the coupled model III, the coupled model II underestimates the cumulative gas production even when the connection between matrix and hydraulic fractures is included. The reason is that the shale matrix becomes one of the main media for gas flow; however, it is blocked by the coupled model II. The coupled model III allows the gas flow between the neighboring matrix sub-blocks so more gas can transfer from other matrix blocks to the matrix blocks where hydraulic fractures are located, which leads to more gas produced from the pathway of matrix-HF. Since the coupled model III takes the mass transfer between hydraulic fractures and matrix as well as the mass transfer between the neighboring matrix sub-blocks into consideration, the coupled model III can more accurately describe the fluid flow in fractured reservoirs. From the test cases, we can see that each of the connections considered by the coupled model III is essential and may result in significantly different simulation results both on the gas distribution and cumulative production. Therefore, it can be concluded that the coupled model III is the most accurate and reasonable model to simulate fractured reservoirs among the three models developed. 5.2. Sensitivity Studies. To study the effects of different parameters on the coupled models, we ran a series of simulation cases in this section. These cases are based on the first case in 5.1. We first examine the geomechanics effects with different material compressibility and see how much difference the coupled models may generate with their variations. Then the effects of the organic/inorganic matter ratio in the shale matrix, the quantities of matrix sub-blocks, and the number of hydraulic fractures are also tested. At last, the cumulative gas production is compared to see the behavior of the coupled models on different periods of time during an oilfield development. 5.2.1. Geomechanics Effects. Two sets of geomechanics parameters shown in Table 5, which have larger fracture compressibility than the geomechanics parameters I, are tested. The cumulative gas production is plotted in Figure 23a−f. From Figure 23a,b, we can see that the increasing of the hydraulic fracture compressibility does not significantly increase the difference of the prediction of cumulative gas production. However, when we increase the natural fracture compressibility from 3e-4 to 4.5e-4 and 6e-4, the natural fracture permeability becomes more sensitive to a stress change and the gas conductivity in fractures decreases dramatically, which makes the connection of matrix-HF the main pathway for production

two models. The reason is that, for the coupled model I, the water saturation as well as the water conductivity increase during the decreasing of the gas saturation, which leads to a large amount of water production, while the gas saturation predicted by the coupled models II and III does not drop as the coupled model I and hence the water production is much less. In the second case, we consider a larger reservoir with complex fracture networks. The reservoir is divided into 165 × 165 × 1 grid blocks. The dimension of each grid block is (12.5 ft, 12.5 ft, 25 ft). There is one horizontal well with six stages of hydraulic fracturing. The main hydraulic fractures connect with the secondary fractures and form irregular geometric patterns, shown in Figure 19a−b. In this case, natural fractures and hydraulic fractures exist. The natural fractures are of small scale and it is impractical to explicitly describe them by the embedded discrete fracture model. The geometry of largescale hydraulic fractures is complex and irregular; therefore it is difficult for the multiple continua methods to accurately simulate them. However, the coupled models developed in this paper can naturally handle this kind of fractured reservoir. The production schedule of the horizontal well is listed in Table 4. Other parameters are from Table 1 and the geomechanics parameters I are used in this case. Table 4. Case II: Well constraints time 0 The 50th day The 100th day The 150th day The 200th day The 250th day The 300th day

constraint Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP: Maximum STG: Minimum BHP:

6000 MMScf 7000 psia 4500 MMScf 5500 psia 4000 MMScf 3000 psia 3500 MMScf 2000 psia 3000 MMScf 1000 psia 2500 MMScf 500 psia 2000 MMScf 200 psia

First, we compare the gas production rate with and without considering the geomechanics effects on the fractures. In Figure 20a−c, we can clearly see that the gas production rate is much higher if the geomechanics effects are not considered. The reason is that it is impossible for the simulation to capture the closure of the fractures if the geomechanics effects are ignored. This phenomenon leads to a dramatic decrease in fracture permeability, so the gas production rate will sharply drop, which is accurately simulated only when the geomechanics effects are taken into consideration. The profile of gas saturation in natural fractures is shown in Figure 21a−c; the cumulative gas and water production is shown in Figure 22a,b; the gas production rate is shown in Figure 22c; the gas saturation at the grid block (43, 63, 1) in the natural fractures is shown in Figure 22d; and the mass transfer rate and the cumulative mass transfer between different media are shown in Figure 22e,f. From these figures, we can see the similar phenomena: (1) the gas saturation in the natural fractures at the grid blocks where the hydraulic fractures are located is higher than that in the surrounding grid blocks if the coupled models II and III are used; (2) the gas production 7769


Article

Energy & Fuels

Figure 20. (a−c) Gas production rate of the case with complex fracture network.

Figure 21. (a−c) Gas saturation profile in the case with a complex fracture network.

II, respectively. Although the predicted production by the coupled models I and II is close, their gas saturation profiles are significantly different; see Figure 24, which is similar to Figure 16a−c and 21a−c. From Figure 25a−b, we can see that, instead

and leads to a huge difference in the predicted gas production; see Figure 23c−f. Taking the results with the parameters VI and VII as examples, the coupled model III predicts 30% and 60% more cumulative gas production than the coupled models I and 7770


Article

Energy & Fuels

Figure 22. (a−f) Case with a complex fracture network.

coupled model I, the connection of matrix-HF is neglected, which may lead to much underestimated gas production. As we discussed, although the connection of matrix-HF is considered in the coupled model II, the matrix blocks disconnect with their neighboring matrix blocks, which means that the main pathway of gas is blocked. Hence the coupled model II also causes much underestimated production. As we can see in Figure 26, for the coupled model III, 65% and 80% of the gas in the hydraulic fractures comes from the connection of matrix-HF when the larger compressibility in parameters VI and VII is used. 5.2.2. Organic and Inorganic Matter Ratio. We use the organic/inorganic matter volume ratio shown in Table 6 to study the behavior of the coupled models with different volume percentages of organic and inorganic matter. The geomechanics parameters I are used in this test. Taking the cumulative gas production predicted by the coupled model I as the base value, we show the difference of production given by the coupled models II and III for the period of time from the 200th day to the 365th day in Figure 27. The coupled model II only has a minor difference from the coupled model I, while the coupled model III has a bigger difference on production with an increase in the content of inorganic matter in the matrix. That is because the fluid is allowed to flow through the inorganic

Table 5. Geomechanics Parameters Parameters II

Parameters III

Parameters IV

Parameters V

Parameters VI

Parameters VII

material

value

unit

matrix natural fracture hydraulic fracture matrix natural fracture hydraulic fracture matrix natural fracture hydraulic fracture matrix natural fracture hydraulic fracture matrix natural fracture hydraulic fracture matrix natural fracture hydraulic fracture

1e-7 3e-4 4.5e-5 1e-7 3e-4 6e-5 1e-7 4.5e-4 3e-5 1e-7 6e-4 3e-5 1e-7 4.5e-4 4.5e-5 1e-7 6e-4 6e-5

1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi 1/psi

of the connection of NF-HF, the connection of matrix-HF contributes the most gas to the hydraulic fractures. For the 7771


Article

Energy & Fuels

Figure 23. (a−f) Cumulative gas production in the cases with different geomechanics parameters.

5.2.4. Hydraulic Fracturing Stages. In this section, we examine the effects of the number of hydraulic fractures on the behavior of the coupled models. We place one, two, four, and six hydraulic fractures in the reservoir. The geomechanics parameters I are employed. The cumulative gas production from the 200th day to the 365th day is plotted in Figure 29. With an increase in the number of the hydraulic fractures, for the coupled model III, more gas is allowed to directly transfer from the shale matrix to the hydraulic fractures and then to the wellbore. Consequently, the production difference between the coupled models I and III increases from 8% to 12.5%. 5.2.5. Gas Adsorption/Desorption and Nonviscous Flow. In this section, we test the impacts of the gas adsorption/ desorption and nonviscous flow. The Langmuir volume VL reflects the capacity of gas supply from the adsorbed gas to the shale matrix, and therefore we set VL = 0, 0.08, 0.16, 0.32 ft3/ lbm and test its impacts. For the coupled model I, with the increasing of the Langmuir volume, more adsorbed gas can flow to shale matrix pores in an equal pressure condition, which means that more gas is capable of flowing into natural fractures and then into hydraulic fractures. As a result, the cumulative gas production predicted by the coupled model I becomes closer to that predicted by the coupled models II and III. We can see this phenomenon from Figure 30a. In Figure 30b, we compare the results with and without considering apparent permeability. When we consider the nonviscous gas flow, the permeability

medium by the coupled model III, and the connection of the neighboring sub-blocks for the inorganic matter and the connection of inorganic-HF are regarded as the main pathways for the gas flow. 5.2.3. Quantities of Matrix Nested Sub-Blocks. In order to determine the effects of the number of sub-blocks in the MINC model, we refine the matrix sub-blocks into (i) one organic subblock and one inorganic sub-block, (ii) two organic sub-blocks and two inorganic sub-blocks, (iii) two organic sub-blocks and three inorganic sub-blocks, (iv) three organic sub-blocks and two inorganic sub-blocks, and (v) three organic sub-blocks and three inorganic sub-blocks. The cumulative gas production is shown in Figure 28. From the result we can see that the difference in the cumulative gas production has a slight increase when we use more matrix sub-blocks. With an increase in the quantities of matrix sub-blocks, a higher resolution of pressure and saturation inside the shale matrix can be obtained. The high-resolution solutions may make the mechanisms discussed in the previous sections more obvious, and thus the error in production that is predicted by the coupled models I and II becomes larger. In our experience, 4−6 levels of matrix subblocks are recommended for the shale matrix and other ultratight matrix in order to capture reasonable physical phenomena, such as the transient flow behavior and water imbibition mechanism. 7772


Article

Energy & Fuels

Figure 24. (a−f) Gas saturation profiles in the cases with geomechanics parameters VI and VII.

Figure 25. (a, b) Gas transfer between hydraulic fractures and natural fractures, hydraulic fractures and organic matter, and hydraulic fractures and inorganic matter in the cases with geomechanics parameters VI and VII.

will increase during the depletion of pressure. The larger the permeability is, the more gas will diffuse into natural fractures. Therefore, there is more gas supply to natural fractures, and the result from the coupled model I will be closer to that from the coupled models II and III. 5.2.6. Behavior with Different Periods of Time. To study the behavior of the coupled models with different periods of time, we extend the simulation time to 1400 days and partition the time interval into three periods of time: (1) from the

beginning to the 300th day; (2) from the 300th day to the 700th day; and (3) from the 700th day to the 1400th day. The geomechanics parameters I are also used in this test. In Figure 31, we can observe that (1) the coupled model II has almost the same production as the coupled model I and (2) the coupled model III has about 10% more production than the other two models for each period of time. From this sensitivity test, it can be concluded that the difference in the production given by the coupled models does not reduce as time goes on. 7773


Article

Energy & Fuels

Figure 26. Percentage of gas transfer between hydraulic fractures and natural fractures, hydraulic fractures and organic matter, and hydraulic fractures and inorganic matter in the cases with geomechanics parameters VI and VII.

Table 6. Organic and Inorganic Volume Percentages Ratio I Ratio II Ratio III

organic matter (%)

inorganic material (%)

20 50 80

80 50 20

Figure 30. (a, b) Gas adsorption/desorption effects on the cumulative gas production difference from the coupled model I.

Figure 27. Cumulative gas production difference from the coupled model I with different organic/inorganic ratios. Figure 31. Cumulative gas production difference from the coupled model I with different periods of time.

6. CONCLUSIONS We have developed a comprehensive model, the coupled model III, for shale gas reservoirs. This model is capable of handling multiphase flow, complicated mechanisms of shale gas, geomechanics effects, and multiscale fractures. It is implemented in our in-house simulator. The coupled model III employs a MINC model to describe the shale matrix and natural fractures, and uses a discrete fracture model to explicitly represent large-scale hydraulic fractures. The mass transfer directly from organic and inorganic matter to hydraulic fractures and the mass transfer between the neighboring matrix blocks are considered by the model. According to the simulation results, it is essential to incorporate the geomechanics effects into fractured shale gas reservoir simulations, and a dramatic drop in the fracture permeability significantly affects the gas production. When the permeability of natural fractures is close to or smaller than that of matrix, rather than from the pathway of matrix-NF-HF, the gas directly from the matrix to hydraulic fractures makes a non-negligible contribution to the production. The coupled model III can capture the reasonable physical phenomena, and avoid underestimating the production as the coupled models I and II do. From the sensitivity studies, we can conclude that the difference in the production prediction by the coupled models increases when

Figure 28. Cumulative gas production difference from the coupled model I with different quantities of matrix sub-blocks.

Figure 29. Cumulative gas production difference from the coupled model I with different hydraulic fracture stages.

7774


Article

Energy & Fuels

homogenized media. SPE Reservoir Evaluation & Engineering 2008, 11, 750−758. (14) Moinfar, A.; Varavei, A.; Sepehrnoori, K.; Johns, R. T. Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional reservoirs. SPE Reservoir Simulation Symposium, 2013. (15) Moinfar, A.; Varavei, A.; Sepehrnoori, K.; Johns, R. T. Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs. SPE Journal 2014, 19, 289−303. (16) Jiang, J.; Younis, R. Hybrid Coupled Discrete-Fracture/Matrix and Multicontinuum Models for Unconventional-Reservoir Simulation. SPE J. 2015, 21, 1−009. (17) Jiang, J.; Younis, R. M. A multimechanistic multicontinuum model for simulating shale gas reservoir with complex fractured system. Fuel 2015, 161, 333−344. (18) Jiang, J.; Younis, R. M. Numerical study of complex fracture geometries for unconventional gas reservoirs using a discrete fracturematrix model. J. Nat. Gas Sci. Eng. 2015, 26, 1174−1186. (19) Zhang, R.-h.; Zhang, L.-h.; Wang, R.-h.; Zhao, Y.-l.; Huang, R. Simulation of a Multistage Fractured Horizontal Well with Finite Conductivity in Composite Shale Gas Reservoir through FiniteElement Method. Energy Fuels 2016, 30, 9036−9049. (20) Ding, D.; Farah, N.; Bourbiaux, B.; Wu, Y.; Mestiri, I. Simulation of Matrix-Fracture Interaction in Low-Permeability Fractured Unconventional Reservoirs. SPE Reservoir Simulation Conference, 2017. (21) Wang, L.; Wang, S.; Zhang, R.; Wang, C.; Xiong, Y.; Zheng, X.; Li, S.; Jin, K.; Rui, Z. Review of multi-scale and multi-physical simulation technologies for shale and tight gas reservoirs. J. Nat. Gas Sci. Eng. 2017, 37, 560−578. (22) Wu, K.; Li, X.; Wang, C.; Yu, W.; Chen, Z. Model for surface diffusion of adsorbed gas in nanopores of shale gas reservoirs. Ind. Eng. Chem. Res. 2015, 54, 3225−3236. (23) Wu, K.; Chen, Z.; Li, X.; Guo, C.; Wei, M. A model for multiple transport mechanisms through nanopores of shale gas reservoirs with real gas effect-adsorption-mechanic coupling. Int. J. Heat Mass Transfer 2016, 93, 408−426. (24) Chen, D.; Pan, Z.; Ye, Z. Dependence of gas shale fracture permeability on effective stress and reservoir pressure: model match and insights. Fuel 2015, 139, 383−392. (25) Wang, J.; Liu, H.; Wang, L.; Zhang, H.; Luo, H.; Gao, Y. Apparent permeability for gas transport in nanopores of organic shale reservoirs including multiple effects. Int. J. Coal Geol. 2015, 152, 50− 62. (26) Ambrose, R. J.; Hartman, R. C.; Diaz-Campos, M.; Akkutlu, I. Y.; Sondergeld, C. H. Shale gas-in-place calculations part I: new porescale considerations. SPE Journal 2012, 17, 219−229. (27) Gilman, J. R.; Kazemi, H. Improvements in simulation of naturally fractured reservoirs. SPEJ, Soc. Pet. Eng. J. 1983, 23, 695−707. (28) Wang, K.; Liu, H.; Chen, Z. A scalable parallel black oil simulator on distributed memory parallel computers. J. Comput. Phys. 2015, 301, 19−34. (29) Wang, K.; Liu, H.; Luo, J.; Chen, Z. A multi-continuum multiphase parallel simulator for large-scale conventional and unconventional reservoirs. J. Nat. Gas Sci. Eng. 2016, 33, 483−496. (30) Wang, K.; Liu, H.; Luo, J.; Yu, S.; Chen, Z.; Zhang, P. Parallel Simulation of Full-Field Polymer Flooding. Big Data Security on Cloud (BigDataSecurity), IEEE International Conference on High Performance and Smart Computing (HPSC), and IEEE International Conference on Intelligent Data and Security (IDS), 2016 IEEE 2nd International Conference, 2016; pp 220−225. (31) Liu, H.; Wang, K.; Chen, Z.; Jordan, K. E.; Luo, J.; Deng, H. A Parallel Framework for Reservoir Simulators on Distributed-memory Supercomputers. SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 2015. (32) Wallis, J.; Kendall, R.; Little, T. Constrained residual acceleration of conjugate residual methods. SPE Reservoir Simulation Symposium. 1985.

fractures become more sensitive to a stress change. An increase in inorganic matter in the shale matrix, the quantities of matrix sub-blocks, the number of hydraulic fracturing stages, and gas desorption also leads to a larger difference in the simulation results.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; [email protected]. ORCID

Kun Wang: 0000-0001-9442-5349 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The support of Department of Chemical and Petroleum Engineering, University of Calgary and Reservoir Simulation Group is gratefully acknowledged. The research is partly supported by NSERC/AIEES/Foundation CMG, AITF (iCore), IBM Thomas J. Watson Research Center, and the Frank and Sarah Meyer FCMG Collaboration Centre for Visualization and Simulation. The research is also enabled in part by support provided by WestGrid (www.westgrid.ca), SciNet (www.scinethpc.ca), and Compute Canada Calcul Canada (www.computecanada.ca).

■

REFERENCES

(1) Zhou, W.; Banerjee, R.; Poe, B. D.; Spath, J.; Thambynayagam, M. Semianalytical production simulation of complex hydraulic-fracture networks. SPE J. 2013, 19, 6−18. (2) Yu, W.; Wu, K.; Sepehrnoori, K.; Xu, W. A Comprehensive Model for Simulation of Gas Transport in Shale Formation With Complex Hydraulic-Fracture Geometry. SPE Reservoir Evaluation Eng. 2016.10.2118/178747-PA (3) Yang, R.; Huang, Z.; Yu, W.; Li, G.; Ren, W.; Zuo, L.; Tan, X.; Sepehrnoori, K.; Tian, S.; Sheng, M. A Comprehensive Model for Real Gas Transport in Shale Formations with Complex Non-planar Fracture Networks. Sci. Rep. 2016, 6.10.1038/srep36673 (4) Warren, J.; Root, P. J. The behavior of naturally fractured reservoirs. SPEJ, Soc. Pet. Eng. J. 1963, 3, 245−255. (5) Pruess, K. A practical method for modeling fluid and heat flow in fractured porous media. SPEJ, Soc. Pet. Eng. J. 1985, 25, 14−26. (6) Wu, Y.-S.; Pruess, K. A multiple-porosity method for simulation of naturally fractured petroleum reservoirs. SPE Reservoir Eng. 1988, 3, 327−336. (7) Wu, Y.-S.; Di, Y.; Kang, Z.; Fakcharoenphol, P. A multiplecontinuum model for simulating single-phase and multiphase flow in naturally fractured vuggy reservoirs. J. Pet. Sci. Eng. 2011, 78, 13−22. (8) Yan, B.; Wang, Y.; Killough, J. E. Beyond dual-porosity modeling for the simulation of complex flow mechanisms in shale reservoirs. Computational Geosciences 2016, 20, 69−91. (9) Wu, Y.-S.; Li, J.; Ding, D.; Wang, C.; Di, Y. A generalized framework model for the simulation of gas production in unconventional gas reservoirs. SPE Journal 2014, 19, 845−857. (10) Karimi-Fard, M.; Firoozabadi, A. Numerical simulation of water injection in 2D fractured media using discrete-fracture model. SPE Annual Technical Conference and Exhibition, 2001. (11) Karimi-Fard, M.; Durlofsky, L. J.; Aziz, K. An efficient discrete fracture model applicable for general purpose reservoir simulators. SPE Reservoir Simulation Symposium, 2003. (12) Xu, C.; Li, P.; Lu, D. Production performance of horizontal wells with dendritic-like hydraulic fractures in tight gas reservoirs. J. Pet. Sci. Eng. 2017, 148, 64−72. (13) Li, L.; Lee, S. H. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and 7775


Article

Energy & Fuels (33) Cao, H.; Tchelepi, H. A.; Wallis, J. R.; Yardumian, H. E. Parallel scalable unstructured CPR-type linear solver for reservoir simulation. SPE Annual Technical Conference and Exhibition, 2005. (34) Liu, H.; Wang, K.; Chen, Z. A family of constrained pressure residual preconditioners for parallel reservoir simulations. Numerical Linear Algebra with Applications 2016, 23, 120−146. (35) Liu, H.; Wang, K.; Chen, Z.; Jordan, K. E. Efficient multi-stage preconditioners for highly heterogeneous reservoir simulations on parallel distributed systems. SPE Reservoir Simulation Symposium, 2015. (36) Cui, T.; Wang, K.; Liu, H.; Luo, J.; Chen, Z. Efficient and Scalable Methods for Dual-Porosity/Permeability Models in Fractured Reservoir Simulations. SPE Asia Pacific Oil & Gas Conference and Exhibition, 2016. (37) Luo, J.; Chen, Z.; Wang, K.; Deng, H.; Liu, H. An Efficient and Parallel Scalable Geomechanics Simulator for Reservoir Simulation. SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 2015. (38) Luo, J.; Wang, K.; Liu, H.; Chen, Z. Coupled Geomechanics and Fluid Flow Modeling in Naturally Fractured Reservoirs. Low Permeability Symposium, 2016.

7776


A Comprehensive Model Coupling Embedded ... - ACS Publications

Recommend Documents