After many numerical experiments to optimize techniques that could be adapted to large scale FEM models, we have chosen and developed the following algorithms and techniques in this study:
As a basis for understanding this model, more detail is given here.
We use the finite-element method to solve the elastodynamics equation, i.e., Navier's equation for an isotropic, heterogeneous medium,
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 and are taken to be piecewise
constant.
The domain of the problem is an elastic halfspace. In order to render the computational domain finite, we introduce absorbing boundaries at the bottom and sides of our finite computational domain that are local in both space and time [5, 22]. 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 [4, 14] 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 the fault, one can obtain the seismic moment tensor which depends on both the seismic moment and fault orientation. The seismic moment tensor must be applied as a set of effective forces applied on nodes adjacent to the fault. 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.
Rayleigh 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 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, this system of ODEs can be integrated in time by using central differences to yield the explicit method
This method exhibits a second-order accuracy in time, and, when it is 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 left-hand 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 right-hand 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 sparse 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 needed for 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. A detailed description of parallel computing techniques and
algorithms [3] have been presented at Supercomputing '96.
