# Dual simulations of fluid flow and seismic wave propagation in a fracture network: effect of pore pressure on seismic signature

S. Vlastos, E. Liu, I.G. Main, M. Schoenberg, C. Narteau, Bertrand Maillot

► To cite this version:

S. Vlastos, E. Liu, I.G. Main, M. Schoenberg, C. Narteau, et al. Dual simulations of fluid flow and seismic wave propagation in a fracture network: effect of pore pressure on seismic signature. Geophysical Journal International, 2006, 166, pp.825-838. hal-00146271

HAL Id: hal-00146271

https://hal.science/hal-00146271v1

Submitted on 14 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Dual simulations of fluid flow and seismic wave propagation in a fractured network: effects of pore pressure on seismic signature

S. Vlastos,<fnref n="1" />,<fnref n="2" /> E. Liu,<fnref n="1" />,* I. G. Main,<fnref n="2" /> M. Schoenberg,<fnref n="3" /> C. Narteau,<fnref n="2" /> X. Y. Li<fnref n="1" /> and B. Maillot<fnref n="4" />

Accepted 2006 May 3. Received 2006 May 2; in original form 2004 April 14

## SUMMARY

Fluid flow in the Earth's crust plays an important role in a number of geological processes. In relatively tight rock formations such flow is usually controlled by open macrofractures, with significant implications for ground water flow and hydrocarbon reservoir management. The movement of fluids in the fractured media will result in changes in the pore pressure and consequently will cause changes to the effective stress, traction and elastic properties. The main purpose of this study is to numerically examine the effect of pore pressure changes on seismic wave propagation (i.e. the effects of pore pressures on amplitude, arrival time, frequency content). This is achieved by using dual simulations of fluid flow and seismic propagation in a common 2-D fracture network. Note that the dual simulations are performed separately as the coupled simulations of fluid flow and seismic wave propagations in such fracture network is not possible because the timescales of fluid flow and wave propagation are considerably different (typically, fluid flows in hours, whereas wave propagation in seconds). The flow simulation updates the pore pressure at consecutive time steps, and thus the elastic properties of the rock, for the seismic modelling. In other words, during each time step of the flow simulations, we compute the elastic response corresponding to the pore pressure distribution. The relationship between pore pressure and fractures is linked via an empirical relationship given by Schoenberg and the elastic response of fractures is computed using the equivalent medium theory described by Hudson and Liu. Therefore, we can evaluate the possibility of inferring the changes of fluid properties directly from seismic data. Our results indicate that P waves are not as sensitive to pore pressure changes as S and coda (or scattered) waves. The increase in pore pressure causes a shift of the energy towards lower frequencies, as shown from the spectrum (as a result of scattering attenuation). Another important observation is that the fluid effects on the wavefield vary significantly with the source–receiver direction, that is, the azimuth relative to the fracture orientation. These results have significant implications for the characterization of naturally fractured reservoirs using seismic methods, and may impact on experimental design to infer such attributes in a real reservoir situation, particularly in acquiring time-lapse seismic data.

**Key words:** cellular automaton, effective medium theory, finite difference methods, fractures, scattering attenuation, wave propagation.

## 1 INTRODUCTION

The idea that seismic waves can be used to identify the presence of fluids and the transport properties of rocks takes its root from theoretical studies that go back at least 50 yr. More recently, many examples in time lapse or 4-D seismic surveys have demonstrated that seismic waves can be used to monitor changes in oil or gas

reservoirs as a function of time (Shapiro & Gurevich 2002; Vasco 2004).

During production from a reservoir, the movement of fluids is accompanied by substantial changes in the pore pressure field. As fluid drains, pore pressure in general decreases, increasing the effective pressure on fractures, grain boundaries and microcracks. Higher static load on such surfaces decreases their compliance non-linearly and decreases fracture opening and/or pore throat size, thus increasing the stiffness of the rock (increasing compressional and shear velocities) and decreasing permeability (Schoenberg 2002).

<footnote n="1">¹British Geological Survey, Murchison House, West Mains Road, Edinburgh EH9 3LA, UK. E-mail: E.Liu@bgs.ac.uk</footnote>

<footnote n="2">²School of GeoSciences, University of Edinburgh, Grant Institute, West Mains Road, Edinburgh EH9 3JW, UK</footnote>

<footnote n="3">³Lawrence Berkeley Laboratory, Berkeley, California, USA</footnote>

<footnote n="4">⁴University de Cergy-Pontoise, Cergy-Pontoise Cedex, France</footnote>

<footnote>*Corresponding author.</footnote>
Conversely, pore pressure build-up due to injection leads to a de-
crease in effective pressure and an increase in rock compliance. Also
there is evidence of relation between stress and pore pressure. In the
overpressured parts of the central North graben in the UK, and the
Sable subbasin of the Scotian Self in Canada, data indicate that
the horizontal stress, S<sub>h</sub>, increases at a rate proportional to but less
than the rate of increase of pore pressure P. That is a condition con-
sistent with a pore pressure induced deformation of the rock, called
poroelastic behaviour. That does not happen in normal pressured
basins where friction governs S<sub>h</sub> (Engelder & Fischer 1994).

The complex nature of the crack/pore structure of rocks, and the
behaviour of fluids occupying and flowing within the pore structure,
has been studied extensively by a number of researchers. Parame-
ters of great interest are stress effects, attenuation and dispersion
mechanisms. A review of the effects of those parameters is given by
Winkler & Murphy (1995). The complex microstructures of most
rocks cause velocities and attenuation to decrease. Increasing con-
fining pressure or decreasing pore pressure, cause velocities to in-
crease and attenuation to decrease. Focusing on the pore space,
we can imagine that very compliant pores (such as thin cracks)
will close under small stresses. Increasing stress will close more
and more pores, thereby stiffening the overall frame of the rock.
Several theoretical models (Cheng & Toksöz 1979; Mavko & Nur
1978; O'Connell & Budiansky 1974; Walsh 1965) have utilized this
approach.

When we examine seismic properties, we use the effective stress,
which is the difference between confining pressure and the product
of pore pressure and effective stress coefficient. Wyllie *et al.* (1958)
showed very clearly that, to first order, velocities are a function of
the effective stress on the rock. It is interesting to note that the ve-
locity is independent of the confining pressure when the effective
stress is held constant (by increasing pore pressure at the same rate
as confining pressure). It is generally assumed that attenuation has
a similar dependence on effective stress. Winkler & Nur (1979)
showed that increasing confining pressure, or decreasing pore pres-
sure, both reduce attenuation in water-saturated rock. When stress
is anisotropic, it causes velocities in rock to vary with direction
(Lockner *et al.* 1977; Nur & Simmons 1969). An example of this
effect is shown by Nur & Simmons (1969). In their experiment, a
granite sample was subjected to uniaxial stress, and velocities were
measured as a function of azimuth, defined as the angle between the
ray path and the uniaxial stress direction. At zero stress, the velocity
is virtually independent of azimuth. As stress increases, a strong
anisotropy develops. *P* and *SH* waves are much more sensitive to
stress when propagating parallel to the stress direction than when
propagating perpendicular to the stress direction. It is likely that
stress anisotropy will also create attenuation anisotropy, but there
are very few experimental data are available to confirm this, except
maybe the studies by Best (1994) and Best *et al.* (1994).

Fractured rock is modelled as a relatively rigid, defect free, back-
ground medium in which sets of linear slip interfaces are embedded.
A linear slip interface is a surface across which anomalously large
strain occurs due to the passage of a wave. In linear slip deformation
theory, the large strain is approximated by a displacement discon-
tinuity across the surface which is linearly related to the dynamic
traction acting on the interface (to the first order). The dynamic elas-
tic properties of the rock are then found by adding the compliance
tensor of the background to an excess fracture compliance tensor
associated with the fractures (Schoenberg & Sayers 1995). The lin-
ear parameters governing the infinitesimal slip on these planes have
been shown to depend on the static stress state (Pyrak-Nolte *et al.*
1990; Hsu & Schoenberg 1993; Hudson 2000) in a highly non-linear

and, most likely, in a very hysteretic behaviour under changing static
stress to changes in the dynamic elastic moduli, which are relevant
for the extremely low strain levels encountered in seismic wave
propagation. If such relations are established, we may begin to be
able to predict static effective stress from knowledge of the dynamic
properties (wave speeds and their associated polarizations), and in
particular, to probe the pore pressure field changes induced by reser-
voir drainage or other fluid movement.

Sedimentary rocks are composed of relatively rigid grains which
are pressed together yielding a mass of grain boundaries and cracked
grains, to which must be added larger cracks, fractures, bedding
planes, joints, and faults, which are all what we refer to as compliant
surfaces. These features permeate the rock at many, if not all ori-
entations, and at many scales. Seismic energy propagating through
a rock mass vibrates the compliant parts with such low energy that
the strains across these surfaces, even though orders of magnitude
larger than the strains in the intact grains, are still very small. Seis-
mic response can be thought of as a non-destructive probe into the
linear properties of these weak surfaces. Yet the linear coefficients
which we want to measure seismically are just exactly the param-
eters that depend in so sensitive a fashion on changes in the static
effective stress in the solid phase, and consequently, to changes in
pore pressure.

The movement of fluids within a cracked solid can have a signif-
icant effect on the properties of seismic waves of long wavelength
propagating through the solid. The cracks may be partially or com-
pletely filled with fluid. The theory of the effective properties of
materials containing cracks has been developed by a number of au-
thors (e.g. O'Connell & Budiansky 1974; Anderson et al. 1974;
Nishizawa 1982). For an effective medium theory it is assumed
that the wavelengths are large compared with the size of the cracks
(<1, generally <0.1). In the studies referred to above the cracks are
isolated and no fluid can flow from one crack to another or within
single cracks during the passage of a seismic wave. The effect of flow
is that fluid pressure is released from the more highly compressed
cracks, or from more highly compressed regions within cracks, to
the less compressed, with consequent effects on the overall wave
speeds and attenuation.

Until quite recently, the analysis of wave propagation in porous
and permeable media was generally based upon the empirical ap-
proach pioneered by Biot (1956), in which no attempt was made to
provide a theoretical connection between the geometry of the porous
microstructure and the parameters appearing in the equations for the
waves. Subsequent attempts to give such a connection (e.g. Burridge
& Keller 1981) are hampered by the complexity of certain averaging
processes. Further progress along these lines could be made if, for
instance, an extremely simple model of the porous microstructure
were employed. One such model is that of circular cracks used by
Pointer et al. (2000). In this model they considered the transfer of
fluid by three different mechanisms:

(a) between cracks through seismically transparent pathways;

(b) within isolated partially saturated cracks and

(c) from the cracks into the background porous matrix. In each
model the cracks may be aligned or randomly oriented.

Theories for (a) interconnected cracks and (b, c) cracks with
equant porosity have already been developed for aligned cracks by
Hudson *et al.* (1996). A theory for isolated cracks that are partially
saturated with a single fluid was put forward by Hudson (1988).
Model (a) was also studied by O'Connell & Budiansky (1977) and
Mavko & Nur (1975). Partially saturated cracks (model b) have been
investigated by Walsh (1965) and Mavko & Nur (1979). Results
for cracks with equant porosity (model c) have also been given by Thomsen (1995) and Chapman (2003). The occurrence of fluid flow either between cracks within cracks, or from cracks into a porous matrix lowers the mechanical strength of rocks, and has a marked effect on seismic properties, even though the volume occupied by the cracks and pores may be small. Seismic velocities exhibit behaviour that is intermediate between that of empty cracks and that of isolated liquid-filled cracks over the range of frequencies in which fluid flow is significant. There is high attenuation and dispersion that becomes anisotropic if the cracks are aligned (Pointer et al. 2000; Chapman 2003).

Pointer et al. (2000) found that for both connected cracks and equant porosity the dependence of the attenuation on frequency is a peaked curve approaching zero at either frequency limit. In all of these cases the peak attenuation rises with decreasing compressibility. Correspondingly, the range of frequencies over which attenuation is significant decreases except for connected aligned cracks. However, for aligned interconnected cracks the behaviour with frequency is reversed because fluid flow is on a wavelength scale rather than on a local scale. They also found that there is negligible dispersion and attenuation of qSH waves, and showed that qP waves travel at a lower speed in the direction normal to the cracks, compared to the direction parallel to the cracks, when there is sufficient movement of crack fluid, as observed in laboratory experiments by Rathore et al. (1995). Finally, they showed that for aligned cracks the attenuation of qP waves is always greater than that of qSV. That is also true for randomly oriented cracks, except when the fluid flow is from crack to crack, when the opposite is true, because there is no local fluid flow between randomly oriented cracks under pure compression and hence no attenuation.

In this study, we concentrate on seismic attributes that respond directly to the presence of fluids, and in particular to quasi-static fluid effects on the compliances of fractured rocks (here by 'static' we mean that fracture patterns will not change as stress and pore pressure changes during fluid injection) and also we have not considered the wave-induced pressure changes, known as local or squirt flow as described by Hudson et al. (1996) [such changes should be considerably smaller than the pressure changes due to fluid flow injection, but their seismic effects may be significant, as demonstrated by Pointer et al. (2000), Chapman (2003) and Yang & Zhang (2002). We simulate fluid movement in a fractured network with a model of pore fluid pressure diffusion. From the simulation we can monitor pore pressure changes at consecutive times. As a result of the changes in pore pressure the effective stress changes. We use an empirical relationship to estimate the changes in compliances of the fractured rock due to the effective stress changes, and use a 2-D finite difference method to model the wave propagation before and during the fluid injection. The main aim of the study is to examine if there is any direct indication of pore pressure changes in synthetic seismograms for a common fracture network. Note that we only consider attenuation due to seismic scattering in the present work and viscoelastic loss and other fluid-related intrinsic attenuation mechanisms are not considered.

## 2 FLUID FLOW SIMULATIONS

Macroscopically, we may consider fluid flow in porous rocks as well as in fracture zones as a diffusive process with anisotropic diffusivities varying spatially and temporally by several orders of magnitude, and with pore pressures also presenting high local variations. When fluids are injected into a porous rock mass at a sufficiently high pressure we can have two possible types of fracture processes. Depending

on the fluid and rock properties and on the local stress field, hydraulic fracturing or induced seismicity may occur. However, any changes in the fracture network will almost certainly complicate the identification of pore pressure changes in seismic signatures, so for the purpose of the present analysis we keep the fracture network unchanged. Such models are called static in terms of the stress field, and the fracture network is used only to account for the porosity and permeability. A major simplification of the model is that the fluid is assumed to have the same bulk modulus $K_f$ as the solid $K_s$. This assumption was made in the original papers of Maillot & Main (1996) and Maillot et al. (1999) to simplify the simulation. Since in the present study, we are only interested in the effects of pore pressure diffusion, not the fluid saturation, on seismic waves, this assumption will not affect our results significantly. However, if the effects of fluid saturation are to be considered, this assumption, that is, $K_f = K_s$, will no longer be valid. We also assume that both fluid and solid phases are chemically inert and at constant temperatures, that the implicit void spaces are fully connected, that the porosity $\phi$ is uniform and constant. Only a single phase of fluid is considered. We combine mass conservation, Darcy's law and a linear equation of state [$\rho_f = \rho_{f0}(1+p/K_s)$], to obtain the time evolution of the fluid pressure $p$:

$$ \frac{\partial p}{\partial t} = \frac{\partial}{\partial x_i} \left( D_{ij} \frac{\partial p}{\partial x_j} \right) - \frac{1}{3} \frac{\partial \sigma_{ij}}{\partial t}, \quad (1) $$

where

$$ D_{ij}(x, t) = \frac{K_s}{\phi v} \kappa_{ij}(x, t). \quad (2) $$

is the diffusion tensor, with $v$ and $\kappa_{ij}(x, t)$, respectively, the viscosity of the fluid and the permeability of the matrix. We add the last term in eq. (1), which acts as a pore pressure source or sink, to include the effect of solid matrix isotropic stress variations (the stress is taken as positive in tension). The eq. (1) was derived by Maillot & Main (1996) and Maillot et al. (1999). We use a lattice Boltzmann method to solve eq. (1), which is valid in the most general media with anisotropic, heterogeneous and time-dependent diffusivity [eq. (2)] (see Maillot & Main 1996).

In the examples presented in this work, we have implemented the fluid flow model in 2-D using a 256 × 256 $d2p9$ lattice for the Bhatnagar, Gross and Krook (BGK) diffusion model. Following the terminology of Quian et al. (1992), a $d2q9$ lattice is a 2-D square lattice where each node is connected to eight neighbours: four horizontally and vertically, four at 45°, and itself. The boundary conditions are periodic, that is, the top side of the grid links to the bottom side, and the left side to the right side. The plane of computation is taken to be horizontal, justifying the absence of a gravity term in eq. (1) and an injection well is inserted in the centre. The dimension of the model is characterized by a length scale $L$ that represents the overall extent of the model, and a timescale $T$ that represents the duration of the fluid injection at the well. The spatial discretization is $\Delta x = 10$ m, and therefore, the dimension of the model is $2560 \times 2560$ m. The time step is $\Delta t_d = 1$ s. The diffusivity in the model was set as follows: in the unfractured (background) rock, it is isotropic and is equal to $D = 10^{-3}$ m$^2$ s$^{-1}$. In the fractures, the principal components are $D_1$ along the local direction of the fractures and $D_2$ normal to the fractures, with $D_1 = 10^4 D = 10$ m$^2$ s$^{-1}$ and $D_2 = 10^2 D = 0.1^{-3}$ m$^2$ s$^{-1}$.

We examine the case of a pre-existing fractured network, which is hydraulically conducting. The fractures are synthetic patterns generated by a stochastic multiscale cellular automaton model (Narteau 2006), also used in our previous publication (Vlastos et al. 2006).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Pore pressure map in a small area around the injection point (1000 × 1000 m). The figures show the pore pressure at four consecutive times, 10, 40, 70 and 100 hr, after the initial injection of the fluid (the unit of pore pressure, i.e. the colour bar, is in MPa). The black lines are the pre-existing fractures of the medium.</figcaption></figure>

The matrix is considered isotropic, but the fractures are aligned. In the fractured medium, fluid is injected in the centre of the model and we present a pore-pressure map at four consecutive time steps. Fig. 1 shows the pore-pressure maps at 10, 40, 70 and 100 hr after the initiation of injection. The area shown in this figure is a small area around the injection point, and not the whole modelled area. The black lines represent pre-existing fractures. The values of pore pressure are presented by a colour code explained on the right side of each map. The emergence of an elliptical pore pressure contour can be seen at 100 hr, with its semi-major axis aligned parallel to the fracture orientation.

are analysed based on changes in fracture compliance and fracture anisotropy as the traction on the fractures varies due to pore pressure changes, while properties of the intact background rock are assumed to be constant. During fluid injection, pore pressure generally increases, resulting in a decrease in the effective pressure on fractures, grain boundaries and microcracks. Lower static load on such surfaces increases the compliance in a non-linear fashion and increases fracture opening and/or pore throat size, decreasing the stiffness of the rock (decreasing compressional and shear velocities) and increasing permeability (Schoenberg 2002). In the following we will explain in detail all the theories and how we come to such a conclusions.

# 3 EFFECTS OF PORE PRESSURE

Fracture surfaces, grain boundaries, microcracks and joint faces are much more compliant and, therefore, sensitive to stress than the intact rock. Based on this assertion, the properties of the fractured rock

## 3.1 Estimation of fracture compliance

For a volume $V$, assumed to be homogeneous except for the presence of compliant surfaces across which displacement discontinuities can occur (Schoenberg & Sayers 1995),
$$ \epsilon_{ij} = [s_{ijkl_b} + s_{ijkl_f}] \sigma_{kl} = s_{ijkl_b} \sigma_{kl} + \frac{1}{2V} \sum_q \int_{S_q} ([u_i]n_j + [u_j]n_i) dS. \quad (3) $$

The strain tensor $\epsilon_{ij}$ and stress tensor $\sigma_{kl}$ are the mean dynamic strain and stress over the region of the surface areas $S_q$ due to the passing of a seismic wave whose wavelength is taken to be much larger than the linear extent of the volume $V$; $s_{ijkl_b}$ is the elastic compliance tensor of the background medium, $s_{ijkl_f}$ is the excess compliance tensor due to the presence of fractures or other surfaces across which displacement discontinuity occurs; $S_q$ is the surface area of $q$th displacement discontinuity; $[u_i]$ is the local displacement discontinuity across the surface; and $n_i$ is the local normal to the fracture surface. For a set of aligned compliant surfaces, or fractures, with orientation defined by unit normal $n_i$, it is assumed in the linear slip deformation theory that there exists a linear relation between the integral of the displacement discontinuity $[u_i]$ due to the passage of a wave across all the fractures and the dynamic traction $\sigma_{jk}n_k$ on a plane parallel to the fractures. This relation may be written:

$$ \frac{1}{V} \sum_q \int_{S_q} [u_i] dS = Z_{ij} \sigma_{jk} n_k, \quad (4) $$

where $Z_{ij}$ is the symmetric non-negative definite fracture compliance tensor (of dimension stress⁻¹) of this set of aligned fractures in the 1, 2, 3 coordinate system. Substituting equation (4) into eq. (3), we have,

$$ s_{ijkl} = s_{ijkl_b} + s_{ijkl_f} = s_{ijkl_b} + \frac{1}{4}(Z_{ik}n_i n_j + Z_{jk}n_i n_j + Z_{il}n_k n_j + Z_{jl}n_k n_i). \quad (5) $$

This expression for $s_{ijkl_f}$ satisfies the symmetry conditions of 4th rank compliance tensor, $s_{ijkl} = s_{ijkl}$, $s_{ijkl} = s_{ijkl}$, $s_{ijkl} = s_{klj}$. For a vertical fracture set with a horizontal unit normal given by

$$ n_i = \delta_{i1} \cos \theta + \delta_{i2} \sin \theta, \quad (6) $$

the excess fracture compliances are found by substituting eq. (6) into eq. (5), where $\theta$ is measured with respect to the fracture normal, that is, $\theta = 0°$ refers to the fracture normal, and $\theta = 90°$ refers to the fracture plane.

## 3.2 Estimation of effective static stress traction on aligned fractures

In the subsurface, first let static stress be of opposite sign than the usual convention so that compression is positive. Then, in a porous medium, the relationship between effective static stress $\sigma_{eff}$ and the applied external static stress $\sigma_{ext}$ is given by

$$ \sigma_{eff} = \sigma_{ext} - p f l, \quad (7) $$

where $p > 0$ is pore pressure, $f$ is a scalar empirical factor or the effective stress coefficient (usually $f \le 1$), and $l$ is the unity matrix. In some cases of anisotropy, $l$ should most likely be replaced by an anisotropic tensor. For isotropy, though, the usual practice is to assume $f = 1$ (so the effective stress is the same as differential stress). As pore pressure is isotropic, symmetry considerations dictate that on any surface, changes in pore pressure can only change the normal component of traction; tangential components of traction are invariant to pore pressure changes. One may think of a pore pressure drop as allowing the grains and other parts of the rock mass to settle into one another. Such compaction generally decreases the acoustic compliance of the compliant surfaces and the permeability of the

fracture network. A pore pressure rise generally does the opposite. The assumption that $Z$ depends on the static effective stress traction on the fracture faces, which will be the working assumption here, contains a far reaching implication—that effective stress changes influence the fracture parameters but do not necessarily cause a great number of new fractures to open suddenly, which would undermine the assumption of a more or less deterministic functional relationship between static stress and fracture parameters. Even when new fractures are included (e.g. as induced seismicity) they are likely to be preceded by a linear static response phase as described here.

We assume an anisotropic state of external static stress with one of the principal directions of the stress tensor being vertical with principal stress $\sigma_3$, and let the horizontal 1- and 2-directions be the other principal directions with principle stresses $\sigma_1$ and $\sigma_2$, respectively. Consider the set of compliant surfaces at a particular orientation with a normal unit vector normal specified in polar coordinates by

$$ \begin{bmatrix} \sin\phi\cos\theta \\ \sin\phi\sin\theta \\ \cos\phi \end{bmatrix}, \quad (8) $$

where $\phi$ is the polar angle and $\theta$ is the azimuthal angle. The effective stress traction, a vector acting on that fracture plane, is,

$$ \begin{bmatrix} \sigma_1 - p & 0 & 0 \\ 0 & \sigma_2 - p & 0 \\ 0 & 0 & \sigma_3 - p \end{bmatrix} \begin{bmatrix} \sin\phi\cos\theta \\ \sin\phi\sin\theta \\ \cos\phi \end{bmatrix} = \begin{bmatrix} (\sigma_1 - p)\sin\phi\cos\theta \\ (\sigma_2 - p)\sin\phi\sin\theta \\ (\sigma_3 - p)\cos\phi \end{bmatrix}. \quad (9) $$

The normal component of this traction, denoted by $\tau_{\perp}$, is

$$ \begin{align} \tau_{\perp} &= [\sin\phi\cos\theta \quad \sin\phi\sin\theta \quad \cos\phi] \begin{bmatrix} (\sigma_1 - p)\sin\phi\cos\theta \\ (\sigma_2 - p)\sin\phi\sin\theta \\ (\sigma_3 - p)\cos\phi \end{bmatrix} \nonumber \\ &= -p + \frac{\sigma_1 + \sigma_2}{4} + \frac{\sigma_3}{2} + \frac{\sigma_1 - \sigma_2}{4}\cos 2\theta \nonumber \\ &\quad - \left[ \frac{\sigma_1 + \sigma_2}{4} - \frac{\sigma_3}{2} + \frac{\sigma_1 - \sigma_2}{4} \right] \cos 2\phi. \tag{10} \end{align} $$

For vertical compliant surfaces, the normal is horizontal, $\phi = \pi/2$ and $\cos 2\phi = -1$. Then, from eq. (10), the normal component of traction becomes

$$ \tau_{\perp} = -p + \frac{\sigma_1 + \sigma_2}{2} + \frac{\sigma_1 - \sigma_2}{2}\cos 2\theta. \quad (11) $$

At any orientation, the tangential components of traction are independent of $p$.

## 3.3 The case of vertical fractures under anisotropic stress

For a medium with a continuous distribution (over orientation) of fracture sets, a fracture compliance density tensor $Z_{ij}(n_k)$ can be defined over a unit hemisphere (u.h.s.). From eq. (5), the excess compliance due to the fracturing may be written as an integral over the solid angle $\Omega$,

$$ s_{ijkl_f} = \frac{1}{4} \int_{u.h.s.} [Z_{ik}(n_q)n_i n_j + Z_{jk}(n_q)n_i n_j + Z_{il}(n_q)n_k n_j + Z_{jl}(n_q)n_k n_i] d\Omega. \quad (12) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Configuration of the model used for seismic modelling at consecutive stages of fluid injection. The source is located at the centre of the model (x = 1280 m, y = 1280 m) and the single receiver where we record the seismic waves is located at x = 1000 m and y = 300 m. The receivers used to examine azimuthal dependence are shown together with the azimuths for each case.</figcaption></figure>

Under anisotropic static stress, slip interfaces will be subjected to traction varying as a function of orientation, and that traction will have a tangential component except on fractures normal to a principal stress direction. For simplicity, we assume that the fracture sets, or sets of compliant surfaces, at all orientations are rotationally symmetric, and that they remain rotationally symmetric. Rotational

symmetry of the fracture set implies that the fracture compliance matrix in the coordinate system convenient to the fractures depends on two compliances, a normal compliance $Z_N$ and a tangential compliance $Z_T$. Following that the fracture compliance density tensor $Z_{ij}$ becomes

$$Z_{ij} = (\delta_{ij} - n_i n_j) Z_T + n_i n_j Z_N = Z_T \delta_{ij} + (Z_N - Z_T) n_i n_j. \quad (13)$$

In the model that we present here we consider only fluid injection into a flat layer cut by vertical fractures. In such a medium, the fracture normals lie in the (1,2)-plane and an arbitrary fracture normal is given in index notation by eq. (6). The integral over the unit hemisphere reduces to an integral over the unit semi-circle, without loss of generality, let $-\pi/2 < \theta < \pi/2$. For rotationally symmetric fractures, $Z_{ij}$ is given by eq. (13), and substitution of eqs (13) and (6) into eq. (12) gives

$$s_{ijklf} = \int_{-\pi/2}^{\pi/2} \frac{Z_T(\theta)}{4} [\delta_{ik}(\delta_{l1}\cos\theta + \delta_{l2}\sin\theta)(\delta_{j1}\cos\theta + \delta_{j2}\sin\theta) \\ + \delta_{jk}(\delta_{l1}\cos\theta + \delta_{l2}\sin\theta)(\delta_{i1}\cos\theta + \delta_{i2}\sin\theta) \\ + \delta_{il}(\delta_{k1}\cos\theta + \delta_{k2}\sin\theta)(\delta_{j1}\cos\theta + \delta_{j2}\sin\theta) \\ + \delta_{jl}(\delta_{k1}\cos\theta + \delta_{k2}\sin\theta)(\delta_{i1}\cos\theta + \delta_{i2}\sin\theta)] d\theta \\ + \int_{-\pi/2}^{\pi/2} [Z_N(\theta) - Z_T(\theta)] (\delta_{i1}\cos\theta + \delta_{i2}\sin\theta)(\delta_{j1}\cos\theta \\ + \delta_{j2}\sin\theta)(\delta_{k1}\cos\theta + \delta_{k2}\sin\theta)(\delta_{l1}\cos\theta + \delta_{l2}\sin\theta) d\theta. \quad (14)$$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Frequency–length distribution of fractures in log–log plots corresponding to fracture distributions in Fig. 2.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Waveforms recorded at a receiver located at x = 1000 m and y = 300 m for consecutive pore pressure stages. The black waveform represents the pre-injection stage, and red, green, blue and orange represent stages of the fluid injection at 10, 40, 70 and 100 hr after the injection, respectively.</figcaption></figure>

For $n_3 = 0$ we have $s_{ij33_f} \equiv 0$ for arbitrary vertical fractures. For rotationally symmetric fractures, we see from eq. (14) that

$$s_{ijk3_f} \equiv 0, i, j, k = 1, 2. \qquad (15)$$

This results in a compliance matrix in 6 × 6 condensed notation of the form:

$$\begin{bmatrix} S_{11_f} & S_{12_f} & 0 & 0 & 0 & S_{16_f} \\ S_{12_f} & S_{22_f} & 0 & 0 & 0 & S_{26_f} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{44_f} & S_{45_f} & 0 \\ 0 & 0 & 0 & S_{45_f} & S_{55_f} & 0 \\ S_{16_f} & S_{26_f} & 0 & 0 & 0 & S_{66_f} \end{bmatrix},$$

which is the form of a monoclinic medium with up-down symmetry (i.e. the (1,2)-plane is a mirror plane of symmetry).

Let us consider the case of such a vertically fractured medium (with rotationally symmetric fractures) subjected to a changing anisotropic external stress field with principal external stresses $\sigma_1$, $\sigma_2$ and $\sigma_3$ in the 1-, 2- and 3-directions, respectively. As a preliminary simplifying assumption let the fracture compliances at any angle be independent of the tangential component of effective stress traction on the fracture faces and depend on just the normal component, $\tau_{\perp}(\theta)$, given by eq. (11).

A reasonable approach is to assume very compliant fractures at low normal stress with fractures approaching low values asymptotically as normal stress becomes large. Approximating such dependence by exponential decay functions, we may write,

$$\begin{align} Z_N(\theta) &= Z_{N\infty}(\theta) + [Z_{N_0}(\theta) - Z_{N\infty}(\theta)] e^{-\tau_{\perp}(\theta)/\tau_N}, \nonumber \\ Z_T(\theta) &= Z_{T\infty}(\theta) + [Z_{T_0}(\theta) - Z_{T\infty}(\theta)] e^{-\tau_{\perp}(\theta)/\tau_T}. \tag{16} \end{align}$$

That is the general case where the parameters governing the exponential decay functions are themselves functions of $\theta$. However, in the model examined here, we assume that $Z_{N\infty}$, $Z_{N_0}$, $\tau_N$, $Z_{T\infty}$,

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Spectra of the waveforms shown in Fig. 3, representing seismic waves recorded at different stages of the fluid injection. The colour code corresponds to that used in Fig. 3.</figcaption></figure>

$Z_{T_0}$, $\tau_T$ are independent of $\theta$ and eq. (16) becomes

$$\begin{align} Z_N(\theta) &= Z_{N\infty} + [Z_{N_0} - Z_{N\infty}] e^{-\tau_{\perp}(\theta)/\tau_N}, \nonumber \\ Z_T(\theta) &= Z_{T\infty} + [Z_{T_0} - Z_{T\infty}] e^{-\tau_{\perp}(\theta)/\tau_T}. \tag{17} \end{align}$$

The $\theta$ dependence arises just through the $\cos 2\theta$ dependence of $\tau_{\perp}$. Note that as $\tau_{\perp}$ is an even function of $\theta$, $Z_N$ and $Z_T$ are also even functions of $\theta$. Now integrating over all values of $\theta$ from $-\pi/2$ to $\pi/2$ will give the excess compliance tensor for this stressed medium
as a function of $\sigma_1, \sigma_2$ and pore pressure $p$. Immediately it may be seen that $S_{16_f} = S_{26_f} = S_{45_f} = 0$ as these involve integrating an odd function times an even function from $-\pi/2$ to $\pi/2$. Thus an immediate consequence of these assumptions and simplifications is that if the background medium is of orthorhombic or higher symmetry, the long wavelength equivalent medium is orthorhombic. So the $6 \times 6$ compliance matrix becomes

$$S_f = \begin{bmatrix} S_{11_f} & S_{12_f} & 0 & 0 & 0 & 0 \\ S_{12_f} & S_{22_f} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{44_f} & 0 & 0 \\ 0 & 0 & 0 & 0 & S_{55_f} & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{66_f} \end{bmatrix}. \quad (18)$$

In the simulations presented in this work, the coefficients $\tau_T$ and $\tau_N$ have been empirically set to 1.35 MPa. Also the compliances at zero stress have been set to $Z_{T_o} = 5.681 \cdot 10^{-10} \text{ Pa}^{-1}$, $Z_{N_o} = 2.8409 \cdot 10^{-11} \text{ Pa}^{-1}$, and the compliances at infinite stress are $Z_{T_\infty} = Z_{T_o}/5 = 1.136 \cdot 10^{-10} \text{ Pa}^{-1}$ and $Z_{N_\infty} = Z_{N_o}/2 = 1.4205 \cdot 10^{-11} \text{ Pa}^{-1}$.

## 4 NUMERICAL SIMULATION OF WAVE PROPAGATION DURING THE FLUID INJECTION

To model the seismic wave propagation accurately in a complicated network of fractures, we use a finite difference technique. For each cell of the finite difference grid the properties of the elastic material is given by an equivalent medium based on the theory used by Vlastos et al. (2003). The readers are directed to the papers by Vlastos et al. (2003) and Vlastos et al. (2006) for the detailed description of how to implement discrete fractures in finite difference grids. By applying a very dense grid we have very high resolution and accuracy in the representation of the rock properties. When the fluid injection starts, the pore pressure changes affect rock properties greatly. As shown above, fracture compliances change with the variation of pore pressure. At each time step of the fluid flow simulation, eq. (17) gives us the corresponding dynamic effective elastic properties of the rock (Schoenberg & Sayers 1995; Liu et al. 2000; Worthington & Hudson 2000). Therefore, we have a continuous feedback from the changing elastic properties due to the pore pressure changes into the seismic simulation. In that manner we can produce consecutive time-lapse data and examine the potential for extracting information about the pore pressure changes directly from seismic waves.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Snapshots at consecutive time steps (150, 250, 350, 450 ms after the source initiation) which show the difference between the simulation of the medium 10 hr after the injection of the fluid and the pre-injection stage.</figcaption></figure>
For the seismic simulation we also assume an isotropic back-
ground medium, so anisotropy is caused only by the presence of
fractures. In addition, energy loss during wave propagation in the
models presented here, is only due to scattering. The modelling tech-
nique can accommodate both single and multiple scattering effects.
The background medium parameters are $V_p = 3300 \text{ m s}^{-1}$, $V_s =
2000 \text{ m s}^{-1}$, and the density $\rho = 2200 \text{ kg m}^{-3}$. We use a $256 \times
256$ grid, with spatial grid-step 10 m and time step 0.001 s. The
source wavelet is a Ricker wavelet with a dominant frequency of
40 Hz, a pulse initial time of 0.1 s, and wavelengths for the P wave
82 m and for the S wave 50 m. The source is located at the centre of
the model ($x = 1280 \text{ m}, y = 1280 \text{ m}$), exactly at the same position
as the fluid injection point, so the waves travelling to the receivers
will be greatly affected by pore pressure changes.

**4.1 Single azimuth**

Initially, we record the seismic waves at a receiver located at x = 1000 m and y = 300 m, at consecutive stages of the injection numerical simulation. Fig. 2 shows the model configuration used for the seismic simulations. The source–receiver direction is at a 16° angle with the y-direction, while on average the direction of the fractures is at a 30° angle with the y-direction. So the source–receiver angle relative to the fracture normal direction on average is 104° (red line on Fig. 2). [We examine the variation for different

azimuths (remaining lines in Fig. 2) in the next section]. The frac-
ture network consists of fractures with a range of sizes from 30 to
300 m. Fig. 3 shows the power spectrum of the fracture size distri-
bution (plotted in log–log scale). The variation can be fitted with a
straight line, so we can say that the size distribution approximates
a power-law distribution. The plot shows that the number of frac-
tures with a certain size is conversely proportional to the size. That
means that there are a small number of large fractures and most of
the fractures have small size.

Fig. 4 shows the x-component of the waveforms recorded at the receiver located at x = 1000 m and y = 300 m. The black trace is recorded at the pre-injection stage when pore pressure is equal throughout the model and is used as a reference. Then we start the fluid injection and we get feedback from the pore pressure changes at 10, 40, 70 and 100 hr after the initiation of the injection as shown in Fig. 1. We simulate wave propagation in the medium at the four consecutive time steps of the fluid injection. The red, green, blue, and orange traces are recorded at the receiver at stages 10, 40, 70 and 100 hr after the injection initiation, respectively. There are clear variations in the features of the waveforms as pore pressure changes. The direct *P* waves recorded at 0.37 s, do not seem to be affected by the pore pressure changes and seem to maintain the same amplitudes. In contrast, the *S* wave and the coda waves exhibit strong amplitude changes which can be uniquely attributed to the pore pressure changes, since all other parameters remain constant. This

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Snapshots at consecutive time steps (150, 250, 350 and 450 ms after the source initiation), which show the difference between the simulation of the medium 40 hr after the injection of the fluid and the pre-injection stage.</figcaption></figure>
is in agreement with the results presented by Tod (2002); Liu *et al.* (2002) and Angerer *et al.* (2002). Increasing pore pressure results in a higher S-wave amplitude. We can see that the orange waveform has the higher amplitude, which represents the state 100 hr after fluid injection, when the pore pressure is higher compared with the remaining cases. On the other hand, in most of the cases of coda waves (after 0.65 s), they seem to be affected by pore pressure in the opposite way. So the coda has higher amplitude in the pre-injection stage (black waveform) and a lower amplitude when the pore pressure is the greater (orange waveform).

To examine further the effects of pore pressure on the waveform, we show in Fig. 5 the spectrum of each of the waveforms presented in Fig. 4. We can see a significant shift of the peak frequency towards lower frequencies as the injection proceeds. Starting from the pre-injection point where the maximum energy is near 40 Hz, we end up in the stage 100 hr after the injection where the maximum energy is near 30 Hz, indicating a systematic increase of scattering attenuation with the increasing pore pressure. The magnitude of the peak frequency shift is about 10 Hz. Another feature seen in the spectrum is that at the initial stages of the injection, at 10 and 40 hr the energy is distributed in two frequency ranges compared to the rest of the cases where the energy is concentrated in a certain frequency range. This may be because the initial stages are a sort of a

transition state of the system between the pre- and the after-injection state.

It is a usual practice in time-lapse seismic monitoring to examine the difference between measurements for two consecutive time steps, to evaluate the effect of pore pressure changes. For the model of Fig. 2, we conduct forward modelling with the same configuration as in Fig. 2, and for the pre-injection state, and the states 10, 40, 70 and 100 hr after the fluid injection. From each simulation, snapshots of the seismic wave are generated at 150 ms, 250 ms, 350 ms and 450 ms, after the initiation of the source. To examine the effect of pore pressure changes, we take the pre-injection stage as a reference stage and find the difference between the snapshots of each stage after injection and the pre-injection stage. Figs 6–9 show the snapshots at consecutive time steps, which are the results of the difference between the simulation stages at 10, 40, 70 and 100 hr after the fluid injection and the pre-injection stage, respectively. From all the figures we can see that the area of strong differences in the seismic signal have an elliptical shape and that its long axis almost follows the direction of the fractures. Also the strong differences are concentrated in the centre of the model, which is exactly where the fluid is injected, the highest differences in pore pressure. Those features indicate that the dual simulation shown here can map with accuracy the effect of pore pressure changes in

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Snapshots at consecutive time steps (150, 250, 350 and 450 ms after the source initiation), which show the difference between the simulation of the medium 70 hr after the injection of the fluid and the pre-injection stage.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Snapshots at consecutive time steps (150, 250, 350 and 450 ms after the source initiation), which show the difference between the simulation of the medium 100 hr after the injection of the fluid and the pre-injection stage.</figcaption></figure>

seismic wave propagation. By taking a closer look, especially at each snapshot taken 450 ms after the initiation of the source, when the seismic wave has covered a significant part of the modelled area, we can see that gradually the area of strong difference spreads gradually from the injection point outwards, following an ellipse which is exactly the shape of the fluid front.

## 4.2 Azimuthal dependence

Another important aspect of the pore pressure effect is azimuthal dependence. We conducted the same simulations as above and recorded the seismic waves at three receivers, at the same distance from the source, so we expect almost the same attenuation due to the distance travelled by the seismic wave, and at 90°, 130° and 180°, from the fracture normal. This is repeated for each of the four states of the fluid flow simulation shown in Fig. 1. Fig. 10 shows the differences in the horizontal (X) components recorded at the three receivers at each azimuth. The differences are computed between: (a) 10 and 40 hr after injection, (b) 10 and 70 hr after injection, and (c) 10 and 100 hr after injection. Fig. 11 shows the corresponding difference in the frequency spectra.

In general, this is a further confirmation that P waves are not greatly affected in contrast to S waves and coda (or scattered) waves. Along the fractures (azimuth 90°) we can see the strongest difference

in S waveforms, while at directions normal to the fractures (azimuth 180°), the strongest difference is in coda (or scattered) waveforms. As the pore pressure increases further, this effect becomes stronger, as expected. It is interesting to note that there is a variation of the frequency content with azimuth. The greatest frequency shift occurs for azimuth of 180° where a significant amount of energy is shifted from 50 to 60 Hz towards 30–40 Hz. At 130° azimuth, energy moves between the same frequency ranges, but in a smaller degree. Finally, at 90° azimuth there is a much more limited shift of energy, and in addition there is significant energy present at the range of 50–60 Hz. In this case the energy is redistributed to both low and high frequencies, in a transition phase, before it shifts to systematically lower frequencies as angle increases from fracture normal. Note that the systematical azimuthal variations of scattered waves in fractured reservoirs have been observed in data and have been successfully modelled by Minsley *et al.* (2004) and Willis *et al.* (2004), and the results from this work provides further support to our claim that analysis of scattered waves may be used to characterize fractured reservoirs.

# 5 DISCUSSION AND CONCLUSIONS

We have conducted systematic dual numerical simulations of fluid flow and seismic wave propagation. We used a realistic model of a
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Differences in horizontal (X) components recorded at receivers having different azimuths from the fracture normals. Parts (a), (b), and (c) correspond to the differences between pore pressure 10 and 40 hr, 10 and 70 hr, and 10 and 100 hr after injection, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Differences in frequency spectra of the signals corresponding to the same cases as in Fig. 10. Parts (a), (b) and (c) represent increasing time difference between the examined stages as in Fig. 9.</figcaption></figure>
fractured network and conduct within the same model both the fluid flow and seismic simulation. For the fluid flow simulation the fluid is injected at the centre of a horizontal fractured layer, and at selected time steps after the injection information about pore pressure is collected. Variations in pore pressure lead to variations in the local effective stress. We used an empirical relationship between the effective stress changes and respective changes in the compliance of the rock. Therefore, at each selected time step of the fluid simulation, we obtain complete information about the updated elastic properties of the medium, and used them to perform seismic simulation. This process gives seismic data at consecutive time steps when pore pressure changes, which is synthetic time lapse seismic data.

Time-lapse (or 4-D) seismic data have proven their value in reservoir management. Time-lapse seismic is the comparison of 3-D seismic surveys at two or more points in time. Successive time-lapse images give better understanding of the movement of fluid phases, the spatial changes in pore fluid saturation and in volumetric pressure. Also they help in the identification of by-passed oil and in-fill drilling opportunities. Finally, 4-D data are used to constrain reservoir models, predict flow units and flow, and improve the performance of enhanced oil-recovery programs. One of the most successful 4-D projects to date was carried out in the Draugen field, offshore Norway for operator Shell and its partners. That involved a 1990 base survey and a 1998 monitor survey to obtain the results. The conclusions were that a very clear time-lapse signal could be observed on the difference data set. Also time lapse yielded information on the location of the waterfronts and seismic history matching of the dynamic reservoir model reduced uncertainties in the forecast for the production profile. Finally, time-lapse results impacted the location of a new production well within 6 months of the completion of the acquisition. This success story was the result of the combination of time-lapse data and a simulator that was updated by the real data. It is clear how important is the use of a simulating tool and how successful it can be when it is combined with real data.

4-D seismic modelling is also used as part of integrated tools for reservoir engineering, when it is combined with time-lapse seismic and production data. In such a study by Huang (2001) an initial static reservoir description is used for simulator model building. The simulator model is used to generate synthetic time-lapse seismic data, which in turn are compared to measured time-lapse data in a model optimization loop. Model optimization is accomplished through numerical minimization of an objective function formed from the errors between the data measured in the field and the ones predicted by the current model state. That is called seismic history matching. The potentials of 4-D time-lapse seismic modelling was also shown by Olden et al. (2001) in a work that examined the effect of fluid flow and stress changes in a hydrocarbon reservoir. They concluded that time-lapse seismic method could solve significant operational problems.

Fluid flows in the rock pores and inside fractures. In this study we assumed that the main path of the fluid is through fractures. An immediate effect of fluid is that it changes the elastic properties of a fracture. As a result the seismic signature of the fracture will change following the variations of pore pressure. That is the basic idea of time-lapse seismic data. Such data are used today to examine fluid flow and relevant pore pressure changes. Our results show a different response between *P* waves and *S* and coda waves to pore pressure changes. *P* waves seem to be less affected or affected in a limited way, while *S* and its coda waves are strongly affected. That means that the exploration techniques for such cases should be based more on *S*-wave records in order to be able to find even small changes in

pore pressure. Such small changes can be very important in cases of CO₂ injection. We also see that the amplitudes increase with increasing pore pressure.

Frequency also is a critical tool, because frequency shift is strong as pore pressure increases. There is an important shift of the peak frequency towards lower frequencies (implying strong attenuation) as pore pressure increases. That means that seismic data need to cover a wide frequency range. There is also azimuthal dependence of the observed variations. Fractures comprise the main path where the fluid moves in the reservoir; that is why the fluid front is an ellipse with its long axis parallel to the fracture orientation. Parallel to the fracture orientations we observe very strong amplitude variation in the *S* waves as pore pressure changes, while normal to that direction the strong difference is observed in the coda (or scattered) waves. Finally, we show that the greater shift of energy in frequency happens when seismic waves travel normal to the flow path.

Our study can help to provide a greater insight in the systematic effects of pore pressure changes on seismic waveforms and attenuation, and to identify potential ways to estimate pore pressure changes from seismic data. The results presented here form part of our strategy to establish a direct link between pore pressure changes and potentially diagnostic variations on seismic waves. In this study we conclude that significant diagnostic interpretation may be made by examining time lapse *S* and *S*-wave coda data.

## ACKNOWLEDGMENTS

We thank our colleagues Russ Evans, Andy Chadwick (both at BGS) and John H. Queen (Hi-Q Geophysical Inc) for useful comments. We would also like to thank the editor C. Ebinger and two reviewers for their very constructive comments, which led to the significant improvement of this paper. This research is supported by the sponsors of the Edinburgh Anisotropy Project. This work is published with the permission of the Executive Director of the British Geological Survey (NERC) and the EAP sponsors: BG Group, BP, Chevron, CNPC, ConocoPhillips, ENI-Agip, GX Technology, Kerr-McGee, Landmark, Marathon, Norsk Hydro, PetroChina, PDVSA, Schlumberger/WesternGeco, Total, and Veritas DGC.

## REFERENCES

Anderson, D.L., Minster, B. & Cole, D.M., 1974. The effect of oriented cracks on seismic velocities, *J. geophys. Res.*, **79**, 4011–4015.

Angerer, E., Crampin, S., Li, X.-Y. & Davis, T.L., 2002. Processing, modelling and predicting time-lapse effects of overpressured fluid-injection in a fractured reservoir, *Geophys. J. Int.*, **149**, 267–280.

Best, A.I., 1994. Seismic attenuation anisotropy in reservoir sedimentary rocks, *64th Ann. Int. Mtg., Soc. Expl. Geophys.*, 822–825.

Best, A.I., McCann, C. & Southcott, J., 1994. The relationships between the velocities, attenuations and petrophysical properties of reservoir sedimentary rocks, *Geophys. Prospect.*, **42**, 151–178.

Biot, M.A., 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid, I. Low frequency range and II. High frequency range, *J. acoust. Soc. Am.*, **28**, 168–191.

Burridge, R. & Keller, J.B., 1981. Poroelasticity equations derived from microstructure, *J. acoust. Soc. Am.*, **70**, 1140–1146.

Chapman, M., 2003. Frequency dependent anisotropy due to meso-scale fractures in the presence of equant porosity, *Geophys. Prospect.*, **51**, 369–379.

Cheng, C.H. & Toksöz, M.N., 1979. Inversion of seismic velocities for the pore aspect ratio spectrum of a rock, *J. geophys. Res.*, **84**, 7533–7543.
Engelder, T. & Fischer, M.P., 1994. Influence of poroelastic behaviour on the magnitude of minimum horizontal stress, $S_h$, in overpressured parts of sedimentary basins, *Geology*, **22**, 949–952.

Huang, X., 2001. Integrating time-lapse seismic with production data: A tool for reservoir engineering, *The Leading Edge*, **20**, 1148–1153.

Hsu, C.-J. & Schoenberg, M., 1993. Elastic waves through a simulated fractured medium, *Geophysics*, **58**, 964–977.

Hudson, J.A., 1988. Seismic wave propagation through material containing partially saturated cracks, *Geophys. J. Int.*, **92**, 33–37.

Hudson, J.A., 2000. The effect of fluid pressure on wave speeds in a cracked solid, *Geophys. J. Int.*, **143**, 302–310.

Hudson, J.A. & Liu, E., 1999. Effective elastic properties of heavily faulted structures, *Geophysics*, **64**, 479–485.

Hudson, J.A., Liu, E. & Crampin, S., 1996. The mechanical properties of materials with interconnected cracks and pores, *Geophys. J. Int.*, **124**, 105–112.

Liu, E., Hudson, J.A. & Pointer, T., 2000. Equivalent medium representation of fractured rock, *J. geophys. Res.*, **105**, 2981–3000.

Liu, E., Tod, S.R. & Li, X.-Y., 2002. Effects of stress and pore fluid pressure on seismic anisotropy in cracked rock, *CSEG Recorder*, **27**, 93–98.

Lockner, D.A., Walsh, J.B. & Byerlee, J.D., 1977. Changes in seismic velocity as an attenuation during deformation of granite, *J. geophys. Res.*, **82**, 5374.

Maillot, B. & Main, I.G., 1996. A lattice BGK model for the diffusion of pore fluid pressure, including anisotropy, heterogeneity, and gravity effects, *Geophys. Res. Lett.*, **23**, 13–16.

Maillot, B., Nielsen, S. & Main, I.G., 1999. Numerical simulation of seismicity due to fluid injection in a brittle poroelastic medium, *Geophys. J. Int.*, **139**, 263–272.

Mavko, G.M. & Nur, A., 1975. Melt squirt in the asthenosphere, *J. geophys. Res.*, **80**, 1444–1448.

Mavko, G.M. & Nur, A., 1978. The effect of nonelliptical cracks on the compressibility of rocks, *J. geophys. Res.*, **83**, 4459–4468.

Mavko, G.M. & Nur, A., 1979. Wave attenuation in partially saturated rocks, *Geophysics*, **44**, 161–178.

Narteau, C., 2006. Formation and evolution of a population of strike-slip faults in a multiscale cellular automaton model, *Geophys. J. Int*, in press.

Minsley, B., Willis, M., Krasovec, M., Burns, D. & Toksoz, N., 2004. Investigation of a fractured reservoir using P-wave AVOA analysis: a case study of Emilio Field with support from synthetic examples, *74th Ann. Internat. Mtg., Soc. of Expl. Geophys.*, 248–251.

Nishizawa, O., 1982. Seismic velocity anisotropy in a medium containing oriented cracks-transversely isotropic case, *J. Phys. Earth*, **30**, 331–347.

Nur, A. & Simmons, G., 1969. Stress induced velocity anisotropy in rock: an experimental study, *J. geophys. Res.*, **74**, 6667–6674.

O'Connel, R.J. & Budiansky, B., 1974. Seismic velocities in dry and saturated cracked solids, *J. geophys. Res.*, **79**, 5412–5426.

O'Connel, R.J. & Budiansky, B., 1977. Visco-elastic properties of fluid saturated cracked solids, *J. geophys. Res.*, **82**, 5719–5735.

Olden, P., Corbett, P., Westerman, R., Somerville, J., Smart, B. & Koutsabelis, N., 2001. Modeling combined fluid and stress change effects in the seismic

response of a producing hydrocarbon reservoir, *The Leading Edge*, **20**, 1154–1163.

Pointer, T., Liu, E. & Hudson, A., 2000. Seismic wave propagation in cracked porous media, *Geophys. J. Int.*, **142**, 199–231.

Pyrak-Nolte, L.J., Myer, L.R. & Cook, N.G.W., 1990. Transmission of seismic waves across single natural fractures, *J. geophys. Res.*, **95**, 8617–8638.

Quian, Y.N., d'Humieres, D. & Lallemand, P., 1992. Lattice BGK models for Navier-Stokes equation, *Europhys. Lett.*, **17**, 479–484.

Rathore, J.S., Fjaer, E., Holt, R.M. & Renlie, L., 1995. Acoustic anisotropy of a synthetic sandstone with controlled crack geometry, *Geophys. Prospect.*, **43**, 805–829.

Schoenberg, M., 2002. Time dependent anisotropy induced by pore pressure variation in fractured rock, *J. seismic Exploration*, **11**, 83–105.

Schoenberg, M. & Sayers, C.M., 1995. Seismic anisotropy of fractured rock, *Geophysics*, **60**, 204–211.

Shapiro, S.A. & Gurevich, B., 2002. Seismic signatures of fluid transport—introduction, *Geophysics*, **67**, 197–198.

Thomsen, L., 1995. Elastic anisotropy due to aligned cracks in porous rock, *Geophys. Prospect.*, **43**, 805–829.

Tod, S.R., 2002. The effects of stress and fluid pressure on the anisotropy of interconnected cracks, *Geophys. J. Int.*, **149**, 149–156.

Vasco, D.W., 2004. Seismic imaging of reservoir flow properties: time-lapse pressure changes, *Geophysics*, **69**, 511–521.

Vlastos, S., Liu, E., Main, I. & Li, X.-Y., 2003. Numerical simulation of wave propagation in media with discrete distributions of fractures: effects of fracture sizes and spatial distributions, *Geophys. J. Int.*, **152**, 649–668.

Vlastos, S., Narteau, C., Liu, E. & Main, I.G., 2006. Numerical simulation of wave propagation in fractured media: scattering attenuation at different stages of the growth of a fracture population, *Geophys. J. Int.*, in press.

Walsh, J.B., 1965. Seismic attenuation in partially saturated rock, *J. geophys. Res.*, **100**, 15407–15424.

Willis, M., Rao, R., Burns, D., Byun, J. & Vetri, L., 2004. Spatial orientation and distribution of reservoir fractures from scattered seismic energy, *74th Ann. Internat. Mtg., Soc. of Expl. Geophys.*, 1535–1538.

Winkler, K. & Nur, A., 1979. Pore fluid and seismic attenuation in rocks, *Geophys. Res. Lett.*, **6**, 1–4.

Winkler, K. & Murphy, W.F., 1995. Acoustic velocity and attenuation in porous rocks, *A Handbook of Physical Constants: Rock Physics and Phase Relations, American Geophysical Union, Washington, AGU Reference Shelf 3*, 20–34.

Worthington, M.H. & Hudson, J.A., 2000. Fault properties from seismic Q, *Geophys. J. Int.*, **143**, 937–944.

Wyllie, M.R.J, Gregory, A.R. & Gardner, G.H.F., 1958. An experimental investigation of factors affecting elastic wave velocities in porous media, *Geophysics*, **23**, 459–493.

Yang, D.H. & Zhang, Z.J., 2002. Poroelastic wave equation including the Biot/Squirt mechanism and the solid/fluid coupling anisotropy, *Wave Motion*, **35**, 223–245.