The Matsuno baroclinic wave test case

The analytic wave solutions obtained by Matsuno (1966) in his seminal work on equatorial waves provide a simple and informative way of assessing the performance of atmospheric models by measuring the accuracy with which they simulate these waves. These solutions approximate the solutions of the shallow-water equations on the sphere for low gravity-wave speeds such as those of the baroclinic modes in the atmosphere. This is in contrast to the solutions of the non-divergent barotropic vorticity equation, used in the Rossby–Haurwitz test case, which are only accurate for high gravity-wave speeds such as those of the barotropic mode. The proposed test case assigns specific values to the wave parameters (gravity-wave speed, zonal wave number, meridional wave mode and wave amplitude) for both planetary and inertia-gravity waves, and suggests simple assessment criteria suitable for zonally propagating wave solutions. The test is successfully applied to a spherical shallow-water model in an equatorial channel and to a global-scale model. By adding a small perturbation to the initial fields it is demonstrated that the chosen initial waves remain stable for at least 100 wave periods. The proposed test case can also be used as a resolution convergence test.

1 Introduction

A cornerstone of global-scale model assessment is the Rossby–Haurwitz test case, originally used by Phillips (1959) as a qualitative way of assessing his shallow-water model. Phillips initialized his model with an analytic wave solution of the non-divergent barotropic vorticity equation obtained by Haurwitz (1940), and examined the spatiotemporal smoothness of the simulated fields at later times. Using this procedure he concluded that the emergent noise in his model was due to a small but significant divergence field missing from the initial fields. Even though the solutions of the non-divergent barotropic vorticity equation are not solutions of the shallow-water equations (SWEs), Phillips' procedure was adopted by Williamson et al. (1992) as a standard test case for shallow-water models and has been extensively used ever since (Jablonowski et al., 2009; Mohammadian and Marshall, 2010; Bosler et al., 2014; Ullrich, 2014; Li et al., 2015, are only five recent examples).

However, there are two known issues with the original Rossby–Haurwitz test case that limit its usefulness (Thuburn and Li, 2000). The first is the generation of small-scale features via a potential enstrophy cascade, which requires adequate dissipation mechanisms to remove enstrophy at the grid scale (in order to mimic a continuous cascade to sub-grid scales). The second is the instability of the initial wave number 4 used in the Rossby–Haurwitz test case. In contrast to Hoskins (1973), who found that wave numbers smaller than or equal to 5 are stable, Thuburn and Li (2000) show that the Rossby–Haurwitz wave number 4 is in fact also unstable.

Recently, Shamir and Paldor (2016) proposed a similar procedure to that of Phillips (1959) where instead of using the solutions of the non-divergent barotropic vorticity equation, the initial fields are the analytic wave solutions of the linearized SWEs on the sphere derived in Paldor et al. (2013). These solutions fully account for the small divergence field and can be computed on any grid given the latitudes and longitudes. In particular, they include the fast propagating inertia-gravity (IG) waves that are completely absent from the non-divergent barotropic vorticity equation. Consequently, the procedure proposed by Shamir and Paldor (2016) provides a more quantitative assessment than the original procedure from Phillips (1959) while it is just as easy to implement.

Both solutions obtained by Haurwitz (1940) and Paldor et al. (2013) approximate the solutions of the SWEs in the asymptotic limit of high gravity-wave speeds. For most practical purposes they are sufficiently accurate for gravity-wave speeds of about 200–300 m s⁻¹ or higher, which are typical of the barotropic mode in Earth's atmosphere and oceans. However, the typical speeds of gravity waves of baroclinic modes in the (tropical) atmosphere are about 20–30 m s⁻¹ (Wheeler and Kiladis, 1999). Thus, the abovementioned procedures are only relevant for assessing the accuracy with which the barotropic wave mode is simulated. In order to assess the accuracy of the baroclinic wave modes we propose, in the present work, to use the analytic wave solutions of the linearized SWEs on the equatorial β-plane obtained by Matsuno (1966) that approximate the solutions of the SWEs on the sphere in the asymptotic limit of low gravity-wave speeds. (De-Leon and Paldor, 2011; Garfinkel et al., 2017).

In addition to being on two opposite ends of the spectrum of gravity-wave speed the solutions obtained by Matsuno (1966) differ from those obtained by both Haurwitz (1940) and Paldor et al. (2013) in their meridional extent. While the former become negligibly small outside a narrow equatorial band, the latter two have non-negligible amplitudes in the vicinity of the poles. Thus, while the Rossby–Haurwitz test case is only relevant to global-scale models, the test case proposed in the present study is applicable to both global-scale and tropical models.

A homonymous, but unrelated, test case is the baroclinic wave test case developed in Jablonowski (2004) and Jablonowski and Williamson (2006) and independently in Polvani et al. (2004), and its variants in Lauritzen et al. (2010) and Ullrich et al. (2014). This test case is concerned with the non-linear generation of synoptic-scale eddies in multilayer models via baroclinic instability. In contrast, the test case proposed here is concerned with linear wave propagation in (non-linear) single-layer models. In particular, while the term baroclinic usually implies the use of multilayer models, here this term is used to denote a single thin layer model of homogeneous density where the gravity-wave speeds are similar to those observed in baroclinic modes in the atmosphere.

The idea of using Matsuno's solutions as a test case in a similar fashion to that of the Rossby–Haurwitz test case is most likely not original, but has never been standardized. Thus, the purpose of the present work is to standardize the Matsuno test case similar to the way that Williamson et al. (1992) standardized the Rossby–Haurwitz test case. We start with a short description of the analytic expressions derived by Matsuno (1966) in Sect. 2. The proposed test procedure, including the choice of wave parameters and assessment criteria, is described in Sect. 3. In Sect. 4 we demonstrate the usefulness of the proposed test case using both an equatorial channel spherical shallow-water model, and a global-scale model. In addition, we examine the smoothness and stability of the initial waves in a similar fashion to that used in Thuburn and Li (2000) and demonstrate the possibility of using the proposed test case as a resolution convergence test. Finally, concluding remarks are given in Sect. 5.

2 The analytic solutions

The proposed test case is based on the analytic solutions of the SWEs on the equatorial β-plane obtained by Matsuno (1966). These solutions have the form of zonally propagating waves, i.e.,

$$\begin{array}{}\text{(1)}& {\displaystyle }\left[\begin{array}{c}u(x,y,t)\\ v(x,y,t)\\ \mathrm{\Phi }(x,y,t)\end{array}\right]=\mathit{\text{Re}}\left\{\left[\begin{array}{c}\widehat{u}\left(y\right)\\ \widehat{v}\left(y\right)\\ \widehat{\mathrm{\Phi }}\left(y\right)\end{array}\right]e^{i(kx-\mathit{\omega }t)}\right\},\end{array}$$

where x and y are the local Cartesian coordinates in the zonal and meridional directions, respectively; t is time; u and v are the velocity components in the zonal and meridional directions, respectively; Φ is the geopotential height; k is the planar zonal wave number (which has dimensions of m⁻¹); ω is the wave frequency; $\widehat{u}\left(y\right),\widehat{v}\left(y\right)$ and $\widehat{\mathrm{\Phi }}\left(y\right)$ are the latitude-dependent amplitudes; and $i=\sqrt{-\mathrm{1}}$ is the imaginary unit. In accordance with the sign convention used in Matsuno (1966) we assume that k is non-negative and let ω take any real value. Note, however, that the sign in front of ω in Eq. (1) is the opposite to that in the theory of Matsuno (1966). The convention chosen here is more intuitive as it implies that positive values of ω correspond to waves that propagate in the positive x direction, i.e., eastward.

The unknown wave frequencies and latitude-dependent amplitudes are derived from the (well-known) energies and eigenfunctions of the (time-independent) Schödinger equation of the quantum harmonic oscillator. The resulting frequencies are given by the solutions of the following cubic equation:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{\omega }}_{n,k}^{\mathrm{3}}-\left[gH{k}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{2}\mathrm{\Omega }\sqrt{gH}}{a}}(\mathrm{2}n+\mathrm{1})\right]{\mathit{\omega }}_{n,k}-{\displaystyle \frac{\mathrm{2}\mathrm{\Omega }gHk}{a}}=\mathrm{0},\end{array}$$

for $n=-\mathrm{1},\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$, where $\mathrm{\Omega }=\mathrm{7.29212}\times {\mathrm{10}}^{-\mathrm{5}}$ rad s⁻¹, $a=\mathrm{6.37122}\times {\mathrm{10}}^{\mathrm{6}}$ m and g=9.80616 m s⁻² are the Earth's angular frequency, mean radius and gravitational acceleration, respectively, and H is the mean layer's depth (thickness).

For n≥1 Eq. (2) has three distinct real roots corresponding to a slowly westward propagating Rossby wave, a fast eastward propagating inertia-gravity (EIG) wave and a fast westward propagating inertia-gravity (WIG) wave. For n=0 one of the three roots, the one corresponding to a westward propagating gravity wave with $\mathit{\omega }=-\sqrt{gH}k$, leads to infinite zonal wind and is thus considered a physically unacceptable solution. The remaining two roots correspond to the lowest (i.e., n=0) EIG wave and the mixed Rossby–gravity (MRG) wave. For $n=-\mathrm{1}$ Eq. (2) has one real root $\mathit{\omega }=\sqrt{gH}k$, which correspond to the equatorial Kelvin wave (see Matsuno, 1966). The existence of the latter two waves on a sphere is discussed in Garfinkel et al. (2017) and Paldor et al. (2018).

For given values of the zonal wave number, k, and meridional mode number, n, the roots of the cubic equation can be obtained in a closed analytic form using the solutions of the general cubic equation as follows (e.g., Abramowitz and Stegun, 1964):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\omega }}_{n,k,j}=\mathit{\text{Re}}\left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left({\mathrm{\Delta }}_j+{\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{0}}}{{\mathrm{\Delta }}_j}}\right)\right\},\text{for}j=\mathrm{1},\mathrm{2},\mathrm{3}\end{array}$$

where j stands for the three roots, and where

$$\begin{array}{}\text{(4a)}& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{0}}& {\displaystyle }=\mathrm{3}\left[gH{k}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{2}\mathrm{\Omega }\sqrt{gH}}{a}}(\mathrm{2}n+\mathrm{1})\right],\text{(4b)}& {\displaystyle }{\mathrm{\Delta }}_j& {\displaystyle }=[{\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{4}}+\sqrt{{\mathrm{\Delta }}_{\mathrm{4}}^{\mathrm{2}}-\mathrm{4}{\mathrm{\Delta }}_{\mathrm{0}}^{\mathrm{3}}}}{\mathrm{2}}}{]}^{\mathrm{1}/\mathrm{3}}\exp \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }j}{\mathrm{3}}}i\right),\text{(4c)}& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{4}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{54}\mathrm{\Omega }gHk}{a}}.\end{array}$$

Given the definitions in Eq. (4), the explicit expressions for the frequencies of the Rossby, WIG and EIG waves are obtained by sorting the values in Eq. (3) as follows:

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }\text{Rossby}:{\mathit{\omega }}_{n,k,\mathrm{R}}=-\underset{j=\mathrm{1},\mathrm{2},\mathrm{3}}{min}\left|{\mathit{\omega }}_{n,k,j}\right|,\text{(5b)}& {\displaystyle }& {\displaystyle }\mathrm{\text{Westward inertia gravity}}:{\mathit{\omega }}_{n,k,\mathrm{WIG}}=\underset{j=\mathrm{1},\mathrm{2},\mathrm{3}}{min}{\mathit{\omega }}_{n,k,j},\text{(5c)}& {\displaystyle }& {\displaystyle }\text{Eastward inertia gravity}:{\mathit{\omega }}_{n,k,\mathrm{EIG}}=\underset{j=\mathrm{1},\mathrm{2},\mathrm{3}}{max}{\mathit{\omega }}_{n,k,j}.\end{array}$$

Having found (one of) the wave frequencies for a given combination of n and k, the corresponding latitude-dependent amplitudes can be written as

$$\begin{array}{ll}\text{(6a)}& {\displaystyle }{\widehat{v}}_n& {\displaystyle }=A{\widehat{H}}_n\left({\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}{\displaystyle \frac{y}{a}}\right)\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}{\displaystyle \frac{y}{a}}\right)}^{\mathrm{2}}\right){\displaystyle }{\widehat{u}}_{n,k}& {\displaystyle }={\displaystyle \frac{gH{\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}}{ia\left({\mathit{\omega }}_{n,k}^{\mathrm{2}}-gH{k}^{\mathrm{2}}\right)}}[-\sqrt{{\displaystyle \frac{n+\mathrm{1}}{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\omega }}_{n,k}}{\sqrt{gH}}}+k\right){\widehat{v}}_{n+\mathrm{1}}\\ \text{(6b)}& {\displaystyle }& {\displaystyle }-\sqrt{{\displaystyle \frac{n}{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\omega }}_{n,k}}{\sqrt{gH}}}-k\right){\widehat{v}}_{n-\mathrm{1}}]{\displaystyle }{\widehat{\mathrm{\Phi }}}_{n,k}& {\displaystyle }={\displaystyle \frac{gH{\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}}{ia\left({\mathit{\omega }}_{n,k}^{\mathrm{2}}-gH{k}^{\mathrm{2}}\right)}}[-\sqrt{{\displaystyle \frac{n+\mathrm{1}}{\mathrm{2}}}}\left({\mathit{\omega }}_{n,k}+\sqrt{gH}k\right){\widehat{v}}_{n+\mathrm{1}}\\ \text{(6c)}& {\displaystyle }& {\displaystyle }+\sqrt{{\displaystyle \frac{n}{\mathrm{2}}}}\left({\mathit{\omega }}_{n,k}-\sqrt{gH}k\right){\widehat{v}}_{n-\mathrm{1}}],\end{array}$$

for $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ (the cases $n=-\mathrm{1},\mathrm{0}$ require special treatment), where $\mathit{ϵ}=(\mathrm{2}\mathrm{\Omega }a{)}^{\mathrm{2}}/gH$ is Lamb's parameter, A is an arbitrary amplitude (that has dimensions of m s−1) and ${\widehat{H}}_n$ are the normalized Hermite polynomials of degree n defined by the following three-term recurrence relation (Press et al., 2007):

$$\begin{array}{}\text{(7a)}& {\displaystyle }{\widehat{H}}_{-\mathrm{1}}\left(x\right)& {\displaystyle }=\mathrm{0},\text{(7b)}& {\displaystyle }{\widehat{H}}_{\mathrm{0}}\left(x\right)& {\displaystyle }={\mathit{\pi }}^{-\mathrm{1}/\mathrm{4}},\text{(7c)}& {\displaystyle }{\widehat{H}}_{n+\mathrm{1}}\left(x\right)& {\displaystyle }=x\sqrt{{\displaystyle \frac{\mathrm{2}}{n+\mathrm{1}}}}{\widehat{H}}_n-\sqrt{{\displaystyle \frac{n}{n+\mathrm{1}}}}{\widehat{H}}_{n-\mathrm{1}}.\end{array}$$

It should be noted that the chosen normalization for the latitude-dependent amplitudes in Eq. (6) is different from the one used in Matsuno (1966). We use the above normalization for convenience, as it guarantees that $\widehat{v}$ is independent of both k or ω. Furthermore, the use of the normalized version of the Hermite polynomials also leads to slightly different pre-factors in front of ${\widehat{v}}_{n+\mathrm{1}}$ and ${\widehat{v}}_{n-\mathrm{1}}$ compared with Matsuno (1966). However, they are generally more computationally stable. Finally, the outer parentheses in Eq. (6a) denote the argument of ${\widehat{H}}_n$ and the exponential function, and not multiplicative factors. In other words, the independent variable in this equation is $({\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}y/a)$, and not simply y.

While the solutions obtained by Matsuno (1966) apply for the equatorial β-plane, the proposed test case is intended for use in spherical models. As is shown in Garfinkel et al. (2017), the SWEs on the equatorial β-plane approximate the SWEs on the sphere to zero-order in powers of $\mathrm{1}/{\mathit{ϵ}}^{\mathrm{1}/\mathrm{4}}$. Thus, the solutions obtained by Matsuno are only accurate in the asymptotic limit ϵ→∞. For the fixed values of Earth's angular frequency and mean radius, this implies that the solutions obtained by Matsuno are only accurate for sufficiently low gravity-wave speeds $\sqrt{gH}$.

In practice, in order to use the solutions of Matsuno (1966) in spherical models, the local Cartesian coordinates x and y in the above Eqs. (1) and (6) have to be replaced by the longitude λ and latitude ϕ of the geographical coordinate system. Recall that the transformation from the Cartesian system to the spherical system is $(x,y)\to a(\cos {\mathit{\varphi }}_{\mathrm{0}}\mathit{\lambda },\mathit{\varphi })$, where ϕ₀ is the central latitude at which the planar approximation is applied. Likewise, the planar wave number k in Eqs. (1)–(6) has to be replaced by its spherical counterpart, k_s, using the transformation $k\to k_{\mathrm{s}}/a\cos {\mathit{\varphi }}_{\mathrm{0}}$. Thus, for the equatorial β-plane where ϕ₀=0, the transformation is simply $(x,y)\to a(\mathit{\lambda },\mathit{\varphi })$ and $k\to k_{\mathrm{s}}/a$. In particular, the reader should note that the planar wave number k has units of m⁻¹, while the spherical wave number k_s is dimensionless.

3 The proposed test procedure

The general procedure of the proposed test case is similar to the Rossby–Haurwitz procedure in that the model in question is initialized with velocity and height fields corresponding to a particular wave solution and the time evolution of that wave is then examined. The initial wave fields in this case are taken from the analytic expressions in Sect. 2. The specific choice of wave parameters and assessment criteria in the present work are discussed below, separately. As is often the case, these choices represent compromises between conflicting factors, e.g., adherence to observations vs. adherence to asymptotic validity of the analytic solutions or rigorous testing vs. simplicity. In any case, these choices may be the subject of discourse as deemed appropriate by the community.

3.1 Wave parameters

The wave parameters consist of the speed of gravity waves, $\sqrt{gH}$, the wave number and wave mode, k and n, the wave amplitude, A, and the wave type. Any given combination of these parameters completely specifies a unique wave using the expressions in Eqs. (1)–(6). We consider all other parameters, including the spatiotemporal resolution and the form of diffusion/viscosity terms, to be modeling choices left to the developers. This approach is aimed at testing the models in their intended mode of operation. However, as noted in Polvani et al. (2004), different choices of the form of diffusion/viscosity terms correspond to different sets of equations and may not converge to the same solutions.

The choice of gravity-wave speed $\sqrt{gH}$ is inspired by the observed speed of gravity waves of the baroclinic modes in the atmosphere. In practice, we keep g fixed to Earth's gravitational acceleration and set the speed of gravity waves by letting H=30 m, which is within the range of observed equivalent depths in the equatorial atmosphere (Wheeler and Kiladis, 1999). As mentioned in Sect. 2, the analytic solutions obtained by Matsuno (1966) on the equatorial β-plane are only accurate approximations of the SWEs on the sphere in the asymptotic limit of low gravity-wave speeds. The above value was found to be sufficiently accurate, using trial and error, in the sense that it yields stable solutions for at least 100 wave periods in the simulations described in Sect. 4.

In addition to the speed of gravity waves, the accuracy of Matsuno's solutions also depends on the wave number and wave mode. For a given value of $\sqrt{gH}$, these solutions become asymptotically accurate in the limits k_s, n→0 (but k_s≠0) (De-Leon and Paldor, 2011). Also, the spatial variability and the required spatial resolution both increase with the wave number or wave mode, so both of these considerations suggest that reasonable choices for the wave number and wave mode consist of small to moderate values. The proposed wave number and wave mode are k_s=5 (i.e., $k=\mathrm{5}/a$ m⁻¹) and n=1, i.e., within the range of dominant values observed in the equatorial atmosphere (Wheeler and Kiladis, 1999), but other choices may work just as well provided k_s and n are not too large.

The proposed test case is based on the solutions of the linear SWEs but is intended to be used in non-linear models. Therefore, the wave-amplitude should be sufficiently small so as to satisfy the linearization condition. The proposed amplitude of $\widehat{v}$ in Eq. (6) is $A={\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, chosen by trial and error so as to enable stable solutions for at least 100 wave periods in the simulations of Sect. 4.

In general, there are two qualitatively different wave types, Rossby and IG, which differ in the magnitude of their divergence and vorticity fields. The former is more solenoidal (non-divergent), whereas the latter is more irrotational. In order to assess the performance of the models in these two qualitatively different limits we suggest using one of each. As Rossby waves are exclusively westward propagating, we choose the EIG wave of the two IG waves as the second case to cover the two directions of longitudinal propagation.

For these chosen values of $\sqrt{gH}$, k and n, the wave periods ($T=\frac{\mathrm{2}\mathit{\pi }}{\mathit{\omega }}$) are T=18.5 d for the Rossby wave and T=1.9 d for the EIG wave.

3.2 Assessment criteria

For sufficiently small wave amplitudes we expect the spatiotemporal structure of the simulated solutions to be that of zonally propagating waves, i.e., $q=\widehat{q}\left(\mathit{\varphi }\right)e^{i(k\mathit{\lambda }-\mathit{\omega }t)}$ (where q stands for any of the dependent variables u, v or Φ), with frequency and latitude-dependent amplitudes corresponding to the initial wave. In this case, it is desirable to assess the accuracy of the zonal and meridional structures of the waves independently. A fast and simple way of doing so is using Hovmöller diagrams, where the temporal change in any direction is isolated by intersecting the fields along a fixed value of the other direction. This results in the following two diagrams:

• i.

  A time–longitude diagram is obtained by intersecting the fields at a certain latitude. The contour lines in the time–longitude plane are the set of points satisfying $k\mathit{\lambda }-\mathit{\omega }t=$ constant (for some real constant). Thus, the expected pattern for this diagram is that of straight lines with slopes that equal the inverse of the wave's phase speed k∕ω. In order to avoid small fluctuations in the vicinity of latitudinal zero crossings, we recommend using latitudinal intersects at or near local extrema.

• ii.

  A latitude–time diagram is obtained by intersecting the fields at a certain longitude. For any two wave fronts with an equal phase $k({\mathit{\lambda }}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}})=\mathit{\omega }(t_{\mathrm{2}}-t_{\mathrm{1}})$. Thus, holding λ fixed while varying t from t₁ to t₂ is equivalent to holding t fixed and varying λ from λ₁ to ${\mathit{\lambda }}_{\mathrm{2}}={\mathit{\lambda }}_{\mathrm{1}}+\mathit{\omega }/k(t_{\mathrm{2}}-t_{\mathrm{1}})$. The resulting pattern is similar to that of a latitude–longitude diagram, but provides an estimate of the time evolution at a particular longitude (as opposed to a latitude–longitude snapshot at a particular time).

Likewise, for zonally propagating waves it is also desirable to isolate the errors in the phase speed and spatial structure. As discussed in Shamir and Paldor (2016), the frequently used spherical l₂ error entangles the two, and is therefore of lesser use for assessing the accuracy with which the model simulates a propagating wave. Thus, for a more quantitative assessment we suggest using the relative difference between the root-mean-square of the analytic solution and the simulated solutions, i.e.,

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\sqrt{I\left[q^{\mathrm{2}}\right]}-\sqrt{I\left[q_{\mathrm{a}}^{\mathrm{2}}\right]}}{\sqrt{I\left[q_{\mathrm{a}}^{\mathrm{2}}\right]}}},\end{array}$$

where the quantities q and q_a (which are scalars in the case of geopotential and vectors in the case of velocity) correspond to the simulated and analytic solutions, respectively, and where

$$\begin{array}{}\text{(9)}& {\displaystyle }I\left[q\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}q(\mathit{\lambda },\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }.\end{array}$$

Henceforth we refer to the quantity in Eq. (8) as the “structure error” as, in contrast to the l₂ error, it is unaffected by phase speed errors (i.e., phase shifts in λ).

4 Results

In this section we demonstrate the usefulness of the Matsuno test case by applying the proposed procedure to both an equatorial channel finite-difference model and a global-scale spectral model. We then examine the stability of the selected waves/modes in a similar fashion to that used in Thuburn and Li (2000) for the wave number 4 Rossby–Haurwitz wave. Finally, we demonstrate the possibility of using the analytic solutions obtained by Matsuno as a resolution convergence test.

4.1 Demonstration using an equatorial channel finite-difference model

The model is a spherical version of the Cartesian model used in Gildor et al. (2016), in which the integration forward in time is carried out using the conservation form of the SWEs

$$\begin{array}{ll}{\displaystyle \frac{\partial U}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}\left({\displaystyle \frac{U^{\mathrm{2}}}{h}}\right)+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left({\displaystyle \frac{UV}{h}}\right)\\ \text{(10a)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}UV\tan \mathit{\varphi }}{ah}}-\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }V=-{\displaystyle \frac{g}{\mathrm{2}a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial h^{\mathrm{2}}}{\partial \mathit{\lambda }}}{\displaystyle }\\ {\displaystyle \frac{\partial V}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}\left({\displaystyle \frac{UV}{h}}\right)+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left({\displaystyle \frac{V^{\mathrm{2}}}{h}}\right)-{\displaystyle \frac{(U^{\mathrm{2}}-V^{\mathrm{2}})\tan \mathit{\varphi }}{ah}}\\ \text{(10b)}& {\displaystyle }& {\displaystyle }+\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }U=-{\displaystyle \frac{g}{\mathrm{2}a}}{\displaystyle \frac{\partial h^{\mathrm{2}}}{\partial \mathit{\varphi }}}{\displaystyle }\\ \text{(10c)}& {\displaystyle \frac{\partial h}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}\left[{\displaystyle \frac{\partial U}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \left(V\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]=\mathrm{0},\end{array}$$

where U=hu, V=hv and h is the total layer thickness. The numerical scheme employs a standard finite-difference shallow-water solver in which the time-differencing follows a leapfrog scheme (center differencing in both time and space). The computations were carried out on an Arakawa C-grid. The model contains provisions for a temporal Robert–Asselin filter, but the filter's coefficient was set to zero in the simulations of the present section. In addition, the model includes no diffusion/viscosity terms.

The computational domain is $-\mathrm{180}{}^{\circ }\le \mathit{\lambda }\le \mathrm{180}{}^{\circ }$ and $-\mathrm{30}{}^{\circ }\le \mathit{\varphi }\le \mathrm{30}{}^{\circ }$. The boundary conditions are periodicity at the zonal boundaries $\mathit{\lambda }=\pm \mathrm{180}{}^{\circ }$ and vanishing meridional velocity at the channel's boundaries $\mathit{\varphi }=\pm \mathrm{30}{}^{\circ }$. The corresponding values of h and U at the boundaries are determined by the differential equations. For the chosen wave parameters the amplitude of the meridional velocity $\widehat{v}$ in Eq. (6) has an e-folding latitude of 11^∘, and its amplitude at $\mathit{\varphi }=\pm \mathrm{30}{}^{\circ }$ decays to $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$ of its maximal value; therefore, the velocity outside the computational domain can be comfortably neglected. The grid spacing and time step are $\mathrm{\Delta }\mathit{\lambda }=\mathrm{\Delta }\mathit{\varphi }=\mathrm{0.5}{}^{\circ }$ and Δt=600 s, which were found to yield stable solutions for at least 100 wave periods.

Figure 1 shows the initial u, v, Φ, ξ, and δ fields (where ξ and δ are the relative vorticity and divergence, respectively) of the chosen Rossby wave mode (Fig. 1a–e), and the resulting latitude–time (Fig. 1f–j) and time–longitude (Fig. 1k–o) Hovmöller diagrams of the simulated solution. The initial fields were obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. The simulated solutions were obtained using the abovementioned equatorial channel model. The chosen intersects used in the calculation of the Hovmöller diagrams are indicated by white dashed lines superimposed on the initial fields, and are also provided in the figure's caption. For the sake of legibility the time domain shown in each panel is only the last wave period of the simulation, i.e., $\mathrm{99}T\le t\le \mathrm{100}T$, where T is the wave period. The fields are normalized on their global maximum at t=0. Thus, white regions correspond to times at which the simulated solution exceeds the initial wave amplitude, momentarily. With this in mind, recall that the patterns in the latitude–time diagrams are similar to those in latitude–longitude diagrams, and can therefore be used for comparison with the initial fields. In general, the initial wave structure is preserved and the dominant slope in the time–longitude diagrams corresponds to the analytic slope indicated by dashed white lines (Fig. 1k–o). There are, however, some noticeable deviations: a slight east–west tilt can be observed in the latitude–time diagrams (Fig. 1f–j), but most egregiously, the divergence field is less regular than the other four. We return to this last point at the end of Sect. 4.3. The phase of the simulated patterns in the latitude–time diagrams fit the expected patterns considering the westward propagation of the Rossby mode at $\mathit{\lambda }=-\mathrm{18}{}^{\circ }$ in 1 wave period after 99 wave periods.

Figure 1 (a–e) The initial u, v, Φ, ξ, and δ Rossby wave fields, obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. (f–j) Latitude–time Hovmöller diagrams of the simulated solutions, obtained by intersecting the fields at $\mathit{\lambda }=-\mathrm{18}{}^{\circ }$ (also indicated by the white vertical dashed lines in a–e). (k–o) Time–longitude Hovmöller diagrams of the simulated solutions, obtained by intersecting u at $\mathit{\varphi }=\mathrm{0}{}^{\circ }$ and all other fields at $\mathit{\varphi }=\mathrm{9}{}^{\circ }$ (also indicated by the white horizontal dashed lines in a–e). The simulated solutions were obtained using the equatorial channel finite-difference model. The fields are normalized on their global maximum at t=0. The wave period for the chosen wave parameters is T=18.5 d. Contour levels range from −1.0 to +1.0 in intervals of 0.2.

Figure 2 (a–e) The initial $u,v,\mathrm{\Phi }$, ξ, and δ EIG wave fields, obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. (f–j) Latitude–time Hovmöller diagrams of the simulated solutions, obtained by intersecting the fields at $\mathit{\lambda }=-\mathrm{18}{}^{\circ }$ (also indicated by the white vertical dashed lines in a–e). (k–o) Time–longitude Hovmöller diagrams of the simulated solutions, obtained by intersecting v at $\mathit{\varphi }=\mathrm{9}{}^{\circ }$, δ at $\mathit{\varphi }=\mathrm{0}{}^{\circ }$ and all other fields at $\mathit{\varphi }=\mathrm{15}{}^{\circ }$ (also indicated by the white horizontal dashed lines in a–e). The simulated solutions were obtained using the equatorial channel finite-difference model. The fields are normalized on their global maximum at t=0. The wave period for the chosen wave parameters is T=1.9 d. Contour levels range from −1.0 to +1.0 in intervals of 0.2.

Similarly, Fig. 2 shows the initial u, v, Φ, ξ, and δ fields of the chosen EIG wave mode (Fig. 2a–e), and the resulting latitude–time (Fig. 2f–j) and time–longitude (Fig. 2k–o) Hovmöller diagrams of the simulated solution. Note that under the normalization used here the initial v field is independent of the wave type and is therefore identical to the one in Fig. 1. In contrast to Fig. 1 the patterns in the latitude–time diagrams of the simulated solutions are noticeably out of phase. However, considering the agreement between the dominant slope in the time–longitude diagrams and the analytic slope indicated by the dashed white lines (Fig. 2k–o), it is reasonable to say that this phase shift only results from a small phase speed error that accumulates over time. In addition, in contrast to the Rossby wave in Fig. 1, the divergence field in this case is just as regular as the other four fields.

The structure error defined in Eq. (8) is shown in Fig. 3 for both Rossby (Fig. 3a) and EIG (Fig. 3b) waves as a function of time. In both cases the structure error fluctuates about a mean value of less than 1 % and there is no visible trend throughout the simulation time of 100 wave periods. Recall that the structure error defined in Eq. (8) is insensitive to phase differences.

Figure 3 The structure error defined in Eq. (8) for both the Rossby (a) and EIG (b) waves as a function of time. Blue: calculated for the velocity vector $\mathit{u}=(u,v)$. Red: calculated for the geopotential Φ.

4.2 Demonstration using a global-scale spectral model

To demonstrate the applicability of the Matsuno wave as a test case for global-scale models we use the Geophysical Fluid Dynamics Laboratory's (GFDL's) spectral transformed shallow-water model which uses spherical harmonics as its basis functions (https://www.gfdl.noaa.gov/idealized-spectral-models-quickstart/). The chosen spectral resolution was T85, i.e., a triangular truncation where both the highest retained wave number and the total wave number equal 85. The chosen time step was Δt=600 s, as in the equatorial channel model. The model contains provisions for hyper-diffusion terms as well as a temporal Robert–Asselin filter, but the coefficients of both were set to zero for the simulations described below.

Figure 4 (a–e) The initial u, v, Φ, ξ, and δ Rossby wave fields, obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. (f–j) Latitude–time Hovmöller diagrams of the simulated solutions, obtained by intersecting the fields at $\mathit{\lambda }=-\mathrm{18}{}^{\circ }$ (also indicated by the white vertical dashed lines in a–e). (k–o) Time–longitude Hovmöller diagrams of the simulated solutions, obtained by intersecting u at $\mathit{\varphi }=\mathrm{0}{}^{\circ }$ and all other fields at $\mathit{\varphi }=\mathrm{9}{}^{\circ }$ (also indicated by the white horizontal dashed lines in a–e). The simulated solutions were obtained using GFDL's global-scale spectral model. The fields are normalized on their global maximum at t=0. The wave period for the chosen wave parameters is T=18.5 d. Contour levels range from −1.0 to +1.0 in intervals of 0.2.

Figure 4 shows the initial u, v, Φ, ξ, and δ fields (where ξ and δ are the relative vorticity and divergence, respectively) of the chosen Rossby wave mode (Fig. 4a–e), and the resulting latitude–time (Fig. 4f–j) and time–longitude (Fig. 4k–o) Hovmöller diagrams of the simulated solution. The initial fields were obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. The simulated solutions were obtained using the abovementioned GFDL global-scale spectral model. The chosen intersects used in the calculation of the Hovmöller diagrams are indicated by white dashed lines superimposed on the initial fields, and are also provided in the figure's caption. For the sake of legibility the time domain shown in each panel is only the last wave period of the simulation, i.e., $\mathrm{99}T\le t\le \mathrm{100}T$, where T is the wave period. The fields are normalized on their global maximum at t=0. Thus, white regions correspond to times at which the simulated solution exceeds the initial wave amplitude, momentarily. With this in mind, recall that the patterns in the latitude–time diagrams are similar to those in latitude–longitude diagrams, and can therefore be used for comparison with the initial fields. Indeed, the patterns in the latitude–time diagrams of the simulated solutions agree quite accurately with those of the initial wave structure, but are noticeably out of phase. Nevertheless, considering the agreement between the dominant slope in the time–longitude diagrams and the analytic slope indicated using the dashed white lines (Fig. 4k–o), it is reasonable to say that this phase shift results from a small phase speed error that accumulates over time. In addition, the divergence field is less regular than the other four fields. We return to this point at the end of Sect. 4.3.

Similarly, Fig. 5 shows the initial u, v, Φ, ξ, and δ fields of the chosen EIG wave mode (Fig. 5a–e), and the resulting latitude–time (Fig. 5f–j) and time–longitude (Fig. 5k–o) Hovmöller diagrams of the simulated solution. Note that under the normalization used in the present paper the initial v field is independent of the wave type and is therefore identical to the one in Fig. 4. As in Fig. 4, the patterns in the latitude–time diagrams of the simulated solutions are noticeably out of phase, but the dominant slope in the time–longitude diagrams (Fig. 5k–o) agrees well with the analytic slope, indicating that the observed phase shift results from a small phase speed error that accumulates over time.

Finally, the structure error in Fig. 6 fluctuates about a mean value of less than 1 % and there are no visible trends throughout the 100 wave period simulations. Recall that the structure error defined in Eq. (8) is insensitive to phase differences.

Figure 5 (a–e) The initial u, v, Φ, ξ, and δ EIG wave fields, obtained using the analytic expressions from Sect. 2 and the wave parameters from Sect. 3.1. (f–j) Latitude–time Hovmöller diagrams of the simulated solutions, obtained by intersecting the fields at $\mathit{\lambda }=-\mathrm{18}{}^{\circ }$ (also indicated by the white vertical dashed lines in a–e). (k–o) Time–longitude Hovmöller diagrams of the simulated solutions, obtained by intersecting v at $\mathit{\varphi }=\mathrm{9}{}^{\circ }$, δ at $\mathit{\varphi }=\mathrm{0}{}^{\circ }$ and all other fields at $\mathit{\varphi }=\mathrm{15}{}^{\circ }$ (also indicated by the white horizontal dashed lines in a–e). The simulated solutions were obtained using the GFDL global-scale spectral model. The fields are normalized on their global maximum at t=0. The wave period for the chosen wave parameters is T=1.9 d. Contour levels range from −1.0 to +1.0 in intervals of 0.2.

Figure 6 The structure error defined in Eq. (8) for both the Rossby (a) and EIG (b) waves as a function of time. Blue: calculated for the velocity vector $\mathit{u}=(u,v)$. Red: calculated for the geopotential Φ.

4.3 Smoothness and stability

In this section we examine the generation of small-scale features and the stability of the proposed wave solutions in a similar fashion to that used in Thuburn and Li (2000) for the original Rossby–Haurwitz wave number 4.

In Thuburn and Li (2000), the generation of small-scale features and the potential enstrophy cascade is observed by examining the potential vorticity field, which generates tongues that wrap up around themselves and break the initial east–west symmetry. For the small wave amplitude $A={\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ used in the present work, the potential vorticity is dominated by the planetary vorticity which is 5–6 orders of magnitude (depending on the wave) larger than the relative vorticity. Thus, instead of the potential vorticity we examine the relative vorticity (as well as the geopotential). Figures 1–2, as well as Figs. 4–5, show the evolution of these two fields between t=99T and t=100T, where T is the wave period in each case. Clearly, both fields remain regular throughout the simulations and do not develop small-scale features such as those observed in Thuburn and Li (2000). Recall that the simulations in the present work were carried out without any diffusion/viscosity terms. Thus, the simulations remain stable for at least 100 wave periods with no need to remove potential enstrophy at the grid scale.

In order to examine the stability of the chosen initial waves we repeat the simulations of the previous section with an added perturbation (white noise) to the initial fields. We demonstrate the stability of the waves using only the global-scale model, which was found to yield more stable results when adding the perturbation.

Figures 7 and 8 show the initial fields of the perturbed Rossby (Fig. 7a–e) and EIG waves (Fig. 8a–e), respectively, and the resulting latitude–time (Figs. 7f–j and 8f–j) and time–longitude (Figs. 7k–o and 8k–o) Hovmöller diagrams of the simulated solution, obtained using the GFDL global-scale spectral model. The initial perturbation in these figures consist of a uniformly distributed random white noise with an amplitude of 5 % of the field's amplitude added to each of the fields u, v, and Φ. Specifically, let q stand for any of the variables u, v or Φ, then the initial perturbation is given by

$$\begin{array}{}\text{(11)}& {\displaystyle }q=q_{\mathrm{a}}+\mathrm{0.05}\underset{\mathit{\lambda },\mathit{\varphi }}{max}\left|q_{\mathrm{a}}\right|\left(\mathrm{2}R-\mathrm{1}\right),\end{array}$$

where q_a is the analytic solutions obtained as in Sect. 2, and R is a uniformly sampled random array with elements in (0,1) whose dimensions are the same as q_a (in the present work a different R was drawn for each of the three variables).

Figure 7 Same as Fig. 4, but for the perturbed Rossby wave.

Overall, the perturbed waves seem to be stable. The u, v and Φ fields are almost as regular as those of the non-perturbed waves, except for the zero contour. The small-scale features in the vorticity field of the perturbed Rossby wave smooth out with time, in contrast to the potential vorticity field of the Rossby–Haurwitz wave number 4. Conversely, the perturbed Rossby wave divergence field is completely eroded. The vorticity and divergence fields of the perturbed EIG wave are not as regular as those of the non-perturbed wave. However, they too become smoother with time and the initial wave remains the most dominant wave throughout the entire 100 wave period simulation. The structure error in Fig. 9 is similar to the previous ones in Figs. 3 and 6. These results are quite surprising, as we would expect a sufficiently large perturbation to excite other modes, regardless of the waves' stability.

Figure 8 Same as Fig. 5, but for the perturbed EIG wave.

Figure 9 Same as Fig. 6, but for the perturbed Rossby (a) and EIG (b) waves.

Both the non-perturbed Rossby wave in Sect. 4.1 and 4.2, and the perturbed Rossby wave in the present section indicate that the divergence field is more sensitive than the other four fields of the Rossby wave. An immediate suspect in this regard is the divergence field amplitude, which is small for the chosen Rossby wave. For reference, the meridional wind amplitude for the chosen waves parameters (of both the Rossby and EIG waves) is $\mathrm{6.4}\times {\mathrm{10}}^{-\mathrm{6}}$, whereas the Rossby wave divergence field amplitude is $\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{12}}$. Conversely, the divergence field amplitude is only 1 order of magnitude smaller than the vorticity field amplitude, which is $\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{11}}$. Regardless of the cause, the fact that all other four fields remain quite regular while the divergence field is completely eroded suggests that the small but significant divergence field described by Phillips (1959) is in fact a small and insignificant divergence field.

4.4 Convergence test for the linear shallow-water models

In addition to the test cases proposed by Williamson et al. (1992) a resolution convergence test of linearized SWEs in which the simulations are compared to higher-order simulations is also useful for ensuring that the errors decrease with the increase in resolution. In this section we demonstrate that Matsuno's analytic wave solutions can be used for this purpose. We use the equatorial channel model which can be easily turned into a linear shallow-water model.

Figure 10 shows the structure error in absolute value as a function of the grid spacing $\mathrm{\Delta }=\mathrm{\Delta }\mathit{\lambda }=\mathrm{\Delta }\mathit{\varphi }$, from $\mathrm{\Delta }=\mathrm{2.5}{}^{\circ }$ to $\mathrm{\Delta }=\mathrm{0.25}{}^{\circ }$ every 0.25^∘. For each resolution, the initial non-perturbed waves were integrated for 100 wave periods. As an estimate of the structure error at each resolution we use the time-series averages (indicated by dots). The error bars were estimated using the standard deviations of the entire time series. As the resolution increases from $\mathrm{\Delta }=\mathrm{2.5}{}^{\circ }$ to $\mathrm{\Delta }=\mathrm{0.25}{}^{\circ }$, the structure error time-series average decreases from about 2 % to less than 1 %, and the standard deviation decreases from about 2 % to about 0.5 %. The time step across all resolutions in this figure was held fixed at Δt=100 s. Note that all of the results of the previous sections were obtained for $\mathrm{\Delta }=\mathrm{0.5}{}^{\circ }$ and Δt=600 s. This time step was found to yield convergent results for $\mathrm{\Delta }=\mathrm{0.5}{}^{\circ }$ in the sense that decreasing the time step by a factor of 2 yields no improvements. Nevertheless, for the convergence test in the present section we have further decreased the time step to Δt=100 s in order to allow a further increase in the spatial resolution by another 0.25^∘. This has also enabled a comparison with the results of the previous sections, thus ensuring that the simulations remain stable. Needless to say, for any time step one can expect to encounter numerical instabilities at some (high) spatial resolution.

Figure 10 Structure error in absolute value as a function of the grid spacing $\mathrm{\Delta }=\mathrm{\Delta }\mathit{\lambda }=\mathrm{\Delta }\mathit{\varphi }$, from $\mathrm{\Delta }=\mathrm{2.5}{}^{\circ }$ to $\mathrm{\Delta }=\mathrm{0.25}{}^{\circ }$ every 0.25^∘. The points correspond to the time averaged structure error over 100 wave periods, and the error bars are determined from the standard deviation. Blue: calculated for the velocity vector $\mathit{u}=(u,v)$. Red: calculated for the geopotential Φ.

5 Concluding remarks

As vertical resolutions in atmospheric and oceanic models increase it is essential to assess the accuracy with which they resolve baroclinic wave modes, typified by low gravity-wave phase speeds, in addition to the barotropic mode. To this end we propose to use a similar procedure to that used in the Rossby–Haurwitz test case but with different initial conditions. Instead of using the analytic solutions obtained by Haurwitz (1940), which are only accurate for high gravity-wave speeds such as those of the barotropic mode, we propose the use of the analytic solutions obtained by Matsuno (1966), which are accurate for lower gravity-wave speeds such as those of the baroclinic modes.

While the solutions from Matsuno (1966) apply for the equatorial β-plane, they approximate the solutions of the SWEs on the sphere for the speeds of gravity waves found in the baroclinic modes in the atmosphere, and as demonstrated in the present work can be accurately simulated in both equatorial channel and global-scale models in spherical coordinates. In addition, unlike the original Rossby–Haurwitz wave number 4, the chosen initial waves of the present test case remain stable for at least 100 wave periods, which for the chosen Rossby wave correspond to about 1850 d.

While the solutions of the SWEs obtained by Matsuno (1966) account for the small divergence field missing from the non-divergent Rossby–Haurwitz waves, the results of the present study suggest that this missing divergence field is insignificant.

Ideally, we expect the proposed test case to stand on an equal footing alongside the Rossby–Haurwitz test case, but in the words of Williamson et al. (1992): “The test will only become standard to the extent that the community finds it useful”.