3D reconstitution of porous media from image processing data using a multiscale percolation system

Journal of Petroleum Science and Engineering 42 (2004) 15 – 28 www.elsevier.com/locate/petrol

3D reconstitution of porous media from image processing data using a multiscale percolation system J.-F. Daı¨an a, C.P. Fernandes b, P.C. Philippi b,*, J.A. Bellini da Cunha Neto b a b

Laboratoire d’e´tude des Transferts en Hydrologie et Environnement (UJF, INPG, CNRS, IRD) BP 53, F38041 Grenoble cedex, France Porous Media and Thermophysical Properties Laboratory, Mechanical Engineering Department, Federal University of Santa Catarina, 88040-900 Floriano´polis, SC, Brazil Received 19 May 2003; accepted 16 October 2003

Abstract The concept of multiscale percolation system (MPS) is proposed as a tool for the reconstitution of a 3D idealized pore space. The statistical properties of the MPS are generated on the basis of the size distributions of pores and solid particles of the porous medium under investigation that were obtained by microphotograph processing. The generation process is based on upscaling by means of renormalization. Renormalization is also used for computing MPS intrinsic permeability. The method is applied to a series of 23 reservoir sandstones which experimental permeability vary in a range 1:400. Permeability values computed according to the MPS model appear to be well correlated with experimental values. D 2003 Elsevier B.V. All rights reserved. Keywords: Permeability; Image analysis; Percolation theory; Renormalization

1. Introduction Two- and three-dimensional representations of porous media and, particularly, of reservoir rocks in petroleum science, have been used for simulating fluid displacements, single phase or multiphase flow, and for the estimation of transport properties. These representations require to transform the microphotographs of sections of a porous medium sample to a 3D structure by means of a given reconstitution process. A first class of methods uses the 2D image for extracting statistical descriptors of the pore space * Corresponding author. Fax: +55-48-234-1519. E-mail address: [email protected] (P.C. Philippi). 0920-4105/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.profnurs.2003.10.002

morphology like autocorrelation functions and the corresponding first- and second-order moments. These descriptors may be directly used to search empirical correlations with the transport properties (Ioannidis et al., 1996). Alternatively, the statistical properties obtained from a cross section may be used to constrain the generation of a 3D pore space (Quiblier, 1984; Thovert et al., 1993; Talukdar et al., 2002). Direct simulation of 3D flow in the reconstructed pore space may be performed using Navier – Stokes equations (Adler et al., 1990; Bentz and Martys, 1994) or lattice gas/lattice Boltzmann simulations (Ferreol and Rothman, 1995; Santos et al., 2002). Other authors previously perform a skeletonization of the reconstructed pore space before modeling transport (Liang et al., 2000a,b).

16

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

The second class of models is based on reconstitution methods in the form of pore networks (Diaz et al., 1987; Ioannidis and Chatzis, 1993; Philippi and Souza, 1995) or random generation of the pore space, with correlation or not (Ioannidis et al., 1997) using the basic laws of Percolation Theory. The present paper proposes a reconstruction method based on the concept of multiscale percolation system (MPS) introduced by Neimark (1989). This concept was used for the reconstruction of porous media with a large pore size distribution in relation with the correlation function of the pore space (Fernandes et al., 1996). The concept of MPS was used elsewhere for reconstituting an idealized pore space structure from mercury intrusion porosimetry data (Xu et al., 1997a,b). For this model, which will be referred as XDQ in the present paper, the original concept of MPS was modified and simplified. While MPS concept involves two size distributions, for the pores and for the solid, XDQ model considers only the distribution of pore sizes, since mercury porosimetry provides no information on the structure of the solid phase. In the present paper, the main interest is to reconstitute the pore space from the data obtained by means of image processing. In this way, considering that image processing enables to characterize the geometrical structure of both phases, it is convenient to come back to the original MPS concept involving two distributions, since this concept provides a more complete description of the geometrical three-dimensional arrangement of pore and solid features.

2. Image processing data In this paper, the data used in the multiscale reconstitution of the porous medium are provided by image analysis. The reconstitution process is applied for the case of microstructures of petroleum reservoir rocks. For this, digital color images of core samples are derived from optical microscopy. Thin plates are prepared by introducing a colored resin into the porous structure, previously evacuated, and the porous phase is identified as the geometrical region related to the particular resin’s color. Thus, solely the pores that are accessible to fluid flow are identified, i.e., the open pore space. An example of a digital color image obtained from optical microscopy is shown in Fig. 1a.

The data needed for the multiscale reconstitution are the pore and solid size distributions. To quantify these distributions, the original color images are transformed to binary images (Fig. 1b) defining the geometric regions related to the pore and solid phases. This transformation is realized, manually, using three adequate thresholds, one for each component histogram in the images coded in the Hue – Saturation – Intensity (HSI) color model (Gonzalez and Wood, 1992). The pore size distribution of the resulting 2D binary sections is obtained by successive opening, derived from mathematical morphology and using balls with increasing radius (Coster and Chermant, 1989). The resulting image can be considered as the union of balls completely enclosed in the porous phase. Thus, after each opening operation, the porous phase that is lost is the union of all the features that are smaller than the opening ball. The cumulative porous distribution is, then, given by: FðrÞ ¼

/  /ðrÞ /

ð1Þ

where / is the total porosity of the original image and /(r) is the volume fraction of the porous phase after the opening operation with a radius-r ball. To reduce time processing, the opening operation is not applied directly to the binary image, but to a transformed image called background distance image (Chassery and Montanvert, 1991). In this image, each pixel is labeled with the smaller distance from it to the neighbor background. This labeling technique uses a sequential algorithm where the Euclidean distance is approximated by a discrete integer distance (Chassery and Montanvert, 1991; Laurent and Frendo-Rosso, 1992). The most commonly used discrete distance is the chamfer distance known as d3 – 4. The solid size distribution is obtained in the same manner, by inverting the binary image.

3. Definition of a multiscale percolation system (MPS) 3.1. Multiscale fragmentation Following Neimark’s concept of multiscale percolation system (MPS) (Neimark, 1989; Fernandes et al.,

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

17

Fig. 1. Microphotograph (a) of a thin plate of a reservoir rock (S23) and the corresponding binarized image (b). The image corresponds to 2.91  2.18 mm. The blue resin appears in soft grey in the original image and in black in the binary image.

1996), we consider n levels of fragmentation of the space, referred as levels 1 to n. Level 1 is associated with the largest blocks of size a1, level n with the finest ones, of size an. In this paper, the size ratio (ai/ ai + 1) is uniformly equal to 2. Fig. 2 illustrates the building process of a MPS.

At the finest rank of fragmentation, the blocks of size an are randomly occupied either by void (pores of size an, or of class n) with probability pn, or by solid (particles of size an, or of class n), with probability sn = 1  pn. At the following ranks i, blocks of size ai = 2n  ian are, successively, superposed. They can be

18

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

networks. The two approaches are developed in this paper. Blocks of the several classes, defined above, are considered as elements (sites or bonds) of constitutive networks of ranks 1, . . ., n. For each constitutive network of mesh size ai, these elements can be decided to be pores with probability pi, solid particles with probability si, or remain available for occupation with finer elements belonging to the classes i + 1, . . ., n. 3.2. Characterization of pore connection by means of renormalization and upscaling Fig. 2. A 2D representation of a MPS.

void features with probability pi, solid blocks with probability si, or remain available for the previous occupation by pores and particles of classes i + 1 to n (Fig. 2). At rank i, the space fraction which was not occupied either by pores or by particles of class i is: qi ¼ 1  ðpi þ si Þ

ð2Þ

The process is continued until rank 1. In the final structure obtained, the visible space fraction occupied by pores of classes i = 1,. . .n is: u1 ¼ p 1

ui ¼ pi qi1 qi2 . . . q1 ð3Þ Pn and the total porosity is i¼1 ui . Similarly, the following visible solid space fractions can be defined: w1 ¼ s1

wi ¼ si qi1 qi2 . . . q1

In classical percolation, networks are randomly occupied by active elements (sites or bonds) and the fraction of elements which are connected at long distance is called the ‘‘strength of the infinite cluster’’ (Stauffer, 1985; Stauffer and Aharony, 1994), which can be denoted Y( p), where p is the concentration (or probability of presence) of the active elements. The infinite cluster can be, also, characterized by the ratio y( p) = Y( p)/p. The main features of function Y( p) are well known in percolation theory, particularly, concerning its behavior close to the percolation threshold. Representations of this function for a cubic bond network and for a cubic site network based on percolation theory and on numerical simulations (Xu et al., 1997a) are shown in Fig. 3.

ð4Þ

In view of studying fluid transport at the macroscopic scale, the problem is to find the portion of the pore space that is connected at long distance, and the contribution vi of each class of pores to this part of the pore space. The study of connection will be, here, considered inside the percolation theory framework. Percolation networks may be of two kinds, because the blocks of the MPS can be considered either as the sites, or as the bonds of a cubic network. For describing space coverage by pores and solid particles, the site approach is more convenient. On the contrary, having in view the transport phenomena occurring in the pore space, it is more convenient to use bond

Fig. 3. The connected fraction (infinite cluster) of a 3D cubic bond or site network. Percolation threshold (critical concentration): pc = 0.307 for site cubic networks, pc = 0.25 for bond cubic networks.

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

In the case of a MPS, the function Y( p) which defines the connected fraction of the active elements cannot be used directly, since the concentration of active elements cannot be immediately defined by a single number like p. All the classes contribute to the pore space and the concentration must be defined as a function of the various contributions pi. Renormalization or rescaling methods may be used to overcome this problem. Following XDQ renormalization (Xu et al., 1997a), we admit that a network of uniform mesh size a with a concentration of active elements p is equivalent, from the viewpoint of connection, to a network of mesh-size 2a with a concentration f( p), provided that the correlation length is the same in the two networks. Like any renormalization method (Stauffer, 1985; Sahimi, 1994), XDQ renormalization does not preserve all the properties of a network. Nevertheless, we have a fixed point at exactly the percolation threshold pc, that is, f( pc) = pc, which is not necessarily the case for traditional small cell renormalization. In this way, XDQ method allows to perform a consistent reconstitution of a porous medium, particularly useful in the case where it is composed of pores and particles of widely varying sizes. Xu et al. (1997a) studied the conditions that the renormalization function f must satisfy to preserve the correlation length. In particular, these conditions must include the existence of a fixed point at the percolation threshold, pc. The shape of the renormalization function f( p) proposed for bond percolation and for site percolation in a cubic network is shown in Fig. 4. The renormalization method enables to build up the connected cluster of pores of a MPS by means of upscaling, in the following way. The network of rank n, of smallest mesh size an, occupied by pores with concentration pn, can be considered as equivalent to a network of mesh size an  1 = 2an occupied by active elements (fictitious pores) with concentration f ( pn). Superposing the pores and the solid particles of class n  1 yields the following concentration of active elements in the network of rank n  1, which accounts for the contributions of pores of the classes n and n  1: pn1 ¼ pn1 þ qn1 f ðpn Þ

ð5Þ

19

Fig. 4. The XDQ renormalization functions for a cubic bond network and for a cubic site network.

The process can be repeated, and a concentration can be defined in each network: pi ¼ pi þ qi f ðpiþ1 Þ

ð6Þ

At rank 1, the pores of all classes are integrated in the structure, yielding the concentration p1, in the network of rank 1. 3.3. Volumetric contributions to the connected porosity at long distance The concentration p1 takes all classes of pores into account. Consequently, the strength of the infinite cluster is Y(p1). This multiscale infinite cluster contains pores of class 1 and active elements representative of the others classes. The respective proportion of these two types of active elements is the same in the concentration p1 and in the infinite cluster. Consequently, the contribution of pores of class 1 to the infinite cluster, which represents the porosity connected at long distance is: v1 ¼

p1 Y ðp1 Þ ¼ p1 yðp1 Þ ¼ u1 yðp1 Þ p1

ð7Þ

To determine the volumetric contribution of anyone of the other classes, i, it must be observed that the renormalization used at each stage of the upscaling process does not preserve either the strength of the

20

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

infinite cluster or the ratio y( p) = Y( p)/p. To account for this, Xu et al. (1997a) introduced a second renormalization function g( p) such that: gðpÞ ¼

yðpÞ y½f ðpÞ

ð8Þ

This function must also be defined for sub-critical concentrations, since the partial concentrations pi may be smaller than the percolation threshold pc. Xu et al. (1997a) studied the requirements for the function g( p) for p < pc in order to preserve the statistical density of finite clusters in a sub-critical network. Using this function at each stage of the upscaling process, the volume contributions of any class of pores to the connected porosity of a MPS is obtained: vi ¼ pi qi1 qi2 . . . q1 gðpi Þgðpi1 Þ . . . gðp2 Þyðp1 Þ vi ¼ ui gðpi Þgðpi1 Þ . . . gðp2 Þyðp1 Þ

ð9Þ

4. Identifying the pore and particle concentrations from image processing data 4.1. From the image of a porous section to a 3D MPS As mentioned in Section 2, the area fraction of the section under analysis identified as pores corresponds to the part of the space which has been intruded by resin, i.e., the connected part of the pore space or open porosity. The material under analysis may contain a closed pore space or not, but in any case, solely the connected pore space is accessible to resin, contributing to the porosity. Image processing yields the size distribution of the connected pore space. This distribution corresponds, in the multiscale concept proposed above, to the size distribution of the connected pore space defined by the specific volumes vi. The part of the image where resin is not visible is identified as solid, but in fact, it corresponds to both solid volumes and unconnected porosity if it exists. The size distribution of this area corresponds therefore, in the multiscale modeling, to the sum: wap i ¼ wi þ ðui  vi Þ

ð10Þ

which will be referred as the apparent solid fraction.

For interpreting a given image in terms of a 3D MPS, we accept the following basic hypotheses: for each species (pores, connected or not, solid, real or apparent), the size distribution with respect to the areas which appear in a section is identical to the size distribution with respect to the volume in the material. This condition is clearly satisfied for the section of a random structure made of embedded cubic blocks such as the one given in Fig. 2. We assume that it can be extended to the case of blocks having no particular shape. Moreover, these size distributions with respect to the volumes are supposed to provide sufficient information on the geometrical structure of the analyzed material, for estimating properties, such as the connection degree of the pores and the role of the various pore sizes in connecting the pore space. These connection properties are essential for the study of fluid invasion and transport. It is clear that these assumptions may be more or less valid, depending on the nature of the material studied. In particular, for materials where the phases present a more or less deterministic spatial organization, including anisotropy or concentration of pores of a given size in particular locations, the multiscale system here proposed, which is, essentially, based on randomness, will probably fail. A given pore size distribution with respect to the volumes may give rise to very different values of the transport properties, depending on the spatial arrangement of the pores. However, this type of problem arises for any reconstitution method that aims some generality. In this way, random generation of porous structures appears to remain the best procedure, when no clear deterministic and quantifiable organization is identifiable in the images. 4.2. Identification process According to the previous assumptions, image processing is used to give, simultaneously, the volumetric pore and the solid particles size distributions. These distributions can be rearranged in order to define n sizes satisfying ai = 2ai + 1. Therefore, the 2n specific volumes obtained in this way, in each distribution, are to be considered as the specific volumes vi and wiap of the representative multiscale structure to be determined.

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

In this manner, the identification of the representative MPS consists in solving an inverse problem involving 2n unknowns, i.e., the values of pi and si. It can be seen that all the (ui + wi ) values (Eq. (10)), are directly known from the disposable data. Adding Eqs. (3) and (4) yields: ðu1 þ w1 Þ ¼ ðp1 þ s1 Þ ðui þ wi Þ ¼ ðpi þ si Þqi1 qi2 . . . q1

21

converges. However, the size distribution of the connected pores, vi, is preserved, since this property is supposed to be the most determinant factor for the permeability.

5. Computing transport coefficients ð11Þ

Therefore, ( pi + si) can be easily determined from i = 1 to i = n, as well as qi (from Eq. (2)). For determining pi, an iterative method is required. The algorithm is outlined in the following. An arbitrary value of p1,p1* is adopted. Eq. (7) directly yields p1. The value of qi being known Eq. (6) for i = 1, yields p2, owing to an inversion routine for function f. Then, Eq. (9), for i = 2, can be used for determining p2. In this way, all values pi can be successively determined, until pn. Using Eq. (6) again, the values of pi can be successively calculated, starting from i = n, until i = 1. This yields a new value for p1,p1** to be compared with the initial value p1*. The solution is reached when p1** is close enough to p1*. It can be easily shown that p1** decreases with p1* and dichotomy can be used to solve the problem. The evidence that the problem involves 2n unknowns for 2n equations is an argument for the uniqueness of the solution, but does not prove its existence. In fact, the algorithm does not converge, sometimes, and it can be concluded that no random MPS has the values of the vi and wiap of the source image. In such cases, we modify the distribution of the apparent solid wiap, performing successive modifications on the initial distribution until the algorithm

The representative multiscale structure being completely identified by the values of pi and si, the procedure for the computation of the various transport coefficients described by Xu et al. (1997b) can be applied without major modification. The difference lies in the values of qi to be used in the upscaling process (Section 3, Eq. (5)), since the version of XDQ model for mercury intrusion data does not account for the presence of solid particles except at rank n. The XDQ method, based on the renormalization method developed by King (1989) for site networks composed of cubic blocks and elsewhere by Hinrichsen et al. (1993) for cubic bond networks, closely follows the upscaling process described in Section 3.2, and uses the values of the concentrations pi. At each step of the upscaling process, it is possible to compute a distribution of conductivities. In the network of rank n, each bond or site representative of a pore has a conductivity rn. A cell of size 2an (Fig. 5) randomly composed of elements of conductivity rn with probability pn, and of elements of zero conductivity with probability 1  pn, can be considered as an element of the renormalized network of mesh size 2an = an  1. The equivalent conductivity of the cells can, then, be computed, leading to the distribution of the conductivities of the renormalized elements. In the renormalized network, the non-zero conductivities of this distribution are ran-

Fig. 5. Two-dimensional representations of the unit cells used in King’s (1989) renormalization and Hinrichsen et al. (1993) renormalization.

22

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

domly distributed over the fraction f( pn) of the sites or bonds. Superposing the pores of class n  1, which have the conductivity rn  1 and probability pn  1 and the solid particles of zero conductivity and probability sn  1, leads to a modified distribution of the conductivities in the network of rank n  1. New cells of size 2an  1 can be randomly defined in this network, and their conductivity distribution can be computed. The upscaling process can be continued until the elements of rank 1 are integrated and the renormalization process can be performed without adding new classes of pores and particles. This produces a conductivity distribution more and more narrow, which is supposed to converge towards a Dirac distribution centered at the macroscopic conductivity when the largest scale of the MPS is reached. Transport coefficients of different nature can be obtained, considering the particular models for the conductivities ri. In the present paper, we are interested in the intrinsic permeability. The size of the pores must be defined in order to use Poiseuille’s law for the definition of ri. Considering the class i of a site MPS, formed with cubic blocks of size ai, laminar flow is considered to occur in each cubic block through a square cross section a i2. This yields the following conductivity, according to Poiseuille law for a square cross section: ri ¼

1 2 a 28:4 i

ð12Þ

For bond networks, the elementary volume a i3 is occupied by three bonds, so the volume of each bond is a i3/3. The bonds are considered as cylindrical ducts of length ai and of cross section a i3/3. This leads to a result that seems contradictory: in a network, in which all the bonds are pores, the porous cross section available for transport in a single direction is only 1/ 3 of the total section of the sample. This paradox has been discussed by Xu et al. (1997b). We may note that Hinrichsen et al. (1993) do not consider the factor 1/3 for the cross section. Hinrichsen’s assumption is more consistent from the point of view of transport, but it is incompatible with the distribution of the pore volumes, considered as network bonds, over three spatial directions.

In spite of these uncertainties, the permeability of bond MPS was computed using the following pore conductivity: ri ¼

1 1 2 d 3 32 i

ð13Þ

The factor 1/32 corresponds to Poiseuille’s law for circular cylinders and di is the mean diameter of the balls of class i identified in the image.

6. Application to reservoir rocks 6.1. Morphologic characterization Table 1 gives a roll of the 23 reservoir-rocks from, which thin plates of about 10 Am thickness have been treated. Ten microphotographs have been taken for each plate. The average properties of each rock have been determined from the analysis of these 10 images. The image shown in Fig. 1a was obtained by transmission of unpolarized light. Three types of regions can be seen. The transparent mineral grains clearly appear in white and the blue resin appears in soft grey. Dark regions between the mineral grains are also visible in the image, but they do not appear in the binary image, since they are not colored in blue. The interpretation of these dark spots is doubtful. They may be due to dead oil or to the presence of pyrites. It is also possible that the penetration of resin into the pores is incomplete. Another problem that may occur when interpreting the images is the difference between the porosity determined from image analysis (optical porosity), and the porosity measured by mercury intrusion or gas expansion (experimental porosity). This difference is sometimes important, particularly for rocks of low porosity and can be due to several factors. Microporous minerals can be present, and these minerals may have been intruded or not by the resin. Even if resin has penetrated into the microporosity, it may remain invisible due to the image resolution. The intrusion of resin may also be incomplete, even in pores of visible size if their entry is controlled by small micropores. In any case, the part of the porosity which does not appear as visible resin in the images, even if rela-

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

23

Table 1 Experimental and computed properties of the 23 reservoir rocks studied Sample

/opt

n

ip

d1 (Am)

dpmax (Am)

dpmin (Am)

ics

dcs (Am)

icb

dcb (Am)

Kexp (mDa)

KXDQs (mDa)

KXDQb (mDa)

S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 S21 S22 S23

0.037 0.059 0.065 0.072 0.076 0.095 0.095 0.100 0.109 0.113 0.128 0.135 0.140 0.152 0.155 0.156 0.159 0.182 0.182 0.182 0.183 0.192 0.206

7 6 6 6 6 6 7 6 6 7 7 6 6 6 6 7 6 7 7 6 6 6 5

4 3 2 3 3 2 3 3 2 3 2 2 2 2 2 3 2 2 2 2 2 2 2

430 341 263 347 379 251 469 302 328 495 301 206 231 244 321 437 334 488 331 230 296 225 170

54 85 132 87 95 126 117 76 164 124 150 103 116 122 161 109 167 244 166 115 148 113 85

6.7 10.7 8.2 10.8 11.8 7.8 7.3 9.4 10.3 7.7 4.7 6.4 7.2 7.6 10.0 6.8 10.4 7.6 5.2 7.2 9.3 7.0 10.6

7 6 6 6 6 6 7 6 6 7 7 5 6 6 6 7 6 7 6 6 6 6 5

6.7 10.7 8.2 10.8 11.8 7.8 7.3 9.4 10.3 7.7 4.7 12.9 7.2 7.6 10.0 6.8 10.4 7.6 10.4 7.2 9.3 7.0 10.6

7 6 6 6 6 6 7 6 6 7 5 5 6 6 6 6 6 6 5 5 6 5 5

6.7 10.7 8.2 10.8 11.8 7.8 7.3 9.4 10.3 7.7 18.8 12.9 7.2 7.6 10.0 13.7 10.4 15.3 20.7 14.4 9.3 14.1 10.6

1.8 3.6 9 23 14.7 26 34.7 37 68 37.2 154 104 34.7 92 68 152 408 197 316 70 540 145 642

0.73 5.79 4.83 9.83 12.8 13.2 14.4 19.2 31.8 20.7 32.1 72 36.9 52.9 79.5 69.4 92.4 100 132 104 140 183 247

0.09 1.02 0.95 1.65 2.28 4.6 5.9 48.8 16.3 10.5 31 36.5 21 32 48 42.5 49 75 126 85 272 132 133

/opt: optical porosity. n: total number of classes. ip: class of the largest pores. d1: typical diameter of class 1. dpmax: typical diameter of largest pores (class ip). dpmin: typical diameter of smallest pores (class n in all cases). ics, (icb): critical class for site (bond) MPS. dcs (dcb): critical pore diameter for site (bond) MPS. Kexp: experimental permeability. KXDQb: permeability computed by XDQ model using a bond MPS. KXDQs: permeability computed by XDQ model using a site MPS.

tively important, represents regions of the pore space of poor accessibility for fluids. Therefore, ignoring this porosity in analyzing the images and in modeling should not induce important errors in permeability estimation. The optical porosity, /opt (Table 1) is widely varying in the range 1:5, with a minimum value of 3.4% and a maximum of 20.6%. The optical porosity is, generally, appreciably smaller than the experimental porosity, directly determined from core samples. The size distribution obtained from image processing was rearranged in n classes, considering that MPS typical diameters vary by a factor 2. A class i is defined in the interval [Di, 2Di] and pffiffiffi the typical diameter associated to this class is Di 2. Two examples of volume size distribution for the pores impregnated by resin (open porosity) and for the apparent solid (particles and, possibly, closed pores) obtained from image processing, can be seen

of Tables 2 and 3. For in the columns vi and wapp i all the samples, global information on the size distribution are given in Table 1. Generally, the number of classes is 6 or 7. The first classes (largest objects) do not present pores, but only solid particles. 6.2. Reconstitution as a MPS and critical diameter Tables 2 and 3 illustrate the reconstruction procedure using sites and bonds MPS, described in Section 4 for, respectively, the less porous sample S1, with porosity 3.4% and for the most porous one S23, with porosity 20.6%. The values of the 2n basic pore and solid concentrations ( pi,si) of the representative MPS are given together with the corresponding pore and solid volumes (ui,wi) and the concentrations of active elements pi involved at each stage of the upscaling process. The concentration pi varies irregularly from rank n to rank 1, since according to Eq. (5), renorm-

24

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

Table 2 Reconstitution of S1 (porosity 3.4%) Image distributions

Site reconstitution

Bond reconstitution

i

vi

w iapp

ui

wi

pi

si

pi

ui

wi

pi

si

pi

1 2 3 4 5 6 7

0.0000 0.0000 0.0000 0.0030 0.0093 0.0072 0.0140

0.2230 0.4085 0.2113 0.0624 0.0265 0.0178 0.0165

0.0000 0.0000 0.0000 0.0056 0.0172 0.0134 0.0260

0.2133 0.3778 0.2030 0.0684 0.0313 0.0253 0.0182

0.0000 0.0000 0.0000 0.0273 0.1307 0.1615 0.5878

0.2133 0.4803 0.4965 0.3326 0.2382 0.3048 0.4121

0.3335 0.3580 0.4741 0.6521 0.7136 0.6303 0.5878

0.0000 0.0000 0.0000 0.0053 0.0164 0.0127 0.0246

0.2230 0.4085 0.2113 0.0601 0.0195 0.0123 0.0058

0.0000 0.0000 0.0000 0.0341 0.1793 0.2294 0.8069

0.2230 0.5258 0.5736 0.3827 0.2131 0.2218 0.1930

0.2731 0.2924 0.3998 0.6142 0.7833 0.7764 0.8069

Percolation threshold (critical concentration), pc = 0.307 for site cubic networks, pc = 0.25 for bond cubic networks.

alization from a value larger than pc as well as pore addition increases the concentration, while addition of solid elements has the opposite effect. For the sample of low porosity, appreciable differences can be observed between the total pore volume of each class in the MPS, ui, and the connected part of this volume, vi. On the contrary, for the sample with higher porosity, the reconstitution process leads to very small differences: almost all the pores of the MPS are connected. For each class, the unconnected pore volume is recovered in the difference between the solid volumes wi and wapp i . The MPS concept associated with XDQ renormalization allows, also, to estimate the critical pore diameter defined by Charlaix et al. (1987) and Thompson et al. (1987). The critical pore diameter is supposed to play an important role in the permeability of a percolation pore network. These authors consider the smallest sub-network of the largest pores connected at a large scale, i.e., which concentration is above pc. The critical diameter dc is the diameter of the smallest pores present in this sub-network. This diameter controls all the main pathways of the fluids and is supposed to play a determinant role in the value

of the intrinsic permeability together with the formation factor F. This is quantified by the following formula (Thompson et al., 1987):

K~

dc2 F

ð14Þ

In a MPS, the XDQ upscaling procedure can be performed successively for sub-networks involving an increasing number of the first classes, until the obtained concentration p1 is larger than pc. The last class i integrated in the connected sub-network gives an estimation of the critical diameter, dc. This estimation is given in Table 1 in the form of the critical class index ic and of its typical diameter dc. We observe that except for some rocks of higher porosity, the connection of the sub-network is only obtained when all classes of pores are integrated (ic = n), i.e., the critical diameter coincides with the minimum pore diameter. It can be observed that the critical diameter is near constant in present rock series, varying by a factor of about 2 for site reconstitution and 3 for bond reconstitution.

Table 3 Reconstitution of S23 (porosity 20.6%) Image distributions

Site reconstitution

Bond reconstitution

i

vi

w iapp

ui

wi

pi

si

pi

ui

wi

pi

si

pi

1 2 3 4 5

0.0000 0.0082 0.0679 0.0666 0.0631

0.1278 0.4127 0.1681 0.0489 0.0362

0.0000 0.0085 0.0701 0.0688 0.0653

0.1278 0.4125 0.1658 0.0467 0.0340

0.0000 0.0097 0.1555 0.3202 0.6573

0.1278 0.4730 0.3677 0.2174 0.3426

0.6219 0.4860 0.6270 0.7573 0.6573

0.0000 0.0083 0.0680 0.0667 0.0632

0.1278 0.4127 0.1679 0.0488 0.0361

0.0000 0.0095 0.1508 0.3106 0.6364

0.1278 0.4732 0.3724 0.2271 0.3635

0.7045 0.4976 0.6231 0.7506 0.6364

Percolation threshold (critical concentration), pc = 0.307 for site reconstitution, pc = 0.25 for bond reconstitution.

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

25

6.3. Intrinsic permeability and XDQ estimation Table 1 gives the experimental values of the intrinsic permeability. Values are widely varying in the series by a factor of about 400. It is remarkable in this series that the permeability increases with the porosity, with a good correlation coefficient, as it can be seen in Fig. 6. In fact, the following empirical relationship can be proposed for this series of rocks, with a correlation coefficient of 0.92: Kexp ¼ 4410 /3:1 opt

ð15Þ

The existence of this type of correlation, without any typical pore size dependence, clearly shows the very particular nature of this rock series. Typical pore sizes, evaluated from optical analysis by dpmin and dpmax or, considering XDQ reconstitution, by the critical pore size dc, do not present any significant variation in the series. This was confirmed by performing a multilinear regression for the logarithmic values of Kexp as a function of the logarithmic values of /opt and dc. The correlation with dc was found to be insignificant. This particular character of the rock-series under consideration must be kept in mind when analyzing the performance of XDQ model in predicting the permeability. Table 1 gives the estimations of the permeability given by XDQ model, using a reconstitution either as

Fig. 6. The empirical relationship between optical porosity and permeability for the series of 23 rocks studied.

Fig. 7. Estimated permeability according to XDQ model using a site MPS.

a bond MPS with Hinrichsen’s renormalization, or as a site MPS with King’s renormalization. Figs. 7 and 8 illustrate the comparison between measured and computed values of permeability presented in a logarithmic scale. For site reconstitution, the points in the logarithmic plot present a general alignment parallel to the diagonal. The computed permeability appears to be statistically proportional to the experimental value. Considering that present model is based on a few statistical properties of the porous microstructure, permeability estimation appears to be fairly good, although the model exhibits a general tendency to underestimate the permeability. In fact, the ratio of

Fig. 8. Estimated permeability according to XDQ model using a bond MPS.

26

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

the computed and experimental values has a logarithmic average of 0.56, with a confidence interval of [0.32, 0.98], estimated from the logarithmic standard deviation of the ratio. This interval represents a standard error of 1.8. Nevertheless, for bond reconstitution, the scattering is more important, and the underestimation of permeability with respect to the experimental value is considerable in the range of low permeability. The ratio of the computed and experimental values is characterized by an average value of 0.28, which is appreciably smaller than for site reconstruction, and the wider confidence interval [0.12, 0.66] indicates a more important scattering. The reasons of this poor performance may be found in the uncertainties discussed in Section 5. This confirms that the use of site MPS is better for the reconstitution from image analysis data. For both site and bond MPS reconstruction, permeability is underestimated. A first reason for this underestimation may be related to the direct use of the size distribution obtained from a 2D image for generating a 3D cubic pavement. As it was explained in Section 4.1, in MPS reconstruction, size distribution in volume is assumed to be the same size distribution in area that was measured on the 2D source image. As it is well known, this is not necessarily true for real porous structures. In this way, consider, e.g., a porous structure with spherical pores randomly located: twodimensional sections yield circles, which diameter distribution is systematically below the diameter distribution of the spheres. In this case, permeability underestimation will be directly related to the underestimation in the block diameters used in the 3D MPS pavement. The permeability underestimation may also be due to the fact that, in the present approach, the pore space is modeled as the percolation cluster extracted from a random spatial distribution of pores over the network matrix. This generation process yields particular geometrical configurations, which present a increasing fractal character and, consequently, a high tortuosity, when the concentration approaches the percolation threshold. The reservoir rocks studied in the present paper are sandstones for which the pore space is, essentially, inter-granular, generated by a more or less compact and deterministic package of solid particles of widely varying size. This process of generation has

no reason to yield fractal geometrical structures, since it is not a percolation process. The importance of the tortuosity effects in relation with the percolation processes can be estimated, for a given XDQ reconstitution, on the basis of the value of the representative concentration of active elements in the final network, p1, compared to the percolation threshold. In the series studied, only a few samples of low porosity present a reconstitution for which concentration p1 is clearly in the region close to pc, where the percolation effects are dominant. An example is given for the sample S1 of porosity 3.4% (Table 2). This analysis suggests that the overestimation of the tortuosity effect due to the use of a method of reconstitution based on percolation should be more and more important when the porosity decreases. For MPS bond reconstitution, we, effectively, observe in Fig. 8 that the deviation of the computed permeability increases in the range of the lower values. But this is not confirmed for the case of MPS site reconstitution (Fig. 7), although the values of p1 also lie in the critical range for the samples of lowest porosity, as it can be seen in Table 2 for the sample S1.

7. Conclusion In this paper, the concept of multiscale percolation system has been used for the reconstruction of a 3D porous medium using the data extracted from the analysis of microphotographs of a cross section of the sample. The process is an extension of XDQ model originally conceived for the treatment of mercury intrusion data. The new model allows to integrate to the MPS parameters, the size distribution of the solid phase provided by image processing, while mercury porosimetry is limited to information on the pore size distribution. An important difference between XDQ model and other reconstruction processes may be emphasized. XDQ process does not effectively generate a numerical sample of the 3D reconstructed medium, but only determines a few statistical properties of the representative MPS in agreement with the corresponding statistical properties of the medium considered. Consequently, the identification process is less sophisticated than for most other models. For the estimation

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

of the permeability, XDQ also avoids the heavy process of solving the flow equations in a complex 3D pore space, using a statistical method based on the treatment of conductivity distributions and upscaling by renormalization. The representative pore space associated to the MPS model is strongly structured by its generation process. Percolation is based on the random occupation of the space previously fragmented by a regular grid. MPS model improves the concept by a multiscale fragmentation of the space, limited, however, in present case, to grids for which the size varies following powers of 2. The reconstituted pore space is obtained by extracting the connected part of the space randomly occupied. It is clear that the structure obtained is submitted to important constraints, particularly when it is close to the percolation threshold. Nevertheless, it can be noticed that the methods based on the properties of the autocorrelation function also induce limitations in the structure of the reconstructed pore space. The application of XDQ model to a series of reservoir rocks of widely varying porosity and permeability shows that the method is able to provide a reasonably accurate estimation of the permeability. Nomenclature ai Mesh size of a network of rank i (m) di Diameter of pores of class i (m) f( p), g( p) Renormalization functions K Intrinsic permeability (m2) n Number of classes in a MPS p Concentration of pores in a percolation network pc Percolation threshold (critical concentration) pi Concentration of pores of class i in a network of rank i qi = 1  ( pi + s i ) si Concentration of solid particles of class i in a network of rank i ui Specific volume of pores of class i in a MPS vi Specific volume of connected pores of class i in a MPS wi Specific volume of solid particles of class i in a MPS Y( p), y( p) connected fraction of the active elements in a percolation network

/ pi ri

27

porosity concentration of active elements in the partial network of rank i conductivity of the elements of class i

Acknowledgements The authors would like to acknowledge CENPES/ PETROBRAS (Centro de Pesquisas e Desenvolvimento Leopoldo A. Miguez de Mello) for providing the images and experimental data for reservoir-rocks. Jean-Francß ois Daı¨an gratefully acknowledges financial support of Brazilian Ministry of Science and Technology (MCT-RHAE). References Adler, P.M., Jacquin, C.G., Quiblier, J.A., 1990. Flow in simulated porous media. Int. J. Multiph. Flow 16, 691 – 712. Bentz, D.P., Martys, N.S., 1994. Hydraulic radius and transport in reconstructed model three-dimensional porous media. Transp. Porous Media 17, 221 – 238. Charlaix, E., Guyon, E., Roux, S., 1987. Permeability of a random array of fractures of widely varying apertures. Transp. Porous Media 2, 31 – 43. Chassery, J.M., Montanvert, A., 1991. Ge´ome´trie Discre`te en Analyse d’Images. Herme`s, Paris. Coster, M., Chermant, J.L., 1989. Precis d’Analyse d’Images. Presses du CNRS, Paris. Diaz, C.E., Chatzis, I., Dullien, F.A.L., 1987. Simulation of capillary pressure curves using bond correlated site percolation on a simple cubic network. Transp. Porous Media 2, 215 – 240. Fernandes, C.P., Magnani, F.S., Philippi, P.C., Daı¨an, J.-F., 1996. Multiscale geometrical reconstruction of porous structures. Phys. Rev., E 54, 1734 – 1741. Ferreol, B., Rothman, D.H., 1995. Lattice – Boltzmann simulations of flow through Fontainebleau sandstone. Transp. Porous Media 20, 3 – 20. Gonzalez, R.C., Wood, R.E., 1992. Digital Image Processing. Addison-Wesley Publishing, USA. Hinrichsen, E.L., Aharony, A., Feder, J., Hansen, A., Jossang, T., 1993. A fast algorithm for estimating large-scale permeability of correlated anisotropic media. Transp. Porous Media 12, 55 – 72. Ioannidis, M.A., Chatzis, I., 1993. Network modeling of pore structure and transport properties of porous media. Chem. Eng. Sci. 48, 951 – 972. Ioannidis, M.A., Kwiecien, M.J., Chatzis, I., 1996. Statistical analysis of the porous microstructure as a method for estimating reservoir permeability. J. Pet. Sci. Eng. 16, 251 – 261. Ioannidis, M.A., Kwiecien, M.J., Chatzis, I., 1997. Electrical conductivity and percolation aspects of statistically homogeneous porous media. Transp. Porous Media 29, 61 – 83.

28

J.-F. Daı¨an et al. / Journal of Petroleum Science and Engineering 42 (2004) 15–28

King, P.R., 1989. The use of renormalization for calculating effective permeability. Transp. Porous Media 4, 37 – 58. Laurent, J.P., Frendo-Rosso, C., 1992. Application of image analysis to the estimation, of AAC thermal conductivity. In: Wittmann, F.H. (Ed.), Advances in Autoclaved Aerated Concrete. Balkema, Rotterdam, pp. 65 – 70. Liang, Z., Ioannidis, M.A., Chatzis, I., 2000a. Geometric and topological analysis of three dimensional porous media: pore space partitioning based on morphological skeletonization. J. Colloid Interface Sci. 221, 13 – 24. Liang, Z., Ioannidis, M.A., Chatzis, I., 2000b. Permeability and electrical conductivity of porous media from 3D replicas of the microstructure. Chem. Eng. Sci. 55, 5247 – 5262. Neimark, A.V., 1989. Multiscale percolation systems. Sov. Phys. JETP 69, 786 – 791. Philippi, P.C., Souza, H.A., 1995. Modeling moisture distribution and isothermal transfer in a heterogeneous porous material. Int. J. Multiph. Flow 21, 667 – 691. Quiblier, J.A., 1984. A new three dimensional modeling technique for studying porous media. J. Colloid Interface Sci. 98, 84 – 102. Sahimi, M., 1994. Applications of Percolation Theory. Taylor & Francis, London. Santos, L.O.E., Philippi, P.C., Damiani, M.C., Fernandes, C.P.,

2002. Using three-dimensional reconstructed microstructures for predicting intrinsic permeability of reservoir rocks based on a Boolean lattice gas method. J. Pet. Sci. Eng. 35 (1 – 2), 109 – 124. Stauffer, D., 1985. Introduction to Percolation Theory. Taylor & Francis, London. Stauffer, D., Aharony, A., 1994. Introduction to Percolation Theory, rev. 2nd ed. Taylor & Francis, London. Talukdar, M.S., Tosaeter, O., Ioannidis, M.A., Howard, J.J., 2002. Stochastic reconstruction of chalk from 2D images. Transp. Porous Media 48, 101 – 123. Thompson, A.H., Katz, A.J., Krohn, C.E., 1987. The microgeometry and transport properties of sedimentary rocks. Adv. Phys. 36, 625 – 694. Thovert, J.F., Salles, J., Adler, P.M., 1993. Computerized characterization of the geometry of real porous media: their discretization, analysis and interpretation. J. Microsc. 170, 65 – 79. Xu, K., Daian, J.-F., Quenard, D., 1997a. Multiscale structures to describe porous media: Part I. Theoretical background and invasion by fluids. Transp. Porous Media 26, 51 – 73. Xu, K., Daian, J.-F., Quenard, D., 1997b. Multiscale structures to describe porous media: Part II. Transport properties and application to test materials. Transp. Porous Media 26, 319 – 338.