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

Domain Reduction Method for Three-Dimensional Earthquake Modelingin Localized Regions, Part II: Verification and Applications
by Chiaki Yoshimura,* Jacobo Bielak, Yoshiaki Hisada, and Antonio Fernández

Abstract Several examples are used to verify the domain reduction method(DRM), a two-step finite-element methodology described in a companion article formodeling earthquake ground motion in highly heterogeneous three-dimensional localized regions. The first set involves a simple, flat-layered system. Verification of
the DRM for this problem is carried out by comparing the results to those calculated
directly by the theoretical Green’s function method. Its applicability to more general
problems is illustrated by two examples: a basin and a hill with and without a weathered surface layer and with the same stratigraphy. We use a Green’s function approach for the first step, which for the examples under consideration needs to be
performed only once. For the second step, the domain of computation is reduced ineach case to a small neighborhood of the geological feature at hand. The second
application considers the ground motion due to a strike-slip double couple buried
14 km below the free surface in an 80-km ⳯ 80-km ⳯ 30-km region that encloses
entirely the Los Angeles basin. This problem is solved first by the finite-element
method using the single-step traditional approach, in which the ground motion is
calculated simultaneously near the seismic source, along the propagation path, andwithin the region of interest with a single model that encompasses the entire geological structure, from the source to the region of interest. The DRM is then used todetermine anew the ground motion over a much smaller (6-km ⳯ 6-km ⳯ 0.6-km)
region contained within the original domain, and the results of the two methods
within this region of interest are compared.
These examples serve to demonstrate that in many applications the DRM can besignificantly more efficient than the traditional approach. The DRM can be particularly advantageous (1) if the source is far from the local structure and the local
structure is much softer than that of the exterior region, (2) if the localized feature
exhibits nonlinear behavior, or (3) if for a prescribed source, one wishes to consider
a sequence of simulations in which the properties of the local feature, which might
include man-made structures, are varied from one simulation to the next.
Introduction
In a companion article (Bielak et al., 2003), hereafterreferred to as Paper I, we have described a two-step, domain
reduction method (DRM) for modeling efficiently source,
path, and site effects during earthquakes. A particular caseof this method, in which the seismic excitation consisted ofa plane-wave incident wave and the examples were restricted
to two dimensions, was presented earlier (e.g., Loukakis,
1988; Loukakis and Bielak, 1995). Alternative two-stepmethods that share the attractive features of the DRM havebeen presented also by other authors (e.g., Aydinoğlu, 1993;

Zahradnı́k and Moczo, 1996; Moczo et al., 1997). A moredetailed bibliography is included in Paper I. Here we assess
the DRM by comparing our results with those from established methods, for two particular three-dimensional problems. We also illustrate the applicability to other problems
and discuss extensions and limitations of the method.
The DRM procedure and definition of regions, boundaries, and variables described in Paper I are summarized in
Figure 1. In step 1 (Fig. 1a), one stores the free-field displacement u0b and u0e within two boundaries C and Ce for a
simpler auxiliary problem without local features (X0 傼 XⳭbounded by CⳭ). These displacements are used for calculateff
ing effective seismic forces Peffb and Pe in step 2 (Fig. 1b).

*Present address: Technology Center, Taisei Corporation, 344-1 Nasecho, Totsuka-ku, Yokohama 245-0051, Japan; yosimura@eng.taisei.co.jp.

825

826

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

These forces are rigorously equivalent within discretization
error to the fault fource Pe. Since they are distributed only
in a layer around X, which contains the local feature, thedomain size for step 2 can be reduced to a smaller region ofinterest. In step 2, the total wave field ui, ub and the residual
wave field we (we ⳱ ue ⳮ u0e , where ue is the total displacement in X̂Ⳮ) for a complex problem with local features (Xeffˆ Ⳮ bounded by Cˆ Ⳮ) are calculated using Peff
傼X
b and Pe .
In this article we consider, as a first particular instance,a flat-layered system for the background structure in step 1(X0 傼 XⳭ). Free-field displacements u0b and u0e are calculated by the Green’s function method (Hisada, 1994, 1995).
Verification of the DRM is carried out by taking the material
in X to be identical to that of the background structure X0
and by comparing the results to those calculated directly by
the Green’s function method. Applicability of the DRM isdemonstrated by replacing X in step 2 with local structures,
including a basin or a hill. Subsequently, we consider an
additional problem involving the Los Angeles basin to illustrate and verify the applicability of the DRM for more

realistic situations, which entail a more complex geometry
and highly heterogeneous material properties.

Flat-Layered System
Model Verification
We consider the two-layer system underlain by an elastic half-space shown in Figure 2; its properties are listed inTable 1. These can, of course, be readily scaled. No material
damping is considered, for simplicity. The seismic source isa dip-slip double couple buried at a point 1 km below thefree surface. The strike, dip, and rake are 0⬚, 90⬚, and 90⬚,
respectively, the seismic moment M0 ⳱ 6 ⳯ 1015 N m, andthe slip function is shown in Figure 3. North is aligned with

Figure 2.

Flat-layered system used to verify the

DRM, with seismic source and region of interest(ROI).

Table 1
Soil Parameters of Layered System
Layer

1
2
3

Thickness(m)Vs(m/sec)Vp(m/sec)Density(g/cm3)

200
400
⬁

250
500
1000

500
1000
2000

2.0
2.0
2.0

Figure 1.

Summary of two-step domain reduction
method (DRM). (a) Step 1 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 2. (b) Step 2, defined overreduced region made up of X and X̂Ⳮ (a truncatedportion of XⳭ). The effective seismic forces Peff areapplied within C and Ce. The unknowns are the total
displacement fields ui in X and ub on C and the residual displacements we in X̂Ⳮ.

Figure 3. Slip function used for seismic source(double couple applied within elastic half-space at
1 km beneath the free surface).

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

the y axis. The region of interest is a 1-km ⳯ 1-km ⳯
1-km cube located 10 km east of the epicenter.
Because of the simplicity of the physical setting and of
the source, in this example we evaluate the displacements
u0b and u0e on C and Ce in step 1 using the three-dimensional
Green’s functions in the computer code by Hisada (1994,
1995). In step 2 the material in X is taken to be identical tothat of the background material, as shown in Figure 4. We
use an elastic wave propagation finite-element simulation
code developed for modeling earthquake ground motion inlarge sedimentary basins (Bao et al., 1998). The wave propagation code is built on top of Archimedes, an environment
for solving unstructured-mesh finite-element problems on
parallel computers (Bao et al., 1998; Archimedes, 1998).
Archimedes includes two- and three-dimensional mesh generators, a mesh partitioner, a parceler, and a parallel code
generator. We use linear tetrahedral elements. We haveadded the capability of automatically determining the surfaces C and Ce once a box that defines the region of interest
has been prescribed and have built in the necessary operations for evaluating the effective forces Peff. These calculations are also performed in parallel since the contribution to
Peff from each element within the layer adjacent to X can beevaluated independently of the other elements within thatlayer. In addition, in this study we have used a lumped mass
matrix approach, in which one-fourth of the mass of eachtetrahedron is assigned to each node. Therefore, the offdiagonal elements of the mass matrix vanish, and consequently, the evaluation of the effective seismic forces involves only a multiplication of the stiffness matrix by afree-field displacement (see equation 8 in Paper I). We use
lumped mass matrices to avoid the need of solving a systemof algebraic equations at each timestep. With this choice, theonly significant algebraic operation at each timestep is amatrix–vector multiplication. In addition, since linear elements have nodes only at the vertices, no displacements need
to be stored inside the layer between C and Ce.
The solution from the DRM for points on a fictitious
borehole that passes through points B and B⬘ (Fig. 4b) isshown in Figure 5. Figure 5a depicts the x (east–west) component of the displacement, and Figure 5b the vertical component, at various depths. (Displacements in the y [north–south] direction are not shown as they essentially vanish,
due to symmetry.) The complete wave field, including bodywaves and surface waves, can be clearly observed in thisfigure. The corresponding results from the Green’s functions
evaluations are also shown in Figure 5, for comparison. Peakvalues, with their signs, are listed on the right columns forboth solutions next to the synthetic seismograms. The agreement between the two sets of waveforms is quite good, withmaximum differences in amplitude on the order of 5%. Thisis consistent with the accuracy we can expect from our finiteelement approximation, which is tailored to 10 points per
wavelength, according to the shear-wave velocity withineach element and a maximum frequency of 1 Hz. The approximate dominant frequency of the surface waves is 0.3

827

Figure 4. Layered system within region of interest. (a) Finite-element mesh tailored to shear-wave
velocity of each layer and the half-space; (b) crosssection on vertical plane through AA⬘. The bolddashed lines show surfaces C and Ce where effective
forces Peff are applied in step 2.

Hz. There are a total of 102,402 elements and 19,143 nodes
in the mesh shown in Figure 4. It took 12 min of CPU timeon eight processors of the T3E computer at the Pittsburgh
Supercomputing Center to solve the corresponding equationof motion, using a second-order central difference methodwith a timestep of 0.01 sec.
Notice that the agreement between the finite-element
solutions and the corresponding Green’s functions remains
quite close down to the interface C, at 700 m. Right below
this point, the finite-element solution almost vanishes. The
same behavior is observed in Figure 6 for the seismograms

828

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 5.

Synthetic seismograms for displacements along downhole line BB⬘ (Fig. 4b). The depth
from free surface and shear-wave velocity of each material is indicated to left of seismograms. (Other properties are listed in Table 1.) The scale, in centimeters,
is shown above the origin of the first seismogram.
Peak displacements from finite-element DRM simulations and corresponding values from Green’s functions calculations are shown to the right of the seismograms. (a) x- (east–west) component; (b) z- (up–down) component.

on the free surface along AA⬘ (Fig. 4a). The difference between the results from the DRM approach and the Green’s
functions does not exceed 5% at these locations, and thedisplacements beyond C also essentially vanish. Recall thatin the outer region X̂Ⳮ, our formulation yields residual displacements; since the material in X is the same as that in the
background structure for this example, we must vanish. The

Figure 6.

Synthetic seismograms for displacements along free-surface (horizontal) line AA⬘ (Fig.
4a). The distance x is measured from the origin of thex axis. Other nomenclature is as in Figure 5.

fact that the numerical values of these residual displacements
are close to zero provides a useful numerical check. An interesting consequence of the vanishing of we for this problem
is that, theoretically, the outer boundary ĈⳭ must play norole in the solution, regardless of the absorbing boundary
conditions one uses there. For the present application theboundary nodes were left unconstrained, thereby implyingthat the outer boundary is traction free. The fact that residual
displacements in X̂Ⳮ are barely visible confirms that for thevalidation problem the boundary condition on ĈⳭ has an
insignificant numerical effect. Moreover, since there are nowaves leaving the region of interest X, one could modify
the material in the exterior region beyond a single-element

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

829

thick layer surrounding the surface Ce, and the results within
X would not change.
Flat-Layered System with Basin and Hill
Idealized Basin. The first example we use to illustrate theapplicability of DRM to more complex situations is one thatinvolves a local structure X with a sedimentary basin embedded into the same two-layer stratigraphic system we considered in the previous section. The basin is in the shape ofa spherical cap and has a maximum depth of 100 m and a
150-m radius at its intersection with the free surface, as
shown in Figure 7. It has a uniform shear-wave velocity of
125 m/sec, P-wave velocity of 250 m/sec, and density of
2 gm/cm3. The seismic source is identical to that for theunperturbed flat-layered system. Thus, there is no need torecalculate the free-field ground motion u0b and u0e on C andCe, as we can reuse the seismograms obtained previously.
On the other hand, in contrast to the flat-layered system forwhich the residual displacement vanished outside the regionof interest, in this problem the basin generates a scatteredwave. Hence, an absorbing boundary must be introduced on
ĈⳭ to limit the occurrence of spurious reflections. We use
a simple dashpot approach (Lysmer and Kuhlemeyer, 1969)
for this purpose, which consists of adding viscous dampers
at each node on ĈⳭ. This gives rise to a diagonal damping
matrix C with nonzero terms associated only with boundary
nodes.
The resulting displacement synthetics along the lineBB⬘ (Fig. 7) are shown in Figure 8, together with the corresponding values for the background (flat-layered) structure. As expected, the basin has the effect of magnifying theamplitude of the free-field ground motion. This amplificationis confined primarily to points within the basin, where it
reaches values of about 50% in the x (east–west) direction.
The vertical amplification is only of the order of 20%.
Figure 9 shows synthetics for locations along AA⬘ (Fig.
7) for the x (east–west) component of the displacement. The
corresponding results for the vertical component are in Figure 10. The top panel in each figure depicts the free-field
ground motion in the background structure with no valley.
There is a discontinuity across the interface C because within
the region of interest X we plot the total displacement, but
in the exterior region only the residual displacement isshown. The middle panel shows the corresponding resultswith the basin present, and the bottom panel the difference
between the previous two. The latter is the residual displacement along the line AA⬘ over the entire region of interest.
In this case, since the stratigraphy in steps 1 and 2 is thesame, the residual wave field corresponds identically to thescattered wave field. The effect of the basin on the groundmotion can be seen explicitly from the bottom panel or bycomparing the top two panels in each figure. The basin effects in the x (east–west) direction are of the same order inthe basin’s interior as the free-field motion, especially in the
deepest part of the basin. This region is affected most be-

Figure 7. Homogeneous basin embedded in flatlayered system. (a) Finite-element mesh; (b) crosssection through AA⬘.

cause the prescribed ground motion excites primarily thebasin’s fundamental mode. Effects for the vertical groundmotion are less pronounced. The residual motion is continuous across C, as expected, because this fictitious interface
has no physical meaning and has been introduced merely forcomputational convenience. It is clear from Figure 9b and c
and Figure 10b and c that the wave motion outside the basin
is purely outgoing and that no visible spurious reflections
are being generated at the absorbing boundary. The peak
amplitude of the scattered wave field at the edges of thecomputational domain (x ⳱ 0, 1000 m) is of the order of
10% of the peak amplitude of the background free-field motion with no valley.

830

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 9. Synthetic seismograms for the east–west component of displacement along free-surface
horizontal line AA⬘ (Fig. 7a). Panels (a) and (b) show
total displacements inside and on C and residual displacements outside of it; panel (c) shows residual displacement at all locations along AA⬘ and thus depicts
directly the effect of the basin on surface groundmotion.
Figure 8.

Synthetic seismograms for displacements along downhole line BB⬘ (Fig. 7b). The solidlines show the response with basin present. The
dashed lines correspond to free-field motion (without
the valley). The right-hand columns show peak values
with and without basin. Traces for points within surface C represent total displacement; those for points
outside this surface show residual displacements withrespect to free-field surface motion of the corresponding points for the flat-layered system. (a) East–west
component; (b) up–down component.

To summarize the results from the surface response of
the basin and its immediate vicinity to the incoming seismic
waves, we have plotted in Figure 11 the distribution of the
maximum value of the total surface displacements overthe entire region of interest. Here we have included also the
response in the y (north–south) direction. Notice the differ-ent scales used for each component; also, the smallest andlargest values of these maxima after spatial smoothing arereported in each panel. The effect of the basin on the freesurface ground can be clearly seen for the various components of the wave field. The peak amplification occurs off
center, especially for the vertical motion. In addition, there
are noticeable backward and forward scattering effects in thevicinity of the basin. The former leads to an increase in the
ground motion with respect to that of the free-field motion,
while the latter has the opposite effect (Fig. 11a,c).
Idealized Hill. To illustrate the applicability of our procedure to the analysis of topographic effects from exposed
geological structures, we consider, as a second example, thecase of a hill supported on the two-layered system, as shownin Figure 12. We model the ground motion for two variationsof the hill problem. In the first instance, the hill is assumed

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

Figure 10. Same as Figure 9, but for the up–down
component of displacement.

to be homogeneous with the same properties as the top layerof the background material; in the second, the hill has aweathered surface layer 25 m thick, with the same properties
as those of the basin in the previous example. The hill has
a square base 350 m ⳯ 350 m, it is 100 m high, and thelateral sides have a slope of 45⬚. The seismic excitation is
the same as before, as is the numerical procedure in almost
every detail. The one difference here is that the localized
topographical feature is located above the free surface of thelayered system, and its lateral and top surfaces are traction
free. It is worth noting that contrary to other methods suchas finite differences, no treatment of any kind is required in
the finite-element method to enforce traction boundary orinterface conditions. These are satisfied automatically as aconsequence of the variational principle that underlies thefinite-element formulation.
Seismograms for several locations along the line AA⬘
on the free surface in Figure 12 are shown in Figure 13 forthe homogeneous hill. The x (east–west) and z (up–down)
components of the displacement are depicted, together witha list of their peak values. A comparison of these results with831

the corresponding free-field displacements in Figures 9a and
10a reveals an amplification on the order of 2.5 at the crestof the hill in the maximum amplitude of the east–west component of displacement. For the vertical component this amplification is only about 1.5. This topographic amplification
can be clearly observed in Figure 14, which shows the residual displacement along the same line. In this case, theresidual motion is the scattered motion from the hill, since
the layered structure beneath the hill is the same as that of
the background structure. The residual motion is significantly greater in the east–west than in the up–down direction. The maximum response of the hill occurs late in the
excitation, and the wave scattering is significantly stronger
than for the basin. The peak value of this scattered motion
is of the order of 35% of the free-field motion at the edgesof the computational domain. This means that the hill’s effect on the free-field ground motion is far from negligible.
The distribution of maximum response of the hill’s freesurface and of its neighboring region is shown in Figure 15.
From this figure it is seen that the prescribed seismic source
excites primarily the fundamental mode of the hill. The peak
amplitudes of the x (east–west) and z (up–down) components
of displacement increase from the base to the top, and the
maximum peak values occur on the eastward side of the crestfor the x component of displacement and on the westward
side and uphill plane for the vertical component. The displacement in the y (north–south) direction is much smaller
than in the previous two cases, but the peak values occur at
the foot of the hill and outside of it. Interestingly, in contrast
to the basin problem, backward scattering here causes deamplification and forward scattering amplification (Fig. 15a).
Even though this example represents an idealized situation,
it suggests that interpreting ground motion in the vicinity of
a topographic feature as free-field ground motion must be
done with caution.
Figures 16–18 show the corresponding results for thehill with the weathered layer. The results are qualitatively
similar to those of the homogeneous hill, except that thesofter layer amplifies significantly the hill response. Compared to the free-field amplitude without the hill, the amplification ratio of the layer east–west displacement component
is about 3.6.
An attractive feature of the DRM methodology exhibitedby the preceding examples is the relative efficiency of theassociated absorbing boundary conditions. We mentioned
earlier that by choosing the residual displacements as theunknown field in the exterior region that surrounds the local
geological features, the residual ground motion in the exterior region is strictly outgoing and corresponds to the deviation of the actual structure from the background structure.
It appears that this perturbation can be small even if theproperties of the local feature differ significantly from thoseof the background structure. In that case, the absorbing
boundary is required to dissipate only a small amount ofenergy compared to that of the free-field motion. This effectis illustrated in Figure 19, which shows snapshots at various

832

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 11. Spatial distribution of maximum absolute value of total displacement
components of free-surface ground motion within the basin and its vicinity. The point
source is located 10 km west (x ⳱ ⳮ10 km). (a) East–west component; (b) north–south component; (c) up–down component. Notice the different scales in each panel.

times taken from an animation of the ground motion for thethree cases considered thus far: the background structure(left-hand column), the basin (middle column), and the homogeneous hill (right-hand column) under the prescribeddouble-couple excitation. The displaced configuration on thevertical plane of symmetry through the line AA⬘ (see Figs.
4a, 7a, and 12a) is superposed on top of the initial configuration for each system. Visible scattered waves emanatefrom the two structures and reach the absorbing boundaries.
These scattered waves, however, are smaller than the freefield ground motion of the background structure. The reason
for this is that some of the input energy is trapped within the
structure and is released only gradually. This implies that theamount of energy that the exterior boundary needs to absorb
at a given instant when the residual wave field is chosen as
the basic unknown in the outer region can be significantly
smaller than that which would need to be absorbed if thetotal displacements were regarded as the unknowns. Thus,
even if the percentage of error is the same for a particular
choice of absorbing boundary condition, its performance canbe expected to be superior for the DRM.
Figure 19 also serves to illustrate the relative responseof the background structures, basin, hill, and their immediate
vicinity. It is clear that at any given instant the response of
the modified systems differs significantly from that of the
background region. Both the basin and the hill exhibit
marked spatial variation over short distances compared tothat of the free-field ground motion. Even though it is stiffer
than the basin in this example, the hill responds morestrongly because the basin is confined within the backgroundstructure, whereas the hill vibrates freely above the freesurface.

Los Angeles Basin
Model Verification
In the preceding examples, because the backgroundstructure consisted of a set of horizontal layers overlying anelastic half-space, we were able to use a theoretical Green’s
function approach to evaluate the free-field motion in step 1of the DRM. If the geometry is complex or the material properties are highly heterogeneous, it becomes necessary to use
a purely numerical procedure, such as finite differences or
finite elements. To test the DRM in a more realistic situation
against the traditional approach, in which the source and theregion of interest are incorporated into a single model, in
this section we consider an 80-km ⳯ 80-km ⳯ 30-km region that encloses entirely the Los Angeles basin and use
the southern California reference three-dimensional seismic
velocity model, version 2 (Magistrale et al., 2000), developed at the Southern California Earthquake Center (SCEC),
to characterize this region. We use the displacement formulation of the finite-element method both for the traditional
approach and for the DRM to determine the ground motion
within a small subdomain (region of interest) of the originalmodel, using as seismic source a buried double couple,
which is located well outside the region of interest.
The SCEC velocity model consists of detailed rule-based
representations of the major southern California basins, embedded in a three-dimensional crust over a variable-depth
Moho. Outside of the basins, the model crust is based onregional tomographic results. Figure 20 presents a plan viewof the shear-wave velocity distribution at the free surface ofthe region considered, together with a vertical cross section
along AA⬘ of the top 15 km of the 30-km deep computa-

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

833

Figure 13.

Synthetic seismograms for displacements of uniform hill along line AA⬘ (Fig. 12a). (a)East–west component; (b) up–down component. Displacements inside C are total; outside they are residual. This is the reason that the seismograms exhibit a
discontinuity across Ce.

Figure 12. Hill on flat-layered system. (a) Finiteelement mesh; (b) cross section through line AA⬘,
which traverses the free surface of the flat-layered
system and that of the hill. Two cases of hill are considered in the simulations: one for a homogeneous
hill, in which its properties are the same as those of
the top surficial layer, and the second in which thehill has a weathered layer with the same properties as
those of the basin in Figure 7.
tional domain across the epicenter of the seismic source. The
scale in the plan view has been capped at 350 m/sec to highlight the large degree of heterogeneity of the model. Thisvelocity is shown to take values beyond 4000 m/sec in the
cross-sectional view.
The small red square shown in Figure 20 represents theregion of interest selected for this demonstration example.
It is 6 km ⳯ 6 km in plan and 0.5 km deep. Due to thecomplexity, heterogeneity, and refinement of the finite-element mesh, we do not sketch the exact location of thetwo bounding surfaces C and Ce on which the effective
forces are to be applied in step 2 (see Fig. 1 for notation).
The lateral limits of this localized region along the crosssection AA⬘ are denoted by vertical yellow lines in Figure20b. Close-up views of the computational domain to be usedin step 2 of the DRM are shown in Figure 21. As in ourprevious models, this domain extends beyond the region ofinterest. The displayed shear-wave velocity exhibits variations between 184 and 274 m/sec at the free surface. Other
regions, outside the limits shown in Figure 21, present shearwave velocities as low as 60 m/sec at some locations.
The source is defined as a strike-slip double couple located at (40 km, 56 km, ⳮ14 km), as shown in Figure 20b,with a strike, slip, and rake of 0⬚, 90⬚, and 0⬚, respectively.
Its seismic moment M0(t) is prescribed as
M0(t) ⳱ M0[1 ⳮ (1 Ⳮ t/T0)eⳮt /T0],
with M0 ⳱ 1 ⳯ 1018 N m, and T0 ⳱ 0.2 sec. The meshgenerated for the simulations is tailored for a maximum frequency of 0.5 Hz; thus, the resulting synthetic ground motions are filtered accordingly.
The verification procedure follows the two steps of the
DRM method: (1) large-scale simulation and calculation of

834

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 14. Synthetic seismograms for residual
displacements along free-surface line AA⬘ in Figure
12. (a) East–west component; the peak value of thisresidual displacement is about twice that of the corresponding value for the flat-layered system. (b) Up–
down component; the peak value is about 70% thatof the corresponding value for the flat-layered system.

the effective forces at the boundaries of the region of interest,and (2) simulation within the reduced domain. To compare
the results of the DRM with those of the traditional procedure, which for this particular problem corresponds to step
1 of the DRM, the velocity response in step 1 is recorded
along two lines of free-surface observation points within theregion of interest, as shown in Figure 21a. The cross section
shown in Figure 21b is taken across the southwest–northeast
diagonal along observation points P. Note that the model
used in step 2 is only 600 m deep. The corresponding synthetics from the traditional approach (step 1) will be compared with those obtained from step 2, in which only the
local region is used in the analysis.
For the first step of the calculations, our elastic wave
propagation finite-element simulation code (Bao et al.,1998) takes as input the original mesh for the entire 80-km
⳯ 80-km ⳯ 30-km model, the DRM limiting surfaces denoted by red lines in the previous figures, the geological andgeometric characteristics of original computational domain,
and the seismic source. Although no material damping is
considered in this model, a simple viscous damping approach is used as in our previous examples to limit the occurrence of spurious reflections at the outer boundaries.
The mesh is partitioned, and a communication graph isdeveloped to distribute the computational load among all
available parallel processors. The finite elements that intersect the limits of the DRM box are tagged as DRM elements
before the beginning of the simulation, and each processor
stores the number and location of its own tagged DRM elements and nodes. The simulation proceeds as a typical wave
propagation analysis, except that in addition to recording
responses at locations of interest, the displacement field isrecorded at the tagged DRM nodes. Parallel synchronization

Figure 15. Spatial distribution of maximum value of displacement components ofground motion on free surface of uniform hill and its vicinity. (a) East–west component;(b) north–south component; (c) up–down component.

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

Figure 16.

Same as Figure 13, but for weathered hill.

is essential for the sequential output procedures since eachprocessor outputs its DRM information to a single output
buffer.
For step 2, the elements that belong to the outer region
and the seismic source are discarded, as shown in Figure 21.
The calculation proceeds with the reduced mesh and the
DRM tagged nodes as multiple seismic sources represented
by the effective forces calculated from the displacement fieldin step 1. This analysis requires much smaller computing

Figure 17.

835

and storage capabilities; it may, therefore, be performed either sequentially or in parallel. In either case, the recorded
DRM node displacements are assigned to the new meshnodes with an interpolation scheme. We have developed an
automated procedure for which the meshes for step 1 andstep 2 need not be identical. Likewise, the simulation timestep may also differ. This represents a clear advantage forcode interaction purposes, as it is common to use different
numerical calculation procedures, meshes, and softwaretools for the large-scale ground-motion simulation (step 1)
and the small-scale ground-motion simulation, soil–structureinteraction, and building response (step 2). In this particular
example, the first and second meshes coincide, and such ascheme is not required. Numerous and repeated numericalsimulations, such as those required by nonlinear analyses orparametric studies, may now be performed with just a fraction of the original computational resources.
The resulting mesh statistics for the present backgroundand local simulations are shown in Table 2. The reduction
in required number of mesh nodes and elements is substantial. This fact translates into considerable computing efficiency for further analysis within the local region. For example, step 1 required 3 hr on 128 parallel processors of the
Cray T3E machine at the Pittsburgh Supercomputing Centerfor 20,000 timesteps of 0.002 sec. In addition, approximately
6 Gb of storage space are used to store mesh and data files.
By contrast, the second step localized calculation may be
performed in a single-processor personal computer with lessthan 2% of the original data storage capabilities. As shownon Table 2, the second step uses only about 1% of the original mesh.
Figure 22 compares the displacement synthetics at theobservation points for the traditional approach (step 1 forthis problem) and the two-step DRM. Clearly, the resulting
synthetics consist both of body and surface waves. The two

Same as Figure 15, but for weathered hill.

836

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 18.

Same as Figure 14, but for weathered
hill. Notice the significant increase of peak response
due to weathering.

sets of results are essentially identical. Notice that eventhough for this example the material properties within thelocalized region are almost uniform in the lateral direction,
the spatial variability of the surface ground motion is quite
strong. The DRM captures accurately the complex groundmotion that is generated as the seismic waves travel fromthe source through the deep and shallow parts of the geological structure within the extended region.
In dealing with the finite-element method, there is one
additional point that should be mentioned here. It is wellknown that the displacement formulation of the finiteelement method fails for incompressible materials (Poisson’s
ratio, ␯ ⳱ 0.5) and leads to inaccurate results as ␯ approaches this limiting value, since the corresponding stiffness matrix becomes nearly singular. For the seismic velocity model of the Los Angeles basin we considered here, there
are many locations where ␯ takes values between 0.4 and
0.44. Independent comparisons of our results with thosefrom finite difference calculations that use an algorithm
based on stress-velocity formulation that is insensitive to
Poisson’s ratio have confirmed the accuracy of our implementation of the traditional finite-element methodology(Day, 2002).

Discussion
The initial motivation for developing the DRM camefrom the desire to deal with problems for which the causative

fault is at some distance from the region of interest. By removing temporarily the local geological feature in step 1,
we showed via the examples in the previous section that it
is possible to greatly simplify the original problem, especially for problems in which some portions of the domain
have very low shear-wave velocity compared to that of the
background geology. This allows one to use coarser meshes
for the background system than would be needed with asingle-step procedure and, therefore, an increased number oftimesteps if the spatially discretized equations of motion aresolved in step 1 by an explicit step-by-step time integration.
Only for step 2 is a finer mesh, and thus smaller timesteps,
required in order to represent accurately the ground motionin the presence of a highly contrasting localized geologicalfeature.
There are other problems for which the DRM may be
advantageous even if the fault is not far from the region ofinterest, for example, for situations in which, because of uncertainty of the geometric and material properties of the local
feature, it may be desirable to repeat the calculations fordifferent combinations of the system parameters, such as inseismic inversion. In that case, step 1 need only be applied
once for a prescribed source. Nonlinear soil behavior extending over a limited region falls in the same category, as
the solution must be determined iteratively. Confining thenonlinearity to step 2 would then be greatly advantageous.
With these applications in mind, we have developed an automated interpolating procedure for which the meshes within
the common domains need not be identical for the two steps.
Similarly, the timesteps in the two simulations may also differ. This allows one to use different codes for the large-scalewave propagation simulation in step 1 and for the smallscale ground-motion simulation in step 2, for which the corresponding code may include provisions for nonlinear material behavior.
For the DRM procedure to be rigorously valid, as we
indicated earlier, the material exterior to X (Fig. 1b) mustbe identical to that of the original problem. However, fromthe numerical results in the preceding section we saw thatfor the basin and hill, the residual wave field in the exterior
region is only a fraction of the complete wave field within
the region of interest. This suggests that one might be able
to simplify considerably step 2 for a general case, yet maintain an acceptable approximation. Nonetheless, we should

Table 2Mesh Statistics for DRM Simulations
Concept

Mesh elements
Mesh nodesInterior boundary nodes
Exterior boundary nodesInterior domain nodes
Exterior domain nodes
DRM elements

Step 1:
Background Calculation

Step 2:
Localized Calculation⬃69,000,000
12,783,986
5,026
5,440
109,439
12,664,081—⬃660,000
119,905
5,026
5,440
109,439—
29,427

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

837

Figure 19.

Snapshots of ground displacement on vertical cross section across AA⬘(plane of symmetry) at various instants, for the
background flat-layered system, the homogeneous basin embedded in the flat-layered system, and the homogeneous hill atop the flatlayered system. Time is measured from theonset of the excitation at the seismic source.
The scale is at the top left of figure. Displacements in the interior region to C are total; thosein the exterior are relative to those corresponding to the background layered system (residual
field). (a) S-wave arrival; (b) multireflection of
S waves and fundamental Rayleigh mode;(c) fundamental Rayleigh mode; (d) first higher
Rayleigh mode. Notice the radiated residual
wave field from basin and hill, which is concentrated primarily on surface layers. Observe,also, the deformation along the boundary of theregion of interest. Full animation, as well as
Figures 9, 10, 13, 14, 16, and 17, shows thenegligible effect of absorbing boundary conditions on ground motion.

point out that if in selecting the region X̂Ⳮ one leaves out
geological features such as a deep layer or a heterogeneitythat is present in the original lithology, the DRM will not be
rigorously equivalent to a single-step procedure that models
the entire region all at once, since any reflections from theheterogeneity or the deep layer due to the residual wave field
will be ignored in step 2. However, provided the backgroundmodel used to determine the free-field motion u0b and u0e is
identical to the original one in the domain XⳭ, then theequivalent seismic forces Peff will be exact to within discretization error. The only approximation will be due to thesecondary reflections generated by the residual wave field,
and these in many cases will be insignificant.

Concluding Remarks
The two-step DRM, described in Paper I and illustrated
by several three-dimensional examples in this article, provides an efficient and reasonably accurate methodology formodeling earthquake ground motion in complex localized
regions with large contrasts in material properties with respect to the background geology. By separating local features with possibly short wavelengths from the backgroundstructure, this methodology can make it possible to model
earthquake ground motion at higher frequencies and withgreater fidelity than has been practical up to now. While thismethod was originally conceived for cases in which the838

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 20. Shear-wave velocity model ofthe Los Angeles basin. (a) Free-surface shearwave velocity distribution, showing a 6-km ⳯
6-km region in which the DRM analysis is performed; (b) cross section AA⬘ shows the shearwave velocity distribution down to 15 km. The
blue zones represent softer soils. Notice thatthe localized region includes the southeastern
portion of the San Fernando Valley. The latitude and longitude of the lower left-handcorner of panel (a) are 33.7275⬚N and
118.9080⬚W.

source is far from the local structure, it can be especially
useful for performing repeated analyses in which the sourceis kept fixed, but the properties of the local feature are varied
from one simulation to the next. This methodology is equally
appropriate if the localized feature exhibits nonlinear behavior or there are engineered structures present within the region of interest. Additional computational savings can begained if the region of interest is restricted to the local geological feature and its vicinity and excludes heterogeneities
or deep layers some distance away. However, errors due tosecondary reflections generated by the residual wave field
will occur if the background region contains heterogeneities
that are not included in the reduced region. While our numerical results indicate that the outgoing waves are small
due to the impedance contrast between the material within
the region of interest and the exterior region, the issue of
secondary reflections deserves further investigation. We be-lieve that the DRM provides a useful tool toward the assessment of seismic hazard and seismic risk reduction in urban
areas within basins with complex topography.

Acknowledgments
This work was supported by the National Science Foundation, underGrant 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. Computing services
were provided by the Pittsburgh Supercomputing Center on the Cray T3E.
This work was also partly supported by a special project of U.S.-Japan
Cooperative Research for Urban Earthquake Disaster Mitigation (No.
11209201), Grant-in-Aid for Scientific Research on Priority Area (CategoryB), Ministry of Education, Culture, Sports, Science, and Technology ofJapan (PI: Tomotaka Iwata), and by Taisei Corporation. We are gratefulfor this support. We thank Peter Moczo and an anonymous reviewer, whose
reviews helped us improve the manuscript. The work was done while oneof the authors (C.Y.) was visiting Carnegie Mellon University as a researchscholar on leave from Taisei Corporation.

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions

Figure 21. Close-up of the DRM region of analysis with the location of the observation points. (a)Plan view with longitudinal (L) and perpendicular (P)cross sections; (b) 600-m deep cross section along P
receiver line. The red line represents the limits of thelocal region. Note that panel (b) shows only the top
600 m of the 30-km deep original model. The displayed values of the shear-wave velocity present variations between 190 and 316 m/sec at the free surface.
However, the model used in the simulations includes
velocities as low as 60 m/sec.

References
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., J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O’Hallaron, J. R.

839

Shewchuk, and J. Xu (1998). Large-scale simulation of elastic wave
propagation in heterogeneous media on parallel computers, Comp.
Methods Appl. Mech. Eng. 152, 85–102.
Bielak, J., K. Loukakis, Y. Hisada, and C. Yoshimura (2003). Domain
reduction method for three-dimensional earthquake modeling in localized regions, Part I: Theory, Bull. Seism. Soc. Am. 93, no. 2, 817–
824.
Day, S. M. (2002). The PEER/SCEC collaboration for validation of 3D wave
propagation codes, Presented at the SCEC Annual Meeting, 8–11 September 2002, Oxnard, California.
Hisada, Y. (1994). An efficient method for computing Green’s function fora 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, II,
Bull. Seism. Soc. Am. 85, 1080–1093.
Loukakis, K. E. (1988). Transient response of shallow layered valleys forinclined incident SV waves calculated by the finite element method,
Master’s Thesis, Carnegie Mellon University, Pittsburgh.
Loukakis, K., and J. Bielak (1995). Seismic response of 2D-valleys: Localsite effects, in Third International Conference on Recent Advances inGeotechnical Earthquake Engineering and Soil Dynamics, St. Louis,
Missouri, 2–7 April 1995, 2, 595–598.
Lysmer, J., and R. L. Kuhlemeyer (1969). Finite dynamic model for infinite
media, J. Eng. Mech. Div. ASCE, 95 (EM4), 859–877.
Magistrale, H., S. Day, R. W. Clayton, and R. Graves (2000). The SCEC
southern California reference three-dimensional seismic velocity
model version 2, Bull. Seism. Soc. Am. 90 (6B), 65–76.
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.
Shewchuk, J. R., and D. R. O’Hallaron (1998). Archimedes, www.cs.
cmu.edu/⬃quake/archimedes.html (last accessed 4 March 2003).
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, Pennsylvania 15213-3890
jbielak@cmu.edu, antoniof@andrew.cmu.edu(C.Y., J.B., A.F.)

Department of Architecture
Kogakuin University
1-24-2 Nishi-shinjuku, Shinjuku-ku
Tokyo 163-8677, Japan
hisada@cc.kogakuin.ac.jp(Y.H.)

Manuscript received 21 September 2001.

840

C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernández

Figure 22. Displacement synthetics (in centimeters). The figure compares the response at the observation points for the traditional single-step and two-step DRM calculations. The results are identical since the effective force nodes of the first- andsecond-step analyses are those from the original mesh; that is, there are no additional
spatial discretization errors between steps 1 and 2. The readings along line L (longitudinal) show more pronounced phase differences than those along line P, consistent
with the location of the hypocenter, the magnitude of the shear-wave velocity, and thedistance between the observation points.