Reconstructed Porous Media

7

Reconstructed Porous Media

7.1

INTRODUCTION

This chapter describes a recent approach to the analysis of macroscopic transport properties of real, and thus random, porous media. It is a mixture of various techniques that are not difficult to handle individually but may require some time to grasp as a whole. The principle of this method is composed of three major steps. The first involves the measurement of any salient geometric features. Different features can be chosen for various materials. In the case studied here (for the Fontainebleau Sandstone), the porosity and autocorrelation function of the pore space are measured. The second step is the reconstruction process. Random samples of porous media are generated in such a way that, on average, they possess the same statistical properties as the real samples that they are assumed to mimic. Once these samples are generated, in the third step, all transports can be studied at least in principle. The method chosen here consists of the full resolution of the field equations inside the reconstructed sample with adequate boundary conditions. Generally, the macroscopic quantity of interest is obtained by a spatial integration of the local field. For instance, the determination of permeability necessitates the resolution of the Stokes equations of motion and the spatial integration of the velocity field. The method is statistical in character, and samples of porous media are generated at will. The permeability of each sample is determined and averaged over all samples. This method has been applied previously to the determination of permeability (Adler et al., 1990) and the formation factor (Adler et al., in prep.) with reasonable success. Many workers have thought of similar methods, but they were less elaborate than the present one, mostly because automatic analysis of the geometry was less developed than now and because computers were much less powerful than today. The idea of examining thin sections of porous material and measuring various quantities (such as the correlation functions of the pore space) was spurred by the development of the variational methods that used the moments of the phase functions explicitly (see Sections 6.6, 6.7, and 6.8). Corson (1974a,b,c,d) was apparently the first to perform such an analysis on powders; this was done manually. Berryman (1985) and Berryman and Blair (1986) performed similar measurements automatically, thanks to the development of digital image

503

504

Reconstructed Porous Media

processing techniques. Other groups used thin sections to determine characteristic geometric quantities. Lin and Cohen (1982) and Koplik et al. (1984) had in mind the reduction of the medium to a capillary network. Doyen (1988) determined the distribution of throat and pore sizes, and Yanuka et al. (1986) reviewed a number of analogous contributions. Various uses have been made of this information, but it should be emphasized that all these authors avoided the full resolution of the field equations. Corson (1974d) used the Beran bounds to estimate the thermal conductivity of solids; improved bounds are given by Equation 6.177. Recall that they used three-points correlations. Berryman and Blair (1986, 1987) determined the permeability by means of a semi-empirical formula from Walsh and Brace (1984), which is presented in Section 7.4. This formula contains the specific surface area estimated from the slope of the two-points correlation function at the origin. Koplik et al. (1984) and Doyen (1988) used effective medium theory (see Section 6.3) to determine the permeability and the conductivity of the equivalent network obtained from image analysis. This chapter is organized as follows. Section 7.2 covers the reconstruction of real media according to the process devised by Joshi (1974) and Quiblier (1984). It shows how one can generate a three-dimensional random porous medium with a given porosity and a given correlation function. The medium is made of elementary cubes that are filled with a solid or liquid. It can be generated in two steps starting from Gaussian and independent variables X(r). Linear combinations of these variables yield a population Y(r) that is still Gaussian but is correlated. The correlation depends upon the set of coefficients â of the linear combinations. An alternate way that uses Fourier transforms is proposed. This population Y(r) is then thresholded, and it becomes a discrete population Z(r) that takes only two values 0 and 1. This transformation can be viewed as a nonlinear function or filter. The correlation of the population Y(r) is modified by this additional filtering. The average value of Z(r) is automatically equal to ε. Section 7.3 addresses the practical simulation of given porous media, with particular application of the Fontainebleau Sandstone. When measurements are taken, one has to solve an inverse problem; this is performed in two steps. Since the experimental correlation is known, one has to determine the correlation function of the population Y(r). The possibility of this determination is illustrated and discussed. Then the coefficients â can be calculated. Since this book is restricted to isotropic media, the determination of these coefficients can be shortened greatly. Some sets of coefficients are given along with the relative precision. Once these coefficients are known for a given sample, artificial porous media can be generated at will. First, the statistical properties of the simulated media are verified to be close to those of the real media they are assumed to match. The influence of various "artificial" parameters is examined, such as the step in which the experimental correlation is sampled. Simulated cross-sections of the media are compared to actual thin sections; the visual aspect of the cross-sections is quite satisfactory. Section 7.4 presents the calculation of the flow field and the permeability of these reconstructed porous media. First, a brief review of the empirical or semi-empirical formulas that are in the literature is given. It should be emphasized that this review is far

Reconstruction of Porous Media

505

from exhaustive and that experimental results are not covered. Next, the permeability of the reconstructed media is determined. The fluid problem, the method of solution, and the possibilities of the numerical program are presented briefly. Then the influence of the three artificial geometric parameters on the permeability is carefully analyzed. A compromise between various constraints is proposed. For a certain choice of values, these parameters do not seem to bias the results too much. The large amplitudes of the statistical fluctuations are then discussed which are encountered during the Monte Carlo simulation of samples. Obtaining precise statistical averages is hindered by the length of the computations. Finally, these results are compared to classic experimental results obtained by Jacquin (1964). The same approach is used for the determination of the formation factor of the Fontainebleau Sandstone.

7.2 R E C O N S T R U C T I O N O F POROUS MEDIA After a general presentation of the objective of reconstruction, Subsections 7.2.2 and 7.2.3 detail the two basic operations of linear and nonlinear filterings.

7.2.1

General Concepts

For completeness, let's recall some definitions (see Section 2.3). The pore space of a porous material can be characterized by the phase function Z(r) as follows: ^

f 1 if r belongs to the pore space [0

^

otherwise,

where r denotes the position with respect to an arbitrary origin. The porosity ε and the correlation function /? z(u) can be defined by the statistical averages (which are denoted by overbars): ε = Z(r), RZW

2

= [ Ζ ( Γ ) - ε ] [ Ζ ( Γ + ϋ ) - ε ] / (ε - ε ) ,

(7.2a) (7.2b)

where ε a positive quantity limited to the [0, 1] interval. A stronger general statement can be made on /? z(u). It can be shown that a function flz(u) is a correlation function if and only if all its Fourier components are nonnegative. A more rigorous version of this property can be found in Adler (1981) (see Equation 2.54). It is also interesting to note that the correlation function of the solid phase is the same as that of the liquid phase (Equation 7.2b). When the material is assumed to be homogeneous, the statistical averages can be replaced by volume averages. When it is assumed to be isotropic, these volume averages

506

Reconstructed Porous Media

can be replaced by surface averages. Hence, the use of thin sections is justified. The surface integrations that correspond to Equation 7.2 can be performed easily by any image analysis software. Section 7.3 illustrates this. The purpose of the present section is to generate a three-dimensional random porous medium with a given porosity ε and a given correlation function Rz(u). The medium is assumed to be translationally invariant at the large scale. Moreover, for simplicity, the medium is assumed to be isotropic, but it will be seen that this requirement is not essential. Equivalently, the idea is to generate a random function of space Z(r) that is equal to zero in the solid phase and to one in the liquid phase. Z(r) must verify the two average properties in Equation 7.2. ε is a given positive number smaller than 1. Rz(u) is a given function of u that verifies the general properties of a correlation function, but is otherwise arbitrary. 3 medium is constructed in a discrete manner. It is For practical purposes, the porous considered to be composed of N small cubes, each of the same size a. These elementary cubes are filled with either liquid or a solid. Examples of such porous media were provided in Figure 6.27. Hence, the spatial variables r and u take only discrete values. The corresponding trios of integers are denoted by r' = r / a = u > = u /, a =

73

(ij,k), <

t\ (r9s,t).

< · >

Note that the correlation function of isotropic media depends only on the norm u of the vector u (see Adler, 1981), in which case it is denoted by Rz(u). An additional condition is imposed by the fact that the sample of generated porous medium has a finite size aNc. It is then standard to consider periodic boundary conditions on the sample for the determination of permeability (see Adler, 1989). The same requirement should be imposed on the generation of the medium itself. Hence, the random field Z(r) must verify the following: ε = Z(r),

(7.4a)

2 Rz(u) = [ Ζ ( Γ ) - ε ] [ Ζ ( Γ , ) - ε ] / ( ε - ε ) ,

(7.4b)

where the translated vector r, is defined modulo aNc for each of its components as follows: r, = r + u

(mod aNc).

(7.5a)

(modAg.

(7.5b)

For example, this equality means that jt = j + s

Because of this spatial periodicity, all physical quantities are independent of the choice of the origin and the external faces of the cells. Several methods exist to generate discrete random variables that verify Equation 7.2.

Reconstruction of Porous Media

507

Here isotropic media are obtained by a simplified version of an algorithm presented by Quiblier (1984) for general three-dimensional porous media. This algorithm was itself an extension of a two-dimensional scheme devised by Joshi (1974). For clarity, this section presents the algorithm briefly and recalls the major properties of the corresponding random functions. It can be shown that a random and discrete field Z(r) can be devised from a Gaussian field X(r) when the latter is successively passed through a linear then a nonlinear filter. Let's summarize the influence of these filters and relate their properties to the statistical properties of the resulting fields.

7.2.2

Linear Filter

First, some general properties of Gaussian variables are given. A real-valued random variable X is Gaussian (or normally distributed) when its probability density fx(x) is given by

2 fx(x)

=

1

Ι(χ-μ)

ζ exp

(7.6)

2πσ

2 The mean μ and the variance σ are finite and expressed by μ=Ε{Χ}9

2

σ

=

(7.7a)

2 Ε{\Χ-μ\ },

(7.7b)

where Ε denotes the expectation or the statistical average. It can be expressed by means of the probability density fx(x) for any function T(X): E{T(X)} =

(7.8)

JT(x)fx(x)dx.

A Gaussian variable is standard when

2 μ=0

(7.9)

σ = 1.

It is convenient to introduce the characteristic function φχ(ή of the random variable X:

11

Φχ(0 = Ε{β

(7.10)

which is simply the Fourier transform of the probability density. It is also convenient to use the following integral (see, e.g., Gradshteyn and Ryshik, 1965): f

/ 2 2 . χ» V^T 74 exp (-ρ χ ± qx)dx = exp \P

^

, j

p>0.

(7.11)

508

Reconstructed Porous Media

It is then easy to show that the characteristic function of a Gaussian variable can be expressed as

22 φχ(ί)

tG^ = exp it μ - -

(7.12)

Characteristic functions are useful when one must handle sums of random variables. If Xl9 . . ., Xn are independent random variables with characteristic functions ίφχ^ί), / = ! , . . . , # } , then the characteristic function φ8ηof the sum

s = x + . . . + x„ n

x

(7.13a)

is equal to the product of the characteristic functions:

(7.13b) This property can be proved easily for two variables (see Doob, 1953). As a direct consequence, the sum of two independent zero mean Gaussian variables Xx and X2 is a zero mean Gaussian variable whose variance is the sum a\ + a\ of the two variances. More generally, the linear combination Ln = axXx+

. . . + *n Xn

of independent zero mean Gaussian variables Xl9 ..., able whose variance is

(7.14a)

Xn is a zero mean Gaussian vari-

(7.14b) i=l

Now let's consider two Gaussian variables {Χλ, X 2} . Because of the previous properties, this couple is a random variable such as any linear combination is a Gaussian. Thus, this couple is a particular case of the multivariate Gaussian variables (Adler, 1981), which can be defined for reference purposes. An Revalued random variable X = {Xl9 . .., Xd] is called multivariate Gaussian if any nontrivial linear combination of its components is Gaussian. The corresponding probability density can be written as

fx(x)

dl2 = {2n)- \\\

-1/2

exp - - ( χ - Ζ Ο - ν - ^ χ - μ /

(7.15a)

where μ is the mean μ=Ε{Χ]

(7.15b)

Reconstruction of Porous Media

509

and V is the d x d covariance matrix \ =

Ε{(Χ-μ)(Χ-μ)}.

(7.15c)

In the specific case of two independent variables (Xl9 X2), the matrix V is diagonal. The general characteristic function φχ is simply as follows (see Equation 7.12): (7.16) Let's now turn to the generation of a random but correlated porous medium of size N*. First, consider the random field X(i, j , k). The random variables X(iJ, k) are assumed to be normally distributed with a mean equal to zero and a variance equal to one. These 3 be defined 3 variables are independent. A linear operator can by an array of coefficients â(u'), where u' belongs to a finite cube (0, L c) in Z . Outside this cube, â(u') is equal to zero. A new random field Y(r) can be expressed as a linear combination of the random variable X(r'): (7.17) u'e[0,^f

where the translated vector r't is defined modulo Nc for each of its components. An overbar is used for the coefficients â to avoid ambiguity with the size a of the elementary cubes. This definition is identical to the definition used by Joshi (1974) and Quiblier (1984), except for the periodic character introduced by the condition modulo Nc. This additional condition is imposed by the transport problems that are solved in Sections 7.4 and 7.5. Periodic boundary conditions are imposed on the sample N*. To better understand the necessity of this requirement, let's describe the result without the condition modulo JVC . Suppose that Nc is larger than L c, as is generally the case. Then, the upstream and downstream faces of the sample are statistically independent. This implies that the calculations on the fields that are done implicitly on an infinite periodic medium whose unit cell is the sample could be completely biased because of this statistical independence. The pores in the upstream and downstream faces are statistically independent and have no reason to be connected. Without any further requirements for the coefficients â(u') of the linear filter, it can be shown that the random variables Y(r) are Gaussian and centered if Nc is larger that L c. This is an immediate consequence of the properties in Equation 7.14, which were just mentioned. Assume further that the variance of Y(r) is equal to one:

2 E{Y (r)}

= h

(7.18)

Hence, the random variables Y(r) have a standard normal distribution although they are no longer statistically independent. Their correlation function RY(u) is easily seen to be

510

Reconstructed Porous Media RY(u)=

£ r,s,tE[0,L ]

c

â r^ , â u ,+ w

(7.19)

where u + r is determined modulo Nc. An alternative way exists to generate a field of correlated variables that uses Fourier transforms. This method is more convenient than the previous one. Let us recall some formulas from Section 4.3. These are specialized here to discrete fields F(I, /, K), which are assumed to be spatially periodic. F can be expanded in the triple Fourier series as follows: F(r') = ^Fm

exp - 2mkm · r',

(7.20a)

where m is an integer vector belonging to the cube

3 # = [ 0 , W c- l ] .

(7.20b)

The reciprocal vector k mis defined as

k

m — K ^ l + 2^2

m =

e (7.21)

+ »*3 3)>

where ( e l5 e 2, e 3) are the three basic unit vectors. The Fourier components Fm can be expressed as

=

2 π

^ 3 Σ

^

P ~

Γ C X

*™ ' '

2 )2 *

This can be proved easily by multiplying Equation 7.20a by exp 2nikm · r', performing the necessary summation, and using the following relationship:

3 Σ

exp 2 7 r / k m· r' = N · <5 k .m

(7.23)

Suppose that one wishes to generate a random and stationary population F(r') with a given spatial correlation RY(u'). The corresponding Fourier components are denoted by

*r,m = NΛ- Σ ^ c

r'etf

( R , )EXP 2 N I K

™'' R

(- ) ? 24A

511

Reconstruction of Porous Media

One should be careful when performing this summation and should remember that a correlation function is an even function of r'. It is usually more convenient to express this property by using the spatially periodic character of the fields: RY(I\J,K)

= RY(Nc-IJ,Kl

Ie[l,Nc-ll

(7.24b)

and other equivalent properties. Let's generate a set of independent, standard, normally distributed values X(r') (r' G €). This field is identical to the field generated previously. The corresponding Fourier transforms are denoted by

*™ 773 ' Σ ' =

X ( R }E X P 2 N I K M

'' R

(- > 7 24C

Now let's consider a population F(r') of random variables defined by its Fourier com: ponents Ym

or equivalently,

1 / 2 F(r') = a £

( % m)

X m exp - 2mkm · r',

(7.25b)

metf

where α is a normalization factor added for convenience. Note that because of the linear character of this transformation, the field F(r') is Gaussian and centered. The correlation function R'Y(u) can be derived from the previous definitions. Since they are centered, it is not necessary to consider the mean. Hence, RÇ(u) = E{Y(r')Y*(r'

+ u)},

(7.26)

where the asterisk denotes the complex conjugate. R'Y(u) can be evaluated easily by using the previous definitions. One obtains the following:

RYW

(7.27)

=

Since RY(u) is an even and periodic function of u, it follows that

= Κγ, m from which it is readily obtained that

Ργ, m '

512

Reconstructed Porous Media * f ( u ) = ^ / ? y( u ) .

(7.28)

This equation determines the normalization coefficient a since the population Y(r') should have a variance equal to one:

12 (7.29)

a = Nl .

This method is equivalent in principle to the previous one, but it is superior numerically since no inversion is required.

7.2.3

Nonlinear Filter

The random field Y(r) is correlated, but is still not satisfactory since it takes its values in IR, while the porous medium must be represented by a discrete value field Z(r) (see Equation 7.1). To extract such a field from 7(r), one applies a nonlinear filter G, that is, the random variable Ζ is a deterministic function of Y: (7.30)

Z = G(Y).

When G is known, the statistical properties of the random field Ζ can be derived from the properties of Y. For completeness, this derivation by Joshi (1974) is repeated here briefly. Generally, one must determine the correlation function Rz(u) when G and RY(u) are known. Rz(u) can be expressed as follows:

£ ( Z M { 2 (} Zr + { } £ { Z l -f > " 2, >i « 2( u ) = t

£{(Z{r}-£{Z}) }

(7.31)

where it is implicitly assumed that the field Z(r) is stationary. Since the random variable Y(r) has a standard normal distribution (with a zero mean and a variance equal to one), its distribution function P(y) is given by the following (see Equation 7.6):

y P(y) = -^==\

y2/2 e- dy.

V27T J-oo

(7.32)

The deterministic function G is defined by the following condition. When the random variable Y is equal to y, Ζ takes the value z: z=l

whenP(y)<£,

ζ=0

otherwise.

(7.33)

Reconstruction of Porous Media 5 1 3

2 Thus, it is evident that the average value o f Z(r) is equal to ε, and that its variance is equal to ε - ε . The most difficult point is the determination o f the correlation function Rz(u) of Z(r) as a function of RY (u). One can start from the fact that the random variable (T(r), Y(r + u)) is a bivariate Gaussian whose probability density is known as a specific case o f Equation 7.15. Let's state this probability density as follows:

p ( y h y i )

=

2π

yl +

exp

:

yl-^Yyiyi

(7.34)

where y , and y 2correspond to Y(r) and Y(r + u ) . Moreover, 7 = *y(lul).

(7.35)

W h e n all this information is introduced into Equation 7 . 3 1 , it yields

l *z(u) = — — f ε —ε J

Ç°°G(yùG(y2 )p(yi,y2)dyldy2.

(7.36)

It should be noted that these calculations do not require any hypothesis of the isotropic character or the porous medium. The standard w a y to perform the double integration o f Equation 7 . 3 6 involves the expansion o f the Gaussian density in terms o f Hermite polynomials. These cumbersome calculations are skipped here; they use classic identities that can be found in Gradshteyn and Ryshik (1965). The results can be summarized as follows. Rz(u) can be expressed as a series in terms

of RY(u):

c R u) RzW =

(7.37a)

^ » y( > m=0

where the coefficients Cmare given by

Cm = - 7 ^ = ;

2l2 \T(y)e-y Hm {y)dy,

(7.37b)

together with

c(y)-

ε-Ι

if P(y)<

ε,

(7.37c)

514

Reconstructed Porous Media

c (y)=

ι ,f

ifP(y)> x

νε(Ι-ε)

(7.37 d)

ε.

Note that the definition of the Hermite polynomials is not standard. The convention used (y) is defined as follows: by Joshi (1974) and Quiblier (1984) has been adopted here. Hm

m Hm (y) = (-l) e>

2

m l2 r1

j^e~y ay

2

12 .

(7.38)

The subsequent use of these polynomials by Joshi (1974) and Quiblier (1984) was found to be consistent with these definitions.

7.3

APPLICATION

This section details a specific application of the general methodology discussed in Section 7.2. Subsection 7.3.1 begins with the measurement of the correlation function on samples of the Fontainebleau Sandstone. Then, Subsection 7.3.2 shows how the relevant information is extracted from the data. Finally, Subsections 7.3.3 and 7.3.4 discuss typical reconstructed media.

7.3.1

Materials

The Fontainebleau Sandstone was selected here because it is known to have remarkable properties for a geological material (Jacquin, 1964). It is a pure sandstone composed of a single mineral, quartz. It does not contain any clay, hence, no physicochemical interaction occurs between the solid phase and aqueous fluids. Permeabilities measured with air and water are usually identical. The geometric structure of this sandstone is simple since it displays only intergranular porosity. No microporosity exists inside the particles themselves. Moreover, it is known to be remarkably homogeneous. Another important advantage is that the porosity can be varied, while the structure is conserved globally. For all these reasons, many measurements have been made on this sandstone. Figure 7.1 shows several thin sections of this sandstone. The pore space was preserved by injecting a dyed glue into the medium; then the sample was cut. The pore space is thus replaced by the dye-hardened glue (Zinszner and Meynot, 1982). Thin sections about 20 μπι thick were obtained by abrasion of the samples. The pore space is clearly defined since the glue is dyed red. Color photographs of these sections were taken with the help of a microscope. The examples displayed in Figure 7.1 show that this material was originally nonconsolidated sand. Subsequently, various geological processes occurred, such as accretion of the sand particles, with the effect that porosity decreased. Usually in such sandstones, porosity ranges from about 0.03 to 0.35. Let's now turn to the details of the image analysis itself. The first serious difficulty is the transformation of the initial picture (such as Figures 7.1a, b, or d) into an image that

Application

515

Figure 7.1 Thin sections of the Fontainebleau Sandstone. The pore space appears black, (a), (b), and (d) are photographs of thin sections, while in (c), the contrast between the solid and "liquid" phases has been enhanced by various treatments, mostly manual. The bar scale in each photograph is 0.5 mm. The names and measured porosity of each photograph are as follows: (a) 2A3, ε = 0.31; (b) GF2, ε = 0.25; (c) CJ, ε = 0.21; (d) 12A13, ε = 0.11.

516

Reconstructed Porous Media

can be binarized readily, such as Figure 7.1c. The alteration of the initial image during this step must be minimized. This point is discussed thoroughly by Adler et al. (1990). The image is then registered and binarized. The length scale is also measured carefully. A length corresponding to one pixel is denoted as a (in μηι). The binarized image S is divided into two halves, SYand S2. Hence, S = SxuS29

SlnS2

=
(7.39)

The porosity is defined simply as the proportion of pore space contained in a given surface. To get an idea of the homogeneity of the sample, this ratio was measured twice, on S and on Sv The corresponding values are denoted by ε and εχ and are given in Figure 7.1 for some of the samples studied. In these samples, the differences are quite small, but this is not always the case. The constant character of porosity was one of the criteria used in the selection of a thin section. To calculate /? z(w), Sx is first translated by a distance u along the jc-axis; It yields Sx(+u). The spatial average indicated in Equation 7.2b is replaced by an intersection of images: Z(JC,y)Z(x + u,y) =

n S.

(7.40)

The other operations indicated in Equation 7.2b are then performed algebraically. Figure 7.2 shows the results of the images shown in Figure 7.1. Several points of practical interest, such as the size of the images, is discussed by Adler et al. (1990). Generally, a compromise must be found between the statistical fluctuations that increase when the image is small and the variability of the average properties of the porous media at the large scale. The four correlations displayed in Figure 7.2 do not show any dramatic change from one sample to another, although porosities are multiplied by a factor of three. Two of them are mostly superimposed, while the correlation for image 12A13 is intermediate between these two and the fourth one. Other representations of these data were tried. For instance, a correlation length can be introduced by the following formula:

uU

Rz(u) = e~ .

(7.41)

Such an expression is well verified for short distances u only for image 12A13, that is, for the lowest porosity. In other cases, the correlation decreases too rapidly. Such an expression is not expected to apply for large u. No further systematic work was done on this point. The concept of a correlation length <£ is used again in Section 7.4 in the loose sense of a length above which the correlation is negligible. Another check of the quality of the data was done as follows. Since the Fourier coefficients of a correlation function should be positive (see Adler, 1981), they were calculated systematically for the experimental correlations. Then, the negative Fourier coefficients were examined and generally found to be small with respect to the largest positive ones. Finally, the correlation functions were reconstructed with the negative coefficients set to zero. The difference between the reconstructed and the original func-

Application

517

1 4

+V ~X

¥ ~

Α - * ¥

30

20

10

J L+- — — >+ 4

Φ~ ~

υ

40

Figure 7.2 Experimental correlation functions # z as functions of the translation u. The length scale a is always equal to a = 3.8 Mm/pixel. Data are for image 2A3 (+), image GF2 (x), image CJ (V), and image 12A13 (o).

tions was negligible since only the third digit was modified, which is much smaller than the other experimental errors. General practical recommendations on the measurements of these correlation functions can be found in Adler et al. (1990).

7.3.2

The Inverse Problem

When one wants to simulate a given porous medium, the first problem is the determination of the correlation function RY(u) and the set of coefficients ά. This is called the inverse problem here. Once these coefficients are calculated, porous media can be simulated and their general properties can be examined critically. Let's first look at the properties of the nonlinear filter given by Equation 7.37. It follows that

/? (/? ,£) = /? (/? ,l-£). z

y

z

r

(7.42)

The limiting case ε = 0.5 is interesting since it can be calculated analytically. By standard manipulations, it can be shown that Rz = 2/ ;rarcsin(/? y), since the coefficients

are given by

ε = 1/2

(7.43)

518

Reconstructed Porous Media 21 π 0

(7.44)

m even,

(21 nm)(m-2)\\l(m-\)\\

m>3,m

odd.

Figure 7.3 represents the function RZ(RY, έ). It was determined numerically as follows. The integral in Equation 7.37b was evaluated for y ranging from - 1 0 to +10 as in Quiblier (1984). Larger integration domains were tried, but without any influence on the final result. The series in Equation 7.37a is limited to a maximum value m = M. Quiblier (1984) chose M to be 12; however, when RY is close to 1, the series converges very slowly. Our data were obtained with M = 50; above this value, numerical difficulties such as overflows occurred. To get an idea of the numerical precision, let's say that only the first three digits were stable with M = 30 for RY < 0.9. Moreover, for various values of ε, RZ(RY = 1) is equal to about 0.9 instead of 1, which illustrates the deterioration of the precision when RY tends toward 1. This last difficulty was solved in the following way. It was assumed that close to 1, Rz can be expanded as

ll

Rz = l-a{(\-RYY

(7.45)

+

This expansion was suggested by the analysis in Equation 7.43. The constant ax depends upon the porosity ε. It was determined by comparison of Equation 7.45 with the numer-

1

05

0

-05 J

-0.5

0

OS

Figure 7.3 The correlation function Rz as a function of the correlation RY(Equation 7.37). Values of porosity ε are (a) 0.1, (b) 0.2, (c) 0.3, (d) 0.4, and (e) 0.5.

Application

519

ical data at RY - 0.9, where the precision is still acceptable. The validity of this expansion was checked by comparing its prediction with the numerical data at RY = 0.8. For instance, for ε = 0.1, Equation 7.45 gives RZ{RY = 0.8) = 0.511 instead of 0.513 obtained numerically. The fit is thus excellent. The inverse problem relative to the nonlinear filter involves the determination of RY when /? zis known. Thus, Equation 7.42 must be solved. The graphical discussion of this equation is evident from Figure 7.3. Remember that the absolute value of a correlation function is always smaller or equal to 1. When Rz is larger or equal to RZ(RY = - 1 ) , Equation 7.42 also has a unique solution. Otherwise, this equation does not have any solution. This imposes an extra condition on the experimental function Rz. RY was determined as a function of the translation vector u for the experimental data shown in Figure 7.2. Figure 7.4 shows the results. The general shape of the correlation functions Rz is conserved. The values that are the most modified are those close to 1. It is gratifying to note that for the data analyzed here, Equation 7.42 always has a solution. A last check was done on RY. Since they are assumed to be correlation functions, they should also have positive Fourier coefficients. As done for the experimental functions Rz, these coefficients were determined systematically for the correlation functions RY. The same conclusions hold as before. Whenever they are negative, the Fourier coefficients are small. If the correlation functions RY are reconstructed with the negative coefficients set to zero, the difference between them and the original functions is negligible. This concludes the discussion of the inverse problem associated with the nonlinear filter. Now let's turn to the linear filter in Equation 7.17. The inverse problem associated with it involves the determination of the coefficients â when RY(u) is known (see Equation 7.19). A few general comments were made by Adler et al. (1990) on Equation 7.19.

14,

0.5*

^x

10

20

30

40

Figure 7.4 Experimental correlation functions RY as functions of the translation u. Conditions are the same as in Figure 7.2.

520

Reconstructed Porous Media

These are not repeated here since the alternate possibility exists of generating porous media by means of the Fourier transform in Equation 7.25. Now let's consider the numerical solution of Equation 7.19. Joshi (1974) and Quiblier (1984) solved3 it directly. This system is simple in its structure, but it leads to a large number (L c+ l ) of unknowns and equations, even for small values of L c. When the 2 2 greatly 2112since the coefficients porous media are isotropic, this number can be reduced a~rst should only be functions of the distance d = (r + s + t ) : (7.46) Hence, for isotropic media, one has to determine the function â(d) in the L c + 1 points a~(d = 0 ) , . . . , â(d = L c) . The points located in between are calculated by a given interpolation scheme. Equations 7.19 are reduced to the L c+ 1 equations as follows: (7.47) u = 0, . . . ,L C. This system was solved by using a combination of two subroutines of the IMSL Library. Whenever necessary, the coefficients a are determined by a subroutine, which interpolates by cubic splines. A single function is defined as the sum of the squares of the L c + 1 Equations 7.47 and another subroutine calculates a minimum of this function by a quasi-Newtonian method. It is verified a posteriori that this minimum corresponds to a solution of Equations 7.47. This combination of subroutines works well; some examples of sets of coefficients â are given in Table 7.1, along with the value of the minimum ^ # o f the function determined by the second subroutine. It is useless to take into account every point of the correlation RY. Usually one point is selected every ri points. For instance, RY(u = 0), RY(u = n'), RY(u = In),..., RY(u = Lcn') are retained. Some additional work was done on these coefficients, but no general order could be found in the variation of these coefficients. For a given set of data, they appear to be erratic functions of distance. In most cases, negative and positive coefficients are obtained. They all display the same sign only in rare instances. When the second subroutine was started from various initial conditions, different solutions were obtained that possessed the same overall precision. Sets α'and α " corresponding to the following conditions were also compared for the same correlation Ry. a~' was obtained for η and L c = 2 / cand a" for In' and lc. However, a~' was not found to be the interpolated value of a" as was expected on intuitive grounds. These precautions are no longer necessary since the Fourier transform of Equation 7.25 does exactly the same thing without the necessity of solving a nonlinear set of equations. The presentation in terms of the linear filter of Equation 7.17 has been retained here, mostly for historical reasons and for its conceptual simplicity. The Fourier transform of Equation 7.25 is far better, especially for large values of L c. For L clarger than 32, it was long, difficult, and sometimes impossible to invert the system of Equations 7.19 numerically.

a Table 7.1

Application

521

Coefficients â(d ) of the Linear Filter Image

Coefficients fl(0) *(D fl(2) e(3) fl(4) â(5) â(6) âÇJ)

m)

2A3 -0.1608 -0.09612 0.04362 -0.02945 -O.02728 0.00449 0.05580 0.07800 0.02928 2 0.15 χ ΙΟ"

GF2 0.1141 0.1163 0.1981 0.09622 -0.03673 -0.03032 -0.01673 -0.03311 -0.3009 3 0.35 χ ΙΟ-

CJ -0.05350 0.00163 0.1547 0.1641 0.04373 -0.00870 -0.00652 -0.02251 -0.01568 3 0.59 χ ΙΟ"

12Α13 0.2130 0.01930 -0.02923 0.03030 0.01642 -0.01928 0.05862 -0.7218 -0.05229 2 0.52 χ ΙΟ"

The images are displayed in Figure 7.1. Data are for L = 8, n' = 4.

c

7.3.3

Generation of Real Porous Media

Once the coefficients of the linear filter are determined, or equivalently, once the Fourier components RYm (see Equation 7.25) are calculated, reconstructed porous media can be generated at will. It should be emphasized that three "artificial" parameters were introduced in these simulations: n\ L c, Nc,

the distance between two experimental data points, the total number of such data, and the size of the generated cube.

Following a general discussion of the average properties of the simulated media, the influence of these parameters is presented briefly. First, let's look at the statistical properties of an example of porous media. Table 7.2 shows that the agreement between the correlation function of the simulated media and the experimental correlation is good. The influence of the total number of simulated con5 is also clear. One needs to generate at least 5 to 10 such configurations figurations (8 x 10 elementary cubes). The agreement is excellent. It is interesting to note that the average difference between the simulated and experimental correlation is on the same order of magnitude as the error obtained for the resolution of Equations 7.47. This has been found to be generally the case. Whenever the resolution of Equations 7.47 is not precise, the average simulated correlation itself is not in good agreement with the experimental correlation it is supposed to fit. For instance, when n'Lc is too small, Rz(n'Lc) cannot be sufficiently small. In this situation, as pointed out by Quiblier (1984), the precision of the solution of Equations 7.47 is poor and the simulated correlations are also poor.

522

Reconstructed Porous Media Table 7.2 RJLu) 1 0.465 0.163 0.006 -0.038 -0.042 -0.045 -0.034 -0.011 ε = 0.25

a

Statistical Properties of Simulated Porous Media δ 0.002 -0.004 0.003 3 - 6 χ 4ΙΟ-ίο4 0.002 -ΙΟ"3 ίο0.002 —

Jf= 10

Jf= 1 1 0.501 0.214 0.064 -0.013 -0.059 -0.062 -0.041 -0.006 ε = 0.27

1 0.443 0.145 0.010 -0.032 -0.034 -0.046 -0.015 -0.001 ε = 0.25

1 0.455 0.165 0.014 -0.032 -0.041 -0.048 -0.040 -0.016 ε = 0.25

. Y= 20 1 0.457 0.168 0.013 -0.033 -0.043 -0.049 -0.038 -0.006 ε = 0.27

JT= 100 1 0.459 0.165 0.005 -0.042 -0.045 -0.047 -0.034 -0.003 ε = 0.25

Data are for image GF2, L = 8, ri = 4, N = 20. The first column recalls the experimental c c correlation, δ is the residue of Equation 7.47. Average correlations and porosities ε are given for ^configurations.

The influence of the parameters n\ Lc and Nc can now be shown. The choice of this set of parameters is partially restricted. The larger Nc is, the better. Sections 7.4 and 7.5 show that Nc is limited by the numerical resolution of the field equations and that L c should be as large as possible with the constraint that it is at least smaller than NJ2. Before studying influence of n\ the statistical variations of the generated media are illustrated in Figure 7.5 for a given set of parameters. The variety of the generated shapes is striking although similarities between the cross-sections in Figures 7.5b, c, and d are apparent. Figures 7.5a, however, looks different Figure 7.6 illustrates the influence of n'. The images presented were chosen deliberately for their average character. The product η x L c was kept constant in Figures 7.6a, b, and c so that the real length of the correlation remains the same from one image to the other. Hence, the change of η is simply equivalent to a change in the unit length scale. An arbitrary unit scale is given in Figure 7.6. For instance, the same real length looks four times larger in Figure 7.6c than in 7.6a. This is consistent with the visual impression given by these images since Figure 7.6c is a plausible magnification of Figure 7.6a. It is also interesting to note that isolated elementary squares are almost absent in Figure 7.6c. Figures 7.6b and d were obtained with η = 4 and 3, respectively, with the other parameters being identical. No difference is apparent between the two. Figure 7.6b can be compared with the real thin section shown in Figure 7.1b. The two figures share a number of features: the randomness, the connection of a few pores to one another, and the existence of several isolated pores. The influence of porosity on the resulting structure can be examined more completely in Figure 7.7 where the images have been simulated with the same set of numerical parameters. In Figure 7.7a (as in Figure 7.1a), the pores are largely connected to one another. As already stated, Figure 7.7b (as Figure 7.6b) contains two kinds of pores; some are connected and some are isolated. The impression given by Figure 7.7c is different since only four large pores appear. This peculiar feature may be caused by the

Application

c.

523

d.

Figure 7.5 Cross-sections of a sample of simulated porous medium. The pores are black; the solid phase is white. The boundaries of the sample are indicated by the dashed lines. This sample has the same characteristics as image 12A13: ε = 0.11, Nc = 27, L c= 8, and η = 4.

c.

'

'

d.

'

1

Figure 7.6 Cross-sections of simulated porous media. Conventions are given in Figure 7.5. Statistical characteristics are for image GF2: ε = 0.25, Nc = 27. Data are for (a) L c= 4, η = 8; (b) L c= 8, η = 4; (c) L c= 16, η = 2; and (d) L c= 8, η = 3.

524

Reconstructed Porous Media

rai

m η

ι Bp

c. Figure 7.7 Cross-sections of simulated porous media of various porosities. Conventions are given in Figure 7.5. The statistical characteristics are as in Figure 7.1, in the same order. Data are for (a), (b), and (d) Nc = 27, L c = 8, and η = 4. Data for (c) are the same except η = 3.

large experimental correlation Rz (see Figure 7.2). Finally, Figure 7.7d is composed mostly of isolated pores as in Figure 7. Id.

7.3.4

Concluding Remarks

It can be safely concluded that the simulation process devised by Joshi (1974) and Quiblier (1984) works quite well. The porosity and correlation function reproduce well. Moreover, the visual features of the simulated cross-sections are not at odds with the features of the real thin sections. The present work was restricted deliberately to isotropic media. On the numerical side, this restriction only modifies the determination of the coefficients â of the linear filter; the rest of the analysis remains unchanged. It should be emphasized that the additional experimental work required for the analysis of anisotropic media is considerably greater than the additional numerical effort required. Section 7.4 shows that the calculation of the flow field imposes extra conditions on the so-called artificial parameters Nc, L c, and n'. More work is being done currently on these reconstructed media. Figure 7.8 shows three-dimensional representations of the reconstructed Fontainebleau Sandstone. Salles et al. (in prep.) have studied their topological properties.

The Permeability of Reconstructed Porous Media

525

b.

a.

Figure 7.8 Three-dimensional reconstructed porous media based on thin sections of Fontainebleau Sandstone corresponding to the pore phase. The porosities of (a) and (b) are 0.11 and 0.14, respectively.

7.4 THE PERMEABILITY O F RECONSTRUCTED POROUS MEDIA Subsections 7.4.1, 7.4.2, and 7.4.3 are a brief summary of Adler et al. (1990). For completeness, the basic equations are recalled, along with some numerical information. Some aspects have already been detailed in Section 6.8. The roles of several necessary parameters that do not convey any physical information are discussed. The results are compared with data obtained by Jacquin (1964) on the Fontainebleau Sandstone. This section concludes with a brief survey of some recent theoretical and experimental developments on permeability and conductivity.

7.4.1

General Concepts

The low Reynolds number flow of an incompressible Newtonian fluid is governed by the usual Stokes equations: ν/7 = μ ν 2ν ,

(7.48a)

V · ν = 0,

(7.48b)

where v, p, and μ are the velocity, pressure, and viscosity of the fluid, respectively. In general, ν satisfies the conditions ν=0

on 5 T

(7.49a)

526 Reconstructed Porous Media ν is spatially periodic.

(7.49b)

3 of the wetted solid inside the unit cell. The volume τ 0 of this cell S denotes the surface is equal to (Nca) . As in Sections 7.2 and 7.3, one considers a finite sample of size Nca. This system of equations and the conditions apply locally at each point R of the interstitial fluid. In addition, it is assumed that either the seepage velocity vector ν is specified, that is, ν = — f Rds · ν = a prescribed constant vector, τ

0 Jdio

(7.50a)

or else that the macroscopic pressure gradient Vp is specified, Vp = a prescribed constant vector.

(7.50b)

As is well known, these two quantities are related by the permeability tensor K, such that ν = -(1/μ)Κ·ν/7.

(7.51)

Here Κ is a symmetric tensor that is positive definite. Usually, only the average permeability Κ is of practical interest. Since the media that have been studied here are isotropic, Κ is a spherical tensor: K =K L

(7.52)

The numerical method is the same as that used in Sections 5.10 and 6.8. It is the artificial compressibility method with a staggered marker and cell (MAC). It is integrated by the alternating direction implicit (ADI) scheme. Further numerical details are given by Adler et al. (1990). The mesh spacing Δ is the same in all three directions. It is a fraction l/n of the smallest length that exists in the medium under consideration. Thus, it is given by Δ=α/η,

(7.53)

where η is an integer. Recall from Section 6.8 that the influence of η on the final results is largely negligible when compared to statistical fluctuations.

7.4.2

Conditions on Numerical Parameters

A number of parameters that do not correspond to any physical reality are necessary to introduce in the course of this analysis. They are discussed here briefly in abstract terms. Examples and detailed results are given in Adler et al. (1990). First, let's recall the various parameters that intervene in the permeability calculation.

The Permeability

of Reconstructed Porous Media

527

The size a of the elementary cubes in the unit cell is related to the image analysis in the following way. When the images are binarized, it is easy to determine the length a (in pm) corresponding to one pixel. The experimental correlation is sampled every η pixels. Hence, the size a can be expressed as a = n'a.

(7.54)

2 Κ is derived simply from the dimensionless permeability The dimensional permeability by multiplying it by (na) . The role of the three artificial geometric parameters n\ L c, and Nc can be summarized as follows. A first condition is implied by the existence of the periodic boundary conditions of Equation 7.5 on the generated unit cell. Because the correlation function is even, to prevent any overlap between the negative and positive values of the translation distance u through the spatial periodicity, L cshould be smaller than NJ2. A second important factor is related to the statistical fluctuations. If NJ2 is much larger than L c, one expects these fluctuations to diminish. Then, a unit cell will contain many "independent" samples of the porous medium. Hence, the first condition can be summarized by (7.55a)

«NJ2 LC

The two other conditions are related to the sampling of the experimental correlation. Recall from Equation 7.41 that <£denotes the correlation length. The correlation should be very close to zero for the maximum sampling distance Lca. In terms of n', this necessitates the following requirement: —

an'

«

L c.

(7.55b)

The second sampling condition is more subtle in its effects on the resulting media. It is natural to require that the correlation length contains many points, that is, 1 « — .

an'

(7.55c)

When this is not the case, the first value Rz(u = a) of the correlation is small, and two adjacent elementary cubes are almost uncorrelated. This corresponds to the classic site percolation that was studied extensively in Chapter 6. When the porosity is smaller than a critical value ec - 0 . 3 1 2 (see Table 6.1), infinite cubic arrays do not percolate, that is, the corresponding permeability is equal to zero. This is exactly what happens here: when the step n' is too large, the correlation between two adjacent cubes is too small. This point has been verified numerically by Adler et al. (1990). All previous conditions can be combined in the following series of inequalities: 1«

«

an

Lc « Nc / 2.

(7.56)

528

Reconstructed Porous Media

The stronger these inequalities are, the better. However, note that the size Nc is limited to 27 on a Cray XMP computer (with a central memory of 3 megawords) and to 100 on a Cray 2 (with a central memory of 256 megawords). The length of the computations should always be kept in mind. The increase in Nc necessitates a large increase in computer time. The reader can find a detailed study of these parameters in Adler et al. (1990), along with numerical illustrations.

7.4.3

Comparison with Experimental Results

The first application of the method of reconstructed porous media was done with the measurements of Jacquin (1964) on the Fontainebleau Sandstone using the samples shown in Figure 7.1. The porosity ε and the permeability Κ were measured on a large number of cylindrical samples of diameter 2.5 cm and length 3 to 4 cm. The global porosity was measured by comparison of three different weights: the weight of the dry sample, the weight of the sample saturated with water, and the apparent weight of the saturated sample immersed in water. Permeability was measured with air and water on a few samples. Except for small permeabilities in which the Klinkenberg effect takes place, no significant difference was observed between air and water. Hence, all other measurements were performed in air only. Figure 7.9 shows permeability data. It is tempting to correlate them by a power law: Koce\

(7.57)

where η is equal to 4.15. This law is in good agreement with the data for intermediate porosities, but it tends to overestimate Κ for low and high porosities. It is interesting to note the existence of large fluctuations in the data: for a given value of the porosity, the permeability can differ by nearly a factor of ten. This is not contradictory with the homogeneous character of the Fontainebleau Sandstone. It was shown by Henriette et al. (1989) that the permeabilities of limestones can vary by two orders of magnitude. The numerical data are also reported in Figure 7.9 and compared to the experimental data. At low porosities, some samples do not percolate, and the corresponding vertical arrows in the graph have no lower ends since the permeability is equal to zero. These nonpercolating samples were taken into account in the determination of the average permeability. It is important to remember that the comparison of Figure 7.9 does not involve any hidden adjusted parameter and that every quantity is measured or calculated. The calculated permeability differs by, at most, a factor of five from the measured permeability. However, the general shape of the experimental curve is predicted in an accurate way as if a systematic error was made in the measurement of the unit scale. It was also shown by Adler et al. (1990) that the Carman equation (see Equation 6.199) is not well verified. Although it is not totally satisfactory, this comparison is still the best to our knowledge. Several tentative ideas have been given by Adler et al. (1990) to explain the discrepancy between the data and the calculations. Three factors may cause it. The first and

The Permeability of Reconstructed Porous Media

529

Κ (mdy)

10

0

0.1

ε 0.2

0.3

Figure 7.9 Permeability (for air in milli Darcy) as a function of the porosity ε. The dots are experimental data. The average numerical permeability Κ was calculated for the four previous samples plus an additional one; they are indicated by crosses. Data are for Nc = 27, L c = 8, and η = 4. The vertical arrows indicate the interval of variation of the individual permeabilities.

simplest cause may be that the samples are too small to make the experimental correlation with small n' and to minimize the statistical fluctuations. This is deferred to future versions of the program. Second, a basic difference exists between the numerical model and the real experiment. In the former, a single elementary cube can block the entire flow without disturbing the average statistical properties of the sample. This situation is never exactly met in experiments, where the fluid can filter through even very small constrictions. Moreover, permeability is measured on large samples (that would correspond roughly to Nc = 170, n' = 4), where such blocking effects are minimized. Such large values of Nc are not totally out of reach in the present development of computer technology. Finite-size effects could also be minimized with large Nc. The third reason is more subtle and may be caused by the simulation process itself. Only two statistical properties of the real media are fitted in the present model. This may not be sufficient to approximate their geometry with acceptable precision for a transport quantity such as permeability, which is very sensitive to connectivity.

7.4.4

Some Recent Developments

Although not main stream to the chapter, it is appropriate here to discuss some experimental work on permeability and its interpretation through adequate correlations. How-

530

Reconstructed Porous Media

ever, this topic is not covered completely here. Rather, a few general references are listed, and some recent contributions are discussed. Some general references have already been cited in Section 6.7, where the Carman-Kozeny equation was given (see Equation 6.199), as well as some specific correlations on beds of spherical particles (see Equation 6.202) and fibers (see Equation 6.203). A recent review on various aspects of sedimentary rocks was done by Thompson et al. (1987). Recent contributions can be classified into two broad categories. In the first, which includes most of the papers, capillary networks are utilized, whereas this underlying model is not used in the second group. It should also be emphasized that many authors address the problem of the formation factor of the permeability simultaneously. When this is done, the contribution are not dissociated into two parts. Walsh and Brace (1984) improved a model originally due to Wyllie and Rose (1950). The pore space is viewed as a series of plane or circular capillaries that are distributed isotropically. By means of elementary arguments that use the average velocity inside a capillary (see Equations 3.78) and 3.79), Walsh and Brace (1984) showed that the permeability Κ and the formation factor F were given by the following:

2 m

ε

(7.58a)

(7.58b)

ε

where m is the hydraulic radius defined as the pore volume divided by the pore area (see Equation 4.281), and THis the tortuosity defined as the square of the ratio of the effective path length to the apparent path length (see Equation 2.13). Also, b is a constant equal to 2 for tubes and 3 for cracks. Berryman and Blair (1987) gave an interesting form to these relationships (see Section 7.1) by eliminating THbetween Κ and F and by introducing the specific area s, defined as the average interface area per unit volume (see Equation 2.9):

2 bs F since it is shown easily that m=-. s

(7.60)

Wong et al. (1984) devised the bond shrinkage model. Consider a capillary network on a ^-dimensional simple cubic lattice in which each capillary corresponds to a cylindrical fluid-filled tube with a radius r,. The model consists of choosing a tube element randomly and reducing its radius by a fixed factor x:

The Permeability of Reconstructed Porous Media

531

(7.61)

η->χη.

The resulting permeability and conductivity of the network can be calculated by application of the standard techniques such as that described in Section 4.13. Analytical calculations can be performed for d = 1, and numerical calculations for higher dimensions. Κ and F are found to follow power laws:

m

K o c

£

F o c

£~

'

m

7.62a) (7.62b)

The coefficients m and m' depend individually upon factor x, but for dimensions higher than two, they are related to one another by m = 2m.

(7.62c)

Hence, the relationship already stated in Equation 6.204 is obtained:

2

K oc F ~ .

(7.62d)

Experiments were performed on artificial rocks made of fused glass beads. Two regimes were observed experimentally for the formation factor: ε>0.2,

m = 3/2,

ε<0.2,

m = 2.3

(7.63)

No explicit permeability data were given by Wong et al. (1984), but Equation 7.62d was well verified for all porosities. Katz and Thompson (1986) developed another model based on an observation by Ambegaokar et al. (1971). Consider a random system with a broad distribution of conductances. The macroscopic conductance is thought to be dominated by the largest conductances. Important simplifications can then be made. One takes into account only the conductances larger than some value g c. This value is chosen in such a way that the set of conductances larger than gc percolates. The conductances larger than gc are replaced by g c, and the smaller ones by zero. Since by definition, conditions are close to a percolation threshold, expansions of the form of Equation 6.56 are proposed for conductivity and permeability. The exponent t is taken to be the same in both cases, that is, t =1.9 (see Equation 6.62). By means of this relationship, the macroscopic transport coefficient is maximized, and this condition determines a characteristic length for each process. The final expressions for conduction and permeability contain unknowns that can be eliminated, and one obtains

532

Reconstructed Porous Media

where / c is a critical length that can be obtained from a mercury injection curve. Let pc be the critical pressure corresponding to the formation of an infinite cluster; / c is expressed by 4cos_^

) _1

where θ is the contact angle. Equation 7.64 was well verified within a factor of 2 or 2 for various sandstones and carbonates (see Katz and Thompson, 1986; Thompson et al., 1987). Johnson et al. (1986) proposed a new characteristic length A defined as follows. First consider a pore space of constant local conductivity <70. The corresponding macroscopic conductivity is σ 0 and the corresponding formation factor is F0. Now consider a disturbed local conductivity σ 0+ <5o(r). It is further assumed that these disturbances affect only the interfacial region and are thus of the form

<5(7(r) = / ( r j ) ,

(7.66)

where η is a local coordinate measured from the solid surface. The corresponding macroscopic conductivity σ can be expanded easily in terms of this disturbance δ σ(τ). One obtains

1

σ = ^ 0- [ σ 0+ 2 Χ 5/ Λ ] ,

(7.67)

where Σ 8and A are expressed by Σs=

(7.68a)

\7(η)άη, Jo

ί |Vyo(r)|''Vr 2

Λ = 2-jr-

[

- j — .

(7.68b)

\V¥o (rfds

ψ0 is the resulting field in the case of a constant local conductivity. All previous studies are global in character, that is, their predictions are comparable with overall measurements. It is only recently that Saleh et al. (1992) undertook systematic local measurements by means of particle image displacement velocimetry or PIDV. This technique has the great advantage over the classic laser doppler anemometry of giving the two-dimensional velocity field in a plane all at once. The information can be recorded on film and then treated after the experiment by automatic image analysis. Measurements can be performed only through transparent media. It was first tried with a medium composed of broken Pyrex, and a mixture of oils was used as the fluid. The composition of the oil is such that its resulting optical index is matched to that of the Pyrex. Figure 7.10 shows an example of streamlines using this technique, along with the corresponding velocity field. At present, this technique is in full development.

The Permeability of Reconstructed Porous Media

533

Figure 7.10 Fields inside an arbitrary slice of porous medium. In (a) and (b), the Pyrex particles are black, (a) Picture for the measurement of the velocity field, (b) Streamlines, (c) The raw velocity field, (d) The filtered velocity field. (Reprinted with permission from Saleh et al. [1992].)

534 Reconstructed Porous Media

7.5 FORMATION FACTOR O F RECONSTRUCTED POROUS MEDIA The method of reconstructed media can be applied to the determination of the macroscopic conductivity a o r , equivalently, its dimensionless inverse, the formation factor F, which is used by practioners in the oil industry. First, media are reconstructed according to the same statistical rules as given previously. Then, the Laplace equation is solved for these media and the results are compared to results obtained for the Fontainebleau Sandstone. This section is based on Adler et al. (in press, 1992), and it concludes with a brief survey of recent developments. First, let's describe briefly the measurements (see Section 7.4 for comments on the use of the Fontainebleau Sandstone). The porosity ε and the macroscopic electrical conductivity σ were measured on a large number of cylindrical samples of diameter 2.5 cm and length 3 to 4 cm. Figure 7.11 shows these results in terms of the formation factor F (Jacquin, 1964):

F = G0/a,

(7.69)

where σ0is the conductivity of the fluid phase alone. This quantity is frequently used in the oil industry. It follows Archie's law (see Section 4.6):

m (7.70)

Foc£~ 9

with a cementation factor of m = 1.64. Since these samples are the same as those used previously, the same thin sections have been used for the reconstruction of numerical media (see Figure 7.1). Let's recall the field equations to be solved in a sample of reconstructed porous medium with periodic boundary conditions at the surface of the unit cell. Electrical and thermal conductions are both governed by a Laplace equation,

2

V 7 = 0,

(7.71)

where Τ is the local field, together with the no-flux boundary condition at the wall Sp when the solid phase is assumed to be insulating: mVr =0

o n S p.

(7.72)

where m is the unit normal vector to 5 p. Vr is assumed to be spatially periodic with a period aNc in all three directions of space. In addition, either the macroscopic temperature gradient Vr or the average heat flux q,

Formation Factor of Reconstructed Porous Media 5 3 5

F 100. 50

10

0J05

0.1

0.2

E

Figure 7.11 The formation factor F as a function of the porosity ε. The dots are experimental data. The average numerical formation factor is indicated by a cross. Numerical data are for η = 2, and Nc = 80. The vertical arrows indicate the interval of variation of the individual formation factors.

is specified. S is the surface of the unit cell. The quantities VT and q are related by the symmetric positive definite conductivity tensor: q = - ( 7 Vf,

(7.74)

which depends only on the geometry of the medium. On average, for an isotropic random medium, σ is a spherical tensor equal to
536

Reconstructed Porous Media

the computational time, a three-dimensional sample with nNc = 1 2 0 requires a memory of about 100 MB and can be handled easily. For the same reasons as in Section 7.4, it is evident that the artificial parameters describing the reconstructed porous media must fulfill the conditions shown in Equation 7.56. Figure 7.11 shows some numerical results. The ratio between the average calculated and measured formation factors is always better than 3. Such an agreement is quite satisfactory when one considers the extreme variability of the transport properties of geological porous media (Henriette et al., 1989) and the fact that no parameter has been fitted in the present comparison. All parameters have been either measured or calculated. At low porosity, some of the samples have a liquid phase that does not percolate, hence, the corresponding formation factors are infinite, as indicated in Figure 7.11. The average formation factor is calculated as the inverse of the average macroscopic conductivities. F{e = 0.21) is slightly larger than F(e = 0.14), simply because the number of calculated configurations is always smaller than 8. Hence, the statistics of the displayed data are poor. This is due to the length of the computations. Note that the numerical results always overestimate the formation factor, while in Section 7.4, the computed permeability always underestimated the conductivity. These results are not contradictory since a permeability is a "conductivity" and a formation factor is the inverse of the conductivity. Hence, in terms of conductivity, both numerical results underestimate the experimental conductivities. It is thus consistent with Section 7.5 that the description and reconstruction of porous media must be improved to obtain better agreement between the calculated and experimental transport coefficients. Overall, it seems that the proposed method of reconstructed porous media is promising.

7.5.1

SOME RECENT DEVELOPMENTS

This section appropriately concludes with a brief survey of the literature of its recent developments. Most relevant references on the formation factor have already been presented throughout this book. Archie's law (Equation 7.70) has been discussed in Section 6.6 (see Equation 6.144). Section 7.1 and Subsection 7.4.4 have given some further information on permeability. Some attention has been given to correlated disorder in networks. Three-dimensional cubic networks were addressed by Webman et al. (1975), as mentioned in Section 6.3. Havlin et al. (1988) studied a one-dimensional network amenable to analytical calculations. The only reference that has not yet been cited is from Roberts and Schwartz (1985). They specifically addressed materials consolidated by diagenesis, that is, materials that were originally packings of particles. They considered ordered or random packings of spheres that are allowed to grow in a uniform manner. The macroscopic conductivity was assumed to be proportional to the minimal cross-sectional area of a typical throat between the fused spheres. This simple model gives surprisingly good results when compared to the experiments of Wong et al. (1984), which were performed on fused glass beads. The comparison is done with a fitted proportionality factor.

References

537

REFERENCES Adler, P.M., 1989. Flow in porous media, in Avnir, D., ed., The fractal approach to heterogeneous Chemistry, John Wiley, New York. Adler, P.M., Jacquin, C.G., and Quiblier, J.A., 1990. Flow in simulated porous media, Int. J. Multiphase Flow, 16,691-712.

Adler, P.M., Jacquin, C.G., and Thovert, J.F., in press, 1992. The formation factor of reconstructed porous media, Water Res. Res. Adler, R.J., 1981. The geometry of random fields, John Wiley, New York. Ambegaokar, V., Halperin, B.I., and Langer, J.S., 1971. Hopping conductivity in disordered systems, Phys. Rev., B4, 2612-2620. Berryman, J.G., 1985. Measurement of spatial correlation functions using image processing techniques, / . Appl. Phys., 57, 2374-2384. Berryman, J.G. and Blair, S.C., 1986. Use of digital image analysis to estimate fluid permeability of porous materials: Application of two-point correlation functions, / . Appl. Phys., 60, 1930-1938. Berryman, J.G. and Blair, S.C., 1987. Kozeny-Carman relations and image processing methods for estimating Darcy's constant, / . Appl. Phys., 62, 2221-2228. Corson, P., 191'4a. Correlation functions for predicting properties of heterogeneous materials, I. Experimental measurement of spatial correlation functions in multiphase solids, / . Appl. Phys., 45, 3159-3164. Corson, P., 1974b. Correlation functions for predicting properties of heterogeneous materials, II. Empirical construction of spatial correlation functions for two-phase solids, / . Appl. Phys., 45,3165-3170. Corson, P., 1974c. Correlation functions for predicting properties of heterogeneous materials, III. Effective elastic moduli of two-phase solids, / . Appl. Phys., 45, 3171-3179. Corson, P., 1974d. Correlation functions for predicting properties of heterogeneous materials, IV. Effective thermal conductivity of two-phase solids, / . Appl. Phys., 45, 3180-3182. Doob, J.L., 1953. Stochastic processes, John Wiley, New York. Doyen, P.M., 1988. Permeability, conductivity and pore geometry of sandstone, / . Geophys. Res., 93,1129-11AO.

Gradshteyn, I.S. and Ryshik, I.M., 1965. Table of integrals series and products, Academic Press, New York. Havlin, S., Blumberg Selinger, R., Schwartz, M., Stanley, H.E., and Bunde Α., 1988. Random multiplicative processes and transport in structures with correlated spatial disorder, Phys. Rev. Lett., 61, 1438-1441. Henriette, Α., Jacquin, C.G., and Adler, P.M., 1989. The effective permeability of heterogeneous porous media, Phys. Chem. Hydrodyn., 11, 63-80. Jacquin, C.G., 1964. Corrélation entre la perméabilité et les caractéristiques géométriques du grès de Fontainebleau, Revue IFP, 19, 921-937. Johnson, D.L., Koplik, J., and Schwartz, L.M., 1986. New pore size parameter characterizing transport in porous media, Phys. Rev. Lett., 57, 2564-2567. Joshi, M.Y., 1974, A class of stochastic models for porous media, Ph.D. Thesis, Univ. of Kansas, Lawrence, KS. Katz, A.J. and Thompson, A.H., 1986. Quantitative prediction of permeability in porous rocks, Phys. Rev., B34, 8179-8181. Koplik, J., Lin, C , and Vermette, M., 1984. Conductivity and permeability from microgeometry, /. Appl. Phys., 56, 3127-3131.

Lin, C. and Cohen, M.H., 1982. Quantitative methods for microgeometric modelling, / . Appl. P/rys., 53, 4152-4165.

538

Reconstructed Porous Media

Quiblier, J.A., 1984. A new three-dimensional modeling technique for studying porous media, /. Coll. Inter/. Sci., 98, 84-102.

Roberts, J.N. and Schwartz, L.M., 1985. Grain consolidation and electrical conductivity in porous media, Phys. Rev., B31, 5990-5997. Saleh, S., Thovert, J.F., and Adler, RM., 1992. Measurement of two-dimensional velocity fields in porous media by particle image displacement velocimetry, Exp. Fluids, 12, 210-212. Salles, J., Thovert, J.F., and Adler, RM., in prep. Automatic characterization of the geometry of porous media. Thompson, A.H., Katz, A.J., and Krohn, C.E., 1987. The microgeometry and transport properties of sedimentary rocks, Adv. Phys., 36, 625-694. Thovert, J.F., Wary, F., and Adler, RM., 1990. Thermal conductivity of random media and regular fractals, / . Appl. Phys., 68, 3872-3883. Walsh, J.B. and Brace, W.F., 1984. The effect of pressure on porosity and the transport properties of rocks, / . Geophys. Res., 89, 9425-9431. Webman, I., Jortner, J., and Cohen, M.H., 1975. Numerical simulation of electrical conductivity in microscopically inhomogeneous materials, Phys. Rev., Β11, 2885-2892. Wong, P.-Z., Koplik, J., and Tomanic, J.P., 1984. Conductivity and permeability of rocks, Phys. Rev., B30, 6606-6614. Wyllie, M.R.J, and Rose, W.D., 1950. Some theoretical considerations related to the quantitative evaluation of the physical characteristics of reservoir rock from electrical log data, Trans. Am. Inst. Mech. Eng., 189, 105-118.

Yanuka, M., Dullien, F.A.L., and Elrick, D.E., 1986. Percolation processes and porous media, I. Geometrical and topological model of porous media using a three-dimensional joint pore size distribution, / . Coll. Interf. Sci., 112, 24-41. Zinszner, B. and Meynot, C , 1982. Visualisation des propriétés capillaires des roches réservoirs, Revue IFP, 39, 337-361.