# Accurate Pseudo-Spectral Acoustic Wave Modelling with Time Dispersion Elimination
## 2 Methodology
\(\partial p/\partial x\) and \(v_{x}\) are calculated at grid points staggered in the \(x\) direction (five-pointed stars), \(\partial p/\partial y\) and \(v_{y}\) are evaluated at grid points staggered in the \(y\) direction (triangles). The remaining variables are evaluated on a regular grid (dots). Equation (4) is the baseline to develop the proposed method. In Equation (4), although the spatial derivative is accurately calculated using the Fourier transform, the temporal derivative is computed with a finite difference approach which typically has larger approximation errors than the Fourier transform. To eliminate such temporal numerical dispersion, [18] proposed using higher-order temporal finite difference approximation which is more time consuming than the commonly used second order approximation in the pseudo-spectral method. [19] present forward and inverse time dispersion transform to reduce temporal dispersion in the non-staggered grid-based modeling method. Here we firstly combine the forward and inverse time dispersion transform with the staggered grid-based pseudo-spectral method. Our present method not only retains the advantages of the staggered grid-based pseudo-spectral method such as coarser spatial sampling, higher spatial derivative approximation accuracy, and non-casual ringing artifacts removal, but also achieves higher temporal derivative approximation accuracy due to the adoption of the additional time-dispersion elimination term. Several numerical examples are demonstrated to verify the effectiveness of our method. Starting from the Taylor series approximation of the second order temporal finite difference operator [3],  
\[\frac{p(t-\Delta t)-2p(t)+p(t+\Delta t)}{\Delta t^{2}}=\sum_{k=1}^{\infty}\frac {2\Delta t^{2(k-1)}}{(2k)!}\frac{\partial^{2}p(t)}{\partial t^{2}} \tag{5}\]  
where \(p(t)\) represents the acoustic wavefield with time index \(t\), here the spatial index of wavefield is omitted for simplicity. For the Fourier transform on the right-hand side of Equation (5), we have  
\[F\left(\sum_{k=1}^{\infty}\frac{2\Delta t^{2(k-1)}}{(2k)!}\frac{\partial^{2}p (t)}{\partial t^{2}}\right)=-\theta^{2}\omega^{2}P(\omega) \tag{6}\]  
with  
\[\theta^{2}=\left(\frac{2\sin\!\left(\frac{\omega\Delta t}{2}\right)}{\omega \Delta t}\right)^{2} \tag{7}\]  
here symbol \(F\) stands for the forward Fourier transform. Reminds that the accurate second order time derivative in the frequency domain can be represented as \(-\omega^{2}\), then \(\theta\) can be the approximation error between the accurate and finite difference approximated temporal derivative. It is also clear that the temporal dispersion factor depends only on the time step and frequency, and is free from the medium parameter or spatial accuracy [25]. Based on factor (7), a transform that removes the undesired temporal dispersion can be defined, and has the form of  
\[p(t)=F^{-1}\left(\int\limits_{-\infty}^{+\infty}\hat{p}(t)e^{-i2\sin^{-1}( \frac{\omega\Delta t}{2})\frac{t}{\Delta t}\,dt}\right) \tag{8}\]  
here \(\hat{p}(t)\) represents the recorded wavefield trace with temporal dispersion errors, and \(p(t)\) is the recorded wavefield trace with temporal dispersion removal. \(F^{-1}\) stands for the inverse Fourier transform. Using a slightly different derivation, [19] proposed an analogous transform of Equation (8), which is referred to as the inverse time dispersion transform (ITDT) to reduce temporal dispersion. They also adopted this transform in a non-staggered grid-based modeling method. Note that to remove the temporal dispersion from the seismic record, an adjoint operator of Equation (8) should be pre-applied on the injection source wavelet for numerical modeling to add temporal dispersion on the seismic record [19]. Here, we firstly combine Equations (4) and (8) to establish a staggered grid-based, temporal dispersion-eliminated, pseudo-spectral method for better performance.