Recursive integral time extrapolation of elastic waves using low-rank symbol approximation

Heterogeneous media

In heterogeneous media, the elastic wave equation can be expressed using the Einstein notation:

$$\rho \ddot{u}_i = (c_{ijkl} u_{k,l})_{,j} = c_{ijkl} u_{k,lj} + c_{ijkl,j} u_{k,l} \;,$$ (23)

where dot on top denotes time derivative, comma in subscript denotes spatial derivative and repeated indices imply summation. Analogously, the Christoffel matrix $\Gamma$ in the case of orthorhombic symmetry can be expanded using the chain rule

$$\Gamma = \mathbf{D}\mathbf{C}(\mathbf{x}) \mathbf{D}^\intercal$$ (24)

$$= \begin{bmatrix} \partial_x & 0 & 0 & 0 & \partial_z & \parti... ...\\ \partial_z & 0 & \partial_x\\ \partial_y & \partial_x & 0 \end{bmatrix}$$

$$= \begin{bmatrix} \partial_x & 0 & 0 & 0 & \partial_z & \partia... ...& C_{55}\partial_x \\ C_{66}\partial_y & C_{66}\partial_x & 0 \end{bmatrix}$$

$$= \begin{bmatrix} C_{11}\partial_x^2+C_{66}\partial_y^2+C_{55}\... ...l_z & C_{55}\partial_x^2+C_{44}\partial_y^2+C_{33}\partial_z^2 \end{bmatrix} +$$

$$\ \hphantom{=} \begin{bmatrix} \partial_xC_{11}\partial_x+\part... ...tial_x+\partial_yC_{44}\partial_y+\partial_zC_{33}\partial_z \end{bmatrix} \;.$$

After spatial Fourier transform, $-\Gamma \xrightarrow{\mathscr{F}} \tilde{\Gamma }$ . The Christoffel matrix $\tilde{\Gamma }$ in the case of orthorhombic anisotropy takes the form:

$$\tilde{\Gamma } = \begin{bmatrix} C_{11}k_x^2+C_{66}k_y^2+C_{55}k_z^2 & (C_{12}... ... & (C_{23}+C_{44})k_yk_z & C_{55}k_x^2+C_{44}k_y^2+C_{33}k_z^2 \end{bmatrix} -$$ (25)

$$\ \hphantom{=} i \begin{bmatrix} \partial_xC_{11}k_x+\partial_y... ... \partial_xC_{55}k_x+\partial_yC_{44}k_y+\partial_zC_{33}k_z \end{bmatrix} \;.$$

When the model is smoothly varying, which is the underlying assumption of wave mode separation (Cheng et al., 2014), the gradients of stiffnesses are insignificant. In such cases, the imaginary matrix can be dropped, which leads to the conventional real-valued Christoffel matrix in Equation 12. However, in the case of strong heterogeneity, such as at medium interfaces, the gradients of stiffnesses become significant and the imaginary part needs to be taken into account. Because $\tilde{\Gamma }$ becomes complex and non-Hermitian, the eigendecomposition of the generic matrix $\mathbf{A}$ is expressed as

$$\mathbf{A}= \frac{\tilde{\Gamma }}{\rho} = \mathbf{Q} \mathbf{V} \mathbf{Q}^{-1}$$ (26)

$$= \begin{bmatrix}\mathbf{a}_1 & \mathbf{a}_2 & \mathbf{a}_3\end{bma... ...athbf{a}}_1^* \\ \hat{\mathbf{a}}_2^* \\ \hat{\mathbf{a}}_3^* \end{bmatrix} \;,$$

where the superscript $*$ denotes conjugate transpose. Note that the eigenvalues $\nu_1^2\,,\nu_2^2$ and $\nu_3^2$ , as well as their corresponding eigenvectors, are complex-valued. Since $\mathbf{Q}$ is square and $\mathbf{Q}^{-1} \mathbf{Q} = \mathbf{I}$ :

$$\begin{cases} & \hat{\mathbf{a}_i}^* \mathbf{a}_j = 1 \;, i ... ...{\mathbf{a}_i}^* \mathbf{a}_j = 0 \;, i \neq j \;. \end{cases}$$

Physically, this means in strongly heterogeneous media, wave mode separation cannot be clearly defined due to wave mode conversion, which is expressed as the imaginary part of eigenvalues. The polarization directions are no longer orthogonal in the real sense - they become orthogonal complex vectors.

The situation stated above does not prevent our framework from performing the correct wavefield extrapolation. The wave extrapolation operator then becomes

$$e^{i\Phi \Delta t} \hat{\mathbf{u}}(k) = \sum\limits_{i=1,2,3} e^{i \nu_i \Delta t} \mathbf{a}_i \hat{\mathbf{a}}_i^* \hat{\mathbf{u}}$$ (27)

$$= \begin{bmatrix}\mathbf{a}_1 & \mathbf{a}_2 & \mathbf{a}_3 \end{bm... ... \\ \hat{\mathbf{a}}_2^* \\ \hat{\mathbf{a}}_3^* \end{bmatrix} \hat{\mathbf{u}}$$

$$= \begin{bmatrix} \hat{s}_{xx} & \hat{s}_{xy} & \hat{s}_{xz} \\ ... ...at{s}_{zx} & \hat{s}_{zy} & \hat{s}_{zz} \\ \end{bmatrix} \hat{\mathbf{u}} \;,$$

where

$$\hat{s}_{xx} = e^{i \nu_1 \Delta t} a_{1_x} \hat{a}_{1_x} + e^{i \nu_2 \Delta t} a_{2_x} \hat{a}_{2_x} + e^{i \nu_3 \Delta t} a_{3_x} \hat{a}_{3_x} \;,$$ (28)

$$\hat{s}_{xy} = e^{i \nu_1 \Delta t} a_{1_x} \hat{a}_{1_y} + e^{i \nu_2 \Delta t} a_{2_x} \hat{a}_{2_y} + e^{i \nu_3 \Delta t} a_{3_x} \hat{a}_{3_y} \;,$$

$$\hat{s}_{xz} = e^{i \nu_1 \Delta t} a_{1_x} \hat{a}_{1_z} + e^{i \nu_2 \Delta t} a_{2_x} \hat{a}_{2_z} + e^{i \nu_3 \Delta t} a_{3_x} \hat{a}_{3_z} \;,$$

$$\hat{s}_{yx} = e^{i \nu_1 \Delta t} a_{1_y} \hat{a}_{1_x} + e^{i \nu_2 \Delta t} a_{2_y} \hat{a}_{2_x} + e^{i \nu_3 \Delta t} a_{3_y} \hat{a}_{3_x} \;,$$

$$\hat{s}_{yy} = e^{i \nu_1 \Delta t} a_{1_y} \hat{a}_{1_y} + e^{i \nu_2 \Delta t} a_{2_y} \hat{a}_{2_y} + e^{i \nu_3 \Delta t} a_{3_y} \hat{a}_{3_y} \;,$$

$$\hat{s}_{yz} = e^{i \nu_1 \Delta t} a_{1_y} \hat{a}_{1_z} + e^{i \nu_2 \Delta t} a_{2_y} \hat{a}_{2_z} + e^{i \nu_3 \Delta t} a_{3_y} \hat{a}_{3_z} \;,$$

$$\hat{s}_{zx} = e^{i \nu_1 \Delta t} a_{1_z} \hat{a}_{1_x} + e^{i \nu_2 \Delta t} a_{2_z} \hat{a}_{2_x} + e^{i \nu_3 \Delta t} a_{3_z} \hat{a}_{3_x} \;,$$

$$\hat{s}_{zy} = e^{i \nu_1 \Delta t} a_{1_z} \hat{a}_{1_y} + e^{i \nu_2 \Delta t} a_{2_z} \hat{a}_{2_y} + e^{i \nu_3 \Delta t} a_{3_z} \hat{a}_{3_y} \;,$$

$$\hat{s}_{zz} = e^{i \nu_1 \Delta t} a_{1_z} \hat{a}_{1_z} + e^{i \nu_2 \Delta t} a_{2_z} \hat{a}_{2_z} + e^{i \nu_3 \Delta t} a_{3_z} \hat{a}_{3_z} \;.$$

Comparing equation 27 with equation 16, we can see that because the eigenvalues that appear in the exponent in equation 27 are no longer real-valued, the corresponding two-step operators to that in equation 20 cannot be constructed using Euler's formula.

Recursive integral time extrapolation of elastic waves using low-rank symbol approximation