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
Comments 0
Log in to post a comment