ImageAnalStereol2009;28:23-34
OriginalResearchPaper
MODELINGOF MULTISCALEPOROUSMEDIA
BIBHUBISWAL1,2,P˚AL-ERICØREN3,RUDOLFJ HELD4,STIGBAKKE3ANDRUDOLF
HILFER1,5
1ICP, Universit¨ atStuttgart, Pfaffenwaldring27,70569St uttgart, Germany,2S.V. College,UniversityofDelhi,
NewDelhi-110021,India,3NumericalRocksAS,N-7041Trondheim,Norway,4StatoilHydro ASA,N-7005
Trondheim,Norway,5Institutf¨ ur Physik,Universit¨ atMainz,55099Mainz,Ger many
e-mail: biswal@icp.uni-stuttgart.de,peoe@numericalro cks.com,rujh@statoilhydro.com,
stig@numericalrocks.com,hilfer@icp.uni-stuttgart.de
(AcceptedJanuary30,2009)
ABSTRACT
A stochastic geometrical modeling method for reconstructi ng three dimensional pore scale microstructures
of multiscale porous media is presented. In this method the p orous medium is represented by a random
but spatially correlated structure of objects placed in the continuum. The model exhibits correlations with
the sedimentary textures, scale dependent intergranular p orosity over many decades, vuggy or dissolution
porosity, a percolating pore space, a fully connected matri x space, strong resolution dependence and wide
variability in the permeabilities and other properties. Th e continuum representation allows discretization at
arbitrary resolutions providing synthetic micro-compute rtomographic images for resolution dependent fluid
flow simulation. Model implementationsfor two different ca rbonaterocksare presented. The method can be
usedtogenerateporescale modelsofa wideclassofmultisca leporousmedia.
Keywords:carbonaterock,imagereconstruction,tomograp hy.
INTRODUCTION
Reservoir rocks are highly heterogeneous and
many of them contain a complex pore structure
with a wide range of length scales. One of the
main objectives of pore scale physics research is
to predict macroscopic petrophysical properties from
the underlying porous microstructure. Apart from the
basic understanding of the physics of flow processes,
thishasalsoenormouspracticalimportanceintermsof
improving uncertain estimates of recovery efficiency
on larger scales. Three dimensional pore scale models
oftherockareessentialforsuchresearch(Adler,1992;
Hilfer, 2000; Thovert et al., 2001; Arns et al., 2003;
OkabeandBlunt, 2007).
Amongthefrequentlyoccurringmultiscaleporous
media in nature, carbonate sediments are of great
economical interest because nearly half of the worlds
petroleum is found in carbonate reservoirs. Carbonate
rocks in general, and dolomites in particular, are
formed by a large number of different physico-
chemical processes and contain a wide variety
of diagenetic textures overprinting their primary
facies.Frequently,theoriginalprimordialdepositional
textures ( e.g., calcite ooids or bioclasts etc.) remain
visible in the final microstructure (Moore, 1989;
Lucia,1999;Moore,2001).Theycontainvuggypores,
defined as pores that are connected only through
interparticle porosity (Lucia, 1999). Dissolution,
nucleation and growth of crystallites are commonphenomenaduring dolomitization processes,and both
porosity and texture can fluctuate in the evolution
of the rocks. Unlike sandstones, permeability in
carbonate rocks can vary widely within a given rock
sample (Fernandes et al., 1996; Anselmetti et al.,
1998).Duetothesereasonsithasbeendifficulttofully
classify and characterizethe pore scalemicrostructure
ofcarbonaterocks(Dunham,1962;Lucia,1999).
The geometrical and petrophysical parameters of
carbonates depend strongly on resolution because
of a very wide range of pore and crystallite sizes
dominatingitsgeometrictexture.Fortypicalcarbonate
rocks,withsmallestcalciteordolomitecrystallitesizes
of 10−7m, one needs at least a voxel size of a≈
10−8mtoresolvethem.Evenforasmallcubicsample
of sizeL=10−2m, the digitized sample requires
prohibitive CPU time and storage space ( ∼1018
digits). To be able to model such a system and at the
sametimecapturethecomplexgeometry,acompletely
different modeling approach and data structure are
needed. We propose here a pore scale model in the
continuum that overcomesthese impediments and can
act as a starting point for more realistic modeling of
thesemultiscaleporousmedia.
Because the diagenetic processes that produce
the large variety of carbonate rock textures are
complex and largely unknown, we attempt a simple
stochastic geometrical model that tries to reproduce
the main features of the pore scale geometry. Any
23
BISWALBET AL:Modelingofmultiscaleporousmedia
modeling method for carbonates requires solving two
key problems: the “pore sizes” and grain diameters
ranging over many decades in length scales and
the carbonate textures showing a strongly correlated
disorder (Song and Sen, 2000). It is also crucial to
correlate pore scale microstructure to the depositional
and diagenetic fabrics and simultaneously distinguish
between intergrain, intragrain, intercrystalline and
vuggy pores. Porosity depends strongly on the shape
and size of the grains and grain packing. The sizes
of intergrain pores depend on the grain size sorting
and decrease as cement occludes the pore space or
as grains are forced closer together by compaction
(Lucia,1999).
The model proposed here aims to incorporate
simultaneously the above diagenetic processes
and the following microstructure features: scale
dependent intergranular porosity over many decades,
vuggy porosity, a percolating pore space, a fully-
connected matrix space, intergranular correlations
from primordial depositional textures and mineral
transformation. It allows discretizations at any
arbitrary resolution. In this paper the feasibility of
this model is tested by reconstructing the generic
microstructure of two different carbonate rocks with
multiple length scales.
THE MODEL
In this continuum growth model the porous
medium is represented by a random sequence of
points decorated with objects (grains or crystallites)
that correlate with the location and properties of the
differentcrystallitephasespresentintheoriginalrock.
This model belongs to the class of so-called germ-
grain models (Stoyan et al., 1995). The deposition of
the points is similar to a random sequentialadsorption
process (Feder, 1980). A stochastic geometer might
call the model a random-field-controlled germ-grain
modelofrandomsequentialadsorptiontype.
CONTINUUMREPRESENTATION
The state space of a rock sample containing N
crystallites occupying a bounded region S⊂R3is the
set
ΩN= (S×[Rmin,Rmax]×E×{1,2,...,g})N(1)
of all sequences ω= (ω0,ω1,...,ωN)∈ΩNwith
[Rmin,Rmax]⊂R1andE={x∈R3:|x|=1}. The
sample contains gdifferent crystallite phases. An
elementωi= (xi,Ri,ai,Ti)of the sequence represents
a crystallite of type Tiat spatial position xi∈Switha radius of the equivalent sphere Riand orientation
ai. These crystallite attributes shall depend on the
originaldepositionaltexturethroughtheirdistribution .
A probability distribution Pr on the space Ωof
sequencesfurtherspecifiesthemodel.
PRIMORDIALFILTER FUNCTION
The crystallite phases and intergranular
correlations from the underlying primordial
depositional textures are first captured by a greyscale
imageGdefined mathematically as a bounded, but
not necessarily continuous function G:S→[0,1].
This crucial input function G(x)is constructed
with information from geological analysis of the
rock, quantitative image analysis of two dimensional
micrographs and three dimensional µ-CT images, if
available.
The location and properties of the crystallites in S
arethendefinedthrough:
Ri=R(G(xi)),ai=A(G(xi)),Ti=T(G(xi)),
(2)
wherethefunctions R:[0,1]→[Rmin,Rmax]∪{0},A:
[0,1]→EandT:[0,1]→ {1,2,...,g}specify their
correlationswiththeprimordialtexture.Thepermitted
crystallite sizes in the interval [Rmin,Rmax]correlate
with the polydispersity in each crystallite phase of
the original sample. The rock sample is represented
as a random sequence of Npoints, each decorated
with crystallites satisfying the primordial correlation
throughEq.2.
DEPOSITIONOFPOINTS
The points xiare chosen randomly in S⊂R3and
added sequentially if the following overlap condition
is satisfied. For each chosen point xia corresponding
sphere radius Riis chosen as defined in Eq. 2 and we
setPr(Ωc) =0for
Ωc={ω∈Ω:∃i,j,O(ωi,ωj)/∈(0,λi]},(3)
where the degree of overlap between the spheres with
radiiRiandRjismeasuredby
O(ωi,ωj) =Ri+Rj−|xi−xj|
Ri+Rj−|Ri−Rj|,(4)
andO(ωi,ωj)is set to zero when RiorRjor both
vanish. The parameter λi=Λ(G(xi))defines the
maximum allowed overlap with the crystallite at xi
andΛ:[0,1]→[λmin,λmax]with 0 <λmin,λmax<
1. In other words, each newly added crystallite has
a finite overlap with an existing crystallite where
the degree of overlap is defined by the primordial
filter function. This ensures full matrix connectivity.
Porosityandporespaceconnectivityofeachcrystallite
24
ImageAnalStereol2009;28:23-34
phase depends on the degree of overlap, the density
of crystallites and the corresponding crystallite size
distribution. The density of crystallites depends on
the density of points, denoted by ρ, deposited in the
corresponding phase. This is a model parameter that
is determined from the image analysis of the original
rock. Porosity and pore space connectivity of the full
sample depends on the original pore phase and vugs
definedin G, inter-crystallite pore spacein eachphase
and the way these different crystallite phasescombine
to form thefull rocksample.
VUGGYPOROSITY
Apart from inter-grain and inter-crystalline pore
space, vuggy pore space is a distinctive feature of
many multiscale porous media. In general, vuggy
poresmaybetouchingornon-touchingandconnected
by other type of porosity (Lucia, 1999). To include
such vuggy pores, deposition of crystallites at vuggy
pore regionsarerestricted bysettingPr (Ωv) =0, for
Ωv={ω∈Ω:∃0≤i≤N;R(G(xi)) =0}.(5)
Thiscorrelatesthevugginesswiththedepositional
texture.
GRAINDECORATION
After depositing random points in a correlated
fashion, the deposited points are then decorated
with crystallites according to the depositional texture
(Eq. 2). A crystallite of type Tiwith orientation
aiis placed at each point xi. For a feasible and
fast convergence of the deposition algorithm, the
overlap condition in Eq. 3 is defined for spherical
volumes associated with each deposited point. Since
the decorated crystallites can be of arbitrary shape,
the spherical volume must either be inscribed within
the crystallite or the size of the crystallite di(Ri)must
be scaled appropriately so as to retain the matrix
connectivity.
DISCRETIZATION
Theporousmediasampleis fully representedbya
list ofNquadruples (xi,di,ai,Ti). It can be discretized
at any arbitrary resolution. In the discretization
procedure,acubicsampleofsidelength ℓissubdivided
into a grid of cubic voxels, each of side length a.
Within each voxel a set of selected collocation points
are chosen and the voxel is then assigned a value
depending on the fraction of these collocation points
that fall inside the deposited crystallites. A larger setof collocation points increases the accuracy of the
discretization,butiscomputationallydemandingwhen
Nis large.
NINE-POINTRULE
In this simple but reasonably effective
discretization procedure a voxel is marked as matrix
and assigned the label iif all of the following nine
pointsfall within the ith crystallite Gi,i=1,2,...,N.
pj=p+a
2(e1+e2+e3)+a
4tj;j=0,1,..,8,(6)
wherepis the position of the first corner of the voxel,
ais the sidelength of the voxel, t0= (0,0,0),t1=
(1,1,1),t2= (−1,1,1),t3= (1,−1,1),t4= (1,1,−1),
t5= (−1,−1,1),t6= (1,−1,−1),t7= (−1,1,−1)
andt8= (−1,−1,−1). These nine points are the
center and the eight symmetric points between the
center and the corners of the voxel. So, if the voxel
is fully inside a single crystallite then its label is the
crystallite number. If all the nine points fall inside
the pore space, i.e., outside all the Ncrystallites,
then the voxel is resolved as pore and assigned a
value 0. In all other cases the voxel status is labeled
asunresolved at the current resolution and assigned
the label −mifmpoints fall inside one or more
crystallites. This label reflects the matrix density
within thevoxel.Suchlabelingcriteria with crystallite
information is usefulfor building network modelsand
computing mechanical properties from the discretized
representation. The above discretization rule with
just nine collocation points provides a reasonably
good discretization for the analysis of the model.
A more accurate discretization method resembling
experimental µ-CT procedureispresentedbelow.
SYNTHETIC µ-CT
The number of collocation points in the
discretization rule is increased from 9 to n3. These
n3collocationpointsarethepointsofa n×n×ncubic
sublattice positionedcentrally inside eachvoxel. Each
voxel is labeled with an integer number mwith the
following meaning
–m=n3:allcollocationpointsarein matrix space.
–m=0:allcollocationpointsareintheporespace.
–0<m<n3:exactlymcollocationpointsareinthe
matrix space.
This procedure produces synthetic µ-CT
discretizations of the computer model of the rock.
Experimental µ-CT images contain noise and a direct
comparison may require addition of noise to the
synthetic µ-CT or filtering of the experimental µ-CT
or both. Digitized representations of the pore scale
can be obtained by choosing an appropriate threshold
25
BISWALBET AL:Modelingofmultiscaleporousmedia
mcand relabeling the voxels with label 0 <m<mc
to 0 (pore) and the voxels with label mc≤m≤n3to
1 (matrix). It is important to note that the CPU time
neededforthe discretizationprocedureis proportional
to the number of collocation points ( n3) and can be
prohibitive if the number of crystallites defining the
microstructure islarge.
COMPUTATION
Thediscretizationprocedurerequiresustocheckif
agivenpoint p∈Gi,i=1,2,...,N.Thecomputational
procedure for dolomite type crystallites, i.e.,Giis
rhombohedral, is as follows: A dolomite crystallite Gi
at position xiis defined by three pairs of intersecting
parallel planes. Each pair is separated by a distance di
and tilted by an angle αdegrees about the coordinate
axesto which they are parallel initially. The equations
ofthesethreepairsofplanesare
nj.x±di
2=0,j=1,2,3, (7)
wheren1= (cosα,0,−sinα),n2= (−sinα,cosα,0),
n3= (0,−sinα,cosα).
The rotations of the crystallites are implemented
throughgeneralizedquaternions.Foreach Gi,withthe
center of the crystallite at the origin of the coordinate
system, the orientation aiis defined by a sequence of
three rotations of θ1,θ2andθ3about the coordinate
axese1= (1,0,0),e2= (0,1,0)ande3= (0,0,1)
respectively.So
ai=A(G(xi)) =Qi3Qi2Qi1, (8)
where the unit quaternion Qijrepresents a rotation of
θjaboutthevector ej.
Thepoint p∈Giif
(nj.p′
r+di
2)(nj.p′
r−di
2)<0,j=1,2,3,(9)
wherep′
r=q−1
i(p−xi)(q−1
i)∗and(.)∗is quaternion
conjugation. The case α=0 corresponds to cubic
Gi. Checking p∈Gifor spherical shaped grains is
straightforward. The above procedure can be readily
generalized to other arbitrary shaped crystallites
definedbycutting planes.
MODELIMPLEMENTATION
Realisticmodelsofmultiscaleporousmediacanbe
generated using the algorithm described in Eqs. 1–5.
Actual implementations may differ from sample
to sample depending on the underlying primordialtextures, crystallite size distributions, type of
crystallites, location and type of vugs etc. Apart from
testing the feasibility of the microscopic modeling for
carbonate rocks, the following description presents
the detailed computational implementation of this
modeling procedure for two different carbonate
microstructures.
Fig.1.Left:Apart( ℓ=0.6mm)oftheoriginaloolithic
dolostone.Right:Apartofthereconstructedsampleof
thesamesize.
SAMPLE 1: DOLOMITIZEDOOID
GRAINSTONEWITHISOPACHOUS
DOLOMITECEMENT
Here we present in detail the specific
reconstruction of an oolithic dolostone (Biswal et al.,
2007). Only a two-dimensional photomicrograph of
the original sample was available as input data in this
case.ApartoftheimageisshowninFig.1.Itcontains
ellipsoidal primordial facies, complete dolomitization
andvuggyporosityoftouching-vugtype.
Imageanalysis
Due to insufficient resolution and small size
of the original image, computerized quantitative
image analysis is not feasible. Therefore, only visual
inspection and manual measurements are carried out.
The rock contains ellipsoidal ooid grains that are
fully dolomitized to two distinct crystallite phases:
a thin isopachous layer (rims of the elliptical ooids)
and the inner core. The diameters of the dolomite
crystallites in these two phases are roughly in the
interval (5 µm to 25 µm) and (1 µm to 10 µm)
respectively.Nodiscerniblecrystallitesizedistributi on
or overlap variation are visible in both phases. The
crystallitepackinginbothphasesshowsaspacefilling
by large number of smaller crystallites around bigger
crystallites. The crystallite packing in the inner cores
of ooids have large overlap. The ooid grains feature
diagenetic replacive dolomitization and dissolution
vugs.Thresholdinganalysisofthe imageindicatesthe
porosityofthesampleto lie in therangeof0.25–0.3.
26
ImageAnalStereol2009;28:23-34
Primordial filter function
In the absence of any three dimensional
information, the primordial filter function G(x)is
constructed from a computer model of polydisperse
sedimentation of elliptical ooids (disc-shaped) in
S⊂R3(Bakke and Øren, 1997). Vuggy pores were
artificially generated by removing a small fraction of
the deposited ooids and enlarging the remaining ooids
by 20% of their size to compensate the inter-ellipsoid
pore space and porosity. The kth ooid represented
by{rk,sk,ek}is centered at rk= (xk,yk,zk)with
semi-axes lengths ak=eksk,bk=eksk,ck=skand
k=1,2,...,1359. The primordial filter function is
definedas
G(x) =/braceleftBigg
1−(dk/sk)ifanydk<sk,
0 otherwise ,(10)
wheredk=/bracketleftbigg
(x1−xk)2
e2
k+(x2−yk)2
e2
k+(x3−zk)2
1/bracketrightbigg1/2
.Asmall
thin section of the primordial geometry is shown in
Fig. 2(top left).
Fig. 2.A thin section of the primordial ooid deposits
(top left), the deposited spheres for the isopachous
layer (top right), deposited spheres for the intraooid
regions(bottomleft)andafterthespheresarereplaced
by the corresponding dolomite crystallites (bottom
right).
Deposition ofpoints
The reconstructed cubic carbonate rock sample S
for Sample 1 has a sidelength ℓ=2 mm,i.e.,S=
[0,ℓ]3. Due to large polydispersity and Rmin/ℓ <<1,the continuum list contains millions of crystallites.
The deposition rule requires a computational scheme
for random deposition of overlapping polydisperse
spheres in Sthat obey a specific size distribution and
eachaddedpoint xi∈SmustsatisfyEq.3.
For computational efficiency we divide S⊂R3
into smaller non-overlapping cubic cells Ueach of
sidelength l,i.e.,S=U1∪U2∪U3...Points are
deposited randomly in each cell Uiand added to the
final list if Eq. 3 is satisfied. This scheme provides
fast convergence in deposition as Eq. 3 is verified
only for a small set of existing points in Uiand the
neighboring cells. Further, to ensure a space filling by
a large number of smaller crystallites around bigger
crystallites as observed in the original image, the
following crystallite size distribution is incorporated
into the deposit scheme. For the jth deposited point at
xjin each cell U, the excluded volume radius Rjis
chosenas
Rj=Rmin+/bracketleftbigg
1−j
n/bracketrightbigg
ξjΔR, (11)
whereΔR=Rmax−Rminandξjis a uniform random
number in (0,1). Although λi=Λ(G(xi)), for this
sample, we have assumed it to be constant, λi=λc.
In each cell n=ρl3number of points are deposited.
The packing density ρis first determined from trial
samples that match the target porosity and show pore
space connectivity. This is checkedusing the Hoshen-
Kopelman algorithm (Stauffer and Aharony, 1992) on
highresolutiondiscretizationsofthetrial samples.
Forthecrystallitesintheisopachouslayers,points
areaddedto Swith Pr (Ω1) =0,for
Ω1={ω∈Ω:∃0≤i≤N,G(xi)/ne}ationslash∈(0,0.125)},
(12)
withRmin=2.5 µm,Rmax=12.5 µm,λc=0.6,l=
62.5 µm and ρ=1.95×10−4. For crystallites in the
innercoreoftheooidgrainspointsareaddedto Swith
Pr(Ω2) =0,for
Ω2={ω∈Ω:∃0≤i≤N,G(xi)∈[0,0.1]},(13)
withRmin=0.5 µm,Rmax=5 µm,λc=0.65,l=
7.5 µm and ρ=3.395×10−2. The reconstruction
procedure for a small thin section of the sample is
shownschematicallyinFig.2.A2Dimageofasection
of the full 3D sample is shown in Fig. 1. The inter-
ooid vuggy porosities correspond to G(x) =0 and
is automatically introduced in Sthrough Eq. 12 and
Eq. 13. Some of the large vugsseen are introduced by
the removal of a small fraction of the ooids from the
primordial deposit.
27
BISWALBET AL:Modelingofmultiscaleporousmedia
Fig. 3.The fullcubicsample Sofsidelength ℓ=2 mm
discretized at different resolutions: a =20 µm (top
left), 10 µm (top right), 5 µm (bottom left) and 2.5 µm
(bottom right). Colorcodes:black(intra-ooidgrains),
grey (grains in isopachous layer), blue (pore space)
andpeach(unresolvedvoxels).
Graindecoration
For this model reconstruction, points
corresponding to the isopachous layer and the inner
cores of ooids are deposited separately and combined
in the following way. First the N1points in the
isopachous layer are decorated with equilateral
rhombohedra ( α=−15). The volume of each
dolomite Giequals the volume of the associated
excluded volume sphere. Although the sphere overlapdoes not imply crystallite overlap, in this sample
the matrix connectivity is retained due to the large
overlap and packing . The orientations θjare chosen
randomly from [−20,20]degrees. From the points
deposited in inner cores of ooids, we only chose the
pointsxi/ne}ationslash∈M1where M1=G1∪G2......∪GN1. The
matrix space Mof the reconstructed rock sample is
fully characterizedby a list of N=N1+N2≃4×107
depositedrhombohedralcrystallites.
Discretizationand property calculation
The reconstructedmodelis discretizedat different
resolutions using the 9-point rule described in
Eq. 6. Four discretizations are shown in Fig. 3.
As in real carbonates, with increasing resolution
of discretization, the fraction of the unresolved
voxels decreases and porosity increases significantly.
Table 1 lists the fraction of the resolved matrix, pore
and unresolved voxels for different resolutions of
discretization. The rough estimate of the porosity at
full resolution obtained by extrapolating these values
is in the targeted range 0.25–0.3. The resolved pore
scale structure shown in Fig. 3 resembles a typical
multiscale microstructure of oolithic carbonate rocks.
The strong resolution dependence of permeabilities
measured on the discretized samples is also listed in
Table 1. Simulation results of resolution dependent
pore space and matrix properties of the discretized
models have been presented in Biswal et al. (2007).
Absolute, single phase permeability and electrical
resistivity, or formation factor resistivity, capture
the pore space connectivity and geometry, whereas
elastic property calculations of bulk and shear moduli
recoverthe connectivityofthe matrix structure. These
property calculations lend themselvesto extrapolation
to infinite resolution, like the porosity determined on
our model at different discretization levels. This was
exemplified in Biswal et al. (2007), and suggests that
Table1.Fractionofresolvedporevoxels(f p), unresolvedvoxels(f u)andmatrix voxels(f m)atdifferentlevelsof
discretization. In column 5and6the absolute permeabilities computed on the discretized sa mples are listed. To
computethe permeability,the unresolvedvoxelsareeither convertedto matrix(K m)orto pore(K p).
Resolution fp fu fm Km Kp
(in µm) (in mD) (in mD)
20 0.12683 0.8463 0.02687 256.32 1.8467 ×107
13.333 0.13868 0.78332 0.078 1643.2 3.7239 ×106
10 0.14595 0.7154 0.13864 34951.0391 ×106
6.6667 0.15562 0.58397 0.2604 6539.5 2.0192 ×105
5 0.1628 0.4756 0.3616 8150.4 7.3075 ×104
4 0.16914 0.3941 0.43677 8520.3 4.0128 ×104
3.333 0.1749 0.3341 0.491 8570.9 2.7383 ×104
2.5 0.18506 0.254 0.561 8558.6 1.834 ×104
0 (extrapolated) 0.253 0.008 0.739
28
ImageAnalStereol2009;28:23-34
the investigated properties on multiscale structures
can be retrieved consistently by means of scale and
resolution dependentsimulationfrom themodel.
SAMPLE2:DOLOMITIZEDGRAINSTONE
WITHDIFFERENT OOIDPHASES
Here we present the detailed modeling procedure
for a second generic carbonate rock microstructure.
It is another dolomitized ooid grainstone for which
three-dimensional µ-CT discretizations are available
at three different resolutions. Separate cross-sections
for each resolution are shown in Fig. 4. The modeling
procedure aims not only to match the morphology but
also to incorporate calibrated porosities at different
levelsofresolution.
Fig. 4.Left window shows the cross-sections from the
experimentalµ-CT images of the original rock sample
at three different resolutions: sidelength of 1.5 mm
at 8.06 µm resolution (top), sidelength of 0.7 mm at
1.4 µm resolution (middle) and sidelength of 0.35 mm
with 0.7 µm resolutions (bottom). On right are the
cross-sections from the synthetic µ-CT images of the
reconstructedrocksampleatthesameresolutionsand
sidelengths.05010015020025000.0050.010.0150.02
n (greyscale value)P(n)
  
µ−CT (8.06 µm)
Model
05010015020025000.20.40.60.81
n (threshold)Porosity
  
µ−CT (8.06 µm)
Model
05010015020025000.0050.010.0150.02
n (greyscale value)P(n)
  
µ−CT (1.4 µm)
Model
05010015020025000.20.40.60.81
n (threshold)Porosity
  
µ−CT (1.4 µm)
Model
05010015020025000.0050.010.0150.02
n (greyscale value)P(n)
  
µ−CT (0.7 µm)
Model
05010015020025000.20.40.60.81
n (threshold)Porosity
  
µ−CT (0.7 µm)
Model
Fig. 5.Greyscale histogram (left, solid lines) and
threshold dependent porosities (right, solid lines)
of the carbonate rock sample computed from
the µ-CT images of size 10003voxels at three
different resolutions. The dotted lines correspond
to the processed synthetic µ-CT images from the
reconstructedrocksampleatthe sameresolutions.
Imageanalysis
Subsample µ-CT images of size 1000 ×1000×
1000 voxels were analyzed for porosity, crystallite
types, crystallite size distribution and packing. The
greyscale histograms and the threshold dependence
of the porosities are shown in Fig. 5 for these three
different resolutions. Although the 8.06 µm image
seems to have a representative volume for porosity
(sidelength of 8.06 mm), the matrix is not well
resolved and one observes a more or less continuous
histogram. Hence a clear cut-off for the porosity is
hard to determine at this resolution. The 1.4 µm
data set also seems to have a representative volume
(sidelengthof1.4mm).Itappearstoresolvethematrix
adequately as confirmed by the clear separation of the
two peaks in the corresponding greyscale histogram
shown in Fig. 5. A similar separation of the two
peaks is also observed in the histograms from 0.7 µm
resolution images. Ideally, the reference threshold
for determining the porosity of the rock should
be ascertained from the highest resolution images
available. However, in spite of the well resolved
matrix, the size of the 0.7 µm images (sidelength of
0.7mm)isclearlynolongerrepresentative.Therefore,
to determine the porosity of the rock at different
29
BISWALBET AL:Modelingofmultiscaleporousmedia
resolutions, we choosea greyscalevalue of 126 as the
reference threshold, determined from the separation
of the peaks in the histogram of 1.4 µm image. At
this threshold, the 8.06 µm data set has a porosity of
26.27%, the 1.4 µm data set has a porosity of 31.21%
and the two data sets with 0.7 µm resolution have
porosities of 30.03% and 31.62%, respectively. The
large fluctuations observed in the porosities of 0.7 µm
imagesareduethesmallvolumesofthesesubsamples.
Moreover, the µ-CT images at different resolutions
analyzed here are cut from different sections of the
rock.
Inabsenceofanyautomatedmethodforestimating
the intraooid porosity, rough areal analysis was done
for a selection of grains. At the 8.06 µm resolution,
no intraooid porosity can be discerned, except in the
hollow ooids. Based on visual inspection at different
resolutions, approximately 20 percent of the ooids
haveahollowinsideandaround10–50%oftheradius
is hollow.Thehollow ooidsaremostly elongated.The
intraooid porosity is found more or less similar at
1.4µm and0.7µm resolutions.
Fig. 6.A section of the primordial ooid deposits used
forthemodeling(left)andthecorrespondingsynthetic
µ-CT image of the reconstructed rock sample (right,
sidelengthof1.5mm atresolution8.06µm).
Primordial filter function
Acomputationalmodelofooidgrain(discshaped)
deposits (Bakke and Øren, 1997) was used as the
primordial geometry. In this continuum primordial
grain deposits, the grain size distribution was made
to match with that extracted from the µ-CT samples.
The ellipsoidal ooid deposits (Fig. 6, left) are all
perfectly oriented. The side length of the cubic box
with primordial deposits is 1 .5 mm. In order to avoid
unwanted border effect it is taken from the center of
the deposited cubic box of sidelength 1 .8 mm with
2152 ooid (disc shaped) deposits. After gridding at
300×300×300 voxels, the porosity of the disc grain
packisroughly34%.Thisestimateissomewhatlower
than the actual porosity of the continuum grainpack.
To compensate the intraooid porosity that will beintroducedlaterandstillbeabletoachievethetargeted
porosity,thesizesoftheprimordialgrainsareenlarged
by2%.Thegreyscaleprimordialfilterfunction G(x)is
thenconstructedusingEq. 10.
Based on the image analysis, five different kinds
of grain parameters (primordial phases) were chosen.
Theprimordialphase 1 correspondsto the isopachous
rim on all the grains. The remaining four primordial
phases are defined as follows. Let Hj
idenote the ith
primordial grain in the jth primordial phase ( {Hj
i,i=
1,2,..,mj,j=2,..,5}).mjis the number of grains
in primordial phase jandm=∑5
i=2mjis the total
number of primordial grains. 30% of the primordial
grainsarerandomlyselectedanddefinedas primordial
phase2.Theyaredenotedby {H2
i,i=1,2,..,m2}and
correspond to the tight intraooid crystallite packing.
20% of the randomly chosen primordial grains are
grouped into the primordial phase 3 corresponding
to the moderate intraooid crystallite packing and
30% of the randomly chosen grains are grouped into
theprimordial phase 4 corresponding to the loose
intraooid crystallite packing. The remaining 20% of
the primordial grains are grouped into the primordial
phase5 that contains grains with tight intraooid
crystallitepackingalongwithhollow,vuggycores.We
define new primordial filter functions separately for
eachoftheprimordial phasesas
Hj(x) =/braceleftBigg
G(x)ifx∈Hj
ifori≤mj,
0 otherwise ,(14)
whereG(x)is definedin Eq.10.
Deposition of centers and grain
parameters
Model parameters for different primordial phases
are chosen as follows and points were deposited
randomly and sequentially in S. For computational
efficiency we divide S⊂R3into smaller non-
overlapping cubic cells Ueach of sidelength l,i.e.,
S=U1∪U2∪U3...Pointsaredepositedrandomlyin
each cell Uiand added to the final list if the overlap
conditiondefinedin Eq.3 is satisfied.
1. Forprimordial phase 1 (isopachous rims of all
the ooids with moderate packing of intraooid
crystallites), we chose Rmin=5 µm,Rmax=
12.5 µm,λ=0.6,l=37.5 µm and ρ=1.8335×
10−3andsetPr (Ω1) =0,for
Ω1={ω∈Ω:∃0≤i≤N,H1(xi)/ne}ationslash∈(0,0.15)}.
(15)
2. Forprimordial phase 2 (30% of the randomly
chosen ooids with tightly packed intraooid
crystallites) we chose Rmin=3.75 µm,Rmax=
30
ImageAnalStereol2009;28:23-34
11.25µm,λ=0.8,l=33.75µmand ρ=5.0301×
10−3andsetPr (Ω2) =0, for
Ω2={ω∈Ω:∃0≤i≤N,H2(xi)∈[0,0.15]}.
(16)
3. Forprimordial phase 3 (20% of the randomly
chosen ooids with intermediate packing, i.e., less
than the moderate packing in the rim,) we chose
Rmin=3.75 µm,Rmax=11.25 µm,λ=0.5,l=
33.75 µm and ρ=2.012×10−3and setPr (Ω3) =
0, for
Ω3={ω∈Ω:∃0≤i≤N,H3(xi)∈[0,0.15]}.
(17)
4. Forprimordial phase 4 (30% of the randomly
chosen ooids with loosely packed intraooid
crystallites) we chose Rmin=3.75µ,Rmax=
11.25µm,λ=0.5,l=33.75µmand ρ=1.6767×
10−3andsetPr (Ω4) =0, for
Ω4={ω∈Ω:∃0≤i≤N,H4(xi)∈[0,0.15]}.
(18)
5. Intraooid vugs are introduced in the remaining
20% of the ooids deposited with parameters
Rmin=3.75 µm,Rmax=11.25 µm,λ=0.8,l=
33.75 µm and ρ=5.0301×10−3. Vug radius
varies randomly (with uniform distribution) from
10% to 50% of the corresponding ooid radius.
With each of the grains in the primordial phase
5, we associate a random number {ηi,i≤m5}
thathasafixedbutrandomvaluein [0.5,1].Weset
Pr(Ω5) =0,for
Ω5={ω∈Ω:∃0≤i≤N,H5(xi)/ne}ationslash∈[η(xi),0.85]},
(19)
whereη(xi) =ηjifxi∈H5
jforj≤m5.
The points are then decorated with equilateral
rhombohedra. Each rhombohedron Gi(before
rotation)centeredat xiistheintersectionofthreepairs
of parallel planes as defined in Eq. 7. The distance
diis chosen such that the volume of Giequals the
volume of the associatedexcluded volume sphere and
α=−15 degrees. Each rotation angle θjis chosen
randomly from [−60,60]degrees. Roughly 6 .4×106
rhombohedral crystallites are sufficient to fill the
whole 1 .5 mm ×1.5 mm ×1.5 mm cubic sample.
A cross-section of the discretized sample with these
decoratedcrystallites is showninFig. 6.
Fig. 7.Cross-sections from the experimental µ-CT
images of the rock sample (left) are compared
with sections from the processed synthetic µ-CT
images (right) of the reconstructed rock samples at
corresponding resolutions: sidelength of 1.5 mm at
8.06 µm resolution (top), sidelength of 0.7 mm at
1.4 µm resolution (middle) and sidelength of 0.35 mm
with 0.7 µm resolutions (bottom). The synthetic µ-
CT images are processed with added white noise of
standarddeviation19,17and15,respectively.
Discretization
The reconstructed samples are then discretized
with 125 collocation points (corresponding to n=5)
asdescribedaboveinthesectionentitled“ SYNTHETIC
µ-CT”. Sevendiscretizeddata setsare generated.One
datasetcontainsthefullsamplediscretizedat8.06µm.
Next,twonon-overlappingsubsamples,halfthesizeof
the full sample, were discretized at 1.4 µm resolution.
Thesetwosubsamplesarethelowerleftcornerandthe
upper right backwards corner of the full sample, i.e.,
the secondsubsample is displaced along the diagonal.
Finally, there are four non-overlapping subsamples
cut along the diagonal, each one quarter the size of
the full sample and discretized at 0.7 µm resolution.
Cross-sections of these greyscale discretizations are
31
BISWALBET AL:Modelingofmultiscaleporousmedia
visualized in Fig. 4 for all three resolutions. Porosity
was measured at different thresholds for all these
discretized samples. Threshold dependent porosities
averagedoverthesubsamplesarepresentedinTable2.
Porosities of the µ-CT images thresholdedat 126 (see
section “ SYNTHETIC µ-CT”) were compared with
model discretizations at corresponding resolutions for
thresholdvaluesof mc=48,52,and56.Theporosities
of the model discretizations correspond well to those
of theµ-CT images at mc=52. Like the 0.7 µm
resolution µ-CT images large porosity fluctuations
are also observed in the discretized samples at this
resolution. Morphologically, the discretized samples
showninFig.4alsomatchwellwiththe µ-CTimages
shownin Fig. 4.
Itisimportanttonotethattheabovediscretizations
represent synthetic µ-CT images from the model
and owing to the continuum representation, can
be generated at arbitrary resolutions. This allows
investigation of the rock properties at both
intermediate and higher resolutions than the available
µ-CT images. However, a direct comparability with
experimental µ-CTimagesrequiresadditionofasmall
amount of Gaussian white noise to these synthetic
images.Forexample,tomatchthegreyscaledensities,
we first shift the peaks and width of the greyscale
density histograms from the model to match with that
of the the experimental µ-CT images. Then we add a
small amount of noise (see Fig. 7) to match the peak
correspondingto the matrix phase.The discretizations
at resolutions 8.06 µm, 1.4 µm and 0.7 µm are added
with white noise of standard deviation 19, 17 and 15
respectively. These images are shown in Fig. 7 and
thethresholddependentporosityrevealsanimpressive
matchwiththatmeasuredfromtheexperimental µ-CT
imagesasshowninFig. 5.
Two-point correlation function
Quantitative characterization and comparison of
themicrostructure iscarriedoutthroughthetwo-point
correlationfunction(Øren etal.,2007)measuredfrom
the discretizations of the model and the original µ-
CT images at different resolutions. These greyscale
representations are thresholded to obtain a digitized
representation I(x). The indicator function I(x) =1 if
x∈matrixspaceand I(x) =0ifx∈porespace,where
xis the position vector of the voxel on the digitized
grid.Assuminghomogeneity,thetwo-pointcorrelation
functionC(r)is definedas
C(r) =/an}bracketle{t(I(x)−φ)(I(x+r)−φ)/an}bracketri}ht
/an}bracketle{t(I(x)−φ)2/an}bracketri}ht,(20)
whereφis the porosity of the sample. For the
following analysis C(r)is computed for raligned in
the X,Y andZdirections.0 50 100 150−0.200.20.40.60.81
lag (microns)Correlation function  
mCT_F
model_F
0 50 100 150−0.200.20.40.60.81
lag (microns)Correlation function  
model_F X
model_F Y
model_F Z
0 50 100 150−0.200.20.40.60.81
lag (microns)Correlation function  
mCT_H av
model_H av (subsample)
0 50 100 150−0.200.20.40.60.81
lag (microns)Correlation function  
mCT_Q av (subsample)
model_Q av (subsample)
Fig.8.Two-pointcorrelationfunctions(averagedover
X, Y and Z directions)calculated from the digitized µ-
CTimagesandthemodeldiscretizationsatresolutions
8.06 µm (top left), 1.4 µm (bottom left) and 0.7 µm
(bottomright).SamplenamesareasdefinedinTable2.
Top right: comparison of the correlation functions in
the individualdirectionsfrom the modeldiscretization
at8.06µm.Bottom left:Twodashedcurvescalculated
forsubsamples.Bottomright:Allcurvescalculatedfor
subsamples.
For this analysis, the µ-CT images (10003voxels,
8–bit greyscale) at resolutions 8.06 µm, 1.4 µm and
0.7 µm are digitized at threshold value of 126. The
model discretizations at resolutions 8.06 µm (full
sample, 1863voxels), 1.4 µm (two subsamples, 5363
voxels) and 0.7 µm (four subsamples, 5363voxels)
are digitized at threshold value mc=52. The two-
pointcorrelationfunctionsmeasuredonthesedigitized
representations are plotted in Fig. 8. At 8.06 µm
resolution, the model measurements agree reasonably
well with the µ-CT data in the short range with
marginal discrepancy in the long range. At 1.4 µm
resolution, the model measurements also agree well
with the µ-CT data. At 0.7 µm resolution, model
measurements from several of the subsamples agree
with the µ-CT data in the short range. The non-
representative sample sizes give rise to large sample
to sample fluctuations making a comparison at this
resolution difficult. The slope of the tangent drawn at
C(0)provides a rough estimate of the specific surface
areaofthesolid/voidandisrelatedtothepermeability
ofthesample(Hilfer,1996).Theaboveresultsindicate
that along with the quantitative agreement in the
porosities (Table 2) and in the morphologies (Fig. 4
and Fig. 7), the transport properties measured from
the model subsamples at different resolutions are
also expected to show reasonable agreement with the
subsamplesoftherockusedfor the µ-CT imaging.
Discrepanciesobservedbetweenthe µ-CT images
and the model discretizations may originate from a
32
ImageAnalStereol2009;28:23-34
Table 2.Variation in threshold dependent porosities of the µ-CT ima ges and the discretized samples (without
noise) at three different resolutions. µ-CT measurement fo r 0.7 µm resolution is averaged over two different
subsamples.Modelmeasurementsfor1.4µmand0.7µmresolut ionsareaveragedovertwoandfoursubsamples
respectively.
Resolution Original φ(n)Model φ(mc)
(in µm) µ-CTn=126 mc=48mc=52mc=56
8.06 mCTF0.2627 modelF0.2537 0.2615 0.27
1.40 mCTH0.3121 modelH0.3096 0.3133 0.317
0.70 mCTQ0.3082 modelQ0.3077 0.3097 0.3116
numberofdifferentsources.Themodelparametersare
onlyapproximate.Thedolomitesizesandtheiroverlap
are measuredfrom two dimensionalmicrographs.The
phases are determined from rough areal analysis of
someselectedgrainsintheexperimental µ-CTimages
shown in Fig. 4. The primordial ooid deposits differ
from the ooids visible in the µ-CT images in a
number of ways. The primordial ooids have circular
cross-sections in the Z-direction and elliptical cross-
section in the other two directions, thus the shape
of the primordial ooids match with the original rock
only in the X and Y directions. As a result of the
depositionprocedure(BakkeandØren,1997),theyare
perfectly aligned with respect to Z-direction, whereas
theorientationoftheooidsintheoriginalrockismore
isotropic. Due to these reasons a large discrepancy is
observed between the correlation functions measured
in the Z-direction and the other two directions for the
model (Fig. 8). This emphasizes the need for using
valid primordial input and model parameters in the
model reconstruction to achieve good correspondence
with the spatial correlations in the microstructure and
transportparameters.
CONCLUSION
In this paper we present a model for multiscale
porous media that provides a continuum description
of the porousrock at arbitrary precision.It reproduces
the enormous variation of “pore sizes” (vugs, pores,
cracks, slits etc) that is typical of multiscale porous
media. The range of length scales that can be
includedis limited onlyby thefloatingpointprecision
of computers. The petrophysical properties of the
reconstructed models can be studied as a function of
resolution without any limit in the magnification, i.e.,
scale dependenceof porosity, connectivity or network
models. The model permits to represent vugginess
(bothseparate,moldicandtouchingsuchasBrecciain
the sense of Lucia (1999). It reproduces intergranular
porosity, intragranular porosity, correlation with
primordialdepositionaltextureandfullconnectivityofpore and matrix space. The model parameters allow
calibration of porosities and physical parameters at
differentresolutions.
A successful reconstruction of a given rock using
this model depends on the accuracy of both the
input parameters obtained from image analysis and
the primordial filter function. The matrix connectivity
is ensured by finite crystallite overlap that is
implementedthrough anoverlapofspheresassociated
with each crystallite. Although it makes the algorithm
simple and fast, this approach may not be suitable
for rocks containing elongatedor complex objects. To
ensure faster and computationally feasible deposition
of millions of spheres, the algorithm presented in
this paper uses a subdivision of the rock into
sub-cells, filled sequentially. A faster algorithm for
depositing hundreds of millions of objects with
a given size distribution, satisfying simultaneously
the overlap rule, would be necessary to overcome
these restrictions. To the best of our knowledge, the
proposed modeling approach is currently the only
feasible method available for generating complex
multiscale carbonate pore-scale models of laboratory
scalerocksamples.
This processinspired continuum model represents
a feasible pore scale modeling technology for
multiscale porous media such as carbonates. Many
elements ofrealistic diageneticprocessesare included
in the modeling procedure. The model was tested
on two examples of oolithic dolostones and can be
easilygeneralizedtomodelawidevarietyofcarbonate
rocks. It can also be used to reconstruct other kinds
of multiscale porous media such as sandstones with
strongheterogeneitiesandfracturedporousmedia.
ACKNOWLEDGEMENTS
BBandRHgratefullyacknowledgefundingofthis
work byStatoilHydro ASA,Norway.
REFERENCES
Adler PM (1992). Porous Media – Geometry and
Transports.Boston:Butterworth-Heinemann.
33
BISWALBET AL:Modelingofmultiscaleporousmedia
Anselmetti FS, Luthi SM, Eberli GP (1998). Quantitative
characterization of carbonate pore systems by digital
imageanalysis.AAPGBulletin 82:1815–36.
Arns CH, Knackstedt MA, Mecke KR (2003).
Reconstructing complex materials via effective
grainshapes.PhysRevLett 91:215506.
Bakke S, Øren PE (1997). 3-d pore-scale modeling of
sandstonesand flow simulationsin pore networks.SPE
J2:136–49.
Biswal B, Øren PE, Bakke S, Held RJ, Hilfer R (2007).
Stochastic multiscale model for carbonate rocks. Phys
RevE 75:061303.
Dunham R (1962). Classification of carbonate rocks
according to depositional texture. In: W. Ham,
eds. Classification of carbonate rocks: American
AssociationofPetroleumGeologists,108–21.
FederJ(1980).Randomsequentialadsorption.JTheorBiol
87:237–54.
Fernandes CP, Magnani FS, Philippi PC, D¨ aian JF
(1996).Multiscalegeometricalreconstructionofporous
structures.PhysRevE 54:1734–41.
Hilfer R (1996). Transport and relaxation phenomena in
porousmedia.AdvChemPhys92:299–424.
Hilfer R (2000). Local porosity theory and stochastic
reconstruction for porous media. In: D. Stoyan,
K. Mecke, eds. R¨ aumliche Statistik und StatistischePhysik.LectNotesPhys254.Berlin:Springer,203–41.
Lucia FJ (1999). Carbonate Reservoir Characterization.
Berlin:Springer.
Moore CH (1989). Carbonate Diagenesis and Porosity.
Amsterdam:Elsevier.
Moore CH (2001). Carboante Reservoirs - Porosity
Evolution and Diagenesis in a Sequence Stratigraphic
Framework.NewYork:Elsevier.
Okabe H, Blunt MJ (2007). Pore space reconstruction of
vuggycarbonatesusingmicrotomographyandmultiple-
pointstatistics. Water ResourRes 43:W12S02.
Øren PE, Bakke S, Held RJ (2007). Direct pore-scale
computation of material and transport properties
for north sea reservoir rocks. Water Resour Res
43:W12S04.
SongYQ,RyeS,SenP(2000).Determiningmultiplelength
scalesinrocks.Nature406:178–81.
Stauffer D, Aharony A (1992). Introduction to Percolation
Theory.London:TaylorandFrancis.
Stoyan D, Kendall WS, Mecke J (1995). Stochastic
GeometryandIts Applications.Chichester: JohnWiley
&Sons.
Thovert JF, Yousefian F, Spanne P, Jacquin CG, Adler
PM (2001). Grain reconstruction of porous media:
Application to a low-porosity fontainebleau sandstone.
PhysRev E63:061307.
34