Geophysical fluid models with simple energy backscatter: explicit flows and unbounded exponential growth

Abstract

Motivated by numerical schemes for large-scale geophysical flow, we consider the rotating shallow water and Boussinesq equations on the whole space with horizontal kinetic energy backscatter source terms built from negative viscosity and stabilising hyperviscosity with constant parameters. We study the impact of this energy input through various explicit flows, which are simultaneously solving the nonlinear equations and the linear equations that arise upon dropping the transport nonlinearity, i.e. the linearisation in the zero state. These include barotropic, parallel and Kolmogorov flows as well as monochromatic inertia gravity waves. With focus on stable stratification, we find that the backscatter generates numerous solutions of this type that grow exponentially and unboundedly, also with vertical structure. This signifies the possibility of undesired energy concentration into specific modes due to the backscatter. Families of steady-state flows of this type arise as well and superposition principles in the nonlinear equations provide explicit sufficient conditions for instability of some of these. For certain steady barotropic flows of this type, we provide numerical evidence of eigenmodes whose growth rates are proportional to the amplitude factor of the flow. For all other arising steady solutions, we prove that this is not possible.

Keywords:

• Explicit flows
• instability
• plane waves
• rotating shallow water equations
• rotating Boussinesq equations

1. Introduction

Spatial resolution of geophysical flows is limited not only in observational data but also in numerical simulations for ocean and climate studies due to lack of computing power. Towards realistic simulations, this is compensated by so-called parameterisations of subgrid effects, i.e. subgrid models that are gauged a priori and simulate the influence of the missing small scale resolution on the resolved large-scale flow. Moreover, these should account for numerical discretisation effects such as over-dissipation and yet admit stable simulations of ocean and climate models. A practical solution to these problems that has come to frequent use are kinetic energy backscatter schemes that effectively introduce negative horizontal viscosity together with hyperviscosity (cf., e.g. Jansen and Held 2014, Zurita-Gotor et al. 2015, Jansen et al. 2019, Juricke et al. 2020, Perezhogin 2020); we refer to Danilov et al. (2019) for a detailed discussion and relations to other approaches. Simulations with backscatter have been found to provide energy “at the right place”, matching the result more closely to observations and high resolution comparisons.

Motivated by this, we consider the rotating shallow water and Boussinesq equations on the whole space with simplified kinetic energy backscatter source terms built from negative horizontal viscosity and stabilising hyperviscosity with constant parameters. The backscatter terms in the numerical scheme have non-constant coefficients from a coupled energy equation that aims to regulate energy consistency and thus these terms vanish for infinite resolution. Our simplified consideration on the continuum level is amenable to an analysis of qualitative features and energy distribution, so that the results point out potential issues and can guide further study of backscatter schemes.

In particular, this idealisation admits a direct analytical study of the influence of backscatter through its impact on various explicit flows. Indeed, explicit flows are frequently used as a tool for benchmarking analytical and numerical studies in this and other contexts (e.g. Weinbaum and O'Brien 1967, Majda 2003, Drazin and Riley 2006, Majda and Wang 2006, Dyck and Straatman 2019, Chai et al. 2020) and also for turbulence studies (e.g. Lelong and Dunkerton 1998, Ghaemsaidi and Mathur 2019, Onuki et al. 2021).

We start our investigations with the rotating shallow equations and then move to the rotating Boussinesq equations. In the spirit of Prugger and Rademacher (2021), we consider flows and waves, which simultaneously solve the nonlinear equations and the linear equations that arise upon dropping the transport nonlinearity. These include analogues of barotropic, parallel and Kolmogorov flows as well as monochromatic inertia gravity waves (MGWs) (cf., Yau et al. 2004, Balmforth and Young 2005, Achatz 2006, Prugger and Rademacher 2021). We aim to identify the occurrence and properties of steady, growing and decaying flows of these types. Our approach is essentially based on direct computations, though these are at times somewhat tedious.

More specifically, we first consider the rotating shallow water equations with backscatter and flat bottom topography in the f-plane approximation given by (1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) (1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) where $\mathit{v}=\mathit{v}(t,\mathit{x})\in {\mathbb{R}}^2$ is the velocity field on the whole space $\mathit{x}\in {\mathbb{R}}^2$ at time $t\ge 0$ and $\eta =\eta (t,\mathit{x})\in \mathbb{R}$ is the deviation of the fluid layer from the characteristic fluid depth $H_0>0$, giving $H_0+\eta$ as the fluid layer thickness. In addition, $f\ne 0$ is the constant Coriolis parameter, g>0 is the gravity acceleration and the backscatter parameters are $b_1,b_2,d_1,d_2>0$.

We then turn to the rotating Boussinesq equations augmented with backscatter (2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) (2b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(2b) (2c) $\begin{array}{rl}\frac{\mathrm{\partial }\mathrm{b}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathrm{b}+N^2{v}_3& =\mu \mathrm{\Delta }\mathrm{b},\end{array}$(2c) with the horizontal backscatter parameters $b_1,b_2,d_1,d_2>0$ and vertical viscosity $\nu \ge 0$ in the diagonal matrix operator. Other quantities in (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ) are the velocity field $\mathit{v}(t,\mathit{x})\in {\mathbb{R}}^3$ for $\mathit{x}\in {\mathbb{R}}^3$, $t\ge 0$, pressure and buoyancy $p(t,\mathit{x}),\mathrm{b}(t,\mathit{x})\in \mathbb{R}$, the vertical unit vector ${\mathit{e}}_{\mathbf{3}}$ and thermal diffusivity $\mu \ge 0$. As usual, the buoyancy considered here is of the form $\mathrm{b}(t,\mathit{x})=-g\left(\rho \right(t,\mathit{x})-\overline{\rho }(z\left)\right)/{\rho }_0\in \mathbb{R}$ with fluid density $\rho (t,\mathit{x})\in \mathbb{R}$ and reference density field $\overline{\rho }\left(z\right)$ depending on the vertical space direction z only, and the characteristic density ${\rho }_0$. Then $N^2=-\left(g/{\rho }_0\right)\mathrm{d}\overline{\rho }/\mathrm{d}z$ is the Brunt–Väisälä frequency for stable stratification $\mathrm{d}\overline{\rho }/\mathrm{d}z<0$.

While the physically natural situation is horizontally isotropic backscatter, $d_1=d_2$, $b_1=b_2$, we allow for (weak) anisotropy. Our motivation for this is to account for an effective anisotropy from the numerical backscatter scheme due to the combination of an anisotropic grid, spatially anisotropic coefficients and the application of filtering (see, e.g. Danilov et al. 2019). It turns out that among the phenomena we find, the only one which requires anisotropic backscatter are explicit growing and certain decaying flows in the shallow water case (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) as discussed below. However, here the anisotropy can be arbitrarily weak, i.e. $|d_1-d_2|+|b_1-b_2|$ can be arbitrarily small. As expected, the additional parameters of anisotropic backscatter often provide more flexibility, but in the Boussinesq case (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ), all general phenomena that we consider occur also in the isotropic case.

The explicit solutions that we consider in the shallow water case (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) are derived from the plane wave ansatz (3) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\eta =\frac{f}{g}\varphi (\mathit{k}\mathbf{\cdot }\mathit{x}),$(3) which implies vanishing nonlinear terms, and in the simplest case are single Fourier modes for ψ and ϕ. For the Boussinesq equations (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ), we find barotropic horizontal flows of similar plane wave type, as well as flows with vertical structure, which relate to parallel flow, Kolmogorov flow and MGWs.

As mentioned, all these flows simultaneously solve the nonlinear equation (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) or (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ), and the respective linearisation when dropping the transport nonlinearity. Consequently, these solutions come as families with a free amplitude factor, but superpositions are heavily constrained in the full nonlinear setting. We identify superposition principles based on ideas from Prugger and Rademacher (2021) and thus find invariant subspaces with linear dynamics in the nonlinear system, in particular those with explicit solutions that grow exponentially and unboundedly. For the shallow water case, this requires at least weakly anisotropic backscatter as mentioned above. We refer to perturbations of a given state, which are realised by possible superpositions with such unboundedly growing explicit solutions of linear type in the nonlinear system, as unbounded instability. Notably, this does not occur in generic damped-driven evolution equations, where unstable manifolds are nonlinear.

The occurrence of unbounded instability highlights the possibility of undesired concentration of energy due to backscatter and is in contrast to the targeted energy redistribution. Solutions with negative growth rate likewise illustrate the possibility of ineffective energy input into certain scales. However, while solutions to consistent discretisations shadow continuum effects over finite time intervals, the specific implications for numerical backscatter schemes require additional investigations beyond the present study.

Concerning the steady states among these explicit flows, we study aspects of their stability in terms of unbounded instability and the linear stability with respect to the aforementioned amplitude factor. For certain steady barotropic flows, we provide numerical evidence for eigenmodes whose growth rates are proportional to the amplitude factor, thus featuring arbitrarily strong local instability. For all other steady solutions of the type considered, we prove that growth rates scale sublinearly with the amplitude factor.

We note that the above analytical realisations of backscatter replace the usual molecular viscosity operator $\nu \mathrm{\Delta }$ by operators of the form $-(d{\mathrm{\Delta }}^2+b\mathrm{\Delta })$, familiar from the scalar Kuramoto–Sivashinsky (KS) equations: $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\frac{1}{2}|\mathbf{\nabla }u{|}^2=-{\mathrm{\Delta }}^{2}u-\mathrm{\Delta }u,\mathit{x}\in {\mathbb{R}}^n,n\le 3.$ These have been derived in various contexts and in particular, the one-dimensional case appears broadly, e.g. for interfacial layers (Wei 2006). Posed on tori, for solutions with globally bounded gradient, the deviations from the spatial mean admit a finite dimensional global attractor (Nicolaenko et al. 1985) and the one-dimensional KS is a paradigm for chaos in a partial differential equation (cf., e.g. Nicolaenko et al. 1985, Smyrlis and Papageorgiou 1991, Kalogirou et al. 2015, and references therein). Differentiating KS yields a system for $\mathit{v}=\mathbf{\nabla }u$ with the fluid transport nonlinearity, thus relating more closely to (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) and (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ) in case all backscatter coefficients are equal, although this relation is clearly far from complete. Moreover, the solutions that we are focussing on, in particular the unboundedly growing ones, are all non-gradient and do not exist on one-dimensional space, therefore they are unrelated even on these levels.

This paper is organised as follows. In section 2, we discuss the horizontal flows of the rotating shallow water equations with backscatter. Section 3 is devoted to the rotating Boussinesq equations with backscatter and the analysis of existence, growth and unbounded instability properties of the aforementioned different types of flows.

2. Rotating shallow water with backscatter

In this section, we consider the rotating shallow water equations with backscatter (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) and first identify certain explicit flows. The inviscid rotating shallow water equations without backscatter, i.e. $b_1=b_2=d_1=d_2=0$, possess the explicit plane wave steady solutions $\mathit{v}={\varphi }^{\prime }(\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\eta =\frac{f}{g}\varphi (\mathit{k}\mathbf{\cdot }\mathit{x}),$ for any wave vector $\mathit{k}\in {\mathbb{R}}^2$ and sufficiently smooth wave shape ϕ (e.g. Prugger and Rademacher 2021). These are also in geostrophic balance, corresponding to Rossby waves.

For the case of backscatter in (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ), we seek solutions of the similar form (3(3) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\eta =\frac{f}{g}\varphi (\mathit{k}\mathbf{\cdot }\mathit{x}),$(3) ) for any wave vector $\mathit{k}=(k_1,k_2{)}^{\mathsf{T}}\in {\mathbb{R}}^2$ and sufficiently smooth wave shapes ψ and ϕ. The time-independence of η results from equation (1b(1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) ), since $\mathit{v}$ is divergence free and the nonlinear terms vanish in this case. Inserting (3(3) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\eta =\frac{f}{g}\varphi (\mathit{k}\mathbf{\cdot }\mathit{x}),$(3) ) into (1a(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) yields the linear equation $\begin{array}{r}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}{\mathit{k}}^{\perp }+f(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\xi }-\psi )\mathit{k}=-|\mathit{k}{|}^2\left(\begin{array}{cc}{d}_1|\mathit{k}{|}^2{\displaystyle \frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}}+b_1{\displaystyle \frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2}}& 0\\ 0& d_2|\mathit{k}{|}^2{\displaystyle \frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}}+b_2{\displaystyle \frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2}}\end{array}\right){\mathit{k}}^{\perp }.\end{array}$ Every vector in ${\mathbb{R}}^2$ on the right-hand side has a unique representation by the orthogonal basis vectors on the left-hand side. The scalar product with ${\mathit{k}}^{\perp }$ and $\mathit{k}$, respectively, gives (4a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-|\mathit{k}{|}^2(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(4a) (4b) $\begin{array}{rl}f(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\xi }-\psi )& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}+(b_1-b_2\left)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2}\right).\end{array}$(4b) We focus on monochromatic solutions, i.e. that contain a single Fourier mode, and later investigate possible superpositions. Equations (4(4a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-|\mathit{k}{|}^2(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(4a) ) restrict such solutions to the form (5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) with arbitrary shifts $\tau ,s\in \mathbb{R}$ and the rest of the real parameters must satisfy (6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) (6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) (6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) We note that (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) is a dispersion relation of growth/decay and wave vector, and (6b(6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) ) an amplitude relation of the two amplitudes ${\alpha }_1$ and ${\alpha }_2$ in terms of the wave vector; (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ) is an auxiliary compatibility condition. Specifically, (1a(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) possesses explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) if the parameters satisfy (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ,b), and the time-independence of η coming from (1b(1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) ) requires condition (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ), which means ${\alpha }_2$ or λ is zero. In particular, (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ) means that these explicit solutions with non-trivial depth variation η are steady. Notably, solutions with $\lambda >0$, whose existence is studied in section 2.1, grow exponentially and unboundedly. As mentioned before, we refer to this as unbounded instability of the zero state, and more generally of any other solution which admits superposition with such growing explicit solutions. In figure 1, we plot the loci of these different solutions. The blue and red regions show the sign of the growth rate λ, which characterise the exponentially decaying and growing explicit solutions. In fact, we show in section 2.2.1 that the red region describes a subset of real unstable eigenmodes of the linearisation of the zero state.

Figure 1. We plot the occurrence of explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) in the plane of wave vectors with fixed parameters $d_1=1.0,b_1=1.5,b_2=2.2,f=0.3,g=9.8,H_0=0.1,{\alpha }_1=1.0$ and $d_2,{\alpha }_2$ as in the subcaptions. Red regions: $\lambda >0$, i.e. unbounded growth; blue regions: $\lambda <0$; black curves: $\lambda =0$, i.e. steady states. The white curves mark solutions with ${\alpha }_2=0$, and the white bullets mark steady solutions with $\lambda =0$. The black dashed curves mark solutions of condition (6b(6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) ) only, for ${\alpha }_2\ne 0$ with values as in the subcaptions, which also solve (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) at the white bullets ($\lambda =0$). In (c) we mark a line of wave vectors in a fixed direction (grey), and the growing or decaying solutions (white circles) on it, whose superpositions with or without the steady state on the grey line again yield explicit solutions (Colour online).

We proceed as follows: in section 2.1, we will discuss the sets of solutions in terms of their wave vectors, the organising parameters and possible superpositions. In section 2.2, we then analyse the unbounded instability and the linear stability of explicit steady solutions. As mentioned in the introduction, we will see that in this case, unbounded instability requires anisotropy $b_1\ne b_2$, $d_1\ne d_2$, although the differences can be arbitrarily small.

2.1. Sets of solutions and superpositions

In order to analyse the existence of solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) of equations (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) in more detail, it is convenient to write the dispersion relation (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) and the amplitude relation (6b(6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) ) in the form (7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) (7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) with real parameter $\sigma =f({\alpha }_1-{\alpha }_2)/{\alpha }_1$ describing the relative difference between the amplitudes of the velocity vector $\mathit{v}$ and the fluid depth variation η. Steady solutions satisfy (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$ and it is then natural to view σ as an adjustment, defined by (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ), of the relation between the amplitudes ${\alpha }_1,{\alpha }_2$ depending in particular on the wave vector $\mathit{k}$. For the time-dependent case $\lambda \ne 0$, we have $\sigma =f$, since ${\alpha }_2=0$ is required due to (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ), and viewing (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) as a definition for λ. The natural free parameter is the wave vector $\mathit{k}$. We will see that the existence and growth or decay properties of solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ), as well as the locations of unboundedly unstable steady states of this kind are strongly connected with the values of σ, which we therefore consider as an organising parameter.

2.1.1. Superpositions of explicit flows

Before discussing existence conditions, we briefly note that superpositions of solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) are also solutions, if all wave vectors $\mathit{k}$ lie on the same line through the origin in the wave vector plane. We plot examples in figure 1(c). The reason is that for these superpositions, the nonlinear terms in (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) still vanish due to the orthogonality of wave vectors and flow directions (cf. Prugger and Rademacher 2021), and the remaining linear equations are satisfied by each superposed explicit solution. This radial superposition principle of wave vectors gives non-trivial subspaces of initial data to (1a(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) in which the dynamics are linear. In the example of figure 1(c), this space is three dimensional, since the negated wave vectors give linearly dependent solutions, and this is the maximum possible as shown below.

2.1.2. Steady explicit solutions

For steady states, we only need to investigate the wave vectors $\mathit{k}$ satisfying the dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$. These form a simple closed curve around the origin in wave vector space that is symmetric with respect to axis reflection, and whose interior is star shaped, i.e. all points of the set are connected with the origin through a direct line contained in the set, but it need not be convex. We plot an example in figure 1. To see this, note that for wave vectors $\mathit{k}=r{\mathit{k}}^{\prime }$, ${\mathit{k}}^{\prime }$ with $|{\mathit{k}}^{\prime }|=1$ fixed, the right-hand side of equation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) is linear in the squared wave vector length $|\mathit{k}{|}^2=r^2$ (after using $\mathit{k}=r{\mathit{k}}^{\prime }$ with $|{\mathit{k}}^{\prime }|=1$ and division by $r^2$). Furthermore, for any fixed ${\mathit{k}}^{\prime }$, there is exactly one $r_0>0$ so that $\lambda =0$ for $r=r_0$, and $\lambda >0$ for $0<r<r_0$ as well as $\lambda <0$ for $r_0<r$. This means that λ is positive in the interior of the closed curve of steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) (red regions in figure 1), except for the origin, where $\lambda =0$, and λ is negative outside (blue regions in figure 1).

In polar coordinates, $\mathit{k}=r\left(\begin{array}{c}\cos \phi \\ \sin \phi \end{array}\right)$ with angle $\phi \in [0,2\pi )$ and wave number $r\ge 0$; the curve for explicit non-trivial steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ), i.e. the wave vectors with $\lambda =0$, is parameterised by the angle φ with the wave number given by (8) $r=\sqrt{\frac{b_1{\sin }^2\phi +b_2{\cos }^2\phi }{d_1{\sin }^2\phi +d_2{\cos }^2\phi }}.$(8) Generally, these steady solutions have different values of σ, the relative difference of amplitudes ${\alpha }_1$ and ${\alpha }_2$; recall that time-dependent explicit solutions ($\lambda \ne 0$) all have the same value $\sigma =f$. In either case, the explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) form a linear space since their amplitudes only enter into the ratio $({\alpha }_2-{\alpha }_1)/{\alpha }_1$ (so into σ) and are therefore naturally parameterised by an arbitrary amplitude parameter $a\ge 0$ that is a common factor of both ${\alpha }_1$ and ${\alpha }_2$, and thus does not change the value of σ.

In the following, we further investigate the conditions (7(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) for the existence of explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ). First, we analyse the occurrence and shapes of the curves defined by the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) and use this to determine the time-dependent explicit solutions, i.e. ${\alpha }_2=0$, as well as the steady explicit solutions for which (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) is satisfied with $\lambda =0$. Second, we discuss the values of σ for which steady solutions exist; clearly any steady solution of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) has a corresponding value of σ. But not every σ admits such a steady solution and the value of σ for the time-dependent solutions ($\lambda \ne 0$) is fixed at $\sigma =f$, since these solutions require ${\alpha }_2=0$.

2.1.3. Set of solutions

In order to investigate the set of explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ), primarily of the time-dependent ones with ${\alpha }_2=0$ and $\lambda \ne 0$, we analyse the shapes of the curves defined by the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ). We start with two special cases:

1. In the isotropic case $b_1=b_2$ and $d_1=d_2$, equation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) requires $\sigma =0$, i.e. ${\alpha }_1={\alpha }_2$, so that in this case all non-trivial solutions are steady, i.e. $\lambda =0$, and have $\mathit{k}=0$ or $\mathit{k}$ on the circle with radius $\sqrt{b_1/d_1}$ defined by dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) (with $\lambda =0$), see figure 2(a). Thus, non-steady solutions ($\lambda \ne 0$) of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) arise from anisotropy in the backscatter.

2. There is also a special anisotropic case. If $d_1\ne d_2$, then the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) is satisfied in the origin and on the circle $|\mathit{k}{|}^2=(b_1-b_2)/(d_1-d_2)$ with $\sigma =0$, and λ defined by the dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) is always constant on that circle. If additionally $b_1/d_1=b_2/d_2$, then all solutions of (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) on the circle $|\mathit{k}{|}^2=(b_1-b_2)/(d_1-d_2)$ also solve (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$, so all of these give explicit steady solutions, which have $\sigma =0$. In case $b_1/d_1\ne b_2/d_2$, the value of σ for the steady states is not constant, as mentioned above.

Figure 2. Possible structures of solution curves of the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) analogous to figure 1. Fixed parameters: $d_1=1.0,b_1=1.5,f=0.5,g=9.8,H_0=0.1,{\alpha }_1=1.0$. It is $\sigma =0$ in (a) and $\sigma =-1$ in (b) and (c) (Colour online).

It remains to discuss the general anisotropic case, for which we consider the wave vectors in polar coordinates $\mathit{k}=r(\cos \phi ,\sin \phi {)}^{\mathsf{T}}$ as above; by symmetry of the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) it suffices to take $\phi \in [0,\pi ]$. The special cases $\phi \in \{0,\pi /2,\pi \}$ requires $\sigma =0$, i.e. steady solutions (since then ${\alpha }_2\ne 0$), and the corresponding wave vectors are $\mathit{k}\in \{\left(\begin{array}{c}0\\ 0\end{array}\right),\pm \sqrt{\frac{b_1}{d_1}}\left(\begin{array}{c}0\\ 1\end{array}\right),\pm \sqrt{\frac{b_2}{d_2}}\left(\begin{array}{c}1\\ 0\end{array}\right)\}.$ We now consider $\phi \in (0,\pi )\mathrm{\setminus }\left\{\pi /2\right\}$ only. In the case $d_1=d_2$, solutions of (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) are (9) $r=\pm \sqrt{\frac{2\sigma }{(b_1-b_2)\sin \left(2\phi \right)}}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{g}\mathrm{n}\left(\sigma \right)=\mathrm{s}\mathrm{g}\mathrm{n}\left(\right(b_1-b_2)\sin (2\phi \left)\right),$(9) with $\mathrm{s}\mathrm{g}\mathrm{n}(\cdot )$ the sign function (see, e.g. figure 2(b)). For $d_1\ne d_2$, we get (10) $r=\sqrt{\frac{b_1-b_2}{2(d_1-d_2)}\pm \sqrt{\frac{(b_1-b_2{)}^2}{4(d_1-d_2{)}^2}-\frac{2\sigma }{(d_1-d_2)\sin \left(2\phi \right)}}},$(10) which gives real solutions to (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ), if and only if the expressions in the square roots of (10(10) $r=\sqrt{\frac{b_1-b_2}{2(d_1-d_2)}\pm \sqrt{\frac{(b_1-b_2{)}^2}{4(d_1-d_2{)}^2}-\frac{2\sigma }{(d_1-d_2)\sin \left(2\phi \right)}}},$(10) ) are non-negative. This means that for fixed angle φ, we have two cases:

1. $b_1=b_2$ or $\mathrm{s}\mathrm{g}\mathrm{n}(b_1-b_2)=-\mathrm{s}\mathrm{g}\mathrm{n}(d_1-d_2)$ requires for at least one solution that (11) $\mathrm{s}\mathrm{g}\mathrm{n}\left(\sigma \right)=-\mathrm{s}\mathrm{g}\mathrm{n}((d_1-d_2)\sin \left(2\phi \right)).$(11)

2. $b_1\ne b_2$ and $\mathrm{s}\mathrm{g}\mathrm{n}(b_1-b_2)=\mathrm{s}\mathrm{g}\mathrm{n}(d_1-d_2)$ requires for at least one solution that (12a) $\begin{array}{rl}\frac{(b_1-b_2{)}^2}{8(d_1-d_2)}\sin \left(2\phi \right)& \ge \sigma \mathrm{f}\mathrm{o}\mathrm{r}\phi \in (0,\frac{1}{2}\pi )\\ & (\mathrm{i}\mathrm{f}{d}_1>d_2,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\phi \in (\frac{1}{2}\pi ,\pi )),\end{array}$(12a) (12b) $\begin{array}{rl}\frac{(b_1-b_2{)}^2}{8(d_1-d_2)}\sin \left(2\phi \right)& \le \sigma \mathrm{f}\mathrm{o}\mathrm{r}\phi \in (\frac{1}{2}\pi ,\pi )\\ & (\mathrm{i}\mathrm{f}{d}_1>d_2,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\phi \in (0,\frac{1}{2}\pi )).\end{array}$(12b)

Two solutions for a fixed angle φ occur if and only if the following three conditions are satisfied: (13a) $\begin{array}{rl}\mathrm{s}\mathrm{g}\mathrm{n}(b_1-b_2)& =\mathrm{s}\mathrm{g}\mathrm{n}(d_1-d_2)\mathrm{s}\mathrm{o}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}{b}_1\ne b_2,\end{array}$(13a) (13b) $\begin{array}{rl}\mathrm{s}\mathrm{g}\mathrm{n}\left(\sigma \right)& =\mathrm{s}\mathrm{g}\mathrm{n}((d_1-d_2)\sin \left(2\phi \right)),\end{array}$(13b) (13c) $\begin{array}{rl}\frac{(b_1-b_2{)}^2}{8|d_1-d_2|}|\sin \left(2\phi \right)|& \ge |\sigma |.\end{array}$(13c) In figure 2 we plot examples, where up to one ( figure 2(a,b)) or up to two ( figure 2(c)) solutions of (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) for certain angles φ arise for fixed σ. The conditions (9(9) $r=\pm \sqrt{\frac{2\sigma }{(b_1-b_2)\sin \left(2\phi \right)}}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{g}\mathrm{n}\left(\sigma \right)=\mathrm{s}\mathrm{g}\mathrm{n}\left(\right(b_1-b_2)\sin (2\phi \left)\right),$(9) )–(13(13a) $\begin{array}{rl}\mathrm{s}\mathrm{g}\mathrm{n}(b_1-b_2)& =\mathrm{s}\mathrm{g}\mathrm{n}(d_1-d_2)\mathrm{s}\mathrm{o}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}{b}_1\ne b_2,\end{array}$(13a) ) thus determine the occurrence and structures of the solution curves of the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) depending on the parameter settings. For instance, changing the sign of certain expressions, but not their absolute values, merely rotates the structures by $\pi /2$.

2.1.4. Structure and values of σ

We next discuss the occurrence of explicit steady solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) in the anisotropic case in more detail, in particular the values of σ, for which explicit steady solutions exist. These are the values of σ for which the curves defined by the dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$ and the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) intersect, see figure 3. Recall that time-dependent explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) all have $\sigma =f$.

Figure 3. Solution curves (dashed lines) of the amplitude relation (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) in terms of σ, as well as their intersections with solution curves of the dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$; the grey lines mark $k_2=m{k}_1$. Denotations as in figure 1 with focus on upper left quadrant. In (a) there are no explicit solutions. In (b) and (c) steady solution occur at intersection point (white dot). Fixed parameters: $d_1=1.0,d_2=1.0,b_1=1.5,b_2=2.2,f=0.3,g=9.8,H_0=0.1,{\alpha }_1=1.0$ (Colour online).

To ease computations, we consider a line $k_2=m{k}_1$ with slope $m\in \mathbb{R}$, i.e. $m=\tan \phi$ in the polar coordinates for wave vector used before. Inserting this into (7(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) gives the values of $k_1$ for which line and curves intersect (14a) $\begin{array}{rl}{k}_1^2& =\frac{b_1{m}^2+b_2}{(d_1{m}^2+d_2)(1+m^2)}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(14a) (14b) $\begin{array}{rl}{k}_1^2& =\frac{\sigma }{m(b_1-b_2)}\mathrm{f}\mathrm{o}\mathrm{r}{d}_1=d_2,\end{array}$(14b) (14c) $\begin{array}{rl}{k}_1^2& =\frac{b_1-b_2}{2(d_1-d_2)(1+m^2)}\pm \sqrt{\frac{m(b_1-b_2{)}^2-\sigma (d_1-d_2\left)\right(1+m^2)}{4m(d_1-d_2{)}^2(1+m^2{)}^2}}\mathrm{f}\mathrm{o}\mathrm{r}{d}_1\ne d_2.\end{array}$(14c) Here, (14a(14a) $\begin{array}{rl}{k}_1^2& =\frac{b_1{m}^2+b_2}{(d_1{m}^2+d_2)(1+m^2)}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(14a) ) is the intersection of the line with the curve defined by (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ), while (14b(14b) $\begin{array}{rl}{k}_1^2& =\frac{\sigma }{m(b_1-b_2)}\mathrm{f}\mathrm{o}\mathrm{r}{d}_1=d_2,\end{array}$(14b) ,c) provides the intersection of the line with the curve defined by (7b(7b) $\begin{array}{rl}0& =((d_1-d_2)|\mathit{k}{|}^2+b_2-b_1)k_1{k}_2+\sigma ,\end{array}$(7b) ) in the two cases (see figure 3(a)).

We next choose σ such that both intersection points are at the same position on the ray (see figure 3(b)). This occurs when the right-hand side of (14a(14a) $\begin{array}{rl}{k}_1^2& =\frac{b_1{m}^2+b_2}{(d_1{m}^2+d_2)(1+m^2)}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(14a) ) equals the right-hand side of (14b(14b) $\begin{array}{rl}{k}_1^2& =\frac{\sigma }{m(b_1-b_2)}\mathrm{f}\mathrm{o}\mathrm{r}{d}_1=d_2,\end{array}$(14b) ) or (14c(14c) $\begin{array}{rl}{k}_1^2& =\frac{b_1-b_2}{2(d_1-d_2)(1+m^2)}\pm \sqrt{\frac{m(b_1-b_2{)}^2-\sigma (d_1-d_2\left)\right(1+m^2)}{4m(d_1-d_2{)}^2(1+m^2{)}^2}}\mathrm{f}\mathrm{o}\mathrm{r}{d}_1\ne d_2.\end{array}$(14c) ). In both cases, we find $\sigma =\sigma \left(m\right)$ is (15) $\sigma \left(m\right)=\frac{(b_1{d}_2-b_2{d}_1)(b_1{m}^2+b_2)m}{(d_1{m}^2+d_2{)}^2},$(15) which satisfies $\sigma (-m)=-\sigma \left(m\right)$ and is zero for $b_1{d}_2-b_2{d}_1=0$, the aforementioned special anisotropic case $b_1/d_1=b_2/d_2$. In the remaining case $b_1{d}_2\ne b_2{d}_1$, we note that $\sigma \left(m\right)$ is differentiable and $\sigma \left(m\right)\to 0$ for $m\to \pm \mathrm{\infty }$, so that it suffices to determine the extrema. The derivative of $\sigma \left(m\right)$ is given by (16) ${\sigma }^{\prime }\left(m\right)=\frac{b_1{d}_2-b_2{d}_1}{(d_1{m}^2+d_2{)}^3}(-b_1{d}_1{m}^4+3(b_1{d}_2-b_2{d}_1)m^2+b_2{d}_2),$(16) whose roots, and therefore the location of the extrema, are (17) $m_{\pm }=\pm \sqrt{\frac{3(b_1{d}_2-b_2{d}_1)}{2b_1{d}_1}+\sqrt{\frac{9(b_1{d}_2-b_2{d}_1{)}^2}{4b_1^2{d}_1^2}+\frac{b_2{d}_2}{b_1{d}_1}}}.$(17) Thus, steady explicit solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) exist for $\sigma \in \left[\sigma \right(m_{-}),\sigma (m_{+}\left)\right]$, if $b_1{d}_2-b_2{d}_1>0$, and $\sigma \in \left[\sigma \right(m_{+}),\sigma (m_{-}\left)\right]$ for the case $b_1{d}_2-b_2{d}_1<0$. We may interpret the endpoints $\sigma \left(m_{-}\right)$ and $\sigma \left(m_{+}\right)$, where the solution curves of (7(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ) with $\lambda =0$ touch each other (see figure 3(c)), as bifurcation points of explicit steady solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ).

Using $m=\tan \phi$, we equivalently obtain σ as a function of the wave vector angle φ that was used above. We plot an example of the resulting function $\sigma \left(\phi \right)$ in figure 5(b).

2.2. Stability analysis of steady solutions

We study the stability of a steady state $({\mathit{v}}_s,{\eta }_s)$ of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) via the linear operator $\mathbf{L}=\mathbf{L}({\mathit{v}}_s,{\eta }_s)$, which results from linearising (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) in $({\mathit{v}}_s,{\eta }_s)$. A spectrum of $\mathbf{L}$ with positive real part then implies that the steady solution $({\mathit{v}}_s,{\eta }_s)$ is linearly unstable. In the following, we will show that the trivial steady state $({\mathit{v}}_s,{\eta }_s)\equiv (0,0,0)$ is linearly unstable, and in certain cases even unboundedly unstable. Afterwards, we focus on the stability of non-trivial steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ). We first analyse the unbounded instability of these in the full nonlinear equations (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) with respect to solutions of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ). Moreover, we study a certain long-wavelength instability in this case and briefly consider the energy of the explicit solutions. Finally, we investigate the linear stability of all steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) with small and large amplitudes.

2.2.1. Linear stability of trivial steady state

For the trivial homogeneous steady solution $({\mathit{v}}_s,{\eta }_s)\equiv (0,0,0)$, the linearisation of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) is exactly (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) without the nonlinear terms, and the corresponding linear operator $\mathbf{L}$ is then the remaining right-hand side of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ). The spectrum of $\mathbf{L}$ in this case can be determined by the dispersion relation $d(\lambda ,\mathit{k}):=det(\lambda \mathsf{I}-\hat{\mathsf{L}})=0,$ with wave vectors $\mathit{k}=(k_1,k_2{)}^{\mathsf{T}}\in {\mathbb{R}}^2$, temporal rates $\lambda =\lambda \left(\mathit{k}\right)\in \mathbb{C}$ and $\hat{\mathsf{L}}=\hat{\mathsf{L}}\left(\mathit{k}\right)$ the Fourier transform of $\mathbf{L}$ given by $\hat{\mathsf{L}}\left(\mathit{k}\right)=\left(\begin{array}{ccc}-d_1|\mathit{k}{|}^4+b_1|\mathit{k}{|}^2& f& -\mathrm{i}g{k}_1\\ -f& -d_2|\mathit{k}{|}^4+b_2|\mathit{k}{|}^2& -\mathrm{i}g{k}_2\\ -\mathrm{i}{H}_0{k}_1& -\mathrm{i}{H}_0{k}_2& 0\end{array}\right).$ The dispersion relation is thus explicitly (18) $d(\lambda ,\mathit{k})={\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(18) with the coefficients $\begin{array}{rl}{c}_2& :=(d_1+d_2)|\mathit{k}{|}^4-(b_1+b_2)|\mathit{k}{|}^2,\\ c_1& :=(d_1|\mathit{k}{|}^4-b_1|\mathit{k}{|}^2)(d_2|\mathit{k}{|}^4-b_2|\mathit{k}{|}^2)+g{H}_0|\mathit{k}{|}^2+f^2,\\ c_0& :=-g{H}_0|\mathit{k}{|}^2\left(\right(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2\left)k_1^2\right).\end{array}$ Recall that the explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) with (6(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) solve both (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) with and without the nonlinear terms, since these terms vanish by construction of these solutions. Therefore, the wave vectors $\mathit{k}$ and growth rates λ of the explicit solutions are in fact real solutions of the dispersion relation (18(18) $d(\lambda ,\mathit{k})={\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(18) ). In other words, all these explicit solutions are real eigenmodes of $\mathbf{L}$ and the values λ defined by the dispersion relation (6(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) are the corresponding real elements in the spectrum of $\mathbf{L}$. Thus, the possible values for λ of the explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) with (6(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) directly provide part of the spectrum of $\mathbf{L}$, for instance all values of λ on the white and black curves in figure 1(c). In particular, the occurrence of positive growth rates λ in (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) implies that the trivial steady solution $({\mathit{v}}_s,{\eta }_s)\equiv (0,0,0)$ is linearly unstable with respect to these exponentially growing explicit solutions. For instance, in figure 1(c), this happens for the wave vectors on the part of the white curves within the red region. More generally, even if the white curves do not intersect the red region, we next show that the red region is filled with unstable real modes of $\mathbf{L}$.

We first note that in $\mathit{k}=(0,0{)}^{\mathsf{T}}$, the dispersion relation reduces to $d(\lambda ,(0,0{)}^{\mathsf{T}})={\lambda }^3+f^2\lambda =0,$ which gives $\lambda =0$ and $\lambda =\pm \mathrm{i}f$, all having zero real part. A subset of the unstable spectrum can be determined by the sign of $c_0=-g{H}_0|\mathit{k}{|}^2\text{β}\left(\mathit{k}\right)$, where $\text{β}\left(\mathit{k}\right)$ is the expression in brackets in the definition of $c_0$ above. The coefficient $c_0$ of the dispersion relation (18(18) $d(\lambda ,\mathit{k})={\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(18) ) is zero if and only if $\mathit{k}\in {\mathbb{R}}^2$ satisfies $\text{β}\left(\mathit{k}\right)=0$, which means that $\lambda =0$ is in the spectrum of $\mathbf{L}$ with corresponding eigenmodes having such wave vectors $\mathit{k}$. Furthermore, $c_0$ is negative if and only if $\text{β}\left(\mathit{k}\right)>0$, so according to the dispersion relation (18(18) $d(\lambda ,\mathit{k})={\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(18) ) there is at least one positive real value $\lambda >0$ for each of these wave vectors $\mathit{k}$. We notice that $\text{β}\left(\mathit{k}\right)$ is also exactly the same expression as on the right-hand side of the dispersion relation (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) or (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ), whose sign we have already analysed above. In other words, the red regions plotted, e.g. in figure 1(c), correspond to a part of the unstable spectrum of $\mathbf{L}$ which are positive real. In particular, we conclude that $({\mathit{v}}_s,{\eta }_s)\equiv (0,0,0)$ is linearly unstable for any choice of parameters with horizontal backscatter. However, the spectrum of $\mathbf{L}$ may also contain non-real unstable parts. We plot an example in figure 4(a), where the unstable region extends into the blue region of figure 1(c). This can be further studied based on the dispersion relation (18(18) $d(\lambda ,\mathit{k})={\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(18) ), but we will not do this here.

Figure 4. Information on the spectrum of the linearisation of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) in the trivial steady solution $({\mathit{v}}_s,{\eta }_s)\equiv (0,0,0)$ for parameters $d_1=1.0,d_2=1.04,b_1=1.5,b_2=2.2,f=-0.3,g=9.8,H_0=0.1$. (a) Signs of the most unstable real part of elements in the spectrum in terms of the wave vector $\mathit{k}$ of the associated eigenmodes – real part positive (red region), negative (blue region); black dot and curve correspond to steady states of (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ). Note the unstable non-real spectrum in addition to, e.g. figure 1. (b) The real part of the spectrum for $k_2=0$, showing unstable spectrum in the vicinity of the origin (Colour online).

Figure 5. Examples for unboundedly unstable explicit steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) for the same parameter values as in figure 1(c). (a) Red arcs mark unboundedly unstable steady solutions. Grey curves mark explicit time-dependent solutions with ${\alpha }_2=0$. (b) Graph of $\sigma \left(\phi \right)$. Red parts mark unboundedly unstable cases; dashed grey line marks the value of the Coriolis parameter f (Colour online).

The previous investigation regarding the instability of the trivial flow in fact shows that the explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) of the full nonlinear equations (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) with $\lambda >0$ are real unstable eigenmodes. Here, we see a specific case of what we refer to as unbounded instability: perturbations of the zero state by one such mode not only leads to infinitesimal or local growth, but to globally in time unbounded growth in the nonlinear system.

2.2.2. Unbounded and long-wavelength instability of non-trivial steady states

In the following we show that some of the steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) can be unboundedly unstable as well. We consider parameter values such that some time-dependent explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) have positive growth rate λ, as in the example of figure 1; recall that this requires at least weak anisotropy of the horizontal backscatter. As already shown, steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) exist on the whole curve defined by the dispersion relation (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) with $\lambda =0$ (see figure 5(a) as well as the black curve in figure 1(c)). Now superpositions of explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ), which have the same wave vector direction (e.g. the intersections of white or black curves with the grey line in figure 1(c)), are also explicit solutions of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ). In the case of figure 1(c), these are in particular a non-trivial steady solution $({\mathit{v}}_s,{\eta }_s)$ (white dot) and an exponentially growing solution $({\mathit{v}}_g,0)$ (white circle in red area). Any superposition $\alpha ({\mathit{v}}_s,{\eta }_s)+ϵ({\mathit{v}}_g,0)$, with arbitrary $\alpha ,ϵ\in \mathbb{R}$, is also an explicit solution of (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ); in particular, ε can be arbitrarily close to zero. For any $ϵ\ne 0$, the resulting solution is exponentially and unboundedly growing. Thus, $({\mathit{v}}_s,{\eta }_s)$ is an unboundedly unstable steady solution.

This implies the unbounded instability of the explicit steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) corresponding to wave vectors on the red arcs in figure 5(a); in figure 1(c), these are between the intersections of black and white curves. These arcs connect intersection points of the curve defined by the dispersion relation (6a(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) for $\lambda =0$, with that for time-dependent explicit solutions defined by the amplitude relation (6b(6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) ) with ${\alpha }_2=0$. The instability of the other explicit steady solutions (black regions in figure 5(a)) is not determined in this way; we discuss some cases later. However, the transition from the black to the red arcs can be associated with a long-wavelength instability (also called sideband or modulational instability). The numerical result plotted in figure 7 shows that this instability should be expected on top of already unstable spectrum.

In order to study the long-wavelength instability, we consider the wave vector angle φ and first discuss the values of $\sigma \left(\phi \right)$, for which the corresponding explicit steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) are unboundedly unstable. Recall that since time-dependent explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) require ${\alpha }_2=0$, according to (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ), these have $\sigma =f$. Thus, steady solutions with $\sigma \left(\phi \right)=f$ lie at the intersections with the curve of time-dependent solutions and all steady solutions “between” those with $\sigma \left(\phi \right)=f$ are unboundedly unstable, since those can be superposed with growing explicit solutions (cf. red regions in figure 5(a)). More precisely, there are at most four angles ${\phi }_j\in [0,2\pi )$, ordered by size, so that $\sigma \left({\phi }_j\right)=f$ (cf. figure 5(b)), and steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) whose wave vectors have angles between ${\phi }_1$ and ${\phi }_2$, or ${\phi }_3$ and ${\phi }_4$, are unboundedly unstable. Hence, a steady solution (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) is unboundedly unstable if and only if $\mathrm{s}\mathrm{g}\mathrm{n}\left(f\right)\sigma \left(\phi \right)>|f|$, with its corresponding value $\sigma \left(\phi \right)$ (see figure 5(b)).

Towards the long-wavelength instability, we parameterise the set of steady solutions by the angle φ of their wave vectors. In figure 6(a) we plot for each φ the wave vector lengths for which an explicit solution of the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) exists and whether it is steady, exponentially decaying or growing. This also readily shows admissible superpositions of explicit solutions, since these must have the same angle φ; with respect to these exponentially growing explicit solutions, we thus have stable and unstable steady solutions. The stability change occurs at the intersections of the curves of the steady and time-dependent solutions, thus providing a long-wavelength instability character at these points, since the difference of wavelength between the steady solution and the unstable mode (the Floquet–Bloch parameter) crosses zero here (see figure 6(b)). Conversely, given any small Floquet–Bloch parameter, one can find a value of φ, so that the corresponding steady solution is unstable with respect to it. See figure 6(c), where the self-intersection point at the origin shows two such points along φ.

Figure 6. Illustration of long-wavelength instabilities of explicit flows. (a and b) Wave vectors $\mathit{k}=|\mathit{k}|(\cos \phi ,\sin \phi {)}^{\mathsf{T}}$ of explicit solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) as functions of φ. Black: steady solutions ${\mathit{k}}_s\left(\phi \right)$; blue: exponentially decaying; red: exponentially growing. (c) Growth rate λ as a function of the difference of wave vector lengths from steady solution. Parameters for all three cases as in figure 1(c) (Colour online).

Figure 7. Shown are the eigenvalues with real part larger than $-0.1$ of an approximation of ${\mathbf{L}}_1$ with N = 10 wave modes, i.e. $3(2N+1{)}^2$ Fourier modes on the periodic domain $[0,2\pi /k_1]\times [0,2\pi /k_2]$, and Bloch modes from the grid with distance $\pi /4$. Parameters are as in figure 1(c) and s = 0. Amplitudes are ${\alpha }_1=1$ and ${\alpha }_2=0$, so $\sigma =f=0.3$ and the selected steady solution corresponds to the point between the red and black arcs in figure 5 with $\mathit{k}\approx (-1.4,0.35)$. In particular, the solution is already unstable at marginal instability with respect to the explicit modes.

We briefly consider some energetic aspects of the explicit solutions. The kinetic energy density of solutions to (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) is given by $\mathrm{K}\mathrm{E}(\mathit{v},\eta )=(H_0+\eta )|\mathit{v}{|}^2/2$ and the potential energy density by $\mathrm{P}\mathrm{E}(\mathit{v},\eta )=g(H_0+\eta {)}^2/2$. The superposed explicit solutions are generally of the form $\begin{array}{rl}\mathit{v}& ={\alpha }_s\cos \left({\xi }_s\right){\mathit{k}}^{\perp }+{\alpha }_n{\mathrm{e}}^{{\lambda }_{n}t}\cos \left({\xi }_n\right){\mathit{k}}^{\perp }+{\alpha }_p{\mathrm{e}}^{{\lambda }_{p}t}\cos \left({\xi }_p\right){\mathit{k}}^{\perp },\\ \eta & =\gamma \frac{f}{g}\sin {\xi }_s+c,\end{array}$ with $|\mathit{k}|=1$, phase variables ${\xi }_a={\kappa }_a\mathit{k}\mathbf{\cdot }\mathit{x}+{\tau }_a$ for $a\in \{s,n,p\}$, growth rates ${\lambda }_n<0$ and ${\lambda }_p>0$, wave numbers ${\kappa }_s,{\kappa }_n,{\kappa }_p\in \mathbb{R}$ as well as arbitrary amplitudes ${\alpha }_a\in \mathbb{R}$ and shifts ${\tau }_a,c\in \mathbb{R}$ for any $a\in \{s,n,p\}$. The corresponding terms are steady, decaying and growing explicit solutions, determined by (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) and (6(6a) $\begin{array}{rl}\lambda & =(b_1-d_1|\mathit{k}{|}^2)k_2^2+(b_2-d_2|\mathit{k}{|}^2)k_1^2,\end{array}$(6a) ) (compare with steady solutions in red region in figure 5(a) and the possible superpositions with solutions which are decaying or growing in time). The energy densities of these explicit solutions explicitly read $\begin{array}{rl}\mathrm{K}\mathrm{E}(\mathit{v},\eta )& =\frac{1}{2}(\gamma \frac{f}{g}\sin {\xi }_s+H_0+c)({\alpha }_s\cos {\xi }_s+{\alpha }_n{\mathrm{e}}^{{\lambda }_{n}t}\cos {\xi }_n+{\alpha }_p{\mathrm{e}}^{{\lambda }_{p}t}\cos {\xi }_p{)}^2,\\ \mathrm{P}\mathrm{E}(\mathit{v},\eta )& =\frac{{\gamma }^2{f}^2}{2g}{\sin }^2{\xi }_s+\gamma f(H_0+c)\sin {\xi }_s+\frac{g(H_0+c{)}^2}{2}.\end{array}$ Notably, being cubic in the sine/cosine terms, the kinetic energy in Fourier space features various diadic and triadic combinations of the wave vectors of the corresponding velocity components of the explicit solution. On the temporal side, the squared linear combination of time-independent, decaying and growing parts yields doubling and adding of the individual rates. The potential energy is ignorant to the dynamic terms, but we note the constant and $2{\kappa }_s\mathit{k}$ Fourier modes from the quadratic term.

2.2.3. Linear stability of non-trivial steady states with small and large amplitudes

We now study stability properties of steady solutions (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) for (asymptotically) small and large amplitudes, the natural asymptotic regimes for families of solutions with a free amplitude parameter. Since such linear spaces of solutions arise more broadly in incompressible fluid equations with transport nonlinearity (cf. Prugger and Rademacher 2021), and for later use in section 3, we set up the notation for the more general setting of an evolution equation with linear term $\mathbf{L}$ and bilinear nonlinearity $\mathit{B}$ given by $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathit{u}=\mathbf{L}\mathit{u}+\mathit{B}(\mathit{u},\mathit{u})+\mathbf{\nabla }p,$ with a pressure p, that is trivial for the rotating shallow water equations (1(1a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}& =-f{\mathit{v}}^{\perp }-g\mathbf{\nabla }\eta -\left(\begin{array}{cc}{d}_1{\mathrm{\Delta }}^2+b_1\mathrm{\Delta }& 0\\ 0& d_2{\mathrm{\Delta }}^2+b_2\mathrm{\Delta }\end{array}\right)\mathit{v},\end{array}$(1a) ) with $\mathit{u}=(\mathit{v},\eta )$, and otherwise will derive from the incompressibility constraint $\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0$.

We assume that there exists a family of steady-state solutions $\mathit{u}=a{\mathit{u}}_0$ with amplitude parameter $a\in \mathbb{R}$ and associated (possibly trivial) pressure $p=a{p}_0$. The spectral stability of the steady state $a{\mathit{u}}_0$ is determined by the linearised right-hand side in $a{\mathit{u}}_0$ and thus the solutions to the generalised eigenvalue problem $\lambda \mathit{U}={\mathbf{L}}_a\mathit{U}+\mathbf{\nabla }P,{\mathbf{L}}_a:=\mathbf{L}+a{\mathbf{L}}_0,{\mathbf{L}}_0:=\mathit{B}(\cdot ,{\mathit{u}}_0)+\mathit{B}({\mathit{u}}_0,\cdot ),$ with eigenvalue parameter $\lambda \in \mathbb{C}$, eigenmode $\mathit{U}=(\mathit{v},\eta )$ and P either trivial or determined by the linearised constraint $\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0$.

Since the resulting spectrum is locally uniformly continuous with respect to the parameter a, we immediately note that for $|a|\ll 1$ it is close to that for a = 0 associated to $\mathbf{L}$. Since its spectrum is unstable for the backscatter setting, as shown in the linear stability analysis of the zero state above, it follows that all the discussed explicit flows for small amplitudes inherit unstable modes of the trivial state, more so for smaller amplitudes.

Regarding large amplitudes, $|a|\gg 1$, we consider eigenvalue parameters that scale with the amplitude, i.e. $\lambda =a\tilde{\lambda }$, and set $P=a\tilde{P}$. This gives the (generalised) eigenvalue problem (19) $\tilde{\lambda }\mathit{U}=(a^{-1}\mathbf{L}+{\mathbf{L}}_0)\mathit{U}+\mathbf{\nabla }\tilde{P}.$(19) The operator of the limiting problem, as $|a|\to \mathrm{\infty }$, is ${\mathbf{L}}_0$, and again by continuity of the spectrum, its stability properties partially predict those of ${\mathbf{L}}_a$ for $|a|\gg 1$. In particular, an unstable eigenmode of ${\mathbf{L}}_0$ implies strongly unstable eigenmodes of ${\mathbf{L}}_a$ for $|a|\gg 1$, for which the growth rate $\mathrm{R}\mathrm{e}\left(\lambda \right)$ is proportional to the amplitude a of the steady solution. However, eigenmodes of ${\mathbf{L}}_a$ for which λ is not proportional to a will move to the origin in the scaled operator as $|a|\to \mathrm{\infty }$, and thus contribute to the kernel of ${\mathbf{L}}_0$. In particular, $a{\mathit{u}}_0$ may be unstable for all a even though ${\mathbf{L}}_0$ does not possess unstable spectrum. Indeed, this turns out to be the case in the present setting. This is consistent with the unstable rates λ of the explicit flows from the analysis of unbounded instability above, which are associated with unbounded growth, as these are constant with respect to a so that in the scaling of ${\mathbf{L}}_a$ satisfy $\tilde{\lambda }=\lambda /a\to 0$ as $a\to \mathrm{\infty }$.

Hence, we consider the limiting problem $\tilde{\lambda }\mathit{U}={\mathbf{L}}_0\mathit{U}+\mathbf{\nabla }\tilde{P},$ whose spectral properties do not seem to be known analytically for the explicit flows we are concerned with. Specifically, for (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) we have $\mathit{u}=(\mathit{v},\eta )$, the bilinear form is $\mathit{B}(({\mathit{v}}_1,{\eta }_1),({\mathit{v}}_2,{\eta }_2))=\left(\begin{array}{c}-({\mathit{v}}_1\mathbf{\cdot }\mathbf{\nabla }){\mathit{v}}_2\\ -({\mathit{v}}_1\mathbf{\cdot }\mathbf{\nabla }){\eta }_2-{\eta }_1\mathbf{\nabla }\mathbf{\cdot }{\mathit{v}}_2\end{array}\right),$ and the steady-state family is generated by ${\mathit{u}}_0=({\mathit{v}}_0,{\eta }_0{)}^{\mathsf{T}}$ from (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ), i.e. with $\xi =\mathit{k}\mathbf{\cdot }\mathit{x}$, ${\mathit{v}}_0=\cos \left(\xi \right){\mathit{k}}^{\perp },{\eta }_0={\alpha }_2\sin \xi +s,$ where ${\alpha }_2$ is chosen so that (6b(6b) $\begin{array}{rl}\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_1}f& =k_1{k}_2\left(\right(d_1-d_2)|\mathit{k}{|}^2+b_2-b_1),\end{array}$(6b) ) holds. We are then interested in the spectrum of the operator ${\mathbf{L}}_0$ defined by $\begin{array}{rl}{\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \eta \end{array}\right)& =-\left(\begin{array}{l}({\mathit{v}}_0\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla }){\mathit{v}}_0\\ ({\mathit{v}}_0\mathbf{\cdot }\mathbf{\nabla })\eta +(\mathit{v}\mathbf{\cdot }\mathbf{\nabla }){\eta }_0+{\eta }_0\mathbf{\nabla }\mathbf{\cdot }\mathit{v}+\eta \mathbf{\nabla }\mathbf{\cdot }{\mathit{v}}_0\end{array}\right)\\ & =-\left(\begin{array}{l}\cos \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\mathit{v}-(\mathit{v}\mathbf{\cdot }\mathit{k})\sin \left(\xi \right){\mathit{k}}^{\perp }\\ \cos \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\eta +{\alpha }_2(\mathit{v}\mathbf{\cdot }\mathit{k})\cos \xi +({\alpha }_2\sin \xi +s)(\mathbf{\nabla }\mathbf{\cdot }\mathit{v})\end{array}\right).\end{array}$ We immediately note that the kernel of ${\mathbf{L}}_0$ is infinite dimensional: any perturbation $\mathit{v},\eta$ of the same form as the steady flow ${\mathit{u}}_0$, i.e. $\mathit{v}\left(\mathit{x}\right)={\varphi }_1\left(\xi \right){\mathit{k}}^{\perp }$, $\eta \left(\mathit{x}\right)={\varphi }_2\left(\xi \right)$ with arbitrary ${\varphi }_1$ and ${\varphi }_2$, lies in the kernel, since $({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })(\mathit{v},\eta )=0$ as well as $\mathit{v}\mathbf{\cdot }\mathit{k}=0$ and $\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0$ (in fact, one can show that here $\mathit{v}$ can be an arbitrary function of ξ). Next, we show that the spectrum of ${\mathbf{L}}_0$ is purely imaginary.

First let ${\alpha }_2=s=0$, i.e. the steady state ${\mathit{u}}_0$ is in the intersection of the sets of steady and time-dependent solutions as in figure 1, so that ${\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \eta \end{array}\right)=-\left(\begin{array}{l}\cos \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\mathit{v}-(\mathit{v}\mathbf{\cdot }\mathit{k})\sin \left(\xi \right){\mathit{k}}^{\perp }\\ \cos \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\eta \end{array}\right).$ It is a diagonal operator where $\mathit{v}$ and η are decoupled. Let us change coordinates to $\xi =k_{1}x+k_{2}y$, $\zeta =-k_{2}x+k_{1}y$. Then $\mathbf{\nabla }$ becomes $(k_1{\mathrm{\partial }}_{\xi }-k_2{\mathrm{\partial }}_{\zeta },k_2{\mathrm{\partial }}_{\xi }+k_1{\mathrm{\partial }}_{\zeta })$ so that ${\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla }$ turns into $K{\mathrm{\partial }}_{\zeta }$, where $K:=k_1^2+k_2^2$. We thus obtain the operator ${\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \eta \end{array}\right)=-\left(\begin{array}{l}\cos \left(\xi \right)K{\mathrm{\partial }}_{\zeta }\mathit{v}-(\mathit{v}\mathbf{\cdot }\mathit{k})\sin \left(\xi \right){\mathit{k}}^{\perp }\\ \cos \left(\xi \right)K{\mathrm{\partial }}_{\zeta }\eta \end{array}\right),$ whose Fourier transform with respect to ζ with wave number parameter ϑ reads $\begin{array}{rl}{\hat{\mathbf{L}}}_0\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right)& =-\left(\begin{array}{l}\mathrm{i}\cos \left(\xi \right)K\vartheta \hat{\mathit{v}}-(\hat{\mathit{v}}\mathbf{\cdot }\mathit{k})\sin \left(\xi \right){\mathit{k}}^{\perp }\\ \mathrm{i}\cos \left(\xi \right)K\vartheta \hat{\eta }\end{array}\right)\\ & =-\left(\begin{array}{cc}\mathrm{i}\cos \left(\xi \right)K\vartheta \mathsf{I}-\sin \left(\xi \right)\mathsf{A}\left(\mathit{k}\right)& 0\\ 0& \mathrm{i}\cos \left(\xi \right)K\vartheta \end{array}\right)\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right),\end{array}$ where $\mathsf{A}\left(\mathit{k}\right)=\left(\begin{array}{cc}-k_1{k}_2& -k_2^2\\ k_1^2& k_1{k}_2\end{array}\right)$.

The lower right entry, which corresponds to η, is a multiplication operator by $\mathrm{i}\cos \left(\xi \right)K\vartheta$ and so its spectrum is the range of this function, which is $\mathrm{i}K\vartheta [-1,1]\subset \mathrm{i}\mathbb{R}$. Since this multiplication operator appears in the upper left entry as multiplying the identity, which commutes with any matrix, $\mathsf{A}\left(\mathit{k}\right)$ can be brought to normal form. This features a double zero eigenvalue so that the operator on the upper left block possesses purely imaginary spectrum. In particular, the spectrum is neutrally stable.

For ${\alpha }_2,s\ne 0$ and writing $\mathit{v}=(u,v)$, we analogously obtain the transformed operator ${\hat{\mathbf{L}}}_0\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right)=-\left(\begin{array}{l}\mathrm{i}\cos \left(\xi \right)K\vartheta \hat{\mathit{v}}-(\hat{\mathit{v}}\mathbf{\cdot }\mathit{k})\sin \left(\xi \right){\mathit{k}}^{\perp }\\ \mathrm{i}\cos \left(\xi \right)K\vartheta \hat{\eta }+\mathbf{B}\hat{\mathit{v}}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{A}}_1& 0\\ \mathbf{B}& {\mathbf{A}}_2\end{array}\right)\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right),$ with $\begin{array}{rl}{\mathbf{A}}_1& :=-\mathrm{i}\cos \left(\xi \right)K\vartheta \mathsf{I}+\sin \left(\xi \right)\mathsf{A}\left(\mathit{k}\right),\\ {\mathbf{A}}_2& :=-\mathrm{i}\cos \left(\xi \right)K\vartheta ,\\ \mathbf{B}\hat{\mathit{v}}& :=-{\alpha }_2(\hat{\mathit{v}}\mathbf{\cdot }\mathit{k})\cos \xi +({\alpha }_2\sin \xi +s)(k_1{\mathrm{\partial }}_{\xi }\hat{u}+k_2{\mathrm{\partial }}_{\xi }\hat{v}+\mathrm{i}\vartheta (k_1\hat{v}-k_2\hat{u})).\end{array}$ If λ lies in the resolvent set of both operators on the diagonal ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$ (by the above this includes any non-purely imaginary value), then λ is also in the resolvent set of the present ${\hat{\mathbf{L}}}_0$, since $({\hat{\mathbf{L}}}_0-\lambda {)}^{-1}\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right)=\left(\begin{array}{cc}({\mathbf{A}}_1-\lambda \mathsf{I}{)}^{-1}& 0\\ -({\mathbf{A}}_2-\lambda {)}^{-1}\mathbf{B}({\mathbf{A}}_1-\lambda \mathsf{I}{)}^{-1}& ({\mathbf{A}}_2-\lambda {)}^{-1}\end{array}\right)\left(\begin{array}{c}\hat{\mathit{v}}\\ \hat{\eta }\end{array}\right).$ Hence, as claimed, the asymptotic operator possesses marginally stable spectrum and we cannot immediately infer in/stability information for large amplitudes. However, numerical computations based on truncated Fourier series suggest that the spectrum is in fact rather strongly unstable, cf. figure 7.

3. Rotating Boussinesq equations with backscatter

We now turn to the study of various explicit solutions in the rotating Boussinesq equations augmented with backscatter (2(2a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2,-\nu )\mathrm{\Delta }\mathit{v},\end{array}$(2a) ). To ease notation, we write these in the form (20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) (20b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(20b) (20c) $\begin{array}{rl}\frac{\mathrm{\partial }\mathrm{b}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathrm{b}+N^2{v}_3& =\mu \mathrm{\Delta }\mathrm{b},\end{array}$(20c) where we focus on horizontal backscatter $d_1,d_2,b_1,b_2>0$ with usual viscosity vertically, $d_3=0,b_3=-\nu \le 0$, and stable stratification $N^2>0$. We emphasise that in contrast to the shallow water equations, the backscatter can be horizontally isotropic for the phenomena we find here, in particular for unbounded instability in the classes of flows we consider. Not surprisingly, the additional parameters of anisotropy often give more freedom of choice. A subtlety is that for one type of MGWs discussed in section 3.2.3, unbounded growth occurs under stable stratification only in the presence of at least some anisotropy. For comparison and illustration, we also discuss briefly the usual horizontal viscous or inviscid cases $d_1=d_2=0$, $b_1,b_2\le 0$, unstable stratification $N^2<0$, and artificial vertical backscatter $d_3,b_3>0$. As in section 2, we are especially interested in the parameter relations and stability properties of steady solutions, in particular unbounded instability, as well as in unboundedly growing explicit solutions. We first investigate the horizontal flows, which are comparable with the explicit solutions of the shallow water equations in section 2, but are less restricted and have additional properties in this case here. Afterwards, we analyse other explicit solutions with vertical structure and coupled buoyancy.

3.1. Horizontal flow and decoupled system

In order to compare with the results of the rotating shallow water equations with backscatter, we consider here the barotropic case with horizontal velocity field. We therefore choose a velocity field $\mathit{v}$ that is independent of the vertical coordinate z and has $v_3\equiv 0$, as well as a horizontally independent buoyancy $\mathrm{b}=\mathrm{b}(t,z)$. This ansatz yields the reduced equations (21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) (21b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(21b) (21c) $\begin{array}{rl}\frac{\mathrm{\partial }\mathrm{b}}{\mathrm{\partial }t}& =\mu \frac{{\mathrm{\partial }}^2\mathrm{b}}{\mathrm{\partial }{z}^2},\end{array}$(21c) with gradient and Laplacian for the horizontal directions $\mathit{x}=(x,y{)}^{\mathsf{T}}$, $\mathit{v}=\mathit{v}(t,\mathit{x})\in {\mathbb{R}}^2$, $\mathrm{b}=\mathrm{b}(t,z)$ and $p=\tilde{p}+B(t,z)$, where $\mathrm{\partial }/\mathrm{\partial }z\left[B\right(t,z\left)\right]=\mathrm{b}(t,z)$. Here the buoyancy is decoupled from the velocity field and determined by the linear heat equation (21c(21c) $\begin{array}{rl}\frac{\mathrm{\partial }\mathrm{b}}{\mathrm{\partial }t}& =\mu \frac{{\mathrm{\partial }}^2\mathrm{b}}{\mathrm{\partial }{z}^2},\end{array}$(21c) ); on the idealised whole space this can be readily solved by Fourier transform.

Regarding the momentum equations, compared with the shallow water equations we may view the equation for the fluid depth (1b(1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) ) to be replaced by the incompressibility condition (21b(21b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(21b) ). This is less restrictive and admits a larger set of explicit flow solutions as discussed in Prugger and Rademacher (2021) for the setting without backscatter. In particular, the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) for $\mathit{v}$ satisfies (21b(21b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(21b) ) and can readily be adjusted to solutions of (21a(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ). However, there are no additional a priori constraints for the pressure akin to condition (1b(1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) ) or (6c(6c) $\begin{array}{rl}{\alpha }_2\cdot \lambda & =0.\end{array}$(6c) ), so that the resulting pressure can be exponentially decaying or growing along with the velocity field. Moreover, linear combinations of any of these solutions with the same wave vector direction but different wavelength (any wave vector on a whole ray like in figure 1(c)) also yield explicit solutions of (21(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ).

The resulting set of explicit solutions of (21a(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ) for which the nonlinear terms vanish can be identified by the following ansatz for wave shape ψ and pressure profile ϕ (22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) where $\mathit{k}\in {\mathbb{R}}^2$ and without loss of generality $|\mathit{k}|=1$ by the freedom in choosing ψ and ϕ. Substitution into (21a(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ) gives the linear equations for ψ and ϕ, (23a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(23a) (23b) $\begin{array}{rl}\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\xi }& =\frac{k_1{k}_2}{f}\left(\right(d_1-d_2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}+(b_1-b_2\left)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2}\right)+\psi .\end{array}$(23b) For the Boussinesq equations with viscosity, instead of the backscatter terms similar equations arise (cf. Prugger and Rademacher 2021). However, in that case the pressure gradient fully compensates the buoyancy and the Coriolis term in the equations, which makes the velocity field geostrophically balanced. In contrast, in the present case, equation (23b(23b) $\begin{array}{rl}\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\xi }& =\frac{k_1{k}_2}{f}\left(\right(d_1-d_2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}+(b_1-b_2\left)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2}\right)+\psi .\end{array}$(23b) ) for the pressure shape ϕ allows the pressure gradient to not only compensate the full Coriolis term but also part of the backscatter terms. In particular, the velocity field (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) with (23(23a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(23a) ) is in general not geostrophically balanced.

3.1.1. Superposition principles

The general wave shape ψ in (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) also contains the superpositions of arbitrary many sinusoidal waves in the same wave vector direction $\mathit{k}$ and any wave number $|\mathit{k}|$. This is possible in the Boussinesq equations, since $\mathbf{\nabla }p$ in the momentum equation (20a(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ) is not further constrained, unlike (1b(1b) $\begin{array}{rl}\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\eta & =-(H_0+\eta )\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{v}\right),\end{array}$(1b) ) for $\mathbf{\nabla }\eta$. It is also possible to superpose by integrating over the wave numbers in the same wave vector direction.

The structure of the Boussinesq equations admits superposing $\mathit{v}$ in the form (5(5) $\mathit{v}={\alpha }_1{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau ){\mathit{k}}^{\perp },\eta ={\alpha }_2\frac{f}{g}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau )+s,$(5) ) in another way, namely with different wave vector directions, but the same wave number (cf. Prugger and Rademacher 2021). For the decoupled momentum equation (21a(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ), and finite superposition of N waves with arbitary $N\in \mathbb{N}$, the resulting superposed sinusoidal explicit solutions of a form similar to (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) are given by (24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) (24b) $\begin{array}{rl}\tilde{p}(t,\mathit{x})& =-\sum _{i=1}^N\sum _{j=i+1}^N{\mathrm{e}}^{({\lambda }_i+{\lambda }_j)t}(({\mathit{k}}_i\mathbf{\cdot }{\mathit{k}}_j){\psi }_i{\psi }_j+s^2\frac{\mathrm{\partial }{\psi }_i}{\mathrm{\partial }\xi }\frac{\mathrm{\partial }{\psi }_j}{\mathrm{\partial }\xi })-f\sum _{i=1}^N\frac{{\gamma }_i}{{\alpha }_i}{\mathrm{e}}^{{\lambda }_{i}t}\frac{\mathrm{\partial }{\psi }_i}{\mathrm{\partial }\xi },\end{array}$(24b) for any fixed s>0 and arbitrary ${\alpha }_i\in \mathbb{R}\setminus \left\{0\right\}$, ${\tau }_i\in \mathbb{R}$ and ${\mathit{k}}_i=(k_{i,1},k_{i,2}{)}^{\mathsf{T}}\in {\mathbb{R}}^2$ with $|{\mathit{k}}_i|=s$ for any $1\le i\le N$. Here, each ${\lambda }_i$ and ${\gamma }_i$ is defined by the dispersion and amplitude relations (25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) (25b) $\begin{array}{rl}f\frac{{\gamma }_i-{\alpha }_i}{{\alpha }_i}& =((d_1-d_2)s^2+b_2-b_1)k_{i,1}{k}_{i,2},\end{array}$(25b) in order to solve (21a(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ). Since each wave in (24a(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) is divergence free, the whole superposed velocity $\mathit{v}$ solves (21b(21b) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{v}& =0,\end{array}$(21b) ) and thus is an explicit solution of (21(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ).

It is also possible to superpose explicit solutions of the form (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) by integrating over the whole circle $S_s:=\{\mathit{k}\in {\mathbb{R}}^2\mid |\mathit{k}|=s\}$ for any fixed s>0. The exact form is then (26a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\int }_{S_s}\alpha \left(\mathit{k}\right){\mathrm{e}}^{\lambda \left(\mathit{k}\right)t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau \left(\mathit{k}\right)){\mathit{k}}^{\perp }\mathrm{d}\mathit{k},\end{array}$(26a) (26b) $\begin{array}{rl}\tilde{p}(t,\mathit{x})& =-s^4{\int }_0^{2\pi }{\int }_{{\phi }_1}^{2\pi }{\alpha }_1{\alpha }_2{\mathrm{e}}^{({\lambda }_1+{\lambda }_2)t}(\cos {\xi }_1\cos {\xi }_2+\cos ({\phi }_1-{\phi }_2)\sin {\xi }_1\sin {\xi }_2)\mathrm{d}{\phi }_2\mathrm{d}{\phi }_1\\ & -f{\int }_{S_s}\gamma \left(\mathit{k}\right){\mathrm{e}}^{\lambda \left(\mathit{k}\right)t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau \left(\mathit{k}\right))\mathrm{d}\mathit{k},\end{array}$(26b) where, for $i=1,2$, we set ${\alpha }_i:=\alpha \left({\mathit{k}}_i\right)$ with ${\mathit{k}}_i:=s(\cos {\phi }_i,\sin {\phi }_i{)}^{\mathsf{T}}$, ${\lambda }_i:=\lambda \left({\mathit{k}}_i\right)$ and ${\xi }_i:={\mathit{k}}_i\mathbf{\cdot }\mathit{x}+\tau \left({\mathit{k}}_i\right)$, for all $0\le {\phi }_i<2\pi$. Sufficient for the convergence of the integrals is $\alpha \in L^1\left(S_s\right)$, $\tau \in L^{\mathrm{\infty }}\left(S_s\right)$ and for almost all $\mathit{k}\in S_s$ we require (27a) $\begin{array}{rl}\alpha \left(\mathit{k}\right)\lambda \left(\mathit{k}\right)& =\alpha \left(\mathit{k}\right)((b_1-d_1{s}^2)k_2^2+(b_2-d_2{s}^2)k_1^2),\end{array}$(27a) (27b) $\begin{array}{rl}f(\gamma \left(\mathit{k}\right)-\alpha \left(\mathit{k}\right))& =\alpha \left(\mathit{k}\right)k_1{k}_2((d_1-d_2)s^2+b_2-b_1),\end{array}$(27b) corresponding to (25(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) if $\alpha \left(\mathit{k}\right)\ne 0$.

The explicit solutions (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with (25(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) differ from (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) with (23(23a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(23a) ), as well as the other explicit solutions before, not only by the structure of the superposition but also by the resulting nonlinear terms, which are not vanishing. Due to the special structure of $\mathit{v}$ in (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) and $|{\mathit{k}}_i|=s$ for all $1\le i\le N$, the nonlinear terms form a gradient that can be fully compensated by the pressure gradient in the momentum equation; this gives the first sum in (24b(24b) $\begin{array}{rl}\tilde{p}(t,\mathit{x})& =-\sum _{i=1}^N\sum _{j=i+1}^N{\mathrm{e}}^{({\lambda }_i+{\lambda }_j)t}(({\mathit{k}}_i\mathbf{\cdot }{\mathit{k}}_j){\psi }_i{\psi }_j+s^2\frac{\mathrm{\partial }{\psi }_i}{\mathrm{\partial }\xi }\frac{\mathrm{\partial }{\psi }_j}{\mathrm{\partial }\xi })-f\sum _{i=1}^N\frac{{\gamma }_i}{{\alpha }_i}{\mathrm{e}}^{{\lambda }_{i}t}\frac{\mathrm{\partial }{\psi }_i}{\mathrm{\partial }\xi },\end{array}$(24b) ). The same holds for the solutions (26(26a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\int }_{S_s}\alpha \left(\mathit{k}\right){\mathrm{e}}^{\lambda \left(\mathit{k}\right)t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau \left(\mathit{k}\right)){\mathit{k}}^{\perp }\mathrm{d}\mathit{k},\end{array}$(26a) ) with (27(27a) $\begin{array}{rl}\alpha \left(\mathit{k}\right)\lambda \left(\mathit{k}\right)& =\alpha \left(\mathit{k}\right)((b_1-d_1{s}^2)k_2^2+(b_2-d_2{s}^2)k_1^2),\end{array}$(27a) ) correspondingly. We refer to Prugger and Rademacher (2021) for further discussion and literature references for explicit solutions with gradient nonlinearities without backscatter.

Comparing the explicit solutions (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with (25(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ), as well as (26(26a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\int }_{S_s}\alpha \left(\mathit{k}\right){\mathrm{e}}^{\lambda \left(\mathit{k}\right)t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}+\tau \left(\mathit{k}\right)){\mathit{k}}^{\perp }\mathrm{d}\mathit{k},\end{array}$(26a) ) with (27(27a) $\begin{array}{rl}\alpha \left(\mathit{k}\right)\lambda \left(\mathit{k}\right)& =\alpha \left(\mathit{k}\right)((b_1-d_1{s}^2)k_2^2+(b_2-d_2{s}^2)k_1^2),\end{array}$(27a) ), with those without backscatter we notice two major differences: first, in the present case the amplitudes of the pressure can be different from those of the velocity; the conditions on the amplitudes are given in (25b(25b) $\begin{array}{rl}f\frac{{\gamma }_i-{\alpha }_i}{{\alpha }_i}& =((d_1-d_2)s^2+b_2-b_1)k_{i,1}{k}_{i,2},\end{array}$(25b) ) and (27b), respectively. The reason is that the pressure gradient in the momentum equation can additionally compensate a part of the backscatter term as well. Second, the growth rates ${\lambda }_i$ and $\lambda \left(\mathit{k}\right)$ can be different, so that each wave is decaying or growing differently. Both differences require anisotropy in the backscatter of the momentum equation, i.e. $d_1\ne d_2$ or $b_1\ne b_2$. Indeed, the usual viscosity is isotropic in this sense.

In conclusion, the above constructions of explicit solutions can be viewed as superposition principles for the (nonlinear) Boussinesq equations in this setting: (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) expresses a radial superposition principle of flows in the same wave vector direction, and (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) an angular superposition principle of flows on the same scale. Superpositions of plane waves with different wave vector directions and scales are, in general, not giving solutions.

3.1.2. Unbounded instability of steady states

Analogous to section 2.2.2, the possible superpositions of explicit solutions, as contained in (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ), imply linear subspaces with linear dynamics, in particular unbounded exponential growth of perturbations. Compared with the rotating shallow water equations, restrictions on wave vectors are absent, and in this section we discuss implications for (in)stability of steady solutions, i.e. those with zero growth rate ${\lambda }_i$ in the dispersion relation (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ).

Since $b_1,b_2,d_1,d_2>0$ in (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ), we have ${\lambda }_i>0$ for any sufficiently small wave numbers s and ${\lambda }_i<0$ for any sufficiently large ones. In particular, the trivial flow with $\mathit{v}\equiv 0$ is unstable with exponential unbounded growth with respect to any small s. More importantly, the above radial superposition principle immediately implies that the same holds for every single-wave steady solution (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ): such steady explicit solutions also arise from (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with N = 1 and ${\lambda }_1=0$ in (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ). Since the definition of ${\lambda }_i$ in (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) is exactly the same as in the dispersion relation (7a(7a) $\begin{array}{rl}\lambda & =-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2,\end{array}$(7a) ), we can use here the results from section 2.1.2 about the growth rate λ as well. Hence, for any $b_1,b_2,d_1,d_2>0$, the set of solutions with ${\lambda }_1=0$ forms a simple closed curve around the origin in the wave vector space that is symmetric with respect to axis reflections and whose interior is star shaped. Furthermore, λ is positive in the interior of this closed curve, except $\lambda =0$ at the origin $\mathit{k}=(0,0{)}^{\mathsf{T}}$, and λ is negative outside the closed curve. Thus, single-wave steady solutions (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) exist in any direction and can be radially superposed with explicit solutions with any smaller wave numbers s, which makes them unboundedly unstable.

However, in the case $b_1/d_1\ne b_2/d_2$, for fixed s there are up to four wave vectors for which ${\lambda }_i=0$; to see this note that for ${\lambda }_i=0$ the dispersion relation (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) is linear in ${\cos }^2\theta$ for ${\mathit{k}}_i=s(\cos \theta ,\sqrt{1-{\cos }^2\theta }{)}^{\mathsf{T}}$ in the first quadrant. Thus, due to the reflection symmetry, there is at most one solution in each quadrant of the wave vector plane. Because of the symmetry, the steady states of the form (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) can consist of (at most) two different wave vector directions and for those we cannot infer instability by the radial superposition principle. In that case, if $b_1/d_1\ne b_2/d_2$ we next appeal to the angular superposition principle. First, we note that for $b_1/d_1\ne b_2/d_2$, there is a unique (up to reversing orientation) longest wave vector ${\mathit{k}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ with length $s_{\mathrm{m}\mathrm{a}\mathrm{x}}$, so that (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) is satisfied with ${\lambda }_i=0$. See figure 8(a) for a typical example. Indeed, due to the aformentioned structures, ${\mathit{k}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is along an axis, though the set of $\mathit{k}$ with $\lambda =0$ need not be convex. Hence, any steady superposed solution (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with $s=s_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is built from $\pm {\mathit{k}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, which lie on the same line in wave vector space. Thus, we can use the radial superposition principle with the exponentially growing explicit solutions of smaller wave numbers, which leads to unbounded growth with respect to modes on any larger scale. The same applies for steady solutions with minimal wave vector length $s_{\mathrm{m}\mathrm{i}\mathrm{n}}$.

Figure 8. Signs of λ as defined in the dispersion relation (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) (red: $\lambda >0$, blue: $\lambda <0$, black: $\lambda =0$) and possible wave vectors for (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) for a fixed wave number s>0 (white circles). White dots mark wave vectors of corresponding steady solutions. Fixed parameters are $d_1=1.0,d_2=1.04,b_1=1.1,b_2=2.2$, i.e. $b_1/d_1\ne b_2/d_2$ (Colour online).

Second, we consider steady superposed solutions (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with $s_{\mathrm{m}\mathrm{i}\mathrm{n}}<s<s_{\mathrm{m}\mathrm{a}\mathrm{x}}$ that can be built from two different directions, cf. the white dots in figure 8(b). Here, we apply the angular superposition principle and superpose with any explicit solution (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) on the same scale, i.e. whose wave vector has the same length s, cf. the white circle in figure 8(b). Since for some wave vectors the corresponding λ defined by the dispersion relation (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) is positive for the length $s<s_{\mathrm{m}\mathrm{a}\mathrm{x}}$, at least for the wave vector $s/s_{\mathrm{m}\mathrm{a}\mathrm{x}}{\mathit{k}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, we again have unbounded instability with respect to a range of modes, here on the same scale. We note that for $s=s_{\mathrm{m}\mathrm{i}\mathrm{n}}$, radial and angular superposition both give ranges of modes with unbounded growth.

In case $b_1/d_1=b_2/d_2$, and thus for the isotropic case as well, we cannot infer the instability of steady superposed solutions (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ), i.e. steady solutions not of single-wave type, using the above superposition principles, since the wave vectors of non-trivial steady solutions $\{\mathit{k}\in {\mathbb{R}}^2\mid \lambda (\mathit{k})=0\}$, with $\lambda \left(\mathit{k}\right)$ from (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ), form a circle with radius $\sqrt{b_1/d_1}$, which means they all have the same length. Thus, explicit solutions (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ) with $s=\sqrt{b_1/d_1}$ consist only of steady solutions; there are also no steady solutions of the form (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) for other wave numbers. However, in some cases, unbounded instability still follows from angular superposition with unboundedly growing parallel flows, cf. section 3.2.1.

3.1.3. Linear stability of steady states with small and large amplitudes

After the investigation of unbounded instability of steady states, we now turn to the linear stability of steady solutions with small and large amplitudes, analogous to section 2.2.3. Concerning small amplitudes, we are naturally led to linear and spectral stability of the zero state, as for the rotating shallow water equations. Instead of analysing the full spectrum of the linearisation of (21(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ) in the trivial steady flow $\mathit{v}\equiv 0$, here we restrict attention to instability with respect to eigenmodes of the form of the horizontal flows (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ), which means solutions to (23a(23a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(23a) ). The Fourier transform of (23a(23a) $\begin{array}{rl}\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}& =-(d_1{k}_2^2+d_2{k}_1^2)\frac{{\mathrm{\partial }}^4\psi }{\mathrm{\partial }{\xi }^4}-(b_1{k}_2^2+b_2{k}_1^2)\frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{\xi }^2},\end{array}$(23a) ) yields the (linear) dispersion relation for perturbation wave vector $\mathit{k}\in {\mathbb{R}}^2$ and temporal rate $\lambda \in \mathbb{C}$, $d_{\psi }(\lambda ,\mathit{k})=-(d_1{k}_2^2+d_2{k}_1^2)|\mathit{k}{|}^2+b_1{k}_2^2+b_2{k}_1^2-\lambda =0,$ which is of course equivalent to (25a(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) with $\lambda ={\lambda }_i$ and $\mathit{k}={\mathit{k}}_i$. The above discussion for steady states built from single direction wave vectors $\pm \mathit{k}$ implies that the spectrum of the linearisation in such horizontal flows is at least as unstable as that of the zero state in this wave vector direction, since the spectra contain the growth rates of the corresponding explicit solutions in this direction. Of course the result is much stronger in that these modes actually grow unboundedly in the nonlinear Boussinesq equations. In contrast, for the linearisation of (21(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ) in steady multi-mode horizontal flow, similar comparison of its spectrum with that of the zero state holds, but with different wave vector directions and the same wave number $|\mathit{k}|=s$.

However, analogous to section 2.2.3, any (superposed) steady horizontal flow inherit the instability of any unstable mode in the dispersion relation $d_{\psi }(\lambda ,\mathit{k})=0$ of the zero state for sufficiently small amplitudes $0<|{\alpha }_i|\ll 1$, $1\le i\le N$, though the growth induced by such modes may be bounded. Recall that the explicit solutions of (21(21a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{v}}^{\perp }+\mathbf{\nabla }\tilde{p}& =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(d_1\mathrm{\Delta }+b_1,d_2\mathrm{\Delta }+b_2)\mathrm{\Delta }\mathit{v},\end{array}$(21a) ) also satisfy the full rotating Boussinesq equations (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ), which admit modes that have vertical structures and are coupled with buoyancy. In section 3.2, we will discuss such explicit solutions of (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ), which also satisfy these equations without the nonlinear terms. In particular, unboundedly growing flows of this type provide additional explicit unstable modes of the zero state $\mathit{v}\equiv 0$, which – in contrast to the horizontal flows – are also influenced by the Brunt–Väisälä frequency $N^2$, thermal diffusivity μ and the vertical viscosity. In addition, these imply linear instability modes of horizontal flows with sufficiently small amplitudes $0<|{\alpha }_i|\ll 1$, $1\le i\le N$.

As to large amplitude(s), where $1\ll |{\alpha }_i|$ for at least one i, we first note that, in the notation of section 2.2.3 and with $\mathit{u}=(\mathit{v},\mathrm{b})$, the bilinear form for the present case reads $\mathit{B}(({\mathit{v}}_1,{\mathrm{b}}_1),({\mathit{v}}_2,{\mathrm{b}}_2))=-\left(\begin{array}{c}({\mathit{v}}_1\mathbf{\cdot }\mathbf{\nabla }){\mathit{v}}_2\\ ({\mathit{v}}_1\mathbf{\cdot }\mathbf{\nabla }){\mathrm{b}}_2\end{array}\right).$ The steady-state family $\mathit{u}=a{\mathit{u}}_0$ in this case has ${\mathit{u}}_0=({\mathit{v}}_0,{\mathrm{b}}_0{)}^{\mathsf{T}}$ and ${\mathit{v}}_0=\sum _{i=1}^N{\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i){\mathit{k}}_i^{\perp },$ where we consider in this section only horizontal (wave) vectors $\mathit{k}:=(\hat{\mathit{k}},0{)}^{\mathsf{T}}$, ${\mathit{k}}^{\perp }:=({\hat{\mathit{k}}}^{\perp },0{)}^{\mathsf{T}}$ for any $\hat{\mathit{k}}\in {\mathbb{R}}^2$. Here, the third component of ${\mathit{v}}_0$ vanishes, and ${\mathrm{b}}_0$ is an arbitrary constant solving (21c(21c) $\begin{array}{rl}\frac{\mathrm{\partial }\mathrm{b}}{\mathrm{\partial }t}& =\mu \frac{{\mathrm{\partial }}^2\mathrm{b}}{\mathrm{\partial }{z}^2},\end{array}$(21c) ). In order to locate strongly unstable modes, whose growth rates $\lambda =a\tilde{\lambda }$ are proportional to the amplitude parameter a, we are concerned with the generalised eigenvalue problem (19(19) $\tilde{\lambda }\mathit{U}=(a^{-1}\mathbf{L}+{\mathbf{L}}_0)\mathit{U}+\mathbf{\nabla }\tilde{P}.$(19) ), as $|a|\to \mathrm{\infty }$, which here reads (28) $\tilde{\lambda }\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)={\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)+\left(\begin{array}{c}\mathbf{\nabla }\overline{p}\\ 0\end{array}\right),\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0,$(28) with $p=a\overline{p}$ and for the operator ${\mathbf{L}}_0$ defined by ${\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)=-\left(\begin{array}{c}({\mathit{v}}_0\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla }){\mathit{v}}_0\\ ({\mathit{v}}_0\mathbf{\cdot }\mathbf{\nabla })\mathrm{b}\end{array}\right).$ In order to simplify and illustrate the main finding, we investigate the stability of a certain superposed steady solution and reduce to only two modes, i.e. ${\alpha }_i=0$ for i>2, translate so that ${\tau }_1={\tau }_2=0$, and set $\xi =\mathit{k}\mathbf{\cdot }\mathit{x}$, $\zeta ={\mathit{k}}^{\perp }\mathbf{\cdot }\mathit{x}$ so that with a certain wave vector $\mathit{k}$ ${\mathit{v}}_0={\alpha }_1\sin \left(\xi \right){\mathit{k}}^{\perp }+{\alpha }_2\sin \left(\zeta \right)\mathit{k}.$ We note that even in the anisotropic case, there is at least one wave vector $\mathit{k}$, such that both $\mathit{k}$ and ${\mathit{k}}^{\perp }$ correspond to a steady mode. We omit the full proof and instead explain the existence of such $\mathit{k}$ based on figure 8. We start with superposed steady solutions with wave number $s=s_{max}$ as in figure 8(a). Reducing s towards the wave number $s_{min}$ as in figure 8(b), there is an intermediate value of s such that a wave vector $\mathit{k}$ with $|\mathit{k}|=s$ exists, for which $\mathit{k}$ and ${\mathit{k}}^{\perp }$ each correspond to a single mode steady solution. This is ensured by the symmetry of the curve defined by (25(25a) $\begin{array}{rl}{\lambda }_i& =(b_1-d_1{s}^2)k_{i,2}^2+(b_2-d_2{s}^2)k_{i,1}^2,\end{array}$(25a) ) with ${\lambda }_i=0$. We then obtain $\begin{array}{rl}{\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)& =-\left(\begin{array}{c}{\alpha }_1\sin \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+{\alpha }_2\sin \left(\zeta \right)(\mathit{k}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+{\alpha }_1\cos \left(\xi \right)(\mathit{v}\mathbf{\cdot }\mathit{k}){\mathit{k}}^{\perp }\\ +{\alpha }_2\cos \left(\zeta \right)(\mathit{v}\mathbf{\cdot }{\mathit{k}}^{\perp })\mathit{k}\\ {\alpha }_1\sin \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })\mathrm{b}+{\alpha }_2\sin \left(\zeta \right)(\mathit{k}\mathbf{\cdot }\mathbf{\nabla })\mathrm{b}\end{array}\right)\\ & =-\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}({\mathbf{L}}_1,{\mathbf{L}}_2)\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right),\end{array}$ with a $4\times 4$ block diagonal matrix operator in which ${\mathbf{L}}_1:=\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}({\mathbf{L}}_3,{\mathbf{L}}_2)$ is a $3\times 3$ block matrix operator build from the $2\times 2$ matrix operator ${\mathbf{L}}_3={\mathbf{L}}_2\mathsf{I}+{\alpha }_1\cos \left(\xi \right)\mathsf{A}\left(\mathit{k}\right)-{\alpha }_2\cos \left(\zeta \right){\mathsf{A}}^{\mathsf{T}}\left(\mathit{k}\right),$ with $\mathsf{A}\left(\mathit{k}\right)$ as in section 2.2.3 (for which ${\mathsf{A}}^{\mathsf{T}}\left(\mathit{k}\right)=-\mathsf{A}\left({\mathit{k}}^{\perp }\right)$) and the linear operator ${\mathbf{L}}_2={\alpha }_1\sin \left(\xi \right)({\mathit{k}}^{\perp }\mathbf{\cdot }\mathbf{\nabla })+{\alpha }_2\sin \left(\zeta \right)(\mathit{k}\mathbf{\cdot }\mathbf{\nabla }).$ Taking the divergence of (28(28) $\tilde{\lambda }\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)={\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)+\left(\begin{array}{c}\mathbf{\nabla }\overline{p}\\ 0\end{array}\right),\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0,$(28) ) gives, using $\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0$, the linear pressure Poisson equation $\mathrm{\Delta }\overline{p}=\mathbf{\nabla }\mathbf{\cdot }\left({\mathbf{L}}_1\mathit{v}\right)$. We denote the solution as $\overline{p}={\mathrm{\Delta }}^{-1}\mathbf{\nabla }\mathbf{\cdot }\left({\mathbf{L}}_1\mathit{v}\right)$ and substitution into (28(28) $\tilde{\lambda }\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)={\mathbf{L}}_0\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)+\left(\begin{array}{c}\mathbf{\nabla }\overline{p}\\ 0\end{array}\right),\mathbf{\nabla }\mathbf{\cdot }\mathit{v}=0,$(28) ) yields the eigenvalue problem (29) $\tilde{\lambda }\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)=({\mathbf{L}}_0+\left(\begin{array}{cc}\mathbf{\nabla }{\mathrm{\Delta }}^{-1}\mathbf{\nabla }\mathbf{\cdot }({\mathbf{L}}_1\circ )& 0\\ 0& 0\end{array}\right))\left(\begin{array}{c}\mathit{v}\\ \mathrm{b}\end{array}\right)$(29) in which $\circ$ denotes the slot for $\mathit{v}$. The resulting operator on the right-hand side features a diagonal block structure such that the spectrum is the union of the spectra of $-{\mathbf{L}}_2$ and $-{\tilde{\mathit{L}}}_1$ defined by ${\tilde{\mathit{L}}}_1:={\mathbf{L}}_1-\mathbf{\nabla }{\mathrm{\Delta }}^{-1}\mathbf{\nabla }\mathbf{\cdot }({\mathbf{L}}_1\circ ).$ Since ${\mathbf{L}}_2$ is skew-adjoint ($\mathrm{i}{\mathbf{L}}_2$ is self-adjoint on suitable spaces), its spectrum is purely imaginary. In case the steady state is a single mode flow, i.e. ${\alpha }_i=0$ for $i>1$, we readily infer as in section 2.2.3 that the spectrum of $-{\tilde{\mathit{L}}}_1$ is also purely imaginary, so that the spectrum of ${\mathbf{L}}_0$ is purely imaginary.

Otherwise, if ${\alpha }_1{\alpha }_2\ne 0$, it appears difficult to determine the spectrum of ${\tilde{\mathit{L}}}_1$ analytically and we resort to numerical computations. For this let $(\cdot {)}_m$ denote the projection onto the mode $\exp (\mathrm{i}{\mathit{k}}_m\mathbf{\cdot }\mathit{x})$. Then, (30) $({\tilde{\mathit{L}}}_1\mathit{v}{)}_m=({\mathbf{L}}_1\mathit{v}{)}_m-\frac{{\mathit{k}}_m\mathbf{\cdot }({\mathbf{L}}_1\mathit{v}{)}_m}{|{\mathit{k}}_m{|}^2}{\mathit{k}}_m=(\mathsf{I}-\mathsf{B}\left({\mathit{k}}_m\right))({\mathbf{L}}_1\mathit{v}{)}_m,$(30) with suitable matrix $\mathsf{B}\left({\mathit{k}}_m\right)$. This admits straightforward numerical computation of spectra on truncated Fourier series. We plot results for an example in figure 9, which gives unstable spectrum and thus strong evidence for unstable spectrum of ${\mathbf{L}}_0$.

Figure 9. Shown is an approximation of part of the spectrum of ${\tilde{\mathit{L}}}_1$. Using (30(30) $({\tilde{\mathit{L}}}_1\mathit{v}{)}_m=({\mathbf{L}}_1\mathit{v}{)}_m-\frac{{\mathit{k}}_m\mathbf{\cdot }({\mathbf{L}}_1\mathit{v}{)}_m}{|{\mathit{k}}_m{|}^2}{\mathit{k}}_m=(\mathsf{I}-\mathsf{B}\left({\mathit{k}}_m\right))({\mathbf{L}}_1\mathit{v}{)}_m,$(30) ) we reduced to two-dimensional wave vectors by fixing the third component of ${\mathit{k}}_m$ at zero. Here, $\mathit{k}=(1,1)$, ${\alpha }_1=0.1,{\alpha }_2=1$. As in figure 7 we use N = 10 wave modes, i.e. $3(2N+1{)}^2$ Fourier modes on the periodic domain $[0,2\pi /k_1]\times [0,2\pi /k_2]$, and Bloch modes in the first component from the grid with distance $\pi /8$. In particular, the spectrum is unstable, so that large amplitude solutions are linear unstable with growth rates proportional to the amplitude.

Notably, this means that for steady states that are mixed mode flows (24(24a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& =\sum _{i=1}^N{\mathrm{e}}^{{\lambda }_{i}t}{\psi }_i{\mathit{k}}_i^{\perp }\text{with}{\psi }_i={\alpha }_i\sin ({\mathit{k}}_i\mathbf{\cdot }\mathit{x}+{\tau }_i),1\le i\le N,\end{array}$(24a) ), it is possible that linear growth rates are proportional to the amplitude parameter a. In contrast, such modes do not exist for steady single mode flows since in this case the spectrum of ${\mathbf{L}}_0$ is purely imaginary, as in section 2.2.3. Again we remark that we expect the growth induced by these modes in the nonlinear system is bounded.

3.2. Flows with vertical structure and coupled buoyancy

The rotating Boussinesq equations with backscatter (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ) also admit explicit solutions of different form in which the velocity and the buoyancy are coupled, and in which the vertical dependence and velocity component are non-trivial. Here we investigate parallel flows, Kolmogorov flows and MGWs. As before, we are particularly interested in the occurrence of unboundedly growing explicit solutions as well as the existence of such steady solutions and their stability properties.

3.2.1. Parallel flow

This class of explicit flows is well known in the inviscid and viscous case (e.g. Wang 1990). It possesses only a vertical velocity component and is thus different from the horizontal flows and admits more general dependence on the horizontal space variables. Specifically, (31) $\mathit{v}(t,\mathit{x})=w(t,x,y){\mathit{e}}_{\mathbf{3}},\mathrm{b}(t,\mathit{x})=\tilde{\mathrm{b}}(t,x,y),p(t,\mathit{x})=\tilde{p}\left(t\right)z,$(31) where w and $\tilde{\mathrm{b}}$ satisfy (with horizontal Laplacian and bi-Laplacian) (32a) $\begin{array}{rl}\frac{\mathrm{\partial }w}{\mathrm{\partial }t}+(d_3{\mathrm{\Delta }}^2+b_3\mathrm{\Delta })w+\tilde{p}& =\tilde{\mathrm{b}},\end{array}$(32a) (32b) $\begin{array}{rl}\frac{\mathrm{\partial }\tilde{\mathrm{b}}}{\mathrm{\partial }t}-\mu \mathrm{\Delta }\tilde{\mathrm{b}}& =-N^{2}w.\end{array}$(32b) Recall that kinetic energy backscatter, which has $d_3=0$ and $b_3\le 0$, has no vertical impact so that parallel flows are in fact independent of backscatter. Plane wave parallel flows can be superposed with the horizontal flows (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) that have zero buoyancy, if their wave vector directions $\mathit{k}$ are the same. In this case, the orthogonality conditions for wave vectors and wave directions are satisfied and the nonlinear terms vanish, so that the superposition of both solutions is also an explicit solution due to the remaining linear system. Thus, a priori, any parallel flow of this form is unboundedly unstable concerning perturbations (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) with the same wave vector $\mathit{k}$ and small enough wave number $|\mathit{k}|$.

Existence and dynamics of parallel flows can be inferred from the dispersion relation of the linear equations (32(32a) $\begin{array}{rl}\frac{\mathrm{\partial }w}{\mathrm{\partial }t}+(d_3{\mathrm{\Delta }}^2+b_3\mathrm{\Delta })w+\tilde{p}& =\tilde{\mathrm{b}},\end{array}$(32a) ). By Fourier transformation with wave vector $\mathit{k}$ and growth rate λ, this is given by $det(\lambda \mathsf{I}-\hat{\mathsf{L}})=0,\hat{\mathsf{L}}:=\left(\begin{array}{cc}-d_3|\mathit{k}{|}^4+b_3|\mathit{k}{|}^2& 1\\ -N^2& -\mu |\mathit{k}{|}^2\end{array}\right)$ or equivalently as the characteristic polynomial (33) ${\lambda }^2+c_1\lambda +c_0=0,$(33) where $c_1:=d_3{K}^2+(\mu -b_3)K$ and $c_0:=\mu (d_{3}K-b_3)K^2+N^2$ with $K:=|\mathit{k}{|}^2\ge 0$. Steady solutions to (32(32a) $\begin{array}{rl}\frac{\mathrm{\partial }w}{\mathrm{\partial }t}+(d_3{\mathrm{\Delta }}^2+b_3\mathrm{\Delta })w+\tilde{p}& =\tilde{\mathrm{b}},\end{array}$(32a) ) require constant $\tilde{p}$ and consist of Fourier modes with $\mathit{k}\ne 0$ that solve (33(33) ${\lambda }^2+c_1\lambda +c_0=0,$(33) ) with $\lambda =0$, i.e. $c_0=0$. Growing spatially non-constant solutions to (32(32a) $\begin{array}{rl}\frac{\mathrm{\partial }w}{\mathrm{\partial }t}+(d_3{\mathrm{\Delta }}^2+b_3\mathrm{\Delta })w+\tilde{p}& =\tilde{\mathrm{b}},\end{array}$(32a) ) exist if and only if (33(33) ${\lambda }^2+c_1\lambda +c_0=0,$(33) ) possesses a root with positive real part and $\mathit{k}\ne 0$, and then do so exponentially and unboundedly. Note that both roots have negative real parts only for $c_1,c_0>0$ and complex conjugate solutions can be superposed to form a real parallel flow solution.

With vertical viscosity, $d_3=0$ and $b_3<0$, and focusing on non-constant solutions $\mathit{k}\ne 0$, we have $c_1>0$ and $c_0=-\mu b_3{K}^2+N^2$. So steady states require $N^2<0$ and then $K^2=N^2/\left(\mu b_3\right)$, or $\mu =N^2=0$ and any K. For $\mu >0$ also, the unstable case $c_0<0$ requires unstable stratification $N^2<0$, and then $c_0<0$ occurs on a disc of wave vectors.

Regarding small amplitudes, analogous to section 3.1.3, any steady (or decaying) parallel flow with small amplitude is unstable, though typically not unboundedly, with respect to unstable modes of the trivial steady state that are exhausted for decreasing amplitude. In the large amplitude scaling, the resulting operator ${\mathbf{L}}_0$ for steady parallel flow $w_0$ is a lower triangular matrix operator with diagonal entries $L_1:=w_0(x,y){\mathrm{\partial }}_z$. Hence, as for the spectrum of ${\mathbf{L}}_0$ in section 3.1.3, the spectrum is given by the diagonal entries. For any (smooth) $w_0$, the operator $L_1$ is skew self-adjoint, similar to ${\mathbf{L}}_2$ in section 3.1.3, and thus the spectrum of ${\mathbf{L}}_0$ is purely imaginary. Hence, no real parts of the spectrum of the steady parallel flow are proportional to its amplitude a.

Finally, in order to illustrate the abstract structure and in preparation of the flows discussed below, next we briefly consider the artificial case $d_3,b_3>0$.

Without thermal diffusion ($\mu =0$), we have $c_0=N^2$ and $c_1={\delta }_3\left(K\right):=(d_{3}K-b_3)K$. For stable stratification, $N^2>0$, growing Fourier modes occur if and only if $c_1<0$ and $c_1^2\ge 4c_0$, which is equivalent to $K<b_3/d_3$ and $N^2\le {\delta }_3^2\left(K\right)/4$, respectively. For $K\in (0,b_3/d_3)$, the global maximum of ${\delta }_3^2\left(K\right)/4$ is $b_3^4/\left(64d_3^2\right)$ at $K=K_1:=b_3/\left(2d_3\right)$; its global minimum is zero at K = 0 and $K=b_3/d_3$. Specifically, if $N^2\le b_3^4/\left(64d_3^2\right)$, i.e. the stability of the stratification is sufficiently weak compared with the backscatter destabilisation, then $c_1^2\ge 4c_0$ for K in a positive interval $I_1\subset (0,b_3/d_3)$; in particular, $I_1=\left\{K_1\right\}$ if $N^2=b_3^4/\left(64d_3^2\right)$. Hence, a parallel flow (31(31) $\mathit{v}(t,\mathit{x})=w(t,x,y){\mathit{e}}_{\mathbf{3}},\mathrm{b}(t,\mathit{x})=\tilde{\mathrm{b}}(t,x,y),p(t,\mathit{x})=\tilde{p}\left(t\right)z,$(31) ) grows exponentially and unboundedly if it contains a Fourier mode with wave vector $\mathit{k}$ in the annulus $\{\mathit{k}\in {\mathbb{R}}^2||\mathit{k}{|}^2\in I_1\}$.

In the presence of thermal diffusion ($\mu >0$), we first note that $c_0$ is a cubic polynomial in K, so $c_0$ has a local maximum at K = 0 and there is a global minimum at some K>0. Specifically, if $0<N^2<4\mu b_3^3/\left(27d_3^2\right)$, then $c_0<0$ in a positive interval $I_2\subset (0,b_3/d_3)$; if $N^2<0$, then $c_0<0$ in an interval $I_3:=[0,K_2)$ for some $K_2>b_3/d_3$. Hence, in this case, a parallel flow (31(31) $\mathit{v}(t,\mathit{x})=w(t,x,y){\mathit{e}}_{\mathbf{3}},\mathrm{b}(t,\mathit{x})=\tilde{\mathrm{b}}(t,x,y),p(t,\mathit{x})=\tilde{p}\left(t\right)z,$(31) ) grows exponentially and unboundedly if it contains a Fourier mode with wave vector $\mathit{k}$ in the annulus $\{\mathit{k}\in {\mathbb{R}}^2||\mathit{k}{|}^2\in I_2\}$ for (not too strongly) stable stratification, or in the disc $\{\mathit{k}\in {\mathbb{R}}^2||\mathit{k}{|}^2\in I_3\}$ for unstable stratification, and its Fourier coefficient vector is not an eigenvector for a possible negative root. Other modes that have $c_0\ge 0$ also yield such growth if $c_1<0$ and $c_1^2\ge 4c_0$. We omit details but note that $c_1<0$ occurs for $K<(b_3-\mu )/d_3$, possibly containing $I_2$.

3.2.2. Kolmogorov flow

Another well-known class of explicit solutions for the Boussinesq equations in the absence of backscatter are the so-called Kolmogorov flows (see, e.g. Balmforth and Young 2005) with wave vectors of the form $\mathit{k}=(k,0,-m{)}^{\mathsf{T}}$, where $k,m\in \mathbb{R}$. Here, we study their occurrence in the case of backscatter and start with the ansatz (34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) and the flow direction $\mathit{a}=\alpha (0,1,0{)}^{\mathsf{T}}+\text{β}(m,0,k{)}^{\mathsf{T}}$.

Compared with the Kolmogorov flows without backscatter and rotation from Prugger and Rademacher (2021), here we have rotation ($f\ne 0$), time dependence ($\lambda \ne 0$) and a nonzero second component of the velocity direction ($\alpha \ne 0$). A superposition of these Kolmogorov flows with the horizontal flow solutions (22(22) $\mathit{v}=\psi (t,\mathit{k}\mathbf{\cdot }\mathit{x}){\mathit{k}}^{\perp },\tilde{p}=f\varphi (t,\mathit{k}\mathbf{\cdot }\mathit{x}),$(22) ) is not possible, since the orthogonality conditions of wave vectors and velocity directions are not satisfied, thus leading to non-gradient terms from the nonlinearity. However, superposition of different Kolmogorov flows is possible, as long as all wave vectors $\mathit{k}$ have the same direction, as is the superposition with such MGWs as discussed in section 3.2.3.

We next determine the necessary relations for the coefficients of (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) in order to solve the Boussinesq equations. For better readability, we define the following terms resulting from the backscatter and thermal diffusion: ${\delta }_{\mu }(k,m)=\mu |\mathit{k}{|}^2,{\delta }_i(k,m)=d_i|\mathit{k}{|}^4-b_i|\mathit{k}{|}^2\text{ for any }1\le i\le 3,$ where $|\mathit{k}{|}^2=k^2+m^2$. Upon inserting (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) into (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ), we find that the coefficients have to satisfy (35) $\left(\begin{array}{cccc}-f& m(\lambda +{\delta }_1)& 0& k\\ \lambda +{\delta }_2& fm& 0& 0\\ 0& k(\lambda +{\delta }_3)& -1& -m\\ 0& N^{2}k& \lambda +{\delta }_{\mu }& 0\end{array}\right)\left(\begin{array}{c}\alpha \\ \text{β}\\ c\\ \gamma \end{array}\right)=0.$(35) For $(k,m)=(0,0)$, these require $\alpha =0$ and c = 0, which is the zero state. From the second row of the $4\times 4$ matrix in (35(35) $\left(\begin{array}{cccc}-f& m(\lambda +{\delta }_1)& 0& k\\ \lambda +{\delta }_2& fm& 0& 0\\ 0& k(\lambda +{\delta }_3)& -1& -m\\ 0& N^{2}k& \lambda +{\delta }_{\mu }& 0\end{array}\right)\left(\begin{array}{c}\alpha \\ \text{β}\\ c\\ \gamma \end{array}\right)=0.$(35) ) and $f\ne 0$, we immediately find that $\alpha =0$ implies $\text{β}m=0$. In case $\text{β}=0$, the Kolmogorov flow (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) is the trivial zero solution, and in case m = 0, (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) is a parallel flow. Hence, we may assume $\alpha \ne 0$. Since (35(35) $\left(\begin{array}{cccc}-f& m(\lambda +{\delta }_1)& 0& k\\ \lambda +{\delta }_2& fm& 0& 0\\ 0& k(\lambda +{\delta }_3)& -1& -m\\ 0& N^{2}k& \lambda +{\delta }_{\mu }& 0\end{array}\right)\left(\begin{array}{c}\alpha \\ \text{β}\\ c\\ \gamma \end{array}\right)=0.$(35) ) is a homogeneous linear system in $(\alpha ,\text{β},c,\gamma )$, non-trivial solutions require a kernel of the associated matrix. Hence, either there is no non-trivial Kolmogorov flow or a linear space of these, which requires vanishing determinant of this matrix. Assuming $(k,m)\ne (0,0)$ and dividing by $-|\mathit{k}{|}^2$, this gives (36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) with coefficients $\begin{array}{rl}{c}_2& :=({\delta }_2+{\delta }_{\mu })+|\mathit{k}{|}^{-2}({\delta }_3{k}^2+{\delta }_1{m}^2),\\ c_1& :={\delta }_2{\delta }_{\mu }+|\mathit{k}{|}^{-2}\left[\right({\delta }_2+{\delta }_{\mu }\left)\right({\delta }_3{k}^2+{\delta }_1{m}^2)+N^2{k}^2+f^2{m}^2],\\ c_0& :=|\mathit{k}{|}^{-2}\left[{\delta }_2{\delta }_{\mu }\right({\delta }_3{k}^2+{\delta }_1{m}^2)+{\delta }_2{N}^2{k}^2+{\delta }_{\mu }{f}^2{m}^2].\end{array}$

Steady Kolmogorov flows and linear stability

The condition $\lambda =0$ for steady Kolmogorov flow reduces (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) to (37) ${\delta }_2{\delta }_{\mu }({\delta }_3{k}^2+{\delta }_1{m}^2)+{\delta }_2{N}^2{k}^2+{\delta }_{\mu }{f}^2{m}^2=0.$(37) For comparison, note that the left-hand side identically equals to zero in the absence of backscatter and viscosity ($d_j,b_j=0$ for $j=1,2,3$) and thermal diffusion ($\mu =0$) so that in this case, steady and non-trivial Kolmogorov flow exists for all $(k,m)\in {\mathbb{R}}^2\setminus \left\{\right(0,0\left)\right\}$. In contrast, in the presence of backscatter $d_2,b_2>0$, but still without thermal diffusion ($\mu =0$), only ${\delta }_2{N}^2{k}^2=0$ remains so that either k = 0 or ${\delta }_2=0$. The latter means $|\mathit{k}{|}^2=b_2/d_2$, so that non-trivial and steady Kolmogorov flow occurs on the m-axis and the circle in the $(k,m)$-plane with radius $\sqrt{b_2/d_2}$, cf. figure 10(b,d). Conversely, we can create steady Kolmogorov flows for any wave vector $(k,m)\ne (0,0)$ by suitable choice of $d_2,b_2$ such that ${\delta }_2=0$.

Figure 10. Unboundedly growing Kolmogorov flows occur for wave vectors $(k,m)$ in the red regions for parameter sets $\mathrm{\Lambda }:=\{N^2,d_1,b_1,d_2,b_2,d_3,b_3,\mu \}$ (with isotropic backscatter), where the Coriolis parameter is fixed at f = 1. Red regions: one of conditions (i), (ii) and (iii) is satisfied; blue regions: none of conditions (i), (ii) and (iii) is satisfied; black curves in (a): loci of $c_0=c_1=c_2=0$; black curves in (b–e): loci of $c_0=0$; white dots: the zero state at $(k,m)=(0,0)$ (Colour online).

Concerning stability, analogous to section 3.1.3, small amplitude steady Kolmogorov flows are unstable, though typically not unboundedly, due to the instability of the zero state under backscatter. In the large amplitude scaling, the resulting operator ${\mathbf{L}}_0$ is a triangular block matrix operator, similar to section 3.2.1, with skew-adjoint parts that imply purely imaginary spectrum. Hence, there are again no growth rates that are proportional to the amplitude of the steady Kolmogorov flow. Regarding stability of Kolmogorov flows without backscatter, we refer to Balmforth and Young (2005).

A source of unbounded instability of steady Kolmogorov flows are possible superpositions with MGWs discussed in section 3.2.3 below. In the following, we will examine the existence of exponentially and unboundedly growing Kolmogorov flows, which then also proves unbounded instability of steady Kolmogorov flows due to possible superpositions.

Unboundedly growing Kolmogorov flows

Such flows with $\text{β},k\ne 0$ transfer the horizontal backscatter to growing vertical velocity component. They correspond to positive roots of (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ), which occur as follows in terms of the sign of $c_0$:

1. If $c_0<0$, then (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) has a positive root.

2. For $c_0=0$, (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) has a positive root if and only if $-c_2+\sqrt{c_2^2-4c_1}>0$.

3. For $c_0>0$, (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) has a positive root if and only if $-c_2+\sqrt{c_2^2-3c_1}>0$ and $2c_2^3-9c_1{c}_2+27c_0+(6c_1-2c_2^2)\sqrt{c_2^2-3c_1}\le 0$.

The conditions in (iii) imply that the local minimum of the cubic polynomial in (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) lies at a positive value and the value on the local minimum is non-positive.

For comparison, we start with the common situation without backscatter, viscosity and thermal diffusion ($\mu ,b_j,d_j=0$ for $j=1,2,3$), where a growing Kolmogorov flow (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) requires unstable stratification $N^2<0$. Indeed, in this case (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) has $c_0=c_2=0$ and $c_1=|\mathit{k}{|}^{-2}(N^2{k}^2+f^2{m}^2)$, so that a positive root occurs if and only if $c_1<0$, which requires $N^2<0$, and thus growing solutions occur for $m^2/k^2<-N^2/f^2$, cf. figure 10(a). More precisely, for such wave vectors, a steady flow co-exists with a growing and a decaying flow (on the red regions in figure 10(a)), which turn into a triple steady flow on the boundary, where $c_0=c_1=c_2=0$ (black curves).

In the presence of (horizontal) backscatters ($b_j,d_j>0$ for $j=1,2$), growth of Kolmogorov flows (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) is possible also for stable stratification $N^2>0$. We next focus on case (i) with negative $c_0$ and omit details of cases (ii) and (iii) with non-negative $c_0$. Some examples are plotted in figure 10.

First, we note that, for $N^2>0$, the coefficient $c_0$ is negative for sufficiently small k and m. Its sign is that of the left-hand side of (37(37) ${\delta }_2{\delta }_{\mu }({\delta }_3{k}^2+{\delta }_1{m}^2)+{\delta }_2{N}^2{k}^2+{\delta }_{\mu }{f}^2{m}^2=0.$(37) ), whose leading order term as $(k,m)\to (0,0)$ is $|\mathit{k}{|}^2(-b_2{N}^2{k}^2+\mu f^2{m}^2)$ and can be negative only if $N^2>0$ and $k\ne 0$. Then, to leading order, growing Kolmogorov flows occur for $m^2/k^2<b_2{N}^2/\left(\mu f^2\right)$ with $\mu >0$ (cf. figure 10(e)) and increase of $N^2$ or $b_2$ enlarges this region in the $(k,m)$-plane near the origin, while increase of μ or f shrinks it. Rewriting the condition as $N^2>\mu f^2{m}^2/\left(b_2{k}^2\right)$, for a given wave vector, then on the one hand increase of μ or f requires sufficiently stable stratification, and on the other hand increase of $b_2$ allows for less stable stratification. The leading order term of the left-hand side of (37(37) ${\delta }_2{\delta }_{\mu }({\delta }_3{k}^2+{\delta }_1{m}^2)+{\delta }_2{N}^2{k}^2+{\delta }_{\mu }{f}^2{m}^2=0.$(37) ) is always positive if $N^2<0$, so $c_0>0$ near the origin, which refers to case (iii) above; we omit details and just plot an example in figure 10(c).

Second, for general $(k,m)$ and vanishing thermal diffusion ($\mu =0$), only ${\delta }_2{N}^2{k}^2$ remains. Hence, $c_0<0$ for $|\mathit{k}{|}^2<b_2/d_2$ if $N^2>0$ and $k\ne 0$ (cf. figure 10(d)), and $c_0<0$ for $|\mathit{k}{|}^2>b_2/d_2$ if $N^2<0$ and $k\ne 0$ (cf. figure 10(b)). For k = 0, $c_0$ is zero so we refer to the case (ii) above for the growing Kolmogorov flows. The numerical computations show that in the examples of figure 10(b,d), the condition in case (ii) is not satisfied for k = 0 so that here a Kolmogorov flow (34(34) $\mathit{v}(t,\mathit{x})={\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x})\mathit{a},\mathrm{b}(t,\mathit{x})=c{\mathrm{e}}^{\lambda t}\cos (\mathit{k}\mathbf{\cdot }\mathit{x}),p(t,\mathit{x})=\gamma {\mathrm{e}}^{\lambda t}\sin (\mathit{k}\mathbf{\cdot }\mathit{x}),$(34) ) is not growing. Near $\mathit{k}=0$, the examples of figure 10(b,c) can be regarded as perturbations of the example of figure 10(a). As before, the red region indicates existence of growing flows; note that decaying or oscillating ones may co-exist. Notably, in figure 10(b,c), the part of the red region's boundary without black marking stems from a bifurcation of saddle-node-type in terms of the wave vector, where (36(36) ${\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0,$(36) ) possesses a positive double root. In contrast, the black curves mark loci of steady flows, where $c_0=0$.

Similar to parallel flow, we briefly consider artificial vertical backscatter since isotropic backscatter $d_j=d>0,b_j=b>0$, $j=1,2,3$, admits an analytical consideration. In this case $\delta :=d|\mathit{k}{|}^4-b|\mathit{k}{|}^2={\delta }_j$ for $j=1,2,3$. With thermal diffusion ($\mu >0$), we then have $c_0<0$ for small wave number $|\mathit{k}{|}^2<b/d$ (i.e. $\delta <0$) if stratification is sufficiently stable $N^2>-\mu |\mathit{k}{|}^2(f^2{m}^2+{\delta }^2|\mathit{k}{|}^2)/\left(k^2\delta \right)>0$, and for large wave number $|\mathit{k}{|}^2>b/d$ (i.e. $\delta >0$) only for sufficiently unstable stratification $N^2<-\mu |\mathit{k}{|}^2(f^2{m}^2+{\delta }^2|\mathit{k}{|}^2)/\left(k^2\delta \right)<0$. Unlike the growing parallel flow, which requires sufficiently weak stable stratification, here the increase of positive $N^2$ increases the set of wave vectors – within the small wave number region – for growing Kolmogorov flow. Another observation is that for the fixed stratification – stable or unstable – larger μ leads to smaller region in wave vector space of growing solutions.

The latter can also be observed (numerically) without vertical backscatter $d_3=0,b_3<0$ as shown in figure 10(d,e) for $N^2>0$ and figure 10(b,c) for $N^2<0$.

3.2.3. MGWs

The last kind of explicit plane wave-type solutions we are aware of are the so-called MGWs, which are for example discussed in Achatz (2006) for the inviscid Boussinesq equations. Here, we study the occurrence in the rotating Boussinesq equations with backscatter. These solutions again form an invariant subspace of linear dynamics since the nonlinear terms vanish, but structurally differ from the aforementioned flows. In particular, the wave profile of a MGW is a time-dependent travelling wave with phase variable $\xi =kx+mz-\omega t$ and takes the form (38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) (38b) $\begin{array}{rl}\mathrm{b}(t,\mathit{x})& ={\text{β}}_1{\mathrm{e}}^{\lambda t}\sin \xi +{\text{β}}_2\omega {\mathrm{e}}^{\lambda t}\cos \xi ,\end{array}$(38b) (38c) $\begin{array}{rl}p(t,\mathit{x})& ={\gamma }_1{\mathrm{e}}^{\lambda t}\cos \xi +{\gamma }_2\omega {\mathrm{e}}^{\lambda t}\sin \xi .\end{array}$(38c) The conditions for these to be (non-trivial) explicit solutions of (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ), with $\mathit{k}=(k,m{)}^{\mathsf{T}}\in {\mathbb{R}}^2\mathrm{\setminus }\{(0,0{)}^{\mathsf{T}}\}$, in particular depend on ω. In fact, for $\omega =0$ these are Kolmogorov flows from section 3.2.2 with $\text{β}=0$, and the existence conditions read, after inserting (38(38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) ) into (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ), (39) $\left(\begin{array}{ccc}\lambda +{\delta }_2& 0& 0\\ mf& -k& 0\\ kf& m& |\mathit{k}{|}^2\\ 0& \lambda +{\delta }_{\mu }& 0\end{array}\right)\left(\begin{array}{c}{\alpha }_1\\ {\text{β}}_1\\ {\gamma }_1\end{array}\right)=0,$(39) with ${\delta }_j={\delta }_j\left(\mathit{k}\right):=d_j|\mathit{k}{|}^4-b_j|\mathit{k}{|}^2$ for any $1\le j\le 3$ and ${\delta }_{\mu }={\delta }_{\mu }\left(\mathit{k}\right):=\mu |\mathit{k}{|}^2$ as above.

For $\omega \ne 0$, the existence conditions are at first the following eight: (40) $\left(\begin{array}{cccccc}\lambda +{\delta }_2& 0& 0& 0& 0& 0\\ mf& -k& 0& {\omega }^2|\mathit{k}{|}^2& 0& 0\\ kf& m& |\mathit{k}{|}^2& 0& 0& 0\\ 0& \lambda +{\delta }_{\mu }& 0& 0& {\omega }^2& 0\\ 1& 0& 0& mf& 0& 0\\ 0& -1& 0& k{N}^2& \lambda +{\delta }_{\mu }& 0\\ 0& 0& 0& {\delta }_3{k}^2+{\delta }_1{m}^2+\lambda |\mathit{k}{|}^2& -k& 0\\ 0& 0& 0& ({\delta }_3-{\delta }_1)km& -m& |\mathit{k}{|}^2\end{array}\right)\left(\begin{array}{c}{\alpha }_1\\ {\text{β}}_1\\ {\gamma }_1\\ {\alpha }_2\\ {\text{β}}_2\\ {\gamma }_2\end{array}\right)=0.$(40) However, most of these can be readily solved directly in terms of the coefficients. In the following, we explicitly determine all non-trivial solutions (38(38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) ). Afterwards we shortly discuss the stability of MGWs and the possible superpositions with other types of explicit solutions.

MGWs with steady phase ($\mathit{\omega }=\mathbf{\text{0}}$)

We start with $\omega =0$, which gives a certain type of Kolmogorov flows. There is no propagation of the travelling wave profile and the second terms of $\mathit{v}$, b and p in (38(38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) ) vanish. Solutions with $k=0,{\alpha }_1=0,{\text{β}}_1=-{\gamma }_{1}m,\lambda =-{\delta }_{\mu }$ have $\mathit{v}\equiv 0$, the pressure depends on the buoyancy only, and there are no further conditions on the wave vector $\mathit{k}=(0,m{)}^{\mathsf{T}}$.

Other solutions with $\omega =0$ satisfy $k\ne 0,{\text{β}}_1={\alpha }_1\frac{m}{k}f,{\gamma }_1=-{\alpha }_1\frac{f}{k},\lambda =-{\delta }_{\mu },{\delta }_2={\delta }_{\mu },$ so $\mathit{v}\ne 0$ for ${\alpha }_1\ne 0$. The last of these equations gives a condition on the wave vector $\mathit{k}$, which is equivalent to $d_2|\mathit{k}{|}^2-b_2-\mu =0$ and constrains the wave number to $|\mathit{k}{|}^2=(b_2+\mu )/d_2$.

In both of these cases, the growth rate is defined by the thermal diffusion as $\lambda =-\mu |\mathit{k}{|}^2$, so that all these MGWs are exponentially decaying for $\mu >0$ and steady for $\mu =0$. Furthermore, both solutions only depend on the parameter of the second momentum or the buoyancy equation but are independent of the Brunt–Väisälä frequency $N^2$.

MGWs with oscillating phase ($\omega \ne 0$)

We turn to non-trivial MGW solutions (38(38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) ) with $\omega \ne 0$. The simplest class are “vertically varying” MGWs with $\begin{array}{c}k=0,{\alpha }_1=-{\alpha }_{2}mf,{\text{β}}_1={\text{β}}_2={\gamma }_1={\gamma }_2=0,\\ \lambda =-{\delta }_2,\omega =\pm f,{\delta }_1={\delta }_2,\end{array}$ so that $\mathrm{b}\equiv p\equiv 0$, while $\mathit{v}\ne 0$ – in contrast to the case $k=\omega =0$ above. Notably, these solutions depend on the parameters of the horizontal momentum equations only, while the solutions with $k=\omega =0$ only depend on those from the buoyancy equation. The last of the above equations is a condition on the wave vector $\mathit{k}=(0,m{)}^{\mathsf{T}}$, which is equivalent to $(d_1-d_2)m^2+b_2-b_1=0$. Thus, these solutions exist for all $m\ne 0$ in the isotropic case, while for $d_1\ne d_2$ they are restricted to $m^2=(b_1-b_2)/(d_1-d_2)>0$.

Since $\lambda =-{\delta }_2$, these MGWs are exponentially and unboundedly growing for wave numbers $m^2<b_2/d_2$, and steady for $m^2=b_2/d_2$, thus transferring the horizontal backscatter to growing vertical dependence. In particular, the thermal diffusion μ and Brunt–Väisälä frequency $N^2$ have no impact, since these solutions trivially satisfy the buoyancy equation by $v_3\equiv \mathrm{b}\equiv 0$.

Another type are “zonally varying” MGWs with the coefficients satisfying $\begin{array}{c}{\alpha }_1=m=0,{\text{β}}_1={\alpha }_2(N^2-\frac{({\delta }_3-{\delta }_{\mu }{)}^2}{4})k,{\text{β}}_2={\alpha }_2\frac{{\delta }_3-{\delta }_{\mu }}{2}k,\\ {\gamma }_1={\gamma }_2=0,\lambda =-\frac{{\delta }_3+{\delta }_{\mu }}{2},{\omega }^2=N^2-\frac{({\delta }_3-{\delta }_{\mu }{)}^2}{4}.\end{array}$ These have vanishing pressure and depend only on the parameters of the buoyancy and vertical momentum equations. In contrast to the MGW solutions before, these here do also depend on the Brunt–Väisälä frequency $N^2$. Furthermore, different from the solutions before, the growth rate λ and phase frequency ω depend on both the vertical term ${\delta }_3$ and the thermal diffusion. Since ${\omega }^2>0$, the last equation is a condition on the wave vector $\mathit{k}=(k,0{)}^{\mathsf{T}}$, which requires stable stratification $N^2>0$ satisfying $({\delta }_3-{\delta }_{\mu }{)}^2<4N^2$. In case $N^2>0$, these MGW solutions exist at least for sufficiently small $k^2$, since ${\delta }_3,{\delta }_{\mu }\to 0$ as $|\mathit{k}|\to 0$. Due to the equation for the growth rate λ, these MGWs are exponentially decaying for kinetic energy backscatter, where ${\delta }_3=-b_3|\mathit{k}{|}^2>0$. Hence, this kind of MGWs are steady or exhibit growth only in the artificial case $d_3,b_3>0$ for wave numbers $k^2\le (b_3-\mu )/d_3$, i.e. only in the case $b_3>\mu$.

The existence analysis for the remaining solutions to (38(38a) $\begin{array}{rl}\mathit{v}(t,\mathit{x})& ={\alpha }_1{\mathrm{e}}^{\lambda t}\sin \left(\xi \right)(0,1,0{)}^{\mathsf{T}}+{\alpha }_2\omega {\mathrm{e}}^{\lambda t}\cos (\xi \left)\right(-m,0,k{)}^{\mathsf{T}},\end{array}$(38a) ) with $\omega ,k,{\alpha }_1\ne 0$ is more involved. These solutions have the coefficients $\begin{array}{c}{\alpha }_1=-{\alpha }_{2}mf,{\text{β}}_1={\alpha }_1\frac{m}{k}f+{\alpha }_2\frac{{\omega }^2}{k}|\mathit{k}{|}^2,{\text{β}}_2=\frac{{\alpha }_2}{k}({\delta }_3{k}^2+{\delta }_1{m}^2-{\delta }_2|\mathit{k}{|}^2),\\ {\gamma }_2={\alpha }_2({\delta }_1-{\delta }_3)\frac{km}{|\mathit{k}{|}^2}+{\text{β}}_2\frac{m}{|\mathit{k}{|}^2},{\gamma }_1=-{\alpha }_1\frac{kf}{|\mathit{k}{|}^2}-{\text{β}}_1\frac{m}{|\mathit{k}{|}^2},\lambda =-{\delta }_2,\end{array}$ with additional conditions, that also define the phase frequency ω, given by

(41a) $\begin{array}{rl}{\omega }^2=\left(\right({\delta }_3-{\delta }_2\left)\right({\delta }_{\mu }-{\delta }_2)+N^2)\frac{k^2}{|\mathit{k}{|}^2}+\left(\right({\delta }_1-{\delta }_2\left)\right({\delta }_{\mu }-{\delta }_2)+f^2)\frac{m^2}{|\mathit{k}{|}^2}& >0,\end{array}$(41a) (41b) $\begin{array}{rl}\left(\right({\delta }_{\mu }-{\delta }_2)N^2+({\delta }_3-{\delta }_2\left)(\right({\delta }_{\mu }-{\delta }_2{)}^2+{\omega }^2))k^2+({\delta }_1-{\delta }_2)(({\delta }_{\mu }-{\delta }_2{)}^2+{\omega }^2)m^2& =0.\end{array}$(41b)

The first condition (41a) has solutions $(k,m)$ for any parameter: for stable stratification $N^2>0$, the factors of $k^2$ and $m^2$ are both positive for $|\mathit{k}|$ sufficiently small, since ${\delta }_j,{\delta }_{\mu }\to 0$ for $|\mathit{k}|\to 0$ for all $1\le j\le 3$. For unstable stratification $N^2<0$, one can first choose $|\mathit{k}|$ sufficiently small, so that the factor of $m^2$ is positive. Then, for a fixed $|\mathit{k}{|}^2=s^2$ so that ${\delta }_j$ and ${\delta }_{\mu }$ are constant for these wave vectors, one can choose $k^2=s^2-m^2$ small enough, so that the term with $k^2$ does not make the whole expression negative.

We omit a complete analysis of the more complicated conditions (41) here, but note that the horizontally isotropic case ${\delta }_1={\delta }_2$ requires unstable stratification for the existence of unboundedly growing MGWs of this kind. Both conditions together generate a set of solutions with rather complex structure, as plotted in figure 11 for the cases with stable stratification and (weakly) anisotropic backscatter as well as unstable stratification with isotropic backscatter, which generate unboundedly growing explicit solutions.

Figure 11. Existence of MGW solutions on the wave vector plane $(k,m)$ in the case $\omega ,{\alpha }_1\ne 0$ (white curves) with growth rate λ positive (red), zero (black) or negative (blue). Fixed parameter: $d_1=0.12,d_3=0.0,b_1=0.14,b_3=-0.1,\mu =0.12,f=0.3$. (a) Stable stratification $N^2=1$ with anisotropic backscatter $d_2=0.11,b_2=0.12$. (b) Unstable stratification $N^2=-0.05$ with isotropic backscatter $d_2=0.12,b_2=0.14$ (Colour online).

Due to the equation for the growth rate $\lambda =-{\delta }_2$, these MGWs are exponentially and unboundedly growing for wave numbers $|\mathit{k}{|}^2<b_2/d_2$ and steady states for $|\mathit{k}{|}^2=b_2/d_2$ (see figure 11). Notably, since $k\ne 0$, these flows can transfer the horizontal backscatter (in fact only the meridional component) to growing vertical velocity.

Superpositions and stability

It is possible to superpose MGWs and Kolmogorov flows as a solution to (20(20a) $\begin{array}{rl}\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}+(\mathit{v}\mathbf{\cdot }\mathbf{\nabla })\mathit{v}+f{\mathit{e}}_{\mathbf{3}}\times \mathit{v}+\mathbf{\nabla }p-{\mathit{e}}_{\mathbf{3}}\mathrm{b}& =-(\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(d_1,d_2,d_3){\mathrm{\Delta }}^2+\mathsf{d}\mathsf{i}\mathsf{a}\mathsf{g}(b_1,b_2,b_3)\mathrm{\Delta })\mathit{v},\end{array}$(20a) ) if these have the same direction of wave vectors $\mathit{k}$, and this can also be in the form of an integral. Indeed, the similar structure of wave vector and velocity direction yields vanishing nonlinear terms and a remaining system of linear equations, which each solution satisfies. Note that Kolmogorov flows exist on the whole wave vector space $(k,m{)}^{\mathsf{T}}\in {\mathbb{R}}^2$, while MGW in general not. This means that superposition in one wave vector direction is possible for arbitrary wave number of the Kolmogorov flow, but the wave number of the MGW is in general restricted. Depending on the wave vectors, the Kolmogorov flows and MGWs in such a superposition can be steady, exponentially growing or decaying.

With this superposition, we can prove in certain cases the unbounded instability of steady MGW solutions due to perturbations with exponentially growing Kolmogorov flows, and vice versa. In section 3.2.2, we found that growing Kolmogorov flows in case $\mu >0$ may occur only in certain directions. In this case, steady MGW solutions are unboundedly unstable a priori only in these certain wave vector directions (see investigation of $c_0$ in section 3.2.2). Without thermal diffusion ($\mu =0$), there are always growing Kolmogorov flows in any direction with $k\ne 0$, for stable stratification $N^2>0$ at least for all $|\mathit{k}{|}^2<b_2/d_2$, and for unstable stratification $N^2<0$ for $|\mathit{k}{|}^2>b_2/d_2$. Thus, in these cases, the steady MGWs with $k\ne 0$ are always unboundedly unstable with respect to certain Kolmogorov flows. However, steady MGWs with k = 0, as well as Kolmogorov flows, are also unboundedly unstable, since there are exponentially and unboundedly growing MGWs with k = 0 and $m^2<b_2/d_2$, with which they can be superposed (see the case of “vertically varying” MGWs with $\omega \ne 0$ and k = 0).

As in section 3.1.3, steady small amplitude MGWs are unstable due to the unstable zero state under backscatter. In the large amplitude scaling, the resulting operator ${\mathbf{L}}_0$ is a triangular block matrix operator with skew-adjoint parts that imply purely imaginary spectrum as for Kolmogorov flows. Hence, there are again no unstable eigenvalues that scale with the amplitude of the steady MGW.

4. Discussion

Motivated by the numerical backscatter scheme (Jansen and Held 2014, Danilov et al. 2019, Juricke et al. 2020, Perezhogin 2020), we have studied the impact of simplified kinetic energy backscatter via classes of explicit flows in the shallow water and rotating Boussinesq equations on the whole space. Here, we have found that backscatter induces unbounded instability of the zero state, as well as of certain non-trivial steady solutions, in the sense of unboundedly and exponentially growing flows, also with vertical structure and also for stable stratification. This highlights the possibility of concentration of energy due to backscatter and is in contrast to the desired energy redistribution. Since all flows we have considered simultaneously solve the linear equations from dropping the transport nonlinearity and the full nonlinear equations, these are directly linked to spectrum and linear stability of the trivial flow. Indeed, in wave vector space, the condition to solve the nonlinear equations emposes additional constraints, which are stronger in the shallow water case than in the purely horizontal Boussinesq equations. We have discussed the corresponding linear and nonlinear constraints on coefficients for these and also for flows with vertical structure and coupled buoyancy that relate to parallel flows, Kolmogorov flows and MGWs.

We have identified superposition principles of these flows in the nonlinear equations and have discussed the resulting unbounded instability of the flows themselves, in particular of steady flows induced by the backscatter. Due to the linear–nonlinear structure, these steady flows come as a family with an amplitude scaling parameter and we have discussed linear stability in the regimes of small and large amplitudes. Flows of small amplitudes inherit the instability of the zero state, but the treatment of large amplitudes is more subtle. Here we have considered a renormalised eigenvalue problem and have found that the resulting spectrum is purely imaginary in all except one class of the aforementioned flows. Based on numerical computations, in the exceptional case of steady multi-mode horizontal flows, unstable rates are proportional to the amplitude of the steady flow, thus leading to arbitrarily strong growth rates.

The approach that we have presented can be applied to operators with other constant coefficient linear or derivative terms. An interesting case that we will pursue further is the inclusion of bottom drag, which admits an onset of instability with finite wave number. Also forcing of plane wave form can be treated, similar to Prugger and Rademacher (2021).

It would be interesting to study analytically and numerically in what way the undesired growth occurs in numerical discretisations. For fixed backscatter coefficients under decreasing grid size, the solutions of the discrete system are expected to track those of the continuum over increasingly long time intervals. Hence, the presence of unbounded instability in the limit is expected to readily imply large, but possibly finite growth for any fixed small grid size. The actual backscatter scheme scales the strength of backscatter with the grid size so that it would be interesting to identify, for a given discretisation scheme, the relation between the growth of flows and the grid-scaling of backscatter coefficients. This may provide another approach to energetic consistency of backscatter schemes.

Acknowledgments

The authors thank Marcel Oliver, Gualtiero Badin, Stephan Juricke and Ulrich Achatz for fruitful discussions, and the anonymous reviewers for comments that helped improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This paper is a contribution to the project M2 (Systematic multi-scale modelling and analysis for geophysical flow) of the Collaborative Research Centre TRR 181 “Energy Transfers in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project number 274762653.