2.1. The governing equations

The basis for our investigation are the compressible Euler equations in a rotating coordinate system on the sphere, reproduced below,

(2.1a)$$\begin{gather} \frac{{\rm D}\boldsymbol{v}}{{\rm D}t}+c_p \theta\boldsymbol{\nabla}{\rm \pi}+g\,\boldsymbol{e}_z +2 \boldsymbol{\varOmega}\times\boldsymbol{v}=0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{{\rm D}\theta}{{\rm D}t}=0, \end{gather}$$

(2.1c)$$\begin{gather}\frac{{\rm D}{\rm \pi}}{{\rm D}t}+\frac{R}{c_v} {\rm \pi}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}=0, \end{gather}$$

where the velocity vector $\boldsymbol {v}=u \boldsymbol {e}_x+v \boldsymbol {e}_y+w \boldsymbol {e}_z$ is composed of the zonal, meridional and vertical wind. The variable $R=c_p-c_v$ denotes the specific gas constant for dry air, where $c_p$ and $c_v$ are the heat capacities at constant pressure and constant volume, respectively. Furthermore, $\boldsymbol {e}_j, j\in \{x,y,z\}$ represents the unit vectors, and $g$ is the gravitational acceleration. The Earth's angular velocity vector in the direction of the Earth's rotation axis is $\boldsymbol {\varOmega } = \varOmega _x \boldsymbol {e}_x + \varOmega _y \boldsymbol {e}_y + \varOmega _z \boldsymbol {e}_z$. The Exner pressure ${\rm \pi}$ and potential temperature $\theta$ are defined by the canonical thermodynamic quantities, temperature $T$ and pressure $p$. These quantities are related by the equations of state

(2.2a)$$\begin{gather} {\rm \pi}=(\,p/p_0)^\kappa, \end{gather}$$

(2.2b)$$\begin{gather}\theta=T/{\rm \pi}, \end{gather}$$

with $p_0$ a reference pressure at $z=z_0$ and $\kappa =R/c_p=2/7$ for two-atomic gases. By means of the ideal gas law, the thermodynamical variables are linked to the density $\rho$ via

(2.3)\begin{equation} \rho=\frac{p_0}{R}\frac{{\rm \pi}^{(1-\kappa)/\kappa}}{\theta}. \end{equation}

The material derivative is given as

(2.4)\begin{equation} \frac{{\rm D}}{{\rm D}t} =\frac{\partial}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}, \end{equation}

where $\boldsymbol {\nabla } = \boldsymbol {e}_x\,\partial / \partial x + \boldsymbol {e}_y\,\partial / \partial y + \boldsymbol {e}_z\,\partial / \partial z$ denotes the nabla operator. The governing equations introduced in (2.1) are based on mass, momentum and energy conservation.

2.2. The non-traditional setting

The effect of the non-traditional setting is most obvious close to the equator, and so we assume $f$–$F$-planes at the equator (Thuburn, Wood & Staniforth Reference Thuburn, Wood and Staniforth2002), i.e. $f=0$ and $F=2\varOmega$ with $\varOmega =|\boldsymbol {\varOmega }|$. More details on the $f$–$F$-planes notation may be found in Appendix A. The $\beta$-plane effect is additionally discussed by means of scale analysis in § 2.8. From here on, we consider only the flow in the $x$–$z$-plane, i.e. we restrict the analysis to meridionally homogeneous fields. The horizontal domain of interest is $x\in [0,x_{MAX}]$ with $z\in [z_0,z_{TOA}]$ for the vertical extent and $z_0$ the reference height, where TOA stands for a fictitious demarcation representing the ‘top of the atmosphere’, and a clearer elaboration on the boundary conditions will be provided in § 2.4. This problem set-up is illustrated in figure 1.

Figure 1. Sketch of the Coriolis acceleration at the equator.

Writing out the compact form of (2.1) explicitly yields the following:

(2.5a)$$\begin{gather} \frac{\partial u}{\partial t} +u\frac{\partial u}{\partial x} +w\frac{\partial u}{\partial z} +c_p\theta\frac{\partial{\rm \pi}}{\partial x}+Fw=0, \end{gather}$$

(2.5b)$$\begin{gather}\frac{\partial w}{\partial t} +u\frac{\partial w}{\partial x} +w\frac{\partial w}{\partial z} +c_p\theta\frac{\partial{\rm \pi}}{\partial z} -Fu+g=0, \end{gather}$$

(2.5c)$$\begin{gather}\frac{\partial \theta}{\partial t} +u\frac{\partial \theta}{\partial x} +w\frac{\partial \theta}{\partial z}=0, \end{gather}$$

(2.5d)$$\begin{gather}\frac{\partial {\rm \pi}}{\partial t} +u\frac{\partial {\rm \pi}}{\partial x} +w\frac{\partial {\rm \pi}}{\partial z} +\frac{R}{c_v}{\rm \pi}\left( \frac{\partial u}{\partial x} +\frac{\partial w}{\partial z} \right)=0. \end{gather}$$

2.3. The hydrostatic background atmosphere

A stationary solution to the governing equations (2.5) is the hydrostatic atmosphere at rest, i.e.

(2.6a)$$\begin{gather} u(x,z,t)=0, \end{gather}$$

(2.6b)$$\begin{gather}w(x,z,t)=0, \end{gather}$$

(2.6c)$$\begin{gather}\theta(x,z,t)=\bar\theta(z), \end{gather}$$

(2.6d)$$\begin{gather}{\rm \pi}(x,z,t)=\bar{\rm \pi}(z). \end{gather}$$

Substituting (2.6) into (2.1), all governing equations but the vertical velocity equation (2.5b) vanish. The remaining equation reduces to the hydrostatic equation

(2.7)\begin{equation} c_p\,\bar\theta\frac{\mathrm{d}\bar{\rm \pi}}{\mathrm{d}z}=-g. \end{equation}

Given a temperature profile $\bar {T}=\bar {T}(z)$, we can integrate (2.7) using the definitions in (2.2) to obtain the hydrostatic background variables, and they may be written as follows:

(2.8)$$\begin{gather} \bar{\rm \pi}(z)=\exp\left[-\int_0^{z}\frac{g}{c_p}\frac{1}{\bar{T}(\zeta)}\ \mathrm{d}\zeta\right], \end{gather}$$

(2.9)$$\begin{gather}\bar{\theta}(z)=T_0\exp\left[\int_0^{z}\frac{1}{\bar{T}(\zeta)}\left(\frac{g}{c_p} +\frac{\mathrm{d}\bar{T}}{\mathrm{d}\zeta}\right)\ \mathrm{d}\zeta\right], \end{gather}$$

(2.10)$$\begin{gather}\bar\rho(z)=\frac{p_0}{RT_0}\exp\left[-\int_0^{z}\frac{1}{\bar{T}(\zeta)}\left(\frac{g}{R} +\frac{\mathrm{d}\bar{T}}{\mathrm{d}\zeta}\right)\ \mathrm{d}\zeta\right]. \end{gather}$$

Subsequently, we consider perturbations from the basic, hydrostatic state by the following ansatz:

(2.11a)$$\begin{gather} u(x,z,t)=u'(x,z,t), \end{gather}$$

(2.11b)$$\begin{gather}w(x,z,t)=w'(x,z,t), \end{gather}$$

(2.11c)$$\begin{gather}\theta(x,z,t)=\bar\theta(z)+\theta'(x,z,t), \end{gather}$$

(2.11d)$$\begin{gather}{\rm \pi}(x,z,t)=\bar{\rm \pi}(z)+{\rm \pi}'(x,z,t). \end{gather}$$

Assuming that the perturbations are infinitesimally small, the governing equations may be linearised around the basic hydrostatic state. Before stating the linearised system, we introduce a transformation of the primed variables concatenated into a vector $\boldsymbol {\chi }=(\chi _u, \chi _w, \chi _\theta, \chi _{\rm \pi} )^\mathrm {T}$ that simplifies the upcoming derivations (cf. Achatz, Klein & Senf Reference Achatz, Klein and Senf2010, equation (2.16))

(2.12)\begin{equation} \boldsymbol{\chi}=\bar{\rho}^{1/2}\left(u', w', \frac{g}{N}\,\frac{{\theta^\prime}}{\bar{\theta}}, \frac{c_p}{C} \bar{\theta}{\rm \pi}'\right)^\mathrm{T}. \end{equation}

Here, the Brunt–Väisälä frequency and the speed of sound are given by

(2.13)$$\begin{gather} N(z)=\sqrt{\frac{g}{\bar{\theta}(z)}\frac{\mathrm{d}\bar{\theta}}{\mathrm{d}z}} =\sqrt{\frac{g}{\bar{T}(z)}\left(\frac{g}{c_p} +\frac{\mathrm{d}\bar{T}}{\mathrm{d}z}\right)}, \end{gather}$$

(2.14)$$\begin{gather}C(z)=\sqrt{\gamma R\bar{T}(z)},\quad \gamma=\frac{c_p}{c_v}, \end{gather}$$

respectively. Inserting (2.12) into (2.11) and the result into (2.5), we may write the linearised system in terms of the transformed perturbation vector components as

(2.15a)$$\begin{gather} \frac{\partial\chi_u}{\partial t}+\frac{\partial}{\partial x}\left(C\chi_{\rm \pi}\right)+F\chi_w=0, \end{gather}$$

(2.15b)$$\begin{gather}\frac{\partial\chi_w}{\partial t}+\left(\frac{\partial}{\partial z}+\varGamma\right)\left(C\chi_{\rm \pi}\right)-N\chi_\theta-F\chi_u=0, \end{gather}$$

(2.15c)$$\begin{gather}\frac{\partial\chi_\theta}{\partial t}+N\chi_w=0, \end{gather}$$

(2.15d)$$\begin{gather}\frac{\partial\chi_{\rm \pi}}{\partial t}+C\left[\frac{\partial\chi_u}{\partial x}+\left(\frac{\partial}{\partial z}-\varGamma\right)\chi_w\right]=0, \end{gather}$$

where we introduced the quantity $\varGamma$, which has dimension of inverse height, as

(2.16)\begin{align} \varGamma(z)&=-\frac{1}{2\bar{\rho}(z)}\frac{\mathrm{d}\bar{\rho}}{\mathrm{d}z} -\frac{1}{\bar{\theta}(z)}\frac{\mathrm{d}\bar{\theta}}{\mathrm{d}z} \nonumber\\ &=\frac{1}{2\bar{T}(z)}\left[\left(\frac{c_v}{R}-1\right)\frac{g}{c_p} -\frac{\mathrm{d}\bar{T}}{\mathrm{d}z}\right] , \end{align}

in accordance with Durran (Reference Durran1989). Notice that by the choice of the variable transformation,

(2.17)\begin{equation} e=\tfrac{1}{2} \left( \chi_u^2+\chi_w^2+\chi_\theta^2+\chi_{\rm \pi}^2 \right), \end{equation}

is the total energy density of the perturbation which evolves in time as

(2.18)\begin{equation} \frac{\partial e}{\partial t}+\frac{\partial}{\partial x}\left(C\chi_u\chi_{\rm \pi}\right) +\frac{\partial}{\partial z}\left(C\chi_w\chi_{\rm \pi}\right)=0. \end{equation}

Hence, Gauss’ theorem ensures that the total energy is locally conserved. The quantity $\chi _u^2+\chi _w^2$ represents the kinetic energy density of the perturbation, and the quantity $\chi _\theta ^2 + \chi _{\rm \pi} ^2$ represents the sum of the potential and internal energy density. The perspective of energy conservation lends an additional physical interpretation to the newly introduced variable $\varGamma$. Considering (2.15), we observe that the quantity $C\varGamma$ represents an energy exchange rate for the conversion between kinetic energy due to vertical motion and internal energy which is associated with the compressibility of the gas.

2.4. Linear stability analysis

Proceeding with our theoretical developments, notice that the solution of the linearised equations (2.15) is of the form

(2.19)\begin{equation} \boldsymbol{\chi}(x,z,t)=\boldsymbol{\psi}(x,z)\exp(-{\rm i}\lambda t),\quad \lambda\in\mathbb{C}. \end{equation}

This solution ansatz yields an eigenvalue problem with eigenvalues ${\rm i}\lambda$, i.e.

(2.20)\begin{equation} (\boldsymbol{\mathsf{T}}-{\rm i}\lambda)\boldsymbol{\psi}=0, \end{equation}

with the differential operator

(2.21) \begin{equation} \boldsymbol{\mathsf{T}}= \begin{pmatrix} 0 & F & 0 & \partial/\partial x(C{\times})\\ -F & 0 & -N & \partial/\partial z(C{\times})+\varGamma C\\ 0 & N & 0 & 0\\ C\partial/\partial x & C\partial/\partial z-C\varGamma & 0 & 0 \end{pmatrix}. \end{equation}

Note that as $\boldsymbol {\psi }$ is complex valued, the solution of (2.20) is a linear combination of the eigenfunctions. This ensures that the solutions in (2.19) are eventually real fields.

Now, let us define an inner product as

(2.22)\begin{equation} \langle\boldsymbol{\psi}|\boldsymbol{\phi}\rangle =\int_{z_0}^{z_{TOA}}\int_0^{x_{MAX}} \boldsymbol{\psi}^{H}\boldsymbol{\phi} \ \mathrm{d}\kern0.06em x\ \mathrm{d}z \end{equation}

in terms of some well-behaved vectors of prognostic variables $\boldsymbol {\psi }=(\psi _u,\psi _w,\psi _\theta,\psi _{\rm \pi} )^\mathrm {T}$ and $\boldsymbol {\phi }=(\phi _u,\phi _w,\phi _\theta,\phi _{\rm \pi} )^\mathrm {T}$ that satisfy the yet to be defined boundary conditions. Furthermore, $\mathrm {H}$ denotes the Hermitian transpose, i.e. the conjugate transpose. This inner product induces a norm that is given as

(2.23)\begin{equation} \|\boldsymbol{\psi}\|^2 =\langle\boldsymbol{\psi}|\boldsymbol{\psi}\rangle =\int_{z_0}^{z_{TOA}}\int_0^{x_{MAX}} \boldsymbol{\psi}^{H}\boldsymbol{\psi} \ \mathrm{d}\kern0.06em x\ \mathrm{d}z =\int_{z_0}^{z_{TOA}}\int_0^{x_{MAX}} 2e \ \mathrm{d}\kern0.06em x\ \mathrm{d}z = 2E, \end{equation}

which corresponds to twice the global energy $E$, and therefore this choice of the inner product is motivated by a physical argument.

At this point, an open question that we still have to consider is the choice of realistic boundary conditions. To answer this question, let us consider the properties of the operator $\boldsymbol{\mathsf{T}}$. If we assume for simplicity periodic boundary conditions in the $x$-direction, such that $\boldsymbol {\psi }(0,z)=\boldsymbol {\psi }({x_{MAX}},z)$ (or vanishing fields at $x\rightarrow \pm \infty$), we then observe that

(2.24)\begin{equation} \langle\boldsymbol{\mathsf{T}}\boldsymbol{\psi}|\boldsymbol{\phi}\rangle =-\langle\boldsymbol{\psi}|\boldsymbol{\mathsf{T}}\boldsymbol{\phi}\rangle +\int_0^{x_{MAX}}\left[ C\psi_{\rm \pi}^\ast\phi_w+C\psi_w^\ast\phi_{\rm \pi} \right]_{z_0}^{z_{TOA}} \ \mathrm{d}\kern0.06em x, \end{equation}

for all $\boldsymbol {\psi }$ and $\boldsymbol {\phi }$. Here, $\ast$ represents the complex conjugate. This observation implies that the operator $\boldsymbol{\mathsf{T}}$ is skew-Hermitian if the following condition holds:

(2.25)\begin{equation} \left[C\psi_{\rm \pi}^\ast\phi_w+C\psi_w^\ast\phi_{\rm \pi} \right]_{z_0}^{z_{TOA}}=0. \end{equation}

Skew-Hermitian operators exhibit purely imaginary eigenvalues ${\rm i}\lambda$ (or rather an imaginary spectrum to be precise), so $\lambda \in \mathbb {R}$. Therefore, the perturbations in (2.15) are bounded in time and ultimately stable if the condition in (2.25) holds.

Whether or not condition (2.25) holds now depends on the choice of the top and bottom boundaries. If we were to assume a free-slip (no-normal-flow) boundary, i.e. $w(x,0,t)=z_0$, for the bottom with no orography, then this boundary condition implies that

(2.26)\begin{equation} \psi_w,\phi_w=0,\quad\text{at }z=z_0. \end{equation}

If we were to also restrict the fields to vanish at the upper boundary, then the condition (2.25) holds, and all perturbations would be stable. A similar result generalised to the equatorial $\beta$-plane was also obtained by Fruman & Shepherd (Reference Fruman and Shepherd2008) who considered zonally symmetric flows. In fact, all trivial boundary conditions would result in stable perturbations. We can illuminate this statement by integrating (2.18) over the spatial domain, which leaves us with the rate of change of the global energy $E$ given in terms of the vertical fluxes of total energy

(2.27)\begin{equation} \frac{\mathrm{d}E}{\mathrm{d}t}=\int_0^{x_{MAX}}\left.C\chi_w\chi_{\rm \pi}\right|_{z=z_0}\ \mathrm{d}\kern0.06em x-\int_0^{x_{MAX}}\left.C\chi_w\chi_{\rm \pi}\right|_{z=z_{TOA}}\ \mathrm{d}x, \end{equation}

at $z_0$ and the top of the atmosphere, respectively. Therefore, for an instability to grow in energy, there must be an energy flux through the boundaries, but for trivial boundary conditions, the flux would vanish. As energy cannot come from space, the only reasonable energy source must be at $z_0$.

2.5. The isothermal background atmosphere

Let us now simplify our model even further by assuming an isothermal background state, $\partial \bar {T}/\partial z=0$. As a consequence, all coefficients in $\boldsymbol{\mathsf{T}}$ become constant and, taking into account the horizontal boundary conditions, the perturbations must now be of the form

(2.28)\begin{equation} \boldsymbol{\psi}(x,z)=\tilde{\boldsymbol{\psi}}\exp({\rm i}kx+\mu z),\quad k\in\mathbb{R}, \mu\in\mathbb{C}, \end{equation}

where $\tilde {\boldsymbol {\psi }}$ is a vector with constant components. The operator $\boldsymbol{\mathsf{T}}$ may be written as a matrix, that is,

(2.29)\begin{equation} \tilde{\boldsymbol{\mathsf{T}}}= \begin{pmatrix} 0 & F & 0 & {\rm i}Ck\\ -F & 0 & -N & C(\mu+\varGamma)\\ 0 & N & 0 & 0\\ {\rm i}Ck & C(\mu-\varGamma) & 0 & 0 \end{pmatrix}. \end{equation}

Dividing (2.29) by $N$, the system may be non-dimensionalised by defining new dimensionless variables that are given as follows:

(2.30a–d)\begin{equation} K=\frac{C}{N}\,k,\quad M=\frac{C}{N}\,\mu,\quad \varLambda=\frac{1}{N}\,\lambda,\quad \varepsilon=\frac{F}{N}. \end{equation}

Writing the eigenvalue problem of (2.20) in terms of the dimensionless variables and incorporating the horizontal periodic boundary conditions leads to

(2.31a,b)\begin{equation} \left(\boldsymbol{\mathsf{S}}-{\rm i}\varLambda\right)\,\tilde{\boldsymbol{\psi}}=0,\quad \boldsymbol{\mathsf{S}}(K,M;\varepsilon)= \begin{pmatrix} 0 & \varepsilon & 0 & {\rm i}K\\ -\varepsilon & 0 & -1 & M+G\\ 0 & 1 & 0 & 0\\ {\rm i}K & M-G & 0 & 0 \end{pmatrix}, \end{equation}

where we have

(2.32)\begin{equation} G=\frac{C\varGamma}{N}=\frac{1-\gamma/2}{\sqrt{\gamma-1}}=\sqrt{\frac{9}{40}}. \end{equation}

The characteristic polynomial of $\boldsymbol{\mathsf{S}}$ is obtained, and it is

(2.33)\begin{equation} \mathcal{P}_{\boldsymbol{\mathsf{S}}(K,M;\varepsilon)}({\rm i}\varLambda)=\varLambda^4 -\left(1+\varepsilon^2+G^2+K^2-M^2\right)\varLambda^2 +2\varepsilon G K\varLambda +K^2. \end{equation}

As this polynomial is a depressed quartic function in $\varLambda$, its roots may be expressed by radicals. However, these expressions would be tedious and provide no further elucidation. Instead, we note that $\varepsilon =O(10^{-2})$ is typically a very small number, and so it may be sensible to employ perturbation theory for the subsequent developments. If we assume that $M,\,K=O(1)$ as $\varepsilon \rightarrow 0$, then we have, to leading order, the following roots for the characteristic polynomial:

(2.34)$$\begin{gather} {\varLambda^{(0)}}^2=\frac{1}{2}L^2\left(1\pm\sqrt{1-4\frac{K^2}{L^4}}\right), \end{gather}$$

(2.35)$$\begin{gather}L^2=K^2-M^2+1+G^2=K^2-M^2+\frac{49}{40}, \end{gather}$$

which is essentially the dispersion relation for acoustic-gravity waves.

2.6. The leading-order eigenvector

Computing the leading-order eigenvector of $\boldsymbol{\mathsf{S}}(K,M;0)$ gives us

(2.36)\begin{equation} \tilde{\boldsymbol{\psi}}^{(0)}= \left(1, \frac{{\rm i}}{K}\frac{{\varLambda^{(0)}}^2-K^2}{M-G}, \frac{1}{{\varLambda^{(0)}} K}\frac{{\varLambda^{(0)}}^2-K^2}{M-G}, \frac{{\varLambda^{(0)}}}{K}\right)^\mathrm{T}. \end{equation}

A particularly interesting situation arises when we assume a free-slip boundary condition at the ground to the leading order, such that $\psi _w^{(0)}=0$ at $z=z_0$, because the leading-order total energy flux then also vanishes according to (2.27). If we were able to show that an instability occurs in the next-order correction, then the energy flux would be, at most, next-to-leading order. Or, in other words, a minuscule energy flux through the bottom boundary that would otherwise be ignored might lead to an exponentially growing perturbation. We consider this situation worthwhile to investigate in more detail. The leading-order eigenvector (2.36) fulfils the free-slip boundary condition non-trivially if

(2.37)\begin{equation} {\varLambda^{(0)}}^2=K^2. \end{equation}

Inserting (2.37) into (2.33), we can solve the leading-order characteristic polynomial to obtain

(2.38)\begin{equation} M=-G, \end{equation}

which ensures a physical solution that decays as $z\rightarrow \infty$. Hence, the leading-order eigenvector becomes

(2.39)\begin{equation} \tilde{\boldsymbol{\psi}}^{(0)}= \left(1, 0, 0, \pm 1\right)^\mathrm{T}, \end{equation}

and this solution describes a Lamb wave (Vallis Reference Vallis2017). Since the characteristic polynomial is a quartic, we expect to find four roots. Yet for the free-slip lower boundary condition we imposed, only two roots were obtained in (2.37). We note that the leading-order polynomial has two more roots for $M=-G$ that do not satisfy the leading-order boundary conditions for all $K$, and these are

(2.40)\begin{equation} {\varLambda^{(0)}}^2=1. \end{equation}

These roots would describe Brunt waves (Walterscheid Reference Walterscheid2003).

2.7. The instability growth rate

Figure 2 depicts all four leading-order roots (grey lines) for varying horizontal wavenumber $K$. Notice that, to leading order, we obtain only real-valued roots, and these translate to purely imaginary and hence stable eigenvalues (spectrum). However, the numerical solution for $\varepsilon \neq 0$ reveals imaginary parts for $\varLambda$ (the dashed blue lines) and, consequentially, the perturbations might be unstable in these regions. Referring to figure 2, we furthermore observe that the potentially unstable eigenvalues appear around $K=kC/N=\pm 1$ where the dispersion curves of Lamb and Brunt waves collide. This point is distinct as the leading-order operator degenerates; its eigenvalues have algebraic multiplicity of two instead of one at this point. Perez et al. (Reference Perez, Delplace and Venaille2022) identified and investigated almost the same points of degeneracy. Remarkably, they found a novel unidirectional mode with this approach. The main difference to this study is induced by the dissimilar boundary conditions investigated here.

Figure 2. Leading-order roots of the characteristic polynomial $\mathcal {P}_{\boldsymbol{\mathsf{S}}(K,G;0)}({\rm i}\varLambda ^{(0)})$ (grey lines, see (2.33)). Real part of the numerically computed roots of $\mathcal {P}_{\boldsymbol{\mathsf{S}}(K,G;0.1)}({\rm i}\varLambda )$ (solid blue line) and its imaginary part (dashed blue line) for $\varepsilon =0.1$. We chose an unrealistically large $\varepsilon$ in order to exemplify its effect.

We now determine the maximum imaginary part of $\varLambda$, and this would give us the instability growth rate of the perturbation. To do so entails solving for an asymptotic solution at the point of degeneracy, $K=\pm 1$. However, perturbation theory for non-Hermitian operators that have leading-order eigenvalues with algebraic multiplicity greater than one is in fact cumbersome (Kato Reference Kato1982). Its rigorous mathematical treatment is beyond the scope of this paper as standard techniques from, say, quantum mechanics are inapplicable. That being the case, we will nevertheless provide an asymptotic analysis of the next-order correction to the eigenvalue, while refraining from the more involved computations for the eigenvectors. For this purpose, we expand $\varLambda$ in terms of powers of $\varepsilon$, i.e.

(2.41)\begin{equation} \varLambda=\varLambda^{(0)}+\varepsilon^{1/2}\varLambda^{(1/2)}+o(\varepsilon^{1/2}). \end{equation}

The choice for the exponent of $1/2$ for the next-to-leading-order term is due to the multiplicity of the eigenvalue. We refer the interested reader to a complete reasoning on this matter in Lidskii (Reference Lidskii1966). Also, Kato (Reference Kato1982) contains a proof on the existence and form of the expansion. Substituting the ansatz together with the leading-order result in the characteristic polynomial (2.33) and collecting terms of the same power in $\varepsilon$, we find that the $O(\varepsilon ^{1/2})$ terms cancel each other. At $O(\varepsilon )$, we obtain

(2.42)\begin{equation} \varLambda^{(1/2)}=\pm {\rm i}\sqrt{\frac{G}{2}}. \end{equation}

Rewriting (2.42) together with the leading-order result in a dimensional form yields

(2.43)\begin{equation} \lambda\approx\pm N\pm {\rm i}\sqrt{\frac{F}{2}C\varGamma} ,\end{equation}

and hence the instability growth rate close to the equator is given as $\sqrt {\varOmega C\varGamma }$.

2.8. Influence of the $\beta$-plane effect

We have assumed the $f$–$F$-planes up to this point of our theoretical investigations. To conclude our theoretical investigations, we now discuss a possible influence by the $\beta$-plane effect which may be estimated by means of scale analysis.

Given an isothermal atmosphere such that $\bar {T}(z) \equiv T_0 = 300$ K, $\gamma = 1.4$, $g \approx 9.81$ m s$^{-1}$, $\varOmega \equiv \varOmega _y = 7.292 \times 10^{-5}$ s$^{-1}$ and $R=287.4$ J kg$^{-1}$ K$^{-1}$, we obtain $N\approx 0.02$ s$^{-1}$ and $C\approx 347$ m s$^{-1}$. With $k=N/C$, this leads to a wavelength of the unstable mode of approximately 120 km, and a theoretical instability growth rate of

(2.44)\begin{equation} \sqrt{\varOmega C \varGamma} \approx 7.9 \times 10^{-4}\,\text{s}^{-1}. \end{equation}

This instability growth rate corresponds to a doubling period of the perturbation's amplitude approximately every 15 min.

By some Wentzel–Kramers–Brillouin-type argument, we may assume that the perturbation is meridionally homogeneous if it remains relatively constant within a latitude range of 10 wavelengths. Let us now compute the rate of the $\beta$-plane effect due to the term $f^\prime =\beta y$ at the northern edge of a meridional domain of size $10\times 120\ {\rm km}=1200\ {\rm km}$. Taking $\beta =2\varOmega /a$ at the equator where $a$ represents the Earth's radius, this amounts to $f'\approx 1\times 10^{-5}$ s$^{-1}$. As the $\beta$-plane effect is linear in the governing equations, we may compare it with the theoretical growth rate obtained from linear stability analysis and from (2.44), we observe that the theoretical growth rate is approximately 80 times greater than $f^\prime$. Hence, while damping of the perturbation might be possible due to the $\beta$-plane effect, the damping rate would be almost two orders of magnitude smaller than the instability growth rate. We therefore anticipate that the inclusion of the $\beta$-plane approximation will have only a small effect on the theory and results presented in this paper.

Summarising this section, we found a linearly unstable meridionally homogeneous mode of the hydrostatic atmosphere at rest close to the equator by means of asymptotic perturbation theory. The mode is, to leading order, identical to Lamb waves which are a known horizontally propagating external mode of the atmosphere. Its phase velocity is the speed of sound. To leading order, a free-slip boundary condition was assumed which prohibits any vertical energy flux. We derived a higher-order correction to the eigenvalues of the linearised system and provided justifications that an analytical determination of the corresponding correction to the eigenvector is beyond the scope of this paper. Yet, in the next section we do obtain numerical approximations to the eigenmode structure and we use them to initiate simulations based on a full nonlinear compressible solver (see the discussion of (3.1) and (3.2) below).

The validity of our asymptotic theory is therefore tested against the numerical solution of the nonlinear governing equations in the next section.

3.1. Numerical model

Numerical corroboration of the Lamb-wave-like unstable mode that we derived in the previous section necessitates a model that is capable of solving the compressible non-hydrostatic Euler equations. The second-order accurate and semi-implicit numerical model by Benacchio & Klein (Reference Benacchio and Klein2019) is well tested and is especially suitable to study flow instabilities numerically due to its conservation characteristics and numerical stability properties. A brief summary of the model is given in Appendix B. Details on the extension of the numerical model to support the non-traditional setting are documented in Appendix C.

3.2. The initial condition

The initial condition for the numerical experiments is obtained by solving the eigenvalue problem (2.20) with the isothermal background assumption

(3.1)\begin{equation} (\tilde{\boldsymbol{\mathsf{T}}} - {\rm i} {\lambda})\tilde{\boldsymbol{\psi}} = 0. \end{equation}

This is done numerically, utilising the LAPACK library's eigenvalue problem solver. The transformed perturbation variables $\boldsymbol {\chi }$ are then obtained from the solution $\tilde {\boldsymbol {\psi }}=(\tilde {\psi }_u,\tilde {\psi }_w,\tilde {\psi }_\theta,\tilde {\psi }_{\rm \pi} )^\mathrm {T}$ by applying (2.19) and (2.28). Specifically, the initial perturbation fields are given as follows:

(3.2a)$$\begin{gather} u^\prime(x,z,0) = A \tilde{\psi}_u \sqrt{\frac{\rho_0}{\bar{\rho}(z)}} \exp({-\varGamma z}) \cos\left( \frac{N}{C} x \right), \end{gather}$$

(3.2b)$$\begin{gather}w^\prime(x,z,0) = A \tilde{\psi}_w \sqrt{\frac{\rho_0}{\bar{\rho}(z)}} \exp({-\varGamma z}) \cos\left( \frac{N}{C} x \right), \end{gather}$$

(3.2c)$$\begin{gather}\theta^\prime(x,z,0) = A \tilde{\psi}_\theta \frac{N}{g} \sqrt{\frac{\rho_0}{\bar{\rho}(z)}} \bar{\theta} \exp({-\varGamma z}) \cos\left( \frac{N}{C} x \right), \end{gather}$$

(3.2d)$$\begin{gather}{\rm \pi}^\prime(x,z,0) = A \tilde{\psi}_{\rm \pi} \frac{C}{c_p} \sqrt{\frac{\rho_0}{\bar{\rho}(z)}} \bar{\theta}^{-1}(z)\exp({-\varGamma z}) \cos\left(\frac{N}{C} x \right), \end{gather}$$

where $A = 10^{-1}$ m s$^{-1}$ is a prescribed small amplitude. Furthermore, the background density is given as

(3.3)\begin{equation} \bar{\rho}(z) = \rho_0 \exp(-z / H_{\rho}), \end{equation}

with $H_{\rho } = RT_0 / g$ representing the density scale height. The expression for $\varGamma$ in (2.16) simplifies to

(3.4)\begin{equation} \varGamma = \frac{1}{H_{\rho}} \left( \frac{1}{\gamma} - \frac{1}{2} \right), \end{equation}

under the assumption of an isothermal atmosphere, and the background potential temperature is

(3.5)\begin{equation} \bar{\theta}(z) = \theta_0 \exp(z/H_{\theta}), \end{equation}

where $H_{\theta } = H_{\rho } / \kappa$ denotes the potential temperature scale height. The quantities $\rho _0$, $T_0$ and $\theta _0$ are reference density, temperature and potential temperature at $z=z_0$, respectively. The initial density field is obtained via the equation of state (2.3).

The horizontal extent of the domain is chosen such that it is equivalent to four wavelengths, i.e. ${x_{MAX}} = 2{\rm \pi} j C / N$ with $j=4$, and this corresponds to $x \in [-244.5\ \mathrm {km},244.5\ \mathrm {km}]$. The full vertical extent of the domain corresponds to $z\in [0.0,80.0\ \mathrm {km}]$.

3.3. The total energy flux

From the initial condition, we additionally computed the vertical flux of total energy as given in (2.27). Figure 3 shows the profile of the horizontal mean of the vertical energy flux. The positive vertical energy flux indicates an upward flow of energy and decays quickly for increasing $z$. In line with physical expectations and as discussed in § 2.4, energy enters the domain through the bottom boundary but cannot escape into space, resulting in an accumulation of energy in the domain and generating a growing instability. Notice that the initial flux at the bottom boundary is very small, as expected from the asymptotic theory. This initial flux is also small compared with typical sensible or latent heat fluxes in the atmosphere especially considering the initial amplitude of 0.1 m s$^{-1}$, which presents a recognisable breeze. Since our simulation domain is finite and the flux at the top is non-zero, an absorbing non-reflecting boundary condition at the TOA is necessary.

Figure 3. Initial horizontally averaged vertical flux of the total energy.

3.4. Numerical representation of the boundary conditions

To emulate a free surface at the top of the boundary, a wave-absorbing layer is applied to the top 20 km to approximate the non-reflecting boundary condition (Durran Reference Durran2010). Specifically, the Rayleigh damping by Durran & Klemp (Reference Durran and Klemp1983), as reproduced below, has been used. The following terms are added to the right-hand side of the governing equations in (2.1):

(3.6a)$$\begin{gather} R_u = \tau(z) u^\prime, \end{gather}$$

(3.6b)$$\begin{gather}R_w = \tau(z) w^\prime, \end{gather}$$

(3.6c)$$\begin{gather}R_\theta = \tau(z) \theta^\prime, \end{gather}$$

(3.6d)$$\begin{gather}R_{\rm \pi} = \tau(z) {\rm \pi}^{\prime}, \end{gather}$$

with

(3.7)\begin{equation} \tau(z) = \begin{cases} 0, & \text{for } z\leq z_D, \\ -\dfrac{\alpha}{2} \left[ 1 - \cos \left(\dfrac{z-z_D}{z_{TOA} - z_D} {\rm \pi}\right) \right], & \text{for}\ 0 \leq \dfrac{z-z_D}{z_{TOA} - z_D} \leq \dfrac{1}{2},\\ -\dfrac{\alpha}{2} \left[ 1 + \left( \dfrac{z-z_D}{z_{TOA}-z_D} - \dfrac{1}{2} \right) {\rm \pi}\right], & \text{for}\ \dfrac{1}{2} \leq \dfrac{z-z_D}{z_{TOA} - z_D} \leq 1, \end{cases} \end{equation}

where $z_{TOA}=80.0$ km, $z_D=60.0$ km, and $\alpha =0.5$. An inverse of the Rayleigh damping is applied to the bottom 3 km of the domain. Instead of relaxing the flow to the ambient state, this bottom forcing relaxes fully to the perturbation variables obtained from the numerical solution of the eigenvalue problem in (3.2), which will be referred to as the ‘semi-analytical eigenmode forcing’ (SA) below. Such a bottom forcing mimics an atmospheric boundary layer that allows a vertical energy transport into the free atmosphere.

Both the top and bottom boundary conditions are approximations of the physical boundaries of the atmosphere, and these may be sources of error. See § 2.4 for more details on these choices of the top and bottom boundaries.

Furthermore, the domain is periodic in the horizontal direction. The simulations in this section are conducted on a grid with $(N_x \times N_z)= (301 \times 120)$ grid points.

For the initial perturbation in (3.2), the prognostic variables in the form of (2.12) are depicted in figure 4 for the Lamb wave and figure 5 for the unstable mode. Note that the colour scale employed in these and the subsequent contour plot figures are nonlinear.

Figure 4. Initial perturbation: distribution of the prognostic variables with the transformation given in (2.12) for the stable Lamb wave case. (a) Transformed horizontal velocity perturbation $\chi _u$, (b) transformed vertical velocity perturbation $\chi _w$, (c) transformed potential temperature perturbation $\chi _\theta$ and (d) transformed Exner pressure perturbation $\chi _{\rm \pi}$. The contours have units of kg$^{1/2}$ m$^{-1/2}$ s$^{-1}$, i.e. the units of the square root of energy density.

Figure 5. Output of the initial transformed prognostic variable fields for the Lamb wave-like instability depicting the initial non-traditional Coriolis effect, i.e. $\varOmega _y = 7.292 \times 10^{-5}$ s$^{-1}$ at $t=0$ s, corresponding to the initial conditions of the LWLI-SA and LWLI-SO cases. Contour units are as in figure 4.

Two further points must be addressed before we present our results in the subsequent subsection. In Appendix D, we present a careful numerical treatment of the thermodynamical quantities that ensures a well-balanced Lamb wave solution in a stably stratified atmosphere. With this development, the linear solution of the Lamb wave in the traditional setting maintains a vertical velocity field up to machine accuracy for 10 000 time steps. In Appendix E, we investigate the energy conservation of the stable Lamb wave solutions with the numerical scheme. In particular, we demonstrate that energy is conserved up to ${\sim }0.4\,\%$ over 36 000 s for a run with Lamb wave initial condition in the traditional setting, i.e. with $F=0$ s$^{-1}$ in (3.1), and top and bottom rigid wall boundaries. This is in line with the theory developed in § 2. Furthermore, energy is similarly well – although not exactly – conserved for a run that is initialised with a stable Lamb wave but features non-zero Coriolis strength.

3.5. Instability growth rate under the non-traditional setting

Three experimental configurations are investigated below: a run with the Lamb wave-like unstable eigenmode as the initial condition, with Coriolis acceleration in the non-traditional setting, and the semi-analytical eigenmode bottom forcing (LWLI-SA); a run with the Lamb wave-like unstable eigenmode as the initial condition, with Coriolis acceleration in the non-traditional setting, and sub-optimal bottom forcing (LWLI-SO–an explanation on the sub-optimal forcing is given below); and a stable Lamb wave as the initial condition without Coriolis acceleration and bottom forcing (LW). The time-step size $\Delta t=10$ s is restricted by an advective Courant–Friedrichs–Lewy number of ${\rm CFL}_{adv}=0.9$.

Having established the stability of the Lamb wave with the traditional approximation in Appendix E (and in particular figure 12), the focus of our study below will be on the Lamb wave-like instability in the non-traditional setting with a Rayleigh damping at the top and a bottom forcing corresponding to the semi-analytical eigenmode that can be obtained by the numerical solution to the eigenvalue problem.

Figure 6 depicts the fields at $3\,600\,{\rm s}=1\,{\rm h}$ for the LWLI-SA case with the non-traditional setting. Here, we observe a significant increase of the amplitude of the perturbation due to the instability. The structure and amplitude of these fields are similar to the corresponding solution of the eigenvalue problem in (3.1), i.e. at $t=1$ h and $F=2\varOmega _y=1.458\times 10^{-4}$ s$^{-1}$. These results therefore closely match the structure of the theoretically predicted instability from § 2.

Figure 6. Output of the fields for the transformed prognostic variables at $t=1$ h for the LWLI-SA case, i.e. a Lamb wave-like instability in the non-traditional setting and the semi-analytical bottom forcing. Contour units are as in figure 4.

For completeness, we include the LWLI-SO result with a sub-optimal bottom forcing in figure 7. This is done by only applying a forcing to the quantities $u^\prime$ and ${\rm \pi} ^\prime$. This test aims to mimic more physical, real-world scenarios where there might be an injection of energy from the bottom atmospheric boundary layer, but that may not necessarily correspond to the non-traditional semi-analytical eigenmode structure. Outside of the bottom 3 km wherein the forcing is applied, the solution structure after 1 h looks almost identical to the structure of the instability in figure 6. The results indicate that a similar instability develops even if the bottom forcing deviates from the optimal structure.

Figure 7. Output of the fields for the transformed prognostic variables at $t=1$ h for the LWLI-SO case, i.e. a Lamb wave-like instability in the non-traditional setting and sub-optimal bottom forcing. Contour units are as in figure 4.

In order to quantify our observations of the simulation results, we compute the relative norm as a measure of the energy of the perturbation. The relative norm is given by

(3.8)\begin{equation} \text{rel. norm} = \frac{\Vert {\boldsymbol{\chi}} \Vert}{\Vert{\boldsymbol{\chi}}_0 \Vert}, \end{equation}

where $\Vert {\cdot } \Vert$ represents the $L^2$-norm in (2.23), and $\boldsymbol {\chi }_0$ denotes the initial $\boldsymbol {\chi }$.

Figure 8 depicts the logarithm of the relative norm over time for the three configurations. The relative norm is computed over the vertical domain of $z \in [3.0, 25.0]$ km for all three runs, i.e. we effectively exclude the lowest 3 km from our analysis. This is to exclude energy contributions from the bottom forcing.

Figure 8. Semi-logarithmic plot of the relative norm against time: a run with stable Lamb wave initial conditions in the traditional approximation (LW; blue dashed curve); a run with unstable Lamb wave-like initial conditions in the non-traditional setting and semi-analytical eigenmode as the bottom forcing (LWLI-SA; red solid curve); a run with the unstable Lamb wave-like initial conditions in the non-traditional setting and sub-optimal bottom forcing (LWLI-SO; green dotted curve); and a best-fit curve (dash-dotted in black) of the increment in the LWLI-SA run. The relative norm for all runs is taken in the limited $z\in [3.0,\,25.0]$ km domain.

As demonstrated in Appendix E, the relative norm of the traditional approximation LW run (blue dashed curve) remains close to the initial value, indicating that the energy of the system is almost exactly conserved over the time scale considered here. For the non-traditional setting LWLI-SA run with the semi-analytical eigenmode structure applied as the bottom forcing (red curve), a best-fit curve is plotted to quantify the instability growth rate (dash-dotted curve in black). Specifically, the best-fit curve is given by

(3.9)\begin{equation} \mathscr{F}(t) = \mathscr{F}_0 \exp({\zeta}t), \end{equation}

where $\mathscr {F}_0$ and $\zeta$ are obtained from the curve-fitting method. From this best-fit curve, we find an instability growth rate of about $7.83 \times 10^{-4}$ s$^{-1}$, which deviates from the theoretically predicted value in (2.44) by merely $0.8\,\%$. Finally, a similar relative norm curve for the non-traditional setting LWLI-SO run with the sub-optimal bottom forcing is depicted (dotted green curve), and the instability growth rate is equally close to the theoretical value as it yields the same growth rate of $7.83 \times 10^{-4}$ s$^{-1}$ up to the third digit (best-fit curve is not shown here).

All experiments presented in this section were performed with the same spatial and temporal resolutions. However, additional runs were conducted with spatial resolutions of $(151 \times 60)$ and $(301 \times 120)$. Furthermore, for each spatial grid resolution, runs with the following time-step sizes $\Delta t$ were carried out (results apart from the configuration studied above are not shown),

(3.10)\begin{equation} \Delta t = (1, 2, 4, 8, 10, 12, 14, 16)\,\text{s}. \end{equation}

With values in the range of $(7.83\pm 0.01) \times 10^{-4}$ s$^{-1}$, the growth rates obtained from all of these runs do not differ significantly from those presented above. The task of estimating the instability growth rate is difficult and may be sensitive to the experimental parameters (McNally, Lyra & Passy Reference McNally, Lyra and Passy2012; Zhelyazkov Reference Zhelyazkov2015). Therefore, the consistent experimental growth rate obtained here across varying temporal and spatial resolutions strongly implies a dynamical nature of the unstable mode.

The numerical results presented in this section resolutely suggest that the unstable mode described in § 2 may be expected to arise in atmospheric flow simulations that implement the non-traditional Coriolis acceleration. In particular, any problem formulation with a boundary layer model that does not per design suppress an energy influx into the bulk of the troposphere from below will permit these growing modes.