This paper describes the formulation and application of a fully implicit, three-dimensional, two-phase (water and oil) mathematical model for the simulation of naturally fractured reservoirs. The model equations are discretized using the control-volume finite element (CVFE) scheme where pressure terms are treated in a finite element manner while mobility terms are upstream-weighted as in the usual finite difference (FD) method. The form of the discretized equations is similar to the usual FD method, thus standard techniques are used to solve the Jacobian. The equations are solved fully implicit for pressure, saturation, and wellbore variables using the iterative Newton-Raphson method of solution.

Two field scale examples, that were first presented by Kazemi et al., were used to test the model. These include a quadrant of a five-spot pattern, and a five-well fractured reservoir having dip and natural water influx.


Performance of naturally fractured reservoirs is complicated and more difficult to predict than conventional reservoirs due to the different mechanisms involved in the production of oil and gas. The complexity arises from the existence of two flow paths (matrix and fracture which have totally different properties) that communicate with each other.

Because of the irregularity of the fracture morphology, distribution, orientation, and extension, some sort of idealization of these reservoirs is necessary for the simulation purposes.

Kazemi et al. (1976) set the foundation for field-scale, dual-porosity, multi-phase simulators. In their idealization, an elemental reservoir volume in a naturally fractured reservoir (Fig. 1) is represented by a set of grid cells. Each grid cell may contain one or several matrix blocks. Accordingly, within any grid cell, all matrix blocks will have the same pressures and saturations. Similarly, the fractures within a grid cell will have identical pressures and saturations which differ from those of matrix blocks. A typical grid cell with its pressure distribution is shown in Fig. 2. This idealization was based on the following assumptions:

  1. the matrix blocks are homogeneous and isotropic and form a uniform assemblage of identical, rectangular, parallelopipeds,

  2. the fractures are contained within an orthogonal system of continuous, uniform, and constant-width paths where the fractures are parallel to the principal axes of permeability.,

  3. the flow takes place from fracture to fracture, from fracture to matrix, and from matrix to fracture,

  4. the matrix blocks have high storage volume and low flow capacity whereas fractures have low storage volume and high flow capacity, and

  5. the wellbore intersects the fracture network and hence the fractures provide the main path for fluid flow into the wellbore.

In many practical problems, this idealization is sufficient. In cases where a matrix block contains several grid nodes, the gravity segregation within the matrix block is calculated. This description points to the critical engineering decision in the selection of the number of grid nodes with respect to the number of matrix blocks.

This content is only available via PDF.
You can access this article if you purchase or spend a download.