NUMERICAL SIMULATION OF A NATURALLY FRACTURED RESERVOIR .

This research is concerned with the development of a numerical model for stratified normally fractured reservoirs. Three dimensional three phase flow black oil simulation model is adopted. The dual porosity-dual permeability model is used. The IMPES (Implicit Pressure Explicit Saturation) method is used to solve the difference equations. The Tertiary trap in K oil field (an Iraqi oil field) was simulated by the numerical model. The trap consists of six layers having different properties. Equally spaced Cartesian grids were used to divide each layer into 1600 cells in the x-y plane with the thickness as the dimension of each grid block in z direction. Applying the two IMPES pressure equations to each grid of the simulated domain resulting in a block seven diagonal coefficient matrix. Gauss-Seidel iterative method was used to solve the system of equations. The time steps are controlled through a maximum saturation difference and a material balance error limits. The actual production histories of the 15 wells in K oil field are used to get the past performance of the field for the production period. The calculated and measured average reservoir pressures, produced gas/oil ratio through the production periods had acceptable match.


Introduction
The development of a simulator for naturally fractured reservoirs (NFR) is a real challenge from both the reservoir description and numerical solution point of view.Fluid flow behavior in fractured reservoir through highpermeability low effective porosity fracture surrounding low-permeability high porosity matrix block has been described extensively in the oil literature during the last thirty-five years.The fluid exchange between the fracture network and the individual matrix blocks is the most physical aspect of the fluid flow problem to characterize.
4 developed 3D, multiple well numerical simulator for simulating single or two-phase flow of water and oil in a fractured reservoir.The simulator equations are a two-phase flow extension of single -phase flow equation derived by 6 The simulator account for relative fluid mobility, gravity forces, imbibition, and variation in reservoir properties.It handles uniformly and non-uniformly distributed fractures and for no fractures at all.A semi-implicit finite difference expansion had been used to solve the original dual porosity equations.Different methods of solving the system of equations were proposed depending on the number of nodes.The results showed the significance of imbibtion on recovery of oil from the rock in reservoirs with interconnected fracture network.
3 modeled the flow in the fracture system by representing fluid transfer from the matrix into the fracture by a "source" term and fluid transfer from the fracture to the matrix by a "sink"(or negative source) term.
8 presented a stable, flexible, fully implicit, finite difference simulator in heterogeneous, dual-porosity reservoir.They used the flow equation proposed by 7.
Cartesian and radial coordinates are included in the model.The conventional five point finite difference in the x-y plane was extended to a special nine-point formulation to account for the directional flow other than the x-y directions.Each node in the model has two properties one for the matrix and the other for the fracture.Good agreement was noticed between the analytical and numerical solution for pressure build up prediction.
Thomas, 5 used the dual porosity model to develop 3D, three phase simulator for NFRs.The same flow equations of 7 were used.The formulation was implicit in pressure, water saturation and gas saturation for both matrixfracture flow and fracture flow.The gravity and capillary effect were incorporated.After expanding the matrix-fracture flow equation in totally implicit form the matrix unknown were eliminated in terms of fracture unknowns to reduce the total number of unknowns.The time steps were controlled automatically using a maximum saturation change of 0.1.
3 described a 3D three phase compositional simulator.A dual permeability and/or a dual porosity system were used to describe complex porous media including highly fractured, micro fractured or non-fractured regions.In addition to the viscous and capillary pressure forces, the matrix-fracture exchange term can handle gravity effects.The conservation equations were expressed in compositional form and equilibrium K-values were used.The fully implicit equations are linearized by Newton-Raphson iteration scheme.Because of the multi-purpose nature of the model, several different choices of discretized time-solution techniques are available.9 presented an empirical formulation for the transfer function representing the matrix-fracture interaction in the dual-porosity model.Depending on the assumption that when water imbibtion is the dominant force for displacing oil from the matrix.The aim of the study is to develop a numerical simulation model for a naturally fractured stratified Iraqi reservoir (the ministry of oil does not give permission to state the name or publish the map of the reservoir) and check the match between the actual reservoir history with that predicted using the simulator.

Field Description
K oil field is an Iraqi oil field located in the north of Iraq.The field is a simple asymmetrical doubly plunging anticline.Its main axes strike NW-SE.The slope of the NE flank is between 9o to 13o while the slope of SW flank reaches 20o in some locations (geological study 1992).The structure is about 17 Km long and 6 Km wide.This work is concentrated on the Tertiary trap which consists of six carbonate units having different thickness and different petrophysical properties.The structural map of the trap is shown in Fig. 1.Core sample studies showed that the fractures are homogeneously distributed in the field and can be divided into open, closed, completely filled and partially filled fractures.The degree of fracturing in the units ranges from rare in the top and bottom units to open fractures and vugs in the middle units.
The reservoir has a large gas dome, medium oil column and water at the flanks.Fifteen producing wells were drilled in the trap during the 1980ths having different production history (field measurements 2007).
Grid Construction Each layer in the reservoir was subdivided using a grid system having equal spacing in the x and y directions (200 meter) spacing.Layer thickness was considered as the spacing in the z direction.
Block centered grid and row ordering methods were used.The grid network with the location of the producing wells is shown in Fig. 2.
Transmissibility Evaluation Single point upstream weighing is used to evaluate transmissibility at the block boundaries.The fluid potential is used to recognize the upstream cell from which fluids are flow to the adjacent one.Each grid block, not in the boundary of the simulated area, is communicated with six blocks.So for each flowing phase the transmissibility at the six blocks boundaries are calculated after assigning the upstream one in each direction.

Flow Equations
The dual porosity dual permeability model is used in which the flow in the reservoir occurs in both fracture and matrix system in addition to fluid exchange between those systems.The equations describing three phase three dimensions flow in the fracture system written in finite difference form are:- The matrix equations are same as eqs.(1, 2 and 3) but using the parameters of the matrix system.The transfer function, τ, governing fluids transfer between the fracture and matrix system is defined as: Where σ is the shape factor, λ is the mobility of phase α (water, oil or gas).The shape factor σ in eq. ( 4) is calculated from the relation (Kazemi 1976) but a factor of 2 is used instead of 4 in the original equation.
The difference equations 1, 2 and 3 are solved using the IMPES method resulting in two main pressure equations in terms of fracture and matrix pressure (Ahmed 2007).
(7) Equations 6 and 7 are applied to the grids in the simulated domain, therefore 2N equations resulted from N grid blocks.The 2N equations are overlapped together to get a final equation having a matrix form of: where A is the coefficient matrix, P is the pressure vector and D is the right hand side vector.
Matrix A is a seven diagonal matrix consisting of six individual matrices (blocks), each individual matrix represents one of the six layers forming the field.Each non-zero entry in matrix A is a 2x2 matrix.The components of matrix A and D are initially calculated using the pressure and saturation values for water, oil and gas resulted from the initialization process.
To solve this equation Gauss-Seidel iteration method is used in the following manner: 1-Take a certain time step and use the pressure distribution values results from initialization process as initial pressure values.
2-Calculate the new pressure values in the 1st layer for the fracture and the matrix using the Gauss-Seidel iterative method.
3-Check the maximum difference between the initial and the calculated values, if the difference <= .01psi then proceed to step 4 other wise use the newly calculated pressure values in step 2 as an initial values and repeat calculations starting from step 2. 4-Repeat steps 2 and 3 for the layers 2 to 6. 5-Use the final pressure distribution resulted from the above calculations as initial pressure values and repeat calculations from step 2 to step 5 checking for maximum pressure difference of .01psi between all initial and calculated values to stop the calculations.
After calculating the pressure distribution in the field at the new time level, calculation of oil and water saturations at the new time level in the fracture are calculated from eq. ( 1) and (2) also these saturations in the matrix are calculated using the same equations but using the matrix properties.
Check is made for maximum saturation difference in all the cells, if the difference is

Fig. 3 .
Fig.3.Comparison between the measured and calculated pressure.