Irreducible water distribution in sandstone rock: Two phase flow simulations in CT-based pore network

Phys. Chem. Earth (A), Vol. 25, No. 2, pp. 169-174, 2000 0 2000 Elsevier Science Ltd All rights reserved 1464.1895/00/$ - see front matter

Pergamon

PII:

Irreducible Water Distribution CT- based Pore Network

S1464-1895(00)00027-2

in Sandstone Rock: l3vo Phase Flow Simulations

Minyao Zhou, Dalian Lu, John Dunsmuir

and Hans Thomann

Exxon Research and Engineering

Clinton Township,

Company,

Received 23 April 1999; accepted 30 September

Route 22 East, Annandale,

in

New Jersey 0880 1, USA

1999

Abstract. A simple cutoff approach based on the capillary bundle model has become an industrial standard method to quickly obtain the irreducible water saturation from NMR T, distribution. However this approach is not always valid. To overcome the shortcoming of the capillary bundle model that ignores pore-pore connectivity, we have conducted a two-phase flow simulation in a CT-based pore network. The CT -based pore network is a representation of a real rock pore structure that is described by a binary X-ray tomographic data set. Simulation on such network mimics a process of the porous plate measurement. The generated capillary curve is quite reasonable. The oil accession distributions at different water saturations plotted as a function of pore size provide an insight for the immiscible displacement process in the real rock pore structure. Water is trapped not only in de4 ends but also in the wellconnected pores due to a pore level by-pass mechanism. At the capillary end point pressure, a plot of the trapped water distribution as a function of pore size has the same lognormal distribution as the pore size distribution, which is much different from what the simple capillary bundle model suggests. 0 2000 Elsevier Science Ltd. All rights reserved.

for the estimation of Swirr. The physical basis underlying the current industry approach for deriving Swirr from NMR logging data is It is assumed that rock pores can be illustrated in Fig.la. represented by a distribution of isolated pores and that the NMR magnetization signal from fluid in a given pore has a single relaxation time Tz. For most reservoir rocks the NMR relaxation surface relaxation is a dominant mechanism. Therefore, the relaxation time is proportional to the ratio of pore volume to pore surface, V, IS, which is a good measurement of the pore size. The observed NMR signal decay is non-exponential because the signal is a sum of contributions from all pores with a wide pore size distribution. Mathematically, the signal with nonexponential decay can be represented by a distribution (or spectrum) of relaxation times. Due to the proportionality between the relaxation time T, and the pore size, the profile of T, distribution can be approximately understood as a profile of the pore size distribution. In order to equate a capillary pressure used in the immiscible displacement process with the pore size derived from the NMR relaxation time, a capillary bundle model has been employed in the industry. In the model, the capillary pressure can be related to the tube radius by:

1.

p _ 2ocose c--=-’ R

Introduction

One parameter of particular interest during the field appraisal of a newly discovered reservoir is the irreducible water saturation, Swim. The irreducible water saturation is For a defined by the capillary end point saturation. capillary pressure curve in which oil displaces water in a water saturated rock, the capillary end point saturation of the water is defined by the pressure at which the slope of the capillary pressure curve approaches infinity. Standard experiments to determine Swirr include the centrifugation method and the porous plate flow method. In situ, Swirr is estimated conventionally through the resistivity logging measurement that however is not a direct Recently NMR measurement of the pore fluid properties. logging, a direct measurement of pore fluid, has been used Correspondence

ocose

(1)

(V/S)

where d is a surface tension between two fluids, R radius of tube and 8 a contact angle. In this simple model, the nonwetting fluid (oil) is assumed to enter the capillary tube with the largest diameter, followed by successively smaller diameter tubes, as the capillary pressure P,, which is defined as the pressure difference between two fluid phases, increases. The distribution for two phases is simply divided by a vertical line as shown in Fig.lb. The position of the vertical line shifts to the left as the oil saturation increases. Following this intuitive picture, the irreducible water saturation is estimated by a simple cutoff, which is calibrated using laboratory data of rock samples from the same well. This simple capillary bundle model is not physically realistic since it does not take into account the pore

to: Zhou, M. 169

M. Zhou et al.: irreducible Water Distribution in Sandstone Rock: Two Phase Flow Simulations in CT-based Pore Network

170

NMR T2 (msec) Fig.1 (a) Model for analysis of NMR relaxation time data: the relaxation time is proportional to the pore size. (b) Bundle pore filling model gives a vertical line as the phase boundary that shifts to left when oil saturation increases by assuming same surface relaxation rate for two phases and a simple cutoff determines the irreducible water saturation.

connectivity

in the three dimensional pore-networkstructure. The model therefore cannot explain the capillary end point saturation. The bundle pore-filling model leaves no water in the rock if the capillary pressure increases with no limit. To derive the irreducible water saturation from NMR logging data, we need an advanced model that takes into account pore-pore connectivity. In principle, all irreducible water is capillary-bound, but a detailed analysis tells us that the bound water can be classified into four types. These are (1) water trapped on high surface area clays; (2) water trapped in pendular rings due to pore nontube geometry and surface roughness; (3) water trapped in dead ends of the pore-network; and (4) water trapped in pores due to the pore level by-pass. The first three types of water trapping are easily understood. The last phenomenon has been reported in the model studies, such as the invasion percolation model (Lenormand & Bories, 1980; Chandler et. al.,1982). However, we have found no reports in the literature about such investigations in a real rock pore structure. Motivated by the NMR data interpretation, we have recently conducted the two-phase flow simulations to study the by-pass phenomena in clean sandstone and found a common feature in the irreducible water distribution that leads to a new pore-filling model. In the section 2 of this paper we first describe the microtomography image (or CT image) of the sandstone and then present a simple method that maps the CT data to a pore network. In section 3, we discuss the detailed physics of the simulation. We report the results of numerical simulation in the section 4. Some important aspects related to the simulation are discussed in the last section. 2.

Tomographic Node/Channel

Data Analysis Network

to

Define

The sandstones used in this study were 5mm diameter by 1Omm long right cylinders cored from epoxy backfilled samples of Fontainebleau sandstone originally part of a comprehensive study (T. Bourbie, and B. Zinszner, 1985). Micro CT data for these specimens were obtained using our imaging facility at beamline X2-B at the National Synchrotron Light Source, Brookhaven National Laboratory (B.P. Flannery et. al., 1987). The scanner is capable of acquiring a 1024 cube of image data with 14 bits of dynamic range in about 1 hour. Although the scanner is capable of resolution approaching 1 micron, we scanned the specimens at lower magnification to capture a representative sample of the rock. For these rocks, higher resolution images would contain only a few grains and would not

provide meaningful results (F. M. Auzerais et. al., 1996). The projection data were reconstructed using direct Fourier inversion. A 500 x 500 x 500 portion of the reconstructed volume was selected for analysis. The smaller cube with a scale of about 7pmlpixel is completely filled by the 3D rock image. We convert the grayscale CT image into a set of binary rock-pore data using an anti-diffusion segmentation algorithm (P. Perona, T. Shiota, J. Malik, 1994). The resulting binary volume represents a real pore space of a 3.5mm cube of sandstone. Calculation of absolute permeability using both the network model described below and Lattice-Boltzmann (N.S. Martys et. al., 1999) give very good agreement with available experimental data. From this we conclude that the microCT for these samples has successfully captured the important aspects of geometry and topology that control fluid flow. In principle, we can do a first principle calculation in real pore space and obtain an accurate solution of two-phase flow but it is hard to implement such calculation in an ordinary workstation. Even the LB simulation, which is faster than the first principle calculation, required a parallel supercomputer. Instead, we conduct the simulation in a pore network that is constructed from the real pore structure. A mapping from CT data to a pore network reduces the data set sizes from tens of millions to tens of thousands of digital words, a size manageable by a workstation. The critical features of the pore structure that determine the flow properties, such as the pore size distribution and the pore connectivity, ine preserved in the mapping process. The mapping algorithm has following steps. (1) The CT data set of L3 elements is divided into L slices with each slice containing a subset data of L’. (2) In every slice there are pore areas disconnected from each other. Such a pore area, referred to as a poral, is labeled with an identifying number. (3) Each poral may connect with porals in the two neighboring slices. A connecting number Z defines the number of these connections. (4) A poral is defined as a node if Z is larger than 2. (5) A pore space between two nodes is defined as a channel. In this mapping procedure, the distribution of coordination numbers Z has a sharp peak at 3. For one sample of Fontainebleau sandstone with porosity 0.128, 2-3 for 84% of the nodes, Z=4 for 12.8%, Z-5 for 2.6% and a few nodes were observed with Z from 6 to 9. The mapping process generates a list of nodes and channels. Each channel is divided into a set of sections. The volume, surface area, and position in the 3D image of each section are also recoded in the mapping procedure.

M. Zhou el al. : Irreducible Water Distribution in Sandstone Rock: Two Phase Flow Simulations in CT-based Pore Network

500 1

SW=052

1

n

171

0)

1 Pore Size (V/S)

Pore Size (V/S)

1

1

Pore Size (V/S)

Pore size (V/S)

poresize(in unit of pixel)is shown by the black filling area. The distribution is plot in terms of channel (pore) number at the water saturation 0.88.0.52.0.3 and 0.11. The last one is the capillary end point.

Fig.2 Oil accession distribution vs.

Many tiny channels are dead ends. It is well understood that water in dead ends is irreducible so that, in order to save the memory and CPU time, we cut off all dead ends from the pore network before we do two-phase flow simulations. It should be understood that our simulation results do not include the part of irreducible water that is trapped in dead ends. 3.

Physics of Two Phase Pore Network

Flow

in a CT-based

It is widely accepted that reservoir was initially filled with brine and water-wet. The oil from source invaded in and displaced brine. We use the same initial condition in our simulations. When one immiscible fluid (oil) displaces another immiscible fluid(water) in a channel, the driving force is the pressure difference AP between two nodes that are connected through this channel. When AP is greater than the capillary pressure associated with the oil-water interface, the interface moves. The capillary pressure is well formulated for a cylindrical tube (see Eq.l), but for a noncylindrical tube, the capillary pressure is not well defined. To get a reasonable expression of capillary pressure, we assume a piston-like displacement (that means that water trapped in pendular rings is excluded in our consideration) and employ a virtual work argument. Under a small displacement Ax, the interfacial area between fluid l(tbe displaced fluid) and the solid wall is reduced by an amount

Ax*$, where $ is the length of the contact line, and the interfacial area between fluid 2(the displacing fluid) and the solid wall increases the same amount. The free energy of system increases as: AF=(y2-y1).hx.d,

(2)

where yi is an interfacial tension between i’th fluid and solid. The energy increment AF equals the work AW done by pressure difference between two fluids, which is AF=AW=P,.h,cross_section.

(3)

Combining Eq.(2) and Eq.(3), we have a formula for capillary pressure, p, =

(

y2 - y,)

perimeter .

cross _ set tion

For a cylindrical tube, the general expression of capillary pressure Eq(4) is just Eq.(l) because of a force balance OCOS~= y2 - y,. We assume that the ratio of perimeter to cross-section is approximately equal to the ratio of the section volume to the section surface so that the capillary pressure, as the interface sits in the k’th section, is P,(k) = ocos8(S/v),

(5)

M. Zhou et al.: Irreducible Water Distribution in Sandstone Rock: Two Phase Flow Simulations in CT-based Pore Network

172

In the reservoir, the flow displacement is so small, the pressure drop to overcome viscous force is neglected in simulations so that all the nodes with same fluid have the same pressure. However, viscous force is not totally ignored. Quantitative calculation of the trapped water requires a formula for flow resistance that accounts for the relative speed of interface movement between channels. For a cylindrical tube, flow flux is proportional to the pressure difference and the inverse of proportional coefficient, the flow resistance, is R = 2 viscosity . (perimeter)2 . length H (cross_section)3

carefully. That requires a smaller time step. However, too short time step consumes CPU time enormously. For practical reason, the time step in our calculation is automatically adjusted such that the maximum position change among all moving interfaces is not larger than a fraction of one section. For every increment of the oil pressure, the interface movements are calculated through many time steps until the whole system reaches a force balance such that no interface can move.

(6)

If the shape of cross-section were different from a disk, numerical factor in Eq.(6) would be larger than 2 but than 3. A numerical factor of 3 is associated with extreme case of infinite thin rectangle. Assuming sections in all channels have a common shape-factor C, channel flow resistance is a sum of section resistance can be formulated by using Eq.(6),

the less an all the that

=

0

:

0.0

I

I

I

I

I

0.2

0.4

0.6

0.6

1.0

Water

Saturation

Fig.3 In capillary Curve plot, the capillary pressure is normalized by acos6Vpixel.

where ?J, AL, S and V are fluid thickness, surface and volume: S = perimeter

viscosity,

AL.

(8)

??

V = cross -section

section’s

AL

??

The cubic sample has six faces. Two opposite faces serve as inlet and outlet respectively and the remaining four faces are sealed in the simulations. The inlet is connected to an oil reservoir and the outlet is connected to a water reservoir. The simulation starts from the zero capillary pressure: the oil pressure equals the water pressure. We fix the water pressure and gradually increase the oil pressure. The oilwater interfaces all sit at the inlet before the oil pressure reaches a number such that one or more inlet-connected channels satisfy the criterion of interface moving, P,il

-

pwat~r

-

P,(k) ’ 09

with section-number condition is

k = 1. The channel

(9) flux

under this

Then the interface moves to a new position during time duration of At and the new position can be located using a simple mass balance calculation. In every time step, all channels are checked by the criterion (9). Water would be trapped in a channel as the two nodes connected by this channel both are accessed and filled by oil. The interface position in a channel at the moment when its end-node is filled up by oil coming from other channels determines the water-trapped-percentage in that channel. To avoid an overestimation of the trapped water, the relative moving speed of interfaces in channels should be calculated

4.

Results

We first run the simulation on a pore network that is constructed from the CT image of Fontainebleau sandstone with porosity 0.128 and permeability around 0.37 Darcy. We also run the simulations on the other 6 samples of Fontainebleau sandstone with porosity varying from 0.072 to 0.22. All the samples gave similar results. As an example, we just report in this article the detail information we obtain from the first sample. (1) Oil accession distribution

at different water saturation

Figure 2. shows the initial channel size distribution (V/S) and the size distribution of oil filled channels at four capillary pressures (a-d). Fig.2a and 2b illustrate that oil does not access the largest pore first as the pore-filling bundle model claims, which indicates that there are some small pores blocking the connection between the large pores and the inlet. Fig.2b, S,=O.52 is a breakthrough point where the displacing fluid oil reaches the outlet. Fig.2d, S,=O.ll represents an end point where all remaining water is trapped no matter how large the inlet pressure is. Fig.2 shows that at the end point water is left not only in small pores but also in some large pores. It is obvious that only a pore level by-pass mechanism can result in such distribution. Examining images of the trapped water distribution in the pore space, we found that most trapped water does not form clusters but rather is trapped in individual pores. (2) Capillary

Curve

One of the industrial standard measurements for rock properties is the capillary pressure curve. The centrifugation and porous plate techniques are two ways to obtain the capillary curve. The centrifugation method,

M. Zhou et al.: Irreducible Water Distribution

in Sandstone

Rock: Two Phase Flow Simulations in CT-based Pore Network

173

0.20 .-E ?

0.15

A

5 (I) 5 Z

A A

0.10

A

i

0.00 I-

4567

’

2

3

4567

’

2

3

1

0.1 Permeability Fig.4 Permeability dependence permeability is laboratory data from the simulation.

1 Pore Size (V/S)

(Darcy)

of the irreducible and the irreducible

water saturation: water saturation is

which requires less time, applies pressure to the sample as a body force. In principle, the capillary curve thus obtained does not give a real endpoint. All water will get out of the sample if the applied capillary pressure is increased. The underlining physics of the porous plate method are consistent with our simulation. The curve goes up vertically in the porous plate measurement when the water This is often saturation reduces to some number. misinterpreted as an operational error, capillary disconnection. The simulation we conducted gives the same behavior (see Fig.3) which might be helpful for resolving the confusion associated with the measurement. This result together with the simulations on the other five samples indicate that the irreducible water saturation decreases with increasing permeability as observed in the core measurements. Fig.4 shows how the calculated Swirr varies with the measured permeability.

(3) Irreducible

water distribution

Fig.5 shows the irreducible water distribution. The volume distribution of the trapped water appears lognoxmal, as the pore size distribution itself, rather than the one described by a simple cutoff. The irreducible water distributions from other seven samples have similar lognormal profiles. 1.

Discussion

The simulation we have conducted is different from the invasion percolation model although both are capillarycontrolled displacements. The difference first is in the aspect of the pore network where the calculations are implemented: the pore network invasion percolation model used is a lattice, where the bonds (pores) have same length but with randomly distributed radii, but our pore network is constructed from the real pore structure. Secondly, our simulation is deterministic but invasion percolation model involves a stochastic process: choosing a site next to the interface by a relative probability that is determined by the relative pore size among the next all available pores. The invasion percolation model ignores the long length correlation in the displacement and therefore results in a

10

Fig.5 Irreducible water distribution vs. the pore size (in unit of pixel) is represented by the filling area that appears a lognormal distribution. As a reference, the pore size distribution itself is plotted too.

trapped water distribution different from the deterministic displacement. The traditional picture of irreducible water distribution does not include the water trapped in pores whose size is larger than expected from the relationship between the capillary pressure and pore size. Typical S,,,, for sandstone rock with permeability about 0.5 D is around 0.2. The irreducible water saturation, 0.11, which we obtained from the simulations considering only the pore level by-pass mechanism, is a large fraction of the typical Swirr , which suggests that we should not ignore the possible phenomenon that the trapped water might exist in large pores. The displacement we simulated is a quasi-equilibrium process. For a dynamic invasion with finite flow rate, the trapped water saturation is expected higher than what we obtained. In principle, the end point saturation should increase with increasing flow rate (Lenormand et. al., 1988). There is a question about the irreducible water distribution: whether the water trapped in large pores would redistribute through the connected water path that is neglected by our assumption of the piston-like interface movement. The answer to the question depends on the time scale. The time scale of redistribution must be long because the relative permeability of the water path near the end point is close to zero. The redistribution of the trapped water may therefore not happen in the time scale of a laboratory measurement. This is supported by the similarity between the profile of NMR T2 distribution after centrifugation and the trapped water distribution from the simulation. If redistribution occurs in the geological time scale then the reservoir water saturation near the capillary end point must be much less than the core measurements and the calibration of the logging data with the core data would be questionable. To date, we have found little evidence against the effectiveness of the core calibration. We believe that there is a mechanism keeping the trapped water from being redistributed, for example the water path in the reservoir becomes disconnected due to a wettability change that is much faster than the water redistribution.

174

M. Zhou ef al.: Irreducible Water Distribution in Sandstone Rock: Two Phase Flow Simulations in CT-based Pore Network

Acknowledgements

The authors wish to thank Teng-fong Wong and Wenlu Zhu of the State University of New York (Stoney Brook) for providing the specimens in this study and for their useful discussions. Research carried out (in part) at the National Synchrotron Light Source, Brookhaven National Laboratorv.which is suooorted bv the U.S. Department of Energy, Division of Materials Sciences and- Division of Chemical Sciences under contract number DE-AC02-98CH10886.

Auzerais, F. M., Dunsmuir, J.H., Ferreol, B.B., Martys, N., Olson, J., Ramakrishnan, T. S., Rothman, D. H., and Schwartz, L. M., 1996, Geophysical Research Letters 23, No.7, 705

Lenormand, R., and Bories, S., 1980, C. R. Acad. Sci. B 291,279 Lenormand, R., Toubol, E., and Zarcone, C., 1988, J. Fluid Mech. 189,

165

Perona, P., Shiota, T., and Malik, J., “Geometry-Driven D@iiion in Computer Vision” B.M. ter Haar Romeny (Ed),Kluwer Acad. Pub.,73-92

Reference

Boston (1994). ISBN 0-7923-3087-O Bourbie, T. and Zinszner, B., 1985, J. Geophys. Res. 90, 11524 Martys, N.S., Hagedom, J.G., Goujon, G., and Devaney, J.E., Large

Chandler, R., J. Koplik, K. Lerman, and J. F. Willemsen, 1982, J. Fluid

scale simulations of single and multi component Jlow in porous media, in Developments in X-ray microtomography II, Ulrich Bonse, Editor.

Mech. 119,249.

Flannery, B.P., Deckman, H.W., Roberge, W.G., D’Amico, K.L., 1987, Science 237, 1389.

Proceedings of SPIE Vol. 3772, 205213 (1999).