finite difference solutions of the wave equation

For simplicity, I use the acoustic wave equation (11), which is second order in t, as opposed to equation (7) [fourth order in derivates of t]. Adding a force function, f(x,y,z,t), to equation (11) yields

   $$\begin{eqnarray} \frac{\partial^2 P}{\partial t^2} = (1+2 \eta) v^2 \left(\frac{... ...\frac{\partial^4 F}{\partial y^2 \partial z^2} \right) +f.\,\,\,\,\end{eqnarray}$$

In addition, we must solve

   $$\begin{eqnarray} P=\frac{\partial^2 F}{\partial t^2}\end{eqnarray}$$

for F along with solving equation (20).

Using the second-order finite-difference approach, we can approximate partial derivatives with difference equations, as follows

$$\frac{\partial^2 P}{\partial t^2} \rightarrow \frac{P_{i-1}- 2P_{i}+P_{i+1}}{{\Delta t}^2},$$

and

$$\frac{\partial^4 F}{\partial x^2 \partial z^2} \rightarrow ... ..._{i+1,j-1}+F_{i+1,j+1}+F_{i-1,j+1}}{{\Delta x}^2 {\Delta z}^2},$$

where i and j are indexes corresponding to steps in position or time. The later equation is derived from the 2-D convolution of the second-order approximation for each of x and z derivatives, and it is conveniently of a fourth-order accuracy.

Therefore, F in equation (21) is calculated using the following recursive formula:

$$F_{i+1}(x,y,z) = 2 F_{i}(x,y,z)-F_{i-1}(x,y,z) + {\Delta t}^2 P_{i}(x,y,z),$$

and P is calculated using

$$P_{i+1}(x,y,z) = 2 P_{i}(x,y,z)-P_{i-1}(x,y,z) + {\Delta t}^2 \left(\frac{\partial^2 P}{\partial t^2} \right)_i,$$

where $\frac{\partial^2 P}{\partial t^2}$ is the finite-difference approximation of the new Laplacian equation (20).

Higher-order finite-difference approximations (4th order) can be used with equation (7) directly. The set of finite-difference equations (20-21) are subjected to the same constraints and rules used in the isotropic case (such as the CFL condition) to avoid numerical dispersion and instability.

 

fddata
Figure 2 The wavefield at 1 s resulting from a source at center of the section for a VTI medium with $\eta$=0.4 (left), and for an isotropic medium (right). In both cases, the velocity (vertical velocity in VTI media) is 1000 m/s. The solid black curves are the solutions of the eikonal equation for both media with the shear wave velocity set to zero, and the dashed curves (not apparent) are the solutions when the shear wave velocity equals half the P-wave velocity. The two curves coincide in the isotropic case and practically coincide in the VTI case.

Figure 2 shows the wavefield at time 1 second caused by an impulse force excited at time 0. On the left side, the medium is homogeneous and VTI with v_v=2 km/s, v=2 km/s, and $\eta$=0.4. On the right side, the medium is isotropic with v_v=2 km/s, v=2 km/s, and $\eta$=0. Both wavefields are calculated using the second-order finite difference applied to the new acoustic wave equation. In the case of VTI, an additional wave type appears in the section and travels at a speed that is lower than the P-wave velocity. This artifact is the additional solution, mentioned earlier, that behaves like a wave for positive $\eta$ and exponentially decays or grows for negative $\eta$.Such an artifact does not appear in the solution for an isotropic medium. Therefore, we may want to place the source in an isotropic layer and take advantage of the evanescent nature of this wave in isotropic media. The black curves in Figure 2 correspond to solutions of the eikonal equation. The solid curves correspond to the solutions that use the acoustic assumption, in equation (16), in which the shear wave velocity equals zero, while the dashed curves correspond to a shear wave velocity equal to half the P-wave velocity. As expected, in the isotropic case both curves exactly coincide and are therefore indistinguishable. In the VTI case, differences between the two curves exist, but are hardly noticeable. The independence of the eikonal equation on the shear wave velocity in VTI medium is in agreement with the results I obtained in an earlier study 1997.

 

fddatatiel Figure 3 The z-component of the elastic wavefield at 1 s caused by a source at time 0 located at the center. The VTI model is the same as in Figure 2, with v=1000 m/s and $\eta$=0.4. The solid curve corresponds to the solution of the acoustic eikonal equation for the same medium.

Figure 3 shows the z-component of the elastic wavefield (computed using the elastic wave equation) for the same model used in Figure 2. The solid curve is the solution of the acoustic eikonal equation. Kinematically, for P-waves, the acoustic and elastic wavefields are similar. Dynamically, they differ considerably; the elastic wavefield includes S-waves (the slower wave), here with triplication, a phenomenon common to S-waves in strongly anisotropic media.