Simulating propagation of decoupled elastic waves using low-rank approximate mixed-domain integral operators for anisotropic media

Pseudo-spectral solution of the elastic wave equation

Wave propagation in general anisotropic elastic media is governed by the linearized momentum balance law and a linear constitutive relation between the stress and strain tensors. These governing equations can be combined to write the displacement formalism as,

$$\rho\partial^2_{tt}\mathbf{u} = \bigtriangledown{\cdot{[\mathbf{C}\cdot(\bigtriangledown^{T}\cdot\mathbf{u})]}} + \mathbf{f},$$ (1)

with $\mathbf{u}=(u_x,u_y,u_z)^{T}$ represents the vector wavefields, $\mathbf{f}$ is the body force vector per unit volume, $\mathbf{C}$ is the $6\times6$ elasticity matrix representing the stiffness tensor with the Voigt's menu, and the spatial differential operator $\bigtriangledown$ has the following matrix representation:

$$\bigtriangledown = \begin{pmatrix}{\partial}_x &0 &0 &0& {\parti... ...tial}_x \cr 0& 0& {\partial}_z & {\partial}_y & {\partial}_x &0 \end{pmatrix}.$$ (2)

The pseudo-spectral method calculates the spatial derivatives using the fast Fourier transform (FFT), while approximating the temporal derivative with a finite-difference. Neglecting the source term, equation 1 is rewritten in the spatial Fourier-domain for a homogeneous medium as,

$$\partial^2_{tt}\hat{\mathbf{u}} + \mathbf{\Gamma}\hat{\mathbf{u}} = \mathbf{0},$$ (3)

in which $\hat{\mathbf{u}}$ is the wavefields in the wavenumber-domain, $\mathbf{k}=(k_x,k_y,k_z)$ is the wavenumber vector, and $\mathbf{\Gamma}=1/\rho\mathbf{L}\cdot\bar{\mathbf{C}}\cdot\mathbf{L}^{T}$ represents the $3\times3$ density normalized Christoffel matrix with the wavenumber-domain counterpart of the space differential operator (removing the imaginary unit $i$ ) satisfies,

$$\mathbf{L}= \begin{pmatrix}k_x & 0 &0 &0 & k_z & k_y \cr 0 & k_y & 0 & k_z &0 & k_x \cr 0 & 0 & k_z & k_y & k_x &0\end{pmatrix}.$$ (4)

To calculate the 2nd-order temporal derivatives, we use the standard leapfrog scheme, i.e.,

$$\partial^2_{tt}{u^{(n)}_{i}} = \frac{u^{(n+1)}_i - 2u^{(n)}_i +u^{(n-1)}_i}{\Delta{t^2}},$$ (5)

in which $\Delta{t}=t^{n+1}-t^n$ is the time-step. For constant density homogeneous media, applying the two-step time-marching scheme leads to the pseudo-spectral formula:

$$\partial^2_{tt}{\mathbf{u}^{(n)}} = \mathbf{\Psi} \mathbf{u}^{(n)},$$ (6)

with the spectral operator defined with the following kernel:

$$\mathbf{\Psi}:= (2\pi)^{-3}\int{\int{\mathbf{\Gamma}(\mathbf{k})e^{i\mathbf{k}\cdot(\mathbf{x}-\mathbf{y})}d\mathbf{y}d\mathbf{k}}}.$$ (7)

Phase terms in the integral operator can be absorbed into forward and inverse Fourier transforms. This implies that the wavefields are first transformed into wavenumber-domain using forward FFTs, then multiplied with the corresponding components of the Christoffel matrix, and finally transformed back into space-domain using inverse FFTs. For locally smooth media, we use a spatially varying Christoffel matrix to tackle the heterogeneity, i.e.,

$$\mathbf{\Psi}:= (2\pi)^{-3}\int{\int{\mathbf{\Gamma}(\mathbf{x},\mathbf{k})e^{i\mathbf{k}\cdot(\mathbf{x}-\mathbf{y})}d\mathbf{y}d\mathbf{k}}}.$$ (8)

The extended formation of this pseudo-spectral elastic wave propagator is shown in Appendix B. Spectral methods are charaterized by the use of Fourier basis functions to describe the field variables and have the advantages over finite-difference schemes that the mesh requirements are more relaxed (Kosloff et al., 1989; Liu and Li, 2000).

Simulating propagation of decoupled elastic waves using low-rank approximate mixed-domain integral operators for anisotropic media