Elastic notes, 1993: W.S. Harlan

Raypath extrapolation of vector wavefields in arbitrary elastic media

William S. Harlan 1993

Preliminaries

Begin with conservation of momentum for a continuous medium:

ρ(x̃) · ü_{j}(x̃, t) − τ_{jk,k}(x̃, t) = F_{j}(x̃, t),  (1)

where t is time, and x̃ (with elements x_{j}) the vector of spatial coordinates. ρ is the density, u_{j} the displacement vector, F_{j} the body force vector, and τ_{jk} the stress tensor. Each dot indicates a time derivative.

Tensor notation will be used. Ordinary subscripts index the spatial components of vectors and tensors. Subscripts appearing after a comma indicate a spatial derivative in the indexed direction. The summation convention requires that repeated subscripts, like k in equation (1), be implicitly summed and eliminated.

Next, Hooke’s Law assumes a linear relationship between the stress tensor and a spatially differentiated displacement vector, called the strain tensor:

τ_{jk}(x̃, t) = c_{jklm}(x̃) · u_{l,m}(x̃, t).  (2)

The tensor c_{jklm} is the elastic stiffness, with a maximum of 21 independent components in a three-dimensional coordinate system. The following symmetries always hold:

\begin{aligned}c_{jklm} & = c_{kjlm} = c_{lmjk}; so \\ c_{jklm} & = c_{kjlm} = c_{jkml} = c_{kjml} = c_{lmjk} = c_{mljk} = c_{lmkj} = c_{mlkj}. \end{aligned}  (3)

Combining the equations (1) and (2) to eliminate the stress tensor produces the elastic wave equation:

ρ(x̃) · ü_{j}(x̃, t) − [c_{jklm}(x̃)u_{l,m}(x̃, t)]_{,k} = F_{j}(x̃, t).  (4)

We use the Fourier kernal exp(i2πft) to transform from frequency f to time t. The inverse transformation uses the complex conjugate.

4π^{2}f^{2}ρ(x̃)ũ_{j}(x̃, f) + [c_{jklm}(x̃)ũ_{l,m}(x̃, f)]_{,k} + F̃_{j}(x̃, f) = 0.  (5)

Tildes will mark Fourier transformed functions.

Energy Flow

A raypath should indicate the flow and direction of energy in a wavefield. Only in isotropic media can we assume that this flow of energy is perpendicular to wavefronts.

Elastic notes, 1993: W.S. Harlan Page 2

To describe energy flow, we do not even need to know how to decompose the wavefield into modes and wavefronts.

Take the dot product of each side of the elastic wave equation (4) with u̇_{j}.

\begin{aligned}F_{j}u̇_{j} & = ρü_{j}u̇_{j} − (c_{jklm}u_{l,m})_{,k}u̇_{j} \\ & = \frac{d}{dt}(\frac{1}{2}ρu̇_{j}u̇_{j}) − (c_{jklm}u_{l,m}u̇_{j})_{,k} + c_{jklm}u_{l,m}u̇_{j,k} \\ & = \frac{d}{dt}(\frac{1}{2}ρu̇_{j}u̇_{j} + \frac{1}{2}c_{jklm}u_{j,k}u_{l,m}) − (c_{jklm}u_{l,m}u̇_{k})_{,j} . \end{aligned}  (6)

(The symmetries of equation (3) were used.) This equation is in the form of a conservation law:

\frac{dσ(x̃, t)}{dt} + ∇ · J̃(x̃, t) = \frac{dσ(x̃, t)}{dt} + J_{j,j}(x̃, t) = 0,  (7)

where

σ(x̃, t) ≡ \frac{1}{2}ρ(x̃)| ˙ũ(x̃, t)|^{2} + \frac{1}{2}c_{jklm}(x̃)u_{j,k}(x̃, t)u_{l,m}(x̃, t)  (8)

is the total energy density, and

J_{j}(x̃, t) ≡ −c_{jklm}(x̃)u_{l,m}(x̃, t)u̇_{j}(x̃, t) = −τ_{jk}(x̃, t)u̇_{k}(x̃, t)  (9)

is the energy flow, or Poynting vector.

Define the group velocity (or energy velocity) g_{j} by

J_{j}(x̃, t) ≡ σ(x̃, t)g_{j}(x̃, t).  (10)

The energy flow equals the energy density times the energy (group) velocity. Raypaths are generally defined to follow energy. Let us define a raypath as a curve which is tangent to the group velocity vector at every point. In the anisotropic case, raypaths will not necessarily be perpendicular to wavefronts.

Wave Mode Decomposition

Assume that a solution for the displacement vector u_{j}(x̃, t) is well approximated locally by the approximation:

u_{j}(x̃, t) ≈ a_{j}(x̃)ϕ[t − T (x̃)]; and ũ_{j}(x̃, f) ≈ ã_{j}(x̃)ϕ̃(f)exp[−i2πfT (x̃)].  (11)

where a_{j} gives the local amplitude and ϕ modulates the phase according to the time delay T.

Substituting equation (11) into (4) gives the following (suppressing arguments):

\begin{aligned}0 = & − ρa_{j}ϕ̈ + [c_{jklm}(a_{l,m}ϕ + a_{l}ϕ̇s_{m})]_{,k} + F_{j} \\ = & − ρa_{j}ϕ̈ + c_{jklm,k}(a_{l,m}ϕ + a_{l}ϕ̇s_{m}) \\ & + c_{jklm}(a_{l,km}ϕ + a_{l,m}ϕ̇s_{k} + a_{l,k}ϕ̇s_{m} + a_{l}ϕ̈s_{k}s_{m} + a_{l}ϕ̇s_{m,k}) + F_{j} \\ = & − (ρa_{j} − c_{jklm}a_{l}s_{k}s_{m})ϕ̈ \\ & + [c_{jklm,k}a_{l}s_{m} + c_{jklm}(a_{l,m}s_{k} + a_{l,k}s_{m} + a_{l}s_{m,k})]ϕ̇ \\ & + (c_{jklm,k}a_{l,m} + c_{jklm}a_{l,km})ϕ + F_{j}, \end{aligned}  (12)

where s_{j}(x̃) = T_{,j}(x̃) is the phase slowness vector, pointing in the direction of a local plane wave, with a magnitude equal to the reciprocal of the wave velocity.

In the frequency domain, we find a similar expression, grouped by powers of f.

\begin{matrix}0 = \{ & 4π^{2}f^{2}(ρa_{j} − c_{jklm}a_{l}s_{k}s_{m}) \\ + & i2πf[c_{jklm,k}a_{l}s_{m} + c_{jklm}(a_{l,m}s_{k} + a_{l,k}s_{m} + a_{l}s_{m,k})] \\ + & (c_{jklm,k}a_{l,m} + c_{jklm}a_{l,km})\}ϕ̃ + F̃_{j}. \end{matrix}  (13)

To allow non-vanishing ϕ, each of these three scaled terms should vanish in areas with a vanishing body force F_{j}.

The first term produces the equation appropriate for a high-frequency limit—a vector version of the Eikonal equation:

a_{j}(x̃) = [c_{jklm}(x̃)s_{k}(x̃)s_{m}(x̃)/ρ(x̃)] · a_{l}(x̃), ∀j.  (14)

Notice that this equation is independent of time or frequency. If we scale equation (14) with a phase velocity v(x̃) = |s̃(x̃)|^{−1}, the reciprocal of the magnitude of the phase slowness, then

v(x̃)^{2}a_{j}(x̃) = [c_{jklm}(x̃)ŝ_{k}(x̃)ŝ_{m}(x̃)/ρ(x̃)] · a_{l}(x̃), ∀j.  (15)

Choose a unit normal ŝ_{j} = vs_{j} to a particular plane wave. The eigenvalues of equation (15) equal the squared velocities, and the eigenvectors give the polarity of different wave modes (i.e. compressional and shear waves).

Equations (14) and (15) resemble the Christoffel equation, which Fourier transforms the spatial dimensions and assumes a homogeneous material. Rather than assume global homogeneity, I prefer a high-frequency approximation that assumes stiffness to change slowly over the spatial wavelengths of the propagating wavefront.

A plot of |s̃| as a function of ˆ s̃ is called a slowness surface. A plot of v = |s̃|^{−1} as a function of ŝ is called a normal surface.

The remaining terms of (13) give “Transport equations” in regions of vanishing body sources:

[c_{jklm}(x̃)s_{m}(x̃)a_{l}(x̃)]_{,k} + c_{jklm}(x̃)s_{k}(x̃)a_{l,m}(x̃) = 0, and
[c_{jklm}(x̃)a_{l,m}(x̃)]_{,k} = 0, ∀j. (16)

These functions are also independent of time and frequency.

Mode energy flow

Substituting the approximation (11) into the total energy density (8), we find

σ(x̃, t) = \frac{1}{2}ρ(x̃)|ã(x̃)|^{2}ϕ̇[t −T (x̃)]^{2} + \frac{1}{2}c_{jklm}(x̃)a_{j}(x̃)s_{k}(x̃)a_{l}(x̃)s_{m}(x̃)ϕ̇[t −T (x̃)]^{2} . (17)

Because of equation (14), these two terms are equal—i.e., the kinetic equals the potential energy density.

σ(x̃, t) = c_{jklm}(x̃)a_{j}(x̃)a_{l}(x̃)s_{k}(x̃)s_{m}(x̃)ϕ̇[t − T (x̃)]^{2} .  (18)

Similarly, the energy flow (9) becomes

J_{i}(x̃, t) = c_{jklm}(x̃)a_{k}(x̃)a_{m}(x̃)s_{l}(x̃)ϕ̇[t − T (x̃)]^{2} .  (19)

Substituting energy density (18) and energy flow (19) into the definition of group velocity (10), we find the following equality:

c_{jklm}(x̃)a_{k}(x̃)a_{m}(x̃)s_{l}(x̃) ≡ [c_{nklm}(x̃)a_{n}(x̃)a_{l}(x̃)s_{k}(x̃)s_{m}(x̃)] · g_{j}(x̃); and

g_{j}(x̃) = \frac{c_{jklm}(x̃)a_{k}(x̃)a_{m}(x̃)s_{l}(x̃)}{c_{nklm}(x̃)a_{n}(x̃)a_{l}(x̃)s_{k}(x̃)s_{m}(x̃)} = [c_{jklm}(x̃)â_{k}(x̃)â_{m}(x̃)/ρ(x̃)]s_{l}(x̃). (20)

Equation (14) has allowed the last equality. Notice that, under this approximation, the group velocity is function only of position, not time or frequency. A plot of |g̃| as a function of direction ˆ g̃ is called a ray surface.

If both sides of equation (20) are dotted with the slowness vector s_{j}(x̃), then we find that

g_{j}(x̃)s_{j}(x̃) = 1.  (21)

This important result states that the dot product of the group velocity vector with the slowness vector is unity. The slowness vector is perpendicular to the wavefront, by construction, but the group velocity vector is not, unless the group and phase velocities are equal.

To propagate a wavefront one small step Δt in time, we first calculate the normal ˆ s̃ to each point on the wavefront. The phase velocity v = |ˆs̃|^{−1} along the wavefront is given by the appropriate eigenvalues of equation (14). With this value of s̃ = ˆ s̃/v, we can calculate the group velocity g̃ by substitution into equation (20). Each point on the wavefront can be extrapolated in the direction of g̃ by a perturbation g̃ · Δt. Finally, we connect the revised points on the wavefront and recalculate the normal vectors ˆ s̃.

It is possible to calculate a group velocity surface from a slowness surface numerically, and vice versa. Assuming the elastic material desribed by the stiffness tensor ˜ C̃ and the density ρ to be fixed, let us perturb the slowness vector s̃+δs̃ and particle motion ˆã + δˆã in a way that continues to satisfy equation (14). Expanded to first order,

δa_{j} = c_{jklm}/ρδs_{k}s_{m}a_{l} + c_{jklm}/ρs_{k}δs_{m}a_{l} + c_{jklm}/ρs_{k}s_{m}δa_{l}  (22)

Take the dot product of both sides with a_{j} and regroup,

a_{j}δa_{j} = 2[c_{jklm}/ρs_{m}a_{j}a_{l}]δs_{k} + [c_{jklm}/ρs_{k}s_{m}a_{j}]δa_{l} = 2g_{k}δs_{k} + a_{l}δa_{l}.  (23)

We have used the definition of group velocity (20), the Eikonal equation (14), and the symmetries of stiffness (3) to simplify the terms in brackets. Subtracting the identical terms we find the simple result

g_{j}δs_{j} = 0.  (24)

And because of (21), we also have

s_{j}δg_{j} = −g_{j}δs_{j} = 0.  (25)

The perturbation s̃+δs̃ must lie along the slowness surface to be a valid perturbation. Thus, s̃ is tangent to the slowness surface, and g̃ is perpendicular to the slowness surface. Similarly, the phase slowness is perpendicular to the ray surface.

The normal relations implied by (25) allow us to calculate the angle between s̃ and g̃. If we know the magnitude of one vector, then equation (21) allows us to calculate the magnitude of the other vector. If we have already constructed a slowness surface, then group velocity directions are calculated as normals to this surface. The magnitudes of the group velocities derive from the cosine relation (21), providing all information necessary to draw the ray surface. Similarly, the slowness surface can be constructed from a ray surface.

Impulsive source

Assume a body force with an arbitrary wavelet at an impulsive location: x̃^{0}:

F_{j}(x̃, t) = w(t)b_{j}δ(x̃ − x̃^{0}), and F̃_{j}(x̃, f) = w̃(f)δ(x̃ − x̃^{0}).  (26)

The only term of equation (13) that could cancel this impulsive term would involve the second spatial derivative of a_{j}(x̃). Thus,

c_{jklm}(x̃)a_{l,km}(x̃)ϕ(t) = −F_{j}(x̃, t) = −w(t)b_{j}δ(x̃ − x̃^{0}).  (27)

Since the approximation (11) allows for an arbitrary scaling of the amplitude and phase terms, we can assume that our solution sets

\begin{aligned}ϕ(t) & = w(t), and \\ c_{jklm}(x̃)a_{l,km}(x̃) & = b_{j}δ(x̃ − x̃^{0}). \end{aligned}  (28)