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

Numerical examples

In this section, we use various synthetic models to test the one-step version of the proposed low-rank RITE method.

ORTw-fd-1,ORTw-ps-1
Figure 1. Wavefield snapshot of wavefield propagation in a two-layer orthorhombic model using the finite-difference method (a) and the pseudospectral method (b) with a step size of $1\;ms$ .

ORTw-lr-1,ORTw-lr-2,ORTw-lr-4,ORTw-lr-8
Figure 2. Wavefield snapshot of wavefield propagation in a two-layer orthorhombic model using the proposed low-rank RITE method with a step size of $1\;ms$ (a), $2\;ms$ (b), $4\;ms$ (c) and $8\;ms$ (d).

We first test the accuracy and stability of the proposed method in a two-layer orthorhombic model, and compare it with the conventional finite-difference and pseudospectral methods. We construct a two-layer orthorhombic model on a $100\times100\times100$ grid with the density-normalized stiffness tensor (in $km^2/s^2$ ) of the first layer given in equation 41 (Schoenberg and Helbig, 1997), and the second layer being a scaled version of the first layer by a factor of $1.8$ . We apply triangle smoothing in the vertical direction with a radius of $4$ . The spatial sampling rate in all directions of the grid is $10\;m$ . A displacement source is oriented at $45^{\circ}$ tilt and $45^{\circ}$ azimuth and injected at $x=0.5\;km,\;y=0.5\;km\;,z=0.4\;km$ . The source wavelet has a peak frequency of $35\;Hz$ . A wavefield snapshot shown in Figure 3 was taken at $0.18\;s$ . The wavefield that was modeled by eighth-order finite-difference in space and second-order finite-difference in time using a time-step size of $1\;ms$ suffers from strong dispersion artifacts (Figure 1a), while the wavefield modeled by pseudospectral method in space and second-order finite-difference in time is free of spatial dispersions. However, increasing the time step size to $2\;ms$ leads to numerical instability of both methods, and thus the results are not shown. In comparison, the proposed low-rank RITE method shows an accurate result free from dispersion and instability using increasing time step sizes of $1\;ms$ , $2\;ms$ , $4\;ms$ and $8\;ms$ (Figure 3).

$$\begin{bmatrix}9 & 3.6 & 2.25 & 0 & 0 & 0 \\ 3.6 & 9.84 & 2.4 & 0... ...2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1.6 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2.182 \end{bmatrix}$$ (41)

ORTw-lr-1,ORTw-lr0-1
Figure 3. Wavefield snapshot of wavefield propagation in a two-layer orthorhombic model using the low-rank RITE method with (a) and without (b) accounting for the gradient of stiffnesses.

inter-lr,inter-ps,inter-pslr
Figure 4. Interleaved shot gathers extracted at $Z=100\;m$ between wavefield modeled by (a) the low-rank method with and without including the stiffness gradient terms; (b) the pseudo-spectral method with and without including the stiffness gradient terms; and (c) the pseudo-spectral method and the low-rank method, both including the stiffness gradient terms. The odd-numbered traces correspond to the methods mentioned first.

Next, to test wave extrapolation accuracy at a medium interface, where strong heterogeneity occurs, we use the same two-layer orthorhombic model without smoothing in the vertical direction. Figure 3 shows wavefield snapshots taken at $t=0.18\;s$ . Because of the strong contrast at the medium interface, the transmitted and reflected waves calculated by the proposed method including the stiffness gradient terms demonstrates noticeable amplitude and phase differences compared with the waves calculated by ignoring the gradient terms. Figure 4a shows interleaved shot gathers extracted at $Z=100\;m$ , which compares the modeled data by the low-rank method and the pseudo-spectral method, both with and without accounting for the stiffness gradient terms. We can observe that for either the low-rank or the pseudo-spectral methods, modeled direct P- and S-arrivals are identical while reflected P- and S-arrivals show obvious difference in both amplitude and phase due to the missing of stiffness gradient terms. To validate our observation, we also carry out wave extrapolation using the pseudo-spectral method by solving first- and second-order elastic wave equations. Their shot gather, shown in Figure 4b, demonstrates the same phenomena. Finally, in Figure 4c, we interleave the amplitude-normalized shot gathers from low-rank RITE solution including the stiffness gradients, togethr with a pseudo-spectral solution of the first-order equations. They demonstrate almost identical amplitude and phase behavior. In terms of computational cost, the numerical ranks of the wave extrapolation matrices are $2$ when the stiffness gradient terms are not included, while the ranks increase to $5$ when the gradient terms are included.

TRIw-lr-p,TRIw-lr-s
Figure 5. Wavefield snapshot of (a) P-wave and (b) coupled S-waves propagation in a perturbed triclinic model using the proposed low-rank RITE method with a step size of $1\;ms$ .

In the third example, we test the ability of the proposed method to model decomposed wave mode propagation in a more general anisotropic medium, a triclinic medium. We use the lab measurements from Mah and Schmitt (2003) as the background model, with the stiffness tensors (in GPa) of the triclinic model shown in equation 42, which is then normalized by a mass density of $1.395\;kg/m^3$ . We then add spherical perturbations centered at $Z=400\;m, X=500\;m, Y=500\;m$ to all the density normalized stiffness tensors to make the models mildly heterogeneous. Using the same mesh size as the previous example, and a time step size of $1\;ms$ , the proposed method is capable of accurately and separately propagating P- and S-waves in the triclinic medium (Figure 5) using equation 16. The two S-wave modes are propagated together since they do not decouple easily. Both the P- and coupled S-waves are free of numerical dispersion and instablity. The numercial ranks of the wave extrapolation matrices for this test are $3$ .

$$\begin{bmatrix}14.9 & 6.3 & 5.2 & 0.7 & 0.9 & -0.5 \\ 6.3 & 14.9 ... ...5 & 0.8 & -0.1 & 3.0 & 0.2 \\ -0.5 & -0.4 & 0.1 & 0.1 & 0.2 & 3.7 \end{bmatrix}$$ (42)

ORTc-11-f
Figure 6. (a) $C11$ of an orthorhommic model based on the French model.

ORTw-lr-p-x,snap-p-x-cut
Figure 7. (a) The x-component of decoupled P-wavefield at $t=0.39s$ . (b) A portion of the wavefield centered at $Z=0.64\;km, X=0.89\;m, Y=0.93\;m$ .

freq-p-x,freq-p-x2
Figure 8. Wavenumber amplitude spectrum of (a) the analytical wavefield, and (b) its real part.

snap-p-x-up,snap-s-x-up,snap-p-x-dn,snap-s-x-dn
Figure 9. The decomposition of wavefield into up- and down-going directions. (a) The up-going P-wave, and (b) the up-going S-wave. (c) The down-going P-wave. (d) The down-going S-wave.

In the last example, we demonstrate one advantage of the proposed method for seismic imaging applications, the accessibility of directional information of wavefield. For a real-valued wavefield, the propagation direction cannot be uniquely defined in the wavenumber domain because of its symmetric wavenumber spectrum about zero (Hu et al., 2016). This is due to the fact that the real-valued wavefield contains both positive and negative frequency components. The analytical wavefield, on the other hand, has decoupled positive and negative frequency components, so its wavenumber spectrum reveals the direction of wavefield propagation. To demonstrate this fact, we construct an orthorhombic model based on the French model (French, 1974), with its $C11$ component shown in Figure 6. We use equation 16 to calculate a decoupled P-wavefield, whose x-component at $t=0.39\;s$ is shown in Figure 7a. We select a portion of the wavefield centered at $Z=0.4\;km, X=0.58\;m, Y=0.58\;m$ propagating in the direction negative along all three axes, as shown in Figure 7b. Figure 8a shows the wavenumber amplitude spectrum of the wavefield, and its energy are all distributed in the positive octant. One the other hand, the wavenumber amplitude spectrum of the real part of the wavefield shows symmetric distribution of energy about the origin, as demonstrated by Figure 8b. Finally, we decouple the upward-traveling and downward-traveling wavefield based on the sign $k_z$ , and the downward-traveling P- and S-waves are shown in Figure 9c and 9d.

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