On spurious instabilities on the β-planes with no mean flows

Recent numerical calculations that suggest the existence of instabilities of the shallow water equations on the β-planes (both mid-latitudes and on the equator) are shown to be associated with spurious instabilities rather than genuine unstable modes. The spurious unstable modes arise in an analysis of the necessary condition for the existence of complex roots to the cubic dispersion relation but it is shown here that these complex roots are not associated with solutions (eigenvalues and eigenfunctions) of the corresponding eigenvalue equation.


Motivation
In Mckenzie (2009) a general condition for the existence of complex phase speeds was derived for the Linearized Shallow Water Equations (LSWE) on the β-plane in the absence of mean flow.The condition was also calculated (using very similar formulation to that of the LSWE) to stratified flows without rotation (where Brunt-Vaisala frequency, N (z), is substituted for Coriolis frequency, f (y)) and also on the equator (where f 0 = 0).This note attempts to cast the formulation of the equations in a more traditional way and examine the range of parameter values where the calculation of the instabilities reported in Mckenzie (2009) is valid.These instabilities defy the necessary conditions for instability link derived in the last 70 years between the existence of unstable perturbations and an energetic source for their growth (see e.g.Kuo, 1949;Pedlosky, 1982).The most common of these sources is a mean flow that can transfer its kinetic energy under certain conditions (e.g. an inflection point of the vorticity) to the unstable perturbation modes.The aim of this note is to clarify the source of the apparent contradiction between the classical GFD instability theory and the numerical calculations of McKenzie (2009).

Energy integral for the Linearized Shallow Water Equations (LSWE) with rotation
For arbitrary latitude-dependent Coriolis frequency, f (y), the LSWE with rotation are: A trivial multiplication of the x-momentum equation (first equation of system 1) by Hu, the y-momentum equation (second equation of system 1) by Hv and the continuity equation (third equation of system 1) by gη (note that g and H stand for the reduced gravity and reduced height, respectively, so the speed of gravity waves, (gH ) 1 / 2 , is a parameter that can vary from 2-3 m s −1 in a baroclinic ocean to 200-300 m s −1 in a barotropic ocean/atmosphere).
Adding the resulting equations one gets: When this equation is integrated over the entire domain (including infinity) the RHS vanishes for suitable boundary conditions (integration of the divergence operator retains only values of the integrand at the boundaries) and the result is: This straightforward calculation (which can be found in any textbook) shows that regardless of the form assumed for f (e.g.f-plane, β-plane) in the SWE with no mean flow, linear instabilities can develop, only if perturbation energy is advected inward from the domain's boundaries.The way by which this advection through the boundaries generates the unstable modes is the subject of many classical instability studies the resulted in a wealth of examples associated, for example, with potential vorticity jumps at the boundaries such as in Rayliegh problem (see Paldor et al., 2009) or the inflection-point conditions (Kuo, 1949).
Published by Copernicus Publications on behalf of the European Geosciences Union.

Instability on the β-planes without mean flow
Solutions system (1) on the β-plane (where f = f 0 + βy = 2 sinφ 0 + 2 cosφ 0 y/R with -earth's frequency of rotation, R -earth's radius and φ 0 -the latitude where the tangential plane touches the spherical earth) are usually assumed to represent zonally propagating waves i.e. u, v and η are assumed to vary in x and t as e ik(x−Ct) .The case R −1 = 0 (i.e. when f (y) = f 0 so β = 0) is called the f-plane and the case φ 0 = 0 is the equatorial β-plane.In all these cases a cubic equation determines the relationship between the phase speed C and the wavenumber, k and the remaining parameters: , R, g (gravity or reduced gravity) and H appear in the coefficients of this relationship, known as the dispersion relation.When the equations are nondimensionalized the four-dimensional parameters g, H , and R combine into a single non-dimensional parameter α = gH /(2 R) 2 .The dispersion relation, which is usually written as either C(k) or kC(k), is obtained when the equations for the y-dependent amplitudes of the wavelike solutions are transformed into a second order (the equation is always second order in y because there is no y-derivative in the linear x-momentum equation) eigenvalue equation whose eigenvalues (which are determined by the boundary conditions) are combinations of all the parameters that appear in the equations.Denoting the eigenvalues by E (in the case of harmonic waves in a channel of non-dimensional width φ the eigenvalues are (nπ/ φ) 2 ).
On the mid-latitude β-plane theory the cubic dispersion relation takes the form (see Eq. 3.1 of Paldor et al., 2007, PRM hereafter): On the equatorial β-plane the cubic dispersion relation is give by the roots of: where ω = kC (see Matsuno, 1966, andErlick et al., 2007).
A similar cubic dispersion relation was derived in LeBlond and Mysak (1978) but solutions for the eigenvalue, E, were only derived on the f-plane.
The lower bound of m = 2/3 3/2 , above which complex values of C are encountered, probably corresponds to the M > 1 condition derived in McKenzie (2009).
Transforming the condition on m(= ba −3/2 ) into a condition on a and b yields: In a mid-latitude channel (a = α + E sin 2 (φ 0 )/k 2 ; b = cos(φ 0 )α/k 2 ) this condition translates to: Defining x = k 2 and rearranging, one gets the following necessary condition for instability: The RHS of Eq. ( 7) becomes infinitely large at both small and large x(= k 2 ) so this condition can only be satisfied (i.e.complex roots exist) at O(1) x-values.The minimal value of the RHS of Eq. ( 7) occurs at: Substituting this minimal value of x into the RHS of Eq. ( 7) one gets the necessary condition for a complex root to occur: Dividing through by α and rearranging, yields: In interpreting the implications of condition (10) one should remember that the eigenvalue, E, depends on α.In the simple case of harmonic waves in a mid-latitude channel (where the term βy set equal to 0; see PRM) the E n in a channel of width 2δφ (2L = 2δφR) are (see Eqs. 2.10a and 4.1 in PRM): Substituting this relationship into Eq.( 10) one obtains the condition: Rearranging, this relationship yields the condition: Viewing this last relation as a quadratic inequality in √ α and requiring that the inequality be satisfied at least at the minimum of the parabola on the LHS yields: In mid-latitudes this condition cannot be satisfied since δφ <min(φ 0 , π/4) for a channel located in one hemisphere.
In the n = 0 case (where the RHS attain their lowest values) the physical interpretation of Eq. ( 14), is that the ratio between the change in the value of Coriolis parameter across the domain (∝ cosφ 0 δφ) and its mean value (∝ sinφ 0 ) is larger than π.Thus, the necessary condition for instability on the mid-latitude β-plane is satisfied only when the very essence of the β-plane approximation (where sinφ is expanded to first order only near some φ 0 > 0) loses its validity.
For waves in a channel on the equatorial β-plane the eigenfunctions are Hermite Functions and the eigenvalues are E n = (2n + 1)α −1/2 .As noted below Eq. ( 4) The equatorial counterpart of condition ( 10) is obtained by substituting α for sin 2 (φ 0 ) in the expression for a and setting cos(φ 0 ) = 1 in the expression for b.Carrying out these substitutions in Eq. ( 10) one obtains: which cannot be satisfied for any integer, n = 0, 1, 2, . . . .Thus, in both equatorial and mid-latitude channels the necessary condition for the existence of (the spurious) instabilities, Eq. ( 10), is satisfied only for physically unacceptable parameter values or for a non-existing mode (i.e., n < 0).

Concluding remarks
The conclusion that can be drawn from the trivial calculations presented above is that inferences based solely on the (cubic) dispersion relation can be irrelevant to the instability problem under study.Only a complete solution of the eigenvalue problem associated with the dispersion relation (such as determining the values of E) can yield genuine instabilities and in the absence of such a solution the analysis is merely an algebraic calculation of the complex roots of a cubic C(k) relation.In the case addressed in McKenzie (2009) an attempt to evaluate the u y (v in the standard GFD notation) eigenfunction that solves his Eq. ( 19) subject to acceptable boundary conditions would have clarified that no instability exist in the shallow water equations in the absence of a mean flow.
From another perspective, the shallow water equation can be formulated as a Hamiltonian system in which the temporal evolution is governed by properly defined Poisson bracket operator (Weinstein, 1983;Piterbarg and Schulman, 1989).Thus, the conservation of energy is an inherent and fundamental property of this system that governs its dynamics via the Hamiltonian and Poisson bracket operator and should not be viewed as a trivial result of the straightforward algebraic calculation presented in Sect. 2. Clearly, energy conservation and Hamiltonian form are not consistent with exponential temporal growth of infinitesimally small initial perturbations i.e., the satisfaction of the necessary conditions for the occurrence of complex roots of a cubic relation in C does not guarantee the existence of unstable modes.