Whereas the mesh generation, partitioning, and parceling steps are currently performed sequentially, the wave propagation equations are solved on a parallel machine. This section describes the numerical techniques we use. Navier's equations of elastodynamics for an isotropic, heterogeneous medium are
where u is the displacement vector field, 
is the density,
and 
and 
are elastic material constants, which depend
on the shear and dilatational wave velocities.
The domain of the problem is an elastic halfspace, i.e.\ semi-infinite. In order to render the computational domain finite, we introduce absorbing boundaries at the bottom and sides of the basin that are local in both space and time [4, 10]. These boundaries allow the passage of outgoing waves with minimum reflection.
Since, in many cases, the earthquake source can be outside the computational domain, its effect must be introduced into the region. This is carried out as described in [3, 6] by means of effective forces. In short, for an arbitrary earthquake excitation these forces are determined in terms of the free-field motion by introducing a fictitious auxiliary surface that surrounds the basin. Across this auxiliary surface one imposes the conditions of continuity of displacement and traction. By selecting the total displacement vector field as the unknown in the resulting interior region and the scattered displacement field in the exterior region, the free-field displacement and traction now appear explicitly in the continuity conditions, which become jump conditions, with the free-field displacement and traction on the righthand side. These non-homogeneous terms on the righthand side are the ones that give rise to the effective forces upon spatial discretization. If, on the other hand, the seismic source is located inside the computational domain, say as a kinematic dislocation across the fault, one can select the fault itself as the auxiliary surface. The procedure is similar, but now one uses the total displacement everywhere as the unknown field; thus, the displacement field again experiences a jump across the interface, but the traction remains continuous. Notice that with this technique, whether the source is originally located inside or outside the computational domain, only outgoing waves will impinge upon the absorbing boundary. Both types of source are implemented in our code.
We also model material damping in the basin via viscous damping. With these modifications, standard Galerkin discretization in space by finite elements produces a system of ordinary differential equations of the form
where M is the mass matrix, C is the damping matrix
associated with the absorbing boundary and material damping, K
is the stiffness matrix, and f is the effective force
vector. Here, M, C, and K are block matrices; the
th block of M is a 
matrix given by
and the 
th block of K is given by
where 
is the finite element global basis function associated
with the i-th node.
Damping is introduced through a proportional damping approximation at the element level, i.e. we take
where 
and 
are scalar constants and the superscript
e indicates an element matrix. The first term leads to as damping
factor that is inversely proportional to frequency, and the second to
one that is linear in frequency. The constants 
and 
,
which may vary within the basin according to the type of material, are
chosen to best fit a prescribed attenuation law.
Given appropriate initial conditions, this system of ODEs can be integrated in time using central differences, yielding the explicit method
This method exhibits second-order accuracy in time, and when coupled
with linear finite elements, we obtain second-order accuracy in space
as well. We use a lumped mass approximation to M, which amounts
to numerically integrating (3) with integration points at
element vertices. This results in a diagonal mass matrix. To render
the lefthand side operator of (6) diagonal, we further
evaluate the off-diagonal components of 
at 
rather than 
. Inversion of the time stepping
operator thus requires only a scaling of the righthand side of
(6), which is carried out just once prior to time
stepping. Forming the products of K and C with vectors
comprises the major computational effort associated with iterating on
(6). The sparsity structure of K is dictated by the
underlying finite element mesh, and is thus very irregular. If shear
waves are not over-resolved, the time step necessitated by stability
is of the order of the time step dictated by accuracy, which is what
an implicit method would take. By choosing an explicit method we avoid
solving linear systems at each time step. Thus, overall, the explicit
method is superior for our application, especially on a parallel
computer.
We have tried several different choices of basis function order, and
have concluded that piecewise-linear functions are the most efficient
for problems requiring engineering accuracy. Our conclusion is based
on numerical experimentation using plane Ricker wavelets on
unstructured homogeneous meshes (in which case we know what the exact
solution should be), but a simple argument can be given as follows. We
recognize first that (spatial) approximation errors are bounded from
above by interpolation errors. We then ask, for a given order of basis
function and a given acceptable level of infinity-norm error, how many
nodal points are required to produce a piecewise-polynomial
interpolant of a simple harmonic wave. Next, we convert the required
number of nodes per wavelength to an estimate of the storage and work
required for an iteration of the explicit method (6). For
example, if N is the total number of nodes, one can show that
trilinear hexahedra require 
words of storage and 498N
flops/time step, while triquadratic hexahedra necessitate 367N words
and 1164N flops/time step. So, for example, triquadratic elements
should require at least 2.2 times fewer nodes in order to be preferred
(for storage reasons) over trilinear elements. However,
one-dimensional interpolation theory tells us that 5% error requires
10 nodes per wavelength using linear elements or 9.4 nodes using
quadratics. Thus, in three dimensions, triquadratics only allow
times fewer nodes than trilinears, and are thus not
warranted. An opposite conclusion is reached if one demands 99%
accuracy. Our confidence in the values of material properties and in
the fidelity of the source models for this problem does not warrant
solution accuracies greater than 95%. Thus, we conclude that for this
level of accuracy, the powers of higher-order interpolation are offset
by their increased cost, both in storage and in increased work.
