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.