Tarje Nissen-Meyer

Abstract. We develop a spectral-element method for computing the full 3-D moment-tensor and point-force response of a spherically symmetric earth model in a 2-D semicircular computational domain. The full elastodynamic response to a six-component moment tensor at an earthquake hypocenter and a three-component point force at a seismic station can be determined by solving six independent 2-D problems, three for a monopole source, two for a dipole source, and one for a quadrupole source. This divide-and-conquer, 3-D to 2-D reduction strategy provides a basis for the efficient computation of exact Fréchet sensitivity kernels in a spherically symmetric earth, with all wavefield features accounted for. To focus on the novel inclusion of the full source in a cylindrical coordinate system, we describe the 2-D weak formulation of the set of elastodynamic equations, its discretization using spectral elements, and the associated axial boundary conditions and source representations for each of the excitation types in the case of a homogeneous, solid elastic sphere. The method is numerically validated against both analytical solutions and normal-mode summation.

Moment-tensor source gallery. Due to the polynomial representation of functions in the spectral-element method, a point-like moment-tensor source, containing w, spreads out over the bearing element in its discretized form. The figure shows all non-zero source vector components f_β as defined in eqs (24)–(30) for the element containing the source at polynomial order N = 8. Arrow lengths represent the relative values of the magnitude of f_β. The physical source location on the axis s = 0 is denoted by a large black dot. Neighbouring elements (not shown) share the same values as the edge of the shown element, thus a maximum of four elements may have non-zero source terms. Note that the location and relative spacing of grid points within the element is different for directions parallel and perpendicular to the axis, reflecting the axial discretization in terms of GLJ points for the ξ (and s) directions, and GLL points for the η (and z) directions.

Global wave propagation upon an (M_xx + M_yy) monopole source located at 344 km depth. The top (bottom) left panel depicts the r (θ) component of both SEM-based (solid lines) and normal-mode displacements (dashed lines) along the surface in the phi = 0 plane for a dominant period of T₀ = 100 s. The amplitudes have been scaled with a factor of θ_rec^1/3 to identify most phases at all distances. Some prominent direct, multiple and surface wave phases are labelled in the top left figure. The panels on the right depict two sample traces from the sections on the left for both SEM (solid lines) and normal-mode simulations (thicker, dashed lines) in absolute displacement values. Below each seismogram is the residual waveform between SEM and normal-mode trace, magnified by a factor of 10. Note the different absolute displacement scales for each of the single trace panels.

Work in collaboration with Tony Dahlen and Alexandre Fournier, accepted in Geophys. J. Int. under
"A 2-D spectral-element method for computing spherical-earth seismograms--II. Background models", 2008.

Abstract. We portray a dedicated spectral-element method to solve the elastodynamic wave equation upon spherically symmetric earth models at the expense of a 2-D domain, computing full 3-D wavefields of arbitrary resolution to obtain Fréchet sensitivity kernels, especially for diffracted arrivals. The meshing process is presented for varying frequencies in terms of its efficiency as measured by the total number of elements, their spacing variations and stability criteria. We assess the mesh quantitatively by defining these numerical parameters in a general non-dimensional form such that comparisons to other grid-based methods are trivial. Efficient-mesh generation for the PREM example and a minimum-messaging domain decomposition and parallelization strategy lay foundations for waveforms up to frequencies of 1 Hz on moderate Beowulf clusters. The parallelization yields extremely beneficial scaling, which is highly desirable especially for many repeated calculations as necessary for a comprehensive kernel database. The discretization of fluid, solid and respective boundary regions is similar to previous spectral-element implementations but instead using a potential formulation in the fluid that incorporates the density, thus obtaining equal boundary terms on both fluid and solid sides. We compare the conventional second-order Newmark time extrapolation scheme with a newly implemented fourth-order symplectic scheme and argue for the latter option in cases of propagation over many wavelengths due to drastic accuracy improvements. Various validation examples such as full moment-tensor seismograms, wavefield snapshots, energy conservation, and some analytical on-the-fly accuracy tests illustrate the overall behavior of the final package in terms of computational cost, accuracy, and dispersion.

Left: The semicircular, solid-fluid domain Ω=Ω^s + Ω^f discretized for the PREM background model using quadrilateral elements Ω_e for dominant source period T₀=20 s. Note that all discontinuities are honored and several conforming coarsening levels are included to maintain a relatively constant resolution throughout the domain. Top right: Enlargement of the crust for one (right) and two (left) crustal layers, and the upper mantle, including one mesh coarsening region. Note the variable vertical spacing due to discontinuities. Bottom right: To circumvent the singularity at the center, we apply the following analytical expressions to reshape rectangular elements: |x|^p+|y|^p=|r|^p, x=s+z, y=s-z, 1≤p≤2. This guarantees an easy handle on grid spacing which varies maximally at the outermost, deformed elements of this central region and hence controls stability and resolution.

Domain decomposition for T₀=10 s and four (left), eight (middle), and sixteen (right two panels) processors. Note that each processor has the same number of elements ("efficiency"), and maximally two neighbors (minimal latency) with small shared boundaries (i.e., small bandwidth). Each processor touches the axis, albeit to varying amounts. Axial terms are however 1-D, and as such do not add any countable CPU time to the overall scheme. The central region is subdivided such that each processor has exactly the same amount of elements.

Selected SEM and normal-mode seismograms for the Vanuatu event at 15 km depth and the PREM model and T₀=30 s. The vertical axes denote the respective components of the actual ground displacement, and the titles the epicentral distances.

Wavefield snapshots for a monopole M_zz source at 650 km depth for T₀=10 s, including magnified panels of the crust-upper mantle and core-mantle boundary regions; labels refer to major seismic phases.