Finite Element Model to Simulate Two-Phase Fluid Flow in Naturally Fractured Oil Reservoirs: Part I.

Naturally fractured reservoirs host more than 20% of the world's total oil and gas reserves. To produce from such reservoirs efficiently, a good understanding of the reservoir behavior at various conditions is essential. This allows us to predict the reservoir performance in advance and assess its economic feasibility. However, the production from such reservoirs is challenging due to (a) uncertainty associated with the fracture map, (b) complex physics phenomena of fluid and rock interaction, and (c) lack of comprehensive knowledge of the extent, orientation, and permeability sensitivity of the fracture network. This paper addresses the abovementioned challenges by presenting a three-dimenisonal (3-D) two-phase fluid flow model in a poroelastic environment. The model is based on a hybrid methodology by combining single continuum and discrete fracture network approaches. Also, the capillary pressure effect, saturation, and relative permeability variations are considered. A mathematical formulation for three-dimensional, two immiscible fluid flows including rock deformation for the fracture network and the rock matrix is presented. A standard Galerkin-based finite element method is applied to discretize the poroelastic governing equations in space and time. The characteristic Galerkin discretization method is used to stabilize the solution of the convection equation in a finite element approach. The 3-D model is validated against IMEX commercial software and finite element package to test its robustness. The results show that the developed model has the ability to predict the two-phase flow behavior precisely, which can be used to assess the production performance of naturally fractured reservoirs. Numerical results of fluid flow profiles for single and multiple fractures with different orientations that match experimental oil drainage tests and field case studies are presented that proves the reliability of the developed multiphase flow model. Results of simulated well production data show an excellent match with the field production data.

Introduction

Multiphase fluid flow simulation in a three-dimenisonal (3-D) naturally fractured reservoir to estimate the production potential remains a notable challenge. The key reason is the complex subsurface fracture network that affects the oil production process and the extremely large amount of produced water which leads to the damage of oil zones.[1] In this study, equations of two-phase flow are derived and presented for the sake of understanding of fluid flow behavior. Different models are used for the simulation of immiscible flow through such fractured reservoirs. These models used the dual porosity/dual permeability and discrete fracture approaches. A discrete fracture network (DFN) approach integrates different mathematical approaches such as finite volume and Galerkin’s finite element.[2,3] Hotteit and Firoozabadi[4] (2006) developed a compositional simulator that used the DFN approach and mixed finite element for two-phase flow simulation in naturally fractured reservoirs in a two-dimensional (2-D) space matrix and one-dimensional (1-D) space fracture. Mathai et al. (2005)[5] developed a 3-D simulation model for flow simulation through fractured porous media using the CVFE method. Xu et al. (2005)[6] extended Mathai et al. (2005)[5] for multiphase flow simulation. A mixed finite element (MFE)[7,8] technique is used to estimate the velocity fields in highly heterogeneous especially for single-phase fluid flow.[9] Hoteit and Firoozabadi[10] extended MFE for multiphase flow taking into account during the simulation process and gravity and capillarity effects. Azim et al.[11,12] used the concept of a discrete fracture approach to generate the subsurface fracture map and for relative permeability upscaling from the core scale to the field scale. In addition, Azim et al.[13] created a 3-D numerical model to evaluate water conning phenomena to estimate the production potential of fractured basement reservoirs. Azim et al. (2017)[14] developed a 3-D fractured model to study the low salinity waterflooding in naturally fractured oil reservoirs. The model is based on a discrete fracture approach and finite element technique. The fluid flow through long fractures is simulated using a cubic law, assuming that the fracture is a continuous porous media without any embedded granular materials.

The fluid flow equations for fractured porous media are solved using analytical and numerical methods, e.g., finite element method;[15] finite volume method,[16] mixed finite element,[17] and boundary element method.[18] These approaches have several drawbacks including the long computational time, and a hybrid methodology of incorporating the permeability tensors is used to avoid such limitations.[19] Park et al. (2000)[20] simulated the flow through the discrete fracture approach and ignored the flow at the fracture matrix interface during the calculation of permeability tensors in 3-D space. Gupta et al. (2001)[21] and Pride et al. (2003)[22] ignored the effect of the matrix in permeability tensor estimation and flow simulation. Fluid production from a reservoir or injection into a reservoir alters the pressure state and leads to deformation of the rock by generating seismic activities. Furthermore, porosity, permeability, and oil recovery are consequently affected. For the above reasons, this study focuses on presenting a novel poroelastic numerical model to assess the oil recovery of fractured reservoirs. The model is for two-phase flow under geomechanics effects to comprise the dynamic behavior of the fracture system. The simulation workflow is based on the upstream flux weighted technique. The water saturation equation is discretized using the standard finite element method to get a stable solution.

To overcome such phenomena, characteristics of the Galerkin’s discretization method is employed to stabilize the equation solution in a finite element approach. Synthetic cases using a commercial black oil numerical reservoir simulator (CMG-IMEX) for the homogenous and fractures system are created to check the reliability of the developed multiphase numerical model. Moreover, the aim of the developed multiphase flow numerical model is to understand the influence of the fractures on a multiphase fluid flow at the fracture matrix interactions. For that purpose, the current study used a two-phase flow experiments of fractured porous media.

In this study, the workflow for the simulation of two phases is based on a hybrid methodology that combines single continuum and DFN approaches. The workflow of this study is to show the derivation of two-phase flow equations in a poroelastic environment. Next, the finite element approach is used to discretize flow equations in a three-dimensional space. In addition, this study presents a developed solution for the water saturation equation in space and time. The model is validated against synthetic cases created using a commercial numerical reservoir simulator (CMG-IMEX) and a finite element toolbox.

Derivation of Multiphase Flow Equations

Mass Conversation Equation

The mass of the medium solid component can be expressed aswhere ϕ is the porosity and ρs is the density of the fluid.

The mass conversation is given byEquation is rearranged asThe mass of the medium is given byThe fluid mass conversation iswhere Sψ is the saturation, ρψ is the density, and qψ is the rate exchange between the matrix and the fracture system.

Equation is discretely written for oil and water fluid phases aswhere U the intrinsic velocity.where k1 is the matrix permeability and p1ψ and p2ψ are the matrix and fracture pressures, respectively.

Darcy velocities are defined asFrom eq 9, intrinsic velocities can be described asSubstituting eq into eq Expanding the derivatives of eqs and 11Equations and 14 can be reformulated and written as followsUsing the total derivative ofEquations and 16 can be rewritten using eq asFrom eq , the porosity changing with respect to time can be given asSubstituting eq into eq orwhere K is the fluid bulk modulus and Ks and Kns are the bulk moduli of solid rock.

Substituting eq into eq

The two-phase fluid flow equation is given as follows:

Combine eq with Darcy’s lawThe water-phase flow governing equation will be as followsThe oil-phase fluid flow governing equation can be written aswhere Pwm and Pwf are the water pressures inside the matrix and fractures.

Momentum Balance

The relationship between applied stresses σ and intergranular (effective) stresses σ is given bywhere α is the Biot’s constant, and δ is the Kronecker delta.

The linear relationship iswhere D is the elasticity matrix.

The motion equation for a solid iswhere F is the traction vector. The relationship between the body strain and its displacement is defined asA Navier-type equation for the displacement u is given aswhere G is the shear modulus, v is the Poisson’s ratio, and is the average pore pressure.

Finite Element Discretization

The 3-D four node tetrahedral element (see Figure ) and the shape functions are defined as followswhere (ξ, η, ζ) are the element coordinate system. The element geometry fan be transformed into the local coordinates (ξ, η, ζ) as followswhere is the Jacobin matrix and can be defined aswhere N is the corresponding shape function, and u̅, P̅w, and P̅o are the nodal unknown variables. is the strain displacement matrix and can be defined asThe equilibrium equation iswhere ∂f is the load.

Figure 1

Tetrahedral element used in the generated mesh for simulation.

Substituting eq into eq and by dividing by ∂t, the equation will be as followsThe average pore pressure (P) can be defined asSubstituting eq into eq , we getEquation is modified by incorporating the capillary pressure as followsDiscretization form for the water phase is given as followsThe 2-D space flow iswhere Ω̅f represents the fracture part of the domain as a 2-D entity, and Ωm represents the matrix domain and Ω is the entire domain.

Discretization in Time

The coupled equations are discretized in time using the fully implicit scheme as followswhere Δt is the time step size.

In this study, an iterative approach is used for nodal unknown estimation at each time step. The iterative technique used in this study to calculate the nonlinear coefficient includes capillary pressure and fluid saturation and relative permeability. A convergence criterion is used for eqs and 45 to attain a stable solution of fluid pressure and displacement aswhere Ri is the number of nodal unknowns, and ε is the convergence limit.

Model Validation

The model is validated in an uncoupled way, in which the elastic and two-phase flow problems are validated separately. For the elastic problem validation, one patch test is used that includes a simple fixed traction test. However, the flow problem is tested and validated against commercial black oil reservoir simulators (CMG-IMEX).

Model Validation Using Elasticity Problems

The test is performed under fixed traction. A 3-D tetrahedral element model with dimensions of 5 m x 5 m x 5 m is used (see Figure ). The model’s Young’s modulus and Poisson’s ratio are 200 GPa and 0.3, respectively. A fixed traction of 500 N/m2 is applied on the top face of the cube and a fixed displacement on the bottom face. The following analytical solution is used to assess the accuracy of the developed elastic model resultswhere u is the displacement, z is the model height, E is the Young’s modulus, and σ is the applied stress.

Figure 2

3-D tetrahedral model used to validate the elastic model of the designed 3-D simulator.

From eq 50, it can be drawn that the results of the displacement due to the application of uniform traction is linear in the z-direction with zero value at the bottom surface and the maximum value of 1.25E-8 m at the top surface. The developed elastic numerical model provided a displacement of 9.5E-09 m at the top surface as shown in Figure . The results of displacement from the analytical solution and the developed elastic model are very close with error of 14%, and it can be improved by increasing the number of tetrahedral elements used in the tested model.

Poroelastic Model Validation

A 2-D circular reservoir with a drainage radius of 1000 and 0.1 m wellbore radius is used (see Figure ) to validate the poroelastic model developed in this study. The drained condition is used, which is obtained using the Kirsch’s[23] problem. The simulation results for the poroelastic model are plotted against the analytical solutions in Figures –6. The results of pressure, tangential stress, displacements, effective stresses, and stress profiles presented in those figures are in good agreement with analytical solutions. The data[24] used in this validation section are presented in Table .

Figure 3

2-D circular reservoir for poroelasticity.

Figure 4

Pore pressure distribution (a) with time and (b) after 1 h of injection under σH = 5800 psi and σh = 5500 psi, Pr = 5500 psi, Pw = 1000 psi, K = 0.01 md, Ky = 0.01 md.

Figure 6

Radial stress along x axis distribution under σH = 5800 psi and σh = 5500 psi, Pr = 5500 psi, Pw = 1000 psi, K = 0.01 md, K = 0.01 md.

Table 1

Parameters Used in the Verification of Poroelasticity Solutions

E-modulus	40 GPa
poisson’s ratio	0.2
ϕ	0.1
Cw	1.0E-4 psi–1
μw	0.1 cp
Biot’s coefficient	1.0
max stress	5800 psi
min stress	5500 psi
initial reservoir pressure	5500 psi
wellbore pressure	1000 psi
formation permeability Kx	0.01 md
formation permeability Ky	0.01 md
rw	0.1 m
re	1000 m

Tangential stress distribution under σH = 5800 psi and σh = 5500 psi, Pr = 5500 psi, Pw = 1000 psi, K = 0.01 md, Ky = 0.01 md.

Validation of the Unfractured Model Using CMG Software

Validation of two-phase fluid flow is performed by creating a 2-D mesh reservoir of 10,000 elements with one injector and one producer as shown in Figure . Minimum horizontal stress of 5000 psi is applied in the y-direction while maximum horizontal stress of 5500 psi is applied in the x-direction. The reservoir pressure initially is 4500 psi. The permeability and porosity of the model are 100 md and 0.2, respectively. The 2-D model is created as well using CMG commercial software. The same inputs are used to compare the water saturation profile between the two models.

Figure 7

2-D simulation model used in the validation process with σmin = 5000 psi and σmax = 5500 psi.

As it can be seen from Figure a,b, the water saturation profile is propagated after 5 min of water injection (Figure a) and 1 h of injection (Figure b). Figure shows the 2-D CMG model after one hour of water injection. Figure shows the comparison of water saturation profile along the injection–production direction for the developed numerical model based on the finite element technique used in this study and CMG software based on the finite difference (FD) approximation method. It can be seen from Figure that a good match was achieved and this proves the ability of the developed numerical model to predict the water saturation profile precisely.

Figure 8

Water saturation propagation (a) after 5 min of water injection and (b) after 1 h of water injection under σmin = 5000 psi, σmax = 5500 psi, and Pr = 4500 psi.

Figure 9

Water saturation propagation after 5 min of water injection using CMG-IMEX software.

Figure 10

Comparison of the water saturation profile along the injection–production direction for the developed numerical model based on the finite element modeling (FEM) in this study and CMG-IMEX software based on the finite difference (FD).

Validation Using the Laboratory Drainage Test

In this section, the developed two-phase fluid flow model is validated against the experimental data. Glass bead experimental data are collected for a drainage test[25] (produced volume of water). The glass bead is simulated using the developed methodology in this study. Thr matrix used in the generated mesh is represented by the tetrahedral element (see Figure ), and triangle elements are used for meshing the 0° vertical fracture. Matrix permeability is 3.4 D and permeability of the vertical fracture is 104 D. The viscosity of soltrol-130 used during the experiment 2 cp and the density is 0.8 gm/cm3, while water viscosity is 1.002 cp and density is 0.998 gm/cm3. The experiment is started by injecting soltrol-130 at the top end of the glass bead at a constant rate of 3 cm3/min, and fluids that are produced due to the displacement process are collected at the bottom end of the glass bead model (see Figure ). The results of this drainage test show worthy agreement between the simulated and experimental data (see Figure ) regarding the cumulative produced volume and pore volume (PV) injected, which proves the robustness of the developed simulation model.

Figure 11

Two-phase flow profile for vertical fracture at 0.3 pore volume injected PVI.

Figure 12

Simulated and experimental volume of the produced fluid for vertical fracture at 0o.

3-D Single-Phase Fractured Reservoir Case Study

In this case study, a block with 100 m x 200 m x 100 m is created using Gmsh with 100 fractures as shown in Figure . The created fractures have a radius of 30–120 m. The reservoir and fluid properties used in this model are shown in Table collected from the Soultz geothermal reservoir. An injector is used with 6200 psi and a producer is used with a constant bottom hole flowing pressure of 2000 psi. The finite element mathematical model developed in this study is used to calculate the pressure and velocity distribution along the fractures and matrix within the fractured block. Figure shows the calculated permeability tensors in 3-D space using the periodic boundary conditions. Figure shows the pressure distribution after one month and 3 months of injection, while Figure shows the velocity distribution. It can be seen from Figure that the pressure is high around the injector and propagates in the direction of fractures and production well. This is due to the low permeability value used for the matrix and high permeability for the fractures.

Figure 13

Reservoir block with dimensions of 100 m × 200 m × 100 m.

Table 2

Model Rock and Fluid Properties

fluid properties	value
fluid compressibility	1.0E-4 psi–1
fluid viscosity	1 cp
producer Pwf	1000 psi
injector Pwf	6200 psi
initial reservoir pressure	2000 psi
formation and fractures properties
formation porosity	0.14
fractures permeability	100 *103 md
formation permeability Kx	1 md
formation permeability Ky	1 md
formation permeability Kz	1 md
horizontal stress	6000 psi
vertical stress	3000 psi

Figure 14

Calculated permeability tensor in 3-D space (a) Kxx, (b) Kyy, and (c) Kzz.

Figure 15

Pressure distribution with Pinj = 6200 psi, Pprod = 2000 psi, σH = 6000 psi, and σv = 2000 psi. (a) After 1 month of water injection and (b) after 3 months of water injection.

Figure 16

Fluid velocity distribution with Pinj = 6200 psi, Pprod = 2000 psi, σH = 6000 psi and σv = 2000 psi. (a) After 1 month of water injection and (b) after 3 months of water injection.

Figure shows the fluid velocity distribution after 1 month and 3 months. It can be seen form Figure that the fluid velocity is high inside the fractures compared to the fluid injected velocity in the adjacent matrix.

3-D Case Study of Multiphase Fluid Flow

This case is for a full field study of a typical fractured basement reservoir in Vietnam. The field dimensions are 10 km × 25 km × 300 m with one well at the reservoir center penetrating the basement. A drill stem test is conducted in this reservoir to test the productivity of the formation. The tested formation thickness was 71 m. Well testing results show that the oil production rate is 850 bbls/d with a tubing head pressure of 1344 kPa. The aim of this section is to test the production potential of the basement under depletion and waterflooding mechanisms. The field dimensions are immense to simulate. Thus, the field is divided into different sectors based on the calculated 3-D permeability tensors (see Figure ). 3-D permeability tensors are calculated for 10,000 fractures with length (l < 800 m) by dividing the whole field into a number of grid blocks with a size of 100 m × 100 m × 50 m. In addition, 4000 fractures (l > 800 m) are explicitly discretized in the reservoir domain. As a result of dividing a full field model into different sectors to make the simulation process simple, flux boundary conditions are applied on each sector established from the full field run.

Figure 17

3-D block permeability tensors for short fractures (K) for the whole field.

The porosity of the basement reservoir is calculated by relating the log interpreted porosity and permeability tensors (see Figure ). Based on the calculated 3-D porosity, the original oil in place (OOIP) of the whole field is calculated using the volumetric method.

Figure 18

3-D calculated porosity map for the whole field using the available well log data.

Results and Discussion

Fluid flow is simulated using the developed multiphase numerical model in this study to test the production potential of the studied reservoir. At first, the numerical simulation results show that the oil production rate obtained at the opening of the tested well is 960 bbls/d and this is very close to the test data (850 bbls/d). The difference in rates between the simulated data and well test data is related to ignoring the effect of wellhead pressure in the simulation model. Then, a sector model is selected to test the reservoir productivity under different drive mechanisms. This sector is located in the eastern south part of the reservoir with dimensions of 5000 m × 5000 km × 300 m (see Figure ). The OOIP for this sector is (1075 MM barrel) calculated using the initial water saturation of 0.3 and the oil formation volume factor of 1.33 bbl/STB.

In this sector, 1100 discrete fractures longer than 800 m were used in the fluid flow simulation model (see Figures and 19). The reservoir parameters used are presented in Table . The reservoir productivity is tested first under a depletion drive mechanism. The fluid flow is simulated in this sector for 14 years using 16 wells. These wells are placed based on the areas that have high values of permeability tensors. Figure shows the oil production rate under the depletion drive mechanism. Figure shows that the oil-producing rate decreases rapidly. The calculated cumulative oil production in this case is 32.7 × 106 bbls after simulation of 14 years. The simulated oil recovery is about 3%. These results show that even by using 16 wells in this sector, the productivity of the reservoir is still very low. Thus, to increase the recovery factor, the waterflooding technique is used with injectors placed in the reservoir center to increase the reservoir sweep efficiency. Four injectors are selected in this sector around the producing wells and water is injected under high constant pressure (30 MPa). The water is injected into the bottom section of the reservoir at depth (1600–1700 m) while the oil producing interval depth is 1400–1500 m. There is 100 m interval separating the injection and production zones. Figure shows the produced water cut at the production wells. As can be seen from Figure , some of these wells have high water cut (#W2, #W12, and #W15) and the others have low produced water cut (#W4, #W10, and #W14). This is related to the high or low permeability regions intersected with these wells. The water production increasing rate is very high and is like a vertical trend in the first 3 years. Contrarily, the oil decreasing rate is high within the first 3 years of production and the decreasing rate becomes stable later (see Figure ). This behavior of water and oil production rates is related to the sudden pressure drop around areas of producing wells. This high drop in pressure resulted from the highly connected fractures surrounding the producing wells. However, later, the oil production rate becomes stable and this is due to the formation of large water buffers around the injected wells at later stages. The calculated cumulative oil production under the waterflooding mechanism is 5.7 × 10+7 bbls. The estimated recovery factor is 5.3%. The results show that under the waterflooding mechanism, the oil recovery is increased using 14 producing wells. Though the increasing rate of oil recovery is still low, it is recommended to increase the number of injection wells to increase the sweep efficiency and to compensate for the drop in pressure at the beginning of the production process.

Figure 19

3-D tetrahedral mesh generated for the selected sector. The well locations shown in this mesh include 4 injectors and 12 producers. The generated mesh includes 1100 discrete fractures.

Table 3

Reservoir Parameters for the Selected Sector

parameter	value
sector 3-D	5000 m × 5000 m × 300 m
aperture of fracture (m)	1.5 × 10–4
permeability of matrix (mD)	10–5
initial pressure (MPa)	15.16
bubble point pressure (MPa)	6.2
Pwf (MPa) (production)	6.21
Pinj (MPa)	30
Swi	0.37
μo (cp)	0.82
ρo (kg/m3)	818
βo	1.4
Cf (MPa–1)	20.3 × 10–5
horizontal stresses (MPa)	33.1
vertical stress (MPa)	41.3
Young’s modulus (GPa)	42
Poisson ratio	0.28
ρs (kg/m3)	2800
wellbore radius (m)	0.1

Figure 20

Rate of produced oil for the depletion drive mechanism scenario with Pinit = 15.6 MPa and Pprod = 6.21 MPa.

Figure 21

Wells’ water cut with Pinj = 30 MPa, Pinit = 15.6 MPa, and Pprod = 6.21 MPa.

Figure 22

Rate of produced oil for the waterflooding mechanism scenario with Pinj = 30 MPa, Pinit = 15.6 MPa, and Pprod = 6.21 MPa.

Conclusions

This study presented the comprehensive derivation of two-phase flow equations through naturally fractured reservoirs in a poroelastic environment. The developed poroelastic numerical model is validated in a decoupled way. First, the numerical model is validated using only elastic problems. Then, the model is validated against commercial black oil reservoir simulators (CMG). A case study of the fractured block is used to validate the robustness of the developed finite element model in predicting the pressure and velocity distribution along the intersected fractures. The results of the validation process indicate that the developed numerical poroelastic model is able to predict the behavior of multiphase fluid flow through the porous fractured reservoirs.

The results of the field case study and experimental drainage test show that the model can be used to assess the production potential of the naturally fractured reservoirs in terms of the in place e and wells production. Additionally, a new methodology will be proposed to upscale two-phase relative permeability from the laboratory scale to field conditions under a poroelastic environment in part II.

The different steps of the presented innovative approach in this study are used together to evaluate the production potential of these reservoirs under waterflooding mechanisms. The results of flow simulation in this case show that the injection of water in the bottom section of the basement assisted in generating a large water buffer and consequently forming oil/water contact. The formation of oil/water contact supports the reservoir pressure, and the oil recovery is increased. In addition, using five spot injection patterns assisted in increasing the reservoir sweep efficiency over peripheral water injection.