Numerical Implementation

Equations 3.11 and those in Appendix A are written in the wavenumber domain and are easily evaluated there. It would be difficult to derive finite-difference (FD) operators that correspond to all of the wavenumber terms. Rather than using a FD implementation on regular or staggered grids, I simply evaluate equations as written. Each of the seven terms required in 2D (just the wavenumber operator multiplied by the wavefield) is individually evaluated, inverse Fourier transformed back to space and then multiplied by the spatially variable coefficient term before summation. I note that these seven terms (twenty-one terms in 3D) can be evaluated in parallel to reduce the elapsed computation time. Once the spatial operators have been evaluated, equations 3.11 can be advanced in time using the rapid expansion method (REM). The REM for reverse time migration problem was proposed by Pestana and Stoffa (2010) for the isotropic case and has been successfully applied to the VTI case as well (Pestana et al., 2011). It is based on a Chebyshev orthogonal polynomial expansion (Tal-Ezer et al., 1987) and was applied to seismic modeling by Kosloff et al. (1989). REM is a good approach that in practice is numerically stable and does not show divergent behavior even in models with strong contrasts (Carcione et al., 2002), and when combined with the pseudospectral spatial operators, the wave propagation will be free of numerical dispersion.

Figure 3.1 shows a wavelet comparison after traveling 8 s in a homogeneous isotropic media using REM and FD method (2nd-order in time, 8th-order in space). It is obvious that 8th-order FD still causes visible dispersion while REM preserves the wavelet quite well. In addition, REM remains stable for any $\Delta t$ since REM is an algorithm that closely approximates the analytical solution.

Figure: Wavefield snapshots at time $t=1~s$ in a homogeneous TTI medium ( $V_{p_0}=3000~m/s$ , $\epsilon =0.24$ , $\delta =0.1$ , $\theta =45^\circ$ ) computed using different TTI equations. (a) is generated from the TTI coupled equations with $V_{s_0}=0$ . (b) is the same as (a) but with a nonzero SV wave velocity. (c) and (d) are the results from the solution of the TTI decoupled P and SV wave equations, respectively.