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

A. Prugger${{\dagger}}$ and J.D.M. Rademacher${{\dagger}}$ and J.
Yang${{\ddagger}}$
${{\dagger}}$University of Bremen, Department 3 - Mathematics, 28359 Bremen,
Germany
${{\ddagger}}$Sun Yat-sen University, School of Mathematics (Zhuhai), 519082
Zhuhai, China Corresponding author. Email<EMAIL_ADDRESS>
###### 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 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. (Chai et al., 2020; Drazin
and Riley, 2006; Dyck and Straatman, 2019; Majda, 2003; Majda and Wang, 2006;
Weinbaum and O’Brien, 1967), 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 are simultaneously
solving the nonlinear equations and the linear equation that arise from
dropping the transport nonlinearity, i.e. the linearisation in the zero state.
These include flows that correspond to barotropic, parallel and Kolmogorov
flows as well as monochromatic inertia gravity waves, cf. (Achatz, 2006;
Balmforth and Young, 2005; Prugger and Rademacher, 2021; Yau et al., 2004) for
vanishing backscatter. More specifically, we identify these solutions from
superpositions of plane waves with suitably selected wave vectors and wave
directions. These are therefore specific eigenmodes of the linear operator
$\mathcal{L}$ arising from infinitesimal perturbations of the zero state and
thus correspond to elements in the spectrum of $\mathcal{L}$. Since these are
equally solutions of the nonlinear equations we find subspaces of linear
dynamics in the nonlinear system, where all these solutions come as a family
with a free scaling factor. We confirm that the presence of backscatter
implies unstable spectrum of $\mathcal{L}$, i.e. spectrum with positive growth
rates, also for stable stratification. Moreover, we identify explicit
solutions of the nonlinear problem that grow exponentially and unboundedly,
including solutions whose vertical velocity is growing unboundedly, despite
the horizontal nature of the backscatter. This highlights the possibility of
undesired concentration of energy due to backscatter, and is in contrast to
the targeted energy redistribution. The solutions which have negative real
spectrum likewise illustrate the possibility of ineffective energy input into
certain scales. However, while solutions to consistent discretisations shadow
continuum effects over finite times, the specific implications for numerical
backscatter schemes require additional investigations beyond the present
study.

We refer to this unbounded exponential growth of linear type in the nonlinear
system as unbounded instability and note that it does not occur in generic
damped-driven evolution equations, where unstable manifolds are nonlinear.
Among these explicit flows there are numerous steady states and we identify
admissible superpositions in the nonlinear equations, which provide explicit
conditions for unbounded instability in many cases. For certain steady
barotropic flows of this type we provide numerical evidence of eigenmodes
whose growth rates are proportional to the amplitude factor of these steady
flows, thus featuring arbitrarily strong local instability. For all other
steady solutions of the type considered, we analyse the resulting large
amplitude regime and prove that growth rates must scale sublinear with the
amplitude. We next summarise and describe the results in more detail.

We first consider the rotating shallow water equations with backscatter and
flat bottom topography in the $f$-plane approximation given by

(1.1c) $\displaystyle\frac{\partial\boldsymbol{v}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}$
$\displaystyle=-f\boldsymbol{v}^{\perp}-g\nabla\eta-\left(\begin{array}[]{cc}d_{1}\Delta^{2}+b_{1}\Delta&0\\\
0&d_{2}\Delta^{2}+b_{2}\Delta\end{array}\right)\boldsymbol{v}$ (1.1d)
$\displaystyle\frac{\partial\eta}{\partial t}+(\boldsymbol{v}\cdot\nabla)\eta$
$\displaystyle=-(H_{0}+\eta)\text{div}(\boldsymbol{v})\,,$

where $\boldsymbol{v}=\boldsymbol{v}(t,\boldsymbol{x})\in\mathbb{R}^{2}$ is
the velocity field on the whole space $\boldsymbol{x}\in\mathbb{R}^{2}$ at
time $t\geq 0$ and $\eta=\eta(t,\boldsymbol{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\neq 0$ is the
constant Coriolis parameter, $g>0$ gravity acceleration and the backscatter
parameters are $b_{1},\,b_{2},\,d_{1},\,d_{2}>0$.
The explicit solutions that we consider in this case derive from the plane
wave ansatz

(1.2)
$\displaystyle\boldsymbol{v}=\psi(t,\boldsymbol{k}\cdot\boldsymbol{x})\boldsymbol{k}^{\perp}\,,\quad\eta=\frac{f}{g}\phi(\boldsymbol{k}\cdot\boldsymbol{x})\,,$

which implies vanishing nonlinear terms and, for single Fourier modes in
$\psi$ and $\phi$, readily yields certain nonlinear relations of parameters in
order to obtain solutions to (1.1). These turn into geostrophically balanced
Rossby waves when the backscatter parameters tend to zero. We consider the
wave vectors $\boldsymbol{k}$ as the primary parameters, which (in this case)
admits a direct comparison to the spectrum of $\mathcal{L}$, where each
$\boldsymbol{k}\in\mathbb{R}^{2}$ corresponds to up to three eigenmodes. From
this viewpoint, the set of $\boldsymbol{k}$ that admit plane waves as
solutions to the nonlinear equations are more constrained, forming curves
$\boldsymbol{k}(s)$, with $s\in\mathbb{R}$, in the $\boldsymbol{k}$-plane with
corresponding growth rates $\lambda(s)\in\mathbb{R}$. We study the geometry of
$\boldsymbol{k}(s)$ and find that, up to sign, at most three lie on the same
line in the $\boldsymbol{k}$-plane through the origin. Due to a radial
superposition principle that allows to superpose plane wave flows whose wave
vectors lie in the same line through the origin, this gives rise to several
invariant subspace with linear dynamics of dimensions $n=1,\,2,\,3$. Regarding
the backscatter coefficients, the latter requires anisotropy in the sense that
$n>1$ requires $b_{1}/d_{1}\neq b_{2}/d_{2}$. In addition to $\boldsymbol{k}$,
an organising parameter in this analysis is the relative size of the
amplitudes of $\boldsymbol{v}$ and $\eta$, which lies in a certain interval
for existing steady solutions (1.2) and determines those, which are
unboundedly unstable.
Moreover, we show that when $\lambda(s)$ changes sign, the transition to
unstable spectrum related to plane waves (1.2) is a long-wavelength,
modulational instability; this occurs in addition to unstable spectrum from
other modes that we numerically find. The free amplitude parameter $a$ of
steady plane flows naturally leads to the asymptotic regimes $|a|\ll 1$ and
$|a|\gg 1$ of small and large amplitudes. In the small amplitude regime the
unstable spectrum of the zero state generates unstable spectrum of the plane
wave flows, implying positive growth rates in addition to the explicit ones,
but these are not expected to create unbounded growth. In order to study the
large amplitude regime, we consider a rescaled problem and show that the
resulting operator possesses purely imaginary spectrum. Hence, for these cases
the unstable growth rates must scale sublinear with respect to the amplitude
factor $a$. Numerical computations suggest such unstable spectrum indeed
occurs.

We then turn to the rotating Boussinesq equations augmented with backscatter

(1.3a) $\displaystyle\frac{\partial\boldsymbol{v}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}+f\boldsymbol{e_{3}}\times\boldsymbol{v}+\nabla
p-\boldsymbol{e_{3}}\mathtt{b}$
$\displaystyle=-\mathrm{diag}\bigl{(}d_{1}\Delta+b_{1},d_{2}\Delta+b_{2},-\nu\bigr{)}\Delta\boldsymbol{v}$
(1.3b) $\displaystyle\nabla\cdot\boldsymbol{v}$ $\displaystyle=0$ (1.3c)
$\displaystyle\frac{\partial\mathtt{b}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\mathtt{b}+N^{2}v_{3}$
$\displaystyle=\mu\Delta\mathtt{b}\,,$

with the horizontal backscatter parameters $b_{i},\,d_{i}>0$, $i=1,\,2$ and
vertical viscosity $\nu\geq 0$ in the diagonal matrix operator, following the
idealised kinetic energy backscatter. Other quantities in (1.3) are the
velocity field $\boldsymbol{v}(t,\boldsymbol{x})\in\mathbb{R}^{3}$ for
$\boldsymbol{x}\in\mathbb{R}^{3}$, $t\geq 0$, pressure and buoyancy
$p(t,\boldsymbol{x}),\,\mathtt{b}(t,\boldsymbol{x})\in\mathbb{R}$, the nonzero
Coriolis parameter $f\in\mathbb{R}\backslash\\{0\\}$, the vertical unit vector
$\boldsymbol{e_{3}}$ and thermal diffusivity $\mu\geq 0$. As usual, the
buoyancy considered here is of the form
$\mathtt{b}(t,\boldsymbol{x})=-g(\rho(t,\boldsymbol{x})-\overline{\rho}(z))/\rho_{0}\in\mathbb{R}$
with fluid density $\rho(t,\boldsymbol{x})\in\mathbb{R}$ and reference density
field $\overline{\rho}(z)$ depending on the vertical space direction $z$ only,
characteristic density $\rho_{0}$ and gravitational acceleration $g$. Then
$N^{2}=-(g/\rho_{0})d\overline{\rho}/dz$ is the Brunt-Väisälä frequency with
stable stratification $d\overline{\rho}/dz<0$.

We focus on stable stratification so that destabilisation develops from
backscatter in the horizontal directions only. Guided by (Prugger and
Rademacher, 2021), we first find barotropic horizontal flows of similar plane
wave type as in the shallow water case, which arise from sinusoidal wave
shapes by means of a radial superposition principle, akin to the shallow water
case, and in addition from an angular superposition principle of plane waves
with the same wave length. Each of these allows for infinite dimensional
invariant subspaces with linear dynamics. In contrast to the shallow water
case, here the nonlinear constraints admit a horizontal flow for any
horizontal wave vector; the pressure can be determined explicitly in all
cases. These flows are in general not fully geostrophically balanced since the
pressure not only compensates the Coriolis force, but also the gradient part
of the nonlinear term, which may occur with angular superposition. In
addition, for anisotropic backscatter (similar to the shallow water case) the
backscatter partially compensates the Coriolis force.
We determine the loci and stability properties of steady states related to
these flows – in particular the unboundedly unstable ones – and also the large
amplitude regime. Unlike the shallow water plane wave flows, we numerically
find that certain superposed horizontal flows have spectral growth rates that
are proportional to the amplitude factor, i.e. arbitrarily strong
instabilities.
In the same spirit we investigate steady and unboundedly growing explicit
flows with vertical structure and coupled buoyancy. These relate to known
flows in absence of backscatter: parallel flow, Kolmogorov flow and
monochromatic inertia gravity waves (MGWs). While the existence of parallel
flows is unaffected by the (horizontal) backscatter, superposition with small
wave number horizontal flows implies unbounded instability of any parallel
flow. However, the situation is more subtle for the other flows. Notably, we
find that in Kolmogorov flows and MGWs, the purely horizontal backscatter
triggers unbounded growth in the vertical velocity. We also identify possible
superpositions of Kolmogorov flows and MGWs. Regarding the large amplitude
regime, for these flows the growth rates cannot be proportional to their
amplitude factors.

We note that the above analytical realisations of backscatter replace the
usual molecular viscosity operator $\nu\Delta$ by operators of the form
$-(d\Delta^{2}+b\Delta)$, familiar from the scalar Kuramoto-Sivashinsky
equations (KS),

$\frac{\partial u}{\partial t}+\frac{1}{2}|\nabla u|^{2}=-\Delta^{2}u-\Delta
u\,,\quad\boldsymbol{x}\in\mathbb{R}^{n}\,,\quad n\leq 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 the references therein.
Differentiating KS yields a system for $\boldsymbol{v}=\nabla u$ with the
fluid transport nonlinearity, thus relating more closely to (1.1) and (1.3) 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 §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.1) 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

$\displaystyle\boldsymbol{v}=\phi^{\prime}(\boldsymbol{k}\cdot\boldsymbol{x})\boldsymbol{k}^{\perp}\,,\quad\eta=\frac{f}{g}\phi(\boldsymbol{k}\cdot\boldsymbol{x})\,,$

for any wave vector $\boldsymbol{k}\in\mathbb{R}^{2}$ and sufficiently smooth
wave shape $\phi$, e.g. (Prugger and Rademacher, 2021). These are also in
geostrophic balance, corresponding to Rossby waves.

For the case of backscatter in (1.1) we seek solutions of the similar form
(1.2) for any wave vector
$\boldsymbol{k}=(k_{1},k_{2})^{\mathsf{T}}\in\mathbb{R}^{2}$ and sufficiently
smooth wave shapes $\psi$ and $\phi$. The time-independence of $\eta$ results
from equation (1.1d), since $\boldsymbol{v}$ is divergence free and the
nonlinear terms vanish in this case. Inserting (1.2) into (1.1c) yields the
linear equation

$\displaystyle\frac{\partial\psi}{\partial
t}\boldsymbol{k}^{\perp}+f\left(\frac{\partial\phi}{\partial\xi}-\psi\right)\boldsymbol{k}=-|\boldsymbol{k}|^{2}\left(\begin{array}[]{cc}d_{1}|\boldsymbol{k}|^{2}\frac{\partial^{4}\psi}{\partial\xi^{4}}+b_{1}\frac{\partial^{2}\psi}{\partial\xi^{2}}&0\\\
0&d_{2}|\boldsymbol{k}|^{2}\frac{\partial^{4}\psi}{\partial\xi^{4}}+b_{2}\frac{\partial^{2}\psi}{\partial\xi^{2}}\end{array}\right)\boldsymbol{k}^{\perp}\,.$

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 $\boldsymbol{k}^{\perp}$ and $\boldsymbol{k}$,
respectively, gives

(2.1a) $\displaystyle\frac{\partial\psi}{\partial t}$
$\displaystyle=-|\boldsymbol{k}|^{2}\left(d_{1}k_{2}^{2}+d_{2}k_{1}^{2}\right)\frac{\partial^{4}\psi}{\partial\xi^{4}}-\left(b_{1}k_{2}^{2}+b_{2}k_{1}^{2}\right)\frac{\partial^{2}\psi}{\partial\xi^{2}}$
(2.1b) $\displaystyle f\left(\frac{\partial\phi}{\partial\xi}-\psi\right)$
$\displaystyle=k_{1}k_{2}\left((d_{1}-d_{2})|\boldsymbol{k}|^{2}\frac{\partial^{4}\psi}{\partial\xi^{4}}+(b_{1}-b_{2})\frac{\partial^{2}\psi}{\partial\xi^{2}}\right)\,.$

We focus on monochromatic solutions, i.e. that contain a single Fourier mode,
and later investigate possible superpositions. Equations (2.1) restrict such
solutions to the form

(2.2) $\displaystyle\boldsymbol{v}=\alpha_{1}e^{\lambda
t}\cos(\boldsymbol{k}\cdot\boldsymbol{x}+\tau)\boldsymbol{k}^{\perp}\,,\quad\eta=\alpha_{2}\frac{f}{g}\sin(\boldsymbol{k}\cdot\boldsymbol{x}+\tau)+s\,,$

with arbitrary shifts $\tau,\,s\in\mathbb{R}$ and the rest of the real
parameters must satisfy

(2.3a) $\displaystyle\lambda$
$\displaystyle=(b_{1}-d_{1}|\boldsymbol{k}|^{2})k_{2}^{2}+(b_{2}-d_{2}|\boldsymbol{k}|^{2})k_{1}^{2}$
(2.3b) $\displaystyle\frac{\alpha_{2}-\alpha_{1}}{\alpha_{1}}f$
$\displaystyle=k_{1}k_{2}\bigl{(}(d_{1}-d_{2})|\boldsymbol{k}|^{2}+b_{2}-b_{1}\bigr{)}$
(2.3c) $\displaystyle\alpha_{2}\cdot\lambda$ $\displaystyle=0\,.$

Specifically, (1.1c) possesses explicit solutions (2.2) with parameters
satisfying (2.3a) and (2.3b), and the time-independence of $\eta$ coming from
(1.1d) requires condition (2.3c), which means $\alpha_{2}$ or $\lambda$ is
zero. In particular, (2.3c) means that these explicit solutions with non-
trivial depth variation $\eta$ are steady. Notably, solutions with
$\lambda>0$, whose existence is studied in §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 $\lambda$, which characterise the exponentially decaying and
growing explicit solutions. In fact, we show in §2.2.1 that the red region
describes a subset of real unstable eigenmodes of the linearisation of the
zero state.

We proceed as follows: in §2.1 we will discuss the sets of solutions in terms
of their wave vectors, the organising parameters and possible superpositions.
In §2.2 we then analyse the unbounded instability and the linear stability of
explicit steady solutions.

(a) $d_{2}=1.0$ and $\alpha_{2}=0.5$

(b) $d_{2}=1.04$ and $\alpha_{2}=-0.5$

(c) $d_{2}=1.04$

Figure 1: We plot the occurrence of explicit solutions (2.2) 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$, the white
bullets mark steady solutions with $\lambda=0$. The black dashed curves mark
solutions of condition (2.3b) only, for $\alpha_{2}\neq 0$ with values as in
the subcaptions, which also solve (2.3a) at the white bullets ($\lambda=0$).
In (c) we mark a line of wave vectors in a fixed direction (gray), and the
growing or decaying solutions (white circles) on it, whose superpositions with
or without the steady state on the gray line again yield explicit solutions.

###  2.1 Sets of solutions and superpositions

In order to analyse the existence of solutions (2.2) to (1.1) in more detail,
it is convenient to write the conditions (2.3a) and (2.3b) in the form

(2.4a) $\displaystyle\lambda$
$\displaystyle=-(d_{1}k_{2}^{2}+d_{2}k_{1}^{2})|\boldsymbol{k}|^{2}+b_{1}k_{2}^{2}+b_{2}k_{1}^{2}$
(2.4b) $\displaystyle 0$
$\displaystyle=\bigl{(}(d_{1}-d_{2})|\boldsymbol{k}|^{2}+b_{2}-b_{1}\bigr{)}k_{1}k_{2}+\sigma\,,$

with real parameter $\sigma=f\frac{\alpha_{1}-\alpha_{2}}{\alpha_{1}}$
describing the relative difference between the amplitudes of the velocity
vector $\boldsymbol{v}$ and the fluid depth variation $\eta$. Steady solutions
satisfy (2.4a) with $\lambda=0$ and it is then natural to view $\sigma$ as an
adjustment, defined by (2.4b), of the relation between the amplitudes
$\alpha_{1},\,\alpha_{2}$ depending in particular on the wave vector
$\boldsymbol{k}$. For the time-dependent case $\lambda\neq 0$ we have
$\sigma=f$, since $\alpha_{2}=0$ is required due to (2.3c), and viewing (2.4a)
as a definition for $\lambda$. The natural free parameter is the wave vector
$\boldsymbol{k}$. The existence and growth or decay properties of solutions of
the form (2.2), as well as the locations of unboundedly unstable steady states
of this kind, are strongly connected with the values of $\sigma$, 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 (2.2) are also solutions, if all wave vectors
$\boldsymbol{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.1) still vanish due to the
orthogonality of wave vectors and flow directions, 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 (1.1c) 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
$\boldsymbol{k}$ satisfying (2.4a) 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 $\boldsymbol{k}=r\boldsymbol{k}^{\prime}$,
$\boldsymbol{k}^{\prime}$ with $|\boldsymbol{k}^{\prime}|=1$ fixed, the right-
hand side of equation (2.4a) is linear in the squared wave vector length
$|\boldsymbol{k}|^{2}=r^{2}$ (after using
$\boldsymbol{k}=r\boldsymbol{k}^{\prime}$ with $|\boldsymbol{k}^{\prime}|=1$
and division by $r^{2}$). Furthermore, for any fixed $\boldsymbol{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 $\lambda$ is positive in the interior of the closed curve of steady
solutions (2.2) (red regions in Figure 1), except for the origin, where
$\lambda=0$, and $\lambda$ is negative outside (blue regions in Figure 1).
In polar coordinates

$\displaystyle\boldsymbol{k}=r\left(\begin{array}[]{c}\cos(\varphi)\\\
\sin(\varphi)\end{array}\right)\,,$

with angle $\varphi\in[0,2\pi)$ and wave number $r\geq 0$, the curve for
explicit non-trivial steady solutions (2.2), i.e. the wave vectors with
$\lambda=0$, is parameterised by the angle $\varphi$ with the wave number
given by

(2.5) $\displaystyle
r=\sqrt{\frac{b_{1}\sin(\varphi)^{2}+b_{2}\cos(\varphi)^{2}}{d_{1}\sin(\varphi)^{2}+d_{2}\cos(\varphi)^{2}}}\,.$

Generally, these steady solutions have different values of $\sigma$, the
relative difference of amplitudes $\alpha_{1}$ and $\alpha_{2}$; recall that
time-dependent explicit solutions ($\lambda\neq 0$) all have the same value
$\sigma=f$. In either case, the explicit solutions (2.2) form a linear space
since their amplitudes only enter into the ratio
$(\alpha_{2}-\alpha_{1})/\alpha_{1}$ (so into $\sigma$), and are therefore
naturally parameterised by an arbitrary amplitude parameter $a\geq 0$ that is
a common factor of both $\alpha_{1}$ and $\alpha_{2}$, and thus does not
change the value of $\sigma$.

In the following we further investigate the conditions (2.4) for the existence
of explicit solutions (2.2). First, we analyse the occurrence and shapes of
the curves defined by (2.4b) and use this to determine the time-dependent
explicit solutions, i.e. $\alpha_{2}=0$, as well as the steady explicit
solutions for which (2.4a) is satisfied with $\lambda=0$. Second, we discuss
the values of $\sigma$, for which steady solutions exist; clearly any steady
solution of the form (2.2) has a corresponding value of $\sigma$. But not
every $\sigma$ admits such a steady solution and the value of $\sigma$ for the
time-dependent solutions ($\lambda\neq 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 (2.2), primarily of the
time-dependent ones with $\alpha_{2}=0$ and $\lambda\neq 0$, we analyse the
shapes of the curves defined by (2.4b). We start with two special cases:
In the isotropic case $b_{1}=b_{2}$ and $d_{1}=d_{2}$, equation (2.4b)
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
$\boldsymbol{k}=0$ or $\boldsymbol{k}$ on the circle with radius
$\sqrt{b_{1}/d_{1}}$ defined by (2.4a) (with $\lambda=0$), see Figure 2(a).
Thus, non-steady solutions ($\lambda\neq 0$) of the form (2.2) arise from
anisotropy in the backscatter.
There is also a special anisotropic case. If $d_{1}\neq d_{2}$, then (2.4b) is
satisfied in the origin and on the circle
$|\boldsymbol{k}|^{2}=\frac{b_{1}-b_{2}}{d_{1}-d_{2}}$ with $\sigma=0$, and
$\lambda$ defined by (2.4a) is always constant on that circle. If additionally
$b_{1}/d_{1}=b_{2}/d_{2}$, then all solutions of (2.4b) on the circle
$|\boldsymbol{k}|^{2}=\frac{b_{1}-b_{2}}{d_{1}-d_{2}}$ also solve (2.4a) with
$\lambda=0$, so all of these give explicit steady solutions, which have
$\sigma=0$. In case $b_{1}/d_{1}\neq b_{2}/d_{2}$ the value of $\sigma$ for
the steady states is not constant, as mentioned above.

(a) $d_{2}=1.0,\,b_{2}=1.5,\,\alpha_{2}=1.0$

(b) $d_{2}=1.0,\,b_{2}=2.2,\,\alpha_{2}=3.0$

(c) $d_{2}=1.04,\,b_{2}=2.2,\,\alpha_{2}=3.0$

Figure 2: Possible structures of solution curves of (2.4b) 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).

It remains to discuss the general anisotropic case, for which we consider the
wave vectors in polar coordinates
$\boldsymbol{k}=r(\cos(\varphi),\sin(\varphi))^{\mathsf{T}}$ as above; by
symmetry of (2.4b) it suffices to take $\varphi\in[0,\pi]$. The special cases
$\varphi\in\\{0,\frac{\pi}{2},\pi\\}$ requires $\sigma=0$, i.e. steady
solutions (since then $\alpha_{2}\neq 0$), and the corresponding wave vectors
are

$\displaystyle\boldsymbol{k}\in\left\\{\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)\right\\}\,.$

We now consider $\varphi\in(0,\pi)\backslash\\{\frac{\pi}{2}\\}$ only. In the
case $d_{1}=d_{2}$ solutions of (2.4b) are

(2.6) $\displaystyle
r=\pm\sqrt{\frac{2\sigma}{(b_{1}-b_{2})\sin(2\varphi)}}\,,\quad\text{for}\quad\mathrm{sgn}(\sigma)=\mathrm{sgn}\bigl{(}(b_{1}-b_{2})\sin(2\varphi)\bigr{)}\,,$

with $\mathrm{sgn}(\cdot)$ the sign function (see, e.g. Figure 2(b)). For
$d_{1}\neq d_{2}$ we get

(2.7) $\displaystyle
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(2\varphi)}}}\,,$

which gives real solutions to (2.4b), if and only if the expressions in the
square roots of (2.7) are non-negative. This means that for fixed angle
$\varphi$ we have two cases:
(1) $b_{1}=b_{2}$ or $\mathrm{sgn}(b_{1}-b_{2})=-\mathrm{sgn}(d_{1}-d_{2})$
requires for at least one solution that

(2.8)
$\displaystyle\mathrm{sgn}(\sigma)=-\mathrm{sgn}\bigl{(}(d_{1}-d_{2})\sin(2\varphi)\bigr{)}\,.$

(2) $b_{1}\neq b_{2}$ and
$\mathrm{sgn}(b_{1}-b_{2})=\mathrm{sgn}(d_{1}-d_{2})$ requires for at least
one solution that

(2.9a) $\displaystyle\frac{(b_{1}-b_{2})^{2}}{8(d_{1}-d_{2})}\sin(2\varphi)$
$\displaystyle\geq\sigma\quad\text{for
}\varphi\in\left(0,\frac{\pi}{2}\right)\left(\text{ if }d_{1}>d_{2}\text{,
otherwise }\varphi\in\left(\frac{\pi}{2},\pi\right)\right)\,,$ (2.9b)
$\displaystyle\frac{(b_{1}-b_{2})^{2}}{8(d_{1}-d_{2})}\sin(2\varphi)$
$\displaystyle\leq\sigma\quad\text{for
}\varphi\in\left(\frac{\pi}{2},\pi\right)\left(\text{ if }d_{1}>d_{2}\text{,
otherwise }\varphi\in\left(0,\frac{\pi}{2}\right)\right)\,.$

Two solutions for a fixed angle $\varphi$ occur if and only if the following
three conditions are satisfied:

(2.10a)
$\displaystyle\mathrm{sgn}(b_{1}-b_{2})=\mathrm{sgn}(d_{1}-d_{2})\,,\quad\text{so
also }b_{1}\neq b_{2}\,,$ (2.10b)
$\displaystyle\mathrm{sgn}(\sigma)=\mathrm{sgn}\bigl{(}(d_{1}-d_{2})\sin(2\varphi)\bigr{)}\,,$
(2.10c)
$\displaystyle\frac{(b_{1}-b_{2})^{2}}{8|d_{1}-d_{2}|}|\sin(2\varphi)|\geq|\sigma|\,.$

In Figure 2 we plot examples, where up to one (Fig. 2(a) and Fig. 2(b)) or up
to two (Fig. 2(c)) solutions of (2.4b) for certain angles $\varphi$ arise for
fixed $\sigma$. The conditions (2.6)-(2.10) thus determine the occurrence and
structures of the solution curves of (2.4b) 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 $\sigma$

We next discuss the occurrence of explicit steady solutions of the form (2.2)
in the anisotropic case in more detail, in particular the values of $\sigma$,
for which explicit steady solutions exist. These are the values of $\sigma$
for which the curves defined by (2.4a) with $\lambda=0$ and (2.4b) intersect,
see Figure 3. Recall that time-dependent explicit solutions (2.2) all have
$\sigma=f$.
To ease computations, we consider a line $k_{2}=mk_{1}$ with slope
$m\in\mathbb{R}$, i.e. $m=\tan(\varphi)$ in the polar coordinates for wave
vector used before. Inserting this into (2.4) gives the values of $k_{1}$ for
which line and curves intersect

(2.11a) $\displaystyle k_{1}^{2}$
$\displaystyle=\frac{b_{1}m^{2}+b_{2}}{(d_{1}m^{2}+d_{2})(1+m^{2})}\qquad\text{and}$
(2.11b) $\displaystyle k_{1}^{2}$
$\displaystyle=\frac{\sigma}{m(b_{1}-b_{2})}\qquad\text{for}\quad d_{1}=d_{2}$
(2.11c) $\displaystyle k_{1}^{2}$
$\displaystyle=\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})(1+m^{2})}{4m(d_{1}-d_{2})^{2}(1+m^{2})^{2}}}\qquad\text{for}\quad
d_{1}\neq d_{2}\,.$

Here (2.11a) is the intersection of the line with the curve defined by (2.4a),
while (2.11b) and (2.11c) the intersection of the line with the curve defined
by (2.4b) in the two cases (see Figure 3(a)).

(a) $\alpha_{2}=-3.5,\,m=-0.4$ $(\sigma=1.35)$

(b) $\alpha_{2}\approx-0.69,\,m=-0.4$ $(\sigma\approx 0.51)$

(c) $\alpha_{2}\approx-1.20,\,m\approx-0.84\qquad\qquad~{}$ $~{}\hskip
14.22636pt(\sigma=\sigma(m_{\pm})\approx 0.66)$

Figure 3:  Solution curves (dashed lines) of (2.4b) in terms of $\sigma$, as
well as their intersections with solution curve of (2.4a) with $\lambda=0$ and
with gray line marking $k_{2}=mk_{1}$. Here we focus on the upper left
quadrant. Denotations as in Figure 1. 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$.

We next choose $\sigma$ 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
(2.11a) equals the right hand side of (2.11b) or (2.11c). In both cases we
find $\sigma=\sigma(m)$ is

(2.12)
$\displaystyle\sigma(m)=\frac{(b_{1}d_{2}-b_{2}d_{1})(b_{1}m^{2}+b_{2})m}{(d_{1}m^{2}+d_{2})^{2}}\,,$

which satisfies $\sigma(-m)=-\sigma(m)$ 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}\neq b_{2}d_{1}$,
we note that $\sigma(m)$ is differentiable and $\sigma(m)\rightarrow 0$ for
$m\rightarrow\pm\infty$, so that it suffices to determine the extrema. The
derivative of $\sigma(m)$ is given by

(2.13)
$\displaystyle\sigma^{\prime}(m)=\frac{b_{1}d_{2}-b_{2}d_{1}}{(d_{1}m^{2}+d_{2})^{3}}\left(-b_{1}d_{1}m^{4}+3(b_{1}d_{2}-b_{2}d_{1})m^{2}+b_{2}d_{2}\right)\,,$

whose roots, and therefore the location of the extrema, are

(2.14) $\displaystyle
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}}}}\,.$

Thus, steady explicit solutions of the form (2.2) exist for
$\sigma\in[\sigma(m_{-}),\sigma(m_{+})]$, if $b_{1}d_{2}-b_{2}d_{1}>0$, and
$\sigma\in[\sigma(m_{+}),\sigma(m_{-})]$ for the case
$b_{1}d_{2}-b_{2}d_{1}<0$. We may interpret the endpoints $\sigma(m_{-})$ and
$\sigma(m_{+})$, where the solution curves of (2.4) with $\lambda=0$ touch
each other (see Figure 3(c)), as bifurcation points of explicit steady
solutions of the form (2.2).

Using $m=\tan(\varphi)$ we equivalently obtain $\sigma$ as a function of the
wave vector angle $\varphi$ that was used above. We plot an example of the
resulting function $\sigma(\varphi)$ in Figure 5.

###  2.2 Stability analysis of steady solutions

We study stability of a steady state $(\boldsymbol{v}_{s},\eta_{s})$ of (1.1)
via the linear operator
$\mathcal{L}=\mathcal{L}(\boldsymbol{v}_{s},\eta_{s})$, which results from
linearising (1.1) in $(\boldsymbol{v}_{s},\eta_{s})$. A spectrum of
$\mathcal{L}$ with positive real part then implies that the steady solution
$(\boldsymbol{v}_{s},\eta_{s})$ is linearly unstable. In the following we will
show that the trivial steady state
$(\boldsymbol{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 (2.2). We first analyse the unbounded
instability of these in the full nonlinear equations (1.1) with respect to
solutions of the form (2.2). 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 (2.2) with small and large amplitudes.

####  2.2.1 Linear stability of trivial steady state

For the trivial homogeneous steady solution
$(\boldsymbol{v}_{s},\eta_{s})\equiv(0,0,0)$ the linearisation of (1.1) is
exactly (1.1) without the nonlinear terms, and the corresponding linear
operator $\mathcal{L}$ is then the remaining right-hand side of (1.1). The
spectrum of $\mathcal{L}$ in this case can be determined by the dispersion
relation

$\displaystyle
d(\lambda,\boldsymbol{k}):=\det(\lambda\mathrm{Id}-\widehat{\mathcal{L}})=0\,,$

with wave vectors
$\boldsymbol{k}=(k_{1},k_{2})^{\mathsf{T}}\in\mathbb{R}^{2}$, temporal rates
$\lambda=\lambda(\boldsymbol{k})\in\mathbb{C}$ and
$\widehat{\mathcal{L}}=\widehat{\mathcal{L}}(\boldsymbol{k})$ the Fourier
transform of $\mathcal{L}$ given by

$\displaystyle\widehat{\mathcal{L}}(\boldsymbol{k})=\begin{pmatrix}-d_{1}|\boldsymbol{k}|^{4}+b_{1}|\boldsymbol{k}|^{2}&f&-\mathrm{i}gk_{1}\\\
-f&-d_{2}|\boldsymbol{k}|^{4}+b_{2}|\boldsymbol{k}|^{2}&-\mathrm{i}gk_{2}\\\
-\mathrm{i}H_{0}k_{1}&-\mathrm{i}H_{0}k_{2}&0\end{pmatrix}\,.$

The dispersion relation is thus explicitly

(2.15) $\displaystyle
d(\lambda,\boldsymbol{k})=\lambda^{3}+c_{2}\lambda^{2}+c_{1}\lambda+c_{0}=0\,,$

with the coefficients

$\displaystyle c_{2}\ $ $\displaystyle:=\
(d_{1}+d_{2})|\boldsymbol{k}|^{4}-(b_{1}+b_{2})|\boldsymbol{k}|^{2}\,,$
$\displaystyle c_{1}\ $ $\displaystyle:=\
(d_{1}|\boldsymbol{k}|^{4}-b_{1}|\boldsymbol{k}|^{2})(d_{2}|\boldsymbol{k}|^{4}-b_{2}|\boldsymbol{k}|^{2})+gH_{0}|\boldsymbol{k}|^{2}+f^{2}\,,$
$\displaystyle c_{0}\ $ $\displaystyle:=\
-gH_{0}|\boldsymbol{k}|^{2}\left((b_{1}-d_{1}|\boldsymbol{k}|^{2})k_{2}^{2}+(b_{2}-d_{2}|\boldsymbol{k}|^{2})k_{1}^{2}\right)\,.$

Recall that the explicit solutions (2.2) with (2.3) solve both (1.1) with and
without the nonlinear terms, since these terms vanish by construction of these
solutions. Therefore, the wave vectors $\boldsymbol{k}$ and growth rates
$\lambda$ of the explicit solutions are in fact real solutions of the
dispersion relation (2.15). In other words, all these explicit solutions are
real eigenmodes of $\mathcal{L}$ and the values $\lambda$ defined by (2.3a)
are the corresponding real elements in the spectrum of $\mathcal{L}$. Thus,
the possible values for $\lambda$ of the explicit solutions (2.2) with (2.3)
directly provide part of the spectrum of $\mathcal{L}$, for instance all
values of $\lambda$ on the white and black curves in Figure 1(c). In
particular, the occurrence of positive growth rates $\lambda$ in (2.3a)
implies that the trivial steady solution
$(\boldsymbol{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 $\mathcal{L}$.

We first note that in $\boldsymbol{k}=(0,0)^{\mathsf{T}}$ the dispersion
relation reduces to

$\displaystyle 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}=-gH_{0}|\boldsymbol{k}|^{2}\beta(\boldsymbol{k})$, where
$\beta(\boldsymbol{k})$ is the expression in brackets in the definition of
$c_{0}$ above. The coefficient $c_{0}$ of the dispersion relation (2.15) is
zero if and only if $\boldsymbol{k}\in\mathbb{R}^{2}$ satisfies
$\beta(\boldsymbol{k})=0$, which means that $\lambda=0$ is in the spectrum of
$\mathcal{L}$ with corresponding eigenmodes having such wave vectors
$\boldsymbol{k}$. Furthermore, $c_{0}$ is negative if and only if
$\beta(\boldsymbol{k})>0$, so according to the dispersion relation (2.15)
there is at least one positive real value $\lambda>0$ for each of these wave
vectors $\boldsymbol{k}$. We notice, that $\beta(\boldsymbol{k})$ is also
exactly the same expression as on the right-hand side of (2.3a) or (2.4a),
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
$\mathcal{L}$ which are positive real. In particular, we conclude that
$(\boldsymbol{v}_{s},\eta_{s})\equiv(0,0,0)$ is linearly unstable for any
choice of parameters with horizontal backscatter. However, the spectrum of
$\mathcal{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 (2.15), but
we will not do this here.

The previous investigation regarding the instability of the trivial flow in
fact shows that the explicit solutions (2.2) of the full nonlinear equations
(1.1) 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.

(a)

(b)

Figure 4:  Information on the spectrum of the linearisation of (1.1) in the
trivial steady solution $(\boldsymbol{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 $\boldsymbol{k}$ of the associated eigenmodes – real part
positive (red region); negative (blue region); black dot and curve correspond
to steady states of (2.2). 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.

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

In the following we show that some of the steady solutions (2.2) can be
unboundedly unstable as well. We consider parameter values such that some
time-dependent explicit solutions (2.2) have positive growth rate $\lambda$,
as in the example of Figure 1. As already shown, steady solutions (2.2) exist
on the whole curve defined by (2.3a) with $\lambda=0$ (see Figure 5 as well as
the black curve in Figure 1(c)). Now superpositions of explicit solutions
(2.2), which have the same wave vector direction (e.g. the intersections of
white or black curves with the gray line in Figure 1(c)), are also explicit
solutions of (1.1). In the case of Figure 1(c) these are in particular a non-
trivial steady solution $(\boldsymbol{v}_{s},\eta_{s})$ (white dot) and an
exponentially growing solution $(\boldsymbol{v}_{g},0)$ (white circle in red
area). Any superposition
$\alpha(\boldsymbol{v}_{s},\eta_{s})+\varepsilon(\boldsymbol{v}_{g},0)$, with
arbitrary $\alpha,\,\varepsilon\in\mathbb{R}$, is also an explicit solution of
(1.1); in particular, $\varepsilon$ can be arbitrarily close to zero. For any
$\varepsilon\neq 0$, the resulting solution is exponentially and unboundedly
growing. Thus, $(\boldsymbol{v}_{s},\eta_{s})$ is an unboundedly unstable
steady solution.

Figure 5:  Examples for unboundedly unstable explicit steady solutions (2.2)
for the same parameter values as in Fig. 1(c). (a) Red arcs mark unboundedly
unstable steady solutions. Gray curves mark explicit time-dependent solutions
with $\alpha_{2}=0$. (b) Graph of $\sigma(\varphi)$. Red parts mark
unboundedly unstable cases; dashed gray line marks the value of the Coriolis
parameter $f$.

This implies the unbounded instability of the explicit steady solutions (2.2)
corresponding to wave vectors on the red arcs in Figure 5; in Figure 1(c)
these are between the intersections of black and white curves. These arcs
connect intersection points of the curve defined by (2.3a) for $\lambda=0$,
with that for time-dependent explicit solutions defined by (2.3b) with
$\alpha_{2}=0$. The instability of the other explicit steady solutions (black
regions in Figure 5) 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 $\varphi$ and first discuss the values of $\sigma(\varphi)$, for which
the corresponding explicit steady solutions (2.2) are unboundedly unstable.
Recall that since time-dependent explicit solutions (2.2) require
$\alpha_{2}=0$, according to (2.3c), these have $\sigma=f$. Thus, steady
solutions with $\sigma(\varphi)=f$ lie at the intersections with the curve of
time-dependent solutions and all steady solutions “between” those with
$\sigma(\varphi)=f$ are unboundedly unstable, since those can be superposed
with growing explicit solutions (cf. red regions in Figure 5). More precisely,
there are at most four angles $\varphi_{j}\in[0,2\pi)$, ordered by size, so
that $\sigma(\varphi_{j})=f$ (cf. Figure 5), and steady solutions (2.2) whose
wave vectors have angles between $\varphi_{1}$ and $\varphi_{2}$, or
$\varphi_{3}$ and $\varphi_{4}$, are unboundedly unstable. Hence, a steady
solution (2.2) is unboundedly unstable if and only if
$\mathrm{sgn}(f)\sigma(\varphi)>|f|$, with its corresponding value
$\sigma(\varphi)$ (see Figure 5).

Figure 6:  Illustration of long-wavelength instabilities of explicit flows.
(a)-(b): Wave vectors
$\boldsymbol{k}=|\boldsymbol{k}|\bigl{(}\cos(\varphi),\sin(\varphi)\bigr{)}^{\mathsf{T}}$
of explicit solutions (2.2) as functions of $\varphi$. Black: steady solutions
$\boldsymbol{k}_{s}(\varphi)$; blue: exponentially decaying; red:
exponentially growing. (c): Growth rate $\lambda$ as a function of the
difference of wave vector lengths from steady solution. Parameters for all
three cases as in Figure 1(c).

Towards the long-wavelength instability, we parameterise the set of steady
solutions by the angle $\varphi$ of their wave vectors. In Figure 6 we plot
for each $\varphi$ the wave vector lengths for which an explicit solution of
the form (2.2) 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 $\varphi$; 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). Conversely, given any small
Floquet-Bloch parameter one can find a value of $\varphi$, so that the
corresponding steady solution is unstable with respect to it. See Figure 6,
where the self-intersection point at the origin shows two such points along
$\varphi$.

We briefly consider some energetic aspects of the explicit solutions. The
kinetic energy density of solutions to (1.1) is given by
$\rm{KE}(\boldsymbol{v},\eta)=\frac{1}{2}(H_{0}+\eta)|\boldsymbol{v}|^{2}$ and
the potential energy density by
$\rm{PE}(\boldsymbol{v},\eta)=\frac{g}{2}(H_{0}+\eta)^{2}$. The superposed
explicit solutions are generally of the form

$\displaystyle\boldsymbol{v}$
$\displaystyle=\alpha_{s}\cos_{s}\boldsymbol{k}^{\perp}+\alpha_{n}e^{\lambda_{n}t}\cos_{n}\boldsymbol{k}^{\perp}+\alpha_{p}e^{\lambda_{p}t}\cos_{p}\boldsymbol{k}^{\perp}\,,$
$\displaystyle\eta$
$\displaystyle=\gamma\frac{f}{g}\sin(\kappa_{s}\boldsymbol{k}\cdot\boldsymbol{x}+\tau_{s})+c\,,$

with $|\boldsymbol{k}|=1$, wave shapes
$\cos_{a}=\cos(\kappa_{a}\boldsymbol{k}\cdot\boldsymbol{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 (2.2) and
(2.3) (compare with steady solutions in red region in Figure 5 and the
possible superpositions with solutions which are decaying or growing in time).
The energy densities of these explicit solutions explicitly read

$\displaystyle\rm{KE}(\boldsymbol{v},\eta)$
$\displaystyle=\frac{1}{2}\bigl{(}\gamma\frac{f}{g}\sin(\kappa_{s}\boldsymbol{k}\cdot\boldsymbol{x}+\tau_{s})+H_{0}+c\bigr{)}\bigl{(}\alpha_{s}\cos_{s}+\alpha_{n}e^{\lambda_{n}t}\cos_{n}+\alpha_{p}e^{\lambda_{p}t}\cos_{p}\bigr{)}^{2}\,,$
$\displaystyle\rm{PE}(\boldsymbol{v},\eta)$
$\displaystyle=\frac{\gamma^{2}f^{2}}{2g}\sin^{2}(\kappa_{s}\boldsymbol{k}\cdot\boldsymbol{x}+\tau_{s})+\gamma
f(H_{0}+c)\sin(\kappa_{s}\boldsymbol{k}\cdot\boldsymbol{x}+\tau_{s})+\frac{g(H_{0}+c)^{2}}{2}\,.$

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}\boldsymbol{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 (2.2) 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 §3, we set up the notation for the more general setting of an evolution
equation with linear term $\boldsymbol{L}$ and bilinear nonlinearity
$\boldsymbol{B}$ given by

$\frac{\partial}{\partial
t}\boldsymbol{u}=\boldsymbol{L}\boldsymbol{u}+\boldsymbol{B}(\boldsymbol{u},\boldsymbol{u})+\nabla
p\,,$

with a pressure $p$, that is trivial for the rotating shallow water equations
(1.1) with $\boldsymbol{u}=(\boldsymbol{v},\eta)$, and otherwise will derive
from the incompressibility constraint $\nabla\cdot\boldsymbol{v}=0$.

We assume there exists a family of steady state solutions
$\boldsymbol{u}=a\boldsymbol{u}_{0}$ with amplitude parameter $a\in\mathbb{R}$
and associated (possibly trivial) pressure $p=ap_{0}$. The spectral stability
of the steady state $a\boldsymbol{u}_{0}$ is determined by the linearised
right-hand side in $a\boldsymbol{u}_{0}$ and thus the solutions to the
generalised eigenvalue problem

$\lambda\boldsymbol{U}=\mathcal{L}_{a}\boldsymbol{U}+\nabla
P\,,\quad\mathcal{L}_{a}:=\boldsymbol{L}+a\boldsymbol{L}_{0}\,,\quad\boldsymbol{L}_{0}:=\boldsymbol{B}(\cdot,\boldsymbol{u}_{0})+\boldsymbol{B}(\boldsymbol{u}_{0},\cdot)\,,$

with eigenvalue parameter $\lambda\in\mathbb{C}$, eigenmode
$\boldsymbol{U}=(\boldsymbol{v},\eta)$ and $P$ either trivial or determined by
the linearised constraint $\nabla\cdot\boldsymbol{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 $\boldsymbol{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\widetilde{\lambda}$, and set
$P=a\widetilde{P}$. This gives the (generalised) eigenvalue problem

(2.16)
$\displaystyle\widetilde{\lambda}\boldsymbol{U}=(a^{-1}\boldsymbol{L}+\boldsymbol{L}_{0})\boldsymbol{U}+\nabla\widetilde{P}\,.$

The operator of the limiting problem, as $|a|\to\infty$, is
$\boldsymbol{L}_{0}$, and again by continuity of the spectrum, its stability
properties partially predict those of $\mathcal{L}_{a}$ for $|a|\gg 1$. In
particular, an unstable eigenmode of $\boldsymbol{L}_{0}$ implies strongly
unstable eigenmodes of $\mathcal{L}_{a}$ for $|a|\gg 1$, for which the growth
rate $\mathrm{Re}(\lambda)$ is proportional to the amplitude $a$ of the steady
solution. However, eigenmodes of $\mathcal{L}_{a}$ for which $\lambda$ is not
proportional to $a$ will move to the origin in the scaled operator as
$|a|\to\infty$, and thus contribute to the kernel of $\boldsymbol{L}_{0}$. In
particular, $a\boldsymbol{u}_{0}$ may be unstable for all $a$ even though
$\boldsymbol{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 $\lambda$ 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 $\mathcal{L}_{a}$
satisfy $\widetilde{\lambda}=\lambda/a\to 0$ as $a\to\infty$.

Hence, we consider the limiting problem

$\widetilde{\lambda}\boldsymbol{U}=\boldsymbol{L}_{0}\boldsymbol{U}+\nabla\widetilde{P}\,,$

whose spectral properties do not seem to be known analytically for the
explicit flows we are concerned with. Specifically, for (2.2) we have
$\boldsymbol{u}=(\boldsymbol{v},\eta)$, the bilinear form is

$\displaystyle\boldsymbol{B}\Big{(}(\boldsymbol{v}_{1},\eta_{1}),(\boldsymbol{v}_{2},\eta_{2})\Big{)}=\begin{pmatrix}-(\boldsymbol{v}_{1}\cdot\nabla)\boldsymbol{v}_{2}\\\
-(\boldsymbol{v}_{1}\cdot\nabla)\eta_{2}-\eta_{1}\nabla\cdot\boldsymbol{v}_{2}\end{pmatrix}\,,$

and the steady state family is generated by
$\boldsymbol{u}_{0}=(\boldsymbol{v}_{0},\eta_{0})^{\mathsf{T}}$ from (2.2),
i.e. with $\xi=\boldsymbol{k}\cdot\boldsymbol{x}$,

$\displaystyle\boldsymbol{v}_{0}=\cos(\xi)\boldsymbol{k}^{\perp}\,,\quad\eta_{0}=\alpha_{2}\sin(\xi)+s\,,$

where $\alpha_{2}$ is chosen so that (2.3b) holds. We are then interested in
the spectrum of the operator $\boldsymbol{L}_{0}$ defined by

$\displaystyle\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\eta\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}(\boldsymbol{v}_{0}\cdot\nabla)\boldsymbol{v}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}_{0}\hfill\\\
(\boldsymbol{v}_{0}\cdot\nabla)\eta+(\boldsymbol{v}\cdot\nabla)\eta_{0}+\eta_{0}\nabla\cdot\boldsymbol{v}+\eta\nabla\cdot\boldsymbol{v}_{0}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\cos(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\boldsymbol{v}-(\boldsymbol{v}\cdot\boldsymbol{k})\sin(\xi)\boldsymbol{k}^{\perp}\hfill\\\
\cos(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\eta+\alpha_{2}(\boldsymbol{v}\cdot\boldsymbol{k})\cos(\xi)+(\alpha_{2}\sin(\xi)+s)(\nabla\cdot\boldsymbol{v})\end{pmatrix}\,.$

We immediately note that the kernel of $\boldsymbol{L}_{0}$ is infinite
dimensional: any perturbation $\boldsymbol{v},\,\eta$ of the same form as the
steady flow $\boldsymbol{u}_{0}$, i.e.
$\boldsymbol{v}(\boldsymbol{x})=\phi_{1}(\xi)\boldsymbol{k}^{\perp}$,
$\eta(\boldsymbol{x})=\phi_{2}(\xi)$ with arbitrary $\phi_{1}$ and $\phi_{2}$,
lies in the kernel, since
$(\boldsymbol{k}^{\perp}\cdot\nabla)(\boldsymbol{v},\eta)=0$ as well as
$\boldsymbol{v}\cdot\boldsymbol{k}=0$ and $\nabla\cdot\boldsymbol{v}=0$ (in
fact, one can show that here $\boldsymbol{v}$ can be an arbitrary function of
$\xi$). Next, we show that the spectrum of $\boldsymbol{L}_{0}$ is purely
imaginary.

First let $\alpha_{2}=s=0$, i.e. the steady state $\boldsymbol{u}_{0}$ is in
the intersection of the sets of steady and time-dependent solutions as in
Figure 1, so that

$\displaystyle\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\eta\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\cos(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\boldsymbol{v}-(\boldsymbol{v}\cdot\boldsymbol{k})\sin(\xi)\boldsymbol{k}^{\perp}\\\
\cos(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\eta\hfill\end{pmatrix}\,.$

It is a diagonal operator where $\boldsymbol{v}$ and $\eta$ are decoupled. Let
us change coordinates to $\xi=k_{1}x+k_{2}y$, $\zeta=-k_{2}x+k_{1}y$. Then
$\nabla$ becomes
$(k_{1}\partial_{\xi}-k_{2}\partial_{\zeta},k_{2}\partial_{\xi}+k_{1}\partial_{\zeta})$
so that $\boldsymbol{k}^{\perp}\cdot\nabla$ turns into $K\partial_{\zeta}$,
where $K:=k_{1}^{2}+k_{2}^{2}$. We thus obtain the operator

$\displaystyle\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\eta\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\cos(\xi)K\partial_{\zeta}\boldsymbol{v}-(\boldsymbol{v}\cdot\boldsymbol{k})\sin(\xi)\boldsymbol{k}^{\perp}\\\
\cos(\xi)K\partial_{\zeta}\eta\hfill\end{pmatrix}\,,$

whose Fourier transform with respect to $\zeta$ with wave number parameter
$\vartheta$ read

$\displaystyle\widehat{\boldsymbol{L}}_{0}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}i\cos(\xi)K\vartheta\hat{\boldsymbol{v}}-(\hat{\boldsymbol{v}}\cdot\boldsymbol{k})\sin(\xi)\boldsymbol{k}^{\perp}\\\
i\cos(\xi)K\vartheta\hat{\eta}\hfill\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}i\cos(\xi)K\vartheta\mathrm{Id}-\sin(\xi)A(\boldsymbol{k})&0\\\
0&i\cos(\xi)K\vartheta\end{pmatrix}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}\,,$

where $A(\boldsymbol{k})=\begin{pmatrix}-k_{1}k_{2}&-k_{2}^{2}\\\
k_{1}^{2}&k_{1}k_{2}\end{pmatrix}$.

The lower right entry, which corresponds to $\eta$, is a multiplication
operator by $i\cos(\xi)K\vartheta$ and so its spectrum is the range of this
function, which is $iK\vartheta[-1,1]\subset i\mathbb{R}$. Since this
multiplication operator appears in the upper left entry as multiplying the
identity, which commutes with any matrix, $A(\boldsymbol{k})$ 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\neq 0$ and writing $\boldsymbol{v}=(u,v)$ we analogously
obtain the transformed operator

$\displaystyle\widehat{\boldsymbol{L}}_{0}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}i\cos(\xi)K\vartheta\hat{\boldsymbol{v}}-(\hat{\boldsymbol{v}}\cdot\boldsymbol{k})\sin(\xi)\boldsymbol{k}^{\perp}\hfill\\\
i\cos(\xi)K\vartheta\hat{\eta}+\mathcal{B}\hat{\boldsymbol{v}}\hfill\end{pmatrix}=\begin{pmatrix}\mathcal{A}_{1}&0\\\
\mathcal{B}&\mathcal{A}_{2}\end{pmatrix}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}\,,$

with

$\displaystyle\mathcal{A}_{1}$
$\displaystyle:=-i\cos(\xi)K\vartheta\mathrm{Id}+\sin(\xi)A(\boldsymbol{k})$
$\displaystyle\mathcal{A}_{2}$ $\displaystyle:=-i\cos(\xi)K\vartheta$
$\displaystyle\mathcal{B}\hat{\boldsymbol{v}}$
$\displaystyle:=-\alpha_{2}(\hat{\boldsymbol{v}}\cdot\boldsymbol{k})\cos(\xi)+(\alpha_{2}\sin(\xi)+s)(k_{1}\partial_{\xi}\hat{u}+k_{2}\partial_{\xi}\hat{v}+i\vartheta(k_{1}\hat{v}-k_{2}\hat{u}))\,.$

If $\lambda$ lies in the resolvent set of both operators on the diagonal
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ (by the above this includes any non-
purely imaginary value), then $\lambda$ is also in the resolvent set of the
present $\widehat{\boldsymbol{L}}_{0}$, since

$\displaystyle(\widehat{\boldsymbol{L}}_{0}-\lambda)^{-1}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}(\mathcal{A}_{1}-\lambda)^{-1}&0\\\
-(\mathcal{A}_{2}-\lambda)^{-1}\mathcal{B}(\mathcal{A}_{1}-\lambda)^{-1}&(\mathcal{A}_{2}-\lambda)^{-1}\end{pmatrix}\begin{pmatrix}\hat{\boldsymbol{v}}\\\
\hat{\eta}\end{pmatrix}\,.$

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.

Figure 7:  Shown are the eigenvalues with real part larger than $-0.1$ of an
approximation of $\mathcal{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
$\boldsymbol{k}\approx(-1.4,0.35)$. In particular, the solution is already
unstable at marginal instability with respect to the explicit modes.

##  3 Rotating Boussinesq equations with backscatter

We now turn to the study of various explicit solutions in the rotating
Boussinesq equations augmented with backscatter (1.3). To ease notation, we
write these in the form

(3.1a) $\displaystyle\frac{\partial\boldsymbol{v}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}+f\boldsymbol{e_{3}}\times\boldsymbol{v}+\nabla
p-\boldsymbol{e_{3}}\mathtt{b}$
$\displaystyle=-\bigl{(}\mathrm{diag}(d_{1},d_{2},d_{3})\Delta^{2}+\mathrm{diag}(b_{1},b_{2},b_{3})\Delta\bigr{)}\boldsymbol{v}$
(3.1b) $\displaystyle\nabla\cdot\boldsymbol{v}$ $\displaystyle=0$ (3.1c)
$\displaystyle\frac{\partial\mathtt{b}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\mathtt{b}+N^{2}v_{3}$
$\displaystyle=\mu\Delta\mathtt{b}\,,$

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\leq 0$, and stable
stratification $N^{2}>0$. For comparison and illustration we also discuss
briefly the usual horizontal viscous or inviscid cases $d_{1}=d_{2}=0$,
$b_{1},\,b_{2}\leq 0$, unstable stratification $N^{2}<0$, and artificial
vertical backscatter $d_{3},\,b_{3}>0$. As in §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 §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 $\boldsymbol{v}$ that is
independent of the vertical coordinate $z$ and has $v_{3}\equiv 0$, as well as
a horizontally independent buoyancy $\mathtt{b}=\mathtt{b}(t,z)$. This ansatz
yields the reduced equations

(3.2a) $\displaystyle\frac{\partial\boldsymbol{v}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}+f\boldsymbol{v}^{\perp}+\nabla\tilde{p}$
$\displaystyle=-\mathrm{diag}\bigl{(}d_{1}\Delta+b_{1},d_{2}\Delta+b_{2}\bigr{)}\Delta\boldsymbol{v}$
(3.2b) $\displaystyle\nabla\cdot\boldsymbol{v}$ $\displaystyle=0$ (3.2c)
$\displaystyle\frac{\partial\mathtt{b}}{\partial t}$
$\displaystyle=\mu\frac{\partial^{2}}{\partial z^{2}}\mathtt{b}\,,$

with gradient and Laplacian for the horizontal directions
$\boldsymbol{x}=(x,y)^{\mathsf{T}}$,
$\boldsymbol{v}=\boldsymbol{v}(t,\boldsymbol{x})\in\mathbb{R}^{2}$,
$\mathtt{b}=\mathtt{b}(t,z)$ and $p=\tilde{p}+B(t,z)$, where $\frac{\partial
B(t,z)}{\partial z}=\mathtt{b}(t,z)$. Here the buoyancy is decoupled from the
velocity field and determined by the linear heat equation (3.2c); on the
idealized 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 (1.1d) to be replaced by the
incompressibility condition (3.2b). 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 (2.2) for
$\boldsymbol{v}$ satisfies (3.2b) and can readily be adjusted to solutions of
(3.2a). However, there are no additional a priori constraints for the pressure
akin to condition (1.1d) or (2.3c), 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 (3.2).

The resulting set of explicit solutions of (3.2a) for which the nonlinear
terms vanish can be identified by the following ansatz for wave shape $\psi$
and pressure profile $\phi$

(3.3)
$\displaystyle\boldsymbol{v}=\psi(t,\boldsymbol{k}\cdot\boldsymbol{x})\boldsymbol{k}^{\perp}\,,\quad\tilde{p}=f\phi(t,\boldsymbol{k}\cdot\boldsymbol{x})\,,$

where $\boldsymbol{k}\in\mathbb{R}^{2}$ and without loss of generality
$|\boldsymbol{k}|=1$ by the freedom in choosing $\psi$ and $\phi$.
Substitution into (3.2a) gives the linear equations for $\psi$ and $\phi$,

(3.4a) $\displaystyle\frac{\partial\psi}{\partial t}$
$\displaystyle=-(d_{1}k_{2}^{2}+d_{2}k_{1}^{2})\frac{\partial^{4}\psi}{\partial\xi^{4}}-(b_{1}k_{2}^{2}+b_{2}k_{1}^{2})\frac{\partial^{2}\psi}{\partial\xi^{2}}$
(3.4b) $\displaystyle\frac{\partial\phi}{\partial\xi}$
$\displaystyle=\frac{k_{1}k_{2}}{f}\left((d_{1}-d_{2})\frac{\partial^{4}\psi}{\partial\xi^{4}}+(b_{1}-b_{2})\frac{\partial^{2}\psi}{\partial\xi^{2}}\right)+\psi\,.$

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 (3.4b) for the pressure
shape $\phi$ allows the pressure gradient to not only compensate the full
Coriolis term, but also part of the backscatter terms. In particular, the
velocity field (3.3) with (3.4) is in general not geostrophically balanced.

####  3.1.1 Superposition principles

The general wave shape $\psi$ in (3.3) also contains the superpositions of
arbitrary many sinusoidal waves in the same wave vector direction
$\boldsymbol{k}$ and any wave number $|\boldsymbol{k}|$. This is possible in
the Boussinesq equations, since $\nabla p$ in the momentum equation (3.1a) is
not further constrained, unlike (1.1d) for $\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 $\boldsymbol{v}$
in the form (2.2) in another way, namely with different wave vector
directions, but the same wave number, cf. (Prugger and Rademacher, 2021). For
the decoupled momentum equation (3.2a), and finite superposition of $N$ waves
with arbitary $N\in\mathbb{N}$, the resulting superposed sinusoidal explicit
solutions of a form similar to (3.3) are given by

(3.5a) $\displaystyle\boldsymbol{v}(t,\boldsymbol{x})$
$\displaystyle=\sum_{i=1}^{N}e^{\lambda_{i}t}\psi_{i}\boldsymbol{k}^{\perp}_{i}\qquad\mbox{with}\quad\psi_{i}=\alpha_{i}\sin(\boldsymbol{k}_{i}\cdot\boldsymbol{x}+\tau_{i})\,,\quad
1\leq i\leq N\,,$ (3.5b) $\displaystyle\tilde{p}(t,\boldsymbol{x})$
$\displaystyle=-\sum_{i=1}^{N}\sum_{j=i+1}^{N}e^{(\lambda_{i}+\lambda_{j})t}\Bigl{(}(\boldsymbol{k}_{i}\cdot\boldsymbol{k}_{j})\psi_{i}\psi_{j}+s^{2}\frac{\partial\psi_{i}}{\partial\xi}\frac{\partial\psi_{j}}{\partial\xi}\Bigr{)}-f\sum_{i=1}^{N}\frac{\gamma_{i}}{\alpha_{i}}e^{\lambda_{i}t}\frac{\partial\psi_{i}}{\partial\xi}\,,$

for any fixed $s>0$ and arbitrary $\alpha_{i}\in\mathbb{R}\setminus\\{0\\}$,
$\tau_{i}\in\mathbb{R}$ and
$\boldsymbol{k}_{i}=(k_{i,1},k_{i,2})^{\mathsf{T}}\in\mathbb{R}^{2}$ with
$|\boldsymbol{k}_{i}|=s$ for any $1\leq i\leq N$. Here, each $\lambda_{i}$ and
$\gamma_{i}$ is defined by

(3.6a) $\displaystyle\lambda_{i}$
$\displaystyle=(b_{1}-d_{1}s^{2})k_{i,2}^{2}+(b_{2}-d_{2}s^{2})k_{i,1}^{2}\,,$
(3.6b) $\displaystyle f\frac{\gamma_{i}-\alpha_{i}}{\alpha_{i}}$
$\displaystyle=\bigl{(}(d_{1}-d_{2})s^{2}+b_{2}-b_{1}\bigr{)}k_{i,1}k_{i,2}\,,$

in order to solve (3.2a). Since each wave in (3.5a) is divergence free, the
whole superposed velocity $\boldsymbol{v}$ solves (3.2b) and thus is an
explicit solution of (3.2).
It is also possible to superpose explicit solutions of the form (3.5) by
integrating over the whole circle
$S_{s}:=\\{\boldsymbol{k}\in\mathbb{R}^{2}\mid|\boldsymbol{k}|=s\\}$ for any
fixed $s>0$. The exact form is then

(3.7a) $\displaystyle\boldsymbol{v}(t,\boldsymbol{x})$
$\displaystyle=\int\displaylimits_{S_{s}}\alpha(\boldsymbol{k})e^{\lambda(\boldsymbol{k})t}\sin\bigl{(}\boldsymbol{k}\cdot\boldsymbol{x}+\tau(\boldsymbol{k})\bigr{)}\boldsymbol{k}^{\perp}d\boldsymbol{k}$
(3.7b)
$\displaystyle\begin{split}\tilde{p}(t,\boldsymbol{x})&=-s^{4}\int_{0}^{2\pi}\int_{\varphi_{1}}^{2\pi}\Bigl{(}\cos\bigl{(}\xi_{1}\bigr{)}\cos\bigl{(}\xi_{2}\bigr{)}+\cos(\varphi_{1}-\varphi_{2})\sin\bigl{(}\xi_{1}\bigr{)}\sin\bigl{(}\xi_{2}\bigr{)}\Bigr{)}\\\
&\hskip
12.80373pt\cdot\alpha_{1}\alpha_{2}e^{(\lambda_{1}+\lambda_{2})t}d\varphi_{2}d\varphi_{1}-f\int\displaylimits_{S_{s}}\gamma(\boldsymbol{k})e^{\lambda(\boldsymbol{k})t}\cos\bigl{(}\boldsymbol{k}\cdot\boldsymbol{x}+\tau(\boldsymbol{k})\bigr{)}d\boldsymbol{k}\,,\end{split}$

where, for $i=1,\,2$, we set $\alpha_{i}:=\alpha(\boldsymbol{k}_{i})$ with
$\boldsymbol{k}_{i}:=s\bigl{(}\cos(\varphi_{i}),\sin(\varphi_{i})\bigr{)}^{\mathsf{T}}$,
$\lambda_{i}:=\lambda(\boldsymbol{k}_{i})$ and
$\xi_{i}:=\boldsymbol{k}_{i}\cdot\boldsymbol{x}+\tau(\boldsymbol{k}_{i})$, for
all $0\leq\varphi_{i}<2\pi$. Sufficient for the convergence of the integrals
is $\alpha\in L^{1}(S_{s})$, $\tau\in L^{\infty}(S_{s})$ and for almost all
$\boldsymbol{k}\in S_{s}$ we require

(3.8a) $\displaystyle\alpha(\boldsymbol{k})\lambda(\boldsymbol{k})$
$\displaystyle=\alpha(\boldsymbol{k})\bigl{(}(b_{1}-d_{1}s^{2})k_{2}^{2}+(b_{2}-d_{2}s^{2})k_{1}^{2}\bigr{)}\,,$
(3.8b) $\displaystyle
f\bigl{(}\gamma(\boldsymbol{k})-\alpha(\boldsymbol{k})\bigr{)}$
$\displaystyle=\alpha(\boldsymbol{k})k_{1}k_{2}\bigl{(}(d_{1}-d_{2})s^{2}+b_{2}-b_{1}\bigr{)}\,,$

corresponding to (3.6) if $\alpha(\boldsymbol{k})\neq 0$.

The explicit solutions (3.5) with (3.6) differ from (3.3) with (3.4), 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 $\boldsymbol{v}$ in (3.5) and
$|\boldsymbol{k}_{i}|=s$ for all $1\leq i\leq 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 (3.5b). The same holds for the
solutions (3.7) with (3.8) 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 (3.5) with (3.6), as well as (3.7) with
(3.8), 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 (3.6b) and
(3.8b) respectively. The reason is that the pressure gradient in the momentum
equation can additionally compensate a part of the backscatter terms as well.
Second, the growth rates $\lambda_{i}$ and $\lambda(\boldsymbol{k})$ 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}\neq d_{2}$ or $b_{1}\neq 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: (3.3) expresses a radial superposition principle of flows in the same
wave vector direction, and (3.5) 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 §2.2.2, the possible superpositions of explicit solutions, as
contained in (3.3), 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 (3.6a).

Since $b_{1},\,b_{2},\,d_{1},\,d_{2}>0$ in (3.6a) 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
$\boldsymbol{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 (3.3): Such steady explicit solutions also arise from (3.5) with
$N=1$ and $\lambda_{1}=0$ in (3.6a). Since the definition of $\lambda_{i}$ in
(3.6a) is exactly the same as in (2.4a), we can use here the results from
§2.1.2 about the growth rate $\lambda$ 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, $\lambda$ is positive in the interior of this closed curve,
except $\lambda=0$ at the origin $\boldsymbol{k}=(0,0)^{\mathsf{T}}$, and
$\lambda$ is negative outside the closed curve. Thus, single-wave steady
solutions (3.3) 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}\neq 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$ (3.6a) is linear in $\cos(\theta)^{2}$ for
$\boldsymbol{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 (3.5) 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}\neq b_{2}/d_{2}$ we
next appeal to the angular superposition principle. First, we note that for
$b_{1}/d_{1}\neq b_{2}/d_{2}$ there is a unique (up to reversing orientation)
longest wave vector $\boldsymbol{k}_{\mathrm{max}}$ with length
$s_{\mathrm{max}}$, so that (3.6a) is satisfied with $\lambda_{i}=0$. See
Figure 8(a) for a typical example. Indeed, due to the aformentioned
structures, $k_{\mathrm{max}}$ is along an axis, though the set of
$\boldsymbol{k}$ with $\lambda=0$ need not be convex. Hence, any steady
superposed solution (3.5) with $s=s_{\mathrm{max}}$ is built from
$\pm\boldsymbol{k}_{\mathrm{max}}$, 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{min}}$.
Second, we consider steady superposed solutions (3.5) with
$s_{\mathrm{min}}<s<s_{\mathrm{max}}$ 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 (3.5) 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 $\lambda$
defined by (3.6a) is positive for the length $s<s_{\mathrm{max}}$, at least
for the wave vector $s/s_{\mathrm{max}}\boldsymbol{k}_{\mathrm{max}}$, we
again have unbounded instability with respect to a range of modes, here on the
same scale. We note that for $s=s_{\mathrm{min}}$ radial and angular
superposition both give ranges of modes with unbounded growth.

In case $b_{1}/d_{1}=b_{2}/d_{2}$ we cannot infer the instability of steady
superposed solutions (3.5) using the above superposition principles, since the
wave vectors of non-trivial steady solutions
$\\{\boldsymbol{k}\in\mathbb{R}^{2}\mid\lambda(\boldsymbol{k})=0\\}$, with
$\lambda(\boldsymbol{k})$ from (3.6a), form a circle with radius
$\sqrt{b_{1}/d_{1}}$, which means they all have the same length. Thus,
explicit solutions (3.5) with $s=\sqrt{b_{1}/d_{1}}$ consist only of steady
solutions; there are also no steady solutions of the form (3.3) for other wave
numbers. However, in some cases unbounded instability still follows from
angular superposition with unboundedly growing parallel flows, cf. §3.2.1.

(a) Case $s=s_{\mathrm{max}}=\sqrt{b_{2}/d_{2}}$ for steady solutions with
largest possible wave number.

(b) Case $s\approx 1.20<s_{\mathrm{max}}$ for steady solutions with different
wave vector directions.

Figure 8:  Signs of $\lambda$ as defined in (3.6a) (red: $\lambda>0$, blue:
$\lambda<0$, black: $\lambda=0$) and possible wave vectors for (3.5) 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}\neq
b_{2}/d_{2}$.

####  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 §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 (3.2) in the trivial steady flow $\boldsymbol{v}\equiv 0$, here we restrict
attention to instability with respect to eigenmodes of the form of the
horizontal flows (3.3), which means solutions to (3.4a). The Fourier transform
of (3.4a) yields the dispersion relation for perturbation wave vector
$\boldsymbol{k}\in\mathbb{R}^{2}$ and temporal rate $\lambda\in\mathbb{C}$,

$d_{\psi}(\lambda,\boldsymbol{k})=-(d_{1}k_{2}^{2}+d_{2}k_{1}^{2})|\boldsymbol{k}|^{2}+b_{1}k_{2}^{2}+b_{2}k_{1}^{2}-\lambda=0\,,$

which is of course equivalent to (3.6a) with $\lambda=\lambda_{i}$ and
$\boldsymbol{k}=\boldsymbol{k}_{i}$. The above discussion for steady states
built from single direction wave vectors $\pm\boldsymbol{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 (3.2) 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
$|\boldsymbol{k}|=s$.

However, analogous to §2.2.3, any (superposed) steady horizontal flow inherit
the instability of any unstable mode in the dispersion relation
$d_{\psi}(\lambda,\boldsymbol{k})=0$ of the zero state for sufficiently small
amplitudes $0<|\alpha_{i}|\ll 1$, $1\leq i\leq N$, though the growth induced
by such modes may be bounded. Recall that the explicit solutions of (3.2) also
satisfy the full rotating Boussinesq equations (3.1), which admit modes that
have vertical structures and are coupled with buoyancy. In §3.2 we will
discuss such explicit solutions of (3.1), 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
$\boldsymbol{v}\equiv 0$, which – in contrast to the horizontal flows – are
also influenced by the Brunt-Väisälä frequency $N^{2}$, thermal diffusivity
$\mu$ and the vertical viscosity. In addition, these imply linear instability
modes of horizontal flows with sufficiently small amplitudes
$0<|\alpha_{i}|\ll 1$, $1\leq i\leq N$.

As to large amplitude(s), where $1\ll|\alpha_{i}|$ for at least one $i$, we
first note that, in the notation of §2.2.3 and with
$\boldsymbol{u}=(\boldsymbol{v},\mathtt{b})$, the bilinear form for the
present case reads

$\displaystyle\boldsymbol{B}\Big{(}(\boldsymbol{v}_{1},\mathtt{b}_{1}),(\boldsymbol{v}_{2},\mathtt{b}_{2})\Big{)}=-\begin{pmatrix}(\boldsymbol{v}_{1}\cdot\nabla)\boldsymbol{v}_{2}\\\
(\boldsymbol{v}_{1}\cdot\nabla)\mathtt{b}_{2}\end{pmatrix}\,.$

The steady state family $\boldsymbol{u}=a\boldsymbol{u}_{0}$ in this case has
$\boldsymbol{u}_{0}=(\boldsymbol{v}_{0},\mathtt{b}_{0})^{\mathsf{T}}$ and

$\displaystyle\boldsymbol{v}_{0}=\sum_{i=1}^{N}\alpha_{i}\sin(\boldsymbol{k}_{i}\cdot
x+\tau_{i})\boldsymbol{k}_{i}^{\perp}\,,$

where we consider in this section only horizontal (wave) vectors
$\boldsymbol{k}:=(\hat{\boldsymbol{k}},0)^{\mathsf{T}}$,
$\boldsymbol{k}^{\perp}:=(\hat{\boldsymbol{k}}^{\perp},0)^{\mathsf{T}}$ for
any $\hat{\boldsymbol{k}}\in\mathbb{R}^{2}$. Here the third component of
$\boldsymbol{v}_{0}$ vanishes, and $\mathtt{b}_{0}$ is an arbitrary constant
solving (3.2c). In order to locate strongly unstable modes, whose growth rates
$\lambda=a\widetilde{\lambda}$ are proportional to the amplitude parameter
$a$, we are concerned with the generalised eigenvalue problem (2.16), as
$|a|\to\infty$, which here reads

(3.9) $\displaystyle\widetilde{\lambda}\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}=\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}+\begin{pmatrix}\nabla\bar{p}\\\
0\end{pmatrix}\,,\quad\nabla\cdot\boldsymbol{v}=0\,,$

with $p=a\bar{p}$ and for the operator $\boldsymbol{L}_{0}$ defined by

$\displaystyle\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}(\boldsymbol{v}_{0}\cdot\nabla)\boldsymbol{v}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}_{0}\hfill\\\
(\boldsymbol{v}_{0}\cdot\nabla)\mathtt{b}\end{pmatrix}\,.$

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=\boldsymbol{k}\cdot\boldsymbol{x}$,
$\zeta=\boldsymbol{k}^{\perp}\cdot\boldsymbol{x}$ so that with a certain wave
vector $\boldsymbol{k}$

$\displaystyle\boldsymbol{v}_{0}=\alpha_{1}\sin(\xi)\boldsymbol{k}^{\perp}+\alpha_{2}\sin(\zeta)\boldsymbol{k}\,.$

We note, that even in the anisotropic case there is at least one wave vector
$\boldsymbol{k}$, such that both $\boldsymbol{k}$ and $\boldsymbol{k}^{\perp}$
correspond to a steady mode. We omit the full proof, and instead explain the
existence of such $\boldsymbol{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 $\boldsymbol{k}$ with
$|\boldsymbol{k}|=s$ exists, for which $\boldsymbol{k}$ and
$\boldsymbol{k}^{\perp}$ each correspond to a single mode steady solution.
This is ensured by the symmetry of the curve defined by (3.6) with
$\lambda_{i}=0$. We then obtain

$\displaystyle\boldsymbol{L}_{0}\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\alpha_{1}\sin(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\boldsymbol{v}+\alpha_{2}\sin(\zeta)(\boldsymbol{k}\cdot\nabla)\boldsymbol{v}+\alpha_{1}\cos(\xi)(\boldsymbol{v}\cdot\boldsymbol{k})\boldsymbol{k}^{\perp}+\alpha_{2}\cos(\zeta)(\boldsymbol{v}\cdot\boldsymbol{k}^{\perp})\boldsymbol{k}\hfill\\\
\alpha_{1}\sin(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)\mathtt{b}+\alpha_{2}\sin(\zeta)(\boldsymbol{k}\cdot\nabla)\mathtt{b}\end{pmatrix}$
$\displaystyle=-\mathrm{diag}(L_{11},L_{22},L_{22})\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}\,,$

with block diagonal matrix operator in which $L_{11}$ is the 2-by-2 matrix
operator

$\displaystyle
L_{11}=L_{22}\mathrm{Id}+\alpha_{1}\cos(\xi)A(\boldsymbol{k})-\alpha_{2}\cos(\zeta)A^{\mathsf{T}}(\boldsymbol{k})\,,$

with $A(\boldsymbol{k})$ as in §2.2.3 (for which
$A^{\mathsf{T}}(\boldsymbol{k})=-A(\boldsymbol{k}^{\perp})$) and

$\displaystyle L_{22}$
$\displaystyle=\alpha_{1}\sin(\xi)(\boldsymbol{k}^{\perp}\cdot\nabla)+\alpha_{2}\sin(\zeta)(\boldsymbol{k}\cdot\nabla)\,.$

Taking the divergence of (3.9) gives, using $\nabla\cdot\boldsymbol{v}=0$, the
linear pressure Poisson equation
$\Delta\bar{p}=\nabla\cdot(L_{1}\boldsymbol{v})$, with block matrix operator
$L_{1}:=\mathrm{diag}(L_{11},L_{22})$. We denote the solution as
$\bar{p}=\Delta^{-1}\nabla\cdot(L_{1}\boldsymbol{v})$ and substitution into
(3.9) yields the eigenvalue problem

(3.10) $\displaystyle\widetilde{\lambda}\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}=\left(\boldsymbol{L}_{0}+\begin{pmatrix}\nabla\Delta^{-1}\nabla\cdot(L_{1}\circ)&0\\\
0&0\end{pmatrix}\right)\begin{pmatrix}\boldsymbol{v}\\\
\mathtt{b}\end{pmatrix}\,,$

in which $\circ$ denotes the slot for $\boldsymbol{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 $-L_{22}$ and $-\widetilde{L}_{1}$
defined by

$\widetilde{L}_{1}:=L_{1}-\nabla\Delta^{-1}\nabla\cdot(L_{1}\circ)\,.$

Since $L_{22}$ is skew-adjoint ($\mathrm{i}L_{22}$ 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 §2.2.3
that the spectrum of $-\widetilde{L}_{1}$ is also purely imaginary, so that
the spectrum of $\boldsymbol{L}_{0}$ is purely imaginary.

Otherwise, if $\alpha_{1}\alpha_{2}\neq 0$, it appears difficult to determine
the spectrum of $\widetilde{L}_{1}$ analytically and we resort to numerical
computations. For this let $(\cdot)_{m}$ denote the projection onto the mode
$\exp(\mathrm{i}\boldsymbol{k}_{m}\cdot\boldsymbol{x})$. Then

(3.11)
$\displaystyle(\widetilde{L}_{1}\boldsymbol{v})_{m}=(L_{1}\boldsymbol{v})_{m}-\frac{\boldsymbol{k}_{m}\cdot(L_{1}\boldsymbol{v})_{m}}{|\boldsymbol{k}_{m}|^{2}}\boldsymbol{k}_{m}=(\mathrm{Id}-B(\boldsymbol{k}_{m}))(L_{1}\boldsymbol{v})_{m}\,,$

with suitable matrix $B(\boldsymbol{k}_{m})$. 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 $\boldsymbol{L}_{0}$.

Figure 9:  Shown is an approximation of part of the spectrum of
$\widetilde{L}_{1}$. Using (3.11) we reduced to two-dimensional wave vectors
by fixing the third component of $\boldsymbol{k}_{m}$ at zero. Here
$\boldsymbol{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 (3.5) 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 $\boldsymbol{L}_{0}$ is purely
imaginary, as in §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 (3.1) 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 is non-
trivial. Here we investigate parallel flows, Kolmogorov flows and
monochromatic inertia gravity waves. 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,

(3.12)
$\displaystyle\boldsymbol{v}(t,\boldsymbol{x})=w(t,x,y)\boldsymbol{e_{3}}\,,\quad\mathtt{b}(t,\boldsymbol{x})=\tilde{\mathtt{b}}(t,x,y)\,,\quad
p(t,\boldsymbol{x})=\tilde{p}(t)z\,,$

where $w$ and $\tilde{\mathtt{b}}$ satisfy (with horizontal Laplacian and bi-
Laplacian)

(3.13a) $\displaystyle\frac{\partial w}{\partial
t}+\left(d_{3}\Delta^{2}+b_{3}\Delta\right)w+\tilde{p}$
$\displaystyle=\tilde{\mathtt{b}}$ (3.13b)
$\displaystyle\frac{\partial\tilde{\mathtt{b}}}{\partial
t}-\mu\Delta\tilde{\mathtt{b}}$ $\displaystyle=-N^{2}w\,.$

Recall that kinetic energy backscatter, which has $d_{3}=0$ and $b_{3}\leq 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 (3.3) that have zero buoyancy, if their wave vector directions
$\boldsymbol{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 (3.3) with the same wave
vector $\boldsymbol{k}$ and small enough wave number $|\boldsymbol{k}|$.

Existence and dynamics of parallel flows can be inferred from the dispersion
relation of the linear equations (3.13). By Fourier transformation with wave
vector $\boldsymbol{k}$ and growth rate $\lambda$ this is given by

$\displaystyle\det(\lambda\mathrm{Id}-\widehat{\boldsymbol{L}})=0\,,\quad\widehat{\boldsymbol{L}}:=\begin{pmatrix}-d_{3}|\boldsymbol{k}|^{4}+b_{3}|\boldsymbol{k}|^{2}&1\\\
-N^{2}&-\mu|\boldsymbol{k}|^{2}\end{pmatrix}\,,$

or equivalently as the characteristic polynomial

(3.14) $\displaystyle\lambda^{2}+c_{1}\lambda+c_{0}=0\,,$

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:=|\boldsymbol{k}|^{2}\geq 0$.
Steady solutions to (3.13) require constant $\tilde{p}$ and consist of Fourier
modes with $\boldsymbol{k}\neq 0$ that solve (3.14) with $\lambda=0$, i.e.
$c_{0}=0$. Growing spatially non-constant solutions to (3.13) exist if and
only if (3.14) possesses a root with positive real part and
$\boldsymbol{k}\neq 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 $\boldsymbol{k}\neq 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}/(\mu b_{3})$, 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 §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 $\boldsymbol{L}_{0}$ for steady parallel flow $w_{0}$ is a
lower triangular matrix operator with diagonal entries
$L_{1}:=w_{0}(x,y)\partial_{z}$. Hence, as for the spectrum of
$\boldsymbol{L}_{0}$ in §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
$L_{22}$ in §3.1.3, and thus the spectrum of $\boldsymbol{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}(K):=(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}\geq
4c_{0}$, which is equivalent to $K<b_{3}/d_{3}$ and
$N^{2}\leq\delta_{3}^{2}(K)/4$, respectively. For $K\in(0,b_{3}/d_{3})$ the
global maximum of $\delta_{3}^{2}(K)/4$ is $b_{3}^{4}/(64d_{3}^{2})$ at
$K=K_{1}:=b_{3}/(2d_{3})$; its global minimum is zero at $K=0$ and
$K=b_{3}/d_{3}$. Specifically, if $N^{2}\leq b_{3}^{4}/(64d_{3}^{2})$, i.e.
the stability of the stratification is sufficiently weak compared with the
backscatter destabilisation, then $c_{1}^{2}\geq 4c_{0}$ for $K$ in a positive
interval $I_{1}\subset(0,b_{3}/d_{3})$; in particular, $I_{1}=\\{K_{1}\\}$ if
$N^{2}=b_{3}^{4}/(64d_{3}^{2})$. Hence, a parallel flow (3.12) grows
exponentially and unboundedly if it contains a Fourier mode with wave vector
$\boldsymbol{k}$ in the annulus
$\\{\boldsymbol{k}\in\mathbb{R}^{2}\,|\,|\boldsymbol{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}/(27d_{3}^{2})$, 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 (3.12) grows exponentially and unboundedly if it contains a
Fourier mode with wave vector $\boldsymbol{k}$ in the annulus
$\\{\boldsymbol{k}\in\mathbb{R}^{2}\,|\,|\boldsymbol{k}|^{2}\in I_{2}\\}$ for
(not too strongly) stable stratification, or in the disc
$\\{\boldsymbol{k}\in\mathbb{R}^{2}\,|\,|\boldsymbol{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}\geq 0$
also yield such growth if $c_{1}<0$ and $c_{1}^{2}\geq 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
absence of backscatter are the so-called Kolmogorov flows, see e.g. (Balmforth
and Young, 2005), with wave vectors of the form
$\boldsymbol{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

(3.15) $\displaystyle\boldsymbol{v}(t,\boldsymbol{x})=e^{\lambda
t}\cos(\boldsymbol{k}\cdot\boldsymbol{x})\boldsymbol{a}\,,\quad\mathtt{b}(t,\boldsymbol{x})=ce^{\lambda
t}\cos(\boldsymbol{k}\cdot\boldsymbol{x})\,,\quad p(t,\boldsymbol{x})=\gamma
e^{\lambda t}\sin(\boldsymbol{k}\cdot\boldsymbol{x})\,,$

and the flow direction
$\boldsymbol{a}=\alpha(0,1,0)^{\mathsf{T}}+\beta(m,0,k)^{\mathsf{T}}.$

Compared with the Kolmogorov flows without backscatter and rotation from
(Prugger and Rademacher, 2021), here we have rotation ($f\neq 0$), time
dependence ($\lambda\neq 0$) and a nonzero second component of the velocity
direction ($\alpha\neq 0$). A superposition of these Kolmogorov flows with the
horizontal flow solutions (3.3) 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
$\boldsymbol{k}$ have the same direction, as is the superposition with such
monochromatic inertia gravity waves as discussed in §3.2.3.

We next determine the necessary relations for the coefficients of (3.15) in
order to solve the Boussinesq equations. For better readability we define the
following terms resulting from the backscatter and thermal diffusion:

$\displaystyle\delta_{\mu}(k,m)=\mu|\boldsymbol{k}|^{2}\,,\quad\delta_{i}(k,m)=d_{i}|\boldsymbol{k}|^{4}-b_{i}|\boldsymbol{k}|^{2}\,\mbox{
for any }\,1\leq i\leq 3\,,$

where $|\boldsymbol{k}|^{2}=k^{2}+m^{2}$. Upon inserting (3.15) into (3.1) we
find that the coefficients have to satisfy

(3.16) $\displaystyle\begin{pmatrix}-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{pmatrix}\begin{pmatrix}\alpha\\\ \beta\\\
c\\\ \gamma\end{pmatrix}=0\,.$

For $(k,m)=(0,0)$ these require $\alpha=0$ and $c=0$, which is the zero state.
From the second row of the 4-by-4 matrix in (3.16) and $f\neq 0$ we
immediately find that $\alpha=0$ implies $\beta m=0$. In case $\beta=0$, the
Kolmogorov flow (3.15) is the trivial zero solution, and in case $m=0$, (3.15)
is a parallel flow. Hence, we may assume $\alpha\neq 0$. Since (3.16) is a
homogeneous linear system in $(\alpha,\beta,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)\neq(0,0)$ and dividing by
$-|\boldsymbol{k}|^{2}$, this gives

(3.17) $\displaystyle\lambda^{3}+c_{2}\lambda^{2}+c_{1}\lambda+c_{0}=0\,,$

with coefficients

$\displaystyle c_{2}$
$\displaystyle:=(\delta_{2}+\delta_{\mu})+|\boldsymbol{k}|^{-2}(\delta_{3}k^{2}+\delta_{1}m^{2})\,,$
$\displaystyle c_{1}$
$\displaystyle:=\delta_{2}\delta_{\mu}+|\boldsymbol{k}|^{-2}[(\delta_{2}+\delta_{\mu})(\delta_{3}k^{2}+\delta_{1}m^{2})+N^{2}k^{2}+f^{2}m^{2}]\,,$
$\displaystyle c_{0}$
$\displaystyle:=|\boldsymbol{k}|^{-2}[\delta_{2}\delta_{\mu}(\delta_{3}k^{2}+\delta_{1}m^{2})+\delta_{2}N^{2}k^{2}+\delta_{\mu}f^{2}m^{2}]\,.$

Steady Kolmogorov flows and linear stability. The condition $\lambda=0$ for
steady Kolmogorov flow reduces (3.17) to

(3.18)
$\displaystyle\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\,.$

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\\{(0,0)\\}$.
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
$|\boldsymbol{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. Fig 10(b), 10(d). Conversely, we can create steady
Kolmogorov flows for any wave vector $(k,m)\neq(0,0)$ by suitable choice of
$d_{2},\,b_{2}$ such that $\delta_{2}=0$.

Concerning stability, analogous to §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 $\boldsymbol{L}_{0}$ is a triangular block matrix operator,
similar to §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 monochromatic inertia gravity waves discussed in §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 $\beta,\,k\neq 0$
transfer the horizontal backscatter to growing vertical velocity component.
They correspond to positive roots of (3.17), which occur as follows in terms
of the sign of $c_{0}$:

  * (1)

If $c_{0}<0$, then (3.17) has a positive root.

  * (2)

For $c_{0}=0$, (3.17) has a positive root if and only if
$-c_{2}+\sqrt{c_{2}^{2}-4c_{1}}>0$.

  * (3)

For $c_{0}>0$, (3.17) 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}}\leq
0$.

The conditions in (3) imply that the local minimum of the cubic polynomial in
(3.17) 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 (3.15) requires unstable stratification
$N^{2}<0$. Indeed, in this case (3.17) has $c_{0}=c_{2}=0$ and
$c_{1}=|\boldsymbol{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. Fig. 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 Fig. 10(a)), which turn into a triple
steady flow on the boundary, where $c_{0}=c_{1}=c_{2}=0$ (black curves).

(a) $\Lambda=\\{-4,0,0,0,0,0,0,0\\}$

(b) $\Lambda=\\{-4,1,2,0.5,1,0,-1,0\\}$

(c) $\Lambda=\\{-4,1,2,0.5,1,0,-1,3\\}$

(d) $\Lambda=\\{4,1,2,0.5,1,0,-1,0\\}$

(e) $\Lambda=\\{4,1,2,0.5,1,0,-1,3\\}$

Figure 10:  Unboundedly growing Kolmogorov flows occur for wave vectors
$(k,m)$ in the red regions for parameters sets
$\Lambda:=\\{N^{2},d_{1},b_{1},d_{2},b_{2},d_{3},b_{3},\mu\\}$, where the
Coriolis parameter is fixed at $f=1$. Red regions: one of conditions (1), (2)
and (3) is satisfied; blue regions: none of conditions (1), (2) and (3) 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)$.

In the presence of (horizontal) backscatters ($b_{j},\,d_{j}>0$ for $j=1,\,2$)
growth of Kolmogorov flows (3.15) is possible also for stable stratification
$N^{2}>0$. We next focus on case (1) with negative $c_{0}$ and omit details of
cases (2) and (3) with non-negative $c_{0}$. Some examples are plotted in Fig.
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
(3.18), whose leading order term as $(k,m)\to(0,0)$ is
$|\boldsymbol{k}|^{2}(-b_{2}N^{2}k^{2}+\mu f^{2}m^{2})$ and can be negative
only if $N^{2}>0$ and $k\neq 0$. Then, to leading order, growing Kolmogorov
flows occur for $m^{2}/k^{2}<b_{2}N^{2}/(\mu f^{2})$ with $\mu>0$ (cf. Fig
10(e)) and increase of $N^{2}$ or $b_{2}$ enlarges this region in the
$(k,m)$-plane near the origin, while increase of $\mu$ or $f$ shrinks it.
Rewriting the condition as $N^{2}>\mu f^{2}m^{2}/(b_{2}k^{2})$, for a given
wave vector, then on the one hand increase of $\mu$ 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 (3.18) is always positive if $N^{2}<0$, so $c_{0}>0$ near the origin,
which refers to case (3) above; we omit details and just plot an example in
Fig. 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
$|\boldsymbol{k}|^{2}<b_{2}/d_{2}$ if $N^{2}>0$ and $k\neq 0$ (cf. Fig.
10(d)), and $c_{0}<0$ for $|\boldsymbol{k}|^{2}>b_{2}/d_{2}$ if $N^{2}<0$ and
$k\neq 0$ (cf. Fig. 10(b)). For $k=0$, $c_{0}$ is zero so we refer to the case
(2) above for the growing Kolmogorov flows. The numerical computations show
that in the examples of Fig. 10(b) and 10(d) the condition in case (2) is not
satisfied for $k=0$ so that here a Kolmogorov flow (3.15) is not growing. Near
$\boldsymbol{k}=0$, the examples of Fig. 10(b) and 10(c) can be regarded as
perturbations of the example of Fig. 10(a). As before, the red region
indicates existence of growing flows; note that decaying or oscillating ones
may co-exist. Notably, in Fig. 10(b) and 10(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 (3.17) 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|\boldsymbol{k}|^{4}-b|\boldsymbol{k}|^{2}=\delta_{j}$ for
$j=1,\,2,\,3$. With thermal diffusion ($\mu>0$), we then have $c_{0}<0$ for
small wave number $|\boldsymbol{k}|^{2}<b/d$ (i.e. $\delta<0$) if
stratification is sufficiently stable
$N^{2}>-\mu|\boldsymbol{k}|^{2}(f^{2}m^{2}+\delta^{2}|\boldsymbol{k}|^{2})/(k^{2}\delta)>0$,
and for large wave number $|\boldsymbol{k}|^{2}>b/d$ (i.e. $\delta>0$) only
for sufficiently unstable stratification
$N^{2}<-\mu|\boldsymbol{k}|^{2}(f^{2}m^{2}+\delta^{2}|\boldsymbol{k}|^{2})/(k^{2}\delta)<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 $\mu$ 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 Figs. 10(d), 10(e) for $N^{2}>0$, and Figs.
10(b), 10(c) for $N^{2}<0$.

####  3.2.3 Monochromatic inertia gravity waves

The last kind of explicit plane wave-type solutions we are aware of are the
so-called monochromatic inertia gravity waves (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

(3.19a) $\displaystyle\boldsymbol{v}(t,\boldsymbol{x})$
$\displaystyle=\alpha_{1}e^{\lambda
t}\sin(\xi)(0,1,0)^{\mathsf{T}}+\alpha_{2}\omega e^{\lambda
t}\cos(\xi)(-m,0,k)^{\mathsf{T}}\,,$ (3.19b)
$\displaystyle\mathtt{b}(t,\boldsymbol{x})$ $\displaystyle=\beta_{1}e^{\lambda
t}\sin(\xi)+\beta_{2}\omega e^{\lambda t}\cos(\xi)\,,$ (3.19c) $\displaystyle
p(t,\boldsymbol{x})$ $\displaystyle=\gamma_{1}e^{\lambda
t}\cos(\xi)+\gamma_{2}\omega e^{\lambda t}\sin(\xi)\,.$

The conditions for these to be (non-trivial) explicit solutions of (3.1), with
$\boldsymbol{k}=(k,m)^{\mathsf{T}}\in\mathbb{R}^{2}\backslash\\{(0,0)^{\mathsf{T}}\\}$,
in particular depend on $\omega$. In fact, for $\omega=0$ these are Kolmogorov
flows from §3.2.2 with $\beta=0$, and the existence conditions read, after
inserting (3.19) into (3.1),

(3.20) $\displaystyle\begin{pmatrix}\lambda+\delta_{2}&0&0\\\ mf&-k&0\\\
kf&m&|\boldsymbol{k}|^{2}\\\
0&\lambda+\delta_{\mu}&0\end{pmatrix}\begin{pmatrix}\alpha_{1}\\\ \beta_{1}\\\
\gamma_{1}\end{pmatrix}=0\,,$

with
$\delta_{j}=\delta_{j}(\boldsymbol{k}):=d_{j}|\boldsymbol{k}|^{4}-b_{j}|\boldsymbol{k}|^{2}$
for any $1\leq j\leq 3$ and
$\delta_{\mu}=\delta_{\mu}(\boldsymbol{k}):=\mu|\boldsymbol{k}|^{2}$ as above.
For $\omega\neq 0$ the existence conditions are at first the following eight:

(3.21) $\displaystyle\begin{pmatrix}\lambda+\delta_{2}&0&0&0&0&0\\\
mf&-k&0&\omega^{2}|\boldsymbol{k}|^{2}&0&0\\\
kf&m&|\boldsymbol{k}|^{2}&0&0&0\\\ 0&\lambda+\delta_{\mu}&0&0&\omega^{2}&0\\\
1&0&0&mf&0&0\\\ 0&-1&0&kN^{2}&\lambda+\delta_{\mu}&0\\\
0&0&0&\delta_{3}k^{2}+\delta_{1}m^{2}+\lambda|\boldsymbol{k}|^{2}&-k&0\\\
0&0&0&(\delta_{3}-\delta_{1})km&-m&|\boldsymbol{k}|^{2}\end{pmatrix}\begin{pmatrix}\alpha_{1}\\\
\beta_{1}\\\ \gamma_{1}\\\ \alpha_{2}\\\ \beta_{2}\\\
\gamma_{2}\end{pmatrix}=0\,.$

However, most of these can be readily solved directly in terms of the
coefficients. In the following we explicitly determine all non-trivial
solutions (3.19). Afterwards we shortly discuss the stability of MGWs and the
possible superpositions with other types of explicit solutions.

MGWs with steady phase ($\omega=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 $\boldsymbol{v}$, $\mathtt{b}$ and $p$ in
(3.19) vanish. Solutions with

$\displaystyle
k=0\,,\quad\alpha_{1}=0\,,\quad\beta_{1}=-\gamma_{1}m\,,\quad\lambda=-\delta_{\mu}$

have $\boldsymbol{v}\equiv 0$, the pressure depends on the buoyancy only, and
there are no further conditions on the wave vector
$\boldsymbol{k}=(0,m)^{\mathsf{T}}$.

Other solutions with $\omega=0$ satisfy

$\displaystyle k\neq
0\,,\quad\beta_{1}=\alpha_{1}\frac{m}{k}f\,,\quad\gamma_{1}=-\alpha_{1}\frac{f}{k}\,,\quad\lambda=-\delta_{\mu}\,,\quad\delta_{2}=\delta_{\mu}\,,$

so $\boldsymbol{v}\neq 0$ for $\alpha_{1}\neq 0$. The last of these equations
gives a condition on the wave vector $\boldsymbol{k}$, which is equivalent to
$d_{2}|\boldsymbol{k}|^{2}-b_{2}-\mu=0$ and constrains the wave number to
$|\boldsymbol{k}|^{2}=\frac{b_{2}+\mu}{d_{2}}$.

In both of these cases the growth rate is defined by the thermal diffusion as
$\lambda=-\mu|\boldsymbol{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\neq 0$). We turn to non-trivial MGW
solutions (3.19) with $\omega\neq 0$. The simplest class are ‘vertically
varying’ MGWs with

$\displaystyle
k=0\,,\quad\alpha_{1}=-\alpha_{2}mf\,,\quad\beta_{1}=\beta_{2}=\gamma_{1}=\gamma_{2}=0\,,\quad\lambda=-\delta_{2}\,,\quad\omega=\pm
f\,,\quad\delta_{1}=\delta_{2}\,,$

so that $\mathtt{b}\equiv p\equiv 0$, while $\boldsymbol{v}\neq 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
$\boldsymbol{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\neq
0$ in the isotropic case, while for $d_{1}\neq d_{2}$ they are restricted to
$m^{2}=\frac{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 $\mu$ and Brunt-
Väisälä frequency $N^{2}$ have no impact, since these solutions trivially
satisfy the buoyancy equation by $v_{3}\equiv\mathtt{b}\equiv 0$.

Another type are ‘zonally varying’ MGWs with the coefficients satisfying

$\displaystyle\alpha_{1}=m=0\,,\quad\beta_{1}=\alpha_{2}\Bigl{(}N^{2}-\frac{(\delta_{3}-\delta_{\mu})^{2}}{4}\Bigr{)}k\,,\quad\beta_{2}=\alpha_{2}\frac{\delta_{3}-\delta_{\mu}}{2}k\,,$
$\displaystyle\gamma_{1}=\gamma_{2}=0\,,\quad\lambda=-\frac{\delta_{3}+\delta_{\mu}}{2}\,,\quad\omega^{2}=N^{2}-\frac{(\delta_{3}-\delta_{\mu})^{2}}{4}\,.$

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 $\lambda$
and phase frequency $\omega$ 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 $\boldsymbol{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}\rightarrow 0$ as $|\boldsymbol{k}|\rightarrow 0$.
Due to the equation for the growth rate $\lambda$, these MGWs are
exponentially decaying for kinetic energy backscatter, where
$\delta_{3}=-b_{3}|\boldsymbol{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}\leq\frac{b_{3}-\mu}{d_{3}}$, i.e. only in the case $b_{3}>\mu$.

The existence analysis for the remaining solutions to (3.19) with
$\omega,\,k,\,\alpha_{1}\neq 0$ is more involved. These solutions have the
coefficients

$\displaystyle\alpha_{1}=-\alpha_{2}mf\,,\quad\beta_{1}=\alpha_{1}\frac{m}{k}f+\alpha_{2}\frac{\omega^{2}}{k}|\boldsymbol{k}|^{2}\,,\quad\beta_{2}=\frac{\alpha_{2}}{k}(\delta_{3}k^{2}+\delta_{1}m^{2}-\delta_{2}|\boldsymbol{k}|^{2})\,,$
$\displaystyle\gamma_{2}=\alpha_{2}(\delta_{1}-\delta_{3})\frac{km}{|\boldsymbol{k}|^{2}}+\beta_{2}\frac{m}{|\boldsymbol{k}|^{2}}\,,\quad\gamma_{1}=-\alpha_{1}\frac{kf}{|\boldsymbol{k}|^{2}}-\beta_{1}\frac{m}{|\boldsymbol{k}|^{2}}\,,\quad\lambda=-\delta_{2}\,,$

with additional conditions, that also define the phase frequency $\omega$,
given by

$\displaystyle\omega^{2}=\Bigl{(}(\delta_{3}-\delta_{2})(\delta_{\mu}-\delta_{2})+N^{2}\Bigr{)}\frac{k^{2}}{|\boldsymbol{k}|^{2}}+\Bigl{(}(\delta_{1}-\delta_{2})(\delta_{\mu}-\delta_{2})+f^{2}\Bigr{)}\frac{m^{2}}{|\boldsymbol{k}|^{2}}>0\,,$
$\displaystyle\Bigl{(}(\delta_{\mu}-\delta_{2})N^{2}+(\delta_{3}-\delta_{2})\bigl{(}(\delta_{\mu}-\delta_{2})^{2}+\omega^{2}\bigr{)}\Bigr{)}k^{2}+(\delta_{1}-\delta_{2})\bigl{(}(\delta_{\mu}-\delta_{2})^{2}+\omega^{2}\bigr{)}m^{2}=0\,.$

The first condition 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 $|\boldsymbol{k}|$ sufficiently small, since
$\delta_{j},\,\delta_{\mu}\rightarrow 0$ for $|\boldsymbol{k}|\rightarrow 0$
for all $1\leq j\leq 3$. For unstable stratification $N^{2}<0$ one can first
choose $|\boldsymbol{k}|$ sufficiently small, so that the factor of $m^{2}$ is
positive. Then, for a fixed $|\boldsymbol{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 second condition here.
Both conditions together generate a set of solutions with rather complex
structure, as plotted in Figure 11 for the cases $\mu=0$ and $\mu>0$.
Due to the equation for the growth rate $\lambda=-\delta_{2}$, these MGWs are
exponentially and unboundedly growing for wave numbers
$|\boldsymbol{k}|^{2}<b_{2}/d_{2}$ and steady states for
$|\boldsymbol{k}|^{2}=b_{2}/d_{2}$ (see Figure 11). Notably, since $k\neq 0$,
these flows can transfer the horizontal backscatter (in fact only the
meridional component) to growing vertical velocity.

(a) $\mu=0.0$

(b) $\mu=0.12$

Figure 11:  Existence of MGW solutions on the wave vector plane $(k,m)$ in the
case $\omega,\alpha_{1}\neq 0$ (white curves) with growth rate $\lambda$
positive (red), zero (black) or negative (blue). Fixed parameter:
$d_{1}=0.12,\,d_{2}=0.11,\,d_{3}=0.0,\,b_{1}=0.14,\,b_{2}=0.12,\,b_{3}=-0.1,\,f=0.3,\,N^{2}=1$.

Superpositions and stability. It is possible to superpose MGWs and Kolmogorov
flows as a solution to (3.1) if these have the same direction of wave vectors
$\boldsymbol{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 §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 §3.2.2).
Without thermal diffusion ($\mu=0$) there are always growing Kolmogorov flows
in any direction with $k\neq 0$, for stable stratification $N^{2}>0$ at least
for all $|\boldsymbol{k}|^{2}<b_{2}/d_{2}$, and for unstable stratification
$N^{2}<0$ for $|\boldsymbol{k}|^{2}>b_{2}/d_{2}$. Thus, in these cases, the
steady MGWs with $k\neq 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\neq 0$ and $k=0$).

As in §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 $\boldsymbol{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 monochromatic inertia gravity waves.

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 unbounded instability in the
limit is expected to readily imply arbitrarily large growth that may however
be bounded 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.

#### Disclosure statement

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

#### 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.

## References

  * Achatz (2006) Achatz, U. (2006).  Gravity-Wave Breakdown in a Rotating Boussinesq Fluid: Linear and Nonlinear Dynamics.  Habilitation Thesis. University of Rostock.
  * Balmforth and Young (2005) Balmforth, N. J. and Young, Y.-N. (2005).  Stratified kolmogorov flow. ii.  J. Fluid Mech., 528:23–42.
  * Chai et al. (2020) Chai, J., Wu, T., and Fang, L. (2020).  Single-scale two-dimensional-three-component generalized-beltrami-flow solutions of incompressible navier-stokes equations.  Physics Letters A, 384(34):126857.
  * Danilov et al. (2019) Danilov, S., Juricke, S., Kutsenko, A., and Oliver, M. (2019).  Toward consistent subgrid momentum closures in ocean models.  In Eden, C. and Iske, A., editors, Energy Transfers in Atmosphere and Ocean, pages 145–192. Springer-Verlag, Cham.
  * Drazin and Riley (2006) Drazin, P. G. and Riley, N. (2006).  The Navier-Stokes equations: a classification of flows and exact solutions, volume 334 of London Mathematical Society Lecture Note Series.  Cambridge University Press, Cambridge.
  * Dyck and Straatman (2019) Dyck, N. J. and Straatman, A. G. (2019).  Exact solutions to the three-dimensional navier–stokes equations using the extended beltrami method.  Journal of Applied Mechanics, 87(1).  011004\.
  * Ghaemsaidi and Mathur (2019) Ghaemsaidi, S. J. and Mathur, M. (2019).  Three-dimensional small-scale instabilities of plane internal gravity waves.  Journal of Fluid Mechanics, 863:702–729.
  * Jansen et al. (2019) Jansen, M. F., Adcroft, A., Khani, S., and Kong, H. (2019).  Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies.  Journal of Advances in Modeling Earth Systems, 11(8):2844–2860.
  * Jansen and Held (2014) Jansen, M. F. and Held, I. M. (2014).  Parameterizing subgrid-scale eddy effects using energetically consistent backscatter.  Ocean Modelling, 80:36–48.
  * Juricke et al. (2020) Juricke, S., Danilov, S., Koldunov, N., Oliver, M., Sein, D., Sidorenko, D., and Wang, Q. (2020).  A kinematic kinetic energy backscatter parametrization: From implementation to global ocean simulations.  Journal of Advances in Modeling Earth Systems, 12(12):e2020MS002175.
  * Kalogirou et al. (2015) Kalogirou, A., Keaveny, E. E., and Papageorgiou, D. T. (2015).  An in-depth numerical study of the two-dimensional Kuramoto-Sivashinsky equation.  Proc. R. Soc. A., 471(2179):20140932, 20.
  * Lelong and Dunkerton (1998) Lelong, M.-P. and Dunkerton, T. J. (1998).  Inertia–gravity wave breaking in three dimensions. part i: Convectively stable waves.  Journal of the Atmospheric Sciences, 55(15):2473 – 2488.
  * Majda (2003) Majda, A. (2003).  Introduction to PDEs and waves for the atmosphere and ocean, volume 9 of Courant Lecture Notes in Mathematics.  New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI.
  * Majda and Wang (2006) Majda, A. J. and Wang, X. (2006).  Non-linear dynamics and statistical theories for basic geophysical flows.  Cambridge University Press, Cambridge.