Sage Journals: Discover world-class research

Abstract

Due to the discontinuity and multi-scale of fracture porous media (FPM), fluid transport in the naturally fractured reservoir (NFR) is more complex than in the conventional reservoirs. Based on the embedded discrete fracture model (EDFM), the effect of fracture pattern (i.e. fracture orientation, intersection relationship, numbers and distribution) on fluid transport is investigated, and the paper also conducts research on pressure distribution and flow behavior. A series of fracture models are built to get a detailed analysis of fracture pattern influences. Results show that when the fracture orientation tends to parallel the main stream direction, the seepage flow through the fracture is more distinct and the pressure drop at two sides of fracture becomes larger. As the intersecting fracture cannot serve as part of main-stream channels, the intersection is almost invalid. In the multi-scale fracture model, the pressure drop in the large-scale fractures network is higher than in the small-scale fractures. Thus, the simulation results show that the fracture pattern has a significant impact on the seepage flow and pressure distribution. Also, this model is simulated for M field which has random natural fracture network to study the flow and transport characteristics in fractured porous media. The comprehensive simulation results the fracture patterns have impact effect on production.

Keywords

Naturally fractured reservoir fluid transport embedded discrete fracture model fracture pattern production prediction

Introduction

Fractured reservoir represents over 40% of the world’s oil and gas reserves (Allan and Sun, 2003). Unlike the conventional reservoir, the fractured reservoir exists extensive fracture networks, with multi-scales and orientated fractures in a specific region, where permeability anisotropy is high. Therefore, the numerical simulation for fractured reservoir is still a challenge. There are extensive studies over the past twenty years, which have proposed different modeling approaches and numerical strategies to simulate different types of fractures.

In 1960s, it was proposed the dual-porosity model and its extensions for naturally fractured porous media with many naturally small-scale fractures (Kazemi, 1969; Warren and Root, 1963). This method is only suitable for small-scale fractures in the reservoir. Thus, it fails to provide accurate depictions and simulations for various scale and arbitrary orientation fractures.

To obtain more accurate simulation for the multi-scale fractures reservoir, the discrete fracture model (DFM) approach was devised (Hoteit and Firoozabadi, 2008; Karimi-Fard and Firoozabadi, 2003; Noorishad and Mehran, 1982). In this method, geometry and locations of fractures are captured accurately by using complex unstructured gridding techniques. However, numbers of small grids is required to capture fracture geometries which makes the computational cost high. Therefore, the DFM has limited applicability for realistic simulation due to the complex grids.

A new approach EDFM is introduced and the schematic as shown in Figure 1 which can provide accurate solutions for multi-scale fractures reservoir simulation (Hajibeygi et al., 2011; Lee et al., 2001; Li and Lee, 2008). In this model, small-scale fractures are homogenized and treated as a continuous and damaged matrix with effective coefficients. On the other side, large-scale fractures are represented by a coupled discrete fracture model. Moreover, the EFDM model has less degrees of freedom than the DFM, thus, it has sufficient coupling to achieve good convergence rates.

Figure 1.

2 D matrix grids with 1 D fracture network of 5 fracture lines connection (The shaded cells are embedded fracture.)

The aim of this research is to study the influence of fracture pattern on fluid transport in the FPM using the EDFM and solved by the FVM. Several fracture models with different pattern (i.e. fracture orientation, intersection, numbers and distribution) are established to investigate the seepage flow and pressure distribution in the FPM. The preliminary simulation results show that the fracture with different patterns have a significant impact on fluid transport in fractured reservoir.

Mathematical model in EDFM

The EDFM uses a structured grid to represent matrix system and introduces additional fractures by computing intersection of fractures with the matrix grids as shown in Figure 1. The modeling approach is organized by governing equations, fracture-matrix transfer equations, fracture-fracture intersection equations and fully coupled simulation process.

Governing equations

The flow in the FPM is often described by the Darcy-law equation for incompressible single-phase flow. We consider the Darcy’s law for incompressible single-phase flow without capillary pressure and gravity effects in this study.

Matrix mass balance equation is written as:

- \nabla \cdot (\frac{K_{m}}{μ} \cdot \nabla p_{m}) + ψ_{m f} = q_{m}

(1)

The fluid velocity can be obtained using Darcy-law:

v_{m} = - \frac{K_{m}}{μ} \cdot \nabla p_{m}

The transport equation in the matrix system is derived to the continuity equation by balancing the components of the transported mass, as:

- \nabla (κ \nabla c) + v_{m} \cdot \nabla c = ϕ_{m} \frac{\partial c}{\partial t}

(2)

The fracture mass balance equation is given by:

- \nabla \cdot (\frac{K_{f}}{μ} \cdot \nabla p_{f}) + ψ_{f m} + T_{f f} = q_{f}

(3)

The fracture permeability is assumed by the cubic law which relates to fracture aperture:

K_{f} = \frac{b_{f}^{2}}{12}

(4)

Here, the transport equation in fracture system is similar to the transport equation (2) of the matrix system.

Fracture-matrix coupling equations

Transport function is treated as a source term in the mass balance equation between matrix and fracture, which is similar as in classical well models (Peaceman, 1978):

ψ_{m f} = C I \cdot \frac{K_{m}}{μ} (p_{f} - p_{m})

(5)

The connectivity index CI is computed based on the position and the properties of matrix grids and fracture segments. It is defined as the area of fracture segment bounded inside grid block over the average distance between matrix and fracture control volumes (Li and Lee, 2008):

C I = \frac{A_{m - f}}{d_{m - f}}

(6)

d_{m - f} = \frac{\int_{V} n \cdot xdv}{V}

(7)

Fracture intersection equations

In fact, fractures often intersect other fractures which has a big impact on the flow dynamics in the FPM. In this model, the effect of fracture-fracture coupling was considered. The transmissivity between intersection fractures can be obtained by a similar approach used in electrical engineering for known as star-delta transformation in circuits (Karimi-Fard and Firoozabadi, 2003):

T_{f f} = \frac{T_{i} \cdot T_{j}}{T_{i} + T_{j}}

(8)

T_{i (j)} = \frac{A_{f} λ_{i (j)}}{0.5 \cdot d x_{f}}

(9)

Discretizing the pressure and transport equations

The FVM is used to discretize the pressure and equations for matrix and fracture system. The domain Ω for elliptic pressure equation can be discretized as the integration over finite control volume Ω _ij . The divergence integral for matrix system can be written as:

- \int_{Ω_{i j}} [\nabla \cdot (λ_{m} \cdot \nabla p_{m}) + ψ_{m f}] d V = - \int_{Ω_{i j}} q_{m} d V

(10)

The divergence integral also can be expressed as the surface integral to the boundary of the volume by using Gauss’s theorem:

- \int_{Ω_{i j}} [\nabla \cdot (λ_{m} \cdot \nabla p_{m}) + ψ_{m f}] \cdot ndS = - \int_{Ω_{i j}} q_{m} d V

(11)

The pressure gradient over the cell boundary is approximated by two points flux approximation scheme with the FVM. The Figure 2 shows the schematic of matrix and fracture grid. The discretized expression for the matrix pressure is given by:

\begin{array}{l} \frac{Δ A_{x} λ_{i - \frac{1}{2}, j}}{Δ x} (p_{i, j} - p_{i - 1, j}) + \frac{Δ A_{x} λ_{i + \frac{1}{2}, j}}{Δ x} (p_{i, j} - p_{i + 1, j}) + \\ \frac{Δ A_{y} λ_{i, j - \frac{1}{2}}}{Δ y} (p_{i, j} - p_{i, j - 1}) + \frac{Δ A_{y} λ_{i, j + \frac{1}{2}}}{Δ y} (p_{i, j} - p_{i, j + 1}) = q_{i, j} + \sum_{Ω_{i j} \cap Ω_{k}} C I_{k} λ_{i, j} (p_{k} - p_{i, j}) \end{array}

(12)

Figure 2.

The schematic of matrix grid cell and fracture grid. (The blue line is a fracture. The shaded cells are embedded fracture.).

Similarly, the pressure for the fracture system can be expressed as:

\frac{b_{f} λ_{k - \frac{1}{2}}}{Δ x_{f}} (p_{k} - p_{k - 1}) + \frac{b_{f} λ_{k + \frac{1}{2}}}{Δ x_{f}} (p_{k} - p_{k + 1}) = q_{k} - \sum_{Ω_{i j} \cap Ω_{k}} C I_{k} λ_{i j} (p_{k} - p_{i j})

(13)

Here, the discretized approach for transport equation is similar to the discretization of pressure equation.

The solving strategy

The FV discretization scheme is used to solve for pressure and transport equations, so they can be assembled respectively as a linear system. The block matrix expressions for both of pressure and transport systems are given by:

(\begin{matrix} T_{m m} & T_{m f} \\ T_{f m} & T_{f f} \end{matrix}) \cdot (\begin{matrix} p_{m} \\ p_{f} \end{matrix}) = (\begin{matrix} Q_{m} \\ Q_{f} \end{matrix})

(14)

(\begin{matrix} A_{m m} & A_{m f} \\ A_{f m} & A_{f f} \end{matrix}) \cdot (\begin{matrix} c_{m} \\ c_{f} \end{matrix}) = (\begin{matrix} M_{m} \\ M_{f} \end{matrix})

(15)

Figure 3 shows a flowchart of calculation procedures. At the initial time (T=t₀), the pressure is defined as p₀, then start calculating the pressure at time step n. next, the pressure solution is used to solve velocity. Then the concentration at time step n is calculated. If the time step reaches the maximum value, the calculation is end and the results are output; if not, the calculation will be continuing.

Figure 3.

The flowchart of calculation procedures.

Model validation and introduction

2D geometry of two vertical crossing fractures in base model is considered as shown in Figure 4. The domain of base case is 9 by 9 m and discretized by 99 by 99 cells. The pressure at the inlet (left-hand-side) is 10 × 10⁵ Pa and the outlet pressure (right-hand-side) is zero Pa. The top and bottom boundaries are no flow boundaries. A concentration boundary of 1 kg/m³ is applied at the inlet. The base model parameters are shown in Table 1. The new model is compared with the reference model of DFM (Pluimers, 2015). Figure 5 shows that the comparison of pressure in two systems between the models and the agreement is quite good, despite the different fracture apertures lead to a slight deviation in the comparison. The error of comparison is less than 4%, which indicates the new EDFM is validated. Besides, the new EDFM contains two order of magnitude less degree of freedom (DOF) than the reference model.

Figure 4.

The base model (black lines are fractures).

Table 1.

Parameters of the verified model.

Variable	Description	Value
φ_m (%)	Matrix porosity	0.3
K_m (mD)	Matrix permeability	1.0
K_f (mD)	Fracture permeability	1 × 10³
μ (cP)	viscosity	1.0
P_in (Pa)	Inlet pressure	10 × 10⁵
P_out (Pa)	Outlet pressure	0.0
c (kg/m³)	Inlet concentration	1.0

Figure 5.

The comparison of pressure between the EDFM and the reference model. (a) Matrix pressure system. (b) Fracture pressure system.

Results and discussion

Based on the previous verified model, the effect of fracture patterns (i.e. Fracture orientation, intersection, numbers and distribution) are systematically studied.

Effect of fracture orientation

Figures 6 and 7 show the fluid transport and pressure distribution with different fracture orientations. The effect of single fracture orientation (0, 45 and 90 degrees) on seepage field is studied. The changes of fracture orientation have a big influence on fluid transport and pressure distribution. When the fracture orientation is 90 degree, the fracture is perpendicular to the main streamline, the fluid gathers the left side of fracture and flow through the short distance of the fracture then getting into matrix. The distribution of pressure field is more homogenous, which indicates that a fracture perpendicular to the mainstream has less contribution to the fluid transport. While the fracture orientation tends to parallel the mainstream line direction (45 degree), the seepage flow through fracture is more distinct and the pressure drop at two sides of fracture becomes larger. As the fracture is parallel to main streamline direction (0 degree), the fracture will contribute largest conductivity the fluid transport.

Figure 6.

Pressure distribution with different fracture orientations.

Figure 7.

Fluid transport with different fracture orientations.

Effect of fracture intersection

Figures 8 and 9 show that the effect of fracture intersections on fluid transport and pressure distribution. The interference of fracture intersection on pressure field and seepage flow is studied. The fluid transport with different fracture intersections is described as follows. When the two perpendicular intersected fractures are in the middle of the domain, only one fracture which is parallel to streamline has a contribution to the fluid flow and the pressure conduction. It illustrates that the intersection is invalid if one intersecting fracture not served as a part of main-stream channels. As fractures rotate 45 degree intersection in the domain, the pressure distribution becomes more heterogeneous and each fracture begins to contribute the flowing of fluid. With the intersection closes to 15 degree, both fractures gradually parallel to main streamline, the more fluid transport is, the faster pressure conduction rate is.

Figure 8.

Pressure distribution with different fracture intersection degrees.

Figure 9.

Fluid transport with different fracture intersection degrees.

Multi-scale fractures simulation

Based on the previous verified model, a numerical simulation of the multi-scale complex fracture model is studied. The number of fractures increases to 10 and the domain in a 500 by 500 m domain. The inlet pressure increases to 50 × 10⁵ Pa. The fractures are distributed in arbitrary angles and different scales as shown in Figure 10.

Figure 10.

The simulation study of the multi-scale fracture model (500 by 500 m domain): (a) Pressure distribution and (b) Fluid transport.

Figure 10(b) shows that the flow front distributes unevenly, which is caused by the permeability contrast between fracture and matrix. The higher permeable fracture and the effective intersection (which tend to main streamline direction) lead to larger pressure drops and higher flow rate. This is happened because the fracture allows the fluid flow through the domain earlier than the matrix. When the longer scale fracture penetrates into matrix, more fluid gathers into fractures and the flow is enhanced more than the smaller scale fracture. Also, the pressure distribution is more heterogeneity at the two ends of the larger scale fracture as shown in Figure 10(a).

The fracture pattern effect on production in M field

M field, located in West Africa, is a granite fractured reservoir. The net pay is 100 m, drainage radius is 50 m. Fluid viscosity is 4.66 cp. Matrix permeability is about 10∼50 mD and porosity is 3∼11%. Fracture density is 0.7∼2 per meter. Fracture width is 0.16∼0.22 mm. fracture orientation ranges from 20 to 75 degree. The well model and complex fracture network as shown in Figure 11(a) and pressure distribution as shown in Figure 11(b). In this model, the fluid is compressible. Single phase flow. Darcy flow in the matrix, while flow in the fracture system is related to fracture width, so fracture permeability changes accordingly. Fluid flow into the wellbore only through fractures while the matrix contribution is not considered. The simulation result is verified by actual production data, the relative error is less than 5%, as shown in Figure 12.

Figure 11.

The simulation for a well in M field: (a) Geometry of discrete fracture network near the wellbore and (b) Pressure distribution near the wellbore.

Figure 12.

The simulation results related to the comparison of EDFM, D-P model and practical production data.

The effect of the fracture orientation

Figure 13 shows the effect of fracture orientation on well production. As the angle between fracture and wellbore increases, the higher the well production is. When the fracture is perpendicular to the wellbore, which fracture orientation is zero degree, the max production rate will be achieved. This is due to the direction of fluid flow is that of maximum fracture permeability.

Figure 13.

The fracture orientation effect on production in M field.

The effect of the fracture numbers

Figure 14 indicates the effect of fracture numbers on well production. As fracture numbers increase, the well production increases accordingly. Since more fractures were drilled and communicated with wellbore, the flow resistance reduced and pressure sweep area increased so that the production enhanced. With fracture increase continuously, the well production increases slowly. This is due to fracture spacing becoming more compact, the increase of flow interference between fractures.

Figure 14.

The fracture numbers effect on production in M field.

The effect of the fracture permeability

Figure 15 shows the effect of fracture permeability on well production. As fracture permeability increases, more fluid flows into the wellbore, then the production rate increased by 40%. Furthermore, if the fracture permeability continues to increase from 800 to 1000 mD, the production increase has slowed down. With the fracture permeability increase to 1500 mD, the production rate reaches a certain value. This is due to the capacity of fluid flowing from the matrix to the fracture is limited at the location of the fracture reservoir.

Figure 15.

The fracture permeability effect on production in M field.

The effect of fracture to matrix permeability ratio

From the vertical zone in M field, the reservoir divided into weathering and fracture zone. When drilling in different reservoir, the production rates are varying. Based on fracture permeability ranges in 10∼50 mD and matrix permeability are 500 ∼1500 mD in vertical zone. The effect of fracture to matrix permeability ratio on well production is studied, as shown in Figure 16.

Figure 16.

The ratio of fracture to matrix permeability effect on production in M field.

As fracture to matrix permeability ratio is 10 times, the permeability difference between fracture and matrix is little, permeability and pressure distribution near the wellbore is homogenous, the reservoir with good communication, so as to obtain high well production. When the ratio of fracture permeability to matrix permeability gradually increases to 30 times, the more uneven the pressure distribution near the wellbore, the faster the pressure decline, lower well production achieved. Furthermore, if the difference between fracture permeability and matrix permeability increases to 100 times, indicating the fractures are main flow channel. Pressure distribution near the wellbore is extremely uneven, the well pressure and production declined dramatically.

Conclusions

In this study, a new productivity model based on EDFM has been established and the mathematic seepage flow model of both matrix and fracture systems are given. The implementation of new EDFM is validated with the reference model (Pluimers, 2015) and a good agreement is observed. The new model is more efficient than the reference model.

When fracture orientation tends to parallel the mainstream line, the pressure drop at two sides of fracture becomes larger and more fluid transport through the fracture. If the intersecting fracture cannot serve as part of main-stream channels, the intersection is meaningless.

The simulation of multi-scale fracture model indicates that the heterogeneity of the flow front happens due to the large permeability contrast between fracture and matrix. The pressure drop in large-scale fractures network is higher than in the small-scale fractures. Thus, the fracture pattern has a big impact on fluid transport in the FPM.

Several simulation cases are established for production study in M field to study the fractures effect on production. The results show that with increases of fracture to matrix permeability ratio, production decreases and pressure distributes unevenly; With fracture numbers and permeability increase, the growth of production is limited.

The study has a potential limitation. The current implementation of the EDFM uses a first order upwind approximation for the advection term, which introduces some artificial numerical diffusion. Multiple approaches will be used to treat the advection term of transport equation.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

We thank the support from Research Institute of Petroleum Exploration & Development, CNPC.

ORCID iD

Jinlin Zhang

References

Allan

Sun

(2003) Controls on recovery factor in fractured reservoirs: Lessons learned from 100 fractured fields. In: The SPE annual technical conference and exhibition. Society of Petroleum Engineers; SPE-84590-MS.

Hajibeygi

Karvounis

Jenny

(2011) Hierarchical fracture model for the iterative multi-scale finite volume method. Journal of Computational Physics 230(24): 8729–8743.

He Y, Cheng S, Sun Z, et al. (2020) Improving oil recovery through fracture injection and production of multiple fractured horizontal wells. Journal of Energy Resources Technology 142(5): 053002.

Hoteit

Firoozabadi

(2008) An efficient numerical model for incompressible two-phase flow in fractured media. Advances in Water Resources 31(6): 891–905.

Karimi-Fard

Firoozabadi

(2003) Numerical simulation of water injection in fractured media using the discrete fracture model and the Galerkin method. SPE Reservoir Evaluation & Engineering 6(2): 117–126.

Kazemi

(1969) Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution. Society of Petroleum Engineers Journal 9(4): 451–462.

Lee

Lough

Jensen

(2001) Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research 37(3): 443–455.

Lee

(2008) Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation & Engineering 11(4): 750–758.

Nandlal K and Weijermars R (2019) Drained rock volume around hydraulic fractures in porous media: Planar fractures versus fractal networks. Petroleum Science 16(11): 1064–1085.

10.

Noorishad

Mehran

(1982) An upstream finite element method for solution of transient transport equation in fractured porous media. Water Resources Research 18(3): 588–596.

11.

Peaceman

(1978) Interpretation of well-block pressures in numerical reservoir simulation. Society of Petroleum Engineers Journal 18(3): 183–194.

12.

Pluimers

(2015) Hierarchical Fracture Modeling Approach. Netherlands: TU Delft.

13.

Warren

Root

(1963) The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal 3(3): 245–255.