Spectral interpolation in semi-implicit semi-Lagrangian methods for shallow water equations on the sphere

Technical Summary: Spectral Interpolation in Semi-Implicit Semi-Lagrangian Methods for Shallow Water Equations on the Sphere

Problem Statement
Semi-implicit semi-Lagrangian (SISL) methods are the standard for operational global weather and climate models because they allow for time steps larger than those permitted by the Courant-Friedrichs-Lewy (CFL) condition. However, a significant numerical inconsistency exists in current operational schemes. While the spatial discretization of prognostic variable derivatives typically employs spectrally accurate methods (such as Spherical Harmonic Transforms [SHT] or the Double Fourier Sphere [DFS] method), the interpolation step required to evaluate variables at off-grid departure points relies on low-order techniques, such as cubic or quintic tensor-product Lagrange polynomials.

This mismatch creates a discrepancy in numerical accuracy. While low-order interpolation is often justified by computational cost concerns (with a complexity of $O(n^2)$ versus $O(n^4)$ for spectral interpolation) and theoretical error bounds for simple advection, there is limited numerical evidence in the literature justifying this trade-off for complex systems like the Shallow Water Equations (SWE). Furthermore, low-order interpolation is known to introduce excessive numerical diffusion and phase errors, potentially degrading accuracy and conservation properties over long integration times.

Methodology
This work presents the first numerical investigation of spectrally accurate interpolation within SISL schemes for the SWE on a sphere. The authors build upon the DFS-based SISL model developed by Yoshimura [56], which utilizes a stable Fourier basis and a two-time-level SETTLS scheme.

The core methodological innovation is the replacement of the standard Lagrange interpolation with a spectral interpolation scheme based on the same DFS Fourier basis used for spatial derivatives. To address the computational cost barrier, the authors employ the Nonuniform Fast Fourier Transform (NUFFT).

• Interpolation Mechanism: The prognostic variables (specifically the terms $R_v$ and $R_h$ in the semi-implicit system) are expanded using the partial regularity DFS basis. These expansions are then evaluated at the non-grid departure points.
• Complexity Reduction: Direct evaluation of bivariate Fourier series at $O(n^2)$ departure points would require $O(n^4)$ operations. By recasting the series into complex exponential form and utilizing the NUFFT (specifically the FINUFFT library), the computational complexity is reduced to $O(n^2 \log^2(1/\epsilon) + n^2 \log n)$, where $\epsilon$ is the accuracy tolerance. This renders the interpolation step computationally consistent with the rest of the spectral SISL scheme.
• Implementation: The authors implemented this spectral interpolation in two contexts:
  1. A DFS-based SISL model (using Yoshimura's stable Fourier basis).
  2. An SHT-based SISL model (using Gauss-Legendre grids), where the DFS method is adapted to generate bivariate Fourier interpolants from the non-equispaced latitude points.

Key Contributions and Results
The authors evaluated the new spectral interpolation schemes against standard low-order Lagrange interpolation (LAG) using a suite of standard SWE test cases (Williamson TC1, TC2, TC5, TC6, and the Galewsky test case). The results demonstrate several distinct advantages:

1. Accuracy and Convergence:

   • TC1 (Advection of a cosine bell): The spectral DFS interpolation maintained the same spatial convergence rate as LAG but achieved errors approximately two orders of magnitude lower. Over long integration times (10 revolutions), the DFS method reduced errors by three orders of magnitude compared to LAG.
   • TC2 (Geostrophic Flow): While LAG interpolation failed to show expected second-order temporal convergence (likely due to dominant spatial errors), the DFS interpolation schemes achieved the theoretical second-order accuracy.
   • TC5 & TC6 (Mountain Flow and Rossby-Haurwitz Waves): The DFS method consistently produced lower errors in height fields compared to LAG, particularly at finer resolutions and smaller time steps. In TC6, the DFS method better preserved the shape of the wave troughs, matching high-order Eulerian reference solutions more closely.

2. Conservation Properties:

   • The spectral interpolation significantly improved the conservation of mass and energy. In TC1, the absolute mass error after 120 days was $1.21 \times 10^{-11}$ for DFS versus $4.34 \times 10^{-4}$ for LAG.
   • In TC2 and TC5, the DFS method maintained stable mass and energy conservation errors over time, whereas LAG errors tended to increase or oscillate more significantly.

3. Numerical Diffusion and Small-Scale Dynamics:

   • Galewsky Test Case: This test involves the evolution of a nonlinear jet into turbulence. Spectral analysis of the kinetic energy spectrum revealed that LAG interpolation caused rapid dissipation of high-wavenumber energy. In contrast, the DFS interpolation schemes preserved the energy cascade, allowing the high-wavenumber spectrum to approach the theoretical $k^{-5/3}$ scaling observed in turbulent rotating flows. LAG interpolation failed to capture this behavior, indicating excessive numerical damping.

4. Computational Efficiency:

   • Wall-clock time comparisons showed that for low-accuracy requirements, LAG interpolation is faster. However, as the required accuracy increases (stricter error tolerances), the DFS interpolation method becomes more efficient, reaching the target accuracy in less time than the LAG approach. This confirms that the NUFFT-accelerated spectral interpolation is viable for high-precision applications.

Significance
The paper claims that incorporating spectrally accurate interpolation into SISL schemes resolves the numerical mismatch between spatial discretization and advection interpolation. The study demonstrates that this approach is not merely a theoretical improvement but offers tangible benefits:

• It eliminates the "accuracy gap" between the spectral derivative calculations and the interpolation steps.
• It reduces numerical diffusion, thereby preserving fine-scale structures and wave-like features over long integration periods.
• It enhances the conservation of mass and energy, which is critical for long-term climate simulations.
• It achieves these improvements without incurring prohibitive computational costs, thanks to the NUFFT acceleration.

The authors conclude that while low-order interpolation has been a practical compromise, the new spectral approach offers a superior alternative for high-accuracy, long-term simulations of atmospheric dynamics, motivating its extension to more complex three-dimensional dynamical cores.