Efficient temporal high-order staggered-grid scheme with a dispersion-relation-preserving method for the scalar wave modeling

Staggered-grid finite-difference (FD) method is widely used to solve the wave equation for the numerical seismic modeling, and it is a key step of the advanced seismic imaging and inversion problem. However, the conventional FD method is prone to instability and dispersion error due to the insufficient approximation accuracy. In this work, we propose an efficient temporal high-order finite-difference (FD) scheme with the cross-rhombus stencil. Compared with the standard cross-rhombus method, the new method has less computational cost due to we simplify the FD scheme. Moreover, the dispersion relation of the new method is easy to obtain the dispersion-relation-preserving (DRP) FD coefficients, which can significantly alleviate the spatial and temporal dispersion errors. Dispersion and stability analyses indicate that the new scheme has better performance in seismic modeling than the conventional method, and numerical experiments also indicate that the new scheme can effectively mitigate dispersion error and improve the numerical accuracy.


Introduction
Staggered-grid finite-difference methods have been extensively applied in the seismic wave simulations due to their straightforward implementation and high computing efficiency (Kindelan et al., 1990;Moczo et al., 2000;Etgen and O¡¯Brien, 2007;Moczo et al., 2011;Moczo et al., 2014;Etemadsaeed et al., 2016;Liu et al., 2019;Zhang et al., 2022). Highorder approximation for temporal derivatives in the staggered-grid FD scheme contributes to suppressing the temporal dispersion errors, and enhancing the stability with a large time step. However, the explicit high-order temporal derivative approximation in the FD scheme is always unstable (Chen, 2007;2011). Generally, we use a second-order temporal approximation and an arbitrary even-order spatial approximation to solve the scalar wave equation.
To improve the temporal accuracy, Dablain (1986) proposed a new FD scheme based on the Lax-Wendroff approach, which can reach fourth-order accuracy in the temporal derivative approximation but has a limitation in the case of the large velocity contrast.
Frontiers in Earth Science 01 frontiersin.org Chen (2007Chen ( , 2011 further developed the fourth-and sixth-order schemes and analyzed the stability condition in the high-order cases. Alternatively, Liu and Sen (2009) proposed the time-space domain FD coefficients by incorporating the dispersion relation of the temporal and spatial terms. The time-space domain method can reach arbitrary even-order accuracy along with some specific propagation angles. However, it is still second-order accuracy along with other angles (Liu and Sen, 2009). To further improve the accuracy, Liu and Sen (2013) proposed a novel rhombus stencil. This new stencil with the time-space domain FD coefficients can reach arbitrary-order accuracy in both temporal and spatial approximations. However, the standard rhombus stencil is not a computationally friendly method for large-scale modeling, and it will increase the computational cost exponentially for the high-order cases. Afterwards, Tan and Huang (2014a), Tan and Huang (2014b) proposed an effective FD stencil with the sixth-order accuracy in the time approximation. Tan's stencil is similar to the rhombus stencil, but it involves fewer grid nodes outside the cross axis compared to the standard rhombus stencil, thus reducing the computational cost significantly. Wang E. et al. (2016) generally defined this stencil as the cross-rhombus stencil with arbitrary even-order temporal accuracy. The cross-rhombus stencil contains a large cross stencil and a small rhombus stencil. Among them, the small rhombus stencil increases the temporal accuracy and ensures computational efficiency, while the large cross stencil has a high-order spatial accuracy. Then, Ren et al. (2017) developed the cross-rhombus stencil in the staggered-grid FD scheme, and presented two methods for solving the FD coefficients. Wang et al. (2019) further developed the cross-rhombus stencil in the 3D case with the general cuboid grid.
To mitigate the dispersion error, the dispersion relation of the FD scheme should require many wavenumbers, because the spatial dispersion error usually comes from the high-wavenumber component. However, the conventional Taylor series expansion (TE) method for solving FD coefficient satisfies the dispersion relation near the zero wavenumbers, so it is prone to dispersion. The optimization method is a feasible way to obtain the FD coefficients (Liu, 2013;Zhang and Yao, 2013;Tan and Huang, 2014b;Chen et al., 2020) for mitigating dispersion, where the dispersionrelation-preserving method (Wang and Teixeira, 2003;Ye and Chu, 2005;Liang et al., 2015;Etemadsaeed et al., 2016;Chen et al., 2020) has been widely concerned because of its simplicity and easy implementation. The DRP-based method expands dispersion relation into an over-determined system associated with a series of wavenumbers and propagation angles, and then solves this overdetermined system numerically to obtain the FD coefficients in the least square sense. The DRP-based FD coefficients satisfy a series of wavenumbers from low to high in the sense of the least square, thus the DRP-based method can effectively mitigate the temporal and spatial dispersion error.
The DRP-based FD coefficients have been successfully applied to the temporal high-order scheme (cross-rhombus stencil) in the regular grid Chen et al. (2020), in which the DRP-based coefficients can significantly mitigate dispersion error, while the cross-rhombus  1141220   TABLE 1 Abbreviation table of different FD methods used for dispersion analyses.   Abbreviations  FD coefficients  FD stencils   TE-C-S  TE-based space domain coefficients  Cross stencil   TE-C-TS  TE-based time-space domain coefficients  Cross stencil   TE-CR-TS  TE- stencil can effectively improve the temporal approximation accuracy. However, for the staggered-grid FD scheme, the DRPbased coefficients cannot be directly obtained because the dispersion relation is difficult to be extended into an over-determined system. Liang et al. (2018) has presented a special FD scheme with a high computational efficiency, in which the second-order FD operator is used to approximate some partial derivatives rather than the global high-order FD operator. And such replacement simplifies the dispersion relation into a form of the linear equation. Motivated by (Liang et al., 2018;Zhou et al., 2022), we propose a general simplified FD scheme for the temporal high-order modeling with a cross-rhombus stencil. The new scheme contains a cross stencil Frontiers in Earth Science 03 frontiersin.org with the analytic FD coefficients and a cross-rhombus stencil with the DRP-based coefficients for different partial derivatives. The cross stencil can simplify the dispersion relation, which makes it easy to construct the over-determined system. The cross-rhombus stencil can make the FD scheme maintain a high-order temporal approximation. The dispersion relation of our new FD scheme can be expanded to an over-determined system with a series of wavenumbers and angles. Solving this over-determined system by the numerical methods (Wang et al., 2014;Wang et al., 2016 Y.;Chen et al., 2020;Wu et al., 2020;Li et al., 2022), we obtain the DRP-based FD coefficients. Therefore, the new FD scheme has three advantages: 1. The DRP-based FD coefficients can effectively mitigate the dispersion error; 2. It still has a temporal higher-order approximation accuracy; 3. The computational cost of the proposed method is significantly reduced compared with the standard crossrhombus scheme.

Review of the standard cross-rhombus scheme
The 2D first-order acoustic wave equations with the constant density are where K = ρυ 2 is the bulk modulus, ρ(x, z) is the density and v(x, z) is the velocity. p(x, z, t) is the pressure, and v = [υ x , υ z ] T is the particle velocity vector. The staggered-grid FD scheme with the cross-rhombus stencil for above equations is Here, p l i,j = p(x + ih, z + jh, t + lτ), h is the grid spacing and τ is the time step. The FD operators H CR x and D CR x are and The superscript CR represents the cross-rhombus stencil composed of a standard cross stencil and a rhombus stencil (as shown in Figure 1). a m,n and b m,n represent the FD coefficients of the operators H CR x and D CR x , respectively. M and N are the spatial and temporal operator length parameters, respectively. Generally, when temporal operator length N > 3, the accuracy increases far less than the increase of the calculation cost. Thus, we recommend that N = 3 is enough. M and N represents the (2M)th-order accuracy in space and (2N)th-order accuracy in time respectively. And the FD operators along the z-axis (H CR z and D CR z ) can be defined in the same way.
Assuming plane wave propagating in the grid, we let p l m,n = p 0 0,0 e i(k x mh+k z nh−ωlτ) , υ l+1/2 x m+1/2,n = υ where k x = k cos(θ) and k z = k sin(θ) are the wavenumbers in x-and z-axes, respectively. θ is a propagation angle of the plane wave, ω is the angular frequency and i = √ −1. Substituting Eq. 5 into Eq. 2, we obtain Here, r = vτ/h is the Courant number. Generally, the operators H CR and D CR adopt the same FD coefficients, i.e. a m.n = b m,n .
Then, Equation 6 can be rewritten as Eq. 8 represents the time-space domain dispersion relation of the staggered-grid FD scheme with the cross-rhombus stencil. Using the Taylor series to expand the trigonometric functions with respect to the propagation angle θ, we can obtain the time-space domain TEbased FD coefficients (Ren et al., 2017). The cross-rhombus stencil with the TE-based FD coefficients can achieve arbitrary even-order temporal accuracy, thus mitigating the temporal dispersion error significantly.
Frontiers in Earth Science 05 frontiersin.org

A new simplified staggered-grid FD scheme with the cross-rhombus stencil
The TE-based FD coefficients satisfy the dispersion relation within a limited wavenumber bandwidth, resulting in the highwavenumber components of seismic wavefield are prone to the spatial dispersion. However, the dispersion relation of the standard staggered-grid scheme is a quadratic equation, which is difficult to expand into the over-determined system for solving DRP-based FD coefficients. In this part, we develop a new simplified staggered-grid FD scheme. The new scheme can not only easily obtain the overdetermined system for the DRP-based coefficients, but also greatly reduce the computational cost.
The FD operators D CR and H CR use same coefficients (a m,n = b m,n ), which causes the dispersion relation to be a secondorder non-linear equation. To obtain a simple dispersion relation, we propose a new simplified FD scheme as Here, the superscript C represents the cross stencil. For example, the FD operator D C x can be defined as We use the cross-stencil-based operators D C to replace part of the operator D CR in the new scheme. Thus, the new scheme contains a cross-rhombus stencil and a cross stencil for different partial derivatives.  For convenience, D C adopts the analytic time-space domain FD coefficients Here, L represents the length of the analytic FD operator. And the analytical cross stencil simplifies the original second-order dispersion relation to a linear form. The new dispersion relation can be easily extended to an over-determined linear system for solving the wide-bandwidth FD coefficients. Moreover, the analytical cross stencil has less computational cost compared with the standard cross-rhombus scheme, especially in the high-order cases.

Determining FD coefficients of the new stencil by the DRP-based method
In this part, we present the method for solving the DRP-based FD coefficients of the new scheme. We assume that the plane wave propagating in the grid, and then substitute Eq. 5 into Eq. 9, yield the new dispersion relation Clearly, the new dispersion relation can be easily extended to the linear system satisfying a series of wavenumbers and propagation angles. It is worth noting that if the FD operator H uses cross stencil, for example: H C and D CR are applied to Equation 9, the dispersion relation does not change. And the cross stencil is applied to the FD operator H or D, the corresponding FD schemes are equivalent. Following our previous work , we define a new function ψ m,β,θ to represent the weights of FD coefficients a m,0 in Eq. 12. Then ψ m,β,θ can be denoted as ψ m,β,θ = sin ((m − 0.5) β cos (θ)) L ∑ l=1 b l,0 sin ((l − 0.5) β cos (θ)) Here, β = kh and θ represents the propagation angle. Similarly, we define another function χ m,n,β,θ to represent the weights of a m,n . The function χ m,n,β,θ is defined as χ m,n,β,θ = 2 sin ((m − 0.5) β cos (θ)) cos (nβ sin (θ)) Since a m,n = a n,m , we define a new function Therefore, the dispersion relation (Eq. 12) of the new FD scheme can be rewritten as Then we extend ψ m,β,θ to a matrix involving a series of β and a fixed angle θ, and the matrix is ] .
Where β i = β max /ξ*i, β max = 2πf max /v with respect to the maximum frequency of the seismic wavefield Chen et al. (2022) and π is the circular constant. Similarly, we extend the function φ n,β,θ to the matrix.
This over-determined system has ζ × ξ rows and 1 + M + N 2 /4 columns with even number N and 1 + M + (N − 1) 2 /4 columns with odd number N. We can easily obtain the DRP-based FD coefficients by solving this over-determined system. We also introduce the simplified FD scheme in the 3D case, the details are shown in the Supplementary Material S1.

Numerical dispersion and stability analyses 3.1 Numerical dispersion analysis
In this part, we analyse the dispersion characteristics of our new scheme. The phase velocity of the new FD scheme can be expressed as where q is Thus, the dispersion parameter δ of the new FD scheme is defined as Frontiers in Earth Science 08 frontiersin.org If δ is not equal to 1, the FD scheme suffers from the numerical dispersion, i.e., has the spatial dispersion error (δ < 1) or temporal dispersion error (δ > 1). We analyze and compare the dispersion parameter δ of the new FD scheme with the other three methods, and the abbreviations of these methods are listed in Table 1. The dispersion curves of δ varying with the kh are shown in Figure 2.
It can be seen that the cross-stencil-based FD schemes (TE-C-S and TE-C-TS methods) have obviously temporal dispersion error (δ > 1). The corresponding temporal dispersion of the crossrhombus stencil (TE-CR-TS and DRP-CR-TS methods) is alleviated (Figures 2C, D) due to the temporal high-order approximation. It is worth noting that the proposed method (DRP-CR-TS method)   satisfies the widest range of the wavenumber (kh), which can mitigate the spatial dispersion error considerably.

Stability analysis
According to the dispersion relation of the new FD scheme, we obtain It is clear that Then, we obtain and r ≤ √1/q.
We consider the Nyquist wavenumber, that is Substituting Eq. 29 into (28), we obtain the stability condition We denote the right-hand side of inequation (30) as the stability factor where the stability factor s is related to the FD coefficients a m,n and b l,0 , and these FD coefficients are determined by the Courant number r. In the following, we analyze the stability factors s varying with the Courant number r, the stability curves of our new scheme and other three methods are shown in Figure 3. It can be seen that the stability factors of the proposed method are slightly less than that of the TE-CR-TS method. Although the DRP-CR-TS method adopts an analytical cross stencil, its stability does not sharp decrease, and it is much larger than that of the conventional methods (TE-C-S and TE-C-TS). Figure 4 shows the maximum value of stability factors satisfying r ≤ s for different orders M. The stability curve of DRP-CR-TS method is volatile due to the use of numerical method to solve the over-determined system, and it is sensitive to order M, but this does not affect the overall stability. It can be seen that the stability curve of the proposed method has the same level with the TE-CR-TS method. Stability analyses in Figures 3, 4 reveals that the proposed scheme has the same level stability as the standard temporal high-order scheme (TE-CR-TS method), and far better than the conventional scheme.

Homogeneous model
In this section, we use a 2-D homogeneous velocity model to examine our new scheme. The 2D homogeneous model has 512 × 512 grid nodes with the grid spacing h = 6m and velocity υ = 1500 m/s. A Ricker-wavelet source with a main frequency of 40Hz is located at the spatial point (1536m, 1536m). Two receivers at the spatial points (768m, 768m) and (768m, 1536m) are used to record the waveforms. Figure 5 shows the snapshots with a time step τ = 0.001s. The TE-C-S and TE-C-TS methods have serious temporal and spatial dispersion errors ( Figure 5A). The temporal dispersion error (red arrow) of the TE-CR-TS method is smaller than that of the TE-C-TS method. However, the TE-CR-TS method still has obvious spatial dispersion error (white arrow), because the TE-based FD coefficients preserve the dispersion relation in a limited range. Figure 5D shows that the corresponding spatial dispersion error is reduced considerably in the proposed method (DRP-CR-TS). Figure 6 shows the seismic records of the four methods at Receiver 1 and Receiver 2. The reference traces represented by the green curves are obtained by the high-order FD scheme under the fine grid. The Receiver 1 (Figure 6A) shows that the temporal dispersion errors of the TE-CR-TS and DRP-CR-TS methods are smaller than that of the TE-C-S and TE-C-TS methods, and the Receiver 2 shows that the spatial dispersion error are serious in the TE-based methods (TE-C-S, TE-C-TS, and TE-CR-TS). But the proposed method (DRP-CR-TS) still has a small level of the dispersion error. Table 2 lists the relative errors of the four methods compared to the reference trace at Receiver 1 and Receiver 2. The relative errors of the crossrhombus stencil are smaller than that of the conventional cross stencil, especially the DRP-CR-TS method reduces the relative error significantly.
We also analyze the snapshots for different orders M, the snapshots are shown in Figure 7. In this case, simply increasing M can not effectively reduce the dispersion error in the TE-C-S method. The temporal dispersion errors of the TE-C-TS and TE-CR-TS methods are gradually reduced, but when M = 12, there are still obvious spatial dispersion error (white arrow). The corresponding spatial dispersion error is mitigated in the proposed method (DRP-CR-TS) even at low order (M = 6). Then, we study the snapshots of the proposed method for different N and L, the results are shown in Figure 8. The dispersion error of the proposed method is small for different L and N, thus we can select an appropriate low-order L or N to improve the computational efficiency.

Inhomogeneous model 4.2.1 2D BP model
We use a widely referred 2D BP velocity model (Figure 9) to test the four methods in the inhomogeneous model. The 2D BP model has 560 × 1000 grid nodes with the variation of velocities from 1500 m/s to 4500 m/s. In this case, we set time step τ = 0.0007 s, M = 8, N = 3, L = M, grid spacing h = 7.5m and main frequency f m = 30Hz for numerical simulation. A total of 1,000 receivers are located on the surface of the model. Figure 10 shows the snapshots of the four methods. The TE-C-S method has obvious temporal dispersion error (red arrow), and the corresponding error in the TE-C-TS and TE-CR-TS methods is reduced (Figures 10B, C). However, the spatial dispersion error (white arrow) is still serious. Figure 10D shows the spatial dispersion error of the proposed method (DRP-CR-TS) is significantly reduced in the low-and high-velocity layers. Figure 11 shows the seismic records of the four methods, and Figure 12 shows the corresponding seismic records at (1920m, 0m). It can be seen that the TE-based methods have serious spatial dispersion error from the first arrivals, but the corresponding error is mitigated in the proposed method (DRP-CR-TS). The reflected wave from the high-velocity layers contain both the spatial and temporal dispersion errors, and the corresponding error in the proposed method (DRP-CR-TS) is smaller than that of the other three methods.

Discussion
In this section, we discuss the computational cost and accuracy simultaneously. Taking the above BP model as an example, we design a series of FD parameters for seismic modeling. The numerical experiments are executed on the same computer (Intel Core I7-700 with 3.6 GHz and 8 GB memory). Table 3 shows the CPU execution times and the relative errors at spatial point (1920 m, 0 m) of the four methods. It is clear that the TE-C-S and TE-C-TS methods have fewer execution times in the numerical experiments (Cases 1 and 2 in Table 3), but their relative errors are larger than that of the TE-CR-TS and DRP-CR-TS methods. The relative error of the DRP-CR-TS method is significantly reduced, and their execution time is less than that of the TE-CR-TS method. It is worth noting that when we reduce the length of the analytic FD operator in the DRP-CR-TS method (Cases 5 and 6), the execution time is reduced considerably, and the relative error is almost unaffected. Besides, the DRP-CR-TS method can select a relatively larger time step to reduce the computational cost due to the temporal high-order approximation.

Conclusion
We propose a new staggered-grid DRP-based FD scheme with a cross-rhombus stencil for solving the scalar wave equation. The new scheme has a simplified dispersion relation, which is convenient for solving the dispersion-relation-preserving FD coefficients. Besides, the simplified scheme uses the cross stencil instead of the cross-rhombus stencil in some FD operators, thus reducing the computational cost considerably. Dispersion analyses reveals that the proposed FD scheme can effectively mitigate the dispersion error, and it still has a temporal higher-order approximation accuracy. The proposed scheme also has better stability compared with the conventional scheme. Numerical experiments show that the proposed scheme has smaller temporal and spatial dispersion errors while ensuring the computational efficiency, and it is an economical way for the large-scale seismic modeling.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions
Conceptualization is contributed by CZ and GC; investigation is contributed by CZ, GC, LF, and XZ; data contributed by CZ, GC; formal analysis is contributed by CZ and GC; writing-original draft preparation is contributed by CZ, GC; writing-review and editing is contributed by CZ, LF, and XZ.
Funding their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.