Bulletin of the Seismological Society of America, Vol. 93, No. 2, pp. 817–824, April 2003

Domain Reduction Method for Three-Dimensional Earthquake Modelingin Localized Regions, Part I: Theory
by Jacobo Bielak, Kostas Loukakis, Yoshiaki Hisada, and Chiaki Yoshimura

Abstract This article reports on the development of a modular two-step, finiteelement methodology for modeling earthquake ground motion in highly heterogeneous localized regions with large contrasts in wavelengths. We target complex geological structures such as sedimentary basins and ridges that are some distance awayfrom the earthquake source. We overcome the problem of multiple physical scales
by subdividing the original problem into two simpler ones. The first is an auxiliary
problem that simulates the earthquake source and propagation path effects with amodel that encompasses the source and a background structure from which the localized feature has been removed. The second problem models local site effects. Its
input is a set of equivalent localized forces derived from the first step. These forces
act only within a single layer of elements adjacent to the interface between theexterior region and the geological feature of interest. This enables us to reduce thedomain size in the second step. If the background subsurface structure is simple, one
can replace the finite-element method in the first step with an alternative efficientmethod. The methodology is illustrated in a companion paper (Yoshimura et al.,2003) for several 3D problems of increasing physical and computational complexity.
We consider first a flat-layered, stratigraphic system. For this simple case, the firststep can be carried out by means of 3D Green’s function evaluations. The extension
to more general problems is illustrated by two examples: a basin and a hill, with thesame background stratigraphy. To verify the two-step procedure with a problem forwhich the finite-element method is used throughout, we model ground motion in asmall region of the Los Angeles Basin, using both the two-step domain-reduction
method and the traditional approach in which the computational domain contains
both the source and the geological region of interest.
Introduction
lems that involve realistic basin models with highly heterogeneous materials, because of their flexibility and simplicity.
Computer codes based on these methods have been usedsuccessfully to model earthquake ground motion in a variety
of applications; for instance, near-source ground motion(e.g., Wald and Heaton, 1994), basin structure and directivity effects (e.g., Olsen and Archuleta, 1996; Pitarka et al.,1998), and edge effects (e.g., Kawase, 1996; Hisada et al.,
1998; Aagaard et al., 2001).
Despite the recent advances in ground-motion simulation capabilities for earthquake excitation (Fig. 1), researchers nowadays are still forced to make restrictive simplifications and approximations in “3D simulations,” such as
limiting the maximum frequency or lowest wave velocities
that can be considered. One reason is that most methods
currently in use for large-size problems are based on uniform
structured grids. The grid size, which is proportional to thelowest shear wave velocity in the model and inversely pro-In the past 10 years tremendous growth has occurred in
the development of physics-based 3D models for simulating
earthquake ground motion in seismic regions. During thisperiod, numerical modeling methods for anelastic wave
propagation that take into consideration the earthquake
source, propagation path, and local site effects have become
increasingly available. There are several types of such methods. Boundary element and discrete wavenumber methods
have been popular for moderate-sized problems with relatively simple geometry and geological conditions (e.g.,
Mossessian and Dravinski, 1987; Kawase and Aki, 1990;
Hisada et al., 1993; Sánchez-Sesma and Luzón, 1995; Bouchon and Barker, 1996). Finite differences (e.g., Frankel andVidale, 1992; Frankel, 1993; Graves, 1993, 1996; Olsen etal., 1995; Pitarka, 1999; Stidham et al., 1999; Sato et al.,1999) and finite elements (e.g., Lysmer and Drake, 1971;
Toshinawa and Ohmachi, 1992; Bao, 1998; Bao et al., 1998;
Aagaard et al., 2001) are better suited for larger-sized prob817

818

J. Bielak, K. Loukakis, Y. Hisada, and C. Yoshimura

Figure 1. Schematic of semi-infinite seismic region, including causative fault, geological structure
and local features.
portional to the highest frequency of interest, is held constantthroughout the computational domain even if the softest soils
occupy only a small region. Finite-element and other methods with irregular grids (Bao et al., 1998; Pitarka, 1999; Aoi
and Fujiwara, 1999; Kristek et al., 1999; Oprsal and Zahradnı́k, 1999, 2002) are more flexible, because they can better tailor the mesh size to the local wavelength of the propagating waves. Regardless of their differences, one featurethat is common to traditional finite-difference and finiteelement methods is that the ground motions near the causative fault, and those along the propagation path and within
the region of interest, are all calculated simultaneously, using a single model that encompasses the whole geologicalstructure, from the fault region to the region of interest.
Thus, source, propagation path, and local site effects are
determined all at once. This single-pass approach works wellfor many applications. However, if the source is far from theregion of interest, the problem size becomes exceedingly
large and the methods become ineffective.
An alternative formulation, which avoids the need torepresent accurately the geometric and material properties of
the whole region within a single model, consists in subdividing the problem into two sequential parts. First, one
considers a background structure from which the localized
geological features have been deleted, and calculates the corresponding ground motion. This computation requires a gridor mesh that is only as fine as dictated by the softest material
in the background model; it needs to be performed only once
for a specified earthquake source. In a second step, only areduced region of interest which contains the localized feature is modeled to the desired accuracy. The ground motion
obtained in the first step is used to determine a set of localized equivalent forces, which are then applied as input overa computational domain that is only slightly bigger than thegeologic feature. Only the second part of the computationneeds to be repeated if any system parameters within theregion of interest need to be varied. Such a two-step procedure was developed by Bielak and his coworkers in theframework of the finite-element method, originally forbuilding–soil–foundation interaction problems (Bielak andChristiano, 1984; Cremonini et al., 1988), and later applied
to ground motion modeling of 2D sedimentary valleys in a
half-space due to incident plane SV waves (Loukakis, 1988;

Loukakis and Bielak, 1994a,b). Similar procedures havebeen presented by Clough and Penzien (1975), Kausel et al.
(1978), and Aydinoğlu (1980, 1993), under some restrictive
assumptions. Other researchers have developed alternative
two-step or hybrid procedures that make use of combinationsof different computational methods; for example, the wavenumber method and finite differences (Zahradnı́k andMoczo, 1996); modal summation method and finite differences (Regan and Harkrider, 1989; Fäh et al., 1993, 1994);
finite-element and boundary integral methods (Mita andLuco, 1987; Bielak et al., 1991); and wavenumber method,
finite differences, and finite elements introduced to represent
irregular geometries (Moczo et al., 1997). A review of various hybrid methods dating back to 1980 can be found in
Moczo et al. (1997). All these methods were concerned with
2D applications. An extension of the finite-difference approach to three dimensions has been developed recently by
Oprsal and Zahradnı́k (2002).
In this set of two articles we extend the modular twostep procedure developed by Bielak and his co-workers to
3D problems. We use a finite-element formulation in which
the primary unknowns are the total wave field within thedomain that contains the localized structure and a scatteredwave field in the exterior domain. This requires that we store
the free-field displacements from the background structure
only in a single layer of elements at the interface of the twodomains in the first step; it enables us to reduce the size ofthe computational domain in the second step. This methodology, which we call the domain reduction method (DRM),
is capable of efficiently modeling 3D wave fields for an arbitrary earthquake source in highly heterogeneous geological systems with large localized impedance contrasts andarbitrary shapes. To validate this procedure, in a companion
article (Yoshimura et al., 2003) we consider a flat-layered
system for which a solution can be obtained readily by evaluating the corresponding Green’s functions (Hisada, 1994,
1995) and illustrate its applicability to more general problems with two examples, one involving an idealized basin
and the other, a hill. To verify the two-step procedure witha problem for which the finite-element method is used inboth steps, we model ground motion in a small region of theLos Angeles Basin, using both the two-step DRM, and the
traditional approach in which the computational domain
contains both the source and the geological region of interest. Whereas the methodology is applicable to elastic, anelastic, and inelastic problems, only the elastic case is considered explicitly in this set of articles for clarity ofpresentation. In the concluding remarks, however, we discuss briefly the extension to the more general situations. It
will be seen that the use of the DRM can be advantageouseven in situations in which the causative fault is not far fromthe region of interest.

Formulation of Domain Reduction Method
We treat the problem of a semi-infinite seismic region
that contains localized geological features such as sedimen-

819

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory

tary valleys and ridges as well as seismically active faults,
as shown in Figure 1, under earthquake excitation. The geometry is arbitrary; the material is linearly elastic, and theearthquake excitation is prescribed as a kinematic source
along a predetermined fault.
Because the causative fault can be far from the geological features, we wish to define a new problem in which theexcitation is brought closer to the region of interest. Suchtransfer, of course, needs to be performed in a way that theresulting ground motion within the region of interest is identical with that due to the original source. To fix ideas, suppose the new excitation is to be specified on the fictitious
surface C shown in Figure 2a. This interface divides theseismic region into two parts: X, which contains the geological features of interest, and XⳭ, the semi-infinite exterior
subdomain, which includes the fault. We will determine theappropriate expressions for the equivalent excitation using afinite-element formulation. First, the original semi-infinite
region needs to be truncated for computational reasons. This
is indicated in Figure 2 by the inclusion of the outer boundary CⳭ. We assume, for the time being, that it is far enoughfrom the fault that no waves reflected from CⳭ reach X
within the time under consideration. This assumption will
be removed later.
Let the vector field of nodal displacements in the interior
domain X, the exterior domain XⳭ, and the boundary between them, C, be denoted, respectively, by ui (interior), ue(exterior), and ub (boundary), as shown in Figure 2a. The
seismic excitation is prescribed as a kinematic source, defined by the jump of the tangential displacements across thefault; the normal displacements and tractions remain continuous. This excitation can be equivalently specified by meansof body forces (e.g., Aki and Richards, 1980). With thisrepresentation, a standard application of Galerkin’s ideas
with finite-element spatial discretization yields a set of nodal

forces, Pe, which act in the vicinity of the fault (e.g., Bao,
1998).
Rather than analyzing simultaneously the entire domain,which includes both the fault and the localized structure, wewould like to focus on the response of a smaller region restricted to a neighborhood of the local structure. To this end,
we partition the total domain into two separate subdomains
as depicted in Figure 2b. One contains the fault and the otherthe localized geological feature. The displacements ub arecontinuous across C, and Pb are the nodal forces transmitted
by XⳭ onto X.
The ground motion within the entire computational domain is governed by Navier’s equations of elastodynamics.
When these equations are discretized spatially by finite elements in X and XⳭ, they can be expressed in partitioned
form as:

冤

MX
MX
ii
ib
X
X
Mbi Mbb

冥冦 冧 冤

üi
KX
KX
ii
ib
Ⳮ
X
X
üb
Kbi Kbb

冥 冦uu 冧 ⳱ 冦P0 冧 , in X
i

b

and

冤

Ⳮ

MX
bb
Ⳮ
MX
eb

Ⳮ

X
Mbe
XⳭ
Mee

冥 冦üü 冧
K
Ⳮ 冤
K
b
e

XⳭ
bb
XⳭ
eb

Ⳮ

X
Kbe
XⳭ
Kee

, in X冥 冦uu 冧 ⳱ 冦ⳮP
P 冧
b
e

b

Ⳮ

(2)

e

In these equations, the matrices M and K denote mass andstiffness matrices, the subscripts i, e, and b refer to nodes ineither the interior or the exterior domain or on their common
boundary, and the superscripts X and XⳭ refer to the domains over which the various matrices are defined.
The traditional form of the equation of motion for the
total domain is obtained by adding (1) and (2):

Truncated seismic region. (a) Outer boundary CⳭ restricts computations
to a finite domain; fictitious interface C divides region into two subdomains: XⳭ, which
includes the seismic source, represented by nodal forces Pe, and X, which contains the
localized geological features. (b) Regions partitioned explicitly into two substructures
across interface C; Pb are nodal forces transmitted from XⳭ onto X; ⳮPb are corresponding reactions from X onto XⳭ; nodal displacements ub are required to be continuous across C.

Figure 2.

(1)

b

820

J. Bielak, K. Loukakis, Y. Hisada, and C. Yoshimura

冤

MX
MX
0
ii
ib
Ⳮ
X
X
XⳭ
Mbi Mbb Ⳮ MX
M
bb
be
Ⳮ
Ⳮ
0
MX
MX
eb
ee

Ⳮ

冥冦 冧
üi
üb
üe

冤

KiiX
KX
0
ib
X
X
XⳭ
XⳭ
Ⳮ Kbi Kbb Ⳮ Kbb
Kbe
Ⳮ
Ⳮ
0
KX
KX
eb
ee

冥冦 冧 冦 冧ui
0ub ⳱ 0uePe(3)ue ⳱ u0e Ⳮ we

Now, to transfer the seismic excitation from the fault to C,
we consider an auxiliary problem in which the exterior region and the material therein, as well as the causative fault,
are identical with those of the original problem. The interior
domain, now denoted as X0, is, however, a simpler, background structure, that does not include the localized geological features. This is illustrated in Figure 3a. X0 is chosensuch that the new problem defined over the total domain X0
傼 XⳭ is easier to solve than the original problem. We denoteby u0i , u0b, ue0, and Pb0 the corresponding nodal displacements and the interface forces, as shown in Figure 3b. The
subscripts i, b, and e, have the same meaning as before. After
spatial discretization, the equations of motion in XⳭ for theauxiliary problem can be written as:

冤

Ⳮ

Ⳮ

冥冦 冧 冤XMX
bb MbeⳭXⳭ
Meb MXeeⳭⳭ冥冦 冧 冦 冧X
ü0b
KX
bb KbeⳭ
0 Ⳮ
XⳭüe
Keb KX
ee

ub0
ⳮPb0
0 ⳱
ue
Pe

Pe ⳱ MebXⳭü0b Ⳮ MeeXⳭ

üe0 Ⳮ KebXⳭ

u0b Ⳮ Keeu0e

Then, by substituting (5) into (3), we can solve for the
displacements ui, ub, and ue for the complete domain. This
formulation by itself, however, offers no advantage over
the traditional approach because (5) includes the terms

Figure 3.

冤

冥冦 冧

X
MX
Mib
0
ii
X
X
XⳭ
XⳭ
Mbi Mbb Ⳮ Mbb Mbe
Ⳮ
Ⳮ
0
MX
MX
eb
ee

üi
üb
ẅe

冤

KX
KX
0
ii
ib
Ⳮ
X
X
XⳭ
Ⳮ Kbi Kbb Ⳮ Kbb KX
be
Ⳮ
Ⳮ
X
0
KX
Kee
eb

(4)

(5)

(6)

That is, the residual field we is the relative displacement fieldwith respect to the reference free field u0e .
Then, substituting (6) into (3), and writing the terms that
contain the free field on the right side, results in:⳱The partitioned mass and stiffness matrices, as well as Pe,
are the same as in (2) because the material properties in XⳭ
and the earthquake source are identical in both cases.
From the second equation in (4), we can now express
the nodal forces Pe in terms of the free field, as follows:
XⳭ

Ⳮ

X
X
Mee
ü0e and Kee
u0e , which require that the free field u0e be
stored throughout the domain XⳭ. This entails an undue
computational effort.
To simplify the analysis, we use a transformation ofvariables, by which we express the total displacement ue as
the sum of the free field due to the background structure andthe residual field due to the localized geological feature:

冥冦 冧
ui
ub
we

冦

0
Ⳮ
0
XⳭ 0
ⳮMX
ü
be
e ⳮ Kbe ueⳭXXⳭ 0Pe ⳮ Meeüe0 ⳮ Kee
ue

(7)

冧

Finally, after substituting for Pe from (5) into (7), we obtain
the desired equation:

冤

冥冦 冧

X
MX
Mib
0
ii
X
X
XⳭ
XⳭ
Mbi Mbb Ⳮ Mbb
Mbe
Ⳮ
Ⳮ
0
MX
MX
eb
ee

üi
üb
ẅe

冤

KX
KX
0
ii
ib
Ⳮ
X
X
XⳭ
Ⳮ Kbi Kbb Ⳮ KX
K
bb
be
Ⳮ
XⳭ
0
KX
Kee
eb

冦

冥冦 冧
冧
ui
ub
we

0
Ⳮ
0
XⳭ 0
⳱ ⳮMX
ü
be
e ⳮ Kbe ue
Ⳮ
0
XⳭ 0
MX
eb üb Ⳮ Keb ub

Auxiliary seismic region. Localized geological features of actual problem
in X have been replaced by a simpler background structure over domain X0. (a) Entire
region; (b) region partitioned into two substructures.

(8)

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory

821

The mass matrix and stiffness matrix in the left hand side of(8) are identical with those of (3). However, the seismic
forces Pe on the fault have been replaced by the effective
nodal forces Peff, given by:

冦 冧 冦

冧

Peff
0
i
Ⳮ
eff
eff
XⳭ 0
0
P ⳱ Pb ⳱ ⳮMbe üe ⳮ KX
be ue
eff
XⳭ 0
XⳭ 0
Pe
Meb üb Ⳮ Keb ub(9)These forces have the key property that they involve only
the submatrices Mbe, Kbe, Meb, and Keb, which vanish everywhere except in a single layer of finite elements in XⳭ adjacent to C. This small domain lies between C and its adjacent surface Ce, as shown in Figure 4. Therefore, the forces
Peff act exclusively within that layer. Also, the only wave
field needed to determine Peff is that obtained from the auxiliary problem at the nodes that lie on C, Ce, and between
these surfaces. This localization of the equivalent seismic
forces around the geologic feature is the key advantage ofthe transformation (6).
Another important consequence of (9) is that all the
waves in the exterior region XⳭ will be outgoing. This suggests that for solving (8), the size of the region XⳭ can be
drastically reduced if one is interested only in the groundmotion near the localized features, provided suitable absorbing boundaries are used to limit the occurrence of spurious
waves. Because of this attractive feature, we name ourmethod Domain Reduction Method (DRM). To emphasize
this reduction in size, we will denote the reduced exteriorˆ Ⳮ and its corresponding outer boundary by Cˆ Ⳮ.
region by X
These results were derived originally in the context of a halfspace and plane wave excitation in a slightly different form
(Bielak and Christiano, 1984; Loukakis, 1988; Loukakis andBielak, 1994a). The present derivation is more rigorous andconcise, and it incorporates explicitly the effect of an extended source on a finite fault. A procedure similar to thatin Loukakis (1988) was developed subsequently by Aydinŏglu (1993), in the context of soil–structure interaction
without explicit treatment of the earthquake source. Insteadof using a finite-element formulation throughout as in Loukakis (1988) and Loukakis and Bielak (1994a), Aydinŏglu(1993) used a boundary integral representation for the tractions at the interface between the interior and exterior domains; to make the equations local at the interface, the traction was approximated in the form of a mass-dashpot-spring,and the material outside the interface C was excluded fromthe computations. It is not clear that this simplification willlead to acceptable approximations. In addition, contrary to
our formulation, in which the effective forces depend only
on the properties of the material within the region exterior
to X, the effective forces are defined on a layer within theinterior region. This could pose some difficulty if one isinterested in considering interior regions that behave nonlinearly, as discussed in the next section. Similar results were
derived also by Zahradnı́k and Moczo (1996) for the finite
difference method in two dimensions using a rectangular

Figure 4. Seismic region with two neighboring
surfaces C and Ce on which effective nodal forces Peffdefined by equation (9) are to be applied. These forces
are equivalent to and replace the original seismic
forces Pe, which act in the vicinity of the causative
fault.

excitation box, and an algorithm similar to that of Alterman
and Karal (1968) to do the coupling.

Discussion and Concluding Remarks
The results described in the previous paragraphs can besummarized as a two-step procedure for analyzing the earthquake response of localized geological features, as follows.
In step I, as shown in Figure 5a, one starts with a backgroundgeological model that is defined over the domains XⳭ and
X0, which include the original earthquake source, and thatdefines the boundary C of what will be the region of interestin step II. Then one calculates the free-field ground motion
u0b and u0e and stores it at all the nodes on the adjacent surfaces C and Ce as well as at any interior nodes within thefinite element layer that lies between them. Suitable absorbing boundary conditions on CⳭ, such as those proposed by
Clayton and Engquist (1977) and Stacey (1988), must be
used to keep spurious wave reflections within acceptable
limits. Spurious reflections from the absorbing boundaries
are generally unavoidable in practice, and can lead to inaccuracies in the numerical results.
If desired, any appropriate method can be used in step
I to determine these wave fields in place of the finite-element
method. The second step is performed on a reduced region
X 傼 X̂Ⳮ, which contains the geological features of interest,
but not the causative fault, as shown in Figure 5b. The effective seismic forces Peff are evaluated first, from (9), using
u0b and u0e as input. Once these forces have been established,
the total wave fields ui and ub, and the residual wave field
we, defined respectively over X, C, and X̂Ⳮ, are obtained bysolving (8). Actually, this equation must be modified slightly
to incorporate the absorbing boundary conditions on ĈⳭ.
The complete solution ue within the domain X̂Ⳮ can be evaluated simply as ue ⳱ we Ⳮ ue0.
It is important to emphasize that the DRM is exact towithin finite element spatial and time discretization errors.
In our implementation we have used piecewise linear finite
elements in space, and second-order central differences intime. Thus, our results are second-order accurate in space

822

J. Bielak, K. Loukakis, Y. Hisada, and C. Yoshimura

Figure 5. Summary of two-step DRM. (a) Step I
defines the auxiliary problem over background geological model. Resulting nodal displacements within
C, Ce, and the region between them are used to evaluate effective seismic forces Peff required for step II.
(b) Step II, defined over reduced region made up of
X and X̂Ⳮ (a truncated portion of XⳭ). The effective
seismic forces Peff are applied within C and Ce. The
unknowns are the total displacement fields ui in X andub on C, and the residual displacements we in X̂Ⳮ.
and time. The main computational efficiency is gained fromthe use of the finite-element method, which allows one totailor the mesh size to the local wavelength of the propagating wave, and to the reduction of the overall computational domain in the second step, which might be quite dramatic if the source is far, the material within the region ofinterest is nonlinear, or repeated solutions with minorchanges in the parameters are required, as in inverse problems. With the goal of increasing the computational efficiency of 3D simulations, other authors have developed
powerful memory optimization methods, especially in thecontext of the finite-difference method, to be able to deal
with more realistic material properties and frequencies closer
to those of engineering interest. In particular, Moczo et al.
(1999) use a lossy compression scheme based on a discrete
wavelet transform to optimize RAM and disk storage requirements. Combining the compression scheme with thememory optimization of Graves (1996) and the memory
variable economization of Day (1998) and Day and Bradley(2001), Moczo et al. (2000, 2001) in some cases achieve a
total disk plus RAM storage reduction approaching or exceeding an order of magnitude. Such optimization tech-niques could be applied also to the DRM to attain further
storage reductions.
It should also be emphasized that the strict validity ofour two-step methodology hinges on the linearity of the material outside the region of interest X and on the requirementthat the material properties in the exterior region be the samefor the auxiliary problem as for the original problem. The
material within X, on the other hand, can be arbitrary. Noticethat its properties do not enter into the computation of Peff.
In fact, the two-step procedure we just described would remain valid even if the material in X were nonlinear. All thatwould change in the formulation are the stiffness terms in(3) and (7) pertaining to the interior domain, which would
then depend on the solution. More generally, the second termin (1) and in subsequent equations would be a nonlinear
function of ui and ub, which would depend on the material
constitutive equations. Another assumption in our derivationis that the material is nondissipative. Clearly, adding linear
dissipation would not change the essence of the procedure,
but merely modify the particular form of the equations. For
instance, the addition of viscous damping would result inadded terms proportional to velocity in the equations (Loukakis, 1988; Loukakis and Bielak, 1994b). More generally,
linear viscoelastic behavior would, rigorously, introduce
convolutions into the semidiscretized formulation. Such
convolutions, however, can be avoided if one approximates
the viscoelastic dependence by a linear combination of multiple relaxation mechanisms that can be expressed in differential form with the aid of auxiliary memory variables (e.g.,
Day and Minster, 1984; Emmerich and Korn, 1987; Carcione et al. 1988; Moczo et al., 1997). One important drawback encountered originally with these methods was theadded computational and storage cost associated with theauxiliary memory variables. This difficulty has been partially removed by the coarse-graining methodology developed by Day (1998), in which only one individual relation
mechanism per node produces highly accurate results.

Acknowledgments
This work was supported by the National Science Foundation, under
grant CMS-9980063 of the KDI program, and under grant EEC-0121989
of the Engineering Research Centers. The cognizant program directors are
Clifford J. Astill and Joy M. Pauschke, respectively. This work was also
partly supported by a special project of U.S.–Japan Cooperative Researchfor Urban Earthquake Disaster Mitigation (11209201), Grant-in-Aid for
Scientific Research on Priority Area (Category B), Ministry of Education,
Culture, Sports, Science, and Technology of Japan (PI: Tomotaka Iwata),
and by Taisei Corporation. We are grateful for this support. We thank Peter
Moczo and an anonymous reviewer, whose reviews helped us improve themanuscript. The work was done while one of the authors (C. Y.) was visiting Carnegie Mellon University as a research scholar on leave from Taisei
Corporation.

References
Aagaard, B. T., J. F. Hall, and T. Heaton (2001). Characterization of nearsource ground motions with earthquake simulations, Earthquake
Spectra 17, 177–207.

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory
Aki, K., and P. G. Richards (1980). Quantitative Seismology: Theory andMethods, W. H. Freeman, New York.
Alterman, Z., and F. C. Karal, Jr. (1968). Propagation of elastic waves inlayered media by finite difference methods, Bull. Seism. Soc. Am. 58,
367–398.
Aoi, S., and H. Fujiwara (1999). 3-D finite difference method using discontinuous grids, Bull. Seism. Soc. Am. 89, 918–930.
Aydinoğlu, M. N. (1980). Unified formulations for soil-structure interaction, in Proc. 7th World Conf. Earthquake Engineering, Istanbul, Turkey, 8–13 September 1980, Vol. 7, 121–128.
Aydinoğlu, M. N. (1993). Consistent formulation of direct and substructure
methods in nonlinear soil-structure interaction, Soil Dyn. Earthquake
Eng. 12, 403–410.
Bao, H. (1998). Finite element simulation of earthquake ground motion inrealistic basins, Ph.D. Thesis, Carnegie Mellon University, Pittsburgh.
Bao, H., J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O’Hallaron, J. R.
Shewchuk, and J. Xu (1998). Large-scale simulation of elastic wave
propagation in heterogeneous media on parallel computers, Comput.
Methods Appl. Mech. Eng. 152, 85–102.
Bielak, J., and P. Christiano (1984). On the effective seismic input for nonlinear soil-structure interaction systems, Earthquake Eng. Struct. Dyn.
12, 107–119.
Bielak, J., R. C. MacCamy, D. S. McGhee, and A. Barry (1991). Unified
symmetric BEM-FEM for site effects on ground motion—SH waves,
J. Eng. Mech. 117, 2265–2285.
Bouchon, M., and J. S. Barker (1996). Seismic response of a hill: The
example of Tarzana, California, Bull. Seism. Soc. Am. 86, 66–72.
Carcione, J. M., D. Kosloff, and R. Kosloff (1988). Wave propagation ina linear viscoacoustic medium, J. R. Astr. Soc. 78, 105–118.
Clayton, R. W., and B. Engquist (1977). Absorbing boundary conditions
for acoustic and elastic wave equations, Bull. Seism. Soc. Am. 67,
1529–1540.
Clough, R. W., and J. Penzien (1975). Dynamics of Structures, McGrawHill, New York, 584.
Cremonini, M. G., P. Christiano, and J. Bielak (1988). Implementation ofeffective seismic input for soil-structure interaction systems, Earthquake Eng. Struct. Dyn. 16, 615–625.
Day, S. M. (1988). Efficient simulation of constant Q using coarse-grained
memory variables, Bull. Seism. Soc. Am. 88, 105–1062.
Day, S. M., and C. R. Bradley (2001). Memory-efficient simulation of
anealstic wave propagation, Bull. Seism. Soc. Am. 91, 520–531.
Day, S. M., and J. B. Minster (1984). Numerical simulation of attenuatedwavefields using a Padé approximant method, Geophys. J. R. Astr.
Soc. 78, 105–118.
Emmerich, H., and M. Korn (1987). Incorporation of attenuation into time
domain computations of seismic wave fields, Geophysics 52, 1252–
1264.
Fäh, D., P. Suhadolc, and G. F. Panza (1993). Variability of seismic groundmotion in complex media: the case of a sedimentary basin in the Friuli(Italy) area, J. Appl. Geophys. 30, 131–148.
Fäh, D., P. Suhadolc, St. Mueller, and G. F. Panza (1994). A hybrid method
for the estimation of ground motion in sedimentary basins: Quantitative modeling for Mexico City, Bull. Seism. Soc. Am. 84, 383–399.
Frankel, A., and J. E. Vidale (1992). A three-dimensional simulation ofseismic waves in the Santa Clara Valley, California from a Loma
Prieta aftershock, Bull. Seism. Soc. Am. 82, 2045–2074.
Frankel, A. (1993). Three-dimensional simulations of ground motions in
the San Bernardino Valley, California, for hypothetical earthquakes
on the San Andreas fault, Bull. Seism. Soc. Am. 83, 1024–1041.
Graves, R. W. (1993). Modeling three-dimensional site response effects in
the Marina District Basin, San Francisco, California, Bull. Seism. Soc.
Am. 83, 1042–1063.
Graves, R. W. (1996). Simulating seismic wave propagation in 3D elastic
media using staggered-grid finte-differences, Bull. Seism. Soc. Am.
86, 1091–1106.
Hisada, Y. (1994). An efficient method for computing Green’s function for823

a layered half-space with sources and receivers at close depths, Bull.
Seism. Soc. Am. 84, 1456–1472.
Hisada, Y. (1995). An efficient method for computing Green’s function fora layered half-space with sources and receivers at close depths (Part2), Bull. Seism. Soc. Am. 85, 1080–1093.
Hisada, Y., K. Aki, and T. L. Teng (1993). 3-D simulations of surface wave
propagation in the Kanto sedimentary basin, Japan. Part 2: Application of the surface wave BEM, Bull. Seism. Soc. Am. 83, 1700–1720.
Hisada, Y., H. Bao, J. Bielak, O. Gattas, and D. R. O’Hallaron (1998).
Simulation of long-period ground motions during the 1995 Hyogoken-Nanbu (Kobe) earthquake using 3-D finite element method, inProc. 2nd International Symposium on the Effect of Surface Geology
on Seismic Motion, 1–3 December 1998, Yokohama, Japan, 1353–
1360.
Kausel, E., R. V. Whitman, J. P. Morray, and F. Elsabee (1978). The spring
method for embedded foundations, Nucl. Eng. Des. 48, 377–392.
Kawase, H. (1996). The cause of the damage belt in Kobe: “the basin-edge
effect”, constructive interference of the direct S-wave with the basininduced diffracted/Rayleigh wave, Seism. Res. Lett. 67, 25–34.
Kawase, H., and K. Aki (1990). Topography effect at the critical SV-wave
incidence: Possible explanation of damage pattern by the Whittier
Narrows, California, earthquake of 1 October 1987, Bull. Seism. Soc.
Am. 80, 1–22.
Kristek, J., P. Moczo, K. Irikura, T. Iwata, and H. Sekiguchi (1999). The
1995 Kobe mainshock simulated by the 3D finite differences, in The
Effects of Surface Geology on Seismic Motion, K. Irikura, et al. (Editors), Vol. 3, Balkema, Rotterdam, 1361–1368.
Loukakis, K. E. (1988). Transient response of shallow layered valleys forinclined incident SV waves calculated by the finite element method,
M.Sc. Thesis, Carnegie Mellon University, Pittsburgh.
Loukakis, K., and J. Bielak (1994a). Seismic response of two-dimensional
sediment-filled valleys to oblique incident SV-waves calculated by
the finite element method, in Proceeding of the Fifth U.S. National
Conference on Earthquake Engineering, 10–14 July 1994, Chicago,
Illinois, Vol. III, 25–34.
Loukakis, K., and J. Bielak (1994b). Layering and damping effects on seismic response of sedimentary valleys to oblique excitation, in Proceedings of the Second International Conference on Earthquake Resistant Construction and Design, 15–17 June 1994, Berlin, Germany,
93–100.
Lysmer, J., and L. A. Drake (1971). The propagation of Love waves acrossnonhorizontally layered structures, Bull. Seism. Soc. Am. 61, 1233–
1252.
Mita, A., and J. E. Luco (1987). Dynamic response of embedded foundations: a hybrid approach, Comput. Methods Appl. Mech. Eng. 63, 233–
259.
Moczo, P., E. Bystrický, J. Kristek, J. M. Carcione, and M. Bouchon (1997).
Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures, Bull. Seism. Soc. Am. 87, 1305–1323.
Moczo, P., J. Kristek, and E. Bystrické (2001). Efficiency and optimizationof the 3D finite-difference modeling of seismic ground motion, J.
Comput. Acoust. 9, 593–609.
Moczo, P., M. Lucká, J. Kristek, and M. Kristekové (1999). 3D displacement finite differences and a combined memory optimization, Bull.
Seism. Soc. Am. 89, 69–79.
Mossessian, T., and M. Dravinski (1987). Application of a hybrid method
for scattering of P, SV, and Rayleigh waves by near-surface irregularities, Bull. Seism. Soc. Am. 77, 1784–1803.
Olsen K. B., J. C. Pechmann, and G. T. Schuster (1995). Simulation of 3D
elastic wave propagation in the Salt Lake basin, Bull. Seism. Soc. Am.
85, 1688–1710.
Olsen K. B., and R. J. Archuleta (1996). 3D simulation of earthquakes on
the Los Angeles fault system, Bull. Seism. Soc. Am., 86, 575–596.
Oprsal, I., and J. Zahradnı́k (1999). Elastic finite-difference method for
irregular grids, Geophysics 64, 240–250.
Oprsal, I., and J. Zahradnı́k (2002). Three-dimensional finite difference

824
method and hybrid modeling of earthquake ground motion, J. Geophys. Res. 107 (B8), 2161, doi:10.1029/2000JB000082.
Pitarka A. (1999). 3D elastic finite-difference modeling of seismic motion
using staggered grids with non-uniform spacing, Bull. Seism. Soc. Am.
89, 54–68.
Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi (1998). Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogoken Nanbu (Kobe), Japan, earthquake, Bull. Seism. Soc. Am. 88, 428–
440.
Regan, J., and D. G. Harkrider (1989). Numerical modelling of SH Lg waves
in and near continental margins, Geophys. J. Int. 98, 107–130.
Sánchez-Sesma, F. J., and F. Luzón (1995). Seismic response of threedimensional valleys for incident P, S, and Rayleigh waves, Bull.
Seism. Soc. Am. 85, 269–284.
Sato, T., R. W. Graves, and P. G. Somerville (1999). Three-dimensional
finite-difference simulations of long-period strong motions in the Tokyo metropolitan area during the 1990 Odawara earthquake (MJ 5.1)
and the great 1923 Kanto earthquake (MS 8.2) in Japan, Bull. Seism.
Soc. Am. 89, 579–607.
Stacey, R. (1988). Improved transparent boundary formulations for theelastic-wave equation, Short Note, Bull. Seism. Soc. Am. 78, 2089–
2097.
Stidham, C., M. Antolik, D. Dreger, S. Larsen, and B. Romanowicz (1999).
Three-dimensional structure influences on the strong motion wavefield of the 1989 Loma Prieta earthquake, Bull. Seism. Soc. Am. 89,
1184–1202.
Toshinawa, T., and T. Ohmachi (1992). Love wave propagation in a threedimensional sedimentary basin, Bull. Seism. Soc. Am. 82, 1661–1667.
Wald, D. J., and T. H. Heaton (1994). Spatial and temporal distribution ofslip for the 1992 Landers, California, earthquake, Bull. Seism. Soc.
Am. 84, 668–691.
Yoshimura, C., J. Bielak, Y. Hisada, and A. Fernández (2002). Domain
reduction method for three-dimensional earthquake modeling in localized regions, Part II: Verification and applications, Bull. Seism.
Soc. Am. 93, no. 2, 825–840.

J. Bielak, K. Loukakis, Y. Hisada, and C. Yoshimura
Zahradnı́k, J. (1995). Simple elastic finite-difference scheme, Bull. Seism.
Soc. Am. 85, 1879–1887.
Zahradnı́k, J., and P. Moczo (1996). Hybrid seismic modeling based on
discrete-wavenumber and finite-difference methods, Pure Appl. Geophys. 148, 21–38.
Department of Civil and Environmental Engineering
Carnegie Mellon University
Pittsburgh, PA 15213-3890
jbielak@cmu.edu(J.B.)
Technical Company of General Construction
Isiodou 22, Athens 10674, Greece
klukakis@compulink.gr(K.L.)
Department of Architecture
Kogakuin University
1-24-2 Nishi-shinjuku, Shinjuku-ku
Tokyo 163-8677, Japan
hisada@cc.kogakuin.ac.jp(Y.H.)
Technology Center
Taisei Corporation
344-1 Nase-cho, Totsuka-ku
Yokohama 245-0051, Japanyosimura@eng.taisei.co.jp(C.Y.)
Manuscript received 21 September 2001.