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 atmospheric and oceanic 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 small speeds of gravity waves such as those of the baroclinic modes in the atmosphere and ocean. 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 large speeds of gravity waves 5 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 amplitude) for both planetary and inertia gravity waves, and confirms the accuracy of the simulation by employing Hovmöller diagrams and temporal and spatial spectra. The proposed test case is successfully applied to a standard finite-difference, equatorial, non-linear, shallow water model in spherical coordinates, which demonstrates that Matsuno’s wave-solutions can be accurately simulated for at least 10 wave-periods, which for oceanic 10 planetary waves is nearly 1300 days. In order to facilitate the use of the proposed test case, we provide Matlab, Python and Fortran codes for computing the analytic solutions at any time on arbitrary latitude-longitude grids. Copyright statement.

with which the barotropic wave mode is resolved. 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), which approximate the solutions of the SWEs on the sphere in the asymptotic limit of small speed of gravity waves (De-Leon and Paldor, 2011;Garfinkel et al., 2017).
In addition to being on two opposite ends of the spectrum in terms of the relevant speeds of gravity wave, the solutions 15 obtained by Matsuno (1966) differ from those obtained by both Haurwitz (1940) and Paldor et al. (2013) in their meridional extent. The former become negligible outside a narrow equatorial band, whereas 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 20 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 multi-layer models via baroclinic instability. In contrast, the proposed test case is concerned with linear wave propagation in (non-linear) single-layer models.
While the customary use of the term baroclinic is associated with density variations in the vertical direction, here the same term is used to denote thin layers of fluid of homogeneous density where the gravity waves speeds are similar to those observed in 25 baroclinic modes in the atmosphere and oceans.
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 in the same spirit that Williamson et al. (1992) standardized the Rossby-Haurwitz one. We start with a short description of the analytic expressions derived by Matsuno (1966) in section 2. The proposed test procedure, including the choice of wave- 30 parameters and assessment criteria, is described in Section 3. We then demonstrate the usefulness of the proposed test case in section 4, using a standard finite-difference, equatorial, (non-linear) shallow water model in spherical coordinates.

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.
where x and y are the local Cartesian coordinates in the zonal and meridional directions, respectively; t is time; u and v are 5 the velocity components in the zonal and meridional directions, respectively; Φ is the geopotential height; k is the zonal wavenumber; ω is the wave-frequency; andû(y),v(y) andΦ(y) are the latitude dependent amplitudes. In accordance with the sign convention used in Matsuno we assume k is non-negative and let ω take any real value. Note, however, that the sign in front of ω in (1) is opposite to that in Matsuno's theory. 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. in the eastward direction.

10
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 for n = −1, 0, 1, 2, . . . , where Ω and a are Earth's angular frequency and mean radius, respectively; g and H are the reduced 15 gravity and equivalent depth.
For n ≥ 1 Equation (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 ω = − √ gHk, leads to infinite zonal wind and is thus discarded as a physically reasonable solution. The remaining two roots correspond to the lowest (i.e. 20 n = 0) EIG wave and the Mixed Rossby-Gravity (MRG) wave. For n = −1 Equation (2) has one real root ω = √ gHk, which correspond to the equatorial Kelvin wave (see Matsuno, 1966).
Having found (one of) the wave-frequencies for a given combination of n and k, the corresponding latitude dependent amplitudes can be written aŝ v n = ψ n (3a) for n = 1, 2, 3, . . . (the cases n = −1, 0 require special treatment), where Here = (2Ωa) 2 /gH is Lamb's parameter, A is an arbitrary amplitude, andĤ n are the normalized Hermite polynomials of 10 degree n. Note: (i) The chosen normalization for the latitude dependent amplitudes in (3) is different from the one used in Matsuno. We use the above normalization for convenience, as it relates the amplitude ofv to that of ψ in a straight forward way and guarantees that it is independent of both k or ω. (ii) The use of the normalized version of the Hermite polynomials also leads to slightly different pre-factors in front of ψ n+1 and ψ n−1 compared to Matsuno. However, they are generally more stable and more convenient for the spectral analyses employed in the following sections.

15
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 1/ 1/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 small speeds of gravity waves √ gH.

20
In practice, to use Matsuno's solutions in spherical models, the local Cartesian coordinates x and y in the above formulae (1) and (4) have to be replaced by the longitude λ and latitude φ of the geographical coordinate system. Recall that the transformation between the two system is (x, y) → a(cos φ 0 λ, φ), where φ 0 is the central latitude to which the β-plane approximation is applied. Thus, for the equatorial β-plane where φ 0 = 0, the transformation is simply (x, y) → a(λ, φ). Likewise, the planar wavenumber k in all of the formulae (1)-(4) has to be replace by its spherical counterpart using the transformation 25 k → k/a cos φ 0 = k/a.
For the small gravity wave speeds (i.e. small values of ) used in the present work, the eigenfunctions in (4) become negligible outside a narrow equatorial band (see Sections 3.1 and 4 below). Thus, the above expression can be used to test both globalscale and tropical models. Finally, using the above formulae to calculate the waves' frequencies and latitude dependent amplitudes requires routines for finding the roots of the cubic equation (2), and for evaluating the normalized Hermite polynomials on arbitrary latitudelongitude grids. The roots of the cubic equation can be obtained in a closed analytic form using the solutions of the general cubic equation as detailed in appendix A. The normalized Hermite polynomials, in turn, can be evaluated using their threeterm recurrence relation, as in Press et al. (2007). In order to facilitate the application of the proposed test case we provide as 5 part of the online supplementary material Matlab, Python and Fortran codes (named matsuno.m, matsuno.py and matsuno.f90, respectively) for computing the analytic fields on arbitrary latitude-longitude grids.

Proposed test procedure
The general procedure of the proposed test case is similar to the Rossby-Haurwitz one 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.

10
The initial wave fields in this case are taken from the analytic expressions in Section 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.

wave-parameters
The wave-parameters consist of the speed of gravity waves √ gH, the wave-number and wave-mode k and n, the waveamplitude A, and the wave-type. Any given combination of these parameters completely specify a unique wave using the expressions in (2)-(4). We consider all other parameters, including the spatio-temporal 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 modus 20 operandi. However, as noted in Polvani et al. (2004), different choices for the form of diffusion/viscosity terms correspond to different sets of equations and may not converge to the same solutions.
We offer two different choices for √ gH inspired by the observed speeds of gravity waves of the baroclinic modes in the atmosphere and oceans. In practice we keep g fixed to Earth's gravitational acceleration, and vary the speed of gravity waves by varying H. For atmospheric models we suggest using H = 30 m, which is within the range of observed equivalent depths in 25 the equatorial atmosphere (Wheeler and Kiladis, 1999). For ocean models we suggest using H = 0.5 m, which corresponds to the observed speed of gravity waves of the first baroclinic mode in the oceans (Chelton et al., 1998). As mentioned in section 2, the analytic solutions obtained by Matsuno for the equatorial β-plane are only accurate approximations of the SWEs on the sphere in the asymptotic limit of small speeds of gravity waves. The above values were found by trial and error to be sufficiently accurate in the sense that they yield stable integrations for at least 10 wave-periods in the simulations demonstrated in Section In addition to the speed of gravity waves, the accuracy of Matsuno's solutions depend on the wave-number and wave-mode as well. For a given value of √ gH, they become asymptotically accurate in the limits k, n → 0 (but k = 0) (De-Leon and Paldor, 2011). In addition, the higher the wave-number or wave-mode are, the greater the spatial variability and the required spatial resolution are. 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 = 5 and n = 1, which are within the range of 5 dominant values observed in the equatorial atmosphere (Wheeler and Kiladis, 1999), but other choices may work just as well provided k 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 waves-amplitude should be sufficiently small so as to satisfy the linearization condition. The proposed amplitude of Equation (4) is A = 10 −8 m s -1 , chosen by trial and error so as to enable stable integrations for at least 10 wave periods in the 10 simulations in Section 4.
In general, there are two qualitatively different wave types, the Rossby and IG wave types, that differ in the magnitude of their divergence field. In order to assess the models' performances in these two qualitatively different limits we suggest using one of each. Since Rossby waves are exclusively westward propagating, we choose the EIG wave from the two IG waves as the second one in order to eliminate potential longitudinal biases.

Assessment criteria
For sufficiently small wave-amplitudes we expect the spatio-temporal structure of the simulated solutions to be that of zonally propagating waves, i.e.ξ(φ)e i(kλ−ωt) (where ξ 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 20 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 longitude-time diagram obtained by intersecting the fields at a certain latitude. The contour lines in the longitude-time plane are the set of points satisfying kλ − ωt = const (for some real const). Thus, the expected pattern for this diagram is that of straight lines whose slopes equal the wave's phase speed k/ω. In order to avoid small fluctuations in the vicinity of 25 latitudinal zero-crossings, we recommend using latitudinal intersects at or near local extrema.
(ii) A latitude-time diagram obtained by intersecting the fields at a certain longitude. For any two wave-fronts with equal phase k(λ 2 − λ 1 ) = ω(t 2 − t 1 ). Thus, holding λ fixed while varying t from t 1 to t 2 is equivalent to holding t fixed and varying λ from λ 1 to λ 2 = λ 1 + ω/k(t 2 − t 1 ). The resulting pattern is similar to that of a latitude-longitude diagram, but provides a testament of the time evolution as opposed to a momentary snapshot.  the difference between the two provides a quantitative error-measure. The Fourier analysis can be trivially found using readily available Fast-Fourier-Transform libraries. The Hermite analysis was found in the present work using the procedure described in Appendix B.
In order to attain reasonable accuracy in terms of the spectral analysis, it is recommended to integrate the initial fields forward in time for at least 10 wave-periods, with a sampling frequency of about 10 samples per period. The proposed integration and 5 sampling times (denoted by T f and T s ) for each of the four simulations, along with the wave-period (denoted by T ), are given in Table 1.

Demonstration
In order to demonstrate the usefulness of the Matsuno test case we run the proposed procedure using a simple finite-difference shallow water model. The model is a spherical version of the Cartesian model used in Gildor et al. (2016), in which the EIG, H = 0.5 m: 5.7 60 12 Table 1. The wave-period T , integration time T f , and sampling spacing Ts. Note, T is derived from the chosen test wave-parameters and is a characteristic of the waves, whereas T f and Ts where chosen to attain reasonable accuracy in terms of the spectral analysis. integration forward in time is carried out using the transport form of the SWEs where U = hu, V = hv and h is the total layer thickness. The numerical scheme employs a standard finite difference shallowwater solver in which the time-differencing follows a leapfrog scheme (centre differencing in both time and space). The computations were done on an Arakawa C-grid. The original Matlab code used in the following simulations is available as part of the on-line supporting material.  The resulting Hovmöller diagrams of the simulated solutions for each of the initial waves are shown in Figure 2. In all cases, the latitude-time diagrams were obtained by intersecting v at λ = −18 • (also indicated by vertical dashed white lines in Figure   1). The longitude-time diagrams were obtained by intersecting v at φ = 9 • for H = 30 m, and φ = 4 • for H = 0.5 m (the latter case is indicated by horizontal dashed white lines in Figure 1). Similar results are obtained using u or Φ, provided the meridional intersects are taken in the vicinity of their local extrema. For the sake of legibility the shown time domain in each 5 panel is 9T ≤ t ≤ 10T , where T is the corresponding wave-period provided in Table 1. Note that the latitude-time diagrams are similar to that of a latitude-longitude diagram and can, therefore, be used to compare with the initial fields. In all cases the initial wave-structure is clearly discernible after 10 wave-periods and the dominant slope corresponds to the analytic slope indicated with dashed white lines. For a more quantitative assessment we use Fourier and Hermite analyses in order to locate the dominant wave-numbers, wave-modes, and wave-frequencies, and compare their amplitudes with the initial wave amplitudes. Figure 3 contains the resulting power-spectra of the simulated solutions obtained using the same latitudinal and longitudinal intersects used for the Hovmöller diagrams in Figure 2, and t = 10T as the time intersect, where T in each case is the corresponding wave-period listed in Table 1. In all cases the dominant wave-number and wave-mode clearly match the initial values k = 5 and n = 1.

5
Likewise the dominant frequencies match the initial ones indicated on the figure using black dashed lines.
|F k |, |Fn| and |Fω| denote the amplitudes of the wave component in the k, n and ω spaces, respectively.
The overall scenario that emerges from the spectral analyses is similar to that of the Hovmöller diagrams, namely that the initial waves remain dominant after 10 wave-periods. However, it provides a quantitative measure of the errors in the simulated solutions, since for an ideal model the amplitudes of the first most dominant wave components should be the same as the analytic ones. Thus, the closer the amplitudes of the first most dominant wave components to the analytic ones, and the smaller the amplitudes of the second most dominant wave components, the better.

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, in addition to the barotropic mode. To this end we propose to use a similar procedure to the one used in the Rossby-Haurwitz test case but replace the initial conditions. Instead of using the analytic solutions obtained by Haurwitz (1940), which are only accurate for large gravity wave speeds such as those of the barotropic mode, we propose to use the analytic solutions obtained by Matsuno (1966), which are accurate for smaller gravity wave speeds such as those of the baroclinic modes.
Unlike the solutions obtained by Haurwitz (1940) for the non-divergent barotropic vorticity equation, the solutions of the SWEs obtained by Matsuno (1966) fully account for the velocity divergence. As a result an initial wave can be accurately simulated for long times and the simulated solution can be compared to the analytic solution to obtain a quantitative assessment. 5 While Matsuno's solutions 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 oceans, and as demonstrated in the present work can be accurately simulated in spherical coordinate models. Finally, in contrast to the barotropic mode, the baroclinic modes are confined to narrow equatorial band. Therefore, the proposed test can also be used to test tropical models.
Ideally, we expect the proposed test case to stand on an equal footing alongside the Rossby-Haurwitz one, but in the words 10 of Williamson et al. (1992): "The test will only become standard to the extent that the community finds it useful".
Code availability. The following files are available online as part of the supplementary material of this work: matsuno.m, matsuno.py, and matsuno.f90: In order to facilitate the application of the proposed test case we provide Matlab, Python, and Fortran codes for computing the analytic solutions obtained by Matsuno as described in section 2. The code can be used to evaluate the horizontal velocity fields u and v, in m s -1 , and the geopotential field Φ, in m 2 s -2 (or h in m), on arbitrary latitude-longitude grids for all t ≥ 0.
shallow_water_model.m: A Matalb code containing the shallow water model used in section 4.

Appendix A: The roots of the cubic equation
For any combination of wave-number and wave-mode, k and n, the wave-frequencies are given by the three roots of the cubic equation in (2). These roots can be written in a closed analytic form using the general solution of the cubic equation (e.g. 20 Abramowitz and Stegun, 1964): where j stands for the three roots, and where ∆ 0 = 3 gHk 2 + 2Ω √ gH a (2n + 1) , Given the definitions in (A2), the explicit expressions for the frequencies of the Rossby, WIG and EIG waves are obtained by sorting the values in (A1) as follows: Rossby : ω n,k,R = − min j=1,2,3 |ω n,k,j |, Westward Inertia-Gravity : ω n,k,WIG = min j=1,2,3 ω n,k,j , Eastward Inertia-Gravity : ω n,k,EIG = max j=1,2,3 ω n,k,j . 5

Appendix B: The Hermite transform
Let ξ(φ) denote a latitudinal series (i.e. an intersect along specific longitude and time) of any one of the dependent variables.

20
Finally, the integral in (B2) was evaluated using the trapezoidal rule.