# Large-scale self-organization in dry turbulent atmospheres

Alexandros Ens Alexakis, Raffaele Marino, Pablo D Mininni, Adrian van Kan, Raffaello Foldes, Fabio Feraco

► To cite this version:

Alexandros Ens Alexakis, Raffaele Marino, Pablo D Mininni, Adrian van Kan, Raffaello Foldes, et al. Large-scale self-organization in dry turbulent atmospheres. Science, 2024, 383 (6686), pp.1005-1009. 10.1126/science.adg8269. hal-04779207

HAL Id: hal-04779207

https://hal.science/hal-04779207v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Large-scale self-organisation in dry turbulent
atmospheres

Alexandros Alexakis¹, Raffaele Marino², Pablo D. Mininni³,
Adrian van Kan⁴, Raffaello Foldes², Fabio Feraco²,⁵

1* LPENS, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris-Diderot

²CNRS, École Centrale de Lyon, INSA de Lyon, Université Claude Bernard Lyon 1,
LMFA, UMR5509, F-69134 Écully, France

³Universidad de Buenos Aires, Departamento de Física,
(INFINA), Ciudad Universitaria, 1428 Buenos Aires, Argentina

⁴Department of Physics, University of California, Berkeley,
Berkeley, 94720, California, USA

⁵Leibniz-Institute of Atmospheric Physics, University of Rostock,
Schloßstraße 6, Kühlungsborn, 18225, Germany

* E-mail: alexakis@phys.ens.fr

This work has been published in:

Science 383, 1005–1009 (2024). DOI:10.1126/science.adg8269

How turbulent convective fluctuations organise to form large-scale
structures in planetary atmospheres remains a question that eludes
quantitative answers. The assumption that this process is the re-
sult of an inverse cascade was suggested half a century ago in
two-dimensional fluids, but its applicability to atmospheric and
oceanic flows remains heavily debated, hampering our understand-
ing of the energy balance in planetary systems. We show with
direct numerical simulations of spatial resolutions of 12288² ×
384 points that rotating and stratified flows can support a bi-
directional cascade of energy, in three dimensions, with a ratio
of Rossby to Froude numbers comparable to that of the Earth's
atmosphere. Our results establish that in dry atmospheres spon-
taneous order can arise via an inverse cascade to the largest spa-
tial scales.

# 1 Introduction

Flow structures thousands of kilometers wide are not uncommon in the atmo-
sphere of the Earth and of other planets. The energy of these structures could
originate from processes associated with the global atmospheric circulation, but
also from smaller scale convective turbulence. In the latter case, small-scale ed-
dies conspire to self-organise into larger structures. Such a process goes against
our daily life experience where turbulence generates smaller-scale erratic struc-
tures as one observes when pouring milk in a coffee cup. One thus needs to come
up with convincing mechanisms for how such large-scale organization can take
place in planetary atmospheres.

One of the most important theoretical discoveries in the twentieth century
in the field of non-equilibrium physics is the phenomenon of self-organization,
which spontaneously creates large-scale order out of small-scale disorder. One of the first examples of this process was given by Onsager with the statistical mechanics of a gas of point vortices [1] that was later generalized to two-dimensional turbulent flows [2, 3, 4, 5]. A two-dimensional flow conserves an additional invariant, the enstrophy, given by the mean squared vorticity. The relation between energy and enstrophy leads to an incompatibility for the simultaneous bulk transfer of both quantities to the small scales. As a result vortices self-interact transferring enstrophy to smaller scales while energy is transferred to larger scales. This process takes place on a continuum of scales, forming a constant flux of energy from small to large scales in what is known as an inverse energy cascade, as opposed to the disordered forward energy cascade observed in three-dimensional turbulence that is directed to small scales.

Although planetary atmospheres are often very thin (the Earth's atmosphere has horizontal synoptic scales of the order of 1000 km, and a pressure scale height of 15 km), the corresponding flows are far from being two-dimensional. Nonetheless, two-dimensionality is not imperative for the appearance of self-organisation. Three-dimensional rotating and stratified flows (two key ingredients of atmospheric dynamics) conserve a different invariant, the potential vorticity, that can also lead to an inverse cascade. This happens in the quasi-geostrophic limit where rotation and stratification are asymptotically strong [6]. In this regime the flow is restricted to modes that satisfy hydrostatic balance, and as a result gravito-inertial waves are filtered out. Inverse cascades can also be present in rotating Rayleigh-Bénard convection [7, 8, 9, 10, 11], where in this case a generalized
quasi-geostrophy limit can be considered that partially preserves gravito-inertial
modes. However for most planetary flows, the quasi-geostrophic limit is at best
a crude approximation, with gravito-inertial waves composing a significant part
of the energy budget cascading energy forward [12, 13, 14, 15]. Thus, an inverse
cascade in planetary atmospheres caused either by two-dimensionality or quasi-
geostrophy remains conjectural.

Could atmospheric dynamics display an inverse cascade away from these lim-
its? In recent years it has been demonstrated that a hybrid state can be reached
such that large scales cascade energy inversely while smaller scales cascade en-
ergy forward in what is now known as a bidirectional cascade [16]. Bidirectional
cascades have been observed early on with direct numerical simulations (DNS)
[17, 18, 19, 20, 21, 22]. In rotating and stratified flows, simulations also indicate
the presence of bidirectional cascades [23, 24, 25], though in a regime where ro-
tation and stratification are comparable in strength, which is typical for the ocean
but not for the atmosphere.

Nonetheless, the existence of self-organisation processes through a bidirec-
tional cascade in planetary atmospheres becomes a compelling possibility as re-
cent research using satellite images with cloud tracking analysis, and in-situ air-
craft measurements, estimated the flux (and thus also the direction) of the energy
cascades in planetary flows in the Earth's atmosphere [26, 27], the ocean [28],
and the Jovian atmosphere [29, 11]. These studies affirmed the presence of both
inverse and forward energy cascades depending on the scale examined or on the al-
titude. However, satellite images constrain the measurements to two-dimensional
slices ignoring thus any processes occurring in the third direction.

Up to today there is no definite evidence of whether planetary atmospheric flows satisfy the necessary conditions for a bidirectional cascade to establish itself. The difficulty in answering such questions lies on the one hand in the fact that information from satellite images is limited, and on the other hand in the extreme parameter values that are met in planetary atmospheres which are hard to obtain in DNS. However, not only the technology to perform high-cadence high-resolution observations of the atmosphere has just started to come along, but also the computational power to perform DNS of stratified atmospheres in a realistic parameter space has become available. Here, with the use of DNS in a large grid using 40,000,000 CPU hours, we establish that the fluid model of a rotating and stably stratified dry atmosphere described by the non-hydrostatic Boussinesq equations can generate a bidirectional cascade leading to large-scale organisation of the flow.

## 2 Results

We consider a fluid in a Cartesian, triply periodic domain of vertical height $H$ and horizontal dimension $L = 32H$, in the presence of gravity, a stable mean density gradient and solid body rotation in the vertical direction [30]. The dynamics of the system are described by the incompressible velocity field $\mathbf{u}$ and the normalized
density variation $\phi$, governed by the Boussinesq equations [31, 32]:

$$ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\Omega \times \mathbf{u} = -\nabla P - \mathbf{e}_z N \phi + \nu \nabla^2 \mathbf{u} + \mathbf{f}, \quad (1) $$

$$ \partial_t \phi + \mathbf{u} \cdot \nabla \phi = N \mathbf{e}_z \cdot \mathbf{u} + \kappa \nabla^2 \phi, \quad (2) $$

where $\Omega$ is the solid body rotation rate, $N$ the Brunt-Väisälä frequency, $P$ is the pressure, $\nu$ the viscosity, $\kappa$ the density diffusivity, and $\mathbf{f}$ an external forcing acting at scales $\ell_F \sim H$ injecting energy at a rate $\epsilon$. Although this model has some strong simplifications, like periodicity or a simplified forcing mechanism, it is the most elementary model capturing the necessary physics to reproduce atmospheric dynamics.

This system has four independent non-dimensional control parameters: (i) the Reynolds number $Re_\epsilon = \epsilon^{1/3}k_H^{-4/3}/\nu$, (ii) the Prandtl number $Pr = \nu/\kappa$, that here is set to unity, (iii) the Rossby number, $Ro_\epsilon = \epsilon^{1/3}k_H^{2/3}/\Omega$, and (iv) the Froude number $Fr_\epsilon = \epsilon^{1/3}k_H^{2/3}/N$ (with $k_H = 2\pi/H$). We can also define dimensionless parameters based on the domain size $L$ and the flow r.m.s. velocity $U$, as, e.g., $Re = UL/\nu$, $Ro = U/(H\Omega)$, and $Fr = U/(NH)$, which are closer to the definitions used in atmospheric measurements.

Simulations were performed at resolutions of 12288² × 384 grid points [33, 30]. As a reference, for $H = 15$ km this corresponds to a domain length of 480 km, and a vertical and horizontal resolution of 39 m. Our simulations are characterised by $Re_\epsilon = 2000$, $Ro_\epsilon = 1$, and $Fr_\epsilon = 0.025$, or alternatively $Re \approx 2 \times 10^6$, $Ro \approx 0.4$, and $Fr \approx 0.01$. These values are compatible, e.g., with that of the Mesosphere-Lower Thermosphere (MLT) [34].
<figure><img src="image_2.png" /><figcaption>Figure 1: **Visualisation of density fluctuations $\phi$ and of the velocity field in the computational domain.** Structures at scales much larger than the forcing (i.e., at the scale of the domain height, with wavenumber $k_f = 2\pi/H$) are abundant in the visible horizontal plane (left), indicative of an efficient transfer of the energy towards the lowest modes along the perpendicular direction in Fourier space. The large-scale patterns are visible in the flow visualisation shown by arrows in a zoom in the top right. At the same time, three-dimensional instability patterns and small-scale features are detectable in both horizontal and vertical cuts of the zoomed simulation domain (bottom right), suggesting the action of a forward turbulent cascade. See the Supplementary Information for a movie of the density fluctuations in the entire domain. Visualisations were done with VAPOR [35].</figcaption></figure>

Fig. 1 shows visualisations of the flow and the density field. Large-scale structures of horizontal width ten times larger than $H$ can already be seen by visual inspection. At the same time, looking at the zoomed-in cross-sections one realizes that these structures are far from being two-dimensional. In the large scales, pancake structures of alternating sign of $\phi$ along the vertical, and the emergence of macroscopic cyclones and anti-cyclones, are visible in Fig. 1 (top right panel). These features are observed in weather maps and are a landmark of the large-
scale energy-containing structures in the Earth’s atmosphere. At the same time
small-scale overturning events can be seen in the zoomed in section that are ten
times smaller than *H* (see the bottom right panel in Fig. 1), and are also detectable
sometimes in the sky as Kelvin-Helmholtz billows. Thus, already at this qualita-
tive level one can testify the presence of a bidirectional cascade.

To become more quantitative we note that the inviscid Boussinesq equations
conserve the total energy $\mathcal{E}_T$ given by the sum of kinetic energy $\mathcal{E}_K$ and poten-
tial energy $\mathcal{E}_P$. Alternatively we decomposed it into the energy of gravito-inertial
modes $\mathcal{E}_{GW}$ and energy of quasi-geostrophic modes $\mathcal{E}_{QG}$ where $\mathcal{E}_T = \mathcal{E}_{GW} + \mathcal{E}_{QG}$.
Gravito-inertial modes are dispersive waves-modes due to the combined restor-
ing force of gravity and Coriolis, while quasi-geostrophic modes are modes that
satisfy hydrostatic balance (see Methods for their exact definitions). These ener-
gies are distributed differently in the Fourier space, among the wavenumbers $k_\|$,
identifying horizontal planes, and $k_\perp$, identifying vertically axisymmetric cylin-
ders. We define three different energy spectra averaged over fixed $k_\|$, $k_\perp$, and
$k = \sqrt{k_\|^2 + k_\perp^2}$. We do not define a new symbol for each spectrum, but distin-
guish between them by their argument, i.e., $E_i(k)$, $E_i(k_\perp)$, and $E_i(k_\|)$ where $i$
is $T$, $K$, $P$, $GW$, or $QG$, which stand for total, kinetic, potential, gravito-inertia
wave, and quasi-geostrophic respectively. In addition we define $E_i(k_\perp, k_\|)$ that
shows the spectral energy density for a given pair $k_\perp, k_\|$.

The left panels in Fig. 2 show the energy spectra $E_i(k)$, $E_i(k_\perp)$ and $E_i(k_\|)$
with the energy component $i$ as indicated in the legend. The inset also shows the
ratios $R_{GW} = E_{GW}/E_T$ and $R_{QG} = E_{QG}/E_T$. In the top panel, for $k > k_H$, the
<figure><img src="image_3.png" /></figure>

Figure 2: **Energy spectra and fluxes.** Left panels: Spherically (top), Cylindrically (middle), and plane (bottom) averaged energy spectra, for all energy components. Insets show the ratios of energy components $R_{GW} = E_{GW}/E_T$ (pink) and $R_{QG} = E_{QG}/E_T$ (green). Right panels: Energy fluxes across spheres (top), cylinders (middle), planes (bottom) in spectral space. Total energy flux (black line), energy flux of kinetic energy (blue line), energy flux of potential energy (red line). The forcing wavenumber $k_{\perp} \approx k_H$ (where energy is injected), the buoyancy wavenumber $k_B = N/U$, and the dissipation wavenumber $k_{\nu}$ (where energy is dissipated) are indicated by vertical dashed lines. The inset in the top right panel shows the inverse flux as a function of Re obtained from different runs.

spectra have been averaged over shells of width $k_H$ since otherwise large peaks of period $k_H$ are observed due to the strong domain anisotropy, shown by the light grey lines for the total energy spectrum. For wavenumbers larger than $k_H$
and smaller than the viscous wavenumber $k_ν$ the spherically averaged spectrum displays a power-law behavior with an exponent very close to Kolmogorov's prediction $k^{-5/3}$ for homogeneous isotropic turbulence. This power-law behavior, composed 70% by gravito-inertial waves, is indicative of a forward energy cascade. At $k$ smaller than $k_H$ a similar power-law is observed (albeit with a smaller prefactor). This indicates the presence of an inverse cascade. This energy at small $k$ is almost exclusively kinetic, dominated by 2D quasi-geostrophic modes.

For $E_i(k_⊥)$ three different power-laws can be observed. First, in the range $k_L < k_⊥ < k_H$ a $k_⊥^{-5/3}$ law is observed, where $k_L = 2\pi/L$. Second, in the range $k_H < k_⊥ < k_B = N/U$ a steeper power-law close to $k_⊥^{-3}$ is observed, where $k_B$ is the buoyancy wavenumber. Finally, at larger $k_⊥$, a shallower power-law slope starts to appear with exponent close to $-5/3$. The spectra for $k_⊥ > k_H$ resemble the so-called Nastrom-Gage spectra [36] observed in the Earth's troposphere, although we should note here that, in the atmospheric case, the transition from the $k^{-3}$ regime to the $k^{-5/3}$ occurs at a scale $k' < k_H$ while here it occurs for $k > k_H$. Finally, the last panel of Fig. 2 shows $E_i(k_‖)$, which displays a $k^{-5/3}$ power-law behavior.

The right panels of Fig. 2 show the energy fluxes across different surfaces in wavenumber space: across constant $k$ spheres $\Pi_i(k)$, constant $k_⊥$ cylinders $\Pi_i(k_⊥)$ and constant $k_‖$ planes $\Pi_i(k_‖)$. As with the spectra, we distinguish between fluxes based on their arguments. Here $i$ is $T$, $K$ or $P$ that stand for total, kinetic, and potential energies respectively. Positive values imply a flux of energy towards large wavenumbers, while negative values indicate a flux towards
small wavenumbers. $\Pi_T(k)$ flux is positive for $k > k_H$ indicating a forward cascade. However, a small fraction, corresponding to 5% of total energy injection rate, cascades towards the large scales. This is seen in the negative flux observed at $k < k_H$. This flux is also constant for more than a decade of wavenumbers almost up to $k_L$. The inset in Fig. 2 shows the amplitude of this negative flux measured from different simulations varying only $Re$. The flux increases with $Re$ and saturates at the largest $Re$. $\Pi_T(k_\perp)$ is also positive for $k_\perp > k_H$ and negative for $k_\perp < k_H$. However, in this case the fraction of energy that cascades towards smaller $k_\perp$ is five times larger than $\Pi_T(k)$. $\Pi_T(k_\|)$ is positive everywhere.

<figure><img src="image_4.png" /><figcaption>Figure 3: **Two-dimensional energy spectra in log-log scale.** White dashed lines indicate isotropic contours (i.e., modes with constant wavenumber $k$). The solid white line indicates the maximum resolved wavenumbers. Left: Total energy spectrum. Right: Gravity wave energy spectrum ratio. The black dotted lines mark $2\Omega k_\| = Nk_\perp$ where inertial wave frequency matches gravity wave frequency. The arrows indicate the direction of the flux of energy.</figcaption></figure>

While in one-dimensional spectra and fluxes it is easier to identify power laws, the energy distribution depends on $k_\|$ and $k_\perp$ independently. In Fig. 3 we show colour shaded plots of $E_T(k_\perp, k_\|$) and $R_{GW} = E_{GW}(k_\perp, k_\|)/E_T(k_\perp, k_\|)$. The
arrows indicate the direction of the energy transfers based on the fluxes in Fig. 2. A part of the injected energy is transferred to larger $k_{\perp}$, producing the $k_{\perp}^{-3}$ spectrum observed in the $k_H < k_{\perp} < k_B$ range. $R_{GW}(k_{\perp}, k_{\|})$ indicates that this forward transfer takes place through GW modes. The peak of $E(k_{\perp}, k_{\|})$ is observed at $k_{\perp} \simeq 2k_L$ and $k_{\|} \simeq 2k_H$ formed by an inverse transfer indicated by the black arrow and is dominated by QG modes. This energy is responsible for the formation of the $k_{\perp}^{-5/3}$ scaling at small $k_{\perp}$ in Fig. 2. When rotation dominated scales are reached at $2\Omega k_{\|} \simeq Nk_{\perp}$ (black dotted line), QG modes transfer their energy to GW modes that cascade it back to small scales. Of the energy that has moved to smaller $k_{\perp}$, a finite amount is transferred (cyan arrow) below the smallest dashed white line ($k_{\perp}^2 + k_{\|}^2 = k_H^2$). The energy in these modes forms the $E(k) \propto k^{-5/3}$ spectrum for $k < k_H$. This energy is the only that escapes to the largest scales $k \to 0$ and corresponds to a true inverse cascade.

### 3 Conclusions

This work has shown that dry turbulent atmospheres modeled by the non-hydrostatic Boussinesq equations can lead to a bidirectional energy cascade. The results showed that indeed there is a flux of energy directed to the small wavenumbers $k$, corresponding to 5% of the total energy injection rate at the largest Reynolds number. This flux, albeit small, is shown to persist up to the largest scales of the system and is $Re$-independent as large values of $Re$ are reached.

Our analysis provides a detailed description of how energy is transferred across scales and between different modes. These transfers, indicated by the arrows
sketched in Fig. 3, summarise the results in this work. Stratification, rotation and the geometric constraint of finite *H* all play a role in the formation of this inverse cascade. In physical terms, at the scale of the forcing stratification is dominant, constraining a large fraction of the energy to hydrostatic (QG) modes. This leads to the formation of pancake structures, known in stratified turbulence [37], that move energy to smaller *k*⊥ and larger *k*∥. This process ceases at wavenumbers where stratification is comparable to rotation *Nk*⊥ ∝ 2Ω*k*∥. Rotation, that tends to bi-dimentionalise the flow [38], prevents larger *k*∥ modes from appearing and energy is converted to gravito-inertial mode energy that cascades back to larger *k*. This is true for all *k* except for the *k*∥ = 0 modes that are unaffected by rotation. These modes that follow two-dimensional dynamics cascade the energy to ever smaller *k*⊥. Their stability against three-dimensional perturbations is assured by rotation and the finiteness of *H* that leads to the *k*∥ = 0 modes being isolated [16]. As a result they cascade energy to ever smaller *k*⊥ with no channel to return this energy back to the small scales.

An important outcome from this picture is that energy fluxes obtained from
horizontal averages can significantly overestimate the inverse energy flux. In the
presently examined simulation Π<sub>T</sub>(k<sub>⊥</sub>) was five times larger than the true inverse
flux Π<sub>T</sub>(k). This important result limits the observational estimates Π<sub>T</sub>(k). Most
present estimates of the inverse flux are based on averages of two-dimensional
slices obtained by satellite images. They thus contain no information on the fields
variation in the third direction and as a result it is Π<sub>T</sub>(k<sub>⊥</sub>) that is measured and not
Π<sub>T</sub>(k) that represents the true inverse flux. As a reference, and for comparisons
with observations, we provide as supplementary material (see Fig. S1 in [30])
spectra and third-order structure functions along horizontal tracks, such as those
resulting from airplane or satellite measurements. The structure functions display
a change in sign indicative of an inverse cascade, but overestimate the inverse flux
just as horizontal averages do. Thus these estimates of the inverse flux could be
significantly larger than their true value.
# 4 Supplementary material

## Materials and Methods

The Boussinesq equations (1) and (2) in the main text describe a rotating and stably stratified non-hydrostatic incompressible fluid. The pressure $P$ provides the correction to the hydrostatic pressure. The field $\mathbf{u} = (u, v, w)$ is the full three-dimensional velocity field, and the pointwise fluid density is $\rho = \rho_0 + z(d\bar{\rho}/dz) + \rho'$, where $\rho_0$ is the mean fluid density, $d\bar{\rho}/dz$ is the uniform background density gradient, and the density fluctuations $\rho'$ are written in terms of the normalised density variation $\phi$ as $\rho' = [-(\rho_0/g)(d\bar{\rho}/dz)]^{1/2}\phi$, where $g$ is the acceleration of gravity. Here we consider the case of stable stratification, $d\bar{\rho}/dz < 0$, which is the case for a large part of the Earth's atmosphere. The squared Brunt-Väisälä frequency is then $N^2 = -(g/\rho_0)(d\bar{\rho}/dz)$. For large $N$ and in the presence of rotation we expect to approach the quasi-geostrophic limit, where gravity is balanced by pressure variations and vertical velocities and gravity waves are suppressed. In some regions of the atmosphere the stratification can be unstable, $d\bar{\rho}/dz > 0$, which has been the focus of many studies [7, 8, 9, 10, 11]. In this case and for strong rotation a quasi-geostrophic limit can also be defined [9, 10, 11] and even develop inverse energy cascades. It is worth noting that in the unstably stratified quasi-geostrophic regime gravitational modes are not entirely filtered out.

Simulations were done using the GHOST code, which is publicly available at bluehttps://doi.org/10.5281/zenodo.8015308 [39]. Equations are integrated at double precision, using a parallel pseudo-spectral method to compute spatial derivatives.
tives, and a second-order Runge-Kutta method to evolve the fields in time [33].
Aliasing instabilities were controlled using the 2/3 rule. No subgrid model or ef-
fective transport parameters were used, and all relevant spatial and temporal scales
in the system were explicitly resolved.

The forcing was delta-correlated in time, acting only on the horizontal com-
ponents of the velocity (u and v), and on all modes in the Fourier shell k ∈
[0.9, 1.1]k_H. To initialise the flow, simulations were done at increasing spatial
resolutions of 1536² × 48, 3072² × 96, 6144² × 192, and 12288² × 384 grid points.
In the simulations with 1536² × 48 and with 3072² × 96 grid points the flow was
started from rest, and the simulations were spun up for ≈ 70 forcing-scale turnover
times, τ_f = (εk_H²)⁻¹/³, which was sufficient to display a well developed inverse
cascade. The simulation at 6144² × 192 resolution was then done starting from the
re-scaled last output of the previous simulation, running for at least ≈ 25 turnover
times. Finally, the simulation using 12288² × 384 grid points was integrated for
≈ 6 turnover times from the last re-scaled state of the 6144² × 192 run. This time
was sufficient for any transients to die out, and for a clear increase of energy due
to the inverse cascade to be observed. To ensure the flows were spatially resolved
we also performed other simulations at lower resolution using different Reynolds
numbers. To remove initial transients, all analyses were performed using the last
turnover times of each simulation, averaging over time.

We also performed a simulation with a resolution of 1536² × 48 grid points
with no-slip boundary conditions at z = 0 and z = H using a Fourier continu-
ation method [40], to verify that the inverse cascade also develops with different
boundary conditions in the vertical direction. Finally, note that all simulations were stopped at a time such that energy did not accumulate at the lowest horizontal wave numbers, thus avoiding the formation of domain-scale features that would be impacted by the finite domain size.

Three-dimensional visualisations were done with the software VAPOR [35], publicly available at [41].

## Supplementary Text

### Energy decomposition

The total energy per mass unit of the system, $\mathcal{E}_T$, is given by the sum of the kinetic and potential energies, given respectively by $\mathcal{E}_K = \langle |\mathbf{u}|^2 \rangle / 2$ and $\mathcal{E}_P = \langle \phi^2 \rangle / 2$, where $\langle \cdot \rangle$ indicates the volume average over the domain. Alternatively, the energy can be split into the energy contained in gravity wave modes and in quasi-geostrophic modes. Following [42, 43], we write for each wave number $\mathbf{k}$ the modes in terms of the vertical vorticity $\hat{w}_{\parallel}(\mathbf{k})$, the vertical velocity $\hat{u}_{\parallel}(\mathbf{k})$, and the density variation $\hat{\phi}(\mathbf{k})$ as

$$ \mathbf{Z}(\mathbf{k}) = [\hat{w}_{\parallel}(\mathbf{k}), -ik\hat{u}_{\parallel}(\mathbf{k}), -k_{\perp}\hat{\phi}(\mathbf{k})]. $$

Assuming a periodic time dependence $e^{i\sigma t}$ and diagonalizing the linear part of the Boussinesq equations, we obtain three eigenfrequencies

$$ \sigma_{QG}(\mathbf{k}) = 0, \quad \text{and} \quad \sigma_{GW}^{\pm}(\mathbf{k}) = \pm\sqrt{(4\Omega^2 k_{\parallel} + N^2 k_{\perp}^2)/k^2}, $$

with $\sigma_{QG}$ corresponding to the quasi-geostrophic mode and the other two to the
gravity wave modes. The eigenvectors are given by

$$ \mathbf{Z}_{QG}(\mathbf{k}) = [iNk_{\perp}, 0, 2\Omega k_{\parallel}], \quad \text{and} \quad \mathbf{Z}_{GW}^{\pm}(\mathbf{k}) = [i2\Omega k_{\parallel}, -k\sigma^{\pm}(\mathbf{k}), Nk_{\perp}]. $$

The projection $\mathcal{P}(\mathbf{Z})$ of a mode $\mathbf{Z}$ to the quasi-geostrophic modes is given by

$$ \mathcal{P}_{QG}(\mathbf{Z}) = \mathbf{Z}_{QG}[\mathbf{Z}_{QG}^* \cdot \mathbf{Z}] / |\mathbf{Z}_{QG}|^2, $$

and for the gravito-inertial modes by

$$ \mathcal{P}_{GW}(\mathbf{Z}) = \mathbf{Z}_{GW}^{+} \frac{(\mathbf{Z}_{GW}^{+})^{*} \cdot \mathbf{Z}}{|\mathbf{Z}_{GW}^{+}|^2} + \mathbf{Z}_{GW}^{-} \frac{(\mathbf{Z}_{GW}^{-})^{*} \cdot \mathbf{Z}}{|\mathbf{Z}_{GW}^{-}|^2}. $$

Energies and energy spectra can then be calculated using these projections.

## Energy spectra and fluxes

Energy is distributed across scales by nonlinear interactions. The simplest way to quantify the distribution of energy among scales is by the Fourier transformed fields $\hat{\mathbf{u}}(\mathbf{k})$ and $\hat{\phi}(\mathbf{k})$, where the notion of scale is given by the inverse wavenumber $k^{-1}$. The kinetic and potential energy spectra are then defined respectively as

$$ E_K(k_i) = \frac{1}{k_L} \sum_{\mathbf{q} \in S_k} |\hat{\mathbf{u}}(\mathbf{q})|^2, \quad E_P(k_i) = \frac{1}{k_L} \sum_{\mathbf{q} \in S_k} |\hat{\phi}(\mathbf{q})|^2, $$

with $k_L = 2\pi/L$, where $S_k$ is the set of wavenumbers $\mathbf{q}$ that satisfy $k \le |\mathbf{q}| < k + k_L$ for the spherically averaged spectra, $k_{\perp} \le |\mathbf{q}_{\perp}| < k_{\perp} + k_L$ for the cylindrically averaged spectra, and $k_{\parallel} \le |q_{\parallel}| < k_{\parallel} + k_L$ for the plane averaged spectra. For the two dimensional spectra, $S_k$ is the set that satisfies both $k_{\perp} \le |\mathbf{q}_{\perp}| < k_{\perp} + k_L$ and $k_{\parallel} \le |q_{\parallel}| < k_{\parallel} + k_L$. The total energy spectrum is always given by $E_T(k) = E_K(k) + E_P(k)$.
To construct the spectra of gravito-inertial modes and of quasi-geostrophic modes, the same process is repeated after the modes $\hat{u}(\mathbf{k})$ and $\hat{\phi}(\mathbf{k})$ are projected to the relevant manifolds as described before.

The kinetic and potential energy fluxes $\Pi_K(k)$ and $\Pi_P(k)$ through a closed set of wavenumbers $S_k$ are defined as

$$ \Pi_K(k) = \langle \mathbf{u}_{S_k} \cdot (\mathbf{u} \cdot \nabla \mathbf{u}) \rangle, \quad \Pi_P(k) = \langle \phi_{S_k} (\mathbf{u} \cdot \nabla \phi) \rangle, $$

where $\mathbf{u}_{S_k}$ and $\phi_{S_k}$ are obtained by filtering the fields $\mathbf{u}$ and $\phi$ so that only wavenumbers in the set $S_k$ are kept. The total energy flux is then given by $\Pi_T = \Pi_K + \Pi_P$. For the flux through a sphere the set $S_k$ corresponds to the set of wavenumbers $\mathbf{q}$ that satisfy $|\mathbf{q}| \le k$, for the flux through cylinders it corresponds to the set of wavenumbers $\mathbf{q}$ that satisfy $|\mathbf{q}_{\perp}| \le k_{\perp}$, and for flux through planes it corresponds to the set of wavenumbers $\mathbf{q}$ that satisfy $|\mathbf{q}_{\parallel}| \le k_{\parallel}$.

Note that the averages over all these manifolds are defined in spectral space. Thus, spectra and fluxes that depend on $k_{\parallel}$ quantify vertical gradients and structures in real space. Spectra and fluxes that depend on $k_{\perp}$ quantify horizontal gradients (as those obtained in atmospheric measurement from flights). Isotropic spectra and fluxes (i.e., dependent on $k$) mix gradients in all directions, and should be interpreted with care in the presence of anisotropies. However, the latter are the relevant spectral quantities to uniquely quantify cascades from a theoretical point of view.

## Velocity spectra and structure functions along vertical and horizontal tracks

We extracted 384 random horizontal tracks from the numerical data (to mimic,
e.g., observations done along straight flight legs by airplanes, or reconstructions of velocity fields using velocimetry techniques from satellite images). We computed wind velocity (u component) power spectra along each individual track, and third-order structure functions

$$S_3(r) = \langle [u(\mathbf{x} + r\mathbf{e}_x) - u(\mathbf{x})]^3 \rangle,$$

where $e_x$ is the unit vector along the track, and the average is done over all points $\mathbf{x}$ in the horizontal track. This structure function is expected to follow a power law and to be negative at scales in which energy is transfered to smaller scales, and positive when the energy flux is reversed [44, 45].

Figure S1 shows the resulting spectra and structure functions along the synthetic tracks. Note the change in sign of the average of $S_3$ over all tracks at $k_H r \approx 1$, as well as the strong fluctuations in the individual tracks. Assuming $S_3(r) = -(4/5)\epsilon r$ [45], the analysis indicates that the energy flux at large scales is negative (i.e., inverse), and results in an estimation of the inverse flux similar to the one obtained from $\Pi_T(k_\perp)$ (i.e., on an overestimation of the inverse flux). The estimation can be slightly improved by using scaling laws that take into account the flow stratification and anisotropy [45], but as shown in [46], in the anisotropic case proper determination of energy fluxes from $S_3$ requires averaging along multiple directions.

Comparisons with observations should take into account the difficulty in obtaining airplane data along long legs, and the dependence of satellite data on indirect methods to reconstruct velocity fields. In spite of these limitations, changes
in the scaling law in airplane data [47] and positive third-order structure func-
tions in airplane and satellite data [44, 48, 49] have been reported. The scale at
which scaling laws or the sign of $S_3$ change in these observations depends on
atmospheric conditions, varying also across different samples and types of obser-
vations. While the results presented here can provide a possible explanation for
this behavior, care must be taken to also consider fluctuations and noise in the
data, the already mentioned overestimation of the inverse energy flux resulting
solely from horizontal measurements, and sources of energy at large atmospheric
scales that can obfuscate inverse fluxes.
<figure><img src="image_5.png" /><figcaption>Figure 4: **Power spectrum and third-order structure function along synthetic tracks.** Top: Power spectrum of the *u* component of the velocity field averaged along horizontal tracks in the *x* direction (red line). The power spectrum obtained considering all directions and points in the simulation is shown in blue line. Gray lines show several spectra obtained for individual tracks, to illustrate variability. Characteristic scales and power laws are shown as a reference. Bottom: Third-order structure function averaged over horizontal tracks (solid and dashed, depending on the sign). Several structure functions for individual tracks are also shown in gray to illustrate variability. A power law is indicated as a reference in the range with negative energy flux.</figcaption></figure>
**Movie S1. (separate file)**

A movie showing density fluctuations at different heights in the entire domain is provided as Supplementary Material. The movie features the development of large-scale structures as well as the formation of small scales. Note the emergence of macroscopic cyclones and anti-cyclones with relative lower and higher densities, as well as small-scale overturning events reminiscent of Kelvin-Helmholtz billows.

**References**

[1] L. Onsager, *Il Nuovo Cimento* (1943-1954) **6**, 279 (1949).

[2] R. H. Kraichnan, *Physics of Fluids* **10**, 1417 (1967).

[3] R. H. Kraichnan, D. Montgomery, *Reports on Progress in Physics* **43**, 547 (1980).

[4] C. E. Leith, *Physics of Fluids* **11**, 671 (1968).

[5] G. K. Batchelor, *Physics of Fluids* **12**, II (1969).

[6] J. G. Charney, *Journal of the Atmospheric Sciences* **28**, 1087 (1971).

[7] P. H. Roberts, *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences* **263**, 93.

[8] F. H. Busse, *Journal of Fluid Mechanics* **44**, 441 (1970).
[9] A. M. Rubio, K. Julien, E. Knobloch, J. B. Weiss, *Physical Review Letters* **112**, 144501 (2014).

[10] B. Favier, L. J. Silvers, M. R. E. Proctor, *Physics of Fluids* **26**, 096605 (2014).

[11] L. Siegelman, et al., *Nature Physics* **18**, 357 (2022).

[12] J. Vanneste, I. Yavneh, *Journal of the Atmospheric Sciences* **61**, 211 (2004).

[13] J. Vanneste, *Annual Review of Fluid Mechanics* **45**, 147 (2013).

[14] J. Thomas, D. Daniel, *Journal of Fluid Mechanics* **911**, A60 (2021).

[15] J. Thomas, D. Daniel, *Journal of Fluid Mechanics* **902**, A7 (2020).

[16] A. Alexakis, L. Biferale, *Physics Reports* **767**, 1 (2018).

[17] L. M. Smith, J. R. Chasnov, F. Waleffe, *Physical Review Letters* **77**, 2467 (1996).

[18] A. Celani, S. Musacchio, D. Vincenzi, *Physical Review Letters* **104**, 184506 (2010).

[19] K. Seshasayanan, S. J. Benavides, A. Alexakis, *Physical Review E* **90**, 051003 (2014).

[20] S. J. Benavides, A. Alexakis, *Journal of Fluid Mechanics* **822**, 364 (2017).

[21] A. Sozza, G. Boffetta, P. Muratore-Ginanneschi, S. Musacchio, *Physics of Fluids* **27**, 035112 (2015).
[22] A. van Kan, A. Alexakis, *Journal of Fluid Mechanics* **899**, A33 (2020).

[23] R. Marino, P. D. Mininni, D. Rosenberg, A. Pouquet, *EPL (Europhysics Letters)* **102**, 44006 (2013).

[24] A. Pouquet, R. Marino, *Physical Review Letters* **111**, 234501 (2013).

[25] R. Marino, A. Pouquet, D. Rosenberg, *Physical Review Letters* **114**, 114504 (2015).

[26] D. Byrne, J. A. Zhang, *Geophysical Research Letters* **40**, 1439 (2013).

[27] G. P. King, J. Vogelzang, A. Stoffelen, *Journal of Geophysical Research: Oceans* **120**, 346 (2015).

[28] D. Balwada, J.-H. Xie, R. Marino, F. Feraco, *Science Advances* **8**, eabq2566 (2022).

[29] R. Young, P. L. Read, *Nature Physics* **13**, 1135 (2017).

[30] Materials and methods are available as supplementary materials.

[31] J. Pedlosky, *Geophysical Fluid Dynamics* (Springer, 1987).

[32] G. K. Vallis, *Atmospheric and Oceanic Fluid Dynamics* (Cambridge University Press, 2017).

[33] P. D. Mininni, D. Rosenberg, R. Reddy, A. Pouquet, *Parallel Computing* **37**, 316 (2011).
[34] H.-L. Liu, P. B. Hays, R. G. Roble, *Journal of the Atmospheric Sciences* **56**, 2152 (1999).

[35] J. Clyne, P. Mininni, A. Norton, M. Rast, *New Journal of Physics* **9**, 301 (2007).

[36] G. Nastrom, K. Gage, W. Jasperson, *Nature* **310**, 36 (1984).

[37] P. Billant, J.-M. Chomaz, *Physics of Fluids* **13**, 1645 (2001).

[38] H. P. Greenspan, *The theory of rotating fluids* (Cambridge University Press, 1968).

[39] D. Rosenberg, P. Mininni, P. C. D. Leoni, pmininni/ghost: Ghost master, version 2.0.0, zenodo (2023); https://doi.org/10.5281/zenodo.8015308.

[40] M. Fontana, O. P. Bruno, P. D. Mininni, P. Dmitruk, *Computer Physics Communications* **256**, 107482 (2020).

[41] S. e. a. Pearse, Ncar/vapor: Vapor, version 3.8.1, zenodo (2023); https://doi.org/10.5281/zenodo.7779648.

[42] P. Bartello, *Journal of Atmospheric Sciences* **52**, 4410 (1995).

[43] C. Herbert, A. Pouquet, R. Marino, *Journal of Fluid Mechanics* **758**, 374 (2014).

[44] J. Y. Cho, E. Lindborg, *Journal of Geophysical Research: Atmospheres* **106**, 10223 (2001).
[45] E. Lindborg, J. Y. Cho, *Journal of Geophysical Research: Atmospheres* **106**, 10233 (2001).

[46] Q. Nie, S. Tanveer, *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* **455**, 1615 (1999).

[47] P. Trent Vonich, G. J. Hakim, *Journal of the Atmospheric Sciences* **75**, 2523 (2018).

[48] P. Héas, E. Mémin, D. Heitz, P. D. Mininni, *Tellus A: Dynamic Meteorology and Oceanography* **64**, 10962 (2011).

[49] P. Rodriguez Imazio, A. Dörnbrack, R. D. Urzua, N. Rivaben, A. Godoy, *Journal of Geophysical Research: Atmospheres* **127**, e2021JD035908 (2022).

## Acknowledgements

**Funding:** Computer resources in Joliot-Curie at CEA were provided by PRACE (research project ID 2020235566), and by GENCI (allocation No. A0110506421). 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 Studienstiftung des deutschen Volkes, by the National Science Foundation (Grants DMS-2009563, DMS-2308337) and by the German Research Foundation (Projektnummer: 522026592). The simulation output was analyzed on HPC facilities at the École Normale Superieure in Paris
(France), at École Centrale de Lyon (PMCS2I) in Ecully (France), and at Departamento de Fisica (FCEN, UBA) in Argentina.