3D Finite-Difference Method Using
Discontinuous Grids
Shin Aoi and Hiroyuki Fujiwara

Abstract

We have formulated a 3D finite-difference method (FDM) using
discontinuous grids, which is a kind of multigrid method. As long as uniform grids are used,
the grid size is determined by the shortest wavelength to be calculated, and thisconstitutes a significant constraint on the introduction of low velocity layers. We use
staggered grids which consist of, on the one hand, grids with fine spacing near the surface
where the wave velocity is low, and, on the other hand, grids whose spacing is three times
coarser in the deeper region. In each region, we calculated the wavefield using velocitystress formulation of second order accuracy and connected these two regions using linear
interpolations. The second order finite-difference (FD) approximation was used for
updating. Since we did not use interpolations for updating, the time increments were thesame in both regions. The use of discontinuous grids adapted to the velocity structure
resulted in a significant reduction of computational requirements, which is model dependent
but typically one fifth to one tenth, without a marked loss of accuracy.

Helmberger and Vidale, 1988; Frankel, 1993;
Graves, 1996), which is a particular issue for
applying finite-difference method (FDM) in
seismology. The FDM is one of the mostpractical waveform simulation methods in use
today.
The improvement in computer capacities
has made it possible to carry out simulations ofa three dimensional wavefield with realistic
velocity models on a large scale, such as theones for the Kanto Plane （e.g. Sato et al., 1998;
Sugawara et al., 1997) and the Los Angeles
Basin （e.g. Yomogida and Etgen, 1993; Olsen
and Archuleta, 1996; Graves, 1998). However,despite its considerable influence on waveforms,
the low velocity layer near the surface cannot
be incorporated into such models. For example,
in the Kobe area where extensive damage
occurred in the 1995 Hyogoken-Nanbu
Earthquake, numerous geophysical explorations
such as a reflection survey, refraction survey,
and microtremor observation were performed,
and detailed 3D seismic wave velocity
structures have been proposed (e.g. Huzita,

Introduction
As a method of seismic wave simulation,
the finite-difference (FD) approximation has
frequently been used to solve equations of
motion numerically for a couple of decades (e.g.
Boore, 1972; Kelly et al., 1976), and theformulation using staggered grids is commonly
employed at present (e.g. Virieux, 1984, 1986;
Levander, 1988; Graves, 1996).
Many
researches have been carried out, such as theresearch of free boundary conditions on thesurface (e.g. Vidale and Clayton, 1986; Stacey,
1994; Pitarka and Irikura, 1996; Ohminato andChouet, 1997), elastic and liquid medium
boundary conditions as boundary conditions on
the seabed (e.g. Okamoto, 1996), non-reflecting(e.g. Cerjan et al., 1985) and absorbing (e.g.
Clayton and Engquist, 1977; Stacey, 1988;
Higdon, 1991) boundary conditions to avoid thereflected waves from the boundary of a finite
computational region, as well as the introductionof a double couple point source (e.g. Alterman
and Karal, 1968; Vidale and Helmberger, 1987;1

region into small grids unnecessarily. In thismethod, the waveforms are first calculated by a
2D FDM using the structure model without thesurface layer of low velocity. Then the effectsof this layer on the waveforms are consideredusing the convolution of 1D transfer functions,evaluated by the propagator-matrix technique(Haskell, 1953). However, this method merely
evaluates the influence of the shallow layer in
an approximate way. Another approach is totake a coarser grid spacing by enhancing theaccuracy of FD approximation, using suchmethods as an FDM with spatial difference of ahigher order (e.g. Yomogida and Etgen, 1993) ora pseudospectral method (e.g. Furumura, 1992).
However, the coarse grid spacing does notenable the modeling of detailed parts of thestructure, and the computation accuracy is not
sufficient in structures having discontinuities
with a high contrast. In order to achieve a highlevel of accuracy in computation, it is necessary
to use sufficiently small grids. Thus, if we wish
to calculate waveforms by the FDM using
models that include near-surface layers withlow velocity, we need to employ non-uniform
grids that are adapted to the velocity structure.
There are two types of non-uniform grids,
continuous and discontinuous.
Continuous grids are the grids with the
optimal distribution of grid spacing achieved bycontinuous mapping.
Examples of methods
using continuous grids include refining the grid
spacing in the vicinity of the free surface (e.g.
Moczo, 1989; Carcione, 1992), reducing the grid
spacing in the vicinity of the fault plane (e.g.
Mikumo et al., 1987; Mikumo and Miyatake,1993), refining the grid spacing within a given
region (e.g. Pitarka, 1999) and making grids thatgenerally follow the interfaces of media (e.g.
Fornberg, 1988; Nielsen et al., 1994). Thesemethods are free of artificial computational
errors resulting from sudden changes in the grid
spacing, since they allow for a continuous
reduction in grid spacing. On the other hand,their shortcoming is that the number of grid
points can be changed only along the coordinate
axis.
With regard to discontinuous grids, thegrid system consists of several regions, each ofthem having a uniform grid. This is a kind of

1996). Ground motion simulations performed
by Iwata et al. (1998) and Pitarka et al. (1998)
using a 3D FDM and models of the undergroundstructure in the Kobe region successfully
reproduced an extension of the severely
damaged band of the earthquake. Meanwhile,
due to the high values of the adopted shear
wave velocity of the near-surface sediments,
the amplitude of the simulated strong motions
were smaller than those observed.
As long as uniform grids are employed,
their size is determined by the shortest
wavelength to be calculated. Thus, the entire
region must be divided into small grids evenwhen the layer of low velocity only occupies asmall part. This considerably increases thecomputational requirements in terms of timeand memory. Using the 3D FDM, in order tocalculate up to the frequency N times higher, or
in order to introduce low velocity layers inwhich the wave velocity is N times smaller (i.e.
reduce the grid size to 1/N), N3 more memory,
and N4 more computation time are required.
Therefore, we cannot depend on the progress ofcomputingcapacities
exclusively
forcomputations on a large scale. Moreover, inorder to perform the underground structure
inversion (e.g. Aoi et al., 1995, 1997), it isnecessary to calculate the waveforms manytimes, so a method enabling quick calculationsof waveforms is required.
When the wavefield is calculated by the
FDM, the grid size near the surface should be
as small as possible for the following reasons:- the wave velocity of the near-surface
sediment is relatively low;
- the underground structure is extremely
heterogeneous;- free surface boundary conditions that must be
imposed on the free surface tend to be unstable,
and in many cases we are interested inwaveforms at the surface;- when the grid is staggered the wave-field
variables are not defined at the same position;- the energy of the surface waves isconcentrated near the surface, and its groupvelocity is slower than S-wave velocity.
Pitarka et al. (1996) represents anattempt to evaluate the influence of the surface
layer on the wavefield without dividing the2

of grids in the horizontal direction cannot be
reduced, and that accordingly, the grid spacingin the horizontal direction cannot becomecoarser. In the present paper, we present an
FD technique that is based on a discontinuous
grid.
We also analyze its accuracy bycomparisons with waveforms produced by thediscrete-wavenumber
method(DWNM)(Bouchon, 1981; Schmidt and Tango, 1986) andby the FDM using uniform grids.

multigrid technique which is already in use in
the field of fluid mechanics (e.g. McBryan et al,1991). At the boundary of each region, theFDM has to be formulated in a way thatmaintains the continuity of the wavefield.
Examples of waveform simulations by the FDM
using discontinuous grids include reducing thegrid spacing in the vicinity of the free surface(e.g. Moczo et al., 1996), refining the grid
spacing in the vicinity of the borehole (e.g. Falket al., 1996; Kessler and Kosloff, 1991) andavoiding grid spacing which is too small in thecentral part of the cylindrical coordinate (e.g.
Furumura et al., 1998). However, these are all
examples of 2D problems. Examples of hybrid
methods using both grid systems include
Jastram and Tessmer (1994).
Numerous issues of seismology deal withstructures in which the wave velocity is lower inthe shallower part and higher in the deeper part.
In such cases, grids that are discontinuous in
the vertical direction are often advantageous.
This is due to the fact that as long as
continuous grids are used, even in the deep part
where the velocity is much higher, the number

MethodFormulation of the FDM with Discontinuous
Grids
We used a discontinuous grid thatconsists of two regions with different grid
spacing. Figure 1 shows the unit cell of thegrid and the 3D discontinuous grid together with
its cross-sections. The grid spacing of Region
I is small (the grid spacings in directions x, y
and z are ∆x , ∆y and ∆z , respectively), whereas
the grid spacing of Region II is three times

Fig. 1: (center）3D discontinuous grid system and a unit cell for staggered grids (inside the circle).
(left) Two transections on the top and at the bottom of the overlapping region of Regions I and II,
where the elimination or the insertion of grid points are necessary. (right) Two profiles of the
discontinuous grid. The arrows A-E show the overlapping region of Regions I and II, and thedetails of the interpolation are given in Table 1.

3

coarser (the grid spacings in directions x, y andz are 3 ∆x , 3 ∆y and 3 ∆z , respectively). In eachregion, a 3D staggered grid FDM of secondorder accuracy in time and space is employed.
The discretized equations of motion are
given by

Table 1
How to up-date the time step of variables oneach planevx in++11//22, j , k = vx in+−11//22, j , knτ nτ xy in, j +1, k − τ xy in, j , k
xx i +1, j , k − τ xx i , j , k+ ∆t × b 
+

∆x
∆y


+

τ xz in, j , k +1 − τ xz in, j , k 
∆z





v y in,+j1+/12/ 2, k = v y in,−j1+/12/ 2, k

+





n
τ n
τ yz in, j +1, k − τ yz in, j , k
xz i +1, j , k − τ xz i , j , k+ ∆t × b 
+

∆x
∆y


τ zz in, j , k +1 − τ zz in, j , k 
∆z





and the discretized stress-strain relations are
represented asvx in++11//22, j , k − vx in−+11//22, j , k




∆x

τ xxin,+j1,k = τ xxin, j , k + ∆t × (λ + 2µ ) v n +1 / 2 − v n +1 / 2
vz in,+j 1, k/+21 / 2 − vz in,+j 1, k/−21 / 2 
y i , j −1 / 2, k
 y i , j +1 / 2, k
+ λ
+

∆y
∆z






τ yy in,+j1,k = τ yy in, j , k + ∆t × (λ + 2µ )v y in,+j 1+1/ 2/ 2, k − v y in,+j 1−1/ 2/ 2, k∆y vx n +1 / 2 − v x n +1 / 2
i −1 / 2, j , k i +1 / 2, j , k

+ λ



∆x

+

vz in,+j 1, k/+21 / 2 − vz in,+j 1, k/−21 / 2 


∆z



vz in,+j 1, k/+21 / 2 − vz in,+j 1, k/−21 / 2
τ zz in,+j1,k = τ zz in, j , k + ∆t × (λ + 2µ )∆z

 v n +1 / 2 − v n +1 / 2v y in, +j 1+1/ 2/ 2, k − v y in,+j1−1/ 2/ 2, k  x i +1 / 2, j , k
x i −1 / 2, j , k

+
+ λ

x
∆
∆y





τ xy in++11/ 2, j +1 / 2, k = τ xy in+1 / 2, j +1 / 2, k + ∆t × µ
 v n +1 / 2v y in++11, /j +21 / 2, k − v y in, +j 1+1/ 2/ 2, k − v x in++11// 22, j , k x i +1 / 2, j +1, k

+


∆y
∆x



τ xz in++11/ 2, j , k +1 / 2 = τ xz in+1 / 2, j , k +1 / 2 + ∆t × µ
 v n +1 / 2− v x in++11// 22, j , kv z in++11, /j ,2k +1 / 2 − v z in, +j 1, k/+21 / 2  x i +1 / 2, j , k +1

+


∆z
∆x



τ yz in, +j 1+1 / 2, k +1 / 2 = τ yz in, j +1 / 2, k +1 / 2 + ∆t × µ
 v n +1 / 2− v y in, +j 1+1/ 2/ 2, kv z in, +j 1+1/ ,2k +1 / 2 − v z in, +j ,1k/+21 / 2  y i , j + 1 / 2 , k +1

+


∆z
∆y



vx , v y , vz

representt xx , t yy , t zz , t xy , t xz , t yz

theparticle

FDM

----

A

FDM

Interpolation

B and C

FDM

----

D

Interpolation

FDM

E

----

FDM

Region II

----

FDM

are the body force components.

and ∆z represent the grid spacing in
the x , y and z directions, respectively. ∆tbis thedenotes the time increment.
buoyancy (inverse of density), and λ and µ
are Lame constants. We used effective media
parameters which were calculated using Grave's
formulation (Graves, 1996).
An elastic
attenuation is introduced in the same way as inGraves (1996).
Only the field variables (velocity andstress components) that are adjacent to thevariable to be updated are required to update
the wavefield from the time level t n to t n+1 .
As it is clear from equations (1) and (2), when
the FD approximation of second-order accuracy
is employed for the spatial derivatives, the field
variables that are within the distance of halfgrid spacing in the x , y and z directions are
required. Accordingly, only the field variablesat the bottom plane of Region I and the top
plane of Region II cannot be calculated bystaggered grid FD operators (Fig. 1 and Table 1).
Therefore, the field variables of these twoplanes must be calculated by inserting or
eliminating the grids of the other region and byusing interpolation of the wavefield across thetwo regions.

vz in,+j1, k/ +21 / 2 = v z in,−j1,k/ +21 / 2+

Region I∆x ,τ yz in, j , k +1 − τ yz in, j , k ∆zRegion IIfx, fy, fz

(1)

n
τ n
τ yyin, j +1, k − τ yyin, j , k
xy i +1, j , k − τ xy i , j , k+ ∆t × b 
+

∆x
∆y


Region I(2).
velocity,

are the stress components, and4∆yInterpolation
In this section we explain the techniqueused to interpolate the wavefield at theboundary between Regions I and II.
Regions I and II overlap in the verticaldirection, covering the distance of 3∆z / 2 . The
field variables at the top plane of Region II
cannot be obtained through the FD solutions.
However, since the locations of the grid points
at the top plane of Region II are identical tothose of Region I, the field variables of thelatter can be employed as those of the former(Fig. 1).
The field variables at the bottom ofRegion I are obtained by an interpolation
scheme, using the field variables in Region II
obtained by the FDM. What is interesting hereis that the interpolations of all field variables
are carried out within one horizontal plane, andthat apart from these interpolations, the time
updates of variables are carried out by the FD
calculations (Table 1).
The linear functions
a 0 ( x) = 1 − x
a1 ( x ) = x(0 ≤ x ≤ 1)

Fig. 2: Grid points location on the plane forinterpolation, where (I, J) and (i, j) are local
numberings for the interpolation.
techniques used to avoid boundary reflections
are absorbing (e.g. Clayton and Engquist, 1977;
Stacey, 1988; Higdon, 1991), or non-reflecting(e.g. Cerjan et al., 1985) boundary conditions.
The former is a method of realizing theboundary conditions that make the reflectionsof the body wave with a specific wavenumber
vanish at a grid or a few grids in the vicinity ofthe boundary. The latter is a method of
eliminating the reflected waves through theirgradual attenuation by setting an absorbing
region outside the boundary.
Absorbingboundary conditions do have certain advantages
in terms of computational requirements.
However, apart from the case of the body wave
with a specific wavenumber (normally a verticalincident wave), of which the reflections vanish,the method cannot realize a perfect absorption.
On the contrary, approximately 10 to 20 %additional memory and computation are required
to realize the non-reflecting boundary
conditions.
This method is capable ofabsorbing the waves almost completely
regardless of whether it is a body wave (of any
wavenumber) or a surface wave. Here we used
the non-reflecting boundary conditions of
Cerjan et al. (1985).
In Cerjan et al. (1985), Gaussian functions
given by
W = exp(−α ( J 0 − j ) 2 ), ( j = 1,2,⋅ ⋅ ⋅ ⋅ ⋅⋅, J 0 )
(5)
and

(3)

are used for the interpolation.
Table 2indicates the weights for the interpolation
obtained from these functions (equation (3)).
These weights correspond to the points, x=0,
1/3, 2/3 and 1, when the grid reduction factor
is 3.
The variables must be interpolated on thex-y plane at the bottom of Region I (Fig. 1,
bottom left), and the grids are positioned asshown in Figure 2, where (I, J) and (i, j) are local
indexings for the interpolation.
The field
variables obtained by the interpolation scheme
are
1

1

u ij = ∑ ∑ 1α iI, ,jJ U I , J(i, j = 0,1,2,3; I , J = 0,1)

I =0 J =0(4)

α iI, ,jJ = aiI ⋅ a Jj

where U I , J and ui , j indicate the field variables
in Regions I and II, respectively.

v np+1 / 2 = W ⋅ v np+1 / 2

Boundary Conditions

τ npq = W ⋅τ npq ( p, q = x, y, z )(6)

are used to attenuate the wavefield close to theboundary. According to Cerjan et al. (1985),
α = 0.015 and J 0 = 20 are the most appropriate
values. Therefore these values are employed

Boundary conditions are imposed to avoid
artificial reflections from the boundary of thefinite computational region. The most popular5

in Region II. In Region I, α = 0.005 and J 0 = 60
are employed because the process ofattenuation must be continuous with Region II.
Though the number of grids J 0 = 60 inRegion I is relatively large, the memory andcomputation time required to realize the nonreflecting boundary condition are negligible
compared to those saved from the use of the
discontinuous grids.
Moreover, as thecomputation scale increases, the ratio ofmemory and computation used for the absorbing
region to those used for the entire calculation
decreases accordingly.
For example, in amodel with 2000*2000 grids in two horizontal
directions, the ratio that the absorbing region
occupies is less than 13 % of the entire region,and consequently, the increase of thecomputational requirement is hardly significant.
However, in a case where the computational
region is extremely flat, this ratio may become
too significant to be neglected. One approach
to solve this problem is to make the grid
spacing in the absorbing region in Region I

coarser, so that it will be identical to the grid
spacing of Region II （ 3∆x , 3∆y ）, thus reducing
the number of grid points J 0 to 20. However,this approach is limited to structures in which
the wave velocity is high near the absorbing
boundary of Region I.
Stability Conditions
The stability condition for the constantgrid spacing FD technique isV p ∆t

1

∆x 2

+

1

∆y 2

+

1

∆z 2

<1

(7).

In the present method, this condition must be
satisfied in both Regions I and II, because weuse the constant grid spacing FD technique ineach region.
Since we do not use
interpolations for updating the wavefield, thetime increments are the same in both Regions I
and II. This means that we use the minimum
values of the time increments ( ∆t ) determinedby equation (7) for both regions.

Table 3
Physical parameters of the structure model

Structure

Shape

Case 1Case 2

Horizontally stratified structure

3-D basin structure----------------------

Vp

2.4 km/s

Vs

0.8 km/s

Density

1.8 km/s

Thickness/Max.
depth

1.0 kmVp

4.3 km/s

Vs

2.5 km/s

Density

2.5 g/cm3

Sediment

Rock

6

Paraboloid of revolution
Diameter 20 km
Center x = -0.05 km, y = 0.0kmTable 4
Source parameters
Case 1

0.46257
0.434450.351930.24104
0.19161
0.26249
0.29423
0.27459
0.273470.264380.25167
0.246370.23867
0.228840.21741

Case 20.19199

Sourcelocation

0.16037x=-1.2 km, y=-1.2 km, z= 9.4 km

0.15207
0.162040.17688

Source timefunction

Ricker wavelet ( Tc = 3 sec )Source type

Single force
point source( x-direction )0.18951
0.198270.204200.20929

Double couple point
source (strike=60o,
rake=30o, dip=60o)

0.21498
0.221800.229700.46830
0.439210.35424
0.23978
0.18933
0.259260.29063
0.27144
0.27155
0.25987
0.253880.248450.24055
0.23056
0.21887
0.20609
0.192940.180290.169160.160390.154630.151960.15201
0.154200.15786
0.162450.167470.17261
0.17760
0.182290.18653
0.19027
0.19352
0.19635
0.198780.20087
0.202710.204390.205990.20757
0.209190.21089
0.212700.21464
0.21671
0.21892
0.221250.223710.226270.22894

Examples of computation
The method proposed in this paper wasused to calculate the waveforms in a 1D
structure (horizontally stratified structure) andin a 3D basin structure. The results arecompared with those obtained with a staggered
grid FDM using uniform grids, in order todemonstrate the accuracy of the proposed
method. The results obtained with the 1D
structure are also compared to those of DWNM.

Fig. 3: Velocity waveform of the x-component at theobservation points on the z-axis ( 0 km ≤ z ≤ 4.9km , withspatial interval of 100 m. This interval is 300 m in
Region II) in the 1D structure, calculated (a) by the
FDM using discontinuous grids and (b) by the FDM
using a uniform grid. The figures on the right are themaximum absolute value of amplitude (hereinafter
maximum amplitude) of each waveform. On the left,
the schematic of the structure employed is shown.
The arrow indicates the boundary between the tworegions.

1D structure case
The structure is horizontally stratified andconsists of two homogeneous layers theparameters of which are shown in Table 3. We
use a single force point source acting in the xdirection which is located at (-1.2 km, -1.2 km,
9.4 km). Its source time function is a Ricker
wavelet with a characteristic period of 3 sec.
(Table 4).
Regardingthediscontinuous
grids
employed in this calculation, their grid spacing is
100 m and 300 m in Regions I and II,
respectively. The depths of Regions I and II
are 1.5 km (15 grid points) and 18 km (60 grid
points), respectively.
The time increments derived from thestability condition, equation (7), in Regions I and
II are ∆t < 0.01343 sec. and ∆t < 0.04049 sec.,
respectively.
Though the grid spacing inRegion II is three times coarser than in Region I
and the time increment can be three times

longer accordingly, in our study the same timeincrement ∆t = 0.0125 sec. is employed in bothregions. The use of a longer time increment inRegion II and the interpolations of time sampling
can further reduce the computation time.
Using a Ricker wavelet with acharacteristic period of 3 sec. the minimum
number of grids per wavelength in the entire
computational region is 12.
In the computation using a continuous
grid, the grid size is 100 m (155 grids in verticaldirection) in the entire region of computation,and the time increment is identical to that of
the computation with the discontinuous grid.

7

Vp=4.3km
Vs=2.5km
ρ =2.5g/cm3
1.0km

0.04772
0.05315
0.058360.06291
0.066710.070110.073710.07786
0.082200.08630
0.090100.093850.09758
0.10101
0.10389
0.106220.10815
0.10973
0.11689
0.12555
0.135030.14528
0.156260.167810.17976
0.192030.204420.216800.229210.24149
0.25383
0.266210.27877
0.291410.304210.321190.34173
0.36192
0.38103
0.39874
0.41462
0.428460.440150.44966
0.457140.46257
0.466140.467810.46762
0.46559

0

4

8

12

16

20

24

28

Vp=4.3km
Vs=2.5km
ρ =2.5g/cm3
1.0km▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼

Vp=2.4km
Vs=0.8km
ρ =1.8g/cm3
Vp=2.4km
Vs=0.8km
ρ =1.8g/cm3(b)▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼

(a)

32

Fig. 4: Velocity waveforms of the x-component at
the observation points on the x-axis( −2.0km ≤ x ≤ 22.5km , with an interval of 500 m) in
the 1D structure, calculated (a) by the FDM using
discontinuous grids and (b) by the FDM using
uniform grids. The figures on the right are themaximum amplitude of each waveform. On theleft, the schematic of the structure employed isshown.

Figure 3 compares synthetic seismograms
calculated at a vertical receivers array with theproposed scheme using a discontinuous grid andthe uniform grid scheme. The array covers thedepth range between 0.0km and 4.9km and theinterval between stations is 100m (300minterval for Regions II). There is an excellent
match between the waveforms computed withthe two grids. This indicates a high level ofaccuracy of the linear interpolation at theboundary between the two Regions.
The
accuracy is also confirmed by the continuous
propagation of the wave passing the boundary(shown by the arrow) which is characterized by

T [sec]0.05087
0.056420.062440.06791
0.073240.077930.08075
0.08358
0.08863
0.09394
0.09685
0.09806
0.100210.10409
0.108480.11215
0.11487
0.116960.121490.129180.13759
0.146770.15668
0.167230.178290.189700.20135
0.213110.22487
0.23658
0.248300.260140.272160.28434
0.29901
0.319930.34106
0.36185
0.381720.400170.41678
0.431390.443840.454170.46225
0.46830
0.472200.474020.47379
0.47152

0

4

8

12

16

20

24

28

32

T [sec]Fig. 5: The results of computation, by the FDM using discontinuous grids (dashed line), the FDM using uniform
grid (thin line) and the DWNM (thick line). The waveforms in the lower part of the figure (dot line) show thedifference in the results by the FDMs using the uniform grids and those using discontinuous grids. On theright of each waveform, the maximum and minimum amplitudes are indicated.

8

accuracy of the former one.

the absence of artificial reflected waves or
disturbed waves (Fig. 3 (a)).
Figure 4 shows the comparison ofsynthetic seismograms calculated at receiverson the free surface. The receivers are aligned
along the x-axis (-2.0km to 22.5km) and their
spatial interval is 500m. As in the previous
case, the comparison between the waveforms
computed with the two techniques is very good.
In order to examine in detail the errors of
the proposed scheme, the waveforms obtainedby the FDM using a discontinuous grid, the FDM
using uniform grids and the DWNM, at theobservation points shown by circles in Figures 3
and 4, are shown in Figure 5. The bottom
traces of each figure are the differences
between the result from the FDM using a
discontinuous grid and the one using uniform
grid. As shown in Figure 5, all three methods
resulted in amplitude and arrival times for allphases that are almost identical.
At allreceivers, the differences in amplitude is within
15 % of that calculated with the FDM using
uniform grids except for the observation points
which are far from the epicenter and whose
amplitude
issmall.
Moreover,thecorresponding results of the FDM and the
DWNM indicate a sufficiently high level of

3D basin structure caseHere, synthetic seismograms calculatedwith the proposed FDM using discontinuous gridare compared to those of an FDM using uniform
grid. In order to check the accuracy of thepresent method in the case of a more complex
structure we used a basin model (Table 3).
The basin model consists of a homogenous half
space (bedrock) and a homogenous sedimentarylayer. The basin geometry is a half paraboloidwith a diameter of 20km and maximum depth of
1.0 km. The center of the paraboloid is at x=0.05 km and y=0.0 km. The grids and source
location employed for the computation are
identical to those used in the 1D case. The
source is a double couple point source with thestrike, rake and dip angles of 60o, 30o and 60o,
respectively (Table 4).
Figures 6 (a) and (b) show the comparisonof the synthetic seismograms calculated with
the FDM using both discontinuous grids anduniform grid at observation points having thesame locations as in the 1D case.

0.573220.61369
0.64475
0.65272
0.63218
0.58256
0.50524
0.423310.305170.141770.12941
0.12398
0.11738
0.10991
0.101880.08625
0.079640.08764
0.106230.122440.13406
0.142400.14995
0.158780.171450.18657
0.20358

0.54492
0.581170.608530.61424
0.59358
0.54587
0.47267
0.39215
0.28309
0.13445
0.12195
0.11673
0.110420.10328
0.095610.087900.080770.077700.07656
0.076700.077960.080150.08304
0.088080.094230.100240.10592
0.111110.11575
0.11982
0.123390.126480.129190.131640.133900.136100.138320.14066
0.14383
0.14757
0.151410.155400.15960
0.163960.168510.173340.17834
0.18352
0.18894
0.19452

Fig. 6: Velocity waveform of the x-component at
the observation points on the z-axis( 0 km ≤ z ≤ 4.9km , with an interval of 100 m. This
interval is 300 m in Region II) in a 3D basin
structure, calculated (a) by the FDM using
discontinuous grid and (b) by the FDM using
uniform grid. The figures on the right are themaximum amplitude of each waveform. On theleft, the schematic of the structure employed isshown. The arrow indicates the boundary

9Fig. 7: The results of computation by the FDM using discontinuous grids (dashed line) and the FDM using
uniform grid (thin line). The waveforms in the lower part of the figure (dot line) show the differences between
them. The maximum amplitude is indicated on the right of each waveform.

The basin induced surface waves are
dominant.
In spite of the complicatedwavefield, both methods resulted in amplitude
and arrival times for all of the phases that arealmost identical, as in the 1D case. The
absence of reflected waves or disturbed waves
at the artificial boundary between Regions I and
II indicates that the wave is continuously
propagated. The fact that surface waves aretrapped by the sedimentary layer indicates theimportance of high computational accuracy inthe surface layer, which means the importanceof using small grid spacing in the surface layer.
Figure 7 shows the waveforms of the xcomponent at five observation points.
Waveforms show an excellent agreement andthe differences between the result from the
FDM with two kinds of grids are within 10 % at
all of the observation points.
In both 1D and 3D model cases, the useof the discontinuous grid enables us to savecomputation time and memory to approximately
1/4.5 of what was needed when the uniform gridis used.

Discussions
In synthesizing waveforms by the FDM
using discontinuous grids, the most important
point is the accuracy of interpolations. In
order to assist the quantitative evaluation of
the accuracy of the linear interpolation scheme
employed in this study, the interpolation of the
1D problem is considered here. One period of
the cosine function, (1 + cos 2πx / λ ) / 2f ( x) = 0(−λ / 2 ≤ x ≤ λ / 2)
( x < −λ / 2, λ / 2 < x)

(9),is discretized by the spatial interval ∆X , andeach discretized point is called X j . The cosine
function discretized with 10 grid points per
wavelength (i.e. λ / ∆X = 10 ) is shown in Figure 8(a). This function is resampled, with an intervalof ∆x (which is ∆X / 3 ） , using linear
interpolations.
The discretized points aredenoted by xi .
Figure 8 (b) shows theresampled function, and Figure 8 (c) indicatesthe difference between the discretized function

10

Fig. 8: Evaluation of error in the linear
interpolation. (a) Cosine function
discretized with 10 grid points per single
wavelength (i.e. λ / ∆X = 10 ). (b) Resampled
function with ane interval of (which is∆X / 3 ）, using linear interpolations. (c)Difference between the discretized function
and the exact value of the original function.
(d) Evaluation of differences, repeated 10
times, shifting the sampling points X j by
∆X / 10 each time. In this case, the error ofinterpolation is 0.022 (2.2 %).

FD approximation of second-order accuracy is
employed.
We have discussed the relation between
the wavelength of the cosine function and theerrors by linear interpolations. Although it isimportant to examine how the interpolation
errors propagate through the wave equation
which is sensitive to discretized partial
derivatives, it is difficult to separate the errors
due to the interpolation scheme alone fromthose caused by the grid dispersion and nonreflecting boundaries. Compared with such
errors that are thought to be small enough inthe applications, interpolation errors are of thesame order or smaller. This means that theinterpolation method we used has a sufficientlevel of accuracy.
The problem of error
propagation has to be solved in the future.
We also evaluated the ratio κ , which

and the exact value of the original function.
The maximum absolute value of this difference
divided by the maximum absolute value of thefunction, which is 1.0, is defined as theinterpolation error:
error = max f i − f ( xi ) / max f ( x)(10).
fi shows the values of the discretized function
at resampling points xi . In this case, theinterpolation error is 2.1 % (0.021). The valuesof interpolation error show some changes when
the sampling points are shifted in the space.
The evaluation of errors is repeated 10 times,
shifting the sampling points X j by X j /10 eachtime. These interpolation errors are shown in
Figure 8 (d), and as a result, in the case of
λ / ∆X = 10 , the interpolation error of theinterpolation is 2.2 %. The evaluation of theinterpolation error is performed for a casewhere λ / ∆X is 3 to 15; the outcome is shown in
Figure 9.
According to this figure, theaccuracy of the interpolation increases
monotonously when the grid spacing becomes
smaller. Based on the sampling criterion of theFD approximation of second-order accuracy,
the grid spacing should satisfy the condition
λ / ∆X > 10 (e.g. Virieux, 1986), and accordingly,
the accuracy of the interpolation is more than
2.2 %. As we explained earlier, the interpolationof each field variable is performed on ahorizontal plane, which means that we apply a
2D interpolation. However, its accuracy isalmost the same as that of the 1D interpolation,
because the interpolations in the x- and ydirections are carried out independently. This
analysis shows that a sufficient accuracy is
obtained by using the linear interpolation, in an

Fig. 9: Errors of linear interpolation. The errors ofinterpolation become smaller when the number ofgrid per single wavelength is smaller. When it isless than 10, the error of interpolation is less than
2.2 %.

11

occupies 10 % of entire region, we can reduce
thecomputational
requirements
toapproximately 1/10 and this means that the grid
spacing can be reduced to less than half
withoutincreasing
thecomputational
requirements.
The examination of the accuracy of linear
interpolation employed in this study by a simple
cosine function showed that the accuracy
increases monotonously when the grid spacing∆X decreases. It was also shown that when
the number of grid points per wavelength ismore than 10, the interpolation error is lessthan 2.2 %. Accordingly, as long as an FD
approximation of second-order accuracy is
employed, which means λ / ∆X must be morethan 10, a sufficient level of accuracy is
obtained by using the linear interpolation.
A more flexible grid system will be further
realized by increasing the number of the regions,
and/or by using the hybrid methods thatcombine discontinuous grids with variablespacing. Also the interpolation of the timesampling can reduce the number of timeincrements, rendering the computation timeever shorter.

represents the memory saving by the present
method in comparison with the constant grid
spacing FD technique. When N ⋅ ∆x × M ⋅ ∆y × L ⋅ ∆z
represents the entire region to be calculated,
and the depth of Region I is given by l ⋅ ∆z ,
κ=1 27 ( N + 6 J 0 )( M + 6 J 0 )(26l + L + 3J 0 )( N + 2 J 0 )(M + 2 J 0 )( L + J 0 )If M = N = L andJ 0 << N(11).

are assumed, thenl   10 J 0 κ = 1 27 1 + 26  × 1 +
L 
N 


(12).

(1 + 26 l L ) 27 is a term representing how muchmemory is saved by using a grid which is threetimes coarser in Region II, and 1 + 10 J 0 N is aterm representing the disadvantage resulting
from the larger absorbing region of Region I.
The typical value of l L being approximately
1/10 to 1/20, (1 + 26 l L ) 27 is approximately 1/10.
When N is 2000, 1+ 10 J 0 N is 1.1, which is aloss that can be neglected. The computation
time for interpolations is less than a few
percent of the entire computation time. Thus,a ratio similar to the case of the memory is
saved in computation time as well, that results
in an approximately tenfold saving of both
memory and time.
This signifies that the
introduction of the present method can reduce
the grid spacing to less than half without
imposing an additional load on the computers.

Acknowledgments
Discussions with Peter Moczo and Ralph
J. Archuleta provided valuable insights that
improved the paper. The reviews by Jean
Virieux and Arben Pitarka were helpful inimproving the quality and clarity of themanuscript. Yuko Kase pointed out many of
the references on FDM. Computation time wasprovided by the Supercomputer Center of the
National Research Institute for Earth Scienceand Disaster Prevention.

Conclusions
One of the biggest problems in 3D wave
propagation modeling using the FDM with auniform grid is the high computational
requirements. These requirements are related
to the oversampling of the models due to theconstant grid spacing in the region with highvelocity. In order to avoid oversamping it isnecessary to use non-uniform grids that are
adapted to the velocity structure. In this study
we used a discontinuous grid that consists of
two regions of different grid spacing. The grid
spacing ratio between the two regions was afactor of three. By connecting these tworegions by eliminating or inserting grids, we havesucceeded in dramatically reducing thecomputational
requirements,
without
significantly diminishing the accuracy of thecalculation.
In the case where Region I

References
Aoi, S., T. Iwata, K. Irikura and F. J. SanchezSesma (1995).
Waveform inversion for
determining the boundary shape of a basin
structure, Bull. Seism. Soc. Am. 85, 1445-

12

(1998). Numerical modeling of seismic wave
propagation in laterally heterogeneous
whole Earth using the pseudospectral
method, abstracts of 1998 Japan Earth andPlanetary Science Joint Meeting, 282.
Graves, R. W. (1996). Simulating seismic wave
propagation in 3D elastic media using
staggered-grid finite differences, Bull.
Seism. Soc. Am. 86, 1091-1106.
Graves, R. W. (1998).
Three-dimensional
finite-difference modeling of the San
Andreas fault: Source parameterization andground-motion Levels, Bull. Seism. Soc. Am.
88, 881-897.
Haskell, N. A. (1953).
The dispersion of
surface waves in multilayered media, Bull.
Seism. Soc. Am. 43, 17-34.
Helmberger, D. V. and J. E. Vidale (1988).
Modeling strong motions produced byearthquakes
withtwo-dimensional
numerical codes, Bull. Seism. Soc. Am. 78,
109-121.
Higdon, R. L. (1991). Absorbing boundary
conditions for elastic waves, Geophysics 56,
231-241.
Huzita, K. (1996). On survey results of active
faults in Osaka-Kobe area, Proc. 9th
seminar on studying active faults ‘on deep
structure of Osaka-bay area’, 1-11, June(in Japanese).
Iwata, T., H. Sekiguchi, A. Pitarka, K. Kamaeand K. Irikura (1998). Evaluation of strong
ground motions in the source area during
the1995
Hyogoken-Nanbu(Kobe)
earthquake, submitted to Proc. of the 10th
Japan Earthquake Engineering Symposium.
Jastram, C. and E. Tessmer (1994). Elastic
modelling on a grid with vertically varying
spacing, Geophysical Prospecting 42, 357370.
Kelly, K. R., R. W. Ward, S. Treitel, and R. M.
Alford (1976). Synthetic seismograms: afinite-difference approach, Geophysics 41,
2-27.
Kessler, D. and D. Kosloff (1991). Elastic
wave
propagationusingcylindricalcoordinates, Geophysics 56, 2080-2089.
Levander, A. R. (1988). Fourth-order finitedifference P-SV seismograms, Geophysics
53, 1425-1436.

1455.
Aoi, S., T. Iwata, H. Fujiwara, and K. Irikura,(1997). Boundary Shape Waveform inversion
for Two-Dimensional Basin Structure Using
Three-Component
Array
Data
with
Obliquely Azimuthal Plane Incident Wave,
Bull. Seism. Soc. Am. 87, 222-233
Alterman, Z. and F. C. Karal, Jr. (1968).
Propagation of elastic waves in layered
media by finite difference methods, Bull.
Seism. Soc. Am. 58, 367-398.
Boore, D. M. (1972).
Finite difference
methods for seismic wave propagation inheterogeneous materials, in Methods in
Computational Physics, Vol.11, B. A. Bolt
(Editor), Academic Press, New York.
Bouchon, M., (1981). A simple method to
calculate Green's functions for elastic
layered media, Bull. Seism. Soc. Am. 71,
959-971.
Carcione, J. M. (1992), Modeling anelastic
singular surface wave in the earth,Geophysics 57, 781-792.
Cerjan, C., D. Kosloff, R. Kosloff, and M. Reshef(1985). A nonreflecting boundary condition
for discrete acoustic and elastic wave
equations, Geophysics 50, 705-708.
Clayton, R. and B. Engquist (1977). Absorbing
boundary conditions for acoustic andelastic wave equations, Bull. Seism. Soc.
Am. 67, 1529-1540.
Falk, J., E. Tessmer, and D. Gajewski (1996).
Tube wave modeling by the finite-difference
method with varying grid spacing,
PAGEOPH 148, 77-93.
Fornberg, B. (1988).
The pseudospectral
method: Accurate representation ofinterfaces in elastic wave calculations,Geophysics 53, 625-637.
Frankel, A. (1993).
Three-dimensional
simulation of ground motion in the San
Bernardino
valley,California,
forhypothetical earthquakes on the San
Andreas fault, Bull. Seism. Soc. Am, 83,
1020-1041
Furumura, T., (1992). Studies on thepseudospectral method for the synthetic
seismograms, doctor thesis of Hokkaido
Univ.（in Japanese).
Furumura, T., B. L. N. Kennett and M. Furumura

13

McBryan, O. A., P. O. Frederickson, J. Linden,
A. Schuller, K. Solchenbach, K. Stuben, C. A.
Thole and U. Trottenberg (1991), Multigrid
methods on parallel computers ----- asurvey of recent developments, Impact
comput. sci. engrg. 3, 1-75
Mikumo, T., K. Hirahara, and T. Miyatake (1987).
Dynamical fault rupture processes inheterogeneous media, Tectonophysics 144,
19-36.
Mikumo, T. and T. Miyatake (1993). Dynamic
rupture processes on a dipping fault, andestimates of stress drop and strength
excess from the results of waveform
inversion, Geophys. J. Int. 112, 481-496.
Moczo, P. (1989). Finite-difference techniquefor SH-waves in 2-D media using irregulargrids – application to the seismic response
problem, Geophys. J. Int. 99, 321-329.
Moczo, P., P. Labák, J. Kristek, and F. Hron(1996).
Amplification and differential
motion due to an antiplane 2D resonance inthe sediment valleys embedded in a layerover the half-space, Bull. Seism. Soc. Am.
86, 1434-1446.
Nielsen, P., F. If, P. Berg, and O. Skovgaard(1994).
Using the pseudospectral
technique on curved grids for 2D acoustic
forward modelling, Geophysical Prospecting
42, 321-341.
Ohminato, T. and B. A. Chouet (1997). A
free-surface boundary condition forincluding 3D topography in the finitedifference method, Bull. Seism. Soc. Am. 87,
494-515.
Okamoto, T. (1996). Complete near-field 2.5D
finite difference seismograms for shallow
subduction zone earthquakes, abstracts ofmethod: Kobe-JMA Station Site, Japan,
case study, Geophys. Res. Lett. 23, 27292732.
Pitarka, A., K. Irikura, T. Iwata, and T. Kagawa(1996). Basin structure effects in the
Kobe area inferred from the modeling ofground motions from two aftershocks of theJanuary 17, 1995, Hyogo-ken Nanbu
Earthquake, J. Phys. Earth, 44, 563-576.
Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi(1998). Three-dimensional simulation of
the near-fault ground motion for the 1995
Hyogo-ken
Nanbu(Kobe),
Japan,
earthquake, Bull. Seism. Soc. Am., 88, 428440
Sato, T., R. W. Graves, P.G. Somerville and S.
Kataoka (1998). Estimates of regional andlocal strong motions during the Great 1923
Kanto, Japan, Earthquake (Ms 8.2). Part 2:
Forward simulation of seismograms using
variable-slip rupture models and estimationof near-fault long-period ground motions,
Bull. Seism. Soc. Am. 88, 206-227.
Schmidt, H. and G. Tango (1986). Efficient
global matrix approach to the computationof synthetic seismograms, Geophys. J. R.
astr. Soc., 84, 331-359.
Stacey, R. (1988).
Improved transparent
boundary formulations for the elastic-wave
equation, Bull. Seism. Soc. Am. 78, 20892097.
Stacey, R. (1994).
New finite-difference
methods for free surfaces with a stability
analysis, Bull. Seism. Soc. Am. 84, 171-184.
Sugawara, O., T. Yamanaka and M. Kamata(1997).
Simulationof
3D-wavepropagation in the Kanto Basin during the1923 Kanto Earthquake, abstracts of 1997

1996 Fall Meeting of the SeismologicalSociety of Japan, P09.

Japan Earth and Planetary Science Joint
Meeting, 108 (in Japanese).

Olsen B. K., Ralph J. Archuleta (1996).
Three-dimensional
simulationofearthquakes on the Los Angeles fault
system, Bull. Seism. Soc. Am, 86, 575-596.
Pitarka, A. (1999). 3D elastic finite-difference
modeling of seismic motion using staggered
grids with nonuniform spacing, Bull. Seism.
Soc. Am., 89, 54-68.
Pitarka, A. and K. Irikura (1996) Modeling 3D
surface topography by finite-difference

Vidale, J. E. and Clayton R. W. (1986). A
stable free-surface boundary condition for
two-dimensional elastic finite-difference
wave simulation, Geophysics 51, 2247-2249.
Vidale, J. E. and D. V. Helmberger (1987).
Path effects in strong motion seismology, in
Seismic Strong Motion Synthetics, (Editor B.
Bolt), Academic Press, New York.
Virieux, J. (1984). SH-wave propagation inheterogeneous
media:
Velocity-stress14

finite-difference method, Geophysics 49,
1933-1957.
Virieux, J. (1986). P-SV wave propagation inheterogeneous
media:
Velocity-stressfinite-difference method, Geophysics 51,
889-901.
Yomogida, K. and J.T. Etgen (1993). 3-D wave
propagation in the Los Angeles Basin for
the Whittier-Narrows Earthquake, Bull.
Seism. Soc. Am. 83, 1325-1344.

National Research Institute for
Earth Science and Disaster
Prevention
3-1 Tenoudai, Tsukuba-city, Ibaraki
305 Japan
aoi@geo.bosai.go.jp(S. A., H. F.)

15