# An arbitrary-order staggered time integrator for the linear acoustic wave equation

An arbitrary-order staggered time integrator for the linear acoustic wave equation Summary We suggest a staggered time integrator whose order of accuracy can arbitrarily be extended to solve the linear acoustic wave equation. A strategy to select the appropriate order of accuracy is also proposed based on the error analysis that quantitatively predicts the truncation error of the numerical solution. This strategy not only reduces the computational cost several times, but also allows us to flexibly set the modelling parameters such as the time step length, grid interval and P-wave speed. It is demonstrated that the proposed method can almost eliminate temporal dispersive errors during long term simulations regardless of the heterogeneity of the media and time step lengths. The method can also be successfully applied to the source problem with an absorbing boundary condition, which is frequently encountered in the practical usage for the imaging algorithms or the inverse problems. Numerical approximations and analysis, Computational seismology, Wave propagation 1 INTRODUCTION Numerical simulation of wave propagation is not only meaningful for forward problem itself but also plays an important role in several applications, such as reverse time migration for subsurface imaging (e.g. Baysal et al. 1983; Whitmore 1983) and full waveform inversion for reconstruction of the subsurface parameters (e.g. Tarantola 1984; Pratt et al. 1998; Shin & Cha 2008, 2009). The more accurate the modelling scheme, the higher the resolution of the obtained image. However, since discretization of wave equations inevitably leads to numerical errors, that is, phase and amplitude errors, various methodologies have been developed to improve the accuracy of numerical modelling of wave phenomena. In the geophysics literature, Virieux (1984) introduced the staggered-grid finite-difference method (FDM) by adopting the method of Yee (1966) to simulate the seismic wave equation represented by a first-order system of the stress and the particle velocity fields. This method defines the wavefields and material properties in the staggered grids separately and updates the stress and the particle velocity fields alternatively. Such gridding suppresses the even-odd decoupling and keeps the numerical solution from high-frequency oscillation. Increasing the spatial order of accuracy have been usually attempted to reduce the spatial dispersion error in the context of FDMs (e.g. Levander 1988). Such methods enlarge the grid stencil and incur additional treatments to impose boundary conditions. The compact-stencil FDM (Lele 1992) can also be used to efficiently yield an accurate solution. The method is a generalization of Padé scheme implicitly calculating the high-order spatial derivatives by solving tridiagonal or pentadiagonal matrix (three-point or five-point stencil), which can simply be solved with the complexity of O(N). This approach can yield more accurate solutions with a smaller finite difference stencil than that of the standard FDM. By optimizing the finite difference coefficients to conform the physical property of the wave phenomena, the dispersive and dissipative errors can be minimized (Ashcroft & Zhang 2003). Similar dispersion-relation-preserving (DRP) techniques for the explicit FDM have been also introduced using standard stencil (e.g. Tam & Webb 1993; Bogey & Bailly 2004; Liu & Sen 2009). The DRP techniques significantly reduce the dispersion error of the solution compared to the standard FDM with the equivalent stencil and the order of accuracy. For a deeper exposition on seismic modelling with FDM, look for Moczo et al. (2014). Schemes based on the finite-element method (FEM) increase the accuracy of the scheme by using a number of high-order basis polynomials locally defined in an element to approximate a solution. In the field of seismology, the spectral element method (SEM), a type of FEM technique, is used for the global seismic modelling (Komatitsch & Vilotte 1998). This method usually uses the high-order Gauss–Lobatto collocation nodes to define nodal (or modal) basis functions, which generate a diagonal global mass matrix naturally, allowing effective massive time-domain modelling. De Basabe et al. (2008) compared the dispersion characteristics of the discontinuous Galerkin (DG) method with that of the SEM. DG schemes are non-conforming methods that evaluate the wave solution of each element separately and compensate the discontinuity of the solution with numerical flux terms, such as Godnov or Lax–Friedrich flux (Arnold et al. 2002; Hesthaven & Warburton 2007). The methods can be used for special purposes when conducting simulation that contains complex geometry, that is, models with complex topography or interfaces of medium of different phases (Hermann et al. 2011), which is almost impossible to solve with conventional FDMs. Pseudospectral methods usually offer spectral accuracy in the spatial domain, which implies that the scheme is dispersion-free within the expressible band of wavenumber (Trefethen 2000). A type of basis function for the method depends on the boundary conditions of the problem. The Fourier basis functions are generally used which intrinsically assume the periodic boundary condition for each spatial dimension (Kosloff & Baysal 1982; Kosloff & Kosloff 1983). In order to impose general boundaries, such as Dirichlet or Neumann boundaries, the Chebyshev polynomial basis are required (Carcione 1996) and research has been conducted to make the stability condition comparable to that of the Fourier pseudospectral method (Renaut & Fröhlich 1996). Recently, Klin et al. (2010) suggested simple but effective method to impose a free surface boundary condition when using the Fourier pseudospectral method, which was verified by simulation of 3-D realistic models. Either Fourier or Chebyshev pseudospectral method incurs a pair of forward and inverse Fourier transforms to calculate spatial derivatives implicitly likewise the compact-stencil FDM approach. This can be accelerated by using the fast Fourier transform algorithm, which has the complexity of O(N log N). In general, the higher the spatial accuracy of a numerical method, the more strict the stability condition. The stability condition can be relaxed by matching the temporal order of accuracy with the spatial one. The Lax–Wendroff method can increase the temporal order of accuracy by converting the higher order time derivatives obtained from the Taylor series to the spatial derivatives. Temporal fourth-order scheme (Cohen 1986; Dablain 1986) was implemented to solve the acoustic wave equation by the Lax–Wendroff method, which has also been applied to derivation of ADER-DG schemes (Dumbser & Käser 2006; Käser & Dumbser 2006). Tal-Ezer (1986) suggested another approach to increase the temporal accuracy. This method expands the exponential map of the characteristic matrix of the hyperbolic system to the series of product of the Bessel function and the Chebyshev polynomial by using the Jacobi-Anger expansion. In the context of geophysics, the methods using the Jacobi-Anger expansion as the time stepping method are referred to as the rapid expansion method. This method was originally used to solve a single source problem (Kosloff et al. 1989) and then was extended to multisource problems for the practical applications (Pestana & Stoffa 2010; Tessmer 2011). In this study, we suggest a staggered time integrator which can extend the order of accuracy arbitrarily using the Lax–Wendroff method to solve the linear acoustic wave system. This method incorporates several low-order schemes such as Long et al. (2013) and Tan & Huang (2014). The proposed numerical method can be implemented easily by using a recurrence relation which enables the numerical solution to be almost non-dispersive in temporal dimension. In addition to the temporal accuracy, spectral accuracy in the spatial domain can also be achieved by spatial differentiations in the wavenumber domain via spatial Fourier transform as in the pseudospectral method. We faithfully carry the derivation process and implementation of the method in Theory part. The proposed method is analysed in the consecutive section. Simple strategy on selecting a proper temporal order of accuracy is also stated based on the qualitative expression of the error between the numerical and the exact wave solution of the initial value problem. Then, we verify the proposed modelling method following the strategy with the homogeneous P-wave velocity model. By numerical simulations, we also examine the applicability for the practical usage such as a problem with heterogeneity or a source problem with the absorbing boundary condition. 2 THEORY The acoustic wave equation in the first-order system consists of the linearized conservation relations with respect to the pressure perturbation p and the particle velocity vector v as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ \mathbf {v} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}0 & -\rho _0 c^2 \nabla \cdot \\ -{\rho _0}^{-1}\nabla & 0 \end{array}\right]}}{\left[\begin{array}{c}p \\ \mathbf {v} \end{array}\right]}, \end{eqnarray} (1)where ρ0 and c are the density and the P-wave propagation speed, respectively. Temporal discretization of such first-order system is usually done with the leap-frog or the two-stage Störmer–Verlet method [Maxwell’s equations (Yee 1966), the elastic wave equation (Virieux 1986), etc.], whose structural preserving property bounds the error of Hamiltonian in some level (Chen 2009). In this study, we impose each physical value on the staggered time axis; p on the integer and v on the half integer of the temporal grid. Consider v(t + Δt/2) with respect to v(t) using the Taylor expansion as   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}(t) + \frac{\Delta t}{2}\frac{\partial \mathbf {v}}{\partial t}\bigg \vert _{t} + \frac{\Delta t^2}{8}\frac{\partial ^2\mathbf {v}}{\partial t^2}\bigg \vert _{t} + \cdots + \frac{\Delta t^n}{n! \, 2^n}\frac{\partial ^n\mathbf {v}}{\partial t^n}\bigg \vert _{t} + \cdots , \end{eqnarray} (2)and v(t − Δt/2) as   \begin{eqnarray} \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) = \mathbf {v}(t) - \frac{\Delta t}{2}\frac{\partial \mathbf {v}}{\partial t}\bigg \vert _{t} + \frac{\Delta t^2}{8}\frac{\partial ^2\mathbf {v}}{\partial t^2}\bigg \vert _{t} + \cdots + \frac{(-\Delta t)^n}{n! \, 2^n}\frac{\partial ^n\mathbf {v}}{\partial t^n}\bigg \vert _{t} + \cdots . \end{eqnarray} (3) We can obtain the equation of v(t + Δt/2) by subtracting (3) from (2) as follows:   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) + \sum \limits _{n=0}^\infty \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\frac{\partial ^{(2n+1)}\mathbf {v}}{\partial t^{(2n+1)}}\bigg \vert _{t}{\, }. \end{eqnarray} (4) In a similar manner, we can derive the equation of p(t + Δt) with respect to p(t) and $$\partial _t^{(2n+1)}p(t+dt/2)$$ as follows:   \begin{eqnarray} p(t+\Delta t) = p(t) + \sum \limits _{n=0}^\infty \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\frac{\partial ^{(2n+1)} p}{\partial t^{(2n+1)}}\bigg \vert _{t+{\Delta t}/{2}}{\, }. \end{eqnarray} (5) Following the Lax–Wendroff method, the series of time derivatives of (4) and (5) can be replaced by the spatial derivatives using (1). Then, the time stepping procedure can be expressed as follows with an arbitrary order of 2l + 2 as   \begin{eqnarray} \begin{array}{lcl} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t) + O(\Delta t^{2l+2}), \\ p(t+\Delta t) = p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + O(\Delta t^{2l+2}){\, }. \end{array} \end{eqnarray} (6) It is noted that (6) is the extended formulation of the leap-frog or the two-stage Störmer–Verlet method. If l is zero, the method is equivalent to the conventional second-order staggered-time method (Yee 1966; Virieux 1986). The fourth-order method of Long et al. (2013) (l = 1) and the sixth-order method of Tan & Huang (2014) (l = 2) are subsets of (6). The structural preserving property is also preserved for (6) that no dissipative error occurs during the simulation, intrinsically. Only dispersive error may debase the numerical solution of (6), which can effectively be reduced by increasing the temporal order of accuracy of the scheme. The spatial derivative operators can be calculated via several methods. In this study, we adopt the Fourier pseudospectral method which gives the almost exact spatial derivatives for sufficiently smooth wavefields (Trefethen 2000) in order to correspond to the temporal accuracy of (6). The x-directional spatial derivative of f via the Fourier pseudospectral method is implemented by using the discrete Fourier transform as   \begin{eqnarray} \frac{\partial f}{\partial x} = \mathcal {F}_x^{-1} ( -ik_x{\, }{\, } \mathcal {F}_x (f) ), \end{eqnarray} (7)where $$\mathcal {F}_x$$ and kx are the Fourier transform operator and the wavenumber in the x-direction. The proposed method (6) of desired order of accuracy can be effectively calculated by using a simple recurrence relation. If we define scalar fields Al and an as   \begin{eqnarray} A_l = \sum \limits _{n=0}^l a_n = \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t), \end{eqnarray} (8)the recurrence relation can be simply expressed as   \begin{eqnarray} a_n = \frac{\Delta t^2}{8n{\, }(2n+1)}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )a_{n-1}, \quad n=1,2, \ldots ,l \end{eqnarray} (9)where   \begin{eqnarray} a_0 = \Delta t {\, } p(t){\, }. \end{eqnarray} (10) Then v(t + Δt/2) is calculated as follows:   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \frac{1}{\rho _0}\nabla A_l . \end{eqnarray} (11)In the same way, we can define scalar fields Bl and bn as follows:   \begin{eqnarray} B_l = \sum \limits _{n=0}^l b_n = \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ). \end{eqnarray} (12)The recurrence relation for bn is equivalent to (9) as   \begin{eqnarray} b_n = \frac{\Delta t^2}{8n{\, }(2n+1)}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )b_{n-1}, \quad n=1,2, \ldots ,l , \end{eqnarray} (13)where   \begin{eqnarray} b_0 = {\Delta t}{\, }{\rho _0 c^2}\nabla \cdot \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ). \end{eqnarray} (14) Then p(t + Δt) can be calculated as follows:   \begin{eqnarray} p(t+{\Delta t}) = p(t) - B_l . \end{eqnarray} (15) We can increase the temporal order of accuracy of (6) by two by adding the subsequent term of the series of the spatial derivatives by using (9) and (13), which only requires to apply the second-order spatial derivative operator ρ0c2∇ · ρ0−1∇ to update the sequence. The operator appears in the scalar wave equation considering the heterogeneity of ρ0 in the domain. The operator is applied to a scalar field f via the Fourier pseudospectral method in 2-D domain as follows:   \begin{eqnarray} {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla f = {\rho _0 c^2} \Big ( \mathcal {F}_x^{-1} \Big ( -ik_x{\, } \mathcal {F}_x \Big ( \frac{1}{\rho _0} \mathcal {F}_x^{-1}\Big ({-}ik_x{\, } \mathcal {F}_x (f) \Big ) \Big ) \Big ) + \mathcal {F}_z^{-1} \Big ({-}ik_z{\, } \mathcal {F}_z \Big ( \frac{1}{\rho _0} \mathcal {F}_z^{-1}\Big (-ik_z{\, } \mathcal {F}_z (f) \Big ) \Big ) \Big ) \Big ). \end{eqnarray} (16) To increase the index of an or bn by one, two pairs of Fourier and inverse Fourier transform are required for each direction to consider the heterogeneity of ρ0. If ρ0 is constant elsewhere, the operator is simplified to the Laplacian, that is, c2∇2, which is calculated by the Fourier pseudospectral method as   \begin{eqnarray} {c^2}\nabla ^2 f = {c^2} \left ( \mathcal {F}_x^{-1} \left ( -k_x^2{\, } \mathcal {F}_x (f) \right ) + \mathcal {F}_z^{-1} \left ( -k_z^2{\, } \mathcal {F}_z (f) \right ) \right ). \end{eqnarray} (17) In this case, the required number of pair of the Fourier and inverse Fourier transform is reduced by half compared to (16). 3 ANALYSIS In this section, we manipulate (6) to demonstrate how we can control the accuracy of the numerical simulation result under given conditions. This section is related to both dispersion error and the stability analysis of the scheme. If we rewrite the scheme (6) with respect to p, we can obtain the following relation.   \begin{eqnarray} p(t+\Delta t) + p(t-\Delta t) = \bigg ( 2 + {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \bigg ( \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\bigg ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \bigg )^{n} \bigg )^2 \bigg ) p(t). \end{eqnarray} (18) For the analysis, we consider the initial value problem with p(0) = p0 and v(0) = 0. As the first step of the time marching, the appropriate calculation of v(Δt/2) is needed. The exact value of v(Δt/2) cannot be obtained by using (6) but directly discretizing (1) using the Lax–Wendroff method as   \begin{eqnarray} \mathbf {y}\Big (\frac{\Delta t}{2}\Big ) = \mathbf {y}(0) + \sum \limits _{n=1}^{2l+2} \frac{\Delta t^n\mathbf {M}^n}{n!{}\, 2^n}\frac{\partial ^{n}\mathbf {y}}{\partial t^n}\bigg \vert _{t=0} + O(\Delta t^{2l+3}), \end{eqnarray} (19)where y(t) = [p(t), v(t) ]⊤ and M is the characteristic matrix of (1). Since v(0) = 0, v(Δt/2) is expressed using p0 with an equivalent temporal order of accuracy as follows:   \begin{eqnarray} \mathbf {v}\left ( \frac{\Delta t}{2} \right ) = \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p_0 + O(\Delta t^{2l+2}). \end{eqnarray} (20)Then we can calculate p(Δt) using (6) and (20) as follows:   \begin{eqnarray} \displaystyle p\left(\Delta t \right) = p_0 + 2{\, }{\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \left( \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \left( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n} \right)^2 p_0. \end{eqnarray} (21)As a result, we now have the relation of p with respect to the number of time steps m as   \begin{eqnarray} \begin{array}{lcl}p\left( 0 \right) = p_0, \\ \displaystyle p\left(\Delta t \right) = {\rm C}(l){\, }p_0, \\ \displaystyle p(m\Delta t + \Delta t) = 2 {\rm C}(l){\, }p(m \Delta t) - p(m\Delta t - \Delta t), \quad m=1, 2, \ldots , \end{array} \end{eqnarray} (22)where   \begin{eqnarray} {\rm C}(l) = 1 + 2{\, }{\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \left( \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \left( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n} \right)^2. \end{eqnarray} (23) Relation (22) is equivalent to the three-term recurrence formulation of the Chebyshev polynomials of the first kind, denoted Tn(x) with the degree of the polynomials of n. Therefore, we can express p at mth time step as follows:   \begin{eqnarray} p(m \Delta t) = T_m( {\rm C}(l) ){\, }p_0, \quad m=0, 1, 2, \ldots . \end{eqnarray} (24) For the further analysis, we set c and ρ0 to be constant elsewhere. Then ρ0c2∇ · ρ0−1∇ is expressed as c2∇2. Since we use the pseudospectral method in this study, C(l) is implemented in the wavenumber domain as   \begin{eqnarray} {\rm \widetilde{C}}(l,\theta ) = 1 - 2\left( \sum \limits _{n=0}^{l} \frac{ (-1)^n }{(2n+1)!{\, }2^{2n+1}} \theta ^{2n+1} \right)^2, \end{eqnarray} (25) where $${\rm \widetilde{C}}(l,\theta )$$ is C(l) to which the spatial Fourier transform is applied and θ = ckΔt where $$k=\sqrt{k_x^2+k_z^2}$$ in the 2-D domain. If we consider the relation   \begin{eqnarray} \cos (\theta ) = 1 - \displaystyle 2\sin ^2\left( \frac{\theta }{2} \right) \approx 1 - 2\left( \displaystyle \sum \limits _{n=0}^{l} \frac{ (-1)^n }{(2n+1)!{\, }2^{2n+1}} \theta ^{2n+1} \right)^2, \end{eqnarray} (26)we can rewrite $${\rm \widetilde{C}}(l,\theta )$$ as the Taylor expansion of cos (θ) as follows:   \begin{eqnarray} {\rm \widetilde{C}}(l,\theta ) = \sum \limits _{n=0}^{2l+1} \frac{ (-1)^n }{(2n)!} \theta ^{2n}. \end{eqnarray} (27)Then p at the mth time step in the wavenumber domain can be expressed by applying spatial Fourier transform to (24) as   \begin{eqnarray} \widetilde{p}(m \Delta t) = T_m( {\rm \widetilde{C}}(l,\theta ) ){\, }\widetilde{p}_0, \quad m=0, 1, 2, \ldots ,\end{eqnarray} (28)where $$\widetilde{p}$$ is p defined in the wavenumber domain. Eq. (28) explicitly shows that the stability condition of the proposed scheme is $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$, otherwise, $$T_m( {\rm \widetilde{C}}(l,\theta ) )$$ diverges as m increases. Such condition is equivalent to that derived by the Von Neumann or Fourier analysis which sets the pressure field as the plane wave expressed in the harmonic function (Strikwerda 2004). Fig. 1 illustrates the region of θ marked in black where the stability criterion is violated with respect to the temporal order of accuracy. The stability region gets enlarged as the order of accuracy in time increases, which is referred as the relaxation of the Courant–Friedrichs–Lewy (CFL) condition. The small regions in the middle are usually located at the odd multiple of π, which are marked because the Taylor series inaccurately approximate the cosine function and are slightly larger than one at these points. Although such a small error causes divergence of the solution, it can be acceptable depending on the modelling conditions. In order to determine the acceptability of the modelling condition or to estimate the quality of the result of the numerical simulation, more informative condition than the stability criterion is needed. Figure 1. View largeDownload slide Regions violating $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$ are marked in black. The wave components of these regions get diverging as m increases. Figure 1. View largeDownload slide Regions violating $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$ are marked in black. The wave components of these regions get diverging as m increases. If l of (27) goes to infinity, $${\rm \widetilde{C}}(l,\theta )$$ converges to cos (θ). Then, the amount of the truncation error of the numerical solution of a wave component θ at the mth time step can be written as follows:   \begin{eqnarray} {\rm error}(m,l,\theta ) = \left( T_m({\rm \widetilde{C}}(l,\theta )) - T_m(\cos (\theta )) \right)\widetilde{p}_0. \end{eqnarray} (29) Relation (29) is a general expression of the misfit between the numerical result and the analytic solution expressing not only dispersion error due to the time discretization, but also the amplification error when the stability condition is violated. If we set a condition as   \begin{eqnarray} |T_m({\rm \widetilde{C}}(l,\theta )) - T_m(\cos (\theta ))| \le \epsilon , \end{eqnarray} (30)we can obtain θmax, the maximum allowable θ with respect to l satisfying the criterion (30) after m time steps. In this study, we set the maximum m and ε as 50 000 and 0.01, respectively, which implies that the error of the solution in the wavenumber domain where θmax is less than or equal to 1 per cent during 50 000 time steps. Table 1 represents θmax of our error criterion with respect to the temporal order of accuracy from 2 to 28. As the order of accuracy increases, θmax is extended so that more relaxed modelling conditions, such as larger Δt, faster wave propagation speed or smaller grid size, are allowed. In Fig. 2, θmax is represented by grey dots and shown together with the instability region marked in black in Fig. 1 for exact comparison. It is shown that θmax tends to increase linearly with respect to the order of accuracy. The error conditions of Table 1 is much stricter than the stability conditions because most grey dots are located at far more left positions to the instability regions. If less restrictive error condition, such as much smaller m or larger ε, is used, θmax will increase and may move to the right position of the small black region which causes the instability in Fig. 2. Figure 2. View largeDownload slide The grey dots are the maximum allowable θ satisfying eq. (30) with respect to the order of accuracy in temporal dimension when m and ε are 50 000 and 0.01, respectively. Figure 2. View largeDownload slide The grey dots are the maximum allowable θ satisfying eq. (30) with respect to the order of accuracy in temporal dimension when m and ε are 50 000 and 0.01, respectively. Figure 3. View largeDownload slide Computational domain of the numerical simulation whose size is 4π × 4π km2. Homogeneous velocity of 5 km s−1 and density of 1 $$\text{g {}cm}^{-3}$$ is used with a grid size of 0.0491 km. Four black squares illustrate receivers to record time traces located at π km for z and 2π (R1), 2.5π (R2), 3π (R3) and 3.5π km (R4) for x, respectively. Figure 3. View largeDownload slide Computational domain of the numerical simulation whose size is 4π × 4π km2. Homogeneous velocity of 5 km s−1 and density of 1 $$\text{g {}cm}^{-3}$$ is used with a grid size of 0.0491 km. Four black squares illustrate receivers to record time traces located at π km for z and 2π (R1), 2.5π (R2), 3π (R3) and 3.5π km (R4) for x, respectively. Table 1. θmax with respect to the temporal order of accuracy (2l + 2) from 2 to 28 when m and ε are 50 000 and 0.01, respectively. 2l + 2  2  4  6  8  10  12  14  16  18  20  22  24  26  28  θmax  0.02  0.5  1.45  2.52  3.14  5.14  6.11  6.28  9.04  9.42  12  12.54  14.94  15.65  2l + 2  2  4  6  8  10  12  14  16  18  20  22  24  26  28  θmax  0.02  0.5  1.45  2.52  3.14  5.14  6.11  6.28  9.04  9.42  12  12.54  14.94  15.65  View Large Now we briefly discuss on how the proposed scheme is applied to real situations. If we assume that h = Δx = Δz, the maximum expressible k by the spatial discretization is $$\sqrt{2}\pi /h$$. Then, we can obtain Δt for given maximum velocity cmax, kmax and θmax as follows:   \begin{eqnarray} \Delta t = \frac{\theta _{\text{max}}h}{\sqrt{2}\pi c_{\text{max}}}. \end{eqnarray} (31) If we set the modelling time as tmax, the number of time steps Nt can be calculated as tmax/Δt. Then, Nop, the total number of operation of c2∇2, is obtained as (2l + 1)Nt using (11) and (15) as follows:   \begin{eqnarray} N_{\text{op}}=\frac{\sqrt{2}\pi c_{\text{max}} t_{\text{max}}}{\theta _{\text{max}} h}(2l+1). \end{eqnarray} (32) Since Nop is proportional to the computational cost of the modelling, we can examine the efficiency using (32) beforehand. Table 2 represents Δt, Nt and Nop with respect to the order of accuracy of the scheme to satisfy accuracy condition listed in Table 1 when cmax is 5 km s−1, h is 10 m and tmax is 10 s. The error of the numerical solution for θ less than θmax are bounded in 1 per cent since Nt is less than 50 000 for all cases. We can confirm that Nop tends to decrease as the higher-order method is used. Speed up factor for each Δt in Table 2 is the inversed value of normalized Nop by that of the fourth-order method, which shows that the simulation can be performed two or three times faster by higher-order methods than by fourth-order method. This tendency does not continue; the speed up factor appears to converge to a value between three and four in this case. If we assume that the maximum required frequency is given as 100 Hz, the optimal order of accuracy of the scheme is 18 or 20 because Δt should be less than 5 ms. Then we can perform the simulation almost three times faster than that by the fourth-order scheme with tiny time step length. In case that the target Δt is given, it is convenient to select the temporal order of accuracy of the scheme corresponding to θmax larger than but closest to $$\theta =\sqrt{2}\pi c_{\text{max}}\Delta t/h$$ because the larger the time steps are, the more efficient the simulation, according to the tendency of the accuracy analysis. For the example above, when Δt is given as 5 ms, we should select the order of accuracy as 22 because θ = 11.107 is less than and closest to θmax = 12. Then Nop of the case is 42 000, which corresponds to 3.17 if the value is converted to speed up factor with respect to Nop of the fourth-order scheme. Table 2. Parameters to estimate the efficiency of the modelling with respect to the order of accuracy from 2 to 28 satisfying the accuracy criterion listed in Table 1 when cmax is 5 km s−1, h is 10 m and tmax is 10 s. 2l + 2  4  6  8  10  12  14  16  18  20  22  24  26  28  Δt (ms)  0.23  0.65  1.13  1.41  2.31  2.75  2.83  4.07  4.24  5.40  5.65  6.73  7.04  Nt × 10−2  444  153  88  71  43  36  35  25  23  19  17  15  14  Nop × 10−3  1332  766  617  636  475  473  530  418  488  389  408  372  383  Speed up  1  1.74  2.16  2.09  2.80  2.82  2.51  3.19  2.98  3.43  3.27  3.59  3.48  2l + 2  4  6  8  10  12  14  16  18  20  22  24  26  28  Δt (ms)  0.23  0.65  1.13  1.41  2.31  2.75  2.83  4.07  4.24  5.40  5.65  6.73  7.04  Nt × 10−2  444  153  88  71  43  36  35  25  23  19  17  15  14  Nop × 10−3  1332  766  617  636  475  473  530  418  488  389  408  372  383  Speed up  1  1.74  2.16  2.09  2.80  2.82  2.51  3.19  2.98  3.43  3.27  3.59  3.48  View Large 4 NUMERICAL EXPERIMENTS In this section, numerical examples are presented to verify the properties of the proposed method explained in the analysis section. 4.1 Initial value problem 4.1.1 Homogeneous model First, the numerical modelling of the initial value problem in 2-D homogeneous media is performed. As shown in Fig. 3, the area of computational domain is 4π × 4π km2 which composed of 512 equally spaced discrete points with a grid interval of 0.0491 km for each spatial dimension and with the homogeneous P-wave speed c of 5 km s−1 and density ρ0 of 1 $$\text{g {}cm}^{-3}$$. The periodic boundary condition is naturally imposed because the pseudospectral method uses the harmonic basis functions to calculate the spatial derivatives. We use a Gaussian profile $$p_0=e^{-240\left((x-2\pi )^2+(z-2\pi )^2\right)}$$. For v at the initial step, we set it as zero for the simplicity. Four different Δt values, such as 1, 5, 10 and 20 ms, are used for the simulation over 50 s. Before performing the modelling, we can predict the performance as in the previous section. Table 3 represents the maximum expressible θ, the required order of accuracy of the scheme, Nop of the simulation under each condition and the relative speed up factor with respect to Nop of the simulation using 1 ms. As Δt increases, a higher-order scheme is required to accommodate the maximum expressible θ during the modelling. The larger the Δt, the smaller the Nop is needed; the computational cost for the numerical simulation with Δt = 20 ms is less than one-third of that with Δt = 1 ms as shown in Table 3. Wavefields at 1, 10, 20, 30, 40 and 50 s are illustrated in Fig. 4, with the results of different Δt values shown at each quadrant of the snapshots. All the wavefields are notably equivalent in terms of the scale and location of the events. Time traces of p at the receivers for each Δt from 49 to 50 s are shown by Fig. 5. We can confirm that all traces of four different Δt values are nearly identical without dispersion or dissipation phenomena. For quantitative comparison, the L2 errors of the numerical solution in Fig. 4 and the analytic solutions are obtained and represented in Table 4. The error values are negligibly small throughout the simulation for every Δt which, ascertain the feasibility of the strategy to select the order of accuracy. Figure 4. View largeDownload slide Combined wavefields obtained by numerical modelling at (a) 1, (b) 10, (c) 20, (d) 30, (e) 40 and (f) 50 s for each Δt at each quadrant. Figure 4. View largeDownload slide Combined wavefields obtained by numerical modelling at (a) 1, (b) 10, (c) 20, (d) 30, (e) 40 and (f) 50 s for each Δt at each quadrant. Figure 5. View largeDownload slide Plots of time traces recorded at (a) R1, (b) R2, (c) R3 and (d) R4 for each Δt values. Figure 5. View largeDownload slide Plots of time traces recorded at (a) R1, (b) R2, (c) R3 and (d) R4 for each Δt values. Table 3. Maximum expressible θ, the required order of accuracy, Nt and Nop with respect to Δt for the homogeneous velocity model with $$c=5{\, }\text{km s}^{-1}$$, h = 0.0491 km. Speed up factor of each case is the inversed value of normalized Nop by that of Δt = 1 ms. Δt (ms)  1  5  10  20  θ  0.4525  2.2627  4.5255  9.05  2l + 2  4  8  12  18  Nt  50 000  10 000  5000  2500  Nop  300 000  140 000  110 000  85 000  Speed up  1  2.143  2.727  3.529  Δt (ms)  1  5  10  20  θ  0.4525  2.2627  4.5255  9.05  2l + 2  4  8  12  18  Nt  50 000  10 000  5000  2500  Nop  300 000  140 000  110 000  85 000  Speed up  1  2.143  2.727  3.529  View Large Table 4. L2 errors between numerical and analytic wavefields at 1, 10, 20, 30, 40 and 50 s for each Δt. t (s)  1  10  20  30  40  50  L2 error (1 ms)  8.37E − 9  6.81E − 8  9.55E − 8  1.14E − 7  1.30E − 7  2.11E − 7  L2 error (5 ms)  5.72E − 10  2.19E − 9  9.09E − 10  2.36E − 8  2.24E − 8  3.91E − 8  L2 error (10 ms)  8.08E − 10  1.19E − 9  4.03E − 9  8.72E − 9  8.15E − 9  9.92E − 9  L2 error (20 ms)  3.57E − 10  6.71E − 10  7.13E − 10  2.28E − 9  3.60E − 9  4.34E − 9  t (s)  1  10  20  30  40  50  L2 error (1 ms)  8.37E − 9  6.81E − 8  9.55E − 8  1.14E − 7  1.30E − 7  2.11E − 7  L2 error (5 ms)  5.72E − 10  2.19E − 9  9.09E − 10  2.36E − 8  2.24E − 8  3.91E − 8  L2 error (10 ms)  8.08E − 10  1.19E − 9  4.03E − 9  8.72E − 9  8.15E − 9  9.92E − 9  L2 error (20 ms)  3.57E − 10  6.71E − 10  7.13E − 10  2.28E − 9  3.60E − 9  4.34E − 9  View Large 4.1.2 Heterogeneous model Now, we perform numerical experiments on whether the proposed modelling scheme can accurately consider the heterogeneity of the material models. The Marmousi-2 model (Martin et al. 2006) illustrated in Fig. 6 are used for the simulation. The computational domain is 17 km for x-axis and 3.5 km for z-axis with a grid size of 12.5 m. The maximum and minimum values of c are 4.7 and 1.028 km s−1, respectively. As with the homogeneous model example, periodic boundary condition is applied to all edges. The initial value for p is the Gaussian profile $$p_0=e^{-240\left((x-8.5)^2\,+\,(z-1.75)^2\right)}$$ and zero for the particle velocity. We perform the numerical modelling over 50 s using four different Δt values, 1, 5, 10, and 20 ms and record the time traces for each Δt at 0.875 km for z-axis and 8.5 (R1), 10.625 (R2), 12.75 (R3) and 14.875 km (R4) for x-axis. Figure 6. View largeDownload slide P-wave velocity (upper) and density (lower) of the Marmousi-2 model. Figure 6. View largeDownload slide P-wave velocity (upper) and density (lower) of the Marmousi-2 model. Table 5 represents the maximum expressible θ, the required order of accuracy of the scheme, Nop of the simulation under each condition and the relative speed up factor with respect to Nop of the simulation using 1 ms. Similar to the previous results using the homogeneous model, it can be seen that the Nop decreases as the Δt increases. The speed up factor shows that computational cost of the modelling using 5, 10, and 20 ms are less than half of that using 1 ms. Time traces for each Δt are compared in Fig. 7. In order to see the effect of accumulated errors over time, we compared the traces from 49 to 50 s. We can find that all traces of four different Δt values are nearly identical without dispersion or dissipation phenomena. The proposed method is shown to yield consistent and stable numerical results irrespective of time step lengths when using the heterogeneous model. Figure 7. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Figure 7. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Table 5. Maximum expressible θ, the required order of accuracy, Nt, Nop with respect to Δt for the simulation using the Marmousi-2 model with $$c_{\text{max}}=4.7{\, }\text{km s}^{-1}$$, h = 0.0125 km. Speed up factor of each case is inversed value of normalized Nop by that of Δt = 1 ms. Δt (ms)  1  5  10  20  θ  1.6705  8.3526  16.7052  33.4105  2l + 2  8  18  30  52  Nt  50 000  10 000  5000  2500  Nop  700 000  340 000  290 000  255 000  Speed up  1  2.059  2.414  2.745  Δt (ms)  1  5  10  20  θ  1.6705  8.3526  16.7052  33.4105  2l + 2  8  18  30  52  Nt  50 000  10 000  5000  2500  Nop  700 000  340 000  290 000  255 000  Speed up  1  2.059  2.414  2.745  View Large 4.2 Source problem with absorbing boundary condition In this section, we examine the feasibility of the proposed method in the context of source problems with the absorbing boundary condition for practical usages. Most conditions are identical to the previous example using the Marmousi-2 model. The Ricker wavelet $$q(t)=(1-2\pi ^2(f_0t-\beta )^2)e^{-\pi ^2(f_0t-\beta )^2}$$ is used as a source wavelet with f0 = 10 Hz and β = 1. The range of frequency band of the wavelet is 30 Hz with a dominant frequency component at 10 Hz. The wavelet is imposed to p with a distributed monopole source $$f(x,z)=e^{-500((x-x_s)^2+(z-z_s)^2)}$$ where xs = 8.5 km, zs = 0.0625 km. The following expression is used to consider the source term:   \begin{eqnarray} p(t+\Delta t) = p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + f(x,z)\frac{q(t) + q(t+\Delta t)}{2}\Delta t. \end{eqnarray} (33) Time traces are recorded by a line of receivers located at the depth of 0.0875 km. The numerical simulation with Δt = 20 ms is excluded because the Nyquist frequency of the numerical solution is 25 Hz which is lower than maximum frequency component of the source wavelet. Instead, we perform the simulation with Δt = 16 ms, whose Nyquist frequency is slightly less than 30 Hz. The required order of accuracy of the scheme and Nop are equal to those listed in Table 5 for Δt of 1, 5 and 10 ms. The required order of accuracy for Δt = 16 ms is 44 because the maximum expressible θ is 26.728. The parameter Nop for the case is 268 750 and speed up factor is 2.604. We apply the absorbing boundary condition suggested by Lavelle & Thacker (2008), a kind of sponge-type damping layer mimicking the perfectly matched layer (PML) because of the instability problem for large Δt. Although this method is not as efficient as PML, performs much better than the simple-sponge. From our experience, the simple-sponge requires a layer of 1.5 km thickness, whereas this method requires layer of 0.5 km to damp out the incident waves perfectly. In the damping area, the quadratic damping profile $$d_x(x)=\xi (x-x_d)^2/L_d^2$$ is used where xd is the starting position. In this study, we set ξ to 100, Ld to 0.5 km as stated above, and xd to 0 and 17 km for left and right boundaries, respectively. In the same way, we set quadratic damping profile in the z-direction, dz(z), with equivalent ξ and Ld where zd are 0 and 3.5 km for top and bottom layers. Fig. 8 shows that kinematics of recorded waves are in good agreement. It is also confirmed that the waves are effectively dampened by 0.5-km-thick absorbing layers. Time traces are plotted in Fig. 9 recorded from the seismogram at 8.5 (R1), 10.625 (R2), 12.75 (R3) and 14.875 km (R4) for x-axis and 0.0625 for z-axis. We can confirm that the arrival time and amplitude of the first arrival and subsequent signals are equivalent to each other regardless of Δt. Figure 8. View largeDownload slide Combination of seismograms for four different Δt values recorded at a line receiver located at 0.0625 km for z. Figure 8. View largeDownload slide Combination of seismograms for four different Δt values recorded at a line receiver located at 0.0625 km for z. Figure 9. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Figure 9. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. 5 CONCLUSIONS In this study, we suggested the staggered-time method whose order of accuracy in temporal dimension can be extended arbitrarily using the simple recurrence relations. The method is applied to the linear acoustic wave system of pressure and particle velocity in a stagnant and stationary media. We analysed the method quantifying the error of the numerical solution and proposed the error criterion that guarantees 1 per cent error during 50 000 time steps. We also propose a simple strategy how to use the result of the analysis; rather than adjusting the modelling conditions to satisfy the stability and accuracy of a given modelling scheme, we set temporal order of accuracy of the scheme satisfying the error condition for given time step length, grid intervals and velocity model. Such property allows us to flexibly determine modelling conditions. For instance, when the scheme is used in the reverse time migration algorithm, there is no need to reduce the time step length more than necessary to match the desired spatial resolution because of the stability issue. The proposed modelling scheme allows not only the flexibility, but also the efficiency; the calculation efficiency is improved significantly with larger time step length and higher order of accuracy, which is proven by the analysis and numerical examples. By virtue of this property, it is efficient for the full waveform inversion algorithms when the method is applied. Because we can accurately perform the simulation with a sufficiently large time step length corresponding to the low frequencies generally used during the inversion process. We also confirmed by the numerical examples that the method perfectly considers the heterogeneity of the media. We proposed how to properly impose the source wavelet, which yield equivalent solutions regardless of the time step lengths whose Nyquist frequency is less than that of the maximum frequency of the wavelet. As stated above, we expect that the proposed method is effective for the geophysical applications. This method can also be applied to simulate other wave phenomena such as electromagnetic or elastodynamic waves which also have many practical applications, that is, ground penetrating radar or medical imaging algorithms. ACKNOWLEDGEMENTS This work was supported by the Energy Efficiency and Resources Core Technology Program of the Korea Institute of Energy Technology Evaluation and Planning (KETEP), granted financial resource from the Ministry of Trade, Industry and Energy, Republic of Korea (20132510100060). And the Institute of Engineering Research at Seoul National University provided research facilities for this work. REFERENCES Arnold D.N., Brezzi F., Cockburn B., Marini L.D., 2002. Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. , 39( 5), 1749– 1779. Google Scholar CrossRef Search ADS   Ashcroft G., Zhang X., 2003. Optimized prefactored compact schemes, J. Comput. Phys. , 190( 2), 459– 477. Google Scholar CrossRef Search ADS   Baysal E., Kosloff D.D., Sherwood J.W., 1983. Reverse time migration, Geophysics , 48( 11), 1514– 1524. Google Scholar CrossRef Search ADS   Bogey C., Bailly C., 2004. A family of low dispersive and low dissipative explicit schemes for flow and noise computations, J. Comput. Phys. , 194( 1), 194– 214. Google Scholar CrossRef Search ADS   Carcione J., 1996. A 2D Chebyshev differential operator for the elastic wave equation, Comput. Methods Appl. Mech. Eng. , 130( 1–2), 33– 45. Google Scholar CrossRef Search ADS   Chen J.-B., 2009. Lax–Wendroff and Nyström methods for seismic modelling, Geophys. Prospect. , 57( 6), 931– 941. Google Scholar CrossRef Search ADS   Cohen G., 1986. Fourth-order schemes for the 2-d wave equation in a homogeneous medium, in SEG Technical Program Expanded Abstracts 1986 , pp. 632– 635, Society of Exploration Geophysicists. Dablain M., 1986. The application of high-order differencing to the scalar wave equation, Geophysics , 51( 1), 54– 66. Google Scholar CrossRef Search ADS   De Basabe J.D., Sen M.K., Wheeler M.F., 2008. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion, Geophys. J. Int. , 175( 1), 83– 93. Google Scholar CrossRef Search ADS   Dumbser M., Käser M., 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—II. The three-dimensional isotropic case, Geophys. J. Int. , 167( 1), 319– 336. Google Scholar CrossRef Search ADS   Hall B., 2015. Lie Groups, Lie Algebras, and Representations: An Elementary Introduction , vol. 222, Springer. Google Scholar CrossRef Search ADS   Hermann V., Käser M., Castro C.E., 2011. Non-conforming hybrid meshes for efficient 2-D wave propagation using the discontinuous Galerkin method, Geophys. J. Int. , 184( 2), 746– 758. Google Scholar CrossRef Search ADS   Hesthaven J.S., Warburton T., 2007. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications , Springer Science & Business Media. Käser M., Dumbser M., 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—I. The two-dimensional isotropic case with external source terms, Geophys. J. Int. , 166( 2), 855– 877. Google Scholar CrossRef Search ADS   Klin P., Priolo E., Seriani G., 2010. Numerical simulation of seismic wave propagation in realistic 3-D geo-models with a Fourier pseudo-spectral method, Geophys. J. Int. , 183( 2), 905– 922. Google Scholar CrossRef Search ADS   Komatitsch D., Vilotte J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Kosloff D.D., Baysal E., 1982. Forward modeling by a Fourier method, Geophysics , 47( 10), 1402– 1412. Google Scholar CrossRef Search ADS   Kosloff D., Kosloff R., 1983. A Fourier method solution for the time dependent Schrödinger equation as a tool in molecular dynamics, J. Comput. Phys. , 52( 1), 35– 53. Google Scholar CrossRef Search ADS   Kosloff D., Queiroz Filho A., Tessmer E., Behle A., 1989. Numerical solution of the acoustic and elastic wave equations by a new rapid expansion method, Geophys. Prospect. , 37( 4), 383– 394. Google Scholar CrossRef Search ADS   Lavelle J., Thacker W., 2008. A pretty good sponge: dealing with open boundaries in limited-area ocean models, Ocean Modelling , 20( 3), 270– 292. Google Scholar CrossRef Search ADS   Lele S.K., 1992. Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. , 103( 1), 16– 42. Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P –SV seismograms, Geophysics, 53( 11), 1425– 1436. Liu Y., Sen M.K., 2009. A new time–space domain high-order finite-difference method for the acoustic wave equation, J. Comput. Phys. , 228( 23), 8779– 8806. Google Scholar CrossRef Search ADS   Long G., Zhao Y., Zou J., 2013. A temporal fourth-order scheme for the first-order acoustic wave equations, Geophys. J. Int. , 194( 3), 1473– 1485. Google Scholar CrossRef Search ADS   Martin G.S., Wiley R., Marfurt K.J., 2006. Marmousi2: an elastic upgrade for Marmousi, Leading Edge , 25( 2), 156– 166. Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Gális M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Pestana R.C., Stoffa P.L., 2010. Time evolution of the wave equation using rapid expansion method, Geophysics , 75( 4), T121– T131. Google Scholar CrossRef Search ADS   Pratt R.G., Shin C., Hick G., 1998. Gauss–Newton and full newton methods in frequency–space seismic waveform inversion, Geophys. J. Int. , 133( 2), 341– 362. Google Scholar CrossRef Search ADS   Renaut R., Fröhlich J., 1996. A pseudospectral Chebychev method for the 2D wave equation with domain stretching and absorbing boundary conditions, J. Comput. Phys. , 124( 2), 324– 336. Google Scholar CrossRef Search ADS   Shin C., Cha Y.H., 2008. Waveform inversion in the Laplace domain, Geophys. J. Int. , 173( 3), 922– 931. Google Scholar CrossRef Search ADS   Shin C., Cha Y.H., 2009. Waveform inversion in the Laplace–Fourier domain, Geophys. J. Int. , 177( 3), 1067– 1079. Google Scholar CrossRef Search ADS   Strikwerda J.C., 2004. Finite Difference Schemes and Partial Differential Equations , SIAM. Google Scholar CrossRef Search ADS   Tal-Ezer H., 1986. Spectral methods in time for hyperbolic equations, SIAM J. Numer. Anal. , 23( 1), 11– 26. Google Scholar CrossRef Search ADS   Tam C.K., Webb J.C., 1993. Dispersion-relation-preserving finite difference schemes for computational acoustics, J. Comput. Phys. , 107( 2), 262– 281. Google Scholar CrossRef Search ADS   Tan S., Huang L., 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation, Geophys. J. Int. , 197( 2), 1250– 1267. Google Scholar CrossRef Search ADS   Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics , 49( 8), 1259– 1266. Google Scholar CrossRef Search ADS   Tessmer E., 2011. Using the rapid expansion method for accurate time-stepping in modeling and reverse-time migration, Geophysics , 76( 4), S177– S185. Google Scholar CrossRef Search ADS   Trefethen L.N., 2000. Spectral Methods in MATLAB , SIAM. Google Scholar CrossRef Search ADS   Virieux J., 1984. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 49( 11), 1933– 1942. Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. Google Scholar CrossRef Search ADS   Whitmore N., 1983. Iterative depth migration by backward time propagation, in SEG Technical Program Expanded Abstracts 1983 , pp. 382– 385, Society of Exploration Geophysicists. Yee K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. , 14( 3), 302– 307. Google Scholar CrossRef Search ADS   APPENDIX A: ABSORBING BOUNDARY CONDITION We implemented the pretty-good-sponge boundary condition (Lavelle & Thacker 2008) because the instability occurs in the damping layer when applying the PML. Although less effective than the PML, it has comparable damping ability with moderate thickness of damping layer. The 2-D acoustic wave system defined in the damping layer is written as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ \mathbf {v} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}-{\rm tr}(\mathbf {D}) & -\rho _0 c^2 \nabla \cdot \\ -{\rho _0}^{-1}\nabla & -\mathbf {D} \end{array}\right]}}{\left[\begin{array}{c}p \\ \mathbf {v} \end{array}\right]}, \end{eqnarray} (A1)where D is a diagonal 2 × 2 matrix composed of damping parameters, that is, D11 = dx(x) and D22 = dz(z) for 2-D problem where dx and dz are the damping profiles for x- and z-directions, respectively. The proposed method defined in the damping layer can be derived by applying (A1) to (4) and (5) omitting the cross-terms between damping parameter and the spatial derivatives as   \begin{eqnarray} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) {}& - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \Big ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \Big )^{n} p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}\mathbf {v}(t), \end{eqnarray} (A2)and   \begin{eqnarray} p(t+\Delta t) \approx p(t) {}& - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1}p\bigg (t+\frac{\Delta t}{2}\bigg ) . \end{eqnarray} (A3) If we set v(t) and p(t + Δt/2) as an averaged value as   $$\mathbf {v}(t) \approx \frac{1}{2}\bigg ( \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) \bigg ),\quad p\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \frac{1}{2}\big ( p(t+\Delta t) + p(t) \big ),$$ (A4)we can rewrite eqs (A2) and (A3) as   \begin{eqnarray} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \mathbf {B}^{-1} \left( \mathbf {A}{\, }\mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t) \right) + O(\Delta t^2), \end{eqnarray} (A5) and   \begin{eqnarray} p(t+\Delta t) \approx \beta ^{-1} \left( \alpha {\, }p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \right) + O(\Delta t^2), \end{eqnarray} (A6)where   \begin{eqnarray} \mathbf {A} &=& I - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}, \nonumber \\ \mathbf {B} &=& I + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}, \nonumber \\ \alpha &=& 1 - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1},\nonumber \\ \beta &=& 1 + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1}. \end{eqnarray} (A7) The temporal order in the damping layer is only two because of the approximation and the time averaging of the variables, which may aggravate the accuracy of the solution in the whole domain if reflections are unignorable. If there is almost no reflection, the wavefield inside the computational domain maintains the desired order of accuracy in time. APPENDIX B: ANALYTIC SOLUTION We can derive the analytic solution of the acoustic wave equation in a homogeneous model using the one-step method, which imposes the periodic boundary condition by default if the system is solved in the wavenumber domain by using the Fourier basis (Trefethen 2000). The linear acoustic wave system (1) defined in the 2-D domain can be rewritten as follows:   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ u \\ w \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c@{\quad}c}0 & -\rho _0 c^2 \frac{\partial }{\partial x} & -\rho _0 c^2\frac{\partial }{\partial z} \\ -\frac{1}{\rho _0}\frac{\partial }{\partial x} & 0 & 0 \\ -\frac{1}{\rho _0}\frac{\partial }{\partial z} & 0 & 0 \end{array}\right]}}{\left[\begin{array}{c}p \\ u \\ w \end{array}\right]}. \end{eqnarray} (B1) This equation can be converted into system ordinary differential equation form via spatial Fourier transformation of the system variables as follows:   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}\widetilde{p} \\ \widetilde{u} \\ \widetilde{w} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c@{\quad}c}0 & i\rho _0 c^2 k_x & i\rho _0 c^2 k_z \\ \frac{i}{\rho _0} k_x & 0 & 0 \\ \frac{i}{\rho _0} k_z & 0 & 0 \end{array}\right]}}{\left[\begin{array}{c}\widetilde{p} \\ \widetilde{u} \\ \widetilde{w} \end{array}\right]}. \end{eqnarray} (B2) Then, the acoustic wave solution in the time-wavenumber domain is expressed as   \begin{eqnarray} \widetilde{y}(t)=e^{t\widetilde{\mathbf {M}}}\widetilde{y}_0, \end{eqnarray} (B3)where $$\widetilde{y}_0$$ and $$\widetilde{y}(t)$$ is the transformed system variable vector at the initial and t, respectively, and $$\widetilde{\mathbf {M}}$$ is the characteristic matrix of (B2) whose exponential map can be formulated as follows (Hall 2015):   \begin{eqnarray} e^{t\widetilde{\mathbf {M}}}= {{\left[\begin{array}{c@{\quad}c@{\quad}c}\cos (ckt) & \frac{ick_x\sin (ckt)}{k} & \frac{ick_z\sin (ckt)}{k} \\ \frac{ick_x\sin (ckt)}{k} & \frac{k_z^2+k_x^2\cos (ckt)}{k^2} & \frac{k_xk_z\left(\cos (ckt)-1\right)}{k^2} \\ \frac{ick_z\sin (ckt)}{k} & \frac{k_xk_z\left(\cos (ckt)-1\right)}{k^2} & \frac{k_z^2+k_x^2\cos (ckt)}{k^2} \end{array}\right]}}. \end{eqnarray} (B4) Therefore, time-space domain solution can be easily obtained by applying the inverse Fourier transform for both sides of equation (B3) as   \begin{eqnarray} y(t) = \mathcal {F}_{xz}^{-1}\left(e^{t\widetilde{\mathbf {M}}} \mathcal {F}_{xz}\left(y_0\right)\right), \end{eqnarray} (B5)where $$\mathcal {F}_{xz}$$ is the 2-D spatial Fourier transform operator for the xz-plane. APPENDIX C: GENERALIZATION The arbitrary-order staggered time integrator can be applied to other linear wave equations, that is, first-order elastodynamic or electromagnetic wave system. If we define system variables as $$\mathbf {p}\in \mathbb {R}^m$$ and $$\mathbf {q}\in \mathbb {R}^n$$, we can express the first-order wave equations in a general form as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}\mathbf {p} \\ \mathbf {q} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}0 & \mathbf {M_1} \\ \mathbf {M_2} & 0 \end{array}\right]}}{\left[\begin{array}{c}\mathbf {p} \\ \mathbf {q} \end{array}\right]}, \end{eqnarray} (C1)where M1 and M2 are the submatrices of the characteristic matrix M. M1 and M2 have dimensions of m × n and n × m, respectively. As can be found in (1), components of these matrices consist of products of spatial differential operators and material properties. For the 2-D linear elastic wave equation in the isotropic media (Virieux 1986), we can set p and q as the stress σ = [σxx, σzz, σxz]⊤ and the particle velocity v = [u, w]⊤, respectively. M1 and M2 of the system can be written as   \begin{eqnarray} \mathbf {M}_1 = {\left[\begin{array}{c@{\quad}c}(\lambda + 2\mu )\frac{\partial }{\partial x} & \mu \frac{\partial }{\partial z} \\ \mu \frac{\partial }{\partial x} & (\lambda + 2\mu )\frac{\partial }{\partial z} \\ \mu \frac{\partial }{\partial z} & \mu \frac{\partial }{\partial x} \end{array}\right]},\quad \mathbf {M}_2 = \frac{1}{\rho _0}{\left[\begin{array}{c@{\quad}c@{\quad}c}\frac{\partial }{\partial x} & 0 & \frac{\partial }{\partial z} \\ 0 & \frac{\partial }{\partial z} & \frac{\partial }{\partial x} \end{array}\right]}, \end{eqnarray} (C2)where λ and μ are Láme parameters. For Maxwell’s equations in the isotropic media (Yee 1966), we can set p and q as the electric field E = [Ex, Ey, Ez]⊤ and the magnetic field H = [Hx, Hy, Hz]⊤, respectively. M1 and M2 of the system can be written as   \begin{eqnarray} \mathbf {M}_1 = \frac{1}{\varepsilon } {\left[\begin{array}{c@{\quad}c@{\quad}c}0 & -\frac{\partial }{\partial z} & \frac{\partial }{\partial y} \\ \frac{\partial }{\partial z} & 0 & -\frac{\partial }{\partial x} \\ -\frac{\partial }{\partial y} & \frac{\partial }{\partial x} & 0 \end{array}\right]},\quad \mathbf {M}_2 = \frac{1}{\mu } {\left[\begin{array}{c@{\quad}c@{\quad}c}0 & \frac{\partial }{\partial z} & -\frac{\partial }{\partial y} \\ -\frac{\partial }{\partial z} & 0 & \frac{\partial }{\partial x} \\ \frac{\partial }{\partial y} & -\frac{\partial }{\partial x} & 0 \end{array}\right]}, \end{eqnarray} (C3)where ε and μ are the permittivity and the permeability of a medium, respectively. Following the process from (2) to (5), we can apply the arbitrary-order staggered time integrator to system (C1). If we impose p on the integer and q on the half integer of the temporal grid, we can obtain the semi-discretized formulation of (C1) as follows:   \begin{eqnarray} &&{\displaystyle \mathbf {q}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {q}\bigg (t-\frac{\Delta t}{2}\bigg ) + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \mathbf {M_2} \Big ( \mathbf {M_1}\mathbf {M_2} \Big )^{n} \mathbf {p}(t) + O(\Delta t^{2l+2}),} \nonumber\\ &&{\mathbf {p}(t+\Delta t) = \mathbf {p}(t) + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \mathbf {M_1}\Big ( \mathbf {M_2}\mathbf {M_1}\Big )^{n} \displaystyle \mathbf {q}\bigg (t+\frac{\Delta t}{2}\bigg ) + O(\Delta t^{2l+2}).} \end{eqnarray} (C4) © 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

# An arbitrary-order staggered time integrator for the linear acoustic wave equation

, Volume 212 (2) – Feb 1, 2018
15 pages

/lp/ou_press/an-arbitrary-order-staggered-time-integrator-for-the-linear-acoustic-DK557eAfhE
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx457
Publisher site
See Article on Publisher Site

### Abstract

Summary We suggest a staggered time integrator whose order of accuracy can arbitrarily be extended to solve the linear acoustic wave equation. A strategy to select the appropriate order of accuracy is also proposed based on the error analysis that quantitatively predicts the truncation error of the numerical solution. This strategy not only reduces the computational cost several times, but also allows us to flexibly set the modelling parameters such as the time step length, grid interval and P-wave speed. It is demonstrated that the proposed method can almost eliminate temporal dispersive errors during long term simulations regardless of the heterogeneity of the media and time step lengths. The method can also be successfully applied to the source problem with an absorbing boundary condition, which is frequently encountered in the practical usage for the imaging algorithms or the inverse problems. Numerical approximations and analysis, Computational seismology, Wave propagation 1 INTRODUCTION Numerical simulation of wave propagation is not only meaningful for forward problem itself but also plays an important role in several applications, such as reverse time migration for subsurface imaging (e.g. Baysal et al. 1983; Whitmore 1983) and full waveform inversion for reconstruction of the subsurface parameters (e.g. Tarantola 1984; Pratt et al. 1998; Shin & Cha 2008, 2009). The more accurate the modelling scheme, the higher the resolution of the obtained image. However, since discretization of wave equations inevitably leads to numerical errors, that is, phase and amplitude errors, various methodologies have been developed to improve the accuracy of numerical modelling of wave phenomena. In the geophysics literature, Virieux (1984) introduced the staggered-grid finite-difference method (FDM) by adopting the method of Yee (1966) to simulate the seismic wave equation represented by a first-order system of the stress and the particle velocity fields. This method defines the wavefields and material properties in the staggered grids separately and updates the stress and the particle velocity fields alternatively. Such gridding suppresses the even-odd decoupling and keeps the numerical solution from high-frequency oscillation. Increasing the spatial order of accuracy have been usually attempted to reduce the spatial dispersion error in the context of FDMs (e.g. Levander 1988). Such methods enlarge the grid stencil and incur additional treatments to impose boundary conditions. The compact-stencil FDM (Lele 1992) can also be used to efficiently yield an accurate solution. The method is a generalization of Padé scheme implicitly calculating the high-order spatial derivatives by solving tridiagonal or pentadiagonal matrix (three-point or five-point stencil), which can simply be solved with the complexity of O(N). This approach can yield more accurate solutions with a smaller finite difference stencil than that of the standard FDM. By optimizing the finite difference coefficients to conform the physical property of the wave phenomena, the dispersive and dissipative errors can be minimized (Ashcroft & Zhang 2003). Similar dispersion-relation-preserving (DRP) techniques for the explicit FDM have been also introduced using standard stencil (e.g. Tam & Webb 1993; Bogey & Bailly 2004; Liu & Sen 2009). The DRP techniques significantly reduce the dispersion error of the solution compared to the standard FDM with the equivalent stencil and the order of accuracy. For a deeper exposition on seismic modelling with FDM, look for Moczo et al. (2014). Schemes based on the finite-element method (FEM) increase the accuracy of the scheme by using a number of high-order basis polynomials locally defined in an element to approximate a solution. In the field of seismology, the spectral element method (SEM), a type of FEM technique, is used for the global seismic modelling (Komatitsch & Vilotte 1998). This method usually uses the high-order Gauss–Lobatto collocation nodes to define nodal (or modal) basis functions, which generate a diagonal global mass matrix naturally, allowing effective massive time-domain modelling. De Basabe et al. (2008) compared the dispersion characteristics of the discontinuous Galerkin (DG) method with that of the SEM. DG schemes are non-conforming methods that evaluate the wave solution of each element separately and compensate the discontinuity of the solution with numerical flux terms, such as Godnov or Lax–Friedrich flux (Arnold et al. 2002; Hesthaven & Warburton 2007). The methods can be used for special purposes when conducting simulation that contains complex geometry, that is, models with complex topography or interfaces of medium of different phases (Hermann et al. 2011), which is almost impossible to solve with conventional FDMs. Pseudospectral methods usually offer spectral accuracy in the spatial domain, which implies that the scheme is dispersion-free within the expressible band of wavenumber (Trefethen 2000). A type of basis function for the method depends on the boundary conditions of the problem. The Fourier basis functions are generally used which intrinsically assume the periodic boundary condition for each spatial dimension (Kosloff & Baysal 1982; Kosloff & Kosloff 1983). In order to impose general boundaries, such as Dirichlet or Neumann boundaries, the Chebyshev polynomial basis are required (Carcione 1996) and research has been conducted to make the stability condition comparable to that of the Fourier pseudospectral method (Renaut & Fröhlich 1996). Recently, Klin et al. (2010) suggested simple but effective method to impose a free surface boundary condition when using the Fourier pseudospectral method, which was verified by simulation of 3-D realistic models. Either Fourier or Chebyshev pseudospectral method incurs a pair of forward and inverse Fourier transforms to calculate spatial derivatives implicitly likewise the compact-stencil FDM approach. This can be accelerated by using the fast Fourier transform algorithm, which has the complexity of O(N log N). In general, the higher the spatial accuracy of a numerical method, the more strict the stability condition. The stability condition can be relaxed by matching the temporal order of accuracy with the spatial one. The Lax–Wendroff method can increase the temporal order of accuracy by converting the higher order time derivatives obtained from the Taylor series to the spatial derivatives. Temporal fourth-order scheme (Cohen 1986; Dablain 1986) was implemented to solve the acoustic wave equation by the Lax–Wendroff method, which has also been applied to derivation of ADER-DG schemes (Dumbser & Käser 2006; Käser & Dumbser 2006). Tal-Ezer (1986) suggested another approach to increase the temporal accuracy. This method expands the exponential map of the characteristic matrix of the hyperbolic system to the series of product of the Bessel function and the Chebyshev polynomial by using the Jacobi-Anger expansion. In the context of geophysics, the methods using the Jacobi-Anger expansion as the time stepping method are referred to as the rapid expansion method. This method was originally used to solve a single source problem (Kosloff et al. 1989) and then was extended to multisource problems for the practical applications (Pestana & Stoffa 2010; Tessmer 2011). In this study, we suggest a staggered time integrator which can extend the order of accuracy arbitrarily using the Lax–Wendroff method to solve the linear acoustic wave system. This method incorporates several low-order schemes such as Long et al. (2013) and Tan & Huang (2014). The proposed numerical method can be implemented easily by using a recurrence relation which enables the numerical solution to be almost non-dispersive in temporal dimension. In addition to the temporal accuracy, spectral accuracy in the spatial domain can also be achieved by spatial differentiations in the wavenumber domain via spatial Fourier transform as in the pseudospectral method. We faithfully carry the derivation process and implementation of the method in Theory part. The proposed method is analysed in the consecutive section. Simple strategy on selecting a proper temporal order of accuracy is also stated based on the qualitative expression of the error between the numerical and the exact wave solution of the initial value problem. Then, we verify the proposed modelling method following the strategy with the homogeneous P-wave velocity model. By numerical simulations, we also examine the applicability for the practical usage such as a problem with heterogeneity or a source problem with the absorbing boundary condition. 2 THEORY The acoustic wave equation in the first-order system consists of the linearized conservation relations with respect to the pressure perturbation p and the particle velocity vector v as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ \mathbf {v} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}0 & -\rho _0 c^2 \nabla \cdot \\ -{\rho _0}^{-1}\nabla & 0 \end{array}\right]}}{\left[\begin{array}{c}p \\ \mathbf {v} \end{array}\right]}, \end{eqnarray} (1)where ρ0 and c are the density and the P-wave propagation speed, respectively. Temporal discretization of such first-order system is usually done with the leap-frog or the two-stage Störmer–Verlet method [Maxwell’s equations (Yee 1966), the elastic wave equation (Virieux 1986), etc.], whose structural preserving property bounds the error of Hamiltonian in some level (Chen 2009). In this study, we impose each physical value on the staggered time axis; p on the integer and v on the half integer of the temporal grid. Consider v(t + Δt/2) with respect to v(t) using the Taylor expansion as   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}(t) + \frac{\Delta t}{2}\frac{\partial \mathbf {v}}{\partial t}\bigg \vert _{t} + \frac{\Delta t^2}{8}\frac{\partial ^2\mathbf {v}}{\partial t^2}\bigg \vert _{t} + \cdots + \frac{\Delta t^n}{n! \, 2^n}\frac{\partial ^n\mathbf {v}}{\partial t^n}\bigg \vert _{t} + \cdots , \end{eqnarray} (2)and v(t − Δt/2) as   \begin{eqnarray} \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) = \mathbf {v}(t) - \frac{\Delta t}{2}\frac{\partial \mathbf {v}}{\partial t}\bigg \vert _{t} + \frac{\Delta t^2}{8}\frac{\partial ^2\mathbf {v}}{\partial t^2}\bigg \vert _{t} + \cdots + \frac{(-\Delta t)^n}{n! \, 2^n}\frac{\partial ^n\mathbf {v}}{\partial t^n}\bigg \vert _{t} + \cdots . \end{eqnarray} (3) We can obtain the equation of v(t + Δt/2) by subtracting (3) from (2) as follows:   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) + \sum \limits _{n=0}^\infty \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\frac{\partial ^{(2n+1)}\mathbf {v}}{\partial t^{(2n+1)}}\bigg \vert _{t}{\, }. \end{eqnarray} (4) In a similar manner, we can derive the equation of p(t + Δt) with respect to p(t) and $$\partial _t^{(2n+1)}p(t+dt/2)$$ as follows:   \begin{eqnarray} p(t+\Delta t) = p(t) + \sum \limits _{n=0}^\infty \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\frac{\partial ^{(2n+1)} p}{\partial t^{(2n+1)}}\bigg \vert _{t+{\Delta t}/{2}}{\, }. \end{eqnarray} (5) Following the Lax–Wendroff method, the series of time derivatives of (4) and (5) can be replaced by the spatial derivatives using (1). Then, the time stepping procedure can be expressed as follows with an arbitrary order of 2l + 2 as   \begin{eqnarray} \begin{array}{lcl} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t) + O(\Delta t^{2l+2}), \\ p(t+\Delta t) = p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + O(\Delta t^{2l+2}){\, }. \end{array} \end{eqnarray} (6) It is noted that (6) is the extended formulation of the leap-frog or the two-stage Störmer–Verlet method. If l is zero, the method is equivalent to the conventional second-order staggered-time method (Yee 1966; Virieux 1986). The fourth-order method of Long et al. (2013) (l = 1) and the sixth-order method of Tan & Huang (2014) (l = 2) are subsets of (6). The structural preserving property is also preserved for (6) that no dissipative error occurs during the simulation, intrinsically. Only dispersive error may debase the numerical solution of (6), which can effectively be reduced by increasing the temporal order of accuracy of the scheme. The spatial derivative operators can be calculated via several methods. In this study, we adopt the Fourier pseudospectral method which gives the almost exact spatial derivatives for sufficiently smooth wavefields (Trefethen 2000) in order to correspond to the temporal accuracy of (6). The x-directional spatial derivative of f via the Fourier pseudospectral method is implemented by using the discrete Fourier transform as   \begin{eqnarray} \frac{\partial f}{\partial x} = \mathcal {F}_x^{-1} ( -ik_x{\, }{\, } \mathcal {F}_x (f) ), \end{eqnarray} (7)where $$\mathcal {F}_x$$ and kx are the Fourier transform operator and the wavenumber in the x-direction. The proposed method (6) of desired order of accuracy can be effectively calculated by using a simple recurrence relation. If we define scalar fields Al and an as   \begin{eqnarray} A_l = \sum \limits _{n=0}^l a_n = \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t), \end{eqnarray} (8)the recurrence relation can be simply expressed as   \begin{eqnarray} a_n = \frac{\Delta t^2}{8n{\, }(2n+1)}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )a_{n-1}, \quad n=1,2, \ldots ,l \end{eqnarray} (9)where   \begin{eqnarray} a_0 = \Delta t {\, } p(t){\, }. \end{eqnarray} (10) Then v(t + Δt/2) is calculated as follows:   \begin{eqnarray} \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \frac{1}{\rho _0}\nabla A_l . \end{eqnarray} (11)In the same way, we can define scalar fields Bl and bn as follows:   \begin{eqnarray} B_l = \sum \limits _{n=0}^l b_n = \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ). \end{eqnarray} (12)The recurrence relation for bn is equivalent to (9) as   \begin{eqnarray} b_n = \frac{\Delta t^2}{8n{\, }(2n+1)}\left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )b_{n-1}, \quad n=1,2, \ldots ,l , \end{eqnarray} (13)where   \begin{eqnarray} b_0 = {\Delta t}{\, }{\rho _0 c^2}\nabla \cdot \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ). \end{eqnarray} (14) Then p(t + Δt) can be calculated as follows:   \begin{eqnarray} p(t+{\Delta t}) = p(t) - B_l . \end{eqnarray} (15) We can increase the temporal order of accuracy of (6) by two by adding the subsequent term of the series of the spatial derivatives by using (9) and (13), which only requires to apply the second-order spatial derivative operator ρ0c2∇ · ρ0−1∇ to update the sequence. The operator appears in the scalar wave equation considering the heterogeneity of ρ0 in the domain. The operator is applied to a scalar field f via the Fourier pseudospectral method in 2-D domain as follows:   \begin{eqnarray} {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla f = {\rho _0 c^2} \Big ( \mathcal {F}_x^{-1} \Big ( -ik_x{\, } \mathcal {F}_x \Big ( \frac{1}{\rho _0} \mathcal {F}_x^{-1}\Big ({-}ik_x{\, } \mathcal {F}_x (f) \Big ) \Big ) \Big ) + \mathcal {F}_z^{-1} \Big ({-}ik_z{\, } \mathcal {F}_z \Big ( \frac{1}{\rho _0} \mathcal {F}_z^{-1}\Big (-ik_z{\, } \mathcal {F}_z (f) \Big ) \Big ) \Big ) \Big ). \end{eqnarray} (16) To increase the index of an or bn by one, two pairs of Fourier and inverse Fourier transform are required for each direction to consider the heterogeneity of ρ0. If ρ0 is constant elsewhere, the operator is simplified to the Laplacian, that is, c2∇2, which is calculated by the Fourier pseudospectral method as   \begin{eqnarray} {c^2}\nabla ^2 f = {c^2} \left ( \mathcal {F}_x^{-1} \left ( -k_x^2{\, } \mathcal {F}_x (f) \right ) + \mathcal {F}_z^{-1} \left ( -k_z^2{\, } \mathcal {F}_z (f) \right ) \right ). \end{eqnarray} (17) In this case, the required number of pair of the Fourier and inverse Fourier transform is reduced by half compared to (16). 3 ANALYSIS In this section, we manipulate (6) to demonstrate how we can control the accuracy of the numerical simulation result under given conditions. This section is related to both dispersion error and the stability analysis of the scheme. If we rewrite the scheme (6) with respect to p, we can obtain the following relation.   \begin{eqnarray} p(t+\Delta t) + p(t-\Delta t) = \bigg ( 2 + {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \bigg ( \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\bigg ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \bigg )^{n} \bigg )^2 \bigg ) p(t). \end{eqnarray} (18) For the analysis, we consider the initial value problem with p(0) = p0 and v(0) = 0. As the first step of the time marching, the appropriate calculation of v(Δt/2) is needed. The exact value of v(Δt/2) cannot be obtained by using (6) but directly discretizing (1) using the Lax–Wendroff method as   \begin{eqnarray} \mathbf {y}\Big (\frac{\Delta t}{2}\Big ) = \mathbf {y}(0) + \sum \limits _{n=1}^{2l+2} \frac{\Delta t^n\mathbf {M}^n}{n!{}\, 2^n}\frac{\partial ^{n}\mathbf {y}}{\partial t^n}\bigg \vert _{t=0} + O(\Delta t^{2l+3}), \end{eqnarray} (19)where y(t) = [p(t), v(t) ]⊤ and M is the characteristic matrix of (1). Since v(0) = 0, v(Δt/2) is expressed using p0 with an equivalent temporal order of accuracy as follows:   \begin{eqnarray} \mathbf {v}\left ( \frac{\Delta t}{2} \right ) = \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p_0 + O(\Delta t^{2l+2}). \end{eqnarray} (20)Then we can calculate p(Δt) using (6) and (20) as follows:   \begin{eqnarray} \displaystyle p\left(\Delta t \right) = p_0 + 2{\, }{\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \left( \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \left( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n} \right)^2 p_0. \end{eqnarray} (21)As a result, we now have the relation of p with respect to the number of time steps m as   \begin{eqnarray} \begin{array}{lcl}p\left( 0 \right) = p_0, \\ \displaystyle p\left(\Delta t \right) = {\rm C}(l){\, }p_0, \\ \displaystyle p(m\Delta t + \Delta t) = 2 {\rm C}(l){\, }p(m \Delta t) - p(m\Delta t - \Delta t), \quad m=1, 2, \ldots , \end{array} \end{eqnarray} (22)where   \begin{eqnarray} {\rm C}(l) = 1 + 2{\, }{\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \left( \sum \limits _{n=0}^{l} \frac{\Delta t^{2n+1}}{(2n+1)!{\, }2^{2n+1}} \left( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right)^{n} \right)^2. \end{eqnarray} (23) Relation (22) is equivalent to the three-term recurrence formulation of the Chebyshev polynomials of the first kind, denoted Tn(x) with the degree of the polynomials of n. Therefore, we can express p at mth time step as follows:   \begin{eqnarray} p(m \Delta t) = T_m( {\rm C}(l) ){\, }p_0, \quad m=0, 1, 2, \ldots . \end{eqnarray} (24) For the further analysis, we set c and ρ0 to be constant elsewhere. Then ρ0c2∇ · ρ0−1∇ is expressed as c2∇2. Since we use the pseudospectral method in this study, C(l) is implemented in the wavenumber domain as   \begin{eqnarray} {\rm \widetilde{C}}(l,\theta ) = 1 - 2\left( \sum \limits _{n=0}^{l} \frac{ (-1)^n }{(2n+1)!{\, }2^{2n+1}} \theta ^{2n+1} \right)^2, \end{eqnarray} (25) where $${\rm \widetilde{C}}(l,\theta )$$ is C(l) to which the spatial Fourier transform is applied and θ = ckΔt where $$k=\sqrt{k_x^2+k_z^2}$$ in the 2-D domain. If we consider the relation   \begin{eqnarray} \cos (\theta ) = 1 - \displaystyle 2\sin ^2\left( \frac{\theta }{2} \right) \approx 1 - 2\left( \displaystyle \sum \limits _{n=0}^{l} \frac{ (-1)^n }{(2n+1)!{\, }2^{2n+1}} \theta ^{2n+1} \right)^2, \end{eqnarray} (26)we can rewrite $${\rm \widetilde{C}}(l,\theta )$$ as the Taylor expansion of cos (θ) as follows:   \begin{eqnarray} {\rm \widetilde{C}}(l,\theta ) = \sum \limits _{n=0}^{2l+1} \frac{ (-1)^n }{(2n)!} \theta ^{2n}. \end{eqnarray} (27)Then p at the mth time step in the wavenumber domain can be expressed by applying spatial Fourier transform to (24) as   \begin{eqnarray} \widetilde{p}(m \Delta t) = T_m( {\rm \widetilde{C}}(l,\theta ) ){\, }\widetilde{p}_0, \quad m=0, 1, 2, \ldots ,\end{eqnarray} (28)where $$\widetilde{p}$$ is p defined in the wavenumber domain. Eq. (28) explicitly shows that the stability condition of the proposed scheme is $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$, otherwise, $$T_m( {\rm \widetilde{C}}(l,\theta ) )$$ diverges as m increases. Such condition is equivalent to that derived by the Von Neumann or Fourier analysis which sets the pressure field as the plane wave expressed in the harmonic function (Strikwerda 2004). Fig. 1 illustrates the region of θ marked in black where the stability criterion is violated with respect to the temporal order of accuracy. The stability region gets enlarged as the order of accuracy in time increases, which is referred as the relaxation of the Courant–Friedrichs–Lewy (CFL) condition. The small regions in the middle are usually located at the odd multiple of π, which are marked because the Taylor series inaccurately approximate the cosine function and are slightly larger than one at these points. Although such a small error causes divergence of the solution, it can be acceptable depending on the modelling conditions. In order to determine the acceptability of the modelling condition or to estimate the quality of the result of the numerical simulation, more informative condition than the stability criterion is needed. Figure 1. View largeDownload slide Regions violating $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$ are marked in black. The wave components of these regions get diverging as m increases. Figure 1. View largeDownload slide Regions violating $$|{\rm \widetilde{C}}(l,\theta )| \le 1$$ are marked in black. The wave components of these regions get diverging as m increases. If l of (27) goes to infinity, $${\rm \widetilde{C}}(l,\theta )$$ converges to cos (θ). Then, the amount of the truncation error of the numerical solution of a wave component θ at the mth time step can be written as follows:   \begin{eqnarray} {\rm error}(m,l,\theta ) = \left( T_m({\rm \widetilde{C}}(l,\theta )) - T_m(\cos (\theta )) \right)\widetilde{p}_0. \end{eqnarray} (29) Relation (29) is a general expression of the misfit between the numerical result and the analytic solution expressing not only dispersion error due to the time discretization, but also the amplification error when the stability condition is violated. If we set a condition as   \begin{eqnarray} |T_m({\rm \widetilde{C}}(l,\theta )) - T_m(\cos (\theta ))| \le \epsilon , \end{eqnarray} (30)we can obtain θmax, the maximum allowable θ with respect to l satisfying the criterion (30) after m time steps. In this study, we set the maximum m and ε as 50 000 and 0.01, respectively, which implies that the error of the solution in the wavenumber domain where θmax is less than or equal to 1 per cent during 50 000 time steps. Table 1 represents θmax of our error criterion with respect to the temporal order of accuracy from 2 to 28. As the order of accuracy increases, θmax is extended so that more relaxed modelling conditions, such as larger Δt, faster wave propagation speed or smaller grid size, are allowed. In Fig. 2, θmax is represented by grey dots and shown together with the instability region marked in black in Fig. 1 for exact comparison. It is shown that θmax tends to increase linearly with respect to the order of accuracy. The error conditions of Table 1 is much stricter than the stability conditions because most grey dots are located at far more left positions to the instability regions. If less restrictive error condition, such as much smaller m or larger ε, is used, θmax will increase and may move to the right position of the small black region which causes the instability in Fig. 2. Figure 2. View largeDownload slide The grey dots are the maximum allowable θ satisfying eq. (30) with respect to the order of accuracy in temporal dimension when m and ε are 50 000 and 0.01, respectively. Figure 2. View largeDownload slide The grey dots are the maximum allowable θ satisfying eq. (30) with respect to the order of accuracy in temporal dimension when m and ε are 50 000 and 0.01, respectively. Figure 3. View largeDownload slide Computational domain of the numerical simulation whose size is 4π × 4π km2. Homogeneous velocity of 5 km s−1 and density of 1 $$\text{g {}cm}^{-3}$$ is used with a grid size of 0.0491 km. Four black squares illustrate receivers to record time traces located at π km for z and 2π (R1), 2.5π (R2), 3π (R3) and 3.5π km (R4) for x, respectively. Figure 3. View largeDownload slide Computational domain of the numerical simulation whose size is 4π × 4π km2. Homogeneous velocity of 5 km s−1 and density of 1 $$\text{g {}cm}^{-3}$$ is used with a grid size of 0.0491 km. Four black squares illustrate receivers to record time traces located at π km for z and 2π (R1), 2.5π (R2), 3π (R3) and 3.5π km (R4) for x, respectively. Table 1. θmax with respect to the temporal order of accuracy (2l + 2) from 2 to 28 when m and ε are 50 000 and 0.01, respectively. 2l + 2  2  4  6  8  10  12  14  16  18  20  22  24  26  28  θmax  0.02  0.5  1.45  2.52  3.14  5.14  6.11  6.28  9.04  9.42  12  12.54  14.94  15.65  2l + 2  2  4  6  8  10  12  14  16  18  20  22  24  26  28  θmax  0.02  0.5  1.45  2.52  3.14  5.14  6.11  6.28  9.04  9.42  12  12.54  14.94  15.65  View Large Now we briefly discuss on how the proposed scheme is applied to real situations. If we assume that h = Δx = Δz, the maximum expressible k by the spatial discretization is $$\sqrt{2}\pi /h$$. Then, we can obtain Δt for given maximum velocity cmax, kmax and θmax as follows:   \begin{eqnarray} \Delta t = \frac{\theta _{\text{max}}h}{\sqrt{2}\pi c_{\text{max}}}. \end{eqnarray} (31) If we set the modelling time as tmax, the number of time steps Nt can be calculated as tmax/Δt. Then, Nop, the total number of operation of c2∇2, is obtained as (2l + 1)Nt using (11) and (15) as follows:   \begin{eqnarray} N_{\text{op}}=\frac{\sqrt{2}\pi c_{\text{max}} t_{\text{max}}}{\theta _{\text{max}} h}(2l+1). \end{eqnarray} (32) Since Nop is proportional to the computational cost of the modelling, we can examine the efficiency using (32) beforehand. Table 2 represents Δt, Nt and Nop with respect to the order of accuracy of the scheme to satisfy accuracy condition listed in Table 1 when cmax is 5 km s−1, h is 10 m and tmax is 10 s. The error of the numerical solution for θ less than θmax are bounded in 1 per cent since Nt is less than 50 000 for all cases. We can confirm that Nop tends to decrease as the higher-order method is used. Speed up factor for each Δt in Table 2 is the inversed value of normalized Nop by that of the fourth-order method, which shows that the simulation can be performed two or three times faster by higher-order methods than by fourth-order method. This tendency does not continue; the speed up factor appears to converge to a value between three and four in this case. If we assume that the maximum required frequency is given as 100 Hz, the optimal order of accuracy of the scheme is 18 or 20 because Δt should be less than 5 ms. Then we can perform the simulation almost three times faster than that by the fourth-order scheme with tiny time step length. In case that the target Δt is given, it is convenient to select the temporal order of accuracy of the scheme corresponding to θmax larger than but closest to $$\theta =\sqrt{2}\pi c_{\text{max}}\Delta t/h$$ because the larger the time steps are, the more efficient the simulation, according to the tendency of the accuracy analysis. For the example above, when Δt is given as 5 ms, we should select the order of accuracy as 22 because θ = 11.107 is less than and closest to θmax = 12. Then Nop of the case is 42 000, which corresponds to 3.17 if the value is converted to speed up factor with respect to Nop of the fourth-order scheme. Table 2. Parameters to estimate the efficiency of the modelling with respect to the order of accuracy from 2 to 28 satisfying the accuracy criterion listed in Table 1 when cmax is 5 km s−1, h is 10 m and tmax is 10 s. 2l + 2  4  6  8  10  12  14  16  18  20  22  24  26  28  Δt (ms)  0.23  0.65  1.13  1.41  2.31  2.75  2.83  4.07  4.24  5.40  5.65  6.73  7.04  Nt × 10−2  444  153  88  71  43  36  35  25  23  19  17  15  14  Nop × 10−3  1332  766  617  636  475  473  530  418  488  389  408  372  383  Speed up  1  1.74  2.16  2.09  2.80  2.82  2.51  3.19  2.98  3.43  3.27  3.59  3.48  2l + 2  4  6  8  10  12  14  16  18  20  22  24  26  28  Δt (ms)  0.23  0.65  1.13  1.41  2.31  2.75  2.83  4.07  4.24  5.40  5.65  6.73  7.04  Nt × 10−2  444  153  88  71  43  36  35  25  23  19  17  15  14  Nop × 10−3  1332  766  617  636  475  473  530  418  488  389  408  372  383  Speed up  1  1.74  2.16  2.09  2.80  2.82  2.51  3.19  2.98  3.43  3.27  3.59  3.48  View Large 4 NUMERICAL EXPERIMENTS In this section, numerical examples are presented to verify the properties of the proposed method explained in the analysis section. 4.1 Initial value problem 4.1.1 Homogeneous model First, the numerical modelling of the initial value problem in 2-D homogeneous media is performed. As shown in Fig. 3, the area of computational domain is 4π × 4π km2 which composed of 512 equally spaced discrete points with a grid interval of 0.0491 km for each spatial dimension and with the homogeneous P-wave speed c of 5 km s−1 and density ρ0 of 1 $$\text{g {}cm}^{-3}$$. The periodic boundary condition is naturally imposed because the pseudospectral method uses the harmonic basis functions to calculate the spatial derivatives. We use a Gaussian profile $$p_0=e^{-240\left((x-2\pi )^2+(z-2\pi )^2\right)}$$. For v at the initial step, we set it as zero for the simplicity. Four different Δt values, such as 1, 5, 10 and 20 ms, are used for the simulation over 50 s. Before performing the modelling, we can predict the performance as in the previous section. Table 3 represents the maximum expressible θ, the required order of accuracy of the scheme, Nop of the simulation under each condition and the relative speed up factor with respect to Nop of the simulation using 1 ms. As Δt increases, a higher-order scheme is required to accommodate the maximum expressible θ during the modelling. The larger the Δt, the smaller the Nop is needed; the computational cost for the numerical simulation with Δt = 20 ms is less than one-third of that with Δt = 1 ms as shown in Table 3. Wavefields at 1, 10, 20, 30, 40 and 50 s are illustrated in Fig. 4, with the results of different Δt values shown at each quadrant of the snapshots. All the wavefields are notably equivalent in terms of the scale and location of the events. Time traces of p at the receivers for each Δt from 49 to 50 s are shown by Fig. 5. We can confirm that all traces of four different Δt values are nearly identical without dispersion or dissipation phenomena. For quantitative comparison, the L2 errors of the numerical solution in Fig. 4 and the analytic solutions are obtained and represented in Table 4. The error values are negligibly small throughout the simulation for every Δt which, ascertain the feasibility of the strategy to select the order of accuracy. Figure 4. View largeDownload slide Combined wavefields obtained by numerical modelling at (a) 1, (b) 10, (c) 20, (d) 30, (e) 40 and (f) 50 s for each Δt at each quadrant. Figure 4. View largeDownload slide Combined wavefields obtained by numerical modelling at (a) 1, (b) 10, (c) 20, (d) 30, (e) 40 and (f) 50 s for each Δt at each quadrant. Figure 5. View largeDownload slide Plots of time traces recorded at (a) R1, (b) R2, (c) R3 and (d) R4 for each Δt values. Figure 5. View largeDownload slide Plots of time traces recorded at (a) R1, (b) R2, (c) R3 and (d) R4 for each Δt values. Table 3. Maximum expressible θ, the required order of accuracy, Nt and Nop with respect to Δt for the homogeneous velocity model with $$c=5{\, }\text{km s}^{-1}$$, h = 0.0491 km. Speed up factor of each case is the inversed value of normalized Nop by that of Δt = 1 ms. Δt (ms)  1  5  10  20  θ  0.4525  2.2627  4.5255  9.05  2l + 2  4  8  12  18  Nt  50 000  10 000  5000  2500  Nop  300 000  140 000  110 000  85 000  Speed up  1  2.143  2.727  3.529  Δt (ms)  1  5  10  20  θ  0.4525  2.2627  4.5255  9.05  2l + 2  4  8  12  18  Nt  50 000  10 000  5000  2500  Nop  300 000  140 000  110 000  85 000  Speed up  1  2.143  2.727  3.529  View Large Table 4. L2 errors between numerical and analytic wavefields at 1, 10, 20, 30, 40 and 50 s for each Δt. t (s)  1  10  20  30  40  50  L2 error (1 ms)  8.37E − 9  6.81E − 8  9.55E − 8  1.14E − 7  1.30E − 7  2.11E − 7  L2 error (5 ms)  5.72E − 10  2.19E − 9  9.09E − 10  2.36E − 8  2.24E − 8  3.91E − 8  L2 error (10 ms)  8.08E − 10  1.19E − 9  4.03E − 9  8.72E − 9  8.15E − 9  9.92E − 9  L2 error (20 ms)  3.57E − 10  6.71E − 10  7.13E − 10  2.28E − 9  3.60E − 9  4.34E − 9  t (s)  1  10  20  30  40  50  L2 error (1 ms)  8.37E − 9  6.81E − 8  9.55E − 8  1.14E − 7  1.30E − 7  2.11E − 7  L2 error (5 ms)  5.72E − 10  2.19E − 9  9.09E − 10  2.36E − 8  2.24E − 8  3.91E − 8  L2 error (10 ms)  8.08E − 10  1.19E − 9  4.03E − 9  8.72E − 9  8.15E − 9  9.92E − 9  L2 error (20 ms)  3.57E − 10  6.71E − 10  7.13E − 10  2.28E − 9  3.60E − 9  4.34E − 9  View Large 4.1.2 Heterogeneous model Now, we perform numerical experiments on whether the proposed modelling scheme can accurately consider the heterogeneity of the material models. The Marmousi-2 model (Martin et al. 2006) illustrated in Fig. 6 are used for the simulation. The computational domain is 17 km for x-axis and 3.5 km for z-axis with a grid size of 12.5 m. The maximum and minimum values of c are 4.7 and 1.028 km s−1, respectively. As with the homogeneous model example, periodic boundary condition is applied to all edges. The initial value for p is the Gaussian profile $$p_0=e^{-240\left((x-8.5)^2\,+\,(z-1.75)^2\right)}$$ and zero for the particle velocity. We perform the numerical modelling over 50 s using four different Δt values, 1, 5, 10, and 20 ms and record the time traces for each Δt at 0.875 km for z-axis and 8.5 (R1), 10.625 (R2), 12.75 (R3) and 14.875 km (R4) for x-axis. Figure 6. View largeDownload slide P-wave velocity (upper) and density (lower) of the Marmousi-2 model. Figure 6. View largeDownload slide P-wave velocity (upper) and density (lower) of the Marmousi-2 model. Table 5 represents the maximum expressible θ, the required order of accuracy of the scheme, Nop of the simulation under each condition and the relative speed up factor with respect to Nop of the simulation using 1 ms. Similar to the previous results using the homogeneous model, it can be seen that the Nop decreases as the Δt increases. The speed up factor shows that computational cost of the modelling using 5, 10, and 20 ms are less than half of that using 1 ms. Time traces for each Δt are compared in Fig. 7. In order to see the effect of accumulated errors over time, we compared the traces from 49 to 50 s. We can find that all traces of four different Δt values are nearly identical without dispersion or dissipation phenomena. The proposed method is shown to yield consistent and stable numerical results irrespective of time step lengths when using the heterogeneous model. Figure 7. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Figure 7. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Table 5. Maximum expressible θ, the required order of accuracy, Nt, Nop with respect to Δt for the simulation using the Marmousi-2 model with $$c_{\text{max}}=4.7{\, }\text{km s}^{-1}$$, h = 0.0125 km. Speed up factor of each case is inversed value of normalized Nop by that of Δt = 1 ms. Δt (ms)  1  5  10  20  θ  1.6705  8.3526  16.7052  33.4105  2l + 2  8  18  30  52  Nt  50 000  10 000  5000  2500  Nop  700 000  340 000  290 000  255 000  Speed up  1  2.059  2.414  2.745  Δt (ms)  1  5  10  20  θ  1.6705  8.3526  16.7052  33.4105  2l + 2  8  18  30  52  Nt  50 000  10 000  5000  2500  Nop  700 000  340 000  290 000  255 000  Speed up  1  2.059  2.414  2.745  View Large 4.2 Source problem with absorbing boundary condition In this section, we examine the feasibility of the proposed method in the context of source problems with the absorbing boundary condition for practical usages. Most conditions are identical to the previous example using the Marmousi-2 model. The Ricker wavelet $$q(t)=(1-2\pi ^2(f_0t-\beta )^2)e^{-\pi ^2(f_0t-\beta )^2}$$ is used as a source wavelet with f0 = 10 Hz and β = 1. The range of frequency band of the wavelet is 30 Hz with a dominant frequency component at 10 Hz. The wavelet is imposed to p with a distributed monopole source $$f(x,z)=e^{-500((x-x_s)^2+(z-z_s)^2)}$$ where xs = 8.5 km, zs = 0.0625 km. The following expression is used to consider the source term:   \begin{eqnarray} p(t+\Delta t) = p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + f(x,z)\frac{q(t) + q(t+\Delta t)}{2}\Delta t. \end{eqnarray} (33) Time traces are recorded by a line of receivers located at the depth of 0.0875 km. The numerical simulation with Δt = 20 ms is excluded because the Nyquist frequency of the numerical solution is 25 Hz which is lower than maximum frequency component of the source wavelet. Instead, we perform the simulation with Δt = 16 ms, whose Nyquist frequency is slightly less than 30 Hz. The required order of accuracy of the scheme and Nop are equal to those listed in Table 5 for Δt of 1, 5 and 10 ms. The required order of accuracy for Δt = 16 ms is 44 because the maximum expressible θ is 26.728. The parameter Nop for the case is 268 750 and speed up factor is 2.604. We apply the absorbing boundary condition suggested by Lavelle & Thacker (2008), a kind of sponge-type damping layer mimicking the perfectly matched layer (PML) because of the instability problem for large Δt. Although this method is not as efficient as PML, performs much better than the simple-sponge. From our experience, the simple-sponge requires a layer of 1.5 km thickness, whereas this method requires layer of 0.5 km to damp out the incident waves perfectly. In the damping area, the quadratic damping profile $$d_x(x)=\xi (x-x_d)^2/L_d^2$$ is used where xd is the starting position. In this study, we set ξ to 100, Ld to 0.5 km as stated above, and xd to 0 and 17 km for left and right boundaries, respectively. In the same way, we set quadratic damping profile in the z-direction, dz(z), with equivalent ξ and Ld where zd are 0 and 3.5 km for top and bottom layers. Fig. 8 shows that kinematics of recorded waves are in good agreement. It is also confirmed that the waves are effectively dampened by 0.5-km-thick absorbing layers. Time traces are plotted in Fig. 9 recorded from the seismogram at 8.5 (R1), 10.625 (R2), 12.75 (R3) and 14.875 km (R4) for x-axis and 0.0625 for z-axis. We can confirm that the arrival time and amplitude of the first arrival and subsequent signals are equivalent to each other regardless of Δt. Figure 8. View largeDownload slide Combination of seismograms for four different Δt values recorded at a line receiver located at 0.0625 km for z. Figure 8. View largeDownload slide Combination of seismograms for four different Δt values recorded at a line receiver located at 0.0625 km for z. Figure 9. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. Figure 9. View largeDownload slide Plots of time traces recorded at R1 (a), R2 (b), R3 (c) and R4 (d) for each Δt values. 5 CONCLUSIONS In this study, we suggested the staggered-time method whose order of accuracy in temporal dimension can be extended arbitrarily using the simple recurrence relations. The method is applied to the linear acoustic wave system of pressure and particle velocity in a stagnant and stationary media. We analysed the method quantifying the error of the numerical solution and proposed the error criterion that guarantees 1 per cent error during 50 000 time steps. We also propose a simple strategy how to use the result of the analysis; rather than adjusting the modelling conditions to satisfy the stability and accuracy of a given modelling scheme, we set temporal order of accuracy of the scheme satisfying the error condition for given time step length, grid intervals and velocity model. Such property allows us to flexibly determine modelling conditions. For instance, when the scheme is used in the reverse time migration algorithm, there is no need to reduce the time step length more than necessary to match the desired spatial resolution because of the stability issue. The proposed modelling scheme allows not only the flexibility, but also the efficiency; the calculation efficiency is improved significantly with larger time step length and higher order of accuracy, which is proven by the analysis and numerical examples. By virtue of this property, it is efficient for the full waveform inversion algorithms when the method is applied. Because we can accurately perform the simulation with a sufficiently large time step length corresponding to the low frequencies generally used during the inversion process. We also confirmed by the numerical examples that the method perfectly considers the heterogeneity of the media. We proposed how to properly impose the source wavelet, which yield equivalent solutions regardless of the time step lengths whose Nyquist frequency is less than that of the maximum frequency of the wavelet. As stated above, we expect that the proposed method is effective for the geophysical applications. This method can also be applied to simulate other wave phenomena such as electromagnetic or elastodynamic waves which also have many practical applications, that is, ground penetrating radar or medical imaging algorithms. ACKNOWLEDGEMENTS This work was supported by the Energy Efficiency and Resources Core Technology Program of the Korea Institute of Energy Technology Evaluation and Planning (KETEP), granted financial resource from the Ministry of Trade, Industry and Energy, Republic of Korea (20132510100060). And the Institute of Engineering Research at Seoul National University provided research facilities for this work. REFERENCES Arnold D.N., Brezzi F., Cockburn B., Marini L.D., 2002. Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. , 39( 5), 1749– 1779. Google Scholar CrossRef Search ADS   Ashcroft G., Zhang X., 2003. Optimized prefactored compact schemes, J. Comput. Phys. , 190( 2), 459– 477. Google Scholar CrossRef Search ADS   Baysal E., Kosloff D.D., Sherwood J.W., 1983. Reverse time migration, Geophysics , 48( 11), 1514– 1524. Google Scholar CrossRef Search ADS   Bogey C., Bailly C., 2004. A family of low dispersive and low dissipative explicit schemes for flow and noise computations, J. Comput. Phys. , 194( 1), 194– 214. Google Scholar CrossRef Search ADS   Carcione J., 1996. A 2D Chebyshev differential operator for the elastic wave equation, Comput. Methods Appl. Mech. Eng. , 130( 1–2), 33– 45. Google Scholar CrossRef Search ADS   Chen J.-B., 2009. Lax–Wendroff and Nyström methods for seismic modelling, Geophys. Prospect. , 57( 6), 931– 941. Google Scholar CrossRef Search ADS   Cohen G., 1986. Fourth-order schemes for the 2-d wave equation in a homogeneous medium, in SEG Technical Program Expanded Abstracts 1986 , pp. 632– 635, Society of Exploration Geophysicists. Dablain M., 1986. The application of high-order differencing to the scalar wave equation, Geophysics , 51( 1), 54– 66. Google Scholar CrossRef Search ADS   De Basabe J.D., Sen M.K., Wheeler M.F., 2008. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion, Geophys. J. Int. , 175( 1), 83– 93. Google Scholar CrossRef Search ADS   Dumbser M., Käser M., 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—II. The three-dimensional isotropic case, Geophys. J. Int. , 167( 1), 319– 336. Google Scholar CrossRef Search ADS   Hall B., 2015. Lie Groups, Lie Algebras, and Representations: An Elementary Introduction , vol. 222, Springer. Google Scholar CrossRef Search ADS   Hermann V., Käser M., Castro C.E., 2011. Non-conforming hybrid meshes for efficient 2-D wave propagation using the discontinuous Galerkin method, Geophys. J. Int. , 184( 2), 746– 758. Google Scholar CrossRef Search ADS   Hesthaven J.S., Warburton T., 2007. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications , Springer Science & Business Media. Käser M., Dumbser M., 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—I. The two-dimensional isotropic case with external source terms, Geophys. J. Int. , 166( 2), 855– 877. Google Scholar CrossRef Search ADS   Klin P., Priolo E., Seriani G., 2010. Numerical simulation of seismic wave propagation in realistic 3-D geo-models with a Fourier pseudo-spectral method, Geophys. J. Int. , 183( 2), 905– 922. Google Scholar CrossRef Search ADS   Komatitsch D., Vilotte J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Kosloff D.D., Baysal E., 1982. Forward modeling by a Fourier method, Geophysics , 47( 10), 1402– 1412. Google Scholar CrossRef Search ADS   Kosloff D., Kosloff R., 1983. A Fourier method solution for the time dependent Schrödinger equation as a tool in molecular dynamics, J. Comput. Phys. , 52( 1), 35– 53. Google Scholar CrossRef Search ADS   Kosloff D., Queiroz Filho A., Tessmer E., Behle A., 1989. Numerical solution of the acoustic and elastic wave equations by a new rapid expansion method, Geophys. Prospect. , 37( 4), 383– 394. Google Scholar CrossRef Search ADS   Lavelle J., Thacker W., 2008. A pretty good sponge: dealing with open boundaries in limited-area ocean models, Ocean Modelling , 20( 3), 270– 292. Google Scholar CrossRef Search ADS   Lele S.K., 1992. Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. , 103( 1), 16– 42. Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P –SV seismograms, Geophysics, 53( 11), 1425– 1436. Liu Y., Sen M.K., 2009. A new time–space domain high-order finite-difference method for the acoustic wave equation, J. Comput. Phys. , 228( 23), 8779– 8806. Google Scholar CrossRef Search ADS   Long G., Zhao Y., Zou J., 2013. A temporal fourth-order scheme for the first-order acoustic wave equations, Geophys. J. Int. , 194( 3), 1473– 1485. Google Scholar CrossRef Search ADS   Martin G.S., Wiley R., Marfurt K.J., 2006. Marmousi2: an elastic upgrade for Marmousi, Leading Edge , 25( 2), 156– 166. Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Gális M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Pestana R.C., Stoffa P.L., 2010. Time evolution of the wave equation using rapid expansion method, Geophysics , 75( 4), T121– T131. Google Scholar CrossRef Search ADS   Pratt R.G., Shin C., Hick G., 1998. Gauss–Newton and full newton methods in frequency–space seismic waveform inversion, Geophys. J. Int. , 133( 2), 341– 362. Google Scholar CrossRef Search ADS   Renaut R., Fröhlich J., 1996. A pseudospectral Chebychev method for the 2D wave equation with domain stretching and absorbing boundary conditions, J. Comput. Phys. , 124( 2), 324– 336. Google Scholar CrossRef Search ADS   Shin C., Cha Y.H., 2008. Waveform inversion in the Laplace domain, Geophys. J. Int. , 173( 3), 922– 931. Google Scholar CrossRef Search ADS   Shin C., Cha Y.H., 2009. Waveform inversion in the Laplace–Fourier domain, Geophys. J. Int. , 177( 3), 1067– 1079. Google Scholar CrossRef Search ADS   Strikwerda J.C., 2004. Finite Difference Schemes and Partial Differential Equations , SIAM. Google Scholar CrossRef Search ADS   Tal-Ezer H., 1986. Spectral methods in time for hyperbolic equations, SIAM J. Numer. Anal. , 23( 1), 11– 26. Google Scholar CrossRef Search ADS   Tam C.K., Webb J.C., 1993. Dispersion-relation-preserving finite difference schemes for computational acoustics, J. Comput. Phys. , 107( 2), 262– 281. Google Scholar CrossRef Search ADS   Tan S., Huang L., 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation, Geophys. J. Int. , 197( 2), 1250– 1267. Google Scholar CrossRef Search ADS   Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics , 49( 8), 1259– 1266. Google Scholar CrossRef Search ADS   Tessmer E., 2011. Using the rapid expansion method for accurate time-stepping in modeling and reverse-time migration, Geophysics , 76( 4), S177– S185. Google Scholar CrossRef Search ADS   Trefethen L.N., 2000. Spectral Methods in MATLAB , SIAM. Google Scholar CrossRef Search ADS   Virieux J., 1984. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 49( 11), 1933– 1942. Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. Google Scholar CrossRef Search ADS   Whitmore N., 1983. Iterative depth migration by backward time propagation, in SEG Technical Program Expanded Abstracts 1983 , pp. 382– 385, Society of Exploration Geophysicists. Yee K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. , 14( 3), 302– 307. Google Scholar CrossRef Search ADS   APPENDIX A: ABSORBING BOUNDARY CONDITION We implemented the pretty-good-sponge boundary condition (Lavelle & Thacker 2008) because the instability occurs in the damping layer when applying the PML. Although less effective than the PML, it has comparable damping ability with moderate thickness of damping layer. The 2-D acoustic wave system defined in the damping layer is written as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ \mathbf {v} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}-{\rm tr}(\mathbf {D}) & -\rho _0 c^2 \nabla \cdot \\ -{\rho _0}^{-1}\nabla & -\mathbf {D} \end{array}\right]}}{\left[\begin{array}{c}p \\ \mathbf {v} \end{array}\right]}, \end{eqnarray} (A1)where D is a diagonal 2 × 2 matrix composed of damping parameters, that is, D11 = dx(x) and D22 = dz(z) for 2-D problem where dx and dz are the damping profiles for x- and z-directions, respectively. The proposed method defined in the damping layer can be derived by applying (A1) to (4) and (5) omitting the cross-terms between damping parameter and the spatial derivatives as   \begin{eqnarray} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) {}& - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \Big ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \Big )^{n} p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}\mathbf {v}(t), \end{eqnarray} (A2)and   \begin{eqnarray} p(t+\Delta t) \approx p(t) {}& - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1}p\bigg (t+\frac{\Delta t}{2}\bigg ) . \end{eqnarray} (A3) If we set v(t) and p(t + Δt/2) as an averaged value as   $$\mathbf {v}(t) \approx \frac{1}{2}\bigg ( \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) + \mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) \bigg ),\quad p\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \frac{1}{2}\big ( p(t+\Delta t) + p(t) \big ),$$ (A4)we can rewrite eqs (A2) and (A3) as   \begin{eqnarray} \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \approx \mathbf {B}^{-1} \left( \mathbf {A}{\, }\mathbf {v}\bigg (t-\frac{\Delta t}{2}\bigg ) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \frac{1}{\rho _0}\nabla \left ( {\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n} p(t) \right) + O(\Delta t^2), \end{eqnarray} (A5) and   \begin{eqnarray} p(t+\Delta t) \approx \beta ^{-1} \left( \alpha {\, }p(t) - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\left ({\rho _0 c^2}\nabla \cdot \frac{1}{\rho _0}\nabla \right )^{n}{\rho _0 c^2}\nabla \cdot \displaystyle \mathbf {v}\bigg (t+\frac{\Delta t}{2}\bigg ) \right) + O(\Delta t^2), \end{eqnarray} (A6)where   \begin{eqnarray} \mathbf {A} &=& I - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}, \nonumber \\ \mathbf {B} &=& I + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}\mathbf {D}^{2n+1}, \nonumber \\ \alpha &=& 1 - \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1},\nonumber \\ \beta &=& 1 + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n}(d_x+d_z)^{2n+1}. \end{eqnarray} (A7) The temporal order in the damping layer is only two because of the approximation and the time averaging of the variables, which may aggravate the accuracy of the solution in the whole domain if reflections are unignorable. If there is almost no reflection, the wavefield inside the computational domain maintains the desired order of accuracy in time. APPENDIX B: ANALYTIC SOLUTION We can derive the analytic solution of the acoustic wave equation in a homogeneous model using the one-step method, which imposes the periodic boundary condition by default if the system is solved in the wavenumber domain by using the Fourier basis (Trefethen 2000). The linear acoustic wave system (1) defined in the 2-D domain can be rewritten as follows:   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}p \\ u \\ w \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c@{\quad}c}0 & -\rho _0 c^2 \frac{\partial }{\partial x} & -\rho _0 c^2\frac{\partial }{\partial z} \\ -\frac{1}{\rho _0}\frac{\partial }{\partial x} & 0 & 0 \\ -\frac{1}{\rho _0}\frac{\partial }{\partial z} & 0 & 0 \end{array}\right]}}{\left[\begin{array}{c}p \\ u \\ w \end{array}\right]}. \end{eqnarray} (B1) This equation can be converted into system ordinary differential equation form via spatial Fourier transformation of the system variables as follows:   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}\widetilde{p} \\ \widetilde{u} \\ \widetilde{w} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c@{\quad}c}0 & i\rho _0 c^2 k_x & i\rho _0 c^2 k_z \\ \frac{i}{\rho _0} k_x & 0 & 0 \\ \frac{i}{\rho _0} k_z & 0 & 0 \end{array}\right]}}{\left[\begin{array}{c}\widetilde{p} \\ \widetilde{u} \\ \widetilde{w} \end{array}\right]}. \end{eqnarray} (B2) Then, the acoustic wave solution in the time-wavenumber domain is expressed as   \begin{eqnarray} \widetilde{y}(t)=e^{t\widetilde{\mathbf {M}}}\widetilde{y}_0, \end{eqnarray} (B3)where $$\widetilde{y}_0$$ and $$\widetilde{y}(t)$$ is the transformed system variable vector at the initial and t, respectively, and $$\widetilde{\mathbf {M}}$$ is the characteristic matrix of (B2) whose exponential map can be formulated as follows (Hall 2015):   \begin{eqnarray} e^{t\widetilde{\mathbf {M}}}= {{\left[\begin{array}{c@{\quad}c@{\quad}c}\cos (ckt) & \frac{ick_x\sin (ckt)}{k} & \frac{ick_z\sin (ckt)}{k} \\ \frac{ick_x\sin (ckt)}{k} & \frac{k_z^2+k_x^2\cos (ckt)}{k^2} & \frac{k_xk_z\left(\cos (ckt)-1\right)}{k^2} \\ \frac{ick_z\sin (ckt)}{k} & \frac{k_xk_z\left(\cos (ckt)-1\right)}{k^2} & \frac{k_z^2+k_x^2\cos (ckt)}{k^2} \end{array}\right]}}. \end{eqnarray} (B4) Therefore, time-space domain solution can be easily obtained by applying the inverse Fourier transform for both sides of equation (B3) as   \begin{eqnarray} y(t) = \mathcal {F}_{xz}^{-1}\left(e^{t\widetilde{\mathbf {M}}} \mathcal {F}_{xz}\left(y_0\right)\right), \end{eqnarray} (B5)where $$\mathcal {F}_{xz}$$ is the 2-D spatial Fourier transform operator for the xz-plane. APPENDIX C: GENERALIZATION The arbitrary-order staggered time integrator can be applied to other linear wave equations, that is, first-order elastodynamic or electromagnetic wave system. If we define system variables as $$\mathbf {p}\in \mathbb {R}^m$$ and $$\mathbf {q}\in \mathbb {R}^n$$, we can express the first-order wave equations in a general form as   \begin{eqnarray} \frac{\partial }{\partial t}{\left[ \begin{array}{c}\mathbf {p} \\ \mathbf {q} \end{array} \right]} = {{\left[\begin{array}{c@{\quad}c}0 & \mathbf {M_1} \\ \mathbf {M_2} & 0 \end{array}\right]}}{\left[\begin{array}{c}\mathbf {p} \\ \mathbf {q} \end{array}\right]}, \end{eqnarray} (C1)where M1 and M2 are the submatrices of the characteristic matrix M. M1 and M2 have dimensions of m × n and n × m, respectively. As can be found in (1), components of these matrices consist of products of spatial differential operators and material properties. For the 2-D linear elastic wave equation in the isotropic media (Virieux 1986), we can set p and q as the stress σ = [σxx, σzz, σxz]⊤ and the particle velocity v = [u, w]⊤, respectively. M1 and M2 of the system can be written as   \begin{eqnarray} \mathbf {M}_1 = {\left[\begin{array}{c@{\quad}c}(\lambda + 2\mu )\frac{\partial }{\partial x} & \mu \frac{\partial }{\partial z} \\ \mu \frac{\partial }{\partial x} & (\lambda + 2\mu )\frac{\partial }{\partial z} \\ \mu \frac{\partial }{\partial z} & \mu \frac{\partial }{\partial x} \end{array}\right]},\quad \mathbf {M}_2 = \frac{1}{\rho _0}{\left[\begin{array}{c@{\quad}c@{\quad}c}\frac{\partial }{\partial x} & 0 & \frac{\partial }{\partial z} \\ 0 & \frac{\partial }{\partial z} & \frac{\partial }{\partial x} \end{array}\right]}, \end{eqnarray} (C2)where λ and μ are Láme parameters. For Maxwell’s equations in the isotropic media (Yee 1966), we can set p and q as the electric field E = [Ex, Ey, Ez]⊤ and the magnetic field H = [Hx, Hy, Hz]⊤, respectively. M1 and M2 of the system can be written as   \begin{eqnarray} \mathbf {M}_1 = \frac{1}{\varepsilon } {\left[\begin{array}{c@{\quad}c@{\quad}c}0 & -\frac{\partial }{\partial z} & \frac{\partial }{\partial y} \\ \frac{\partial }{\partial z} & 0 & -\frac{\partial }{\partial x} \\ -\frac{\partial }{\partial y} & \frac{\partial }{\partial x} & 0 \end{array}\right]},\quad \mathbf {M}_2 = \frac{1}{\mu } {\left[\begin{array}{c@{\quad}c@{\quad}c}0 & \frac{\partial }{\partial z} & -\frac{\partial }{\partial y} \\ -\frac{\partial }{\partial z} & 0 & \frac{\partial }{\partial x} \\ \frac{\partial }{\partial y} & -\frac{\partial }{\partial x} & 0 \end{array}\right]}, \end{eqnarray} (C3)where ε and μ are the permittivity and the permeability of a medium, respectively. Following the process from (2) to (5), we can apply the arbitrary-order staggered time integrator to system (C1). If we impose p on the integer and q on the half integer of the temporal grid, we can obtain the semi-discretized formulation of (C1) as follows:   \begin{eqnarray} &&{\displaystyle \mathbf {q}\bigg (t+\frac{\Delta t}{2}\bigg ) = \mathbf {q}\bigg (t-\frac{\Delta t}{2}\bigg ) + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \mathbf {M_2} \Big ( \mathbf {M_1}\mathbf {M_2} \Big )^{n} \mathbf {p}(t) + O(\Delta t^{2l+2}),} \nonumber\\ &&{\mathbf {p}(t+\Delta t) = \mathbf {p}(t) + \displaystyle \sum \limits _{n=0}^l \frac{\Delta t^{2n+1}}{(2n+1)! \, 4^n} \mathbf {M_1}\Big ( \mathbf {M_2}\mathbf {M_1}\Big )^{n} \displaystyle \mathbf {q}\bigg (t+\frac{\Delta t}{2}\bigg ) + O(\Delta t^{2l+2}).} \end{eqnarray} (C4) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 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