Seismic monitoring of reservoir fluid flow:
Fundamental theory and examples

David E. Lumley

ABSTRACT

Time-lapse 3-D seismic monitoring of subsurface rock property changes incurred during reservoir fluid-flow processes is an emerging new diagnostic technology for optimizing hydrocarbon production. I discuss the physical theory relevant for three-phase fluid flow in a producing oil reservoir, and rock physics transformations of fluid-flow pressure, temperature and pore-fluid saturation values to seismic P-wave and S-wave velocity. I link fluid-flow physical parameters to seismic reflection data amplitudes and traveltimes through elastic wave equation modeling and imaging theory. I demonstrate with synthetic and field data examples that changes in fluid flow can be monitored and imaged from repeated seismic surveys acquired at varying production calendar times.

INTRODUCTION Hydrocarbon reservoirs are increasingly recognized as spatially heterogeneous entities, in terms of pore-fluid content, pore-fluid saturation, porosity, permeability, lithology, and structural control. Knowledge of these reservoir parameters and their spatial variation is critical in the evaluation of the total volume of hydrocarbon reserves in place, in understanding and predicting physical processes in the reservoir such as fluid flow and heat transfer, and in projecting and monitoring reservoir fluid production and recovery as a function of time. A new diagnostic role for seismic has recently been proposed in which several repeat 3-D seismic surveys are acquired in time-lapse mode to monitor reservoir production fluid-flow processes, e.g, Nur . This 4-D monitoring concept offers a method for better characterization of reservoir complexity by monitoring the flow of fluids with time in a producing reservoir.

In this paper, I attempt to link fluid-flow physics, rock physics, and the physics of seismic wave-theoretic modeling and imaging. Figure schematically shows how these three disciplines are coupled in the seismic fluid-flow monitoring problem, and how the critical physical parameters in each discipline are related. The equations of fluid-flow describe changes in pore pressure p, temperature T, and multi-phase pore-fluid saturation S due to fluid flow in porous media. Rock physics transformations relate these fluid-flow parameters to seismic compressional-wave and shear-wave propagation velocities $\alpha$ and $\beta$.Elastic wave theory demonstrates that scattered wave amplitudes A and traveltimes $\tau$ of reflected seismic waves ${\bf u}^{^{P\!P}}$ contain information about the fluid-flow parameters, and more importantly, that time-varying aspects of fluid-flow can be isolated from static background geology and identified separately by images constructed from multiple time-lapse seismic data sets.

 

forward
Figure 1 A schematic showing the critical physical parameters in each discipline of fluid flow, rock physics, and seismic wave theory, and how they are naturally coupled together in the seismic fluid-flow monitoring problem. Fluid-flow parameters are permeability $\kappa$, porosity $\phi$, viscosity $\eta$, pressure p, temperature T, and saturation S_i. Rock physics parameters are bulk and shear moduli K and $\mu$, density $\rho$, compressional and shear velocities $\alpha$ and $\b$. Seismic parameters are traveltime $\tau$, amplitude A, compressional and shear reflectivities R_p and R_s.

FLUID-FLOW THEORY A simple view of a hydrocarbon reservoir under production can be modeled as an isothermal (constant reservoir temperature), immiscible (no chemical fluid mixing) three-phase fluid flow, e.g., Dake . The three phases present in the pore space of the reservoir rock are oil, gas and water. During production, pore pressure decreases near oil producing wells, and increases near gas/water/steam injection wells, stimulating three-phase fluid flow in the three spatial dimensions ${\bf x}$ of the reservoir as a function of production time t. This fluid flow can be approximated by coupling fluid-mass conservation with Darcy's Law, which relates a gradient in pore pressure $p({\bf x},t)$ to the rate of fluid-flow $q({\bf x},t)$, given the permeability $\kappa({\bf x})$ and porosity $\phi({\bf x})$ of the rock, and the viscosity $\eta({\bf x})$of the fluid. Assuming the porosity varies slowly in space, the three phases of fluid-flow are coupled by the following immiscible fluid displacement equations:

 

$$\nabla\left(\frac{\kappa_o}{\eta_o}\nabla p_o\right) - \phi\partial_tS_o = Q_o$$ (1)

 

$$\nabla\left(\frac{\kappa_w}{\eta_w}\nabla p_w\right) - \phi\partial_tS_w = Q_w$$ (2)

 

$$\nabla\left(\rho_g\frac{\kappa_g}{\eta_g}\nabla p_g\right) - \phi\partial_t(\rho_g S_g) = \rho_g Q_g \;\; \;.$$ (3)

The subscripts o, w, and g refer to oil, water and gas respectively. S_i is the saturation of the ith fluid component in the pore space on a scale from zero to unity, and Q_i is a fluid source or sink term which can represent fluid withdrawal from a producing well, or fluid addition from an injection well. p_i is the partial pressure for each phase of oil, gas or water. Equations (2) and (1) assume incompressible fluid, whereas (3) incorporates the significant expansion and compression effects of gas under variable pressure conditions by including the gradient terms of gas fluid density $\rho_g$.These flow equations are coupled with the statements that the total pore saturation is complete and conserved with time:

 

S_w + S_o + S_g = 1 (4)

 

$$\partial_tS_w + \partial_tS_o + \partial_tS_g = 0 \;.$$ (5)

Equations (1)-(5) describe simple three-phase fluid flow in the hydrocarbon reservoir over production time. They are typically solved on a variable 3-D reservoir mesh by finite-difference or finite-element methods for the pressure and saturation spatial distributions at several time steps. An example of a fluid-flow simulation mesh is shown in Figure for a faulted and uplifted reservoir in the Troll field, offshore Norway. The important parameters to monitor are the evolution of the pore pressure and saturation changes in space and time. More complex equations are required to model thermal effects from steam injection processes, miscible floods in which the fluid phases are allowed to mix by chemical reaction, and other complicated phenomena such as oil fractionation, emulsions, fluid interfingering and gravity override.

 

res-grid
Figure 2 An expanded section of a reservoir fluid-flow simulation grid.

ROCK PHYSICS TRANSFORMATIONS

Given the fluid-flow equations of the previous section, and some description of the reservoir geology, we can use rock physics analysis to transform reservoir pressure, temperature and fluid saturation data into seismic parameters. The most important parameters are the particle displacement velocities of the elastic waves that may propagate and scatter through the reservoir. These seismic wave velocities are denoted $\alpha({\bf x})$ and $\b({\bf x})$ for the compressional (P) and shear (S) waves respectively.

Typically, dry rock properties for the reservoir are measured in the lab from core samples as a function of mineralogy, porosity, pressure and temperature. Then, effective fluid bulk moduli are computed for three-phase fluid mixtures of oil, gas and water, including the effects of temperature and pressure. Finally, saturated rock properties are calculated using Gassmann's equation to combine the dry-rock data and effective fluid moduli as a function of pressure, temperature, porosity, and fluid saturation.

Dry rock properties

Rock cores obtained from boreholes in the reservoir are cleaned and oven-dried prior to dry rock measurements. Dry rock porosity $\phi$and density $\rho$ measurements are performed. Compressional and shear wave velocities are measured in the dry core samples with an ultrasonic wave generator, oscilloscope, and computer controlled measurement apparatus. Traveltimes for the P and S waves to propagate in the core sample are measured at the 100 kHz to 1 MHz frequency range, and dry rock values for $\alpha$ and $\b$ are computed. Figure shows an example of compressional waveforms measured across different core samples in a lab experiment, Lumley et al. . These measurements may be repeated under varying lab conditions of confining pressure and temperature to map out the response of a dry rock sample to reservoir pore pressure and temperature. Based on the dry $\alpha$, $\b$ and $\rho$ data, the dry bulk moduli K_dry and dry shear moduli G_dry of the core samples can be obtained using the relation:

 

$$K_{dry}= \rho(\alpha^2 - \frac{4}{3} \b^2) \;\;\;\; ; \;\;\;\; G_{dry}= \rho \b^2 \;,$$ (6)

where $\rho$ is the dry density.

 

pwaves
Figure 3 High-frequency compressional waveforms measured in various core samples under rock physics lab conditions. The lower two waveforms were measured in dry and saturated Massillon sandstone cores respectively.

Saturated rock properties

Unfortunately, ultrasonic lab measurements of saturated rock properties are not representative of field saturated rock properties in the surface-seismic frequency bandwidth (10-200 Hz). This is due to dispersive wave effects caused by frequency-dependent fluid oscillations in the core sample pore space at ultrasonic frequencies. However, the velocity effect of saturation can be calculated at seismic frequencies by using Gassmann's formulas, e.g., Bourbié et al. . This formula relates the effective elastic moduli of a dry rock to the effective moduli of the same rock containing fluid at low frequencies:

 

$$K_{sat}= K_{solid}\; \frac{\phi K_{dry}-(1+\phi)K_{fluid}K_... ...{(1-\phi)K_{fluid}+ \phi K_{solid}- K_{fluid}K_{dry}/K_{solid}}$$ (7)

and  

$$G_{sat}= G_{dry}\;,$$ (8)

where K_sat and G_sat are the effective bulk and shear moduli of the saturated rock. Gassmann's relations require knowledge of the effective shear and bulk moduli of the dry rock (G_dry and K_dry), the bulk modulus of the mineral material making up the rock (K_solid), the effective bulk modulus of the saturating pore fluid (K_fluid), and the porosity $\phi$. Equation (7) is used to compute the low-frequency bulk modulus of saturated rock from high-frequency dry rock data.

For partially saturated rocks at sufficiently low frequencies, we can use an effective modulus K_fluid for the pore fluid that is an isostress average of the moduli of the liquid and gaseous phases:

 

$$\frac{1}{K_{fluid}} = \frac{S}{K_{liquid}} + \frac{(1-S)}{K_{gas}} \;.$$ (9)

This requires knowledge of the bulk modulus of the liquid phase (K_liquid), the bulk modulus of the gas phase (K_gas), and the saturation values (S). In general, if the pore fluid includes more than two phases, we can calculate the mixture's effective bulk modulus K_fluid based on the the number of fluid components N, the volumetric concentrations c_i of the ith component, and their bulk moduli K_i:

 

$$\frac{1}{K_{fluid}} = \left. \sum_{i=1}^{N} \frac{c_i}{K_i} \right. \;.$$ (10)

Finally, we can use the following formulas to find seismic velocities in saturated rocks:

 

$$\alpha_{sat}= \sqrt{ \frac{(K_{sat}+ \frac{4}{3}G_{sat}) }{\rho_{sat}} }$$ (11)

and  

$$\beta_{sat}= \sqrt{ \frac{G_{sat}}{\rho_{sat}} } \;,$$ (12)

where $\rho_{sat}$ is the density of the saturated rock:

 

$$\rho_{sat}= (1-\phi)\rho_{solid}+ \phi\rho_{fluid}\;,$$ (13)

$\rho_{solid}$ is the density of the solid phase, and $\rho_{fluid}$ is the density of the fluid mixture obtained as an arithmetic mean of the volume-concentration-weighted fluid density components $\rho_i$:

 

$$\rho_{fluid}= \sum_{i=1}^N c_i\rho_i \;.$$ (14)

This process allows us to calculate fluid properties that depend on fluid-flow saturation values. These properties will also depend on pressure and temperature through the dry rock measurements, and through the variation in bulk modulus and density of reservoir gas. Pressure- and temperature-dependent gas properties may be calculated as discussed by Batzle and Wang . This total combined rock physics analysis allows us to calculate seismic velocities in saturated reservoir rock as a function of mineralogy, fluid type and saturation value, pore pressure and temperature. Therefore, we can map P-wave and S-wave velocity and density as a function of the reservoir grid directly from the fluid-flow simulation values of pressure, temperature, and oil, gas and water saturation.

SEISMIC WAVE THEORY

Given the fluid-flow pressure, temperature and saturation data, mapped to seismic P-wave and S-wave velocity and density, the response of these fluid-flow changes can be modeled and imaged in seismic data by considering basic elastic wave theory.

Elastic wave modeling

Consider the elastodynamic wave equation for a seismic particle displacement vector wavefield $\u({\bf x},\omega)$ and a second order tensor stress field ${\mbox{$\boldmath\sigma$}}({\bf x},\omega)$ due to a body force vector excitation ${\bf f}({\bf x},\omega)$:

 

$$\nabla{\bf \cdot}{\mbox{$\boldmath\sigma$}}-\rho \omega^2 \u = {\bf f}\;,$$ (15)

e.g., Aki and Richards . Assume further a linear elastic stress-strain relationship in the material continuum such that

 

$${\mbox{$\boldmath\sigma$}}= {\bf C}{\bf :}\nabla\u \;,$$ (16)

where ${\bf C}({\bf x})$ is the fourth-order elastic stiffness tensor C_ijkl, and the ``${\bf :}$'' symbol means a second order inner contraction. A volume integral representation for the $\grave{P}\!\acute{P}$ ``P-wave to P-wave'' scattered wavefield ${\bf u}^{^{P\!P}}$ can be expressed as:

 

$${\bf \hat{a}}_r{\bf \cdot}{\bf u}^{^{P\!P}}= \int_{{\cal V}} {\bf f}^{^{P}}{\bf \cdot}{\bf u}^{^{P}}\, d{\cal V}\;.$$ (17)

Equation (17) is the volume integral representation of the reflected $\grave{P}\!\acute{P}$ wavefield ${\bf u}^{^{P\!P}}({\bf x}_r,\omega;{\bf x}_s)$ measured at a receiver ${\bf x}_r$ along an arbitrary vector component direction ${\bf \hat{a}}_r$, due to the excitation of a body force reflection-diffraction scattering potential ${\bf f}^{^{P}}({\bf x};{\bf x}_s,{\bf x}_r)$ at each subsurface point ${\bf x}$, excited by the incident source wavefield ${\bf u}^{^{P}}({\bf x},\omega;{\bf x}_s)$ generated by a seismic source located at ${\bf x}_s$.

The assumption of isotropic elastic WKBJ (ray-valid) Green's tensors for P waves leads to:

 

$${\bf u}^{^{P}}({\bf x},\omega;{\bf x}^{^{\prime}}) = A^{^{P}... ...^{\prime}})} = A^{^{P}}{\bf \hat{t}}^{^{P}}e^{i\omega\tau} \;,$$ (18)

where A^<<701>>P and $\tau$ are the ray-valid P-wave amplitude and traveltime from the ``source'' location ${\bf x}^{^{\prime}}$ to the ``observation'' point ${\bf x}$, and are related by the eikonal and transport equations respectively:

 

$$\vert\nabla\tau\vert^2 = \nabla\tau{\bf \cdot}\nabla\tau = \frac{1}{\alpha^2}$$ (19)

 

$$A\;\nabla^2\tau + 2\nabla\tau{\bf \cdot}\nabla A = 0 \;.$$ (20)

The unit vector ${\bf \hat{t}}^{^{P}}$ is the direction parallel to P-wave propagation, as shown in Figure , and is perpendicular to the wavefronts $\tau = \mbox{constant}$.

For a linear isotropic elastic solid, the stress-strain relationship is

 

$${\mbox{$\boldmath\sigma$}^{^{P}}}= \left[ \lambda(\nabla{\bf... ...{P}}){\bf I}+ 2\mu{\bf \hat{t}}^{^{P}}{\bf u}^{^{P}}\right] \;,$$ (21)

where $\lambda({\bf x})$ and $\mu({\bf x})$ are the Lamé parameters such that

 

$$\lambda= \rho\alpha^2 \;\;\;\;\; \mbox{and} \;\;\;\;\; \mu=\rho\b^2 \;,$$ (22)

and ${\bf I}$ is the second-order identity matrix $\delta_{ij}$.Lumley and Beydoun showed that the P-wave reflection-diffraction scattering potential ${\bf f}^{^{P}}$ is:

 

$${\bf f}^{^{P}}= -\frac{i\omega}{\alpha^2}\vert{\bf u}^{^{P}}... ...bf \cdot}({\bf \hat{t}}^{^{P}}{\bf \hat{t}}^{^{P}}) \right] \;.$$ (23)

Equation (23) is a body force equivalent for scattering-surface reflectivity excitations, and is clearly dependent on material property contrasts (gradients): $\nabla\alpha$, $\nabla\lambda$ and $\nabla\mu$.

After some algebraic manipulation, (17) can be represented as:

 

$${\bf \hat{a}}_r{\bf \cdot}{\bf u}^{^{P\!P}}({\bf x}_r) = \... ..._{sr}} \grave{P}\!\acute{P}\cos\phi_r \,d{\cal V}({\bf x}) \;.$$ (24)

To recap (24), the density at a subsurface point is denoted $\rho({\bf x})$, and the geometric reflection coefficient at that point is $\grave{P}\!\acute{P}({\bf x})$. The amplitude terms A_s and A_r represent the cumulative geometric spreading, transmission loss, Q-attenuation, etc., from the source and receiver to the subsurface point ${\bf x}$ respectively. The factor $\cos a_r$ involves the vector component projection at the surface location ${\bf x}_r$ of ${\bf u}^{^{P\!P}}$ onto the arbitrary direction ${\bf \hat{a}}_r$.The term $\tau_{sr}=\tau_s+\tau_r$ is the total traveltime from the source at ${\bf x}_s$ to the subsurface point ${\bf x}$ and back up to the receiver at ${\bf x}_r$.Finally, the diffraction weight $\cos\phi_r$represents the angle between the anticipated geometric specular reflection direction ${\bf \hat{t}}^{^{P\!P}}_s$ and the non-geometric diffraction direction ${\bf \hat{t}}^{^{P}}_r$. In the case of specular reflection when ${\bf \hat{t}}^{^{P\!P}}_s = {\bf \hat{t}}^{^{P}}_r$, $\phi_r=0$ and so $\cos\phi_r=1$.The generalized reflection ray and angle geometries are shown in Figure .

A linearized version of the nonlinear $\grave{P}\!\acute{P}$ reflection coefficient can be parameterized as:

 

$$\grave{P}\!\acute{P}({\bf x},\cos\bar{\theta}) \approx 0.5... ...heta}-0.5\tan^2\bar{\theta}\right) \frac{\Delta\rho}{\rho} \;,$$ (25)

where $\gamma=\b/\alpha$ and $\bar{\theta}$ is the reflection angle between the incident wave direction ${\bf \hat{t}}^{^{P}}$ and the local gradient of the compressional P-impedance structure $\nabla(\rho\alpha)$.This linearization is also given in a somewhat different parameterization by Aki and Richards , and is accurate when the relative perturbations in impedance and density are small, and the reflection angles $\bar{\theta}$ are less than the critical angle at which conical ``head waves'' emerge.

The seismic modeling equations (19), (20), (24) and (25) show that changes in fluid-flow pressure, temperature and saturation, mapped to $\alpha$, $\b$ and $\rho$ changes through rock physics transformations, will have an effect on the traveltimes $\tau$ and reflection amplitudes $\grave{P}\!\acute{P}$ in the seismic data ${\bf u}^{^{P\!P}}$. If several seismic surveys are recorded at different phases of production fluid-flow, the seismic response will change with calendar time due to the coupled equations in fluid-flow, rock physics and elastic wave theory shown above. For example, Figure shows modeled CMP gather seismograms, after moveout correction, before and after oil production from a horizontal well. The next section addresses the topic of imaging changes in fluid-flow directly from multiple seismic data sets recorded in ``monitoring'' mode.

 

anglegeom
Figure 4 Generalized reflection ray and angle geometries.

 

noise13-ann
Figure 5 Modeled CMP gather seismograms, before and after oil production from a horizontal well. Note near offset waveform traveltime and amplitude changes, and far offset AVO tuning.

Seismic wavefield imaging

Given a seismic data set recorded at some calendar time T₁, we would like to be able to image the subsurface reflectivity structure R₁ which generated the reflected waves observed in that seismic data. Furthermore, we would like to obtain several reflectivity estimates R₁, R₂, R₃,... corresponding to surveys over a producing reservoir at calendar times T₁, T₂, T₃..., and infer something about the change in subsurface fluid flow from the changes in the R_i maps. The required imaging procedure is called ``seismic migration'' in the seismic exploration industry.

I briefly derive a ``kinematic'' Kirchhoff prestack depth migration equation which is suitable for either 2-D or 3-D data acquisition, and incorporates single-arrival traveltime and phase estimates. This migration equation yields accurate estimates of reflectivity amplitudes for near-offset data, and provides an efficient and robust structural imaging condition for far-offset data, Lumley .

Given the Helmholtz variable-velocity $\alpha$ scalar wave equation

 

$$\left\{ \nabla^2 + (\omega/\alpha)^2 \right\} P({\bf x},\omega) = S({\bf x},\omega) \;,$$ (26)

the ``downgoing wavefield'' D generated by a single source at location ${\bf x}_s$ can be evaluated at any subsurface location ${\bf x}$ within a volume ${\cal V}$ from the frequency-domain integral representation:

 

$$D({\bf x},\omega;{\bf x}_s) = \int_{{\cal V}} G({\bf x},\omega;{\bf x}') \; S({\bf x}',\omega;{\bf x}_s) \;d{\bf x}' \;,$$ (27)

where $G({\bf x},\omega;{\bf x}')$ is the Green's function solution to (26) associated with the source location, and S is the source wave function. If we neglect the absolute amplitude of the source and consider only relative amplitudes in the migrated section, and assume the source has a compact delta function shape in both space and time: $\delta(t)\delta({\bf x}'-{\bf x}_s)$, then the downgoing wavefield can be approximated by the source Green's function alone: $D({\bf x},\omega;{\bf x}_s) \approx G({\bf x},\omega;{\bf x}_s)$.

The ``upgoing wavefield'' U reflected from the subsurface ${\bf x}$ due to a source at ${\bf x}_s$ can be reconstructed from the seismic (scalar) data $P={\bf \hat{a}}_r{\bf \cdot}{\bf u}^{^{P\!P}}$ recorded at receivers ${\bf x}_r$ using a Kirchhoff integral representation:

 

$$U({\bf x},\omega;{\bf x}_s) = \int_{\S} {\bf \hat{n}}{\bf \c... ...a;{\bf x}_r) \; P({\bf x}_r,\omega;{\bf x}_s) \; d{\bf x}_r\;,$$ (28)

where $G({\bf x},\omega;{\bf x}_r)$ is the receiver Green's function and ${\bf \hat{n}}$ is the unit vector normal to the recording surface $\S$ that bounds the image volume ${\cal V}$ of interest. The gradient operator $\nabla$ is taken with respect to the subsurface coordinate ${\bf x}$ along the recording surface at ${\bf x}={\bf x}_r$.

Given that the reflected wavefield U can be modeled as a convolution of the subsurface reflectivity R with the source wavefield D, a local least-squares estimate of R can be obtained as the weighted zero-lag correlation of the source and reflected wavefields: $R \approx \sum_{\omega} W\; U\,D^{\ast}$, where W are as yet unspecified weights, and $D^{\ast}$ is the complex conjugate of D. If this weighted zero-lag correlation is further averaged for all such single shot-profile migrations, the frequency-domain Kirchhoff migration equation becomes:

 

$$R({\bf x}) \approx \int_{\omega} \int_{{\bf x}_s} \int_{{\bf... ...bf x}_r,\omega;{\bf x}_s) \; d{\bf x}_r\,d{\bf x}_s\,d\omega\;.$$ (29)

It can be shown that the reflectivity image R is proportional to a reflection-angle averaged version of the $\grave{P}\!\acute{P}$ coefficient in the modeling equation (24), and is a first order estimate of the relative P-wave impedance contrast in the earth, $\Delta(\rho\alpha)/(\rho\alpha)$, as inferred from equation (25).

Assume a parametric form for the Green's functions G such that:

 

$$G({\bf x},\omega;{\bf x}_a) \approx A_a({\bf x};{\bf x}_a) \; e^{\pm i(\omega\tau_a + \phi_a)} \;,$$ (30)

where opposite signs are chosen in the exponential for the source (outgoing) and receiver (reverse-time extrapolated) Green's functions respectively. The parameters A_a, $\tau_a$ and $\phi_a$ are the single-arrival Green's function amplitudes, traveltimes and phase rotations from location ${\bf x}_a$ to ${\bf x}$. These parameters are often estimated by conventional high-frequency asymptotic ray methods.

Given the parametric form (30), an efficient time-domain version of (29) can be obtained as:

 

$$R({\bf x}) \approx \int_{{\bf x}_s} \int_{{\bf x}_r} \cos\th... ...P}({\bf x}_r,{\bf x}_s;t=\tau_{sr}) \;d{\bf x}_r\,d{\bf x}_s\;,$$ (31)

where the ``obliquity factor'' $\cos\theta_r$ is a function of the incident angle at each receiver with respect to the surface normal, and is obtained as the dot product $(\alpha_r \nabla\tau_{r}{\bf \cdot}{\bf \hat{n}})$.The Kirchhoff space-time migration equation (31) is a weighted diffraction stack of the preprocessed, deconvolved (but not divergence-corrected) data $\hat{P}$, after phase-rotation by the Green's function parameters $\phi_{sr}=\phi_s+ \phi_r$, evaluated along the diffraction trajectories given by the Green's function traveltimes $\tau_{sr}=\tau_s+\tau_r$.

I define kinematic migration by setting the migration weights $\hat{W}$ to unity. I note that (31) is suitable for 2-D migration if all spatial coordinates are 2-vectors, e.g. ${\bf x}=(x,z)$, and $\hat{P}$ is preprocessed by the half-time derivative operator $\sqrt{i\omega}$. However, (31) is equally suitable for 3-D migration if all spatial coordinates are 3-vectors, e.g. ${\bf x}=(x,y,z)$, and $\hat{P}$ is preprocessed by the full time derivative $\d_t$.

Seismic velocity analysis

The seismic migration equation (31) images subsurface reflectivity structure, which is proportional to short-wavelength impedance variation. However, (31) is a nonlinear function of wave propagation velocity $\alpha$ to first order in traveltimes $\tau$, and to second order in amplitudes A. The coherency of any reflectivity image is therefore dependent upon the accuracy to which the long-wavelength velocity structure $\alpha$ is known. Velocity estimation is nonlinear, and reflectivity estimation is linear, given a smooth estimate of the background velocity field. In general, the problem of estimating the short-wavelength reflectivity and the long-wavelength velocity are nonlinearly coupled, and need to be solved simultaneously. In practice, the problem is usually assumed to be separable and solved separately; first for the velocity and then for the reflectivity.

A general nonlinear inverse problem can be posed to solve for long-wavelength velocity and short-wavelength impedance as follows. A measure of coherency can be defined as some function $\cal F$ of the reflectivity image R, which is itself a function of the velocity $\alpha$,

 

$$\mbox{Coherency} \;\; = \;\; {\cal F} \{R(\alpha)\} \;.$$ (32)

Since the image is assumed to be most coherent at the correct velocity model, a nonlinear optimization problem ensues, e.g., Symes and Carazzone . An optimal velocity solution is obtained when the coherency does not improve with slight adjustments to the velocity model:

 

$$\frac{ \d {\cal F} \{R(\alpha)\} } {\d\alpha} \;\; = \;\; 0 \;\;\;\rightarrow \;\;\; \alpha({\bf x}) \;.$$ (33)

Equation (33) is the gradient of the image coherency with respect to velocity perturbation, and can be used in a nonlinear steepest-descent or conjugate-gradient method to iterate towards a final solution.

Seismic fluid-flow monitoring

The seismic migration equation (31) can be used to obtain an estimate of the subsurface reflectivity R for any seismic data set. This reflectivity image is an angle-average estimate of the elastic $\grave{P}\!\acute{P}$ scattering coefficient, and is a first order estimate of short-wavelength P-wave impedance contrasts, $\Delta(\rho\alpha)/(\rho\alpha)$, in the subsurface. Information on the long-wavelength velocity structure $\alpha({\bf x})$ is available in the seismic data from the traveltime information $\tau_{sr}$ in the wavefield ${\bf u}^{^{P\!P}}$, by solution of the nonlinear velocity analysis system (33). The velocity structure $\alpha$ and reflectivity image R estimated from a single seismic survey will be comprised of coupled contributions from the reservoir geology and the fluid-flow states in pore space. The estimation and interpretation of this information from a single seismic data set is known as seismic reservoir characterization:

$$\mbox{Characterization} \;\;=\;\; \alpha({\bf x}),\; R_p({\bf x}), \; R_s({\bf x}) \;\;=\;\; \mbox{rock + fluid} \;.$$

Seismic reservoir characterization is a very difficult task because of the ambiguity in trying to separate geology effects from fluid-flow effects in a single seismic data set.

However, when multiple seismic surveys are conducted at separate calendar times, it is expected that the reservoir geology will not change from survey to survey, but the state of fluid flow will change. Therefore, differencing a series of reflectivity images R_i and velocity model estimates $\alpha_i$ will remove the static geologic contribution to the seismic data, and isolate time-varying seismic changes in the reservoir which are due to time-varying fluid-flow changes alone. The process of estimating and comparing reflectivity images and velocity estimates from multiple seismic data sets recorded at different calendar times is known as seismic reservoir monitoring,

$$\mbox{Monitoring} \;\;=\;\; \d_t \;\;\left\{ \;\alpha({\bf x... ..., \; R_s({\bf x},t)\; \right\} \;\;=\;\; \mbox{fluid flow} \;.$$

Seismic reservoir monitoring is potentially a much less ambiguous task than characterization, because the effects of geology and fluid-flow may be uncoupled by comparing time-varying seismic data sets.

SYNTHETIC EXAMPLE

In this section, I present a simulated example of seismic reservoir monitoring on a North Sea reservoir, Lumley et al. , in order to determine if it would be feasible to observe changes in seismic monitor data given a specific fluid flow and reservoir geology scenario. Fluid-flow simulations were computed by Norsk Hydro for the case of solution-gas-drive oil production from a horizontal well in the reservoir after 56 and 113 days of production. Reservoir geology and core data measurements were combined to make rock physics transformations of the fluid-flow saturation, pressure and temperature data to seismic velocities and densities. Multi-offset elastic reflection seismogram surveys were simulated at each production phase. Prestack migration images show gas coning during production in reasonable seismic noise levels and frequency bandwidth.

Fluid-flow simulation

The oil zone is located at a depth of approximately 1550 meters, and it varies in thickness from a few meters to more than 20 meters. Horizontal drilling technology makes it possible to commercially produce oil from this thin oil zone. The reservoir fluid-flow grid, shown in Figure , consisted of 312 individual blocks, and each grid cell was given with a fluid-flow simulation data value of porosity, pore pressure, gas saturation, oil saturation, and water saturation. As an example of the simulation data, Figure shows the change in oil saturation after 113 days of production. The thin oil zone is compressed due to gascap expansion from above and upward water coning from below.

 

oil-sat
Figure 6 The change in oil saturation after 113 days of oil production. The thin oil zone at 1.55 km depth is compressed from above due to downward gascap coning (black), and upward water coning from below (dark gray). Black implies a 100% decrease in oil saturation.

Rock physics

Given the fluid-flow simulation data and reservoir geology, rock physics analysis is used to transform pressure, temperature and saturation data into $\alpha$, $\b$ and $\rho$ distributions within the reservoir. Estimates of dry rock properties in the reservoir are based on the results in a similar reservoir reported by Blangy and Strandenes and Blangy , which involved laboratory experiments conducted on 38 core samples. Figure shows an example of the $\alpha_{dry}$ core measurements made as a function of reservoir sand porosity at 20 MPa differential pressure.

 

vp-phi
Figure 7 P-wave velocity $\alpha$ measured versus porosity in 38 dry core samples from the Troll reservoir at 20 MPa differential pressure.

Figure shows an example of the change in P impedance mapped from the fluid-flow simulation at 113 days of oil production compared to the the initial reservoir state prior to production. There is a significant decrease in P impedance above the oil zone where the gascap has expanded downward, and a smaller increase in P impedance below the oil zone where the aquifer has coned upward during production. The maximum relative decrease in P impedance is about 15% which is strong enough to cause a significant seismic reflection.

 

ip-change
Figure 8 The change in P impedance after 113 days of production in units of (m/s g/cc). The gas cone above the oil zone shows a 15% decrease in impedance (black), and the aquifer uplift below shows a 5% increase in impedance (white), compared to initial conditions.

Seismic monitoring

Given the $\alpha$, $\b$ and $\rho$ distributions in the reservoir corresponding to the fluid-flow simulation data and rock physics transformations, I simulated a multi-offset surface seismic survey at each of three separate monitoring phases of oil production. These three prestack datasets were then processed to produce stacked and prestack-migrated reflectivity images. Difference images obtained by subtracting combinations of stacked and migrated sections clearly show that reservoir fluid production is visible in the seismic data and can be monitored in the presence of reasonable levels of seismic noise. Figure shows the difference section obtained by a simple subtraction of the Base Survey migration image, prior to oil production, from the Monitor 2 migration image, after 113 days of oil production. The gas coning clearly stands out from the background seismic noise as a bright spot at 2 km and 1.6 s, and accurately defines the spatial extent of the true impedance model anomaly.

 

mig13
Figure 9 A difference section of the Base Survey prestack migration subtracted from the Monitor 2 prestack migration.

FIELD DATA EXAMPLE

In this section I briefly show a field data example from a steamflood in the Duri field, Indonesia, first reported by Bee et al.  and Lumley .

Field data

Figure shows two inline profile-view sections cut from a 3-D migration cube. The left panel is before steam injection, the right panel is after 5 months of steam injection. After steam injection, the heat, pressure and desaturation of the steamed zone show up clearly as differences in the time-lapse seismic sections. The shallow portion of the borehole is even visible in the seismic data due to heating of the casing and surrounding rock.

Figure shows two time-slice map-view sections cut from a 3-D migration cube at about the depth that steam should be penetrating the reservoir. The lower panel shows a circular disk of steam ``S'' about 50 m in diameter (white), and a large outer annulus interpreted to be a transient pressure front ``P'' (dark gray). These examples clearly demonstrate that reservoir fluid-flow changes can be seen in some field conditions with time-lapse seismic monitor data.

CONCLUSION

I have discussed the physical theory relevant for three-phase fluid flow in a producing oil reservoir, and rock physics transformations of fluid-flow pressure, temperature and pore-fluid saturation values to seismic P-wave and S-wave velocity. I have linked fluid-flow physical parameters to seismic reflection data amplitudes and traveltimes through elastic wave equation modeling and imaging theory. I have demonstrated with both synthetic and field data examples that changes in fluid-flow can be monitored and imaged in certain conditions from repeated seismic surveys acquired at varying production calendar times.

ACKNOWLEDGMENTS

I thank Amos Nur, Jack Dvorkin and James Packwood of the Stanford Rock Physics Laboratory for their collaboration on the synthetic data example. Sverre Strandenes and Norsk Hydro were very helpful in providing the reservoir geology information and fluid-flow simulation data. The steamflood data were kindly provided by Steve Jenkins (CalTex), and Fred Herkenhoff and Ray Ergas (Chevron).

 

inline-migs-ann
Figure 10 Inline profile-view sections cut from a 3-D migration cube. The left panel is before steam injection, the right panel is after 5 months of steam injection. Steam near the central injection well causes bright reflection amplitudes and time delay of underlaying events. The polarity reversal on the base of reservoir reflection at 200 ms is probably due to the presence of a transient pressure front.

 

tslice-migs-ann
Figure 11 Time slice map-view sections cut from a 3-D migration cube. The top panel is before steam injection, the lower panel is after 5 months of steam injection. The steam zone ``S'' is located within the bright amplitude disk with a diameter of 50 m, and a pressure transient ``P'' is outlined by the black polarity reversal.

[MISC,GEOTLE,GEOPHYSICS,SEGCON,SEP]