An unstable mode of the stratified atmosphere under the non-traditional Coriolis acceleration

Abstract

The traditional approximation neglects the cosine components of the Coriolis acceleration, and this approximation has been widely used in the study of geophysical phenomena. However, the justification of the traditional approximation is questionable under a few circumstances. In particular, dynamics with substantial vertical velocities or geophysical phenomena in the tropics have non-negligible cosine Coriolis terms. Such cases warrant investigations with the non-traditional setting, i.e. the full Coriolis acceleration. In this manuscript, we study the effect of the non-traditional setting on an isothermal, hydrostatic and compressible atmosphere assuming a meridionally homogeneous flow. Employing linear stability analysis, we show that, given appropriate boundary conditions, i.e. a bottom boundary condition that allows for a vertical energy flux and non-reflecting boundary at the top, the atmosphere at rest becomes prone to a novel unstable mode. The validity of assuming a meridionally homogeneous flow is investigated via scale analysis. Numerical experiments were conducted, and Rayleigh damping was used as a numerical approximation for the non-reflecting top boundary. Our three main results are as follows: (i) experiments involving the full Coriolis terms exhibit an exponentially growing instability, yet experiments subjected to the traditional approximation remain stable; (ii) the experimental instability growth rate is close to the theoretical value; (iii) a perturbed version of the unstable mode arises even under sub-optimal bottom boundary conditions. Finally, we conclude our study by discussing the limitations, implications and remaining open questions. Specifically, the influence on numerical deep-atmosphere models and possible physical interpretations of the unstable mode are discussed.

1. Introduction

The Coriolis force arises from the Earth's rotation. Named after Coriolis (1835), the effect of the Coriolis force is the deflection of a particle moving along the surface of the globe. Such an effect is perpendicular to the axis of the globe's rotation and the direction of the particle's velocity. Atmospheric and oceanic phenomena are subjected to the effects of the Coriolis force. Yet, as the atmosphere and ocean are shallow relative to the Earth's radius, the dynamics we are interested in is largely horizontal.

Based on the argument of the difference in the horizontal and vertical scales of the dynamics, Laplace (1835) concluded that, in studying the Coriolis effects, we may neglect the vertical acceleration arising from the Coriolis force, keeping only the horizontal. Specifically, the full Coriolis acceleration results in a tilted outward-pointing vector along the surface of the globe. This vector may be decomposed into sine and cosine components. The former are perpendicular to the surface of the globe and interact with only horizontal motion, while the latter are parallel and associated with vertical motion. Appendix A provides an elaboration on the full Coriolis acceleration. Following the arguments by Laplace, one may drop the cosine terms.

The approximation introduced by Laplace has since been widely used. This led Eckart (1960) to coin the term ‘traditional approximation’ when describing this approximation. The traditional approximation is also introduced in other standard textbooks on atmospheric and oceanic phenomena, see e.g. Pedlosky (2013), Vallis (2017) and Achatz (2022). However, the validity of the traditional approximation may be questioned under certain circumstances. For instance, at the equator, where the contributions from the cosine terms are the largest, or for phenomena involving substantial vertical motion, wherein the vertical Coriolis acceleration should not be ignored.

The limitations of the traditional approximation are gaining increasing attention, and the review by Gerkema et al. (2008) may be a good starting point for a study beyond the traditional approximation. A few other relevant references to the literature are mentioned below. White & Bromley (1995) proposed a set of dynamically consistent quasi-hydrostatic equations with full Coriolis support. They also demonstrated via scale analysis that, in the tropics, the cosine Coriolis terms may be as large as 10 % of the important terms in the equations, and thus they are non-negligible. The condition for linear stability has been derived by Fruman & Shepherd (2008) for a zonally symmetric, compressible, and non-traditional set-up. They noticed that an instability arises for a particular angular momentum profile near the equator. Fruman (2009) computed the analytical solution for meridionally confined zonally propagating waves for the linearised hydrostatic Boussinesq equations including the full Coriolis terms. Studies involving non-traditional effects were also conducted by de Verdière & Schopp (1994), Maas & Harlander (2007), Igel & Biello (2020), Rodal & Schlutow (2021) and Perez, Delplace & Venaille (2022).

Challenges to the assumptions made in the traditional approximation are not limited to the analysis of physical phenomena. In the example of a numerical model, Borchert et al. (2019) introduced an upper-atmosphere extension to the ICOsahedral Non-hydrostatic (ICON) model (Zängl et al. 2015) that includes the full Coriolis term. Another deep-atmosphere numerical model with full Coriolis support has been published by Smolarkiewicz et al. (2016). From here on, we use the term ‘non-traditional setting’ to refer to a setting where the full effect of the Coriolis force is considered.

In this manuscript, we study the stability of an isothermal hydrostatic atmosphere under the non-traditional setting and assuming a meridionally homogeneous flow. Via linear stability analysis of the compressible inviscid Euler equations, we show that, given a non-reflective boundary condition at the top and a forcing boundary at the bottom, a small perturbation of the isothermal hydrostatic atmosphere at rest leads to an exponentially growing instability. We identify the structure of this unstable mode and quantify its theoretical growth rate. In contrast to this choice of boundary conditions, it can be shown by an argument from functional analysis in combination with energy conservation that the atmosphere is unconditionally linearly stable if a free-slip boundary condition is assumed at the flat top and bottom boundaries. Section 2 contains these theoretical developments. The assumption of meridionally homogeneous flow is also discussed in detail there.

For the case of a non-reflecting upper boundary and a bottom forcing that allows for a vertical flux of energy into the domain, the unstable mode is investigated through numerical experiments. Section 3.1 and Appendix B introduce the numerical model by Benacchio & Klein (2019), and Appendix C extends the numerical model to support the non-traditional setting. Appendix D details the necessary numerical extension to ensure a well-balanced Lamb wave solution, and Appendix E establishes the energy-conserving properties of the stable Lamb wave solution by the numerical model. The non-reflective top boundary is realised numerically via a Rayleigh damping layer.

Experiments made in § 3 demonstrate that, in the presence of the full Coriolis acceleration, the isothermal, hydrostatic atmosphere at rest becomes unstable with an experimental instability growth rate that is close to the theoretically predicted value. The instability growth is also observed in an experiment where the bottom forcing applied deviates from the one that was used to closely represent the analytical theory within the numerical simulation framework. Results with this sub-optimal forcing suggest that related unstable modes exist for a variety of bottom boundary conditions. In nature, the atmospheric boundary layer flow will provide effective bottom boundary conditions, for instance, and in conjunction with non-homogeneous surface conditions, these may induce an effective bottom forcing akin to the sub-optimal forcing employed here.

The existence of this unstable mode due to the full Coriolis effect has two main implications. First, this unstable mode may provide a theoretical understanding of hitherto unexplained physical phenomena. Second, the presence of the unstable mode is relevant to numerical weather and climate prediction models that solve a full representation of the Coriolis force, in particular numerical models of the deep atmosphere. Along with a summary of the results of this manuscript, § 4 provides further elaborations and discussions of these implications.

2. Theory

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 2002), 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 2010, 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 (1989). 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 (2008) 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 2017). 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 2003).

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. (2022) 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 1982). 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 (1966). Also, Kato (1982) 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. Numerical experiment and results

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 (2019) 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 2010). Specifically, the Rayleigh damping by Durran & Klemp (1983), 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 2012; Zhelyazkov 2015). 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.

4. Conclusions

This paper presents a novel unstable mode that arises from considering the non-traditional terms of the Coriolis force. This unstable mode lends itself to physical interpretations and may have substantial implications for numerical schemes that employ the non-traditional setting.

Considering the linearised compressible Euler equations in the non-traditional setting, an inner product of the state vector of the prognostic variables for perturbations is derived inducing a norm. We showed that this norm is proportional to the perturbation's total energy. We also found that the differential operator emerging from the linear system, after the usual separation of variables, is not skew-Hermitian with respect to the inner product when the boundaries are non-trivial. By some version of the spectral theorem, this observation gives rise to a potentially unstable spectrum. And indeed, by means of asymptotic theory, we showed that the isothermal, hydrostatic atmosphere at rest subject to the non-traditional Coriolis acceleration becomes prone to an unstable Lamb wave-like mode under infinitesimal initial perturbations. The theoretical growth rate of the unstable mode is quantified.

With a small perturbation of the isothermal hydrostatic atmosphere at rest, numerical experiments demonstrate that this initial configuration supports a linearly unstable mode when the full Coriolis acceleration is considered. The experimental instability growth rate is close to the theoretical value. In contrast, the simulation remains stable if the traditional approximation is applied. Even if the injection of energy at bottom boundary is sub-optimal, i.e. if it deviates from the theoretical structure, a similar unstable mode develops nevertheless. Therefore, this unstable mode may present itself under conditions that are less stringent than the theory implies.

Physical evidence from observations for this novel type of unstable mode can be found in Nishida, Kobayashi & Fukao (2014). They observed Lamb waves at around $30^\circ$ North with a vast array of high-resolution microbarometers. In contrast to intermittent waves due to large-scale events such as volcanic eruptions, they found omnipresent background waves. The authors discussed several excitation mechanisms as a cause for the phenomenon observed. The results presented in this paper imply that these Lamb waves might be generated by the instability due to the Coriolis force, as the frequency spectrum of the observations fits well with our theoretical predictions. Further corroboration is that most waves the authors observed travel along the East–West axis as predicted by our theory.

Furthermore, our investigation presents relevant implications for numerical models that use the non-traditional setting. One would encounter this unstable mode in numerical simulations if there is energy flux from the atmospheric boundary layer. This effect would be especially significant for the numerical simulations of the deep atmosphere, e.g. the papers by Borchert et al. (2019) and Smolarkiewicz et al. (2016) as mentioned in the introduction, or in tropical dynamics with the non-traditional setting. Therefore, the characterisation of the unstable mode presented here is a worthwhile pursuit that may broaden our understanding of the simulation results from these numerical models.

Our derivation of the unstable mode and the numerical experiments have two limitations that may be investigated in a future study. (i) The effect of the Earth's curvature has not been considered, and (ii) meridional homogeneity was assumed. On point (i), it is possible to account for the Earth's curvature by incorporating its effect in the $\beta$-plane term following (Maas & Harlander 2007). However, in contrast to these authors, we have considered the compressible rather than the incompressible Euler equations in this study, and adapting their approach to the present application will require substantial work, which we leave to a potential future study. On point (ii), we plan to carry out three-dimensional numerical experiments that will explore the influence of the $\beta$-plane effect on the eigenmodes found in the present study. An exhaustive examination of the physical interpretations of the unstable mode could also be carried out.

Finally, it may also be fruitful to investigate the unstable mode in the more operationally ready numerical models at the German Weather Service and the European Centre for Medium-Range Weather Forecasts. Two questions are particularly worth answering. Can the unstable mode described in this paper be reproduced by the numerical models? If the answer is yes, what influences does the unstable mode have on atmospheric dynamics and circulation?