Analysis and improvement of a semi-Lagrangian exponential scheme for the shallow-water equations on the rotating sphere

TL;DR

The work addresses stiff linear dynamics coupled with nonlinear advection in the shallow-water equations on the rotating sphere by analyzing and improving semi-Lagrangian exponential time integrators. It identifies the integration-factor discretization as the source of first-order accuracy in prior SL-ETDRK schemes and introduces a second-order SE22 method with a trapezoidal-like integration factor, supported by truncation-error and stability analyses and extensive numerical tests. Results show SE22 delivers improved stability and accuracy relative to Eulerian ETDRK and SL-SI-SETTLS, with competitive computational cost, suggesting potential suitability for operational atmospheric applications while highlighting ongoing challenges related to commutativity and viscosity. Overall, the paper provides a rigorous pathway to higher-order, stable time stepping for advection-dominated, stiff systems in geophysical flows, with clear guidance for practical deployment and future refinements.

Abstract

In this work, we study and extend a class of semi-Lagrangian exponential methods, which combine exponential time integration techniques, suitable for integrating stiff linear terms, with a semi-Lagrangian treatment of nonlinear advection terms. Partial differential equations involving both processes arise for instance in atmospheric circulation models. Through a truncation error analysis, we show that previously formulated semi-Lagrangian exponential schemes are limited to first-order accuracy due to the approximation of the integration factor acting on the discretization of the linear term; we then formulate a new discretization leading to second-order accuracy. Also, a detailed stability study is conducted to compare several Eulerian and semi-Lagrangian exponential schemes, as well as a well-established semi-Lagrangian semi-implicit method, which is used in operational atmospheric models. Numerical simulations of the shallow-water equations on the rotating sphere are performed to assess the orders of convergence, stability properties, and computational cost of each method. The proposed second-order semi-Lagrangian exponential method was shown to be more stable and accurate than the previously formulated schemes of the same class at the expense of larger wall-clock times; however, the method is more stable and has a similar cost compared to the well-established semi-Lagrangian semi-implicit method; therefore, it is a competitive candidate for potential operational applications in atmospheric circulation modeling.

Figures (12)

• Figure 1: Exact Lagrangian trajectory $\Gamma^{n+1}_j$ (in blue) in $[t_n, t_{n+1}]$, and estimated trajectory provided by SETTLS (in red), where $V^{n+1/2}_m$ is an estimation of the velocity field on its midpoint $(t_{n+1/2}, x_m)$ (see Eq. \ref{['eq:SETTLS']}).
• Figure 2: Stability regions (in red) of ETD1RK and ETD2RK applied to \ref{['eq:pde_eulerian']} (considering a nonlinear treatment of the advection term), and SL-SI-SETTLS, SE11 and SE12 applied to \ref{['eq:pde_SL']}, in the [id=R2-new]($\Im(\xiN),\Im(\xiL))$$\Im(\xiN)-\Im(\xiL)$ plane, considering both $\xiN$ and $\xiL$ to be purely imaginary.
• Figure 3: Stability functions of ETD1RK and ETD2RK applied to \ref{['eq:pde_eulerian']} (considering a nonlinear treatment of the advection term), and SL-SI-SETTLS, SE11 and SE12 applied to \ref{['eq:pde_SL']}, as a function of $\Im(\xiN)$ for fixed values of $\xiL$ and $\Re(\xiN) = 0$. The stability regions of ETD1RK and ETD2RK are empty in the case $\xiL = 0$. The horizontal, dashed line indicates the boundary of the stability region requiring $|A| \leq 1$. The same legend holds for both plots.
• Figure 4: Numerical verification of the convergence order of the semi-Lagrangian schemes considered in this work (relative $L_{2}$ error vs time step size), considering governing equations with only advective and linear terms.
• Figure 5: Numerical verification of the convergence order of the schemes considered in this work (normalized $\xLtwo$ norm of the error on the geopotential field perturbation vs the time step size). Similar results are obtained with the $\xLinfty$ norm.