# Eliminating time dispersion from seismic wave modeling

Eliminating time dispersion from seismic wave modeling Summary We derive an expression for the error introduced by the second-order accurate temporal finite-difference (FD) operator, as present in the FD, pseudospectral and spectral element methods for seismic wave modeling applied to time-invariant media. The ‘time-dispersion’ error speeds up the signal as a function of frequency and time step only. Time dispersion is thus independent of the propagation path, medium or spatial modeling error. We derive two transforms to either add or remove time dispersion from synthetic seismograms after a simulation. The transforms are compared to previous related work and demonstrated on wave modeling in acoustic as well as elastic media. In addition, an application to imaging is shown. The transforms enable accurate computation of synthetic seismograms at reduced cost, benefitting modeling applications in both exploration and global seismology. Numerical modelling, Computational seismology, Wave propagation 1 INTRODUCTION The understanding of the earth's subsurface from seismology observations is increasingly held back more by computing power rather than available data. In particular, propagating seismic data through synthetic models dominates the cost of wave equation-based imaging methods such as reverse time migration (RTM, Baysal et al.1983) and full-waveform inversion (Tarantola 1984). Coarsely discretizing the wavefield typically curbs the modeling costs, but introduces ‘numerical dispersion’ errors. The numerical dispersion distorts the modeled wavefield, affecting particularly the high-frequency content of the data. In practice, thus, only the low-frequency end of the full data bandwidth can be imaged accurately at reasonable cost. Breaking the barrier between modeling costs and accuracy is therefore important for many applications. In particular, modeling accurate wavefields with coarse discretization would enable accurate imaging at low cost. Accurate wavefield modeling implies that no errors arise from the space and time discretization; a topic that has been addressed in various ways by researchers in the past. For instance, the spatial and temporal error may be combined such that they partially cancel each other (see e.g. Etgen & Brandsberg-Dahl 2009; Wang et al.2014). The drawback of such approaches is the need for medium- and simulation-dependent stencils (Stork 2013). Instead, it is also possible to minimize the spatial and temporal error separately, as discussed below. Modeling seismic waves involve the computation of spatial gradients of the wavefield and of the medium. In exploration seismology, these spatial gradients are often computed with finite-difference (FD) approximations. The FD method (Moczo et al.2014) replaces derivatives in the wave equation with FD approximations of limited accuracy. A fine spatial discretization improves the accuracy of the FD method at significant computing cost. For example, reducing the spacing in 3-D scales to the power four with computing cost, for equal signal bandwidth. Instead, it is common to use higher order FD operators to reduce errors, even at a coarse discretization (see e.g. Song et al.2013; Liu 2013, 2014; Yao et al.2016). In the limit of using the entire wavefield to compute derivatives, the pseudospectral (PS) method achieves very high accuracy (Kosloff & Baysal 1982). Another high-order accurate spatial operator, derived from the weak formulation of the wave equation, is the spectral element method (SEM; Tromp et al.2008). SEM is used in global seismology to model the phases of seismograms with high accuracy, but the temporal error remains, and its removal would greatly benefit seismologists. Within FD, PS and SEM modeling, the time integration is carried out with explicit time stepping, an FD approximation. Accurate time stepping is then achieved using short steps in time, often at large computing cost. Alternatively, higher order Lax–Wendroff corrections can be used to reduce the errors (e.g. Dablain 1986; Blanch & Robertsson 1997; Dumbser et al.2007). Such correction schemes allow long time steps, although the Courant–Friedrichs–Lewy (CFL) condition defines an upper time-step limit for stable simulations. If short time steps are required for a stable simulation, large temporal errors are naturally prevented. However, temporal errors accumulate and may still become significant for a given modeling purpose, as the errors increase with number of propagated wavelengths and the associated frequency, which we will further address in the discussion section (Stork 2013). Additionally, optimized Lax–Wendroff schemes increase the stable time-step upper limit two- or threefold (Soubaras & Zhang 2008; Chu et al.2009; Liu et al.2014; Amundsen & Pedersen 2017), putting the temporal error at the forefront again. Stork (2013) took a radically different approach by suggesting that the temporal error can be filtered from seismograms at modest additional cost. A 1-D homogeneous example to illustrate the concept put forward by Stork is shown in Fig. 1. The procedure consists of three steps. First, a train of wavelets is modeled in an arbitrary model (e.g. a homogeneous model) to obtain a trace which gains temporal errors with increasing time (Fig. 1b). Second, using a simulation with very small steps in time, or with the use of the source wavelet itself (both options contain negligible or no temporal errors), a deconvolution is applied to the dispersed wavelets to create a ‘bank’ of filters that can remove the temporal error at various times (Fig. 1c). Finally, the desired wavefield simulation is carried out in the actual seismic domain of interest (which may be of arbitrarily high complexity and spatial dimension), and the previously computed filters are applied to remove the temporal error. Details of this approach can be found in Stork (2013), Li et al. (2013, 2016), Dai et al. (2014) and Anderson et al. (2015). A correct implementation of the filter is tedious, requiring many dispersed wavelets corresponding to the particular time-step characteristics of that simulation, and interpolation to correct wavelets at intermediate times. A simpler alternative using analytical filters was proposed by Wang & Xu (2015a,b), and later rederived in Qin et al. (2017), Mittet (2017) and Xu et al. (2017). The suggested analytical filters apply an altered discrete Fourier transform to seismograms which eliminates the time-stepping error in one operation, thus do not require the creation a filter bank ahead of the simulation or interpolation of the filter operation. However, adding and correspondingly removing the error with the proposed filters does not reproduce the input signal fully (see Section 2.5 below), suggesting that the filters are incomplete. Figure 1. View largeDownload slide The concept of a time-dispersion filter, cf. Stork (2013), illustrated in a 1-D homogeneous medium. (a) A simulation is created with a given source wavelet and time step. (b) In the absence of spatial effects distorting the wavelet, for example by carrying out a pseudospectral simulation in a 1-D homogeneous medium, a trace is gained with dispersion due to the time-stepping error only. (c) A filter bank is created which undoes the time–frequency-varying phase shift due to the temporal error. The created filter bank can then also be applied to traces modeled in complex media of higher dimensions. Figure 1. View largeDownload slide The concept of a time-dispersion filter, cf. Stork (2013), illustrated in a 1-D homogeneous medium. (a) A simulation is created with a given source wavelet and time step. (b) In the absence of spatial effects distorting the wavelet, for example by carrying out a pseudospectral simulation in a 1-D homogeneous medium, a trace is gained with dispersion due to the time-stepping error only. (c) A filter bank is created which undoes the time–frequency-varying phase shift due to the temporal error. The created filter bank can then also be applied to traces modeled in complex media of higher dimensions. To this date, the time-stepping error has only been derived for constant-velocity acoustic media, and it has only been conjectured that the error is independent of the spatial medium (e.g. Stork 2013). In this paper, we solve these outstanding questions by deriving an expression for the time-stepping error that is independent of the spatial medium and we refer to it as ‘time dispersion’. We show that time dispersion must be added to the injection wavelet, to model the correct kinematics and to correctly remove time dispersion afterwards. We derive a new exact filter pair to add and remove time dispersion from traces, which is fully invertible. We demonstrate the proposed method on acoustic and anisotropic elastic media, using FD, PS and SEM modeling. The results are valid for solving the second-order wave equation, or the system of first-order partial differential equations (Virieux 1986). Finally, the method is demonstrated on a simple RTM example, to show that the results apply to generic seismic modeling problems. 2 THEORY In a general form, the wave equation can be written as:   $$\frac{\partial ^2 \mathbf {u}(\mathbf {x},t)}{\partial t^2} = L(\mathbf {x})\mathbf {u}(\mathbf {x},t) + \mathbf {f}(\mathbf {x},t),$$ (1)where u(x, t) represents the particle displacement at position x and time t, L represents the elastic potential term and f the source function. We assume that the medium is time invariant, but may be heterogenous and anisotropic. We will discuss the implications of time dispersion for viscoelasticity in Section 5.2. 2.1 Finite-difference operator Numerical wave propagation methods such as FD, PS and SEM, typically compute the second-order time derivative in the wave equation using an FD approximation. In particular, a second-order accurate FD operator derived from a Taylor series is often used:   $$\frac{\mathbf {u}({t\!-\!\Delta t})-2\mathbf {u}({t}) + \mathbf {u}({t\!+\!\Delta t})}{\Delta t^2} = \sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}},$$ (2)where Δt denotes the time step. Disregarding O(Δt2) and higher order terms from eq. (2), thus assuming that the operator represents just the second-order time derivative, introduces a numerical error. Accurate simulations are typically achieved using short steps in time or by using higher order time-stepping schemes (Dablain 1986). Note that, although eq. (2) corresponds to an approximation of a second-order time derivative, the results described in the following sections apply equally to a system of first-order partial differential equations solved with a leap-frog time-stepping scheme (Virieux 1986) as shown in Appendix A. The theory in this paper may also be used to derive similar correction filters if higher order time-stepping schemes are used, for example to allow longer steps in time. Appendix B demonstrates this procedure for a higher order time-stepping scheme (Soubaras & Zhang 2008; Amundsen & Pedersen 2017). 2.2 Describing time dispersion Denoting the Fourier transform operator by $$\mathcal {F}[\cdot ]$$ and its inverse by $$\mathcal {F}^{-1}[\cdot ]$$, we begin by transforming the right-hand side of eq. (2) to the frequency domain:   \begin{eqnarray} \mathcal {F}\left[\sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}}\right]= -\theta ^2\omega ^2 {\tilde{\mathbf {u}}}(\omega ). \end{eqnarray} (3)The right-hand side of eq. (3) corresponds to the second-order time derivative multiplied with an additional factor θ2 corresponding to the undesired remaining Taylor expansion terms:   $$\theta ^2\!=\!\sum _{k=1}^\infty \!\frac{(-i\omega \!\Delta t)^{2(k-1)}}{\frac{1}{2}(2k)!}\!=\!1-\frac{\Delta t^2\omega ^2}{\frac{1}{2}4!}+\ldots \!=\!\left(\!\frac{\sin (\!\frac{\omega \Delta t}{2}\!)}{\frac{1}{2}\omega \Delta t}\!\right)^2\!.$$ (4) We insert the second-order accurate FD operator in time in eq. (1), and find that:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},{t\!-\!\Delta t})\!-\!2\mathbf {u}(\mathbf {x},{t})\!+\!\mathbf {u}(\mathbf {x},{t\!+\!\Delta t})}{\Delta t^2}\!\approx \!L(\mathbf {x}){\mathbf {u}(\mathbf {x},t)}\!+\mathbf {f}( \mathbf {x},\!t), \end{eqnarray} (5)thus represents a modified wave equation, solving exactly:   $$\frac{\partial ^2 \mathbf {\hat{u}}(\mathbf {x},t)}{\partial t^2} = \mathcal {F}^{-1}\Bigg [\frac{1}{\theta ^2}\mathcal {F}\left [L(\mathbf {x}){\mathbf {\hat{u}}(\mathbf {x},t)}+\mathbf {f}(\mathbf {x},t)\right ]\Bigg ].$$ (6)The formulation in eq. (6) follows from applying a forward and inverse Fourier transform on eq. (5), and eliminating θ from the left-hand side. We explicitly note that $$\mathbf {\hat{u}}(\mathbf {x},t)$$ is the time-dispersed solution to the FD system, without specifying anything about the spatial operator L or its accuracy. Solutions that satisfy the FD system are closely related to solutions to the wave equation, as is detailed in the section below. In Appendix C, we present a proof that the relation holds in the continuous case. We will show that eq. (6) describes an acceleration of the original wavefield. Consider first an illustrative example where we substitute t → 2t in the general wave equation (1):   \begin{eqnarray} \frac{\partial ^2 \mathbf {u}(\mathbf {x},2t)}{\partial (2t)^2} = L(\mathbf {x})\mathbf {u}(\mathbf {x},2t) + \mathbf {f}(\mathbf {x},2t), \end{eqnarray} (7)  \begin{eqnarray} \frac{\partial ^2 \mathbf {u}(\mathbf {x},2t)}{\partial t^2} = 2^2\left (L(\mathbf {x})\mathbf {u}(\mathbf {x},2t) + \mathbf {f}(\mathbf {x},2t)\right ). \end{eqnarray} (8)Eq. (8) computes the wavefield u(x, 2t) using the original time differential ∂t, by multiplying the spatial terms with factor 22. By analogy, eq. (5) similarly represents an accelerated wavefield if we denote $$\mathbf {u}(\mathbf {x},\hat{t})$$ as the accelerated solution and modify the source function $$\mathbf {f}(\mathbf {x},t)\rightarrow \mathbf {f}(\mathbf {x},\hat{t})$$:   $$\frac{\partial ^2 \mathbf {u}(\mathbf {x},\hat{t})}{\partial t^2} = \mathcal {F}^{-1}\Bigg [\frac{1}{\theta ^2}\mathcal {F}\left [L(\mathbf {x}){\mathbf {u}(\mathbf {x},\hat{t})}+\mathbf {f}(\mathbf {x},\hat{t})\right ]\Bigg ].$$ (9)Eq. (9) implies that the second-order temporal FD operator solves an accelerated version of the wave equation perfectly after proper pre-processing of the source wavelet. Would we use the original injection wavelet, the accelerated wavefield would correspondingly contain erroneous wavelengths, changing the kinematic and dynamic properties of the simulation. To our knowledge, this important observation is absent from previously published literature on time dispersion. As an explicit example, consider the heterogeneous acoustic wave equation:   $$\frac{\partial ^2 p(\mathbf {x},t)}{\partial t^2} = c^2(\mathbf {x})\rho (\mathbf {x})\nabla \cdot \rho ^{-1}(\mathbf {x})\nabla p(\mathbf {x},t),$$ (10)where ∇ represents the spatial derivative operator. Using second-order differencing in time solves the modified equation:   $$\frac{\partial ^2 p(\mathbf {x},\hat{t})}{\partial t^2}=\mathcal {F}^{-1}\Bigg [\frac{c(\mathbf {x})^2}{\theta ^2}\mathcal {F}\left [\rho (\mathbf {x})\nabla \cdot \rho ^{-1}(\mathbf {x})\nabla p(\mathbf {x},\hat{t})\right ]\Bigg ],$$ (11)The factor c(x)/θ takes the role of the effective velocity, that is, higher frequencies propagate faster than lower frequencies. This effect, caused by the second-order finite differencing in time, is the time dispersion shown in Fig. 2. Note that Dablain (1986) derived a similar expression for the time dispersion, assuming a homogenous acoustic medium and using plane-wave theory. However, this severely restrictive limitation is not required. The theory derived in this section remains valid for mode conversions, anisotropy or complex medium responses: signals are simply accelerated as a function of their frequency. Figure 2. View largeDownload slide The function 1/θ originating from the second-order FD operator represents a frequency-dependent velocity increase. The horizontal axis represents the range from 0 frequency up to f = 1/2Δt. Figure 2. View largeDownload slide The function 1/θ originating from the second-order FD operator represents a frequency-dependent velocity increase. The horizontal axis represents the range from 0 frequency up to f = 1/2Δt. 2.3 Time-dispersion transforms By recognizing time dispersion as an acceleration of the wavefield, we can consider it as a problem that only affects the time axis. Accelerating the source $$\mathbf {f}(\mathbf {x},t)\rightarrow \mathbf {f}(\mathbf {x},\hat{t})$$, thus adding time dispersion, or decelerating a recorded trace $$\mathbf {u}(\mathbf {x},\hat{t})\rightarrow \mathbf {u}(\mathbf {x},t)$$, thus removing time dispersion, can therefore be carried out using a filter operating on a single time-series. Adding dispersion to a signal of frequencies ω0 gives rise to the following spectrum, due to the acceleration factor 1/θ:   $$\omega = \frac{\omega _0}{\theta }.$$ (12)Note that this relation between the dispersed and original frequencies is independent of the medium. The factor θ only depends on the time step and frequency of the wavelet, describing time dispersion with no knowledge of the propagation path, medium, or even spatial accuracy. Substituting eq. (4) in eq. (12) provides a relation between the original and time-dispersed frequencies:   \begin{eqnarray} \omega _0 = \frac{2}{\Delta t}\sin \left(\frac{\omega \Delta t}{2} \right), \end{eqnarray} (13a)  \begin{eqnarray} \omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t}{2} \right) . \end{eqnarray} (13b) The spectral content of the recorded signal is thus identical to the correct signal, it is merely recorded at the wrong frequency. However, we can map the data back to the desired frequency. We propose a map between the correct and dispersed phase rotations of a recording by inserting eq. (13a) in a forward Fourier transform. We thus add dispersion to a trace u(t) as follows:   $$\hat{u}(t)= \mathcal {F}^{-1}\left[ \int _{-\infty }^\infty u(t) e^{-i2\sin \left(\frac{\omega \Delta t}{2}\right)\frac{t}{\Delta t}} \textrm{d}t\right].$$ (14)We will refer to this operation as the Forward Time Dispersion Transform (FTDT). Eq. (13b) allows us to map a temporally dispersed signal back to the correct phase rotations, in a similar fashion as the FTDT. We thus remove time dispersion from trace $$\hat{u}(t)$$ as follows:   \begin{eqnarray} {u}(t) = \mathcal {F}^{-1}\left[ \int _{-\infty }^\infty \hat{u}(t) e^{-i2\sin ^{-1}\left(\frac{\omega _0\Delta t}{2}\right)\frac{t}{\Delta t}} \textrm{d}t\right]. \end{eqnarray} (15)We will refer to this operation as the Inverse Time Dispersion Transform (ITDT). 2.4 Sampling a dispersed wavefield As we record the accelerated wavefield on the seismograms, we must carefully consider the sampling requirements. The Nyquist frequency fN = 1/(2Δt) represents the maximum frequency that can be recorded without aliasing. Considering a maximum signal acceleration of π/2 ≈ 1.6 (Fig. 2), we must sample the simulated field with an increased sampling rate fN, d = π/2 fN = π/(4Δt) to capture signals at their original Nyquist frequency. The support of the sin −1 function in eq. (13b) coincides exactly with this frequency range. 2.5 Comparison to filters derived in Wang & Xu (2015a) While the ITDT, as denoted by eq. (15), is equal to that found in Wang & Xu (2015a), the FTDT as given in eq. (14) is not. Wang & Xu (2015a) use for their FTDT, the adjoint of the ITDT as given in eq. (15). The adjoint of this modified Fourier transform does not provide a complete inverse, which leads to an incorrect phase and amplitude response. We demonstrate this from a pure signal-processing perspective. Applying time dispersion and correspondingly removing it, should return the same data as we merely map one frequency into another and then reverse the process. In Fig. 3, the filter pair from Wang & Xu (2015a) has been applied on a Ricker wavelet. The central frequency is 8 Hz with time step Δt = 15 ms. We note that the filter pair from Wang & Xu (2015a) does not reproduce the input wavelet, while the filter pair presented in this paper does. The effect becomes more notable at higher frequencies. Figure 3. View largeDownload slide Adding and correspondingly removing time dispersion should return the input data. The filters from Wang & Xu (2015a) do not achieve this, whereas the filters in this paper do. Wang & Xu (2015a) use for the FTDT an adjoint version of the ITDT, which is an incorrect filter. Figure 3. View largeDownload slide Adding and correspondingly removing time dispersion should return the input data. The filters from Wang & Xu (2015a) do not achieve this, whereas the filters in this paper do. Wang & Xu (2015a) use for the FTDT an adjoint version of the ITDT, which is an incorrect filter. 2.6 Implementation of time-dispersion filters Eqs (14) and (15) are discretized for sampled data. We implement the FTDT using a discrete Fourier transform and a standard inverse fast Fourier transform (IFFT), for a discrete signal u[n] of N samples recorded uniformly with a sampling interval Δt. The FTDT becomes:   $$\hat{u}[n] = {\rm IFFT}\!\left[ \sum _{j=0}^{N\!-\!1} u[j]e^{\!-\!i2\!\sin \left(\!\frac{\pi f}{N}\!\right) j } \right],\quad \, f\!\in \!\left\lbrace 0,..., \frac{N\!\!-\!1}{2}\right\rbrace .$$ (16)Similarly for the ITDT, we obtain:   $${u}[n] = {\rm IFFT}\!\left[ \sum _{j=0}^{N\!-\!1}\!\hat{u}[j]e^{\!-\!i2\!\sin ^{\!-\!1}\left(\!\frac{\pi f}{N}\!\right) j } \right]\!,\quad \, f\!\in \!\left\lbrace 0,..., \frac{N\!\!-\!1}{\pi }\right\rbrace .$$ (17)The frequency vector f is limited to the support of the sin −1 function, as noted in Section 2.4. Note that discretizing the transforms in this formation makes the operations independent of Δt. We can thus apply the transforms without any knowledge of simulation parameters. We zero-pad the signal to twice the length to avoid adverse edge effects from the Fourier transforms. A MATLAB implementation of the proposed FTDT and ITDT is provided with the digital supplements to this paper. 2.7 Application of time-dispersion filters Eqs (16) and (17) add and remove time dispersion from signals. This result holds for any time-invariant medium L(x), spatial error O(Δxn), and is valid for both solving the second-order wave equation as well as the corresponding system of first-order partial differential equations (Virieux 1986). To remove time dispersion from seismic traces, we thus follow the workflow from Fig. 4: Apply the FTDT to the source function f, as given by eq. (9). The pre-processing of the wavelet ensures that the kinematics of the wavefield is retained. From a numerical perspective, it means that we inject a wavelet with an identical amount of time dispersion as that accumulating in the wavefield. Run the standard seismic wave simulation application and record traces wherever desired. Apply the ITDT to recorded traces after the simulation is completed. This recovers the original, time-dispersion-free wavefield. Figure 4. View largeDownload slide The proposed workflow to eliminate time dispersion from synthetic seismograms. No changes are required to the standard seismic modeling application. Furthermore, the filters are independent of the medium, and depend only on the time stepping and frequency of the signals. Figure 4. View largeDownload slide The proposed workflow to eliminate time dispersion from synthetic seismograms. No changes are required to the standard seismic modeling application. Furthermore, the filters are independent of the medium, and depend only on the time stepping and frequency of the signals. 3 EXAMPLES In this section, we apply the workflow of Fig. 4 on FD, PS and SEM modeling results for various dimensions and models. The results are compared to analytical 1-D solutions, or against simulations using short steps in time, to show that equivalent accuracy is reached at reduced cost. 3.1 1-D acoustic examples We demonstrate the FTDT and ITDT with two simple 1-D examples. The first is a PS computed solution to the acoustic wave equation, free of spatial errors. In addition, a staggered second-order spatial FD operator is used to generate a solution that illustrates that the transforms are independent of spatial errors. We use a spatial step of Δx = 5 m, a temporal step of Δt = 15 ms and a constant velocity c = 2000 m s−1, with source and receiver separated by 1 km. We inject a Ricker wavelet of 35 Hz central frequency, with FTDT applied, thus the recorded wavelet has propagated about 20 wavelengths. The code to perform these simulations is provided with the digital supplements to this paper. 3.1.1 Pseudospectral spatial operator The PS method is spectrally accurate in space, and thus only accumulates time dispersion. Fig. 5 compares the PS solution with an analytical solution of the expected wavelet arrival. We observe that applying the FTDT to the analytical solution reproduces the time-dispersed solution, with a maximum error on the order of 10−10. The error is likely a combination of minor artefacts of the discrete Fourier transform in both the PS spatial operator and in the temporal dispersion transforms. Conversely, we construct a dispersion-free recording by applying the ITDT to the PS solution. Figure 5. View largeDownload slide Time dispersion owing to the pseudospectral method (blue) compared to the analytical solution (green). The time-dispersion transform workflow (Fig. 4) recreates the analytical solution (black dots), while applying the FTDT to the analytical solution recreates the signal computed with the PS method (red dots). The maximum difference between the analytical solution and the computed, corrected recording is on the order of 10−10. Figure 5. View largeDownload slide Time dispersion owing to the pseudospectral method (blue) compared to the analytical solution (green). The time-dispersion transform workflow (Fig. 4) recreates the analytical solution (black dots), while applying the FTDT to the analytical solution recreates the signal computed with the PS method (red dots). The maximum difference between the analytical solution and the computed, corrected recording is on the order of 10−10. 3.1.2 Second-order accurate staggered spatial operator The staggered second-order spatial FD operator introduces spatial errors in addition to time dispersion. Fig. 6 compares two seismograms using a long and short time step, differing by a factor 10 to create a ‘course’ and a ‘fine’ simulation in time. It can be concluded that the FTDT and ITDT correctly model and remove time dispersion, independent of the spatial errors. We observe that the error between the two solutions is on the order of 10−3, which predominantly can be explained by remaining time dispersion in the short time-step solution. Figure 6. View largeDownload slide Time dispersion in the second-order accurate spatial FD method (blue), compared to a short time-step solution (green). The time-dispersion transform workflow (Fig. 4) recreates the solution from the short time steps (black dots), while applying the FTDT to the short time-step solution recreates the signal computed with long time steps. All remaining ‘ringing’ from the spatial errors is left untouched by the transforms. Figure 6. View largeDownload slide Time dispersion in the second-order accurate spatial FD method (blue), compared to a short time-step solution (green). The time-dispersion transform workflow (Fig. 4) recreates the solution from the short time steps (black dots), while applying the FTDT to the short time-step solution recreates the signal computed with long time steps. All remaining ‘ringing’ from the spatial errors is left untouched by the transforms. 3.2 2-D elastic examples 3.2.1 2-D anisotropic elastic example using fourth-order accurate FD The time-dispersion transforms are not limited or influenced by the number of spatial dimensions. We demonstrate their performance on a 2-D elastic anisotropic model, using a staggered grid fourth-order spatial FD operator, to solve the first-order system of equations of the wave equation in a leap-frog fashion (Virieux 1986). The model is shown in Fig. 7, with a vertical transverse isotropy given by Thomsen's parameters δ = ε = 0.1 (Thomsen 1986). The model is discretized with Δx = Δz = 1 m, and an explosion-type Ricker wavelet of 80 Hz central frequency is used. Fig. 8 shows snapshots of the wavefield displaying a large variety of different wave types. As shown on a trace recorded close to the free surface, we fully recreate the short time-step solution from the coarse simulation. The difference between the two solutions is small for early arrivals, but becomes more apparent for later times. The code to perform this simulation is provided with the digital supplements of this paper. Figure 7. View largeDownload slide Elastic model used for the simulation shown in Fig. 8. A vertically transverse isotropy is additionally used, with Thomsen's parameters ε = δ = 0.1. The source and receiver positions are shown by the red star and blue triangle, respectively. Figure 7. View largeDownload slide Elastic model used for the simulation shown in Fig. 8. A vertically transverse isotropy is additionally used, with Thomsen's parameters ε = δ = 0.1. The source and receiver positions are shown by the red star and blue triangle, respectively. Figure 8. View largeDownload slide Snapshots of a 2-D elastic anisotropic simulation with a source at the red star and receiver at the blue triangle. The horizontal particle velocity recorded at the receiver is shown on the right. The long (Δt = 0.32 ms) and short (Δt = 0.032 ms) time-step simulations are compared, as well as the solution after applying the time-dispersion transforms. The solution after applying the time-dispersion transforms matches the short time-step solution for all times. Figure 8. View largeDownload slide Snapshots of a 2-D elastic anisotropic simulation with a source at the red star and receiver at the blue triangle. The horizontal particle velocity recorded at the receiver is shown on the right. The long (Δt = 0.32 ms) and short (Δt = 0.032 ms) time-step simulations are compared, as well as the solution after applying the time-dispersion transforms. The solution after applying the time-dispersion transforms matches the short time-step solution for all times. 3.2.2 2-D isotropic elastic example using the spectral element method We demonstrate the time-dispersion transforms on the SEM, using SPECFEM2D with its default isotropic elastic model shown in Fig. 9 (Tromp et al.2008). We inject a Ricker wavelet of 25 Hz central frequency, with FTDT applied. We compare solutions using a long time step (Δt = 1 ms) and a short time step (Δt = 0.1 ms). Fig. 10 shows that applying the ITDT to the long time-step solution correctly reproduces the solution from the short time step, where the presence of time dispersion is negligible. Figure 9. View largeDownload slide The SPECFEM 2-D (Tromp et al.2008) default four layer model as shown with four different colours, with as overlay the used mesh, as well as the vertical particle displacement at 0.38 s. The receivers are shown in black along the surface and the blue rhombus represents the trace that is zoomed-in on in Fig. 10. Figure 9. View largeDownload slide The SPECFEM 2-D (Tromp et al.2008) default four layer model as shown with four different colours, with as overlay the used mesh, as well as the vertical particle displacement at 0.38 s. The receivers are shown in black along the surface and the blue rhombus represents the trace that is zoomed-in on in Fig. 10. Figure 10. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1 ms) and a short (Δt = 0.1 ms) time-step 2-D elastic spectral element simulation using SPECFEM2D (Tromp et al.2008). The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. Figure 10. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1 ms) and a short (Δt = 0.1 ms) time-step 2-D elastic spectral element simulation using SPECFEM2D (Tromp et al.2008). The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. 3.3 3-D elastic example using high-order FD Spatial numerical errors can be limited using a high-order spatial FD operator. Here, we demonstrate the time-dispersion transforms for such a case on the 3-D elastic SEG/EAGE overthrust model (Aminzadeh 1996), shown in Fig. 11. We use a shear wave velocity of 0.6 times the P-wave velocity and derive the density model through Gardner's relation (Gardner et al.1974). We use a 20 point spatial operator with optimized coefficients (Liu 2014) and a second-order accurate FD operator in time to solve the system of first-order equations in the conventional leap-frog fashion (Virieux 1986). We inject a Ricker wavelet of 25 Hz central frequency, with FTDT applied. Fig. 12 shows that the ITDT applied to a long time-step seismogram removes time dispersion, thereby reproducing the short time-step solution where the presence of time dispersion is negligible. Figure 11. View largeDownload slide The P-wave velocities for the 3-D SEG/EAGE overthrust model (Aminzadeh 1996), the black lines denote the location of the slices projected on the 3-D cube. The receivers are denoted with blue (overlapping the black line), the black rhombus represents the trace that is zoomed-in on Fig. 12, the source is located at 100 m depth below the black star. Figure 11. View largeDownload slide The P-wave velocities for the 3-D SEG/EAGE overthrust model (Aminzadeh 1996), the black lines denote the location of the slices projected on the 3-D cube. The receivers are denoted with blue (overlapping the black line), the black rhombus represents the trace that is zoomed-in on Fig. 12, the source is located at 100 m depth below the black star. Figure 12. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1.6 ms) and a short (Δt = 0.16 ms) time-step 3-D elastic FD simulation. The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths, extracted from the full 3-D gather. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. Figure 12. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1.6 ms) and a short (Δt = 0.16 ms) time-step 3-D elastic FD simulation. The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths, extracted from the full 3-D gather. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. 4 APPLICATION TO REVERSE TIME MIGRATION RTM recovers subsurface reflectors by evaluating the zero-lag time correlation between a forward propagated source wavefield and a time-reversed backpropagated receiver wavefield. The wavefields accumulate time dispersion in opposing time directions, which blurs the migration. Stork (2013) notes that this blurring is mitigated by adding time dispersion to the receiver traces before time reversal – the time-reversal process removes the added time dispersion as the simulation runs backwards in time, and as a result the two wavefields contain equal time dispersion at all times. We will refer to such an RTM as a temporally spectral accurate (TSA) RTM. Effectively, we pre-process the data by adding time dispersion in accordance with the modified wave equation (9) that we compute numerically. The workflow is shown in Fig. 13. Figure 13. View largeDownload slide Workflow to eliminate time dispersion from reverse time migration. The ITDT and FTDT are each other's inverse, thus the dispersed synthetic traces may be used straight as input to the reverse run. Conversely, the recorded data are pre-processed to prepare it for the reverse-time simulation. Effectively, thus, only two processing steps are required. Figure 13. View largeDownload slide Workflow to eliminate time dispersion from reverse time migration. The ITDT and FTDT are each other's inverse, thus the dispersed synthetic traces may be used straight as input to the reverse run. Conversely, the recorded data are pre-processed to prepare it for the reverse-time simulation. Effectively, thus, only two processing steps are required. We show a simple example of the TSA RTM in an acoustic case with a second-order accurate in time and eight-order accurate in space-staggered FD simulation (Virieux 1986). Our model contains a shallow as well as a deep box-shaped target (velocity anomalies). The simplicity of the model makes the size of the errors more visible, as blurring of phases in complex media is less apparent when multiple reflectors blur. Fig. 14(a) corresponds to a case where we inject ‘recorded traces’ (numerical dispersion-free) and a forward modeled source wavefield (numerically dispersed) in a smooth velocity model, and reverse propagate the two simulations in tandem to produce an image by cross-correlation of snapshots. The image of the shallow target is well recovered (due to the short propagation paths of reflections in the shallow part of the model), but the deep target is imaged at the incorrect depth and with distorted phases. In Fig. 14(b), we have injected the time-reversed receiver data after applying the FTDT. The procedure images both the shallow and deep targets accurately. This result comes at no extra computing cost other than pre-processing the recorded traces with the altered forward Fourier transform and standard inverse (fast) Fourier transform constituting the FTDT. An important consequence is that relatively long time steps may be used in the RTM simulation without compromising on accuracy. Figure 14. View largeDownload slide Comparing a conventional RTM against a temporally spectral accurate (TSA) RTM. The boxes in red display the location of the targets (velocity anomalies). Top: zoom on shallow target. Bottom: zoom on deep target. (a) Injecting the data without applying dispersion results in a blurred image of the deep target at incorrect depth. (b) Injecting the data after applying the FTDT, creates a sharp image of the deep target at the correct depth. The effect becomes more pronounced at greater depths or longer propagation paths. Figure 14. View largeDownload slide Comparing a conventional RTM against a temporally spectral accurate (TSA) RTM. The boxes in red display the location of the targets (velocity anomalies). Top: zoom on shallow target. Bottom: zoom on deep target. (a) Injecting the data without applying dispersion results in a blurred image of the deep target at incorrect depth. (b) Injecting the data after applying the FTDT, creates a sharp image of the deep target at the correct depth. The effect becomes more pronounced at greater depths or longer propagation paths. 5 DISCUSSION 5.1 Computing cost The computational cost of computing the response of the wave equation (1) on a model using for instance an FD or SEM solution scales proportionally to NxNyNzNt, where Nx, y, z represent the number of gridpoints in x-, y- and z-directions respectively, and Nt represents the number of time steps. The time-dispersion transforms, conversely, apply a discrete Fourier transform and IFFT on M seismograms, which scales with $$M(N_t^2 + N_t\log _2(N_t))$$. The computing cost of the transforms is therefore modest compared to the full simulation when the number of recorded seismograms M is limited to a subset of the full grid, that is, M ≪ NxNyNz. The methodology proposed in this paper enables the use of large steps in time to reduce Nt, in contrast to the conventional solution to achieve higher temporal accuracy by increasing Nt. The CFL condition defines an upper time-step limit that may be used. However, we have shown that time dispersion accumulates with time, and that it distorts the kinematics of the wavefield propagating in the model, particularly at higher frequencies. Therefore, even with small time-step sizes, the time dispersion may play a significant role for long simulations. The conventional CFL limit may furthermore be increased with the use of higher order time-stepping schemes (Soubaras & Zhang 2008; Chu et al.2009; Amundsen & Pedersen 2017), whose time dispersion can be similarly removed with the time-dispersion transforms as shown in Appendix B. As time dispersion accumulates with time, and thus is more severe for longer simulations, it is difficult to categorically state where and when its use is advantageous. The most straightforward suggestion is to always execute simulations at maximum stable Δt, using the FTDT on the source wavelet to ensure correct kinematics of waves propagating on the model. One may apply the ITDT on a synthetic seismogram to test if the correction significantly alters the phases for a given modeling purpose. If the correction is judged to be significant, one may then apply the ITDT to all traces. 5.2 Viscoelastic media The derived workflow provides a ‘black box’ method to eliminate time dispersion without requiring any changes to the seismic modeling engine. The method then returns equivalent results regardless of the chosen time-step length. However, the method does not directly generalize to viscoelastic media. Consider a model with constant attenuation, used for a simulation with large time steps compared to a simulation with small time steps. The wavefield in the coarse simulation is accelerated with respect to the simulation with the small time steps. However, both simulations contain identical amounts of attenuation at any given time. The simulation with large time dispersion, thus, attenuates waves at a lower rate along its propagation path than the simulation with the small time step. The workflow presented in this paper cannot correct for this reduced amplitude attenuation effect. We speculate that the solution to this problem may lie in the use of a frequency-varying attenuation model which applies a larger than nominal amount of attenuation to higher frequencies. In such a case, the workflow presented in this paper cannot be used as a ‘black box’ method any longer, and requires modification on the side of the modeling engine. The exact form of the required changes to the attenuation model are the subject of future research. Finally, the noted effects of time dispersion on kinematics (larger than expected wavelengths due to acceleration) and viscoelasticity (less than expected attenuation) are the result of second-order accurate time stepping, and not of the time-dispersion transforms. The time-dispersion transforms provide a solution to undo the distorted kinematics of the simulation. The attenuation model will be distorted regardless, particularly notable when simulating long time-series with large steps in time. 6 CONCLUSIONS Coarse modeling of the wave equation enables the generation of synthetic seismograms at low computing cost, but may introduce ‘numerical dispersion’ errors that affect particularly the high-frequency range of the signal bandwidth. In this paper, we showed that the error introduced by the temporal FD operator, which we refer to as ‘time disperison’, can be analysed separately from that introduced by the spatial operator. In our paper, we present a novel description of time dispersion which is equally valid for the second-order system as well as for a first-order system of partial differential equations solved through a leap-frog time-stepping scheme. The time-dispersion error is described as an acceleration of the wavefield, as a function of frequency and the time step only, independent of the propagation path. We derived two simple transforms that add and remove time dispersion from time-series, which are both required to remove the adverse effects of time dispersion from seismic wave modeling. The two transforms were demonstrated on 1-D, 2-D and 3-D acoustic and elastic examples, showing perfect removal of time dispersion in FD, PS and SEM modeling. A simple application to RTM showed how pre-processing the data removes the detrimental effects of time dispersion in the resulting image. The time-dispersion transforms have a modest computing cost compared to the full simulation, when only a limited number of seismograms is required. Furthermore, the time-dispersion transforms enable the use of as large time steps as possible while honouring the CFL stability condition. We thus showed that we can model the wave equation coarsely in time, without the adverse effects of time dispersion. A note is made that the time-dispersion transforms do not apply to viscoelastic media without making alterations to the attenuation model. Acknowledgements This work was supported by SNF grant 2-77220-15. We wish to thank Andreas Fichtner and Marlies Vasmel for insightful discussions and suggestions. We thank four anonymous reviewers, Tarje-Nissen Meyer, Alexander Droujine, Jozef Kristek, the Editor Jörg Renner and the Associate Editor Ludovic Métivier for their detailed proofreading and comments that helped to significantly improve the manuscript. Footnotes 1Choosing coefficients c−1 = c1 = 1/Δt2 and c0 = −2/Δt2, and c∀|j| > 1 = 0 for the second-order accurate FD operator:   \begin{equation*} \theta _{\rm 2nd}(\omega ) = \sqrt{\frac{2-2\cos (\omega \Delta t)}{\omega ^2\Delta t^2}}= \frac{\sin (\frac{1}{2}\omega \Delta t)}{\frac{1}{2}\omega \Delta t}.\qquad\qquad\qquad\qquad\qquad(C6) \end{equation*} REFERENCES Aminzadeh F., 1996. “3-D salt and overthrust models”, in Applications of 3-D Seismic Data to Exploration and Production  chap. 28, pp. 247– 256, eds Weimer P., David T.L., AAPG/SEG, Tulsa. Amundsen L., Pedersen Ø., 2017. Time step n-tupling for wave equations, Geophysics  82 T249– T254. https://doi.org/10.1190/geo2017-0377.1 Google Scholar CrossRef Search ADS   Anderson J.E., Brytik V., Ayeni G., 2015. Numerical temporal dispersion corrections for broadband temporal simulation, RTM and FWI, in 85th Annual International Meeting  SEG, Expanded Abstracts, pp. 1096– 1100. Baysal E., Kosloff D.D., Sherwood J.W., 1983. Reverse time migration, Geophysics  48 1514– 1524. https://doi.org/10.1190/1.1441434 Google Scholar CrossRef Search ADS   Blanch J.O., Robertsson J. O.A., 1997. A modified Lax-Wendroff correction for wave propagation in media described by Zener elements, Geophys. J. Int.  131 381– 386. https://doi.org/10.1111/j.1365-246X.1997.tb01229.x Google Scholar CrossRef Search ADS   Chu C., Stoffa P.L., Seif R., 2009. 3D elastic wave modeling using modified high-order time stepping schemes with improved stability conditions, in 79th Annual International Meeting  SEG, Expanded Abstracts, pp. 2662– 2666. Dablain M., 1986. The application of high-order differencing to the scalar wave equation, Geophysics  51 54– 66. https://doi.org/10.1190/1.1442040 Google Scholar CrossRef Search ADS   Dai N., Liu H., Wu W., 2014. Solutions to numerical dispersion error of time FD in RTM, in 84th Annual International Meeting  SEG, Expanded Abstracts, pp. 4027– 4031. Dumbser M., Käser M., De La Puente J., 2007. Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D, Geophys. J. Int.  171 665– 694. https://doi.org/10.1111/j.1365-246X.2007.03421.x Google Scholar CrossRef Search ADS   Etgen J.T., Brandsberg-Dahl S., 2009. The pseudo-analytical method: application of pseudo-laplacians to acoustic and acoustic anisotropic wave propagation, in 79th Annual International Meeting  SEG, Expanded Abstracts, pp. 2552– 2556. Gardner G., Gardner L., Gregory A., 1974. Formation velocity and density—the diagnostic basics for stratigraphic traps, Geophysics  39 770– 780. https://doi.org/10.1190/1.1440465 Google Scholar CrossRef Search ADS   Kosloff D.D., Baysal E., 1982. Forward modeling by a Fourier method, Geophysics  47 1402– 1412. https://doi.org/10.1190/1.1441288 Google Scholar CrossRef Search ADS   Li Y., Wong M., Clapp R., 2013. Equivalent accuracy at a fraction of the cost: overcoming temporal dispersion, Tech. Rep. 150, Stanford Exploration Project. Li Y.E., Wong M., Clapp R., 2016. Equivalent accuracy at a fraction of the cost: overcoming temporal dispersion, Geophysics  81 T189– T196. https://doi.org/10.1190/geo2015-0398.1 Google Scholar CrossRef Search ADS   Liu H., Dai N., Niu F., Wu W., 2014. An explicit time evolution method for acoustic wave propagation, Geophysics  79 T117– T124. https://doi.org/10.1190/geo2013-0073.1 Google Scholar CrossRef Search ADS   Liu Y., 2013. Globally optimal finite-difference schemes based on least squares, Geophysics  78 T113– T132. https://doi.org/10.1190/geo2012-0480.1 Google Scholar CrossRef Search ADS   Liu Y., 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling, Geophys. J. Int.  197 1033– 1047. https://doi.org/10.1093/gji/ggu032 Google Scholar CrossRef Search ADS   Mittet R., 2017. On the internal interfaces in finite-difference schemes, Geophysics  82 T159– T182. https://doi.org/10.1190/geo2016-0477.1 Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Gális M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures  Cambridge University Press. Google Scholar CrossRef Search ADS   Qin Y., Quiring S., Nauta M., 2017. Temporal dispersion correction and prediction by using spectral mapping, in 79th Annual International Conference & Exhibition 2017  EAGE, Extended Abstracts, doi:10.3997/2214-4609.201700677. Song X., Fomel S., Ying L., 2013. Lowrank finite-differences and lowrank fourier finite-differences for seismic wave extrapolation in the acoustic approximation, Geophys. J. Int.  193 960– 969. https://doi.org/10.1093/gji/ggt017 Google Scholar CrossRef Search ADS   Soubaras R., Zhang Y., 2008. Two-step explicit marching method for reverse time migration, in 78th Annual International Meeting  SEG, Expanded Abstracts, pp. 2272– 2276. Stork C., 2013. Eliminating nearly all dispersion error from FD modeling and RTM with minimal cost increase, in 75th Annual International Conference & Exhibition  EAGE, Extended Abstracts, doi:10.3997/2214-4609.20130478. Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics  49 1259– 1266. https://doi.org/10.1190/1.1441754 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics  51 1954– 1966. https://doi.org/10.1190/1.1442051 Google Scholar CrossRef Search ADS   Tromp J., Komatitsch D., Liu Q., 2008. Spectral-element and adjoint methods in seismology, Commun. Comput. Phys.  3 1– 32. Virieux J., 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method, Geophysics  51 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Wang M., Xu S., 2015a. Finite-difference time dispersion transforms for wave propagation, Geophysics  80 WD19– WD25. https://doi.org/10.1190/geo2015-0059.1 Google Scholar CrossRef Search ADS   Wang M., Xu S., 2015b. Time dispersion prediction and correction for wave propagation, in 85th Annual International Meeting  SEG, Expanded Abstracts, pp. 3677– 3681. Wang Y., Liang W., Nashed Z., Li X., Liang G., Yang C., 2014. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method, Geophysics  79 T277– T285. https://doi.org/10.1190/geo2014-0078.1 Google Scholar CrossRef Search ADS   Xu Z., Jiao K., Cheng X., Sun D., King R., Nichols D., Vigh D., 2017. Time-dispersion filter for finite-difference modeling and reverse time migration, in 87th Annual International Meeting  SEG, Expanded Abstracts, pp. 4448– 4452. Yao G., Wu D., Debens H.A., 2016. Adaptive finite difference for seismic wavefield modelling in acoustic media, Sci. Rep.  6, doi:10.1038/srep30302. SUPPORTING INFORMATION Supplementary data are available at GJI online. Time_dispersion_1D_examples.m ITDT.m FTDT.m FD_elastic_anisotropic.tar.gz Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: FIRST-ORDER SYSTEM OF EQUATIONS We will first show how the factor θ can be recovered from the finite-difference operator itself. Note that in general, an nth-order derivative is written as $$\mathcal {F}\left[{\partial ^n {u}(t)}/{\partial t^n}\right] = (-i\omega )^n{\tilde{u}}(\omega )$$. Phase shifts are denoted as $$\mathcal {F}[u(t+\Delta t)]=e^{-i\omega \Delta t}\tilde{u}(\omega )$$. The temporal FD operator in the Fourier domain is thus:   \begin{eqnarray} \mathcal {F}\left[ \frac{{u}(t\!-\!\Delta t)\!-\!2{u}(t)\!+\!{u}(t\!+\!\Delta t)}{\Delta t^2} \right]\!\!=\!\left(\!\frac{e^{-i\omega \Delta t}\!-2\!+\!e^{i\omega \Delta t}}{\Delta t^2}\!\right)\tilde{u}(\omega ),\nonumber\\ \end{eqnarray} (A1)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\qquad= -\left[\frac{2-2\cos (\omega \Delta t)}{\Delta t^2}\right]\tilde{u}(\omega ),\nonumber\\ \end{eqnarray} (A2)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\quad= -\left[\frac{\sin (\frac{\omega \Delta t}{2})}{\frac{1}{2}\omega \Delta t} \right]^2\omega ^2 \tilde{u}(\omega ), \end{eqnarray} (A3)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\,\,\,\,\,\,= -\theta ^2\omega ^2 \tilde{u}(\omega ). \end{eqnarray} (A4)We note that −ω2 represents the second-order time derivative, and the additional factor θ2 corresponds to the undesired remaining Taylor expansion terms. We use this method to show the equivalence in the first-order system of equations solved in a leap-frog manner (Virieux 1986):   \begin{eqnarray} \mathcal {F}\left[\frac{u\left(t+\frac{\Delta t}{2}\right)-u\left(t-\frac{\Delta t}{2}\right)}{\Delta t} \right] &=& \left(\frac{e^{-i\omega \Delta t/2}-e^{i\omega \Delta t/2}}{\Delta t}\right)\tilde{u}(\omega ), \end{eqnarray} (A5)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=& -\left( \frac{2i\sin \left(\frac{\omega \Delta t}{2}\right)}{\Delta t} \right)\tilde{u}(\omega ), \end{eqnarray} (A6)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=& -\left( \frac{\sin \left(\frac{\omega \Delta t}{2}\right)}{\frac{1}{2}\omega \Delta t} \right)i\omega \tilde{u}(\omega ), \end{eqnarray} (A7)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=&-\theta i\omega \tilde{u}(\omega ). \end{eqnarray} (A8)We note that −iω represents the first-order time derivative, and the additional factor θ corresponds to the undesired remaining Taylor expansion terms. The factor θ is identical in both derivations, representing the same acceleration of the wavefield in every calculation. In other words, the time-dispersion error is identical when solving either the first- or second-order system of equations. The results derived in this paper apply thus equally to both cases. APPENDIX B: FOURTH-ORDER ACCURATE TIME DISPERSION B1 Taylor-series fourth-order accurate scheme The method developed in this paper may be extended to any order time stepping. We demonstrate this for a fourth-order accurate scheme in time, where we solve a system based on eq. (2), the Taylor series of the FD operator:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},t\!-\!\Delta t)\!-\!2\mathbf {u}(\mathbf {x},t)\!+\!\mathbf {u}(\mathbf {x},t\!+\!\Delta t)}{\Delta t^2} = \sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}}, \end{eqnarray} (B1)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\approx \left(L+ \frac{\Delta t^2L^2}{\frac{1}{2}4!}\right)\mathbf {u}(\mathbf {x},t).\nonumber\\ \end{eqnarray} (B2)In other words, the Taylor series of the temporal FD operator is truncated, and we substitute spatial operators using the wave equation to eliminate higher order temporal error terms. For this time-stepping scheme, we get the following relation between the time-dispersed (ω) and original frequencies (ω0):   $$\omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t\sqrt{1 - \frac{(\omega _0\Delta t)^2}{\frac{1}{2}4!}}}{2} \right) .$$ (B3)The square-root term corresponds to θ from eq. (3), truncated to the two terms solved in the fourth-order accurate scheme. Kinematically speaking, the second-order FD operator, that is, left-hand side of eq. (B2), accelerates the wavefield as discussed in the paper, while the fourth-order accurate scheme, that is, right-hand side of eq. (B2), decelerates the original wavefield, leading to a reduced overall time dispersion. Solving the system for ω0 gives:   \begin{eqnarray} \omega _0 = \frac{1}{\Delta t}\sqrt{6-\sqrt{12+24\cos (\omega \Delta t)}},\\ \ \ \ \ \ {\rm for }\quad\omega \Delta t\le \frac{2\pi }{3}\approx 2.1 \nonumber . \end{eqnarray} (B4)One can substitute identities (B3) and (B4) in the ITDT and FTDT, replacing eqs (13a) and (13b), to obtain an ITDT and FTDT for the fourth-order scheme. The ITDT based on eq. (B3) is equal to that found in Wang & Xu (2015a,b), who do not report the limitation ωΔt ≤ 2.1, which may cause unintended leaking of low frequencies into the higher signal spectrum. B2 Optimized-series fourth-order accurate scheme The advantage of a fourth-order accurate scheme in time is reduced time dispersion and a larger CFL condition. The Taylor series in eq. (B1) may be replaced with a more efficient series expansion that achieves an even larger CFL condition, at the cost of accumulating more time dispersion. Soubaras & Zhang (2008) and Amundsen & Pedersen (2017) describe such a procedure that achieves twice the CFL condition using the following series expansion:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},t\!-\!\Delta t)\!-\!2\mathbf {u}(\mathbf {x},t)\!+\!\mathbf {u}(\mathbf {x},t\!+\!\Delta t)}{\Delta t^2} \approx \left(L+ \frac{\Delta t^2L^2}{16}\right)\mathbf {u}(\mathbf {x},t).\nonumber\\ \end{eqnarray} (B5)Comparing eqs (B5) with (B2), we observe that the single difference between the two schemes is the factor 16. For this scheme, we may similarly propose the FTDT and ITDT, using identities:   \begin{eqnarray} \omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t\sqrt{1 - \frac{(\omega _0\Delta t)^2}{16}}}{2} \right) , \end{eqnarray} (B6a)  \begin{eqnarray} \omega _0= \frac{4}{\Delta t}\sin \left( \frac{\omega \Delta t}{4} \right). \end{eqnarray} (B6b) Note that eq. (B6a) is valid for $$\omega _0\Delta t \le 4\sin (\pi /4)=2\sqrt{2}\approx 2.8$$, related to the Nyquist sampling explained in Section 2.4. We show a simple example, using the pseudospectral method. The standard CFL limit is 2/π, whereas this fourth-order scheme allows stable computations using a CFL limit of 4/π. Fig. A1 shows a recording of a 25 Hz central frequency Ricker wavelet, using the same model as in the 1-D examples in this paper: Δx = 5 m, c = 2000 m s−1 and the distance between source and receiver is 1 km, but we use time steps at double the CFL limit of 1.2. As the PS method is spectrally accurate in space, only time dispersion accumulates. Using the proposed time-dispersion workflow proposed in this paper, we recover the Ricker wavelet without dispersion. Figure A1. View largeDownload slide Recording due to a pseudospectral, modified fourth-order accurate in time scheme, eq. (B5) (blue), and the derived FTDT and ITDT identities of eq. (B6) restore the original input Ricker wavelet (red). Figure A1. View largeDownload slide Recording due to a pseudospectral, modified fourth-order accurate in time scheme, eq. (B5) (blue), and the derived FTDT and ITDT identities of eq. (B6) restore the original input Ricker wavelet (red). APPENDIX C: CORRESPONDENCE BETWEEN SOLUTIONS SATISFYING THE FINITE-DIFFERENCE SYSTEM AND THE WAVE EQUATION In this section, we will show the relation between continuous solutions satisfying the wave equation and those satisfying the finite-difference (FD) system. Let u be a solution to the wave equation:   $$\mathcal {A}u(x,t) \stackrel{{\mbox{def}}}{=}\left( \frac{\partial ^2}{\partial t^2} - \mathcal {L} \right)u(x,t) = f(x,t),$$ (C1)where t > 0 denotes time, $$x\in \mathbb {R}^d$$, f describes the source term and $$\mathcal {L}$$ is a spatial operator independent of time. Now consider a related solution v to the FD system:   $$\mathsf {A}v(x,t) \stackrel{{\mbox{def}}}{=}(\mathsf {D}-\mathcal {L})v(x,t) = g(x,t),$$ (C2)where $$\mathsf {D}v(x,t) = \sum _n c_n v(x,t-n\Delta t)$$, with coefficients cn chosen such that $$\mathsf {D}$$ approximates the second-order time derivative. We will show that we can write v(x, t) as a function of u(x, t), to give a relation between solutions satisfying the FD system and those satisfying the wave equation. Applying a Fourier transform to eq. (C2) yields:   \begin{eqnarray} \int _{-\infty }^\infty \mathsf {A}v(x,t)e^{-i\omega t}\textrm{d}t = \left(\sum _{n}c_n e^{i \omega \Delta t n} - \mathcal {L}\right)\widehat{v}(x,\omega ), \end{eqnarray} (C3)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \left( -\omega ^2\theta ^2 - \mathcal {L}\right)\widehat{v}(x,\omega ), \end{eqnarray} (C4)using:   $$\theta (\omega ) = \sqrt{\frac{\sum _n c_n e^{i\omega \Delta t n}}{-\omega ^2}},$$ (C5)where it is assumed that Δt is small enough such that the quantity in the square root is positive.1 Now try as ansatz the following relation between the solution v of the FD system and u of the actual wave equation:   \begin{eqnarray} v(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{u}(x,\omega \theta ) e^{i\omega t}\textrm{d}\omega , \end{eqnarray} (C7)with Fourier transform:   \begin{eqnarray} \widehat{v}(x,\omega ) = \widehat{u}(x,{\omega }{\theta }). \end{eqnarray} (C8) Continuing eq. (C4):   \begin{eqnarray} \int _{-\infty }^\infty \mathsf {A}v(x,t) e^{-i\omega t}\textrm{d}t = ( -\omega ^2\theta ^2 - \mathcal {L})\widehat{v}(x,\omega ), \end{eqnarray} (C9)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \left( -\omega ^2\theta ^2 - \mathcal {L} \right)\widehat{u}(x,\omega \theta ), \end{eqnarray} (C10)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \widehat{f}(x,\omega \theta ). \end{eqnarray} (C11)The last step is due to the definition that eq. (C10) is a solution to the wave equation (C1). Thus, if f and $$\mathcal {L}$$ are regular enough such that the Fourier transform u is well defined, we can use definitions:   $$v(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{u}(x,\omega \theta ) e^{{i\omega t}}{\rm d}\omega ,$$ (C12)and   $$g(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{f}(x,\omega \theta )e^{{i\omega t}}\textrm{d}\omega ,$$ (C13)that satisfy the FD system for all t:   $$\mathsf {A}v(x,t) = g(x,t).$$ (C14)The solution in the FD system can thus be written as phase-shifted solutions to the true wave equation, independent of the spatial location x or the exact form of operator $$\mathcal {L}$$. Numerically solving the FD system may cause aliasing, in which case this continuous relation does not hold for the discrete situation. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Eliminating time dispersion from seismic wave modeling

, Volume 213 (1) – Apr 1, 2018
12 pages

/lp/ou_press/eliminating-time-dispersion-from-seismic-wave-modeling-R4qAJWqeGb
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx563
Publisher site
See Article on Publisher Site

### Abstract

Summary We derive an expression for the error introduced by the second-order accurate temporal finite-difference (FD) operator, as present in the FD, pseudospectral and spectral element methods for seismic wave modeling applied to time-invariant media. The ‘time-dispersion’ error speeds up the signal as a function of frequency and time step only. Time dispersion is thus independent of the propagation path, medium or spatial modeling error. We derive two transforms to either add or remove time dispersion from synthetic seismograms after a simulation. The transforms are compared to previous related work and demonstrated on wave modeling in acoustic as well as elastic media. In addition, an application to imaging is shown. The transforms enable accurate computation of synthetic seismograms at reduced cost, benefitting modeling applications in both exploration and global seismology. Numerical modelling, Computational seismology, Wave propagation 1 INTRODUCTION The understanding of the earth's subsurface from seismology observations is increasingly held back more by computing power rather than available data. In particular, propagating seismic data through synthetic models dominates the cost of wave equation-based imaging methods such as reverse time migration (RTM, Baysal et al.1983) and full-waveform inversion (Tarantola 1984). Coarsely discretizing the wavefield typically curbs the modeling costs, but introduces ‘numerical dispersion’ errors. The numerical dispersion distorts the modeled wavefield, affecting particularly the high-frequency content of the data. In practice, thus, only the low-frequency end of the full data bandwidth can be imaged accurately at reasonable cost. Breaking the barrier between modeling costs and accuracy is therefore important for many applications. In particular, modeling accurate wavefields with coarse discretization would enable accurate imaging at low cost. Accurate wavefield modeling implies that no errors arise from the space and time discretization; a topic that has been addressed in various ways by researchers in the past. For instance, the spatial and temporal error may be combined such that they partially cancel each other (see e.g. Etgen & Brandsberg-Dahl 2009; Wang et al.2014). The drawback of such approaches is the need for medium- and simulation-dependent stencils (Stork 2013). Instead, it is also possible to minimize the spatial and temporal error separately, as discussed below. Modeling seismic waves involve the computation of spatial gradients of the wavefield and of the medium. In exploration seismology, these spatial gradients are often computed with finite-difference (FD) approximations. The FD method (Moczo et al.2014) replaces derivatives in the wave equation with FD approximations of limited accuracy. A fine spatial discretization improves the accuracy of the FD method at significant computing cost. For example, reducing the spacing in 3-D scales to the power four with computing cost, for equal signal bandwidth. Instead, it is common to use higher order FD operators to reduce errors, even at a coarse discretization (see e.g. Song et al.2013; Liu 2013, 2014; Yao et al.2016). In the limit of using the entire wavefield to compute derivatives, the pseudospectral (PS) method achieves very high accuracy (Kosloff & Baysal 1982). Another high-order accurate spatial operator, derived from the weak formulation of the wave equation, is the spectral element method (SEM; Tromp et al.2008). SEM is used in global seismology to model the phases of seismograms with high accuracy, but the temporal error remains, and its removal would greatly benefit seismologists. Within FD, PS and SEM modeling, the time integration is carried out with explicit time stepping, an FD approximation. Accurate time stepping is then achieved using short steps in time, often at large computing cost. Alternatively, higher order Lax–Wendroff corrections can be used to reduce the errors (e.g. Dablain 1986; Blanch & Robertsson 1997; Dumbser et al.2007). Such correction schemes allow long time steps, although the Courant–Friedrichs–Lewy (CFL) condition defines an upper time-step limit for stable simulations. If short time steps are required for a stable simulation, large temporal errors are naturally prevented. However, temporal errors accumulate and may still become significant for a given modeling purpose, as the errors increase with number of propagated wavelengths and the associated frequency, which we will further address in the discussion section (Stork 2013). Additionally, optimized Lax–Wendroff schemes increase the stable time-step upper limit two- or threefold (Soubaras & Zhang 2008; Chu et al.2009; Liu et al.2014; Amundsen & Pedersen 2017), putting the temporal error at the forefront again. Stork (2013) took a radically different approach by suggesting that the temporal error can be filtered from seismograms at modest additional cost. A 1-D homogeneous example to illustrate the concept put forward by Stork is shown in Fig. 1. The procedure consists of three steps. First, a train of wavelets is modeled in an arbitrary model (e.g. a homogeneous model) to obtain a trace which gains temporal errors with increasing time (Fig. 1b). Second, using a simulation with very small steps in time, or with the use of the source wavelet itself (both options contain negligible or no temporal errors), a deconvolution is applied to the dispersed wavelets to create a ‘bank’ of filters that can remove the temporal error at various times (Fig. 1c). Finally, the desired wavefield simulation is carried out in the actual seismic domain of interest (which may be of arbitrarily high complexity and spatial dimension), and the previously computed filters are applied to remove the temporal error. Details of this approach can be found in Stork (2013), Li et al. (2013, 2016), Dai et al. (2014) and Anderson et al. (2015). A correct implementation of the filter is tedious, requiring many dispersed wavelets corresponding to the particular time-step characteristics of that simulation, and interpolation to correct wavelets at intermediate times. A simpler alternative using analytical filters was proposed by Wang & Xu (2015a,b), and later rederived in Qin et al. (2017), Mittet (2017) and Xu et al. (2017). The suggested analytical filters apply an altered discrete Fourier transform to seismograms which eliminates the time-stepping error in one operation, thus do not require the creation a filter bank ahead of the simulation or interpolation of the filter operation. However, adding and correspondingly removing the error with the proposed filters does not reproduce the input signal fully (see Section 2.5 below), suggesting that the filters are incomplete. Figure 1. View largeDownload slide The concept of a time-dispersion filter, cf. Stork (2013), illustrated in a 1-D homogeneous medium. (a) A simulation is created with a given source wavelet and time step. (b) In the absence of spatial effects distorting the wavelet, for example by carrying out a pseudospectral simulation in a 1-D homogeneous medium, a trace is gained with dispersion due to the time-stepping error only. (c) A filter bank is created which undoes the time–frequency-varying phase shift due to the temporal error. The created filter bank can then also be applied to traces modeled in complex media of higher dimensions. Figure 1. View largeDownload slide The concept of a time-dispersion filter, cf. Stork (2013), illustrated in a 1-D homogeneous medium. (a) A simulation is created with a given source wavelet and time step. (b) In the absence of spatial effects distorting the wavelet, for example by carrying out a pseudospectral simulation in a 1-D homogeneous medium, a trace is gained with dispersion due to the time-stepping error only. (c) A filter bank is created which undoes the time–frequency-varying phase shift due to the temporal error. The created filter bank can then also be applied to traces modeled in complex media of higher dimensions. To this date, the time-stepping error has only been derived for constant-velocity acoustic media, and it has only been conjectured that the error is independent of the spatial medium (e.g. Stork 2013). In this paper, we solve these outstanding questions by deriving an expression for the time-stepping error that is independent of the spatial medium and we refer to it as ‘time dispersion’. We show that time dispersion must be added to the injection wavelet, to model the correct kinematics and to correctly remove time dispersion afterwards. We derive a new exact filter pair to add and remove time dispersion from traces, which is fully invertible. We demonstrate the proposed method on acoustic and anisotropic elastic media, using FD, PS and SEM modeling. The results are valid for solving the second-order wave equation, or the system of first-order partial differential equations (Virieux 1986). Finally, the method is demonstrated on a simple RTM example, to show that the results apply to generic seismic modeling problems. 2 THEORY In a general form, the wave equation can be written as:   $$\frac{\partial ^2 \mathbf {u}(\mathbf {x},t)}{\partial t^2} = L(\mathbf {x})\mathbf {u}(\mathbf {x},t) + \mathbf {f}(\mathbf {x},t),$$ (1)where u(x, t) represents the particle displacement at position x and time t, L represents the elastic potential term and f the source function. We assume that the medium is time invariant, but may be heterogenous and anisotropic. We will discuss the implications of time dispersion for viscoelasticity in Section 5.2. 2.1 Finite-difference operator Numerical wave propagation methods such as FD, PS and SEM, typically compute the second-order time derivative in the wave equation using an FD approximation. In particular, a second-order accurate FD operator derived from a Taylor series is often used:   $$\frac{\mathbf {u}({t\!-\!\Delta t})-2\mathbf {u}({t}) + \mathbf {u}({t\!+\!\Delta t})}{\Delta t^2} = \sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}},$$ (2)where Δt denotes the time step. Disregarding O(Δt2) and higher order terms from eq. (2), thus assuming that the operator represents just the second-order time derivative, introduces a numerical error. Accurate simulations are typically achieved using short steps in time or by using higher order time-stepping schemes (Dablain 1986). Note that, although eq. (2) corresponds to an approximation of a second-order time derivative, the results described in the following sections apply equally to a system of first-order partial differential equations solved with a leap-frog time-stepping scheme (Virieux 1986) as shown in Appendix A. The theory in this paper may also be used to derive similar correction filters if higher order time-stepping schemes are used, for example to allow longer steps in time. Appendix B demonstrates this procedure for a higher order time-stepping scheme (Soubaras & Zhang 2008; Amundsen & Pedersen 2017). 2.2 Describing time dispersion Denoting the Fourier transform operator by $$\mathcal {F}[\cdot ]$$ and its inverse by $$\mathcal {F}^{-1}[\cdot ]$$, we begin by transforming the right-hand side of eq. (2) to the frequency domain:   \begin{eqnarray} \mathcal {F}\left[\sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}}\right]= -\theta ^2\omega ^2 {\tilde{\mathbf {u}}}(\omega ). \end{eqnarray} (3)The right-hand side of eq. (3) corresponds to the second-order time derivative multiplied with an additional factor θ2 corresponding to the undesired remaining Taylor expansion terms:   $$\theta ^2\!=\!\sum _{k=1}^\infty \!\frac{(-i\omega \!\Delta t)^{2(k-1)}}{\frac{1}{2}(2k)!}\!=\!1-\frac{\Delta t^2\omega ^2}{\frac{1}{2}4!}+\ldots \!=\!\left(\!\frac{\sin (\!\frac{\omega \Delta t}{2}\!)}{\frac{1}{2}\omega \Delta t}\!\right)^2\!.$$ (4) We insert the second-order accurate FD operator in time in eq. (1), and find that:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},{t\!-\!\Delta t})\!-\!2\mathbf {u}(\mathbf {x},{t})\!+\!\mathbf {u}(\mathbf {x},{t\!+\!\Delta t})}{\Delta t^2}\!\approx \!L(\mathbf {x}){\mathbf {u}(\mathbf {x},t)}\!+\mathbf {f}( \mathbf {x},\!t), \end{eqnarray} (5)thus represents a modified wave equation, solving exactly:   $$\frac{\partial ^2 \mathbf {\hat{u}}(\mathbf {x},t)}{\partial t^2} = \mathcal {F}^{-1}\Bigg [\frac{1}{\theta ^2}\mathcal {F}\left [L(\mathbf {x}){\mathbf {\hat{u}}(\mathbf {x},t)}+\mathbf {f}(\mathbf {x},t)\right ]\Bigg ].$$ (6)The formulation in eq. (6) follows from applying a forward and inverse Fourier transform on eq. (5), and eliminating θ from the left-hand side. We explicitly note that $$\mathbf {\hat{u}}(\mathbf {x},t)$$ is the time-dispersed solution to the FD system, without specifying anything about the spatial operator L or its accuracy. Solutions that satisfy the FD system are closely related to solutions to the wave equation, as is detailed in the section below. In Appendix C, we present a proof that the relation holds in the continuous case. We will show that eq. (6) describes an acceleration of the original wavefield. Consider first an illustrative example where we substitute t → 2t in the general wave equation (1):   \begin{eqnarray} \frac{\partial ^2 \mathbf {u}(\mathbf {x},2t)}{\partial (2t)^2} = L(\mathbf {x})\mathbf {u}(\mathbf {x},2t) + \mathbf {f}(\mathbf {x},2t), \end{eqnarray} (7)  \begin{eqnarray} \frac{\partial ^2 \mathbf {u}(\mathbf {x},2t)}{\partial t^2} = 2^2\left (L(\mathbf {x})\mathbf {u}(\mathbf {x},2t) + \mathbf {f}(\mathbf {x},2t)\right ). \end{eqnarray} (8)Eq. (8) computes the wavefield u(x, 2t) using the original time differential ∂t, by multiplying the spatial terms with factor 22. By analogy, eq. (5) similarly represents an accelerated wavefield if we denote $$\mathbf {u}(\mathbf {x},\hat{t})$$ as the accelerated solution and modify the source function $$\mathbf {f}(\mathbf {x},t)\rightarrow \mathbf {f}(\mathbf {x},\hat{t})$$:   $$\frac{\partial ^2 \mathbf {u}(\mathbf {x},\hat{t})}{\partial t^2} = \mathcal {F}^{-1}\Bigg [\frac{1}{\theta ^2}\mathcal {F}\left [L(\mathbf {x}){\mathbf {u}(\mathbf {x},\hat{t})}+\mathbf {f}(\mathbf {x},\hat{t})\right ]\Bigg ].$$ (9)Eq. (9) implies that the second-order temporal FD operator solves an accelerated version of the wave equation perfectly after proper pre-processing of the source wavelet. Would we use the original injection wavelet, the accelerated wavefield would correspondingly contain erroneous wavelengths, changing the kinematic and dynamic properties of the simulation. To our knowledge, this important observation is absent from previously published literature on time dispersion. As an explicit example, consider the heterogeneous acoustic wave equation:   $$\frac{\partial ^2 p(\mathbf {x},t)}{\partial t^2} = c^2(\mathbf {x})\rho (\mathbf {x})\nabla \cdot \rho ^{-1}(\mathbf {x})\nabla p(\mathbf {x},t),$$ (10)where ∇ represents the spatial derivative operator. Using second-order differencing in time solves the modified equation:   $$\frac{\partial ^2 p(\mathbf {x},\hat{t})}{\partial t^2}=\mathcal {F}^{-1}\Bigg [\frac{c(\mathbf {x})^2}{\theta ^2}\mathcal {F}\left [\rho (\mathbf {x})\nabla \cdot \rho ^{-1}(\mathbf {x})\nabla p(\mathbf {x},\hat{t})\right ]\Bigg ],$$ (11)The factor c(x)/θ takes the role of the effective velocity, that is, higher frequencies propagate faster than lower frequencies. This effect, caused by the second-order finite differencing in time, is the time dispersion shown in Fig. 2. Note that Dablain (1986) derived a similar expression for the time dispersion, assuming a homogenous acoustic medium and using plane-wave theory. However, this severely restrictive limitation is not required. The theory derived in this section remains valid for mode conversions, anisotropy or complex medium responses: signals are simply accelerated as a function of their frequency. Figure 2. View largeDownload slide The function 1/θ originating from the second-order FD operator represents a frequency-dependent velocity increase. The horizontal axis represents the range from 0 frequency up to f = 1/2Δt. Figure 2. View largeDownload slide The function 1/θ originating from the second-order FD operator represents a frequency-dependent velocity increase. The horizontal axis represents the range from 0 frequency up to f = 1/2Δt. 2.3 Time-dispersion transforms By recognizing time dispersion as an acceleration of the wavefield, we can consider it as a problem that only affects the time axis. Accelerating the source $$\mathbf {f}(\mathbf {x},t)\rightarrow \mathbf {f}(\mathbf {x},\hat{t})$$, thus adding time dispersion, or decelerating a recorded trace $$\mathbf {u}(\mathbf {x},\hat{t})\rightarrow \mathbf {u}(\mathbf {x},t)$$, thus removing time dispersion, can therefore be carried out using a filter operating on a single time-series. Adding dispersion to a signal of frequencies ω0 gives rise to the following spectrum, due to the acceleration factor 1/θ:   $$\omega = \frac{\omega _0}{\theta }.$$ (12)Note that this relation between the dispersed and original frequencies is independent of the medium. The factor θ only depends on the time step and frequency of the wavelet, describing time dispersion with no knowledge of the propagation path, medium, or even spatial accuracy. Substituting eq. (4) in eq. (12) provides a relation between the original and time-dispersed frequencies:   \begin{eqnarray} \omega _0 = \frac{2}{\Delta t}\sin \left(\frac{\omega \Delta t}{2} \right), \end{eqnarray} (13a)  \begin{eqnarray} \omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t}{2} \right) . \end{eqnarray} (13b) The spectral content of the recorded signal is thus identical to the correct signal, it is merely recorded at the wrong frequency. However, we can map the data back to the desired frequency. We propose a map between the correct and dispersed phase rotations of a recording by inserting eq. (13a) in a forward Fourier transform. We thus add dispersion to a trace u(t) as follows:   $$\hat{u}(t)= \mathcal {F}^{-1}\left[ \int _{-\infty }^\infty u(t) e^{-i2\sin \left(\frac{\omega \Delta t}{2}\right)\frac{t}{\Delta t}} \textrm{d}t\right].$$ (14)We will refer to this operation as the Forward Time Dispersion Transform (FTDT). Eq. (13b) allows us to map a temporally dispersed signal back to the correct phase rotations, in a similar fashion as the FTDT. We thus remove time dispersion from trace $$\hat{u}(t)$$ as follows:   \begin{eqnarray} {u}(t) = \mathcal {F}^{-1}\left[ \int _{-\infty }^\infty \hat{u}(t) e^{-i2\sin ^{-1}\left(\frac{\omega _0\Delta t}{2}\right)\frac{t}{\Delta t}} \textrm{d}t\right]. \end{eqnarray} (15)We will refer to this operation as the Inverse Time Dispersion Transform (ITDT). 2.4 Sampling a dispersed wavefield As we record the accelerated wavefield on the seismograms, we must carefully consider the sampling requirements. The Nyquist frequency fN = 1/(2Δt) represents the maximum frequency that can be recorded without aliasing. Considering a maximum signal acceleration of π/2 ≈ 1.6 (Fig. 2), we must sample the simulated field with an increased sampling rate fN, d = π/2 fN = π/(4Δt) to capture signals at their original Nyquist frequency. The support of the sin −1 function in eq. (13b) coincides exactly with this frequency range. 2.5 Comparison to filters derived in Wang & Xu (2015a) While the ITDT, as denoted by eq. (15), is equal to that found in Wang & Xu (2015a), the FTDT as given in eq. (14) is not. Wang & Xu (2015a) use for their FTDT, the adjoint of the ITDT as given in eq. (15). The adjoint of this modified Fourier transform does not provide a complete inverse, which leads to an incorrect phase and amplitude response. We demonstrate this from a pure signal-processing perspective. Applying time dispersion and correspondingly removing it, should return the same data as we merely map one frequency into another and then reverse the process. In Fig. 3, the filter pair from Wang & Xu (2015a) has been applied on a Ricker wavelet. The central frequency is 8 Hz with time step Δt = 15 ms. We note that the filter pair from Wang & Xu (2015a) does not reproduce the input wavelet, while the filter pair presented in this paper does. The effect becomes more notable at higher frequencies. Figure 3. View largeDownload slide Adding and correspondingly removing time dispersion should return the input data. The filters from Wang & Xu (2015a) do not achieve this, whereas the filters in this paper do. Wang & Xu (2015a) use for the FTDT an adjoint version of the ITDT, which is an incorrect filter. Figure 3. View largeDownload slide Adding and correspondingly removing time dispersion should return the input data. The filters from Wang & Xu (2015a) do not achieve this, whereas the filters in this paper do. Wang & Xu (2015a) use for the FTDT an adjoint version of the ITDT, which is an incorrect filter. 2.6 Implementation of time-dispersion filters Eqs (14) and (15) are discretized for sampled data. We implement the FTDT using a discrete Fourier transform and a standard inverse fast Fourier transform (IFFT), for a discrete signal u[n] of N samples recorded uniformly with a sampling interval Δt. The FTDT becomes:   $$\hat{u}[n] = {\rm IFFT}\!\left[ \sum _{j=0}^{N\!-\!1} u[j]e^{\!-\!i2\!\sin \left(\!\frac{\pi f}{N}\!\right) j } \right],\quad \, f\!\in \!\left\lbrace 0,..., \frac{N\!\!-\!1}{2}\right\rbrace .$$ (16)Similarly for the ITDT, we obtain:   $${u}[n] = {\rm IFFT}\!\left[ \sum _{j=0}^{N\!-\!1}\!\hat{u}[j]e^{\!-\!i2\!\sin ^{\!-\!1}\left(\!\frac{\pi f}{N}\!\right) j } \right]\!,\quad \, f\!\in \!\left\lbrace 0,..., \frac{N\!\!-\!1}{\pi }\right\rbrace .$$ (17)The frequency vector f is limited to the support of the sin −1 function, as noted in Section 2.4. Note that discretizing the transforms in this formation makes the operations independent of Δt. We can thus apply the transforms without any knowledge of simulation parameters. We zero-pad the signal to twice the length to avoid adverse edge effects from the Fourier transforms. A MATLAB implementation of the proposed FTDT and ITDT is provided with the digital supplements to this paper. 2.7 Application of time-dispersion filters Eqs (16) and (17) add and remove time dispersion from signals. This result holds for any time-invariant medium L(x), spatial error O(Δxn), and is valid for both solving the second-order wave equation as well as the corresponding system of first-order partial differential equations (Virieux 1986). To remove time dispersion from seismic traces, we thus follow the workflow from Fig. 4: Apply the FTDT to the source function f, as given by eq. (9). The pre-processing of the wavelet ensures that the kinematics of the wavefield is retained. From a numerical perspective, it means that we inject a wavelet with an identical amount of time dispersion as that accumulating in the wavefield. Run the standard seismic wave simulation application and record traces wherever desired. Apply the ITDT to recorded traces after the simulation is completed. This recovers the original, time-dispersion-free wavefield. Figure 4. View largeDownload slide The proposed workflow to eliminate time dispersion from synthetic seismograms. No changes are required to the standard seismic modeling application. Furthermore, the filters are independent of the medium, and depend only on the time stepping and frequency of the signals. Figure 4. View largeDownload slide The proposed workflow to eliminate time dispersion from synthetic seismograms. No changes are required to the standard seismic modeling application. Furthermore, the filters are independent of the medium, and depend only on the time stepping and frequency of the signals. 3 EXAMPLES In this section, we apply the workflow of Fig. 4 on FD, PS and SEM modeling results for various dimensions and models. The results are compared to analytical 1-D solutions, or against simulations using short steps in time, to show that equivalent accuracy is reached at reduced cost. 3.1 1-D acoustic examples We demonstrate the FTDT and ITDT with two simple 1-D examples. The first is a PS computed solution to the acoustic wave equation, free of spatial errors. In addition, a staggered second-order spatial FD operator is used to generate a solution that illustrates that the transforms are independent of spatial errors. We use a spatial step of Δx = 5 m, a temporal step of Δt = 15 ms and a constant velocity c = 2000 m s−1, with source and receiver separated by 1 km. We inject a Ricker wavelet of 35 Hz central frequency, with FTDT applied, thus the recorded wavelet has propagated about 20 wavelengths. The code to perform these simulations is provided with the digital supplements to this paper. 3.1.1 Pseudospectral spatial operator The PS method is spectrally accurate in space, and thus only accumulates time dispersion. Fig. 5 compares the PS solution with an analytical solution of the expected wavelet arrival. We observe that applying the FTDT to the analytical solution reproduces the time-dispersed solution, with a maximum error on the order of 10−10. The error is likely a combination of minor artefacts of the discrete Fourier transform in both the PS spatial operator and in the temporal dispersion transforms. Conversely, we construct a dispersion-free recording by applying the ITDT to the PS solution. Figure 5. View largeDownload slide Time dispersion owing to the pseudospectral method (blue) compared to the analytical solution (green). The time-dispersion transform workflow (Fig. 4) recreates the analytical solution (black dots), while applying the FTDT to the analytical solution recreates the signal computed with the PS method (red dots). The maximum difference between the analytical solution and the computed, corrected recording is on the order of 10−10. Figure 5. View largeDownload slide Time dispersion owing to the pseudospectral method (blue) compared to the analytical solution (green). The time-dispersion transform workflow (Fig. 4) recreates the analytical solution (black dots), while applying the FTDT to the analytical solution recreates the signal computed with the PS method (red dots). The maximum difference between the analytical solution and the computed, corrected recording is on the order of 10−10. 3.1.2 Second-order accurate staggered spatial operator The staggered second-order spatial FD operator introduces spatial errors in addition to time dispersion. Fig. 6 compares two seismograms using a long and short time step, differing by a factor 10 to create a ‘course’ and a ‘fine’ simulation in time. It can be concluded that the FTDT and ITDT correctly model and remove time dispersion, independent of the spatial errors. We observe that the error between the two solutions is on the order of 10−3, which predominantly can be explained by remaining time dispersion in the short time-step solution. Figure 6. View largeDownload slide Time dispersion in the second-order accurate spatial FD method (blue), compared to a short time-step solution (green). The time-dispersion transform workflow (Fig. 4) recreates the solution from the short time steps (black dots), while applying the FTDT to the short time-step solution recreates the signal computed with long time steps. All remaining ‘ringing’ from the spatial errors is left untouched by the transforms. Figure 6. View largeDownload slide Time dispersion in the second-order accurate spatial FD method (blue), compared to a short time-step solution (green). The time-dispersion transform workflow (Fig. 4) recreates the solution from the short time steps (black dots), while applying the FTDT to the short time-step solution recreates the signal computed with long time steps. All remaining ‘ringing’ from the spatial errors is left untouched by the transforms. 3.2 2-D elastic examples 3.2.1 2-D anisotropic elastic example using fourth-order accurate FD The time-dispersion transforms are not limited or influenced by the number of spatial dimensions. We demonstrate their performance on a 2-D elastic anisotropic model, using a staggered grid fourth-order spatial FD operator, to solve the first-order system of equations of the wave equation in a leap-frog fashion (Virieux 1986). The model is shown in Fig. 7, with a vertical transverse isotropy given by Thomsen's parameters δ = ε = 0.1 (Thomsen 1986). The model is discretized with Δx = Δz = 1 m, and an explosion-type Ricker wavelet of 80 Hz central frequency is used. Fig. 8 shows snapshots of the wavefield displaying a large variety of different wave types. As shown on a trace recorded close to the free surface, we fully recreate the short time-step solution from the coarse simulation. The difference between the two solutions is small for early arrivals, but becomes more apparent for later times. The code to perform this simulation is provided with the digital supplements of this paper. Figure 7. View largeDownload slide Elastic model used for the simulation shown in Fig. 8. A vertically transverse isotropy is additionally used, with Thomsen's parameters ε = δ = 0.1. The source and receiver positions are shown by the red star and blue triangle, respectively. Figure 7. View largeDownload slide Elastic model used for the simulation shown in Fig. 8. A vertically transverse isotropy is additionally used, with Thomsen's parameters ε = δ = 0.1. The source and receiver positions are shown by the red star and blue triangle, respectively. Figure 8. View largeDownload slide Snapshots of a 2-D elastic anisotropic simulation with a source at the red star and receiver at the blue triangle. The horizontal particle velocity recorded at the receiver is shown on the right. The long (Δt = 0.32 ms) and short (Δt = 0.032 ms) time-step simulations are compared, as well as the solution after applying the time-dispersion transforms. The solution after applying the time-dispersion transforms matches the short time-step solution for all times. Figure 8. View largeDownload slide Snapshots of a 2-D elastic anisotropic simulation with a source at the red star and receiver at the blue triangle. The horizontal particle velocity recorded at the receiver is shown on the right. The long (Δt = 0.32 ms) and short (Δt = 0.032 ms) time-step simulations are compared, as well as the solution after applying the time-dispersion transforms. The solution after applying the time-dispersion transforms matches the short time-step solution for all times. 3.2.2 2-D isotropic elastic example using the spectral element method We demonstrate the time-dispersion transforms on the SEM, using SPECFEM2D with its default isotropic elastic model shown in Fig. 9 (Tromp et al.2008). We inject a Ricker wavelet of 25 Hz central frequency, with FTDT applied. We compare solutions using a long time step (Δt = 1 ms) and a short time step (Δt = 0.1 ms). Fig. 10 shows that applying the ITDT to the long time-step solution correctly reproduces the solution from the short time step, where the presence of time dispersion is negligible. Figure 9. View largeDownload slide The SPECFEM 2-D (Tromp et al.2008) default four layer model as shown with four different colours, with as overlay the used mesh, as well as the vertical particle displacement at 0.38 s. The receivers are shown in black along the surface and the blue rhombus represents the trace that is zoomed-in on in Fig. 10. Figure 9. View largeDownload slide The SPECFEM 2-D (Tromp et al.2008) default four layer model as shown with four different colours, with as overlay the used mesh, as well as the vertical particle displacement at 0.38 s. The receivers are shown in black along the surface and the blue rhombus represents the trace that is zoomed-in on in Fig. 10. Figure 10. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1 ms) and a short (Δt = 0.1 ms) time-step 2-D elastic spectral element simulation using SPECFEM2D (Tromp et al.2008). The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. Figure 10. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1 ms) and a short (Δt = 0.1 ms) time-step 2-D elastic spectral element simulation using SPECFEM2D (Tromp et al.2008). The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. 3.3 3-D elastic example using high-order FD Spatial numerical errors can be limited using a high-order spatial FD operator. Here, we demonstrate the time-dispersion transforms for such a case on the 3-D elastic SEG/EAGE overthrust model (Aminzadeh 1996), shown in Fig. 11. We use a shear wave velocity of 0.6 times the P-wave velocity and derive the density model through Gardner's relation (Gardner et al.1974). We use a 20 point spatial operator with optimized coefficients (Liu 2014) and a second-order accurate FD operator in time to solve the system of first-order equations in the conventional leap-frog fashion (Virieux 1986). We inject a Ricker wavelet of 25 Hz central frequency, with FTDT applied. Fig. 12 shows that the ITDT applied to a long time-step seismogram removes time dispersion, thereby reproducing the short time-step solution where the presence of time dispersion is negligible. Figure 11. View largeDownload slide The P-wave velocities for the 3-D SEG/EAGE overthrust model (Aminzadeh 1996), the black lines denote the location of the slices projected on the 3-D cube. The receivers are denoted with blue (overlapping the black line), the black rhombus represents the trace that is zoomed-in on Fig. 12, the source is located at 100 m depth below the black star. Figure 11. View largeDownload slide The P-wave velocities for the 3-D SEG/EAGE overthrust model (Aminzadeh 1996), the black lines denote the location of the slices projected on the 3-D cube. The receivers are denoted with blue (overlapping the black line), the black rhombus represents the trace that is zoomed-in on Fig. 12, the source is located at 100 m depth below the black star. Figure 12. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1.6 ms) and a short (Δt = 0.16 ms) time-step 3-D elastic FD simulation. The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths, extracted from the full 3-D gather. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. Figure 12. View largeDownload slide Comparison between the recorded horizontal particle velocity from a long (Δt = 1.6 ms) and a short (Δt = 0.16 ms) time-step 3-D elastic FD simulation. The top panel compares traces from the two simulations, in addition to the ITDT applied to the long time-step solution. The left bottom panel shows the difference between the two shot gathers of different time-step lengths, extracted from the full 3-D gather. The right bottom panel shows the difference after applying the ITDT to the long time-step solution. 4 APPLICATION TO REVERSE TIME MIGRATION RTM recovers subsurface reflectors by evaluating the zero-lag time correlation between a forward propagated source wavefield and a time-reversed backpropagated receiver wavefield. The wavefields accumulate time dispersion in opposing time directions, which blurs the migration. Stork (2013) notes that this blurring is mitigated by adding time dispersion to the receiver traces before time reversal – the time-reversal process removes the added time dispersion as the simulation runs backwards in time, and as a result the two wavefields contain equal time dispersion at all times. We will refer to such an RTM as a temporally spectral accurate (TSA) RTM. Effectively, we pre-process the data by adding time dispersion in accordance with the modified wave equation (9) that we compute numerically. The workflow is shown in Fig. 13. Figure 13. View largeDownload slide Workflow to eliminate time dispersion from reverse time migration. The ITDT and FTDT are each other's inverse, thus the dispersed synthetic traces may be used straight as input to the reverse run. Conversely, the recorded data are pre-processed to prepare it for the reverse-time simulation. Effectively, thus, only two processing steps are required. Figure 13. View largeDownload slide Workflow to eliminate time dispersion from reverse time migration. The ITDT and FTDT are each other's inverse, thus the dispersed synthetic traces may be used straight as input to the reverse run. Conversely, the recorded data are pre-processed to prepare it for the reverse-time simulation. Effectively, thus, only two processing steps are required. We show a simple example of the TSA RTM in an acoustic case with a second-order accurate in time and eight-order accurate in space-staggered FD simulation (Virieux 1986). Our model contains a shallow as well as a deep box-shaped target (velocity anomalies). The simplicity of the model makes the size of the errors more visible, as blurring of phases in complex media is less apparent when multiple reflectors blur. Fig. 14(a) corresponds to a case where we inject ‘recorded traces’ (numerical dispersion-free) and a forward modeled source wavefield (numerically dispersed) in a smooth velocity model, and reverse propagate the two simulations in tandem to produce an image by cross-correlation of snapshots. The image of the shallow target is well recovered (due to the short propagation paths of reflections in the shallow part of the model), but the deep target is imaged at the incorrect depth and with distorted phases. In Fig. 14(b), we have injected the time-reversed receiver data after applying the FTDT. The procedure images both the shallow and deep targets accurately. This result comes at no extra computing cost other than pre-processing the recorded traces with the altered forward Fourier transform and standard inverse (fast) Fourier transform constituting the FTDT. An important consequence is that relatively long time steps may be used in the RTM simulation without compromising on accuracy. Figure 14. View largeDownload slide Comparing a conventional RTM against a temporally spectral accurate (TSA) RTM. The boxes in red display the location of the targets (velocity anomalies). Top: zoom on shallow target. Bottom: zoom on deep target. (a) Injecting the data without applying dispersion results in a blurred image of the deep target at incorrect depth. (b) Injecting the data after applying the FTDT, creates a sharp image of the deep target at the correct depth. The effect becomes more pronounced at greater depths or longer propagation paths. Figure 14. View largeDownload slide Comparing a conventional RTM against a temporally spectral accurate (TSA) RTM. The boxes in red display the location of the targets (velocity anomalies). Top: zoom on shallow target. Bottom: zoom on deep target. (a) Injecting the data without applying dispersion results in a blurred image of the deep target at incorrect depth. (b) Injecting the data after applying the FTDT, creates a sharp image of the deep target at the correct depth. The effect becomes more pronounced at greater depths or longer propagation paths. 5 DISCUSSION 5.1 Computing cost The computational cost of computing the response of the wave equation (1) on a model using for instance an FD or SEM solution scales proportionally to NxNyNzNt, where Nx, y, z represent the number of gridpoints in x-, y- and z-directions respectively, and Nt represents the number of time steps. The time-dispersion transforms, conversely, apply a discrete Fourier transform and IFFT on M seismograms, which scales with $$M(N_t^2 + N_t\log _2(N_t))$$. The computing cost of the transforms is therefore modest compared to the full simulation when the number of recorded seismograms M is limited to a subset of the full grid, that is, M ≪ NxNyNz. The methodology proposed in this paper enables the use of large steps in time to reduce Nt, in contrast to the conventional solution to achieve higher temporal accuracy by increasing Nt. The CFL condition defines an upper time-step limit that may be used. However, we have shown that time dispersion accumulates with time, and that it distorts the kinematics of the wavefield propagating in the model, particularly at higher frequencies. Therefore, even with small time-step sizes, the time dispersion may play a significant role for long simulations. The conventional CFL limit may furthermore be increased with the use of higher order time-stepping schemes (Soubaras & Zhang 2008; Chu et al.2009; Amundsen & Pedersen 2017), whose time dispersion can be similarly removed with the time-dispersion transforms as shown in Appendix B. As time dispersion accumulates with time, and thus is more severe for longer simulations, it is difficult to categorically state where and when its use is advantageous. The most straightforward suggestion is to always execute simulations at maximum stable Δt, using the FTDT on the source wavelet to ensure correct kinematics of waves propagating on the model. One may apply the ITDT on a synthetic seismogram to test if the correction significantly alters the phases for a given modeling purpose. If the correction is judged to be significant, one may then apply the ITDT to all traces. 5.2 Viscoelastic media The derived workflow provides a ‘black box’ method to eliminate time dispersion without requiring any changes to the seismic modeling engine. The method then returns equivalent results regardless of the chosen time-step length. However, the method does not directly generalize to viscoelastic media. Consider a model with constant attenuation, used for a simulation with large time steps compared to a simulation with small time steps. The wavefield in the coarse simulation is accelerated with respect to the simulation with the small time steps. However, both simulations contain identical amounts of attenuation at any given time. The simulation with large time dispersion, thus, attenuates waves at a lower rate along its propagation path than the simulation with the small time step. The workflow presented in this paper cannot correct for this reduced amplitude attenuation effect. We speculate that the solution to this problem may lie in the use of a frequency-varying attenuation model which applies a larger than nominal amount of attenuation to higher frequencies. In such a case, the workflow presented in this paper cannot be used as a ‘black box’ method any longer, and requires modification on the side of the modeling engine. The exact form of the required changes to the attenuation model are the subject of future research. Finally, the noted effects of time dispersion on kinematics (larger than expected wavelengths due to acceleration) and viscoelasticity (less than expected attenuation) are the result of second-order accurate time stepping, and not of the time-dispersion transforms. The time-dispersion transforms provide a solution to undo the distorted kinematics of the simulation. The attenuation model will be distorted regardless, particularly notable when simulating long time-series with large steps in time. 6 CONCLUSIONS Coarse modeling of the wave equation enables the generation of synthetic seismograms at low computing cost, but may introduce ‘numerical dispersion’ errors that affect particularly the high-frequency range of the signal bandwidth. In this paper, we showed that the error introduced by the temporal FD operator, which we refer to as ‘time disperison’, can be analysed separately from that introduced by the spatial operator. In our paper, we present a novel description of time dispersion which is equally valid for the second-order system as well as for a first-order system of partial differential equations solved through a leap-frog time-stepping scheme. The time-dispersion error is described as an acceleration of the wavefield, as a function of frequency and the time step only, independent of the propagation path. We derived two simple transforms that add and remove time dispersion from time-series, which are both required to remove the adverse effects of time dispersion from seismic wave modeling. The two transforms were demonstrated on 1-D, 2-D and 3-D acoustic and elastic examples, showing perfect removal of time dispersion in FD, PS and SEM modeling. A simple application to RTM showed how pre-processing the data removes the detrimental effects of time dispersion in the resulting image. The time-dispersion transforms have a modest computing cost compared to the full simulation, when only a limited number of seismograms is required. Furthermore, the time-dispersion transforms enable the use of as large time steps as possible while honouring the CFL stability condition. We thus showed that we can model the wave equation coarsely in time, without the adverse effects of time dispersion. A note is made that the time-dispersion transforms do not apply to viscoelastic media without making alterations to the attenuation model. Acknowledgements This work was supported by SNF grant 2-77220-15. We wish to thank Andreas Fichtner and Marlies Vasmel for insightful discussions and suggestions. We thank four anonymous reviewers, Tarje-Nissen Meyer, Alexander Droujine, Jozef Kristek, the Editor Jörg Renner and the Associate Editor Ludovic Métivier for their detailed proofreading and comments that helped to significantly improve the manuscript. Footnotes 1Choosing coefficients c−1 = c1 = 1/Δt2 and c0 = −2/Δt2, and c∀|j| > 1 = 0 for the second-order accurate FD operator:   \begin{equation*} \theta _{\rm 2nd}(\omega ) = \sqrt{\frac{2-2\cos (\omega \Delta t)}{\omega ^2\Delta t^2}}= \frac{\sin (\frac{1}{2}\omega \Delta t)}{\frac{1}{2}\omega \Delta t}.\qquad\qquad\qquad\qquad\qquad(C6) \end{equation*} REFERENCES Aminzadeh F., 1996. “3-D salt and overthrust models”, in Applications of 3-D Seismic Data to Exploration and Production  chap. 28, pp. 247– 256, eds Weimer P., David T.L., AAPG/SEG, Tulsa. Amundsen L., Pedersen Ø., 2017. Time step n-tupling for wave equations, Geophysics  82 T249– T254. https://doi.org/10.1190/geo2017-0377.1 Google Scholar CrossRef Search ADS   Anderson J.E., Brytik V., Ayeni G., 2015. Numerical temporal dispersion corrections for broadband temporal simulation, RTM and FWI, in 85th Annual International Meeting  SEG, Expanded Abstracts, pp. 1096– 1100. Baysal E., Kosloff D.D., Sherwood J.W., 1983. Reverse time migration, Geophysics  48 1514– 1524. https://doi.org/10.1190/1.1441434 Google Scholar CrossRef Search ADS   Blanch J.O., Robertsson J. O.A., 1997. A modified Lax-Wendroff correction for wave propagation in media described by Zener elements, Geophys. J. Int.  131 381– 386. https://doi.org/10.1111/j.1365-246X.1997.tb01229.x Google Scholar CrossRef Search ADS   Chu C., Stoffa P.L., Seif R., 2009. 3D elastic wave modeling using modified high-order time stepping schemes with improved stability conditions, in 79th Annual International Meeting  SEG, Expanded Abstracts, pp. 2662– 2666. Dablain M., 1986. The application of high-order differencing to the scalar wave equation, Geophysics  51 54– 66. https://doi.org/10.1190/1.1442040 Google Scholar CrossRef Search ADS   Dai N., Liu H., Wu W., 2014. Solutions to numerical dispersion error of time FD in RTM, in 84th Annual International Meeting  SEG, Expanded Abstracts, pp. 4027– 4031. Dumbser M., Käser M., De La Puente J., 2007. Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D, Geophys. J. Int.  171 665– 694. https://doi.org/10.1111/j.1365-246X.2007.03421.x Google Scholar CrossRef Search ADS   Etgen J.T., Brandsberg-Dahl S., 2009. The pseudo-analytical method: application of pseudo-laplacians to acoustic and acoustic anisotropic wave propagation, in 79th Annual International Meeting  SEG, Expanded Abstracts, pp. 2552– 2556. Gardner G., Gardner L., Gregory A., 1974. Formation velocity and density—the diagnostic basics for stratigraphic traps, Geophysics  39 770– 780. https://doi.org/10.1190/1.1440465 Google Scholar CrossRef Search ADS   Kosloff D.D., Baysal E., 1982. Forward modeling by a Fourier method, Geophysics  47 1402– 1412. https://doi.org/10.1190/1.1441288 Google Scholar CrossRef Search ADS   Li Y., Wong M., Clapp R., 2013. Equivalent accuracy at a fraction of the cost: overcoming temporal dispersion, Tech. Rep. 150, Stanford Exploration Project. Li Y.E., Wong M., Clapp R., 2016. Equivalent accuracy at a fraction of the cost: overcoming temporal dispersion, Geophysics  81 T189– T196. https://doi.org/10.1190/geo2015-0398.1 Google Scholar CrossRef Search ADS   Liu H., Dai N., Niu F., Wu W., 2014. An explicit time evolution method for acoustic wave propagation, Geophysics  79 T117– T124. https://doi.org/10.1190/geo2013-0073.1 Google Scholar CrossRef Search ADS   Liu Y., 2013. Globally optimal finite-difference schemes based on least squares, Geophysics  78 T113– T132. https://doi.org/10.1190/geo2012-0480.1 Google Scholar CrossRef Search ADS   Liu Y., 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling, Geophys. J. Int.  197 1033– 1047. https://doi.org/10.1093/gji/ggu032 Google Scholar CrossRef Search ADS   Mittet R., 2017. On the internal interfaces in finite-difference schemes, Geophysics  82 T159– T182. https://doi.org/10.1190/geo2016-0477.1 Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Gális M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures  Cambridge University Press. Google Scholar CrossRef Search ADS   Qin Y., Quiring S., Nauta M., 2017. Temporal dispersion correction and prediction by using spectral mapping, in 79th Annual International Conference & Exhibition 2017  EAGE, Extended Abstracts, doi:10.3997/2214-4609.201700677. Song X., Fomel S., Ying L., 2013. Lowrank finite-differences and lowrank fourier finite-differences for seismic wave extrapolation in the acoustic approximation, Geophys. J. Int.  193 960– 969. https://doi.org/10.1093/gji/ggt017 Google Scholar CrossRef Search ADS   Soubaras R., Zhang Y., 2008. Two-step explicit marching method for reverse time migration, in 78th Annual International Meeting  SEG, Expanded Abstracts, pp. 2272– 2276. Stork C., 2013. Eliminating nearly all dispersion error from FD modeling and RTM with minimal cost increase, in 75th Annual International Conference & Exhibition  EAGE, Extended Abstracts, doi:10.3997/2214-4609.20130478. Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics  49 1259– 1266. https://doi.org/10.1190/1.1441754 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics  51 1954– 1966. https://doi.org/10.1190/1.1442051 Google Scholar CrossRef Search ADS   Tromp J., Komatitsch D., Liu Q., 2008. Spectral-element and adjoint methods in seismology, Commun. Comput. Phys.  3 1– 32. Virieux J., 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method, Geophysics  51 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Wang M., Xu S., 2015a. Finite-difference time dispersion transforms for wave propagation, Geophysics  80 WD19– WD25. https://doi.org/10.1190/geo2015-0059.1 Google Scholar CrossRef Search ADS   Wang M., Xu S., 2015b. Time dispersion prediction and correction for wave propagation, in 85th Annual International Meeting  SEG, Expanded Abstracts, pp. 3677– 3681. Wang Y., Liang W., Nashed Z., Li X., Liang G., Yang C., 2014. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method, Geophysics  79 T277– T285. https://doi.org/10.1190/geo2014-0078.1 Google Scholar CrossRef Search ADS   Xu Z., Jiao K., Cheng X., Sun D., King R., Nichols D., Vigh D., 2017. Time-dispersion filter for finite-difference modeling and reverse time migration, in 87th Annual International Meeting  SEG, Expanded Abstracts, pp. 4448– 4452. Yao G., Wu D., Debens H.A., 2016. Adaptive finite difference for seismic wavefield modelling in acoustic media, Sci. Rep.  6, doi:10.1038/srep30302. SUPPORTING INFORMATION Supplementary data are available at GJI online. Time_dispersion_1D_examples.m ITDT.m FTDT.m FD_elastic_anisotropic.tar.gz Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: FIRST-ORDER SYSTEM OF EQUATIONS We will first show how the factor θ can be recovered from the finite-difference operator itself. Note that in general, an nth-order derivative is written as $$\mathcal {F}\left[{\partial ^n {u}(t)}/{\partial t^n}\right] = (-i\omega )^n{\tilde{u}}(\omega )$$. Phase shifts are denoted as $$\mathcal {F}[u(t+\Delta t)]=e^{-i\omega \Delta t}\tilde{u}(\omega )$$. The temporal FD operator in the Fourier domain is thus:   \begin{eqnarray} \mathcal {F}\left[ \frac{{u}(t\!-\!\Delta t)\!-\!2{u}(t)\!+\!{u}(t\!+\!\Delta t)}{\Delta t^2} \right]\!\!=\!\left(\!\frac{e^{-i\omega \Delta t}\!-2\!+\!e^{i\omega \Delta t}}{\Delta t^2}\!\right)\tilde{u}(\omega ),\nonumber\\ \end{eqnarray} (A1)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\qquad= -\left[\frac{2-2\cos (\omega \Delta t)}{\Delta t^2}\right]\tilde{u}(\omega ),\nonumber\\ \end{eqnarray} (A2)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\quad= -\left[\frac{\sin (\frac{\omega \Delta t}{2})}{\frac{1}{2}\omega \Delta t} \right]^2\omega ^2 \tilde{u}(\omega ), \end{eqnarray} (A3)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\,\,\,\,\,\,= -\theta ^2\omega ^2 \tilde{u}(\omega ). \end{eqnarray} (A4)We note that −ω2 represents the second-order time derivative, and the additional factor θ2 corresponds to the undesired remaining Taylor expansion terms. We use this method to show the equivalence in the first-order system of equations solved in a leap-frog manner (Virieux 1986):   \begin{eqnarray} \mathcal {F}\left[\frac{u\left(t+\frac{\Delta t}{2}\right)-u\left(t-\frac{\Delta t}{2}\right)}{\Delta t} \right] &=& \left(\frac{e^{-i\omega \Delta t/2}-e^{i\omega \Delta t/2}}{\Delta t}\right)\tilde{u}(\omega ), \end{eqnarray} (A5)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=& -\left( \frac{2i\sin \left(\frac{\omega \Delta t}{2}\right)}{\Delta t} \right)\tilde{u}(\omega ), \end{eqnarray} (A6)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=& -\left( \frac{\sin \left(\frac{\omega \Delta t}{2}\right)}{\frac{1}{2}\omega \Delta t} \right)i\omega \tilde{u}(\omega ), \end{eqnarray} (A7)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,&=&-\theta i\omega \tilde{u}(\omega ). \end{eqnarray} (A8)We note that −iω represents the first-order time derivative, and the additional factor θ corresponds to the undesired remaining Taylor expansion terms. The factor θ is identical in both derivations, representing the same acceleration of the wavefield in every calculation. In other words, the time-dispersion error is identical when solving either the first- or second-order system of equations. The results derived in this paper apply thus equally to both cases. APPENDIX B: FOURTH-ORDER ACCURATE TIME DISPERSION B1 Taylor-series fourth-order accurate scheme The method developed in this paper may be extended to any order time stepping. We demonstrate this for a fourth-order accurate scheme in time, where we solve a system based on eq. (2), the Taylor series of the FD operator:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},t\!-\!\Delta t)\!-\!2\mathbf {u}(\mathbf {x},t)\!+\!\mathbf {u}(\mathbf {x},t\!+\!\Delta t)}{\Delta t^2} = \sum _{k=1}^\infty \frac{\Delta t^{2(k-1)}}{\frac{1}{2}(2k)!}\frac{\partial ^{2k}\mathbf {u}(t)}{\partial t^{2k}}, \end{eqnarray} (B1)  \begin{eqnarray} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\approx \left(L+ \frac{\Delta t^2L^2}{\frac{1}{2}4!}\right)\mathbf {u}(\mathbf {x},t).\nonumber\\ \end{eqnarray} (B2)In other words, the Taylor series of the temporal FD operator is truncated, and we substitute spatial operators using the wave equation to eliminate higher order temporal error terms. For this time-stepping scheme, we get the following relation between the time-dispersed (ω) and original frequencies (ω0):   $$\omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t\sqrt{1 - \frac{(\omega _0\Delta t)^2}{\frac{1}{2}4!}}}{2} \right) .$$ (B3)The square-root term corresponds to θ from eq. (3), truncated to the two terms solved in the fourth-order accurate scheme. Kinematically speaking, the second-order FD operator, that is, left-hand side of eq. (B2), accelerates the wavefield as discussed in the paper, while the fourth-order accurate scheme, that is, right-hand side of eq. (B2), decelerates the original wavefield, leading to a reduced overall time dispersion. Solving the system for ω0 gives:   \begin{eqnarray} \omega _0 = \frac{1}{\Delta t}\sqrt{6-\sqrt{12+24\cos (\omega \Delta t)}},\\ \ \ \ \ \ {\rm for }\quad\omega \Delta t\le \frac{2\pi }{3}\approx 2.1 \nonumber . \end{eqnarray} (B4)One can substitute identities (B3) and (B4) in the ITDT and FTDT, replacing eqs (13a) and (13b), to obtain an ITDT and FTDT for the fourth-order scheme. The ITDT based on eq. (B3) is equal to that found in Wang & Xu (2015a,b), who do not report the limitation ωΔt ≤ 2.1, which may cause unintended leaking of low frequencies into the higher signal spectrum. B2 Optimized-series fourth-order accurate scheme The advantage of a fourth-order accurate scheme in time is reduced time dispersion and a larger CFL condition. The Taylor series in eq. (B1) may be replaced with a more efficient series expansion that achieves an even larger CFL condition, at the cost of accumulating more time dispersion. Soubaras & Zhang (2008) and Amundsen & Pedersen (2017) describe such a procedure that achieves twice the CFL condition using the following series expansion:   \begin{eqnarray} \frac{\mathbf {u}(\mathbf {x},t\!-\!\Delta t)\!-\!2\mathbf {u}(\mathbf {x},t)\!+\!\mathbf {u}(\mathbf {x},t\!+\!\Delta t)}{\Delta t^2} \approx \left(L+ \frac{\Delta t^2L^2}{16}\right)\mathbf {u}(\mathbf {x},t).\nonumber\\ \end{eqnarray} (B5)Comparing eqs (B5) with (B2), we observe that the single difference between the two schemes is the factor 16. For this scheme, we may similarly propose the FTDT and ITDT, using identities:   \begin{eqnarray} \omega = \frac{2}{\Delta t}\sin ^{-1}\left(\frac{\omega _0\Delta t\sqrt{1 - \frac{(\omega _0\Delta t)^2}{16}}}{2} \right) , \end{eqnarray} (B6a)  \begin{eqnarray} \omega _0= \frac{4}{\Delta t}\sin \left( \frac{\omega \Delta t}{4} \right). \end{eqnarray} (B6b) Note that eq. (B6a) is valid for $$\omega _0\Delta t \le 4\sin (\pi /4)=2\sqrt{2}\approx 2.8$$, related to the Nyquist sampling explained in Section 2.4. We show a simple example, using the pseudospectral method. The standard CFL limit is 2/π, whereas this fourth-order scheme allows stable computations using a CFL limit of 4/π. Fig. A1 shows a recording of a 25 Hz central frequency Ricker wavelet, using the same model as in the 1-D examples in this paper: Δx = 5 m, c = 2000 m s−1 and the distance between source and receiver is 1 km, but we use time steps at double the CFL limit of 1.2. As the PS method is spectrally accurate in space, only time dispersion accumulates. Using the proposed time-dispersion workflow proposed in this paper, we recover the Ricker wavelet without dispersion. Figure A1. View largeDownload slide Recording due to a pseudospectral, modified fourth-order accurate in time scheme, eq. (B5) (blue), and the derived FTDT and ITDT identities of eq. (B6) restore the original input Ricker wavelet (red). Figure A1. View largeDownload slide Recording due to a pseudospectral, modified fourth-order accurate in time scheme, eq. (B5) (blue), and the derived FTDT and ITDT identities of eq. (B6) restore the original input Ricker wavelet (red). APPENDIX C: CORRESPONDENCE BETWEEN SOLUTIONS SATISFYING THE FINITE-DIFFERENCE SYSTEM AND THE WAVE EQUATION In this section, we will show the relation between continuous solutions satisfying the wave equation and those satisfying the finite-difference (FD) system. Let u be a solution to the wave equation:   $$\mathcal {A}u(x,t) \stackrel{{\mbox{def}}}{=}\left( \frac{\partial ^2}{\partial t^2} - \mathcal {L} \right)u(x,t) = f(x,t),$$ (C1)where t > 0 denotes time, $$x\in \mathbb {R}^d$$, f describes the source term and $$\mathcal {L}$$ is a spatial operator independent of time. Now consider a related solution v to the FD system:   $$\mathsf {A}v(x,t) \stackrel{{\mbox{def}}}{=}(\mathsf {D}-\mathcal {L})v(x,t) = g(x,t),$$ (C2)where $$\mathsf {D}v(x,t) = \sum _n c_n v(x,t-n\Delta t)$$, with coefficients cn chosen such that $$\mathsf {D}$$ approximates the second-order time derivative. We will show that we can write v(x, t) as a function of u(x, t), to give a relation between solutions satisfying the FD system and those satisfying the wave equation. Applying a Fourier transform to eq. (C2) yields:   \begin{eqnarray} \int _{-\infty }^\infty \mathsf {A}v(x,t)e^{-i\omega t}\textrm{d}t = \left(\sum _{n}c_n e^{i \omega \Delta t n} - \mathcal {L}\right)\widehat{v}(x,\omega ), \end{eqnarray} (C3)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \left( -\omega ^2\theta ^2 - \mathcal {L}\right)\widehat{v}(x,\omega ), \end{eqnarray} (C4)using:   $$\theta (\omega ) = \sqrt{\frac{\sum _n c_n e^{i\omega \Delta t n}}{-\omega ^2}},$$ (C5)where it is assumed that Δt is small enough such that the quantity in the square root is positive.1 Now try as ansatz the following relation between the solution v of the FD system and u of the actual wave equation:   \begin{eqnarray} v(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{u}(x,\omega \theta ) e^{i\omega t}\textrm{d}\omega , \end{eqnarray} (C7)with Fourier transform:   \begin{eqnarray} \widehat{v}(x,\omega ) = \widehat{u}(x,{\omega }{\theta }). \end{eqnarray} (C8) Continuing eq. (C4):   \begin{eqnarray} \int _{-\infty }^\infty \mathsf {A}v(x,t) e^{-i\omega t}\textrm{d}t = ( -\omega ^2\theta ^2 - \mathcal {L})\widehat{v}(x,\omega ), \end{eqnarray} (C9)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \left( -\omega ^2\theta ^2 - \mathcal {L} \right)\widehat{u}(x,\omega \theta ), \end{eqnarray} (C10)  \begin{eqnarray} \qquad\qquad\qquad\qquad\,\,= \widehat{f}(x,\omega \theta ). \end{eqnarray} (C11)The last step is due to the definition that eq. (C10) is a solution to the wave equation (C1). Thus, if f and $$\mathcal {L}$$ are regular enough such that the Fourier transform u is well defined, we can use definitions:   $$v(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{u}(x,\omega \theta ) e^{{i\omega t}}{\rm d}\omega ,$$ (C12)and   $$g(x,t) = \frac{1}{2\pi }\int _{-\infty }^\infty \widehat{f}(x,\omega \theta )e^{{i\omega t}}\textrm{d}\omega ,$$ (C13)that satisfy the FD system for all t:   $$\mathsf {A}v(x,t) = g(x,t).$$ (C14)The solution in the FD system can thus be written as phase-shifted solutions to the true wave equation, independent of the spatial location x or the exact form of operator $$\mathcal {L}$$. Numerically solving the FD system may cause aliasing, in which case this continuous relation does not hold for the discrete situation. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations