Instability of combined gravity-inertial-Rossby waves in atmospheres and oceans

The properties of the instability of combined gravity-inertial-Rossby waves on a β-plane are investigated. The wave-energy exchange equation shows that there is an exchange of energy with the background stratified medium. The energy source driving the instability lies in the background enthalpy released by the gravitational buoyancy force. It is shown that if the phase speed of the westward propagating low frequency-long wavelength Rossby wave exceeds the Poincaŕ e-Kelvin (or “equivalent” shallow water) wave speed, instability arises from the merging of Rossby and Poincaŕ e modes. There are two key parameters in this instability condition; namely, the equatorial/rotational Mach (or Froude) numberM and the latitudeθ0 of the β-plane. In general waves equatorward of a critical latitude for given M can be driven unstable, with corresponding growth rates of the order of a day or so. Although these conclusions may only be safely drawn for short wavelengths corresponding to a JWKB wave packet propagating internally and located far from boundaries, nevertheless such a local instability may play a significant role in atmosphere-ocean dynamics.


Introduction
In this paper we extend some aspects of an instability of combined gravity-inertial-Rossby waves on a β-plane (McKenzie, 2009) to include the equatorial waveguide and the case of prescribed meridional wave numbers.The instability arises from the merging of westward propagating Rossby and Poincaré modes.The two key parameters in the problem, namely the equatorial rotational Mach (or Froude) number M and the latitude θ 0 on which the β-plane is constructed, define, through the marginal stability condition, the unstable region as lying equatorward of a critical latitude for any given M.
The governing equations given in the next section reduce to a single partial differential wave equation which, for Fourier type plane wave modes, becomes an ordinary second order differential equation describing the latitudinal structure of the perturbations.It is shown that the wave energy equation may be cast into a wave energy-exchange form from which it is evident that there is indeed an energy exchange with the inhomogeneous, stratified background.The energy source which may drive an instability is identified as the background enthalpy released by gravitational buoyancy.
In Sect. 4 we present a local stability analysis based on JWKB solutions of the wave structure equation.These solutions describe the propagation of "short wavelength" wave packets located well within boundaries and therefore boundary terms (conditions) are surely irrelevant to the properties of these waves.The corresponding local dispersion equation shows that a wave coupling instability sets in when the low frequency-long wavelength Rossby wave speed exceeds the Poincaré-Kelvin (or "equivalent" shallow water) wave speed.If the β-plane is centred on the equator a waveguide system can be formed, and analytic solutions in terms of Hermite polynomials are available which yield eigenvalues which do not exhibit instability.However these solutions are good approximations only if the effective rotational Mach number is large, to ensure that the first reflection point lies well inside the Tropics.For moderate Mach numbers of O(1) the local JWKB solutions, appropriate to "short wavelength" wave packets, may exhibit instability.The mid-latitude case reiterates the previous analysis (McKenzie, 2009), however using a different normalization, which now also includes the equatorial case, as well as prescribed meridional wave numbers.In the latter case it is shown that the unstable region is further diminished to a fairly narrow belt around the equator for moderate values of that wave number.On the other hand for waves evanescent about a given latitude the unstable region is extended.
The work is summarized in Sect. 5 in which we refer to pertinent recent work (Maas and Harlander, 2007;Paldor et al., 2007).Further work on full wave solutions of the wave equation for freely propagating waves will be important in the possible extension of the present conclusions based on local JWKB solutions.It may well be that the instability is of a local nature and that globally the system is stable.However, that does not imply that the instability should be dismissed as insignificant.After all, the Schwarzschild criterion for atmospheric stability, N 2 > 0, predicts local instability (N 2 < 0), the onset of overturning and the development of the convection zone of the Sun which is itself globally stable.Hence, the non-linear evolution of this instability may play an important, hitherto unrecognised, role in atmosphere-ocean dynamics.

Governing equations
The linearized equations of motion for small amplitude perturbations about a background atmosphere stratified hydrostatically and rotating with frequency ẑ may be written (e.g.Eckart, 1960;Gill, 1982;Pedlosky, 1987) as in which q = ρ 0 u is the perturbation mass flux or momentum density (the background density ρ 0 multiplied by the perturbation fluid velocity u), p e is the pressure perturbation and ρ e g is the buoyancy force.In the β-plane approximation the Coriolis parameter (or frequency) f , at any latitude θ, is given by in which: where y = aδθ, θ = θ 0 + δθ and a is the radius of the planet.The spherical geometry is replaced by a β-plane constructed tangential to the planet at a given latitude θ 0 in which the local Cartesian co-ordinates are x directed eastward, y northward and z vertically.In the Boussinesq approximation the continuity equation assumes the incompressible form div q = 0. (4) The small amplitude density perturbation, ρ e , however evolves according to the buoyancy-adiabatic condition (Eckart, 1960;Lighthill, 1980), and where c 2 0 = γp 0 /ρ 0 .We shall assume that the atmosphere is stably stratified so that the square of the Brunt-Vaisälä frequency N 2 > 0. In effect the Boussinesq approximation ignores the variations in the perturbation density ρ e , unless multiplied by g as in the buoyancy force, and therefore filters out higher frequency acoustic waves.The background density ρ 0 (z) has a scale height H and is stratified according to where p 0 (z) is the background pressure distribution.Equations (1), ( 4) and (5) readily lead to the following wave equation for the system (e.g.McKenzie, 2009) This linear partial differential equation, which displays both dispersive (through the frequencies f and N ) and anisotropic (through the preferred directions ŷ and ẑ) properties, lends itself to Fourier plane wave analysis.In the simplified case with N assumed constant we may choose zonally and vertically propagating solutions of the form in which the latitudinal structure in y is governed by the second order differential equation which follows from substituting Eq. ( 9) into Eq.( 8).The wave number is given by Equation ( 10) admits classical JWKB solutions, appropriate to a slowly varying medium in the short-wave limit, of the form in which the local wave number k y is given by the local dispersion equation which follows from Eq. ( 11) and may be written in the wave normal surface form The wave normal surface is an ellipsoid or hyperboloid, whose axis of revolution is parallel to the k z axis, but displaced −β/2ω units along the k x axis, according as ω ≤ f and ω < N .In this form it is particularly useful in constructing ray trajectories for given ω, k x and k z , and with k y thus determined at each latitude, the ray direction is normal to the wave normal surface (Lighthill, 1980).On the other hand the local dispersion equation may also be written in the diagnostic form in which This is the appropriate form for any investigation of possible instabilities characterized by complex conjugate values for the frequency ω (McKenzie, 2009).This form is also appropriate to shallow water theory of an ocean of depth h in which V is then replaced by √ gh, the shallow water speed.We shall refer to V as the Poincaré-Kelvin speed (Pedlosky, 1987).
Before considering energy arguments (next section) and the stability analysis (Sect.4) we briefly recapitulate the waves given by Eq. ( 14).For wave frequencies ω ω i (the Poincaré frequency) Eq. ( 14) yields the Rossby wave dispersion equation On the other hand for higher frequencies Eq. ( 14) gives the Poincaré-inertial modes (modified by the β-effect), Note that Eq. ( 16) shows that the long wavelength (k 2 → 0) Rossby wave propagates westward at the speed ω/k x = βV 2 /f 2 , whereas the short wavelength -(k 2 → ∞) -high frequency Poincaré-inertial waves propagate at the Poincaré-Kelvin speed V .Intuitively it is clear that if βV /f 2 > 1 the Rossby wave may couple with a westward Poincaré mode over some intermediate band of frequencies and wave numbers.Analysis shows that this is indeed the criterion for the onset of instability.

Wave energy equation and the energy reservoir
The system of Eqs. ( 1), ( 4) and ( 5) possess a wave energy equation which follows in a straightforward fashion by taking the scalar product of the equation of motion (1) with the mass flux q and using Eqs.( 4) and ( 5) to guide the right hand side terms q • ∇p e and gρ e q z into yielding the conservation form (18) in which E is a measure of the wave energy density consisting of the kinetic energy (the first term) and the thermobaric energy (the second term) (Eckart, 1960), whilst F is a measure of the wave energy flux.The Coriolis term makes no contribution to the wave energy equation because the Coriolis "force" is perpendicular to q.Although the wave energy equation is in fact a redundant equation (since it follows from the system's equations), we quote and examine it here because of its obvious physical interest.Equation ( 18) when integrated over a volume V enclosed by a surface S takes the form (using Gauss's theorem on the divergence term) where dS is directed along the inward pointing normal to the surface.If the surface (or surfaces) enclosing the volume are rigid, requiring q n = 0 on the boundary so that the pressure p e from without can do no work on the volume of the fluid contained within, then the surface integral on the right hand side is indeed zero and E =const.However, in the case of an ocean, although the bottom may be regarded as rigid, the surface is free at which pressure must balance across it, and through which there is no flow; that is to say, the surface is a streamline.Hence there is a free surface contribution to the surface integral on the right hand side of Eq. ( 20) which cannot be put to zero.Therefore the wave energy measure E may change with time.Similarly in the case of an atmosphere bounded below by an ocean the surface integral cannot, a priori, be set to zero.In the case of an atmosphere over ground q n is indeed zero at ground, but conditions at large heights are subject to radiation type conditions with the result that again the right hand side cannot automatically be put to zero.Therefore in principle the wave energy Eq. ( 18) and its integral form (20) do permit the wave energy E to change with time and the wave system, described by Eq. ( 8), may admit unstable solutions which grow in time without violating a wave energy theorem.Apologies to "old school" fluid dynamicists who will regard these conclusions as merely the wave perturbation form of the general energy theorem which states that "the total energy of a volume of fluid increases at a rate equal to that at which work is being done on the boundary by pressure from without" (e.g.Lamb, 1932).Nevertheless, the further question arises as to what is the possible energy source which may drive unstable solutions of Eq. ( 8), if they exist.There is no obvious mean flow source to drive Kelvin-Helmholtz-baroclinic type instabilities.However observe that the background hydrostatic Eq. ( 7) may be cast in the form in which w 0 is the enthalpy given by The particular form for w 0 depends on the equation of state and whatever heating/cooling processes are active in the background state.For example, in polytropic (adiabatic) processes, p 0 ∝ ρ γ 0 , the enthalpy is given by In general Eq. ( 21) may be regarded as the energy equation for the background state which exhibits a vertical (negative) gradient of the enthalpy in a fashion analogous to the temperature gradient which may drive a thermal-convective instability.Hence enthalpy (the total thermodynamic heat content of a system) released by gravitational buoyancy is available.This energy source may be "tapped" to drive an instability provided the system admits of wave modes capable of releasing this latent state of background energy (Chandrasekhar, 1968).
Finally in this section we observe that the wave energy Eq. ( 18) may be cast in the more "physical" or wave energyexchange form, where the quantity (ρ 0 u 2 + V 2 s ρ 2 e /ρ 0 )/2 is the wave energy density and p e u is the actual wave energy flux.In this form it is evident that there is indeed a wave energy exchange with the inhomogeneous stratified background through the term p e u z /H on the right hand side of Eq. ( 24).A similar exchange term of the form ρ 0 u x u y dV 0 /dy is present in the case of a sheared mean flow where the background zonal flow speed V 0 (y) varies latitudinally, and this in turn may give rise to Kelvin-Helmholtz type instability which feeds off the mean flow energy.

JWKB dispersion equation and local stability analysis
It has been shown (McKenzie, 2009) that the diagnostic form of the dispersion equation, a cubic in ω, yields complex conjugate roots for ω, characteristic of a convective instability (Akhiezer et al., 1967), when the dimensionless parameter m = |cosφ|βV /f 2 0 > 1, where (k x ,k y ) = k(cosφ,sinφ).Here the propagation angle φ lies in the second or third quadrant indicating westward propagation.The gravity-inertial wave coalesces with a westward propagating Rossby wave.A similar instability arises in wave hierarchies in which a main wave interacts with a higher order wave of approximately the same speed and instability arises as "a consequence of an unresolvable competition between the two sets of waves" (Whitham, 1974).Here we summarize the results of the instability calculation using a different normalization from (McKenzie, 2009) which enables us to investigate both the equatorial β-plane and the mid-latitudinal case, as well as the possibility that k 2 y may be a given quantity, set, for example, by boundary conditions.This last case, as we shall see, leads to a somewhat modified form of the unstable region of Mach number-latitude space.On normalizing ω to 2 and k to a the diagnostic form ( 14) may be written where: Here we assume that ω N so that V .= N/k z .Recall that these equations also apply to shallow water theory in an ocean of depth h so that V = √ gh and M = a/ √ gh.

The equatorial waveguide
At the equator, θ 0 = 0, f 0 = 0 and f ≈ βy with the result that the wave frequency ω exceeds f close to the equator.However f increases with y so that at some point ω < f may be achieved, or more importantly there may exist a latitude y = y r at which reflection (given by k 2 y = 0 in the structure Eq. ( 10) takes place where Near y = y r the wave structure Eq. ( 10) assumes the classical Airy form characteristic of wave reflection processes.For |y| < y r waves bounce between these two latitudes to form a waveguide whilst propagating zonally.This case may be solved explicity in terms of Hermite polynomials which yield eigenvalues (arising from evanescent requirements as |y| → ∞) given by (Cane and Sarachik, 1976;Moore and Philander, 1977) where n = 0,1,2,..., is a non-negative integer.We note that even in the most favourable case of the fundamental n = 0, Eq. ( 28) does not posess complex conjugate roots for ω and there is no instability.However, this analysis can only claim to be approximately valid if the first reflection point given by lies well inside the Tropics, y 4a, or equivalently that the equatorial Rossby deformation radius is very much less than the Earth's radius.If the equatorial rotational Mach (or Froude) number M 1 the description in terms of Hermite polynomials is indeed a good approximation, and would also apply to atmospheric waves with k z H 1 because the wave speed In this case M is indeed large provided k z H 1 even if a/c is order unity (actually for the Earth ∼ 1.5).Similarly the "baroclinic" mode in a two-layer system of the ocean can yield wave speeds V of order 5 m s −1 once again making a/c large.
It is interesting to note that the dispersion Eq. ( 28), for the fundamental (n = 0) is of precisely the same form as the JWKB Eq. (25) in which ω 2 i is replaced by βV + V 2 k 2 x .It is as if the stability parameter m were set equal to unity which corresponds to the Rossby and Poincaré-Kelvin speeds being equal.Therefore strictly speaking, although the dispersion Eq. ( 28), with n = 0, does not exhibit complex ω, the system is marginally stable.
However if M is not large but is in fact O(1) the solution in terms of Hermite polynomials is not a good approximation.In this case the local JWKB approximation for freely propagating waves located far from boundaries is appropriate.In the zeroth order we may set f = 0 (y a) and with ω 2 i = V 2 k 2 in Eq. (15a) (and its normalization counterpart in Eq. 24) which becomes For westward propagating waves (cosφ < 0) the double root condition shows this equation has complex roots for ω provided Near marginal stability the growth rate γ is approximately and is centred around the frequency Recent work (Maas and Harlander, 2007) on equatorial wave attractors and inertial oscillators in which they find a localized instability, not requiring shear flow driving, is particularly important and possibly relevent to the present work.

The mid-latitudinal case
In this case the instability condition (complex conjugate roots for ω (Eq.25)) becomes Since the right hand side has a minimum value of sin 2 θ 0 occuring at k2 /4M 2 = sin 2 θ 0 , the instability requirement (Eq.35) may be written This is an alternative and equivalent form of the instability condition already given by McKenzie (2009) in which the stability parameter is > 1.This a formal, analytic way of expressing the requirement βV /f 2 > 1 already obtained on intuitive grounds.Hence for a given M the unstable region lies at latitudes less than the critical latitude θ c given by where the cosφ has been absorbed into M.Note that the condition (Eq.36) is not a large β requirement.β is given by the planet's characteristics.The crucial quantity is in fact the wave speed V .We note that in fact the instability condition (36) lends itself to the more easily understood interpretation that instability sets in if the Rossby wave speed exceeds the Poincaré-Kelvin speed.From Eq. ( 37) it follows that waves with M > 3.4, the corresponding θ c lies inside the Tropics, which is therefore better described by the foregoing equatorial approximation.However the fast barotropic wave (V ∼ 200 m s −1 ), for which M is as "low" as M ∼ 2.33, gives θ c = 27 • and therefore a β-plane centred on the tropics may provide a good approximation.In such a case there is a band of wave numbers exhibiting instability with growth rate γ approximately given by (39) Alternatively one may prescribe a real frequency ω and solve Eq. ( 14) for k as a function of ω of the form This shows that if cosθ 0 |cosφ| the expression inside the radical is negative in the frequency band which separates the low frequency Rossby wave from the high frequency inertial-gravity mode.The wavenumber k is then complex and represents (stable) solutions decaying away from the meridional boundary.The "evanescent" condition (Eq.40b), pertaining to the frequency gap between inertial/gravity waves and Rossby waves is in fact just the opposite of the instability condition (Eq.36).

The instability condition and growth rate for prescribed meridional wavenumbers
It is some interest to examine the instability condition (Eq.36) for prescribed (given) values of the N-S wave number k y .This can arise if zonal walls are erected yielding wave-guide values for k y inversely proportional to the width between the walls.There is also the interesting case of evanescent solutions, k 2 y < 0, which arise for values of k x lying outside the Rossby wave normal circle, that is if then k 2 y < 0. The wave structure equation would then yield exponential decay about that latitude.In these cases the instability condition (Eq.36) takes the form where we have used the fact that the right side of Eq. ( 35) exhibits a minumum value at Hence the instability region in the (M,θ 0 ) plane lies beneath the curve given by cosθ In the case k 2 y > 0 the unstable region is confined to latitudes less than For moderate to large values of k y (> 1) the unstable region is therefore confined to a fairly narrow belt around the equator.
On the other hand for evanescent waves, k 2 y < 0, the unstable region is extended in the parameter space.These properties are depicted in Fig. 1.

Summary and discussion
The latitudinal wave structure governing combined inertialgravity-Rossby waves on a β-plane in the Boussinesq approximation admits JWKB wave packet solutions which are valid in the "short wavelength" limit.The corresponding local dispersion equation shows that if the "coupling" between a westward propagating inertial-gravity wave and Rossby wave is sufficiently strong these two modes coalesce and give rise to a convective instability.A necessary and sufficient condition for the unstable mode coupling to take place is that the low frequency Rossby wave speed exceeds the high frequency Poincaré-Kelvin speed.The implication of this phenomenon is outlined in Sect. 4 for both the equatorial and mid-latitudinal cases, as well as for the situation in which the meridional wave number may be assigned prescribed values.
The set of papers by Paldor et al. (2007); Paldor and Sigalov (2008); Leon and Paldor (2009) is pertinent and of special interest since it develops numerical and analytic solutions to the "full-wave" structure equation.Unfortunately Ann.Geophys., 29,[997][998][999][1000][1001][1002][1003]2011 www.ann-geophys.net/29/997/2011/these solutions appear to apply only to cases where the coupling parameter m < 1 and hence the "local" instability condition was not put to a proper test.Moreover zonal wall boundary conditions, although providing the structure equation with a fixed eigenvalue problem, are not appropriate to freely propagating waves.Instead it would be preferable to impose radiation-type conditions or finiteness requirements over the poles as is done for full wave solutions covering the sphere (e.g.Longuet-Higgins, 1968).Indeed recent work (Maas and Harlander, 2007) shows that a local instability, not requiring mean shear flow, may arise in equatorial regions as a result of wave pile-up.
The wave energy-exchange Eq. ( 24) shows that there is an energy exchange through the stratification term on the right hand side of that equation.This provides the channel between horizontal and vertical momentum and access to the enthalpy latent in the background state to drive the instability arising from the coalescence of the two westward propagating modes.The non-linear evolution of this instability may play an important, hitherto unrecognised, role in atmoshereocean dymanics.

Fig. 1 .
Fig. 1.The marginal stability curves (given by Eq. 42) in the (θ 0 ,M) plane for various values of the meridional wave number k 2 y (including an evanescent case k 2 y = −4).The regions beneath the curves are unstable.