This section presents the parallel numerical techniques we use to approximate the solution of the wave propagation equations. 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
according to:
The domain of the problem is semi-infinite, yet our computational domain must be rendered finite. Here we chose a rectangular parallelepiped whose size is determined by the extent of the valley. We prescribe absorbing boundary conditions that are local in both time and space on all sides of the computational domain except for the top surface, which is traction free. In this application, we use the simplest possible absorbing boundary--a damper. While a viscous damper is a sufficient choice for high frequencies and normal wave incidences, it is, in general, a poor approximant of the exact condition at the truncation interface. However, we choose the truncation surfaces so that they lie outside the sedimentary valley, that is, mostly in rock. In this way, we ensure that the amplitudes of the waves impinging upon the truncation boundaries will have been greatly reduced, thus minimizing spurious reflections.
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 [4, 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 right-hand side. These non-homogeneous terms on the right-hand 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 a 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 Rayleigh 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 a 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, the system of ODEs (3) can be integrated in time using central differences, yielding the explicit method
This method exhibits second-order accuracy in time; 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 (4) with integration points at
element vertices. This results in a diagonal mass matrix. To render
the left-hand side operator of (7) 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 right-hand side of
(7), 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
(7). 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 (7). For
example, on a regular grid, 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 suggests 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.
