Simulating propagation of separated wave modes in general anisotropic media, Part I: qP-wave propagators

Pseudo-pure-mode qP-wave equation

To describe propagation of separated qP-waves in anisotropic media, we first revisit the classical wave mode separation theory. In isotropic media, scalar P-wave can be separated from the extrapolated vector wavefield $\mathbf{u}$ by applying a divergence operation: $P = \bigtriangledown\cdot{\mathbf{u}}$ . In the wavenumber domain, this can be equivalently expressed as a dot product that essentially projects the wavefield $\widetilde{\mathbf{u}}$ onto the wave vector $\mathbf {k}$ , i.e.,

$$\widetilde{P} = i\mathbf{k}\cdot{\widetilde{\mathbf{u}}},$$ (8)

Similarly, for an anisotropic medium, scalar qP-waves can be separated by projecting the vector wavefields onto the true polarization directions of qP-waves by (Dellinger and Etgen, 1990),

$$q\widetilde{P} = i\mathbf{a_{p}}\cdot{\widetilde{\mathbf{u}}},$$ (9)

where $\mathbf{a}_{p}=(a_{px},a_{py},a_{pz})^{T}$ represents the polarization vector for qP-waves. For heterogeneous models, this scalar projection can be performed using nonstationary spatial filtering depending on local material parameters (Yan and Sava, 2009).

To provide more flexibility for characterizing wave propagation in anisotropic media, we suggest to split the one-step projection into two steps, of which the first step implicitly implements partail wave mode separation (like in equation 8) during wavefield extrapolation with a transformed wave equation, while the second step is designed to correct the projection deviation implied by equations 8 and 9. We achieve this on the base of the following observations: the difference of the polarization between an ordinary anisotropic medium and its isotropic reference at a given wave vector direction is usually small, though exceptions are possible (Tsvankin and Chesnokov, 1990; Thomsen, 1986); The wave vector can be taken as the isotropic reference of the polarization vector for qP-waves; It is a material-independent operation to project the elastic wavefield onto the wave vector.

Therefore, we introduce a similarity transformation to the Christoffel matrix, i.e.,

$$\widetilde{\overline{\mathbf{\Gamma}}} = \mathbf{M_p}\widetilde{\mathbf{\Gamma}}\mathbf{M_p}^{-1},$$ (10)

with a invertible $3\times3$ matrix $\mathbf{M}_{p}$ related to the wave vector:

$$\mathbf{M_p}= \begin{pmatrix}i{k_x} & 0 &0 \cr 0 & i{k_y} &0 \cr 0 & 0 & i{k_z}\end{pmatrix}.$$ (11)

Accordingly, we derive an equivalent Christoffel equation,

$$\widetilde{\overline{\mathbf{\Gamma}}}\widetilde{\overline{\mathbf{u}}} = \rho{\omega}^2\widetilde{\overline{\mathbf{u}}},$$ (12)

for a transformed wavefield:

$$\widetilde{\overline{\mathbf{u}}} = \mathbf{M_p}\widetilde{\mathbf{u}}.$$ (13)

The above similarity transformation does not change the eigeinvalues of the Christoffel matrix and thus introduces no kinematic errors for the wavefields. By the way, we can obtain the same transformed Christoffel equation if matrix $\mathbf{M}_{p}$ is constructed using the normalized wavenumbers to ensure all spatial frequencies are uniformly scaled. For a locally smooth medium, applying an inverse Fourier transform to equation 12, we obtain a coupled linear second-order system kinematically equivalent to the original elastic wave equation:

$$\rho\frac{\partial^2\overline{\mathbf{u}}}{\partial t^2} = \overline{\mathbf{\Gamma}}\overline{\mathbf{u}},$$ (14)

where $\overline{\mathbf{u}}$ represents the time-space domain wavefields, and $\overline{\mathbf{\Gamma}}$ represents the Christoffel differential-operator matrix after the similarity transformation.

For the transformed elastic wavefield in the wavenumber-domain, we have

$$\widetilde{\overline{u}} = \widetilde{\overline{u}}_x + \widetild... ...{u}}_y + \widetilde{\overline{u}}_z = i\mathbf{k}\cdot{\widetilde{\mathbf{u}}}.$$ (15)

And in space-domain, we also have

$$\overline{u} = \overline{u}_x + \overline{u}_y + \overline{u}_z = \bigtriangledown\cdot{\mathbf{u}},$$ (16)

with

$$\overline{u}_x = \frac{\partial u_x}{\partial x},\qquad \overline{u}_y = \frac{\partial u_y}{\partial y}, \qquad \overline{u}_z=\frac{\partial u_z}{\partial z}.$$ (17)

These imply that the new wavefield components essentially represent the spatial derivatives of the original components of the displacement wavefield, and the transformation (equation 13) plus the summation of the transformed wavefield components (like in equation 15 or 16) essentially finishes a scalar projection of the displacement wavefield onto the wave vector. For isotropic media, such a projection directly produces scalar P-wave data. In an anisotropic medium, however, only a partial wave-mode separation is achieved becuase there is usually a direction deviation between the wave vector and the polarization vector of qP-wave. Generally, this deviation turns out to be small and its maximum value rarely exceeds $20^\circ$ for typical anisotropic earth media(Psencik and Gajewski, 1998). Because of the orthognality of qP- and qS-wave polarizations, the projection deviations of qP-waves are generally far less than those of the qSV-waves when the elastic wavefields are projected onto the isotropic references of the qP-wave's polarization vectors. As demonstrated in the synthetic examples of various symmetry and strength of anisotropy, the scalar wavefield $\overline{u}$ represents dominantly the energy of qP-waves but contains some weak residual qS-waves. This is why we call the coupled system (equation 14) a pseudo-pure-mode wave equation for qP-wave in anisotropic media.

Substituting the corresponding stiffness matrix into the above derivations, we get the extended expression of pseudo-pure-mode qP-wave equation for any anisotropic media. As demonstrated in Appendix A, pseudo-pure-mode qP-wave equation in vertical TI and orthorhombic media can be expressed as

$$\begin{split}\rho\frac{\partial^2\overline{u}_x}{\partial t^2... ...4})\frac{\partial^2{\overline{u}_y}}{\partial z^2}. \end{split}$$ (18)

Note that, unlike the original elastic wave equation, pseudo-pure-mode wave equation dose not contain mixed partial derivatives. This is a good news because it takes more computational cost to compute the mixed partial derivatives using a finite-difference algorithm with required accuracy. In the forthcoming text, we focus on demonstration of pseudo-pure-mode qP-wave equations for TI media while briefly supplement similar derivation for orthorhombic media in Appendix B.

Simulating propagation of separated wave modes in general anisotropic media, Part I: qP-wave propagators