Angular and Frequency-Dependent Wave Velocity and Attenuation in Fractured Porous Media
JOSE
´
M. CARCIONE,
1
BORIS GUREVICH,
2,3
JUAN E. SANTOS,
4,5,6
and STEFANO PICOTTI
1
Abstract—Wave-induced fluid flow generates a dominant
attenuation mechanism in porous media. It consists of energy loss
due to P-wave conversion to Biot (diffusive) modes at mesoscopic-
scale inhomogeneities. Fractured poroelastic media show signifi-
cant attenuation and velocity dispersion due to this mechanism.
The theory has first been developed for the symmetry axis of the
equivalent transversely isotropic (TI) medium corresponding to a
poroelastic medium containing planar fractures. In this work, we
consider the theory for all propagation angles by obtaining the five
complex and frequency-dependent stiffnesses of the equivalent TI
medium as a function of frequency. We assume that the flow
direction is perpendicular to the layering plane and is independent
of the loading direction. As a consequence, the behaviour of the
medium can be described by a single relaxation function. We first
consider the limiting case of an open (highly permeable) fracture of
negligible thickness. We then compute the associated wave
velocities and quality factors as a function of the propagation
direction (phase and ray angles) and frequency. The location of the
relaxation peak depends on the distance between fractures (the
mesoscopic distance), viscosity, permeability and fractures com-
pliances. The flow induced by wave propagation affects the quasi-
shear (qS) wave with levels of attenuation similar to those of the
quasi-compressional (qP) wave. On the other hand, a general
fracture can be modeled as a sequence of poroelastic layers, where
one of the layers is very thin. Modeling fractures of different
thickness filled with CO
2
embedded in a background medium
saturated with a stiffer fluid also shows considerable attenuation
and velocity dispersion. If the fracture and background frames are
the same, the equivalent medium is isotropic, but strong wave
anisotropy occurs in the case of a frameless and highly permeable
fracture material, for instance a suspension of solid particles in the
fluid.
Key words: Fractures, anisotropy, attenuation,
boundary conditions.
1. Introduction
The acoustic characterization of fractures and
cracks is important from the point of view of reservoir
development (e.g., G
UREVICH et al., 2009). Moreover, it
is also important in CO
2
storage to monitor the injected
plumes as faults and fractures are generated, where the
gas can leak to the surface (P
ICOTTI et al., 2012). A
dense set of planar fractures behaves as an effective
long-wavelength TI medium, leading to azimuthally
varying velocity and attenuation of seismic waves. One
of the important mechanisms of seismic attenuation in
porous media is wave-induced fluid flow, by which the
fast P wave is converted to slow (diffusive) modes of the
Biot type (W
HITE et al., 1975;PRIDE et al., 2004;CAR-
CIONE
and PICOTTI; 2006). The phenomenon has been
studied for alternating thin poroelastic layers, along the
direction perpendicular to the layer planes (W
HITE
et al., 1975) at all frequencies and propagation angles
(K
RZIKALLA and MU
¨
LLER, 2011;CARCIONE et al., 2011).
The case of planar fractures embedded in a poroelastic
background medium is a particular case of the thin layer
problem, where one of layers becomes extremely thin
and compliant. Alternative models, based on distribu-
tions of finite cracks, are given in C
HAPMAN (2003)and
G
ALVIN and GUREVICH (2009). A comprehensive review
of the wave-induced fluid-flow loss mechanisms, with a
discussion of the related relaxation frequencies, is given
in the review paper by M
U
¨
LLER et al.(2010). Readers
can also refer to M
AVKO et al.(2004) for detailed
mathematical expressions of the basic loss mechanisms.
In the case of a solid (non-porous) background,
C
HICHININA et al.(2009) and CARCIONE et al.(2012)
obtained analytical solutions for the TI and more
general anisotropic cases, respectively, i.e., the
1
Istituto Nazionale di Oceanografia e di Geofisica Speri-
mentale (OGS), Borgo Grotta Gigante 42c, 34010 Sgonico, Trieste,
Italy. E-mail: jcarcione@inogs.it
2
Department of Explorations Geophysics, Curtin University,
GPO Box U1987, Perth, WA 6845, Australia. E-mail:
B.Gurevich@curtin.edu.au
3
CSIRO Earth Science and Resource Engineering, ARRC,
26 Dick Perry Avenue, Kensington, Perth, WA 6151, Australia.
4
CONICET, Instituto del Gas y del Petro
´
leo, Facultad de
Ingenierı
´
a, Universidad de Buenos Aires, C1127AAR Buenos
Aires, Argentina. E-mail: santos@math.purdue.edu
5
Department of Mathematics, Purdue University, West
Lafayette, IN 47907-2067, USA.
6
Universidad Nacional de La Plata, La Plata, Argentina.
Pure Appl. Geophys. 170 (2013), 1673–1683
Ó 2013 Springer Basel
DOI 10.1007/s00024-012-0636-8
Pure and Applied Geophysics
complex and frequency-dependent stiffness compo-
nents and corresponding wave velocities and quality
factors. Regarding poroelastic media, G
UREVICH
(2003) and BRAJANOVSKI et al.(2005) found the low-
and high-frequency limit elasticities of the equivalent
TI medium under the long-wavelength assumption.
Moreover, they obtained the expression of the P-wave
modulus for waves propagating normal to the fracture
planes as a function of frequency. L
AMBERT et al.
(2006) validated the theory by numerical simulations
using a poroelastic reflectivity algorithm. D
U et al.
(2011) have developed a similar stress–strain relation
for fractured porous media in the low-frequency limit.
To obtain the wave properties as a function of
frequency and propagation (phase and ray) angle, we
need all the stiffness components of the equivalent TI
medium. An heuristic approach has been outlined by
L
AMBERT et al.(2005), who analyzed the qP wave. A
physical interpretation has been given by K
RZIKALLA
and MU
¨
LLER (2011) for layered media, who verified the
correctness of the method by comparison to numerical
solutions. Further tests for layered media have been
performed by C
ARCIONE et al.(2011). Basically, the
approach assumes a 1D character of the fluid pressure
equilibration process which generates diffusive modes
from the fast P wave, i.e., the fluid-flow direction is
perpendicular to the fracture plane. In the presence of
horizontal plane layers, the initial fluid pressure field is
independent of the type of excitation, i.e., P waves
traveling horizontally or vertically (mainly compres-
sion), or S waves (shearing) will excite the fluid
pressure in such a way as to maintain its distribution.
For more general layer geometries, the process may
depend on the direction of wave propagation and these
assumptions could not be strictly valid. As a conse-
quence, the model considers one relaxation function,
corresponding to the symmetry-axis P-wave stiffness.
From this relaxation function and the high- and low-
frequency elastic limits of the stiffness tensor, we
obtain the five stiffnesses of the equivalent medium.
We then obtain the quality factors and wave velocities
as a function of frequency and propagation angle.
The examples consider a fractured-sandstone
model with parameters mainly taken from B
RAJANOV-
SKI
et al.(2005), i.e., the properties of the background
porous medium and the normal and tangential fracture
compliances. Generally, in this work, fractures are
planar, rotationally symmetric, highly permeable and
have infinite extent. The degree of attenuation and
velocity dispersion caused by different types of het-
erogeneities in the rock properties, namely, porosity,
grain and frame moduli, permeability, and fluid
properties were studied in detail in C
ARCIONE and PI-
COTTI
(2006). Here, we focus on the anisotropy
properties of wave velocity and attenuation.
2. Highly Permeable Fracture
Let us consider a set of fractures parallel to the
(x, y)-plane in an isotropic porous medium of porosity
/, separated by an average distance L, and denote by
Z
N
and Z
T
the normal and shear excess fracture com-
pliances. Here, we assume that the fracture thickness is
nil (extremely small) and that the porosity within the
fracture is one, i.e., there is no fracture frame. More-
over, let us denote by K
m
and l the dry-rock bulk and
shear moduli of the background porous medium,
respectively, with E
m
= K
m
? 4l/3 and k
m
= K
m
-
2l/3. Let K
s
and K
f
be the grain and fluid moduli,
respectively. Biot’s effective stress coefficient is
a ¼ 1
K
m
K
s
: ð1Þ
G
UREVICH (2003) and BRAJANOVSKI et al.(2005)
obtained the low-frequency limit (relaxed) elasticities
of the equivalent TI medium,
c
r
11
¼
c
r
13
2
c
r
33
þ 4 ð1 cÞ l þ
a
2
c
2
b
1 þ Z
N
b
;
c
r
13
¼ c
r
33
1 2c þ 2ac
M
E
G
a þ Z
N
E
m
1 þ Z
N
b
;
c
r
33
¼ E
G
1 þ
Z
N
ðaM E
G
Þ
2
E
G
ð1 þ Z
N
bÞ
"#
1
;
c
r
55
¼ðl
1
þ Z
T
Þ
1
;
c
r
66
¼ l;
ð2Þ
where
E
G
¼K
m
þ a
2
M þ
4
3
l; b ¼ ME
m
=E
G
;
c ¼ l=E
m
: ð3Þ
(The first Eq. (35) in B
RAJANOVSKI et al.(2005)
should be (l
b
-1
? Z
T
)
-1
instead of l
b
-1
? Z
T
, using
their notation).
1674 J. M. Carcione et al. Pure Appl. Geophys.
The high-frequency limit elasticities (unrelaxed)
are given by (B
RAJANOVSKI et al., 2005),
c
11
¼ c
33
¼ E
G
;
c
13
¼ c
11
2l;
c
55
¼ c
r
55
;
c
66
¼ l;
ð4Þ
where
M ¼
K
s
1 / K
m
=K
s
þ /K
s
=K
f
: ð5Þ
Following the approach used by K
RZIKALLA and
M
U
¨
LLER (2011) for a layered porous medium (see also
C
ARCIONE et al., 2011), we obtain the five stiffnesses
of the equivalent TI medium as
p
IJ
ðxÞ¼c
IJ
þ
c
IJ
c
r
IJ
c
33
c
r
33
½p
33
ðxÞc
33
; ð6Þ
where x is the angular frequency, and
p
33
ðxÞ¼
1
E
G
þ
ððaM=E
G
Þ1Þ
2
b# cot # þ Z
N
1
"#
1
;
# ¼
ffiffiffiffiffi
x
iD
r
L
2
; D ¼
jb
g
ð7Þ
(B
RAJANOVSKI et al, 2005), where g is the fluid vis-
cosity, j is the permeability and i ¼
ffiffiffiffiffiffiffi
1
p
: (Note an
omission in Eq. (15) of B
RAJANOVSKI et al.(2005):
Using their notation, M
b
/C
b
should be replaced by
L
b
M
b
/C
b
in the denominator. Eq. (8)inGUREVICH
et al.(2009) is correct). We have used the opposite
sign convention to express the properties in the
Fourier domain, i.e., x has been replaced by -x in
Eq. (15)ofB
RAJANOVSKI et al.(2005) and Eq. (8)of
G
UREVICH et al.(2009). The quantity x L
2
/(4D)is
the normalized frequency. Equations (6) and (7) are
deduced from Eqs. (25), (41) and (42) of K
RZIKALLA
and MU
¨
LLER (2011). Equation (6) is obtained by
considering the 1D character of the fluid pressure
equilibration process between the background por-
ous medium and the fractures, assuming that the
fluid-flow direction is perpendicular to the fracture
plane. As a consequence, the model considers one
relaxation function, corresponding to the symmetry-
axis P-wave stiffness. Therefore, knowing this
relaxation function and the high- and low-frequency
elastic limits of the stiffness tensor, one can obtain
the five complex and frequency-dependent stiff-
nesses of the equivalent viscoelastic medium. The
equations hold for frequencies x x
B
= gp /
(jq
f
), where q
f
is the fluid density, so that the fluid
flow in the pores is of Poiseuille type and the
effective medium approximation (Backus averaging)
is valid.
The low-frequency regime occurs when pressure
has enough time to equilibrate between the back-
ground medium and the fractures within a wave
period. This happens when the diffusion length
ffiffiffiffiffiffiffiffiffiffi
D=x
p
(or wavelength of Biot’s slow wave) is much
larger than the period L. High frequency means that
the frequency is much less than x
B
but much larger
than D/L
2
. In both cases, the wavelength of the pulse
must be much larger than the spacing L (e.g., C
AR-
CIONE
, 2007).
The bulk density is given by
q ¼ð1 /Þq
s
þ /q
f
; ð8Þ
where q
s
is the grain density.
3. General Fractures
The equivalent medium corresponding to a set of
highly permeable fractures can be obtained as a limit
case of a layered medium, where one of the layers is
very soft and thin. Equation (7) has been obtained in
this manner by B
RAJANOVSKI et al.(2005). In this
section, we do not consider any approximation, the
possibility of having a different fluid saturating the
fracture, a finite fracture thickness and a poroelastic
medium forming the fracture material. Consider a
stack of two thin alternating porous layers of thick-
ness l
1
and l
2
, such that the period of the stratification
is L = l
1
? l
2
. The complex and frequency depen-
dent stiffness p
33
is given by
p
33
¼
1
E
G
þ
2ðr
2
r
1
Þ
2
ixðl
1
þ l
2
ÞðI
1
þ I
2
Þ
"#
1
; ð9Þ
where
r ¼
aM
E
G
; I ¼
g
ja
coth
aL
2
; a ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffi
ixgE
G
jME
m
r
;
ð10Þ
Vol. 170, (2013) Angular and Frequency-Dependent Wave Velocity 1675
for each single layer (WHITE et al., 1975;NORRIS,
1993;C
ARCIONE and PICOTTI 2006) [see also CARCIONE
(2007), Eq. (7.400)]. If /
2
! 1 and Z
N
¼ lim
l
2
!0
½l
2
=ðLE
m2
Þ; we obtain Eq. (7). A similar limit for the
shear compliance Z
T
is obtained by replacing E
m2
with l
2
. To use Eq. (9) for a highly permeable frac-
ture, one must consider l
2
/L 1, /
2
& 1, E
m2
=
l
2
/(LZ
N
), l
2
= l
2
/(LZ
T
), a
2
& 1 and j
2
1D.
In the general case (an arbitrarily permeable
fracture), we consider l
2
/L 1 and general poro-
elastic parameters. The equivalent medium is
obtained from Eq. (6), where p
33
is given by Eq. (9),
with
c
r
66
¼ B
1
¼hli;
c
r
11
2c
r
66
¼ c
r
12
¼ B
2
¼ 2
k
m
l
E
m
þ
k
m
E
m
2
1
E
m
1
þ
B
6
2
B
8
;
c
r
13
¼ B
3
¼
k
m
E
m
1
E
m
1
þ
B
6
B
7
B
8
;
c
r
33
¼ B
4
¼
1
E
m
1
þ
B
7
2
B
8
¼
1
E
m
a
E
m
2
E
G
ME
m
1
"#
1
;
c
r
55
¼ B
5
¼hl
1
i
1
;
B
6
¼B
8
2
al
E
m
þ
a
E
m
k
m
E
m
1
E
m
1
!
;
B
7
¼B
8
a
E
m
1
E
m
1
;
B
8
¼
1
M
þ
a
2
E
m
a
E
m
2
1
E
m
1
"#
1
; ð11Þ
and
c
66
¼ c
r
66
;
c
11
2c
66
¼ c
12
¼ 2
ðE
G
2lÞl
E
G
þ
E
G
2l
E
G
2
1
E
G
1
;
c
13
¼
E
G
2l
E
G
1
E
G
1
;
c
33
¼
1
E
G
1
;
c
55
¼ c
r
55
ð12Þ
where we kept the notation of G
ELINSKY and SHAPIRO
(1997) for reference (see also CARCIONE et al. 2011),
and where hui¼ðl
1
u
1
þ l
2
u
2
Þ=L:
Appendix shows how to obtain the energy
velocities and dissipation factors as a function of the
ray and phase angles from the complex stiffnesses (6)
and the composite density (8).
4. Example
The examples consider brine- and oil-saturated
sandstones, with grain properties: K
s
= 37 GPa, l
s
=
44 GPa and q
s
= 2,650 kg/m
3
. Brine has the proper-
ties: K
f
= 2.25 GPa, q
f
= 1,040 kg/m
3
, g = 0.0018 cP,
while those of oil are K
f
= 2 GPa, q
f
= 870 kg/m
3
, g =
0.3 cP. The dry-rock bulk and shear moduli are given
by the Krief model,
K
m
K
s
¼
l
l
s
¼ð1 /Þ
3=ð1/Þ
ð13Þ
(K
RIEF et al., 1990;BRAJANOVSKI et al., 2005).
Porosity and permeability are related by an equation
derived by C
ARCIONE et al.(2000),
j ¼
r
2
g
/
3
45ð1 /Þ
2
; ð14Þ
where r
g
=20lm denotes the average radius of the
grains.
We consider two sets of compliances taken from
B
RAJANOVSKI et al.(2005), obtained as
Z
N
¼
1
E
m
1
DN
1
1
and Z
T
¼
1
l
1
DT
1
1
;
where ðDN; DTÞ = (0.2, 0.5), giving 1/Z
N
= 121 GPa
and 1/Z
T
= 14 GPa, and ðDN; DTÞ= (0.02, 0.05), giving
1/Z
N
= 1,483 GPa and 1/Z
T
= 264 GPa. In the second
set, the two faces of a single fracture are in better
contact, since perfect bonding or absence of fractures
occurs when Z
N
! 0 and Z
T
! 0: The quantities DN
and DT are dimensionless fracture weaknesses (see
H
SU and SCHOENBERG, 1993;BAKULIN et al., 2000).
First, we assume / = 0.25, L = 0.2 m and DN =
0.2, and compute the P-wave phase velocity and
dissipation factor along the direction perpendicular to
the fracture plane (see Fig. 1). As can be seen, the
relaxation peak for oil is located at the seismic fre-
quency band, while that of brine is located at the
sonic band. The peak relaxation frequency is pro-
portional to j/g [see Eq. (7.401) in C
ARCIONE (2007)]
1676 J. M. Carcione et al. Pure Appl. Geophys.
and can substantially be affected by the normal
fracture compliance, as shown in Fig. 2, where the P-
wave peak frequency and peak dissipation factor for
oil saturation are represented as a function of DN:
The peak moves from the sonic band for DN = 0.02 (4
kHz and Q = 286) to the seismic band for DN = 0.2
(35 kHz and Q = 26).
Next, we consider oil saturation, DT = 0.5 and use Eq.
(6) to obtain the complex stiffnesses at 35 Hz, where
the peak quality factor has a value of 26. We obtain
p
11
¼ð33:68; 0:001Þ;
p
13
¼ð5:75; 0:034Þ;
p
33
¼ð30:71; 1:17Þ;
p
55
¼ 6:96;
p
66
¼ 13:92
ð15Þ
in GPa. The energy velocities and dissipation factors
as a function of the ray angle w and propagation
(phase) angle h for oil saturation are represented in
Fig. 3a and b, respectively. There is a noticeable
shear-wave splitting, also called seismic birefrin-
gence, and the coupling between the qP and qS waves
generates strong shear attenuation (Q & 40) at
approximately 50°, with magnitudes comparable to
those of the qP wave. The SH wave is lossless.
Energy velocity rather than phase velocity is shown,
since it represents the wavefront (e.g., C
ARCIONE,
2007).
The same plots are repeated for oil saturation with
ðDN; DTÞ = (0.02, 0.05) (Figs. 4, 5). In this case, the
fracture is stiffer than above, i.e., the two faces of the
fracture are in better contact. As can be seen, the peak
(a)
(b)
Figure 1
P-wave phase velocity (a) and dissipation factor (b) along the
direction perpendicular to the fracture planes for brine and oil
filling the pore space with DN = 0.2
(a)
(b)
Figure 2
P-wave peak relaxation frequency (a) and peak dissipation factor
(b) as a function of the normal fracture parameter DN for oil
saturation
Vol. 170, (2013) Angular and Frequency-Dependent Wave Velocity 1677
