CFL Optimized Forward-Backward Runge-Kutta Schemes for the Shallow Water Equations

We present the formulation and optimization of a Runge–Kutta-type time-stepping scheme for solving the shallow-water equations, aimed at substantially increasing the effective allowable time step over that of comparable methods. This scheme, called FB-RK(3,2), uses weighted forward–backward averaging of thickness data to advance the momentum equation. The weights for this averaging are chosen with an optimization process that employs a von Neumann–type analysis, ensuring that the weights maximize the admittable Courant number. Through a simplified local truncation error analysis and numerical experiments, we show that the method is at least second-order in time for any choice of weights and exhibits low dispersion and dissipation errors for well-resolved waves. Further, we show that an optimized FB-RK(3,2) can take time steps up to 2.8 times as large as a popular three-stage, third-order strong stability-preserving Runge–Kutta method in a quasi-linear test case. In fully nonlinear shallow-water test cases relevant to oceanic and atmospheric flows, FB-RK(3,2) outperforms SSPRK3 in admittable time step by factors roughly between 1.6 and 2.2, making the scheme approximately twice as computationally efficient with little to no effect on solution quality. Significance Statement: The purpose of this work is to develop and optimize time-stepping schemes for models relevant to oceanic and atmospheric flows. Specifically, for the shallow-water equations we optimize for schemes that can take time steps as large as possible while retaining solution quality. We find that our optimized schemes can take time steps between 1.6 and 2.2 times larger than schemes that cost the same number of floating point operations, translating directly to a corresponding speedup. Our ultimate goal is to use these schemes in climate-scale simulations.

Keywords: Shallow-water equations; Numerical analysis/modeling; Model evaluation/performance

Large-scale climate models, such as the U.S. Department of Energy's Energy Exascale Earth System Model (E3SM) ([13]), grow in complexity every year as new numerical techniques are developed and as high-performance computing (HPC) systems become more powerful. Complex, climate-scale simulations can take days or weeks or more of real time to complete; it is vital that these run times are controlled without making undue sacrifices to overall model quality. In the present work, we develop a way to formulate and optimize a class of time-stepping schemes for the shallow-water equations (SWEs). Our ultimate goal is to improve the efficiency of simulations at the climate scale, but here we focus on the shallow-water equations as a starting point. Although this work is directed to models of ocean and atmosphere circulation, the ideas and methods developed here may generalize to other coupled systems.

The central issue here is the following. The Courant–Friedrichs–Lewy (CFL) condition is a well-known necessary condition for the stability of a given model under an explicit time-stepping scheme. In essence, the CFL condition places a restriction on the size of the model's time step, where this restriction depends on the model problem itself and the chosen time and space discretizations. Formally, it is necessary for stability that

(1) $\nu =c\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\le {\nu }^{\text{max}},$

where c is a speed such as a gravity wave speed. We refer to ν as the Courant number and νmax as the maximal Courant number for the given situation. In this work, we develop time-stepping schemes optimized to have maximal νmax. We do this by starting with a three-stage second-order Runge–Kutta method from [36], presently used in the Model for Prediction Across Scales-Atmosphere (MPAS-A) ([33]). This scheme is augmented by adding forward–backward averaging of the thickness tendencies used to update the fluid velocity. Resulting schemes have large νmax and good dispersion and dissipation relations, and they compare well against other schemes in both solution accuracy and effective maximal Courant number.

The SWEs serve as a standard test problem for algorithms and methods designed for applications to ocean and atmosphere modeling; in both cases, the model equations used for climate-scale simulations are generalizations of the SWEs, in which additional vertical tendencies and forcing terms are included, and variable density is allowed. In recent years, strong stability preserving Runge–Kutta (SSPRK) methods ([14]) for solving the SWEs have become quite popular; they have the property that at each time level, a spatial norm of the solution is less than or equal to the norm of the solution at the preceding level. These methods are used to solve the SWEs in a variety of numerical paradigms and applications, such as with finite element spatial discretizations ([2]; [7]; [11]), finite volume discretizations ([3]; [4]; [15]; [17]; [23]; [30]; [34]), and layered models that use a time splitting of barotropic and baroclinic motions ([16]; [19]). Despite this widespread use, strong stability preserving (SSP) time-integration methods may not be the best choice for most large-scale geophysical SWE applications. The SSP property is useful in systems where strong shocks and discontinuities can develop and propagate, but such dynamics are rare in climate-scale oceanic and atmospheric flows, in which barotropic evolution is typically characterized by Rossby- and Kelvin-type wave modes. Numerous contemporary large-scale ocean and atmosphere models employ time-stepping methods that do not have the SSP property, motivating an investigation of its optimality for this class of flow. For example, to solve the barotropic equations in split barotropic/baroclinic systems, the Model for Prediction Across Scales-Ocean (MPAS-O) ([29]; [26]) and Model for Prediction Across Scales-Atmosphere (MPAS-A) ([33]) use a two-stage predictor-corrector method and a three-stage second-order Runge–Kutta method, respectively. The Community Atmosphere Model (CAM) ([6]), the Parallel Ocean Program (POP) ([18]), the Nucleus for European Modeling of the Ocean (NEMO) ([24]), and the Hybrid Coordinate Ocean Model (HYCOM) ([5]) all use variations of semi-implicit leapfrog time integration schemes for their barotropic systems. The Regional Ocean Modeling System (ROMS) ([31]) uses a multistep forward–backward Adams–Bashforth/Adams–Moulton method.

With the potential performance benefits of the ability to take larger time steps in mind, we turn our attention toward methods that are optimal in this sense. Taking advantage of the structure of the SWEs, we employ a strategy similar to that used by [25], wherein the thickness equation is advanced in time with a simple forward Euler step, then the momentum equation is advanced with an explicit backward Euler step using the recently obtained data for the thickness. The method, referred to by Mesinger as the forward–backward method, was more computationally efficient than popular leapfrog methods of the time. The idea of advancing the variables of a system of ODE separately, using the most recently computed data from previously advanced variables, has been employed by the developers of ROMS. Specifically, [31], [32]) present an extension of a two-stage second-order Runge–Kutta scheme that, within each RK stage, first advances the thickness variable, then uses a weighted average of the most recently computed data with previously computed data to advance the momentum equation in a shallow-water system. The coefficients of this scheme are selected to optimize the allowable time-step length for a one-dimensional gravity wave system based on a von Neumann–type stability analysis. In this work, we develop a new, forward–backward weighted three-stage Runge–Kutta method with favorable CFL limits using a similar approach. The optimal weights for our forward–backward scheme are obtained by von Neumann–type analysis on the two-dimensional system, including a linear Coriolis term and linearization of the advection operators about a nonzero mean state. Using a numerical optimization algorithm, we obtain a scheme that is very efficient in terms of its maximal Courant number, as well as its dissipation and dispersion characteristics.

The paper is structured as follows. We begin by introducing our new time-stepping scheme, referred to as FB-RK(3,2), for a general system of ODEs. Then, we describe the process by which we obtain optimal weights for the FB averaging within each RK stage. This process includes the derivation of a von Neumann formulation of the linearized SWEs with a nonzero mean flow and the formulation of an optimization problem. Next, we present the results of this optimization and give a brief discussion of the dispersion and dissipation errors of the optimal FB-RK(3,2) schemes. Finally, we perform a number of numerical experiments that demonstrate the computational efficiency and efficacy of the optimal FB-RK(3,2) schemes as compared to the three-stage third-order SSPRK scheme (SSPRK3) from [14]. These experiments include an investigation of CFL performance across five shallow-water test cases, convergence tests, and comparisons of solution quality. We show that an optimized FB-RK(3,2) scheme can take time steps up to 2.2 times as large as SSPRK3 in nonlinear test cases while retaining solution quality, and because both schemes consist of three RK stages, this translates directly to a corresponding increase in efficiency.

Here, we introduce a time-stepping scheme for solving the shallow-water equations (SWEs) that uses weighted forward–backward averaging in the layer thickness data to advance the momentum equation. We then use a von Neumann analysis and a numerical optimization scheme to tune the forward–backward weights so that the scheme has a large Courant number when applied to the linearized SWEs.

The time-stepping scheme presented here is an extension of the three-stage, second-order Runge–Kutta time-stepping scheme RK(3,2) from [36]. RK(3,2) was originally presented as a third-order method and referred to as RK3, but it is shown by [27] that the method is third-order only on linear problems, reducing to second-order on nonlinear problems. When applied to the SWEs, our extension of RK(3,2) allows us to use the most recently obtained data for the layer thickness to update the momentum data within each Runge–Kutta stage of the method. This is done by taking a weighted average of layer thickness data at the old time level tn and the most recent RK stage, then applying this to the momentum equation. For the remainder of the manuscript, we will call our extension of RK(3,2) that uses forward–backward averaging FB-RK(3,2), as it is a three-stage, second-order scheme.

Consider a general system of ODEs in independent variables u = u(t) and h = h(t) of the following form:

(2) $\begin{array}{l}\frac{du}{dt}=\mathrm{\Phi }(u,h),\\ \frac{dh}{dt}=\mathrm{\Psi }(u,h),\end{array}$

where t is the time coordinate. Let un ≈ u(tn) and hn ≈ h(tn) $h^n\approx h(t^n)$ be the numerical approximations to u and h at time t = tn. Let Δt be a time step such that tn+1 = tn + Δt. Also, let $t^{n+1/m}=t^n+(\mathrm{\Delta }t/m)$ for any positive integer m. Then, FB-RK(3,2) is given by

(3a) $\begin{array}{ll}{\overline{h}}^{n+1/3}\hfill & =h^n+\frac{\mathrm{\Delta }t}{3}\mathrm{\Psi }(u^n,h^n),\hfill \\ {\overline{u}}^{n+1/3}\hfill & =u^n+\frac{\mathrm{\Delta }t}{3}\mathrm{\Phi }(u^n,h^{*}),\hfill \\ h^{*}\hfill & ={\beta }_1{\overline{h}}^{n+1/3}+(1-{\beta }_1)h^n,\hfill \end{array}$

(3b) $\begin{array}{ll}{\overline{h}}^{n+1/2}\hfill & =h^n+\frac{\mathrm{\Delta }t}{2}\mathrm{\Psi }({\overline{u}}^{n+1/3},{\overline{h}}^{n+1/3}),\hfill \\ {\overline{u}}^{n+1/2}\hfill & =u^n+\frac{\mathrm{\Delta }t}{2}\mathrm{\Phi }({\overline{u}}^{n+1/3},h^{**}),\hfill \\ h^{**}\hfill & ={\beta }_2{\overline{h}}^{n+1/2}+(1-{\beta }_2)h^n,\hfill \end{array}$

(3c) $\begin{array}{ll}{h}^{n+1}\hfill & =h^n+\mathrm{\Delta }t\mathrm{\Psi }({\overline{u}}^{n+1/2},{\overline{h}}^{n+1/2}),\hfill \\ u^{n+1}\hfill & =u^n+\mathrm{\Delta }t\mathrm{\Phi }({\overline{u}}^{n+1/2},h^{***}),\hfill \\ h^{***}\hfill & ={\beta }_3{h}^{n+1}+(1-2{\beta }_3){\overline{h}}^{n+1/2}+{\beta }_3{h}^n\cdot \hfill \end{array}$

The weights β1, β2, and β3 are called the forward–backward (FB) weights. The best choice for the FB weights βi is a primary concern of this work.

The choice to calculate data at time tn+1/3 and tn+1/2 in the first and second stages comes from the original RK(3,2) scheme; the treatment of the thickness variable h in (3) shows the formulation of RK(3,2). FB-RK(3,2) can be approximately reduced to RK(3,2) by taking β1 = 0, ${\beta }_2=2/3$ , and β3 = 0. This produces an approximation to RK(3,2) as, while the first and third stages are exact, the FB averaged $h^{**}=(2/3){\overline{h}}^{n+1/2}+(1/3)h^n$ data in the second stage is an approximation to ${\overline{h}}^{n+1/3}$ . The FB averages are chosen in this way so that in each stage, data are symmetrically distributed in time as preliminary experiments showed that this provided the best CFL performance. In the first two RK stages, the FB averaging is done on data from the beginning and the end of the current time interval, i.e., in the first stage, we compute values for h and u at time tn+1/3, so the FB averaging is done with data at time tn and time tn+1/3. In the third RK stage, we use data from the beginning, middle, and end of the current time interval.

In the context of this work, h would be the fluid layer thickness and u the fluid velocity. We have chosen to use an FB average of thickness to advance the velocity variable. However, one could reasonably define a scheme similar to that given by (3) by taking FB averages of u to advance h instead. The subsequent analysis would remain largely the same, but the optimal choices for the FB weights could change substantially, along with the performance of the scheme as a whole. What happens in this case is an open question.

As stated above, the original RK(3,2) scheme is $\mathcal{O}[{(\mathrm{\Delta }t)}^2]$ . For any choice of the FB weights, FB-RK(3,2) remains $\mathcal{O}[{(\mathrm{\Delta }t)}^2]$ while greatly increasing the maximum allowable time step for most problems. A generalized local truncation error analysis is outside the scope of this paper, so we show this holds for a particular problem of interest in section a of the appendix.

A primary concern of this work is the application of FB-RK(3,2) to the nonlinear SWEs on a rotating sphere. Application of a von Neumann–type analysis requires a linear set of PDEs, and we consider the two-dimensional system with linear Coriolis and linearized advection operators to ensure the time-step optimization is conducted on a system that is a close approximation to the fully nonlinear equations. The linearized SWEs centered about a nonzero mean flow are given by

(4) $\begin{array}{l}\frac{\partial \tilde{u}}{\partial t}+f{(\tilde{U}+\tilde{u})}^{\perp }+\tilde{\zeta }{\tilde{U}}^{\perp }=-\tilde{U}\cdot \nabla \tilde{u}-g\nabla \tilde{h},\\ \frac{\partial \tilde{h}}{\partial t}+H(\nabla \cdot \tilde{u})+\nabla \cdot (\tilde{h}\tilde{U})=0,\end{array}$

where $\tilde{u}=[\tilde{u}(x,y,t),\tilde{\upsilon }(x,y,t)]$ is a perturbation of the horizontal fluid velocity, x and y are the spatial coordinates, t is the time coordinate, f = 2Ω sinφ is the Coriolis parameter, where Ω is the angular velocity of the rotating sphere and φ denotes latitude, $\tilde{U}=(\tilde{U},\tilde{V})$ is the constant horizontal mean fluid velocity, ${\tilde{u}}^{\perp }=(-\tilde{\upsilon },\tilde{u})$ , $\tilde{\zeta }=(\partial \tilde{\upsilon }/\partial x)-(\partial \tilde{u}/\partial y)$ is the vertical component of vorticity $\nabla \times \tilde{u}$ , g is the gravitational constant, $\tilde{h}=\tilde{h}(x,y,t)$ is a perturbation of the layer thickness, and H is the (constant in space) fluid layer thickness when at rest. Note that while $\tilde{u}$ and $\tilde{h}$ depend on x and y throughout this work, we will often suppress this dependence in the notation. The full derivation of these equations is given in section b of the appendix.

We can write these in terms of dimensionless quantities by writing $c=\sqrt{gH},\mathbf{u}=(u,\upsilon )=(\tilde{u}/c),\mathbf{U}=(\tilde{U}/c)$ , and $\eta =\tilde{h}/H$ . Then, the dimensionless, linearized SWEs centered about a nonzero mean flow are given by

(5) $\begin{array}{l}\frac{\partial \mathbf{u}}{\partial t}+f{(\mathbf{U}+\mathbf{u})}^{\perp }+c\zeta {\mathbf{U}}^{\perp }=-c\mathbf{U}\cdot \nabla \mathbf{u}-c\nabla \eta ,\\ \frac{\partial \eta }{\partial t}+c(\nabla \cdot \mathbf{u})+c[\nabla \cdot (\eta \mathbf{U})]=0.\end{array}$

In the analysis that follows, we assume that the system (5) is discretized on a rectangular, staggered Arakawa C-grid, where η is computed at cell centers, u is computed at the left and right edges of cells, and υ is computed at the top and bottom edges of cells ([1]). For brevity, we skip some details of calculations relating to the spatial discretization in the main text and instead provide them in section c of the appendix.

With the goal of using a von Neumann–type analysis to obtain an optimal choice of FB weights, we apply a Fourier transform in x and y to the linearized SWEs (5). This is equivalent to seeking solutions of the following form:

$\begin{array}{l}u(x,y,t)=\widehat{u}(k,\mathcal{l},t)e^{\mathit{ikx}+i\mathcal{l}y},\\ \upsilon (x,y,t)=\widehat{\upsilon }(k,\mathcal{l},t)e^{\mathit{ikx}+i\mathcal{l}y},\\ \eta (x,y,t)=\widehat{\eta }(k,\mathcal{l},t)e^{\mathit{ikx}+i\mathcal{l}y},\end{array}$

where k is the wavenumber with respect to x, and $\mathcal{l}$ is the wavenumber with respect to y. In doing this, spatial derivatives are now all of the form $(\partial /\partial x)(e^{\mathit{ikx}+i\mathcal{l}y})$ or $(\partial /\partial y)(e^{\mathit{ikx}+i\mathcal{l}y})$ . Assuming a rectangular C-grid, using centered differences these spatial derivatives become

(6) $\frac{\partial }{\partial x}(e^{\mathit{ikx}+i\mathcal{l}y})=\frac{i2\sin \left(k\frac{\mathrm{\Delta }x}{2}\right)}{\mathrm{\Delta }x}{e}^{\mathit{ikx}+i\mathcal{l}y},$

(7) $\frac{\partial }{\partial y}(e^{\mathit{ikx}+i\mathcal{l}y})=\frac{i2\sin \left(\mathcal{l}\frac{\mathrm{\Delta }y}{2}\right)}{\mathrm{\Delta }y}{e}^{\mathit{ikx}+i\mathcal{l}y}.$

To compute the Coriolis terms, we need the values of u at υ points and the values of υ at u points; call these uυ and υu, respectively. On a rectangular C-grid, we achieve this with unweighted four-point averages of the nearest data, which gives

(8) $u_{\upsilon }(x,y,t)=\cos \left(k\frac{\mathrm{\Delta }x}{2}\right)\cos \left(\mathcal{l}\frac{\mathrm{\Delta }y}{2}\right)\widehat{u}(k,\mathcal{l},t)e^{\mathit{ikx}+i\mathcal{l}y},$

(9) ${\upsilon }_u(x,y,t)=\cos \left(k\frac{\mathrm{\Delta }x}{2}\right)\cos \left(\mathcal{l}\frac{\mathrm{\Delta }y}{2}\right)\widehat{\upsilon }(k,\mathcal{l},t)e^{\mathit{ikx}+i\mathcal{l}y}.$

The details of these calculations are given in section c of the appendix. Set $K=2\sin [k(\mathrm{\Delta }x/2)],L=2\sin [\mathcal{l}(\mathrm{\Delta }y/2)]$ , and $\psi =f\cos [k(\mathrm{\Delta }x/2)]\cos [\mathcal{l}(\mathrm{\Delta }y/2)]$ . Then, accounting for the Fourier transform and the C-grid spatial discretization, (5) becomes

(10) $\begin{array}{l}\frac{\partial \widehat{u}}{\partial t}=fV+\psi \widehat{\upsilon }-\left(iU\frac{cK}{\mathrm{\Delta }x}+iV\frac{cL}{\mathrm{\Delta }y}\right)\widehat{u}-i\frac{cK}{\mathrm{\Delta }x}\widehat{\eta },\\ \frac{\partial \widehat{\upsilon }}{\partial t}=-fU-\psi \widehat{u}-\left(iU\frac{cK}{\mathrm{\Delta }x}+iV\frac{cL}{\mathrm{\Delta }y}\right)\widehat{\upsilon }-i\frac{cL}{\mathrm{\Delta }y}\widehat{\eta },\\ \frac{\partial \widehat{\eta }}{\partial t}=-\left(i\frac{cK}{\mathrm{\Delta }x}\widehat{u}+i\frac{cL}{\mathrm{\Delta }y}\widehat{\upsilon }\right)-\left(iU\frac{cK}{\mathrm{\Delta }x}+iV\frac{cL}{\mathrm{\Delta }y}\right)\widehat{\eta }\cdot \end{array}$

Let ${\nu }_x=c(\mathrm{\Delta }t/\mathrm{\Delta }x)$ and ${\nu }_y=c(\mathrm{\Delta }t/\mathrm{\Delta }y)$ be the Courant numbers with respect to x and y, respectively. Also, set φ = Δtψ. Then, we apply FB-RK(3,2) to (10):

(11a) $\begin{array}{l}{\overline{\widehat{\eta }}}^{n+1/3}={\widehat{\eta }}^n-\frac{1}{3}\left[iK{\nu }_x{\widehat{u}}^n+iL{\nu }_y{\widehat{\upsilon }}^n+(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\widehat{\eta }}^n\right],\\ {\overline{\widehat{u}}}^{n+1/3}={\widehat{u}}^n+\frac{1}{3}[\mathrm{\Delta }\mathit{tfV}+\phi {\widehat{\upsilon }}^n-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\widehat{u}}^n-iK{\nu }_x{\widehat{\eta }}^{*}],\\ {\overline{\widehat{\upsilon }}}^{n+1/3}={\widehat{\upsilon }}^n+\frac{1}{3}[-\mathrm{\Delta }\mathit{tfU}-\phi {\widehat{u}}^n-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\widehat{\upsilon }}^n-iK{\nu }_x{\widehat{\eta }}^{*}],\\ {\widehat{\eta }}^{*}={\beta }_1{\overline{\widehat{\eta }}}^{n+1/3}+(1-{\beta }_1){\widehat{\eta }}^n,\end{array}$

(11b) $\begin{array}{l}{\overline{\widehat{\eta }}}^{n+1/2}={\widehat{\eta }}^n-\frac{1}{2}[iK{\nu }_x{\overline{\widehat{u}}}^{n+1/3}+iL{\nu }_y{\overline{\widehat{\upsilon }}}^{n+1/3}+(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{\eta }}}^{n+1/3}],\\ {\overline{\widehat{u}}}^{n+1/2}={\widehat{u}}^n+\frac{1}{2}[\mathrm{\Delta }\mathit{tfV}+\phi {\overline{\widehat{\upsilon }}}^{n+1/3}-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{u}}}^{n+1/3}-iK{\nu }_x{\widehat{\eta }}^{**}],\\ {\overline{\widehat{\upsilon }}}^{n+1/2}={\widehat{\upsilon }}^n+\frac{1}{2}[-\mathrm{\Delta }\mathit{tfU}-\phi {\overline{\widehat{u}}}^{n+1/3}-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{\upsilon }}}^{n+1/3}-iK{\nu }_x{\widehat{\eta }}^{**}],\\ {\widehat{\eta }}^{**}={\beta }_2{\overline{\widehat{\eta }}}^{n+1/2}+(1-{\beta }_2){\widehat{\eta }}^n,\end{array}$

(11c) $\begin{array}{l}{\widehat{\eta }}^{n+1}={\widehat{\eta }}^n-\left[iK{\nu }_x{\overline{\widehat{u}}}^{n+1/2}+iL{\nu }_y{\overline{\widehat{\upsilon }}}^{n+1/2}+(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{\eta }}}^{n+1/2}\right],\\ {\widehat{u}}^{n+1}={\widehat{u}}^n+\left[\mathrm{\Delta }\mathit{tfV}+\phi {\overline{\widehat{\upsilon }}}^{n+1/2}-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{u}}}^{n+1/2}-iK{\nu }_x{\widehat{\eta }}^{***}\right],\\ {\widehat{\upsilon }}^{n+1}={\widehat{\upsilon }}^n+[-\mathrm{\Delta }\mathit{tfU}-\phi {\overline{\widehat{u}}}^{n+1/2}-(\mathit{iUK}{\nu }_x+\mathit{iVL}{\nu }_y){\overline{\widehat{\upsilon }}}^{n+1/2}-iK{\nu }_x{\widehat{\eta }}^{***}],\\ {\widehat{\eta }}^{***}={\beta }_3{\widehat{\eta }}^{n+1}+(1-2{\beta }_3){\overline{\widehat{\eta }}}^{n+1/2}+{\beta }_3{\widehat{\eta }}^n.\end{array}$

Note that (11) forms a linear system in ${\widehat{u}}^n,{\widehat{\upsilon }}^n$ , and ${\widehat{\eta }}^n$ , that is, there is a matrix $\mathsf{G}$ and column vector b such that

(12) ${\widehat{w}}^{n+1}=\mathsf{G}{\widehat{w}}^n+\mathbf{b},$

where ${\widehat{w}}^n={({\widehat{u}}^n,{\widehat{\upsilon }}^n,{\widehat{\eta }}^n)}^{\mathrm{T}}$ . It follows from (12) that

(13) ${\widehat{w}}^n={\mathsf{G}}^n{\widehat{w}}^0+(\mathsf{I}-{\mathsf{G}}^n){(\mathsf{I}-\mathsf{G})}^{-1}\mathbf{b},$

where $\mathsf{I}$ is the identity matrix. Note that in the above, the superscript on $\widehat{w}$ is an index on the time step, whereas the superscript on $\mathsf{G}$ is an exponent. The matrix $\mathsf{G}$ is called the amplification matrix; properties of von Neumann analysis tell us that a time-stepping method is stable if the eigenvalues of its amplification matrix reside within in the unit circle centered at the origin in the complex plane ([20]). This property is fundamental to the analysis in section 2d.

With FB-RK(3,2) applied to the von Neumann formulation of the linearized SWEs (11) and the corresponding amplification matrix $\mathsf{G}$ from (12), we can more clearly state the goal of the analysis presented in this section. We are interested in selecting the FB weights β1, β2, and β3 in FB-RK(3,2) so that the Courant numbers ${\nu }_x=c(\mathrm{\Delta }t/\mathrm{\Delta }x)$ and ${\nu }_y=c(\mathrm{\Delta }t/\mathrm{\Delta }y)$ can be taken as large as possible while retaining stability of the method.

The eigenvalues of the amplification matrix $\mathsf{G}$ depend on the Courant numbers νx and νy and on β1, β2, and β3. To simplify the analysis, we will assume that νx = νy = ν is the Courant number (this is equivalent to assuming that the spatial discretization is on a grid consisting of sqaure cells). Then, we can write $\mathsf{G}$ = $\mathsf{G}$ (ν, β1, β2, β3). For a given choice of β1, β2, and β3, let ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ be the largest value of the Courant number for which FB-RK(3,2) is stable, i.e., largest value of the Courant number for which the eigenvalues of the amplification matrix $\mathsf{G}$ reside within the unit circle. Symbolically, ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ is the largest number such that

(14) $\left|{\lambda }_{\text{max}}[\mathsf{G}(\nu ,{\beta }_1,{\beta }_2,{\beta }_3)]\right|\le 1,\hspace{0.17em}\text{for}\hspace{0.17em}\text{all}\hspace{0.17em}\nu \in \left(0,{\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}\right],$

where λmax( $\mathsf{G}$ ) is the largest eigenvalue, in absolute value, of the matrix $\mathsf{G}$ .

While the matrix $\mathsf{G}$ is too complex to analyze directly, we can use numerical optimization tools to find optimal FB weights. That is, we can formulate this problem explicitly as an optimization problem with an appropriate cost function and constraints, then use numerical optimization algorithms to search the space of admittable FB weights. We can maximize νmax by minimizing an appropriate cost function:

(15) $C=C({\beta }_1,{\beta }_2,{\beta }_3),$

subject to the constraint that the eigenvalues of $\mathsf{G}$ satisfy (14).

The particular numerical optimization scheme we employ here is the simplicial homology global optimization (SHGO) algorithm developed by [8], implemented in the Python programming language via the open-source Python library SciPy ([35]). SHGO is suitable for so-called black-box optimization problems wherein the cost function can be nonsmooth and discontinuous. Further, SHGO is guaranteed to converge to a global optimum as the number of iterations increases. As such, SHGO is an appropriate choice for our purposes, where the eigenvalues of the matrix $\mathsf{G}$ are sensitive to small changes in the FB weights.

The matrix $\mathsf{G}$ also depends on additional parameters Δtf, kΔx, and $\mathcal{l}\mathrm{\Delta }y$ as can be seen in (11). For the purposes of the numerical optimizations we perform here, we set Δtf = 10−2 since taking f = 10−4 represents the value of the Coriolis parameter at midlatitudes and we expect a FB-RK(3,2) to be capable of a time step on the order of 102 in the nonlinear SWE test cases described in section 3. Preliminary tests showed that the choice of value for Δtf had little to no effect on the optimizations in practice. We take $k\mathrm{\Delta }x=\mathcal{l}\mathrm{\Delta }y=\pi$ , which corresponds to the case of gridscale waves. This choice is informed by the fact that instability often comes from gridscale processes. That is, we perform the optimization on gridscale waves, which represents a worst-case with the hope that this provides the best overall stability.

Given a particular choice of β1, β2, and β3, we can easily estimate the value of ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ to arbitrary precision. Our goal is to maximize this parameter, so we seek to minimize the cost function:

(16) $C_1({\beta }_1,{\beta }_2,{\beta }_3)=\frac{1}{{\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}}.$

In the special case where the mean flow U = (U, V) is zero, we additionally consider a second choice of cost function. When U = 0 the column vector b in (12) is also zero, so (13) simplifies to

(17) ${\widehat{w}}^n={\mathsf{G}}^n{\widehat{w}}^0.$

In this case, it is easy to show that the numerical amplification matrix $\mathsf{G}$ approximates $\tilde{\mathsf{G}}:=e^{\mathsf{A}\mathrm{\Delta }t}$ , where

(18) $\mathsf{A}=\left(\begin{array}{ccc}0& f& -\mathit{ick}\\ -f& 0& -ic\mathcal{l}\\ -\mathit{ick}& -ic\mathcal{l}& 0\end{array}\right).$

The matrix $\tilde{\mathsf{G}}$ is the analytical counterpart to the numerical amplification matrix $\mathsf{G}$ that arises from solving a spatially continuous version of (10). In short, if we assume that our time-stepping method, which is determined by the matrix $\mathsf{G}$ , produces the exact solution at time t = Δt, then $\mathsf{G}=\tilde{\mathsf{G}}$ . We give the details of this in appendix, section d.

Since we wish for $\mathsf{G}$ to approximate $\tilde{\mathsf{G}}$ for a given choice of FB weights, we aim to minimize a second choice of cost function:

(19) $C_2({\beta }_1,{\beta }_2,{\beta }_3)=\frac{1}{{\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}}+{\displaystyle {\int }_0^{\pi /6}\Vert \tilde{\mathsf{G}}(\nu )-\mathsf{G}(\nu ,{\beta }_1,{\beta }_2,{\beta }_3)\Vert }d\nu ,$

where ∥⋅∥ is the Frobenius norm:

(20) $\Vert \mathsf{A}\Vert =\sqrt{{\displaystyle {\sum }_{i,j}{(a_{ij})}^2}}.$

This choice of cost function seeks to maximize νmax while also producing a scheme that has favorable dissipation and dispersion behavior. The choice of norm in the integral term of C2 was chosen as it is quick to compute during the optimization process. There is a certain freedom in the choice of limits of integration for the integral term in C2; this term requires that the difference between $\mathsf{G}$ and $\tilde{\mathsf{G}}$ is minimal as for small values of ν, which correspond to well-resolved waves. In essence, while C1 only considers the CFL performance of the resulting method, C2 does this while also considering that $\mathsf{G}$ should accurately approximate its analytical counterpart $\tilde{\mathsf{G}}$ .

In summary, we consider two choices of cost function given by C1 in (16) and C2 in (19). The function C1 is valid in any case, while C2 applies only when the mean flow U is taken to be zero.

In the context of the dimensionless Eq. (5), the parameter |U| is called the Froude number. We perform the optimization described above five total times, considering different values of the Froude number |U| that appear in oceanic and atmospheric flows.

As noted above, the SHGO algorithm used here is guaranteed to converge to a globally optimal solution if run for enough iterations. As our problem is relatively complex and computing resources finite, it is more likely in practice that the method will convergence to a local extremum. For example, consider the first two rows of Table 1 and notice that the value of ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ given by optimizing C2 is greater than that given by optimizing C1. If these were globally optimal solutions, we would have that ${\nu }_{(0.500,0.500,0.344)}^{\text{max}}\ge {\nu }_{(0.516,0.532,0.331)}^{\text{max}}$ since C1 consists only of a term that asks for a maximal ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ , while C2 is the sum of the same term with a nonnegative term that asks for $\mathsf{G}$ to approximate $\tilde{\mathsf{G}}$ . Despite the problem being sufficiently complex that globally optimal solutions are difficult to obtain, the results from section 3 show that the choices of FB weights from Table 1 are extremely efficient in practice.

Graph: Table 1. Results of optimizing FB-RK(3,2) for the given cost functions, rounded to three decimal places. Note that the reported values of ν(β1,β2,β3)max are the maximal Courant numbers in the context of a problem consisting only of gridscale waves, i.e., the problem for which FB-RK(3,2) was optimized.

A similar Runge–Kutta-based method using forward–backward averaging in the momentum equation is developed by [31] (from here on denoted ShMc). In ShMc, the efficiency of their optimized FB method is reported in terms of the parameter α := ckΔt, by considering the maximal value αmax for which the method is stable. Translating this into the notation used here, we obtain

(21) $\begin{array}{l}\alpha =ck\mathrm{\Delta }t\\ =\nu k\mathrm{\Delta }x.\end{array}$

For their optimal RK2-based FB method, ShMc reports that (νkΔx)max ≈ 2.141; since their scheme uses two RK stages, this translates to what we refer to as an effective CFL number of 2.141/2 ≈ 1.071. Here, we have taken kΔx = π and have νmax = 1.804 in the best case, so (νkΔx)max = 1.804π ≈ 5.667. FB-RK(3,2) uses three RK stages, so we obtain an effective CFL number of 5.667/3 ≈ 1.889, an increase of approximately 1.764 times over ShMc's RK2 FB scheme. Note that ShMc's (νkΔx)max comes from their FB scheme as applied to a one-dimensional wave system analogous to (22), whereas ours is obtained from the two-dimensional linearized SWEs. The comparison between these two values is not one-to-one, but do serve as a general comparison of the efficiency of both methods.

To visualize FB-RK(3,2)'s treatment of gravity waves within the SWEs, we can apply the scheme to a simple wave system and analyze the resulting dissipation and dispersion errors (Fig. 1). Applying a Fourier transform in space to the 1D linearized SWEs (A1) gives

(22) $\begin{array}{l}\frac{\partial \widehat{\eta }}{\partial t}=-\mathit{ick}\widehat{u}\\ \frac{\partial \widehat{u}}{\partial t}=-\mathit{ick}\widehat{\eta }\cdot \end{array}$

Choose a time step Δt and let $\tilde{K}\hspace{0.17em}=k\mathrm{\Delta }x$ and $\nu =c(\mathrm{\Delta }t/\mathrm{\Delta }x)$ . The FB-RK(3,2) scheme for this system is then given by

(23a) $\begin{array}{ll}{\overline{\eta }}^{n+1/3}\hfill & ={\eta }^n-\frac{i\tilde{K}\nu }{3}{u}^n,\hfill \\ {\overline{u}}^{n+1/3}\hfill & =u^n-\frac{i\tilde{K}\nu }{3}[{\beta }_1{\overline{\eta }}^{n+1/3}+(1-{\beta }_1){\eta }^n],\hfill \end{array}$

(23b) $\begin{array}{ll}{\overline{\eta }}^{n+1/2}\hfill & ={\eta }^n-\frac{i\tilde{K}\nu }{2}{\overline{u}}^{n+1/3},\hfill \\ {\overline{u}}^{n+1/2}\hfill & =u^n-\frac{i\tilde{K}\nu }{2}[{\beta }_2{\overline{\eta }}^{n+1/2}+(1-{\beta }_2){\eta }^n],\hfill \end{array}$

(23c) $\begin{array}{ll}{\eta }^{n+1}\hfill & ={\zeta }^n-i\tilde{K}\nu {\overline{u}}^{n+1/2},\hfill \\ u^{n+1}\hfill & =u^n-i\tilde{K}\nu [{\beta }_3{\eta }^{n+1}+(1-2{\beta }_3){\overline{\eta }}^{n+1/2}+{\beta }_3{\eta }^n].\hfill \end{array}$

Equation (23) defines a system of linear equations, which in turn defines a numerical amplification matrix $\mathsf{G}$ 0 similar to that in (12). One can show that the general solution to the system (22) is

(24) $\left[\begin{array}{c}\eta (t)\\ u(t)\end{array}\right]=\exp \left[t\left(\begin{array}{cc}0& -\mathit{ick}\\ -\mathit{ick}& 0\end{array}\right)\right]\left[\begin{array}{c}\eta (0)\\ u(0)\end{array}\right],$

and that eigenvalues of the matrix exponential above are λ = e±ickt. At time t = Δt, these eigenvalues are $\lambda =e^{\pm i\tilde{K}\nu }$ .

Graph: Fig. 1. An illustration of dissipation and dispersion errors for FB-RK(3,2) with each choice of FB weights from Table 1 and SSPRK3. Eigenvalues λ=e±iK˜ν of the matrix in the exact solution (24) are given by orange dots. Drawing a ray from the origin to one of these dots and measuring the angle that ray makes with the positive real axis gives the corresponding value of K˜ν. In other words, K˜ν is the argument of the corresponding complex number. Eigenvalues of the numerical amplification matrix G 0 as K˜ν =kΔxν increases are given by the blue curve and the blue stars. Dashed lines connect exact eigenvalues to numerically determined counterparts, which gives a measure of dissipation and dispersion errors. Values of K˜ν close to zero correspond to well-resolved waves, while values approaching π correspond to temporal gridscale waves.

Note that although we have used the Courant number $\nu =c(\mathrm{\Delta }t/\mathrm{\Delta }x)$ in this discussion in an effort to use notation analogous to the rest of this section, the system considered here is not discretized in space; we have that $\tilde{K}\nu =k\mathrm{\Delta }x\times c(\mathrm{\Delta }t/\mathrm{\Delta }x)=ck\mathrm{\Delta }t$ . Values of $\tilde{K}\nu$ close to zero correspond to well-resolved waves, while values approaching π correspond to temporal gridscale waves. Dissipation error can be read as the distance a numerical eigenvalue is from the unit circle, resulting in the corresponding wave being damped out of the system. Dispersion error can be read as the difference in the angles made with the positive real axis between a numerically determined eigenvalue with its analytical counterpart, resulting in an error in the phase speed of the corresponding wave. For all choices of FB weights in Table 1, dissipation and dispersion errors for well resolved waves are very low (Fig. 1). As waves become less well resolved, each of the five schemes exhibits different behavior. In particular, note the behavior of the eigenvalues in Fig. 1c; well resolved waves are handled well as expected, and as waves move toward grid scale, the dissipation errors become very high, i.e., the blue curve becomes closer to the origin. This means that waves that are not well-resolved are quickly damped out of the system, which may explain the strong performance of FB-RK(3,2) with the corresponding choices of FB weights. Contrast this with the behavior of SSPRK3 in Fig. 1f; as we move toward the grid scale, rather than damping out the potentially problematic motions, the eigenvalues simply leave the unit circle, signaling that the method has become unstable.

Note that the paradigm in which the plots of Fig. 1 exist is different than the paradigm for which the FB weights were optimized. That is, the matrix $\mathsf{G}$ 0 comes from applying FB-RK(3,2) to a one-dimensional gravity wave system (22), whereas the FB weights presented in Table 1 are optimized for the linearized SWEs (10). We present the plots in Fig. 1 to show that the schemes optimized for stability in the linearized SWEs treat well-resolved gravity waves correctly in the best case, where the equations are not discretized in space. This also allows for a quick visual comparison between other methods presented with similar plots in the literature; in particular with the methods presented in ShMc.

Here, we present a number of numerical experiments that demonstrate the efficiency and accuracy of the CFL-optimized FB-RK(3,2) schemes presented in section 2. These experiments compare SSPRK3 and RK(3,2) to FB-RK(3,2) in terms of time-stepping efficiency, temporal convergence, and solution accuracy. We consider a total of five different nonlinear shallow-water test cases described below. The fully nonlinear SWEs are given by (A14) in section b of the appendix.

• Quasi-linear gravity wave (QLW): A simple gravity wave traveling from the north pole to the south, then returning to the north, a process which takes approximately seven days of simulated time. This is the case of an external wave mode, also called the Lamb mode in the context of atmospheric modeling. The momentum advection terms in the momentum equation of (A14) are turned off, making the momentum equation purely linear. The thickness equation (A14) is formally nonlinear, but because the fluid depth is sufficiently large compared to any perturbation in the thickness, this behaves as though it is linear. Therefore, we refer to this test case as quasi-linear. The fluid velocity is initialized to zero, while the thickness is initialized as a Gaussian bell curve centered at the north pole of the form $\exp [-100{(x-\pi )}^2–100y^2]+\hspace{0.17em}500$ .
• Barotropically unstable jet (BUJ): Originally formulated by [12], this case consists of a strong zonal flow at midlatitude with a small perturbation that causes the barotropic jet to become unstable. This flow is driven by strongly nonlinear momentum advection tendencies, physical processes that are of particular relevance to weather prediction and climate modeling.
• Williamson test case 2 (WTC2): A steady state solution to the nonlinear SWEs, one of the standard SWE test cases developed by [37]. The case consists of a geostrophically balanced zonal flow that is constant for all time.
• Williamson test case 4 (WTC4): A forced nonlinear flow consisting of a translating low pressure center superimposed on a jet stream, one of the standard [37] SWE test cases.
• Williamson test case 5 (WTC5): A zonal flow over an isolated mountain, one of the standard [37] SWE test cases. Here, the initial flow is exactly as in WTC2, but the seafloor topography is nonuniform, consisting of an isolated mountain at midlatitude. This perturbation alters the geostrophic balance of the flow, triggering the development of advection-driven instabilities.

All test cases take place on a rotating aquaplanet with angular velocity Ω = 7.292 × 10−5 s−1 and a radius of 6371.22 km. The seafloor topography is uniform unless otherwise specified. Each test case is implemented using a TRiSK spatial discretization ([28]), which is a finite volume-type spatial discretization made for unstructured, variable-resolution polygonal grids. Here, we slightly alter the standard TRiSK discretization by using an upwind potential vorticity flux to remove gridscale noise in the evolution of vorticity. Icosahedral meshes were used throughout, with quasi-uniform 60-km resolution employed in all test cases. Meshes were optimized according to the spherical centroidal Voronoi tessellation (SCVT) metric using the JIGSAW library ([9]), ensuring the TRiSK discretization achieves approximately second-order accuracy in the L2 norm. Balanced initial conditions for the BUJ and WTC test cases are generated following the procedure described in section e of the appendix. The code used to run all the experiments described in this section is open source and available on GitHub and Zenodo; see the data availability statement at the end of this document.

On all of the test cases described above, the optimized FB-RK(3,2) schemes from section 2d admit time steps between 1.59 and 2.81 times as large as SSPRK3 or RK(3,2). This demonstrates the ability for the FB averaging employed by FB-RK(3,2) and the subsequent optimization of the FB weights to greatly increase the stability range of an existing method across multiple problems. SSPRK3 is used as a baseline comparison here since it admits the smallest time steps for each case out of the schemes considered here, and it also consists of three Runge–Kutta stages, costing the same number of floating-point operations as RK(3,2) and FB-RK(3,2).

The largest increases in the admittable time step are in the QLW test case, with the largest increase at 2.81 times in the best case (Table 2). For practical reasons, we used a nonlinear SWE model to run the above test cases, but QLW was intentionally designed so that the model equations, while still formally nonlinear, approximate the behavior of the linear SWEs. That is, QLW provides a test case as close as possible to the paradigm for which FB-RK(3,2) was optimized. As noted, the comparison is not one-to-one. FB-RK(3,2) was optimized for the linearized SWEs on a square C-grid, while QLW uses quasi-linear SWEs and a hexagonal grid with a TRiSK discretization. Nonetheless, the improvement in CFL performance is remarkable and clearly demonstrates the efficacy of the optimization performed in section 2d.

Graph: Table 2. CFL performance of FB-RK(3,2) with optimal FB weights from Table 1 compared to SSPRK3 and RK(3,2) on various test cases. The recorded values for Δ t represent the largest time step for which the given test case could run without becoming unstable for the given duration. Time steps were determined experimentally by running the model for the given duration, increasing the time step until the model became unstable. For practical reasons, time steps were found as multiples of 5.

The remaining four test cases are fully nonlinear, covering a range of standard shallow-water test cases relevant to atmospheric and oceanic flows. In these, FB-RK(3,2) admits time steps between 1.59 and 2.18 times larger than SSPRK3. As SSPRK3 consists of the same number of Runge–Kutta stages as FB-RK(3,2), this directly translates to an increase of overall efficiency of between 1.59 and 2.18. The increases in admittable time step in these nonlinear cases are slightly less than those in the QLW case in part because these problems include fully nonlinear motions that cannot be accounted for in the optimization process used here; because our optimization relies on properties of von Nuemann stability analysis, we are limited to optimizing on a linear set of equations. However, Table 1 suggests that optimizing for the linearized SWEs provides substantial benefits even for the nonlinear SWEs. Indeed, as noted above, our analysis and the subsequent optimization was performed under a number of simplifications versus the paradigm in which these numerical experiments were performed. Nonetheless, our results show that the optimized FB-RK(3,2) schemes are robust when applied to fully nonlinear problems with a full TRiSK discretization. This provides evidence that our schemes are at least partially agnostic to the details of the spatial discretization.

Additionally, one can observe that the optimized FB-RK(3,2) schemes are sensitive to the Froude number; as this increases, the corresponding values of ${\nu }_{({\beta }_1,{\beta }_2,{\beta }_3)}^{\text{max}}$ decrease (Table 1). As the Froude number increases, the flow becomes dominated by advection rather than gravity wave propagation. When the stability is dominated by fast waves, FB-RK(3,2) works best as evidenced by its performance in the QLW case. As the wave and advective phase speeds come closer together, the speed-up from FB is reduced, due to the flow taking on more advective nature as can be seen in the four fully nonlinear test cases.

Turning our attention to the four nonlinear test cases, we see that FB-RK(3,2) with the FB weights obtained by optimizing C1 with |U| = 0.05 admits the largest time step in each case except WTC2 where it instead admits the second largest time step. This points to (β1, β2, β3) = (0.531, 0.531, 0.313) as a robust option for a broad span of problems. However, it is almost certain that the optimal choice of coefficient will be problem dependent.

In section a of the appendix, we perform a simplified local truncation error analysis that shows that FB-RK(3,2) is second-order in time for any choice of FB weights on the one dimensional linearized SWEs. Here, we support the assertion that FB-RK(3,2) is second-order in time with numerical experiments.

Each tested variation of FB-RK(3,2) is second-order in time (Fig. 2). Differences in accuracy between the variations of FB-RK(3,2) are negligible, showing that the choice of optimized FB weights has little effect on accuracy in this case. RK(3,2), the scheme on which FB-RK(3,2) is based, is third-order as expected as this is a linear problem ([27]).

Graph: Fig. 2. Temporal convergence for FB-RK(3,2) for a selection of FB weights. Convergence tests were performed using the QLW case, run for 7 simulated days. The L 2 errors are computed against a reference solution generated by the classical four-stage, fourth-order Runge–Kutta method (RK4) with a time step of 10 s.

We evaluate the solution quality of an optimal FB-RK(3,2) to that of SSPRK3 in two test cases of interest, BUJ and WTC5, comparing the relative vorticity produced by FB-RK(3,2) with (β1, β2, β3) = (0.531, 0.531, 0.313) to that produced by SSPRK3. The strong stability preserving property of SSPRK3 makes it a good choice for a point of comparison; if FB-RK(3,2) can produce a solution of comparable quality, we can be confident that it is a good choice for the SWEs and related problems.

In both cases, shown in Figs. 3 and 4, FB-RK(3,2) produces solutions of comparable quality to SSPRK3 even with a significantly increased time step. In BUJ, the two schemes produce relative velocities that differ by at most 10−7. In WTC5, this difference is on the order of 10−8. We also ran both test cases using RK(3,2) which produced plots, which we omit, that are visually identical to those given here.

Graph: Fig. 3. Relative vorticity produced by (a) SSPRK3 and (b) FB-RK(3,2) with (β 1 , β 2 , β 3) = (0.531, 0.531, 0.313) on the BUJ test case after 6 days of simulated time. SSPRK3 uses a time step of 108 s and FB-RK(3,2) uses a time step of 192 s. (c) The absolute difference between the two solutions. The plots themselves are lat/lon projections of the sphere.

Graph: Fig. 4. Relative vorticity produced by (a) SSPRK3 and (b) FB-RK(3,2) with (β 1 , β 2 , β 3) = (0.531, 0.531, 0.313) on the WTC5 test case after 50 days of simulated time. SSPRK3 uses a time step of 135 s and FB-RK(3,2) uses a time step of 288 s. (c) The absolute difference between the two solutions. The plots themselves are lat/lon projections of the sphere.

FB-RK(3,2) uses a time step 1.78 and 2.13 times larger that that used by SSPRK3 on BUJ and WTC5, respectively. Despite taking about half the wall-clock time and computational resources, FB-RK(3,2) produces vorticities that are qualitatively equivalent to SSPRK3. As noted in section 3a, though FB-RK(3,2) is lower order in time than SSPRK3, the produced solutions are qualitatively equivalent; likely partially due to lower-order spatial errors which dominate. In operational configurations of climate-scale models, such spatial errors are generally dominant, so one can use a more economical time-stepping scheme like FB-RK(3,2) without sacrificing model quality.

With the goal of producing computationally efficient time-stepping schemes for the SWEs, we have presented a class of three-stage second-order Runge–Kutta schemes that use forward–backward averaging to advance the momentum equation. The FB weights in these FB-RK(3,2) schemes were then optimized using a von Neumann–type analysis on the shallow-water equations linearized about a nonzero mean flow. Even though it is impossible to perform a similar analysis and optimization on the fully nonlinear SWEs, the results of our linearized analysis clearly result in significant increases in efficiency when solving nonlinear test problems.

Our optimized FB-RK(3,2) schemes compare favorably in terms of effective CFL number to an existing forward–backward scheme in [31]. In numerical experiments, the schemes are up to 2.8 times more computationally efficient than the RK(3,2) scheme from which they were derived and a popular three-stage third-order strong stability preserving Runge–Kutta scheme (Table 2). Despite being significantly more computationally economical than other schemes, FB-RK(3,2) produces solutions that are of comparable quality to other schemes while using a time step approximately twice as long (Figs. 3 and 4). FB-RK(3,2) schemes are generally applicable to geophysical fluid flows in the ocean and atmosphere, with a range of Froude numbers, as evidenced by their performance across a number of relevant test cases. Their use provides an economical pathway to increasing the computational efficiency of models of any scope, but particularly models at the climate scale wherein time-to-solution can be restrictive.

Moving forward, we are interested in adapting FB-RK(3,2) to a local time-stepping (LTS) scheme similar to that used in [7] or [15], thereby combining the benefits of the CFL performance of FB-RK(3,2) with the obvious performance implications of LTS schemes. Another potentially valuable application of FB-RK(3,2) would be its use as a barotropic solver within a split-explicit, layered model wherein barotropic and baroclinic motions are separated and solved by different methods. The barotropic equations are essentially the SWEs and our work here clearly demonstrates that FB-RK(3,2) is a good choice for solving such a system. As the number of MPI processes increases, the communication costs with in the barotropic solver begin to dominate within a split ocean model. Taking a longer barotropic time step would help to assuage some of these costs.

JRL was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by ORAU under Contract DE-SC0014664. JRL was also supported by the National Science Foundation (NSF) Mathematical Sciences Graduate Internship (MSGI). GC, MRP, and DE were supported by the Earth System Model Development program area of the U.S. DOE, Office of Science, Office of Biological and Environmental Research as part of the multi-program, collaborative Integrated Coastal Modeling (ICoM) project. The authors would also like to thank Pedro S. Peixoto and two anonymous reviewers for their insightful comments and suggestions.

What follows is a simplified local truncation error analysis of FB-RK(3,2) that shows that the method is at least $\mathcal{O}[{(\mathrm{\Delta }t)}^2]$ on a particular problem for any choice of the FB weights βi. To ease notation and simplify symbolic calculations, consider a simplified version of the linearized SWEs (5) with zero mean flow, in one spatial dimension, and where f = 0:

(A1) $\begin{array}{l}\frac{\partial u}{\partial t}=-c\frac{\partial \eta }{\partial x},\\ \frac{\partial \eta }{\partial t}=-c\frac{\partial u}{\partial x}.\end{array}$

Assume that u is sufficiently smooth that derivatives up to sixth order in both t and x exist, and that η is sufficiently smooth that derivatives up to fifth order in both t and x exist. Further, assume that both u and η are sufficiently smooth that $(\partial /\partial t)({\partial }^k/\partial x^k)u=({\partial }^k/\partial x^k)(\partial /\partial t)u$ and $(\partial /\partial t)({\partial }^k/\partial x^k)\eta =({\partial }^k/\partial x^k)(\partial /\partial t)\eta$ for k = 1, 2. This second assumption gives us

(A2) $\frac{{\partial }^{k}u}{\partial t^k}=\{\begin{array}{ll}-c^k\frac{{\partial }^k\eta }{\partial x^k},\hfill & \text{if}\hspace{0.17em}k\hspace{0.17em}\text{is}\hspace{0.17em}\text{odd}\hfill \\ c^k\frac{{\partial }^{k}u}{\partial x^k},\hfill & \text{if}\hspace{0.17em}k\hspace{0.17em}\text{is}\hspace{0.17em}\text{even}\hfill \end{array},$

(A3) $\frac{{\partial }^k\eta }{\partial t^k}=\{\begin{array}{ll}-c^k\frac{{\partial }^{k}u}{\partial x^k},\hfill & \text{if}\hspace{0.17em}k\hspace{0.17em}\text{is}\hspace{0.17em}\text{odd}\hfill \\ c^k\frac{{\partial }^k\eta }{\partial x^k},\hfill & \text{if}\hspace{0.17em}k\hspace{0.17em}\text{is}\hspace{0.17em}\text{even}\hfill \end{array},$

for k = 1, 2. Also, note that in the analysis that follows, the dependence of u and η on x is suppressed to simplify notation.

Now, by applying FB-RK(3,2) to (A1), we can write

(A4) $\begin{array}{l}\frac{u^{n+1}-u^n}{\mathrm{\Delta }t}=-c\frac{\partial {\eta }^n}{\partial x}+\frac{\mathrm{\Delta }t}{2}{c}^2\frac{{\partial }^2{u}^n}{\partial x^2}-\frac{{(\mathrm{\Delta }t)}^2}{6}{c}^3({\beta }_3+1)\frac{{\partial }^3{\eta }^n}{\partial x^3}+\frac{{(\mathrm{\Delta }t)}^3}{36}{c}^4(-4{\beta }_1{\beta }_3+2{\beta }_1+9{\beta }_2{\beta }_3)\frac{{\partial }^4{u}^n}{\partial x^4}\\ -\frac{{(\mathrm{\Delta }t)}^4}{12}{c}^5{\beta }_2{\beta }_3\frac{{\partial }^5{\eta }^n}{\partial x^5}+\frac{{(\mathrm{\Delta }t)}^5}{36}{c}^6{\beta }_1{\beta }_2{\beta }_3\frac{{\partial }^6{u}^n}{\partial x^6}\\ =U(u^n,{\eta }^n),\end{array}$

(A5) $\begin{array}{l}\frac{{\eta }^{n+1}-{\eta }^n}{\mathrm{\Delta }t}=-c\frac{\partial u^n}{\partial x}+\frac{\mathrm{\Delta }t}{2}{c}^2\frac{{\partial }^2{\eta }^n}{\partial x^2}-\frac{{(\mathrm{\Delta }t)}^2}{4}{c}^3{\beta }_2\frac{{\partial }^3{u}^n}{\partial x^3}+\frac{{(\mathrm{\Delta }t)}^3}{12}{c}^4{\beta }_2\frac{{\partial }^4{\eta }^n}{\partial x^4}-\frac{{(\mathrm{\Delta }t)}^4}{36}{c}^5{\beta }_1{\beta }_2\frac{{\partial }^5{u}^n}{\partial x^5}\\ =H(u^n,{\eta }^n).\end{array}$

Notice that the right-hand side of (A4) and (A5) are equations in un and ηn; call them U(un, ηn) and H(un, ηn), respectively.

The local truncation error for a given time-stepping scheme can be thought of as the error introduced by said scheme during one time step, assuming exact data at the beginning of the time step. For this problem, the local truncation errors in u and η are defined to be the quantities ${\tau }_u^n$ and ${\tau }_{\eta }^n$ such that

(A6) $\frac{u(t^{n+1})-u(t^n)}{\mathrm{\Delta }t}=U[u(t^n),\eta (t^n)]+{\tau }_u^n,$

(A7) $\frac{\eta (t^{n+1})-\eta (t^n)}{\mathrm{\Delta }t}=H[u(t^n),\eta (t^n)]+{\tau }_{\eta }^n.$

To calculate ${\tau }_u^n$ and ${\tau }_{\eta }^n$ , rearrange (A6) and (A7), then expand u(tn+1) and η(tn+1) in t as Taylor polynomials centered at t = tn to get

(A8) $\begin{array}{l}\mathrm{\Delta }t{\tau }_u^n=u(t^{n+1})-u(t^n)-\mathrm{\Delta }tU[u(t^n),\eta (t^n)]\\ =\left\{{\displaystyle \sum _{k=0}^6\frac{{(\mathrm{\Delta }t)}^k}{k!}}\frac{{\partial }^{k}u}{\partial t^k}(t^n)+\mathcal{O}[{(\mathrm{\Delta }t)}^7]\right\}-u(t^n)-\mathrm{\Delta }tU[u(t^n),\eta (t^n)],\end{array}$

(A9) $\begin{array}{l}\mathrm{\Delta }t{\tau }_{\eta }^n=\eta (t^{n+1})-\eta (t^n)-\mathrm{\Delta }tH[u(t^n),\eta (t^n)]\\ =\left\{{\displaystyle \sum _{k=0}^5\frac{{(\mathrm{\Delta }t)}^k}{k!}}\frac{{\partial }^k\eta }{\partial t^k}(t^n)+\mathcal{O}[{(\mathrm{\Delta }t)}^6]\right\}-\eta (t^n)-\mathrm{\Delta }tH[u(t^n),\eta (t^n)].\end{array}$

Now, simplify (A8) and (A9). Starting with (A8), we get

(A10) $\begin{array}{c}\Delta t{\tau }_u^n=\left\{{\displaystyle \sum _{k=0}^3\frac{{\left(\Delta t\right)}^k}{k!}}\frac{{\partial }^{k}u}{\partial t^k}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}-u\left(t^n\right)-\Delta tU\left[u\left(t^n\right),\eta \left(t^n\right)\right]\\ =\left\{u\left(t^n\right)+\Delta t\underset{=-c\frac{\partial \eta }{\partial x}\left(t^n\right)}{\underbrace{\frac{\partial u}{\partial t}\left(t^n\right)}}+\frac{{\left(\Delta t\right)}^2}{2}\underset{=c^2\frac{{\partial }^{2}u}{\partial x^2}\left(t^n\right)}{\underbrace{\frac{{\partial }^{2}u}{\partial t^2}\left(t^n\right)}}+\frac{{\left(\Delta t\right)}^3}{6}\frac{{\partial }^{3}u}{\partial t^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}\\ +\left\{-u\left(t^n\right)+\Delta tc\frac{\partial {\eta }^n}{\partial x}\left(t^n\right)-\frac{{\left(\Delta t\right)}^2}{2}{c}^2\frac{{\partial }^2{u}^n}{\partial x^2}\left(t^n\right)+\frac{{\left(\Delta t\right)}^3}{6}{c}^3\left({\beta }_3+1\right)\frac{{\partial }^3{\eta }^n}{\partial x^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}\\ =\left\{\frac{{\left(\Delta t\right)}^3}{6}\frac{{\partial }^{3}u}{\partial t^3}\left(t^n\right)+\frac{{\left(\Delta t\right)}^3}{6}{c}^3\left({\beta }_3+1\right)\frac{{\partial }^3\eta }{\partial x^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}.\end{array}$

From (A10), we see that $\mathrm{\Delta }t{\tau }_u^n=\mathcal{O}[{(\mathrm{\Delta }t)}^3]$ no matter the values of the FB weights, so

(A11) ${\tau }_u^n=\mathcal{O}[{(\mathrm{\Delta }t)}^2].$

Next simplify (A9) in a similar way:

(A12) $\begin{array}{c}\Delta t{\tau }_{\eta }^n=\left\{{\displaystyle \sum _{k=0}^3\frac{{\left(\Delta t\right)}^k}{k!}}\frac{{\partial }^k\eta }{\partial t^k}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}-\eta \left(t^n\right)-\Delta tH\left[u\left(t^n\right),\eta \left(t^n\right)\right]\\ =\left\{\eta \left(t^n\right)+\Delta t\underset{=-c\frac{\partial u}{\partial x}\left(t^n\right)}{\underbrace{\frac{\partial \eta }{\partial t}\left(t^n\right)}}+\frac{{\left(\Delta t\right)}^2}{2}\underset{=c^2\frac{{\partial }^2\eta }{\partial x^2}\left(t^n\right)}{\underbrace{\frac{{\partial }^2\eta }{\partial t^2}\left(t^n\right)}}+\frac{{\left(\Delta t\right)}^3}{6}\frac{{\partial }^3\eta }{\partial t^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}\\ +\left\{-\eta \left(t^n\right)+\Delta tc\frac{\partial u^n}{\partial x}\left(t^n\right)-\frac{{\left(\Delta t\right)}^2}{2}{c}^2\frac{{\partial }^2{\eta }^n}{\partial x^2}\left(t^n\right)+\frac{{\left(\Delta t\right)}^3}{4}{c}^3{\beta }_2\frac{{\partial }^3{u}^n}{\partial x^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}\\ =\left\{\frac{{\left(\Delta t\right)}^3}{6}\frac{{\partial }^3\eta }{\partial t^3}\left(t^n\right)+\frac{{\left(\Delta t\right)}^3}{4}{c}^3{\beta }_2\frac{{\partial }^{3}u}{\partial x^3}\left(t^n\right)+\mathcal{O}\left[{\left(\Delta t\right)}^4\right]\right\}.\end{array}$

From (A12), we see that with no assumptions on the FB weights, we have $\mathrm{\Delta }t{\tau }_{\eta }^n=\mathcal{O}[{(\mathrm{\Delta }t)}^3]$ , so

(A13) ${\tau }_{\eta }^n=\mathcal{O}[{(\mathrm{\Delta }t)}^2].$

Here, we present the derivation of (4), the linearized SWEs centered about a nonzero mean flow. The nonlinear SWEs on a rotating sphere are given by

(A14) $\begin{array}{l}\frac{\partial \stackrel{⌣}{u}}{\partial t}+(\nabla \times \stackrel{⌣}{u}+f\mathbf{k})\times \stackrel{⌣}{u}=-\nabla \frac{{\left|\stackrel{⌣}{u}\right|}^2}{2}-g\nabla \stackrel{⌣}{h},\\ \frac{\partial \stackrel{⌣}{h}}{\partial t}+\nabla \cdot (\stackrel{⌣}{h}\stackrel{⌣}{u})=0,\end{array}$

where $\stackrel{⌣}{u}(x,y,t)=[\stackrel{⌣}{u}(x,y,t),\stackrel{⌣}{\upsilon }(x,y,t)]$ is the horizontal fluid velocity, x and y are the spatial coordinates, t is the time coordinate, f is the Coriolis parameter, g is the gravitational constant, and $\stackrel{⌣}{h}$ is the fluid thickness.

To linearize (A14) about a nonzero mean flow, assume that

(A15) $\stackrel{⌣}{u}(x,y,t)=\tilde{U}+\tilde{u}(x,y,t),$

(A16) $\stackrel{⌣}{h}(x,y,t)=H+\tilde{h}(x,y,t),$

where $\tilde{U}=(\tilde{U},\tilde{V})$ is the constant (with respect to time and space) mean horizontal fluid velocity, $\tilde{u}=[\tilde{u}(x,y,t),\tilde{\upsilon }(x,y,t)]$ is a small perturbation of $\tilde{U}$ , H is the constant (with respect to time and space) fluid thickness when at rest, and $\tilde{h}=\tilde{h}(x,y,t)$ is a small perturbation of H. In both cases, small means that the quantities themselves, as well as their derivatives, are negligible when multiplied together or squared.

Now, we make the substitutions $\stackrel{⌣}{u}=\tilde{U}+\tilde{u}$ and $\stackrel{⌣}{h}=H+\tilde{h}$ and simplify term by term. The kinetic energy term becomes

(A17) $\begin{array}{l}-\nabla \frac{{\stackrel{⌣}{u}}^2}{2}=-\nabla \frac{1}{2}({\stackrel{⌣}{u}}^2+{\stackrel{⌣}{\upsilon }}^2)\\ =-\nabla \frac{1}{2}({\tilde{U}}^2+2\tilde{U}\tilde{u}+{\tilde{u}}^2+{\tilde{V}}^2+2\tilde{V}\tilde{\upsilon }+{\tilde{\upsilon }}^2)\\ =-\tilde{U}\nabla \tilde{u}-\tilde{V}\nabla \tilde{\upsilon }-\frac{1}{2}(\nabla {\tilde{u}}^2+\nabla {\tilde{\upsilon }}^2)\\ \approx -\tilde{U}\nabla \tilde{u}-\tilde{V}\nabla \tilde{\upsilon }\\ =\tilde{U}\cdot \nabla \tilde{\mathbf{u}},\end{array}$

where $\nabla \tilde{u}$ is the Jacobian matrix,

$\left(\begin{array}{cc}{\partial }_x\tilde{u}& {\partial }_x\tilde{\upsilon }\\ {\partial }_y\tilde{u}& {\partial }_y\tilde{\upsilon }\end{array}\right),$

and $\tilde{U}\cdot \mathbf{\nabla }\tilde{\mathbf{u}}$ denotes a 2 × 1 column vector wherein the ith entry is the dot product of $\tilde{U}$ with the ith row of $\nabla \tilde{u}$ . The potential vorticity and Coriolis term becomes

(A18) $\begin{array}{l}\left(\nabla \times \stackrel{⌣}{u}+fk\right)\times \stackrel{⌣}{u}=\left[\left(\frac{\partial \stackrel{⌣}{\upsilon }}{\partial x}-\frac{\partial \stackrel{⌣}{u}}{\partial y}\right)k+fk\right]\times \stackrel{⌣}{u}\\ =i\left[-\stackrel{⌣}{\upsilon }\left(\frac{\partial \stackrel{⌣}{\upsilon }}{\partial x}-\frac{\partial \stackrel{⌣}{u}}{\partial y}+f\right)\right]+j\left[\stackrel{⌣}{u}\left(\frac{\partial \stackrel{⌣}{\upsilon }}{\partial x}-\frac{\partial \stackrel{⌣}{u}}{\partial y}+f\right)\right]\\ =i\left[-f\left(\tilde{V}+\tilde{\upsilon }\right)-\tilde{V}\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)-\tilde{\upsilon }\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)\right]\\ +j\left[f\left(\tilde{U}+\tilde{u}\right)+\tilde{U}\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)+\tilde{u}\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)\right]\\ \approx i\left[-f\left(\tilde{V}+\tilde{\upsilon }\right)-\tilde{V}\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)\right]+j\left[f\left(\tilde{U}+\tilde{u}\right)+\tilde{U}\left(\frac{\partial \tilde{\upsilon }}{\partial x}-\frac{\partial \tilde{u}}{\partial y}\right)\right]\\ =f{\left(\tilde{U}+\tilde{u}\right)}^{\perp }+\tilde{\zeta }{\tilde{U}}^{\perp },\end{array}$

where $\tilde{\zeta }=(\partial \tilde{\upsilon }/\partial x)-(\partial \tilde{u}/\partial y)$ is the vertical component of the vorticity $\nabla \times \tilde{u}$ . The thickness gradient becomes