1. Introduction
The concept of inverse energy cascades has played a central role in the development of turbulence theory, with its first applications in the context of two-dimensional and quasi-two-dimensional flows (Onsager Reference Onsager1949; Kraichnan Reference Kraichnan1967; Kraichnan & Montgomery Reference Kraichnan and Montgomery1980; Herring Reference Herring1988). First introduced as an application of statistical mechanics in fluids (Onsager Reference Onsager1949), it found a possible application to explain the emergence of large-scale structures in planetary atmospheres, and the absence of a spectral gap in the mesoscale range of the Earth’s atmospheric energy spectrum (Lilly Reference Lilly1983; Nastrom, Gage & Jasperson Reference Nastrom, Gage and Jasperson1984; Herring Reference Herring1988; Falkovich Reference Falkovich1992). However, early studies relied heavily on idealised models or on the interpretation of observational data within simplified models, primarily considering two-dimensional (2-D) or quasi-geostrophic (QG) flows (Charney Reference Charney1971). The latter approximation filters gravito-inertial waves and is valid in the hydrostatic limit of infinite rotation and infinite stratification (Julien et al. Reference Julien, Knobloch and Milliff2006; Vallis Reference Vallis2017). These models provide compelling theoretical frameworks, but only within strict asymptotic limits.
While such approximations can capture key aspects of geophysical flows at large scales, they remain crude simplifications. Real planetary atmospheres are strongly stratified, rotating and exhibit fast gravito-inertial (GI) waves and fully three-dimensional (3-D) motions (Riley & deBruynKops Reference Riley and deBruynKops2003; Waite Reference Waite2013; Dritschel & McKiver Reference Dritschel and McKiver2015). Under these conditions, studies showed that inverse energy cascades can develop, albeit their validity remains only in idealised limits, at very large planetary scales, or for parameters that are far from those expected in geophysical scenarios (Métais et al. Reference Métais, Bartello, Garnier, Riley and Lesieur1996; Smith & Waleffe Reference Smith and Waleffe2002; Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2013; Pouquet et al. Reference Pouquet, Marino, Mininni and Rosenberg2017; Oks et al. Reference Oks, Mininni, Marino and Pouquet2017). Nonetheless, the gap between idealised theory and the complexity observed in natural flows has decreased in recent years. It should be noted that planetary atmospheres, being highly dynamic and not necessarily homogeneous, affected by seasonality, and with density decreasing with altitude, exhibit significant variability in their governing parameters. For instance, estimates of the Reynolds number obtained through rocket and radar observations in the upper Earth’s atmosphere (Chau et al. Reference Chau, Urco, Avsarkisov, Vierinen, Latteck, Hall and Tsutsumi2020; Mesquita et al. Reference Mesquita, Larsen, Azeem, Stevens, Williams, Collins and Li2020) were found to range from
$10^3$
to
$10^4$
across the mesosphere and lower thermosphere. These values are consistent with those estimated in large numerical models (Liu, Hays & Roble Reference Liu, Hays and Roble1999; Lund & Fritts Reference Lund and Fritts2012) and compatible with the investigation we presented in a previous study (Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024). Numerical and theoretical studies have shown that when certain conditions are met even fully 3-D flows can exhibit an inverse cascade of energy (Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2010; Xia et al. Reference Xia, Byrne, Falkovich and Shats2011; Pouquet & Marino Reference Pouquet and Marino2013; Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014; Marino, Pouquet & Rosenberg Reference Marino, Pouquet and Rosenberg2015; Benavides & Alexakis Reference Benavides and Alexakis2017; Sahoo, Alexakis & Biferale Reference Sahoo, Alexakis and Biferale2017; Alexakis & Biferale Reference Alexakis and Biferale2018; van Kan & Alexakis Reference van Kan and Alexakis2020, Reference van Kan and Alexakis2022; Alexakis Reference Alexakis2023; Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024). These findings suggest that a bidirectional cascade can exist where larger scales cascade the energy inversely, while smaller scales cascade energy forward, as is observed also in the ocean (Balwada et al. Reference Balwada, Xie, Marino and Feraco2022). In this regime, the development of inverse cascades is not limited to idealised 2-D or QG models, but can also occur in physically realistic settings, provided that key geometrical or dynamical conditions are met, and that a source of turbulent energy exists at those scales (Gage & Nastrom Reference Gage and Nastrom1986).
Observational evidence from planetary atmospheres lends further motivation to consider these questions. In the Earth, observations indicate that energy cascades are mostly direct (Cho & Lindborg Reference Cho and Lindborg2001; Lindborg & Cho Reference Lindborg and Cho2001), albeit upscale energy transfer events have been reported (King, Vogelzang & Stoffelen Reference King, Vogelzang and Stoffelen2015). On Jupiter, high-resolution imaging and spectral analysis reveal the presence of coherent, long-lived structures fed by upscale energy transfer (Young & Read Reference Young and Read2017; Siegelman et al. Reference Siegelman2022). Similarly, Titan’s atmosphere – although much smaller and slower evolving – exhibits large-scale waves, complex wave–eddy interactions and planetary-scale structures in cloud cover (Mitchell & Lora Reference Mitchell and Lora2016). Inverse cascades have also been observed in laboratory experiments of rotating turbulent flows, further supporting the potential for energy to be transferred upscale in three dimensional systems, even at Rossby numbers as high as
$0.33$
(Gallet et al. Reference Gallet, Campagne, Cortet and Moisy2014; Campagne et al. Reference Campagne, Gallet, Moisy and Cortet2014), and in deep rotating flows where energy is injected at small scales (Yarom, Vardi & Sharon Reference Yarom, Vardi and Sharon2013; Yarom & Sharon Reference Yarom and Sharon2014). Inverse energy cascades in thick fluid layers were also observed in Xia et al. (Reference Xia, Byrne, Falkovich and Shats2011).
In this paper we examine the presence of inverse energy cascades in rotating stably stratified turbulent flows constrained to thin domains, in a range of parameters that are relevant for planetary atmospheres. In particular, we focus on regimes with thin aspect ratios, moderate Rossby numbers and small Froude numbers, that correspond to regimes close to those found in planetary atmospheres at scales that span from mesoscales down to the microscales, i.e. scales below geostrophic scales, in which turbulence and energy cascades, either rotation-dominated or stratification-dominated, become relevant. Our results show that, under certain conditions, inverse energy cascades can emerge, suggesting that this process can play a role in atmospheric self-organisation processes. We thus aim here at assessing within our simplified set-up under what conditions rotating stratified flows can develop inverse cascades. Our resolution remains limited, and definitive conclusions about the power-law exponents of the energy spectra cannot be drawn, especially at scales smaller than the forcing. However, one of the cases presented here was studied at larger spatial resolution, (see Alexakis et al. (Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024)). The choice to limit our resolution compared with Alexakis et al. (Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024) was made in order to cover a large parameter space, exploring how the inverse cascade is affected by rotation and stratification. Future work must also consider other effects not addressed here, such as the layer thickness, the effect of boundaries, the large-scale circulation that can feed a direct energy cascade from planetary scales, the modulation of stratification by the diurnal cycle and the influence of moisture.
2. Methodology
2.1. Set-up and numerical simulations
We consider a flattened triple periodic orthogonal box of dimensions
$2\pi L \times 2\pi L \times 2\pi H$
(where
$H=L/32$
corresponds to the
$z$
direction), filled with a linearly stratified fluid of mean density
$\overline {\rho }(z)=\rho _0+zS$
(with stable stratification
$S\lt 0$
) in presence of gravity
$\boldsymbol{-g}$
along the
$z$
direction, satisfying the Boussinesq equations
\begin{align} \partial _t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}+2\varOmega \times \boldsymbol{u} &= -\boldsymbol{\nabla }P - \boldsymbol{e}_z N \phi +\nu \boldsymbol{\nabla} ^2 \boldsymbol{u} + \boldsymbol{f} , \end{align}
\begin{align} \partial _t \phi + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi &= N \boldsymbol{e}_z\boldsymbol{\cdot }\boldsymbol{u} + \kappa \boldsymbol{\nabla} ^2 \phi , \end{align}
where
$\boldsymbol{u}$
is the velocity field,
$\phi$
the mass density fluctuations scaled to have units of velocity,
$\varOmega$
is the solid body rotation rate,
$N=(-gS)^{1/2}$
the Brunt–Väisälä frequency,
$P$
is the correction to hydrostatic pressure,
$\nu$
the viscosity,
$\kappa$
the density fluctuation diffusivity and
$\boldsymbol{f}$
an external forcing. The forcing is three-dimensional (varying along all three directions) but acting only on the horizontal components. It is random, with zero mean, and white in time, acting on wavenumbers
$\boldsymbol{k}$
satisfying
$ 0.9 \leqslant |\boldsymbol{k}|H \leqslant 1.1$
. This makes the forcing scale
$\ell _F$
be of the same order as the layer height
$H$
. The forcing is injecting energy on average at a rate
$\epsilon$
only on the horizontal components of the velocity field. The absence of a vertical component in the forcing implies that it is not exciting GI waves directly, although these can spontaneously appear as the flow evolves (Kafiabad & Bartello Reference Kafiabad and Bartello2016, Reference Kafiabad and Bartello2018; Thomas & Daniel Reference Thomas and Daniel2020). The Reynolds number for the system is defined as
$Re=\epsilon ^{1/3}k_{_H}^{-4/3}/\nu$
, and the Prandtl number as
$Pr=\nu /\kappa$
, and here is set to unity. The Rossby and Froude numbers for the system are defined respectively as
\begin{equation} \textit{Ro} =\epsilon ^{1/3} k_H^{2/3} /\varOmega \quad \mathrm{and} \quad \textit{Fr} =\epsilon ^{1/3} k_H^{2/3} /N \end{equation}
(with
$k_H=1/H$
and
$k_L=1/L$
). Note that all of these dimensionless numbers are defined using the vertical scale,
$H$
, of the computational domain, which is also the forcing scale.
These equations in the absence of forcing and dissipation conserve the total energy
\begin{equation} \mathcal{E} = \frac {1}{2}\bigl\langle |\boldsymbol{u}|^2 \bigr\rangle + \frac {1}{2}\bigl\langle \phi ^2 \bigr\rangle . \end{equation}
The equations are solved numerically using the GHOST code (Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011; Rosenberg et al. Reference Rosenberg, Mininni, Reddy and Pouquet2020), that is a pseudo-spectral code with a second-order Runge–Kutta method for the advancement in time and de-aliased using the 2/3 rule. A large set of 30 different high-resolution simulations was performed, each simulation on a grid of
$6144\times 6144\times 192$
, all with
$Re=500$
, and varying
$Ro$
in the range from
$Ro^{-1}=1/4$
to
$Ro^{-1}=4$
and
$Fr$
in the range
$Fr^{-1}=5$
to
$Fr^{-1}=160$
, both increasing in powers of 2. This range was chosen so that we observed a clear transition from inverse cascading to non-inverse cascading flows. Additionally, it specifically considers the regime characterised by
$Ro \approx 1$
and
$Fr \lt 1$
in planetary atmospheric scales below synoptic scales, where quasi-geostrophy does not hold. The rotation and stratification are categorised as strong or weak based on the specific parameter range examined.
Each simulation was run for several turnover times until small scales (smaller than
$H$
) have reached a steady state, and large scales (larger than
$H$
) have either reached a steady state as well or a sustained increase of their energy was observed with constant rate
$\gamma =\text{d}\mathcal{E}/\text{d}t$
. This constant increase is due to the presence of an inverse cascade where the nonlinearity continuously transfers a fraction of the injected energy to larger scales at which (unlike the energy transferred to small scales) it cannot be rapidly dissipated. This energy is then stored in the large scales and increases linearly with time
$\mathcal{E}=\gamma t$
. The coefficient
$\gamma$
is equal to the amplitude of the inverse cascade. We also confirmed that this energy increase was due to the presence of an inverse cascade by explicitly studying the energy spectrum and flux. This process would continue up until the largest scales of the domain are reached. However, we stop our simulations well before such a state is reached.
2.2. Energy spectra and fluxes
To quantify the scale-by-scale distribution of energy we consider the Fourier transformed fields
$\hat {\boldsymbol{u}}(\boldsymbol{k})$
and
$\hat {\phi }(\boldsymbol{k})$
. The kinetic and potential energy spectra are then given by
\begin{equation} E_K(k_i)= \frac {1}{2k_{_L}} \sum _{\boldsymbol{q}\in S_k} \vert \hat {\boldsymbol{u}}(\boldsymbol{q})\vert ^2, \qquad E_P(k_i)= \frac {1}{2k_{_L}} \sum _{\boldsymbol{q}\in S_k} \vert \hat {\phi }(\boldsymbol{q}) \vert ^2, \end{equation}
where
$S_k$
is the set of wavenumbers
$\boldsymbol{q}$
that satisfy
$k \leqslant \vert \boldsymbol{q} \vert \lt k+k_{_L}$
for spherically averaged (isotropic) spectra,
$k_\perp \leqslant \vert \boldsymbol{q}_\perp \vert \lt k_\perp +k_{_L}$
for cylindrically averaged (axisymmetric) spectra and
$k_\parallel \leqslant \vert {q}_\parallel \vert \lt k_\parallel +k_{_L}$
for plane-averaged (vertical) spectra (with
$k_{_L}=1/L$
). For 2-D spectra,
$S_k$
is the set that satisfies both
$k_\perp \leqslant \vert \boldsymbol{q}_\perp \vert \lt k_\perp +k_{_L}$
and
$k_\parallel \leqslant \vert {q}_\parallel \vert \lt k_\parallel +k_{_L}$
. The total energy spectrum is always given by
$E_T(k)=E_K(k)+E_P(k)$
, such that
$\mathcal{E} = \int E_T(k)\,\text{d}k$
.
We can also decompose these energies into the energy contained in GI waves and in QG modes, following Bartello (Reference Bartello1995) and Herbert, Pouquet & Marino (Reference Herbert, Pouquet and Marino2014). Given the vertical vorticity
$\hat {\omega }_\parallel (\boldsymbol{k})$
, the vertical velocity
$\hat {u}_\parallel (\boldsymbol{k})$
and the density variation
$\hat {\phi }(\boldsymbol{k})$
, a field given by
\begin{equation} \boldsymbol{Z}(\boldsymbol{k}) = [ \hat {\omega }_\parallel (\boldsymbol{k}),\, -ik\hat {u}_\parallel (\boldsymbol{k}),\, -k_\perp \hat {\phi }(\boldsymbol{k})] , \end{equation}
can be projected into GI and QG modes respectively as
\begin{equation} \mathcal{P}_{\textit{GI}}(\boldsymbol{Z})= \boldsymbol{Z}^+_{\textit{GI}} \frac {(\boldsymbol{Z}^+_{\textit{GI}})^*\boldsymbol{\cdot }\boldsymbol{Z}}{\vert \boldsymbol{Z}^+_{\textit{GI}}\vert ^2}+ \boldsymbol{Z}^-_{\textit{GI}} \frac {(\boldsymbol{Z}^-_{\textit{GI}})^*\boldsymbol{\cdot }\boldsymbol{Z}}{\vert \boldsymbol{Z}^-_{\textit{GI}}\vert ^2} , \end{equation}
and
\begin{equation} \mathcal{P}_{\textit{QG}}(\boldsymbol{Z})=\boldsymbol{Z}_{\textit{QG}} [\boldsymbol{Z}_{\textit{QG}}^*\boldsymbol{\cdot }\boldsymbol{Z}]/\vert \boldsymbol{Z}_{\textit{QG}}\vert ^2, \end{equation}
with
\begin{equation} \boldsymbol{Z}^\pm _{\textit{GI}}(\boldsymbol{k}) = [ i 2\varOmega k_\parallel ,\, -k\sigma ^\pm (\boldsymbol{k}),\, {\textit{Nk}}_\perp ], \quad \boldsymbol{Z}_{\textit{QG}}(\boldsymbol{k}) = [ i{\textit{Nk}}_\perp ,\, 0,\, 2\varOmega k_\parallel ] , \end{equation}
and where
$\sigma ^\pm (\boldsymbol{k})=\pm (4\varOmega ^2k_\parallel + N^2 k_\perp ^2)^{1/2}/k$
are the frequencies of the GI waves. The QG modes (also referred to as `slow modes' or `vortical modes' in the literature, in both cases corresponding to the same dynamical modes) have zero frequency. Energies and energy spectra can then also be computed using these projected fields.
The kinetic and potential energy fluxes
$\varPi _K(k)$
and
$\varPi _P(k)$
through a closed set of wavenumbers
$S_k$
are finally defined as
\begin{equation} \varPi _K(k) = \langle \boldsymbol{u}_{S_k} \boldsymbol{\cdot }(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}) \rangle , \qquad \varPi _P(k) = \langle {\phi }_{S_k} (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi ) \rangle , \end{equation}
where
$\boldsymbol{u}_{S_k}$
and
${\phi }_{S_k}$
are obtained by filtering the fields
$\boldsymbol{u}$
and
$\phi$
so that only wavenumbers in the set
$S_k$
are kept. The total energy flux is then given by
$\varPi _T=\varPi _K+\varPi _P$
.
In the simulations, energy spectra and fluxes are time averaged over a short time interval of the order of a few turnover times, in order to reduce naturally occurring fluctuations. In the case where we have an inverse energy cascade and where a steady state cannot be reached, the range of wavenumbers over which the inverse cascade is observed increases with integration time. We have practically chosen to terminate our runs when scales approximately ten times the forcing scale are excited.
3. Results
Table 1 presents the measured growth rate of the energy normalised by the energy injection rate,
$\gamma /\epsilon$
, in the last stages of the computation for all the simulations performed. The system detailed energy balance implies that
$0\leqslant \gamma /\epsilon \leqslant 1$
, with
$\gamma /\epsilon =0$
implying the absence of an inverse cascade, while
$\gamma /\epsilon =1$
implies that all of the injected energy is transferred towards the larger scales. Green boxes mark runs that did not display an inverse cascade, while red boxes indicate cases with a clear increase of energy, and thus the presence of an inverse cascade. Light pink marks simulations that displayed a weak increase of energy, and for which we cannot make a definite statement about the presence or not of an inverse cascade. We note that none of the runs displayed a strict inverse cascade
$\gamma /\epsilon =1$
such that all the energy cascades inversely. This implies that, when an inverse cascade is observed, it is part of a bidirectional or split cascade where part of the energy cascades to larger scales and part to smaller scales (Alexakis & Biferale Reference Alexakis and Biferale2018).
Measured mean inverse energy flux ratios as a function of
$Ro$
and
$Fr$
for all simulations. Green cells indicate no inverse cascade, while red cells indicate the presence of an inverse cascade, with a very weak inverse cascade marked by light pink.
Table 1 can be seen as a phase diagram, distinguishing regions of parameter space that develop or not an inverse cascade. Some clear trends are observed from this phase diagram. First, all runs with
$Ro^{-1}$
smaller than unity do not display an inverse cascade, while simulations with
$Ro^{-1}\geqslant 1$
can develop an inverse cascade depending on
$Fr$
with increasing inverse cascade as rotation increases. Note also that our definition of
$Ro$
measures the strength of rotation at the forcing scale, with
$Ro=1$
usually corresponding to weak rotation (Cambon, Rubinstein & Godeferd Reference Cambon, Rubinstein and Godeferd2004). Stratification on the other hand has the opposite effect. As stratification
$Fr^{-1}$
increases, the amplitude of the inverse cascade decreases. For very large stratification the inverse cascade is even suppressed. At both large stratification and large rotation rates the two effects compete making the boundary between inverse cascading runs and non-inverse cascading runs to be located at
$Ro/Fr \simeq 80$
. We note that the present results hold for a forcing scale
$\ell _F$
of the same order as the domain height
$H$
. The phase diagram in table 1 is expected to change if thinner domains are considered, with a larger range of parameters displaying an inverse cascade.
It is interesting that many cases with dimensionless numbers compatible with those of planetary atmospheres display inverse cascades. The case with
$Fr^{-1}=40$
and
$Ro^{-1}=1$
was examined in Alexakis et al. (Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024) at higher resolution and larger
$Re$
, showing that saturation of the inverse total energy flux occurs as turbulence strengthens, attaining levels compatible with those relevant for real atmospheres (Lilly Reference Lilly1983). In order to illustrate how controlling parameters determine features of the energy transfer and energy distribution in the rotating and stratified flows under study, in what follows we examine in more detail some of the extreme cases of the simulations performed.
3.1. Weak rotation and moderate stratification
Results for
$Ro=4$
and
$Fr=1/5$
: (a) 2-D energy spectrum as a function of
$k_\perp$
and
$k_\parallel$
. Note that the axes are shifted and scaled in such a manner that, in the log–log plot, the values of
$k_\parallel =0$
and
$k_\perp =0$
are visible. (b) Fraction of GI wave energy,
$R_{\textit{GI}}=E_{\textit{GI}}/E_T$
, as a function of
$k_\perp$
and
$k_\parallel$
, where
$E_T$
is the total energy. The black dotted line indicates the
$2\varOmega k_\parallel = {\textit{Nk}}_\perp$
balance where rotation and stratification are equally important. Above this line modes are rotation dominated, while below this line modes are stratification dominated. The white dashed lines indicate
$|\boldsymbol{k}|=$
constant spheres. Most of the energy injection occurs around
$k_\perp \approx k_\parallel \approx k_H$
. (c) Isotropic energy spectra, with several power laws indicated as references. The different energy components are described in the text. The red vertical dashed lines indicate the forcing wavenumbers. The inset shows the ratios of energy in GI and QG modes. (d) Isotropic energy fluxes. (e) Axisymmetric energy spectra, as a function of
$k_\perp$
. (f) Axisymmetric energy fluxes. (g) Plane energy spectra, as a function of
$k_\parallel$
. (h) Plane energy fluxes.
First, we focus in the case of weak rotation with
$Ro^{-1}=1/4$
, and moderate stratification with
$Fr^{-1}=5$
. From the examined runs, this simulation is the closest to the classical isotropic turbulence, and displays only a direct energy cascade. Figure 1 shows all energy spectra and fluxes for this case. Figure 1(a) shows the 2-D energy spectral density as a function of
$k_\parallel$
and
$k_\perp$
. The dashed white lines indicate isotropic contours (i.e. modes with constant wavenumber
$k$
), and the black dotted line marks the modes in which the inertial wave frequency matches the gravity wave frequency, i.e. modes with
$2\varOmega k_\parallel = {\textit{Nk}}_\perp$
. As we will see, in many regimes, energy tends to accumulate in the vicinity of these modes. Figure 1(b) shows the 2-D GI wave energy spectral density, normalised in such a way that wavenumbers with spectral density of 1 have energy only on wave modes, and wavenumbers with spectral density of 0 have all the energy in QG modes. Figures 1(c) and 1(d) display, respectively, isotropic energy spectra and fluxes. The total (T), kinetic (K), potential (P), GI wave and QG components are shown, with several power laws indicated as references. The vertical dashed lines indicate the range of forced wavenumbers. In the inset in figure 1(c) the ratios of the energy components
$R_{\textit{GI}} = E_{\textit{GI}}/E_T$
and
$R_{\textit{QG}} = E_{\textit{QG}}/E_T$
are shown, respectively, in pink and green. Figures 1(e) and 1(f) display the same spectra and fluxes as a function of the perpendicular wavenumber
$k_\perp$
. Finally, figures 1(g) and 1(h) show the same quantities as a function of the parallel wavenumber
$k_\parallel$
. The same quantities will be considered for all other cases in the next sections.
In this regime, the isotropic and perpendicular energy fluxes for
$k$
or
$k_\perp$
smaller than
$k_H$
, are zero or negligible for all energy components, indicating the absence of an inverse cascade, as also testified by
$\gamma =0$
observed in table 1. The cascade is thus strictly forward. In the smaller scales, where the forward cascade develops, energy is distributed almost isotropically (note, however, the decrease of energy for modes above the
$2\varOmega k_\| = {\textit{Nk}}_\perp$
curve indicated with dashes). The energy is dominated by gravity waves and develops isotropic and axisymmetric spectra not far from a
$k^{-5/3}$
power law. The plane-averaged energy spectra (panel g) show, however, a steeper slope, closer to
$k_\|^{-3}$
. Finally, scales larger than the forcing make a negligible contribution to the flow.
3.2. Strong rotation and moderate stratification
Results for
$Ro=1/4$
and
$Fr=1/5$
. References for all panels are as in figure 1.
Figure 2 shows the same quantities for the case with
$Ro^{-1}=4$
and
$Fr^{-1}=5$
. All references are the same as in figure 1. A strong inverse cascade of kinetic energy develops, as is evident from the isotropic and perpendicular fluxes that are negative and close to
$-1$
, with a weak direct energy cascade. It is worth noting that a very small fraction of the energy cascades in the parallel direction. The energy spectrum peaks at wavenumbers with
$k_\perp$
smaller than
$k_H$
while
$k_\|=0$
. It is dominated by the kinetic energy and has a strong contribution of QG modes at those scales. This case is thus very close to the 2-D inverse energy cascade. The small scales display a steep
$\sim k^{-3}$
scaling with some possible recovery of isotropy at very small scales. The overall behaviour is similar to that reported in previous studies of rotating and stratified turbulence with an inverse cascade when
$Ro/Fr$
is
$\mathcal{O}(1)$
(Métais et al. Reference Métais, Bartello, Garnier, Riley and Lesieur1996; Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2013). The emergence of an inverse cascade process is very similar to that observed in purely, strongly rotating flows in which energy goes towards 2-D modes through near-resonant or four-wave interactions (Smith & Waleffe Reference Smith and Waleffe1999; Bourouiba & Bartello Reference Bourouiba and Bartello2007; Deusebio et al. Reference Deusebio, Boffetta, Lindborg and Musacchio2014; Clark Di Leoni & Mininni Reference Clark di Leoni and Mininni2016; van Kan & Alexakis Reference van Kan and Alexakis2020). A thin domain helps this transfer, as less resonant interactions become available; once energy is in the 2-D modes, the decoupling of these modes from 3-D modes caused by rotation allows most of the energy to go towards smaller
$k$
(Clark Di Leoni et al. Reference Clark Di Leoni, Alexakis, Biferale and Buzzicotti2020).
3.3. Weak rotation and strong stratification
Results for
$Ro=4$
and
$Fr=1/160$
. References for all panels are as in figure 1.
Figure 3 shows the same quantities for the case with
$Ro^{-1}=1/4$
and
$Fr^{-1}=160$
. While there is a small peak in the energy for
$k_\perp \lt k_H$
in the axisymmetric kinetic energy spectrum, and the axisymmetric energy flux becomes negative in figure 3(f), the vertical energy flux in figure 3(h) is positive and close to unity, making the isotropic flux zero and the isotropic spectrum display small amplitudes for
$k_\perp \lt k_H$
. This phenomenon has been reported before in purely stratified or weakly rotating and stratified turbulence simulations (Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2014; Sozza et al. Reference Sozza, Boffetta, Muratore-Ginanneschi and Musacchio2015). It corresponds to the development of same density layers and strongly anisotropic horizontal winds with a large-scale correlation in the horizontal direction but small-scale correlation in the vertical (Smith & Waleffe Reference Smith and Waleffe2002; Caulfield Reference Caulfield2021). This process does not correspond to an inverse energy cascade. Energy goes towards smaller
$k_\perp$
but towards larger
$k_\parallel$
as these winds develop, resulting in the formation of pancake-like structures, and in a net transfer of energy towards larger
$k$
even though the energy transfer in
$k_\perp$
may seem inverse.
3.4. Strong rotation and strong stratification
Results for
$Ro=1/4$
and
$Fr=1/160$
. References for all panels are as in figure 1.
Figure 4 displays the case with
$Ro^{-1} = 4$
and
$Fr^{-1}=160$
. This case is the closest to the QG limit where both stratification and rotation are considered large (Charney Reference Charney1971; Vallgren & Lindborg Reference Vallgren2010; van Kan & Alexakis Reference van Kan and Alexakis2022), albeit our values are far from being those needed to be in the QG asymptotic limit. The first difference, compared with the previous cases, is the strong dominance in figure 4(b) of GI waves in the modes with
$2 \varOmega k_\parallel = N k_\perp$
, while all other modes are mostly QG. This feature is present in all simulations with strong rotation and stratification. This, together with the layer thickness
$H/\ell _{_F}\approx 1$
, is enough to drive an inverse energy cascade even under very strong stratification provided
$Ro/Fr = N/\varOmega \lesssim 80$
.
The spectra and fluxes in figure 4 confirm the development of an inverse energy cascade of the kinetic energy. A
$\sim k^{-5/3}$
scaling range can be seen at the isotropic and axisymmetric energy spectra. Both the isotropic and the axisymmetric kinetic energy fluxes display a range of wavenumbers with negative flux for
$k$
and
$k_\perp$
smaller than
$k_H$
. Note, however, that the inverse isotropic flux is significantly smaller than the axisymmetric flux. It is the isotropic flux the one that confirms the development of a net inverse energy cascade, while the axisymmetric flux significantly overestimates the inverse transfer (Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024).
The energy in this regime, which is the more relevant for planetary atmospheres, is self-organised at large scales in a strongly anisotropic process. Rotation, stratification and the constraint of the thickness of the domain all play a role, feeding QG modes in a regime far from the asymptotic limit assumed in that model. Despite the dominant stratification that aids the transfer towards larger
$k_z$
an inverse cascade is present. The panels in figure 4 demonstrate how it can develop in a situation what would be expected to favour the forward cascade. At the forcing scale, stratification constrains energy to QG modes and to the formation of thin layers in the flow. These structures move energy to smaller
$k_\perp$
and larger
$k_\parallel$
, just as in the case in § 3.3. However, as the pure inertial wave frequency increases with
$k_\parallel$
, this process ends at wavenumbers at which both rotation and stratification become comparable, i.e. with
$2\varOmega k_\parallel \approx {\textit{Nk}}_\perp$
. Rotation then becomes dominant and tends to bidimensionalise the flow as in the case in § 3.2, preventing energy from going to modes with larger
$k_\parallel$
. Instead, a fraction of the energy cascades forward as GI waves, while the rest goes to modes with
$k_\parallel =0$
, which are not affected by rotation, and in which energy can inverse cascade towards smaller values of
$k$
as in 2-D turbulence. The stability of these modes is ensured by rotation and by the thin domain geometry (Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024).
The scaling observed at small scales (i.e. wavenumbers larger than
$k_H$
) in this regime is strongly dependent on
$Ro$
,
$Fr$
and
$Re$
, with different power laws developing in the direct cascade range depending on the strength of rotation and of stratification.
4. Conclusions
We studied the conditions under which inverse energy cascades develop in rotating and stably stratified turbulence, using direct numerical simulations in thin domains with parameters of interest for planetary atmospheres. By varying the Rossby and Froude numbers with fixed domain aspect ratio and forcing scale, we constructed a phase diagram of inverse cascade cases. Our results indicate that rotation helps the emergence of the inverse cascade while stratification tends to suppress it. In domains of the size of the forcing scale, such as the ones considered here, rotation such that
$Ro \leqslant 1$
is enough for the system to transfer energy from the forcing to larger scales, leading to the development of coherent structures, for cases with moderate to strong stratification
$Fr^{-1}\lesssim 80$
. In the large rotation (
$Ro^{-1}\gtrsim 1$
) and large stratification (
$Fr^{-1}\gtrsim 80$
) limit the inverse cascade appears to be present when
$Ro/Fr\lesssim 80$
, although further studies at even smaller values of
$Ro$
and
$Fr$
are required to verify this asymptotic behaviour. We note that this critical ratio
$Ro/Fr\approx 80$
is also expected to depend on
$\ell _f/H$
, that here was kept fixed to unity. Larger values of
$\ell _f/H$
that favour the inverse cascade are expected to increase the critical ratio
$Ro/Fr$
, while smaller values of
$\ell _f/H$
should decrease it. We also note that in the simultaneous fast rotating limit
$Ro\to 0$
and
$\ell _f/H\to 0$
, that was examined in van Kan & Alexakis (Reference van Kan and Alexakis2022), the phase diagram was different than what is shown in table 1 here. Thus,
$\ell _f/H$
is another important parameter that remains to be examined.
The resulting flows, for both inverse cascading regimes and strictly forward cascading regimes, were analysed in terms of energy fluxes and 2-D spectra, showing a much more complex distribution of energy in the
$(k_\perp ,k_\parallel )$
plane than what has been typically modelled as power laws. Furthermore, the role of the gravity-inertial modes in driving the forward cascade has been emphasised. Finally, we would like to mention that intermittency, that is known to be important in stratified flows (Rorai, Mininni & Pouquet Reference Rorai, Mininni and Pouquet2014; Feraco et al. Reference Feraco, Marino, Pumir, Primavera, Mininni, Pouquet and Rosenberg2018; Petropoulos et al. Reference Petropoulos, Couchman, Mashayek, de Bruyn Kops and Caulfield2024) could play an important role in the development of these spectra.
We emphasise that the observed inverse cascades are not the result of idealised 2-D or QG approximations. Instead, they arise in 3-D flows constrained by rotation, stratification and geometry, a set of conditions relevant in geophysical and planetary contexts. Notably, several simulations with Rossby and Froude numbers comparable to those found in Earth’s and Jupiter’s atmospheres show inverse cascades. These results provide theoretical support to the idea that such mechanisms can in principle partially contribute to the large-scale organisation of flows in nature.
However, we emphasise that the presence of an inverse cascade in the simulations does not prove that the same process develops in nature. Their occurrence in practice, whether sustained or transient in time, requires a turbulent injection mechanism at small scales. Such mechanisms, that include convection-induced turbulence in the upper layers of planetary atmospheres or gravity wave breaking, have been long hypothesised (Lilly Reference Lilly1983; Nastrom et al. Reference Nastrom, Gage and Jasperson1984). But while they inspire our forcing, these are not realistically captured in our simulations. Real atmospheres are also subject to additional effects not captured in our simplified 3-D rotating-Boussinesq system, such as the effect of boundary conditions and topography, moisture, radiative forcing, the modulation in stratification resulting from the diurnal cycle and the presence of a large-scale circulation at planetary scales with its associated forward energy cascade. Our aim is instead to assess whether inverse cascades can occur in physically plausible regimes of this simplified, yet non-idealised, 3-D systems. This study thus provides a stepping stone between asymptotic models that allow analytical demonstration of inverse cascades and the far more complex physics of real atmospheres. The results suggest that this concept, originally developed in these idealised asymptotic models, can remain relevant for interpreting self-organisation in more realistic natural settings.
Acknowledgements
The authors acknowledge HPC facilities at the École Normale Superieure in Paris (France), at École Centrale de Lyon (PMCS2I) in Écully (France), and at Departamento de Fisica (FCEN, UBA) in Argentina for the support in analysing the data.
Funding
Computer resources in Joliot-Curie at CEA were provided by PRACE (research project ID 2020235566) and by GENCI-TGCC & GENCI-CINES (Project No. A0170506421). This work was supported by the projects ‘DYSTURB’ (no. ANR-17-CE30-0004) and ‘EVENTFUL’ (no. ANR-20-CE30-0011) funded by the French Agence Nationale de la Recherche (ANR), by the project ‘STREAM’ (no. APR-2025) funded by the French Centre National d’Études Spatiales (CNES), and by Proyecto REMATE of Redes Federales de Alto Impacto, Argentina.
Declaration of interests
The authors report no conflicts of interest.
Data availability statement
The GHOST code used for all the simulations is openly available in https://doi.org/10.5281/zenodo.15472741.
Open access
