# A Rapid and Efficient Methodology to Convert Fractured Reservoir Images Into a Dual-Porosity Model

B. Bourbiaux, M. C. Cacas, S. Sarda, J. C. Sabathier

► To cite this version:

B. Bourbiaux, M. C. Cacas, S. Sarda, J. C. Sabathier. A Rapid and Efficient Methodology to Convert Fractured Reservoir Images Into a Dual-Porosity Model. Revue de l'Institut Français du Pétrole, 1998, 53 (6), pp.785-799. 10.2516/ogst:1998069. hal-02079046

HAL Id: hal-02079046

https://ifp.hal.science/hal-02079046

Submitted on 25 Mar 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A RAPID AND EFFICIENT METHODOLOGY TO CONVERT FRACTURED RESERVOIR IMAGES INTO A DUAL-POROSITY MODEL*

B. BOURBIAUX, M.-C. CACAS, S. SARDA
and J.-C. SABATHIER

Institut français du pétrole¹

MÉTHODOLOGIE RAPIDE ET EFFICACE POUR CONVERTIR LES IMAGES DE RÉSERVOIR FRACTURÉ EN MODÈLE À DOUBLE POROSITÉ

La caractérisation et la simulation dynamique des réservoirs naturellement fracturés ont bénéficié d'avancées importantes ces dernières années. Toutefois, l'ingénieur réservoir reste confronté à la difficulté de paramétrer le modèle équivalent à double porosité utilisé pour représenter de tels réservoirs. En particulier, les perméabilités de fracture équivalentes et les dimensions du bloc matriciel équivalent ne peuvent pas être facilement déduites de l'observation des images complexes de réseaux naturels de fractures. Cet article décrit une technique nouvelle et systématique pour calculer ces paramètres équivalents. Les résultats de sa mise en œuvre au moyen d'un logiciel spécifique démontrent sa validité et son efficacité pour l'étude de cas de champ. Un tenseur des perméabilités de fracture équivalentes est déduit de calculs d'écoulements stationnaires et monophasiques dans le réseau réel de fractures, ce dernier étant assimilé à un réseau 3D de résistances et soumis à des conditions aux limites spécifiques. Les dimensions du bloc équivalent dans chaque couche sont obtenues rapidement par identification d'une fonction géométrique représentative d'un déplacement capillaire. La méthodologie a été validée par comparaison avec des simulations à maillage fin réalisées au moyen d'un simulateur de réservoir conventionnel. Puis, une image complexe d'affleurements d'une formation de grès a été traitée à titre de démonstration. L'outil innovant associé à cette méthodologie donne la possibilité à l'ingénieur réservoir de construire un modèle à double porosité qui reproduit de façon plus fidèle le comportement hydraulique du milieu réel fracturé.

## A RAPID AND EFFICIENT METHODOLOGY TO CONVERT FRACTURED RESERVOIR IMAGES INTO A DUAL-POROSITY MODEL

Both characterization and dynamic simulation of naturally-fractured reservoirs have benefited from major advances in recent years. However, the reservoir engineer is still faced with the difficulty of parameterizing the dual-porosity model used to represent such reservoirs. In particular, the equivalent fracture permeabilities and the equivalent matrix block dimensions of such a model cannot be easily derived from observation of the complex images of natural fracture networks. This paper describes a novel and systematic methodology to compute these equivalent parameters. The results of its implementation with specially-designed software demonstrate

<footnote>(1) 1 et 4, Avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex - France</footnote>

<footnote>* This paper was prepared for presentation at the 1997 SPE Annual Technical Conference and Exhibition held in San Antonio, Texas, 5-8 October, 1997.</footnote>
its validity and efficiency in dealing with field situations. A tensor of equivalent fracture permeability is derived from single-phase steady-state flow computations on the actual fracture network using a 3D resistor network method and specific boundary conditions. The equivalent block dimensions in each layer are derived from the rapid identification of a geometrical function based on capillary imbibition. The methodology was validated against fine-grid reference simulations with a conventional reservoir simulator. Then, a complex outcrop image of a sandstone formation was processed for demonstration purposes. This innovative tool enables the reservoir engineer to build a dual-porosity model which best fits the hydraulic behavior of the actual fractured medium.

## METODOLOGÍA RÁPIDA Y EFICAZ PARA CONVERTIR LAS IMÁGENES DE YACIMIENTOS FRACTURADOS EN MODELO DE DOBLE POROSIDAD

La caracterización y la simulación dinámica de los yacimientos naturalmente fracturados han aprovechado las ventajas derivadas de los importantes adelantos que han surgido durante estos últimos años. No obstante, el ingeniero de yacimientos sigue estando acuciado por la dificultad consistente en parametrar el modelo equivalente de doble porosidad utilizado para representar este género de yacimientos. Básicamente, las permeabilidades de fractura equivalentes y las dimensiones del bloque matricial no se pueden deducir fácilmente de la observación de las imágenes complejas de redes naturales de fracturas. Se describe en el presente artículo una técnica, nueva y sistemática, para calcular estos parámetros equivalentes. Los resultados de su implementación por medio de un software específico demuestran su validez y su eficacia para el estudio de casos de campo. Un tensor de las permeabilidades de fractura se deduce de cálculos de flujos estacionarios y monofásicos en la red real de fracturas asimiladas a una red 3D de resistencia y sometido a condiciones límites específicas. Las dimensiones del bloque equivalente en cada estrato se obtienen rápidamente por identificación de una función geométrica representativa de un desplazamiento capilar. La metodología se ha validado por comparación con las simulaciones de mallado fino, realizadas por medio de un simulador de yacimiento convencional. A continuación, se ha procesado a título de demostración una imagen compleja de afloramientos de una formación de arenisca. La herramienta innovadora asociada a esta metodología proporciona la posibilidad al ingeniero de yacimiento de construir un modelo de doble porosidad que reproduce de forma más exacta el comportamiento hidráulico del medio real fracturado.

## INTRODUCTION

Modeling naturally-fractured reservoirs remains a difficult reservoir engineering task, although a great deal of effort has recently been made to better characterize the geometry and distribution of fractures within the reservoir and to increase the capabilities of dual-porosity simulators. Actually, the 3D images of fracture networks derived from field fracturing data integration and stochastic modeling are most often too complex to be used as direct input for a reservoir simulator. Typical horizontal and vertical dimensions of a reservoir simulator cell are respectively of about several hundreds and tens of meters due to limitations in the size of field-scale numerical models. Within such cells, the centi-metric to decametric fractures behave like an equivalent medium, that is, the flow properties of the small-scale fracture network are representative at cell scale and can be described by a permeability tensor [1]. However, faults crossing the reservoir cannot be incorporated into the equivalent fracture medium defined at cell scale but have to be discretized individually.

In the early sixties, Warren and Root [2] introduced an approach to simulate flows within fractured reservoirs modeled as the superposition of two equivalent fracture and matrix media. This equivalent model consists of an array of parallelepiped matrix blocks limited by a set of uniform orthogonal fractures (Fig. 1). Fracture flows are computed within the fracture grid, matrix-fracture transfers are computed at each grid block position, and block-to-block exchanges are also accounted for in the dual-permeability option of

<figure><img src="image_2.png" /><figcaption>Figure 1
Conventional Warren and Root representation of a fractured reservoir.</figcaption></figure>
such simulators. The Warren and Root representation of fractured reservoirs is used in most industrial simulators. It minimizes computation requirements for field-scale simulation, however, it does not account for flow distribution at scales below that of homogenization.

The first difficulty encountered in the use of the Warren and Root equivalent model is to find equivalent fracture permeabilities and equivalent parallelepiped blocks with which to simulate the reservoir flow dynamics.

A methodology is presented in this paper to compute equivalent permeabilities and equivalent blocks from 3D geological images of the fracture network. First, a general format has been defined for processing the 3D images of any fracture characterization model. The methodology itself includes the determination of an equivalent fracture permeability tensor and of an equivalent block section for any layer or stack of layers. It has been validated against reference solutions provided by finely-discretized models. And finally, its implementation with specially-developed software is demonstrated for a naturally-fractured outcrop rock volume.

Before detailing our methodology, it is worthwhile to carry out a quick literature review of existing approaches. All published works concern the determination of equivalent permeabilities. No systematic method is reported for determining equivalent block dimensions.

The approaches detailed [3-7] differ in the flow geometry, the role of the matrix medium, the technique of discretization of the fracture network and the method used for computing flows in fractures. Long *et al.* [3], Cacas *et al.* [4] and Massonnat *et al.* [5] developed 3D flow models. Lough *et al.* [6] developed a 2D fracture flow model including the contribution of 3D matrix flows and based on the boundary element method.

The matrix permeability is taken into account in Massonnat's, Odling's [7] and Lough's models.

The technique used for representing and discretizing the fracture network is probably the main feature that distinguishes each of those models. Discrete fracture models are found in [3, 4, 6]. Long built a 3D model with discs representing fractures and nodes at disc intersections. Cacas developed a similar model but placed nodes at the centers of discs to minimize the number of nodes. Lough developed a 2D fracture flow model with a numerical technique that takes into

account 3D matrix flows without having to discretize matrix blocks. The 2D model developed by Odling is close to the discrete models, however fractures are discretized as the linear elements of a fine square grid superposed on the actual 2D network. Massonnat used a fine regular 3D grid to discretize both fracture and matrix media, with an average permeability of fracture and matrix media assigned to the cells crossed by one fracture. Such a method involves a very high number of cells and cannot be applied to large highly-fractured media.

All discrete models are used to compute steady-state flows on the basis of the resistor network approach, i.e. mass-balance equations written at each node, and of pressure or flow conditions imposed on the boundaries. Fracture flows are formulated using Poiseuille's law in all models except for Long's model involving analytical fracture flow equations.

# 1 THE METHODOLOGY

The starting point for our study is a 3D fractured reservoir image representing a naturally-fractured rock volume that is provided by any fracture characterization software. The equivalent model for this image is determined in the framework of the dual-porosity concept. Therefore, the fracture network alone is taken into account to compute equivalent fracture permeabilities, since fracture flows and matrix flows are computed within two distinct grids in a dual-porosity simulator. Fractures are assumed to be vertical or, more precisely, orthogonal to bed limits, which is acceptable for many reservoirs. Hence, the matrix medium is continuous along the vertical direction and only the dimensions of the horizontal section of an equivalent block need to be determined.

## 1.1 Image Format

The parallelepiped 3D image is discretized vertically as a list of layers. Each layer contains a certain number of fracture plane elements denoted by rectangles (Fig. 2). Each rectangle is defined by a fracture set number, an origin, azimuth and dip, dimensions, neighboring rectangles, and by a fracture aperture or a conductivity. The conductivity is defined as the product of the intrinsic fracture permeability and of the aperture *e* of the fracture represented by parallel
<figure><img src="image_3.png" /><figcaption>Figure 2
Format used to represent fractures.</figcaption></figure>

walls. Poiseuille's law may be used to express the intrinsic fracture permeability and conductivity respectively as $e^2/12$ and $e^3/12$.

This general format makes it possible to process 3D images of any characterization model with our methodology and software.

## 1.2 Equivalent Permeabilities

They are derived from incompressible steady-state flows computed on the actual 3D fracture network. For each direction x, y, z, the flow rate resulting from a given pressure drop $\Delta p$ applied between the two opposite faces of the parallelepiped network is computed with fixed lateral boundary conditions.

Either a zero-flow-rate condition or a linearly-varying pressure condition is applied on the lateral limits, as shown in Figure 3. The flow rate solution is unique and leads, through the use of Darcy's law, to a direct determination of the equivalent permeability in the direction considered. The steps and algorithms involved are described hereafter.

First, the network is discretized as nodes placed at rectangle intersections. A mass-balance equation is written for each node i connected to n nodes j:

$$ \sum_{j=1}^{n} q_{ij} = 0 \qquad (1) $$

Using Darcy's law, each flow rate is expressed as a function of node pressures $p_i$ and $p_j$ and transmissivity $T_{ij}$ between the two nodes i and j:

$$ q_{ij} = T_{ij} (p_i - p_j) \qquad (2) $$

The transmissivity between two nodes is computed as a function of fluid viscosity and of the dimensions and conductivity of the fracture element linking them.

The first equation can be rewritten as:

$$ \left( \sum_{j=1}^{n} T_{ij} \right) p_i - \sum_{j=1}^{n} T_{ij} p_j = 0 \qquad (3) $$

<figure><img src="image_4.png" /><figcaption>Figure 3
The two kinds of boundary conditions for computing equivalent fracture permeabilities.</figcaption></figure>

<figure><img src="image_5.png" /></figure>
If node *i* is also connected with a limit at a fixed known pressure $p_B$, the previous equation becomes:

$$ \left[ \left( \sum_{j=1}^{n} T_{ij} \right) + T_{iB} \right] p_i - \sum_{j=1}^{n} T_{ij} p_j = T_{iB} p_B \quad (4) $$

For the total number of nodes *N* discretizing the fracture network, this leads to a system of *N* equations with *N* pressures as unknowns. The resulting system of equations leads to a symmetric positive definite coefficient matrix which is solved with a conjugate gradient method. The flow rate is then derived from Equation (2) for each couple of nodes as well as the total flow rate through the fracture network in the direction of flow considered. A Darcy fluid velocity *v* is determined on the scale of the entire fractured medium as the ratio between the total flow rate and the cross-section area of the fractured medium. The equivalent permeability is finally calculated by applying Darcy's law between the upstream and downstream faces of the 3D volume. If lateral faces are impermeable, a single equivalent permeability can be determined for each flow direction, such as, for instance, in the *x* direction:

$$ k_{xx} = k_x = \frac{\mu v_x L_x}{\Delta p} \quad (5) $$

If linearly-varying pressures are imposed on lateral faces, then flow rates can be computed in transverse directions. One can then derive the extra-diagonal terms of an equivalent permeability tensor, such as, for instance:

$$ k_{xy} = \frac{\mu v_{xy} L_y}{\Delta p} \quad (6) $$

Considering that the direction orthogonal to bedding planes is generally the direction of least permeability, mixed conditions may be imposed to determine the following partial permeability tensor:

$$ \begin{pmatrix} k_{xx} & k_{yx} & 0 \\ k_{xy} & k_{yy} & 0 \\ 0 & 0 & k_{zz} \end{pmatrix} \quad (7) $$

These mixed conditions consist of imposing a zero flow rate on all lateral limits for vertical flows, and a zero flow rate on horizontal limits and a pressure varying linearly on vertical lateral limits for horizontal flows. The partial permeability tensor derived is sufficient to determine the main directions of flow in the bedding planes with the associated equivalent permeabilities.

## 1.3 Equivalent Block

Following the assumption of vertical fractures, the problem addressed is that of determining the square or rectangular section of an equivalent matrix block for each layer or group of layers having similar fracturing features. To this end, a simple and rapid procedure based on the water-oil capillary imbibition mechanism and denoted as the geometrical method has been developed. The matrix medium is assumed to have constant petrophysical properties throughout the bedding plane. Capillary imbibition of a matrix block involves countercurrent flows of water and oil [8] with a gradual advance of a smoothly-tilted front from the boundaries to the center of the block. The geometrical method is a simplified representation of this flow mechanism, based on the two following assumptions:

- piston-type invasion of the matrix medium by displacing fluid;

- a distance *X* of the front from block boundaries changing with time *t* according to a function which is independent of block shape for a given rock-fluid system.

Because of the piston-type flow assumption, the matrix oil recovery *R* can be expressed as the normalized area *A*, i.e. as the ratio between invaded and total cross-section areas of the block. In addition, with the bi-univocal relationship between *t* and *X*, one replaces the recovery as a function of time, *R(t)*, by the normalized area as a function of distance, *A(X)*. Finally, the geometrical method, schematized in Figure 4, determines the dimensions *a* and *b* of an equivalent block section through a procedure minimizing the difference between the two normalized area versus distance functions referring respectively to the layer, *A(X)*, and to the equivalent square or rectangle, *A<sub>eq</sub>(X)*. *A(X)* is computed using a rapid, patented procedure for image processing. *A<sub>eq</sub>(X)* is a second-order polynomial:

$$ A_{eq}(X) = \frac{2X}{a} + \frac{2X}{b} - \frac{(2X)^2}{ab} \quad (8) $$

Normalized area versus distance curves are discretized with a constant normalized area step to give the same weight to any volume of oil recovered independently of the time at which it is produced. The geometrical method can be applied to a stack of layers with similar petrophysical properties of the matrix medium. For such a case, the normalized area versus distance function is a thickness-weighted sum of the normalized area versus distance functions that refer to the constitutive layers.
<figure><img src="image_6.png" /><figcaption>Figure 4
The geometrical method.</figcaption></figure>

## 1.4 Implementation

Software has been developed to apply the methodology described above to geological images of naturally-fractured media. This software can process the images given by any characterization model as long as the format described above is respected. Among the capabilities of the software, we must mention the possibility to delimit in area and vertically sub-volumes within the 3D fractured medium, making it possible to determine whether the equivalent parameters are representative. The image definition can also be modified and optimized in order to be able to process large images with a sufficient definition of blocks. Equivalent fracture permeabilities are computed with either zero-flow-rate lateral conditions or mixed conditions as described above. Equivalent blocks are determined for each layer and for any pre-defined stack of layers.

# 2 VALIDATION

The previous methods have been validated for fracture networks of limited complexity and for blocks of various shapes.

## 2.1 Equivalent Permeabilities

Reference solutions are computed on fine regular grids discretizing 3D fracture networks containing a limited number of orthogonal fractures oriented along the *x*, *y* and *z* axis. Grid refinement has to be sufficient to avoid the presence of several fractures in a given cell or of disconnected fractures in two contiguous cells. Thus, the permeability input value for each cell of the fine-grid model is the permeability equivalent to matrix and fracture media if one fracture crosses the cell in the direction considered, and zero in the other situations. Validation was performed for different networks with increasing numbers of fractures and connectivity. Three networks were studied, with respectively 3, 10 and 20 fractures. The conductivity was arbitrarily set to 1 mD·m for all fractures. Equivalent permeabilities were computed with the zero-flow-rate lateral boundary condition because this is the only kind of condition available in the conventional flow simulator used to compute fine-grid reference solutions.

The first medium was as simple as possible to be able to find analytical flow solutions. It consists of a cube with a lateral size of 5 m, and divided into five layers. As shown in Figure 5, the three fractures have variable horizontal lengths and are present either in a few layers, or in all of them, making it possible to obtain a 3D flow geometry in the vertical direction. Only nine nodes are required to represent this network with our methodology. To obtain reference values of equivalent permeability, this simple network was discretized with a regular grid of 125 cubic cells, with respectively 25, 15 and 9 cells to represent fractures F1, F2 and F3.

The solutions obtained with the fine-grid model and with the method developed are compared in Table 1. Analytical solutions of the equivalent permeabilities in the *x* and *y* directions could also be calculated since flows were 2D, consisting of either parallel or in-series flow paths, in these directions. The *x* and *y* equivalent permeabilities obtained with the published method are in excellent agreement both with the analytical solutions and with the fine-grid results. This was expected as, in 2D conditions, the
flow problem can be solved exactly with a resistor network approach. A slight difference, of 5%, is obtained in the z direction between the fine grid solution and the method developed. This difference results from the 3D flow geometry and the minimum discretization of fracture planes used in our methodology. Actually, in 3D flow conditions, such a discretization does not account for the detailed distribution of fluids in the vicinity of fracture intersections.

<figure><img src="image_7.png" /><figcaption>Figure 5

Simple fracture network used to validate equivalent fracture permeabilities.</figcaption></figure>

The two other fracture sets were more complex and differed in the number and connectivity of fractures. The same validation approach was followed and led to the same conclusions as before. In addition, the difference in the vertical equivalent permeability values mentioned before was found to decrease when the density and connectivity of fractures increase. The situations where the vertical permeability inaccuracy becomes significant concern media with a small number of large and poorly-connected fractures, that is, situations for which the dual-porosity concept underlying our methodology is not applicable.

<table><caption>TABLE 1
Validation of equivalent permeabilities</caption><thead><tr><th>Solution</th><th>k<sub>x</sub><br>(mD)</th><th>k<sub>y</sub><br>(mD)</th><th>k<sub>z</sub><br>(mD)</th><th>k<sub>y</sub>/k<sub>z</sub></th><th>k<sub>z</sub>/√<span style="text-decoration:overline;">k<sub>x</sub> k<sub>y</sub></span></th></tr></thead><tbody><tr><td>This work (*)</td><td>0.120</td><td>0.227</td><td>0.267</td><td>1.89</td><td>1.62</td></tr><tr><td>Fine grid (**)</td><td>0.119</td><td>0.224</td><td>0.255</td><td>1.88</td><td>1.56</td></tr><tr><td>Analytical</td><td>0.120</td><td>0.226</td><td>-</td><td>1.88</td><td>-</td></tr><tr><td>(*) - (**)<br>difference (%)</td><td>&lt; 1</td><td>+ 1</td><td>+ 5</td><td>&lt; 1</td><td>+ 4</td></tr></tbody></table>

## 2.2 Main Flow Directions

As extra-diagonal terms of equivalent permeability cannot be determined with conventional flow simulators, the procedure for determining the principal flow directions in the layers of a fractured medium was validated indirectly using two networks having exactly the same fracturing features but with different fracture orientations. The two main fracture sets were oriented along the x and y axes in the first network and -30° and 60° from the x axis in the second one. Complex fracture networks could be used since no reference fine-grid solution was computed for this validation step. A nearly diagonal equivalent permeability tensor was obtained for the first network with principal flow directions oriented at 89° and 2° from the x axis. The tensor of the second network was no more diagonal and its principal flow directions were at 57° and -31° from the x axis. These two sets of results are consistent with each other and with the orientation of the two fracture sets in both images.

## 2.3 Equivalent Blocks

Reference solutions were validated for single blocks of various shapes, U-shaped, triangular and comb-shaped blocks. These blocks are shown in Figure 6 with the square or rectangular equivalent block section in gray determined with the geometrical method. The validation of equivalent block dimensions was based on a comparison of the simulated water-oil imbibition recovery curves referring to the 2D finely-discretized real block model and to the equivalent rectangular or square block determined by the geometrical method. Imbibition was simulated with matrix porosity and permeability equal to 0.30 and 1 mD, and with water and oil viscosities respectively equal to 0.35 and 0.43 cp.
<figure><img src="image_8.png" /><figcaption>Figure 6

Single blocks used to validate the geometrical method.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7

Water-oil imbibition in the comb-shaped block and in the equivalent rectangle.</figcaption></figure>

Relative permeability and capillary pressure data are given in Table 2.

<table><caption>TABLE 2

Water-oil relative permeabilities and capillary pressures</caption><thead><tr><td>Water saturation</td><td>Relative permeability to water</td><td>Relative permeability to oil</td><td>Capillary pressure (bar)</td></tr></thead><tbody><tr><td>0.20</td><td>0</td><td>0.80</td><td>2.30</td></tr><tr><td>0.25</td><td>0.0003</td><td>0.68</td><td>1.10</td></tr><tr><td>0.30</td><td>0.002</td><td>0.57</td><td>0.80</td></tr><tr><td>0.35</td><td>0.005</td><td>0.47</td><td>0.57</td></tr><tr><td>0.40</td><td>0.010</td><td>0.37</td><td>0.37</td></tr><tr><td>0.45</td><td>0.018</td><td>0.28</td><td>0.25</td></tr><tr><td>0.50</td><td>0.028</td><td>0.20</td><td>0.17</td></tr><tr><td>0.55</td><td>0.041</td><td>0.13</td><td>0.11</td></tr><tr><td>0.60</td><td>0.057</td><td>0.072</td><td>0.06</td></tr><tr><td>0.65</td><td>0.077</td><td>0.025</td><td>0.03</td></tr><tr><td>0.70</td><td>0.100</td><td>0</td><td>0</td></tr></tbody></table>

The comb-shaped block is a rectangular block of 7.6 m (25 feet) per 4.6 m (15 feet) surrounded by fractures and including several fractures initiated at block outer boundary and ending in the heart of the block. The geometrical method determined an equivalent rectangular block of 6.5 by 2.1 m. Figure 7 shows the very good consistency of the two recovery curves computed for the real block and for the

equivalent rectangle. The reference matrix oil recovery by water imbibition in the comb-shaped block was computed with a model of 1500 cells with very thin fracture cells and a nine-point scheme. Looking in further detail, it may be observed that the two curves are exactly superimposed during the first two weeks of block imbibition and slightly diverge later on. This variation results from the assumptions underlying the geometrical method. During the first 10-20 days, the flow behavior is infinite-acting, therefore the imbibition curve referring to the equivalent block solution determined under the piston-type assumption does not differ from the exact solution. Afterwards, the two recovery curves differ slightly because the effect of limited block size on the exact imbibition solution becomes sensitive. Nevertheless, the difference observed can be neglected if compared to other sources of uncertainty in reservoir engineering.

The geometrical method was also validated for a distribution of square, triangular and U-shaped blocks as shown in Figure 8. The reference imbibition solution for this distribution was computed as the sum of the imbibition recoveries of the constitutive blocks. The geometrical method provided two equivalent block solutions, either a rectangle or a square. Again, Figure 9 shows that both solutions give imbibition recovery curves which are very close to the exact solution of the distribution of blocks.
<figure><img src="image_10.png" /><figcaption>Figure 8
Distribution of square, triangular and U-shaped blocks.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9
Water-oil imbibition in the distribution of blocks and in the equivalent square or rectangular block.</figcaption></figure>

In conclusion, despite simplifying assumptions, the geometrical method provides equivalent block solutions which are reliable in terms of multiphase flow behavior.

## 3 APPLICATION

For demonstration purposes, the same procedures were extensively applied to a large fractured rock volume generated with our in-house stochastic fracture simulator. The geological information underlying this 3D demonstration image is derived from an outcrop study of a fractured Devonian sandstone outcropping in South Morocco. Fracture characterization of this formation revealed two main fracture sets, regional systematic joints and fold-related joints, both of which are associated with compressive tectonic events. Other sets related to later tectonic events were also identified. The 3D image, shown in Figure 10, has x, y and z dimensions equal to 100 m, 100 m and 16.6 m, is made up of 9 layers, and involves only the two main fracture sets referred to above. The first fracture set, oriented along the x axis, corresponds to regional systematic joints and includes about 200 fractures with a horizontal extension of 30 to 40 m. The second set, orthogonal to the previous one, is made up of fold-related joints, among which we can distinguish one subset of some 20 fractures crossing all layers and

extending over the whole domain, and another subset of more than 1000 small fractures. The whole fracture network is input into our software as 7120 rectangles, i.e. fracture plane elements, with 460 rectangles representing the first set, and 180 and 6480 rectangles representing the two subsets of the second set.

Equivalent parameters were computed on centered sub-volumes with areas increasing from 5 m x 5 m to 100 m x 100 m. Thus, the representativeness threshold for either equivalent parameter could be determined. The mixed boundary conditions described above were applied to determine the main horizontal flow directions. The conductivity was assumed equal to 1 mD·m for all fractures of both sets. The discretization of the whole fracture network involves 90 600 nodes. It is important to realize that the regular gridding of such a fractured medium would have required several million cells.

### 3.1 Equivalent Permeabilities

Results along the three axes are shown in Figure 11 as a function of sub-volume dimension. The equivalent permeability along the y axis is about four times the permeability along the x axis which is consistent with the much higher fracture density observed in the y direction. The principal horizontal flow directions are shown in Figure 12. For all sub-volumes, the minimum
<figure><img src="image_12.png" /><figcaption>Figure 10
3D view of the demonstration image.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11
Equivalent fracture permeabilities in the demonstration image. Variation with sub-volume size.</figcaption></figure>

<caption>Figure 12
Principal horizontal flow directions in the demonstration image. Variation with sub-volume size.</caption>
and maximum horizontal flow directions are oriented respectively along the x and y axes as expected, the results differing by less than ±5° from one sub-volume to another. Consequently, the equivalent permeabilities along these principal directions do not significantly differ from those along the x and y axes. The variation of the horizontal permeability anisotropy ratio, $k_y/k_x$, and of the vertical to horizontal permeability ratio, $k_z/\sqrt{k_xk_y}$, is shown in Figure 13. Whereas the permeability anisotropy is significant within layer planes, the vertical to horizontal permeability anisotropy is very small, probably due to the predominant role played by the subset of fractures crossing the entire 3D domain in both the y and z directions.

The variation of equivalent permeability values as a function of sub-volume size (Fig. 11) reveals results that vary significantly for horizontal dimensions of subvolume lower than 50 m, and less so beyond. This lateral dimension of 50 m is an order of magnitude of the threshold scale for equivalent permeability to be representative, which is usually denoted as the Representative Elementary Volume (REV).

<figure><img src="image_14.png" /><figcaption>Figure 13
Equivalent fracture permeability anisotropy in the demonstration image. Variation with sub-volume size.</figcaption></figure>

## 3.2 Equivalent Blocks

The horizontal dimensions of the equivalent blocks of each layer for the whole area studied are given in Table 3 and schematized in Figure 14. Square equivalent blocks were found for each layer. Equivalent block size is more or less in proportion to layer thickness, which is consistent with the correlation between fracture density and layer thickness that was taken into account for fracture generation.

<table><caption>TABLE 3
Equivalent blocks in the demonstration image</caption><thead><tr><th>Layer</th><th>Thickness (m)</th><th>Equivalent block size (m)</th></tr></thead><tbody><tr><td>1</td><td>2.4</td><td>4.8</td></tr><tr><td>2</td><td>1.8</td><td>4.3</td></tr><tr><td>3</td><td>5.8</td><td>8.6</td></tr><tr><td>4</td><td>1.3</td><td>3.3</td></tr><tr><td>5</td><td>2.7</td><td>5.7</td></tr><tr><td>6</td><td>0.8</td><td>2.3</td></tr><tr><td>7</td><td>0.3</td><td>0.8</td></tr><tr><td>8</td><td>0.9</td><td>2.6</td></tr><tr><td>9</td><td>0.6</td><td>1.9</td></tr><tr><td>1 to 9</td><td>16.6</td><td>5.2</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Figure 14
Equivalent blocks in the different layers of the demonstration image.</figcaption></figure>
The variation of the equivalent block dimension as a function of sub-volume size is shown in Figure 15 for layers 3 and 7, and for the nine layers together. Layers 3 and 7 are respectively the thickest and least fractured layer, and the thinnest and most fractured layer, as illustrated by Figures 16 and 17. The equivalent block dimension of layer 7 is between 0.73 and 0.92 m, whatever the sub-volume considered. In this highly-fractured layer, small equivalent blocks are found and solutions are representative even for a very small lateral sub-volume size, down to 15 m. On the contrary, in layer 3, the equivalent block dimension increases from 3 to 8 m when the sub-volume lateral dimension increases in the range of 5 to 60 m. Beyond 60 m, the equivalent block solutions remain between 8.1 and 9.0 m. For all layers taken together, the variation of solutions as a function of sub-volume dimensions is attenuated but remains sensitive for sub-volume horizontal dimensions lower than 60 m. The general trend we observed is that the representative area for determining the equivalent block of a given layer decreases when the density and interconnection of fractures increase, that is, when the size of equivalent blocks decreases.

<figure><img src="image_16.png" /><figcaption>Figure 15
Equivalent blocks in the demonstration image. Variation with sub-volume size.</figcaption></figure>

In conclusion, representative solutions for equivalent permeabilities and equivalent blocks are found for similar minimum areas of the fractured volume. However, a separate determination of the REV's for these two parameters is recommended in the general case, because these two REV's are not directly related and depend on fracture attributes.

The high performance of our methodology in terms of computation times must also be underlined. With an Ultra-Spark 1-140 Sun workstation, the computation of the equivalent permeability tensor for the previous complete 3D network, involving 90 600 nodes, requires only 7 minutes of CPU time. The determination of an equivalent block is still more rapid. For instance, the image of a geological layer discretized with 5 million pixels is processed in 10 to 15 seconds of CPU time.

## 4 RECOMMENDATIONS FOR DUAL-POROSITY MODELING

The methodology presented here and the corresponding software have been designed for fractured reservoir engineering. Within this framework, the following recommendations can be made. First, a careful geological analysis of fracture attributes should be made to delineate reservoir areas with homogeneous fracturing criteria.

The different scales of fracturing also have to be identified. For instance, large faults which obviously cannot be homogenized within a limited reservoir area will not be included in the equivalent model but represented in a discrete way.

For each zone, a study of representativeness must be made both for equivalent fracture permeability and equivalent block size.

Equivalent parameters computed on several representative images of each area will then constitute the cell input of the dual-porosity reservoir model. As for the choice of boundary conditions for equivalent permeability computations, the zero-flow-rate lateral condition should be used because the extra-diagonal terms of the fracture permeability tensor are not taken into account in most simulators. However, additional computations with mixed boundary conditions should be performed to determine the main horizontal flow directions within the reservoir. If these directions are relatively constant throughout the reservoir, then the reservoir grid may be oriented accordingly.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16
Top view of layer 3 in the demonstration image.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 17
Top view of layer 7 in the demonstration image.</figcaption></figure>
For complex reservoirs, with continuously-changing fracture attributes, systematic computations of equivalent parameters may have to be performed on images equally distributed over the entire reservoir.

On an intermediate scale between fractures and large faults, sub-seismic faults remain a difficult issue because their dimensions are of the order of magnitude of reservoir grid blocks, which precludes treating them as part of the fracture network used to define the equivalent fracture medium at grid block scale. Nor can they be represented as discrete objects, because of the stochastic approach used to model them. Different stochastic realizations may actually look equivalent at reservoir scale but may lead to contrasting flow predictions on the scale of a few cells especially in the presence of wells. Therefore, the solution for dealing with subseismic faults might be to constrain their stochastic modeling by using other information derived from wells, for instance. Thus, a discrete representation of these faults could be achieved, at least in those reservoir locations where their effect on flow is drastic, especially in the vicinity of wells. Far from wellbores, it might be possible to integrate them into an equivalent fracture medium defined on a scale larger than cell size.

Another limitation concerns low-density fracture networks. The use of an equivalent fracture medium is not adapted to simulate flows through poorly-connected fractures. First, a scale of representativeness of flow properties is difficult to find near the percolation threshold of a fracture network. Second, if fractures are actually disconnected, a single medium representation should be used with equivalent permeabilities taking into account both matrix and fracture flow contributions.

Finally, one must recall the necessity of modeling the capillary continuity of the matrix medium in reservoirs with fractures sub-orthogonal to layer limits. To this end, the use of a dual-porosity dual-permeability model is recommended, with a proper vertical discretization that simulates vertical multiphase equilibrium with sufficient accuracy.

## CONCLUSION

A rapid and reliable methodology has been designed and implemented to determine systematically the equivalent parameters of dual-porosity models from geological images. This approach and the corresponding software can be applied to large geological

sections, on the order of magnitude of reservoir grid blocks. The accuracy of solutions is quite satisfactory for reservoir engineering applications.

This tool contributes to the integration of disciplines by facilitating the systematic use of fracture characterization data in the downstream activities of reservoir exploitation.

Further developments of this methodology are planned. They include an extension to complex fracture networks with various dips, the integration of dynamic information on the reservoir such as well tests and the development of correlations between fracturing attributes and equivalent parameters.

## NOMENCLATURE

* A normalized cross-section area of block(s)

* a, b equivalent block dimensions

* e fracture aperture

* k equivalent fracture permeability

* L flow length

* p pressure

* q flow rate

* R oil recovery from matrix block(s)

* t time

* T transmissivity

* v fluid velocity

* X distance of water-oil front from block limits

* μ fluid dynamic viscosity.

## SUBSCRIPTS

* i, j nodes of the discretized fracture network

* B boundary of the 3D fractured volume

* x, y, z directions of flow.

## ACKNOWLEDGMENTS

This research program has been sponsored by an industrial consortium that includes *Aramco Overseas Company B. V.*, *Elf Aquitaine Production, Enterprise Oil PLC*, *Qatar General Petroleum Corporation and Total Exploration Production*. The authors wish to thank the members of this consortium for their technical and financial support.
REFERENCES

1 Long J.C.S., Remer J.S., Wilson C.R. and Witherspoon P.A. (1982) Porous media equivalent for networks of discontinuous fractures. *Water Resources Research*, 18, 3, 645-658.

2 Warren J.E. and Root P.J. (1963) The behavior of naturally fractured reservoirs. *SPE Journal*, 245-255.

3 Long J.C.S., Gilmour P. and Witherspoon P.A. (1985) A model for steady fluid flow in random three-dimensional networks of disc-shaped fractures. *Water Resources Research*, 21, 8, 1105-1115.

4 Cacas M.C. et al. (1990) Modeling fracture flow with a stochastic discrete fracture network: calibration and validation. *Water Resources Research*, 26, 3, 479-500.

5 Massonnat G. and Manisse E. (1994) Modélisation des réservoirs fracturés et recherche de paramètres équivalents:

étude de l'anisotropie verticale de perméabilité. *Bull. Centres Rech. Explor.-Prod. Elf Aquitaine*, 18, 171-209.

6 Lough M.F., Lee S.H. and Kamath J.(1996) A new method to calculate the effective permeability of grid blocks used in the simulation of naturally fractured reservoirs. *Paper SPE 36730* presented at the *SPE Annual Technical Conference and Exhibition*, Denver, Co.

7 Odling N.E. (1992) Permeability of natural and simulated fracture patterns, structural and tectonic modelling and its application to petroleum geology. *NPF Special Publication* 1, 365-380, Elsevier, Norwegian Petroleum Society (NPF).

8 Bourbiaux B. and Kalaydjian F. (1990) Experimental study of cocurrent and countercurrent flows in natural porous media, *SPE Reservoir Engineering*, 361-368.

Final manuscript received in October 1998