The reduction of earthquake risk to the general population is a major problem facing countries located in highly seismic regions. Assessing the free-field ground motion to which a structure will be exposed during its lifetime is a critical first step for the design of new structures and retrofit of existing ones.
The main objective of this paper is to describe a methodology for modeling ground motion on parallel computers in large sediment-filled basins during earthquakes, and to illustrate with an application to a real basin. Observations during recent strong earthquakes have shown that three-dimensional local site effects, caused by the waves generated inside a basin, can result in strong ground motion amplification and longer duration of the surface ground motion, with respect to that in rock, and a rapid spatial variation of the ground motion that can cause large differential base motion of extended structures such as bridges or dams [19, 10, 3]. See [1] for a general overview, and [12, 11, 24, 8, 9, 18, 16], for instance, for representative recent work in this field.
Simulating the earthquake response of a large basin is accomplished by solving numerically the partial differential equations of elastic wave propagation, i.e. the Navier equations of elastodynamics. Several numerical methods have been used for approximating the solution to these problems. Boundary element methods have been popular for recovering solutions primarily in the frequency domain and for moderately-sized linear problems. Inhomogeneities, nonlinearities, the large scale of such basins as Los Angeles and Kanto, and the desire to model directly in the time domain, preclude their use here. On the other hand, uniform grid domain methods such as structured finite differences become impractical for the very large problem sizes involved.
To see why uniform grids are impractical, consider the Los Angeles
Basin. For a shear-wave velocity of 0.4 km/s and a frequency of 2 Hz, a
regular discretization of the elasticity operator would place grid
points 0.02 km apart to achieve second order accuracy. The region of
interest has dimensions 140 km 
100 km 
20 km; thus, a
regular discretization, governed by the softest layer, requires 35
billion grid points with three displacement components per grid
point. At least a terabyte of primary memory would be needed, and on
the order of 
operations would be required at each time
step. The stability condition associated with explicit time
integration of the semidiscrete equations of motion imposes a time
increment at least as small as 0.004s. Thus, a computer would have to
perform at a sustained teraflop per second for two days to simulate a
minute of shaking.
Instead, we use unstructured mesh finite element methods that tailor the mesh size to the local wavelength of propagating waves [2]. For a basin such as Los Angeles, the shear-wave velocity varies from 220 m/s to 4500 m/s throughout the basin and its vicinity. Since in three dimensions mesh density varies with the cube of shear-wave velocity, and since the softest soils are concentrated in the top layers, this means that an unstructured mesh method may yield three orders of magnitude fewer equations than with structured grids. Modeling the Los Angeles Basin for values of earthquake period and wave velocity that are desirable for engineering purposes thus becomes practical on the largest of today's parallel supercomputers.
We favor finite element methods for their ability to efficiently
resolve multiscale phenomena, the ease with which they handle traction
boundary conditions, and their firm theoretical foundation. For
temporal approximation, we have studied both explicit and
preconditioned conjugate gradient-based implicit methods. For
hyperbolic problems, explicit methods become unstable if the time step
is greater than the time it takes an elastic wave to cross any
element--the Courant condition. Unconditional stability, on the other
hand, can be achieved if one uses implicit methods. This implies that
larger time steps can be taken. However, the very characteristic that
makes them stable--the fact that the solution at a node at time
requires information from all nodes at time t as
opposed to just the neighbors--renders them unattractive on
distributed memory computers, since this implies global information
exchange. One approach to making implicit methods efficient on
parallel computers is to use iterative methods for their solution,
effectively rendering them explicit. However, we have found that our
mesh generators give us such good control over mesh resolution that
the Courant condition for explicit methods is not onerous. The result
is that the more readily parallelizable explicit methods perform
better for elastic wave propagation problems. In this paper we
consider only a single-step explicit time integration method.
While unstructured mesh methods for simulating wave propagation through heterogeneous media result in many fewer equations, they introduce a number of computational difficulties that must be overcome. First, mesh resolution must closely follow wavelength; too coarse a resolution will introduce error, too fine will result in unnecessary computation as well as excessively small time steps (when explicit integration methods are used). Second, element aspect ratios must remain small; large aspect ratios will eventually result in instability in the time integration scheme. Highly heterogeneous basins, in which wavelengths vary rapidly in space, introduce special difficulties when trying to follow the wavelength change without severely stretching the mesh. Third, unstructured mesh methods are not easy to program on parallel computers; their irregular data structures require nontrivial mappings onto parallel machines, and irregular communication patterns are generated. Thus, we have developed fast, robust computational geometry and mesh generation techniques for highly spatially-variable meshes, and compilers and tools that simplify the programming of unstructured mesh methods on parallel systems.
For an alternative approach to parallel ground motion modeling on distributed memory machines, see e.g., the work of Olsen, Archuleta, and Matarese [17], which employs finite differences on regular grids. See also the references to prior finite difference modeling work on sequential machines contained therein. In addition to the finite element method described in this paper, there have been recent efforts to endow finite difference wave propagation methods with multi-resolution capabilities. See the work described in [14], which uses composites of regular grids to achieve variable resolution.
In the remainder of this paper, we present numerical methods and geometric algorithms for modeling earthquake-induced ground motion in highly heterogeneous basins. We also describe an automatic code generator for solution of unstructured mesh PDE problems on parallel distributed memory computers. As an application, we model the response of the San Fernando Valley in Southern California to an aftershock of the 1994 Northridge earthquake, and give performance results on the Cray T3D.