Spectral Elements Introduction recalling the elastic wave equation The spectral-element method: General concept domain mapping from space-continuous to space-discrete time extrapolation Gauss-Lobatto-Legendre interpolation and integration A special flavour of the spectral-element method: SES3D programme code description computation of synthetic seismograms long-wavelength equivalent models

Scope: Understand the principles of the spectral element method and why it is currently maybe the most important method for wave propagation. This lecture based on notes by Andreas Fichtner. Spectral element method 1 THE ELASTIC WAVE EQUATION Elastic wave equation: L(u, , C) f 2 t

( x, t ) : u( x, t) d L(u, ,C) (x) u( x, t) C Subsidiary conditions: u(x, t) | t t 0 0 t t u(x, t) |t t 0 0 n C ( x, t ) : u( x, ) d |x 0 In this formulation visco-elastic dissipation is included as well as a general anisotropic description of elasticity. Spectral element method 2

SPECTRAL-ELEMENT METHOD: General Concept Subdivision of the computational domain into hexahedral elements: (a) 2D subdivision that honours layer boundaries (b) Subdivision of the globe (cubed sphere) Spectral element method (c) Subdivision with topography 3 SPECTRAL-ELEMENT METHOD: General Concept Mapping to the unit cube: Spectral element method

4 SPECTRAL-ELEMENT METHOD: General Concept Choice of the collocation points: Interpolation of Runges function R(x) using 6th-order polynomials and equidistant collocation points 1 R ( x) 1 ax 2 interpolant Runges phenomenon Spectral element method 5

SPECTRAL-ELEMENT METHOD: General Concept Choice of the collocation points: Interpolation of Runges function R(x) using 6th-order polynomials and Gauss-Lobatto-Legendre collocation points [ roots of (1-x2)LoN-1= completed Lobatto polynomial ] 1 R ( x) 1 ax 2 interpolant We should use the GLL points as collocation points for the Lagrange polynomials. Spectral element method

6 SPECTRAL-ELEMENT METHOD: General Concept Example: GLL Lagrange polynomials of degree 6 collocation points = GLL points global maxima at the collocation points Spectral element method 7 The SE system Diagonal mass matrix M Spectral element method

8 SPECTRAL-ELEMENT METHOD: General Concept Numerical quadrature to determine mass and stiffness matrices: Quadrature node points = GLL points The mass matrix is diagonal, i.e., trivial to invert. This is THE advantage of the spectral-element method. Time extrapolation: (t ) u u(t t ) 2u(t ) u(t t ) t 2 u(t t ) 2u(t ) u(t t ) t 2 M 1 f (t ) Ku (t )

Spectral element method 9 SPECTRAL-ELEMENT METHOD: General Concept Representation in terms of polynomials: N u(x, t ) u i (t ) (iN ) (x) (within the unit interval [-1 1]) i 0 (iN ) (x) : Nth-degree Lagrange polynomials We can transform the partial differential equation into an ordinary differential equation

where we solve for the polynomial coefficients: i K kiu i f k M kiu M ki : mass matrix K ki : stiffness matrix Spectral element method 10 SES3D: General Concept Simulation of elastic wave propagation in a spherical section. Spectral-element discretisation. Computation of Frchet kernels using the adjoint method. Operates in natural spherical coordinates!

3D heterogeneous, radially anisotropic, visco-elastic. PML as absorbing boundaries. Programme philosophy: Puritanism [easy to modify and adapt to different problems, easy

implementation of 3D models, simple code] Spectral element method 11 SES3D: Example Southern Greece 1. Input files [geometric setup, source, receivers, Earth model] 8 June, 2008

2. Forward simulation [wavefield snapshots and seismograms] Mw=6.3 3. Adjoint simulation [adjoint source, Frchet kernels] Spectral element method 12 SES3D: Input files Par: - Numerical simulation parameters - Geometrical setup - Seismic Source - Parallelisation stf: - Source time function

recfile: - Receiver positions Spectral element method 13 SES3D: Parallelisation Spherical section subdivided into equal-sized subsections Each subsection is assigned to one processor. Communication: MPI Spectral element method 14 SES3D: Source time function Source time function

- time step and length agree with the simulation parameters - PMLs work best with bandpass filtered source time functions - Example: bandpass [50 s to 200 s] Spectral element method 15 Simulating delta functions? Spectral element method 16 LONG WAVELENGTH EQUIVALENT MODELS Single-layered crust that coincides with

the upper layer of elements and PREM below boundary between the upper 2 layers of elements lon=142.74 lat=-5.99 d=80 km vertical displacement SA08 lon=150.89 lat=-25.89 Dalkolmo & Friederich, 1995. Complete synthetic seismograms for a spherically symmetric Earth , GJI, 122, 537-550 Spectral element method

17 LONG WAVELENGTH EQUIVALENT MODELS verification 2-layered crust that does not coincide with a layer of elements and PREM below boundary between the upper 2 layers of elements lon=142.74 lat=-5.99

d=80 km vertical displacement SA08 lon=150.89 lat=-25.89 Dalkolmo & Friederich, 1995. Complete synthetic seismograms for a spherically symmetric Earth , GJI, 122, 537-550 Spectral element method 18 LONG WAVELENGTH EQUIVALENT MODELS long wavelength equivalent models Replace original crustral model by a longwavelength equivalent

model which is transversely isotropic [Backus, 1962]. The optimal smooth model is found by dispersion curve matching. Fichtner & Igel, 2008. Efficient numerical surface wave propagation through the optimisation of discrete crustal models, GJI. Spectral element method 19 LONG WAVELENGTH EQUIVALENT MODELS

long wavelength equivalent models Minimisation of the phase velocity differences for the fundamental and higher modes in the frequency range of interest through simulated annealing. Spectral element method 20 LONG WAVELENGTH EQUIVALENT MODELS long wavelength equivalent models crustal thickness map (crust2.0) 3D solution: interpolation of long wavelength equivalent profiles to obtain 3D crustal model. Problem 1: crustal structure not well constrained (receiver function non-uniqueness) Problem 2: abrupt changes in crustal structure (not captured by pointwise RF studies)

Spectral element method 21 SES3D: Calls to caution! 1. Long-term instability of PMLs - All PML variants are long-term unstable! - SES3D monitors the total kinetic energy Etotal. - When Etotal increases quickly, the PMLs are switched off and - absorbing boundaries are replaced by less efficient multiplication by small numbers. 2. The poles and the core - Elements become infinitesimally small at the poles and the core. - SES3D is efficient only when the computational domain is sufficiently far from the poles and the core. 3. Seismic discontinuities and the crust - SEM is very accurate only when discontinuities coincide with element boundaries. - SES3Ds static grid may not always achieve this. - It is up to the user to assess the numerical accuracy in cases where discontinuities run through elements. [Implement long-wavelength equivalent models.]

- Generally no problem for the 410 km and 660 km discontinuities. Spectral element method 22 Spectral elements: summary Spectral elements (SE) are a special form of the finite element method. The key difference is the choice of the basis (form) functions inside the elements, with which the fields are described.

It is the Lagrange polynomials with Gauss-Lobato-Legendre (GLL) collocation points that make the mass matrix diagonal This leads to a fully explicit scheme without the need to perform a (sparse) matrix inverse inversion Material parameters can vary at each point inside the elements SE works primarily on hexahedral grids The hexahedra can be curvilinear and adapt to complex geometries (cubed sphere, reservoir models) Two open-source codes are available here: www.geodynamics.org (specfem3d) regional and global scale www.geophysik.uni-muenchen.de/Members/fichtner (ses3d) - regional scale Spectral element method 23