The direct spectral element method for the calculation of synthetic seismograms in self-gravitating, spherically symmetric planets

This paper introduces $\texttt{DSpecM1D}$, a new implementation of the direct spectral element method using a displacement formulation to calculate synthetic seismograms in self-gravitating, spherically symmetric, anelastic, and transversely isotropic Earth models, which is validated against existing codes with excellent agreement.

The Gist

Imagine the Earth not as a static rock, but as a giant, vibrating bell. When a massive earthquake strikes, it doesn't just shake the ground; it sets the entire planet ringing like a bell struck by a hammer. These "rings" are called free oscillations, and they travel through the Earth's core, mantle, and crust. By listening to these rings, scientists can figure out what the Earth is made of deep down, much like a doctor uses an ultrasound to see inside a body.

However, listening to the Earth is incredibly difficult. The Earth is heavy, it has its own gravity pulling on itself, it has liquid oceans and a molten core, and it isn't perfectly round or uniform. Calculating how these vibrations move through such a complex, self-gravitating object is a mathematical nightmare.

This paper introduces a new, super-smart calculator called DSpecM1D that solves this problem. Here is how it works, explained simply:

1. The Problem: The "Gravity" and "Liquid" Trap

Previous methods for calculating these Earth vibrations had two big flaws:

• The Gravity Issue: Most methods treated the Earth's gravity as a separate, messy side-effect. But in reality, the Earth's own weight (self-gravity) is a huge part of the physics. Ignoring it is like trying to predict how a trampoline bounces without accounting for the weight of the person jumping on it.
• The Liquid Issue: The Earth's outer core is liquid. In liquids, things behave differently than in solids (think of how water sloshes vs. how a rubber ball bounces). Old methods often had to pretend the liquid core was "neutral" or perfectly balanced to make the math work. If the real Earth wasn't perfectly balanced, those old methods failed.

2. The Solution: The "Lego Tower" Approach

The authors built a new method using something called the Direct Spectral Element Method.

Imagine you want to measure the sound of a giant drum. Instead of trying to calculate the sound for the whole drum at once (which is impossible), you break the drum down into tiny, manageable slices, like a stack of Lego bricks.

• The Bricks (Spectral Elements): They divide the Earth into many thin, horizontal layers (elements).
• The Magic (Spectral): Inside each "brick," they don't just use simple straight lines; they use complex, high-order curves (polynomials) that can wiggle and bend to perfectly match the shape of the wave. This is like using a flexible, high-tech ruler instead of a stiff, cheap one.

3. The Big Innovation: One Language for Everyone

The most clever part of this new code is that it speaks the same language for both the solid rock and the liquid core.

• Old Way: It was like having a translator for the solid parts and a completely different translator for the liquid parts. When they met at the boundary, the translation often got messy, especially if the liquid was sloshing around (stratified).
• New Way (DSpecM1D): They use a displacement formulation for the entire planet. Whether it's solid rock or liquid metal, the math treats the movement of particles in the exact same way. This allows the code to handle the liquid core even if it's not perfectly balanced, capturing the "sloshing" effects that other codes miss.

4. Why Does This Matter?

Think of the Earth as a giant, complex musical instrument.

• The Goal: Scientists want to know exactly how this instrument sounds so they can figure out what materials are inside it (like finding a hidden treasure map inside the Earth).
• The Result: The authors tested their new calculator (DSpecM1D) against two other famous calculators (MINEOS and YSpec). The results were nearly identical, proving their new method is accurate.
• The Future: This code is actually a "training wheels" or a "pre-conditioner" for a much bigger project. The authors are building a system to simulate the Earth with all its wobbles, rotations, and uneven bumps (3D models). To do that efficiently, they need a perfect calculator for a simple, round Earth to help speed up the complex calculations. DSpecM1D is that perfect, round-Earth calculator.

In a Nutshell

The authors have built a new, highly accurate digital tool that simulates how the Earth vibrates after an earthquake. By treating the solid and liquid parts of the Earth with the same mathematical rules and fully accounting for the planet's own gravity, they have created a more reliable way to "listen" to the Earth's deep interior. This tool is a crucial stepping stone toward creating the most detailed 3D movies of seismic waves ever made.

Technical Summary

Here is a detailed technical summary of the paper "The direct spectral element method for the calculation of synthetic seismograms in self-gravitating, spherically symmetric planets" by Myhill and Al-Attar.

1. Problem Statement

The paper addresses the challenge of accurately and efficiently calculating synthetic seismograms for self-gravitating, spherically symmetric, non-rotating, anelastic, and transversely isotropic (SNRATI) Earth models.

• Context: Long-period free oscillation studies are critical for understanding Earth's density structure and lateral variations (e.g., Large Low Velocity Provinces). However, modeling these spectra in laterally heterogeneous models is computationally difficult.
• Limitations of Existing Methods:
  • Normal Mode Summation: Requires solving a massive eigenvalue problem first, which is computationally more expensive than solving the forced equations of motion directly. Furthermore, it suffers from a theoretical limitation where eigenfunction bases may lack the necessary degrees of freedom to satisfy boundary conditions in heterogeneous models, potentially preventing convergence to the true solution.
  • Direct Radial Integration: While efficient, previous implementations often struggled with arbitrary fluid stratification or required potential formulations in fluid regions that assumed neutral stratification.
  • 3D Spectral Element Methods (Time Domain): Prohibitively expensive for the hundreds of hours of time-series data required in free oscillation studies.
• Specific Challenge: The presence of the Earth's fluid outer core introduces "undertones" (essential spectrum) near zero frequency. Numerical methods often struggle to resolve these without assuming neutral stratification or using complex potential formulations.

2. Methodology

The authors developed a new code, DSpecM1D (Direct Spectral Element Method 1D), which implements a Direct Solution Method (DSM) using radial spectral elements.

• Formulation:
  • The method uses a displacement formulation throughout the entire model (both solid and fluid regions), unlike previous DSM implementations that switched to a potential formulation in fluids.
  • It solves the frequency-domain equations of motion directly, avoiding the need to solve the associated eigenvalue problem first.
  • Self-gravitation is fully incorporated.
• Numerical Discretization:
  • Galerkin Approach: The wavefield is discretized using a 1D spectral element basis. The domain is divided into elements, and within each, the solution is approximated using high-order Gauss-Lobatto-Legendre (GLL) polynomials.
  • Spherical Symmetry: The problem is reduced to a system of Ordinary Differential Equations (ODEs) for each spherical harmonic degree ($l$) and order ($m$) via generalised spherical harmonic expansions.
  • Fluid Handling: By using a displacement formulation, the code can handle arbitrary fluid stratification (non-neutral) without artificial constraints. It naturally resolves the "undertone" issue by ensuring the mesh is fine enough to resolve the oscillatory nature of the solution in the fluid core.
  • Boundary Conditions: Continuity of displacement is enforced at solid-solid boundaries. At solid-fluid boundaries, tangential slip is allowed, consistent with inviscid fluid dynamics.
• Implementation Details:
  • Anelasticity: Attenuation and dispersion are handled by allowing the Love parameters (elastic moduli) to be frequency-dependent.
  • Truncation: To improve efficiency, the matrix equations are truncated at a "turning depth" where the solution becomes evanescent, effectively applying a free boundary condition.
  • Parallelization: The code is highly parallelized over spherical harmonic degrees ($l$) and frequencies, as these are independent.
  • Solver: Sparse LU decomposition is used to solve the linear system at each frequency.

3. Key Contributions

1. Unified Displacement Formulation: The primary innovation is the use of a displacement formulation across both solid and fluid regions. This allows for the full inclusion of self-gravitation and arbitrary fluid stratification, overcoming the limitations of potential-based methods in the outer core.
2. DSpecM1D Code: The release of a new, open-source code that significantly extends the Direct Solution Method (DSM) to SNRATI models.
3. Preconditioner for 3D Problems: The code is designed to serve as a "spherical Earth preconditioner" for future iterative solutions of 3D, laterally heterogeneous, and rotating Earth models (building on work by Maitra & Al-Attar, 2019).
4. Handling of Undertones: The study demonstrates that with sufficient mesh resolution, the spectral element method can accurately resolve the complex oscillatory behavior of fluid undertones without requiring special potential representations or assuming neutral stratification.

4. Results and Benchmarks

The DSpecM1D code was benchmarked against two established codes:

• MINEOS: A normal mode summation code.
• YSpec: A direct radial integration code.

Key Findings:

• Accuracy: Agreement between DSpecM1D, YSpec, and MINEOS is excellent for both elastic and anelastic models.
  • Spectra: For the 1994 Bolivia earthquake, the average relative misfit between DSpecM1D and YSpec was 0.007%.
  • Waveforms: For displacement and acceleration traces, the mean misfit between DSpecM1D and YSpec was <1% across all components. Misfits against MINEOS were slightly higher (around 1-7% in specific cases) but still within acceptable limits for seismological modeling.
• Convergence: The method exhibits high-order convergence. A 5-point Spectral Element Method (SEM) was identified as a cost-effective balance, requiring roughly one element per wavelength to achieve 0.1% relative error.
• Fluid Core Behavior: In tidal forcing simulations (low frequencies), the code successfully resolved the excitation of undertones in the outer core. The results showed that while coarse meshes produce unphysical responses, sufficiently fine meshes resolve the essential spectrum correctly, confirming that the displacement formulation is robust for non-neutral stratification.
• High-Frequency Application: The code successfully generated record sections up to 2-second periods, demonstrating its versatility for higher-frequency seismological applications beyond just free oscillations.

5. Significance

• Theoretical Advancement: The paper resolves a subtle theoretical issue regarding the convergence of normal mode coupling in heterogeneous models by providing a direct solution method that does not rely on truncated eigenfunction bases.
• Computational Efficiency: By solving the forced equations directly in the frequency domain using high-order spectral elements, the method offers a computationally efficient alternative to time-domain 3D spectral element methods for long-duration free oscillation studies.
• Foundation for Future Work: DSpecM1D is a critical step toward solving the full 3D seismic wave propagation problem in self-gravitating, rotating, and laterally heterogeneous Earth models. It provides the necessary preconditioner to accelerate iterative solvers in these complex scenarios.
• Open Science: The availability of the code (DSpecM1D) on GitHub promotes reproducibility and further development in global seismology.

In summary, this paper presents a robust, high-accuracy numerical method for calculating synthetic seismograms in complex spherical Earth models, bridging the gap between theoretical normal mode theory and practical spectral element implementation, with a specific focus on overcoming challenges related to self-gravitation and fluid cores.