Advanced finite-difference methods for seismic modeling

Yang Liu1,2 and Mrinal K Sen2

1
State Key Laboratory of Petroleum Resource and Prospecting (China University of Petroleum, Beijing), Beijing, 102249, China
2
The Institute for Geophysics, John A. and Katherine G. Jackson School of Geosciences, The University of Texas at Austin, 10100
Burnet Road, R2200 Austin, TX 78758, USA
E-mail: wliuyang@vip.sina.com, mrinal@ig.utexas.edu

Abstract

The finite-difference methods (FDMs) have been widely used in seismic modeling and migration. In this paper,
we review the conventional arbitrary-order explicit FDMs and their recent developments, including arbitrary even-order
implicit FDMs for standard grids and arbitrary even-order time-space domain FDMs for acoustic wave equations. For
a given accuracy, an arbitrary even-order FDM can provide a trade-off between the order number and grid size. These
explicit FDMs are the most popular in seismic modeling community. The finite-difference methods that are implicit in
space are not very common because they generally require more computer resources than the explicit FDMs. Here we
show that a new class of implicit FDMs can be derived that require solving a tri-diagonal system, which makes the
resulting algorithm computationally very efficient. Therefore, some high-order explicit FDM may be replaced by an
implicit FDM of some order to decrease the computation time without affecting the accuracy. We further demonstrate
that we can also derive the FD coefficients in the joint time-space domain. These high-order spatial finite-difference
stencils designed in joint time-space domain, when used in acoustic wave equation modeling, can provide even greater
accuracy than those designed in the space domain alone under the same discretization. We demonstrate performance of
these algorithms using some realistic 2D numerical examples.

Introduction variable time step (Tessmer, 2000) and implicit methods

(Emerman et al, 1982; Zhang and Zhang, 2007).
The finite-difference method (FDM), one of the most
popular methods of numerical solution of partial differential One obvious approach to improve on the accuracy
equations, has been widely used in seismic modeling (e.g., of FD calculation is to design high order operators for
Kelly et. al., 1976; Virieux, 1986; Igel et. al. 1995; Etgen, 2007; computing the derivatives (Dablain, 1986; Fornberg, 1987;
Bansal and Sen, 2008) and migration (Claerbout, 1985; Li, Crase, 1990; Liu et al, 1998; Dong, 2000; Pei, 2004; Chen,
1991; Zhang et. al., 2000; Fei et. al. 2008). The development of 2007; Liu et al, 2008; Liu and Sen, 2009). In this paper, we
an FD operator is based generally on a Taylor series review the developments of arbitrary even-order FDMs for
representation of a function (the wavefield) resulting in a seismic modeling, including EFDMs (explicit FDMs), IFDMs
formulation for numerical evaluation of the derivatives (implicit FDMs) and time-space domain FDMs.
appearing in the wave equation. The accuracy of these
operators is dependent on the order of approximation, namely, First, we briefly review the seismic wave equations,
the number of terms used in the Taylor series representation finite difference, dispersion, stability and the necessity of
of the function. Thus a discrete version of the wave equation high-order finite differences. We restrict our discussion to
is derived where the wavefield is propagated starting from acoustic wave equation because of the ease of algebra. These
the source location (the initial condition). In general, accuracy are followed by description of EFDMs, new IFDMs and time-
of the derivative calculation for a given order depends on space domain FDMs. Finally, we show some numerical
grid spacing. A small grid size helps to increase the accuracy examples to demonstrate the salient features of different
but results in a large number of grids to represent a geologic algorithms.
model. Thus, there is a substantial increase in memory and
floating point operations making these methods prohibitive Background
for realistic modeling of complex geological structures in the
seismic frequency range. Seismic wave equations: finite difference
approximations
To improve the accuracy and stability of a FDM in
numerical modeling, many methods have been developed Seismic wave equations, used to describe seismic
including difference schemes of variable grid (Wang et. al., wave propagation in the subsurface, are typically partial
1996; Hayashi et. al., 1999), irregular grid (Opral et. al., 1999), differential equations containing spatial and temporal
standard staggered grid (Virieux, 1986; Levander, 1988; derivatives. The constant density acoustic wave equation
Robertsson et al., 1994; Graves, 1996), rotated staggered grid has been used most widely in seismic modeling, migration
(Gold et al., 1997; Saenger et. al., 2004; Bohlen et. al., 2006),

GEOHORIZONS December 2009/5

and inversion. The 1D acoustic wave equation in a wavefield values at the current grid and its two neighboring
homogenous medium is given by grids, at a given time step. Such a scheme is called an explicit
scheme in space. However, unlike the temporal derivatives, a
..(1) high-order finite difference on the spatial derivatives is stable.
, A (2N)th-order finite difference formula for the 2nd-order
spatial derivative is given by
where,
, ....(5)
is a scalar wave field, and is the velocity.

From Equation (1), we can see that the wave field is where, are the finite-difference coefficients and their
a function of the spatial coordinate and time. expressions are given in the following section. When ,
we can obtain , , then Equation (5) is similar to
The following 2nd-order finite difference is usually Equation (2). When , then , , .
used for calculating temporal derivatives, Also, a smaller grid size can help attain greater accuracy for
spatial derivatives.
, ..(2)
Substituting Equations (2) and (5) into Equation (1), we have

where,
.
, ..(3) ..(6)

h and are grid size and time step respectively (see Figure 1). Rearranging Equation (6), we obtain,
For the 2nd-order continuous-time wave field,

. ..(4) ,

..(7)

where, r = v/h. Equation (7) is the recursion formula for

solving 1D wave equation by a finite-difference method. The
recursion begins with the wave field values known at two
successive time steps t=0 and t=. The wavefield values at a
future time step at all the spatial locations are computed using
Equation (7).

Dispersion

One obvious observation we can make is that

truncation of terms in a Taylor series expansion causes
inaccuracies in the evaluation of a derivative. It also causes
frequency-dependent numerical errors due to truncation of
Fig.1 Illustration of grids for wave fields in the time-space domain higher order terms, even when the material properties are not
frequency dependent. In fact, we observe a numerical
Thus, a smaller time step leads to greater accuracy for phenomenon called dispersion that is similar to what is
computing temporal derivatives. Higher order approximation observed in dissipative or attenuating media in which
of Equation (2) can also be derived from the Taylor series. propagation velocities are frequency-dependent. Frequency-
Note that Equation (2) requires the wavefield values at the dependent numerical dispersion is observed even in non-
current time step, past time step and a future time step at a attenuating media due to inadequate sampling of wavefields
given spatial point. Such a scheme is called an explicit scheme in space and time.
in time. Generally, an explicit high-order finite difference on
the temporal derivatives requires large computer memory and Here, we discuss the dispersion generated by the
is usually unstable. An equation similar to Equation (2) can recursion formula (7). Using a plane wave theory, we let
also be derived for a spatial derivative which will require the

GEOHORIZONS December 2009/6

 , ..(8) the dispersion curve of this finite-difference modeling. It
follows that high-wavenumber components have lower than
the true velocity. This is well demonstrated with Figure 2(b).
where, k is the wavenumber, is the angular frequency, and In Figure 2(b), because of their different velocities, different-
. Substituting Equation (8) into (7) and simplifying wavenumber components gradually disperse as the wave
it, we have propagates, which is the so-called dispersion phenomenon.

, ..(9)

which can be solved for a frequency ; a group velocity is

given by the ratio of frequency to wavenumber k. The
group or dispersion velocity is

..(10)
(a)
Note that, is the velocity with which the wave
propagates through numerical grids. This expression shows
that the dispersion velocity depends on the medium velocity
, grid size , time step , wavenumber and finite-different
coefficients .

Using Equation (10) and r = v/h, we define a

parameter to describe the dispersion

..(11)

Note that the parameter is the ratio of dispersion (b)

velocity to the true velocity of the medium and in fact it is a
measure of numerical error of finite difference approximation.
If equals 1, there is no dispersion. If is far from 1, a large
dispersion will occur. Thus our primary goal is to make
close to 1 by a proper choice of grid spacing for a given
velocity and frequency. Since is equal to at the Nyquist
frequency, only ranges from 0 to when calculating .

Here, we show a 1D modeling example to examine

this effect. Modeling parameters are shown in Figure 2. Initial
conditions for the modeling are

, ..(12a)
(c)

..(12b) Fig.2 Illustration of dispersion in finite-difference modeling: (a)

Snapshots from 1D analytic solution; (b) Snapshots from 1D
acoustic wave equation modeling by finite-difference recursion
where x0 is the location of the center of the source, and 2 is Formula (7), N=2. Snapshots (1) to (6) are recorded at 0.05s,
0.1s, 0.15s, 0.2s, 0.25s and 0.3sm, respectively. (c) The
an attenuation coefficient. Figures 2(a) and 2(b) show exact dispersion curve of this finite-difference modeling. The model
and modeling results with snapshots. In this 1D case, it can is homogeneous, v = 3000m/s, model size is 1000m. The
be seen from the snapshots that high-wavenumber modeling parameters are: h = 10m, = 0.001s, , .
components propagate with lower velocity. Figure 2(c) shows

GEOHORIZONS December 2009/7

Stability Necessity of high-order finite differences on the
spatial derivatives
Recursion formula, such as the one used in Equation
(7), is commonly used to perform finite-difference modeling.
In the previous section, we demonstrated the presence
If the modeling parameters are inappropriate, for example, if
of grid dispersion in the acoustic wavefield modeled by a 2nd-
time step is too large, the recursion will be unstable. In other
order finite-difference operator. There are several ways to
words, the wavefield will grow and will eventually exceed the
reduce the dispersion. The first way is to decrease the grid
machine precision.
size. Due to the limitation of computer memory, we cannot use
We will now use a conventional eigenvalue method a very small grid size, especially for 3D modeling. The second
(Chen and Guan, 1989) of stability analysis to derive the approach is to decrease the time step. Figure 4(a) shows
stability condition for recursion formula (7). Let dispersion curves of low-order finite-difference modeling for
different time steps, which suggest that decreasing the time
step will decrease the dispersion of low-wavenumber
. components, but will not decrease the dispersion of middle-
and high-wavenumber components. The third approach is to
..(13)
use high-order finite differences. Figure 4(b) shows dispersion
According to Equation (7) and (13), we can obtain curves for high-order finite-difference modeling for different
(see Appendix A) time steps, which suggest that decreasing time step will
decrease the dispersion of low- and middle-wavenumber
components. Comparing Figure 4(a) with 4(b), we conclude
, ..(14) that a high-order finite difference on the spatial derivatives is
necessary to reduce dispersion. One other approach to
suppress the grid dispersion is to use staggered grids (e.g.,
where G is a transition matrix,
Virieux, 1986; Kindelan et al 1990; Graves, 1996; Saenger and
Bohlen, 2004; Bohlen and Saenger; 2006). However, in this
. ..(15) paper we will discuss standard grids only.

When the absolute values of the transition matrix

eigenvalues are less than or equal to 1, the recursion is stable.
If , the roots of the eigenvalue equation
will be less than or equal to 1. Then, the following stability
condition can be obtained

, ..(16)

where, N1 =int [(N+1)/2] and int is a function to get the integer

part of a value. We calculate the variation of s with N, which is
(a)
shown in Figure 3. From this figure, we can see that the area of
r for stable recursion decreases with the increase of N.

(b)
Fig. 4 Dispersion curves for 1D acoustic wave equation modeling by
finite-difference recursion Formula (7) for different time steps.
Fig.3 Illustration of stability condition of 1D acoustic wave equation v = 3000m/s, h = 10m. (a) Low-order finite difference on
modeling by finite-difference recursion Formula (7). The spatial derivatives, N=2; (b) High-order finite difference on
recursion is stable when . spatial derivatives, N=10.

GEOHORIZONS December 2009/8

Explicit finite-difference methods

Centered finite-difference formulas in a standard grid ..(21)

scheme are usually used to numerically solve partial
.
differential equations, which generally involve the first-order
and the second-order derivatives in a seismic modeling
application. In the following, we introduce centered finite-
When ,
difference formulas with arbitrary even-order accuracy for
the first-order and the second-order derivatives.
. . ...(22)
Explicit finite-difference formula for the 1st-order
derivatives To examine the accuracy of these finite-difference
formulas (17) and (20), the following expressions are used
A (2N)th-order finite-difference formula for the first-order (Liu and Sen, 2009a)
derivatives can be expressed as (Dablain, 1986)

, . ...(23a)
, ............(17)

where, p(x) is a function of x, x is a real variable, h is a small , . ...(23b)

constant value of grid size, N is a positive integer; the finite-
difference coefficient cn is an odd discrete sequence. Appendix
B gives a methodology based on plane wave theory for where, , k is wavenumber. When calculating the
deriving finite-difference coefficients. cn can be computed values of the above expressions, only ranges from 0 to /2
by the following formula (Fornberg, 1987; Liu et. al., 1998, because kh is equal to at the Nyquist frequency. The
2008; Liu and Sen, 2009a) calculated values from formulas (23a) and (23b) are compared
with the exact values of for different order numbers
determined by N, shown in Figure 5. We notice that the
. ...(18)
(a)

When ,

. ...(19)

Thus, there exists a limit at which we can attain

infinite accuracy theoretically and at this limit we arrive at the
pseudospectral method (Fornberg, 1987).

Explicit finite-difference formula for the 2nd-order

(b)
derivatives

A (2N)th-order finite-difference formula for the

second-order derivatives can be expressed as

, ...(20)

where, the difference coefficient, cn an even discrete sequence,

can be computed by the following formulas (Liu et. al., 1998, Fig. 5 Accuracy of explicit finite-difference formulas
2008; Liu and Sen, 2009a) (a) The exact and f 1 ()of EFDM with different order
numbers for the first-order derivatives
(b) The exact and f2() of EFDM with different order numbers
for the second-order derivatives

GEOHORIZONS December 2009/9

accuracy increases with the increase of the order number.
, ..(27).
Implicit finite-difference methods

We noted in the previous section that an EFDM where, the difference operator in the denominator is a 2nd-
directly calculates the derivative value at some point using order centered finite difference stencil, , an odd discrete
the function values at that point and at its neighboring points. sequence, and b can be computed by solving the following
However, an IFDM (implicit FDM) expresses the derivative equations (Liu and Sen, 2009a)
value at some point in terms of the function values at that
point and at its neighboring points and the derivative values
at its neighboring points. For example, a compact finite-
difference method (CFDM) is one such IFDM (Lele 1992). ....(28)
However, IFDMs are usually considered expensive due to
the requirement of solving a larger number of equations and
are therefore not very popular. Liu and Sen (2009a) derived .
new implicit formulas for space derivatives with arbitrary
even-order accuracy for arbitrary-order derivatives. Their
approach requires solving tridiagonal matrix equations only. Implicit finite-difference formula for the 2nd-order
They also showed that a high-order explicit method may be derivatives
replaced by an implicit method of some order, which will
increase the accuracy but not the computational cost. An implicit finite-difference formula for the second-
order derivatives can be expressed as (Liu and Sen, 2009a)
We first introduce Claerbouts idea (Claerbout 1985)
upon which the methods of Liu and Sen (2009a) are developed.
The 2nd-order difference operator for a function is expressed as , .. (29).

. ..(24).

where, cn is an even discrete sequence, ,

The accuracy of finite difference can be improved
and b can be computed by solving the
by adding higher-order terms, that is,
following equations (Liu and Sen, 2009a)

, ..(25).

where b is an adjustable constant. This expression is rarely

used it can be changed to the following form (Claerbout 1985)

. ..(26). .
.. (30)

This equation has greater accuracy than equation (24). The following expressions are respectively used to
In order to use the above representation of the 2nd- order examine the accuracy of implicit finite-difference formulas
difference operator, one needs to multiply the FD equation (27) and (29) (Liu and Sen, 2009a)
through out by the denominator and then rearrange terms to
solve for the unknown.

, .. (31a)
Implicit finite-difference formula for the 1st-order
derivatives

Motivated by Formula (26), an implicit finite-

difference formula for the first-order derivatives can be . .. (31b)
expressed as (Liu and Sen, 2009a)

GEOHORIZONS December 2009/10

 The exact values of the above expressions are , respectively. The finite-difference coefficients are derived
and the calculated values are compared with for different independently for time and space derivative operators and
order numbers determined by N, shown in Figure 6. We can the approximations are incorporated into wave equations to
see that the accuracy increases with the increase of the order arrive at a finite difference based recursion operator.
number. Comparing Figure 5 with Figure 6, we note that for Alternatively, we can apply finite difference approximation
the same order number, the accuracy of the implicit method is of time and space derivatives directly to the wave equations
better than that of the corresponding explicit method. It has simultaneously and then derive the coefficients.
been found (Liu and Sen, 2009a) that for standard-grid finite
differences, the accuracy of a (2N+2)th-order implicit formula To demonstrate this, we start with 1D acoustic wave
(27) is nearly equivalent to that of a (6N+2)th-order explicit equation in homogenous media given by Equation (1). A 2nd-
formula (17) for the first-order derivative, and a (2N+2)th- order finite difference on the time derivatives, given by
order implicit formula (29) is nearly equivalent to a (4N+2)th- Equation (2) is usually used. Generally, the modeling accuracy
order explicit formula (20) for the second-order derivative. is improved by using a high-order finite difference on the
Because the computation cost of a (2N+2)th-order IFDM is spatial derivatives; a (2N)th-order finite-difference formula is
equal to that of an EFDM of some order plus that of solving given by Equation (5).
tridiagonal equations, a high-order EFDM may be replaced
by an IFDM of some order, which will increase the accuracy Substituting Equations (2) and (5) into Equation (1)
but not the cost of computation. and rearranging, we have

.
.. (32)

In the conventional method, the finite-difference

coefficients of the space derivatives are determined in the
space domain alone and can be determined by Equation (21)
to satisfy Equation (5). It can be proved that when we use the
(2N)th-order space domain finite-difference and the 2nd-order
time domain finite-difference stencils to solve the 1D acoustic
wave equation, the accuracy order is 2. Obviously, the
conclusion is the same for the 2D and 3D acoustic wave
equations. It is worthwhile to note that increasing N may
decrease the magnitude of the finite-difference error but may
(a) Exact and g1() of IFDM with different order not increases the accuracy order. The main reason is that the
numbers for the first-order derivative finite-difference stencils are designed in the space and time
domains separately, but the wave equation must be solved in
both the space and time domains simultaneously.

To address this issue, i.e., not to satisfy Equation

(5) but to satisfy Equation (32), a unified methodology
(Finkelstein and Kastner, 2007) has been proposed to derive
the finite-difference coefficients in the joint time-space
domain. The key idea of this method is that the dispersion
relation is completely satisfied at some designated
frequencies. Thus several equations are formed and the finite-
difference coefficients are obtained by solving these
equations simultaneously. It is obvious that there is no
dispersion at these designated frequencies when these
coefficients are used in modeling. However, since different
(b) Exact and g2() of IFDM with different order specific frequencies give different results and dispersion at
numbers for the second-order derivative other frequencies is not easily controlled, this method may
Fig. 6 Accuracy of implicit finite-difference formulas not be very useful in practice.
Time-space finite-difference methods Liu and Sen (2009b) developed a new method similar
to Finkelstein and Kastner (2007) which employs a plane wave
In the previous sections, we derived finite-difference
theory and the Taylor series expansion of dispersion relation
stencils of time and space derivatives by using the plane
to derive the finite-difference coefficients in the joint time-
wave theory and the truncated Taylor series in time and space
space domain for acoustic wave equation. Using the approach

GEOHORIZONS December 2009/11

of Liu and Sen (2009b), the following finite-difference
coefficients are obtained

..(33)

where, r = /h. When r = 0 , the finite-difference coefficients

are the same as Equation (21). That is, the conventional (b) New method
method is just a special case of the new method. The 1D Fig. 7 1D dispersion curves of the conventional and new methods
for different space point numbers 2N+1, N = 2, 4, 10, 20, v =
finite-difference modeling method by using these new 3000m/s, = 0.001s, h = 10m.
coefficients from Equation (33) can reach (2N)th-order
accuracy.
12(a) and 12(b) are used to perform modeling. Figure 8 shows
the records for different velocities. The dispersion of the
A parameter is defined as follows to describe the
new method with finite-difference coefficients determined by
dispersion of 1D finite-difference modeling
Equation (33) is smaller than that of the conventional method
finite-difference coefficients determined by Equation (5). The
variation of the velocity affects the results of the conventional
. ..(34)
method more than the new method because the finite-
difference coefficients of the new method depend on velocity.
If equals 1, there is no dispersion. If is far from 1, a
large dispersion will occur. Because kh is equal to at the Nyquist
frequency, kh only ranges from 0 to when calculating .

Figure 7 shows the variation of the dispersion

parameter with kh for different space point numbers. The
involved parameters are listed in the figure. This figure
demonstrates that

Dispersion decreases with the increase of N.

For the conventional method, nearly equals 1 when
. The area where nearly equals 1 does not
extend with the increase of N. That is, increasing N
decreases the magnitude of the dispersion error but (a) v = 4500m/s
does not increase the accuracy order.
For the new method, the area where nearly equals 1
obviously extends with the increase of N.
The accuracy of the new method is greater than that
of the conventional method.
As an example, the initial conditions of Equations

(b) v = 5500m/s
Fig. 8 1D modeling records by conventional EFDM and new time-
space domain EFDM for different velocities. (1), (4) and (7)
are analytic solutions; (2), (5) and (8) are modeling results
from the conventional method, (3), (6) and (9) are modeling
results from the new method. Distances between source center
and these three receivers are 100m, 350m and 600m
(a) Conventional method respectively. 2 = 2, h = 10m, = 0.001s, N = 20.

GEOHORIZONS December 2009/12

 This new time-space domain method has been
extended to solve 2D acoustic wave equation and can reach
(2N)th-order accuracy along 8 directions (Liu and Sen, 2009c).

Numerical modeling examples

To demonstrate the effects of different finite-

difference methods stated in this paper, we perform numerical
modeling of 2D acoustic wave equation for the SEG/EAGE
salt model using the following 2D acoustic wave equation

. ..(35) (b) Snapshot computed by the conventional 14th-order

EFDM
In the modeling algorithm, temporal derivatives are
discretized by second-order finite difference similar to
Equation (2). Spatial derivatives are discretized by following
three finite-difference methods discussed earlier. The first
method is a conventional 14th-order EFDM; its discretization
formula is given in Equation (20) and the finite-difference
coefficients are determined by Equations (21). The second
one is a 14th-order IFDM; its discretization formula is given
in Equation (29) and the coefficients are determined by
Equations (30). The third one is a 28th-order time-space
domain FDM; its discretization formula is given in Equation
(20) and the coefficients are determined in the joint time-
space domain (Liu and Sen, 2009c). It has been reported that (c) Snapshot computed by the 14th-order IFDM (equivalent
the 14th-order IFDM is nearly equivalent to the conventional to the conventional 28th-order EFDM)
28th-order EFDM (Liu and Sen, 2009a).

The SEG/EAGE salt model is shown in Figure 9(a);

some modeling parameters are given in the caption of Figure
9. Here, we simply extend the model spatially to avoid
reflections from the top and other edges of the model. The
modeling snapshots are shown in Figures 9(b), 9(c) and 9(d).
The modeling records are shown in Figures 9(e), 9(f) and
9(g). Partially enlarged views of these three figures are shown
in Figures 9(h), 9(i) and 9(j). Comparing these modeling results,
we can see that the accuracy of the 14th-order IFDM is greater
than that of the conventional 14th EFDM. The 28th-order
time-space domain FDM can retain the waveform better than
the conventional 28th-order EFDM or the 14th-order IFDM. (d) Snapshot computed by the 28th-order time-space
domain FDM

(a) SEG/EAGE salt model

(e) Seismograms by the conventional 14th-order EFDM

GEOHORIZONS December 2009/13

 (f) Seismograms computed by the 14th-order IFDM (j) Zoom of (g)
(equivalent to the conventional 28th-order EFDM)
Fig.9 Numerical modeling results of 2D acoustic wave equation for
the SEG/EAGE salt model by the conventional 14th EFDM,
the 14th-order IFDM (equivalent to conventional 28th-order
EFDM) and the 28th-order time-space domain FDM.
= 0.002s, h = 20m. The coordinate of the source is (6000m,
20m). A one-period sine function with 20Hz frequency is used
as the source wavelet. Time of snapshots is 1600ms.

Conclusions

The finite difference methods for seismic wavefield

modeling and imaging are considered to be the most accurate
since these methods are capable of computing all the wave
modes. In addition there is no restriction on the complexity
(g) Seismograms computed by the 28th-order time-space
of a geologic model as long as we can represent the model in
domain FDM
a discrete form. However, these methods indeed suffer from
numerical inaccuracy resulting from approximating a
derivative operator and discretization of a geologic model.
Novel finite-difference operators can be derived to address
this issue. In this paper, we have reviewed the developments
of arbitrary even-order FDMs with application to acoustic
wave equation. The accuracy of these methods increases
with the increase of order number. Since the accuracy is
determined by the order number and grid size, arbitrary even-
order FDMs can provide a trade-off between the order number
and grid size for a given accuracy. The IFDMs require solving
tridiagonal matrix equations only. Some high-order EFDM
may be replaced by an IFDM of some order to decrease the
(h) Zoom of (e) computation time and retain the accuracy. High-order spatial
FD stencils designed in the joint time-space domain can
improve the accuracy further for acoustic wave equation
modeling. This time-space domain method can be extended
to solve not only the acoustic wave equation but also other
similar partial difference equations.

Acknowledgements

Liu would like to thank China Scholarship Council for

their financial support for this research and UTIG for providing
with the facilities. This research is also partially supported by
NSFC under contract No. 40839901 and the National 863
Program of China under contract No. 2007AA06Z218.

(i) Zoom of (f)

GEOHORIZONS December 2009/14

 International Meeting, SEG, Expanded Abstracts, 523-527.
References Igel, H., P. Mora, and B. Riollet, 1995, Anistotropic wave
propagation through finite-difference grids: Geophysics, 60,
Bansal, R., and M. K. Sen, 2008, Finite-difference modelling of S- 12031216.
wave splitting in anisotropic media: Geophysical
Prospecting, 56, 293-312. Kelly, K. R., R. Ward, W. S. Treitel, and R. M. Alford, 1976,
Synthetic seismograms: A finite-difference approach:
Bohlen, T., and E. H. Saenger, 2006, Accuracy of heterogeneous Geophysics, 41, 227.
staggered-grid finite-difference modeling of Rayleigh waves:
Geophysics, 71, T109T115. Kindelan M., A. Kamel, and P. Sguazzero, 1990, On the construction
and efficiency of staggered numerical differentiators for the
Chen, J., 2007, High-order time discretizations in seismic modeling: wave equation: Geophysics, 55, 107-110.
Geophysics, 72, SM115SM122.
Lele, S. K., 1992, Compact finite difference schemes with spectral-
Chen J. F., and Guan Z., 1989, Numerical solution of partial like resolution, Journal of Computational Physics, 103, 16
differential equations, Tsinghua University Press, Beijing, 42
China
Levander, A., 1988, Fourth-order finite-difference P-SV
Claerbout, J. F., 1985, Imaging the earths interior: Blackwell seismograms: Geophysics, 53, 14251436.
Scientific Publications, Inc., Palo. Alto, CA, USA
Li, Z., 1991, Compensating finite-difference errors in 3-D migration
Crase E., 1990, High-order (space and time) finite-difference and modeling: Geophysics, 56, 1650-1660.
modeling of the elastic wave equation: 60th Annual
International Meeting, SEG, Expanded Abstracts, 987-991. Liu, Y., C. Li, and Y. Mou, 1998, Finite-difference numerical modeling
of any even-order accuracy: Oil Geophysical Prospecting
Dablain, M. A., 1986, The application of high-order differencing to (Abstract in English), 33, 1-10.
the scalar wave equation: Geophysics, 51, 5466.
Liu, Y., and X. Wei, 2008, Finite-difference numerical modeling
Dong L. G., Z T Ma, et. al., 2000, A staggered-grid high-order with even-order accuracy in two-phase anisotropic media:
difference method of one-order elastic wave equation: Chinese Applied Geophysics, 5, 107-114.
Journal of Geophysics, 43: 411-419.
Liu Y., and M. K. Sen, 2009a, A practical implicit finite-difference
Emerman, S., W. Schmidt, and R. Stephen, 1982, An implicit finite- method: examples from seismic modeling: Journal of
difference formulation. of the elastic wave equation: Geophysics and Engineering, 6, 231-249.
Geophysics, 47, 1521-1526.
Liu Y., and M. K. Sen, 2009b, A new time-space domain high-order
Etgen, J. T., and M. J. OBrien, 2007, Computational methods for finite-difference method for acoustic wave equation, CPS/
large-scale 3D acoustic finite-difference modeling: A tutorial: SEG Beijing 2009 International Geophysical Conference,
Geophysics, 72, SM223SM230. Expanded abstract, ID54

Fei, T., and C. L. Liner, 2008, Hybrid fourier finite difference 3D Liu Y., and M. K. Sen, 2009c, 2D acoustic wave equation modeling
depth migration for anisotropic media: Geophysics, 73, S27- with a new high-accuracy time-space domain finite-difference
S34. stencil, 71st EAGE Conference, Extended Abstracts, S011

Fornberg, B., 1987, The pseudospectral method - comparisons Opral, I., and J. Zahradnk, 1999, Elastic finite-difference method
with finite differences for the elastic wave equation: for irregular grids: Geophysics, 64, 240250.
Geophysics, 52, 483-501.
Pei Z., 2004, Numerical modeling using staggered-grid high order
Gold, N., S. A. Shapiro, and E. Burr, 1997, Modeling of high finite-difference of elastic wave equation on arbitrary relief
contrasts in elastic media using a modified finite difference surface: Oil Geophysical Prospecting (Abstract in English),
scheme: 68th Annual International Meeting, SEG, Expanded 39, 629-634.
Abstracts, ST 14.6.
Robertsson, J., J. Blanch, and W. Symes, 1994, Viscoelastic finite-
Finkelstein, B., and R. Kastner, 2007, Finite difference time domain difference modeling: Geophysics, 59, 14441456.
dispersion reduction schemes: Journal of Computational
Physics, 221, 422438. Saenger, E., and T. Bohlen, 2004, Finite-difference modeling of
viscoelastic and anisotropic wave propagation using the
Graves, R., 1996, Simulating seismic wave propagation in 3D elastic rotated staggered grid: Geophysics, 69, 583591.
media using staggered-grid finite differences: Bulletin of the
Seismological Society of America, 86, 10911106. Tessmer, E., 2000, Seismic finite-difference modeling with spatially
varying time steps: Geophysics, 65, 1290-1293.
Hayashi, K., and D. R. Burns, 1999, Variable grid finite-difference Virieux, J., 1986, P-SV wave propagation in heterogeneous media:
modeling including surface topography: 69th Annual Velocity stress finite difference method: Geophysics, 51,
889901.

GEOHORIZONS December 2009/15

 simplifying it, we obtain
Wang, Y. and G. T. Schuster, 1996, Finite-difference variable grid
scheme for acoustic and elastic wave equation modeling:
66th Annual International Meeting, SEG, Expanded
Abstracts, 674-677. . (A7)

Zhang, G., Y. Zhang, and H. Zhou, 2000, Helical finite-difference

schemes for 3-D. depth migration: 69th Annual International Appendix B : A methodology based on plane wave theory for
Meeting, SEG, Expanded Abstracts, 862865.
deriving finite-difference coefficients
Zhang, H., and Y. Zhang, 2007, Implicit Splitting Finite Difference
Scheme for Multi-dimensional Wave Simulation: 75th Annual Finite-difference coefficients are determined by
International Meeting, SEG, Expanded Abstracts, 2011- satisfying finite-difference formulas, such as Equations (17),
2014. (20), (27), (29) and (32). We give a methodology based on
plane wave theory for deriving finite-difference coefficients
Appendix A : Derivation of Equation (14)
with an example for Equation (17). The methodology includes
Equation (7) can be rewritten as follows
5 steps. First, according to a plane wave theory, we express
the wave field as
. (A1)
Substituting Equation (13) into (A1) and simplifying it, we .
have
Second, substitute this function into finite-difference formula
. (A2) and simplify it, this results in

Equation (13) gives .

. (A3)
Third, use the Taylor series expansion for trigonometric
Write (A2) and (A3) as matrix form function,

. (A4) .

Fourth, compare the coefficients and obtain several

Equation (13) also gives
equations,

. (A5) .

Finite-difference coefficients for 2nd-order derivative satisfy Finally, solve these equations to obtain finite-difference
coefficients .
. (A6) The detailed derivations of finite-difference
coefficients for Equations (17), (20), (27), (29) and (32) can be
Substituting Equations (A5) and (A6) into (A4) and found in some papers (Liu and Sen, 2009a, 2009b, 2009c).

GEOHORIZONS December 2009/16