PHYSICAL PROPERTIES - Home Pages of All Faculty at KFUPM

militaryzoologistΤεχνίτη Νοημοσύνη και Ρομποτική

1 Δεκ 2013 (πριν από 3 χρόνια και 4 μήνες)

70 εμφανίσεις

PHYSICAL PROPERTIES

Shipboard measurements of physical properties provide quantitative
information about the composition and lithology of core material and
are used to characterize lithologic units and to correlate core data with
downhole logging and seism
ic reflection data. All physical properties
measurements were taken on cores after they equilibrated to room
temperature (~25°C). Equilibration to room temperature takes 2
-
4 hr.
Magnetic susceptibility, gamma ray attenuation bulk density,
compressional wav
e (
P
-
wave) velocity, and natural gamma radiation
were measured on whole cores using the MST. Thermal conductivity
was measured on each core, using the whole core where possible.
After core splitting, undrained shear strength, index properties, and
addition
al measurements of
P
-
wave velocity were conducted on the
working half.

Multisensor Track Measurements

The MST, which is described in detail by Blum (1997), consists of four
sensors: the magnetic susceptibility logger, gamma ray attenuation
densiometer (GRA
),
P
-
wave logger (PWL), and natural gamma ray
detector (NGR). MST data were sampled at discrete intervals along the
core. The sample interval and the data acquisition period for each
sensor were set to optimize the resolution of data acquired within the
sa
mpling time available for each core. MST data are significantly
degraded if the core liner is only partially filled or if the core is
disturbed. When RCB or XCB drilling was used, the core diameter was
less than the nominal 6.6
-
cm core diameter. The reduce
d core
diameter required corrections of the values measured by the MST. The
values in the database do not reflect these corrections, but the figures
presented in the following chapters show corrected data.

Magnetic Susceptibility Logger

Magnetic susceptibi
lity is the degree to which a material can be
magnetized by an external magnetic field. If the ratio of magnetic
susceptibility is expressed per unit of volume, volume susceptibility is
defined as

=
M
/
H
,


where
M

= the volume magnetization induced in a material of
susceptibility (
) by the applied external field (
H
). Volume

susceptibility is a dimensionless quantity. It can be used to help detect
changes in magnetic properties caused by variations in lithology or by
alteration. Magnetic susceptibility was measured at 5
-
cm intervals
along the core using a Bartington meter (mo
del MS2C) with an 88
-
mm
coil diameter and a 2
-
s integration period. The Bartington meter
operates at a frequency of 0.565 kHz and creates a field intensity of
80 A/m (= 0.1 mT), significantly lower than the field intensity needed
to change the field orient
ation of magnetite grains (~50 mT). The
width of the instrument response to a thin layer of material with a high
magnetic susceptibility is ~10 cm. For this reason, the first and last
measurement of each core section was taken 4 cm from the core
section en
ds.

Gamma Ray Attenuation Densiometer

The GRA densiometer estimates bulk density by measuring the
attenuation of gamma rays traveling through the core from a
137
Cs
source. The gamma rays are attenuated by Compton scattering as
they pass through the sample.

The transmission of gamma rays
through the sample is related to the electron density of the sample by

Y
t

=
Y
i

x e
-
nsd
,


where

Y
t

= the transmitted flux,


Y
i

= the incident flux on a scatterer of thickness
d,


n

= the number of scatterers per unit volume o
r the electron density,
and


s

= the cross
-
sectional area per electron.


The bulk density (
) of the material is related to the electron density
(
n
) by

n

=
x N
AV

x (
Z
/
A
),


where

Z

= the atomic number or number of electrons,


A

= the atomic mass of the material, and


N
AV
= Avogadro's number.


Bulk density estimates are therefore

accurate as long as the ratio
Z
/
A

of the constituent elements is approximately constant and corresponds
to the ratio
Z
/
A

of the calibration standard. The GRA densiometer was
calibrated to a standard consisting of varying amounts of water and
aluminum so t
hat the densities of sediments can be accurately
determined. GRA density was measured using a 2
-
s integration period
at 5
-
cm intervals along the core.

Compressional Wave (
P
-
Wave) Logger

The compressional wave (
P
-
wave) logger (PWL) measures the
ultrasonic t
raveltime of a 500
-
kHz compressional wave pulse through
the core and the core liner. A pair of displacement transducers
monitors the separation between the
P
-
wave transducers, and the
distance is used to convert ultrasonic traveltime into velocity after
co
rrecting for the liner. Good coupling between the liner and the core
is crucial to obtaining reliable measurements. The PWL is calibrated by
placing a water core between the transducers. The PWL was set to
take the mean of 1000 velocity measurements over a

2
-
s period at 5
-
cm intervals along the core.

Natural Gamma Ray Detector

The NGR measures the discrete decay of
40
K,
232
Th, and
238
U, three
long
-
period isotopes that decay at essentially constant rates within
measurable timescales. Minerals that include K,

Th, and U are the
primary source of natural gamma rays. These minerals are found in
clays, arkosic silts and sandstones, potassium salts, bituminous and
alunitic schists, phosphates, certain carbonates, some coals, and acid
or intermediate igneous rocks (
Serra, 1984). The operation of the NGR
is outlined by Hoppie et al. (1994). The NGR system contains four
scintillation counters arranged at 90º angles from each other in a
plane orthogonal to the core track. The counters contain doped sodium
iodide crystal
s and photomultipliers to produce countable pulses. The
total response curve of the instrument is estimated to be ~40 cm and
so integrates a relatively long length of core in comparison to the other
instruments of the MST. Natural gamma ray emissions were
measured
over a 20
-
s period at 10
-
cm intervals. The NGR was calibrated in port
against a thorium source and during Leg 195 by measuring sample
standards at the end of operations at every site.

Thermal Conductivity

Thermal conductivity is the measure of the

rate at which heat flows
through a material. It is dependent on the composition, porosity,
density, and structure of the material. Thermal conductivity profiles of
sediments and rock sections are used, along with temperature
measurements, to estimate heat

flow. Thermal conductivity is
measured through the transient heating of a core sample with a known
geometry using a known heat source and recording the change in
temperature with time, using the TK04 system described by Blum
(1997). For soft sediment, the
rmal conductivity measurements are
made using a needle probe (Von Herzen and Maxwell, 1959) on whole
-
core sections; the reported value is the mean of three repeated
measurements. For materials too hard for the needle probe to
penetrate, thermal conductivit
y measurements are made after core
splitting, using the needle probe in a half
-
space configuration
(Vacquier, 1985); the reported value is the mean of four repeated
measurements. Thermal conductivity measurements were made at an
interval of at least one pe
r core unless variations in lithology required
more frequent sampling.

Undrained Shear Strength

The undrained and residual shear strength of sediments and
serpentinite mud was measured using a Wykeham
-
Farrance motorized
vane shear apparatus following proce
dures described by Boyce (1977).
In making vane shear measurements, it is assumed that a cylinder of
sediment is uniformly sheared around the axis of the vane in an
undrained condition. The vane used for all measurements has a 1:1
length to diameter blade
ratio with a dimension of 1.28 cm. A high
vane rotation rate of 90°/min was used to minimize pore fluid
expulsion while measurements take place. Torque and strain
measurements at the vane shaft were made using a torque transducer
and potentiometer. Undrain
ed shear strength measurements were
made at least once per core section unless variations in lithology
required more frequent sampling.

P
-
Wave Velocity

Discrete
P
-
wave velocity measurements were made in three directions
in the sediments using two pairs of
insertion transducers (PWS1 and
PWS2) with fixed separations of 7 and 3.5 cm, respectively, and a pair
of contact transducers (PWS3) in a modified Hamilton Frame. PWS1,
PWS2, and PWS3 use a 500
-
kHz compressional wave pulse to measure
ultrasonic traveltimes
, which, when combined with transducer
separation data, can be used to determine velocity. PWS1 and PWS2
were only used to measure velocity in soft sediments, where they were
inserted into the face of the split core. PWS1 is aligned with the core
axis (the

z
-
direction), and PWS2 is aligned perpendicular to the core
axis (the y
-
direction). PWS3 is mounted vertically with one transducer
fixed and the other mounted onto a screw, allowing the transducer
separation to be altered. PWS3 measures velocity in the x
-
direction in
split cores but is also used to measure velocity in discrete samples of
hard sediments or crystalline rock. Distilled water is applied to PWS3
to improve the acoustic coupling between the transducers and the
sample.
P
-
wave velocity measurement
s were made at least once per
core section.

Index Properties Measurements

Minicore samples of ~10 cm
3

were collected using a piston sampler in
soft sediment or an electric drill in rocks. Samples were taken at least
once per section. Sediment samples were
placed in a 20
-
mL beaker
and sealed to prevent moisture loss. Rock samples were soaked in
seawater for 24 hr before determining the wet mass. Samples were
then dried in an oven at 105° ± 5°C for 24 hr and allowed to cool in a
desiccator before measuring dr
y weights and volumes (method C in
Blum, 1997). Wet and dry sample masses and dry volumes were
measured and used to calculate wet bulk density, dry density, grain
density, water content, and porosity. Sample mass was determined
using two Scientech electron
ic balances. The balances are equipped
with a computerized averaging system that corrects for ship
accelerations. The sample mass is counterbalanced by a known mass
such that the mass differentials are generally <1 g. Sample volumes
were measured at least
three times, or until a consistent reading was
obtained, using a helium
-
displacement Quantachrome penta
-
pycnometer. A standard reference volume was included with each
group of samples during the measurements and rotated among the
cells to check for instrum
ent drift and systematic error; each time an
error was detected in the measurement of the reference volume, the
offending cell was calibrated. The following relationships can be
computed from the two mass measurements and dry volume
measurements (taken fro
m Blum, 1997, pp. 2
-
2 to 2
-
3). When a
beaker is used, its mass and volume are subtracted from the
measured total mass and volume. This results in the following directly
measured values:

M
b

(bulk mass),


M
d

(dry mass) = mass of solids (
M
s
) + mass of residua
l salt, and


V
d

(dry volume) = volume of solids (
V
s
) + volume of evaporated salt
(
V
salt
).


Variations in pore water salinity (s) and density (
pw
) that typically
occ
ur in marine sediments do not affect the calculations significantly,
and standard seawater values under laboratory conditions are used:

s = 0.035 wt% and


pw
= 1.02
4 g/cm
3
.


Pore water mass (
M
pw
), mass of solids (
M
s
), and pore water volume
(
V
pw
) can then be calculated:

M
pw

= (
M
b

-

M
d
)/(1
-

s),

M
s

=
M
b

-

M
pw

= (
M
d

-

[s x
M
b
])/(1
-

s), and

V
pw

=
M
pw

/
pw

= (
M
b

-

M
d
)/[(1
-

s) x
pw
].


Additional parameters required are the mass and volume of salt (
M
salt

and
V
salt
, respectively) to account for

the phase change of pore water
salt during drying. It should be kept in mind that for practical
purposes, the mass of salt is the same in solution and as a precipitate,
whereas the volume of salt in solution is negligible. Thus,

M
salt

=
M
pw

-

(
M
b

-

M
d
) =
[(
M
b

-

M
d
) x s]/(1
-

s), and

V
salt

=
M
salt

/
salt

= {[(
M
b

-

M
d
) x s]/(1
-

s)}/
salt
,


where the salt density (
salt

= 2.20 g/cm
3
) is a calculated value for
average seawater salt.

Moisture content is the pore water mass expressed either as
percentage of wet bulk mass or as a percentage of the mass of salt
-
corrected solids:

W
b

=
M
pw
/
M
b

= (
M
b

-

M
d
)/[
M
b

x (1
-

s)], and

W
s

=
M
pw
/
M
s

= (
M
b

-

M
d
)/[
M
d

-

(s x
M
b
)].


Calculations of the volume of solids and bulk volume are as follows:

V
s

=
V
d

-

V
salt

and

V
b

=
V
s

+
V
pw
.


Bulk density (
b
), density of solids or grain density (
s
), dry density (
d
), porosity (
P
), and void ratio (
e
) are then calculated according to the
following equations:

b

=
M
b
/
V
b
,

s

=
M
s
/
V
s
,

d

=
M
s
/
V
b
,

P

=
V
pw
/
V
b
, and

e

=
V
pw
/
V
s
.


Electrical Resistivity and Formation Factor

The electrical resistivity of the sediment was measured using a four
-
electrode configuration. The instrument used was modified at the
University of California, Santa Cruz, from th
e design of Andrews and
Bennett (1981) and was built at the University of Hawaii. The
electrodes consisted of four stainless steel pins that are 2 mm in
diameter, 15 mm in length, and spaced 13 mm apart. A 20
-
kHz
square
-
wave current was applied on the oute
r electrodes, and the
difference in potential between the two inner electrodes was
measured. The size of the current (typically 50 mA) was measured
over a resistor in the outer circuit.

The main purpose of measuring sediment resistivity was to determine
th
e formation factor, defined as the ratio of the resistivity of sediment
with included pore water divided by the resistivity of the pore water
alone. In practice, the formation factor is approximated by measuring
the apparent resistivity of the sediment in
the split core liner and
dividing that value by the apparent resistivity of seawater of similar
salinity and the same temperature in a 30
-
cm length of split core liner.
Using the same configuration for the measurement of the apparent
resistivities removes
the effects of geometry from the determination of
the formation factor.

Hydraulic Conductivity and Specific Storage

The hydraulic conductivity and specific storage of the serpentinite mud
was measured during a consolidation test. In this test, an axial sur
face
load is applied to a laterally constrained sample. The axial load
produces an excess pore fluid pressure along the length of the core.
The bottom of the sample is drained so that the excess pore fluid
pressure at that point is zero. The loads and boun
dary conditions are
applied by a Manheim squeezer, and the amount of fluid displaced is
measured as a function of time. Figure
F10

is a cartoon of the
ap
paratus and the boundary conditions. Also shown is the pressure
profile along the length of the sample at various times. The
assumption of incompressible mineral grains and water, common to
soil mechanics (Wang, 2000), allows the volume of water discharged

from the sample to be converted to axial displacement using the cross
-
sectional area of the sample. Because the frame, not the mineral
grains or the water, is compressed, we can calculate the axial
displacement using the following equation:

w

= volume of water discharged/cross
-
sectional area of sample,


where
w

= th
e axial displacement. We then use the relationship for
displacement in an infinite length cylinder as a function of time (Wang,
2000):

,


to determ
ine the lumped product of constants (on the right hand side
of the following equation) by plotting the slope of the displacement
over the square root of time

,


where

c
m

= the vertical compressibility,


= the loading efficiency,


z

= the axial load, and


D

= the hydraulic diffusivity.


Figure
F11

compares experimentally determined displacements with
calculated disp
lacements as a function of time. Only the early time
portion of the plot is used to determine the lumped product of
constants. Early in the experiment, the decrease in pore pressure has
not yet diffused to the end of the sample and so the approximation of
an infinite cylinder is still valid. To determine the hydraulic diffusivity
(
D
) from the lumped product, we need to determine the other
unknown factors,
c
m

and
.

The vertical compressibility is defined by

.


Because all the components of the equation above, axial displacement
(
w
), axial stress (
z
), and sample length (
w
o
) are measured, it is
possible to calculate the vert
ical compressibility. Also, because the
pore pressure throughout the entire length of the sample returns to
zero at very long times, the boundary condition of no change in pore
pressure (
P
pore

= 0) is met.

The loading efficiency is defined as

,


where
P
pore

= the pore fluid pressure. The assumption of
in
compressible grains and pore fluid leads to a value = 1 (Wang,
2000).

The specific storage (
S
s
) is related to the vertical compressibility under
the assumptions of incompressible grains and pore fluid by

S
s

=
c
m

x
f

x g,


where

f

= the fluid density, and


g = the acceleration of gravity.


Finally, we can determine the hydrauli
c conductivity from the hydraulic
diffusivity and specific storage using

K

=
D

x
S
s
.




PREDICTING PHYSICAL PROPERTIES OF RESERVOIR ROCKS
FROM

THE MICROSTRUCTURAL ANALYSIS OF PETROGRAPHIC
THIN

SECTIONS

M.C. Damiani
*1
, C.P. Fernandes
+2
, A. D. Bueno
+3
, L.O. E. Santos
+4
, J.A.B. da
Cunha Neto
+5
,

P.C. Philippi
+
6i

(*) Enginee
ring Simulation and Scientific Software (ESSS)

Parque Tecnológico de Florianópolis
-

Rodovia SC 401 km 001

88030
-
000 Florianópolis, SC, Brazil

(+) Porous Media and Thermophysical Properties Laboratory (LMPT)

Mechanical Engineering Department. Federal Unive
rsity of Santa Catarina BP476

88040
-
900 Florianópolis, SC, Brazil.

Abstract

Important physical properties of petroleum reservoirs are, usually, obtained by applying
standard

experimental tests on rock samples collected along selected depths of the petroleu
m well.
This is,

presently, a laborious and high cost procedure. Present paper describes Imago, which is a
new

petrographic analysis software, developed for predicting the physical properties of
reservoir rocks,

starting from petrographic thin sections rou
tinely produced in petroleum industry for
microscopic

analysis. Imago includes five main modules for i) segmentation of the porous phase from
gray level

and/or color digital images, ii) geometrical characterization of the porous microstructure,
iii) three
-
dimensional

reconstruction, iv) prediction of phase distribution in wetting
-
non wetting fluid

equilibrium displacement and v) prediction of intrinsic permeability. Simulation results
are

compared with experimental results for Berea and some other petroleum

reservoir
sandstones.

1
Marcos Cabral Damiani is a software development engineer from the Engineering Simulation and
Scientific Software

(ESSS) staff. Presently, he is working in the software development of Imago at Porous Media and
Thermophysical

Propert
ies Laboratory (LMPT).
2
Celso Peres Fernandes is a Visiting Professor from Brazilian Petroleum
Association (ANP), working at Porous Media

and Thermophysical Properties Laboratory (LMPT). His fields of interest include image analysis and
processing, 3D

geo
metrical reconstruction and multiphase flow in porous media

3
André Duarte Bueno has BS and MSc degrees in Civil Engineering. He is presently making his Doctoral
studies in

Porous Media and Thermophysical Properties Laboratory (LMPT) which includes equilib
rium phase
distribution and

prediction of relative permeability in porous media.

4
Luis Orlando Emerich dos Santos has a BS degree in Electrical Engineering and a MSc degree in
Physics. He is

presently making his Doctoral studies in Porous Media and Thermo
physical Properties Laboratory
(LMPT) which

includes the use of lattice gas automata (LGA) for the analysis of single and two
-
phase flow displacements
in 3D

porous representations.
5
José A. Bellini da Cunha Neto is Associate Professor at Federal Universit
y of
Santa Catarina. He specializes in the

macroscopic approach to fluid flow in porous media, image analysis and mercury intrusion processes.

6
Paulo C. Philippi is Professor at Federal University of Santa Catarina where he heads the Porous Media
and

Ther
mophysical Properties Laboratory (LMPT). He specializes in image analysis and processing, single
and

multiphase flow in porous media and lattice gas automata (LGA) methods.



Corresponding Author (philippi@lmpt.ufsc.br).
1. Introduction

Important physical properties of petroleum reservoir rocks are, usually, obtained by

applying standard experimental tests on rock samples collected along selected depths of
the

petroleum well.
In addition to sampling costs, laboratory routine includes the
manufacture of thin

plates for petrographic analysis and the measurement of porosity, intrinsic and relative
permeability

and mercury intrusion tests, which are, presently, laborious and high c
ost procedures.

Recently, the advancement of image analysis methods, applied on thin sections of
reservoir

rocks, appears to open a promising fast and low cost technique for predicting these
physical

properties, from the solely knowledge of the porous micr
ostructure of the reservoir rock.
In this

way, by using statistical homogeneity and spatial isotropy, the three
-
dimensional porous
structure is

supposed to be univocally related to two
-
dimensional representation obtained from thin
plates.

Furthermore, the
physical properties of reservoir rocks, in petroleum extraction, are
supposed to be

given from the porous microstructure geometry and the physicochemical interaction
between the

involved fluids and the porous surface.

The main aim of porous media study is
the description of equilibrium and transfer
processes

that are dependent on its porous structure. Common phenomenological methods use
adjustable

idealized models and/or experimental data to describe fluid occupation and the rates of
fluid

transfer. They ar
e not able to explore local geometric information of the porous structure.
Image

analysis techniques, starting from gray
-
level and/or color pictures taken with an electron
scanning

and/or optical microscope of polished porous sections, are suitable to the
analysis of the
geometrical

structure of porous media and could allow to the knowledge of their influence on fluid
equilibrium

and transfer. The main advantage of this approach is that the detailed knowledge of the
porous

structure enable to perform this s
tudy, starting from very fundamental physical laws
which are not

dependent on the porous structur
e.Virtual core analysis is a relatively recent research
field and has

been included in the SPE (Society of Petroleum Engineers) main research lines. At
authors


knowledge, the only existing operational software developed to perform core analysis
was

constructed by Porous Media Research Institute (PMRI)
(http://panda.uwaterloo.ca/pmri/) at

Waterloo University and named VCLab (Virtual Core Laboratory). VCLab start
s from
binary

segmented images and permeability calculations are based on percolation networks,
obtained from

the reconstructed three
-
dimensional representation.

Present paper describes Imago, which is a new petrographic analysis software, developed
for

pr
edicting the physical properties of reservoir rocks, starting from petrographic thin
sections

routinely produced in petroleum industry for microscopic analysis.

Imago includes six main modules for:

i) filtering and segmentation of the porous phase from gra
y level and/or color digital

images,

ii) geometrical characterization of the porous microstructure,

iii) three
-
dimensional (3
-
D) reconstruction,

iv) prediction of phase distribution in wetting non
-
wetting fluid equilibrium

displacement including mercury in
trusion,

v) prediction of intrinsic permeability based on skeleton method ,

vi) prediction of intrinsic permeability using lattice gas automata (LGA), by the

integration of local velocity fields.

The segmentation module includes a neural network method for

pattern recognition.

Geometrical three
-
dimensional reconstruction uses a new algorithm which starts from the

measured porosity and autocorrelation function and uses a Fourier transform method for
reducing

processing time [1].

Phase distribution is predict
ed by supposing quasi
-
static displacement by a wetting and/or
a

non
-
wetting invader fluid into the porous structure [2]..Intrinsic permeability is predicted,
in a first module, from the 3
-
D representation of the

microstructural skeleton [3]. Another permea
bility module include, a boolean lattice gas
automata

(LGA) for predicting permeability by integrating the velocity field inside the three
-
dimensional

representation of the porous structure.

Constructed in C++ object
-
oriented language to be a usable tool i
n petroleum engineering

practice, Imago is multiplatform (Win32 and X
-
Window) and user
-
friendly, with tools
for

generating and visualizing intermediate results: lists, graphics and image properties,
process log and

status. Special tools for three
-
dimension
al scientific visualization such as rendering,
rotation,

translation, zoom and two
-
dimensional cross
-
sections, are also provided.

Simulation results are compared with experimental results for Berea and some other

petroleum reservoir sandstones.

2. Geometri
cal analysis and processing: filtering, segmentation and 3
-
D
reconstruction

Imago was developed for predicting the physical properties of reservoir rocks, based on

petrographic thin sections routinely produced in petroleum industry for microscopic
analysis

(Figure 1). Imago starts from color and/or gray level digital images, obtained after,
respectively,

optical and/or electron
-
scanning microscopy performed on reservoir rocks thin plates.

2.1 Filtering

Small features in the image are frequently found as the

result of polishing procedures
and/or

image acquisition and appear as image noises. These noises are eliminated from the
analysis by

using spatial filters. The most commonly used filters are the
low
-
pass
filter, which
attenuates color

gradients and
median

filter that eliminates small discontinuities without any geometrical

consequence.

2.2 Segmentation

For fluid flow simulation purpose, it is only necessary to segment the original digital
image

in solid and porous phases, resulting in binary black and whit
e images. Thresholding is
performed,

by Imago, in accordance with the following methods.

2.2.1 Color images

Color digital images are derived from optical microscopy. Thin plates are prepared by

introducing a colored resin into the porous structure and the
porous phase is identified as
the

geometrical region related to the particular resin’s color. Digital images are frequently
acquired and

stored as 24 bits files in the RGB color model. In this model, each color is the result of a

combination of components
red, green and blue [4]. Imago presents two packages for
color images

segmentation. The first one works directly with RGB color model and the second one
starts from

HSI model, where the information of color is stored in the hue (H) component, while
compone
nts S

and I give information on saturation and intensity, respectively. Thresholding is
performed using

inter
-
class variance maximization method [5], applied to each component histogram in
both color

models. Taking petroleum industry needs, segmentation me
thod was automated, although
manual

segmentation is also enabled in this module, when, due to polishing and/or acquisition
problems,

original digital image has not a good quality. In this way,
zoom, hand
and
preview
tools
were

included in this module to si
mplify manual segmentation work (Figure 2)..2.2.2 Gray
-
level images

Imago includes three gray
-
level segmentation modules. The first module is manual and

based on the histogram threshold which appears to best separate the gray
-
level classes
associated to

so
lid and porous phase, respectively. The second module is automatic and based on the
inter
-
class

variance maximization method. The third module uses the global entropy method
introduced by [6]

and was, also, automated.

Imago also includes a neural network m
odule for both gray level and color images that is,

presently, in its developing phase. Preliminary results are very promising and the
intention is to

validate this method, which is very suitable for automating.

2.3 Measure of geometrical parameters

2.3.1
Autocorrelation function

In porous materials, one can theoretically distinguish between the solid and pore phase.
The

pore space of porous media can be characterized by the phase function
Z(x)
as follows:




otherwise 0

porespace the to belongs when

ì1
x Z(x)
(1)

where
x
denotes the position with respect to an arbitrary origin.

The porosity e, the autocorrelation function C(u), and the normalized autocovariance

function
R
Z(u) can be, respectively, defined by the following statistical averages (denote
d
by an

overbar):

e= Z(x) (2)

u)

(x)Z(x Z C(u)
(3)

) R
2 Z

e (e

e]




u) [Z(x [Z(x) (u)
(4)

Porosity

is obviously a positive quantity which is limited to the (0, 1)
-
interval. It can be

shown that a function
R
z(u) is a normalized covariance function if all its Fourier
comp
onents are

nonnegative. This work is restricted to homogeneous media, where statistical parameters
are

assumed to be independent of position
x
in space. Thus, the porosity is constant and
Z
R
(u)depends

only upon the vector
u
being independent of position
x. Moreover, when the porous
media is

isotropic,
Z
R is a function of only
u

u , and does not depend on the direction of u.

A simple method to calculate the correlation function and normalized covariance function

was used by [7], [8]. Let S be a section of a porous medium, given by a 2
-
D binary
representation,.porous phase is rep
resented in black and the solid matrix in white. The
binary image S is divided

into two halves
S1
and
S
2. Hence,

S=
S

È
S
2,
S

Ç
S
2=

(5)

In order to calculate
R
z(u),
S1
is first translated by a distance u along the x
-
axis; yielding

S
1(+u). The spatial average

indicated in Eq. (3) is calculated as an intersection of images,

Z(x y Z(x u y , ) , )

=
S
1(+u) ÇS (6)

giving the correlation function.

Since C(u) is related to the
probability
of finding two points separated by u and
belonging to

the same phase, it i
s, however, advantageous to calculate the autocorrelation function
C(u) as a

function of the two
-
dimensional vector
u
=(x,y) and, then, to take its mean value around a
circle

with radius u=
u
. This last procedure produces more reliable C(u) values because

it
increases the

number of realizations needed to calculate this probability. For an image f(x,y), the
Fourier

transform of the autocorrelation function is also the power spectrum of f(x,y) (Wiener
-
Khinchin

theorem). The 2
-
D autocorrelation function is ca
lculated in Imago by using the Fourier
transform.

Fluctuations are drastically reduced, when C(u) is calculated using Fourier transform.

In this way, with the Fourier transform of phase function, the correlation function can be

obtained, rapidly, using Fou
rier transform methods. Figure 3 show, a surface display
representation

of the power spectrum of a porous section binary image. The correlation function is the
inverse

Fourier transform of the power spectrum.

Given the two dimensional estimate ) y , x ( C
, we can obtain the desired one
dimensional

(isotropic) correlation function C(u) by averaging over the ) y , x ( C values at a fixed
radius u.

Except for the cases (0,u) and (u, 0), ) y , x ( C will not generally be known at the points
of interest

(see Fi
gure 4).

Figure 5 shows the correlation function obtained with Fourier transform method,
compared

with the previous displacement method.

2.3.2 Pore Size Distribution

In Imago, pore size distribution of binary 2
-
D sections is obtained by successive
openin
g,

derived from mathematical morphology and using balls with increasing radius [5]. The
resulting

image can be viewed as the union of balls completely enclosed in the porous phase. Thus,
after

opening, porous phase lost all the features smaller than the open
ing ball. The cumulative
porous

distribution is, then, given by:



e

e e
r r F

(7)

where
e
is the total porosity of the original image and

r
e
is the volume fraction of
the porous

phase after opening with a radius r ball..To reduce time processing, opening operation is
not applied directly to the binary im
age, but

to a transformed image called background distance imag
e.
In this image, each pixel is
labeled with

the smaller distance from it to the neighbor background. This labeling technique uses a
sequential

algorithm [9], where the Euclidean distance is ap
proximated by a discrete integer
distance. The

most commonly used discrete distance is the chamfer distance known as
d3
-
4, where each
neighbor

from a given point, taken following the horizontal and vertical coordinated axis, are
considered to

be 3 measurin
g units distant from the starting point. The diagonal neighbors are
considered to be at

4 measuring units from that point. Thus, this discrete distance gives the numerical
approximation of

4/3 = 1.333...to the square root of 2 (1.4142...). The main advanta
ge of using this discrete
distance is

related to the lower computer storage as only integers are stored in resident memory.
Balls

generated with
d3
-
4
distance that are used to perform opening operation are octagonal in
shape.

2.4 Three
-
dimensional reconstr
uction

Three
-
dimensional reconstruction is based on a modification of Adler
et
al. Gaussian

truncated method [7] proposed by [1]. In Adler ’s method, three
-
dimensional stochastic
simulation

of porous structure is based on the generation of a random non
-
cor
related field X(x),
which gives

the correct binary phase function Z(x), after the successive application of a linear and a
non
-
linear

filters. Linear filter gives Y(x), which is a convolution

Y(x)=a*X=



s

s) a(s)X(x (8)

inside a spherical domain with radius lequal to the correlation length, measured at the
original 2
-
D

binary image. Non
-
linear filter transforms Y(x), which is a real variable field to a binary
field

Z(x)Î{0,1}. By assuming homogeneity
and isotropy, 3
-
D pore structure can thus be
constructed

from 2
-
D porous sections, preserving porosity and autocorrelation functio
n.
Another way
to carry

out linear filter is to generate Y(x) from X(x) using Fourier transform [10]. From a
computational

poi
nt of view, the use of the fast Fourier transform algorithm, instead of laborious
solution of

nonlinear equation, makes the Fourier transform superior to linear filter method. Further,
although

Fourier transform
ˆ(p)
a
of a real field a(r) is complex, it c
an be easily show that in
discrete

representations, due to periodicity and conjugate symmetry properties of Fourier
transform,

computer resident memory requirements are the same for both fields. A truncated
Gaussian method

by using Fourier transform was pr
oposed by [1]. The difference between this method and
the

previous ones is that the field Y(x) is directly generated from its autocovariance function
R
Y(u). In

fact, by the Wiener
-
Khinchin theorem, the Fourier transform of the autocovariance of a
function
is

the power spectrum of this function, i.e.,

2 ^ 2

Y

^

Y
Y | R

(u)) (p)
(9)

Therefore, if the autocorrelation function is known for an arbitrary field Y(x), Fourier

transform can be used to generate this field, with the same autocorrelatio
n function. In
fact, the

above equation means that the Fourier transform R

^

Y
(p)
of
RY
(
u) is only related to the magnitude

^

Y of

^

Y

. This means that the phase angle of

^

Y does not depend on the autocovariance.R

^

Y
(p)
, or, in other words, th
at any two functions,

^

Y , with different phase angle may give the

same autocorrelation. This method does not need the linear filter and so avoids solving
the

nonlinear equations, reducing computer running time. Using the fast Fourier transform
makes this

algorithm more efficient. This idea was also used by [11] in reconstructing two
-
dimensional serial

sections in their hybrid approach. However, their method needs linear filter to generate 3
-
D porous

structure representations from these non
-
correlated seri
al sections. Present approach
avoids solving

a set of nonlinear equations associated with linear filter transform. In addition, when
compared with

Adler’s method [10], it reduces the resident memory requirements, because the
independent

gaussian field data

X(x) are not needed. Therefore, both operating time and computer
memory

requirements are improved. These are the advantages of the truncated Gaussian method
by using

Fourier transform used by Imago.

Figure 6 shows a schematic description of Liang
et a
l.’s

method used in the
reconstruction

module of Imago and Figure 7 presents an example of output of this package, showing a

comparison for the cumulative distribution of pore volume fraction between the values
measured on

the original binary image and the mea
n value obtained by measuring on several 2
-
D cross
sections

of the reconstructed microstructure.

Opening operation, section 2.3.2, was used in this computation step. Factor N is related to

the linear size of the 3
-
D spatial representation and factor n is a

sampling factor taken
over
R
Z(u)

data: sampling factor n means that the only points considered for
R
Z(u) data were the
ones separated

by n
-
1 pixels.

Results show that, when sampling factor n is reduced, three
-
dimensional reconstruction

gives a systematic
deviation with a larger amount of small size pores when compared to
the original

binary image. The better agreement was obtained, in present sample case, for n=5 and
n=6. In fact,

it must be remembered that present reconstruction method is based: i) on the

hypothesis
that the

target original two
-
dimensional representation, is a realization of a stochastic process
which is

considered to be ergodic and stationary and ii) on the hypothesis that this process is
Gaussian and

can be, inherently, described by its
only two first moments.

Imago has a three dimensional window for visualizing the three
-
dimensional
reconstructed

microstructure presented by iso
-
surfaces and/or iso
-
lines, with important operational tools
including

orthogonal cuts, zoom, translation and ro
tation. Figure 8 shows an example of this
software’s

feature.

3. Prediction of phase distribution in immiscible displacement

Determination of phase distribution inside the porous structure is an important problem in

oil recovery, giving reliable informatio
n for the prediction of recovery ratio and rate on
secondary

extraction.

Imago’s package for the prediction of phase distribution is based on [2], [12
].
Phase

distribution is simulated by supposing that equilibrium states are attained through quasi
-
static

processes. Method is based on the relationship between capillary pressure and the
curvature radius

of the interface given by Young
-
Laplace’s equation, when this interface is supposed to be
spherical:

i i

o
B A

AB i

P P

(d r



(10).where d is the Eu
clidean dimension of the space,
sAB
is the interfacial tension
between fluids A

and B,
i

o
A
P is the pressure at the domain of fluid A which is connected to a free region L, and
i
B
P
is

the pressure of fluid B, at step i.

In this way, for a given capilla
ry pressure
i

o
A
P
-

i
B
P , non
-
wetting phase will be found in the

geometrical region where it is possible to be place a ball with radius r
i
, given by the
above equation,

and problem is reduced to a geometrical problem that can be solved by image analy
sis
methods.

Method is fully described in [2], [12] and readers are referred to these papers for further
details.

Imago’s 3
-
D visualization window enable the software operator to follow fluid
displacement

in the course of invasion process. An example of su
ch a visualization is given in Figure 9
showing,

wetting fluid location in an imbibition process, at an intermediate pressure step.

3.1 Mercury intrusion simulation

Mercury intrusion on samples of reservoir rocks is a routinely performed experimental
test

in petroleum industry, giving information on porosity, size distribution of pore
constrictions and

thresholding constriction diameter.

Imago’s package for the simulation of mercury intrusion is the same used for the
prediction

of phase distribution explain
ed above considering a quasi
-
static displacement of non
-
wetting fluid

into a reconstructed porous structure, initially filled with an ideally compressible fluid.
Contact

angle with respect to non
-
wetting fluid is further corrected from the 180
o
value, use
d in
the

simulation, to the 140
o
value that is commonly accepted between experimenters.

Figure 10 gives a comparison of simulation results with experimental results for
Berea500

sandstone. Best sampling factor was n=5. For testing statistical homogeneity,

simulation
was

performed for three different linear sizes: N=60, 80 and 100.

4. Flow simulation

4.1 Intrinsic permeability

The simplest approaches for the calculation of intrinsic permeability, based on the idea of

conduit flow, ignore the fact that diffe
rent pores are interconnected with each other.
These are

called capillary permeability models, among which the so
-
called Carman
-
Kozeny model
enjoys

much greatest popularity. Another permeability model is the cut
-
and
-
random
-
rejoin
model. The

sample is secti
oned into two by a plane perpendicular to the flow direction, and the two
parts are

rejoined together in a random fashion. The flow rate in the capillaries is assumed to be
described by

a Hagen
-
Poiseuille relationship. In empirical permeability models perm
eability is usually
correlated

with some characteristic parameters. Network models have been used in the last four
decades. The

simplest network is a regular lattice, which consists of a regular arrangement of sites or
“pores”,

connected by bonds or “throa
ts”. Pores and throats are characterized by their radii, with
an

independent probability distribution function for each, and by the center
-
to
-
center
distance.

Furthermore, the lattice is characterized by its coordination number, the number of
throats meeti
ng

at a pore. Regular networks such as square, hexagonal, kagomé, trigonal, cube, and cross
square

were treated. Network models are all based on the some information of pore structure,
such as pore

size distribution and coordination numbers. But, these dat
a are almost assumed in the
previous

works.

Imago has two packages for permeability calculation. The first package is based on the

Skeleton Method and the second one is based on a Lattice Boolean Model..
4.2 Skeleton
Method

Imago uses skeleton method [3], w
here a detailed analysis and validation results can be

found. The skeleton of a porous structure is a connected line composed by all points p
that are

equally distanced from porous surface. In this way, it can be considered as the line giving
a first

appro
ximation to the
main flow line
through the porous structure. Skeleton is, thus, an
intermediate

level model, between permeability models based on percolation networks and models
based on the

integration of local velocity fields. In this way, considering, i
n addition, the meaningless
of

secondary flows
for the calculation of fluid transfer rate, in low Reynolds number flows,
Imago

calculation method is based on the following assumptions:

i) Skeleton nodes are considered as flow bifurcation points;

ii) Flow b
etween two successive nodes is considered to be one
-
dimensional, through

a variable area cylindrical duct.

The reconstructed 3
-
D porous structure is usually a cube with side length
Lx ´Ly ´
Lz. The

x
-
axis of the medium is chosen to be parallel to the direct
ion of macroscopic flow.
Impervious

boundary conditions are applied to the external surfaces of the cube that are parallel to
axes y and z.

The skeleton of the pore structure is obtained by using Ma’s thinning algorithms [3]. The
3
-
D

skeleton of pore struc
ture provides a real network. The hydraulic conductance of the fluid
in the

network can be calculated from the local hydraulic radius. Therefore, permeability can be
estimated

from the above data.

To compute intrinsic permeability, we need to define the hy
draulic conductance of the
fluid

in the network. Flow is assumed to be sufficiently slow. For each point p in a given link,
there is

associated a pore section in a plane which contains p and is orthogonal to the link. This
section is

obtained by performing

an oblique slice of unity thickness in the reconstructed 3
-
D
structure,

orthogonal to the link. Furthermore, one makes the simplification that the resistance to
fluid flow in

the link point p may be characterized in terms of the hydraulic radius
rH
of thi
s pore
section:

section cross of perimeter of length

section cross of area
rH


(11)

Thus, the conductance of a fluid in a point of pore link,
g
L, is given by Poiseuille’s
equation and

may be written as:

(p) g

4

H

L





(12)

where

is the fluid visco
sity and l is the length of the pore (here l =1).

The overall resistance to flow between the two neighboring nodes i and j, is




j

p=i L ij
(p) g

1

g

1 (13)

The flow rate of fluid between the two nodes,.) P (P g Q
j i ij ij

(14)

where
Pi
and
Pj
are t
he nodal pressures. Since the fluid is incompressible, mass
conservation

requires that

0 Q

j

ij


(15)

where j runs over all the links connected to node i. Equation (15) together with the
appropriate

boundary conditions form a complete solution to the st
eady flow of an incompressible
fluid in the

pore skeleton network. The equations are solved using successive over
-
relaxation:

i

j ij

j j ij

i

g

P g P





(16)

where

is a relaxation parameter.

By imposing a pressure drop DP across the skelet
on network and computing the resulting

single phase flow rate Q, the intrinsic permeability k of the pore network is calculated
from Darcy’s

law:

L L

Q k

z y

x






(17)

where Q is the volumetric flow rate, A is the normal cross
-
sectional area

of the sample, L
is the

length of the sample in the macroscopic flow direction, DP º
P
1
-
P2
is hydrostatic pressure
drop, and


is the viscosity of the fluid.

To illustrate the process to calculate permeability, a 2
-
D skeleton of the pore space is
shown

in
Figure 11, calculated at each porous point, by solving Stokes equations (Present paper
used a

lattice
-
Boltzmann method for the calculation of full velocity field).The links and nodes of
the graph

are, respectively, composed of points with exactly two neigh
bors and three or more
neighbors. Dead

ends are associated with stagnated flows and are eliminated from the skeleton. Figure
illustrates a

detailed porous region showing the skeleton superposed to the local velocity field, and
depicts a

geometrical region
where flow is stagnated, presenting a low speed vortex. It is seen that,
although,

with a strictly geometrical background, skeleton method
can capture
stagnation regions
by

associating dead
-
ends in the geometrical description of these regions.

Table 1 give
s a sample of validation result for Berea sandstone with a nominal
permeability

of 500 mD. Best results where found for sampling factors n=5 and n=6 which
corresponds to the

sampling factors giving the best agreement for pore size distribution (see Figure
7).

As secondary flows associated to flow deviations inside complex porous structures do not

significantly contribute to the main flow resistance when the Reynolds number is very
low, the

present method appears to be very suitable and easier to use when co
mpared with
methods based on

numerical solutions of Stokes equation, giving the full velocity field. In addition, it does
not suffer.from the well
-
known limitations of methods based on percolation networks. In
fact, the skeleton is

constructed trying to pr
eserve the fine details of the pore structure along the flow path
and can, thus,

better describe its influence on main flow.

4.3 Lattice Gas Cellular Automata Model

Imago includes a package for predicting intrinsic permeability of porous media based on

Lat
tice Gas Cellular Automata (LGA) methods. Details about the model and validation
results are

given in [13], a companion paper submitted to present meetin
g.
LGA is a relatively recent
method

developed to perform hydrodynamic calculations being object of con
siderably interest in
the last

years. The method is due originally to Frish, Hasslacher and Pomeau, [14]. In its simplest
form, it

consists of a regular lattice populated with particles that hop from site to site in discrete
time steps

in a process, often,

called
propagatio
n. After propagation, the particles in each site
interact with

each other in a process called
collisio
n, in which the number of particles and momentum
are

preserved. An exclusion principle is imposed in order to achieve better computation
al
efficiency.

In despite of its simplicity, this model evolves in agreement with Stokes equation for low

Mach and low Reynolds numbers. Some aspects make LGA methods very attractive for
simulating

flows through porous media. In fact, the method is explici
t in time and boundary
conditions in

flows through complex geometry structures are very easy to describe in LGA simulation.

Simulations are performed with integers, needing less resident memory capability and
boolean

arithmetic reduces running time. LGA is
, furthermore, very suitable for parallel
processing.

Taking into account the above features, LGA model was introduced as a package of
Imago

software:

i) to increase the precision in permeability calculations, using a model based on the full

velocity field
,

ii) for validation purposes in the construction of simpler (and faster) permeability

models.

Figure 12 shows a typical output from Imago’s LGA package for the prediction of a

Brazilian sandstone intrinsic permeability. Figure also shows pore size distrib
ution
measured as the

mean from 10 original 2
-
D binary image, compared with the mean value obtained by
measuring on

several 2
-
D cross sections taken from the reconstructed microstructures. Sampling factor
that gave

the best agreement for pore size distribu
tion was n=5, which was taken as a pattern value
for the

simulation with LGA model on reconstructed microstructures. Experimental permeability
value was

k=441 mD. Predicted permeability value was k=368 mD for a reconstructed porous
structure with

linear si
ze N=100 and sampling factor n=5 and k=462 mD for N=200 and n=5. In fact, in
addition

to sampling factor, the correct choice of linear size N is important to ensure statistical
homogeneity

in permeability calculations. In general, statistical homogeneity i
n flow calculations is,
only,

ensured, at a linear size N that is, often, greater than the value that assures geometrical
homogeneity

in the calculation of pore size distribution. Further validation results can be found in [13].

5. Conclusions

Prediction o
f physical properties of reservoir rocks needs the correct knowledge of their

three
-
dimensional porous structure and the correct understanding of fluids behavior
inside porous

space and of its physicochemical interaction with solid surfaces. These are, yet
, very
difficult and

actual problems, which call for a very intensive research work.

This paper presented Imago, a new petrographic analysis software that includes very
recent

analysis methods, derived from geometrical stochastic analysis and fluid mechani
cs,
intended to be

a usable tool in petroleum industry..Imago includes as it main new features:

i) a fast algorithm for three
-
dimensional geometrical reconstruction based on Fourier

transform, with optimized computer resident memory storage requirements;

i
i) a fast skeleton method for predicting intrinsic permeability, avoiding the well
-
known

limitations of percolation networks;

iii) a very precise method based on lattice gas cellular automata for the prediction of

intrinsic permeability, calculated from th
e integration of full velocity field inside porous
space.

All these methods were integrated in a single software written in C++ object oriented

language and based on a multiplatform library (Win32 and X
-
Window) which enables
Imago to be

run in a very large

class of computers, ranging from the very diffused personal computers
to high
-
level

Unix
-
based workstations. A very elaborated graphical interface with on
-
line output
windows,

hand tools, lists, export/import modules and 3
-
D visualization make Imago user
-
friendly
and

flexible.

Validation results were also presented for Berea and some brazilian sandstones, including

pore size distribution, mercury intrusion and intrinsic permeability.

Acknowledgements

The authors would like to thanks CENPES/PETROBRAS (Centr
o de Pesquisas e

Desenvolvimento Leopoldo A. Miguez de Mello) for providing the images and
experimental data of

reservoir rocks. C.P. Fernandes gratefully acknowledges financial support of ANP
(Agência

Nacional do Petróleo).

References

[1] Liang, Z. R., Fernandes, C.P., Magnani, F.S. & Philippi, P.C., 1998,
A Reconstruction

Technique of 3
-
D Porous Media by using Image Analysis and Using Fourier Transfor
m,
Journal of

Petroleum Science and Engineering, 2
1, 3
-
4, 273
-
283.

[2] F.S. Magnani P.
C. Philippi, Liang Z.R. & C.P. Fernandes, 2000,
Modelling two
-
phase

equilibrium in three
-
dimensional porous microstructures,
International Journal of
Multiphase

Flow, 26
(1
),
99
-
12
3.

[3] Z.R. Liang, P.C. Philippi, C.P. Fernandes, and F.S. Magnani, 1999,
SP
E Paper
56006:

Prediction of Permeability From the Skeleton of 3D Pore Structur
e,
SPE Reservoir
Evaluation &

Engineering, 2
(2), 161
-
168 .

[4] Gonzalez, R. C., Wood, R.E. (1992)
Digital Image Processin
g, Addison
-
Wesley
Publishing

Company.

[5] Coster, M. an
d Chermant, J.L., 1989,
Precis D’analyse D’image
s. Presses du CNRS,
Paris.

[6] Kapur, J.N., Sahoo D.K., Wong A.K.C., 1981, A New Method for Gray Level Picture

Thresholding using the Entropy of the Histogram
Computer Graphics, Vision and
Image

Processing 2
9
, 273
-
285.

[7] Adler, P.M., Jacquin C.G. and Quiblier J.A., 1990.
Flow in Simulated Porous Medi
a.
Int. J.

Multiphase Flo
w, 16: 691
-
712..[8] Philippi, P.C., Yunes, P.R., Fernandes, C.P.,
Magnani, F.S., 1994,
The Microstructure of Porous

Building Materials:
Study of a Cement and Lime Morta
r,
Transport in Porous Medi
a,
1
4, 219
-
245

[9]. J. M. Chassery, A. Montanvert,
Géometrie Discrète en Analyse d’Image
s, Editions
Hermes,

Paris, 1991.

[10] Adler, P.M., 1992.
Porous Media: Geometry and Transport
s. Butterworth
-
H
einemann, New

York.

[11] Ioannidis, M.A., Kwiecien M. and Chatzis I., 1995. Computer Generation and
Application of

3
-
D Model Porous Media: From Pore
-
Level Geostatistics to the Estimation of Formation
Factor.

Paper SPE 30201 presented at the Petroleum Compu
ter Conference, Houston, TX.

[12] P.C. Philippi, F. S. Magnani and A.D. Bueno
Two Phase Equilibrium Distribution in
Three
-
Dimensional

Porous Microstructure
s. Submitted to Produccion 2000 / Aplicaciones de la

ciencia en la ingeniería de petróleo,
May 08
-
12
/2000, Foz de Iguaçu.

[13] L.O. E. Santos, P.C. Philippi, M.C. Damian
i.
A Boolean Lattice Gas Method for
Predicting

Intrinsic Permeability of Porous Medi
a.
Submitted to Produccion 2000 / Aplicaciones
de la

ciencia en la ingeniería de petróleo,
May 08
-
12 /2
000, Foz de Iguaçu.

[14] Frisch, U., Hasslacher B., Pomeau Y., 1986,
Lattice
-
Gas Automata for the Navier
-
Stokes

Equation,
Physical Review Letters, 5
6, 1505
-
1508..
Figures and Tables

Figure 1. Imago methodology.

Figure 2. HSI segmentation module for color im
ages..Figure 3. Surface display of power
spectrum for a porous section binary image.

Figure 4. Obtaining the isotropic correlation function from discrete values for a particular
image..
0

0.05

0.1

0.15

0.2

0.25

0 50 100 150 200 250

displacement u (pixels)

C
( u)

x direction

y direction

average

Fourier transform

Figure 5. Comparison of correlation function.

R
Z(u)
R
Y(u) R
Y

^

(p)


^

Y(p)

+

Random phase

angle

^

Y(p)

Y(x) Z(x)


-
1

Non linear

filter

Figure 6. Schematic description of Liang et al.’s method for 3D reconstruction [1]..
0

0.2

0.4

0.6

0.8

1

1.2

0.00 20.00 40.00 60.00 80.00 100.00

Pore radius mm

Cumulative volume fraction

n=6, N=200

n=5, N=200

n=4, N=200

n=3, N=200

n=2, N=200

original

Figure 7. Comparison for the cumulative distribution of pore volume fraction of Berea
sandstone

500mD, between values measured on the original binary image and the mean value
obtained by

measu
ring on several 2
-
D cross sections of the reconstructed microstructure.

Figure 8. Imago three
-
dimensional visualization window..Figure 9. Simulation of
wetting
-
fluid invasion (in light) into a reconstructed 3D structure.

0

0.2

0.4

0.6

0.8

1

1.2

0 10 20 30
40 50 60 70 80

Ball radius ( mm )

Hg saturation

Experimental

N=100, n=5

N=80, n=5

N=60, n=5

Figure 10. Comparison of simulation results with experimental data for Berea 500 mD
sandstone..Figure 11. 2
-
D skeleton of the pore space superposed to the local vel
ocity
field, illustrating

a geometrical region where flow is stagnated, presenting a low speed vortex.

0

0.2

0.4

0.6

0.8

1

0 5 10 15 20 25 30 35

radius [pixel]

Cumulative volume fraction

Original, k=441

n=5, N=100, k=368

n=5, N=200, k=462

Figure 12. Pore s
ize distribution and predicted LGA permeability values for a

typical Brazilian sandstone. Experimental value measured at

CENPES/Petrobras was k=441 mD..Table 1. Comparison of the result for Berea

sandstone 500 mD for different sampling factors n

and linear

size N. Table also shows the number of

skeleton nodes where pressure is to be calculated in

the simulation.

n N Node

numbers

K (mD)

60 2,098 169.1

80 5,187 181.4

100 10,957 226.4 4

120 19,323 335.7

40 697 209.4

100 13,125 465.4

110 18,390 478.3

120 24,805 432.9

5

150 48,422 452.2

40 745 362.5

60 2,748 433.7

80 6,897 440.1

100 14,197 439.4

6

120 25,300 459.5

30 327 933.5

40 909 946.8 8

60 3,490 928.3