## 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.

## 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). Conversely, pore pressure build-up due to injection leads to a decrease 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 _{h}*, increases at a rate proportional to but less than the rate of increase of pore pressure

*P*. That is a condition consistent with a pore pressure induced deformation of the rock, called poroelastic behaviour. That does not happen in normal pressured basins where friction governs

*S*(Engelder & Fischer 1994).

_{h}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. Parameters 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 confining pressure or decreasing pore pressure, cause velocities to increase 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 (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 velocity 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 pressure, both reduce attenuation in water-saturated rock. When stress is anisotropic, it causes velocities in rock to vary with direction (Lockner 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 (1994).

Fractured rock is modelled as a relatively rigid, defect free, background 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 discontinuity across the surface which is linearly related to the dynamic traction acting on the interface (to the first order). The dynamic elastic 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 linear parameters governing the infinitesimal slip on these planes have been shown to depend on the static stress state (Pyrak-Nolte 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 reservoir 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 orientations, 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. Seismic 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 parameters 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 significant effect on the properties of seismic waves of long wavelength propagating through the solid. The cracks may be partially or completely filled with fluid. The theory of the effective properties of materials containing cracks has been developed by a number of authors (e.g. O'Connell & Budiansky 1974; Anderson 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 approach 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 (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 (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 2000; Chapman 2003).

Pointer (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 (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 (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 (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.

## 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*. This assumption was made in the original papers of Maillot & Main (1996) and Maillot (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,

_{s}*K*=

_{f}*K*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 φ 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 [ρ

_{s}_{f}=ρ

*f*

_{0}(1 +

*p*/

*K*)], to obtain the time evolution of the fluid pressure p: (1+

_{S}*p*/

*K*

_{S})], to obtain the time evolution of the fluid pressure

*p*:

_{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 (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 *d*2*p*9 lattice for the Bhatnagar, Gross and Krook (BGK) diffusion model. Following the terminology of Quian (1992), a *d*2*q*9 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 Δ*x*= 10 m, and therefore, the dimension of the model is 2560 × 2560 m. The time step is Δ*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 2006). 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.

## 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 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.

### 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),

_{ij}and stress tensor σ

_{kl}are the mean dynamic strain and stress over the region of the surface areas

*S*due to the passing of a seismic wave whose wavelength is taken to be much larger than the linear extent of the volume is the elastic compliance tensor of the background medium, is the excess compliance tensor due to the presence of fractures or other surfaces across which displacement discontinuity occurs;

_{q}*S*is the surface area of

_{q}*q*th displacement discontinuity; [

*u*] is the local displacement discontinuity across the surface; and

_{i}*n*is the local normal to the fracture surface. For a set of aligned compliant surfaces, or fractures, with orientation defined by unit normal

_{i}*n*, it is assumed in the linear slip deformation theory that there exists a linear relation between the integral of the displacement discontinuity [

_{i}*u*] due to the passage of a wave across all the fractures and the dynamic traction σ

_{i}_{jk}

*n*on a plane parallel to the fractures. This relation may be written: where

_{k}*Z*is the symmetric non-negative definite fracture compliance tensor (of dimension stress

_{ij}^{−1}) of this set of aligned fractures in the 1, 2, 3 coordinate system. Substituting equation (4) into eq. (3), we have, This expression for satisfies the symmetry conditions of 4th rank compliance tensor,

*s*=

_{ijkl}*s*,

_{ijlk}*s*=

_{ijkl}*s*,

_{jikl}*s*=

_{ijkl}*s*. For a vertical fracture set with a horizontal unit normal given by the excess fracture compliances are found by substituting eq. (6) into eq. (5), where θ is measured with respect to the fracture normal, that is, θ= 0° refers to the fracture normal, and θ= 90° refers to the fracture plane.

_{klij}### 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 σ_{eff} and the applied external static stress σ_{ext} is given by

*p*> 0 is pore pressure,

*f*is a scalar empirical factor or the effective stress coefficient (usually

*f*≤ 1), and

*I*is the unity matrix. In some cases of anisotropy,

*fI*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 σ_{3}, and let the horizontal 1- and 2-directions be the other principal directions with principle stresses σ_{1} and σ_{2}, respectively. Consider the set of compliant surfaces at a particular orientation with a normal unit vector normal specified in polar coordinates by

_{⊥}, is For vertical compliant surfaces, the normal is horizontal, φ=π/2 and

*cos*2φ=−1. Then, from eq. (10), the normal component of traction becomes At any orientation, the tangential components of traction are independent of

*p*.

### 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*) 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 Ω,

_{k}*Z*and a tangential compliance

_{N}*Z*. Following that the fracture compliance density tensor

_{T}*Z*becomes

_{ij}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 −π/2 < θ < π/2. For rotationally symmetric fractures, *Z _{ij}* is given by eq. (13), and substitution of eqs (13) and (6) into eq. (12) gives

*n*

_{3}= 0 we have ≡ 0 for arbitrary vertical fractures. For rotationally symmetric fractures, we see from eq. (14) that This results in a compliance matrix in 6 × 6 condensed notation of the form: 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 σ_{1}, σ_{2} and σ_{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, τ_{⊥}(θ), 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,

That is the general case where the parameters governing the exponential decay functions are themselves functions of θ. However, in the model examined here, we assume that , , τ_{N}, , , τ

_{T}are independent of θ and eq. (16) becomes The θ dependence arises just through the

*cos*2θ dependence of τ

_{⊥}. Note that as τ

_{⊥}is an even function of θ,

*Z*and

_{N}*Z*are also even functions of θ. Now integrating over all values of θ from −π/2 to π/2 will give the excess compliance tensor for this stressed medium as a function of σ

_{T}_{1}, σ

_{2}and pore pressure

*p*. Immediately it may be seen that = = = 0 as these involve integrating an odd function times an even function from −π/2 to π/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 × 6 compliance matrix becomes In the simulations presented in this work, the coefficients τ

_{T}and τ

_{N}have been empirically set to 1.35 MPa. Also the compliances at zero stress have been set to = 5.681 · 10

^{−10}Pa

^{−1}, = 2.8409 · 10

^{−11}Pa

^{−1}, and the compliances at infinite stress are = /5 = 1.136

^{−10}Pa

^{−1}and = /2 = 1.4205

^{−11}Pa

^{−1}.

## 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 (2003). The readers are directed to the papers by Vlastos (2003) and Vlastos (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 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.

For the seismic simulation we also assume an isotropic background 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 technique can accommodate both single and multiple scattering effects. The background medium parameters are *V _{P}*= 3300 m s

^{−1},

*V*= 2000 m s

_{S}^{−1}, and the density ρ= 2200 kg m

^{−3}. We use a 256 × 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 m,

*y*= 1280 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.

### 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 fracture network consists of fractures with a range of sizes from 30 to 300 m. Fig. 3 shows the power spectrum of the fracture size distribution (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 fractures 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 is in agreement with the results presented by Tod (2002); Liu (2002) and Angerer (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 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.

### 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 (2004) and Willis (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.

## Discussion and Conclusions

We have conducted systematic dual numerical simulations of fluid flow and seismic wave propagation. We used a realistic model of a 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 spacial 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 (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_{2} 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 diagonstic 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

*S*, in overpressured parts of sedimentary basins

_{h}*P*-wave AVOA analysis: a case study of Emilio Field with support from synthetic examples

*A Handbook of Physical Constants: Rock Physics and Phase Relations*

**AGU Reference Shelf 3**,