ADVANCES IN GEOPHYSICS, VOL.

48, CHAPTER 8

THE FINITE-DIFFERENCE TIME-DOMAIN

METHOD FOR MODELING OF SEISMIC
WAVE PROPAGATION
P ETER M OCZO1 , J OHAN O.A. ROBERTSSON2 AND L EO E ISNER3
1 Faculty of Mathematics, Physics and Informatics, Comenius University, Bratislava, Slovak Republic
2 Western Geco Oslo Technology Center, Schlumberger House, Asker, Norway
3 Schlumberger Cambridge Research Ltd., Cambridge, United Kingdom

A BSTRACT

We present a review of the recent development in finite-difference time-domain

modeling of seismic wave propagation and earthquake motion. The finite-difference
method is a robust numerical method applicable to structurally complex media.
Due to its relative accuracy and computational efficiency it is the dominant method
in modeling earthquake motion and it also is becoming increasingly more impor-
tant in the seismic industry and for structural modeling. We first introduce basic
formulations and properties of the finite-difference schemes including promising
recent advances. Then we address important topics as material discontinuities,
realistic attenuation, anisotropy, the planar free surface boundary condition, free-
surface topography, wavefield excitation (including earthquake source dynamics),
non-reflecting boundaries, and memory optimization and parallelization.

Keywords: Anisotropy, Attenuation, Earthquake motion, Earthquake source dy-

namics, Finite-difference method, Free surface, Non-reflecting boundaries, Numer-
ical modeling, Optimally accurate operators, Seismic waves

1. I NTRODUCTION

Faithfully synthesizing observed seismic data requires simulation of seismic

wave propagation in realistic computational models which can include anisotropic
media, non-planar interfaces between layers and blocks, velocity/density/quality-
factor gradients inside layers, and often with free-surface topography. In particu-
lar, the rheology of the medium should allow for realistic attenuation. Anisotropy
increasingly becomes more important particularly in structural studies and is es-
sential in many applications in seismic exploration. The modeling of large earth-
quakes requires at least kinematic modeling of the rupture propagation but the
dynamic modeling is likely to provide more realistic simulations in the near fu-
ture. Since analytical methods do not provide solutions for realistic (structurally

421 © 2007 Elsevier Inc. All rights reserved

ISSN: 0065-2687
DOI 10.1016/S0065-2687(06)48008-0
422 MOCZO ET AL.

complex) models, approximate methods are necessary. Among them, the finite-
difference (FD) method is still the dominant method in modeling earthquake
motion and it also is becoming increasingly more important in the seismic in-
dustry and structural modeling. This is because the FD method can handle rela-
tively complex models, provides the “complete” solution as waves interact with
the model, can be relatively accurate, and, at the same time, relatively compu-
tationally efficient. In addition to this, the FD method can be relatively easily
parallelized. It is important to stress the word relatively—this chapter should ex-
plain why.
In the FD method, a computational domain is covered by a space–time grid, that
is, by a set of discrete grid positions in space and time. The functions describing
a wavefield as well as those describing material properties of the medium are rep-
resented by their values at the grid positions. In principle, the space–time grid
may be arbitrary and usually no assumption is made about the function values in-
between the grid points. Spatial and time derivatives of a function at a given grid
position are approximated by the so-called FD formulae, the derivative being ex-
pressed using the function values at a specified set of the grid positions in a neigh-
borhood of the given position. The original differential equation is thus replaced
by a system of algebraic (FD) equations. The system of FD equations and their
numerical solution have three basic properties—consistency of the FD equations
with the original differential equations, stability and convergence of the numerical
solution. These properties have to be analyzed prior to the numerical calculation.
In principle, the FD method can be applied either in the time or frequency do-
main. While for the forward modeling the time-domain formulation requires less
calculations (at least for simple simulations), the frequency-domain formulation
may be more efficient in inverse problems when simulations for multiple source
locations are required, at least in 2-D (Pratt, 1990; Pratt et al., 1998). As most
of the FD applications to date focus on the forward modeling, this chapter only
addresses the FD time-domain (FDTD) formulation.
The FD method belongs to the domain methods together with, for example,
the finite-element, spectral element or the pseudospectral method (Chaljub et
al., 2006). Therefore, in general, it is less accurate than boundary methods but
much more efficient when applied to complex models (for comparison see, e.g.,
Takenaka et al., 1998).
Because, formally, the FD method almost always can provide some numerical
results, those who are not familiar with the properties of the FD method and par-
ticularly with its inherent limitations, sometimes overestimate the capability and
accuracy of the FD method, and especially some simple, user-friendly looking FD
schemes and codes. An improper application of the FD method can give very in-
accurate results. On the other hand, when properly treated, the FD method shows
accuracy similar to that of other full waveform techniques, and can be a strong
modeling tool applicable to many important problems of recent seismology and
seismic prospecting.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 423

In this chapter we try to provide a limited review of the recent progress in

FDTD modeling of seismic wave propagation and earthquake motion as well as
partial tutorial describing in detail specific FDTD techniques that can be used to
solve real practical problems. Because explicit, heterogeneous, staggered-grid FD
schemes clearly have dominated in the recent FDTD modeling we focus on them.
In the following subsections we briefly introduce basic concepts and properties
of the FD method, some of them directly in relation to solving the elastody-
namic equation. A more detailed introduction to the FD method would require
considerably large space. Basics of the FD method can be found in applied math-
ematical textbooks such as, for example, Forsythe and Wasow (1960), Isaacson
and Keller (1966), Richtmyer and Morton (1967), Marchuk (1982), Anderson
et al. (1984), Mitchell and Griffiths (1994), Morton and Mayers (1994), Durran
(1999), Cohen (2002). Detailed introductory texts on the application of the FD
method to seismic wave propagation and seismic motion modeling can be found
in Boore (1972), Levander (1989), Moczo (1998), Carcione et al. (2002), and
Moczo et al. (2004b). Though focused on the computational electrodynamics, the
extensive book of Taflove and Hagness (2005) is a good reference to the applica-
tion of the FDTD method to partial differential equations in physics.

1.1. The Grid

Consider a Cartesian coordinate system (x, y, z) and a computational do-
main in the four-dimensional space of variables (x, y, z, t) with t meaning time.
Consider a set of discrete points (xI , yJ , zK , tm ) given by xI = x0 + I x,
yJ = y0 + J y, zK = z0 + Kz, tm = t0 + mt; I, J, K, m = 0, 1, 2, . . . .
The spatial increments x, y and z are usually referred to grid spacings, while
t is the time step. The set of points (positions) defines a space–time grid. In
many applications, the regular (uniform) rectangular grid with the grid spacings
x = y = z = h is a natural and reasonable choice. The value of a function u
at a grid position (xI , yJ , zK , tm ), that is u(I, J, K, m) or um
I,J,K , is approximated
m
by a grid function UI,J,K = U (xI , yJ , zK , tm ).
Depending on the particular problem, other than Cartesian coordinate systems
can be used to define a grid. For instance, spherical coordinates can be used for
the whole Earth’s models, and cylindrical for modeling boreholes. The choice
of the grid determines the structure and properties of the FD approximations to
derivatives and consequently the properties of the FD equations. Here we focus
on FD schemes constructed for grids corresponding to the Cartesian coordinate
systems.
In some problems it may be advantageous to define a non-uniform grid. Exam-
ples are grids with irregularly varying size of the grid spacing or discontinuous
(combined) grids with a sudden change in size of the grid spacing. Such grids can
better accommodate geometry of the model or reduce the total number of grid
points covering the computational space. These grids belong to the structured
424 MOCZO ET AL.

grids: at a grid point the neighbor grid points are always known (for example,
they are defined using some mathematical rule). On the other hand, at a grid point
of the unstructured grid some additional information is needed about the neighbor
grid points. Most FD techniques use structured grids because the algorithms are
faster than those on the unstructured grids.
Another and perhaps more important aspect is whether all functions are ap-
proximated at the same grid position or not. In a conventional grid, all functions
are approximated at the same grid positions. In a partly-staggered grid, displace-
ment or particle-velocity components are located at one grid position whereas
the stress-tensor components are located at another grid position. In a staggered
grid, each displacement and/or particle-velocity component and each shear stress-
tensor component has its own grid position. The normal stress-tensor components
share another grid position. The staggered distribution of quantities in space is re-
lated (through the equation of motion) to the staggered distribution of quantities in
time. In all types of grids, an effective density is assigned to a grid position of each
displacement or particle-velocity component while an effective elastic modulus is
assigned to each grid position of the stress-tensor components. The so-called grid
cells of the conventional, partly-staggered and staggered grids are illustrated in
Fig. 1.

1.2. The FD Approximations to Derivatives

Let the function (x) have a continuous first derivative. The forward-difference
formula
d . (x0 + h) − (x0 )
(x0 ) = , (1)
dx h
the backward-difference formula
d . (x0 ) − (x0 − h)
(x0 ) = , (2)
dx h
and the central-difference formula
d . (x0 + h) − (x0 − h)
(x0 ) = (3)
dx 2h
are three different approximations to the 1st derivative of function (x0 ). Sub-
stituting Taylor expansions of functional values (x0 + h) and (x0 − h) about
the point x0 in Eqs. (1) and (2) shows that the difference between the first deriv-
ative and the value of the right-hand side expression, that is, the truncation error,
has the leading term proportional to h1 . The FD formulae (1) and (2) are the
1st-order approximations to the first derivative. Similarly, it is easy to check that
the FD formula (3) is the 2nd-order approximation to the first derivative because
the truncation error is proportional to h2 . For a chosen derivative, set of the grid
points and order of approximation, it is possible to find a FD formula by construct-
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 425

F IG . 1. Spatial grid cells in the conventional, partly-staggered and staggered grids. All displace-
ment-vector components U , V and W are located at each grid position in the conventional grid.
Either displacement or particle-velocity components U , V and W share the same grid positions
whereas stress-tensor components Txx , Tyy , Tzz , Txy , Tyz , and Tzx share other grid positions in the
partly-staggered grid. Displacement and/or particle-velocity components U , V and W are located at
different grid positions as well as stress-tensor components Txx , Tyy , Tzz , Txy , Tyz , and Tzx are in
the staggered grid. Because the normal stress-tensor components are determined by the same spatial
derivatives of the displacement components, they share one grid position.

ing a system of algebraic equations based on Taylor expansions and equating

the coefficients of identical powers of the grid spacing h (e.g., Durran, 1999;
Moczo et al., 2004b). This is also true for approximating a derivative in a plane
or in a volume.
A frequently used 2nd-order approximation to the second derivative is
d2  . (x0 − h) − 2(x0 ) + (x0 + h)
(x0 ) = . (4)
dx 2 h2
The approximation is used in the conventional-grid displacement FD schemes.
An important 4th-order approximation to the 1st derivative is
     
d . 1 3 3
(x0 ) = a  x0 + h −  x0 − h
dx h 2 2
    
1 1
+ b  x0 + h −  x0 − h (5)
2 2
426 MOCZO ET AL.

with a = −1/24, b = 9/8. The approximation is used in the staggered-grid FD

schemes.
For other FD approximations, higher-order FD approximations and FD ap-
proximations on arbitrary spaced grids see also Anderson et al. (1984), Dablain
(1986), Fornberg (1988), Geller and Takeuchi (1995, 1998), Klimeš (1996), and
Cohen (2002).

1.3. Basic Properties of the FD Equations and their Solution

Denote a partial differential equation as PDE and a FD equation(s) as FDE
(a FD scheme may be used instead of FDE). A FDE is consistent with the PDE
if the difference between the FDE and the PDE (the truncation error) vanishes as
the sizes of the time step and spatial grid spacing go to zero independently, that
is, |PDE − FDE| → 0 if t → 0, h → 0. If this is true only when a certain
relationship is satisfied between t and h, the FDE is conditionally consistent.
A FDE is stable if it produces a bounded solution when the exact solution is
bounded, and is unstable if it produces an unbounded solution when the exact
solution is bounded. If the solution of the FDE is bounded for all values of t
and h, the FDE is unconditionally stable. If the solution of the FDE is bounded
only for certain values of t and h, the FDE is conditionally stable. If the solution
of the FDE is unbounded for all values of t and h, the FDE is unconditionally
unstable. The stability analysis can be performed only for linear PDE. A nonlin-
ear PDE must be first linearized locally. The FDE of the linearized PDE can be
analyzed for stability. The most commonly used method for the stability analysis
is the von Neumann method. The basic idea of the von Neumann method is to
represent a discrete solution at a time mt and spatial point I h, that is at one grid
point, by a finite Fourier series, and examine stability of the individual Fourier
components. Thus, the method investigates the local stability. The discrete so-
lution is stable if and only if each Fourier component is stable. Von Neumann
analysis is applicable to linear FDE with constant coefficients. Though a spatial
periodicity is assumed for the finite Fourier series, the analysis can give a useful
result even if this is not the case.
A FDE is convergent if the solution of the FDE approaches the exact solution of
the PDE as the sizes of the time step and spatial grid spacing go to zero. Denoting
the solutions obtained by the PDE and FDE as um m
I,J,K and UI,J,K , respectively,
the convergence means that UI,J,K → uI,J,K if t → 0, h → 0.
m m

Whereas the consistency is the property of the FDE because it relates the FDE
to the PDE, stability and convergence are properties of the numerical solution of
the FDE. In general, while it is easy to analyze the consistency, proving conver-
gence can be a very difficult mathematical problem. Therefore, it is very helpful
that the convergence is related to the consistency and the stability: It follows from
the Lax equivalence theorem that if the FDE is consistent and stable, it is also
convergent.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 427

Due to the discrete nature of the FD solution, the phase and group velocity in
the grid differ from the true velocities in the medium. The grid velocities depend
on the spatial sampling ratio s = h/λ, where λ is the wavelength that is to be prop-
agated in the grid, and also on Courant number ct/ h, where c is the velocity.
The grid dispersion is a very important grid phenomenon. It has a cumulative ef-
fect on the wave propagation—the longer the travel distance, the larger the effect
of the difference between the grid and true velocity. Therefore, the grid dispersion
has to be analyzed prior to the numerical calculations. This is especially important
if the medium is viscoelastic. The viscoelastic medium is intrinsically dispersive
and thus a possible superposition of two dispersion effects has to be prevented by
minimizing the grid dispersion; see Robertsson et al. (1994). The grid-dispersion
relation can be obtained from the stability analysis. The grid dispersion has to be
taken into account in planning the numerical calculation. For a detailed analysis
of stability, grid dispersion and accuracy of the FD schemes solving the equation
of motion on the conventional and staggered-grid schemes in 2-D and 3-D prob-
lems in homogeneous media see, for example, papers by (Marfurt (1984); Crase
et al. (1992); Igel et al. (1995); Geller and Takeuchi (1995, 1998); Klimeš (1996);
Takeuchi and Geller (2000); Mizutani et al. (2000); Moczo et al. (2000)).

1.4. Explicit and Implicit FD Schemes

In an explicit scheme, the motion at any (one) spatial grid point can be updated
for the next time level using an explicit FD formula which uses only values of
motion at previous time levels (and, obviously, using also material grid parame-
ters). In the case of an implicit scheme, there is no explicit formula for updating
motion only in one grid point. In an implicit scheme, the motion at a given time
level is calculated simultaneously at all spatial grid points from the motion val-
ues at previous time levels using the inverse of a matrix. Obviously, the explicit
schemes are computationally simpler. A vast majority of earthquake ground mo-
tion modeling and exploration seismology studies use explicit FD schemes. For
the implicit schemes see, e.g., Emerman et al. (1982), Mufti (1985).

1.5. Homogeneous and Heterogeneous FD Schemes

Motion in a smoothly heterogeneous elastic or viscoelastic continuum is gov-
erned by the equation of motion. The equation can be solved using a proper FD
scheme. If the medium contains a material discontinuity, i.e., an interface between
two homogeneous or smoothly heterogeneous media, at which density and elastic
moduli change discontinuously, the equation of motion still governs motion away
from the discontinuity but boundary conditions apply at the discontinuity. There
are two approaches to deal with this situation—the homogeneous and heteroge-
neous approaches. In the homogeneous approach, a FD scheme for the smoothly
428 MOCZO ET AL.

heterogeneous medium is applied at grid points away from the material discon-
tinuity while a FD scheme obtained by a proper discretization of the boundary
conditions is applied at grid points at or near the material discontinuity. In the al-
ternative heterogeneous approach only one FD scheme is used for all interior grid
points (points not lying on boundaries of a grid) no matter what their positions
are with respect to the material discontinuity. The presence of the material dis-
continuity is accounted for only by assigning appropriate values of elastic moduli
and density. Therefore, except for treating the free surface, the heterogeneous ap-
proach has been much more popular since the beginning of the seventies.
A homogeneous FD scheme is specific for a particular problem. Its application
to complex models with curved material discontinuities is difficult and therefore
impractical. Moreover, finding a stable and sufficiently accurate FD approxima-
tion to the boundary conditions is not a trivial problem, see, e.g., Kummer and
Behle (1982), Slawinski and Krebes (2002).
While widely used, the heterogeneous approach has not been addressed prop-
erly until very recently. The point is that the heterogeneous approach can be
justified only when a heterogeneous formulation of the equation of motion and
Hooke’s law, that is, the same form of the equations for both a point away from
the material discontinuity and a point at the material discontinuity, can be found.
This question will be analyzed in the next section.

1.6. The Equation of Motion, Hooke’s Law, and FD Schemes

Given a model, the equation of motion and Hooke’s law (the stress-strain re-
lation, constitutive law) together with the initial and boundary conditions fully
describe a problem of seismic wave propagation and motion. Consider a Carte-
sian coordinate system (x1 , x2 , x3 ), for example, with the x1 -axis horizontal and
positive to the right, and the x3 -axis positive downward. Let ρ(xi ); i ∈ {1, 2, 3} be
density, cij kl (xq ) tensor of elastic coefficients, κ(xi ) bulk modulus, μ(xi ) shear
modulus, u(xi , t) displacement vector, v(xi , t) particle-velocity vector, f(xi , t)
body force per unit volume, σij (xk , t) and εij (xk , t); i, j, k ∈ {1, 2, 3} stress- and
strain-tensors (from now on, x1 , x2 , x3 and x, y, z may be used interchangeably;
similarly, 1, 2, 3 and x, y, z in the subscripts of the displacement and stress-tensor
components), δij Kronecker delta. A partial time derivative will be denoted by
a dot above the symbol or the operator ∂t ; for example, ∂ui /∂t = u̇i = ∂t ui .
A derivative with respect to the coordinate xj will be denoted by a comma fol-
lowed by xj or the operator ∂j ; for example, ∂σij /∂xj = σij,j = ∂j σij . The
Einstein summation convention for repeated indices will be assumed unless stated
otherwise.
With reference to the FD schemes developed during the last decades, the fol-
lowing alternative formulations of the equation of motion and Hooke’s law for
anisotropic or isotropic media can be used as a starting point for deriving the FD
schemes:
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 429

displacement-stress

ρ üi = σij,j + fi ,
 
1
σij = cij kl εkl or σij = κεkk δij + 2μ εij − εkk δij , (6)
3
displacement-velocity-stress

ρ v̇i = σij,j + fi , vi = u̇i ,

 
1
σij = cij kl εkl or σij = κεkk δij + 2μ εij − εkk δij , (7)
3
velocity-stress

ρ v̇i = σij,j + fi ,
 
1
σ̇ij = cij kl ε̇kl or σ̇ij = κ ε̇kk δij + 2μ ε̇ij − ε̇kk δij , (8)
3
displacement
  
2
ρ üi = κ − μ uk,k + (μui,j ),j + (μuj,i ),j + fi or
3 ,i
ρ üi = (cij kl uk,l ),j + fi . (9)
Because it is not differentiated with respect to the spatial coordinates, the strain
tensor εij = 12 (ui,j +uj,i ) was used here in the first three formulations for brevity.
The above are so-called strong formulations of the equation of motion. How-
ever, a FD method can also be applied to a (integral) weak form of the equa-
tion of motion. For example, a weak form of the Galerkin type is typical
for the standard finite-element method (e.g., Zienkiewicz and Taylor, 1989;
Ottosen and Petersson, 1992).
When solving the strong form of the equation of motion, the boundary con-
ditions at material discontinuities must be explicitly satisfied. The traction con-
tinuity at internal discontinuities and vanishing traction at the free surface are
automatically satisfied by the weak-form solution (this is an advantage of the
weak form). In contrast, continuity of displacement is an essential continuity con-
dition that must be explicitly satisfied by the weak-form solutions.
Whereas most of the recent FD schemes solve one of the strong forms, Geller
and Takeuchi (1995, 1998) developed their optimally accurate FD schemes in
application to the weak form of Strang and Fix (1973) and Geller and Ohminato
(1994).
In principle, any formulation can be used with any of the three types of the
grids (conventional, partly-staggered, staggered). However, it is obvious that the
displacement formulation is naturally connected with the conventional grid. This
is because only displacement values are explicitly present both in the equations
430 MOCZO ET AL.

and grid. Similarly, the velocity-stress formulation is naturally connected with the
staggered grid as all field quantities in the equation of motion and Hooke’s law
are explicitly present in the grid. The particular structure (i.e., relative positions
of the field quantities in space and time) of the staggered grid is unambiguously
determined by the structure of the equation of motion and Hooke’s law, that is, by
the temporal and spatial derivatives of the field quantities.
In the early days of the FD method applied to seismology and seismics,
the displacement formulation and conventional grid were used; for example,
Alterman and Karal (1968), Boore (1970, 1972), Kelly et al. (1976). Because the
conventional-grid displacement FD schemes had problems with instabilities in
models with high-velocity contrasts and with grid dispersion in media with high
Poisson’s ratio, Virieux (1984, 1986) introduced the staggered-grid velocity-stress
FD schemes for modeling seismic wave propagation. Virieux followed Madariaga
(1976) who introduced the staggered-grid formulation in his dynamic modeling
of the earthquake rupture. Bayliss et al. (1986) and Levander (1988) introduced
the 4th-order staggered-grid FD schemes which in 2-D and 3-D need at least four
and eight times less memory, respectively, compared to the 2nd-order schemes.
In terms of CPU the improvement is 5–8 times in 2-D and 10–16 in 3-D. This
is related to the grid dispersion. The staggered-grid FD schemes have become
the dominant type of schemes in the FD modeling of seismic wave propagation
and earthquake motion. In order to further reduce the memory requirements, Luo
and Schuster (1990) suggested a staggered-grid displacement-stress 2-D P-SV FD
scheme which they called a parsimonious scheme. Because the scheme does not
integrate stress in time, the stress-tensor components are only temporary quanti-
ties. Thus, the displacement-stress scheme in 3-D needs only 75% of the memory
needed by the velocity-stress scheme. Rodrigues (1993), and Yomogida and Et-
gen (1993) used the 8th-order 3-D displacement-stress FD schemes, Ohminato
and Chouet (1997) applied the 2nd-order while Moczo et al. (2000, 2002) the
4th-order approximations. Moczo et al. (2000) analyzed the grid dispersion of the
displacement-stress schemes (4th and 2nd order) and pointed out that the stability
and grid dispersion of the displacement-stress, displacement-velocity stress and
velocity-stress schemes are the same. The advantage of the displacement-velocity-
stress scheme is that both displacement and particle-velocity are calculated at no
extra cost.
The partly-staggered grid was probably first used in seismology by Andrews
(1973) who applied it in modeling the fault rupture using his traction-at-split-
node method. Magnier et al. (1994) realized disadvantages of the staggered-grid
schemes in treating the anisotropic media and used the partly-staggered grid.
Zhang (1997) used the partly-staggered grid in his 2-D velocity-stress FD model-
ing. However, the developed schemes have not attracted much attention. Recently,
the use of the partly-staggered grid was promoted by Saenger et al. (2000) and
Saenger and Bohlen (2004). They called the grid rotated staggered grid since they
obtained the spatial FD operator by the rotation of the standard staggered-grid
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 431

operator. The term “rotated staggered grid” is somewhat unfortunate because as-
suming one spatial grid position for the stress tensor and another position for the
displacement vector, it is possible to find a variety of FD schemes—depending
on the order of approximation. Only in one particular case can the spatial FD op-
erator be obtained by the rotation of the standard staggered-grid operator though
it is also easy to obtain it without explicit consideration of the rotation. In fact,
the particular scheme used by Saenger et al. (2000) is a simple consequence of
requirement of the same truncation error with respect to all coordinate axes.
While the reason to use the partly-staggered grids for anisotropic media is obvi-
ous (all stress-tensor components located at the same grid position thus requiring
no interpolation of particle velocities or strains), the application of the partly-
staggered grid to account for the material heterogeneity has not been analyzed.
Therefore, this application still needs to be theoretically more rigorously ana-
lyzed.
A special aspect of the development is the FD schemes formulated for non-
uniform grids. These will be mentioned in the section on the memory optimization
and parallelization.
Apart from the recent dominance of the staggered-grid FD schemes, Geller
and his co-workers developed another important approach to the FD modeling
of seismic wave propagation based on the weak form of the equation of motion.
In their schemes displacement is the sole dependent variable, as opposed to the
staggered-grid schemes. Geller and Takeuchi (1995) derived a general criterion
for optimally accurate numerical operators, and used it to derive an optimally
accurate frequency-domain scheme. Geller and Takeuchi (1998) used the crite-
rion to derive an optimally accurate FDTD scheme for 1-D problems. Takeuchi
and Geller (2000) then developed optimally accurate FDTD operators for 2-D
and 3-D problems. These derivations yield implicit schemes which are solved
using a predictor-corrector algorithm, so that the actual computational schemes
are explicit. Whereas optimally accurate FDTD schemes require at least twice
the CPU time per grid point and time step compared to 2nd-order staggered-grid
FD schemes, they yield accuracy improvements on the order of 10 (for 1-D),
50 (for 2-D), or 100 (for 3-D). From this point of view they are cost-effective.
The optimally accurate schemes have not yet been widely used in practical FD
modeling. The likely reasons are (a) the theory might appear relatively compli-
cated compared to that of the standard staggered-grid schemes, (b) the fact that
quantification and minimization of computational error have not heretofore been
widely viewed as high priorities, (c) inertia with respect to traditional approaches
and the lack of user-friendly codes for optimally accurate schemes. However, it is
likely that optimally accurate FD schemes will be more widely used in the future.
An interesting approach to develop alternative FD scheme was presented by
Holberg (1987). Instead of minimizing the error measured in terms of higher-
order derivatives he minimized the relative error in group velocity caused by
the grid dispersion within a specific frequency band emitted by active sources.
432 MOCZO ET AL.

Therefore he did not develop the FD operator with the predetermined order of the
truncation error using the Taylor expansions. Holberg defined the FD operator as
a differentiator realized by a convolutory (FIR) operator.
Whereas most of the recent FD schemes are 2nd-order accurate in time, it is
possible to increase the order of accuracy by employing the Lax–Wendroff correc-
tion (Lax and Wendroff, 1964; Dablain, 1986). The leading term in the truncation
error for the 2nd-order FD approximation to time derivative is replaced by a term
with only spatial derivatives using the equation of motion. Although the form of
the Lax–Wendroff schemes is quite different from the optimally accurate schemes
of Geller and co-workers, Mizutani et al. (2000) have shown that these two types
of scheme are in fact essentially identical. An efficient implementation of the ap-
proach for viscoelastic media was presented by Blanch and Robertsson (1997).
An alternative approach for the time stepping is to use the Chebychev expansion
method by Tal-Ezer et al. (1990) which combines computational efficiency with
spectral accuracy in time.

1.7. Algorithms for Enhancing Computational Performance

A significant part of the literature has been devoted to enhancing the com-
putational efficiency of FDTD by taking a more holistic view as to how FDTD
are applied to solve a computational problem. Several methods for hybrid mod-
eling were developed in the eighties and nineties where different computational
techniques are used either for temporal/spatial dependences (e.g., Alekseev and
Mikhailenko, 1980) or in different parts of the model appropriate for the local
wave propagation regime (e.g., Shtivelman, 1984, 1985; Kummer et al., 1987;
Stead and Helmberger, 1988; Emmerich, 1989, 1992; Fäh, 1992; Robertsson et
al., 1996; Zahradník and Moczo, 1996; Moczo et al., 1997; Lecomte et al., 2004).
For instance, for a reflection seismic problem, in a smoothly varying overburden
a ray-based solution may be appropriate whereas FDTD are needed only in the
vicinity of a complex target where the energy reflects and is reverberated back to-
wards the surface. Such an approach can enable computational savings of several
orders of magnitude depending on the specific application and model and will be
particularly significant in 3-D. Problems with implementing and using such meth-
ods are related to model generation and interfacing the different methods which
tend to be different from problem to problem making the process very difficult to
implement in an automatic fashion.
A somewhat different approach was taken by, for instance, Robertsson and
Chapman (2000) and Robertsson et al. (2000) who proposed to use a technique,
referred to as FD-injection, fully based on finite differences for regenerating the
FD response in a model following local model alterations. The computational sav-
ings can be very significant, particularly for applications such as full waveform
inversion or time-lapse seismic analyses.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 433

Recently, van Manen et al. (2005) have proposed a method based on concepts
from time-reversal acoustics (Fink, 1997), that may become an important tool in
synthesizing seismic data. The approach relies on a representation theorem of the
wave equation to express the Green’s function between points in the interior as an
integral over the response in those points due to sources on a surface surrounding
the medium. Following a predictable initial computational effort, Green’s func-
tions between arbitrary points in the medium can be computed as needed using a
simple cross-correlation algorithm.

2. O PTIMALLY ACCURATE FD O PERATORS

A linear mechanical system or a finite volume of elastic continuum preferably

supports oscillatory motion at certain frequencies which are called eigenfrequen-
cies or resonant frequencies, i.e. at normal modes. In other words, the oscillatory
motion is naturally most amplified at these frequencies. The same is also true
about a numerical error in a discrete numerical system which is a numerical ap-
proximation to the true physical system. The basic idea of Geller and Takeuchi
(1995) therefore seems quite obvious: to minimize the error of the numerical so-
lution first of all at eigenfrequencies.
Geller and Takeuchi (1995) used first-order Born theory and a normal mode
expansion to obtain formal estimates of the relative error of the numerical solu-
tion and a general criterion for what they called optimally accurate operators. The
criterion requires that the inner product of an eigenfunction and the net error of
the discretized equation of motion should be approximately equal to zero when
the operand is the eigenfunction and the frequency is equal to the corresponding
eigenfrequency. An important aspect of the approach is that it is not necessary to
know the actual values of the eigenfrequencies and eigenfunctions to use the cri-
terion to derive optimally accurate operators. Geller and Takeuchi (1995) showed
that in the case of a heterogeneous medium the criterion is the logical extension
of the criterion to minimize grid dispersion of phase velocity for a homogeneous
medium. Based on this criterion, Geller and Takeuchi (1998) derived optimally
accurate 2nd-order weak-form FDTD scheme for the elastic 1-D case. Takeuchi
and Geller (2000) then generalized their approach to the 2-D and 3-D cases.
Though the optimally accurate FD operators are not applicable to the staggered-
grid schemes (they would lead to apparently intractable implicit schemes—
according to Geller and Takeuchi, 1998), we consider them important and as-
sume their wider applications in future FDTD modeling. Therefore, we briefly
present the basics of Geller and Takeuchi’s approach, closely following Geller
and Takeuchi (1995, 1998).
434 MOCZO ET AL.

2.1. General Criterion for Optimally Accurate FD Operators

In the Direct Solution Method (DSM; Geller and Ohminato, 1994) the equation
of motion for an anelastic solid is transformed into its Galerkin weak form (Strang
and Fix, 1973),
 2
ω T − H c = − g, (10)
where ω is the angular frequency, T mass matrix, H stiffness matrix, c vector of
expansion coefficients for the trial functions, g force vector,
(r) ∗ (s) (r) ∗ (s)
Trs = φi ρφi dV , Hrs = φi ,j cij kl φk ,l dV ,
V V
(r) ∗
gr = φi fi dV , (11)
V
(r)
φi is the ith component of the rth trial function, and ∗ means complex conjugate
quantity. The displacement is represented as
(r)
ui = cr φi . (12)
r
If an infinite trial function expansion were used, Eq. (10) would yield exact so-
lutions. In any practical application the trial function expansion will be finite and
there will be some numerical error. The exact equation of motion can be formally
written as
 2 e
ω T − He ce = − g. (13)
The relation between the numerical and exact quantities is assumed in a form
−
→
T = Te + δT, H = He + δH, c = ce + δc, (14)
−
→
where δT, δH and δc are errors of the numerical operators and solution, respec-
tively.
The normal modes satisfy equation
 2 e
ωp T − He cp = 0, (15)
where ωp is an eigenfrequency of the pth mode and cp is the eigenvector. Ortho-
normalization is assumed in a form

cp∗ He cq = ωp2 cp∗ Te cq = ωp2 δpq . (16)

Substituting Eqs. (14) into the l.h.s. of Eq. (10), replacing the r.h.s. of Eq. (10) by
the l.h.s. of Eq. (13), and neglecting terms with products of errors (the first-order
Born approximation) leads to
 2 e −→ 
ω T − He δc = − ω2 δT − δH ce . (17)
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 435

−
→
Equation (17) enables to determine the error of the numerical solution, δc, if the
exact solution and errors of the operators, i.e., ce , δT and δH are known.
The solution of Eq. (13) can be represented in terms of an eigenfunction expan-
sion

ce = dpe cp . (18)

p

Substituting Eq. (18) into Eq. (13), and using Eq. (16) leads to

dpe = cp∗ g/ ω2 − ωp2 . (19)
The expansion coefficient dpe will be large, when ω is close to ωp . Otherwise, it
will be negligible.
The solution of Eq. (17) can also be represented in terms of an eigenfunction
expansion
−
→
δc = δdp cp . (20)
p

Substituting expansions (18) and (20) into Eq. (17), and using Eqs. (16) leads to
 2 ∗ 
δdp = − cq − cp∗ δH
ω cp δT cq dqe ω2 − ωp2 . (21)
q

The expansion coefficient δdp will be large only when ω is close to ωp . In such
a case obviously only dpe will be large. Therefore, in the vicinity of ω = ωp , the
q = p terms in Eq. (21) can be neglected, i.e., the relative error of the numerical
solution in the vicinity of ωp will approximately be
δdp  2 ∗ ∗
 2
= − ω 
c p δT
cp − 
c p δH
cp ω − ωp2 . (22)
dpe
It follows from Eq. (22) that the relative error will in general greatly increase
with ω → ωp . However, if the numerator of Eq. (22) is also proportional to
ω − ωp , the relative error will remain approximately constant as ω → ωp . Such
proportionality can be achieved if and only if
.
ωp2 cp∗ δT
cp − cp∗ δH
cp = 0 (23)
for each mode. If Eq. (23) is approximately satisfied, then Eq. (22) can be simpli-
fied:
 
 δdp 
  ∗
 d e  ≈ cp δT
cp . (24)
p

This means that the relative error for a given grid can be reliably estimated in
advance of calculation.
436 MOCZO ET AL.

Geller and Takeuchi (1995) defined optimally accurate operators, say T and
H ,as operators that satisfy Eq. (23):
 .
cp∗ ωp2 δT − δH cp = 0. (25)
Substituting first two of Eqs. (14) for operators T and H into Eq. (25) and using
Eq. (15) leads to equivalent equation
 .
cp∗ ωp2 T − H cp = 0. (26)
Equation (25) will be satisfied if the leading term of the truncation error of the dis-
cretized equation is zero when the operand is an eigenfunction and the frequency
is equal to the corresponding eigenfrequency, in other words if
 2  .
ωp δT − δH cp = 0. (27)
On the other hand, however, it is not necessary for Eq. (27) to be satisfied in order
for Eq. (25) to be satisfied, because even if the quantity on the l.h.s. of Eq. (27)
is non-zero, its inner product with cp can still be approximately zero. Geller and
Takeuchi (1995, 1998) and Takeuchi and Geller (2000) take advantage of this
fact to derive optimally accurate operators using a two-step methodology. First,
for interior points, they derived optimally accurate operators that satisfy Eq. (27).
Second, to complete the derivation of the optimally accurate operators, they “fill
in” the few remaining degrees of freedom (corresponding to boundary points or
their neighbors) so that Eq. (25) is approximately satisfied, even though Eq. (27)
is not necessarily satisfied.
In this review, for simplicity, we discuss only the simplest case of optimally
accurate 1-D TDFD operators for interior points only. For a discussion of the
treatment of the boundary operators, see, for example, Sections 3 and 4 of Geller
and Takeuchi (1995).
Consider the equation

exact LHS(ω, u) = f (28)

and such its discretization which gives
discretized LHS(ω, u)
h2 
= exact LHS(ω, u) + exact LHS(ω, u) + · · · , (29)
a
where the primes denote spatial differentiation. The normal modes satisfy the
equation

exact LHS(ωp , up ) = 0, (30)

which implies

exact LHS(ωp , up ) = 0. (31)
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 437

Considering normal modes in Eq. (29), and substituting Eqs. (30) and (31) into
Eq. (29) leads to
.
discretized LHS(ωp , up ) = 0, (32)
which corresponds to condition (26). This shows that the leading term of the trun-
cation error of each FD approximation used to discretize the l.h.s. of Eq. (28) has
the same coefficient, h2 /a, and the displacement 2 times more differentiated than
in the approximated term. Thus, we have an indication for constructing optimally
accurate discretization.

2.2. Optimally Accurate FD Operators for 1-D problem in a Homogeneous

Medium
In order to illustrate the above theory, we closely follow the 1-D problem for a
homogeneous medium in Geller and Takeuchi (1998). Consider 1-D equation
∂ 2u
ρ ü = C + f. (33)
∂x 2
A weak-form FDTD approximation to the equation can be written as

I (M, i) − KI (M, i) U (M, i) = FI ,

Am m m
(34)
where m is the time level, at which Eq. (34) is approximated, I index of the spatial
position at which Eq. (34) is approximated, M time summation index, i spatial
summation index. Matrix Am I has the form
⎡ m ⎤
AI (m + 1, I − 1) Am I (m + 1, I ) AI (m + 1, I + 1)
m

I (M, i) =
⎣ Am (m, I − 1) ⎦,
Am I AmI (m, I ) AmI (m, I + 1)
AI (m − 1, I − 1) AI (m − 1, I ) AI (m − 1, I + 1)
m m m

(35)
m
matrix KI has an analogous structure. The conventional operators, corresponding
to the standard difference formula (4), are
   
ρ 0 1 0 C 0 0 0
AI = 2 0 −2 0 ,
m
KI = 2 1 −2 1 .
m
(36)
 t 0 1 0 h 0 0 0
The operator error at a single point in space and time approximately is (only the
leading term is given)
 
. 2 t ∂ 4 u h2 ∂ 4 u
δAmI (M, i) − δK m
I (M, i) U (M, i) = ρ − C , (37)
12 ∂t 4 12 ∂x 4 m,I
where δAm m
I and δKI are differences between the conventional and exact opera-
tors.
438 MOCZO ET AL.

In order to derive an optimally accurate FDTD scheme, Geller and Takeuchi

(1998) take the Fourier transform of Eq. (34), to obtain the FD equation in the
frequency domain for purposes of error analysis:

ÃI (i) − K̃I (i) Ũ (i) = F̃I . (38)

Because Eq. (38) has essentially the same form as Eq. (10), the above theory of the
error of the numerical solution can be applied. The difference is that the operator
errors here depend on frequency and wavenumber.
The Fourier transform of Eq. (37) is
  2 
. ω 2 2 t 2 h2 d2 d u
δ ÃI (i) − δ K̃I (i) Ũ (i) = ρω u − 2
C 2 . (39)
12 12 dx dx I
It is clear from Eq. (39) that the error is not equal to zero when the operand is an
eigenfunction and frequency is equal to eigenfrequency. This would be the case if
  
. ω 2 2 t h2 d2 d2 u
δ ÃI (i) − δ K̃I (i) Ũ (i) = − ρω 2
u + C (40)
12 12 dx 2 dx 2 I
because the expression inside the brackets is the l.h.s. of the exact homogeneous
equation of motion in the frequency domain, which is zero when the operand is
an eigenfunction and the frequency is the corresponding eigenfrequency. A time-
domain equivalent to Eq. (40) can be obtained by applying the inverse Fourier
transform to Eq. (40):

δAmI (M, i) − δKI (M, i) U (M, i)

m
 2 2 2   2 
.  t ∂ ∂ u ∂ 2u h2 ∂ 2 ∂ u ∂ 2u
= ρ − C + ρ − C . (41)
12 ∂t 2 ∂t 2 ∂x 2 12 ∂x 2 ∂t 2 ∂x 2 m,I
The operators that yield error (41) are
 
ρ 1/12 10/12 1/12
AI = 2 −2/12 −20/12 −2/12 ,
m
 t 1/12 10/12 1/12
 
C 1/12 −2/12 1/12
KmI = 2 10/12 −20/12 10/12 . (42)
h 1/12 −2/12 1/12
Though, obviously, the operators (42) come from the condition (41), which
corresponds to conditions (25) and (26), it is now easy to see relation to the con-
ventional operators. Consider, e.g., approximation to the second time derivative
at time level m and spatial position I . While the conventional approximation, for-
mula (4), is

∂ 2 u  . 1 
2  = 2 UIm−1 − 2UIm + UIm+1 , (43)
∂t m,I  t
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 439

the approximation corresponding to the optimally accurate operator Am

I in
Eq. (42) is
     
∂ 2 u  . 1 ∂ 2 u  ∂ 2 u  ∂ 2 u 
= + 10 2  + 2 (44)
∂t 2 m,I 12 ∂t 2 m,I −1 ∂t m,I ∂t m,I +1
where the conventional approximation (43) is used to approximate derivatives at
time level m at all three spatial positions I −1, I and I +1. Similarly, the optimally
accurate spatial second derivative is approximated over three time levels m − 1,
m and m + 1.
It is obvious that optimally accurate operators yield an implicit FD scheme. In
order to avoid necessity to solve a system of simultaneous linear equations at each
time step, Geller and Takeuchi (1998) applied a predictor-corrector scheme. Note
that the predictor-corrector optimally accurate scheme and the Lax–Wendroff
scheme are essentially the same (Mizutani et al., 2000).
As already mentioned, Takeuchi and Geller (2000) generalized the approach
based on the optimally accurate operators to solve 2-D and 3-D problems.

3. C OMPLEX V ISCOELASTIC M EDIA WITH M ATERIAL D ISCONTINUITIES

As already mentioned, accounting for realistic attenuation and heterogeneity of

the medium, particularly in the presence of material discontinuities in the Earth’s
interior, is of crucial importance in modeling seismic wave propagation and earth-
quake motion. While seismologists tried to model material discontinuities from
the very beginning of the application of the FDM in seismology, the incorpo-
ration of the realistic attenuation was made possible considerably later. We first
briefly review the problem of material heterogeneity and material discontinuities
in the strong formulation of the equation of motion.

3.1. A Material Discontinuity in the Elastic Medium

In one of the pioneering efforts on the application of the FDM to seismic wave
propagation, Alterman and Karal (1968) introduced the concept of fictitious grid
points in order to approximate boundary conditions on material discontinuities
in their displacement FD scheme. Difficulties in application of the homogeneous
approach to curved discontinuities led Boore (1972) to explicitly include a stress-
continuity condition on discontinuities differently from the homogeneous and
heterogeneous approaches. Due to poor numerical properties of his explicit con-
tinuous stress method, Boore (1972) applied the heterogeneous approach in his
SH modeling. In order to follow detailed variation of the torsion modulus and,
at the same time, to avoid derivatives of the modulus, he used the mathematical
440 MOCZO ET AL.

trick of Tikhonov and Samarskii (see, e.g., Mitchell, 1969, p. 23) and calculated
effective grid moduli as integral harmonic averages along a grid line between two
neighboring grid points. Ilan et al. (1975) and Ilan and Loewenthal (1976) applied
the homogeneous approach to the 2-D P-SV problem on horizontal and vertical
planar discontinuities. They used Taylor expansions of displacement to couple the
equation of motion with the boundary conditions. Kelly et al. (1976) presented
their heterogeneous P-SV schemes with simple intuitive averaging of material
parameters. They numerically compared the heterogeneous approach with the ho-
mogeneous one and showed unacceptable difference between the two approaches
in the case of the corner-edge model. Kummer and Behle (1982) followed the
approach of Ilan et al. (1975) and derived the 2nd-order SH schemes for grid
points lying on different types of segments of the step-like polygonal discontinu-
ity. Virieux (1984, 1986) used the velocity-stress formulation and staggered grid
introduced to seismology by Madariaga (1976). His P-SV FD scheme did not suf-
fer from stability problems caused by large values of Poisson’s ratio, which was
the case of all displacement schemes on conventional grids. Since the work by
Virieux, the staggered-grid FD schemes became the dominant schemes applied
to seismic wave propagation in heterogeneous media. Virieux also discussed the
discrepancy between the homogeneous and heterogeneous approaches found by
Kelly et al. (1976). He found it difficult to explain features of the solution obtained
by the homogeneous approach.
An attempt to incorporate internal boundary conditions into a displacement
FD scheme was made by Sochacki et al. (1991) who integrated the equation
of motion over a grid cell that could possibly contain a material discontinuity.
The integrated equation of motion was then discretized. Schoenberg and Muir
(1989) developed a calculus to replace a stack of thin flat elastic anisotropic
homogeneous layers by an equivalent (in the long-wavelength limit) homoge-
neous anisotropic medium. As a result, they found that stress-strain relation for
an averaged medium satisfies the boundary conditions at interfaces. Applying the
Schoenberg–Muir calculus simplifies modeling of wave propagation and, at the
same time, accounts for transversal anisotropy. Muir et al. (1992) applied the
Schoenberg–Muir calculus to a grid cell containing material discontinuity, i.e., in
general, they treated contents of the cell as a stack of thin flat layers that can be av-
eraged by the Schoenberg–Muir calculus. The two papers did not have the impact
on the heterogeneous FD schemes that they deserved—likely because the authors
did not make explicit reference to the question of the heterogeneous FD schemes.
Zahradník and Priolo (1995) explicitly pointed out the problem of heterogeneous
FD schemes. Assuming discontinuous material parameters in the equation of mo-
tion they obtained an expression whose dominant term is equivalent to the traction
continuity condition. As other developers of the heterogeneous FD schemes they
were not aware of the work by Schoenberg and Muir. Graves (1996) intuitively
suggested a formula for determination of effective material grid parameters in
the 3-D 4th-order velocity-stress staggered-grid schemes and numerically demon-
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 441

strated the good level of accuracy. Zhang and Symes (1998) developed a 2-D
4th-order full-stencil immersed interface technique to account for a curved ma-
terial discontinuity. In the first step, all grid points are solved using a standard
FD scheme. In the second step, each grid point whose stencil includes grid points
from both sides of the discontinuity is recalculated using previous time step’s
values with a special 25-point scheme determined using local boundary condi-
tions. Moczo et al. (2002) analyzed the 1-D problem in a medium consisting of
two halfspaces. They demonstrated their method on a simple physical model of
the contact of two media, found a heterogeneous formulation of the equation of
motion and Hooke’s law, and derived a heterogeneous FD scheme. Finally they
analyzed the 3-D problem, suggested a 4th-order heterogeneous staggered-grid
FD scheme, and demonstrated its accuracy compared to the standard staggered-
grid FD schemes. Their analysis will be followed in the next sections.

3.1.1. 1-D problem

Consider two elastic halfspaces with a welded interface in the plane x = 0. The
wave propagation in the halfspaces is described by equations
ρ ± ü± = σ,x
±
+ f ±, σ ± = C ± u±
,x (45)
+ −
where the superscript refers to one halfspace and superscript to the other.
Either u(x, t) is the x-component of the displacement u(ux , 0, 0), σ (x, t) the xx-
component of the stress tensor, f (x, t) the x-component of the body force per unit
volume f(fx , 0, 0) and C(x) = λ(x) + 2μ(x) in the case of P wave, or u(x, t)
is the y-component of the displacement u(0, uy , 0), σ (x, t) the xy-component
of the stress tensor, f (x, t) the y-component of the body force per unit volume
f(0, fy , 0) and C(x) = μ(x) in the case of the SH wave (the coordinate system
can always be rotated so that the S wave could be the SH wave). At the welded
interface the continuity of displacement and traction applies: u− (0) = u+ (0),
σ − (0) = σ + (0). Because a heterogeneous FD scheme should be nothing else than
a discrete approximation to a differential problem, a heterogeneous formulation of
the differential problem is to be found. This means the same form of the equation
of motion and Hooke’s law for a point at a material discontinuity as a point away
from the material discontinuity.
Let ϕ ± (x), c± (x) and g ± (x) be real functions of a real argument x such that
ϕ ± (x) = c± (x)g ± (x) (46)
and
ϕ − (0) = ϕ + (0). (47)
Functions c± (x) and g ± (x) may have discontinuities of the first order at x = 0.
Define
ḡ(0) = 0.5 · g − (0) + g + (0) . (48)
442 MOCZO ET AL.

Then it follows that

ϕ − (0) = ϕ + (0) = c̄(0)ḡ(0), c̄(0) = 2/ 1/c− (0) + 1/c+ (0) . (49)
If c± (x) = 1/r ± (x), then
1
ϕ − (0) = ϕ + (0) = ḡ(0), r̄(0) = 0.5 · r − (0) + r + (0) . (50)
r̄(0)
It follows from Eqs. (46) to (50) that the equation of motion and Hooke’s law for
a point at the material discontinuity have the form
ρ̄(0)ü(0) = σ,x (0) + f¯(0) (51)
and
σ (0) = C̄(0)u,x (0), (52)
respectively, with a density equal to the arithmetic average of the densities in the
two halfspaces, and elastic modulus equal to the harmonic average of the moduli
in the two halfspaces:
ρ̄(0) = 0.5 · ρ − (0) + ρ + (0) ,
C̄(0) = 2/ 1/C − (0) + 1/C + (0) . (53)
The average spatial derivatives of the stress and displacement are
σ,x (0) + f¯(0) = 0.5 · σ,x
− +
(0) + σ,x (0) + f − (0) + f + (0) ,
u,x (0) = 0.5 · u− +
,x (0) + u,x (0) . (54)
It is obvious that Eqs. (51) and (52) for a point at the material discontinuity
have the same form as the equation of motion and Hooke’s law (45), at a point
away from the material discontinuity. This provides a basis for the heterogeneous
1-D FD scheme.
Moczo et al. (2002) also showed that two Hooke elements (elastic springs) con-
nected in series make an appropriate rheological model for considering traction
continuity at the welded interface of two elastic materials.
Consider, for example, the velocity-stress formulation. Let VIm , TIm and FIm be
the discrete approximations to particle velocity vIm = v(I h, mt), stress σIm =
σ (I h, mt) and body force fIm = f (I h, mt). One possible heterogeneous 4th-
order staggered-grid FD scheme for the 1-D problem is
t  m−1/2 m−1/2
TIm+1/2 = TIm−1
+1/2 + CI +1/2 h a VI +2
H
− VI −1
 m−1/2 m−1/2
+ b VI +1 − VI ,
m+1/2 m−1/2 1 t  m 
VI = VI + A a TI +3/2 − TIm−3/2 + b TIm+1/2 − TIm−1/2
ρI h
t m
+ A FI (55)
ρI
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 443

with
xI +1/2  xI +1 −1
1 1 1
ρIA = ρ(x) dx, CIH+1/2 = dx . (56)
h xI −1/2 h xI C(x)
The coefficients in Eqs. (55) are a = −1/24 and b = 9/8. In the scheme, the
averaging of the spatial derivatives of the functions at the material discontinuity is
neglected while the harmonic averaging of the elastic moduli and arithmetic aver-
aging of densities at the material discontinuity is taken into account (see Moczo et
al., 2002 for details). In general, the integrals are evaluated numerically. It is easy
to check that the scheme yields very good accuracy in smoothly and/or discon-
tinuously heterogeneous media. The scheme is capable to sense a true position of
the material discontinuity no matter what the position of the discontinuity is with
respect to the grid points.
While Tikhonov and Samarski (e.g., Mitchell, 1969, p. 23) obtained the har-
monic averaging as a result of the mathematical ‘trick’ to avoid spatial derivatives
of the coefficients in the 2nd-order displacement formulation, the harmonic aver-
age in the heterogeneous formulation (52) is due to traction-continuity condition
at the material discontinuity.

3.1.2. 3-D problem

Define stress and strain vectors

σ = [σxx , σyy , σzz , σxy , σyz , σzx ]T ,
ε = [εxx , εyy , εzz , εxy , εyz , εzx ]T . (57)
Hooke’s law for an isotropic medium can be written as
σ = Eε (58)
where
⎡ ⎤
κ + 43 μ κ − 23 μ κ − 23 μ 0 0 0
⎢κ − 23 μ κ + 43 μ κ − 23 μ 0 ⎥
⎢ 0 0 ⎥
⎢ ⎥
⎢ − 23 μ κ − 23 μ κ + 43 μ 0 0 0 ⎥
E = ⎢κ ⎥ (59)
⎢ 0 ⎥
⎢ 0 0 0 2μ 0 ⎥
⎣ 0 0 0 0 2μ 0 ⎦
0 0 0 0 0 2μ
is the elasticity matrix. Let moduli κ and μ have a discontinuity of the first-order
across a surface S with normal vector n. The surface S defines the geometry of
the material discontinuity (interface). The welded-interface boundary conditions
η) and traction T (
are continuity of displacement u( η, n) across the surface:
u+ (
η) = u− (
η), T + (
η, n) = T − (
η, n). (60)
444 MOCZO ET AL.

For simplicity, consider first the planar surface S parallel to the xy-coordinate
plane with a normal vector n = (0, 0, 1). The conditions (60) imply
+ − + − + −
σzx = σzx , σzy = σzy , σzz = σzz ,
+ − + − + −
εxx = εxx , εyy = εyy , εxy = εxy . (61)
The components σxx , σyy , σxy , εzx , εzy and εzz may be discontinuous across the
material discontinuity. Define averaged stress and strain vectors at the material
discontinuity:
1 + 1
σ A = σ + σ − , εA = ε+ + ε− . (62)
2 2
Due to the boundary conditions,
T
σ A = σxx
A A
, σyy A
, σzz , σxy , σyz , σzx ,
T
εA = εxx , εyy , εzz
A A
, εxy , εyz A
, εzx . (63)
Then Hooke’s law for a point on the material discontinuity is

σ A = ẼεA (64)
with the elasticity matrix
⎡  + 2μA   0 0 0 ⎤
⎢   + 2μA  0 0 0 ⎥
⎢  ⎥
⎢ κ + 43 μ
H
0 ⎥
⎢   0 0 ⎥
Ẽ = ⎢ ⎥ (65)
⎢ 0 0 0 2μA 0 0 ⎥
⎢ ⎥
⎣ 0 0 0 0 2μH 0 ⎦
0 0 0 0 0 2μH
and
 A 2    
κ − 23 μ 4 H (κ − 23 μ)μ A
= · κ + μ + 2 · ,
κ + 43 μ 3 κ + 43 μ
   
κ − 23 μ A 4 H
= · κ + μ . (66)
κ + 43 μ 3
Superscripts A and H denote arithmetic and harmonic averages, respectively.
Relation (64) means that for a point on the material discontinuity it is possi-
ble to find the same form of Hooke’s law as for a point inside a homogeneous
or smoothly heterogeneous medium, given by Eq. (58). Considering the point on
the discontinuity as a point of the averaged medium characterized by matrix Ẽ
assures traction continuity at the point. There is, however, an important differ-
ence between laws (64) and (58). The matrix Ẽ for the averaged medium given
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 445

by Eq. (65) has 5 independent nonzero elements and the averaged medium is
transversely isotropic. Matrix E for any of the two isotropic media in contact
has only 2 independent nonzero elements. This means that the exact hetero-
geneous formulation for a planar material discontinuity parallel with a coordi-
nate plane increases the number of elastic coefficients necessary to describe the
medium.
Next let us consider a planar material discontinuity in a general position in
a Cartesian coordinate system. Let the normal vector be n = (nx , ny , nz ) with
all non-zero elements. Find a Cartesian coordinate system x  y  z in which n is
parallel to the z -axis. Then it is possible to find a matrix Ẽ with 5 independent
non-zero elements. Transforming the matrix Ẽ into a matrix Ẽ in the original co-
ordinate system xyz yields a symmetric elasticity matrix Ẽ which has, in general,
all elements non-zero (5 being independent). This means that all strain-tensor
components are necessary to calculate each stress-tensor component at a point
of the interface (not the case with the standard staggered grid), and 21 non-zero
elastic coefficients are necessary at the point.
If the geometry of a material discontinuity is defined by a non-planar smooth
surface S, the surface may be locally approximated by a planar surface tangential
to surface S at a given point.
It is obvious that finding a heterogeneous formulation of the differential prob-
lem as a basis for a heterogeneous FD scheme for a medium with a material
discontinuity in a general 3-D problem is far more complicated than that in the
1-D problem. A non-simplified treatment and a corresponding FD scheme would
lead to a substantial increase in memory requirement.
Therefore, Moczo et al. (2002) suggested a simplified approach: (a) They
wanted to keep the structure, number of operations and memory requirements
of the standard 4th-order staggered-grid scheme. (b) At the same time, they chose
to determine an effective grid elastic modulus (κ or μ) at each grid position of
the stress-tensor components as volume harmonic average of the modulus within
a volume of the grid cell centered at the grid position. The latter choice was based
on the fact that harmonic averaging is exact in the 1-D case, see Eq. (49), and is a
part of exact averaging in the 3-D case, see Eqs. (64)–(66).
At each position of the displacement or particle-velocity component an effec-
tive grid density is determined as a volume arithmetic average of density within a
volume of the grid cell centered at the grid position. The averaging applies to both
smoothly and discontinuously heterogeneous media. The averages are evaluated
by numerical integration.
Let TIxx,m
+1/2,J +1/2,K+1/2 be the discrete approximation to the stress-tensor
component σxx [(I + 1/2)h, (J + 1/2)h, (K + 1/2)h, mt]. Similarly, let
x,m y,m z,m x,m
VI,J +1/2,K+1/2 , VI +1/2,J,K+1/2 , VI +1/2,J +1/2,K and FI,J +1/2,K+1/2 be discrete
approximations to the particle-velocity and body-force components. Examples of
446 MOCZO ET AL.

the FD schemes for the stress-tensor and particle-velocity components σxx and vx
are
TIxx,m
+1/2,J +1/2,K+1/2

= TIxx,m−1
+1/2,J +1/2,K+1/2
 
t 4
+ κIH+1/2,J +1/2,K+1/2 + μH
h 3 I +1/2,J +1/2,K+1/2
 x,m−1/2 x,m−1/2
× a VI +2,J +1/2,K+1/2 − VI −1,J +1/2,K+1/2
 x,m−1/2 x,m−1/2
+ b VI +1,J +1/2,K+1/2 − VI,J +1/2,K+1/2
 
2 H
+ κI +1/2,J +1/2,K+1/2 − μI +1/2,J +1/2,K+1/2
H
3
 y,m−1/2 y,m−1/2
× a VI +1/2,J +2,K+1/2 − VI +1/2,J −1,K+1/2
 y,m−1/2 y,m−1/2
+ b VI +1/2,J +1,K+1/2 − VI +1/2,J,K+1/2
 z,m−1/2 z,m−1/2
+ a VI +1/2,J +1/2,K+2 − VI +1/2,J +1/2,K−1

 z,m−1/2 z,m−1/2
+ b VI +1/2,J +1/2,K+1 − VI +1/2,J +1/2,K (67)

with
 xI +1 yJ +1 zK+1 −1
1 1
κIH+1/2,J +1/2,K+1/2 = dx dy dz , (68)
h3 xI yJ zK κ
 xI +1 yJ +1 zK+1 −1
1 1
I +1/2,J +1/2,K+1/2 =
μH dx dy dz (69)
h3 xI yJ zK μ
and
x,m+1/2
VI,J +1/2,K+1/2
x,m−1/2  A
= VI,J +1/2,K+1/2 + t/ hρI,J +1/2,K+1/2
 xx,m
× a TI +3/2,J +1/2,K+1/2 − TIxx,m −3/2,J +1/2,K+1/2

+ b TIxx,m
+1/2,J +1/2,K+1/2 − T xx,m
I −1/2,J +1/2,K+1/2
 xy,m xy,m  xy,m xy,m
+ a TI,J +2,K+1/2 − TI,J −1,K+1/2 + b TI,J +1,K+1/2 − TI,J,K+1/2
 zx,m zx,m  zx,m zx,m
+ a TI,J +1/2,K+2 − TI,J +1/2,K−1 + b TI,J +1/2,K+1 − TI,J +1/2,K
 x,m
+ t/ρI,J A
+1/2,K+1/2 FI,J +1/2,K+1/2 (70)
with
1 xI +1/2 yJ +1 zK+1
+1/2,K+1/2 =
A
ρI,J ρ dx dy dz. (71)
h3 xI −1/2 yJ zK
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 447

F IG . 2. The upper part: Positions of the upper and lower layer-halfspace interfaces in five models
of a layer (P-wave and S-wave velocities 1125 m/s and 625 m/s) between two identical halfspaces
(P-wave and S-wave velocities 5468 m/s and 3126 m/s)—shown schematically in one vertical grid
plane. The five models differ from each other by position of the upper layer-halfspace interface in the
spatial grid (the same grid for all models) and thus by the layer thickness. The lower part: Comparison
of our FD and DWN synthetics for the five models. Note very good accuracy of the FD synthetics for
any position of the layer-halfspace interface with respect to the spatial grid. Also note considerable
differences between synthetics due to variations in the layer thickness that is smaller than one grid
spacing. Reproduced from Moczo et al. (2002).
448 MOCZO ET AL.

Moczo et al. (2002) used a set of models to numerically test the scheme. The
tests demonstrated its very good numerical accuracy. Here we illustrate an im-
portant property of the scheme—the capability to sense the position of a material
discontinuity regardless of its position with respect to the spatial grid. Five dif-
ferent models of a single horizontal homogeneous layer located in between two
homogeneous halfspaces differ from each other by the thickness of the layer.
The spatial grid is one and the same in all five models (see the upper part of
Fig. 2). A double-couple point source was located in the lower halfspace. The
FD synthetics (ux component) are compared with those calculated by the discrete
wavenumber (DWN) method (Bouchon, 1981; computer code Axitra by Coutant,
1989), see the lower part of Fig. 2. The FD and DWN synthetics agree very well
regardless of the position of the upper layer-halfspace interface with respect to the
spatial grid. It is also clear from Fig. 2 that differences in thickness of the layer—
smaller than one grid spacing—cause considerable changes in seismic motion.
This often is underestimated by many modelers who consider the size of one grid
cell as “atom of resolution” within which a FD scheme cannot see differences.
This and other examples given by Moczo et al. (2002) clearly show that this is
not the case if the scheme is sufficiently accurate.

3.2. Incorporation of the Realistic Attenuation

3.2.1. Stress-Strain Relation in Viscoelastic Medium—the 1-D Case

The behavior of real Earth’s material can be described as a combination of

elastic solids and viscous fluids. The stress-strain relation therefore also should
depend on time. The rheology of a viscoelastic medium seems appropriate for
quantitative description of seismic wave propagation. Observations show that the
internal friction in the Earth is nearly constant over the seismic frequency range,
e.g. McDonal et al. (1958), Liu et al. (1976), Spencer (1981), Murphy (1982).
The stress-strain relation in a linear isotropic viscoelastic material is given by
the Boltzmann superposition principle. For the 1-D problem this is
t
σ (t) = ψ(t − τ )ε̇(τ ) dτ, (72)
−∞

where σ (t) is the stress, ε̇(t) the time derivative of the strain, and ψ(t) the stress
relaxation function defined as a stress response to a Heaviside unit step function
in strain. According to Eq. (72), the stress at a given time t is determined by the
entire history of the strain until time t. The integral in Eq. (72) represents a time
convolution of the relaxation function and the strain rate. Using the symbol ∗ for
the convolution, Eq. (72) can be written as

σ (t) = ψ(t) ∗ ε̇(t). (73)

 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 449

It follows from the definition of the relaxation function that its time derivative
M(t) = ψ̇(t) (74)
is the stress response to the Dirac δ-function in strain and that
σ (t) = M(t) ∗ ε(t). (75)
Let F and F −1
denote the direct and inverse Fourier transforms
  ∞
F x(t) = x(t) exp(−iωt) dt,
−∞

  1 ∞
F −1 X(ω) = X(ω) exp(iωt) dω,
2π −∞
where ω is the angular frequency. Relation (75) can be Fourier-transformed into
the frequency domain:
σ (ω) = M(ω) · ε(ω). (76)
In general, M(ω) is a complex, frequency-dependent viscoelastic modulus. Due
to properties of the Fourier transform
 
M(ω)
ψ(t) = F −1 . (77)
iω
Equation (76) indicates the correspondence principle in the linear theory of
viscoelasticity—in the frequency-domain, relations for the viscoelastic medium
are obtained by replacing real frequency-independent moduli by complex,
frequency-dependent quantities. Thus, the incorporation of the attenuation into
the frequency-domain computations is much easier than that in the time domain.
The time derivative of the stress is
σ̇ (t) = M(t) ∗ ε̇(t). (78)
An instantaneous elastic response of the viscoelastic material is given by the so-
called unrelaxed modulus MU , a long-term equilibrium response is given by the
relaxed modulus MR :
MU = lim ψ(t) = lim M(ω), MR = lim ψ(t) = lim M(ω).
t→0 ω→∞ t→∞ ω→0
(79)
The modulus defect or relaxation of modulus is
δM = MU − MR . (80)
Given the viscoelastic modulus, the quality factor Q(ω) is defined as
Q(ω) = Re M(ω)/ Im M(ω). (81)
Due to large computer time and memory requirements, the stress-strain relation
(72) in some cases only allowed simplified Q(ω) laws, e.g. linear Q(ω).
450 MOCZO ET AL.

3.2.2. Conversion of the Convolutory Stress-Strain Relation into

a Differential Form

If M(ω) is a rational function, the inverse Fourier transform of Eq. (76) yields
the nth-order differential equation for σ (t), which can be numerically solved
much more easily than the convolution integral. Day and Minster (1984) assumed
that, in general, the viscoelastic modulus is not a rational function. Therefore they
suggested approximating a viscoelastic modulus by an nth-order rational function
and determining its coefficients by the Padé approximant method. They obtained
n ordinary differential equations for n additional internal variables, which replace
the convolution integral. The sum of the internal variables multiplied by the un-
relaxed modulus gives an additional viscoelastic term to the elastic stress. The
work of Day and Minster not only developed one particular approach but, in fact,
indirectly suggested the future evolution—a direct use of the rheological models
whose M(ω) is a rational function of iω. In response to work by Day and Minster
(1984); Emmerich and Korn (1987) realized that an acceptable relaxation func-
tion corresponds to a rheology of what they defined as the generalized Maxwell
body—n Maxwell bodies and one Hooke element (elastic spring) connected in
parallel; Fig. 3. Because in the rheological literature the generalized Maxwell
body is defined without the additional single Hooke element, the abbreviation
GMB-EK will hereafter be used for the model defined by Emmerich and Korn.
Because the viscoelastic modulus of the GMB-EK has a form of a rational func-
tion, Emmerich and Korn (1987) obtained similar differential equations as Day

F IG . 3. Rheological model of the Generalized Maxwell Body (GMB-EK) defined by Emmerich

and Korn (1987). MH and Ml denote elastic moduli, ηl viscosity.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 451

and Minster (1984). In order to fit an arbitrary Q(ω) law they chose the relax-
ation frequencies logarithmically equidistant over a desired frequency range and
used the least-square method to determine weight factors of the relaxation mech-
anisms (classical Maxwell bodies). Emmerich and Korn (1987) demonstrated that
their approach is better than the approach based on the Padé approximant method
both in accuracy and computational efficiency. Independently, Carcione et al.
(1988a, 1988b), in accordance with the approach by Liu et al. (1976), assumed
the generalized Zener body (GZB)—n Zener bodies (ZB, standard linear bodies),
connected in parallel; Fig. 4. Carcione et al. (1988a, 1988b) developed a theory
for the GZB and introduced the term memory variables for the obtained additional
variables.
After publications by Emmerich and Korn (1987) and Carcione et al. (1988a,
1988b) different authors chose to use either the GMB-EK (for example, Emmerich,
1992; Fäh, 1992; Moczo and Bard, 1993; Moczo et al., 1997; Kay and Krebes,
1999) or GZB (for example, Robertsson et al., 1994; Blanch et al., 1995;

F IG . 4. Rheological model of the Generalized Zener Body (GZB). For a classical Zener body
(standard linear body) there are two equivalent models: H-p-M, that is, Hooke element connected in
parallel with Maxwell body, and H-s-KV, that is, Hooke element connected in series with Kelvin–Voigt
body. In the H-p-M model it is easier to recognize the relaxed modulus MRl and modulus defect δMl .
M1l and M2l in the H-s-KV model denote elastic moduli. In both models ηl stands for viscosity.
452 MOCZO ET AL.

Xu and McMechan, 1995; Robertsson, 1996; Hestholm, 1999). In both cases

the authors followed the corresponding mathematical formalisms. Moczo et al.
(1997) applied the GMB-EK approach also in the finite-element method and hy-
brid FD-finite-element method. Emmerich and Korn (1987), Emmerich (1992),
Fäh (1992), and Moczo and Bard (1993) defined one memory variable for one
displacement component. Robertsson et al. (1994) introduced the memory vari-
ables based on the GZB rheology into the staggered-grid velocity-stress FD
scheme. Blanch et al. (1995) suggested an approximate single-parameter method,
τ -method, to approximate the constant Q(ω) law. Xu and McMechan (1998) used
simulated annealing for determining a best combination of relaxation mechanisms
to approximate a desired Q(ω) law.
There appears to have been no or little comments by the authors using the
GZB on the rheology of the GMB-EK and the corresponding algorithms, and vice
versa. Thus, two parallel sets of publications and algorithms had been developed
during years. Therefore, Moczo and Kristek (2005) addressed this development
and showed relation between the two rheologies.

3.2.3. Rheologies of the GMB-EK and GZB

There are simple rules in the time and frequency domains for the mathematical
representation of the linear rheological models consisting of Hooke and Stokes
elements (springs and dashpots) connected in parallel or series. In the frequency
domain, the stress-strain relations for the Hooke and Stokes elements are σ (ω) =
M · ε(ω) and σ (ω) = iωη · ε(ω), respectively, where M is the elastic modulus,
and η viscosity. If two elements are connected in series, stresses are equal while
strains additive. If two elements are connected in parallel, stresses are additive
while strains are equal.
The application of the frequency-domain rules to the GMB-EK yields
n
iMl ω Ml
M(ω) = MH + , ωl = , l = 1, . . . , n, (82)
ωl + iω ηl
l=1

where ωl is a relaxation frequency. The relaxed and unrelaxed moduli are

n
MR ≡ lim M(ω) = MH , MU ≡ lim M(ω) = MR + Ml . (83)
ω→0 ω→∞
l=1

Since MU = MR + δM, we get Ml = δMl , and it is possible to assume

n
δMl = al δM, al = 1 (84)
l=1
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 453

without any simplification. Then

n
ial ω
M(ω) = MR + δM . (85)
ωl + iω
l=1

Using relation (77) it is straightforward to obtain the relaxation function

 n

ψ(t) = MR + δM al e−ωl t · H (t), (86)
l=1

where H (t) is the Heaviside unit step function. The above formulae were pre-
sented by Emmerich and Korn (1987).
From the two equivalent models of the ZB, shown in Fig. 4, we choose the
H-p-M type to obtain M(ω). This is because it is easy to recognize the relaxed
modulus and modulus defect in the ZB. The application of the frequency-domain
rules to the GZB results in
n
1 + iτεl ω
M(ω) = MRl (87)
1 + iτσ l ω
l=1

with relaxation times

η l MU l ηl τεl MU l
τεl = , τσ l = , = (88)
δMl MRl δMl τσ l MRl
and

MU l = MRl + δMl . (89)

The unrelaxed and relaxed moduli are
n n n
τεl
MR = MRl , MU = MRl = MR + δMl . (90)
τσ l
l=1 l=1 l=1

Relations (77) and (87) yield the relaxation function

 n     
τεl t
ψ(t) = MRl 1 − 1 − exp − · H (t). (91)
τσ l τσ l
l=1

Assumption of simplification (Carcione, 2001)

1
MRl = MR (92)
n
leads to
n
MR 1 + iτεl ω
M(ω) = ,
n 1 + iτσ l ω
l=1
454 MOCZO ET AL.
 n    
1 τεl t
ψ(t) = MR 1− 1− exp − · H (t). (93)
n τσ l τσ l
l=1

Formulae (92) and (93) were presented by Carcione (2001). Unfortunately, all pa-
pers dealing with the incorporation of the attenuation based on the GZB, starting
from Liu et al. (1976) until now, despite the book by Carcione (2001), have the
same error—the missing factor 1/n in the viscoelastic modulus and relaxation
function.

3.2.4. The Relation Between the GZB and GMB-EK

Following Moczo and Kristek (2005), consider again the ZB (H-p-M) model.
The application of the frequency-domain rules to the lth ZB yields
   
1 1 MRl MRl
σl (ω) · + = 1+ + · ε(ω). (94)
δMl iηl ω δMl iηl ω
Defining
δMl
ωl = (95)
ηl
and rearranging Eq. (94) gives
iδMl ω
σl (ω) = Ml (ω) · ε(ω), Ml (ω) = MRl + . (96)
ωl + iω
For n ZB connected in parallel, that is, for the GZB (Fig. 4), the stress is
n
 n 
σ (ω) = σl (ω) = Ml (ω) · ε(ω) (97)
l=1 l=1

and thus
n n
iδMl ω
M(ω) = MRl + . (98)
ωl + iω
l=1 l=1

Since
n n
MR = MRl , M U = MR + δMl , MU = MR + δM, (99)
l=1 l=1

it is possible to define
n
δMl = al δM, al = 1 (100)
l=1
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 455

without loss of generality and obtain

n
ial ω
M(ω) = MR + δM . (101)
ωl + iω
l=1

The viscoelastic modulus (101) obtained for the GZB (H-p-M), Fig. 4, is exactly
the same as it has been obtained by Emmerich and Korn (1987) for their GMB-
EK, Fig. 3. Obviously, M(ω) for the GZB (H-s-KV) would be the same. It is also
easy to rewrite the non-simplified ψ(t) for the GZB, Eq. (91), into the form of
ψ(t) for the GMB-EK, Eq. (86), without any simplification. In other words, the
rheology of the GMB-EK and GZB is one and the same. As a consequence, the
GMB-EK will be used in the following.

3.2.5. Introduction of the Anelastic Functions (Memory Variables)—the 1-D Case

In order to focus on the essential aspects of the implementation of the realistic

attenuation in the time-domain computations, we continue by considering the 1-D
case with one stress and one strain component. Using the unrelaxed modulus, the
viscoelastic modulus (101) and relaxation function (86) are rewritten as
n
al ω l
M(ω) = MU − δM ,
ωl + iω
l=1
 n

 −ωl t
ψ(t) = MU − δM al 1 − e · H (t). (102)
l=1

The time derivative of the relaxation function (time-dependent modulus) is

M(t) = ψ̇(t)
n
= −δM al ωl e−ωl t · H (t)
l=1
 n


+ MU − δM al 1 − e−ωl t · δ(t). (103)
l=1

Inserting Eq. (103) into Eq. (75) yields

n t
σ (t) = MU · ε(t) − δM al ω l ε(τ ) · e−ωl (t−τ ) dτ. (104)
l=1 −∞

The convolution integral can be replaced by additional functions (internal vari-

ables, new variables, memory variables, anelastic functions). While Day and
Minster (1984), Emmerich and Korn (1987), and Carcione et al. (1988a, 1988b)
defined the additional functions as material-dependent, Kristek and Moczo (2003)
456 MOCZO ET AL.

defined their anelastic functions as material-independent (the reason will be ex-

plained later):
t
ζl (t) = ωl ε(τ ) · e−ωl (t−τ ) dτ, l = 1, . . . , n. (105)
−∞
The stress-strain relation then becomes
n
σ (t) = MU · ε(t) − δMal ζl (t). (106)
l=1
Equations necessary to solve for the anelastic functions are easily obtained by
taking the time derivative of Eq. (105):
d t
ζ̇l (t) = ωl ε(τ ) · e−ωl (t−τ ) dτ = ωl −ζl (t) + ε(t) (107)
dt −∞
and
ζ̇l (t) + ωl ζl (t) = ωl ε(t), l = 1, . . . , n. (108)
Equations (106) and (108) define the time-domain stress-strain relation for the
viscoelastic medium whose rheology corresponds to the rheology of the GMB-
EK (and to its equivalent—the GZB). If the staggered-grid velocity-stress FD
scheme is to be used, then the time derivative of the stress is needed. It is easy to
obtain
n
σ̇ (t) = MU · ε̇(t) − δMal ξl (t) (109)
l=1
and
ξ̇l (t) + ωl ξl (t) = ωl ε̇(t), l = 1, . . . , n. (110)
As mentioned before, formalism developed specifically for the GZB was used
in many papers. Therefore we give here equations equivalent to those presented
by Robertsson et al. (1994). Using Eqs. (91), (74), (78) and (90) it is easy to obtain
n
σ̇ (t) = MU · ε̇(t) − rl (t), (111)
l=1
  t
MRl τεl
rl (t) = 1− ε̇(τ ) · exp −(t − τ )/τσ l dτ, l = 1, . . . , n,
τσ l τσ l −∞
(112)
 
1 MRl τεl
ṙl (t) + rl (t) = 1− ε̇(t), l = 1, . . . , n. (113)
τσ l τσ l τσ l
Note that anelastic functions (memory variables) rl (t) are material-dependent.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 457

Defining anelastic coefficients (different from those used by Emmerich and

Korn, 1987)

Yl = al δM/MU , l = 1, . . . , n, (114)
the stress-strain relations (106) and (109) become
n
σ (t) = MU · ε(t) − MU Yl ζl (t),
l=1
n
σ̇ (t) = MU · ε̇(t) − MU Yl ξl (t). (115)
l=1

The related Eqs. (108) and (110) remain unchanged.

It is clear that the stress or its time derivative can be calculated if the anelastic
coefficients and unrelaxed modulus are known.
The anelastic coefficients Yl , l = 1, . . . , n have to be determined from Q(ω)-
law. Using the anelastic coefficients, the viscoelastic modulus and quality factor
(81) are
 n

ωl
M(ω) = MU 1 − Yl (116)
ωl + iω
l=1

and
 n
 n

1 ωl ω ωl2
= Yl 2 1− Yl . (117)
l +ω ωl2 + ω2
Q(ω) ω 2
l=1 l=1

Equation (117) yields

ωl ω + ωl2 Q−1 (ω)

n
Q−1 (ω) = Yl . (118)
l=1
ωl2 + ω2

Equation (118) can be used to numerically fit any Q(ω)-law. A sufficiently

accurate approximation to nearly constant Q(ω) is obtained if the relaxation
frequencies ωl are distributed logarithmically equidistant over the frequency
range of interest. If, for example, Q(ω) values are known at frequencies ω̃k ;
k = 1, . . . , 2n − 1, with ω̃1 = ω1 , ω̃2n−1 = ωn , Eq. (118) can be solved for
Yl , l = 1, . . . , n using the least square method. A more detailed discussion of the
frequency range and its sampling at frequencies ω̃k can be found in the papers by
Blanch et al. (1995) and Graves and Day (2003; Eqs. (13) and (14)).
If an elastic P-wave velocity α or S-wave velocity β is known, then, in the
considered 1-D problem, MU = ρα 2 for the P-wave or MU = 2ρβ 2 for the
S-wave. In practice, a phase velocity at a certain reference frequency ωr can be
458 MOCZO ET AL.

measured or estimated. The P- or S-wave phase velocity c(ω) is given by

  
1 M(ω) −1/2
= Re . (119)
c(ω) ρ
It follows (Moczo et al., 1997) from Eqs. (116) and (119) that
R + 1  1/2
MU = ρc2 (ωr ) , R = 21 + 22 , (120)
2R 2

n n
1 ωr /ωl
1 = 1 − Yl , 2 = Yl . (121)
1 + (ωr /ωl )2 1 + (ωr /ωl )2
l=1 l=1

As already pointed out, a constant or almost constant Q is of great importance.

Therefore, Blanch et al. (1995) addressed the question of an efficient and suffi-
ciently accurate curve-fitting procedure in the case of constant Q. Their τ -method
is based on the fact that the level of attenuation caused by a ZB can be determined
by a dimensionless variable τ = (τε − τσ )/τσ . Blanch et al. (1995) derived ex-
plicit closed formula to determine parameters of the GZB for a desired constant Q,
for P- and S-waves, respectively. The GZB obtained by tuning through a single
parameter τ yields a very good constant-Q approximation.

3.2.6. A FD Scheme for the Anelastic Functions in the 1-D Case

The 2nd-order approximations to the anelastic functions ζl and ζ̇l , l = 1, . . . , n

give
. 1
ζl (tm ) = ζl (tm+1/2 ) + ζl (tm−1/2 ) ,
2
. 1 
ζ̇l (tm ) = ζl (tm+1/2 ) − ζl (tm−1/2 ) , (122)
t
where tm denotes the mth time level. Then each of the equations for the anelastic
functions can be solved by
2ωl t 2 − ωl t
ζl (tm+1/2 ) = ε(tm ) + ζl (tm−1/2 ). (123)
2 + ωl t 2 + ωl t
The value of ζl (tm ) needed in the stress-strain relation
n
σ (tm ) = MU · ε(tm ) − MU YlM ζl (tm ), (124)
l=1

is obtained from ζl (tm−1/2 ) and ζl (tm+1/2 ) using Eq. (122). This means that two
values have to be kept in memory for one spatial position at one time. It is, how-
ever, possible (Kristek and Moczo, 2003) to avoid the necessity to keep in memory
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 459

both values. It follows from Eqs. (123) and (122) that

ωl t 2
ζl (tm ) = − ε(tm ) + ζl (tm+1/2 ). (125)
2 − ωl t 2 − ωl t
Then the stress-strain relation (124) can be obtained in the form
n
σ (tm ) = M̃ε(tm ) − ỸlM ζl (tm+1/2 ) (126)
l=1

where
 n

M̃ = MU 1 + G1l YlM , ỸlM = G2l MU YlM ,
l=1
ωl t 2
G1l = , G2l = . (127)
2 − ωl t 2 − ωl t
Using scheme (123) and a proper scheme for Eq. (126) it is sufficient to have only
one variable for one anelastic function at one grid position at one time. In the
case of the staggered-grid velocity-stress FD scheme, the form of equations is the
same; only ζl and ε have to be replaced by ξl and ε̇, respectively.

3.2.7. A Material Discontinuity in the Viscoelastic Medium—the 1-D Case

It is not a trivial task to find a heterogeneous formulation to the differential

problem if the stress is given in the form of Eq. (115). Kristek and Moczo (2003)
suggested an approximate approach which has been shown sufficiently accurate
using numerical tests against the discrete wavenumber method (Bouchon, 1981;
Coutant, 1989).
Consider a contact of two viscoelastic media with the GMB-EK rheology. Each
of the two media is described by a real density and complex frequency-dependent
modulus given by Eq. (116). The question is how to determine density, elastic
(unrelaxed) modulus M̄U , and anelastic coefficients YlM̄ , l = 1, . . . , n for an aver-
aged medium that should represent the contact of two media (that is the boundary
conditions at the interface between the two media) if a material discontinuity goes
through a grid cell. There is no reason to consider other than volume arithmetic
averaging for the density using formula (56). An averaged viscoelastic modu-
lus M̄ can be obtained by numerical averaging in the frequency domain over the
grid cell. From the averaged viscoelastic modulus, the quality factor correspond-
ing to this modulus can be determined, Eq. (81), for example, at frequencies ω̃k ,
k = 1, . . . , 2n − 1; QM̄ (ω̃k ) = Re M̄(ω̃k )/ Im M̄(ω̃k ). Assuming that the rheol-
ogy of the averaged medium can be approximated by the GMB-EK rheology, the
anelastic coefficients YlM̄ , l = 1, . . . , n for the averaged medium can be obtained
using Eq. (118).
460 MOCZO ET AL.

It follows from Eq. (79) that M̄U = limω→∞ M̄(ω). An implication is that,
in the limit, the averaging of the viscoelastic modulus gives the averaging of the
unrelaxed modulus. This means that the unrelaxed (elastic) modulus M̄U for the
averaged viscoelastic medium can be obtained in the same way as in the perfectly
elastic medium.

3.2.8. A Summary of Equations in the 3-D Case

In the 3-D case it is assumed that the rheology of the medium is described by
one GMB-EK (or, equivalently, GZB) for the complex frequency-dependent bulk
modulus and one GMB-EK for the complex frequency-dependent shear modulus.
The stress-strain relation is (Kristek and Moczo, 2003)
 
1
σij = κεkk δij + 2μ εij − εkk δij
3
n   
μ ij 1 kk
− κYl ζl δij + 2μYl ζl − ζl δij
κ kk
(128)
3
l
where i, j, k ∈ {1, 2, 3}, the equal-index summation convention does not apply to
μ
l, κ(xi ) and μ(xi ) are unrelaxed (elastic) bulk and shear moduli, and Ylκ and Yl
are the corresponding anelastic coefficients. Assuming a measured or estimated
Qα (ω) for the P- and Qβ (ω) for the S-waves, the corresponding anelastic coeffi-
β
cients Ylα and Yl are obtained using Eq. (118). Then the anelastic coefficients Ylκ
μ
and Yl are
   
4 2 β  2 4 2
Yl = α Yl − β Yl
κ 2 α
α − β ,
3 3
μ β
Yl = Yl , l = 1, . . . , n. (129)
ij
There are n material-independent anelastic functions ζl for each of 6 strain-
tensor components satisfying equations
ij ij
ζ̇l + ωl ζl = ωl εij , l = 1, . . . , n, (130)
where the equal-index summation convention does not apply to l.
While Eqs. (128) and (130) are applicable to the conventional displacement, or
staggered-grid displacement-stress and displacement-velocity-stress FD schemes,
equations
 
1
σ̇ij = κ ε̇kk δij + 2μ ε̇ij − ε̇kk δij
3
n   
μ ij 1
− κYlκ ξlkk δij + 2μYl ξl − ξlkk δij (131)
3
l
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 461

and
ij ij
ξ̇l + ωl ξl = ωl ε̇ij , l = 1, . . . , n (132)
are needed for the staggered-grid velocity-stress FD schemes.
In analogy to the 1-D case, it is possible to obtain
ij 2ωl t 2 − ωl t ij
ζl (tm+1/2 ) = εij (tm ) + ζ (tm−1/2 ),
2 + ωl t 2 + ωl t l
l = 1, . . . , n (133)
and
 
1
σij (tm ) = κ̃εkk (tm )δij + 2μ̃ εij (tm ) − εkk (tm )δij
3
n 
− Ỹlκ ζlkk (tm+1/2 )δij
l=1
 
μ ij 1
+ 2Ỹl ζl (tm+1/2 ) − ζlkk (tm+1/2 )δij , (134)
3
where
 n
  n

μ
κ̃ = κ 1 + G1l Ylκ , μ̃ = μ 1 + G1l Yl ,
l=1 l=1
μ μ
Ỹlκ = G2l κYlκ , Ỹl = G2l μYl ,
ωl t 2
G1l = , G2l = . (135)
2 − ωl t 2 − ωl t
Equations (133) to (135) are ready for programming. In the case of the staggered-
ij
grid velocity-stress FD scheme, the form of equations is the same; only ζl and
ij
εij are replaced by ξl and ε̇ij .

3.2.9. Coarse Spatial Sampling

The incorporation of realistic attenuation considerably increases the number of

operations and variables/parameters that have to be kept in computer (core) mem-
ory. In order to reduce the increased memory requirements and also computational
time, Zeng (1996), independently Day (1998) and Day and Bradley (2001) intro-
duced coarse spatial sampling of the anelastic functions and coefficients. In Day’s
ij
(1998) approach, one anelastic function ξl for one relaxation frequency ωl is dis-
tributed with a spatial period of 2h, h being a grid spacing. Consequently, n = 8.
Considering, for example, location of the stress-tensor component T zx at 8 cor-
ners of a grid cube h × h × h, only one of the 8 ξlzx anelastic functions is assigned
to one of the 8 corners (say, ξ1zx is assigned to one position, ξ2zx to other position,
462 MOCZO ET AL.

and so on). Consequently, the total number of ξlzx , l = 1, 2, . . . , 8 in the whole

2 · 2 · 2 · 8 = MX · MY · MZ, MX, MY and MZ being the numbers
grid is MX MY MZ

of the grid cells in the three Cartesian directions. Because there are 6 independent
stress-tensor components, the total number of all the anelastic functions in the
μ
whole grid is MX · MY · MZ · 6. Since the anelastic coefficients Ylκ and Yl at the
μxy μyz μzx
grid positions of the normal stress-tensor components, and Yl , Yl and Yl
at the grid positions of the shear stress-tensor components are distributed in the
same coarse manner, the total number of the anelastic coefficients in the grid is
MX·MY ·MZ·5. Thus, the additional memory due to attenuation in Day’s (1998)
approach in the staggered-grid scheme in the case of 8 relaxation frequencies is
equivalent to the case of just one relaxation frequency without coarse sampling,
which is significant.
Graves and Day (2003) analyzed stability and accuracy of the scheme with the
coarse spatial sampling and defined the effective modulus and the quality factor
necessary to achieve sufficient accuracy.
As discussed earlier, Moczo et al. (2002) demonstrated that a position of a
material discontinuity within one grid cell can be sensed by a sufficiently accurate
FD scheme. In a structurally complex model there are material discontinuities
going through grid cells in different orientations with respect to the coordinate
system. In such a case and with the originally suggested spatial sampling (Day,
1998; Day and Bradley, 2001) it can happen that the medium from one side of
the material discontinuity is characterized over one half of the whole considered
frequency range while the medium from the other side of the discontinuity is
characterized over the other half of the considered frequency range. Since the
behavior of the two media in contact is characterized in two disjunctive frequency
sub-intervals, the two media cannot physically interact. Consequently, the two
media cannot be averaged.
In principle, the geometry of the coarse spatial sampling shown in the papers
by Day (1998) and Day and Bradley (2001) is not the only one possible. Keep-
ing the same spatial periodicity of the anelastic quantities, it is possible to avoid
division of a grid cell into two parts characterized in two disjunctive frequency
sub-intervals. Still the best possible alternative situation would be characteriza-
tion of one medium in contact using, for example, relaxation frequencies ω1 , ω3 ,
ω5 , ω7 and characterization of the other medium in contact using ω2 , ω4 , ω6 , ω8 ,
which again is not satisfactory.
In evaluating the sum that makes an anelastic term in the stress-strain relation,
Eqs. (128), (131) or (134), at a given spatial grid position, it would be possi-
ble to account for the anelastic coefficients and functions which are not located
at that grid point by their properly weighted values. Such averaging, however,
poses a problem: Because the anelastic functions (that is, internal variables or
memory variables) introduced by Day and Minster (1984), Emmerich and Korn
(1987), Carcione et al. (1988a, 1988b) and Robertsson et al. (1994) are material-
dependent, any such spatial averaging (accounting for the functions missing at the
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 463

F IG . 5. Coarse spatial distribution of grid cells and anelastic functions. The number on a cell face
indicates the relaxation frequency of the anelastic functions localized in the cell. For example, grid
yy xy yz
cell 1 contains ξ1xx , ξ1 , ξ1zz , ξ1 , ξ1 , ξ1zx . Reproduced from Kristek and Moczo (2003).

considered grid point) would introduce and additional artificial averaging of the
material parameters. There is no reason for such an additional averaging.
There would be no problem with the coarse spatial sampling and at the
same time with weighted spatial averaging of the anelastic functions at a grid
point with only one of the all anelastic functions if the anelastic functions were
material-independent. Therefore Kristek and Moczo (2003) introduced material-
independent anelastic functions (as given above). Moreover, they also suggested
an alternative coarse spatial distribution of the anelastic functions which only re-
quires n = 4 relaxation frequencies, keeping the same memory requirements
as in Day (1998), and Day and Bradley (2001). The distribution is shown in
Fig. 5. Kristek and Moczo (2003) demonstrated accuracy of their FD scheme
with the material-independent anelastic functions and new coarse distribution of
the anelastic functions.

3.2.10. Incorporation of Attenuation in Anisotropic Media

FD modeling in anisotropic media is addressed in the following section. Here

we only point out the main difference between incorporating attenuation in the
isotropic and anisotropic media. In the isotropic medium it is possible to con-
sider that the rheology of the medium is described by two separate viscoelastic
bodies—one for the complex frequency-dependent bulk modulus (corresponding
to the dilatational part of the strain) and one for the complex frequency-dependent
shear modulus (corresponding to the deviatoric part of the strain). In terms of the
quality factors for P- and S-waves, the quality factors can be strictly separated.
This makes the incorporation of the attenuation in the isotropic medium much
464 MOCZO ET AL.

easier compared to the anisotropic medium, where the two quality factors cannot
be simply separated if they are not equal.
Robertsson and Coates (1997) presented a 2-D velocity-stress staggered-grid
FD scheme for modeling qP- and qS-wave propagation in anisotropic media based
on the rheological model described by Carcione and Cavallini (1994). Through
eigenvalue decomposition of the stress and stiffness tensors, relaxation functions
and memory variables are associated with the so-called eigenstiffness and eigen-
stresses. The decomposition of stresses and stiffness in this fashion was first
observed by Lord Kelvin (Thomson, 1856).

4. A NISOTROPIC M EDIA

FD modeling of wave propagation in anisotropic media has attracted limited at-

tention as the main focus of both applied and academic studies has been on wave
propagation in isotropic heterogeneous models. However, a need for modeling of
anisotropic wave propagation arises not only from modeling of wave propaga-
tion in truly anisotropic materials such as for instance aligned anisotropic crystals
in the upper mantle or shales, but also from wave propagation in an equivalent
anisotropic medium, see Eq. (65). An equivalent anisotropic medium is a low-
frequency approximation for wave propagation in heterogeneous isotropic media
(Backus, 1962; Helbig, 1984). An equivalent anisotropic medium can be used to
represent for instance fine-scale layering or heterogeneity on a coarse scale appro-
priate for the scale of discretization of the FD modeling grid. Methods to compute
equivalent anisotropic media are well known in 1-D (Schoenberg and Muir, 1989)
and can applied to FD modeling of smooth interfaces in 2-D and 3-D (Muir et al.,
1992; Moczo et al., 2002). An equivalent anisotropic medium can also be used for
modeling of wave propagation through stress induced cracked materials, where
equivalent homogeneous anisotropic material represents average orientation and
size of cracks (e.g., Crampin and Chastin, 2003). A general equivalent medium
theory for 3-D heterogeneous media is a formidable challenge and still a topic of
research (and likely to remain one for quite some time).
Another reason for the limited use of FD modeling in anisotropic media is
the lack of suitable formulation of anisotropic finite differencing. Whereas the
conventional-grid schemes allow modeling of arbitrary anisotropic propagation,
they may become unstable at fluid/solid boundaries (among other problems). The
staggered-grid schemes are stable at fluid/solid boundaries but wave propagation
in general anisotropic media is significantly more difficult. For example, Hooke’s
law for the general anisotropic medium has the form

σij = cij kl uk,l .

 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 465

F IG . 6. Staggered grid in 2-D. Vx and Vz are the particle displacement/velocity components, Txx ,
Tzz and Tzx are the stress-tensor components.

Tensor of elastic constants cij kl is symmetric in i and j , k and l, ij and kl indexes.

The stress tensor σij requires evaluation of terms such as c1112 u1,2 . Figure 6 il-
lustrates that the derivative of u1 along the 2nd axis is not centered on the stress
component σ11 . In fact, in 3-D anisotropic media all derivatives with coefficients
c1112 , c1113 , c1123 , c2212 , c2213 , c2223 , c3312 , c3313 , c3323 , c2313 ,
c2312 , c1312
are not centered on the corresponding stress components with the standard
staggered-grid formulation. Therefore only wave propagation in an orthorhombic
anisotropic medium with axis of symmetry aligned along the Cartesian coordinate
system does not require any interpolation and is easily implemented. Neglecting
the shift of the derivatives for a generally anisotropic medium introduces signifi-
cant errors. Igel et al. (1995) therefore proposed interpolation operators that shift
the derivatives to the corresponding stress positions on the staggered grid. Anal-
ogously to the differential operators, the higher-order interpolation operators are
more accurate, but require more floating point operations. In fact the interpolation
significantly increases the number of operations at every time-step. In general
anisotropic media we must carry out 15 interpolations (12 for every coefficient
shown above and three of these coefficients require interpolations in two direc-
tions), which leads to a substantial increase in terms of computational cost. For
466 MOCZO ET AL.

example, the 3-D second-order staggered pressure-velocity scheme, with no at-

tenuation and no source representation, requires 16 floating-point operations per
time-step for acoustic media. Analogously, the partly-staggered grid (discussed
below) formulation (stress-velocity) requires 51 operations for an isotropic elastic
medium and 81 operations per time-step for an anisotropic medium (without the
cost interpolation required for staggered grids). The cost of interpolation operators
in anisotropic media leads to an increase in the number of operations by 500 or
more operations per time step (the number of float operations depends on particu-
lar implementation as some of the intermediate variables can be stored and reused
for further calculations; if no intermediate variables are reused, each interpola-
tion requires 55 operations for the fourth-order interpolation operator). The total
number of operations per grid points is therefore on the order of 600 for the 3-D
anisotropic case. In other words, anisotropic finite differences on staggered grids
are approximately 50 times slower (per time step) than acoustic modeling, and
more than 15–20 times slower than isotropic modeling. Note that in addition as
elastic simulations typically require shorter time steps and larger simulation times
to propagate S-waves thus making acoustic modeling comparably even more com-
putationally inexpensive.
A FD formulation on a partly-staggered grid, originally developed for isotropic
wave propagation (Andrews, 1973 and Zhang, 1997), was successfully applied
to anisotropic wave propagation by Saenger and Bohlen (2004) (they used term
rotated staggered). The formulation combines features of the conventional and
staggered-grid schemes. In a grid cell, all stress-tensor components are located at
one grid position, all displacement or particle-velocity components are located at
other grid position. Numerical analysis of the stability condition shows that the
formulation is stable for fluid/solid boundaries. For anisotropic finite-differences
on the partly-staggered grid the calculation is roughly only 60% slower than
isotropic modeling and the relative cost further decreases with the higher-order
space derivative operators that are used as both types of media require the same
number (18) of spatial derivatives (e.g., for the popular 4th-order in space and
2nd-order in time formulation, the anisotropic calculation requires only 35% ad-
ditional computations). This comparison illustrates that there is no significant sav-
ing in computational cost by modeling only restricted degrees of anisotropy on the
partly-staggered grid (e.g., modeling only orthorhombic versus fully anisotropic
media). Although the partly-staggered grid appears to have some quite attractive
properties it is not widely used today and the properties of the partly-staggered-
grid FD schemes, such as stability for various qP-to-qS wave speeds ratio and
interfaces, have to our knowledge not been fully analyzed and tested. Further
research is therefore required before adapting the partly-staggered grid formula-
tion instead of the standard staggered-grid formulation introduced by Madariaga
(1976) and Virieux (1984) that we focus on in our chapter.
Next we shall derive numerical stability conditions for a FD approximation that
is 2nd-order accurate in time and of arbitrary order accuracy in space (von Neu-
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 467

mann condition; O’Brien et al., 1951) for wave propagation in a fully anisotropic
homogeneous medium. We shall follow the derivation of numerical stability by
Crase (1990), Igel et al. (1995), and Saenger and Bohlen (2004). The stability
conditions for anisotropic heterogeneous media are evaluated by stability of the
equivalent homogeneous medium.
The velocity-stress formulation of the equation of motion for a general
anisotropic medium is (without the body-force term), Eqs. (8),
ρ v̇i = σij,j , σ̇ij = cij kl vk,l .
An approximation to the time derivatives yields
   
1 1 . 1
vi I, J, K, m + − vi I, J, K, m − = t ∂j σij (I, J, K, m) ,
2 2 ρ
  
. 1
σij (I, J, K, m) − σij (I, J, K, m − 1) = tcij kl ∂l vk I, J, K, m − .
2
(136)
Here, ∂i = ∂/∂xi . Taking the plane-wave ansatz, i.e., a harmonic wave ui (t, k) =
Ai exp[i(ωt − kj xj )], k = |k| being a wavenumber, which satisfies the elastody-
namic equation, we can rewrite Eq. (136) to contain only particle velocity at one
time:
 
ωt
−4 sin2 vi (I, J, K, m)
2
1  
= −t 2 ∂j cij kl ∂l vk (I, J, K, m) . (137)
ρ
The spatial operators on r.h.s. of Eq. (137) are also discretized. The ansatz solution
for the harmonic wave in a homogeneous anisotropic medium allows analytical
evaluation of the FD approximation to the real wavenumber (note that we ignore
the issue of interpolation here):
N /2  
. 2 S 2n − 1
∂i = k̃i = pn sin ki xi , (138)
xi 2
n=1

where k̃ is the numerical wavenumber, xi grid spacing in the xi -direction, NS

the order of the spatial operator, and pn the coefficients of the spatial operator.
Now we can evaluate complete discretized Eq. (137) for the ansatz solution and
obtain the following expression for the angular frequency:
 1/2
2 1 t 2
ω= arcsin λq (k̃r , cij kl ) . (139)
t 2 2
Here λp are diagonal elements of the matrix
1
Dik = k̃j cij kl k̃l . (140)
ρ
468 MOCZO ET AL.

If the expression under the square root in Eq. (139) is greater than one or less than
zero, the numerical angular frequency is not real and the FD scheme becomes
unstable (because the operator on particle velocity of Eq. (137) is unstable). Thus
a general condition for stability of the FD approximation is

t 2
0≤ λq (k̃r , cij kl ) ≤ 1. (141)
4
To evaluate stability condition (141) we need the eigenvalues of the matrix (140).
To evaluate these we bound the FD wavenumbers of Eq. (138) by
NS /2
2
k̃i ≤ |pn |, (142)
xi
n=1

and evaluate eigenvalues of the FD operator (140). In isotropic media the condi-
tions (139)–(140) can be evaluated analytically and they correspond to the CFL
condition (e.g., Mitchell and Griffiths, 1994, p. 167; Taflove and Hagness, 2005,
p. 42). The analytic solution exists for isotropic homogeneous media where sta-
bility condition (141) is
Nt /2
N /2 2n
2n 2n S
n t α
0 ≤ −2 (−1) |pj | ≤ 1, (143)
n=1
xi2n (2n)! j =1

where α is the P-wave velocity. The condition is valid for both P- and S-wave
velocities, but P-wave velocity imposes a stricter constraint as it is larger than the
S-wave velocity. However, for a general anisotropic medium the stability condi-
tion (141) has to be evaluated numerically.

5. F REE S URFACE

5.1. Traction-Free Boundary Condition

In exploration and earthquake seismology recordings are often made at or close
to the Earth’s surface. In most seismological applications, the air/fluid (ocean)
or air/solid (land) interface can be thought of as a free surface where a fluid or
solid abruptly terminates and is replaced by vacuum, that is, the Earth’s surface
may be approximated by a surface with vanishing traction. Let the surface S with
normal vector n define the geometry of the Earth’s surface. Let T (
u, n) be the
traction vector corresponding to the displacement vector u and normal vector n.
The traction-free condition is

T (
u, n) = 0. (144)
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 469

A real location for the seismic recordings can be affected by at least two signifi-
cant factors causing a large impact on wave propagation that must be modeled
or understood in many applications. Apart from the fact that the near-surface
medium structure may be highly complex with large velocity variations (e.g.,
Eisner and Clayton, 2002), the free surface itself may display significant topo-
graphic variations. At the same time, the FD method inherently has difficulty to
implement traction boundary condition. Therefore the accuracy and efficiency of
the free-surface approximation are key issues in FD modeling of seismic wave
propagation and earthquake motion.
Though, in principle, the same free-surface boundary condition should be sim-
ulated for planar and non-planar (topographic) free surfaces, due to the definition
of the FD method and unlike, for example, the finite-element method, the imple-
mentation of the traction-free condition in the case of topography is a considerably
more difficult problem. This is also reflected in the recent development. Therefore,
the two geometries will be addressed separately.

5.2. Planar Free Surface

Let the surface S be a planar free surface at z = 0 with the unit normal vector
n = (0, 0, −1). Then

σzη = 0, η ∈ {x, y, z} (145)

is the desired boundary condition at the horizontal (flat) free surface.
Recall the 4th-order approximation to the 1st derivative used in the staggered-
grid FD schemes for interior grid points. The 1st derivative of a function (ξ ) at
ξ = ξ0 is approximated by Eq. (5)

 (ξ0 ) = ,ξ (ξ0 )

     
1 3 3
= a  ξ0 + h −  ξ0 − h
h 2 2
    
1 1
+ b  ξ0 + h −  ξ0 − h
2 2
with a = −1/24, b = 9/8 and h being a grid spacing. Evaluation of the
z-derivative of a function at z = 0 then requires the function values at positions
−h/2 and −3h/2 above the free surface. Similarly, evaluation of the z-derivatives
of a function at z = h/2 and z = h requires values at z = −h and z = −h/2,
respectively. This implies two principal possibilities:
• the application of the FD scheme for the interior grid points with the field or
material parameters somehow defined above the free surface,
• the application of a different, say, adjusted, FD scheme which does not re-
quire any values above the free surface.
470 MOCZO ET AL.

Obviously, for the same reason which makes the heterogeneous schemes easier to
use than the homogeneous ones, the first approach is preferable.
The first approach either leads to the so-called vacuum formalism, medium ta-
per or imaging method. The vacuum formalism applies zero moduli above the free
surface. While this approach gives good accuracy in the displacement formula-
tion, e.g., Zahradník and Priolo (1995) in 2-D, Moczo et al. (1999) in 3-D, Graves
(1996) and other authors did not find the approach satisfactory in the staggered-
grid modeling. A density taper was used by Frankel and Leith (1992) in their
conventional-grid displacement FD scheme.

5.2.1. Stress Imaging

The stress imaging was introduced by Levander (1988) in his 2-D P-SV 4th-
order staggered-grid velocity-stress FD scheme. The stress-imaging technique
applies explicit boundary conditions to the stress-tensor component(s) located
at the grid plane coinciding with the free surface, and uses imaged values of
the stress-tensor components above the free surface assuming their antisymme-
try about the free surface. The antisymmetry
σzη (−z) = −σzη (z), η ∈ {x, y, z} (146)
ensures that the boundary condition given by Eqs. (145) are satisfied.
As summarized by Robertsson (1996), there are three possibilities for treating
the displacement or particle-velocity values formally required by the FD scheme:
1. The values are calculated using the 2nd-order approximations to the bound-
ary condition and imaged stress-tensor components. The approach was used
by Levander (1988), Graves (1996), Kristek et al. (2002) and others.
2. The values are mirrored as even values with respect to the free surface. The
approach was used by Crase (1990) and Rodrigues and Mora (1993). As
pointed out by Robertsson (1996), the even values of the particle velocity
values violate the boundary conditions.
3. The values are set to zero. This was proposed by Robertsson (1996).
Given the staggered grid, there are two natural options for locating the free sur-
face. In one, say H formulation, the horizontal displacement or particle-velocity
components, and stress-tensor components Txx , Tyy , Tzz and Txy are located at
the free surface. In the other, say W formulation, the vertical displacement or
particle-velocity component and Tzx and Tzy are at the free surface. Rodrigues
(1993) developed a 3-D 8th-order staggered-grid displacement-stress scheme and
used the stress-imaging technique in the H formulation. He found that it is nec-
essary to use more than twice the number of grid points compared to inside the
medium in order to avoid a significant numerical dispersion. Therefore, he com-
bined the stress-imaging technique with a vertically refined grid near the free
surface and achieved good accuracy. Kristek et al. (2002) numerically tested both
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 471

the H and W formulations of the stress-imaging technique against the discrete-

wavenumber method. They demonstrated that in the 3-D case the stress-imaging
technique in the 4th-order FD modeling requires at least twice as many grid points
per wavelength compared to what is sufficient inside the medium if the Rayleigh
waves are to be propagated without significant grid dispersion even in the case of
the simple homogeneous halfspace. They also tested the 4th-order version of the
Rodrigues (1993) approach. While sufficiently accurate, the approach needs three
times smaller time step (the factor of 3 is due to the most natural refinement of
the staggered-grid).
It is obvious that either at least twice denser spatial sampling or three times
smaller time step degrade the efficiency of the 4th-order staggered-grid modeling
inside the medium.

5.2.2. Adjusted FD Approximations (AFDA)

The principle of the AFDA technique used by Kristek et al. (2002) is simple.
The technique
1. Directly prescribes zero values of σzz at the free surface in the H formula-
tion or σzx and σzy in the W formulation.
2. Applies adjusted FD approximations to calculate the z-derivatives at the grid
points at the free surface and depths h/2 and h; the adjusted approximation
uses only function values in the medium.
As a consequence, no imaged (virtual) values above the free surface are needed.
Kristek et al. (2002) showed that while H-AFDA results in slightly better phases,
W-AFDA results in better amplitudes. They concluded with the recommenda-
tion to use W-AFDA for earthquake ground motion modeling. The calculation of
the stress-tensor and displacement components in W-AFDA can be summarized
as follows (if the velocity-stress formulation is considered, displacement compo-
nents are simply replaced by the particle-velocity components):
1. Direct application of the boundary condition:

Tzx (0) = 0, Tzy (0) = 0.

2. The following 4th-order FD approximations are used to calculate the stress-
tensor and displacement vector components:
    
 1 352 35 h 35 3
 (z0 ) = − (z0 ) +  z0 + −  z0 + h
h 105 8 2 24 2
   
21 5 5 7 
+  z0 + h −  z0 + h + O h4 , (147)
40 2 56 2
472 MOCZO ET AL.
      
1 11 h 17 h 3 3
 (z0 ) = −  z0 − +  z0 + +  z0 + h
h 12 2 24 2 8 2
   
5 5 1 7 
−  z0 + h +  z0 + h + O h4 , (148)
24 2 24 2
    
1 h 577 h 201 h
 (z0 ) = −  (z0 − h) −  z0 − +  z0 +
h 22   528  2 176 2
9 3 1 5  4
−  z0 + h +  z0 + h + O h ,
176 2 528 2
(149)
    
1 16 31 h 29 h
 (z0 ) = (z0 − h) −  z0 − +  z0 +
h 105   24  2 24
 2
3 3 1 5  4
−  z0 + h +  z0 + h + O h .
40 2 168 2
(150)
2.1. The calculation of the stress-tensor components:
Txx (h/2) is obtained from the 4th-order FD approximation to Hooke’s
law for σxx ; derivative uz,z is approximated by formula (148);
Tyy (h/2) and Tzz (h/2)—similar to Txx (h/2);
Tzx (h) is obtained from the 4th-order FD approximation to Hooke’s
law for σzx ; derivative ux,z is approximated by formula (149) in
which ux,z (0) is replaced by uz,x due to condition σzx (0) = 0;
Tzy (h) is obtained from the 4th-order FD approximation to Hooke’s
law for σzy ; derivative uy,z is approximated by formula (149) in
which uy,z (0) is replaced by uz,y due to condition σzy (0) = 0.
2.2. The calculation of the displacement-vector components:
W (0) is obtained from the 4th-order FD approximation to the equation
of motion for uz ; derivative σzz,z is approximated by formula (147)
in which condition σzz (0) = 0 is used;
U (h/2) is obtained from the 4th-order FD approximation to the equa-
tion of motion for ux ; derivative σzx,z is approximated by for-
mula (148);
V (h/2) is obtained from the 4th-order FD approximation to the equa-
tion of motion for uy ; derivative σzy,z is approximated by formula
(148);
W (h) is obtained from the 4th-order FD approximation to the equation
of motion for uz ; derivative σzz,z is approximated by formula (150)
in which condition σzz (0) = 0 is used.
In the W formulation, displacement component W , and stress-tensor components
σzx and σyz are located at the free surface. The corresponding grid material para-
meters are evaluated as integral averages in the half grid-cell volumes, that is, the
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 473

upper half of the volume located above the free surface is not taken into account.
For example, density and unrelaxed moduli are evaluated as
2 xI +1 yJ +1 z1/2
A
ρW = ρIA+1/2,J +1/2,0 = ρ dx dy dz, (151)
h3 xI yJ z0

 xI +1/2 yJ +1 z1/2 −1

2 1
μH
zx = μH
I,J +1/2,0 = 3 dx dy dz , (152)
h xI −1/2 yJ z0 μ

 xI +1 yJ +1/2 z1/2 −1

2 1
yz = μI +1/2,J,0 =
μH H
dx dy dz . (153)
h3 xI yJ −1/2 z0 μ
Using numerical comparisons against the DWN method, Kristek et al. (2002)
demonstrated that with the W-AFDA technique it is possible to apply the same
spatial sampling as inside the medium. Because in many realistic models of the
Earth’s interior lateral material discontinuities reach the free surface, Moczo et al.
(2004a) tested the accuracy of the W-AFDA technique against the finite-element
method for which it is easy and natural to satisfy boundary condition at the free
surface. Detailed numerical tests demonstrated the sufficient accuracy of the W-
AFDA technique in models with near-surface material discontinuities and the
capability of the FD scheme to ‘see’ the true position of the material disconti-
nuities in the spatial grid.

5.3. Free-Surface Topography

5.3.1. Approaches to Model Free-Surface Topography

The classification of approaches is similar to that given in the previous section.

The approaches based on modifying material properties at or in the vicinity of
the free surface to implicitly satisfy the boundary conditions (e.g., Mittet, 2002;
Zahradník and Hron, 1992; Zahradník et al., 1993) lend themselves to be general-
ized to incorporate topographic variations without much difficulty. However, they
tend to require significant spatial oversampling. Frankel and Leith (1992) mod-
eled the topographic effects on seismic waves generated at a Russian test site by
using a technique where a smoothly varying density taper was used to model the
transition from vacuum to the elastic sub-surface. Ohminato and Chouet (1997)
describe one of the first techniques implemented in 3-D where the exact location
of the free surface is chosen such that it follows a staircase approximating the
topographic surface. In their technique normal stresses are never located on the
free surface. Instead, the location of the surface is chosen such that only the shear
stresses which should be zero are at the free surface. The free-surface condition
is simulated by setting the shear modulus to zero at the free surface and all elastic
474 MOCZO ET AL.

moduli to zero in the vacuum above the free surface. In the 2nd-order FD scheme
at least 25 grid spacings per minimum wavelength are needed to achieve an ac-
curate and stable solution. Pitarka and Irikura (1996) also developed a similar
method which was applied in 3-D to study wave propagation site effects at the
Kobe-JMA station.
Robertsson (1996) presented a method of simulating the free surface in
staggered-grid FD modeling. The method is based on stress imaging and also
results in a staircase-shaped surface approximation to the free-surface topogra-
phy. By splitting the update at each time step in two iterations (in 2-D) that
contain spatial derivatives in one direction only, the imaging conditions can be
satisfied also in the case of surface topography. The method is relatively sim-
ple to implement and yields sufficiently accurate results if at least 15–20 grid
points are used per minimum wavelength. The method has been successfully
applied both to land seismic applications (e.g., Robertsson and Holliger, 1997;
Holliger and Robertsson, 1998) as well as for modeling scattering from a rough
sea surface (Laws and Kragh, 2002; Robertsson et al., 2006).
As mentioned above the problem with spatial oversampling the wavefield in
the vicinity of the free surface by roughly a factor of three can be circumvented
by introducing a simple grid-refinement scheme in the vicinity of the free surface
(Rodrigues, 1993; Robertsson and Holliger, 1997). A similar approach was also
taken by Hayashi et al. (2001). The technique by Robertsson (1996) and its exten-
sion by Robertsson and Holliger (1997) will be described in the next subsection.
Several authors developed approaches to avoid a staircase-shaped free surface.
Ilan (1977) considered an arbitrary polygonal free surface. Ilan’s technique did
not address the transition points between the segments of various slopes and re-
quired a non-uniform grid that decreased accuracy. An improved representation
of the arbitrary polygonal free surface was developed by Jih et al. (1988). Fol-
lowing a predefined classification scheme, the different segments were treated
using a one-sided approximation of the free-surface condition. Unfortunately, the
one-sided difference approximations reduce the accuracy of the method such that
a significant spatial oversampling is needed. Another approach to overcome the
“staircase problem”, initially developed by Fornberg (1988) and then further de-
veloped by a number of authors (Tessmer et al., 1992; Carcione and Wang, 1993;
Carcione, 1994; Nielsen et al., 1994; Hestholm and Ruud, 1994; Tessmer and
Kosloff, 1994; Hestholm, 1999; Hestholm and Ruud, 2002) is to solve the wave
equation on a curved grid whose line/surface coincides with the topographic sur-
face (or another internal surface). This is achieved by solving the equation of
motion written in Cartesian coordinates and involves first computing the spatial
derivatives in the new (conformably mapped) coordinate system (curved grid)
and then applying the chain rule to calculate the required Cartesian spatial deriv-
atives. The main drawbacks of the method are stability problems associated with
modeling very rough surface topography as well as the computational overhead
caused by the chain rule (25% in 2-D and 50% in 3-D according to Komatitsch
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 475

et al., 1996). Komatitsch et al. (1996) therefore proposed solving the equation
of motion directly on the curved grid and thus avoiding the chain-rule computa-
tions. Although the approach addresses the computational efficiency, it requires
additional memory.
As the FD method inherently has problems when incorporating the traction-
free condition, particularly in the case of an arbitrary surface topography, Moczo
et al. (1997) instead used a hybrid approach where a finite-element solution was
used in the vicinity of the free-surface topography and a FD solution in the rest
of the model. The method was shown to work very well and to be accurate and
stable. Obviously, in order not to loose the computational efficiency it is desirable
to minimize the region where the finite-element method is applied.

5.3.2. A Method for Modeling Surface Topography in 2-D

To summarize the discussion above, a completely accurate and efficient tech-

nique applicable in general to model free-surface topography does not exist.
Here we describe a relatively robust technique proposed by Robertsson (1996).
The method can be viewed as a generalization of the stress-imaging method of
Levander (1988) (H formulation) with one important modification. Instead of
updating the particle velocities in the vicinity of the free surface such that the
free surface condition is explicitly satisfied by using second-order accurate dif-
ference approximations, the particle velocities are simply set to zero above the
free surface. Stresses on the other hand are imaged such that tractions perpendic-
ular to the free surface always are zero at the free surface. The generalization to
the viscoelastic case is straightforward since no spatial derivatives of the memory
variables occur.
A description of the staggered grid in the vicinity of the free surface is shown
in Fig. 7, for a crest and a trough. The free surface is discretized such that all grey
grid-cells belong to it. It is critical where the boundary is located within the stag-
gered grid-cells (see Robertsson, 1996, for a discussion). Numerically, the free
surface itself is located along the thick black line. The grid-points along the dis-
cretized boundary belong to one and only one of the seven following categories:
(1) horizontal boundary (H);
(2) vertical boundary with vacuum to the left (VL);
(3) inner corner with vacuum above to the left (IL);
(4) outer corner with vacuum below to the left (OL);
(5) vertical boundary with vacuum to the right (VR);
(6) inner corner with vacuum above to the right (IR); and
(7) outer corner with vacuum below to the right (OR).
It is worthwhile to recall the physical meaning of the imaging technique. The
imaging is carried out to ensure that the normal and shear stresses perpendicular to
the boundary under consideration (σzz and σzx for the 2-D case of a flat horizontal
 476
MOCZO ET AL.
F IG . 7. Staggered FD grid in the vicinity of the free-surface boundary where it forms a crest or (left) and a trough (right). The light large squares represent
the locations of the grid-cells in the sub-surface. The grey grid-cells are located along the free-surface boundary, which runs exactly along the thick black line.
All boundary grid-points are classified H-, VL-, IL-, OL-, VR-, IR-, or OR-points as described in the text. Within the grid-cells, the solid squares represent the
σxx and σzz components, the light squares the σzx component, the solid circles the vx component, and the light circles represent the vz component. Reproduced
from Robertsson (1996).
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 477

surface) are zero. Notice that within all seven categories of boundary grid-cells,
the free surface is always parallel to the grid. Imaging, therefore, only takes place
in the directions of the x- and z-coordinate axes.
What happens in the vicinity of inner corners, where grid-points are subject to
imaging from different directions? Briefly, by calculating separately the vertical
and horizontal derivatives of the stress components in the update of the particle
velocities, such problems are readily avoided (see below).
Horizontal boundary grid-point (H-point): An H-point has one neighbor on ei-
ther horizontal side that is either an H-, OL-, OR-, IL- or IR-point. The H-points
are treated identically to the flat free-surface approximation. Particle velocities
are set to zero in the vacuum. Imaging of stresses only takes place in the vertical
direction. Just as for the flat free surface, the σxx component is updated using the
fourth-order accurate central FD approximation along the surface.
Vertical boundary grid-point with vacuum to the left (VL-point): A VL-point
has one neighbor on either vertical side that is either a VL-, IL- or OL-point. The
VL-points are treated similar to the H-points with the exception that here it is the
σxx component that should be zero at the boundary, and the imaging therefore
only takes place in the horizontal direction. The σzz component is updated in the
same way that the σxx component is updated for the H-points.
Inner corner grid-point with vacuum above to the left (IL-point): An IL-point
has one neighbor vertically above that is either an OL- or VL-point, and one
horizontally to the left that is either an H- or OL-point. No imaging takes place
around the IL-points. The σxx and σzz components are located along the boundary,
perpendicular to parts of it, and are therefore set to zero.
Outer corner grid-point with vacuum to the left (OL-point): An OL-point has
either an IL- or H-point immediately to the right and an IL- or VL-point immedi-
ately below. The free surface is located through the vx , vz , and σzx components,
while the normal stresses are located in the vacuum above the free surface (Fig. 7).
The σzx component is therefore set to zero. The σxx and σzx components are im-
aged horizontally with respect to the rightmost vertical part of the free surface,
and the σzz and again the σzx components are imaged vertically with respect to
the lowermost horizontal part of the free-surface (Fig. 7). If the OL-point is ad-
jacent to an IL-point, the particle velocity component in between must be set to
zero to obtain a stable and accurate solution.
Vertical boundary grid-point with vacuum to the right (VR-point): A VR-point
has one neighbor on either vertical side that is either a VR-, IR- or OR-point. The
VR-points are treated analogously to the VL-points, with the exception that the
imaging takes place in the opposite horizontal direction.
Inner corner grid-point with vacuum above to the right (IR-point): An IR-point
has one neighbor vertically above that is either an OR- or VR-point, and one
horizontally to the right that is either an H- or OR-point. Only vertical imaging
of the σzx component takes place for the IR-points. The σxx and σzz components
478 MOCZO ET AL.

are located along the boundary, perpendicular to parts of it, and are therefore set
to zero.
Outer corner grid-point with vacuum to the right (OR-point): An OR-point has
either an IR- or H-point immediately to the left and an IR- or VR-point im-
mediately below. The free-surface makes a step immediately to the left of the
OR-point and only intersects with its vz component, whereas the other velocity-
and stress-fields are located in the vacuum above the free surface (Fig. 7). The
σzx component immediately to the left of the boundary point is set to zero. The
σxx and σzx components are imaged horizontally with respect to the vertical seg-
ment of the free-surface boundary to the left of the grid-point (Fig. 7). The σzz
component is imaged vertically with respect to the step at the vz component in the
OR-point. If the OR-point is adjacent to an IR-point, the particle velocity compo-
nent in between must be set to zero to obtain a stable and accurate solution.
The free-surface boundary points have to be classified prior the FD. The max-
imum computational efficiency is obtained by setting the reciprocal values of the
densities to zero everywhere in the FD grid above the free surface. The particle ve-
locities will then automatically be zero above the free surface after every update.
The imaging algorithm thus does not add a substantial amount of the computa-
tional cost. However, computations are wasted in the region occupied by vacuum
above the topography in the FD grid.
Each of the equations for the particle velocities, consist of two derivatives (one
vertical and one horizontal) of the stress-tensor components. Here we do not show
explicitly equations for the considered 2-D case as they can be easily obtained
from Eqs. (8) for the 3-D case. By first performing all vertical imaging of the
stress-tensor components (H-, OR-, OL-, and IR-points) and then calculating and
adding only the vertical derivatives in equations, the free-surface condition is sat-
isfied completely in all grid-cells along the boundary. Next, the horizontal imaging
of the stress components (VL-, VR-, OR-, and OL-points) is carried out followed
by an update with the remaining horizontal derivatives in equations, again sat-
isfying the free-surface boundary condition. Following the updates of equations
it is necessary to set the vx component to zero in the VR- and OR-points, since
the reciprocal velocity is not zero in these points. Subsequently, the stress-fields
are updated. This action does not involve any imaging of the variable fields. Fol-
lowing the update, the correction of the normal stress parallel to the surface must
be made at the H- (σxx ), VL- and VR-points (σzz ), as described above for the
horizontal free-surface.
There is a modeling limit as to how narrow the “troughs” or “crests” in the
topography can be. A crest cannot be narrower than the number of grid-points
imaged around its horizontal sides. Fig. 7 shows the narrowest crest and trough
that are allowed when using fourth-order accurate spatial central-difference ap-
proximations.
Since the technique requires approximately three times dense spatial sam-
pling than that inside the medium, Robertsson and Holliger (1997) used a grid-
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 479

F IG . 8. Grid-refinement technique (Robertsson and Holliger, 1997) used together with the method
for modeling surface topography (Robertsson, 1996). The upper part of the grid (light grid cells) is
three times more densely sampled than the lower part (dark grid cells). Reproduced from Robertsson
and Holliger (1997).

refinement as illustrated in Fig. 8. For the update at grid points of the finer grid
close to the coarser grid a simple linear interpolation in the vertical direction is
sufficient between the wavefield components in the coarse and fine parts of the
grid. However, in the horizontal direction, some care should be taken since waves
that propagate close to parallel to the transition region will undergo reflection due
to very slight numerical errors in the interpolation between the two regions. An
excellent performance can be obtained by using a more accurate sinc interpolation
(Martin Musil, personal communication).
To illustrate the accuracy of the method for modeling surface topography we
show an acoustic example that has been published elsewhere (Robertsson et al.,
2006). The example addresses the problem of modeling scattering from a rough
sea surface modeled as a self-similar surface using the method by Pierson and
Moskowitz (1964). The so-called Significant Wave Height (SWH) was set to 2
meters (a moderately rough sea as marine reflection seismic data typically are
acquired in seas up to 4 m SWH).
Figure 9 shows the reflected response due to a plane wave incident vertically
from below as well as at 30◦ incidence angle with respect to the vertical, recorded
6 m below the average level of the sea surface. We show solutions computed using
three different methods. First, the FDTD method by Robertsson (1996) described
here. Second, a solution computed using the spectral element method (Chaljub
et al., 2006) and finally a solution computed using the Kirchhoff approximation
(Laws and Kragh, 2002).
480 MOCZO ET AL.

F IG . 9. Rough sea surface reflected response for a plane wave incident vertically from below (top
row) and at an angle of 30◦ to the vertical (bottom row). The reflected/scattered field is recorded at a
point placed at 6 m depth. The plots on the right are enlargements of the scattered coda in the plots
on the left. Black: FDTD response. Grey: Response computed using a Kirchhoff method (Laws and
Kragh, 2002). Dashed: Response computed using a spectral element method (Chaljub et al., 2006).
Figure reproduced from Robertsson et al. (2006).

Careful convergence tests were carried out to ensure that details of the reflected
and scattered response were not contaminated by numerical artifacts in the re-
spective method. The FDTD method determined the spatial discretization of the
sea surface used for all three methods, since it required the densest sampled sea
surface (15–20 grid-points per minimum wavelength).
As can be seen from Fig. 9, there is very good match between the FDTD and
the spectral element solutions both in terms of amplitude and phase far into the
coda. Although the much simpler (and more efficient Kirchhoff method) provides
a reasonably accurate solution, the synthesized coda differs significantly from that
computed using the two other methods.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 481

6. WAVEFIELD E XCITATION

6.1. Direct Modeling of the Point Sources

The average properties of complex seismic sources are usually represented as
point sources in a continuous elastic medium (see, for example, Chapter 3 of Aki
and Richards, 1980). The two simplest point sources are body forces and moment-
tensor sources. A point source representation of a body-force type of source can
be implemented directly as the increment of the corresponding components as
given by the equation of motion, e.g. Eqs. (7):

ρ v̇i = σij,j + fi .

Note, that the body force source-time function is not differentiated with respect
to time. The usual implementation of the body-force point source at the time tm is
then
1
vi (tm ) = vi (tm ) + tfi (tm ).
ρ
Implementation of the moment-tensor source can be included either by stress
(Virieux, 1986; Coutant et al., 1995) or by particle velocity (Frankel, 1993;
Yomogida and Etgen, 1993; Graves, 1996). The implementations are equivalent
due to the body-force equivalent theorem (e.g., pages 40–44 of Aki and Richards,
1980). In the particle-velocity implementation, each component of the moment
tensor is implemented by corresponding couple of the body forces with discrete
arm length between the forces. For example, in the staggered-grid velocity-stress
formulation the Mxx component of the moment tensor is equivalent to couple of
the forces (Mxx /x) acting along the x-axis in the opposite directions. Because
these forces are applied at one grid point, the appropriate volume is one grid cell,
x · y · z. Therefore, the particle-velocity update for the Mxx component of
the moment tensor at grid node (I, J, K) is
   
1 1 1 Mxx (I, J, K, m)
vx I + , J, K, m = vx I + , J, K, m + t ,
2 2 ρ x 2 yz

   
1 1 1 Mxx (I, J, K, m)
vx I − , J, K, m = vx I − , J, K, m − t .
2 2 ρ x 2 yz
However, the equivalent body forces for the representation of the Mxy component
of the moment tensor are not located along the grid line I and they must be aver-
aged from four equivalent body forces (Mxy /2y) acting along the x-axis in the
opposite directions with a force arm of length 2y as illustrated in Fig. 10. There-
fore the particle velocity update for the Mxy component of the moment tensor at
482 MOCZO ET AL.

F IG . 10. An example of equivalent body force implementation of moment tensor source compo-
nent Mxy (top) and Mxx (bottom). The FD grid shows only the XY plane of the 3-D grid as all body
forces are implemented in the same plane (plane K).

grid node (I, J, K) is

 
1
vx I ± , J + 1, K, m
2
 
1 1 Mxy (I, J, K, m)
= vx I ± , J + 1, K, m + t ,
2 ρ 4xy 2 z
 
1
vx I ± , J − 1, K, m
2
 
1 1 Mxy (I, J, K, m)
= vx I ± , J − 1, K, m − t .
2 ρ 4xy 2 z
Analogously the remaining components of the moment tensor can be imple-
mented as equivalent body forces centered at grid node (I, J, K). Note that this
moment tensor implementation allows modeling of the explosive sources (equal
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 483

diagonal elements and vanishing non-diagonal elements of the moment tensor),

pure shear slip (if the moment tensor is determined from the dip, strike and rake
and seismic moment, pages 117–118 of Aki and Richards, 1980), or a compen-
sated linear vector dipole (CLVD) source.

6.2. Introducing the Source Wavefield Along a Boundary Internal to the Grid
Alterman and Karal (1968), Kelly et al. (1976) and Levander (1989) describe
how to introduce a source field into a FD grid by “injecting” an analytical source
solution on an internal artificial surface surrounding the source to for instance
avoid a (point) source singularity in the FD computation. The source wavefield
is introduced so that it radiates from the outside of the surface surrounding the
source. In principle, all we need is to satisfy the principle of superposition and the
continuity of the superimposed wavefields across the surface. In practice because
the FD calculation is discontinuous at the injection surface S, some care must be
taken in the FD calculations where the spatial FD stencil intersects the surface.
Figure 11 shows a 2-D staggered FD grid in the vicinity of the injection sur-
face S. The region outside the surface is referred to as the external region Ve

F IG . 11. Introduction of a source wavefield along an internal artificial surface in the staggered grid
illustrating the update of normal stresses. Normal stresses are located on the solid grid whereas shear
stresses are located on the dashed grid. The region with grey background corresponds to Vi whereas
the remaining part of the grid (white background) corresponds to Ve . The so-called injection surface
S separates the two regions. Two FD stencils are shown in the picture (one updating a point in Vi and
one a point in Ve ). Bold lines in the grey region correspond to grid points where the source wavefield
should be subtracted at the appropriate points during the update whereas bold lines in the white region
show grid points where it should be added. Reproduced from Robertsson and Chapman (2000).
484 MOCZO ET AL.

(white in Fig. 11) whereas the region inside the surface is referred to as the inter-
nal region Vi (grey in Fig. 11). The process of introducing the source wavefield
is straightforward and involves manipulating the update of the wavefield only in
the vicinity of S for the update of points where the spatial extent of the FD stencil
intersects S (where the wavefield in the grid is discontinuous). In Fig. 11 these
points have been marked with by the thick bold lines.
The source field to be injected is known (e.g., analytically) in the vicinity of
the surface S for all times (this is possible at least if the surface is located en-
tirely within a homogeneous part of the model). The source wavefield will interact
with the model throughout the FD grid within Vi and Ve and generate a scattered
wavefield. In Ve both the source wavefield and the scattered wavefield is present
whereas in Vi only the scattered wavefield is present. Since the source wavefield
is known (e.g., analytically) this can be added and subtracted as appropriate for
the update of points along the bold lines for the parts of the spatial stencils that
intersect S (subtract the source wavefield at the appropriate points when updating
points along the bold lines in the grey region and adding the source wavefield
at the appropriate points when updating points along the bold lines in the white
region).
In each FD time step the stresses (and corresponding memory variables for the
viscoelastic case) are first updated in the entire grid. When the update is complete,
we go back and correct the update at the points where the spatial FD stencil in-
tersected the surface S. Here we describe how this is done for the normal stresses
(illustrated in Fig. 11). Outside S in Ve , the wavefield is updated as if the injected
wavefield were propagating through the entire grid. Therefore we must add the
injected source wavefield to the particle-velocity components corresponding to
the parts of the stencil that are inside S in Vi . In Fig. 11 this occurs when normal
stresses along the bold lines in the white region along each side of Ve are updated.
For a 4th-order accurate scheme we therefore need to know the source wavefield
along the two closest grid points inside S for the upper and left edges of the rec-
tangular injection region shown in Fig. 11, and along the closest grid point for the
lower and right edges. Inside S in Vi , the wavefield is updated as if no wavefield
were injected. Therefore we must subtract the injected source wavefield from the
particle-velocity components corresponding to the parts of the stencil that are out-
side S in Ve (for update of normal stresses located along the bold lines in the grey
region along each side in Vi ).
Next, we advance the calculation by half a FD time step and update the particle
velocities in the entire grid using equation. The wavefield injection is performed
using the same procedure by adding and subtracting the stress components of the
injection wavefield at the three grid points around S. In total, we therefore need to
know the values at all times of stresses and particle velocities at three grid points
around S, staggered appropriately both in time and space. By iterating these two
steps of the update, the entire FD simulation is stepped through and the wavefield
is injected along the surface S.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 485

F IG . 12. Example of a hybrid simulation in a homogeneous medium, where the FD computation is

driven by an analytical wavefield for a point source in a homogeneous medium. The source is located
at a fictitious position in the upper corner of the FD grid. The injection surface is outlined as a dashed
box. Snapshots from 0.14 s, 0.24 s and 0.30 s are shown.
486 MOCZO ET AL.

As noted by for instance, Fäh (1992), Zahradník (1995), Robertsson et al.

(1996), Zahradník and Moczo (1996), Moczo et al. (1997), Robertsson and Chap-
man (2000), and Takeuchi and Geller (2003), the same technique for introducing
the source wavefield along an internal boundary can be used for hybrid model-
ing purposes where the wavefield from another computational technique (e.g.,
discrete-wavenumber or ray method) is introduced inside a FD grid. In that case
we reverse the regions Vi and Ve so that the source wavefield is injected inside the
surface S.
Figure 12 shows snapshots from tests of such a hybrid technique. Here the
source wavefield is given by an analytical solution for a point source in a homo-
geneous medium. Since the FD grid also corresponds to a homogeneous medium
we would expect the scattered wavefield in Ve to be zero. In the snapshots from
0.14 s and 0.24 s, we see no evidence of any wavefield being present in Ve . How-
ever, in the snapshot from 0.30 s, part of the wavefield leaks through into Ve . This
is because the wavefront that has propagated in the FD grid has suffered from nu-
merical dispersion and no longer exactly match the wavefield from the analytical
solution which ideally should destructively interfere with the FD wavefront as it
reaches the injection surfaced S.
In the so-called FD-injection technique (Robertsson and Chapman, 2000), the
source wavefield is generated by a FD solution in an unperturbed model to drive
the update on small FD sub-grids surrounding regions of change to compute the
wavefield in a perturbed model. This can provide a powerful method in for in-
stance waveform inversion applications.

6.3. Dynamic Modeling of Earthquake Rupture

6.3.1. Fault Boundary Condition

In many seismological problems an earthquake fault may be represented by a

surface embedded in heterogeneous elastic or viscoelastic pre-stressed medium.
A non-zero initial equilibrium stress is due to tectonic loading and residual stress
after previous earthquakes on the fault. An earthquake itself may be modeled as
spontaneous rupture propagation along the fault. The rupture generates seismic
waves which then propagate from the fault into the embedding medium. In gen-
eral, several ruptures can propagate along the fault at one time. Inside the rupture
displacement and particle-velocity vectors are discontinuous across the fault. At
the same time traction is continuous. Let n(xi ) be a unit normal vector to the fault
surface pointing from the ‘−’ to ‘+’ side of the surface (Fig. 13), D u(xi , t) =
u+ (xi+ , t)− u− (xi− , t) slip (discontinuity in displacement vector across the fault),
D v(xi , t) = v+ (xi+ , t)− v− (xi− , t) slip rate (discontinuity in the particle-velocity
vector across the fault), T ( n; xi , t) = T 0 (
n; xi ) + T ( n; xi , t) total traction on

the fault, T (0 n; xi ) initial traction, and T ( n; xi , t) traction variation. The latter
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 487

F IG . 13. Fault surface and the normal vector n.

is due to the rupture propagation. Inside the rupture the total traction is related to
slip at the same point through the friction law T = T f (D u, D v, θ ), where T f is
frictional traction and θ represents a set of state variables. Given the initial trac-
tion and (visco)elastic material parameters of the fault, it is the friction law which
controls initialization, propagation and healing of the rupture.
Consider further only shear faulting: there is no opening of the fault and no
interpenetrating of the fault materials. Assume a locked fault. If the magnitude
of the shear traction is smaller than the frictional strength (product of the normal
traction and the friction coefficient) at a point, the fault remains locked and slip
rate zero at the point. Should the shear traction exceed the frictional strength, slip
occurs. The shear traction decreases over a finite distance from the static down to
the dynamic frictional level, following the friction law. The slipping is opposed
by the friction.
Let subscripts sh and n denote the shear and normal components. Let S denote
frictional strength. The boundary conditions on the fault can be formulated as
follows (Day, 1982, 2005; Day et al., 2005).
Shear faulting:
D un = 0, D vn = 0, D ush = 0, D vsh = 0. (154)
Shear traction bounded by the frictional strength:
|Tsh | ≤ S. (155)
Colinearity of the shear traction and slip rate:
SD vsh − Tsh (
n)|D vsh | = 0. (156)
The fact that the frictional traction opposes the slipping is consistent with the
colinearity requirement because we consider vector n oriented in the direction
from the ‘−’ to ‘+’ side of the fault and slip as the relative motion of the ‘+’ side
with respect to the ‘−’ side of the fault: both T (
n) and D v are viewed from the
same side of the fault. If slip was defined as the relative motion of the ‘−’ side
with respect to the ‘+’ side of the fault, requirement of the antiparallelism with
the ‘+’ sign in Eq. (156) would be consistent with the frictional traction opposing
the relative motion of the fault faces.
While semi-analytical boundary integral equation (BIE) method is perhaps the
most accurate method to account for the fault boundary conditions, especially
488 MOCZO ET AL.

F IG . 14. Halfspaces H − and H + , partial nodes p.n.+ and p.n.− , and the normal vector n.

on non-planar faults (e.g., Aochi and Fukuyama, 2002), its application is limited
because it cannot include heterogeneity of the medium. Because the FDM can
account for material heterogeneity and is computationally more efficient, it has
been extensively applied to study source dynamics particularly on planar faults
parallel to grid planes.
The FDM has been applied to the dynamic rupture propagation independently
by Andrews (1973, 1976), and Madariaga (1976), and then by many others. Two
main approaches have been developed and applied. In the split-node approach the
fault is represented by a grid surface of split (partial) nodes. At a grid point, each
of the two partial nodes belongs to only one side of the fault and the two nodes
may experience a relative motion (slip) along the fault. In the zone approach the
fault is represented by a “thick” zone and slip is evaluated either using inelastic
strain or as a difference between displacements at grid points separated by the
fault zone. Before we mention other approaches we explain the split-node and
zone methods in some detail.

6.3.2. Traction-at-Split-Nodes (TSN) Method

The TSN method has been developed independently by (Andrews (1973, 1976,
1999); Day (1977, 1982, 2005)) see also Day et al. (2005). Consider a halfspace
H − covered by a FD grid and a grid node p.n.− on the free surface of the half-
space. Similarly, consider halfspace H + and a grid node p.n.+ on its free surface
(Fig. 14). Define an outer normal vector n to the surface of the halfspace H −
pointing to the halfspace H + (i.e., n is in the ‘p.n.− → p.n.+ ’ direction).
Let M − and M + be masses of the two partial nodes. The partial node p.n.−
is accelerated by a force F − which is due to deformation in the halfspace H −
and, possibly, by body forces acting in the halfspace. Similarly, the partial node
p.n.+ is accelerated by a force F + . Thus, the accelerations are a ± = F ± /M ± .
Couple the halfspaces H − and H + along their surfaces in order to simulate a fault.
The coupling can be accomplished by a constraint surface traction acting at the
contact. Consider a traction T c (
n) quantifying a contact force with which material
in H + acts upon material in H − . Let A be an area of the fault surface associated
with each partial node. Then the acceleration a − of the partial node p.n.− is
contributed by the force F − (which is due to deformation in the halfspace H − )
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 489

and by the constraint force F c,− = F c = A · T c ( n) (which is due to the action of

the halfspace H + ): a − = (F − + F c,− )/M − = (F − + A · T c )/M − . Similarly,
the acceleration a + of the partial node p.n.+ is contributed by the force F + (due
to deformation in the halfspace H + ) and by the constraint force F c,+ = −F c =
−A · T c ( n) (due to the action of the halfspace H − ): a + = (F + + F c,+ )/M + =
(F − A · T c )/M + . Consider some initial equilibrium state with traction T 0 (
 + n).
The traction does not contribute to the acceleration of the partial node p.n.− . If
T c is the total traction, only the difference T c − T 0 contributes to the acceleration.
Considering the accelerations at time t,
 
a ± (t) = F ± (t) ∓ A · T c (t) − T 0 /M ± . (157)
Though the initial traction is nonzero, the initial strain is considered zero. Then
forces F ± correspond to deformations caused only by the dynamic changes due
to rupture. The particle velocities and displacements of the partial nodes in the
2nd-order approximation are then
   
dt dt dt{F ± (t) ∓ A · [T c (t) − T 0 ]}
v± t + = v± t − + (158)
2 2 M±
and
 
dt
u± (t + dt) = u± (t) + dt · v± t + . (159)
2
For the slip rate we obtain from Eq. (158)
 
dt
D v t +
2
   − + 
. dt M F (t) − M + F − (t)  c (t) + T 0 ,
= D v t − + dt B − T
2 A · (M − + M + )
(160)
where B = A(M − + M + )/(M − M + ). Find a constraint traction T c (t) = T ct (t)
that assures zero slip rate before two partial nodes start slipping as well as van-
ishing slip rate when the slipping ceases. The question is how to time condition
D v = 0. If D v(t) = 0 is required, the trial traction acts for the interval from
t − dt/2 to t + dt/2 and can reverse the slipping (i.e., produce back-slip) by the
time it is integrated all the way up to t + dt/2. This results in the traction driving
slip rather than opposing it and thus in violating conservation of energy. There-
fore, D v(t +dt/2) = 0 has to be required. Assume D v(t +dt/2) = 0 in Eq. (160)
and obtain the trial traction

. dt −1 M − M + D v(t − dt/2) + M − F + (t) − M + F − (t)

T ct (t) = T 0 + .
A · (M − + M + )
(161)
490 MOCZO ET AL.

Find a constraint traction during the slip, that is, frictional traction Tsh
c (t) = T f (t)
sh
such that D v(t +dt/2) = 0. Assuming first D v(t +dt/2) = 0 for Tsh c (t) = T f (t)
sh
in Eq. (160), and then D v(t + dt/2) = 0 for T c (t) = T ct (t) leads to
 
dt .
= dt B Tsh (t) − Tsh (t) .
f
D vsh t + ct
(162)
2
Recalling the colinearity condition (156)
 
S(t)D vsh (t) − Tsh (t)D vsh (t) = 0.
f
(163)
Using approximation D vsh (t) = 12 [D vsh (t − dt/2) + D vsh (t + dt/2)] and
Eq. (162) we obtain from the colinearity (163)
     
 
D vsh (t) + S(t) dt B T f (t) = . S(t)
D vsh t −
dt 
+ dt B Tsh (t) .
ct
sh
2 2 2
(164)
Define an auxiliary vector γ
 
dt
γ = D vsh t − + dt B Tsh
ct
(t). (165)
2
Equations (164) and (165) imply that Tsh (t) has the direction of the vector
f

ϒ = γ /|γ |. Therefore, enforcing the boundary conditions on the fault can be

formulated as follows:
 ct 
If Tsh (t) ≤ S(t) then T c (t) = T ct (t). (166)
 ct 
If Tsh (t) > S(t) then Tsh (t) = S(t)ϒ,
c  Tn (t) = Tn (t).
c ct
(167)
The above approach, based on finding a trial traction T ct (t) ensuring D v(t +
dt/2) = 0, and colinearity requirement at time t, can cause in some rare cases
large oscillations of rake direction just around the time of rupture arrest (Day,
2005). It is possible to avoid this problem by the modification of the colinearity
condition (Day, 2005):
    
dt  dt  .
− Tsh (t)D vsh t +
f
S(t)D vsh t + = 0. (168)
2 2 
Inserting Eq. (162) into Eq. (168) yields
 ct  f .
S(t) + Tsh (t) − Tsh (t) Tsh (t) = S(t)Tsh
f ct
(t). (169)

Equation (169) means that Tsh (t) has the same direction as T ct (t). Then condition
f

(167) is replaced by
   
If T ct (t) > S(t) then T c (t) = S(t)T ct (t)/T ct (t), T c (t) = T ct (t).
sh sh sh sh n n
(170)
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 491

The modified approach behaves always well.

An assumption of the small displacements is necessary for the TSN method.
The assumption means that the accumulating slip does not change the configu-
ration of the partial nodes adjacent to each other. This means that h  |D ush |,
where h is a spatial grid spacing. Also, a necessary condition is that the time-
stepping algorithm is explicit and a force at a node accelerates only that node.
Accuracy of the TSN implementation heavily depends on the accuracy of cal-
culation of the body forces F ± due to deformations in the halfspaces. Formally,
at each time the surfaces of the halfspaces are the free surfaces. In fact, both
Andrews (1973, 1976, 1999) and Day (1982) used an FD formulation on the
partly-staggered grids in which the spatial differencing is equivalent to the finite-
element method.

6.3.3. Stress Glut (SG) Method

Andrews (1999) describes the SG method that he used in his 1976 paper
(Andrews, 1976). Instead of a surface with an explicit displacement disconti-
nuity in the TSN method, an inelastic zone of a finite thickness is used in the
SG method to represent the fault surface. The method can be implemented in
the partly-staggered or staggered grid. Consider for simplicity a horizontal zone
(perpendicular to the z-axis oriented positive downward) centered at a grid plane
with grid positions of the stress tensor and bounded in the vertical direction by
grid planes with grid positions of the displacement and particle-velocity vectors.
The thickness of the fault zone is thus equal to one grid spacing h. Assume the
velocity-stress formulation. At each time level stresses are updated as if there
were no fault zone. Denote shear stress-tensor components located at the central
horizontal grid plane as trial values (superscript t) and evaluate a magnitude of the
t | = [(τ t )2 +(τ t )2 ]1/2 . If, at a grid point, |T t | ≤ S,
shear trial traction vector |Tsh zx zy sh
S being the frictional strength of the fault, then τzi = τzit , i ∈ {x, y} at that grid
point. If |Tsh
t | > S, then, due to the colinearity, τ = (S/|T t |)τ t , i ∈ {x, y}.
zi sh zi
The latter means an inelastic stress adjustment. Assuming that the corresponding
inelastic strain is distributed over the fault thickness, the offset in stress, i.e., stress
glut, is related to the increment of the seismic moment tensor in the volume of the
grid cell:

Mzi = τzit − τzi · h3 , i ∈ {x, y}. (171)
At the same time, the increment in the seismic moment tensor can be interpreted
under the assumption of a shear slip along the central grid plane of the fault zone.
Assuming a normal vector n(0, 0, 1) and slip rate at a given grid position on the
central grid plane Dvi = vi+ (z+ ) − vi− (z− ),

Mzi = μDvi h2 t, i ∈ {x, y}. (172)

492 MOCZO ET AL.

Equating Eqs. (171) and (172) we obtain for the slip rate

h τzit − τzi
Dvi = , i ∈ {x, y}. (173)
t μ
Andrews (1999) discusses computational details as well as the more complicated
implementation of the method on the staggered grid.

6.3.4. Thick-Fault-Zone (TFZ) Method

Madariaga et al. (1998) described a velocity-stress staggered-grid method to

study dynamic faulting. As in the SG method, no surface with explicit displace-
ment and particle-velocity discontinuity is considered. Instead, an alternative fault
zone of the finite thickness is centered at the grid plane with grid positions of the
normal stress-tensor components. Considering again, e.g., a horizontal fault zone,
frictional boundary conditions are applied on the two nearest grid planes with grid
positions of the shear stress-tensor components, i.e., on the grid planes half grid
spacing below and above the central grid plane. Because the horizontal displace-
ment components are located on the central grid plane, the slip rate is evaluated
as the difference between particle-velocity values at grid planes one grid spacing
above and one grid spacing below the central grid plane, that is, over a thickness
of two grid spacings. Such a configuration preserves symmetry of stresses and
particle velocities about the fault plane. Madariaga et al. (1998) apply standard
4th-order FD formulae to calculate spatial derivatives at all grid positions.

6.3.5. Comparison of the TSN, SG and TFZ Methods

Dalguer and Day (2004, 2006) performed an extensive numerical comparison

of the three methods. Calculations were performed with the 2nd-order TSN and
4th-order staggered-grid SG and TFZ formulations using the same grid spacing
h and uniform grid. The rupture propagation velocity in the fault-zone model is
lower than that in the split model—likely due to the blunting of the stress concen-
tration on the rupture front in the fault-zone models. The TSN solution converges
substantially better than the SG solution, and far better than the TFZ solution.
The SG method reaches convergence with a grid spacing h/2, while grid spacing
smaller than h/4 is necessary for the TFZ method. Dalguer and Day therefore
tested modified grid configurations for the SG and TFZ methods. With the fault
thickness reduced to 0.75h in the SG and 0.5h in the TFZ modeling, the rupture
velocities approach that in the TSN modeling. At the same time, the behaviors
of both fault zones depend on the ratio between h and the fault thickness which
has to be adjusted ad hoc. In any case, the theoretical efficiency of the 4th-order
fault-zone formulations is lost if sufficiently accurate results are to be obtained.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 493

6.3.6. Alternative Approaches

An interesting approach was presented by Ionescu and Campillo (1999) who

studied the 2-D problem of slip instability under slip-dependent friction. They
used a combination of two space–time FD grids. An interior (or body) grid is
used to solve the equation of motion away from the fault. The other grid is used to
implement frictional boundary conditions on the fault. The first-order hyperbolic
system (velocity-stress formulation) is reduced on the fault to a first-order system
in one space variable (orthogonal to the fault) for the shear velocity and shear
stress-tensor components. Using the integration on the characteristic lines and
frictional boundary conditions a nonlinear unstable ordinary differential equation
for slip is obtained. The equation depends on the computed values of the stresses
and particle velocities away from the fault. A small local time step is used to solve
the equation during the larger time step applied in the body grid.
Nielsen modified a standard application of the 4th-order staggered-grid FD
scheme to investigate rupture propagation (e.g., Nielsen et al., 2000; Nielsen and
Carlson, 2000). In order to satisfy the frictional boundary conditions with a slip
rate and tractions evaluated at the same time, he approximates a field variable at
a given time using a 2nd-order interpolation between their values at consecutive
time levels. Thus he reformulates the velocity-stress formulation of the equation
of motion and Hooke’s law with the same time-derivative operator applied to the
particle-velocity and stress-tensor components. In addition to this, Nielsen avoids
application of the spatial-derivative operator across the fault if the operator should
apply to a discontinuous field. He achieves this by reducing the spatial-derivative
operator down to the 2nd-order in evaluating shear stress-tensor components at a
grid point on the grid planes one grid spacing above and one grid spacing below
the fault, if the grid point one grid spacing away is slipping.

6.3.7. Modeling of Non-Planar Faults Using Partly-Staggered Grids

Cruz-Atienza and Virieux (2004) developed a new FD approach to model the

dynamic rupture propagation on non-planar faults. They used a 2-D velocity-
stress formulation on a partly staggered grid. Their formulation is based on a new
application of the fault boundary conditions. No split nodes are considered, that
is, there is no explicit discontinuity at any grid point. A finite source, rupturing
fault, is represented by a set of neighboring numerical cells placed alongside a
fault without sharing any stress grid point. Because this is the first FD approach
enabling both medium heterogeneity and non-planar fault geometry, we describe
its essential features in some detail.
Consider a Cartesian coordinate system and a corresponding spatial rectangu-
lar grid. A square h × h FD grid cell has particle-velocity vector positions at four
corners and stress-tensor
√ position
√ at the center. Then one n-point square numer-
ical cell is a set of n × n neighboring square FD cells. This means that one
494 MOCZO ET AL.

n-point square numerical cell has n stress-tensor grid positions (points). Consid-
ering a fixed spatial support S (square area), the scaling relation for a numerical
cell is S = n · h2n , where n is the number of FD square hn × hn grid cells mak-
ing one n-point numerical cell. The greater n the smoother the fault-geometry
discretization is.
Boundary conditions on the fault are applied locally in each numerical cell. The
application means the following. A planar fault tangential to the discretized fault
geometry is considered inside the cell. Let ϑ be an angle between the x-axis and
the fault plane. The Cartesian stress-tensor components at a stress grid point are
rotated into shear and normal stresses with respect to the fault plane. The same
boundary condition, following a considered friction law, is applied at each stress
grid point within the numerical cell. After imposing the boundary condition, the
stresses are rotated back to the original coordinate system connected with the
grid.
The fault plane always passes through the center of the numerical cell. Ve-
locity grid points in a numerical cell belong to either positive or negative fault
block. At a velocity grid point I , displacement uI parallel to the fault plane is
obtained from the projected particle velocity by simple time integration. Because
no split nodes are considered, displacement uI depends on a distance between
the point and the fault. For any time t, uI will be at its maximum, if the point
lies on a line passing through the center of the cell and perpendicular to the fault.
It will be zero, if the grid point lies on the fault. For a given ϑ, it may happen
that no velocity grid point lies on a line passing through the center of the cell
and perpendicular to the fault. A weight function HI (ϑ) for a grid point I is de-
fined as a ratio between the maximum displacement and displacement, that is,
HI (ϑ) = maxϑ {uI (t, ϑ)}/uI (t, ϑ). HI (ϑ) does not depend on time and has to
be determined for all grid points in a given type of a numerical cell. Then the
weighted displacements ũI (t) = uI (t, ϑ) · HI (ϑ) do not depend on the fault ori-
entation in the spatial grid. The relative displacement of the, for example, positive
fault block, D + is then determined as D + = ( I =1 ũ2I (t))1/2 /p, where p is the
p

number of the velocity grid points. Points within an angular vicinity γ around
the fault (two opposed sectors of a circle centered at the cell’s center) are not
considered because their displacements are too small. The angle is determined as
√
γ = arctan( n − 1).
Cruz-Atienza and Virieux (2004) found that the zigzag discrete shape of the
fault associated with a low number n of grid points leads to unwanted destructive
dynamic interaction of cells. At the same time, as the grid spacing is smaller, the
dynamic interaction is stronger. As a consequence, in order to increase numerical
resolution, it is necessary both to reduce grid spacing and enlarge n. To guar-
antee the low-frequency equivalence between different discretizations of a fault
√
length L, the scaling relationship L ≥ 30hn n has to be followed.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 495

6.4. Non-Reflecting Boundaries

6.4.1. Absorbing Boundary Conditions

For computational reasons it is almost always necessary to limit the computa-

tional domain to the part of the Earth in the vicinity of the seismic source and
receivers as well as the Earth structure contributing to form the seismic response.
This is achieved by applying so-called non-reflecting boundaries or absorbing
boundary conditions (ABC) around the tractable truncated computational domain
such that no energy is transmitted or reflected from the boundary of the compu-
tational domain back into its interior. Thus from the perspective of an observer
located inside the computational domain the signals appear to have been perfectly
absorbed by the boundary.
A significant proportion of the literature on time-domain FD modeling is de-
voted to the design of efficient ABC but most published approaches can be divided
into two groups. The first group, including the popular Clayton and Engquist
(1977); Higdon (1991); Lindman (1975) and Liao et al. (1984), attempt to ex-
trapolate the wavefield beyond the edge of the computational domain and then
use this extrapolated value in the spatial discretization operator used to update lo-
cations inside the computational domain. An interesting approach was suggested
by Peng and Toksöz (1994, 1995).
Collino (1993) and later other authors attempted to apply so-called high-order
ABCs which, in general, is an infinite sequence of ABCs with increasing accu-
racy. An advantage of the high-order ABCs is that they are implementable for an
arbitrary high order (Givoli, 2004).
The second group, including the approach of Kosloff and Kosloff (1986) and
Perfectly Match Layers (PML; Bérenger, 1994), gradually attenuate the amplitude
of the wavefield within a “sponge” layer within the boundary.
For completeness, we should also mention work on complementary bound-
ary conditions and operators (Smith, 1974; Schneider and Ramahi, 1998), which
does not fit within either of the two groups. Dirichlet and Neumann boundary
conditions are complementary as they have opposite reflection coefficients. By
adding the response from two simulations with different complementary boundary
conditions, the first-order boundary reflections will destructively interfere. Unfor-
tunately, higher-order interactions will not necessarily cancel and the methods
have therefore not found widespread application.
Each ABC approach has it own characteristics, advantages and disadvantages.
In the first group mentioned above, for example, the Clayton–Engquist and Hig-
don methods have the advantage of requiring relatively little memory but only
work well within a limited range of angles of incidence. Lindman’s approach
works well over a wide range of angles but only works when the material along
the boundary is homogeneous. Liao’s ABC works well over most angles and for
some types of heterogeneity adjacent to the boundary but requires double preci-
496 MOCZO ET AL.

sion implementation to maintain stability. It also requires more computer memory

than the other methods listed above.
In the second group of ABCs, Kosloff’s method (Kosloff and Kosloff, 1986)
amounts to applying an exponential damping function of the wavefield within
a “sponge” region surrounding the computational domain. Kosloff’s method is
effective at normal incidence and is very easy to implement but requires a rela-
tively large amount of memory and can generate significant reflection artifacts
away from normal incidence. A related method known as perfectly matched
layers (PML) was first introduced for electromagnetic wave propagation by
Bérenger (1994), and later generalized for elastic wave propagation by (Chew
and Liu (1996); Collino and Tsogka (1998, 2001); Komatitsch and Tromp (2003);
Marcinkovich and Olsen (2003); Festa and Nielsen (2003)), and for viscoelastic
anisotropic media by Chen et al. (2000). PMLs have now become established as
the most efficient absorbing boundary condition available, offering simplicity of
implementation, stability, good absorption over a wide range of incident angles
and frequencies and less memory requirements than Kosloff’s ABC. In this sec-
tion we will therefore focus on PMLs, review the theory and discuss issues related
to their implementation.

6.4.2. The Perfectly Matched Layer (PML)

(γ )
Following Chen et al. (2000), we define the split particle velocities vi such
that the equation of motion in a general heterogeneous anisotropic elastic medium
is decomposed into three equations
(γ )
ρ v̇i = σiγ ,γ (174)

with
3
(γ )
vi = vi (175)
γ =1

where the Einstein summation convention is not assumed for the Greek letter γ .
Similarly, the constitutive equation for the split stiffness tensors may be written
as
(γ )
σ̇ij = cij γ q vq,γ (176)

with
3
(γ )
σij = σij . (177)
γ =1
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 497

Following Chew and Liu (1996), we introduce a modified spatial differencing

operator
1
∂˜γ = ∂γ , (178)
eγ
where the denominator eγ has been interpreted as a coordinate stretching variable.
Outside the PML layers, vector e = (1, 1, 1). Inside the PML layers, a complex eα
is needed in order to attenuate the split (v (γ ) , σ (γ ) ) in the γ direction. Assuming
a time dependence of e−iωt we choose
iγ
e γ = aγ + (179)
ω
where aγ and γ are the so-called PML medium profiles. In the time domain,
this leads to the following equations for the partial variables:
(γ )
ρ(aγ ∂t + γ )vi = ∂γ σiγ (180)
and
(γ )
(aγ ∂t + γ )σij = cij γ q ∂γ vq . (181)
The selection of the PML material profile aγ and γ is key to the effective im-
plementation of the PML ABC. Theoretically, these parameters can be arbitrary
in the continuum space and the PML layers would perfectly absorb any incident
wave. This is, however, not true in the discretized model. The discretization error
is proportional to the product of the spatial discretization size and the contrast of
these parameters (Chew and Weedon, 1994; Chew and Liu, 1996). A PML sponge
of width N = 10 with the following profile at each layer i = 1, 2, . . . , N has been
found to result in excellent absorption for many seismic applications for a wide
range of frequencies, Courant numbers and wave types (Chen et al., 2000):
  2.8
1 N/2 − i + 1
aγ (i) = 1 + amax 1 + sin π (182)
2 N
and
  2.8
1 N/2 − i + 1
γ (i) = max 1 + sin π (183)
2 N
with amax = 0.075 and max = 0.91/t.
We note that the PML ABC differs fundamentally from the traditional material
sponge ABCs such as the method by Kosloff and Kosloff (1986) as a result of the
split wavefield formulation. For each set of split wavefields, only the propagation
vector component normal to the boundary is complex. The propagation vector
components parallel to the boundary remain the same across the boundary layers.
As a result, a perfect match can be obtained while the wavefields are attenuated
by the perfectly matched layers, at least in theory in the non-discretized model.
498 MOCZO ET AL.

Finally, the extension of PMLs to viscoelastic media is straightforward (Chen

et al., 2000). The memory variables governing the viscoelastic behavior of the
medium due to propagating waves are updated from the change in stress and
record the history of the stress. The memory variables in the PML region are
therefore split in exactly the same way as the stresses. The equation for updating
the split stress fields becomes, see Eq. (111),
n
= cij γ q ∂˜γ vq +
(γ ) ij (γ )
(aγ ∂t + γ )σij rl (184)
l=1
and the equation for updating the split memory variable fields is, see Eq. (113),
 
1 ij (γ ) 1 τεl
[cij γ q ]R ∂˜γ vq , l = 1, . . . , n,
ij (γ )
ṙl (t) + rl (t) = 1−
τσ l τσ l τσ l
(185)
where [cij γ q ]R means relaxed moduli. The equations for updating the particle
velocities in split form remain the same, Eq. (180).

6.4.3. PML Reflection Coefficients

In Fig. 15 we show reflection and conversion coefficients for incident P- and

S-waves. The curves were obtained from 2-D simulations in a homogeneous
isotropic elastic medium with a P-velocity of 2500 m/s, an S-velocity of 1500 m/s
and a density of 1100 kg/m3 where a point source was positioned close to the
PML boundary (Chen et al., 2000). For incident P-waves the source was explo-
sive, for S-waves the source was a force oriented normal to the PML boundary. In
each case the source was a 50 Hz Ricker wavelet. The reflections were corrected
for propagation distance and radiation pattern to obtain reflection and conversion
coefficients.
Simulations were carried out for PML boundary thicknesses of 10, 20, 30,
40 and 50 cells, and compared to a 40 cell thick Kosloff ABC (shown in the
dashed lines). First, the results show that the choice of a 10 grid-point wide PML
boundary is roughly as efficient as the 40 grid-point wide boundary (justifying
the values we used for cost comparisons in the previous section). However, the
poor performance at low grazing angles may not be acceptable for many appli-
cations. Increasing the thickness of the PML boundaries drastically improves the
effectiveness so that we are close to single-precision machine accuracy for the 50
grid-point wide boundaries. The PS-plot (top right of Fig. 15) reveals very little
P-to-S conversion. Machine precision is achieved with a 20 cell thick PML.
The SS-reflection coefficient (bottom left in Fig. 15) shows a behavior similar
to the PP-results. The SP-plot (bottom right in Fig. 15) does not show the same
remarkable performance as the PS-converted wave reflection although the per-
formance is still satisfactory. It is possible that the SS-and SP-results could be
degraded by P-wave energy generated by the source used.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 499

F IG . 15. Amplitude reflection coefficient as a function of incidence angle for PML boundaries of
width 10, 20, 30, 40 and 50 grid-points (solid) as well as for a 40 grid-points wide Kosloff sponge
(dashed). Isotropic model. Top left: Normalized PP reflection. Top right: Normalized PS reflection.
Bottom left: Normalized SS reflection. Bottom right: Normalized SP reflection. Reproduced from
Chen et al. (2000).

Figure 16 shows the corresponding results for an anisotropic model with the
same density and vertical velocities as the isotropic model but with a vertical to
horizontal P-velocity ratio of 1.12 and an anellipticity of 0.16 as defined by Carrio
et al. (1992). The results do not differ significantly from the ones shown in Fig. 15
for the isotropic case.

6.4.4. PML Implementation Issues and Computational Cost

In general Eq. (180) is a system of nine equations: three partial, or split, veloci-
ties for each velocity component. The split elastic constitutive relation (181) is in
general a system of eighteen equations: three partial stresses for each of six inde-
pendent elements in the stress tensor. For the viscoelastic case, Eq. (181) translates
into Eqs. (184) and (185) which, in general, are systems of eighteen equations. If
all of these variables were required on the complete 3-D FD grid, the resulting
500 MOCZO ET AL.

F IG . 16. Same as Fig. 15 but for an anisotropic model. Reproduced from Chen et al. (2000).

memory requirements would make this approach unattractive. However, because

updating split field variables require only derivative of the full velocities and
stresses, the domain in which the split field variables are defined (and stored)
is restricted to a relatively narrow region around the boundaries.
The PML formulation presented here requires almost the same number of cal-
culations per grid point as the interior of the computational domain since most of
the calculation required for updating each split field is nothing but a constituent of
the conventional time-domain FD calculation which must be evaluated in any case
(e.g., the spatial derivatives). However, because of the splitting of variables, there
is a significant price to pay in storage per grid point (approximately a factor three
depending on dimension, model properties, and other factors). Nevertheless, with
a 10 grid-points wide PML sponge the PML conditions are still computationally
significantly more efficient both in terms of CPU and memory compared to other
ABC (Chen et al., 2000).
Recently Wang and Tang (2003) proposed an alternative PML implementation
which avoids splitting the wavefields. The storage requirements per grid point
are therefore equivalent to those in the interior of the grid. However, the CPU
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 501

requirements increase to solve for a new set of differential equations with similar
characteristics to the viscoelastic equations presented above.
We note that it is straightforward to interface PML boundary conditions with a
free-surface condition.
As we have seen, despite what the name might suggest, PMLs are still far from
perfect. Because of inaccuracy related to discrete numerical implementations of
PMLs, a finite sponge width is required to avoid artificial reflections. The cost in
terms of memory is therefore not insignificant, particularly for hybrid-modeling
scenarios or grids that are narrow in one direction (e.g., cross-line), when the
size of the absorbing boundary regions actually may be comparable to the size
of the computational domain of interest. Moreover, if energy is incident under
very low grazing angles, the PMLs will cause significant reflections. The search
for a truly perfect absorbing boundary condition is therefore likely to continue as
demands for absorbing boundary conditions that are effective in small (deformed)
time-domain FD grids and extreme low angles of incidence is increasing with the
interest in 3-D applications and hybrid modeling scenarios on the rise. One could
then imagine solving the wave equation in 3-D on grids that were only one or a
few grid-points wide in the cross-line direction, at a cost that would be comparable
to that of a 2-D simulation. This would enable synthesizing seismic data with true
3-D amplitudes such that multiples, surface waves, guided waves and body waves
all had the correct relative phases and amplitudes.

7. M EMORY O PTIMIZATION AND PARALLELIZATION

7.1. Memory Optimization

7.1.1. Memory Requirements in the Staggered-Grid Schemes

Many examples of FD modeling of earthquake motion in large sedimentary

basins or seismic wave propagation in complex models in seismic exploration
require very large computer memory and CPU. In order to solve the most chal-
lenging problems of seismology, it is necessary to reduce both memory and CPU
requirements. Thus, memory optimization, sophisticated programming and par-
allelization are practically necessary. Here we briefly review several possibilities
to reduce memory requirements and the number of operations. We also briefly
mention parallelization.
Consider staggered-grid FD schemes. Assume that the material parameters of
the medium can vary between any two grid positions. Such heterogeneity can
be referred to as the point-to-point heterogeneity. Though it is sometimes an im-
proper simplification, a homogeneous medium inside a grid cell is assumed by
some modelers. Such material parameterization can be referred to as homoge-
neous material cells.
502 MOCZO ET AL.

Denote the displacement and particle-velocity components at time level m as

U m , V m , W m , and U̇ m , V̇ m , Ẇ m , respectively. In the case of the point-to-point
heterogeneity, each grid position of the displacement or particle-velocity compo-
nents is assigned its value of density, that is, there are ρ U , ρ V and ρ W assigned
to a grid cell. Similarly, the elastic (unrelaxed) moduli μxy , μyz and μzx are as-
signed to grid positions of the shear stress-tensor components, while κ and μ are
assigned to the joint grid position of the normal stress-tensor components. Cor-
μ μxy μyz
responding to the elastic moduli, the anelastic coefficients Ylκ , Yl , Yl , Yl
μzx
and Yl , l = 1, . . . , n, are assigned to the grid positions of the stress-tensor
components.
An assumption of homogeneous material cells reduces the material parameters
μ
in a grid cell to ρ, κ, μ, Ylκ and Yl , l = 1, . . . , n. Consider displacement-stress
(DS), displacement-velocity-stress (DVS) and velocity-stress (VS) staggered-grid
FD schemes (for the equations and simplified schemes see, e.g., Moczo et al.,
2001) which are 2nd-order accurate in time. The three schemes require the fol-
lowing variables to be stored in primary memory for a grid cell:

DS: U m , V m , W m , U m−1 , V m−1 , W m−1 , (186)

DVS: U m , V m , W m , U̇ m−1/2 , V̇ m−1/2 , Ẇ m−1/2 , (187)

VS: U̇ m−1/2 , V̇ m−1/2 , Ẇ m−1/2 , Txx

m−1 m−1
, Tyy , Tzzm−1 , Txy
m−1
,
m−1
Tyz m−1
, Tzx . (188)

As discussed in the section on incorporation of the realistic attenuation, the coarse

spatial sampling of the anelastic functions (memory variables) can be considered.
In the spatial distribution used by Kristek and Moczo (2003), n = 4 and each
grid cell accommodates all six anelastic functions for just one of the relaxation
yy xy yz
frequencies, that is ξlxx , ξl , ξlzz , ξl , ξl and ξlzx , l ∈ {1, 2, 3, 4}.
Consider a computational region with dimensions XL, Y L and ZL. Let βmin
and αmax be the minimum S-wave and maximum P-wave velocity, respectively.
Let fac be the frequency up to which the FD computation should be sufficiently
accurate. The maximum spatial grid spacing is the h = (βmin /fac ) · s. Here, s is
the spatial sampling ratio s = h/λmin which has to be chosen based on the grid
dispersion in the considered FD scheme and the wave propagation distance. The
time step is t = q · p · (h/α)min , 0 < p ≤ 1, where p is the stability ratio and q
depends on the FD scheme. For√example, for the 2nd-order in time, 4th-order in
space staggered-grid, q = 6/(7 3); see, e.g. Moczo et al. (2000).
Assuming a uniform grid with the grid spacing h in the three Cartesian direc-
tions, the numbers of the grid cells in the three directions are MX = (XL/ h) + 1,
MY = (Y L/ h) + 1, MZ = (ZL/ h) + 1.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 503

In summary, the number of material parameters and field variables in the

staggered-grid FD schemes are:

DS, DVS: MX · MY · MZ · (3 + 5 + 5 · 4) + (6 + 6)
= MX · MY · MZ · 40 (189)

and

VS: MX · MY · MZ · (3 + 5 + 5 · 4) + (9 + 6)
= MX · MY · MZ · 43. (190)

Here, (3 + 5 + 5 · 4) stands for 3 densities, 5 unrelaxed moduli and 5 · 4 anelastic

coefficients, while (6 + 6) for 6 displacement/particle-velocity components and 6
anelastic functions in one grid cell. In the VS scheme, (9 + 6) stand 3 for particle-
velocity components, 6 stress-tensor components and 6 anelastic functions. If the
anelastic coefficients are spatially distributed in the same manner as the anelastic
functions, 5 · 4 is replaced by 5 in formulae (189) and (190). In the case of the
homogeneous grid cells, 3 + 5 + 5 · 4 is replaced by 1 + 2 + 2 · 4.
Clearly, depending on the total number of the grid cells, MX · MY · MZ, the
memory requirements can be very large.
The total number of the grid cells and thus the memory requirements can be
reduced by using a higher-order approximation in space (for example, Dablain,
1986), grid with varying size of the grid spacing or discontinuous grid. The
4th-order staggered-grid schemes were introduced by Bayliss et al. (1986) and
Levander (1988). In recent FD modeling, the 4th-order accuracy in space is al-
most necessary. Yomogida and Etgen (1993) and Rodrigues (1993) used the 8th-
order displacement-stress schemes. The rectangular grid with a varying size of
the grid spacings was first used by Boore (1970) in the 1-D problem. Mikumo
and Miyatake (1987) applied the varying size of the grid spacing in the 3-D
case in a homogeneous medium. Moczo (1989) applied the grid with the vary-
ing size of the grid spacing to the 2-D SH problem in the laterally heterogeneous
medium, Pitarka (1999) presented the 3-D velocity-stress scheme. Jastram and
Behle (1992), Jastram and Tessmer (1994), Falk et al. (1996), Moczo et al. (1996),
Kristek et al. (1999), Aoi and Fujiwara (1999), Hayashi et al. (2001), and Wang
et al. (2001) introduced discontinuous grids.
Given some spatial grid, the core memory can be reduced using core memory
optimization. While simple variants were used before, Graves (1996) described
an optimized procedure. Moczo et al. (1999) presented a combined memory op-
timization which naturally combines core and disk memory optimizations. One
other possibility to reduce the core memory is the use of the material cell types
(e.g., Moczo et al., 2001). Here we briefly characterize the corresponding reduc-
tions.
504 MOCZO ET AL.

7.1.2. Material Cell Types

In most models it is possible to efficiently describe heterogeneity of the medium

by a spatial distribution of integer numbers (a look-up table). In this fashion each
grid cell is assigned an integer number which represents a type of a material grid
cell, that is, a set of material parameters characterizing the medium inside the grid
cell. Such a description is efficient if there are sub-volumes of the computational
model that can be covered by many grid cells of the same material type. Let K be
the total number of different material types necessary to characterize heterogene-
ity of the whole model. Then the number of the material parameters and variables
needed by the schemes is reduced to
DS, DVS: MX · MY · MZ · (6 + 6 + 1) + (3 + 5 + 5 · 4) · K (191)
and
VS: MX · MY · MZ · (9 + 6 + 1) + (3 + 5 + 5 · 4) · K (192)
assuming that K < MX · MY · MZ which is the case even in relatively complex
models. The number one added to the numbers of the field variables represents the
distribution of the integers corresponding to the material cell types. In the case of
the coarsely distributed anelastic coefficients or homogeneous material grid cells,
the same reductions as in formulae (189) and (190) apply.

7.1.3. Core Memory Optimization

In core memory optimization method (Graves, 1996) such that only limited
number of grid planes are kept in core memory all possible time updates for these
planes are carried out. The subset of planes, say, NP planes, repeatedly moves
throughout the whole model space, and the displacement or particle-velocity com-
ponents, stress-tensor components and anelastic functions are successively (plane
by plane) and periodically overwritten in disk. The key aspect in that all possi-
ble time updates are performed for the planes in the core memory means that the
procedure is split into three parts—roll-in, cascade, and roll-out. Consider, for
example, that the subset is made of horizontal grid planes. Then in the formu-
lae (191) and (192) MZ is replaced by NP. The smaller the NP, the larger the
reduction is.

7.1.4. Combined Memory Optimization

It is obvious that whereas the requirement for core memory can be significantly
reduced, the needed amount of disk space can become very large. Although avail-
able disk space becomes larger and access to disk memory becomes faster, still,
in principle, it is possible to reduce also the disk memory needed in the core-
memory-optimization procedure. In the disk memory optimization (Moczo et al.,
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 505

1999), the wavelet transform is applied first to 2-D array of each displacement or
particle-velocity component and each anelastic function. The transform decreases
the information entropy. Then the sets of the wavelet coefficients are compressed
by a standard compression procedure. Thus, for example, instead of the MX · MY
displacement-component values in one grid plane only a relatively short stream
of zeros and ones has to be written in disk. The total number of the material para-
meters and field variables stored in disk is then
DS, DVS: MX · MY · MZ · (6 + 6)/CR + 1 + (3 + 5 + 5 · 4) · K
(193)
and
VS: MX · MY · MZ · (9 + 6)/CR + 1 + (3 + 5 + 5 · 4) · K (194)
where CR is the compression ratio. A reasonable value is about 10. Moczo et al.
(1999) found that the increase of the CPU time due to one passage of the subset
of planes with compression was always smaller than 0.75% of the time for one
passage without compression. This is because the increase due to compression
itself is partly compensated by the smaller number of the I/O operations.

7.1.5. Spatial Discontinuous Grid

In many models the S- and P-wave velocities are lower near the Earth’s sur-
face. In such a case it is advantageous to cover the lower part with a coarser
spatial grid. Due to the structure of the staggered grid, the most natural combined
(discontinuous) grid is the one whose upper part is the h × h × h grid and the
lower part the 3h × 3h × 3h grid (Aoi and Fujiwara, 1999; Kristek et al., 1999;
Robertsson and Holliger, 1997). Let MZH be the number of the grid cells in the
z-direction in the upper grid. Let MZ3H be the number of the grid cells in the
z-direction in the lower grid. Assuming only the material cell types, the number
of material parameters and variables needed by the schemes is
    
MX − 1 MY − 1
DS, DVS: MX · MY · MZH + +1 · +1
3 3

× MZ3H · 6 + 6 + 1 + (3 + 5 + 5 · 4) · K (195)

and
    
MX − 1 MY − 1
VS: MX · MY · MZH + +1 · +1
3 3

× MZ3H · 9 + 6 + 1 + (3 + 5 + 5 · 4) · K . (196)

Formulae (195) and (196) are easy to modify if the core or combined memory
optimization are applied.
506 MOCZO ET AL.

A combination of the material cell types, combined memory optimization and

discontinuous grid typically can reduce the memory requirements by more than
one order of magnitude. An example is given by Moczo et al. (2001) for the
modeling of the 1995 Kobe, Japan, earthquake.

7.1.6. Spatially Varying Time Step

Falk et al. (1998) and Tessmer (2000) introduced a combined grid in which a
smaller time step is applied to the upper part of the grid while a larger time step
is applied to the lower part. Their techniques considerably reduce CPU cost.

7.1.7. Discontinuous Space–Time Grids

Kang and Baag (2004a, 2004b) developed efficient techniques for 2-D and 3-D
4th-order staggered-grid modeling. They combined discontinuous grid in space
with a discontinuous time step. While time integration in a finer grid with grid
spacing h is performed with time step t, time integration in a coarser grid with
grid spacing 3h is performed with time step 3t. The finer grid covers 2-D or 3-D
rectangular subregion which may have a planar free surface. This enables efficient
modeling of localized surface sedimentary structures. Proportionality of the time
step to the grid spacing is due to the fact that the two spatial grids have to overlap
in the medium with a higher speed. The technique considerably reduces both the
number of grid points and the number of operations.

7.2. Parallelization
Because the FD operators are local, the FD algorithms are suitable for paral-
lelization. Over the past decade, several distinct approaches have been applied to
parallelize FD codes.

7.2.1. Message Passing Libraries

Message-passing libraries, such as the Message Passing Interface (MPI, Gropp

et al., 1994), represent an approach suitable for shared or distributed memory
architectures. The MPI typically requires the involvement of the sender-receiver
communication: the source process makes a call to send data and the destination
process makes a call to receive it. While the scaling of the MPI parallelized codes
can be very good, the preparation of the code can require considerable time and
effort. The resulting MPI parallelized code often differs substantially from the
original source code.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 507

7.2.2. Parallelizing Compilers

Some compilers are capable of producing a parallel code that is portable to

shared and distributed memory machines. The compiler analyzes the inside-code
dependences and can produce a parallelized code with or without additional user’s
directives and/or language extensions. The High-performance Fortran language
(HPF, Koelbel et al., 1994) is one of the best known examples. In some cases the
resulting parallel code is quite efficient (e.g., Caserta et al., 2002) but deficiencies
of the approach are also known. In an automatic parallelization regime, compilers
often make conservative assumptions on data dependence, which usually yields
lower efficiency.

7.2.3. Interactive Parallelization Tools

There are parallelization tools which provide a user with the possibility of com-
bining the automatic parallelization analysis with user’s knowledge of the code.
Whereas such interactive tools usually lead to quick parallelization, the efficiency
can be limited due to non-optimized nesting, multiple decompositions and For-
tran 90 constructs. Often additional manual modifications of the source code are
necessary to obtain good performance.

7.2.4. High-Level Library-Based Tools

Library-based tools exist which are designed to help the user with the applica-
tion of the lower-level libraries, such as MPI. Particularly, a user does not need to
handle all details of the MPI parallelization. Some tools enable the use of addi-
tional optimizations specific for the machine architecture. The preparation of the
parallel code may still be very time consuming and invasive.

7.2.5. Directive-Based Parallelization

Some computer manufacturers, for example Cray and SGI, introduced the pos-
sibility to manually supplement parallelization directives in the source code. In the
beginning, the directives were mainly used to support loop-level shared memory
parallelization. Recently SGIs OpenMP has become the best known and stan-
dard tool for relatively efficient and quick parallelization. One problem is that the
OpenMP is most efficient only with shared-memory architectures.

ACKNOWLEDGEMENTS

The authors would like to thank Y.H. Chen, L.A. Dalguer, S.M. Day, R.T. Coates, M. Gális, R.J.
Geller, I. Ionescu, J. Kristek, V. Maupin, S. Nielsen, and J. Virieux for useful discussions and com-
508 MOCZO ET AL.

ments on selected topics. Thanks also go to D. Lafleur and W. Kimman for misprint corrections. The
work was supported in part also by the Geophysical Institute, Slovak Academy of Sciences.

R EFERENCES

Aki, K., Richards, P.G. (1980). Quantitative Seismology. Theory and Methods, vols. I and II. Freeman,
San Francisco.
Alekseev, A.S., Mikhailenko, B.G. (1980). The solution of dynamic problems of elastic wave prop-
agation in inhomogeneous media by a combination of partial separation of variables and finite-
difference methods. J. Geophys. 48, 161–172.
Alterman, Z., Karal, F.C. (1968). Propagation of elastic waves in layered media by finite-difference
methods. Bull. Seismol. Soc. Am. 58, 367–398.
Anderson, D.A., Tannehill, J.C., Pletcher, R.H. (1984). Computational Fluid Mechanics and Heat
Transfer. Hemisphere Publishing Corporation.
Andrews, D.J. (1973). A numerical study of tectonic stress release by underground explosions. Bull.
Seismol. Soc. Am. 63, 1375–1391.
Andrews, D.J. (1976). Rupture propagation with finite stress in antiplane strain. J. Geophys. Res. 81,
3575–3582.
Andrews, D.J. (1999). Test of two methods for faulting in finite-difference calculations. Bull. Seismol.
Soc. Am. 89, 931–937.
Aochi, H., Fukuyama, E. (2002). Three-dimensional non-planar simulation of the 1992 Landers earth-
quake. J. Geophys. Res. 107, doi:10.1029/2000JB000061.
Aoi, S., Fujiwara, H. (1999). 3-D finite-difference method using discontinuous grids. Bull. Seismol.
Soc. Am. 89, 918–930.
Backus, G.E. (1962). Long-wave elastic anisotropy produced by horizontal layering. J. Geophys.
Res. 67, 4427–4440.
Bayliss, A., Jordan, K.E., LeMesurier, B.J., Turkel, E. (1986). A fourth-order accurate finite-difference
scheme for the computation of elastic waves. Bull. Seismol. Soc. Am. 76, 1115–1132.
Bérenger, J.-P. (1994). A perfectly matched layer for the absorption of electromagnetic waves. J.
Comput. Phys. 114, 185–200.
Blanch, J.O., Robertsson, J.O.A. (1997). A modified Lax–Wendroff correction for wave propagation
in media described by Zener elements. Geophys. J. Int. 131, 381–386.
Blanch, J.O., Robertsson, J.O.A., Symes, W.W. (1995). Modeling of a constant Q: Methodology and
algorithm for an efficient and optimally inexpensive viscoelastic technique. Geophysics 60, 176–
184.
Boore, D.M. (1970). Love waves in nonuniform waveguides: Finite difference calculations. J. Geo-
phys. Res. 75, 1512–1527.
Boore, D. (1972). Finite-difference methods for seismic wave propagation in heterogeneous materials.
In: Bolt, B.A. (Ed.), In: Methods in Computational Physics, vol. 11. Academic Press, New York.
Bouchon, M. (1981). A simple method to calculate Green’s functions for elastic layered media. Bull.
Seismol. Soc. Am. 71, 959–971.
Carcione, J.M. (1994). The wave equation in generalized coordinates. Geophysics 59, 1911–1919.
Carcione, J.M. (2001). Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic and
Porous Media. Pergamon.
Carcione, J.M., Cavallini, F. (1994). A rheological model for anelastic anisotropic media with appli-
cations to seismic wave propagation. Geophys. J. Int. 119, 338–348.
Carcione, J.M., Herman, G.C., ten Kroode, A.P.E. (2002). Seismic modeling. Geophysics 67, 1304–
1325.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 509

Carcione, J.M., Kosloff, D., Kosloff, R. (1988a). Wave propagation simulation in a linear viscoacoustic
medium. Geophys. J. 93, 393–407.
Carcione, J.M., Kosloff, D., Kosloff, R. (1988b). Wave propagation simulation in a linear viscoelastic
medium. Geophys. J. 95, 597–611.
Carcione, J.M., Wang, J.P. (1993). A Chebyshev collocation method for the elastodynamic equation
in generalized coordinates. Comput. Fluid Dyn. J. 2, 269–290.
Carrio, P., Costa, J., Ferrer-Pinheiro, J.E., Schoenberg, M. (1992). Cross-borehole tomography in
anisotropic media. Geophysics 57, 1194–1198.
Caserta, A., Ruggiero, V., Lanucara, P. (2002). Numerical modelling of dynamical interaction between
seismic radiation and near-surface geological structures: a parallel approach. Comput. Geosci. 28,
1069–1077.
Chaljub, E., Komatitsch, D., Vilotte, J.P., Capdeville, Y., Valette, B., Festa, G. (2006). Spectral Ele-
ment Analysis in Seismology. In: Wu, R.-S., Maupin, V. (Eds.), Advances in Wave Propagation in
Heterogeneous Earth. In: Advances in Geophysics, vol. 48. Dmowska, R. (Ed.). Elsevier, p. 365
(this book).
Chen, Y.H., Coates, R.T., Robertsson, J.O.A. (2000). Extension of PML ABC to Elastic Wave Prob-
lems in General Anisotropic and Viscoelastic Media: Schlumberger OFSR research note.
Chew, W.C., Liu, Q.H. (1996). Perfectly matched layers for elastodynamics: a new absorbing bound-
ary condition. J. Comput. Acoustics 4, 341–359.
Chew, W.C., Weedon, W.H. (1994). A 3-D perfectly matched medium from modified Maxwell’s equa-
tions with stretched coordinates. Microwave Opt. Tech. Lett. 7, 599–604.
Clayton, R., Engquist, B. (1977). Absorbing boundary conditions for acoustic and elastic wave equa-
tions. Bull. Seismol. Soc. Am. 67, 1529–1540.
Cohen, G.C. (2002). Higher-Order Numerical Methods for Transient Wave Equations. Springer.
Collino, F. (1993). High order absorbing boundary conditions for wave propagation models. In: Klein-
man, R., et al. (Eds.), Proceedings of the Second International Conference on Mathematical and
Numerical Aspects of Wave Propagation. SIAM, Delaware, pp. 161–171.
Collino, F., Tsogka, C. (1998). Application of the PML absorbing layer model to the linear elastody-
namic problem in anisotropic heterogeneous media. INRIA, Rapport de Recherche, No. 3471.
Collino, F., Tsogka, C. (2001). Applications of the PML absorbing layer model to the linear elastody-
namic problem in anisotropic heterogeneous media. Geophysics 66, 294–305.
Coutant, O. (1989). Program of Numerical Simulation AXITRA. Research Reports LGIT. Université
Joseph Fourier, Grenoble (in French).
Coutant, O., Virieux, J., Zollo, A. (1995). Numerical source implementation in a 2-D finite difference
scheme for wave propagation. Bull. Seismol. Soc. Am. 85, 1507–1512.
Crampin, S., Chastin, S. (2003). A review of shear wave splitting in the crack-critical crust. Geophys.
J. Int. 155, 221–240.
Crase, E. (1990). High-order (space and time) finite-difference modeling of the elastic wave equa-
tion. In: 60th Annual International Meeting, Society of Exploration Geophysicits, pp. 987–991.
Expanded Abstracts.
Crase, E., Wideman, Ch., Noble, M., Tarantola, A. (1992). Nonlinear elastic waveform inversion of
land seismic reflection data. J. Geophys. Res. 97, 4685–4703.
Cruz-Atienza, V.M., Virieux, J. (2004). Dynamic rupture simulation of non-planar faults with a finite-
difference approach. Geophys. J. Int. 158, 939–954.
Dablain, M.A. (1986). The application of high-order differencing to the scalar wave equation. Geo-
physics 51, 54–66.
Dalguer, L.A., Day, S.M. (2004). Split nodes and fault zone models for dynamic rupture simulation.
EoS Trans. Am. Geophys. Union 85 (47). Fall Meet. Suppl., Abstract S41A-0944.
Dalguer, L.A., Day, S.M. (2006). Comparison of fault representation methods in finite difference sim-
ulations of dynamic rupture. Bull. Seismol. Soc. Am., in press.
Day, S.M. (1977). Finite element analysis of seismic scattering problems. Ph.D. Dissertation, Univer-
sity of California, San Diego.
510 MOCZO ET AL.

Day, S.M. (1982). Three-dimensional simulation of spontaneous rupture: the effect of nonuniform
prestress. Bull. Seismol. Soc. Am. 72, 1881–1902.
Day, S.M. (1998). Efficient simulation of constant Q using coarse-grained memory variables. Bull.
Seismol. Soc. Am. 88, 1051–1062.
Day, S.M. (2005). Personal communication.
Day, S.M., Bradley, C.R. (2001). Memory-efficient simulation of anelastic wave propagation. Bull.
Seismol. Soc. Am. 91, 520–531.
Day, S.M., Dalguer, L.A., Lapusta, N., Liu, Y. (2005). Comparison of finite difference and bound-
ary integral solutions to three-dimensional spontaneous rupture. J. Geophys. Res. B 110, 12307,
doi:10.1029/2005JB003813.
Day, S.M., Minster, J.B. (1984). Numerical simulation of wavefields using a Padé approximant
method. Geophys. J. R. Astron. Soc. 78, 105–118.
Durran, D.R. (1999). Numerical Methods for Wave Equations in Geophysical Fluid Dynamics.
Springer.
Eisner, L., Clayton, R. (2002). Equivalent medium parameters for numerical modeling in media with
near-surface low velocities. Bull. Seismol. Soc. Am. 92, 711–722.
Emerman, S.H., Schmidt, W., Stephen, R.A. (1982). An implicit finite-difference formulation of the
elastic wave equation. Geophysics 47, 1521–1526.
Emmerich, H. (1989). 2-D wave propagation by a hybrid method. Geophys. J. Int. 99, 307–319.
Emmerich, H. (1992). PSV-wave propagation in a medium with local heterogeneities: A hybrid for-
mulation and its application. Geophys. J. Int. 109, 54–64.
Emmerich, H., Korn, M. (1987). Incorporation of attenuation into time-domain computations of seis-
mic wave fields. Geophysics 52, 1252–1264.
Falk, J., Tessmer, E., Gajewski, D. (1996). Tube wave modelling by the finite-differences method with
varying grid spacing. Pure Appl. Geophys. 148, 77–93.
Falk, J., Tessmer, E., Gajewski, D. (1998). Efficient finite-difference modelling of seismic waves using
locally adjustable time steps. Geophys. Prosp. 46, 603–616.
Fäh, D. (1992). A hybrid technique for the estimation of strong ground motion in sedimentary basins.
Diss. ETH Nr. 9767, Swiss Federal Institute of Technology, Zurich.
Festa, G., Nielsen, S. (2003). PML absorbing boundaries. Bull. Seismol. Soc. Am. 93, 891–903.
Fink, M. (1997). Time reversed acoustics. Phys. Today 50, 34–40.
Fornberg, B. (1988). Generation of finite difference formulas on arbitrary spaced grids. Math. Com-
put. 51, 699–706.
Forsythe, G.E., Wasow, W.R. (1960). Finite Difference Methods for Partial Differential Equations.
Wiley and Sons, New York.
Frankel, A. (1993). Three-dimensional simulations of ground motions in the San Bernardino Valley,
California, for hypothetical earthquakes on the San Andreas fault. Bull. Seismol. Soc. Am. 83,
1020–1041.
Frankel, A., Leith, W. (1992). Evaluation of topographic effects on P- and S-waves of explosions at
the Northern Novaya Zemlya test site using 3-D numerical simulations. Geophys. Res. Lett. 19,
1887–1890.
Geller, R.J., Ohminato, T. (1994). Computation of synthetic seismograms and their partial derivatives
for heterogeneous media with arbitrary natural boundary conditions using the Direct Solution
Method. Geophys. J. Int. 116, 421–446.
Geller, R.J., Takeuchi, N. (1995). A new method for computing highly accurate DSM synthetic seis-
mograms. Geophys. J. Int. 123, 449–470.
Geller, R.J., Takeuchi, N. (1998). Optimally accurate second-order time-domain finite difference
scheme for the elastic equation of motion: One-dimensional case. Geophys. J. Int. 135, 48–62.
Givoli, D. (2004). High-order local non-reflecting boundary conditions: A review. Wave Motion 39,
319–326.
Graves, R.W. (1996). Simulating seismic wave propagation in 3-D elastic media using staggered-grid
finite differences. Bull. Seismol. Soc. Am. 86, 1091–1106.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 511

Graves, R.W., Day, S.M. (2003). Stability and accuracy analysis of coarse-grain viscoelastic simula-
tions. Bull. Seismol. Soc. Am. 93, 283–300.
Gropp, W., Lusk, E., Skjellum, A. (1994). Using MPI, Portable Parallel Programming with the Mes-
sage Passing Interface. MIT Press.
Hayashi, K., Burns, D.R., Toksöz, M.N. (2001). Discontinuous-grid finite-difference seismic modeling
including surface topography. Bull. Seismol. Soc. Am. 91, 1750–1764.
Helbig, K. (1984). Anisotropy and dispersion in periodically layered media. Geophysics 49, 364–373.
Hestholm, S. (1999). Three-dimensional finite difference viscoelastic wave modelling including sur-
face topography. Geophys. J. Int. 139, 852–878.
Hestholm, S., Ruud, B.O. (1994). 2-D finite difference elastic wave modeling including surface topog-
raphy. Geophys. Prosp. 42, 371–390.
Hestholm, S., Ruud, B.O. (2002). 3-D free-boundary conditions for coordinate-transform finite-
difference seismic modeling. Geophys. Prosp. 50, 463–474.
Higdon, R.L. (1991). Absorbing boundary conditions for elastic waves. Geophysics 56, 231–241.
Holberg, O. (1987). Computational aspects of the choice of operator and sampling interval for numer-
ical differentiation in large-scale simulation of wave phenomena. Geophys. Prosp. 35, 629–655.
Holliger, K., Robertsson, J.O.A. (1998). Effects of the shallow subsurface on the upper crustal seismic
reflection images. Tectonophysics 286, 161–169.
Igel, H., Mora, P., Riollet, B. (1995). Anisotropic wave propagation through finite-difference grids.
Geophysics 60, 1203–1216.
Ilan, A. (1977). Finite-difference modeling for P-pulse propagation in elastic media with arbitrary
polygonal surface. J. Geophys. 43, 41–58.
Ilan, A., Loewenthal, D. (1976). Instability of finite-difference schemes due to boundary conditions in
elastic media. Geophys. Prosp. 24, 431–453.
Ilan, A., Ungar, A., Alterman, Z.S. (1975). An improved representation of boundary conditions in
finite difference schemes for seismological problems. Geophys. J. R. Astron. Soc. 43, 727–745.
Ionescu, I.R., Campillo, M. (1999). Influence of the shape of the friction law and fault finiteness on
the duration of initiation. J. Geophys. Res. 104, 3013–3024.
Isaacson, E., Keller, H.B. (1966). Analysis of Numerical Methods. Wiley and Sons, New York.
Jastram, C., Behle, A. (1992). Acoustic modeling on a grid of vertically varying spacing. Geophys.
Prosp. 40, 157–169.
Jastram, C., Tessmer, E. (1994). Elastic modelling on a grid with vertically varying spacing. Geophys.
Prosp. 42, 357–370.
Jih, R.-S., McLaughlin, K.L., Der, Z.A. (1988). Free-boundary conditions of arbitrary polygonal
topography in a two-dimensional explicit elastic finite-difference scheme. Geophysics 53, 1045–
1055.
Kang, T.-S., Baag, Ch.-E. (2004a). Finite-difference seismic simulation combining discontinuous
grids with locally variable timesteps. Bull. Seismol. Soc. Am. 94, 207–219.
Kang, T.-S., Baag, Ch.-E. (2004b). An efficient finite-difference method for simulating 3-D seismic
response of localized basin structures. Bull. Seismol. Soc. Am. 94, 1690–1705.
Kay, I., Krebes, E.S. (1999). Applying finite element analysis to the memory variable formulation of
wave propagation in anelastic media. Geophysics 64, 300–307.
Kelly, K.R., Ward, R.W., Treitel, S., Alford, R.M. (1976). Synthetic seismograms: A finite-difference
approach. Geophysics 41, 2–27.
Klimeš, L. (1996). Accuracy of elastic finite differences in smooth media. Pure Appl. Geophys. 148,
39–76.
Koelbel, C., Loverman, D., Shreiber, R., Steele Jr., G., Zosel, M. (1994). The High-Performance For-
tran Handbook. MIT Press.
Komatitsch, D., Coutel, F., Mora, P. (1996). Tensorial formulation of the wave equation for modeling
curved interfaces. Geophys. J. Int. 127, 156–168.
Komatitsch, D., Tromp, J. (2003). A Perfectly Matched Layer (PML) absorbing condition for the
second-order elastic wave equation. Geophys. J. Int. 154, 146–153.
512 MOCZO ET AL.

Kosloff, R., Kosloff, D. (1986). Absorbing boundaries for wave propagation problems. J. Comput.
Phys. 63, 363–376.
Kristek, J., Moczo, P. (2003). Seismic wave propagation in viscoelastic media with material
discontinuities—a 3-D 4th-order staggered-grid finite-difference modeling. Bull. Seismol. Soc.
Am. 93, 2273–2280.
Kristek, J., Moczo, P., Archuleta, R.J. (2002). Efficient methods to simulate planar free surface in the
3-D 4th-order staggered-grid finite-difference schemes. Studia Geophys. Geod. 46, 355–381.
Kristek, J., Moczo, P., Irikura, I., Iwata, T., Sekiguchi, H. (1999). The 1995 Kobe mainshock simulated
by the 3-D finite differences. In: Irikura, K., Kudo, K., Okada, H., Sasatani, T. (Eds.), The Effects
of Surface Geology on Seismic Motion, vol. 3. Balkema, Rotterdam, pp. 1361–1368.
Kummer, B., Behle, A. (1982). Second-order finite-difference modeling of SH-wave propagation in
laterally inhomogeneous media. Bull. Seismol. Soc. Am. 72, 793–808.
Kummer, B., Behle, A., Dorau, F. (1987). Hybrid modelling of elastic-wave propagation in two-
dimensional laterally inhomogeneous media. Geophysics 52, 765–771.
Laws, R., Kragh, E. (2002). Rough seas and time-lapse seismic. Geophys. Prosp. 50, 195–208.
Lax, P.D., Wendroff, B. (1964). Difference schemes for hyperbolic equations with high order accuracy.
Commun. Pure Appl. Math. 27.
Lecomte, I., Gjøystdal, H., Maaø, F., Bakke, R., Drottning, Å., Johansen, T.-A. (2004). Efficient and
flexible seismic modelling of reservoirs: The HybriSeis concept. The Leading Edge 23, 432–437.
Levander, A.R. (1988). Fourth-order finite-difference P-SV seismograms. Geophysics 53, 1425–1436.
Levander, A.R. (1989). Finite-difference forward modeling in seismology. In: James, D.E. (Ed.), The
Encyclopedia of Solid Earth Geophysics. Van Nostrand Reinhold, pp. 410–431.
Liao, Z., Wong, H.L., Baipo, Y., Yifan, Y. (1984). A transmitting boundary for transient wave analysis.
Sci. Sinica A 27, 1063–1076.
Lindman, E.L. (1975). Free space boundary conditions for the time dependent wave equation. J. Com-
put. Phys. 18, 66–78.
Liu, H.-P., Anderson, D.L., Kanamori, H. (1976). Velocity dispersion due to anelasticity; implications
for seismology and mantle composition. Geophys. J. R. Astron. Soc. 47, 41–58.
Luo, Y., Schuster, G. (1990). Parsimonious staggered grid finite-differencing of the wave equation.
Geophys. Res. Lett. 17, 155–158.
Madariaga, R. (1976). Dynamics of an expanding circular fault. Bull. Seismol. Soc. Am. 67, 163–182.
Madariaga, R., Olsen, K., Archuleta, R. (1998). Modeling dynamics rupture in a 3-D earthquake fault
model. Bull. Seismol. Soc. Am. 88, 1182–1197.
Magnier, S.-A., Mora, P., Tarantola, A. (1994). Finite differences on minimal grids. Geophysics 59,
1435–1443.
Marchuk, G.I. (1982). Methods of Numerical Mathematics. Springer Verlag.
Marcinkovich, C., Olsen, K. (2003). On the implementation of perfectly matched layers in a three-
dimensional fourth-order velocity-stress finite difference scheme. J. Geophys. Res. B 108 (5),
2276, doi:10.1029/2002JB002235.
Marfurt, K.J. (1984). Accuracy of finite-difference and finite-element modeling of the scalar and elas-
tic wave equations. Geophysics 49, 533–549.
McDonal, F.J., Angona, F.A., Mills, L.R., Sengbush, R.L., van Nostrand, R.G., White, J.E. (1958).
Attenuation of shear and compressional waves in Pierre shale. Geophysics 23, 421–439.
Mikumo, T., Miyatake, T. (1987). Numerical modeling of realistic fault rupture processes. In: Bolt,
B.A. (Ed.), Seismic Strong Motion Synthetics. Academic Press, pp. 91–151.
Mitchell, A.R. (1969). Computational Methods in Partial Differential Equations. Wiley and Sons,
London.
Mitchell, A.R., Griffiths, D.F. (1994). The Finite Difference Method in Partial Differential Equations.
Wiley and Sons, New York.
Mittet, R. (2002). Free-surface boundary conditions for elastic staggered-grid modeling schemes. Geo-
physics 67, 1616–1623.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 513

Mizutani, H., Geller, R.J., Takeuchi, N. (2000). Comparison of accuracy and efficiency of time-domain
schemes for calculating synthetic seismograms. Phys. Earth Planet. Int. 119, 75–97.
Moczo, P. (1989). Finite-difference technique for SH-waves in 2-D media using irregular grids—
application to the seismic response problem. Geophys. J. Int. 99, 321–329.
Moczo, P. (1998). Introduction to Modeling Seismic Wave Propagation by the Finite-Difference
Method. Lecture Notes. Kyoto University. Available in pdf format at ftp://ftp.nuquake.sk/pub/
Papers.
Moczo, P., Bard, P.-Y. (1993). Wave diffraction, amplification and differential motion near strong
lateral discontinuities. Bull. Seismol. Soc. Am. 83, 85–106.
Moczo, P., Bystrický, E., Kristek, J., Carcione, J.M., Bouchon, M. (1997). Hybrid modeling of P-SV
seismic motion at inhomogeneous viscoelastic topographic structures. Bull. Seismol. Soc. Am. 87,
1305–1323.
Moczo, P., Kristek, J. (2005). On the rheological models used for time-domain methods of seismic
wave propagation. Geophys. Res. Lett. 32, L01306, doi:10.1029/2004GL021598.
Moczo, P., Kristek, J., Bystrický, E. (2001). Efficiency and optimization of the 3-D finite-difference
modeling of seismic ground motion. J. Comput. Acoust. 9, 593–609.
Moczo, P., Kristek, J., Gális, M. (2004a). Simulation of planar free surface with near-surface lateral
discontinuities in the finite-difference modeling of seismic motion. Bull. Seismol. Soc. Am. 94,
760–768.
Moczo, P., Kristek, J., Halada, L. (2000). 3-D 4th-order staggered-grid finite-difference schemes: Sta-
bility and grid dispersion. Bull. Seismol. Soc. Am. 90, 587–603.
Moczo, P., Kristek, J., Halada, L. (2004b). The Finite-Difference Method for Seismologists. An In-
troduction. Comenius University, Bratislava. Available in pdf format at ftp://ftp.nuquake.sk/pub/
Papers.
Moczo, P., Kristek, J., Vavryčuk, V., Archuleta, R.J., Halada, L. (2002). 3-D heterogeneous staggered-
grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging
of elastic moduli and densities. Bull. Seismol. Soc. Am. 92, 3042–3066.
Moczo, P., Labák, P., Kristek, J., Hron, F. (1996). Amplification and differential motion due to an
antiplane 2-D resonance in the sediment valleys embedded in a layer over the halfspace. Bull.
Seismol. Soc. Am. 86, 1434–1446.
Moczo, P., Lucká, M., Kristek, J., Kristeková, M. (1999). 3-D displacement finite differences and a
combined memory optimization. Bull. Seismol. Soc. Am. 89, 69–79.
Morton, K.W., Mayers, D.F. (1994). Numerical Solution of Partial Differential Equations. Cambridge
Univ. Press.
Mufti, I.R. (1985). Seismic modeling in the implicit mode. Geophys. Prosp. 33, 619–656.
Muir, F., Dellinger, J., Etgen, J., Nichols, D. (1992). Modeling elastic fields across irregular bound-
aries. Geophysics 57, 1189–1193.
Murphy III, W.F. (1982). Effects of partial saturation on attenuation in Massilon sandstone and Vycor
porous glass. J. Acoust. Soc. Am. 71, 1458–1468.
Nielsen, S., Carlson, J.M. (2000). Rupture pulse characterization: Self-healing, self-similar, expanding
solutions in a continuum model of fault dynamics. Bull. Seismol. Soc. Am. 90, 1480–1497.
Nielsen, S., Carlson, J.M., Olsen, K. (2000). Influence of friction and fault geometry on earthquake
rupture. J. Geophys. Res. B 105, 6069–6088.
Nielsen, P., If, F., Berg, P., Skovgaard, O. (1994). Using the pseudospectral technique on a curved grid
for 3-D acoustic forward modeling. Geophys. Prosp. 42, 321–341.
O’Brien, G., Hyman, M., Kaplan, S. (1951). A study of the numerical solution of partial differential
equations. J. Math. Phys. 29, 233–251.
Ohminato, T., Chouet, B.A. (1997). A free-surface boundary condition for including 3-D topography
in the finite-difference method. Bull. Seismol. Soc. Am. 87, 494–515.
Ottosen, N.S., Petersson, H. (1992). Introduction to the Finite Element Method. Prentice Hall.
Peng, C., Toksöz, M.N. (1994). An optimal absorbing boundary condition for finite difference model-
ing of acoustic and elastic wave propagation. J. Acoust. Soc. Am. 95, 733–745.
514 MOCZO ET AL.

Peng, C., Toksöz, M.N. (1995). An optimal absorbing boundary condition for elastic wave modeling.
Geophysics 60, 296–301.
Pierson, W.J., Moskowitz, L. (1964). A proposed spectral form for fully developed wind seas based
on the similarity theory of S.A. Kitaigorodskii. J. Geophys. Res. 69, 5181–5190.
Pitarka, A. (1999). 3-D elastic finite-difference modeling of seismic motion using staggered grids with
nonuniform spacing. Bull. Seismol. Soc. Am. 89, 54–68.
Pitarka, A., Irikura, K. (1996). Modeling 3-D surface topography by finite-difference method: Kobe-
JMA station site, Japan, case study. Geophys. Res. Lett. 23, 2729–2732.
Pratt, R.G. (1990). Inverse theory applied to multi-source cross-hole tomography. Part 2: Elastic wave-
equation method. Geophys. Prosp. 38, 311–329.
Pratt, R.G., Shin, C., Hicks, G.J. (1998). Gauss–Newton and full Newton methods in frequency-space
seismic waveform inversion. Geophys. J. Int. 133, 341–362.
Richtmyer, R.D., Morton, K.W. (1967). Difference Methods for Initial Value Problems. Wiley and
Sons, New York. (reprinted by Kreiger, New York, 1994).
Robertsson, J.O.A. (1996). A numerical free-surface condition for elastic/viscoelastic finite-difference
modeling in the presence of topography. Geophysics 61, 1921–1934.
Robertsson, J.O.A., Blanch, J.O., Symes, W.W. (1994). Viscoelastic finite-difference modeling. Geo-
physics 59, 1444–1456.
Robertsson, J.O.A., Chapman, C.H. (2000). An efficient method for calculating finite-difference seis-
mograms after model alterations. Geophysics 65, 907–918.
Robertsson, J.O.A., Coates, R.T. (1997). Finite-difference modeling of Q for qP- and qS-waves in
anisotropic media. In: 67th Annual International Meeting, Society of Exploration Geophysicists,
pp. 1846–1849. Expanded Abstracts.
Robertsson, J.O.A., Holliger, K. (1997). Modeling of seismic wave propagation near the Earth’s sur-
face. Phys. Earth Planet. Int. 104, 193–211.
Robertsson, J.O.A., Laws, R., Chapman, C.H., Vilotte, J.-P., Delavaud, E. (2006). Modelling of scat-
tering of seismic waves from a corrugated rough sea surface: A comparison of three methods.
Geophys. J. Int., in press.
Robertsson, J.O.A., Levander, A., Holliger, K. (1996). A hybrid wave propagation simulation tech-
nique or ocean acoustic problems. J. Geophys. Res. 101, 11225–11241.
Robertsson, J.O.A., Ryan-Grigor, S., Sayers, C., Chapman, C.H. (2000). A finite-difference injection
approach to modeling of seismic fluid flow monitoring. Geophysics 65, 896–906.
Rodrigues, D. (1993). Large scale modelling of seismic wave propagation. Ph.D. Thesis, Ecole Cen-
trale Paris.
Rodrigues, D., Mora, P. (1993). An efficient implementation of the free-surface boundary condition
in 2-D and 3-D elastic cases. In: 63th Annual International Meeting, Society of Exploration Geo-
physicists, pp. 215–217. Expanded Abstracts.
Saenger, E.H., Gold, N., Shapiro, S.A. (2000). Modeling the propagation of elastic waves using a
modified finite-difference grid. Wave Motion 31, 77–92.
Saenger, E.H., Bohlen, T. (2004). Finite-difference modeling of viscoelastic and anisotropic wave
propagation using the rotated staggered grid. Geophysics 69, 583–591.
Schneider, J.B., Ramahi, O.M. (1998). The complementary operators method applied to acoustic
finite-difference time-domain simulations. J. Acoust. Soc. Am. 104, 686–693.
Schoenberg, M., Muir, F. (1989). A calculus for finely layered anisotropic media. Geophysics 54,
581–589.
Shtivelman, V. (1984). A hybrid method for wave field computation. Geophys. Prosp. 32, 236–257.
Shtivelman, V. (1985). Two-dimensional acoustic modelling by a hybrid method. Geophysics 50,
1273–1284.
Slawinski, R.A., Krebes, E.S. (2002). The homogeneous finite-difference formulation of the P-SV-
wave equation of motion. Studia Geophys. Geod. 46, 731–751.
Sochacki, J.S., George, J.H., Ewing, R.E., Smithson, S.B. (1991). Interface conditions for acoustic
and elastic wave propagation. Geophysics 56, 168–181.
 THE FINITE-DIFFERENCE TIME-DOMAIN METHOD 515

Smith, W.D. (1974). A non-reflecting plane boundary for wave propagation problems. J. Comput.
Phys. 15, 492–503.
Spencer Jr., J.W. (1981). Stress relaxation at low frequencies in fluid-saturated rocks. J. Geophys.
Res. 86, 1803–1812.
Stead, R.J., Helmberger, D.V. (1988). Numerical-analytical interfacing in two dimensions with appli-
cations to modeling NTS seismograms. In: Aki, A., Wu, R.S. (Eds.), Scattering and Attenuation
of Seismic Waves. Birkhauser, Basel, Switzerland, pp. 157–193.
Strang, G., Fix, G.J. (1973). An Analysis of the Finite Element Method. Prentice Hall, Englewood
Cliffs, NJ.
Taflove, A., Hagness, S.C. (2005). Computational Electrodynamics: The Finite-difference Time-
Domain Method. Artech House.
Takenaka, H., Furumura, T., Fujiwara, H. (1998). Recent developments in numerical methods for
ground motion simulation. In: Irikura, K., Kudo, K., Okada, H., Sasatani, T. (Eds.), The Effects of
Surface Geology on Seismic Motion, vol. 2. Balkema, Rotterdam, pp. 91–101.
Takeuchi, N., Geller, R.J. (2000). Optimally accurate second order time-domain finite difference
scheme for computing synthetic seismograms in 2-D and 3-D media. Phys. Earth Planet. Int. 119,
99–131.
Takeuchi, N., Geller, R.J. (2003). Accurate numerical methods for solving the elastic equation of
motion for arbitrary source locations. Geophys. J. Int. 154, 852–866.
Tal-Ezer, H., Carcione, J.M., Kosloff, D. (1990). An accurate and efficient scheme for wave propaga-
tion in linear viscoelastic media. Geophysics 55, 1366–1379.
Tessmer, E. (2000). Seismic finite-difference modeling with spatially varying time steps. Geo-
physics 65, 1290–1293.
Tessmer, E., Kosloff, D. (1994). 3-D elastic modeling with surface topography by a Chebyshev spec-
tral method. Geophysics 59, 464–473.
Tessmer, E., Kosloff, D., Behle, A. (1992). Elastic wave propagation simulation in the presence of
surface topography. Geophys. J. Int. 108, 621–632.
Thomson (Lord Kelvin), W. (1856). Elements of a mathematical theory of elasticity. Philos. Trans. R.
Soc. 166, 481–498.
van Manen, D.J., Robertsson, J.O.A., Curtis, A. (2005). Modeling of wave propagation in inhomoge-
neous media. Phys. Rev. Lett. 94, 164301–164304.
Virieux, J. (1984). SH-wave propagation in heterogeneous media: Velocity-stress finite-difference
method. Geophysics 49, 1933–1957.
Virieux, J. (1986). P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference
method. Geophysics 51, 889–901.
Wang, Y., Xu, J., Schuster, G.T. (2001). Viscoelastic wave simulation in basins by a variable-grid
finite-difference method. Bull. Seismol. Soc. Am. 91, 1741–1749.
Wang, T., Tang, X. (2003). Finite-difference modeling of elastic wave propagation: A nonsplitting
perfectly matched layer approach. Geophysics 68, 1749–1755.
Xu, T., McMechan, G.A. (1995). Composite memory variables for viscoelastic synthetic seismograms.
Geophys. J. Int. 121, 634–639.
Xu, T., McMechan, G.A. (1998). Efficient 3-D viscoelastic modeling with application to near-surface
land seismic data. Geophysics 63, 601–612.
Yomogida, K., Etgen, J.T. (1993). 3-D wave propagation in the Los Angeles basin for the Whittier–
Narrows earthquake. Bull. Seismol. Soc. Am. 83, 1325–1344.
Zahradník, J. (1995). Comment on ‘A hybrid method for estimation of ground motion in sedimentary
basins: Quantitative modeling for Mexico City’ by D. Fäh, P. Suhadolc, St. Mueller and G.F. Panza.
Bull. Seismol. Soc. Am. 85, 1268–1270.
Zahradník, J., Hron, F. (1992). Robust finite-difference scheme for elastic waves on coarse grids.
Studia Geophys. Geod. 36, 1–19.
Zahradník, J., Moczo, P., Hron, F. (1993). Testing four elastic finite-difference schemes for behaviour
at discontinuities. Bull. Seismol. Soc. Am. 83, 107–129.
516 MOCZO ET AL.

Zahradník, J., Moczo, P. (1996). Hybrid seismic modeling based on discrete-wavenumber and finite-
difference methods. Pure Appl. Geophys. 148, 21–38.
Zahradník, J., Priolo, E. (1995). Heterogeneous formulations of elastodynamic equations and finite-
difference schemes. Geophys. J. Int. 120, 663–676.
Zeng, X. (1996). Finite difference modeling of viscoelastic wave propagation in a generally heteroge-
neous medium in the time domain, and a dissection method in the frequency domain. Ph.D. Thesis,
University of Toronto.
Zhang, J. (1997). Quadrangle-grid velocity-stress finite-difference method for elastic-wave-
propagation simulation. Geophys. J. Int. 131, 127–134.
Zhang, Ch., Symes, W.W. (1998). Fourth order, full-stencil immersed interface method for elastic
waves with discontinuous coefficients. 1998 SEG Expanded Abstracts.
Zienkiewicz, O.C., Taylor, R.L. (1989). The Finite Element Method, vol. 1, 4th edition. McGraw–Hill.