SPE 166037

Minimization of Disparity between Well test and Log Driven Permeabilities

through Combination of Rock-physics And Advanced Logs with Proper
Averaging Technique
Sergey Vorobiev, PETRONAS Carigali Sdn Bhd, Michael M. Altunbay, PETRONAS Carigali Sdn Bhd, Vladimir
2
Vorobyev , Lukoil Overseas, Azlan Shah M. Johari, PETRONAS Carigali Sdn Bhd.

Copyright 2013, Society of Petroleum Engineers

This paper was prepared for presentation at the SPE Reservoir Characterisation and Simulation Conference and Exhibitionheld in Abu Dhabi, UAE, 16–18 September 2013.

This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not been
reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect any position of the Society of Petroleum Engineers, its
officers, or members. Electronic reproduction, distribution, or storage of any part of this paper without the written consent of the Society of Petroleum Engineers is prohibited. Permission to
reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright.

Abstract
Minimization of disparities between welltest and log-derived average permeabilities has always been an issue, particularly in
carbonates where complex pore structures add on challenges to permeability estimation from wireline log data. The
disagreement between permeability averages from logs and well tests originates from the combined effects of measurement-
scale of static porosity components for permeability models, dual flow system of fractures and matrix, tensorial nature of
permeability and the averaging techniques used.

The proposed workflow exploits rock-physics templates to identify and to quantify secondary porosity. Rock-physics
templates employ conventionally derived total porosity and shear modulus as inputs. Fracture and vug porosity identified by
the proposed workflow through rock-physics agree with other qualitative and quantitative evidences of non-primary porosity
obtained from NMR, Image logs and core data1. Matrix and connected-vug permeabilities are computed, calibrated and
integrated via “Chen-Jacobi” connectivity-driven model2 by using NMR and acoustic log data. Fracture permeability is
estimated from “fracture aperture” and fracture-porosity by using image log data and rock-physics algorithms. The final
permeability profile is computed with a selective-replacement step. This step ensures that in co-presence of matrix, connected
vugs and fracture permeabilities at a given discrete depth level, the greater one would dominate and replace the lesser one. The
final step in efforts of lessening the disparity between averages of wireline-driven and well test/DST permeabilities for a given
interval is the usage of proposed averaging technique for the integrated wireline-driven permeability profile.

The rock-physics templates used in this study combine Kuster and Toksoz3 “inclusions” theory with the Dvorkin-Nur4
granular media model (1996). We have observed appreciable correlations between secondary porosity driven from shear
velocities against the secondary porosity determined from NMR and Image logs and core data. These correlations further
provide routes for newer permeability models that can be solely based on the rock-physics.

Comparisons of permeability averages computed from wireline-driven permeability profiles against DST or welltest
permeability showed significant improvements toward parity via proposed methodology and averaging technique. The
workflow presented in this study is to guide the reader through numerous steps of the proposed algorithm in detail.

Introduction
Permeability is the key parameter in understanding producibility of a field. Being the coefficient in the Darcy flow equation, it
can be estimated directly either from welltest, RFT or from wireline data by modeling. Disparities between these three
estimations originate from the differences in volumes involved in investigations of wireline, DST and well test as well as
differences in flow geometries. The tensorial nature of permeability is not reflected in the wireline-driven permeability because
it is driven from model equations with no directional /dynamic data involved. However, flow in DST/well test takes place
from a volume of rock with the contributions from directional components to yield the transmissivity that is translated into
permeability with the knowledge of thickness of the interval tested and viscosity of the hydrocarbon produced.
2 SPE 166037

Observations in most of the cases show that horizontal permeability is greater than the vertical one. However, this assumption
may not hold true for specific cases where a fracture network can promote the vertical permeability; therefore, generalizations
may bring surprises that may be too costly to ignore.

Without aiming to cover all the issues (measurement type, scale and volume and tensorial nature of permeability) mentioned
kv
above, let’s assume that permeability anisotropy ( ) is known from core measurements and the disparity between wireline
kh
and DST/well test driven permeability values are negligible. Then, it is possible to predict the outcome of DST/well test in
terms of potential rates with application of petrophysical data (saturations, irreducible saturation and permeability). Ability to
make such predictions on the planning stage of a DST/well test, as fast and as precise as possible, can translate into large
savings by not running it or justification of the cost that is incurred from rig-time and DST/well test.

Derivation of permeability index as a continuous profile is possible only from the wireline logs with a requirement of
calibration. Unfortunately there is no logging tool capable of measuring permeability directly. Numerous models are available
in Petrophysics to generate a continuous permeability profile by tying permeability to log-driven rock properties such as
porosity, clay volume, irreducible water saturation etc. However most of the models are for clastics reservoirs with
predominantly intergranular flow conduits. Carbonates pose additional challenges to permeability prediction due to their very
nature where matrix intergranular porosity may co-exist with secondary porosity, fractures and connected vugs. Earliest
attempts to predict matrix permeability in carbonates were linking permeability with the pore-size either through rock fabric
(Lucia et al 2001)5 or through the porosity and irreducible water saturation (Chilingarian et al, 1996)6. Multiple models were
also published describing and quantifying fracture permeability, either as a stand-alone fracture or as a fracture-system.
Finally, Chen-Jacoby model for the vuggy carbonates was published, where identification of the vuggy zones was tied to the
NMR and sonic log data. Timur-Coates equation was modified accounting for the degree of connectivity between the vugs to
estimate permeability for the matrix and vuggy system. However all the models mentioned above fell short of emphasizing on
the integration of the fracture and vuggy-matrix dual system flow.

The identification of fractures and vuggy intervals and quantification of secondary porosity for vuggy and fractured systems
have also been problematic.The earliest known petrophysical approach is quantifying the positive difference between density-
neutron and acoustic porosity as the secondary porosity based on the postulation of acoustic waves travelling through the
interconnected matrix frame and by-passing the isolated vugs; hence, under calling the porosity as the connectivity between
the vugs diminishes. Chang7 used NMR logs to identify vugs assuming that vugs show up in late time of T2- distributions
because they have large pore-body diameters. NMR T2 cutoff 750ms was introduced to identify vugs in a given T2-
distribution.

Both techniques have merits and demerits. T2-

distribution driven porosity with a cutoff at 750ms and
above represents large pore-sizes but the cutoff is
indifferent to the pore-connectivity and T2-distribution
pattern is sensitive to the presence of hydrocarbons. On
the acoustic side, P-sonic is sensitive to saturation,
especially to the light hydrocarbons that affect results as
well. It is also worth to mention that linear response
1.a 1.b 1.c between porosity and sonic transit travel time,
(represented by Wyllie equation) fits well only into
either moderately-cemented or sands with high effective
pressure trends (over 30mPa). Linear relation between
porosity and sonic travel time deteriorates as the degree
of cementation increases or the rock becomes
unconsolidated. This problem was observed very early
by Raymer et al (1980)14. In addition; neither of the
approaches is conducive to quantifying fracture-
1.d 1.e 1.f occupied pore-space. Thus, rock-physics approach is
Fig. 1 Vugs, macro, micro fractures can be seen on core and image log studied in this paper to overcome the aforementioned
difficulties.
In the last section of this paper, we will explain a better averaging technique for wireline-driven permeability for minimizing
the disparity between the DST/well tests driven permeabilities against the wireline-driven permeability averages. In addition, a
comparison of fluid flow rates from wireline-driven data against DST/well test reported rates is also presented for
completeness and verification of appropriateness of the proposed methodologies and averaging techniques. We used the Power
SPE166037 3

Law averaging8 and its specific case MLE (Most Likely Estimates) 9 based on the specifics of permeability distributions. Both
techniques have been tested and will be discussed in the consequent sections of the paper.

Proposed workflow and methodologies have been tested on two different carbonate fields from two different regions. First
field is located in Indonesia (Madura basin). The studied reservoirs are comprised of grainstones to mudstones, predominantly
wakestones. They are karstified and fractured as the results of unique tectonic events. Vugs, macro and micro scale fractures
can be easily seen on the core and Image logs (Fig. 1). Comprehensive log suites are available in quite a number of wells
(Neutron, Density, GR, NMR, Dipole Sonic, and Image logs) together with good core coverage creating an excellent data set
for the study. Second field is in the Middle East region. Subject reservoirs of the second field comprised of bundstones,
grainstones and wackestones, showing significant vuggy porosity that are easily identified on the core samples and Image logs.
The available data set includes Neutron, Density, GR, NMR, Dipole Sonic and Image logs. It should also be noted that
acoustic properties together with visually measured vuggy porosity are also available on the core samples.

Rock-physics model and secondary porosity evaluation methodology

Carbonate rocks are comprised of grains and skeletal fragments at depositional conditions. At early stages of digenesis with
minimum or no cementation on the grain-contact, carbonate rock still can be treated as a granular medium. Ongoing diagenesis
causes increasing cementation on grain-contacts until pore-isolation or total lithification occurs. At the same time, dissolution
along the pore channels and intra-grain boundaries create either connected or disconnected vugs. Presence of shear or tensile
stress on these rocks may create system of fractures, which then, could be enhanced by dissolutions as well. All the processes
described above can be seen in the following core thin-sections:

2.a 2.b 2.c

2.d 2.e 2.f

Fig. 2 Thin Sections

Fig 2a represents granular carbonate with no or insignificant cementation on the contacts; Fig 2c shows both grains and
skeletal clasts with well-developed intra-skeletal moldic porosity. Fractures can be seen at both micro scale (Fig 2e-2f) and
macro-scale (Fig 1c).

We used Dvorkin-Nur16 model to describe sonic-wave behavior in our un-lithified carbonate clasts. Effective media in this
model is represented by two end-points. First end-point is zero-porosity and represented by bulk and shear modules of the
minerals. Second end-point is the media at the critical/depositional porosity conditions and it is represented by Hertz 27-
Midlin17 theory for well-sorted rocks of the reservoir:
1
2
5 − υ ⎡ 3 ⋅ n 2 ⋅ (1 − φc ) µ 2 ⎤ 3
µhm = ⋅ ⎢ ⋅ Peff ⎥ …… 1
5 ⋅ (2 − υ ) ⎣ 2 ⋅ π 2 (1 − υ )2 ⎦
4 SPE 166037

1
⎡ n 2 ⋅ (1 − φc )2 µ 2 ⎤ 3
K hm = ⎢ 2
P eff ⎥ …… 2
⎣ 18 ⋅ π (1 − υ )
2
⎦

where Khm and µhm are the dry rock bulk and shear modulus at critical porosity (φc), ν is the Poison’s ratio, Peff is the effective
stress and is the difference between overburden stress and pore pressure. While “n” is the coordination number, by Mavko et
al., 199818, characterizing average number of contacts per grain and given by Murphy correlation10:

n = 20 − 34 ⋅ φc + 14 ⋅ φc2 …… 3

Interpolation between these two end points is done using “Hashin-Shtrikman bounds11” that defines the upper and lower
bounds for the elastic moduli of any material reveals the following for the lower bound:

−1
⎡ φ φ ⎤
⎢ φ 1 − ⎥
φc
µdry = ⎢ c
+ ⎥ − z ……4
⎢ µ hm + z µ + z ⎥
⎢⎣ ⎥⎦

µ hm ⎛ 9 ⋅ K hm + 8 ⋅ µ hm ⎞
z= ⋅ ⎜ ⎟ …… 5
6 ⎜⎝ K hm + 2 ⋅ µ hm ⎟⎠

−1
⎡ φ φ ⎤
⎢ 1− ⎥
φc φc 4 …… 6
K dry = ⎢ + ⎥ − ⋅ µ hm
⎢ K + 4 ⋅ µ hm K + 4 ⋅ µ hm ⎥ 3
⎢ hm 3 3 ⎥
⎣ ⎦
Saturated rocks bulk modulus derived from the Gassman12 equation:

K sat K dry Kf
= + …… 7
K − K sat K − K dry φ ⋅ (K − K f )

We assume that non-lithified limestones with intergranular porosity follows Dvorkin-Nur predicted trend. Stiff-end-point is
represented by bulk and shear moduli of calcite mineral (or a mix of calcite and dolomite) at zero porosity. While soft-end-
point is represented by moduli computed from “Hertz-Midlin contact theory” at depositional/critical porosity (approximately
40%). However lithified vuggy or fractured rocks are expected to deviate from this trend.

Assuming that molds and isolated pores are voids in the bulk body of rocks, the best representation of vugs will be spherical
inclusions into the bulk body of solid rocks filled by pore-fluids. Similarly, we believe that the fractures can be approximated
as “prolonged-inclusions filled with pore-fluids”. There are various models that can describe wave scattering on the different
shape inclusions namely “Self-consistent approximation” (Budiansky 1965)20, “Differential Effective medium” (Claery21
1980; Norris22 1985; Zimmerman23 1991), “Hudson’s cracked media” (Hudson24 1980), “Eshelby-Cheng cracked anisotropic
media model” (Cheng25,26 1987, 1993). In this study, we have selected Kuster and Toksoz’s “Long Wavelength First Order
Scattering Theory” (Kuster and Toksoz 1974) for its simplicity. The model describes wave scattering on the inclusions of
different shapes, in particular, spheres, needles, disks and penny-cracks. Spherical inclusions filled with pore-fluids area close
approximation for the vugs, while penny-crack type is a good approximation for the fractures perpendicular to the borehole in
horizontal direction.
SPE166037 5

General form of the Kuster and Toksoz’ formula for the different shapes of inclusions can be written as:

(K kt ) ((KK
− Km ⋅ m + µm ) n
= ∑ X i ⋅ (Ki − K m ) ⋅ Pmi
+ µm ) i =1
…… 8
kt

(µkt − µm )⋅ (µm + ε m ) = ∑ X i ⋅ (µi − µm )⋅ Qm

n
…… 9
(µkt + ε m ) i =1 i

Where Kkt and µkt are the bulk and shear moduli of the effective media with inclusions; Km and µm are the bulk and shear
moduli of the host media that, in our case, is carbonate. Xi is the volume occupied by inclusions and Ki and µi are the bulk and
shear moduli of the inclusions which in our case are fluids, thus µI is zero. Shear modulus (εm) is given by the equation 10.

εm =
µm
∗
(9 ⋅ K m + 8 ⋅ µm )
…… 10
6 (K m + 2 ⋅ µm )
Pmi and Qmi are the coefficients depending on the shape of inclusions. For the spherical inclusions coefficients can be written
as:
4
K m + ⋅ µm
Pmi = 3
4 …… 11
Ki + µm
3

µm + ε m
Qmi = …… 12
µi + ε m

While for the “penny-crack shape”:

4
K m + ⋅ µi
Pmi = 3 …… 13
4
K i + µi + π ⋅ α ⋅ β
3

⎛ 2 ⎞
Ki + ⋅ (µi + µm ) ⎟
1 ⎜ 8 ⋅ µm 3
Qmi = ⋅ ⎜1 + + 2⋅ ⎟ …… 14
5 ⎜ 4 ⋅ µi + π ⋅ α ⋅ (µm + 2 ⋅ β ) 4
⎜ Ki + ⋅ µi + π ⋅ α ⋅ β ⎟⎟
⎝ 3 ⎠

Where:
3 ⋅ K m + µm
β = µm …… 15
3 ⋅ K m + 4 ⋅ µm

And α is the fracture aspect ratio (width to length). Remembering that shear modulus (µi) in our case equals to zero, therefore,
the equations can be simplified further. We assume that carbonates with intergranular porosity will follow the trend predicted
by Dvorkin-Nur model, while carbonates with vuggy or fracture porosity (in the absence of intergranular porosity) would
follow either “spherical inclusions” or “penny-crack” Kuster-Toksoz predicted trends, respectively.

Due to the absence of the fluid effects on shear-modulus, vuggy and fracture-porosity estimations can be made by using shear-
modulus. Rocks are stiffer for vug-dominated systems than intergranular systems; and softer for fracture dominated (Fig. 3).

In our study, we use prediction of shear-modulus for different systems (intergranular, vuggy and fractured) as a basis for the
secondary porosity quantification. Intervals with shear modules higher than predicted by the Dvorkin-Nur model are identified
as “vuggy” and with shear modules lower than predicted labeled as “fractured”. However, we are cognizant of the fact that all
6 SPE 166037

three types of pores may co-exist simultaneously in some cases. Therefore, for any given volume of rocks, total porosity is
made of three components - intergranular, isolated vugs and fracture porosity as formulated below:

φt = φint + φvug + φ frc …… 16

In the first case with our model simplifications, we have only vugs and intergranular
pores, thus intergranular porosity can be estimated as the difference between total and
vuggy porosity:

φint = φt − φvug …… 17

And for the intergranular-fractured system, intergranular porosity can be estimated as

the difference between total and fractured porosity:
3.a 3.b φint = φt − φ frc
Fig. 3 Vugs dominated (3.a), and …… 18
fracture dominated (3.b) systems

We identify a “vuggy-intergranular” system when

shear-modulus is stiffer than what Dvorkin-Nur
model predicts (Fig. 4). While an “intergranular-
fractured” system is detected if it has lower shear
modulus than the one predicted by the Dvorkin-
Nur model. In order to estimate vuggy porosity, we
used the same assumptions as the Dvorkin-Nur
model. The soft member is represented by the
Hertz-Midlin predicted modulus for the
depositional porosity, while the stiff member is
represented by what Kuster-Toksoz’ model
predicts. Interpolation between these two is carried
out with “Hashin-Shtrikman lower bound” using
intergranular porosity as the input.

To compute the value of vuggy porosity, we start

the iterative calculation with zero vug- porosity.
Initially, the model follows “Dvorkin-Nur” trend
and predicts the shear modulus that is lower than
the one computed from the sonic and density logs.
Iteratively increasing vug- porosity causes the
calculated shear-modulus to shift toward to the
Kuster-Toksoz model of sphere predictions that
produces an increasing-trend in shear modulus.
Iterative process continues until a convergence
between calculated and sonic-density driven shear
moduli has been reached. Then, the value of vug-
porosity in the model with minimum disparity
between observed and modeled shear moduli is
reported as the final percentage of the vug-
occupied pore space.
Fig. 4 Proposed workflow

Similarly for the intergranular-fractured system, stiff-end is represented by the Kuster-Toksoz predicted shear modulus for the
penny-crack inclusions, while soft-end represented by Hertz-Midlin estimation. To calculate fracture-porosity, we start with
zero fracture porosity as well. Initially, predicted shear modulus is higher than the one calculated from the logs. By increasing
fracture-porosity percentage, predicted shear modulus starts shifting toward the predictions by the Kuster-Toksoz’s crack
model with a decreasing trend in shear modulus. Iterative calculations lead to convergence between model-predicted versus
log-driven shear moduli, and result in the final value of shear modulus for deriving fracture-porosity. For both cases,
intergranular porosity estimated is the difference between total and either fractured or vuggy porosity.
SPE166037 7

We should note that the “penny-crack

inclusions” model is extremely sensitive to
the aspect ratio and change of the aspect ratio
from 0.01 to 0.002 will lead to fracture-
porosity decrease from 3.5% to 0.5% for the
same inclusion (Fig. 5 a-b). Thus, secondary
porosity estimation for the “intergranular-
fractured” segment of the model is rather
qualitative then quantitative, and requires
careful calibration. The only control on the
aspect ratio value is that in low-porosity-
rocks, estimated fracture porosity cannot
exceed total porosity value. When estimated
porosity turns up exceeding the total porosity,
it indicates that the aspect ratio is
overestimated.
5.a
Prior knowledge of fracture-orientation or
lack of it drives the accuracy/inaccuracy of
the model respectively. The “Penny-crack
inclusions” model assumes that fracture-
planes are orthogonal to the borehole. This
assumption holds for the studied cases based
on available Image log data. However it may
not be correct for the high-angle fractures. In
such cases, “T-matrix inclusion” model15 may
be applied instead of the “Kuster-Toksoz”
model.

Permeability modeling
Total porosity components (intergranular,
fractures and vuggy) computed as discussed
above are used in the “Chen-Jacoby” model
equation (Eqn. 19) to compute permeability.
5.b Outcome of the “Chen-Jacoby” model is a
Fig 5. Fracture porosity 3.5% (5.a) and 0.5% (5.b) showing the effect of changing permeability profile that integrates the matrix
aspect ratios
and connected-vug permeability as a function
of the “level of connectedness” of the vugs. Connectivity factor (p) is derived from acoustic and NMR log data as explained in
publication by Kumar et al., 200919. For verification and integration purposes, total porosity stochastically derived from the
GR, Neutron, Density, PE logs is also compared with NMR-derived total porosity prior to applying the “Chen-Jacobi”
equation.

m n
⎛ φ ⎞ ⎛ p ⋅ BVM ⎞
k , mD = ⎜ ⎟ ⋅ ⎜⎜ ⎟⎟ …… 19
⎝ C ⎠ ⎝ BVI + (1 − p ) ⋅ BVM ⎠

Where: BVM and BVI are bulk volume of movable and irreducible respectively; C, m and n are constant and exponents. P-
factor (p=1-secpor/BVM), density (frequency of fractures) and fracture aperture are computed from image logs. As many
different model equations can be used for fracture permeability, we opted to use the Equation 20 for fracture permeability
estimation.
k frc , mD = 833 ⋅105 ⋅ φ frc ⋅ w2frc
…… 20
Where: wfrc is the average fracture aperture and φfrc is the fracture porosity estimation either from rock-physics model or form
image logs or as follows: φ frc = 0.001 ⋅ w frc ⋅ D frc ⋅ k fi
Where wfrc is the fracture width in mm, Dfrc is number of fractures per m; and kfi is the number of fracture directions.

There are various practices for averaging permeability profiles that are originating from matrix, fracture networks and
connected vugs of dual-porosity systems. In our approach to calculating average permeability for a given interval so that it can
8 SPE 166037

be compared against the DST/well test permeability, we presumed that flow would ultimately be controlled by the
predominant permeability, thus, final permeability should be estimated as the maximum of intergranular-vuggy and/or fracture
systems permeability values for the given discrete-point in depth when all co-exist. The literature survey shows that the most
common averaging method for permeability in dual porosity systems is using geometric mean of the fracture and vuggy-
intergranular matrix permeability values. For the sake of completeness, we have also computed geometric averages of final
permeability profile of the case studies for the selected intervals. In our case studies, the estimations of final average
permeability for the DST intervals were made by using Power Law, Most Likely Estimate (MLE), arithmetic and geometric
averaging techniques for comparison.

Prediction of permeability in fractured and vuggy intervals

The Middle East case is represented by
dolomitized limestone (from 0% up to 60% of
dolomite) reservoirs that are grainstone to
packstone with some wackestone. Rock-physics
model shows (Fig. 6) vuggy porosity
domination (2-15% with 7-10% in average in
main reservoir intervals) with insignificant
presence of fractures. Intervals of predicted
vuggy and fractured porosity are in good
agreement with data from whole core and
Image log. Core sections demonstrated
significant developments of moldic porosity. No
fractures were observed on core samples, but
Image log showed presence of partly healed
fractures (Fig. 1f). The most significant data
Fig 6. Middle East case, model calibration available for this case are bulk and shear
modulus measured on the core samples. Vuggy
porosity was estimated from the core plugs by visual count on the plug surface of the vug occupied area versus total area of the
plug. This count does not include vugs that are in the body of the core plug; hence invisible to the eye. Comparison of visible
and the rock-physics estimates of the vug-porosity showed that visible vug-porosity from core plugs may under call the actual
vug-porosity by 2-3% (Fig. 7). However general shape of the trend produced by the rock-physics model is in agreement with
the trends revealed by the core measurements with 2-3% corrections (Fig. 7).

The South East (SE) Asian case is represented

by Madura basin carbonates. Core XRD
showed predominantly calcite limestones with
insignificant percent of dolomite (less than
10%). Reservoir sections are comprised of
packstone, wackestone and mudstone.
Secondary porosity development for
wackestones is crucial for the reservoir
producibility because the lack of it may lower
the hydrocarbon storage to uneconomical
levels. The available thin sections demonstrate
whole spectrum of reservoirs: intergranular,
intergranular with molds, vugs dominated
reservoirs, intergranular dominated with micro
fractures (Fig. 2a-f). The fracture networks and
vugs are observed at the whole-core scale and
on the Images (Fig. 1a-c). Wide spectrum of
Fig 7. Middle east case model against core calibration the reservoir types makes it a perfect case for
testing the rock-physics model. Rock-physics evaluation have resulted in 16% of vuggy porosity (6-12% in average) and up to
4.5% fracture porosity with predominant range of 0.5-2.5% in the fractured intervals. The vuggy and fractured intervals have
been detected through the rock-physics model are in good agreement with the Image logs-driven trends (Fig. 8, Fig. 9).

It should be noted that the fracture porosity driven from the rock physics model is strongly dependent on the “fracture aspect
ratio”. As demonstrated in the case study, fracture porosity can vary in the range of 0.5% to 4.5% (averaging at 3.5%) to 0.1%
to1.5 averaging at 0.5%) with a five-fold decrease in the aspect ratio (from 0.01 to 0.002). (Fig. 5)
SPE166037 9

In co-existence of fractures and vugs in the

same interval, the rock-physics model will
detect only predominant porosity of these
two, in terms of combined stiffness trends
from fracture and vugs. The estimated
predominant porosity would be lower than
insitu porosity due to superposition of the
soft and stiff pore-inclusions, namely
fractures and vugs (Fig. 10). Fig. 11 shows
comparison of the vug-porosity estimated
from the rock-physics model and secondary
porosity index from P-sonic. However, P-
sonic-driven secondary porosity index
underestimates vug-porosity in most of the
intervals due to combined effect of the fluid
and the stiffness of the connected grains.

Permeability averaging and flow

prediction South East Asian case:
Fig 8. South East Asia, vuggy interval detected by the rock physics template Three sources of data were available for the
calibration of wireline-driven permeability.
These are: permeability values measured on
the core plugs, the whole core probe-
permeability and DST permeability driven
from the mobility. DST interpretation is
beyond the scope of this paper. We would
like to use the DST-driven permeability for
the subject interval as the “ground truth”
for validating both permeability estimates
(from connectivity-model that requires
porosity components) and permeability
averaging techniques used.

Core plug permeabilities agree with probe-

permeabilities at the lowest scale of
comparisons (Fig. 12.a). Minor disparities
among these two can be explained by the
scale of measurements as the depth of
investigation (DOI) of a probe
permeameter being different than the plug
exposed to the flow in a core holder. As
Fig 9. South East Asia, fractured interval detected by the rock physics template well as, vugs and micro-fractures affecting
probe measurements may not be a part of
the plug volume. At this juncture, we should also note that none of the macro-scale fractures seen in the Image logs can be a
part of the minute volumes of rocks (plugs) tested in a core holder.

Combined intergranular and vuggy matrix permeability has been calculated by “Chen-Jacobi” method as discussed earlier. The
porosity components of the Equation 19 were obtained directly from NMR log as well as from the rock physics model.
Fracture permeability that requires the fracture-aperture and porosity is computed by using Image log driven aperture and
fracture porosity from rock physics model and Image logs. Final permeability profile was treated with different averaging
techniques. The origins of the employed permeability models and averaging techniques used are listed below:
1. Geometric mean of the fracture and combined vuggy-intergranular matrix permeability (both fracture and vug-porosity
estimated from the rock-physics model)
2. Dominant permeability estimates of the fracture and combined vuggy-intergranular matrix permeability (both fracture
and vuggy porosity estimated from the rock-physics model)
10 SPE 166037

3. Dominant permeability estimates

of the fracture and combined vuggy-
intergranular matrix permeability (using
Image estimated fracture porosity for the
fracture permeability and P-sonic secondary
porosity index for the vuggy-intergranular
matrix permeability)
4. Dominant permeability estimates
of the fracture and combined vuggy-
intergranular matrix permeability (using
Image estimated fracture porosity for the
fracture permeability and rock-physics
vuggy porosity for the vuggy-intergranular
Fig 10. South East Asia, coexistence of fractures and vugs in the same interval
under call dominant secondary porosity from the rock-physics model matrix).

Comparison of the permeability profiles shows that in well-defined fractured intervals (with no vugs), permeability estimated
from rock-physics agrees with permeability constructed from Image log-driven fracture porosity (Fig. 12b).
Rock-physics estimates only vuggy-porosity in co-presence of fractures and vugs, thus, under-estimating permeability in
comparison to Image log driven one (Fig. 12c). We used auto-detection of fractures from Image logs. Therefore, in the
intervals where fractures are hardly seen on the image (either due to fracture-discontinuity or image quality), rock-physics
template functions better for detecting fractures than various algorithms applied in auto detection routines (Fig. 12d).

Permeability profiles derived as described

above are averaged over DST interval using
Arithmetic, Geometric, Power Law and
MLE averaging techniques. In comparison
with DST permeability, Geometric mean of
fracture combined with connected vug and
intergranular permeability does not reflect
the volumetric scale of the DST; hence
cannot be used as the averaged permeability
for carbonates. Arithmetic mean of
aforementioned permeability can produce
acceptable results for carbonates with
minimal fractures and vugs.

Fig 11. South East Asia, comparison of the rock-physics estimated vuggy porosity
However, the closest permeability average
with secondary porosity index from P-sonic from contributing permeability profiles can
be obtained by using Power Law and/or
MLE averaging techniques (Fig.14). Power
Law will be the closest for the cases when
permeability values vary in orders of
magnitudes as in South East Asian case,
where both fractures with up to few
1000mD permeability and vuggy intervals
with 0.05-100mD are seen over the DST
interval.

The histogram of permeability distribution

with “long tail” fits better to the Power Law
model since the Gaussian distribution of
well-behaving carbonates are replaced with
Paretian behaving distributions due to the
12.a Comparison of geometrical mean of permeability against core plug and probe dominating contributions from fracture
permeabilites networks and connected vugs. MLE
estimate will stay close to the DST results for the cases of disappearing “long tail” of permeability distribution i.e. more
homogenous than heavily-fractured and vuggy, but more complex than well-behaving –not fractured and less vuggy-
carbonates.
SPE166037 11

However, the closest permeability average from contributing permeability profiles can be obtained by using Power Law and/or
MLE averaging techniques (Fig.14) for the case study that is heavily-fractured and vuggy carbonates.

Power Law will be the closest in case when

permeability over averaging interval varies
a few orders of magnitudes. This suits well
to our South East Asian case, where both
fractures with up to 1000mD permeability
and vuggy intervals with 0.05-100mD are
seen over the DST interval.

Rate comparisons and discussion:

The development of productivity
12.b Fractures are easily identified on image log and the rock physics template
estimations and comparison of the results
against DST, provide an assessment of the
most appropriate methodology that can
produce the best porosity, permeability and
saturation and averaging technique for a
specific carbonate type.

We made certain assumptions for the data

that were not available at the time (drainage
radius, viscosity, skin factor) and used the
12.c The rock physics identifies vugs only
drawdown pressure information from DST
in productivity calculations for the
generation of PI (Productivity Index) and
rates. Calibrated and integrated
permeability-profile averages from various
techniques (Power Law, MLE, Geometric
Mean, and Arithmetic Mean) are used
together with petrophysically driven
12.d Fractures invisible to image log are identified by the rock physics saturation data to get the relative
permeability via Corey-Brooks, Pirson and
Koederitz correlations for water-wet
carbonate formations. The relative
permeability at computed saturation and
integrated permeability (closest
approximation to absolute permeability)
produced the effective permeability average
of the subject interval.

The treatment of Darcy’s radial flow

equation with a conceptual tank-model
generated the first estimates of the
productivity index (PI) and the rates. The
gross simplification of conceptual tank-
model can be overcome by computing the PI
for each well by simply creating individual
tanks in the drainage area of the subject well
with the help of the geological model for the
spatial dimensions of the subject carbonate.

This effort can be furthered with more robust

reservoir simulation practices by permitting
flow from each tank to the others in a single
row leading an increase from zero-
Fig 13. South East Asia comparison of the permeability dimensional tank model to one –dimensional
case.
12 SPE 166037

The reservoir simulation practices can also

allow the horizontal flow from one tank to
another at all common faces, then by stacking
these tanks on top of each other and adding
vertical flow element incrementing the
dimensionality of the simulated result (flow
dimensions) from zero to two respectively.

Conclusions and discussion

Our study showed that applied rock-physics
models enable the estimation of secondary
porosity originating from vugs and remain
more on the qualitative side for fractures in
carbonate formation evaluations. Elaborate
evaluation of secondary porosity in
carbonates requires a special suite of wireline
logs comprised of P and S sonic with
Stoneley amplitudes, NMR, Image logs, and
Fig 14 Rate comparisons against DST from different approaches
in some cases, dielectric logs. This
requirement is a must, especially for the exploration wells.

Based on our case studies, we have observed the followings as common traits in required assumptions and also limitations for
the applied rock physics models:

Kuster-Toksoz model is extremely sensitive to “aspect ratio”

Applied model inherits the assumption of “orthogonal fractures to the wellbore”; therefore, the results are only
accurate when the assumption is met. The model will have to be replaced with other models when image log indicate
the presence of “non-orthogonal to borehole” fracture networks

Tested rock physics models cannot address the superimposed vugs and fractured cases leading to underestimation of
porosity; hence, underestimation of permeability

Agreement between wireline and rock-physics-driven porosity, and therefore, permeability-average and the WFT/DST/well
test permeability is the highest when “Power Law model” or “MLE” is used as the averaging technique. This conclusion is
valid especially for carbonates of fractured and vuggy nature with varying degree of connectedness. However, wireline and
rock physics-driven permeability profiles need to be derived and calibrated as discussed and proposed in this study for
lessening the disparity between them.

It has also been observed that geometric mean of the integrated permeability profile (matrix, connected vugs and fracture
permeabilities) of a carbonate formation may agree with whole core and core plug permeability averages while failing to
match with WFT/DST/well test permeability values as heterogeneity in spatial form would affect the WFT/DST/well test
permeability and not being present in the core plugs or the whole core.

Arithmetic mean of the integrated permeability from wireline may be at parity with WFT/DST/well test permeability for the
same interval provided that the carbonate formation neither has fractures nor any connected vugs.

Geometric mean of the wireline-driven combined permeability is at parity with probe, and plug permeabilities when the
carbonate formation neither has fractures nor any connected vugs while failing to match with WFT/DST/well test
permeability.

ACKNOWLEDGEMENTS
We would like to express our gratitude to the managements of PETRONAS & PETRONAS Carigali and LUKOIL Overseas
for support in publication of this study.
SPE166037 13

NOMENCLATURE

REFERENCES
1. M.M. Altunbay, A H A Ramli, R.E. Poit, ”Approaching a Parity Between Well Test and Wireline-driven Average Permeability for
Carbonates”, SPE-166148-MS, SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, Sep-2 Oct 2013
2. S. Chen, D. Jacobi, H. Kwak, M. Altunbay, J. Kloos, ”Pore-Connectivity based Permeability Model for Complex Carbonate
Formations”,SPWLA 2008-E,49th Annual Logging Symposium, May 25 - 28, 2008 , Austin, Texas
14 SPE 166037

3. Kuster, G.T. and Toksoz1, M.N., 1974, “Velocity and attenuation of seismic waves in two phase media” Geophys., 39, 587-618
4. Dvorkin, J. and Nur2, A., 1996, “Elasticity of high-porosity sandstones: theory for two North Sea datasets” Geophys., 61, 1363-
1370
5. Lucia, J., Jennings, J.,Rahnis, M., and Meyer, F., 2001,”Permeability and Rock Fabric from wireline Logs, Arab-D Reservoir,
Ghawar Field, Saudi Arabia”, GeoArabia, v.6, no 4
6. Chilingarian, G.V., Mazzullo, S.J., and Rieke, H.H., 1996,”Carbonate Reservoir Characterization A Geological Engineering
Analysis, Part II.”, Elsevier, Amsterdam, The Netherlands, 994 pp
7. Chang, D., Vinegar, H., Morris, C., and Straley, C. 1997, “Effective porosity, producible fluid, and permeability in carbonates
from NMR logging”, Log Analyst 38, 60-72
8. Deutsch, C., “Calculating effective absolute permeability in sandstone/shale sequences,”, SPE Formation Evaluation, September,
1989
9. J.P. Hoogenboom, W.K. den Otter, H.L. Offerhaus, ”Accurate and unbiased estimation of power-law exponents from single-
emitter blinking data”, The journal of Chemical Physics 125, 204713 (2006)
10. Murphy, W.F., 1982,”Effect of microstructure and pore fluids on the acoustic properties of granular sedimentary materials”, PhD
dissertation, Stanford University
11. Hashin, Z., and Shtrikman3, S., 1963, “A variational approach to the elastic behavior of multiphase minerals”, J. Mech. Phys.
Solids, 11, 127-140.
12. Gassmann, F., 1951, Über die Elastizität poröser Medien: Vier. Der Natur. Gesellschaft in Zurich, 96, 1–23
13. Prasad, M. and Dvorkin, J., “Velocity to porosity transformation in Marine Sediments” Petrophysics, Vol42, No 5, P 429-437
14. Raymer, L.L., Hunt, E.R. and Gardner, J.S., 1980, “An improved sonic transit time-to-porosity transform” Well Log Analyst 21st
Annual Logging Symposium, Paper P
15. Jacobsen, M., Hudson, J.A., and Johansen, T.A., 2003, “T-matrix approach to shale acoustic” Geophysic. J. Int., 154, 533-558
16. Dvorkin, J., and Nur, A., 1996, “Elasticity of high-porosity sandstones, Theory for two Sea data sets”, Geophysics, 61, p 559-564
17. Mindlin, R.D., 1949, “Compliance of Elastic Bodies in Contact”4, J. Appl. Mech., 16, p 259-268
18. Mavko, G., Mukerji, T., 1998, The Rock Physics Handbook: tools for seismic analysis in porous media, Cambridge University
Press, 329.
19. Kumar, M., M. Altunbay, Arina, W., “Permeability Estimation in Vuggy Carbonates”, PGCE, 2009, Kuala Lumpur, Malaysia
20. Budiansky, B., 1965, J. Mech. Phys. Solids 13, 223
21. Cleary, M.P., Chen, I. W., and Lee, S.M., 1980, “Self-consistent Techniques for Heterogeneous Solids,” ASCE Journal of
Engineering Mechanics, Vol. 106, pp 861-887
22. Norris, A.N., 1985, “A Differential Scheme for the Effective Moduli of Composites,”, Mechanics of Materials, Vol. 4, pp. 1-16
23. Zimmerman R.W., 1991, Compressibility of Sandstones, Elsevier Science Publishing Co.
24. Hudson, J. (19080), Overall Properties of Cracked Solid., Math. Proc. Cambridge Philos. Soc. 88, 371-384
25. Preprint. Article published as: Detournay, F., Cheng, A.H.D., “Fundamentals of poroelasticity”, Chapter 5 in Comprehensive Rock
Engineering: Principles, Practice and Projects, Vol. II, Analysis and Design Method, ed. C. Fairhurst, Pergamon Press, pp 113-
171, 1993
26. Eshelby, J.D., 1957, The Determination of the elastic field of an ellipsoidal inclusion and related problems, Proc. R. Soc. Lond. A,
241, 376-396
27. H. Hertz, Über die berührungfester elastischer Körper (On the contact of rigid elastic solids). In: Miscellaneous Papers. Jones and
Schott, Editors, J. reine und angewandte Mathematik 92, Macmillan, London (1896), p. 156 English translation: Hertz, H.
28. Hertz, H. R., 1882, Ueber die Beruehrung elastischer Koerper (On Contact between Elastic Bodies), in GesammelteWerke
(Collected Works), Vol. 1, Leipzig, Germany, 1895.

1
Kuster-Toksoz method calculates the effective moduli for randomly distributed inclusions based on a long-wavelength, first-order scattering theory.
2
Dvorkin-Nur model predicts bulk and shear moduli behavior versus porosity in a granular media. At the soft-end, it uses Hertz-Mindlin contact theory to
predict moduli behavior in un-cemented grain conditions at depositional porosity with confine pressure applied to the grains. Confine pressure in this case
assumed to be difference between overburden and pore- pressures. The stiff end is represented by mineral bulk and shear moduli. Interpolation between these
two points is done through the lower Hashin-Shtrikman bound. Hashin-Shtrikman bounds represent mixing of two different moduli (one stiff and one soft) in
assumption of spheres where material of one stiffness coating another.

3
Methodology for computing the bounds on the elastic moduli and tensors of transversally isotropic and isotropic composites.