The Spectral Element Method: An Efficient Tool To Simulate The Seismic Response of 2D and 3D Geological Structures

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/232707622
The Spectral Element method: an efﬁcient tool to simulate the seismic response
of 2D and 3D geological structures
Article in Bulletin of the Seismological Society of America · April 1998
CITATIONS READS
1,021 3,633
2 authors, including:
Jean-Pierre Vilotte
Institut de Physique du Globe de Paris, Paris, France
185 PUBLICATIONS 6,793 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
Belmont Forum e-Infrastructures and Data Management Project View project
Dynamic rupture modelling View project
All content following this page was uploaded by Jean-Pierre Vilotte on 27 July 2014.
The user has requested enhancement of the downloaded file.

Bulletin of the Seismological Society of America, Vol. 88, No. 2, pp. 368-392, April 1998
The Spectral Element Method: An Efficient Tool to Simulate

the Seismic Response of 2D and 3D Geological Structures
by Dimitri Komatitsch and Jean-Pierre Vilotte
Abstract We present the spectral element method to simulate elastic-wave prop-

agation in realistic geological structures involving complicated free-surface topog-
raphy and material interfaces for two- and three-dimensional geometries. The spectral
element method introduced here is a high-order variational method for the spatial
approximation of elastic-wave equations. The mass matrix is diagonal by construc-
tion in this method, which drastically reduces the computational cost and allows an
efficient parallel implementation. Absorbing boundary conditions are introduced in
variational form to simulate unbounded physical domains. The time discretization is
based on an energy-momentum conserving scheme that can be put into a classical
explicit-implicit predictor/multi-corrector format. Long-term energy conservation
and stability properties are illustrated as well as the efficiency of the absorbing con-
ditions. The associated Courant condition behaves as At c < 0 (n~ lind N-2), with nel
the number of elements, na the spatial dimension, and N the polynomial order. In
practice, a spatial sampling of approximately 5 points per wavelength is found to be
very accurate when working with a polynomial degree of N = 8. The accuracy of
the method is shown by comparing the spectral element solution to analytical solu-
tions of the classical two-dimensional (2D) problems of Lamb and Garvin. The flex-
ibility of the method is then illustrated by studying more realistic 2D models in-
volving realistic geometries and complex free-boundary conditions. Very accurate
modeling of Rayleigh-wave propagation, surface diffraction, and Rayleigh-to-body-
wave mode conversion associated with the free-surface curvature are obtained at low
computational cost. The method is shown to provide an efficient tool to study the
diffraction of elastic waves by three-dimensional (3D) surface topographies and the
associated local effects on strong ground motion. Complex amplification patterns,
both in space and time, are shown to occur even for a gentle hill topography. Exten-
sion to a heterogeneous hill structure is considered. The efficient implementation on
parallel distributed memory architectures will allow to perform real-time visualiza-
tion and interactive physical investigations of 3D amplification phenomena for seis-
mic risk assessment.
Introduction
The use of elastic-wave equations to model the seismic bility of a low-order method with the exponential conver-
response of heterogeneous geophysical media with topog- gence rate associated with spectral techniques and suffers
raphy and internal interfaces is a subject that has been in- from minimal numerical dispersion and diffusion.
tensively investigated by seismologists. The challenge is to Two types of problems have motivated this study. One
develop high-performance methods that are capable of solv- is elastic waveform modeling, in order to understand and
ing the elastic-wave equations accurately and that allow one extract quantitative information from complex seismic data.
to deal with large and complicated computational domains Realistic geological media that present complicated wave
as encountered in realistic 3D applications. phenomena require methods providing solutions of high ac-
This article describes a practical spectral element curacy that correctly simulate boundary conditions, surface
method to solve the 2D and 3D elastic-wave propagation in topography, and irregular interfaces with nonhomogeneous
complex geometry. The method, which stems from a weak properties. The other problem is related to the assessment of
variational formulation, allows a flexible treatment of site effects in earthquake ground motion. In particular, 3D
boundaries, or interfaces, and deals with free-surface bound- surface topography, local velocity variations, and layering
ary conditions naturally. It combines the geometrical flexi- can produce complicated amplification patterns and energy
368
The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures 369
scattering. Such effects can modify the ground shaking to a Alternatively, numerical solutions to the wave equation
large extent and are relevant for the seismic design of struc- have been sought via techniques based on integral represen-
tures. tations of the problem relating quantities on the physical
The need of new high-performance methods for the boundaries. Integral formulations make use of fundamental
elastic-wave propagation can be simply assessed when look- solutions as weighting functions together with Green's the-
ing back to the continuous efforts that have been devoted to orem (Manolis and Beskos, 1988; Bonnet, 1995). In seis-
this subject. mology, these methods can be traced back to the pioneering
Finite-difference methods have been widely imple- work of Aki and Lamer (1970) who used a discrete super-
mented with a varying degree of sophistication. Unfortu- position of plane waves. They have since been extended by
nately, conventional schemes suffer from "grid dispersion," many authors (Bouchon, 1979; S~inchez-Sesma and Rosen-
near large gradient of the wave field, or when too-coarse blueth, 1979; Dravinski and Mossessian, 1987; Horike et al.,
computational grids are used. For realistic applications 1990; Ohori et al., 1992). While direct boundary element
(Frankel, 1993; Olsen and Archuleta, 1996; Pitarka and Ir- methods formulate the problem in terms of the unknown
ikura, 1996), balancing of the trade-off between numerical tractions and displacements (Zhang and Chopra, 1991), in-
dispersion and computational cost turns out to be rather dif- direct methods make use of a formulation in terms of force
ficult. For classical second-order centered finite-difference and moment boundary densities (S~nchez-Sesma and Cam-
methods, at least 15 points must be used for the wavelength pillo, 1993). The combination of discrete wavenumber ex-
corresponding to the upper half-power frequency (Kelly et pansions for Green's functions (Bouchon and Aki, 1977;
al., 1976; Alford et al., 1974). Grid dispersion and aniso- Bouchon, 1979), either with indirect boundary integral rep-
tropy can be reduced when using the staggered-grid for- resentations (Campillo and Bouchon, 1985; Gaffet and Bou-
mulation (Madariaga, 1976; Virieux, 1986; Levander, 1988), chon, 1989) or with direct methods (Kawase and Aki, 1989;
which is based on the symmetric first-order hyperbolic form Kim and Papageorgiou, t993), has lead to successful meth-
of linear elastodynamics (Hughes and Marsden, 1978). This ods with the advantage of seeking solutions over a domain
can also be achieved by using fourth-order centered schemes one dimension lower than the original form of the problem,
both in space and time, based on modified wave-equation with sources at the boundary (removing the uncertainty
techniques (Dablain, 1986; Bayliss et al., 1986). Another about their location), and an a priori satisfaction of the ra-
difficulty with finite differences is their inability to imple- diation condition. On the other hand, methods of this kind
ment free-surface conditions with the same accuracy as in are most often limited to linear and homogeneous problems
the interior regions of the model and their lack of geomet- and are known to encounter difficulties, such as possible
rical flexibility. Even though some techniques have incor- nonuniqueness of the solution of the continuous boundary
porated surface topographies using methods based on grid integral equations at characteristic wavenumbers of the cor-
deformation or vacuum-to-solid taper (Boore, 1972; Jih et responding interior problems, leading to ill-conditioned dis-
al., 1988; Robertsson, 1996; Ohminato and Chouet, 1997) crete equations if left uncorrected. Moreover, the resulting
combined with the staggered grid formulation, they often linear systems of equations in these methods are very large,
remain limited to simple geometrical transformations and nonsymmetric, and dense. The computational advantage in
may affect the stability criterion in the case of grid-defor- processing time and storage requirements that would be ex-
mation techniques, or they require up to 15 grid points per pected intuitively is therefore not always realized in the case
shortest wavelength in the case of vacuum-to-solid tech- of realistic problem sizes.
niques, which puts some limitations for narrow free-surface Spectral methods, introduced in fluid dynamics around
structures. All these ripples make finite-difference methods 20 years ago by S. A. Orszag, have also been proposed for
difficult to use for simulating Rayleigh and interface waves elastodynamics (Gazdag, 1981; Kosloff and Baysal, 1982).
in practical situations. To deal with general boundary conditions, a set of algebraic
Although more suited to heterogeneous media with polynomials (Chebyschev or Legendre in space) replaced
complicated geometries, finite-element methods, based on a the original set of truncated Fourier series. The so-called
variational formulation of the wave equations that allows a global pseudo-spectral method (Kosloff et al., 1990) became
natural treatment of free-boundat-y conditions, have attracted one of the leading numerical techniques in the 1980s in view
somewhat less interest among seismologists (Lysmer and of its accuracy, in terms of the minimum number of grid
Drake, 1972; Toshinawa and Ohmachi, 1992). Apparently, points needed to represent the Nyquist wavelength for non-
the main reason for that is that low-order finite-element dispersive propagation. In these methods, numerical solution
methods exhibit poor dispersion properties (Marfurt, 1984), is derived so as to satisfy the wave equation in differential
while higher-order classical finite elements raise some trou- form at some suitably chosen collocation points. The accu-
blesome problems like the occurrence of spurious waves. racy is shown to depend strongly on this choice. Unfortu-
Recently, space-time finite-element methods based on Ham- nately, global spectral methods suffer from severe limita-
ilton's principle, or on time-discontinuous Galerkin formu- tions: nonuniform spacing of the collocation points for
lation, have been introduced for elastodynamics (Hulbert algebraic polynomials puts stringent constrains on the time
and Hughes, 1990; Richter, 1994). step that cannot be easily removed (Kosloff and Tal-Ezer,
370 D. Komatitsch and J.-P. Vilotte
1993); complicated geometries and heterogeneous material including physical and external boundaries, and t E I =
properties cannot be handled easily nor, when the method is [0,7], with I the time interval of interest.
based on a strong formulation of the differential equations, For elastic-wave propagation, the equations of motion
realistic free-surface boundary conditions. The use of cur- can be written in compact notation form as
vilinear coordinate systems has been proposed to overcome
such a limitation (Fomberg, 1988; Tessmer et al., 1992; Car- p+ = div [o1 + f, (1)
cione and Wang, 1993; Komatitsch et al., 1996) but remains
restricted to smooth global mappings of little use for realistic pu = pv, (2)
geological models. Another idea is to couple domain decom-
position techniques to spectral discretization (Canuto and with the initial conditions
Funaro, 1988; Carcione, 1991); however, this requires a sig-
nificant increase of the computational cost. u(x, t) = Uo(X) (3)
Understanding of the similarity between collocation
methods and variational formulations with consistent quad- and
rature (Gottlieb, 1981; Maday and Quarteroni, 1982) leads
in fluid dynamics to the spectral element method (Patera, v(x, 0 = v0(x) (4)
1984; Maday and Patera, 1989) that may be related to thep
and h - p versions of the finite-element methods (Babugka on part of the internal boundary 1-'~t,
et al., 1981; Babugka and Don', 1981). These methods,
which bring new flexibility to treat complex geometries, tr(x, t) • n(x) = T(x, t), (5)
have been proposed for wave propagation recently by Priolo
et al. (1994) and Faccioli et al. (1996). and on the other part F~nt,
This article describes a practical spectral element
method to solve 2D and 3D elastic-wave propagation in u(x, t) = g(x, t), (6)
complex geometry. The potentialities of the method are
demonstrated on various 2D and 3D problems. In contrast where p = p(x) is the mass density; fix, t) is a generalized
with the method used by Priolo et al. (1994), our formulation body force; u0(x) and Vo(X) are, respectively, the initial dis-
is based on Legendre polynomials and Gauss-Lobatto Le- placement and velocity fields; 6(x, t) is the stress tensor; T(x,
gendre quadrature, leading to fully explicit schemes while t) is the prescribed boundary traction (Neumann condition);
retaining the efficient sum-factorization techniques (Orszag, and g(x, 0 is the prescribed displacement (Dirichlet condi-
1980). Although the particular choice of the sets of algebraic tion). A dot over a symbol indicates partial differentiation
polynomials (Chebyshev or Legendre) and collocation with respect to time. In component forms, div [a] is o-i:J and
points (related to the numerical quadrature) does not gen- 6 - n is aijn j.
erally affect the error estimates significantly, it greatly af- Two simple source terms are considered here: a point
fects the conditioning and sparsity of the resulting set of force,
algebraic equations and is critical for the efficiency of a par-
allel iterative procedure (Fisher, 1990). f(x, t) = f¢~i 6(x - Xo) 9(t - to), (7)
and a body force derived from a seismic moment tensor den-

Formulation of the Problem sity distribution,
Initial Boundary-Value Problem f(x, t) = -div [mo(x)] 9(t - to), (8)

When solutions are assumed to extend to infinity along
some directions, a fundamental obstacle to the direct appli- where m0(x ) is a symmetric tensor and ~(t) is a Ricker wave-
cation of numerical methods is the presence of an unbounded let in time. The stress is determined by the generalized
domain. The boundary-value problem is therefore converted Hooke' s law:
to a formulation that is defined over a bounded region by
introducing an artificial external boundary with appropriate aU (x) = c~jk~(x) u~,t (x, t), (9)
boundary conditions.
We consider an elastic inhomogeneous medium occu- where uk, l = OutfOx1 is the displacement gradient. The elastic
pying an open, bounded region f~ C R he, where nd is the coefficients Cgjk~= Qil, t(x) are positive definite and have all
number of space dimensions (2 or 3 here). The boundary of the required symmetries. By the minor symmetries cii~z ----
f~ is denoted F and can be decomposed into F = F int t_J Cjikl = Colk, ~ depends on the symmetric part of the displace-
F ext, where F e×t is the artificial external boundary. The dis- ment gradient only.
placement and velocity vectors are denoted by u(x, t) and The representation of the radiation condition associated
v(x, t), respectively, where x E ~ , with ~ the closed region with the external boundary is a difficult problem, and nu-
The Spectral Element Methoc# An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures 371
merous approximate schemes have been proposed in the q2 = {w(x) E H I ([~),d : f~ ___>Rne;
geophysical literature, see Bayliss and Turkel (1980) and w(x) = 0 onFg },
int (t2)
Givoli (1991) for a review. Exact nonlocal boundary con-
ditions employing an asymptotic expansion of the far-field the weak form of the governing equations can be obtained
solution to generate a sequence of local boundary operators by multiplying equations (1) through (4) by time-indepen-
(Givoli and Keller, 1990) have now been derived. They are, dent test functions w. Integration by parts and the use of
however, computationally expensive. We assumed here a boundary conditions lead to the abstract formulation
simple local approximation based on the variational formu-
lation of the paraxial condition originally introduced by (w, pi,) + a(w, u) = (w, f)
Engquist and Majda (1977) and Clayton and Engquist
(1977). The local transient impedance of the artificial bound-
+ (w, T)ia¢, + (w, t)rex,, (13)
ary is approximated by use of a limited wavenumber expan-
sion of the elastodynamics equation in the Fourier domain (w, pti) = (w, pv), (14)
along the boundary surface. Such an approximation is ac-
curate for high-frequency waves and for waves impinging with
on the boundary at small angles from the normal direction
only. In the following, a first-order approximation, close to [w, pu(., t)lt=o] = (w, puo), (15)
the one proposed by Stacey (1988), is retained. On the ar-
tificial boundary F ext, the condition is therefore expressed as [w, pv(., t)l,=o] = (w, pVo), (16)
where a(., .) denotes the bilinear form that expresses the
t = Cp p[v" n]n + c,pvr, (10)
virtual work of the internal stresses, defined as
where t is the traction vector on the boundary, n is the unit

outward normal to the surface, vr = v - [v-n]n is the a(w,u) : f o ' V w dV= fnVw'e'Vu dV, (17)
projection of the velocity field on the surface, and cp and c~
are the propagation velocities of P and S waves, respectively.
where in component form ~ : Vw = aoOwi/Oxi with the im-
Such a damping condition was originally proposed by Lys-
plicit summation convention and c is the elastic tensor de-
met and Kuhlemeyer (1969).
fined in equation (9).
The bilinear forms (., .) can be interpreted as the virtual
The Variational Form of the Governing Equations
work of the inertial terms,
While some methods of approximation directly start
from the previous formulation of the initial boundary-value
problem (strong form)--the most notable example being the (w, pi,) = fn pi, • w dV, (18)
finite difference method--a less restrictive approach consists
in considering a weak formulation, or variational formula- and of the external forces,
tion, of the original problem that admits a broader range of
solutions in terms of regularity or smoothness. The most
commonly used formulation is based on the principle of vir- (w, f) = ~n f" w dV;
tual work or virtual displacement (Hughes, 1987; Szab6 and
Babugka, 1991).
(w, T)r~7, = fret T " W dF; (19)
The solution is then searched in the space of the kine-
matically admissible displacements that is defined, accord-
ing to the Dirichlet boundary conditions, as (W, t)roxt = frext t " W dF.
,St = {U(X, t) ~ H 1 (~'~)na : ~'~ Numerical Discretization
× I-~Rne; u(x, 0 = g(x,t) on Fignt × I}, (11) In this section, the Legendre spectral element discreti-
zation of the variational statement of the elastic-wave equa-
tions is outlined.
where the subscript t of St refers to time and H 1 denotes the Like in a standard finite-element method, the original
space of square-integrable functions that possess square- domain is discretized into ne~ nonoverlapping quadrilateral
integrable generalized first derivatives. elements: ~ -- Ue=lne~~e- The restriction of w to the element
Introducing the space q/of the test functions that, in the ~e is denoted Whls~e. Each element ~e is mapped onto a
case of the virtual work formulation, is also called the space reference volume [ ] that is defined, in a local ~ system of
of the virtual displacements, coordinates, as a square or a cube A na with A = [ - 1, 1].
372 D. Komatitsch and J.-P, Vilotte
The invertible element mapping CFeis defined as: element method, to a coupled system of ordinary differential
C/Ze: [ ] ---) ~'r~e s u c h t h a t x ( ~ ) = CFe(~). equations:
With respect to the spatial discretization, we shall re-
quire that the variational statement be satisfied for the piece- M% = 1Text
--t -- Ftint (ut, vt). (23)
wise-polynomial approximation spaces S~ × Vh, where h
denotes the characteristic length scale associated with the
underlying mesh, defined as Let nnodebe the total number of nodes of the global integra-
tion grid EN defined as the assembly of the element domain
~'e.
s)v = { u h ~ s : u h ~ C2(f~) "~ integration grids EN = tOe--N, then ut and v t denote the n~od~
(20) displacement and velocity vectors, respectively, at a given
and uhl~e 0 c12e ~ [PN (/7)] rid} time t; F int is the internal nodal force vector; and F ext is the
external forcing vector that includes the source term and the
and radiation conditions. An attractive property of the method is
that, due to the consistent integration scheme and the use of
= {w h ~ v: w h e L 2 ( ~ ) ~" Gauss-Lobatto Legendre formulas, the mass matrix M is by
(21)
and whlae 0 % E [PN (H)lna}, construction always diagonal.
The spectral-element method therefore combines the
where L2(~) denotes the space of square-integrable func- geometric flexibility of the finite-element method with the
tions defined on ~ and [ p ~ ) ] , a denotes the tensor-product fast convergence associated with spectral techniques. Con-
space of all polynomials of degree N, or less, in each of the sidering only spatial errors, an exponential convergence can
na spatial directions within the reference volume [~. The be ensured (Bernardi and Maday, 1992) for the spectral-
spectral-element spatial discretization can be characterized element approximation when n~t is fixed and N --+ oo. Such
by the discretization pair (nel, IV). Each element integral in- a superconvergence derives from the good stability and ap-
volved in the variational formulation, defined over the do- proximation properties of the polynomial spaces and from
main ~ in the x space, is pulled back, using the local map- the accuracy associated with the Gauss-Lobatto Legendre
ping CFe,on the parent domain [ ] and numerically integrated quadrature and interpolation. The discrete solution suffers
using the numerical quadrature defined as the tensor product from minimal numerical dispersion and diffusion, a fact of
of the 1D Gauss-Lobatto Legendre formulas. In order to take primary importance in the solution of realistic geophysical
advantage of efficient sum-factorization techniques, the problems, see Tordjman (1995) and Komatitsch (1997) for
(N + 1) na basis points for PN are taken to be the same as more details.
the quadrature points on each element ~e and define a col- The semi-discrete momentum equation is then enforced
location grid E~v = {{i, qj, ffk} that is the ha-tensor product in conservative form at t~+ ~, using two parameters fl and 7
of the N + 1 Gauss-Lobatto Legendre integration points. (Simo et aL, 1992):
The piecewise-polynomial approximation w~ of w is
defined using the Lagrange interpolation operator I u on the
1M [Vn+ 1 -- Vn ] = Fen~+ta - F int ( n n + a , V n + a ) , (24)
Ganss-Lobatto grid --N"~e"XN(Wlae) is the unique polynomial At
of PN(D) that coincides with wlh~ at the (N + 1)"a points
of ~N"
"a~ If I~(~) denotes the characteristic Lagrange polyno- u,+l = un + At[(l - ~)Vn + fl-Vn+l]y (25)
mial of degree N associated with the Gauss-Lobatto point i
of the 1D quadrature formula, the approximation of wla~ is
defined as
whla, (X, y, z) = xN (win) (22) a~+l - ?At - an, (26)

N
= ~ If (¢) l~/(tl) l~ (() W~k, where Un+a -- OeUn+l + (1 -- a)Un and Fn+ exta - O/lPn
~,ext
i,j,k = O +1 +
(1 - a)r~. . . . t . Noteworthy properties of this algorithm are
with x = qce(~i, t/j, ~k) and W~k = WhNIf~,O C/~e(~i, r/j, ~). exact conservation of the total angular momentum for a =
A polynomial of degree -<_ 2N - 1 can be integrated fl/? = 1/2 (these values define an acceleration-independent
exactly with N + 1 points. Such a consistent integration, time-marching algorithm) and second-order accuracy if and
where the quadrature points are the same as the basis points, only if a = 1/2; for o~ = fl/? = 1/2, a linear analysis shows
is shown to be sufficient for complex geometries (Maday that the spurious root at zero sampling frequency vanishes
and R0nquist, 1990), an important consideration when wave if and only if 7 = 1 (Simo et al., 1992). This Newmark-
equations are to be solved on significantly deformed ge- type scheme can be generalized to a predictor-multi-correc-
ometry or with variable coefficients. tor format (see Appendix B) that improves its properties and
The procedure outlined above leads, like in the finite- allows an efficient parallelization.
Numerical Implementation. In the simulations presented in the Jacobian matrix becomes nondiagonal, allowing, even
this article, the parameters of the predictor-multi-corrector for small geometrical distortion, one to check the full ac-
algorithm are a = 1/2, fl = 1/2, and ~ = 1. The Courant curacy of the method and the numerical integration. The
number can be defined as n c = max [cAt/Ax], where c is analytical solution (see Pilant, 1979) is calculated using the
the compressional wave speed and Ax is the collocation grid Cagniard-de Hoop technique (de Hoop, 1960) and con-
spacing. The associated Courant condition can be shown to volved numerically with the source-time function.
behave as Atc < O(n~ lind N - 2). The physical problem is a homogeneous elastic half-
The method is shown to work accurately with a low space with a P-wave velocity of 3200 m- sec- 1, an S-wave
number of grid points per minimum wavelength, corre- velocity of 1847.5 m . s e c -1, and a mass density of 2200
sponding to the upper power frequencyfm~, defined in prac- kg.m -3. The numerical model has a width of 4000 m and a
tice as the frequency above which the spectral amplitude of height of 2000 m on the left border. The free surface is tilted
the source becomes less than 5% of the maximum value by 10°, and the source is located exactly on the free surface,
associated with the fundamental frequency fo. For a Ricker at x s = 1720 m (see Fig. 1). Absorbing boundaries are used
wavelet in time, one gets the simple relation fmax = 2.5 f0. on the left, right, and bottom edges of the model. The mesh
In practice, a spatial sampling of approximatively 5 points is composed of 50 X 30 spectral elements, with a polyno-
per minimum wavelength has been found very accurate mial degree of 8 used in each direction, leading to a total of
when working with a polynomial degree of order N = 8. 96,641 points. The line of receivers goes from x = 600 m
Below this value, the solution quickly develops significant to x = 3400 m at the free surface. The variation in time of
numerical oscillations during the propagation. Such an the point-force source is a Ricker wavelet with a central
abrupt transition is characteristic of methods with minimal frequency of 14.5 Hz. The spatial resolution is of the order
numerical dispersion and diffusion. of 5 points per minimum wavelength. In order to compare
Typically, for 2D simulations with a 100,000-point cur- to the analytical solution, a rather restrictive time step has
vilinear grid, the memory occupation is of the order of 30 been used, At = 0.25 m sec, which corresponds to n c =
Megabytes, and the CPU time, for a simulation over 2000 0.25 for this mesh, approximately 40% of the critical Cour-
time steps, is of the order of 15 min on an Ultra Sparc 1 ant number that is determined empirically to be n~ ~ = 0.60.
(140 MHz). For large 3D simulations in an heterogeneous The total number of time steps is 6000.
medium, using a 5,000,000-point curvilinear grid, the mem- A typical snapshot of this simulation (Fig. 1) shows a
ory occupation is of the order of 1.5 Gigabyte, and the CPU strong Rayleigh wave with the expected elliptical polariza-
time, for a simulation over 2000 time steps, is of the order tion and a clear head wave. A very good agreement with the
of 1.5 hr on the CM5 128 nodes of the French Centre Na- exact solution, at receiver 75 located in (x, z) = (2695, 2475)
tional de Calcul Parallble en Sciences de la Terre. In the 3D m and at receiver 100 located in (x, z) = (3400, 2600) m,
case, the memory occupation on the CM5 can be higher is found, with a maximum error of the order of less than 1%,
when using the CMSSL optimized library, which clearly fa- as shown in Figure 2, where the true amplitudes (without
vors CPU efficiency to the detriment of the memory opti- normalization) and the residuals (with a magnification of 10)
mization. For details concerning the parallel implementation between the numerical and analytical solutions are plotted.
of the method, we refer to Fisher (1989), Komatitsch (1997), It can also be checked in this figure that the amplitude of the
and Komatitsch et al. (1997). Rayleigh wave does not decay when it propagates along the
slope, as expected (Garvin, 1956) in plane strain formula-
tion.
Numerical E x a m p l e s
Garvin's Problem
Lamb' s Problem
The second validation example is Garvin's problem
Lamb's problem (Lamb, 1904) is a classical test. The (Garvin, 1956), which involves another type of source than
solution can be regarded as the superposition of three main in Lamb's problem, for example, a line moment-tensor den-
events, the direct P wave, the direct S wave, and the Rayleigh sity that produces a purely explosive source. For a homo-
wave. When the source is located exactly at the free surface, geneous elastic half-space, and a source inside the medium,
the main event is the propagation of a strong Rayleigh wave the main events are a direct P wave, a reflected P wave, and
along the surface, which is nondispersive when the medium a P-to-S mode conversion at the free surface.
is homogeneous and the surface is flat. The model used, the frequency of the source, and all the
Lamb's problem is a challenging test for the accuracy numerical parameters are the same as in Lamb's test, except
of the numerical approximation of the free-surface boundary that the explosive source is now located inside the model at
condition and for the capability of the method to simulate (x, z) = (2236, 1396.5) m. Comparison with the analytical
Rayleigh waves correctly. A detailed numerical analysis of solution, at receiver 55 located in (x, z) = (2307, 2018.5)
this problem can be found in Khun (1985). For the numerical m and at receiver 80 located in (x, z) = (3014.5, 2123) m,
simulation, the free surface has been tilted and the direction is very good with an error than remains less than 1%, as can
of the point force kept normal to it. For this type of mesh, be seen in Figure 3, where the true amplitudes, as well as
RIO(,
--4 ---+ ,.~ ~
--4
X--~ ."4"
Figure 1. Snapshot showing wave propa-
~ "--+ ~-+ gation in tilted Lamb's problem. The line-force
source is located exactly at the free surface.
The cross indicates the source position, and the
----4-------
-ll
t-
--
_
-i-
--t-
diamonds represent the position of the receiv-
ers. The displacement vector field is displayed.
The main event (arrow a) is a strong nondis-
persive Rayleigh wave that propagates along
the free surface with an elliptical polarization.
III The head wave is shown by arrow b and the
{{{ {}{ direct P wave by arrow c.
the residuals between the analytical and numerical solutions, T w o - D i m e n s i o n a l Numerical Simulations
are plotted as in Figure 2. Numerical errors are limited to
More realistic 2D simulations are now considered to
small spurious reflections at the absorbing boundaries.
illustrate the ability of the method to take into account sur-
face topography, deformed interfaces, and complex bound-
Energy-Momentum Conservation and Stability
ary conditions.
of the Absorbing Conditions
To assess some of the properties of the time-marching

algorithm, the propagation of elastic waves in a rectangular Free-Surface and Interface Topography
homogeneous domain, with the same elastic properties as in
Lamb's test, has been simulated. A layered elastic half-space with a prescribed topogra-
The four sides of the model are first assigned to be free phy is considered. The upper layer is characterized by a P-
surfaces, and an explosive source is prescribed within the wave velocity of 2000 m- sec- l, an S-wave velocity of 1300
model. The time step used in this experiment is At = 0.40 m- sec- :, and a mass density of 1000 kg" m -3. The lower-
m sec, corresponding to a Courant number of 0.30, 50% of layer elastic parameters are a P-wave velocity of 2800
the maximal Courant number n~ ~. The time-marching al- m. sec- l, an S-wave velocity of 1473 m. sec- 1, and a mass
gorithm is the predictor-multi-corrector scheme of Appen- density of 1500 kg. m -3. A strong contrast both in velocity
dix B, with a = 1/2, fl = 1/2, and y = 1. The qualitative and in Poisson's ratio is hence modeled, with v = 0.13 for
features of the dynamics of the system are displayed in Fig- the lower layer and v -- 0.38 for the upper layer. The source
ure 4. The long-term stability is apparent from these results. is explosive and located inside the upper layer, just below
After l0 s time steps, the total energy is perfectly conserved the free surface, so as to excite a strong Rayleigh wave along
without oscillations, and correct conversions between kinetic the curved topography.
and potential energy are observed. The numerical model has a width of 2500 m and an
Using the same model, but now enforcing first-order average height of 1700 m. The line of receivers goes from
radiation conditions on the four sides, the efficiency of the x = 7 0 0 m t o x = 2 1 0 0 m a t z = 1380m. Thesource
absorbing boundary conditions can be assessed. For exact position is (x, z) = (620, 1853) m. The mesh is composed
absorbing boundaries, the total energy is expected to decay of 48 × 40 elements, with a polynomial approximation of
monotonically to zero as the spherical waves propagate out- order N = 8, and a total number of collocation grid points
side the domain. The total energy during the simulation is of 123,585. The explosive source is a Ricker wavelet in time,
shown in Figure 5 to decay as expected with no oscillations. with a central frequency of 15 Hz. The number of points per
A small residual energy, typically of the order of 10 - 4 of minimal wavelength is of the order of 5. The time step is At
the total energy of the source, remains in the medium and is -- 0.30 m sec, which corresponds to a Courant number of
related to spurious reflected waves of small displacement 0.45, 70% of the maximal Courant number n~ ax. The simu-
amplitude (a few percents of the amplitude of the outgoing lation is done over 5000 time steps.
waves) produced by the low-order-absorbing boundary con- A typical snapshot of the simulation is shown in Figure
ditions. The spurious reflections are themselves rapidly ab- 6, and the seismograms recorded by the receivers are shown
sorbed, and the energy finally drops to zero. in Figure 7. A strong Rayleigh wave can be observed prop-
Ux residuals trace 75 Uz residuals trace 75

le-05 2.5e-05
Numerical results Numerical results - -
Residuals * 10 ....... 2e-05 Residuals * 10 .......
5e-06
1.5e-05
1e-05
0
g g 5e-06
m
oJ
"0
Q.
-5e-06
"O
0 #p
E
<
E -5e-06
- 1e-05
- 1e-05
-1.5e-05
-1.5e-05
-2e-05
-2e-05 i i i I I f i
-2.5e-05 i i f i i i r
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4
Time (s) Time (s)

1e-05 2.5e-05
2e-05
5e-06
1.5e-05
L~ 1e-05
0
g ~ 5e-06
03 II1
"O
.._= -5e-06 " 0
<E
~" -5e-06
<
-le-05
- 1e-05
-1.5e-05
-1.5e-05
Numerical results ~ Numerical results
Residuals * 10 ....... -2e-05 Residuals * 10 .......
-2e-05 i i i i i
-2.5e-05 i i i i i i i
0 0.2 0.4 0.6 0.8 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4
Time (s) Time (s)
Figure 2. Time evolution of the horizontal and vertical components of the displace-
ment at receivers 75 and 100, along the free surface, for tilted Lamb's problem with a
line-force source as in Figure 1. The amplitude of the Rayleigh wave remains constant
along the free surface, as expected for plane strain. The residuals with respect to the
analytical solution are also displayed, with an amplification factor of 10, on the same
plot. A very good agreement is found, and a minimal dispersion is observed. The
maximum relative error remains less than 1%.
agating along the topography, with a Rayleigh-to-body- mean height of the topography is of the order of 1300 m.
wave mode conversion, particularly clear on the snapshot The source position is (x, z) = (1000, 670) m. The mesh is
after the S wave converted at the free surface, when the composed of 60 X 12 elements, with a polynomial approx-
Rayleigh wave encounters a change of curvature at the free imation of order N = 8 and a total number of collocation
surface. Such a mode conversion is in agreement with the points of 46,657. The explosive source is a Ricker wavelet
theoretical study of Rulf (1969) and previous numerical in time, with a central frequency of 12 Hz. The number of
studies by Jih et al. (1988) and Komatitsch et al. (1996). points per minimal wavelength is of the order of 5. The time
The second model includes a real topography derived step is At = 0.30 m sec, which corresponds to n c = 0.45.
from a structure in the Peruvian Andes (courtesy of Elf Aq- The number of time steps in this simulation is 5000. The
uitaine) as shown in Figure 8 with no vertical exaggeration. line of receivers is placed at the surface between x = 900
We consider a simple homogeneous elastic half-space. The m and x = 5000 m.
elastic parameters are a P-wave velocity of 3200 m- s e c - 1, A typical snapshot of the simulation is shown in Figure
an S-wave velocity of 1848 m . sec-1, and a mass density of 8, and the recorded seismograms are shown in Figure 9. For
2200 kg- m - 3. The explosive source is located near the free a flat free surface, the seismograms should be characterized
surface inside the model. by a single arrival. As expected for a foothill zone, topo-
The numerical model has a width of 5500 m, and the graphical effects are of major importance and induce, in par-

3e-05 8e-05
Numerical results Numerical results - -
Residuals * t0 ....... 6e-05 Residuals * 10 .......
2e-05
1e-05
,¢ 4e-05
2e-05
g 0
0
"O
-2e-05
O_ - 1 e-05
E -4e-05
-2e-05 -6e-05
-8e-05
-3e-05
-0,0001
-4e-05 i i a r i i
-0.00012 I I I I i i i
0 0.2 0.4 0.6 0.8 t 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4
Time (s) Time (s)
U x residuals trace 80 Uz residuals trace 80

6e-05 4e-05
Numerical results - - Numerical results
Residuals * 10 ....... 3e-05 Residuals * 10 .......
4e-05
2e-05
2e-05
~ le-05
0~
"10
R
0
"10
0
-1 e-05
J V
-2e-05
E E
<
< -2e-05
-4e-05
-3e-05
-6e-05
-4e-05
i i i
-8e-05 i i i i i I
-5e-05 i i
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4
Time (s) Time (s)
Figure 3. Time evolution of the horizontal and vertical components of the displace-
ment at receivers 55 and 80 for tilted Garvin's problem with a compressive volumetric
source. The residuals with respect to the analytical solution are also displayed, with an
amplification factor of 10, on the same plot. A very good agreement is observed, and
the maximum error remains less than 1%. A small spurious reflection from the ab-
sorbing boundaries can be observed and is indicated by arrow a.
ticular, some strong diffracted phases and conversions from free surface. The elastic parameters are a P-wave velocity
Rayleigh to body waves. of 3200 m . sec-1, an S-wave velocity of 1847.5 m . sec-1,
and a mass density of 2200 kg- m - 3. The numerical model
Step at the Free Surface has a width of 4000 m and a height of 2400 m. To simulate
The problem considered is the transmission and reflec- a quarter-space, absorbing boundaries are used on the left
tion of a Rayleigh wave at a step of the free surface, for and bottom sides of the mesh. The source position is (x, z)
example, the so-called "comer problem," and can be related, = (2440, 2340) m. The normal and tangential components
for instance, to a cliff or a fault at the free surface (see, for of the displacement are recorded by a line of receivers placed
example, Lapwood, 1961; Hudson and Knopoff, 1964; de around the comer on the free surface, between (x, z) =
Bremaecker, 1958). This problem is very sensitive to the (2600, 2400) m and (x, z) = (4000, 1000) m. The mesh is
accuracy of the free-surface boundary condition, and ad hoc composed of 50 X 30 elements, with a polynomial approx-
numerical treatments of the free-surface boundary condition imation of order N = 8. The total number of collocation
at the comer are sometimes introduced that may perturb the points is 96,641. The line-force source is a Ricker wavelet
relative amplitudes of the different phases. in time with a central frequency of 14.5 Hz. The number of
The physical problem is a homogeneous elastic quarter- points per minimum wavelength is of the order of 5. The
space. The source is a line force located just below the upper time step is At = 0.50 m sec, and the simulation is done
4e+06 over 5000 time steps. The normal and tangential components
Total energy ........
3.5e+06 Potential energy of the displacement are recorded by the receivers on both
surfaces.
3e+06 The time response of the displacement is shown in Fig-
~..~2.5e+06
ure 10. A strong Rayleigh wave propagating along the upper
side is observed. This Rayleigh wave is partly reflected and
2e+06 partly transmitted at the comer. The transmitted Rayleigh
w 1.5e+06 wave can be seen traveling down the right side of the model,
and the comer acts as a diffraction point, triggering spherical
le+06 P and S waves and small Rayleigh waves.
500000
0 i I i i i i i J
5 10 15 20 25 30 35 40 Three-Dimensional Numerical Simulations

Time (s)
To illustrate the capabilities of the spectral element
Figure 4. Long-term stability of the energy-mo- method for 3D problems, a simple parametric study of the
mentum conserving time scheme for the case of a
diffraction of elastic waves by a local surface topography
homogeneous elastic medium, bounded by free-sur-
face boundaries, with an explosive source inside the has been carried out. This study is motivated by the recent
volume. The time evolution of the total and potential results of Bouchon et al. (1996) and Bouchon and Barker
energy is displayed for 105 time steps and a step size (1996) on the effects of a small 3D hill and their implication
of At = 0.40 m sec. The total energy is shown to for strong ground motion, as, for example, the strong accel-
remain constant during the simulation.
erations, of the order of 1.82 g, recently recorded on top of
the Tarzana hill during the 1994 Northridge, California,
earthquake (Spudich et aL, 1996).
Surface topography has long been recognized as re-
sponsible for energy diffraction and scattering effects that
are associated with the observations of late arrivals, mostly
Rayleigh waves [see, for example, the early study of Gilbert
and Knopoff (1960)], and these effects are considered today
3.5e+06 i
as an important source of seismic coda (Clouser and Lang-
Total Energy ........ ston, 1995). Numerous field experiments, performed on hills
3e+06 Potential energy - -
and mountain ridges, have also shown that surface topog-
..........................."........
raphy can induce strong amplifications and variations of the
2.5e+06 ground motion during earthquakes as a result of subsurface
focusing and scattering of the radiated energy (Umeda et al.,
2e+06 ...
1986; Carver et al., 1990; S~inchez-Sesma, 1997). Recently,
1.5e+06
i these effects were evidenced by the well-instrumented dra-
matic examples of the 1994 Northridge, California, earth-
le+06 quake (Spudich et al., 1996; Gao et al., 1996) and the 1995
Kobe, Japan, earthquake (Told et al., 1995; Iwata et al.,
500000 1996), both with epicenters inside urban areas.
Although the need to include the influence of local site
0 i i t
0.2 0.3 0.4 0.5 0.6 0.7 topography in seismic risk assessment and microzonation
Time (s) studies is now well recognized (Bard, 1982; Gaffet and Bou-
chon, •989; S~nchez-Sesma and Campillo, •993; Pedersen
Figure 5. Energy evolution for the case of a ho- et aL, 1994), the calculated levels of amplification are very
mogeneous medium, all bounded by absorbing
often underestimated. Besides possible effects due to the
boundaries, with an explosive source. The total en-
ergy decreases rapidly as the energy is radiated out- nonlinear behavior of soils, or variations in shear-wave ve-
side the domain. Small spurious reflections may ap- locity, comparisons of the observations of Spudich et al.
pear for waves impinging on the sides at small (1996) with theoretical and numerical results indicate that
incidence due to the zeroth-order approximation used the three dimensionality of the topography and its internal
here. These spurious reflections, that keep a small
fraction of the total energy in the system (less than a structure are important factors. In the case of Tarzana hill,
few percent), disappear as soon as they reach a new transverse oscillations of the 3D elongated structure were
absorbing boundary. This effect is hardly visible on found to be of major importance (Spudich et al., 1996; Bou-
this figure due to the scale and can only be observed chon and Barker, 1996), both in the observations and in the
when zooming the area indicated by arrow a. simulations. Surprisingly, very few theoretical investigations
a/
--- iiL
01- . . . .
Figure 6. Snapshot showing wave propa-

gation in the case of a layered elastic half-space
with curved free surface and interface. The
source is explosive and is located within the
volume just below the free surface as indicated
by the cross. The dotted line represents the po-
sition of the receivers. Displayed is the dis-
placement vector field. A Rayleigb wave that
propagates along the topography is observed
(arrow a), with a strong Rayleigh-to-body-
N~
wave mode conversion (arrow b) just after the
converted S wave at the free surface (arrow c),
that occurs when the Rayleigh wave encounters
a change of curvature along the free surface.
x (m) x (m)
1000 1500 2000 1000 1500 2000
O,
t: 01
Ux component Uz component
Figure 7. Time responses of the horizontal and vertical components of the displace-
ment, at the receivers placed inside the model, in the case of the layered elastic half-
space of Figure 6. The strong later event (arrow a), traveling with an S-wave velocity,
corresponds to a Rayleigh-to-body-wave mode conversion when the Rayleigh wave
encounters a strong change of curvature along the free surface. Arrivals at around t =
1.2 sec are reflections from the curved interface, with also some spurious reflections
from the absorbing boundaries.
Amplification by a Three-Dimensional Hill

of 3D surface irregularities have been reported (S~inchez- Following Bouchon et al. (1996), the topography is a
Sesma, 1983; S~nchez-Sesma and Luz6n, 1996; Bouchon et gentle 3D hill and is simply parameterized by
al., 1996; Bouchon and Barker, 1996).
In the following, the spectral-element method is shown z(x, y) = Zo +
to be an effective tool for investigating the diffraction of
( (x - Xo)2~ _ ---Y°)Z~ (27)
elastic waves by 3D structures, in the case of both a homo- ho exp /'
geneous and a stratified medium. l
Figure 8. Snapshots showing wave propa-

gation for a homogeneous elastic section
across the Peruvian Andes (courtesy Elf-Aq-
uitaine). Displayed are the displacement vector
field and the real topography with no vertical
exaggeration. The length of the section is 5500
m, and the average height of the topography is
of the order of 1300 m. The dotted line repre-
sents the position of the receivers, and the
source position is indicated by a cross. Even
for such a simple homogeneous model, the ef-
fect of the strong topography is very important,
with diffracted phases that are superimposed to
the direct wave. Rayleigh-to-body-wave con-
versions can also be seen.
x(m) x (m)
1000 2000 3000 4000 1000 2000 3000 4000
f.-
IIIIJJll]l I
Figure 9. Time responses of the horizontal and vertical components of the displace-
l
LflIJIrlll
ment at the receivers placed along the free surface in the case of the Peruvian cross
section of Figure 8. The section is supposed to be elastic and homogeneous. Diffraction
from the strong topography is particularly clear on both components. Without topog-
raphy, only a single arrival would be observed.
with z0 = 1050 m, Xo = 1040 m, Y0 = 1040 m, ax = 250 Case of a Homogeneous Medium. A homogeneous elastic
m, cry = 125 m, and h 0 = 180 m. The size of the model is medium is first considered, as in Bouchon et aL (1996), with
2080 X 2080 X 1050 m, and the mesh is composed of 26 a P-wave velocity of 3200 m - s e c - 1, an S-wave velocity of
x 26 X 14 elements, with a polynomial order of N = 8, 1847.5 m - sec - 1, and a density of 2200 kg. m - 3. The struc-
corresponding to a total number of points of 4,935,953 (see ture is excited by a vertically incident plane S wave that is
Fig. l l a ) . polarized along either the minor or the major axis of the
x(m) x(m)
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
I--
Figure 10. Time responses of the horizontal and vertical components of the dis-
placement vector at the receivers placed around the corner for the model consisting of
a step at the free surface. A strong Rayleigh wave propagating along the upper free
surface can be observed (arrow a). This surface wave is partly reflected (arrow b) and
partly transmitted (arrow c) by the corner that acts as a diffraction point, triggering P
and S waves (arrows d and e, respectively) with typical diffraction hyperbola.
(a) (b) the lateral edges of the grid to simulate a plane wave in an
qll
infinite half-space.
For a vertically incident S wave of wavelength 2s = h
polarized along the minor axis, Figure 12 shows the surface
ground motion recorded along the minor axis of the topog-
raphy. The main features are a strong amplification of the
amplitude of motion, on the hill, and very clear diffracted
waves, away from the topography, composed of a surface P
wave and a Rayleigh wave. In the vicinity of the hill, the
Figure 11. Three-dimensional models: (a) a 3D amplitude of the Rayleigh wave is of the same order as that
Gaussian shape topography is considered in the case of the surface P wave, but, as it propagates, the P-wave
of an homogeneous elastic half-space and (b) a strat- amplitude decreases faster, with a 1/r decay for the P wave
ified elastic half-space with an nonsymmetrical basin. versus a 1/~/~ decay for the Rayleigh wave, and the Rayleigh
The topography is characterized by a height of h 0 =
180 m and a standard deviation of 250 and 125 m wave becomes the dominant feature. A strong directivity
along the x and y axis, respectively. The size of the effect induced by the topography is observed when compar-
model is 2080 × 2080 × 1050 m. The mesh is com- ing the time response of the y component of the displacement
posed of 26 × 26 × 14 elements, with a polynomial recorded along the minor and major axes, as shown in Figure
order of N = 8 in each direction. The total number 13: The diffracted surface wave clearly propagates prefer-
of collocation points is 4,935,953.
entially along the minor axis.
The response of the elastic structure for two vertically
incident S waves with 2s = h and 2~ = 3h, polarized along
topography. The source is a Ricker wavelet in time, of max- the minor or the major axis of the topography, is shown in
imal amplitude 0.5, and two fundamental wavelengths are Figure 14. The time evolution of the horizontal displacement
investigated: 2s = h or ,~s = 3h, where h is the maximum amplitude recorded at the surface is plotted for the two in-
height of the hill. The corresponding fundamental frequen- cident wavelengths, 2, = h and 2s = 3h, in the case of an
cies of the Ricker arefh = 10.2 Hz andf3h = 3.4 Hz. The incident Ricker wavelet having a maximal amplitude of 0.5.
simulated timescale spans 0.8 sec, with a time step of At = Each frame of these movies is numbered in the upper left
0.5 m sec. Periodic boundary conditions are implemented on comer. This number indicates the time elapsed from the be-
The Spectral Element Method." An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures 38 !
Uy component along Y Uz component along Y
2000 2000-
1500 1500-
~" 1000- ~ 1000-
500. 500-
d 012 014 016

° 0I 012 014 016
Time (s) Time (s)
Figure ]2. Time responses of the y (horizontal) and z (vertical) components of the
displacement vector at the receivers placed on the free surface, along the minor axis
(y-direction) of the topography, for the 3D homogeneous model of Figure 1 l(a) with
a vertical incident S wave of wavelength 2s = h, polarized along the minor axis. The
direct S wave is clearly seen on the y horizontal component (arrow a), with two waves
diffracted by the topography: a 1/r decaying P wave (arrow b) and a Rayleigh wave
(arrow e). Small artificial phases (arrow d) are observed near the edges of the model
due to the lateral periodic boundary conditions used in the simulation.
Uy component along Y Uy component along X
2000
1500-
x I000- x
500 -
0-
6 012 0'.4 016 0 012 014 016

Time (s) Time (s)
Figure 13. Time responses of the y (horizontal) component of the displacement

vector at the receivers placed on the free surface, along the minor axis (y direction)
and the major axis (x direction) of the topography for the 3D homogeneous model. The
vertical incident S wave, of wavelength 25 = h, is polarized along the minor axis. The
direct S wave can be seen on both directions. A strong directivity effect, due to the 3D
geometry of the topography, is observed on the diffracted P (arrow b) and Rayleigh
(arrow c) waves that appear mainly on the recordings done along the minor axis.
Figure 14. Snapshots showing the amplitude of the horizontal projection of the displacement
vector field at the surface for an incident S wave, with a wavelength ,~s = h or )~ = 3h, polarized
along the minor or the major axis, The slides are numbered in the upper left corner, and the time
step between two consecutive numbers is 20 m sec. A strong directivity in the wave propagation
along the topography can be observed. The maximum amplification, of the order of two, is observed
for the case of an incident wavelength )~ = h and a polarization along the minor axis.
ginning of the simulation, in units of 20 m sec. The first few wave with a nice retrograde elliptical polarization is ob-
frames, corresponding to the arrival of the direct wave, have served as well as a strong diffracted P wave with a small
been suppressed. polarization effect due to the curvature of the topography.
A strong directivity in the propagation of the diffracted In contrast, along the major axis, a very small, and rapidly
energy can again be observed: The surface wave propagates attenuating, diffracted S wave can be observed. With the
preferentially along the minor axis of the topography inde- same wavelength but incident polarization along the major
pendently of the polarization of the incident wave, while the axis (see Fig. 14), the amplification gets weaker, of the order
surface P wave propagates in the direction of polarization. of 1.2 on frame 15, and while a wave is still propagating at
The relative amount of diffraction at the top and bottom of the surface along the minor axis, the diffracted P wave now
the hill is difficult to quantify for these wavelengths, even if propagates along the major axis (see frames 26 and 27). This
most of the diffracted energy appears to be generated at the can also be analyzed on the polarization diagrams of Figure
top of the hill. The maximum amplification is of a factor 2 16. Along the minor axis, besides the direct wave, no Ray-
and is recorded at the summit for 2s = h and for an incident leigh wave can be observed. The diffracted energy is com-
polarization along the minor axis (frame 15). In this case, posed of a small body wave that is interpreted as an S wave,
the diffracted surface P and the Rayleigh waves are gener- with an 1/r attenuation, propagating at the surface. Along
ated preferentially along the minor axis (frames 26 and 27). the major axis, the diffracted energy is mainly a P wave. At
This can be clearly seen on the polarization diagram re- the summit, a strong amplification is observed with a nice
corded along the minor axis and the major axis (see Fig. 15). polarization due to the curvature of the topography.
Along the minor axis, besides the direct wave, a Rayleigh In contrast, for a larger incident wavelength, 2 s = 3h,
(a) (b)
Figure 15. Polarization diagrams at differ-

ent receivers along (a) the minor axis (y direc-
tion), and (b) the major axis (x direction) of the
topography for the 3D homogeneous model
and a vertically incident S wave, of wavelength
2s -= h, polarized along the minor axis. Dis-
played is the projection of the displacement
vector onto the plane defined by the y and z
directions. The upper display corresponds to a
receiver located on top of the topography,
while the others are for receivers located at 80,
160, 240, 320, and 400 m from the summit.
Along the minor axis, a strong diffraction (ar-
row a) is observed with, besides the direct S
wave, a clear Rayleigh wave that propagates
down the hill with an elliptical polarization (ar-
row b) and a diffracted P wave with small po-
larization due to the curvature effect of the to-
pography. Along the major axis, very small
diffracted energy is produced, and only a small
surface S wave (arrow c) can be observed.
(a) (b)
Figure 16. Polarization diagrams at differ-

ent receivers along (a) the minor axis (y direc-
tion), and (b) the major axis (x direction) of the
topography for the 3D homogeneous model
and a vertically incident S wave, of wavelength
2, = h, polarized along the major axis. Dis-
played is the projection of the displacement
vector onto the plane defined by the x and z
directions. The positions of the receivers are
the same as in Figure 15. Along the minor axis,
besides the direct S wave (arrow a), no clear
diffracted P wave can now be observed, while
a diffracted S wave (arrow b) still propagates
at the surface along that direction in contrast
with Figure 15. A diffracted P wave (arrow c)
now propagates along the major axis but is of
smaller amplitude than in the case of Figure
15.
a strong amplification, of the order of 1.7, is now observed fication events, like, for example, in the oscillation that ap-
at the summit for an incident polarization along the major pears in the case of 2 s = 3h with a polarization along the
axis (see frame 20 in Fig. 14). For such a large incident major axis.
wavelength, the diffracted waves are difficult to identify (see The transfer function of the structure is shown in Figure
frame 26). An interesting breathing phenomena can be ob- 17. It is computed as the ratio between the spectral amplitude
served over a period of 80 m sec between frames 20 and 25, of the component of the displacement colinear to the polar-
the strong amplification on frame 20 is dimmed (frame 21) ization of the incident wave, and the spectral amplitude of
and disappears on frames 22 and 24 to reappear on frames the incident wave. In Figure 17, the ratio has been calculated
25 and 26. No unique interpretation has yet been obtained. using the results of the simulation with a Ricker wavelet of
It may be related to an oscillation of the structure, for ex- fundamental wavelength '~s = h.
ample, the hill coupled to the elastic half-space, or to a res- For a polarization along the minor axis, a strong ampli-
idence time for constructive phase interference of the dif- fication is observed at the top of the hill for high frequencies
fracted waves within the hill structure. (15 to 30 Hz). At low frequency, the central mode of the
The amplification pattern on the hill exhibits a strong structure can be observed, with a resonance effect. For about
variability. The maximal amplification can occur in a small 5 Hz, two maxima of the amplification are observed at mid-
area around the summit, or at mid-slope, while the summital slope, while a strong deamplification occurs at the top of the
area is exhibiting a deamplification; see frame 16 for 2s = hill. These two maxima can be associated with a kind of
3h with an incident polarization along the minor axis. A resonance due to the beginning of constructive interferences
given area may be submitted several times to strong ampli- inside the hill structure. The transfer function along the mi-
1300 1300 ,
1200 ~ Malor a Is 1200 Minor axis
1100 1100
1000 ' ' ' , I 1000 I , , , ,
400 800 1200 1600 400 800 1200 1600
5'
~,10 ~10
= ="
~15' ~t5
~'20, ~'~
25 25
40O 800 12OO 16OO 400 800 1200 1600

Offset (m) Offset (m)
• c~'°a~ ~ • ~ ~,
(b)
5-
~,10 ~,10-
z
o>'15-
g g
~20 20-
25 25-
30
400 800 1200 1600 30- 4~o ~o 12'oo- 16'oo
Offset (m) Offset (m)
(c) (d)
Figure 17. Transfer function for the response of the 3D homogeneous model to a
vertically incident S wave of wavelength 2, = h: (a) response along the major axis for
an incident polarization along the minor axis; (b) response along the minor axis for an
incident polarization along the minor axis; (c) response along the major axis for an
incident polarization along the major axis; and (d) response along the minor axis for
an incident polarization along the major axis. The transfer function is computed as the
ratio between the spectral amplitude of the component of the displacement vector along
the direction of polarization and the spectral amplitude of the incident wave. A strong
variation of the amplification pattern can be observed.
nor axis is clearly modulated by the Rayleigh and surface P lower layer is characterized by a P - w a v e velocity of 3200
waves that travel downward. m- s e c - 1, an S-wave velocity of 1847.5 m - s e c - 1, and a den-
For a polarization along the major axis, a broad ampli- sity of 2200 kg- m - 3 as before. The upper layer now has a
fication occurs at the top of the hill for low frequencies, up P-wave velocity of 2000 m . s e c - 1, an S-wave velocity of
to 10 Hz, while at higher frequencies (15 to 30 Hz), a small 1155 m . sec -1, and a density of 1500 k g . m -3. The free-
deamplification is now observed. Due to the directivity ef- surface topography and the dimensions of the model remain
fect, modulation by the surface wave propagating down the the same as in the homogeneous case. The geometry of the
minor axis can be observed, while the modulation along the lower surface of the sedimentary basin is given by
major axis is now due to the diffracted P-wave propagation.
The transfer function computed with the results of the z(x, y) = Zb -
simulation done for 2s = h can be convolved with the source
function for 2 s = 3h in order to generate the seismograms h b exp ( - (x - x0)21 - y0)2/ (28)
along the minor and major axes. Comparison of the synthetic
J exp (- -Cr7 -J'
seismograms obtained by convolution to the simulated ones
for 2~ = 3h allows one to check the accuracy, or at least the with zb = 825 m, h b = 200 m, x o = 900 m, Y0 = 1040 m,
consistency, of the different simulations. The result is shown and~r x = 2 4 6 m , i f x = < x 0 ; c r x = 2 3 7 m , i f x > x 0 ; % =
in Figure 18 for an incident polarization along the minor 320 m, if y _-< Y0; and Cry = 246 m, if y > Y0. The basin is
axis and at two receivers along the minor axis only: on top not centered with respect to the free-surface topography, and
and at the base of the topography. The agreement is very the edges of this basin themselves are nonsymmetrical (see
good and attests the consistency of the method over a wide Fig. 1 lb).
range of frequencies. The wiggle observed at 0.8 sec on the Due to the contrast in velocity and Poisson's ratio be-
records is due to the fact that the actual simulated time was tween the upper and lower layers of the model, this simu-
exactly 0.8 sec. lation clearly requires more computer resources than in the
case of the homogeneous example. In particular, the simu-
Case of an Asymmetrical Sedimentary Basin in the Presence lated timescale has to be substantially increased in order to
of Topography. A two-layer elastic half-space having the capture the response of the structure. The results presented
same surface topography has also been considered. The here are just a preliminary attempt, and extended simulations
3 3
Calculation Calculation
Convolution ........ 2 Convolution ........
2
1 1
g
0 e- 0
0
E E
o~
-1 -1
.~_
a -2 a -2
-3 -3
Receiver # 105 Receiver # 65

-4 -4
i i I i i i i
0.15 0.4~ 0.65 0.90 0.15 0.40 0.65 0.90

Time (s) Time (s)
(b)
Figure 18. Comparison between the time response of the y component of the dis-
placement vector, at two receivers, obtained from a direct simulation and from the
convolution of the transfer function by the source function of the simulation. The
receivers are (a) at the top of the hill and (b) at the base of the hill along the minor
axis. The simulation was performed for a vertically incident S wave polarized along
the minor axis with 2s = 3h. The transfer function is the one of Figure 17a computed
from the results of the simulation with an incident wavelength of 2, = h. A very good
agreement is observed up to t = 0.8 sec, which is exactly the time at which the
numerical simulation ended.
Uy component along X Uy component along Y
2000- 2000-
1500 1500
× 1000. x 1000
500 500-
0 0,2 0.4 0.6 0.8 1.0 0 0,2 0.4 0.6 0.8 1.0
Time (s) Time (s)
(~) (b)
Uz component along X Uz component along Y
2000 2000
1500 1500
-~ lOOO i ooo
500 500 -
0 0 ~
0 0.2 0.4 0.6 0.8 1.0 0 012 014 016 018 110
Time (s) Time (s)
(c) (d)
Figure 19. Time responses of the y (horizontal) and z (vertical) components of the
displacement vector at the receivers placed on the free surface, along the major axis (x
direction) and the minor axis (y direction) of the topography, for the 3D heterogeneous
model of Figure 11 a0 with a vertically incident S wave of wavelength 2s = 2h polarized
along the minor axis. The direct S wave is clearly identified on the y component of the
displacement (arrow a), with a deflection due to the shape of the basin. The first multiple
can be observed (arrow b) on the y component in both directions and produces a distortion
of the diffracted Rayleigh (arrow d) and P waves that is particularly clear along the
minor axis. For the vertical z component of the displacement, the amplitudes along the
major axis have been saturated to a third of the maximal amplitude in order to extract
the converted S to P wave at the bottom of the basin (arrow c). Along the minor axis,
the Rayleigh wave has a vertical component that can be clearly seen (arrow d).
are still required. They are, however, interesting enough to nor and the major axes of the topography are shown in Fig-
illustrate the application of the spectral-element method to ure 19.
realistic 3D structures. The main features, for this incident wavelength, are the
The simulation is performed for an incident plane S result of the superposition of two weakly interacting effects:
wave of wavelength 2~ = 2h, polarized along the minor axis. the shape of the surface topography and the shape of the
The seismograms recorded at the free surface along the mi- sedimentary basin. The effects of the basin structure remain
limited in this simulation. The study of more complex ex- length. Moreover, long-term energy-conserving and stability
amples would have required larger models (mainly in terms properties of the method have been shown.
of computer memory) that were beyond our computer re- The capabilities of the method to handle complex free-
sources. surface geometries and deformed internal interfaces have
In Figures 19a and 19b, the distortion of the direct wave been illustrated by solving realistic 2D problems: One in-
due to the structure of the basin can be seen as well as the volves a step at the free surface, while the other includes a
first multiple, which is quite strong here due to the high realistic geological topography based on a cross section
reflection coefficient, of the order of 0.4, that has been as- along the Peruvian Andes. In both cases, Rayleigh waves
sumed in this simulation. In particular, in Figure 19b, the and complex surface-to-body-wave conversions have been
distortion of the diffracted P and Rayleigh waves by the first accurately modeled. The method provides a very flexible
S-wave multiple is quite clear. Due to the limited timescale tool to understand and extract, at low computational cost,
of the simulation, only the first multiple has been simulated. quantitative physical information from complicated wave
In Figure 19c, the amplitudes have been saturated at about phenomena such as diffraction, conversion, and generation
a third of the maximal amplitude in order to extract the con- of Rayleigh or interface waves that occur in geophysical
verted S to P wave at the bottom of the basin that can be applications.
seen here as the first arrival. In Figure 19d, the vertical com- Finally, the spectral method is shown to be an efficient
ponent of a clear Rayleigh wave diffracted by the topogra- tool for studying the diffraction of elastic waves by 3D sur-
phy can be observed and must be considered as the comple- face topographies and its effect on strong ground motion.
ment of the horizontal component observed in Figure 19b, Complex amplification patterns, in space and time, are
in order to give the classical elliptical polarization. These shown to occur even for a gentle 3D hill. The results ob-
limited responses already show some of the interesting phe- tained in this study are in very good agreement with those
nomena produced by the complicated geological structure: obtained by Bouchon et al. (1996) using a boundary integral
smaller Rayleigh-wave velocity in the case of the sedimen- method. The method allows one to handle a heterogeneous
tary layer, presence of an S-wave multiple created by reflec- internal structure below the topography, which leads to in-
tion at the bottom of the basin, and small S-to-P mode con- teresting geophysical applications for seismic risk assess-
version generated on the edges of the basin. The arrival time ment and microzonation studies.
and the waveform of the first multiple, created in the sedi- The method can be efficiently implemented on distrib-
mentary layer, have been checked and shown to be consis- uted memory parallel machines. The typical CPU time for
tent with the thickness and the velocities of this layer. an average 2D simulation, using a mesh of the order of
100,000 points, and simulating 2000 time steps, is 4 rain on
64 nodes of a CM5. The largest 3D simulation treated in this
Conclusions article involves a 26 × 26 × 14 elements mesh, with a
polynomial approximation order of N = 8, leading to a
A practical spectral-element method for calculating the 5,000,000-point curvilinear grid. Such a simulation, with 64
propagation of seismic waves through 2D and 3D geological bits computation and 2000 time steps, requires 1.5 hr on 128
structures has been presented. The method is based on a nodes of a CM5. This will allow real-time visualization and
high-order variational formulation of elastodynamics that al- interactive physical analysis of amplification phenomena
lows the natural treatment of an irregular free surface. The and seismic risk assessment using modern distributed mem-
procedure proposed in this article is based upon semi-dis- ory parallel architectures.
cretization: For the spatial discretization, the spectral-ele-
ment approximation is used, producing a system of ordinary Acknowledgments
differential equations in time, which in turn is discretized
using a finite-difference method for ordinary differential The authors are very grateful to F. J. Sfinchez-Sesma for numerous dis-
cussions of the 3D results and for providing them with Garvin's solution.
equations. More specifically, an energy-momentum conserv-
They would also like to thank Y. Maday and R. Madariaga for fruitful
ing algorithm in time, which can be put into a classical pre- discussions all along this work. Stimulating discussions with G. Seriani and
dictor-multi-corrector format, has been used. The spectral- E. Priolo are also acknowledged. Many thanks to C. Caquineau and P.
element method is shown on various examples to combine Stoclet for their help in the implementation of the code on the CM5. G.
the geometrical flexibility of a low-order finite-element Moguilny provided an invaluable support for the 3D visualizations. The
constructive remarks of the reviewers T. Ohminato and L.-J. Huang are
method with the rapid convergence rate associated with
also acknowledged. This work has been partly supported by the French
spectral techniques, even when dealing with deformed ge- Centre National de Calcul Parall~le en Sciences de la Terre (CNCPST).
ometries or heterogeneous elastic properties.
Classical 2D problems, Lamb and Garvin, for which
analytical solutions exist, have been studied to assess the References
accuracy of the method. The discrete solution is shown to Aki, K. and K. L. Lamer (1970). Surface motion of a layered medium
present minimal numerical dispersion and diffusion. A high having an irregular interface, due to incident plane SH waves, J. Geo-
accuracy is obtained using only 5 points per minimal wave- phys. Res. 75, 933-954.
The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response o f 2D and 3D Geological Structures 389
Alford, R., K. Kelly, and D. Boore (1974). Accuracy of finite difference Fisher, P. (1990). Analysis and application of a parallel spectral element
modeling of the acoustic wave equation, Geophysics 39, 834-842. method for the solution of the Navier-Stokes equations, Comp. Meth.
Babu~ka, I. and M. R. Dorr (1981). Error estimates for the combined h and AppL Mech. Eng. 80, 483--491.
p version of the finite element method, Numer. Math. 37, 257-277. Fomberg, B. (1988). The pseudospectral method: accurate representation
Babu~ka, I., B. A. Szab6, and I. N. Katz (1981). The p version of the finite of interfaces in elastic wave calculations, Geophysics 53, 625-637.
element method, SIAM J. Numer. Anal. 18, 512-545. Frankel, A. (1993). Three-dimensional simulations of ground motions in
Bard, P. Y. (1982). Diffracted waves and displacement field over two di- the San Bernardino valley, California, for hypothetical earthquakes
mensional topographies, Geophys. J. R. Astr. Soc. 71, 731-760. on the San Andreas fault, Bull. Seism. Soc. Am. 83, 1020-1041.
Bayliss, A. and E. Turkel (1980). Radiation boundary conditions for wave- Gaffet, S. and M. Bouchon (1989). Effects of two-dimensional topographies
like equations, Comm. Pure AppL Math. 33, 707-725. using the discrete wavenumber-boundary integral equation method in
Bayliss, A., K. E. Jordan, B. J. LeMesurier, and E. Turkel (1986). A fourth- P-SV cases, J. Acoust. Soc. Am. 85, 2277-2283.
order accurate finite-difference scheme for the computation of elastic Gao, S., H. Liu, P. M. Davis, and L. Knopoff (1996). Localized amplifi-
waves, Bull. Seism. Soc. Am. 76, 1115-1132. cation of seismic waves and correlation with damage due to the North-
Bernardi, C. and Y. Maday (1992). Approximations Spectrales de Probl~- ridge earthquake: evidence for focusing in Santa Monica, Bull Seism.
rues aux Limites Elliptiques, Springer-Verlag, Paris. Soc. Am. 86, $209-$230.
Bonnet, M. (1995). l~quations Intggrales et l~l~ments de Frontibre, Paris, Garvin, W. W. (1956). Exact transient solution of the buried line source
CNRS t~ditions, Eyrolles. problem, Proc. R. Soc. London Ser. A 234, 528-541.
Boore, D. M. (1972). A note on the effect of simple topography on seismic Gazdag, J. (1981). Modeling of the acoustic wave equation with transform
SH waves, Bull. Seism. Soc. Am. 62, 275-284. methods, Geophysics 46, 854-859.
Bouchon, M. (1979). Discrete wavenumber representation of elastic wave Gilbert, F. and L. Knopoff (1960). Seismic scattering from topographic
fields in three space dimensions, J. Geophys. Res. 84, 3609-3614. irregularities, J. Geophys. Res. 65, 3437-3444.
Bouchon, M. and K. Aki (1977). Discrete wavenumber representation of Givoli, D. (1991). Non-reflecting boundary conditions: review article, J.
seismic-source wave fields, Bull. Seism. Soc. Am. 67, 259-277. Comput. Phys. 94, 1-29.
Bouchon, M. and J. S. Barker (1996). Seismic response of a hill: the ex-
Givoli, D. and J. B. Keller (1990). Non-reflecting boundary conditions for
ample of Tarzana, California, Bull. Seism. Soc. Am. 86, 66-72.
elastic waves, Wave Motion 12, 261-279.
Bouchon, M., C. A. Schultz, and M. N. Ttksoz (1996). Effect of three-
Gottlieb, D. (1981). The stability of pseudospectral Chebyshev methods,
dimensional topography on seismic motion, J. Geophys. Res. 101,
Math. Comp. 36, 107-118.
5835-5846.
Horike, M., H. Uebayashi, and T. Takeuchi (1990). Seismic response in
Campillo, M. and M. Bouchon (1985). Synthetic SH seismograms in lat-
three dimensional sedimentary basin due to plane S wave incidence,
erally varying medium by discrete wavenumber method, Geophys. J.
J. Phys. Earth 38, 261-284.
R. Astr. Soc. 83, 307-317.
Canuto, C. and D. Funaro (1988). The Schwarz algorithm for spectral meth- Hudson, J. A. and L. Knopoff (1964). Transmission and reflection of sur-
ods, SIAM J. Numer. Anal. 25, 24-40. face waves at a comer, part 2: Rayleigh waves (theoretical), J. Geo-
C arcione, J. M. (1991). Domain decomposition for wave propagation prob- phys. Res. 69, 281-289.
lems, J. Sci. Comp. 6, 453-472. Hughes, T. J. R. (1987). The Finite Element Method, Linear Static and
Carcione, J. M. and P. J. Wang (1993). A Chebyshev collocation method Dynamic Finite Element Analysis, Prentice-Hall International, Engle-
for the wave equation in generalized coordinates, Comp. Fluid Dyn. wood Cliffs, New Jersey.
J. 2, 269-290. Hughes, T. J. R. and J. E. Marsden (1978). Classical elastodynamics as a
Carver, D. L., K. W. King, E. Cranswick, D. W. Worley, P. Spudich, and linear symmetric hyperbolic system, J. Elasticity 8, 97-110.
C. Mueller (1990). Digital recordings of aftershocks of the October Hulbert, G. M. and T. J. R. Hughes (1990). Space-time finite element meth-
17, 1989, Loma Prieta, California, earthquake: Santa Cruz, Los Gatos, ods for second-order hyperbolic equations, Comp. Meth. Appl. Mech.
and surrounding areas, U.S. Geol. Surv. Open-File Rept. 90-683, Eng. 84, 327-348.
204 pp. Iwata, T., K. Hatayama, H. Kawase, and K. Irikura (1996). Site amplifi-
Clayton, R. and B. Engquist (1977). Absorbing boundary conditions for cation of ground motions during aftershocks of the 1995 Hyogo-ken
acoustic and elastic wave equations, Bull. Seism. Soc. Am. 67, 1529- Nanbu earthquake in severely damaged zone, J. Phys. Earth 44, 563-
1540. 576.
Clouser, R. H. and C. A. Langston (1995). Modeling observed P-Rg con- Jih, R. S., K. L. McLanghlin, and Z. A. Der (1988). Free-boundary con-
versions from isolated topographic features near the NORESS array, ditions of arbitrary polygonal topography in a two-dimensional ex-
Bull. Seism. Soc. Am. 85, 859-873. plicit elastic finite-difference scheme, Geophysics 53, 1045-1055.
Dablain, M. A. (1986). The application of high-order differencing to the Kawase, H. and K. Aid (1989). A study on the response of a soft sedimen-
scalar wave equation, Geophysics 51, 54-456. tary basin for incident S, P and Rayleigh waves, with special reference
de Bremaecker, J. C. (1958). Transmission and reflection of Rayleigh waves to the long duration observed in Mexico City, Bull Seism. Soc. Am.
at comers, Geophysics 23, 253-266. 79, 1361-1382.
de Hoop, A. T. (1960). A modification of Cagniard's method for solving Kelly, K. R., R. W. Ward, S. Treitel, and R. M. Alford (1976). Synthetic
seismic pulse problems, Appl. Sci. Res. B8, 349-356. seismograms: a finite difference approach, Geophysics 41, 2-27.
Dravinski, M. and T. K. Mossessian (1987). Scattering of plane harmonic Khun, M. J. (1985). A numerical study of Lamb's problem, Geophys.
P, SV and Rayleigh waves by dipping layers of arbitrary shape, Bull. Prosp. 33, 1103-1137.
Seism. Soc. Am. 77, 212-235. Kim, J. and A. Papageorgiou (1993). Discrete wavenumber boundary ele-
Engquist, B. and A. Majda (1977). Absorbing boundary conditions for the ment method for 3D scattering problems, J. Eng. Mech. ASCE 119,
numerical simulation of waves, Math. Comp. 31, 629-651. 603-624.
Faccioli, E., F. Maggio, A. Quarteroni, and A. Tagliani (1996). Spectral- Komatitsch, D. (1997). Mtthodes spectrales et 616ments spectraux pour
domain decomposition methods for the solution of acoustic and elastic l'tquation de l'61astodynamique 2D et 3D en milieu htttrog~ne,
wave equations, Geophysics 61, 1160-1174. Ph.D. Thesis, Institut de Physique du Globe de Paris, Paris.
Fisher, P. (1989). Spectral element solution of the Navier-Stokes equations Komatitsch, D., F. Coutel, and P. Mora (1996). Tensorial formulation of
on high performance distributed-memory parallel processors, Ph.D. the wave equation for modelling curved interfaces, Geophys. J. Int.
Thesis, Massachusetts Institute of Technology, Cambridge, Massa- 127, 156-168.
chusetts. Komatitsch, D., J. P. Vilotte, R. Vai, and F. J. SLnchez-Sesma (1998). The
390 D. Kornatitsch and J.-P. Vilotte
Spectral Element method for elastic wave equations: application to S~inchez-Sesma, F. J. (1997). Strong ground motion and site effects, in
2D and 3D seismic problems, Int. J. Num. Meth. Eng., submitted. Computer Analysis and Design of Earthquake Resistant Structures,
Kosloff, D. and E. Baysal (1982). Forward modeling by the Fourier method, D. Beskos and S. Angnostopoulos (Editors), Comp. Mech. Publica-
Geophysics 47, 1402-1412. tions, Southampton, 201-239.
Kosloff, D. and H. Tal-Ezer (1993). A modified Chebyshev pseudospectral Sfinchez-Sesma, F. J. and M. Campillo (1993). Topographic effects for
method with an O(N-l) time step restriction, J. Comput. Phys. 104, incident P, SV and Rayleigh waves, Tectonophysics 218, 113-125.
457-469. S~nchez-Sesma, F. J. and F. Luz6n (1996). Can horizontal P waves be
Kosloff, D., D. Kessler, A. Q. Fitho, E. Tessmer, A. Behle, and R. Strahil- trapped and resonate in a shallow sedimentary basin? Geophys. J. Int.
evitz (1990). Solution of the equations of dynamics elasticity by a 124, 209-214.
Chebyshev spectral method, Geophysics 55, 748-754. S~nchez-Sesma, F. J. and E. Rosenblueth (1979). Ground motion at can-
Lamb, H. (1904). On the propagation of tremors over the surface of an yons of arbitrary shape under incident SH waves, Int. J. Earthquake
elastic solid, Phil Trans. R. Soc. London. Ser. A 203, 1-42. Eng. Struct. Dyn. 7, 441-450.
Lapwood, E. R. (1961). The transmission of a Rayleigh pulse round a cor- Simo, J. C., N. Tarnow, and K. K. Wong (1992). Exact energy-momentum
ner, Geophys. J. 4, 174-196. conserving algorithms and symplectic schemes for nonlinear dynam-
Levander, A. R. (1988). Fourth-order finite-difference P-SV seismograms, ics, Comp. Meth. Appl. Mech. Eng. 100, 63-116.
Geophysics 53, 1425-1436. Spudich, P., M. Hellweg, and W. H. K. Lee (1996). Directional topographic
Lysmer, J. and L. A. Drake (1972). A finite element method for seismology, site response at Tarzana observed in aftershocks of the 1994 North-
in Methods in Computational Physics, Vol. 11, Academic, New York. ridge, California, earthquake: implications for mainshock motions,
Lysmer, J. and R. L. Kuhlemeyer (1969). Finite dynamic model for infinite Bull Seism. Soc. Am. 86, $193-$208.
media, J. Eng. Mech. Div. ASCE 95 EM 4, 859-877. Stacey, R. (1988). Improved transparent boundary formulations for the elas-
Madariaga, R. (1976). Dynamics of an expanding circular fault, Bull Seism. tic wave equation, Bull Seism. Soe. Am. 78, 2089-2097.
Soc. Am. 65, 163-182. Szab6, B. and I. Babugka (1991). Finite Element Analysis, Wiley, New
Maday, Y. and A. T. Patera (1989). Spectral element methods for the in- York.
compressible Navier-Stokes equations, in State of the Art Survey in Tessmer, E., D. Kosloff, and A. Behle (1992). Elastic wave propagation
Computational Mechanics, A. K. Noor and J. T. Oden (Editors), simulation in the presence of surface topography, Geophys. J. Int.
ASME, New York, 71-143. 108, 621-632.
Maday, Y. and A. Quarteroni (1982). Approximation of Burgers equation Toki, K., K. Irikura, and T. Kagawa (1995). Strong motion records in the
by pseudospectral methods, RAIRO Anal Numer. 16, 375-404. source area of the Hyogo-ken Nanbu earthquake, January 17, 1995,
Maday, Y. and E. M. Rcnquist (1990). Optimal error analysis of spectral J. Nat. Disast. Sci. 16, 23-30.
methods with emphasis on non-constant coefficients and deformed Tordjman, N. (1995). t~lrments finis d'ordre 61ev6 avec condensation de
geometries, Comp. Meth. AppL Mech. Eng. 80, 91-115. masse pour l'rquation des ondes, Ph.D. Thesis, Universit6 Paris IX
Manolis, G. and D. E. Beskos (1988). Boundary Element Methods in Elas- Dauphine, Paris.
todynamics, Umwin Hyman, London. Toshinawa, T. and T. Ohmachi (1992). Love wave propagation in three-
Marfurt, K. J. (1984). Accuracy of finite-difference and finite-element mod- dimensional sedimentary basin, Bull. Seism. Soc. Am. 82, 1661-1667.
eling of the scalar wave equations, Geophysics 49, 533-549. Umeda, Y., A. Kuroiso, K. Ito, Y. Ito, and T. Saeki (1986). High acceler-
Ohminato, T. and B. A. Chouet (1997). A free-surface boundary condition ations in the epicentral area of the western Nagana prefecture, Japan,
for including 3D topography in the finite difference method, Bull. earthquake of 1984, J. Seism. Soc. Jpn. 39, 217-228.
Seism. Soc. Am. 87, 494-515. Virieux, J. (1986). P-SV wave propagation in heterogeneous media: veloc-
Ohori, M., K. Koketsu, and T. Minomi (1992). Seismic response of three ity-stress finite-difference method, Geophysics 51, 889-901.
dimensional sediment filled valleys due to incident plane waves, J. Zhang, L. and A. K. Chopra (1991). Three-dimensional analysis of spatially
Phys. Earthquake 40, 209-222. varying ground motions around a uniform canyon in a homogeneous
Olsen, K. B. and R. J. Archuleta (1996). 3-D simulation of earthquakes on half-space, Int. J. Earthquake Eng. Struct. Dyn. 20, 911-926.
the Los Angeles fault system, Bull Seism. Soc. Am. 86, 575-596.
Orszag, S. A. (1980). Spectral methods for problems in complex geome-
tries, J. Comput. Phys. 37, 70-92.
Patera, A. T. (1984). A spectral element method for fluid dynamics: laminar Appendix A
flow in a channel expansion, J. Comput. Phys. 54, 468--488.
Pedersen, H. A., B. LeBrun, D. Hatzfeld, M. Campillo, and P. Bard (1994). The semi-discrete variational form of the momentum
Ground motion amplitude across ridges, Bull. Seism. Soc. Am. 84, equation, corresponding to equation (13), can be given by
1786-1800.
Pilant, W. L. (1979). Elastic Waves in the Earth, Elsevier, Amsterdam.
Pitarka, A. and K. Irikura (1996). Basin structure effects on long period
A p;ch. W h d V + Vwh:e:Vu h d
strong motions in the San Femando valley and the Los Angeles basin e=l e
from the 1994 Northridge earthquake and aftershocks, Bull. Seism.
Soc. Am. 86, S126-S137.
= A f. w h dV
Priolo, E., J. M. Carcione, and G. Seriani (1994). Numerical simulation of
interface waves by high-order spectral modeling techniques, £
e=l e (A1)
Acoust. Soc. Am. 95, 681-693.
Richter, G. R. (1994). An explicit finite element method for the wave equa- "4- hat[ T . dr
•~IT e
tion, AppL Num. Math. 16, 65-80.
Robertsson, J. O. A. (1996). A numerical flee-surface condition for elas-
tic/viscoelastic finite-difference modeling in the presence of topog-
raphy, Geophysics 61, 1921-1934.
Rulf, B. (1969). Rayleigh waves on curved surfaces, J. Acoust. Soc. Am.
45, 493-499.
Sfinchez-Sesma, F. J. (1983). Diffraction of elastic waves by three-dimen- w h e r e A is t h e s t a n d a r d a s s e m b l y o p e r a t o r o v e r t h e t o t a l
sional surface irregularities, Bull Seism. Soc. Am. 73, 1621-1636. n u m b e r o f s p e c t r a l e l e m e n t s , a n d uhu E ShN a n d w h E q/hu.
All the mathematical integrations involved at the ele- belong to the definition of the spatial support of the moment
ment level are approximated by numerical ones. With the density distribution.
help of the local geometrical mapping c/re, each integral is In contrast with classical finite-element methods, the
pulled back from the global coordinate system x to the local mass matrix is diagonal by construction, leading to fully
coordinate system {. On each parent domain, the integral is explicit time schemes. Heterogeneities can be handled nat-
evaluated using the Gauss-Lobatto Legendre quadrature: urally in two different ways: either by prescribing different
material properties for each spectral element or by assigning
different properties at each collocation point, allowing,
V :c:Vu dV =
therefore, for the description of highly variable velocity
e i,j,k
structures.
~7wh ]~e (~i' qj' (k):C(~i' /~j' (k): (A2)
~Tuhln~ (~i, t/j, (k)lJ~(~i, rlj, (k)lC%k, Appendix B
The semi-discrete momentum equation is enforced in

f n f . w ~ dV = ~ (A3) conservative form at t, + a" This energy-momentum conserv-
e
ing scheme can be generalized to a predictor-multi-corrector
Wh~ (~i, r/j, ~k) " fln~ (~i, rlj, ~k)lJe(~i, rlj, ~k)logijk, format that improves its characteristics under some circum-
stances, and it allows for an efficient parallel implementa-
where --N'~e= {(~i, /~j, (k); 0 =< i =
< N; 0 --<j _<--N; 0 =< tion. The implementation sequence is as follows:
k -< N} denotes the integration grid of(N + 1)3 points within
the reference domain; og/jk = co,~o/okwith (Dl > 0 denotes
the 1D Gauss-Legendre Lobatto quadrature weights, and i = 0 (i is the iteration number)
denotes the pull back of the gradient operator in the
coordinate system: Predictor phase:
VxW~/I~'~e ----~ ~W~/[~'~e (~) ---- V~WhN]~[~e( ~ ) F e 1 (~), (A4) d(ni) 1 -- d n + l "v(i)
n+l = 0 a(ni)+1 = fin+ 1"
where Fe(~) = 0gc-/Ce(~) is the gradient of the local geomet- Solution phase:
rical transformation. With the help of the Lagrangian inter-
polation, this leads to the system of ordinary differential
1 MAv<° = Fn+a
~-~ ext _ Fint (d(i)
,un+a, --n+a,
1 ~,~r,,,(i)
..,r(i) a -- ~-~ l,*~--n+ 1
- vn).
equations (23), where now
M = A
ne, . ( I N ( ~ I N ( ~ I N ) f ~ ( 1 N ( ~ I N @ I N ) p dV, Collector phase:
e=l e
(A5) v(i+
n + l 1) ~--- Vn+l
,tr(i) -[- AV,
F int = A n'Fo
=1 ~
V(/N (~ Y (~) / N ) : c : V u ~ dV
(A6)
d(i+l)
n+l :
a.+ 1 + flAt
?
v(ni++~
),
1
- .-Iv°x'l~ (ff @ IN @ IN)t(v~) dF/, a(i+ 1) ~ 1 ..,(i+1)

n+l : an+l -~t--n+l ,
F ext = A ( f f @ l N @ I N ) f dV
e= 1 e (A7) where the predictors are defined as
+
Irintl (IN @ IN @ IN)T d r ]
•d ~ T e
an+, = d n + At 1 -7 vn + At2 -7 an,

For a seismic source described in terms of a distribution of
moment-tensor density fix, 0 = -div[m(x, t)], with a lo- ~ln+l = (1 -- ~ ) an -- ~-~
1 Vn,
calized spatial support, the generalized body force can be
written after integration by part as
with d, +a = ~dn+ 1 + (1 - a)dn and F ~ a = ¢~F~nX~m --[-
(1 -- ot)Fext
n .
eA=l e v(IN @ IN @ lN):m(x' t) dV, (A7)
Exact conservation of the total angular momentum is

where the assembly operation involves all the elements that achieved for the values a = fl/y = 1/2 corresponding to the
conservation form o f the mid-point rule. These values define handle a more complex constitutive behavior, in particular
an acceleration-independent time-marching algorithm. Sec- to incorporate attenuation.
ond-order accuracy is achieved if and only if a = 1/2. For
the parameter values o~ = fl/7 = 1/2, a linear analysis shows D6partement de Sismologie (URA 195)
Institut de Physique du Globe de Paris
that the spurious root at zero sampling frequency vanishes
4, Place Jussieu
if and only if y = 1. This value yields the post-processing 75252--Paris Cedex 05, France
formula an+ 1 = (vn+ 1 - vn)/At. Within this energy-mo-
mentum conserving framework, there is no difficulty to Manuscript received 17 July 1997.
View publication stats

The Spectral Element Method: An Efficient Tool To Simulate The Seismic Response of 2D and 3D Geological Structures

Uploaded by

Copyright:

Available Formats

The Spectral Element Method: An Efficient Tool To Simulate The Seismic Response of 2D and 3D Geological Structures

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

The Spectral Element Method: An Efficient Tool To Simulate The Seismic Response of 2D and 3D Geological Structures

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

Article in Bulletin of the Seismological Society of America · April 1998

Belmont Forum e-Infrastructures and Data Management Project View project

Dynamic rupture modelling View project

The user has requested enhancement of the downloaded file.

The Spectral Element Method: An Efficient Tool to Simulate

Abstract We present the spectral element method to simulate elastic-wave prop-

and a body force derived from a seismic moment tensor den-

Initial Boundary-Value Problem f(x, t) = -div [mo(x)] 9(t - to), (8)

where t is the traction vector on the boundary, n is the unit

,St = {U(X, t) ~ H 1 (~'~)na : ~'~ Numerical Discretization

whla, (X, y, z) = xN (win) (22) a~+l - ?At - an, (26)

--4 ---+ ,.~ ~

To assess some of the properties of the time-marching

Ux residuals trace 75 Uz residuals trace 75

Ux residuals trace 100 Uz residuals trace 100

Ux residuals trace 55 Uz residuals trace 55

U x residuals trace 80 Uz residuals trace 80

5 10 15 20 25 30 35 40 Three-Dimensional Numerical Simulations

Figure 6. Snapshot showing wave propa-

Amplification by a Three-Dimensional Hill

Figure 8. Snapshots showing wave propa-

Uy component along Y Uz component along Y

~" 1000- ~ 1000-

d 012 014 016

Uy component along Y Uy component along X

6 012 0'.4 016 0 012 014 016

Figure 13. Time responses of the y (horizontal) component of the displacement

Figure 15. Polarization diagrams at differ-

Figure 16. Polarization diagrams at differ-

40O 800 12OO 16OO 400 800 1200 1600

Receiver # 105 Receiver # 65

0.15 0.4~ 0.65 0.90 0.15 0.40 0.65 0.90

Uy component along X Uy component along Y

~Tuhln~ (~i, t/j, (k)lJ~(~i, rlj, (k)lC%k, Appendix B

The semi-discrete momentum equation is enforced in

- .-Iv°x'l~ (ff @ IN @ IN)t(v~) dF/, a(i+ 1) ~ 1 ..,(i+1)

an+, = d n + At 1 -7 vn + At2 -7 an,

Exact conservation of the total angular momentum is

View publication stats

You might also like