## 1. Introduction

The propagation of seismic waves through the Earth’s subsurface is described by the wave equation, one of the partial differential equations (PDEs), which describe many of the fundamental natural laws. When the subsurface earth structure is complex, it is difficult to obtain the analytic results. The finite difference (FD) method is one of most widely used numerical methods for wave equation modeling because of its high efficiency, smaller memory requirement, and easy implementation [1–7].

The first application of the FD method to wave equation modeling can be possibly traced back to Alterman and Karal [1]. Alford et al. took the grid dispersion analysis for the second-order and fourth-order FD operators and stated that it is necessary to use high-order FD operators [8]. Kelly et al. further illustrated the grid dispersion, absorbing boundary condition, and other implementation aspects of the FD method [9].

Madariaga developed a staggered grid FD scheme to solve a rupture propagation problem [10]. Virieux adapted this scheme to elastic SH waves and P-SV waves in a 2D Cartesian system [2, 3]. Levander introduced a fourth-order staggered grid FD operator in the space domain to improve accuracy [11].

Grid dispersion is one of the key numerical problems affecting practical usage when utilizing the FD method. Since the traditional FD coefficient obtained in the space domain with the Taylor expansion method is only accurate for a very limited wavenumber range [4], many efforts are paid to reducing the grid dispersion with optimized FD coefficient. Yang et al. proposed the nearly analytic discrete method for wave equation and later improved this method [12, 13]. Chen proposed high-order time discretization method to reduce the dispersion caused by the temporal discretization [14, 15]. The Fourier FD was introduced by Song and Fomel with the combination of fast Fourier transform and finite difference operators [16]. Chu and Stoffa improved the FD methods with a scaled binomial windowed FD scheme that leads to more precise discrete operators [17]. Fomel et al. introduced low-rank approximation of the wave propagator matrix to reduce the cost of wave extrapolation [18].

Generally, the FD coefficients of the spatial derivative are determined only in the spatial domain. However, wave equations are solved in the temporal and spatial domains simultaneously. Finkelstein and Kastner propose a systematic design methodology for obtaining FD coefficients to reduce dispersion, which allows the exact phase velocity or (and) group velocity dispersion relationship to be satisfied at some designated frequencies in the temporal-spatial domain [19, 20]. Etgen proposed minimizing the phase velocity error using the least squares (LS) method [21]. Liu and Sen propose a new time-space domain method to determine the higher order FD coefficients for 1D, 2D, and 3D wave equations [22], and then they use this method to get the staggered grid FD coefficients [23]. Zhang and Yao proposed the use of the simulated annealing algorithm and gave an error limitation for determining the FD coefficients in the space or the time-space domain [24]. Liang et al. proposed utilizing the linear method to determine the FD coefficient in the time-space domain [25]. Ren and Liu developed a novel optimal time-space domain staggered grid FD scheme and used least squares method to get the FD coefficients [26]. Wang et al. proposed the regularized optimization method to get the staggered grid FD coefficient in the time-space domain [27]. Chen et al. used *K* space operator-based high-order staggered grid FD method to improve accuracy [28]. Yong et al. proposed using the optimized equivalent staggered grid FD method with three sets of FD coefficients to improve the simulation accuracy [29]. Compared with the traditional high-order staggered grid FD coefficient obtained by the Taylor expansion method, these methods greatly improved the accuracy with the optimized FD coefficient.

Another way to improve the accuracy and efficiency of the FD methods is using new FD stencil. Liu and Sen studied the rhombus stencil and found that it can reach high-order accuracy along all directions [30]. Liu et al. formulated an explicit time evolution scheme with high temporal accuracy by using a new FD stencil for the second-order wave equation [31]. Tan and Huang propose a staggered grid FD stencil with added points in the diagonal direction for the first-order wave equation [32, 33]. Compared with the traditional staggered grid FD stencil, these methods improved the efficiency by using a larger time step while still preserving high accuracy.

## 2. Acoustic/elastic wave equations

The first-order velocity-stress acoustic wave equation can be described as

where *P* is the acoustic pressure fluctuation, *v* is the wave propagation speed, and *v _{x}
* and

*v*are the particle velocities.

_{z}Substituting Eqs. (2) and (3) into Eq. (1), the second-order acoustic wave equation can be written as

The first-order elastic wave equations in 2D heterogeneous media are [3]

where (*v _{x}
*,

*v*) is the velocity vector,(

_{z}*τ*,

_{xx}*τ*,

_{zz}*τ*) is the stress vector, and

_{xz}*α*and

*β*are the

*P-*and

*S*wave propagation speeds, respectively.

Substituting Eqs. (7)–(9) into Eqs. (5)–(6), the second-order elastic wave equation can be written as

## 3. Finite difference operators

The commonly used staggered grid FD scheme for the first-order acoustic wave equation is as follows:

where *M*
_{1} and *M*
_{2} are the length of the FD operators, *c _{m}
* is the staggered grid FD coefficients, and

*h*is the spatial grid interval.

The second-order FD operator is usually used for the first-order time derivative:

where Δ*t* is the time step.

The commonly used staggered grid FD scheme for the first-order elastic wave equation is as follows:

where *M*
_{1} and *M*
_{2} are the length of the FD operators, *c _{m}
* is the staggered grid FD coefficient to be determined, and

*h*is the spatial grid interval.

## 4. Optimizing finite difference operators

### 4.1. Optimizing finite difference operators for the acoustic wave equation

Using the plane wave theory, let

The following dispersion relation can be obtained by substituting Eqs. (13)–(14) into Eq. (12) [23, 26, 27]:

where *r* = *v Δt*/

*h*,

*M*, and(

_{1}= M_{2}= M*k*,

_{x}*k*) =

_{z}*k*(cos

*θ*, sin

*θ*). It can be observed from Eq. (24) that the dispersion relation is complex and optimized methods are needed to address this problem.

Let *c* be the vector form of the FD coefficients, and denote the left side of Eq. 24 by [27]:

and the right side of Eq. 24 by

The aim is to minimize the dispersion error for a fixed range of wavenumbers:

The upper limit of the wavenumber range used for calculating the FD coefficients is based on the source frequency, the space grid interval, and the wave propagation speed [25]:

The direct minimization of the objective function *Φ* for the FD coefficient may lead to unstable results. Therefore, regularizing technique was applied to restore stability. The regularization model is established as

where *α* > 0 is a user-defined regularization parameter and *D* is a scale operator. The new task is the minimization of *J(c)*, and then the regularized optimization staggered grid FD coefficient can be obtained.

Another way to improve the efficiency and accuracy of the staggered grid FD methods is the utilization of the new staggered grid FD scheme. Different with the previous staggered grid FD scheme, the simplest centered second-order staggered grid FD operator can be used for the spatial derivatives in Eqs. (2) and (3), for example,

The staggered grid FD scheme in Eqs. (30)–(32) can be seen as a new staggered grid FD scheme for the first-order acoustic wave equation. The new staggered grid FD scheme is exactly the same as the traditional staggered grid FD scheme except if the staggered grid FD operator length is shorter for Eqs. (31) and (32). By carefully comparing Eqs. (12) and (14) with Eqs. (30) and (32), we find that the new staggered grid FD scheme is more efficient and can save about 45% of simulation time when *M* equals 7. It looks like the particle velocities *v _{x}
* and

*v*in Eqs. (31) and (32) are inaccurate since only the second-order staggered grid FD operators are used. However, this is not true since the staggered grid FD coefficient in Eq. (30) is optimized with Eqs. (31) and (32) in consideration. In the following, the huge advantage of the new staggered grid FD scheme will be demonstrated because it can reduce the simulation time while still preserving high accuracy compared with the traditional staggered grid FD scheme.

_{z}To get the staggered grid FD coefficient in Eq. (30), we substitute Eqs. (31) and (32) into Eq. (30), using the plane wave theory. Then we get

(33) |

From Eq. (33), the following dispersion relation can be obtained in the frequency-wavenumber domain (it is a special case of Eq. (24)):

Using the basic trigonometric function

we obtain Eq. (36) from Eq. (34):

Similarly, the new dispersion relation for the 3D first-order acoustic wave equation in the frequency-wavenumber domain is

where

Compared with the traditional dispersion relation in Eq. (24), the new dispersion relation in Eqs. (36) and (37) is linear and much simpler.

We assume that there are *M* equally distributed wavenumber points satisfying the dispersion relation within the wavenumber range specified by Eq. (28). Then, we establish the linear equation from Eq. (37) for the 3D case [25]:

(38) |

where
*i*th component of *k _{l}
*(

*l*=

*x*,

*y*,

*z*) is represented as

*k*

_{ i, l },

*k*=

_{x}*k*cos

*θ*cos

*φ*,

*k*=

_{y}*k*cos

*θ*sin

*φ*,

*k*=

_{z}*k*sin

*θ*, and

*k*(

*i*) for each

*i =*1,2,…,

*M +*1 is equally distributed between 0 and

*Ratio*×

*π*/

*h*, where

*Ratio*is determined by Eq. (28). In the following, we will demonstrate that the new staggered grid FD scheme in Eqs. (30)–(32) has similar accuracy compared with the computational intensive traditional staggered grid FD scheme in Eqs. (12)–(14).

The 2D dispersion error *δ* of the new staggered grid FD scheme is defined as

where

The difference between the FD propagation time and the exact propagation time through one grid is defined as [23]

**Figures 1**
and
**2**
show the dispersion error curves of the traditional and the new staggered grid FD schemes for the homogeneous acoustic model in 2D. All the FD coefficients are determined in the time-space domain with M = 7. From
**Figures 1**
and
**2**
, we get the conclusion that the new staggered grid FD scheme can also preserve the dispersion relation in a pretty wider range compared with the traditional staggered grid FD methods. For example, with *r* = 0.0075 in the 2D case, both of them can preserve the dispersion error under 10^{−5} within 80% of *kh* range. However, the new staggered grid FD scheme saves wave equation simulation time because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

Let the left part of Eq. (36) as.

From dispersion relation Eq. (36), it is obvious that.

Then, the stability condition of the new staggered grid FD scheme is (from Eq. (36) with *kh* = *π*).

**Figure 3**
shows the stability condition of the traditional and the new staggered grid FD scheme in 2D. We can see that the stability condition becomes stricter with the increase of the FD operator length. It also shows that the new staggered grid FD scheme’s stability condition is a little bit better than the previous staggered grid FD scheme. For example, the stability conditions are *r* < 0.54 and *r* < 0.57, respectively, for the traditional and the new staggered grid FD scheme with M = 7.

### 4.2. Optimizing finite difference operators for the elastic wave equation

Eqs. (10) and (11) can be written as [11].

The two roots give the following dispersion relation [11]:

Suppose *D _{
xx
}D_{zz} = D_{
xz
}D_{
xz
}
*, then two equations can be obtained from Eq. (46):

Usually, Eq. (48) is used to determine the FD coefficient. For the first-order staggered grid FD scheme, the following dispersion relation can be obtained from Eq. (48):

where *r* = *β Δt*/

*h*,

*M*, and (

_{1}= M_{2}= M*k*,

_{x}*k*) =

_{z}*k*(cos

*θ*, sin

*θ*). It can be observed from Eq. (49) that the dispersion relation is nonlinear and regularized optimized methods can address this problem similarly.

Different with previous staggered grid FD scheme for the first-order elastic wave equation, the simplest centered second-order staggered grid FD operator can be used for the spatial derivatives in Eqs. (7)–(9):

The staggered grid FD scheme in Eqs. (52)–(54) is more efficient than the staggered grid FD scheme in Eqs. (19)–(21). It will be demonstrated later that the staggered grid FD scheme in Eqs. (52)–(54) is accurate for the stress vector (*τ _{xx}
*,

*τ*,

_{zz}*τ*) even when only second-order staggered grid FD operator is used.

_{xz}Then, the new dispersion relation can be obtained from Eq. (49) in the frequency-wavenumber domain:

The staggered grid FD coefficient can be obtained similarly using the linear method.

## 5. Experiments

### 5.1. Acoustic wave equation

#### 5.1.1. Numerical modeling in the layered velocity model

We first consider a layered velocity model. The velocity is 1500 m/s for the first layer and 2500 m/s for the second layer as shown in
**Figure 4**
. The sponge boundary code in CREWES Toolbox is used to reduce artificial reflection waves [34]. A Ricker wavelet with the main frequency as 14.3 Hz was used as the seismic source. The seismic source position is denoted as a asterisk, and the receivers A and B are denoted as a circle and a diamond from top to bottom, respectively. The space grid interval is 20m, the FD operator length *M* is 7 and the time step is 1.5 ms. The staggered grid FD coefficients used in
**Figure 4**
are shown in
**Table 1**
.

The seismograms recorded at positions A and B by different methods are presented in
**Figure 5**
.
**Figure 5(a)**
is obtained with the traditional staggered grid FD scheme with the FD coefficient obtained in the space domain by Taylor expansion method [17]. The grid dispersion is obvious.
**Figure 5(b)**
is obtained with the traditional staggered grid FD scheme with the coefficient obtained in the space domain by the least squares method [35]. The staggered grid FD coefficient provided by Liu is one of the best staggered grid FD coefficient provided in recent years [35].
**Figure 5(c)**
is obtained with the new staggered grid FD scheme with the coefficient determined in the time-space domain by the linear method.
**Figure 5(d)**
is obtained with the pseudo-spectrum method with the second-order acoustic wave equation. We observe that the grid dispersion in
**Figure 5(b)**
and
**(c)**
is similar to each other and is close to the nearly analytic results obtained with the pseudo-spectrum method in
**Figure 5(d)**
. However, the required simulation time is reduced by using the new staggered grid FD scheme because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

#### 5.1.2. Numerical modeling in the salt model

**Figure 6**
shows the salt model from Society of Exploration of geophysicists with variations of velocities from 1486 to 4790 m/s. The seismic source function is the same as the previous example. The spatial sampling interval is 20 m, temporal step is 1 ms, and *M* = 7 for the staggered grid FD operators in
**Figure 7(a)**
and
**(b)**
. In
**Figure 7(c)**
, the parameters are *M* = 7 for the spatial derivatives in Eq. (1), and *M* = 1 for the spatial derivatives in Eqs. (2) and (3). The pseudo-spectrum method is used for the second-order acoustic wave equation as shown
**Figure 7(d)**
.

**Figure 7(a)**
is obtained the with the traditional staggered grid FD scheme with the coefficient obtained in the space domain by Taylor expansion method. The grid dispersion is obvious.
**Figure 7(b)**
is obtained with the traditional staggered grid FD scheme with the coefficient obtained in the time-space domain by the least squares method [27]. Most of the grid dispersion is suppressed.
**Figure 7(c)**
is obtained the with the new staggered grid FD scheme with the coefficient obtained in the time-space domain by the linear method. The grid dispersion in
**Figure 7(c)**
is very similar to the grid dispersion in
**Figure 7(b)**
. However, the simulation time to get
**Figure 7(c)**
is reduced compared with the simulation time to get
**Figure 7(b)**
. Both seismic records in
**Figure 7(b)**
and
**(c)**
are close to seismic record in
**Figure 7(d)**
. We want to mention that the linear method is faster than the LS method to determine the FD coefficients.

**Figure 8**
Further compares the seismograms in
**Figure 7**
at position x/dx = 400. It is also observed that with the coefficient obtained in the space domain by Taylor expansion method, the grid dispersion is serious in the simulation result. The simulation results are almost overlapped for the traditional staggered grid FD scheme and new staggered grid FD scheme with optimized FD coefficient. However, the required simulation time is reduced by using the new staggered grid FD scheme because Eqs. (31) and (32) are much simpler than Eqs. (13) and (14).

**Figure 9**
compares snapshots of particle velocity *v _{x}
* with the different staggered grid FD schemes at 2500 ms. it is also observed that with the coefficient obtained in the space domain by Taylor expansion method, the grid dispersion is most serious. The grid dispersion in

**Figure 9(c)**is very similar to the grid dispersion in

**Figure 9(b)**. It demonstrated that the new staggered grid FD scheme is accurate for the particle velocities in Eqs. (32) and (33) even when only second-order staggered grid FD operator is used.

### 5.2. Elastic wave equation

#### 5.2.1. Numerical modeling in the homogeneous media

We first consider a homogeneous model. The P wave propagation speed is 2598 m/s, and the S wave velocity is 1500 m/s. The seismic source position is at the center of the model. The grid space interval is 20 m, the time step is 1 ms, and the operator length *M* is 7. A Ricker wavelet with the main frequency as 14.3 Hz was used as the seismic source.

The snapshots of the horizontal component obtained by different staggered grid FD methods are presented in
**Figure 10(a)–(c)**
.
**Figure 10(a)**
is obtained with the traditional staggered grid FD scheme with the traditional FD coefficient. The grid dispersion is obvious.
**Figure 10(b)**
is obtained with the traditional staggered grid FD scheme with the new FD coefficient. Compared with
**Figure 10(a)**
, the grid dispersion is suppressed.
**Figure 10(c)**
is obtained with the new staggered grid FD scheme. The grid dispersion curves in
**Figure 10(b)**
and
**(c)**
are very similar, which is further demonstrated in
**Figure 10(d)**
. However, with the new staggered grid FD scheme, we can save about 45% of the modeling time.

#### 5.2.2. Numerical modeling in the homogeneous media

**Figure 11**
shows the salt model from Society of Exploration of geophysicists. The S wave velocity is obtained from the P wave velocity. The seismic source function is plotted as a red asterisk. The spatial sampling interval is 12.5 m, the temporal step is 1 ms, and *M* = 7 for staggered grid FD operators.

**Figure 12**
displays the seismic records of the horizontal component obtained by different staggered grid FD methods.
**Figure 12(a)**
is obtained with the traditional FD scheme with the traditional staggered grid FD coefficient. The grid dispersion is severe.
**Figure 12(b)**
is obtained with the traditional FD scheme with the staggered grid FD coefficient obtained by the least squares method.
**Figure 12(c)**
is obtained the with the new FD scheme with the staggered grid FD coefficient obtained by the linear method. It is observed that the grid dispersion in
**Figure 12(b)**
and
**12(c)**
is smaller than the grid dispersion in
**Figure 12(a)**
.
**Figure 12(d)**
is seismograms obtained from
**Figure 12(a)–(c)**
. It further demonstrated that the grid dispersion in
**Figure 12(b)**
and
**(c)**
is similar to each other and smaller than the grid dispersion in
**Figure 12(a)**
. However, with the new FD scheme, the simulation time is reduced about 45%. In our simulation, there are 525 grids in the z direction and 850 grids in the x direction. With the traditional FD scheme, the simulation time is 920 seconds. With the new staggered FD grid scheme, the simulation time is 530 seconds. The huge reduction in simulation time is due to using the shorter staggered FD operator for the spatial derivatives in Eqs. (7)–(9).
**Figure 13**
is the seismic records of the vertical component obtained by different FD methods. The same pattern can be observed from
**Figure 13(a)–(d)**
.

## 6. Discussion and conclusion

The FD method is the most commonly used numerical method for wave equation modeling. Suppressing the grid dispersion is an important research area. Optimization method is usually used to determine the FD coefficients which could preserve the dispersion relation in a wider range of wavenumber (Zhang and Yao [24]; Ren and Liu [26]; Tan and Huang [32, 33]). We introduced the regularized optimization method to determine the FD coefficient which would be more robust for extreme conditions. The other way to suppress the grid dispersion is the utilization of the new FD scheme for the spatial derivatives. We introduce to use different FD operators for different spatial derivatives in the first-order wave equation. With the new staggered grid FD scheme, the wave equation modeling speed was accelerated while still preserving high accuracy. Through numerical modeling, we conclude that the introduced methods are more efficient while still preserving high accuracy for the first-order acoustic/elastic wave equation modeling. As a result, the introduced methods can be a substitute for the traditional FD methods used in acoustic/elastic wave equation modeling, which are essential in forward seismic wave modeling and reverse-time migration.