Mathematical foundation of elastic wave propagation

Seismic Data Analysis

Series	Investigations in Geophysics
Author	Öz Yilmaz
DOI	http://dx.doi.org/10.1190/1.9781560801580
ISBN	ISBN 978-1-56080-094-1
Store	SEG Online Store

L.1 Stress-strain relation

Seismic waves induce elastic deformation along the propagation path in the subsurface. The equation of wave propagation in elastic solids are derived by using Hooke’s law and Newton’s second law of motion. We shall begin with the stress-strain relation for elastic solids.

Solid bodies, such as rocks are capable of propagating forces that act upon them. Imagine an infinitesimally small volume around a point within a solid body with dimensions (dx, dy, dz) as depicted in Figure L-1. Stress is defined as force per unit area. The stress acting upon one of the surfaces, say dy − dz, can in general be at some arbitrary direction. It can, however, be decomposed into three components: P_xx which is normal to the surface, and P_xy and P_xz which are tangential to the surface. The first subscript refers to the direction of the normal to the surface, and the second subscript refers to the direction of the stress component. The components of the stress that are normal to the surfaces associated with the infinitesimal volume shown in Figure L-1 are called the normal stress components and the components of the stress that are tangential to the surfaces are called the shear stress components. A normal stress component is tensional if it is positive, and compressional if it is negative. To retain the solid volume depicted in Figure L-1 in equilibrium requires nine stress components that make up the stress tensor

$\{P_{ij}\}={\begin{pmatrix}P_{xx}&P_{xy}&P_{xz}\\P_{yx}&P_{yy}&P_{yz}\\P_{zx}&P_{zy}&P_{zz}\end{pmatrix}}.$

(L-1)

The diagonal elements are the normal stress components and the off-diagonal elements are the shear stress components.

As the dimensions of the volume surrounding a point inside a solid body are made infinitesimally smaller, the sum of the moments of all surface forces about any axis must vanish. This requirement causes P_xy = P_yx, P_xz = P_zx, and P_yz = P_zy, making the stress tensor (L-1) symmetrical.

Fluids cannot support shear stress. In a fluid medium, only one independent stress component remains, since P_xx = P_yy = P_zz = −p, the hydrostatic pressure (the negative sign signifies the compressive nature of the hydrostatic pressure), and P_xy = P_xz = P_yz = 0.

Stress induces deformation of solid bodies. Consider two points, P and Q, within a solid body as indicated in Figure L-2. Subject to a stress field, the solid is deformed in some manner and the particles at points P and Q are displaced to new locations P′ and Q′. The displacement components δu, δv, and δw can be expressed by the relation

${\begin{pmatrix}\delta u\\\delta v\\\delta w\end{pmatrix}}={\begin{pmatrix}\partial u/\partial x&\partial u/\partial y&\partial u/\partial z\\\partial v/\partial x&\partial v/\partial y&\partial v/\partial z\\\partial w/\partial x&\partial w/\partial y&\partial w/\partial z\end{pmatrix}}{\begin{pmatrix}\delta x\\\delta y\\\delta z\end{pmatrix}}.$

(L-2)

To understand the physical meaning of the elements of the tensor containing the partial derivatives in equation (L-2), instead of arbitrary displacement as shown in Figure L-2, consider deformations of specific types. The simplest deformation is an extension in one direction as a result of a tensional stress as illustrated in Figure L-3a. (For simplicity, only one surface of the volume in Figure L-1 is considered.) The fractional change in length in the x-direction is δu/δx, which, in the limit as the volume becomes infinitesimally small, is defined as the normal strain component.

• 

  Figure L-1  Stress components acting on an infinitesimally small volume surrounding a point within an elastic solid.

• 

  Figure L-2  Two points, P and Q, within a solid body subject to a stress field. The solid is deformed in some manner and the particles at points P and Q are displaced to new locations P′ and Q′. Here, δu, δv, and δw are the displacement components.

Strain is a dimensionless quantity. There are three normal strain components, e_xx, e_yy, and e_zz; if the deformation is small, these are given by:

$e_{xx}={\frac {\partial u}{\partial x}},$

(L-3a)

$e_{yy}={\frac {\partial v}{\partial y}},$

(L-3b)

$e_{zz}={\frac {\partial w}{\partial z}}.$

(L-3c)

The stress field away from the typical seismic source is so small that it does not cause any permanent deformation on rock particles along the propagation path. Hence, the strain induced by seismic waves is very small, usually around 10⁻⁶. A positive strain refers to an extension and a negative strain refers to a contraction.

Figure L-3  Deformations caused by stress acting on one surface of the volume in Figure L-1: (a) linear deformation that results in extension of the side AB in the x direction by an amount BB′; (b) shearing only; (c) rotation only; (d) combined angular deformation (α) and rotation (β). See text for details.

Other types of deformation are caused by shearing (Figure L-3b), rotation (Figure L-3c), and a combination of the two (Figure L-3d). The angular deformations ξ and ζ, in the limit as the volume becomes infinitesimally small, can be expressed in terms of displacement components δu and δw, by using the relations

$\xi -\zeta ={\frac {\partial w}{\partial x}}$

(L-4a)

and

$\xi +\zeta ={\frac {\partial u}{\partial z}},$

(L-4b)

which, by adding both sides, yield

$\xi =e_{xz}={\frac {1}{2}}\left({\frac {\partial w}{\partial x}}+{\frac {\partial u}{\partial z}}\right)=e_{zx}.$

(L-5a)

Similarly, we obtain

$e_{xy}={\frac {1}{2}}\left({\frac {\partial v}{\partial x}}+{\frac {\partial u}{\partial y}}\right)=e_{yx}$

(L-5b)

and

$e_{zy}={\frac {1}{2}}\left({\frac {\partial w}{\partial y}}+{\frac {\partial v}{\partial z}}\right)=e_{yz}.$

(L-5c)

The quantities described by equations (L-3) and (L-5) make up the strain tensor

$\{e_{ij}\}={\begin{pmatrix}e_{xx}&e_{xy}&e_{xz}\\e_{yx}&e_{yy}&e_{yz}\\e_{zx}&e_{zy}&e_{zz}\end{pmatrix}}.$

(L-6)

Note that, by way of equations (L-5a,L-5b,L-5c), the strain tensor (L-6) is symmetric.

Also from equations (L-4a,L-4b), we obtain by subtraction

$\zeta =\theta _{xz}=-{\frac {1}{2}}\left({\frac {\partial w}{\partial x}}-{\frac {\partial u}{\partial z}}\right)=-\theta _{zx}.$

(L-7a)

Similarly,

$\theta _{xy}=-{\frac {1}{2}}\left({\frac {\partial v}{\partial x}}-{\frac {\partial u}{\partial y}}\right)=-\theta _{yx}$

(L-7b)

and

$\theta _{yz}=-{\frac {1}{2}}\left({\frac {\partial w}{\partial y}}-{\frac {\partial v}{\partial z}}\right)=-\theta _{zy}.$

(L-7c)

Note that by definition

$\theta _{xx}=\theta _{yy}=\theta _{zz}=0.$

The angular deformations given by equations (L-5a,L-5b,L-5c) are called shear strains since they result in a shearing of the volume around a point within a solid body (Figure L-3b). The quantities given by equations (L-7a,L-7b,L-7c) are associated with rotation without deformation (Figure L-3c).

The displacement tensor in equation (L-2) can be expressed as

${\begin{array}{l}{\begin{pmatrix}\partial u/\partial x&\partial u/\partial y&\partial u/\partial z\\\partial v/\partial x&\partial v/\partial y&\partial v/\partial z\\\partial v/\partial x&\partial v/\partial y&\partial v/\partial z\end{pmatrix}}={\frac {1}{2}}{\begin{pmatrix}\partial u/\partial x+\partial u/\partial x&\partial v/\partial x+\partial u/\partial y&\partial w/\partial x+\partial u/\partial z\\\partial u/\partial y+\partial v/\partial x&\partial v/\partial y+\partial v/\partial y&\partial w/\partial y+\partial v/\partial z\\\partial u/\partial z+\partial w/\partial x&\partial v/\partial z+\partial w/\partial y&\partial w/\partial z+\partial w/\partial z\end{pmatrix}}\\\qquad \qquad \qquad \qquad \qquad \qquad \qquad -{\frac {1}{2}}{\begin{pmatrix}0&\partial v/\partial x-\partial u/\partial y&\partial w/\partial x-\partial u/\partial z\\\partial u/\partial y-\partial v/\partial x&0&\partial w/\partial y-\partial v/\partial z\\\partial u/\partial z-\partial w/\partial x&\partial v/\partial z-\partial w/\partial y&0\end{pmatrix}}.\end{array}}$

(L-8)

The elements of the matrices in equation (L-8) can then be expressed in terms of the normal strain components e_xx, e_yy, and e_zz, given by equations (L-3a,L-3b,L-3c), the shear strain components, e_xz, e_xy, and e_yz, given by equations (L-5a,L-5b,L-5c), and the rigid rotational components θ_xz, θ_xy, and θ_yz, given by equations (L-7a,L-7b,L-7c). All of these components are then substituted into equation (L-2) to get

${\begin{pmatrix}\delta u\\\delta v\\\delta w\end{pmatrix}}=\left[{\begin{pmatrix}e_{xx}&e_{xy}&e_{xz}\\e_{yx}&e_{yy}&e_{yz}\\e_{zx}&e_{zy}&e_{zz}\end{pmatrix}}+{\begin{pmatrix}0&\theta _{xy}&\theta _{xz}\\\theta _{yx}&0&\theta _{yz}\\\theta _{zx}&\theta _{zy}&0\end{pmatrix}}\right]{\begin{pmatrix}\delta x\\\delta y\\\delta z\end{pmatrix}}.$

(L-9)

Next, we need to establish the relation between the stress tensor (L-1) and the strain tensor (L-6). For an elastic solid, Hooke’s law states that the strain at any point is directly proportional to the stresses applied at that point. First, consider a volume (Figure L-1) under tensional principle stresses, P_xx, P_yy, and P_zz. The stress-strain relations are written as

$P_{xx}=Ee_{xx},$

(L-10a)

$-\sigma P_{xx}=Ee_{yy},$

(L-10b)

and

$-\sigma P_{xx}=Ee_{zz},$

(L-10c)

where E and σ are proportionality constants, which are specific to the material and are called Young’s modulus and Poisson’s ratio, respectively.

Consider a cylindrical rod that is subjected to a longitudinal extension in the x-direction. All other stresses being zero, this causes a lateral contraction in the y-direction. Note from equation (L-10a) that Young’s modulus is the ratio of the longitudinal stress P_xx to the longitudinal strain e_xx. Substitute P_xx from equation (L-10a) into equation (L-10b) and note that Poisson’s ratio is the ratio of the lateral contraction defined by the strain component −e_yy to longitudinal extension defined by the strain component e_xx. Since strain is a dimensionless quantity, Young’s modulus has the dimensions of stress, and Poisson’s ratio is a pure number.

Similar to equations (L-10a,L-10b,L-10c), we establish the following relations

$P_{yy}=Ee_{yy},$

(L-11a)

$-\sigma P_{yy}=Ee_{xx},$

(L-11b)

$-\sigma P_{yy}=Ee_{zz},$

(L-11c)

and

$P_{zz}=Ee_{zz},$

(L-12a)

$-\sigma P_{zz}=Ee_{xx},$

(L-12b)

$-\sigma P_{zz}=Ee_{yy}.$

(L-12c)

By combining equations (L-10), (L-11), and (L-12), we rewrite the principal stress-strain relations as

$P_{xx}-\sigma P_{yy}-\sigma P_{zz}=3Ee_{xx},$

(L-13a)

$-\sigma P_{xx}-P_{yy}-\sigma P_{zz}=3Ee_{yy},$

(L-13b)

$-\sigma P_{xx}-\sigma P_{yy}-P_{zz}=3Ee_{zz}.$

(L-13c)

By rewriting equations (L-13a,L-13b,L-13c), we get

$(1+\sigma )P_{xx}-\sigma (P_{xx}+P_{yy}+P_{zz})=Ee_{xx},$

(L-14a)

$(1+\sigma )P_{yy}-\sigma (P_{xx}+P_{yy}+P_{zz})=Ee_{yy},$

(L-14b)

$(1+\sigma )P_{zz}-\sigma (P_{xx}+P_{yy}+P_{zz})=Ee_{zz},$

(L-14c)

and adding them, we obtain

$(1-2\sigma )(P_{xx}+P_{yy}+P_{zz})=E(e_{xx}+e_{yy}+e_{zz}).$

(L-15)

Refer to Figure L-1. If the unstrained volume is (δxδyδz), and the strained volume is (δx + δu)(δy + δv)(δz + δw), then the fractional change in volume, Δ, is

$\Delta =\left[(\delta x+\delta u)(\delta y+\delta v)(\delta z+\delta w)-(\delta x\ \delta y\ \delta z)\right]/(\delta x\ \delta y\ \delta z).$

By neglecting the higher-order terms in this ratio, and referring to the definitions of the principal strain components given by equations (L-3a,L-3b,L-3c), we find that the fractional change in volume, or dilatation Δ, in the limiting case of the volume becoming infinitesimally small, is given by

$\Delta =e_{xx}+e_{yy}+e_{zz},$

(L-16)

which, upon substitution into equation (L-15), yields

$P_{xx}+P_{yy}+P_{zz}={\frac {E\Delta }{1-2\sigma }}.$

(L-17)

Finally, by using the relations (L-16) and (L-17), and back substituting into equations (L-13a,L-13b,L-13c), we obtain the relations between the principal stress and principal strain components:

$P_{xx}=\lambda \Delta +2\mu e_{xx},$

(L-18a)

$P_{yy}=\lambda \Delta +2\mu e_{yy},$

(L-18b)

$P_{zz}=\lambda \Delta +2\mu e_{zz},$

(L-18c)

where λ and μ are the elastic moduli for the solid given by

$\lambda ={\frac {\sigma E}{(1+\sigma )(1-2\sigma )}}$

(L-19a)

and

$\mu ={\frac {E}{2(1+\sigma )}}.$

(L-19b)

These elastic moduli have the dimensions of stress.

Now consider the volume in Figure L-1 under shear stress components, P_xy, P_xz, and P_yz. The deformations associated with these stress components are the shear strains, e_xy, e_xz, and e_yz. These stress and strain components also are linearly related for elastic solids:

$P_{xy}=2\mu e_{xy},$

(L-20a)

$P_{xz}=2\mu e_{xz},$

(L-20b)

$P_{yz}=2\mu e_{yz},$

(L-20c)

where the proportionality constant μ is known as the modulus of rigidity. Note, from equation (L-20a), that modulus of rigidity is the ratio of shear stress to shear strain.

By combining equations (L-18a,L-18b,L-18c) and (L-20a,L-20b,L-20c), we obtain the stress-strain relation for elastic solids:

${\begin{pmatrix}P_{xx}&P_{xy}&P_{xz}\\P_{xy}&P_{yy}&P_{yz}\\P_{zx}&P_{zy}&P_{zz}\end{pmatrix}}=\lambda \Delta {\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}}+2\mu {\begin{pmatrix}e_{xx}&e_{xy}&e_{xz}\\e_{xy}&e_{yy}&e_{yz}\\e_{zx}&e_{zy}&e_{zz}\end{pmatrix}}.$

(L-21)

This is the formal expression for Hooke’s law that relates the stress tensor on the left to the strain tensor on the right. Recall that these two tensors are symmetric, and that the diagonal elements represent the normal components and the off-diagonal elements represent the shear components. Equation (L-21) holds for homogeneous, isotropic, elastic solids (for which elastic behavior is independent of direction), and for deformations that are sufficiently small (usual case for seismic waves) to satisfy the linear relationship between stress and strain.

Since the stress tensor (L-1) and the strain tensor (L-6) are symmetric, equation (L-21) can be written in the form

${\begin{pmatrix}P_{xx}\\P_{yy}\\P_{zz}\\P_{xy}\\P_{xz}\\P_{yz}\end{pmatrix}}={\begin{pmatrix}\lambda +2\mu &\lambda &\lambda &0&0&0\\\lambda &\lambda +2\mu &\lambda &0&0&0\\\lambda &\lambda &\lambda +2\mu &0&0&0\\0&0&0&2\mu &0&0\\0&0&0&0&2\mu &0\\0&0&0&0&0&2\mu \\\end{pmatrix}}{\begin{pmatrix}e_{xx}\\e_{yy}\\e_{zz}\\e_{xy}\\e_{xz}\\e_{yz}\end{pmatrix}}.$

(L-22)

This equation is a special form of the generalized Hooke’s law given by

${\begin{pmatrix}P_{xx}\\P_{yy}\\P_{zz}\\P_{xy}\\P_{xz}\\P_{yz}\end{pmatrix}}={\begin{pmatrix}c_{11}&c_{12}&c_{13}&c_{14}&c_{15}&c_{16}\\c_{21}&c_{22}&c_{23}&c_{24}&c_{25}&c_{26}\\c_{31}&c_{32}&c_{33}&c_{34}&c_{35}&c_{36}\\c_{41}&c_{42}&c_{43}&c_{44}&c_{45}&c_{46}\\c_{51}&c_{52}&c_{53}&c_{54}&c_{55}&c_{56}\\c_{61}&c_{62}&c_{63}&c_{64}&c_{65}&c_{66}\\\end{pmatrix}}{\begin{pmatrix}e_{xx}\\e_{yy}\\e_{zz}\\e_{xy}\\e_{xz}\\e_{yz}\end{pmatrix}},$

(L-23)

where c_ij = c_ji are the 21 independent constants for an elastic medium. The stress vector on the left is related to the strain vector on the right by the stiffness matrix {c_ij}. Equation (L-23) states that a stress component is a linear combination of all the strain components. This relation is the foundation of the linear theory of elasticity and its physical basis is the assumption that elastic deformations in solids are infinitesimally small.

For an isotropic solid, the number of independent constants is reduced to two — Lamé’s constants, λ and μ given by equations (L-19a,L-19b); hence equation (L-23) reduces to the special form of equation (L-22). For a transversely isotropic solid, for which elastic behavior is the same in two orthogonal directions but different in the third direction, the number of independent constants is five .

L.2 Elastic wave equation

We now review the equations of motion and examine how stress fields are propagated through an elastic solid. Refer to Figure L-1 and note that, for small displacements and velocities, and neglecting body forces such as gravity acting on the volume surrounding the points, Newton’s second law of motion gives, for the displacement component u,

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}={\frac {\partial P_{xx}}{\partial x}}+{\frac {\partial P_{xy}}{\partial y}}+{\frac {\partial P_{xz}}{\partial z}},$

(L-24)

where ρ is the density. From the expression (L-21) for Hooke’s law, we have

$P_{xx}=\lambda \Delta +2\mu e_{xx},$

(L-25a)

$P_{xy}=2\mu e_{xy},$

(L-25b)

and

$P_{xz}=2\mu e_{xz}.$

(L-25c)

By making these substitutions into equation (L-24), and assuming that the elastic parameters λ and μ are constant, we obtain

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}=\lambda {\frac {\partial \Delta }{\partial x}}+2\mu {\frac {\partial e_{xx}}{\partial x}}+2\mu {\frac {\partial e_{xy}}{\partial y}}+2\mu {\frac {\partial e_{xz}}{\partial z}}.$

(L-26)

By using the definitions for the principal strains given by equations (L-3a,L-3b,L-3c) and that for the dilatation given by equation (L-16), the latter can be expressed in terms of the displacement components as

$\Delta ={\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}.$

(L-27)

Next, by substituting equation (L-27) and the relations for the strain components in terms of the displacement components given by equations (L-3a,L-3b,L-3c) and (L-5a,L-5b,L-5c) into equation (L-26), we obtain

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}=\lambda {\frac {\partial }{\partial x}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)+2\mu {\frac {\partial ^{2}u}{\partial x^{2}}}+\mu \left({\frac {\partial ^{2}v}{\partial y\partial x}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)+\mu \left({\frac {\partial ^{2}w}{\partial z\partial x}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\right).$

(L-28a)

Rearrange the terms on the right-hand side to get the following expression

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}=(\lambda +\mu ){\frac {\partial }{\partial x}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)+\mu \left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\right).$

(L-28b)

Finally, back substitution of equation (L-27) and the definition for the Laplacian operator, ∇² : (∂²/∂x² + ∂²/∂y² + ∂²/∂z²) into equation (L-28b), yields the equation of motion for the displacement component u:

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}=(\lambda +\mu ){\frac {\partial \Delta }{\partial x}}+\mu \nabla ^{2}u.$

(L-29a)

Similarly, we can derive the equations of motion for the displacement components v and w

$\rho {\frac {\partial ^{2}v}{\partial t^{2}}}=(\lambda +\mu ){\frac {\partial \Delta }{\partial y}}+\mu \nabla ^{2}v$

(L-29b)

and

$\rho {\frac {\partial ^{2}w}{\partial t^{2}}}=(\lambda +\mu ){\frac {\partial \Delta }{\partial z}}+\mu \nabla ^{2}w.$

(L-29c)

Define the displacement vector u :(u, v, w) and combine the three equations (L-29a,L-29b,L-29c) to get

$\rho {\frac {\partial ^{2}{\mathbf {u} }}{\partial t^{2}}}=(\lambda +\mu )\nabla \Delta +\mu \nabla ^{2}{\mathbf {u} }.$

(L-30)

This is the equation of wave propagation in homogeneous, isotropic, and elastic solids.

L.3 Seismic wave types — body waves and surface waves

Equation (L-30) can be specialized to describe various wave types that travel within solids and fluids (body waves), and along free surfaces and layer boundaries (surface waves). We shall derive the equations for two types of body waves that are of interest in exploration seismology — the compressional wave for which the displacement is in the direction of propagation and the shear wave for which the displacement is in the direction perpendicular to the direction of propagation. Surface waves are more of an interest in earthquake seismology; as such, we shall not deal with them here.

To begin with, differentiate equation (L-29a) with respect to x

$\rho {\frac {\partial }{\partial x}}\left({\frac {\partial ^{2}u}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial x^{2}}}+\mu {\frac {\partial }{\partial x}}\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\right),$

(L-31a)

equation (L-29b) with respect to y

$\rho {\frac {\partial }{\partial y}}\left({\frac {\partial ^{2}v}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial y^{2}}}+\mu {\frac {\partial }{\partial y}}\left({\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial ^{2}v}{\partial y^{2}}}+{\frac {\partial ^{2}v}{\partial z^{2}}}\right),$

(L-31b)

and equation (L-29c) with respect to z

$\rho {\frac {\partial }{\partial z}}\left({\frac {\partial ^{2}w}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial z^{2}}}+\mu {\frac {\partial }{\partial z}}\left({\frac {\partial ^{2}w}{\partial x^{2}}}+{\frac {\partial ^{2}w}{\partial y^{2}}}+{\frac {\partial ^{2}w}{\partial z^{2}}}\right).$

(L-31c)

Next, reorder the differentiation on the left-hand side of equations (L-31a,L-31b,L-31c) and add them to get

${\begin{array}{l}\rho {\frac {\partial ^{2}}{\partial t^{2}}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)=(\lambda +\mu )\left({\frac {\partial ^{2}\Delta }{\partial x^{2}}}+{\frac {\partial ^{2}\Delta }{\partial y^{2}}}+{\frac {\partial ^{2}\Delta }{\partial z^{2}}}\right)\\\qquad \qquad \quad \quad +\mu {\frac {\partial }{\partial x}}\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\right)+\mu {\frac {\partial }{\partial y}}\left({\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial ^{2}v}{\partial y^{2}}}+{\frac {\partial ^{2}v}{\partial z^{2}}}\right)+\mu {\frac {\partial }{\partial z}}\left({\frac {\partial ^{2}w}{\partial x^{2}}}+{\frac {\partial ^{2}w}{\partial y^{2}}}+{\frac {\partial ^{2}w}{\partial z^{2}}}\right).\end{array}}$

(L-32)

Then, change the order of differentiation and rearrange the terms on the right-hand side to get

${\begin{array}{l}\rho {\frac {\partial ^{2}}{\partial t^{2}}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)=(\lambda +\mu )\left({\frac {\partial ^{2}\Delta }{\partial x^{2}}}+{\frac {\partial ^{2}\Delta }{\partial y^{2}}}+{\frac {\partial ^{2}\Delta }{\partial z^{2}}}\right)\\\qquad \qquad \quad \quad +\mu {\frac {\partial ^{2}}{\partial x^{2}}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)+\mu {\frac {\partial ^{2}}{\partial y^{2}}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right)+\mu {\frac {\partial ^{2}}{\partial z^{2}}}\left({\frac {\partial u}{\partial x}}+{\frac {\partial v}{\partial y}}+{\frac {\partial w}{\partial z}}\right).\end{array}}$

(L-33)

Finally, substitution of equation (L-27) and the definition for the Laplacian operator, ∇² : (∂²/∂x² + ∂²/∂y² + ∂₂/∂z²) into the above expression, yield the equation for dilatational or compressional-wave propagation:

$\rho {\frac {\partial ^{2}\Delta }{\partial t^{2}}}=(\lambda +2\mu )\nabla ^{2}\Delta .$

(L-34)

The velocity with which compressional waves propagate in an elastic solid then is

$\alpha ={\sqrt {\frac {\lambda +2\mu }{\rho }}}.$

(L-35)

Equation (L-34) also can be derived by using the more compact vector algebra. Specifically, by taking the divergence of both sides of equation (L-30) and reordering the differentiation, we have

$\rho {\frac {\partial ^{2}}{\partial t^{2}}}(\nabla \cdot {\mathbf {u} })=(\lambda +\mu )\nabla \cdot (\nabla \Delta )+\mu \nabla ^{2}(\nabla \cdot {\mathbf {u} }).$

(L-36)

By referring to equation (L-27), note that the dilatation is the divergence of the displacement vector: $\Delta =\nabla \cdot \mathbf {u} .$ Substitution of this relation into equation (L-36) and observing that the divergence of a gradient is the Laplacian, $\nabla \cdot (\nabla \Delta )=\nabla ^{2}\Delta$ , yield equation (L-34) that describes the compressional-wave propagation.

Differentiate equation (L-29b) with respect to z

$\rho {\frac {\partial }{\partial z}}\left({\frac {\partial ^{2}v}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial z\partial y}}+\mu {\frac {\partial }{\partial z}}(\nabla ^{2}v),$

(L-37a)

and equation (L-29c) with respect to y

$\rho {\frac {\partial }{\partial y}}\left({\frac {\partial ^{2}w}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial y\partial z}}+\mu {\frac {\partial }{\partial y}}(\nabla ^{2}w).$

(L-37b)

Next, change the order of differentiation and subtract equation (L-37a) from (L-37b)

$\rho {\frac {\partial ^{2}}{\partial t^{2}}}\left({\frac {\partial w}{\partial y}}-{\frac {\partial v}{\partial z}}\right)=\mu \nabla ^{2}\left({\frac {\partial w}{\partial y}}-{\frac {\partial v}{\partial z}}\right).$

(L-38)

Finally, substitute equation (L-7c) to get

$\rho {\frac {\partial ^{2}\theta _{yz}}{\partial t^{2}}}=\mu \nabla ^{2}\theta _{yz}.$

(L-39)

Similarly, differentiate equation (L-29a) with respect to z

$\rho {\frac {\partial }{\partial z}}\left({\frac {\partial ^{2}u}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial z\partial x}}+\mu {\frac {\partial }{\partial z}}(\nabla ^{2}u),$

(L-40a)

and equation (L-29c) with respect to x

$\rho {\frac {\partial }{\partial x}}\left({\frac {\partial ^{2}w}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial x\partial z}}+\mu {\frac {\partial }{\partial x}}(\nabla ^{2}w).$

(L-40b)

Next, change the order of differentiation and subtract equation (L-40a) from (L-40b)

$\rho {\frac {\partial ^{2}}{\partial t^{2}}}\left({\frac {\partial w}{\partial x}}-{\frac {\partial u}{\partial z}}\right)=\mu \nabla ^{2}\left({\frac {\partial w}{\partial x}}-{\frac {\partial u}{\partial z}}\right),$

(L-41)

and substitute equation (L-7a) to get

$\rho {\frac {\partial ^{2}\theta _{xz}}{\partial t^{2}}}=\mu \nabla ^{2}\theta _{xz}.$

(L-42)

Finally, differentiate equation (L-29a) with respect to y

$\rho {\frac {\partial }{\partial y}}\left({\frac {\partial ^{2}u}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial y\partial x}}+\mu {\frac {\partial }{\partial y}}(\nabla ^{2}u),$

(L-43a)

and equation (L-29b) with respect to x

$\rho {\frac {\partial }{\partial x}}\left({\frac {\partial ^{2}v}{\partial t^{2}}}\right)=(\lambda +\mu ){\frac {\partial ^{2}\Delta }{\partial x\partial y}}+\mu {\frac {\partial }{\partial x}}(\nabla ^{2}v).$

(L-43b)

Next, change the order of differentiation and subtract equation (L-43a) from (L-43b)

$\rho {\frac {\partial ^{2}}{\partial t^{2}}}\left({\frac {\partial v}{\partial x}}-{\frac {\partial u}{\partial y}}\right)=\mu \nabla ^{2}\left({\frac {\partial v}{\partial x}}-{\frac {\partial u}{\partial y}}\right),$

(L-44)

and substitute equation (L-7b) to get

$\rho {\frac {\partial ^{2}\theta _{xy}}{\partial t^{2}}}=\mu \nabla ^{2}\theta _{xy}.$

(L-45)

Define the rotational vector Θ : (θ_xy, θ_xz, θ_yz), and combine the three equations (L-39), (L-42), and (L-45) to get the equation for rotational or shear-wave propagation:

$\rho {\frac {\partial ^{2}\Theta }{\partial t^{2}}}=\mu \nabla ^{2}\Theta .$

(L-46)

The velocity with which shear waves propagate in an elastic solid then is

$\beta ={\sqrt {\frac {\mu }{\rho }}}.$

(L-47)

Equation (L-46) also can be derived by using the more compact vector algebra. Specifically, by taking the curl of both sides of equation (L-30) and reordering differentiation, we have

$\rho {\frac {\partial ^{2}}{\partial t^{2}}}(\nabla \times {\mathbf {u} })=(\lambda +\mu )\nabla \times \nabla \Delta +\mu \nabla ^{2}(\nabla \times {\mathbf {u} }).$

(L-48)

By referring to equations (L-7a,L-7b,L-7c) note that the rotational vector Θ : (θ_xy, θ_xz, θ_yz) is the curl of the displacement vector: Θ = (1/2)∇ × u. Substitution of this relation into equation (L-48) and observing that the curl of a gradient is zero, $\nabla \times \nabla \Delta =0$ , yield equation (L-46) for rotational or shear-wave propagation.

By taking the ratio of equations (L-35) and (L-47), and using the relations (L-19a,L-19b), we obtain the expression for the ratio of the shear-wave velocity β to the compressional-wave velocity α

${\frac {\beta }{\alpha }}={\sqrt {\frac {1-2\sigma }{2(1-\sigma )}}},$

(L-49)

which indicates that the velocity ratio β/α, depends on Poisson’s ratio σ, and that the compressional-wave velocity α is always greater than the shear-wave velocity β.

The wave type described by equation (L-34) is commonly referred to as the P-wave (or equivalently, longitudinal wave, dilatational wave, or compressional wave), and the wave type described by equation (L-46) is commonly referred to as the S-wave (or equivalently, rotational wave, transverse wave, or shear wave). For P-waves, the particle motion is in the direction of wave propagation, whereas for S-waves, the particle motion is in the direction perpendicular to the direction of wave propagation.

To understand the nature of P- and S-wave propagation, consider the simple case of plane-wave propagation in the x-direction. This makes the displacement vector u : (u, v, w) a function of x only. Refer to equations (L-29a,L-29b,L-29c), substitute for the cubical dilatation from equation (L-27), and finally, retain only the derivatives with respect to x to obtain

$\rho {\frac {\partial ^{2}u}{\partial t^{2}}}=(\lambda +2\mu ){\frac {\partial ^{2}u}{\partial x^{2}}},$

(L-50a)

$\rho {\frac {\partial ^{2}v}{\partial t^{2}}}=\mu {\frac {\partial ^{2}v}{\partial x^{2}}},$

(L-50b)

and

$\rho {\frac {\partial ^{2}w}{\partial t^{2}}}=\mu {\frac {\partial ^{2}w}{\partial x^{2}}}.$

(L-50c)

From equation (L-50a), note that the particle displacement u (the longitudinal vibration) which is in the direction of propagation (in this case in the x-direction) travels with the velocity of compressional waves (equation L-35). From equations (L-50b,L-50c), the particle displacements v and w (the transverse vibrations), which are in the directions perpendicular to the direction of propagation, travel with the velocity of shear waves (equation L-47).

To a great extent, seismic waves traveling in the earth can be considered as elastic waves. Therefore, we can use the wave equation (L-30) for elastic solids to describe seismic wave propagation. Source types used in seismic prospecting often generate P-waves. Furthermore, in marine exploration, because of the water layer surrounding the source, S-waves cannot propagate. Nevertheless, at non-normal incidence, part of the P-wave energy is converted to S-waves. The P-to-S conversion phenomenon is the basis for the 4-C seismic method (4-C seismic method).

Equations for elastic wave propagation can be specialized to describe acoustic wave propagation induced by hydrostatic stress. Hooke’s law given by equation (L-21) takes the form

$P_{xx}=\lambda \Delta +2\mu e_{xx},$

(L-51a)

$P_{yy}=\lambda \Delta +2\mu e_{yy},$

(L-51b)

and

$P_{zz}=\lambda \Delta +2\mu e_{zz}.$

(L-51c)

Also, note that the normal stress components are equal to the hydrostatic pressure, P_xx = P_yy = P_zz = −P, and that shear stress components vanish, P_xy = P_xz = P_yz = 0. By summing the three equations (L-51a,L-51b,L-51c) and substituting for the dilatation from equation (L-16), we obtain

$\kappa =\lambda +{\frac {2}{3}}\mu ,$

(L-52a)

where κ = −P/Δ is the ratio of the hydrostatic stress P to the volumetric strain Δ and is called the incompressibility or bulk modulus.

Solve equation (L-52a) for λ and substitute into equation (L-35) to get the expression for the compressional-wave velocity in terms of bulk modulus κ and modulus of rigidity μ

$\alpha ={\sqrt {\frac {\kappa +{\frac {4}{3}}\mu }{\rho }}}.$

(L-52b)

For fluids, modulus of rigidity is zero. Hence, by setting μ = 0, κ = λ and using the relation (L-52), for acoustic waves traveling in fluid media, equation (L-34) takes the form

$\rho {\frac {\partial ^{2}P}{\partial t^{2}}}=\kappa \nabla ^{2}P,$

(L-53)

where P represents the pressure wavefield, which travels with a speed ${\sqrt {\kappa \!/\!\rho }}$.

Since the P- and S-waves travel within an elastic solid, they are called body waves. When an elastic solid is bounded, it can also support waves traveling along the boundary, such as the free surface. Waves that travel along the surface are called surface waves. Surface waves are generally characterized as dispersive; that is, each frequency component travels at a different speed.

There are two common types of surface waves: Rayleigh waves and Love waves. Rayleigh waves travel along the free surface and give rise to particle motion that is always in a vertical plane, and elliptical and retrograde with respect to the direction of propagation. Although the wave motion attenuates rapidly with depth, Rayleigh waves attenuate less rapidly with distance along the direction of propagation because of their two-dimensional (2-D) character, compared to P- and S-waves that have three-dimensional (3-D) character. The velocity with which Rayleigh waves travel is a little less than that of shear waves in the same medium.

In many seismograms recorded on land where the near-surface comprises a low-velocity layer, Rayleigh waves are observed in the form of a low-velocity, large-amplitude, and low-frequency wavetrain. These waves are commonly known as ground roll, and they obscure the useful reflection energy on recorded seismograms.

Love waves are surface waves with particle motion that is horizontal and transverse to the direction of propagation. They occur when there is a low-velocity layer overlying a high-velocity medium. Because the corresponding particle motion is horizontal, Love waves are usually not recorded in seismic prospecting by receivers that respond only to vertical displacements.

L.4 Wave propagation phenomena — diffraction, reflection, and refraction

A seismic energy source generates a spherical wave in a 3-D earth. Following an explosion, a spherical cavity is created with its periphery forming a zone of permanent deformation. However, a short distance away from the cavity, the seismic energy induces an elastic deformation. The particle motion associated with this deformation can be described by a time-varying function, commonly known as the source waveform.

It turns out that every point on the spherical wavefront acts like a secondary source that generates its own spherical wave. This is known as Huygens’ principle which is the physical basis for wave propagation. Huygens’ principle can be used to conveniently describe reflection and refraction at layer boundaries, and diffraction from sharp discontinuities in the subsurface.

At sufficiently large distances away from the source location, the spherical wavefront can be represented by plane waves traveling at different angles from the vertical. Each plane wave carries along the propagation path the seismic energy contained in the source waveform with a certain frequency band. The propagation path is known as the raypath which is perpendicular to the wavefront in isotropic media.

As the spherical wave propagates outward and away from its source, two things happen to the source waveform. First, the energy generated by the source is spread over the area of the spherical wavefront, which increases as the square of the radius of the wavefront. Hence, the energy per unit area is inversely proportional to the square of the radius, while the wave amplitude is inversely proportional to the radius. The phenomenon of wave amplitude decaying with distance from the source is called spherical spreading.

Second, part of the seismic energy is lost along the path of propagation because of intrinsic attenuation of rocks. This attenuation is in the form of frictional dissipation of energy into heat. Attenuation mechanisms are still subject to research; nevertheless, a probable cause of loss of seismic energy from absorption is widely believed to be the presence of pore fluids. When set to motion by the passing wave, fluids take out a significant part of the seismic energy. Wave attenuation is generally described by an amplitude function that decays exponentially with distance from the source. The rate of decay depends on the rock type and is described by a constant factor. The constant fractional energy loss per cycle of seismic waves is highest for weathered rocks (3 dB/wavelength), and it is about 0.5 dB/wavelength for normal rock types. Since the energy loss is considered to be constant for a given wavelength, higher frequencies are attenuated more rapidly than lower frequencies. As the source waveform travels down in the earth, it gets broader along the propagation path because of gradual loss of higher frequencies. This means that the earth behaves as a low-pass filter.

Another aspect of wave attenuation is related to the observation made about many of the earth materials that do not behave as perfectly elastic when subjected to extended loading. This nonelastic behavior is called viscoelasticity , and it requires some modifications to Hooke’s law given by equation (L-21). Specifically, stress components are related to strain components and the rate of change of strain components themselves. When viscosity is introduced into the wave propagation in solids, we can still identify P- and S-wave types. However, these waves are now damped and dispersive — they are subjected to frequency-dependent attenuation. It is generally believed that the dispersive nature attributable to viscoelastic wave propagation is negligible within the range of frequencies used in seismic prospecting.

Seismic waves — whether compressional or shear, are subjected to reflection and refraction at layer boundaries with impedance contrast and diffraction at a sharp discontinuity. In this section, we shall briefly review diffraction phenomenon and discuss reflection and refraction of seismic waves within the context of amplitude variation with offset (AVO) analysis and amplitude inversion.

When the incident plane wave encounters a sharp discontinuity with the radius of curvature much less than the seismic wavelength, such as a fault edge, then the energy propagates outward from that point acting as a Huygens’ secondary source in the form of a spherical wave. This phenomenon is known as diffraction.

The amplitude of the diffracted wave along the hyperbolic traveltime trajectory decays with distance away from the diffractor. For an incident plane compressional plane wave, Sommerfeld gives the exact solution for the diffracted wave amplitude as

$A(\omega ,x,z)=A_{1}\exp {\big [}i({\frac {\omega }{\alpha }}r+{\frac {\pi }{4}}){\big ]},$

(L-54a)

where

$A_{1}=-A_{0}{\sqrt {8\pi {\frac {\omega }{\alpha }}}}{\big [}{\rm {sec}}{\frac {1}{2}}(\psi -\varphi )-{\rm {csc}}{\frac {1}{2}}(\psi +\varphi ){\big ]},$

(L-54b)

and A₀ is the incident wave amplitude, α is the wave velocity, ψ is the angle of incidence, $r={\sqrt {x^{2}+z^{2}}},\ \varphi$ is the angle at which energy diffracts, and ω is the angular frequency. Note that the solution given by equation (L-54a) is for 2-D geometry and for sufficiently large distances, thus the term far-field solution.

It is interesting to note that a reflector can be represented by a series of discrete point diffractors. The response of a series of diffractors is simply the superposition of the individual responses. In the limit, we obtain the traveltime response associated with the reflector which is made up of diffractors with an infinitesimal separation between them.

Figure L-4  A two-layered 2-D earth model used in deriving the Zoeppritz equations (L-78).

We now review the reflection and refraction phenomena. For simplicity, consider a monochromatic compressional plane wave that impinges at normal incidence upon a flat layer boundary at z = 0. The incident energy is partitioned between a reflected and transmitted compressional plane wave. This special case is discussed in analysis of amplitude variation with offset. In the next section, we shall treat the general case of a compressional plane wave that is incident at an oblique angle to the interface. Derivation of the reflected and transmitted wave functions gets complicated, and we find that the reflection coefficient changes with angle of incidence. Moreover, at non-normal incidence, the incident compressional-wave energy is partitioned at the interface into four components — reflected compressional, reflected shear, transmitted compressional, and transmitted shear waves.

L.5 The Zoeppritz equations

Consider the case of a 2-D medium in (x, z) coordinates that comprises two solid layers separated by a flat interface as shown in Figure L-4. We shall adopt the derivation of reflected and transmitted wave amplitudes by Officer for this 2-D case between an acoustic layer (water layer) on top and a solid layer (elastic half space) below to derive the Zoeppritz equations that describe the reflected and refracted P- and S-wave amplitudes for the case of two solid layers, again separated by a flat interface at z = 0.

In many of the wave propagation problems, it is convenient to work with displacement potentials in lieu of the displacement vector u : (u, w) itself . We shall define two displacement potentials, a compressional-wave potential Φ and a shear-wave potential Ψ in terms of the displacements u and w as

$u={\frac {\partial \Phi }{\partial x}}+{\frac {\partial \Psi }{\partial z}}$

(L-55a)

and

$w={\frac {\partial \Phi }{\partial z}}-{\frac {\partial \Psi }{\partial x}}.$

(L-55b)

Choice of the displacement potentials is made on the basis that they satisfy the wave equation while enabling us to describe the displacement vector u in terms of a compressional component and a rotational component. To see that the two potentials defined by equations (L-55a,L-55b) satisfy the wave equation, proceed as follows.

Adapt the dilatation Δ given by equation (L-27) to the geometry of Figure L-4 in (x, z) coordinates

$\Delta ={\frac {\partial u}{\partial x}}+{\frac {\partial w}{\partial z}},$

(L-56)

and substitute equations (L-55a,L-55b)

$\Delta ={\frac {\partial }{\partial x}}\left({\frac {\partial \Phi }{\partial x}}+{\frac {\partial \Psi }{\partial z}}\right)+{\frac {\partial }{\partial z}}\left({\frac {\partial \Phi }{\partial z}}-{\frac {\partial \Psi }{\partial x}}\right).$

(L-57)

Then, simplify to get the relation

$\Delta ={\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial z^{2}}}.$

(L-58)

Finally, substitute equation (L-58) into the compressional-wave equation (L-34) to obtain

${\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial z^{2}}}={\frac {1}{\alpha ^{2}}}{\frac {\partial ^{2}\Phi }{\partial t^{2}}},$

(L-59)

where the compressional-wave velocity α is given by equation (L-35). Equation (L-59) indicates that the displacement potential Φ obeys the wave equation.

Now consider the rotational shear component θ_xz given by equation (L-7a). This is the shear strain that is relevant to the wave motion in (x, z) coordinates of Figure L-4. Substitute equations (L-55a,L-55b) into equation (L-7a)

$\theta _{xz}=-{\frac {1}{2}}\left[{\frac {\partial }{\partial x}}\left({\frac {\partial \Phi }{\partial z}}-{\frac {\partial \Psi }{\partial x}}\right)-{\frac {\partial }{\partial z}}\left({\frac {\partial \Phi }{\partial x}}+{\frac {\partial \Psi }{\partial z}}\right)\right],$

(L-60)

and simplify to obtain

$\theta _{xz}={\frac {1}{2}}\left({\frac {\partial ^{2}\Psi }{\partial x^{2}}}+{\frac {\partial ^{2}\Psi }{\partial z^{2}}}\right).$

(L-61)

Finally, substitute equation (L-61) into the shear-wave equation (L-42) to obtain

${\frac {\partial ^{2}\Psi }{\partial x^{2}}}+{\frac {\partial ^{2}\Psi }{\partial z^{2}}}={\frac {1}{\beta ^{2}}}{\frac {\partial ^{2}\Psi }{\partial t^{2}}},$

(L-62)

where the shear-wave velocity β is given by equation (L-47). Equation (L-62) indicates that the displacement potential Ψ also obeys the wave equation.

The plane-wave solutions Φ(x, z, t) and Ψ(x, z, t) to the wave equations (L-59) and (L-62) in (x, z) coordinates of Figure L-4 are expressed as follows:

${\begin{array}{l}\Phi _{1}=A_{0}\exp \left(i{\frac {\omega }{\alpha _{1}}}\sin \varphi _{0}\,x+i{\frac {\omega }{\alpha _{1}}}\cos \varphi _{0}\,z-i\omega t\right)\\\quad \,\,\,+A_{1}\exp \left(i{\frac {\omega }{\alpha _{1}}}\sin \varphi _{1}\,x-i{\frac {\omega }{\alpha _{1}}}\cos \varphi _{1}\,z-i\omega t\right),\end{array}}$

(L-63a)

$\Psi _{1}=B_{1}\exp \left(i{\frac {\omega }{\beta _{1}}}\sin \psi _{1}\,x-i{\frac {\omega }{\beta _{1}}}\cos \psi _{1}\,z-i\omega t\right),$

(L-63b)

$\Phi _{2}=A_{2}\exp \left(i{\frac {\omega }{\alpha _{2}}}\sin \varphi _{2}\,x+i{\frac {\omega }{\alpha _{2}}}\cos \varphi _{2}\,z-i\omega t\right),$

(L-63c)

$\Psi _{2}=B_{2}\exp \left(i{\frac {\omega }{\beta _{2}}}\sin \psi _{2}\,x+i{\frac {\omega }{\beta _{2}}}\cos \psi _{2}\,z-i\omega t\right),$

(L-63d)

where A₀, A₁, B₁, A₂, and B₂ are the incident P-wave, reflected P-wave, reflected S-wave, transmitted P-wave, and transmitted S-wave amplitudes, α₁ and β₁ are the P- and S-wave velocities in layer 1 (top layer), and α₂ and β₂ are the P- and S-wave velocities in layer 2 (bottom layer). The angles φ₀, φ₁, φ₂, ψ₁, and ψ₂ are denoted in Figure L-4. Finally, ω is the angular frequency, and the z variable is positive in the downward direction.

Note that the displacement potential Φ₁ given by equation (L-63a) has two components associated with the incident P-wave and the reflected P-wave. To prove that the wave functions given by equations (L-63a,63b,63c,63d) are solutions to the wave equation, simply substitute these functions into the respective equations (L-59) and (L-62).

Given the incident P-wave amplitude A₀, our objective is to compute the amplitudes of the partitioned wave components — A₁, B₁, A₂, and B₂. We will require the following four boundary conditions to be satisfied at the layer boundary z = 0:

1. The displacement component tangential to the interface is continuous: u₁ = u₂.
2. The displacement component normal to the interface is continuous: w₁ = w₂.
3. The stress component normal to the interface is continuous: (P_zz)₁ = (P_zz)₂.
4. The stress component tangential to the interface is continuous: (P_xz)₁ = (P_xz)₂.

By way of Snell’s law that relates the propagation angles

${\frac {\sin \varphi _{0}}{\alpha _{1}}}={\frac {\sin \varphi _{1}}{\alpha _{1}}}={\frac {\sin \varphi _{2}}{\alpha _{2}}}={\frac {\sin \psi _{1}}{\beta _{1}}}={\frac {\sin \psi _{2}}{\beta _{2}}}.$

(L-64)

Note that the phase of the wave functions given by equations (L-63a,63b,63c,63d) coincide at the layer boundary z = 0. With this observation in mind, refer to equation (L-55a) and consider the first boundary condition expressed as

${\frac {\partial \Phi _{1}}{\partial x}}+{\frac {\partial \Psi _{1}}{\partial z}}={\frac {\partial \Phi _{2}}{\partial x}}+{\frac {\partial \Psi _{2}}{\partial z}}.$

(L-65)

Differentiate the respective wave functions given by equations (L-63a,63b,63c,63d) and substitute into equation (L-65), while noting from Snell’s law that φ₀ = φ₁

${\frac {\sin \varphi _{1}}{\alpha _{1}}}A_{0}+{\frac {\sin \varphi _{1}}{\alpha _{1}}}A_{1}-{\frac {\cos \psi _{1}}{\beta _{1}}}B_{1}={\frac {\sin \varphi _{2}}{\alpha _{2}}}A_{2}+{\frac {\cos \psi _{2}}{\beta _{2}}}B_{2}.$

(L-66)

Next, refer to equation (L-55b) and consider the second boundary condition expressed as

${\frac {\partial \Phi _{1}}{\partial z}}-{\frac {\partial \Psi _{1}}{\partial x}}={\frac {\partial \Phi _{2}}{\partial z}}-{\frac {\partial \Psi _{2}}{\partial x}}.$

(L-67)

Differentiate the respective wave functions given by equations (L-63a,63b,63c,63d) and substitute into equation (L-67)

${\frac {\cos \varphi _{1}}{\alpha _{1}}}A_{0}-{\frac {\cos \varphi _{1}}{\alpha _{1}}}A_{1}-{\frac {\sin \psi _{1}}{\beta _{1}}}B_{1}={\frac {\cos \varphi _{2}}{\alpha _{2}}}A_{2}-{\frac {\sin \psi _{2}}{\beta _{2}}}B_{2}.$

(L-68)

Next, refer to equation (L-18c) rewritten below

$P_{zz}=\lambda \Delta +2\mu e_{zz}$

and substitute equation (L-56) for Δ and equation (L-3c) for e_zz to get

$P_{zz}=\lambda \left({\frac {\partial u}{\partial x}}+{\frac {\partial w}{\partial z}}\right)+2\mu \left({\frac {\partial w}{\partial z}}\right).$

(L-69a)

Combine the terms

$P_{zz}=(\lambda +2\mu ){\frac {\partial w}{\partial z}}+\lambda {\frac {\partial u}{\partial x}},$

(L-69b)

and substitute equations (L-55a,L-55b) for the displacement potentials

$P_{zz}=(\lambda +2\mu ){\frac {\partial }{\partial z}}\left({\frac {\partial \Phi }{\partial z}}-{\frac {\partial \Psi }{\partial x}}\right)+\lambda {\frac {\partial }{\partial x}}\left({\frac {\partial \Phi }{\partial x}}+{\frac {\partial \Psi }{\partial z}}\right),$

(L-69c)

and simplify to get

$P_{zz}=(\lambda +2\mu ){\frac {\partial ^{2}\Phi }{\partial z^{2}}}+\lambda {\frac {\partial ^{2}\Phi }{\partial x^{2}}}-2\mu {\frac {\partial ^{2}\Psi }{\partial z\partial x}}.$

(L-69d)

Rewrite equation (L-69d) as

$P_{zz}=(\lambda +2\mu )\left({\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial z^{2}}}\right)-2\mu \left({\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Psi }{\partial z\partial x}}\right),$

(L-69e)

and substitute equations (L-35), (L-47), and (L-59) to get the final expression for the normal stress P_zz in terms of the displacement potentials Φ and Ψ

$P_{zz}=\rho {\frac {\partial ^{2}\Phi }{\partial t^{2}}}-2\rho \beta ^{2}\left({\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Psi }{\partial z\partial x}}\right).$

(L-69f)

Now, refer to equation (L-69f) and consider the third boundary condition expressed as

$\rho _{1}{\frac {\partial ^{2}\Phi _{1}}{\partial t^{2}}}-2\rho _{1}\beta _{1}^{2}\left({\frac {\partial ^{2}\Phi _{1}}{\partial x^{2}}}+{\frac {\partial ^{2}\Psi _{1}}{\partial z\partial x}}\right)=\rho _{2}{\frac {\partial ^{2}\Phi _{2}}{\partial t^{2}}}-2\rho _{2}\beta _{2}^{2}\left({\frac {\partial ^{2}\Phi _{2}}{\partial x^{2}}}+{\frac {\partial ^{2}\Psi _{2}}{\partial z\partial x}}\right).$

(L-70)

Differentiate the respective wave functions given by equations (L-63a,63b,63c,63d) and substitute into equation (L-70) to get the expression at the boundary z = 0

${\begin{array}{l}-\rho _{1}A_{0}-\rho _{1}A_{1}+2\rho _{1}{\frac {\beta _{1}^{2}}{\alpha _{1}^{2}}}\sin ^{2}\varphi _{0}\,A_{0}+2\rho _{1}{\frac {\beta _{1}^{2}}{\alpha _{1}^{2}}}\sin ^{2}\varphi _{1}\,A_{1}-2\rho _{1}\sin \psi _{1}\cos \psi _{1}B_{1}\\\quad \quad \quad \quad \quad =-\rho _{2}A_{2}+2\rho _{2}{\frac {\beta _{2}^{2}}{\alpha _{2}^{2}}}\sin ^{2}\varphi _{2}\,A_{2}+2\rho _{2}\sin \psi _{2}\cos \psi _{2}B_{2}.\end{array}}$

(L-71a)

Simplify by combining terms and using the trigonometric relation sin 2ψ = 2 sin ψ cos ψ to get

${\begin{array}{l}-\left(1-2{\frac {\beta _{1}^{2}}{\alpha _{1}^{2}}}\sin ^{2}\varphi _{1}\right)A_{0}-\left(1-2{\frac {\beta _{1}^{2}}{\alpha _{1}^{2}}}\sin ^{2}\varphi _{1}\right)A_{1}-\sin 2\psi _{1}B_{1}\\\quad \quad \quad \quad \quad =-{\frac {\rho _{2}}{\rho _{1}}}\left(1-2{\frac {\beta _{2}^{2}}{\alpha _{2}^{2}}}\sin ^{2}\varphi _{2}\right)A_{2}+{\frac {\rho _{2}}{\rho _{1}}}\sin 2\psi _{2}B_{2}.\end{array}}$

(L-71b)

By way of Snell’s law given by equation (L-64), note that

${\frac {\beta _{1}^{2}}{\alpha _{1}^{2}}}\sin ^{2}\varphi _{1}=\sin ^{2}\psi _{1}$

(L-72a)

and

${\frac {\beta _{2}^{2}}{\alpha _{2}^{2}}}\sin ^{2}\varphi _{2}=\sin ^{2}\psi _{2}.$

(L-72b)

Substitute equations (L-72a,L-72b) into equation (L-71b), while noting the trigonometric relation 1 − 2 sin² ψ = cos 2ψ

$-\cos 2\psi _{1}A_{0}-\cos 2\psi _{1}A_{1}-\sin 2\psi _{1}B_{1}=-{\frac {\rho _{2}}{\rho _{1}}}\cos 2\psi _{2}A_{2}+{\frac {\rho _{2}}{\rho _{1}}}\sin 2\psi _{2}B_{2}.$

(L-73)

Finally, refer to equation (L-20b) rewritten below as

$P_{xz}=2\mu e_{xz}$

and substitute equation (L-5a) for e_xz to get

$P_{xz}=\mu \left({\frac {\partial w}{\partial x}}+{\frac {\partial u}{\partial z}}\right).$

(L-74a)

Next, substitute equations (L-55a,L-55b) into equation (L-74a) for the displacement potentials

$P_{xz}=\mu \left[{\frac {\partial }{\partial x}}\left({\frac {\partial \Phi }{\partial z}}-{\frac {\partial \Psi }{\partial x}}\right)+{\frac {\partial }{\partial z}}\left({\frac {\partial \Phi }{\partial x}}+{\frac {\partial \Psi }{\partial z}}\right)\right],$

(L-74b)

and simplify to get

$P_{xz}=\mu \left(2{\frac {\partial ^{2}\Phi }{\partial z\partial x}}+{\frac {\partial ^{2}\Psi }{\partial z^{2}}}-{\frac {\partial ^{2}\Psi }{\partial x^{2}}}\right).$

(L-74c)

Now, refer to equation (L-74c) and consider the fourth boundary condition expressed as

$\mu _{1}\left(2{\frac {\partial ^{2}\Phi _{1}}{\partial z\partial x}}+{\frac {\partial ^{2}\Psi _{1}}{\partial z^{2}}}-{\frac {\partial ^{2}\Psi _{1}}{\partial x^{2}}}\right)=\mu _{2}\left(2{\frac {\partial ^{2}\Phi _{2}}{\partial z\partial x}}+{\frac {\partial ^{2}\Psi _{2}}{\partial z^{2}}}-{\frac {\partial ^{2}\Psi _{2}}{\partial x^{2}}}\right).$

(L-75)

Differentiate the respective wave functions given by equations (L-63a,63b,63c,63d) and substitute into equation (L-75) to get the expression at the boundary z = 0

${\begin{array}{l}\mu _{1}\left(-{\frac {2}{\alpha _{1}^{2}}}\sin \varphi _{0}\cos \varphi _{0}A_{0}+{\frac {2}{\alpha _{1}^{2}}}\sin \varphi _{1}\cos \varphi _{1}A_{1}-{\frac {1}{\beta _{1}^{2}}}\cos ^{2}\psi _{1}B_{1}+{\frac {1}{\beta _{1}^{2}}}\sin ^{2}\psi _{1}B_{1}\right)\\\quad \quad =\mu _{2}\left(-{\frac {2}{\alpha _{2}^{2}}}\sin \varphi _{2}\cos \varphi _{2}A_{2}-{\frac {1}{\beta _{2}^{2}}}\cos ^{2}\psi _{2}B_{2}+{\frac {1}{\beta _{2}^{2}}}\sin ^{2}\psi _{2}B_{2}\right).\end{array}}$

(L-76a)

Simplify by combining terms and using the trigonometric relation sin 2ψ = 2 sin ψ cos ψ to get

${\begin{array}{l}-{\frac {\mu _{1}}{\alpha _{1}^{2}}}\sin 2\varphi _{1}A_{0}+{\frac {\mu _{1}}{\alpha _{1}^{2}}}\sin 2\varphi _{1}A_{1}-{\frac {\mu _{1}}{\beta _{1}^{2}}}(\cos ^{2}\psi _{1}-\sin ^{2}\psi _{1})B_{1}\\\quad \quad =-{\frac {\mu ^{2}}{\alpha _{2}^{2}}}\sin 2\varphi _{2}A_{2}-{\frac {\mu _{2}}{\beta _{2}^{2}}}(\cos ^{2}\psi _{2}-\sin ^{2}\psi _{2})B_{2}.\end{array}}$

(L-76b)

Substitute equation (L-47) into equation (L-76b), while noting the trigonometric relations sin 2ψ = 2 sin ψ cos ψ and cos 2ψ = cos² ψ − sin² ψ to get

$-\sin 2\varphi _{1}A_{0}+\sin 2\varphi _{1}A_{1}-{\frac {\alpha _{1}^{2}}{\beta _{1}^{2}}}\cos 2\psi _{1}B_{1}=-{\frac {\rho _{2}\beta _{2}^{2}\alpha _{1}^{2}}{\rho _{1}\beta _{1}^{2}\alpha _{2}^{2}}}\sin 2\varphi _{2}A_{2}-{\frac {\rho _{2}\alpha _{1}^{2}}{\rho _{1}\beta _{1}^{2}}}\cos 2\psi _{2}B_{2}.$

(L-77)

Equations (L-66), (L-68), (L-73), and (L-77) are the Zoeppritz equations that can be solved for the four unknowns — the reflected compressional-wave amplitude A₁, the reflected shear-wave amplitude B₁, the transmitted compressional-wave amplitude A₂, and the transmitted shear-wave amplitude B₂.

Equations (L-66), (L-68), (L-73), and (L-77) are best displayed in matrix form normalized by the incident-wave amplitude A₀ = 1:

${\begin{pmatrix}\cos \varphi _{1}&{\frac {\alpha _{1}}{\beta _{1}}}\sin \psi _{1}&{\frac {\alpha _{1}}{\alpha _{2}}}\cos \varphi _{2}&-{\frac {\alpha _{1}}{\beta _{2}}}\sin \psi _{2}\\-\sin \varphi _{1}&{\frac {\alpha _{1}}{\beta _{1}}}\cos \psi _{1}&{\frac {\alpha _{1}}{\alpha _{2}}}\sin \varphi _{2}&{\frac {\alpha _{1}}{\beta _{2}}}\cos \psi _{2}\\-\cos 2\psi _{1}&-\sin 2\psi _{1}&{\frac {\rho _{2}}{\rho _{1}}}\cos 2\psi _{2}&-{\frac {\rho _{2}}{\rho _{1}}}\sin 2\psi _{2}\\\sin 2\varphi _{1}&-{\frac {\alpha _{1}^{2}}{\beta _{1}^{2}}}\cos 2\psi _{1}&{\frac {\rho _{2}\beta _{2}^{2}\alpha _{1}^{2}}{\rho _{1}\beta _{1}^{2}\alpha _{2}^{2}}}\sin 2\varphi _{2}&{\frac {\rho _{2}\alpha _{1}^{2}}{\rho _{1}\beta _{1}^{2}}}\cos 2\psi _{2}\end{pmatrix}}{\begin{pmatrix}A_{1}\\B_{1}\\A_{2}\\B_{2}\end{pmatrix}}={\begin{pmatrix}\cos \varphi _{1}\\\sin \varphi _{1}\\\cos 2\psi _{1}\\\sin 2\varphi _{1}\end{pmatrix}}.$

(L-78)

Given the model parameters and angles of propagation as in Figure L-4, equation (L-78) can be used to model angle-dependent reflection amplitudes. In the next section, we shall discuss the practical use of an approximate solution of the Zoeppritz equations for P-to-P reflection amplitudes to derive the AVO attributes.

L.6 Prestack amplitude inversion

Consider the two-layered earth model in Figure L-4. The Zoeppritz equations (L-78) can be solved for the amplitudes of the reflected and refracted P- and S-wave potentials of equations (L-63a,63b,63c,63d) — A₁, B₁, A₂, and B₂. However, our interest in exploration seismology is largely the angle dependency of the P-to-P reflections given by the coefficient A₁. Specifically, we wish to infer or possibly estimate elastic parameters of reservoir rocks from reflection amplitudes and relate these parameters to reservoir fluids.

The exact expression for A₁ derived from the solution of the Zoeppritz equations (L-78) is complicated and not intuitive in terms of its practical use for inferring petrophysical properties of reservoir rocks. Instead, we shall use the approximation provided by Aki and Richards as the starting point for deriving a series of practical AVO equations. Now that we only need to deal with the P-to-P reflection amplitude A₁, we shall switch to the conventional notation by replacing A₁ with R(θ) as the angle-dependent reflection amplitude for AVO analysis.

By assuming that changes in elastic properties of rocks across the layer boundary are small and propagation angles are within the subscritical range, the exact expression for R(θ) given by the Zoeppritz equation can be approximated by

$R(\theta )=\left[{\frac {1}{2}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta \alpha }{\alpha }}-\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \beta }{\beta }}+\left[{\frac {1}{2}}\left(1-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }},$

(L-79)

where α = (α₁ + α₂)/2, average P-wave velocity and Δα = (α₂ − α₁), β = (β₁ + β₂)/2, average S-wave velocity and Δβ = β₂ − β₁, ρ = (ρ₁ + ρ₂)/2, average density and Δρ = ρ₂ − ρ₁, and θ = (ψ₁ + ψ₂)/2, average of the incidence and transmission angles for the P-wave (Figure L-4).

Consider the case of N-fold CMP data represented in the domain of angle of incidence. In discrete form, equation (L-79) can be rewritten as

$R_{i}=a_{i}{\frac {\Delta \alpha }{\alpha }}+b_{i}{\frac {\Delta \beta }{\beta }}+c_{i}{\frac {\Delta \rho }{\rho }},$

(L-80)

where i is the trace index and the coefficients a_i, b_i, and c_i are given by

$a_{i}={\frac {1}{2}}\left(1+\tan ^{2}\theta _{i}\right),$

(L-81a)

$b_{i}=-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i},$

(L-81b)

and

$c_{i}={\frac {1}{2}}\left(1-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i}\right).$

(L-81c)

Note that the reflection amplitude R_i is a linear combination of three elastic parameters — fractional change in P-wave velocity Δα/α, fractional change in S-wave velocity Δβ/β, and fractional change in density Δρ/ρ.

We have, for each CMP location and for a specific reflection event associated with a layer boundary, N such equations (L-80) and three unknowns — Δα/α, Δβ/β, and Δρ/ρ. We have more equations than unknowns; hence, we encounter yet another example of generalized linear inversion (GLI) problem. The objective is to determine the three parameters such that the difference between the modeled reflection amplitudes R_i represented by equation (L-80) and the actual reflection amplitudes X_i is minimum in the least-squares sense . The least-squares error energy S is defined as

$S=\sum \limits _{i}^{N}(R_{i}-X_{i})^{2}.$

(L-82)

Smith and Gidlow argue in favor of solving for only two of the three parameters by making use of the empirical relation between density ρ and P-wave velocity α :

$\rho =k\alpha ^{1/4},$

(L-83a)

where k is a scalar. This relation holds for most water-saturated rocks. Differentiate to get

${\frac {\Delta \rho }{\rho }}={\frac {1}{4}}{\frac {\Delta \alpha }{\alpha }}.$

(L-83b)

Now substitute equation (L-83b) into the original Aki-Richards equation (L-80) and combine the terms with Δα/α

$R(\theta )=\left[{\frac {1}{2}}\left(1+\tan ^{2}\theta \right)+{\frac {1}{8}}\left(1-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right)\right]{\frac {\Delta \alpha }{\alpha }}-\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \beta }{\beta }}.$

(L-84)

Simplify the algebra to obtain the desired two-parameter model equation for prestack amplitude inversion

$R(\theta )=\left[{\frac {5}{8}}-{\frac {1}{2}}{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta +{\frac {1}{2}}\tan ^{2}\theta \right]{\frac {\Delta \alpha }{\alpha }}-\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \beta }{\beta }}.$

(L-85)

Redefine the coefficients a_i and b_i and write the discrete form of equation (L-85)

$R_{i}=a_{i}{\frac {\Delta \alpha }{\alpha }}+b_{i}{\frac {\Delta \beta }{\beta }},$

(L-86)

where

$a_{i}={\frac {5}{8}}-{\frac {1}{2}}{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i}+{\frac {1}{2}}\tan ^{2}\theta _{i}$

(L-87a)

and

$b_{i}=-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i}.$

(L-87b)

Now substitute equation (L-86) into the error-energy equation (L-82) and expand the terms

$S=\sum \limits _{i}^{N}\left[a_{i}^{2}\left({\frac {\Delta \alpha }{\alpha }}\right)^{2}+2a_{i}b_{i}{\frac {\Delta \alpha }{\alpha }}{\frac {\Delta \beta }{\beta }}+b_{i}^{2}\left({\frac {\Delta \beta }{\beta }}\right)^{2}-2\left(a_{i}{\frac {\Delta \alpha }{\alpha }}+b_{i}{\frac {\Delta \beta }{\beta }}\right)X_{i}+X_{i}^{2}\right].$

(L-88)

Minimization of E with respect to the parameters Δα/α and Δβ/β that we wish to estimate from prestack amplitude inversion requires that

${\frac {\partial S}{\partial (\Delta \alpha /\alpha )}}=0$

(L-89a)

and

${\frac {\partial S}{\partial (\Delta \beta /\beta )}}=0.$

(L-89b)

Compute the partial derivatives using equation (L-88) and substitute into equations (L-89a,L-89b) to get two simultaneous equations

$\left(\sum \limits _{i}^{N}a_{i}^{2}\right){\frac {\Delta \alpha }{\alpha }}+\left(\sum \limits _{i}^{N}a_{i}b_{i}\right){\frac {\Delta \beta }{\beta }}=\sum \limits _{i}^{N}a_{i}X_{i}$

(L-90a)

and

$\left(\sum \limits _{i}^{N}b_{i}^{2}\right){\frac {\Delta \beta }{\beta }}+\left(\sum \limits _{i}^{N}a_{i}b_{i}\right){\frac {\Delta \alpha }{\alpha }}=\sum \limits _{i}^{N}b_{i}X_{i}.$

(L-90b)

These equations are put into matrix form

${\begin{pmatrix}\sum \nolimits _{i}^{N}\,a_{i}^{2}&\sum \nolimits _{i}^{N}\,a_{i}b_{i}\\\sum \nolimits _{i}^{N}\,a_{i}b_{i}&\sum \nolimits _{i}^{N}\,b_{i}^{2}\end{pmatrix}}{\begin{pmatrix}{\frac {\Delta \alpha }{\alpha }}\\{\frac {\Delta \beta }{\beta }}\end{pmatrix}}={\begin{pmatrix}\sum \nolimits _{i}^{N}\,a_{i}X_{i}\\\sum \nolimits _{i}^{N}\,b_{i}X_{i}\end{pmatrix}},$

(L-91)

and solved for Δα/α and Δβ/β.

The two parameters Δα/α and Δβ/β represent fractional changes in P- and S-wave velocities; as such, they are related to P- and S-wave reflectivities, ΔI_P/I_P and ΔI_S/I_S, respectively, where I_P and I_S are the P- and S-wave impedances given by

$I_{P}=\alpha \rho$

(L-92a)

and

$I_{S}=\beta \rho .$

(L-92b)

The fractional change in P-wave impedance is given by

$\Delta I_{P}={\frac {\partial I_{P}}{\partial \alpha }}\Delta \alpha +{\frac {\partial I_{P}}{\partial \rho }}\Delta \rho .$

(L-93a)

Apply differentiation to equation (L-92a) and substitute into equation (L-93a), then divide both sides of the resulting equation by I_P to get the expression for P-wave reflectivity

${\frac {\Delta I_{P}}{I_{P}}}={\frac {\Delta \alpha }{\alpha }}+{\frac {\Delta \rho }{\rho }}.$

(L-93b)

Similarly, the fractional change in S-wave impedance is given by

$\Delta I_{S}={\frac {\partial I_{S}}{\partial \beta }}\Delta \beta +{\frac {\partial I_{S}}{\partial \rho }}\Delta \rho .$

(L-94a)

Apply differentiation to equation (L-92b) and substitute into equation (L-94a), then divide both sides of the resulting equation by I_S to get the expression for S-wave reflectivity

${\frac {\Delta I_{S}}{I_{S}}}={\frac {\Delta \beta }{\beta }}+{\frac {\Delta \rho }{\rho }}.$

(L-94b)

Petrophysical interpretation is then made using the crossplot of P- and S-wave reflectivities. Note from the definitions of the coefficients a_i and b_i given by equations (L-87a,L-87b) that you have to choose a value for the ratio β/α to compute the two reflectivities.

An alternative formulation of prestack amplitude inversion is based on transforming the Aki-Richards equation (L-79) to new variables in terms of Lamé’s constants λ and μ . Refer to equation (L-35) and define a new quantity ν = λ + 2μ. The relationship between the two sets of variables (ν, u) and (α, β, ρ), by way of equations (L-35) and (L-47), is given by

$\nu =\rho \alpha ^{2}$

(L-95a)

and

$\mu =\rho \beta ^{2}.$

(L-95b)

The fractional changes in ν and μ are given by

$\Delta \nu ={\frac {\partial \nu }{\partial \alpha }}\Delta \alpha +{\frac {\partial \nu }{\partial \rho }}\Delta \rho$

(L-96a)

and

$\Delta \mu ={\frac {\partial \mu }{\partial \beta }}\Delta \beta +{\frac {\partial \mu }{\partial \rho }}\Delta \rho .$

(L-96b)

Apply differentiation to equation (L-95a) and substitute into equation (L-96a), then divide both sides by ν

${\frac {\Delta \nu }{\nu }}=2{\frac {\Delta \alpha }{\alpha }}+{\frac {\Delta \rho }{\rho }}.$

(L-97a)

Similarly, apply differentiation to equation (L-95b) and substitute into equation (L-96b), then divide both sides by μ

${\frac {\Delta \mu }{\mu }}=2{\frac {\Delta \beta }{\beta }}+{\frac {\Delta \rho }{\rho }}.$

(L-97b)

Now solve equations (L-97a,L-97b) for Δα/α and Δβ/β, respectively, and substitute into the Aki-Richards equation (L-79)

${\begin{array}{l}R(\theta )=\left[{\frac {1}{4}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta \nu }{\nu }}-\left[{\frac {1}{4}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }}\\\quad \quad \ -\left[2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \mu }{\mu }}+\left[2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \rho }{\rho }}+\left[{\frac {1}{2}}\left(1-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }}.\end{array}}$

(L-98)

Finally, simplify to get the Aki-Richards equation transformed to the new variables Δν/ν and Δμ/μ

$R(\theta )=\left[{\frac {1}{4}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta \nu }{\nu }}-\left[2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \mu }{\mu }}+\left[{\frac {1}{4}}\left(1-\tan ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }}.$

(L-99)

As for the original Aki-Richards equation (L-79), consider the case of N-fold CMP data represented in the domain of angle of incidence. In discrete form, equation (L-99) can be rewritten as

$R_{i}=a_{i}{\frac {\Delta \nu }{\nu }}+b_{i}{\frac {\Delta \mu }{\mu }}+c_{i}{\frac {\Delta \rho }{\rho }},$

(L-100)

where i is the trace index and the coefficients a_i, b_i, and c_i are given by

$a_{i}={\frac {1}{4}}\left(1+\tan ^{2}\theta _{i}\right),$

(L-101a)

$b_{i}=-2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i},$

(L-101b)

and

$c_{i}={\frac {1}{4}}\left(1-\tan ^{2}\theta _{i}\right).$

(L-101c)

The reflection amplitude R_i in equation (L-100) is a linear combination of three elastic parameters — fractional changes Δν/ν, Δμ/μ, and Δρ/ρ. Based on the least-squares minimization given by equation (L-82), these three parameters can, in principle, be estimated for a specific event from CMP data.

A variant to the transformation given by equation (L-99) can be formulated in terms of the P- and S-wave reflectivities, ΔI_P/I_P and ΔI_S/I_S, respectively. Solve equation (L-93b) for Δα/α and equation (L-94b) for Δβ/β, and substitute into the Aki-Richards equation (L-79)

${\begin{array}{l}R(\theta )=\left[{\frac {1}{2}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta I_{P}}{I_{P}}}-\left[{\frac {1}{2}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }}\\\quad \quad \,\,-\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta I_{S}}{I_{S}}}+\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \rho }{\rho }}+\left[{\frac {1}{2}}\left(1-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right)\right]{\frac {\Delta \rho }{\rho }}.\end{array}}$

(L-102)

Simplify to get the Aki-Richards equation transformed to the new variables ΔI_P/I_P and ΔI_S/I_S

$R(\theta )=\left[{\frac {1}{2}}\left(1+\tan ^{2}\theta \right)\right]{\frac {\Delta I_{P}}{I_{P}}}-\left[4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta I_{S}}{I_{S}}}-\left[{\frac {1}{2}}\tan ^{2}\theta -2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta \right]{\frac {\Delta \rho }{\rho }}.$

(L-103)

Again, consider the case of N-fold CMP data represented in the domain of angle of incidence. In discrete form, equation (L-103) can be rewritten as

$R_{i}=a_{i}{\frac {\Delta I_{P}}{I_{P}}}+b_{i}{\frac {\Delta I_{S}}{I_{S}}}+c_{i}{\frac {\Delta \rho }{\rho }},$

(L-104)

where i is the trace index and the coefficients a_i, b_i, and c_i are given by

$a_{i}={\frac {1}{2}}\left(1+\tan ^{2}\theta _{i}\right),$

(L-105a)

$b_{i}=-4{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i},$

(L-105b)

and

$c_{i}=-\left({\frac {1}{2}}\tan ^{2}\theta -2{\frac {\beta ^{2}}{\alpha ^{2}}}\sin ^{2}\theta _{i}\right).$

(L-105c)

The reflection amplitude R_i in equation (L-104) is a linear combination of three parameters — ΔI_P/I_P, ΔI_S/I_S, and Δρ/ρ. Again, based on the least-squares minimization given by equation (L-82), these three parameters can be estimated for a specific event from CMP data.

Goodway implemented a specific form of equation (L-104) to derive the AVO attributes ΔI_P/I_P and ΔI_S/I_S. For a specific value of α/β = 2 and small angles of incidence, the third term in equation (L-104) vanishes. The resulting equation then is solved for the P- and S-wave reflectivities, ΔI_P/I_P and ΔI_S/I_S, respectively, using the least-squares solution as in equation (L-91)

${\begin{pmatrix}\sum \nolimits _{i}^{N}\,a_{i}^{2}&\sum \nolimits _{i}^{N}\,a_{i}b_{i}\\\sum \nolimits _{i}^{N}\,a_{i}b_{i}&\sum \nolimits _{i}^{N}\,b_{i}^{2}\end{pmatrix}}{\begin{pmatrix}{\frac {\Delta I_{P}}{I_{P}}}\\{\frac {\Delta I_{S}}{I_{S}}}\end{pmatrix}}={\begin{pmatrix}\sum \nolimits _{i}^{N}\,a_{i}X_{i}\\\sum \nolimits _{i}^{N}\,b_{i}X_{i}\end{pmatrix}}.$

(L-106)

See also

• Introduction to reservoir geophysics
• Seismic resolution
• Analysis of amplitude variation with offset
• Acoustic impedance estimation
• Vertical seismic profiling
• 4-D seismic method
• 4-C seismic method
• Seismic anisotropy
• Exercises

External links

find literature about Mathematical foundation of elastic wave propagation