High-order exponential integration for seismic wave modeling

TL;DR

The paper addresses memory and efficiency challenges in seismic waveform inversion by solving the semi-discrete wave equation $dU/dt = HU + f(t)$ under absorbing boundary conditions. It benchmarks seven time-stepping schemes, including Faber-based exponential FA, Krylov-based KRY, and High-order RK (HORK), against classic low-order methods on 2D acoustic models with PML. Key findings show that Krylov-based exponential integration offers the best convergence among high-order methods, and high-order schemes enable larger time steps with meaningful memory savings, though Krylov can incur higher per-step costs; Leap-Frog remains fastest per step but limited by stability, while RK9-7 shows strong dissipation performance. The results suggest that leveraging large time steps with high-order exponentials can substantially reduce memory requirements in full waveform inversion, guiding method choice and highlighting potential improvements via adaptive stepping and advanced Krylov techniques like KIOPS.

Abstract

Seismic imaging is a major challenge in geophysics with broad applications. It involves solving wave propagation equations with absorbing boundary conditions (ABC) multiple times. This drives the need for accurate and efficient numerical methods. This study examines a collection of exponential integration methods, known for their good numerical properties on wave representation, to investigate their efficacy in solving the wave equation with ABC. The purpose of this research is to assess the performance of these methods. We compare a recently proposed Exponential Integration based on Faber polynomials with well-established Krylov exponential methods alongside a high-order Runge-Kutta scheme and low-order classical methods. Through our analysis, we found that the exponential integrator based on the Krylov subspace exhibits the best convergence results among the high-order methods. We also discovered that high-order methods can achieve computational efficiency similar to lower-order methods while allowing for considerably larger time steps. Most importantly, the possibility of undertaking large time steps could be used for important memory savings in full waveform inversion imaging problems.

Figures (14)

• Figure 1: Uniform staggered grid in 2D with the relative positions of the acoustic wave equations' variables and parameters. $u,\;v$ and $c$ are collocated. The shaded region represents the PML domain.
• Figure 2: Snapshots of the reference solution at times $t=0.6s$ and $T=1.3s$ within the homogeneous medium $\Omega=[0\text{km},6\text{km}]\times[0\text{km},5\text{km}]$. The Ricker signal source position (blue dot) and the receiver location (black square) are highlighted. During the time interval $t\in[0.6,1.3]s$, the front wave propagates through the receiver location.
• Figure 3: Dependence of $\Delta t_{\text{max}}$ on the approximation degree of the numerical scheme. A higher number of stages leads to an increase in the maximum allowable time step without significantly increasing the number of computations.
• Figure 4: Variation of $\Delta t_{\text{max}}$ (left) and $\text{N}^{\text{\tiny disp}}_{\text{op}}$ (right) concerning the numerical scheme and the number of stages utilized, according to the numerical dispersion error for a Ricker source peak frequency of $f_M=15$Hz. Generally, a higher number of stages leads to an increase in the maximum allowable time step size without significantly increasing the number of computations. * Here we neglect the computational complexity of creating the Krylov subspaces.
• Figure 5: Variation of $\Delta t_{\text{max}}$ (left) and $\text{N}^{\text{\tiny disp}}_{\text{op}}$ (right) concerning the numerical scheme and the number of stages utilized, according to the numerical dissipation error for a Ricker source peak frequency of $f_M=15$Hz. Generally, a higher number of stages leads to an increase in the maximum allowable time step size without significantly increasing the number of computations. * Here we neglect the computational complexity of creating the Krylov subspaces.