# Global climate modeling of Saturn's atmosphere. Part III: Global statistical picture of zonostrophic turbulence in high-resolution 3D-turbulent simulations

Simon Cabanes, Aymeric Spiga, Roland M. B. Young

► To cite this version:

Simon Cabanes, Aymeric Spiga, Roland M. B. Young. Global climate modeling of Saturn's atmosphere. Part III: Global statistical picture of zonostrophic turbulence in high-resolution 3D-turbulent simulations. Icarus, 2020, 345, 10.1016/j.icarus.2020.113705. insu-03726979

HAL Id: insu-03726979

https://insu.hal.science/insu-03726979v1

Submitted on 19 Sep 2025

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
# Global climate modeling of Saturn's atmosphere.
## Part III: Global statistical picture of zonostrophic turbulence in high-resolution 3D-turbulent simulations

Simon Cabanes<sup>a,b,*</sup>, Aymeric Spiga<sup>a,c</sup>, Roland M. B. Young<sup>a,d</sup>

## Abstract

We conduct an in-depth analysis of statistical flow properties calculated from the reference high-resolution Saturn simulation obtained by global climate modelling in Part II. In the steady state of this reference simulation, strongly energetic, zonally dominated, large-scale structures emerge, which scale with the Rhines scale. Spectral analysis reveals a strong anisotropy in the kinetic energy spectra, consistent with the zonostrophic turbulent flow regime. By computing spectral energy and enstrophy fluxes we confirm the existence of a double cascade scenario related to 2D-turbulent theory. To diagnose the relevant 3D dynamical mechanisms in Saturn's turbulent atmosphere, we run a set of four simulations using an idealized version of our Global Climate Model devoid of radiative transfer, with a well-defined Taylor-Green forcing and over several rotation rates (4, 1, 0.5, and 0.25 times Saturn's rotation rate). This allows us to identify dynamics in three distinctive inertial ranges: (1) a “residual-dominated” range, in which non-axisymmetric structures dominate with a −5/3 spectral slope; (2) a “zonostrophic inertial” range, dominated by axisymmetric jets and characterized by the pile-up of strong zonal modes with a steeper, nearly −3, spectral slope; and (3) a “large-scale” range, beyond Rhines' typical length scale, in which the reference Saturn simulation

<footnote><sup>a</sup>Laboratoire de Météorologie Dynamique (LMD/IPSL), Sorbonne Université, Centre National de la Recherche Scientifique, École Polytechnique, École Normale Supérieure, Paris, France.</footnote>

<footnote><sup>b</sup>DICEA, Sapienza Università di Roma, Via Eudossiana 18, 00184 Rome, Italy.</footnote>

<footnote><sup>c</sup>Institut Universitaire de France (IUF), Paris, France</footnote>

<footnote><sup>d</sup>Department of Physics & National Space Science and Technology Center, UAE University, Al Ain, United Arab Emirates</footnote>

<footnote>*Corresponding author. Email address: cabanes.simon@gmail.com</footnote>
and our idealized simulations differ. In the latter range, the dynamics is dom-
inated by long-lived zonal modes 2 and 3 when a Saturn-like seasonal forcing
is considered (reference simulation), and a steep energetic decrease with the
idealized Taylor-Green forcing. Finally, instantaneous spectral fluxes show
the coexistence of upscale and downscale enstrophy/energy transfers at large
scales, specific to the regime of zonostrophic turbulence in a 3D atmosphere.

***Keywords:*** Global climate modeling, zonostrophic regime, Saturn's zonal jets, spectral analysis.

**1. Introduction**

The self-organization of large-scale coherent flows is ubiquitous in nature and is an essential, cutting edge research topic in fundamental fluid dynamics, geophysics and planetary sciences (Galperin and Read, 2019). The strong zonal (i.e. east-west) winds observed on the gas giants, i.e. Jupiter, Saturn, Uranus and Neptune (Ingersoll, 1990) are unmistakable markers of powerful, self-organized, atmospheric dynamics, insofar as they have been difficult to access by direct measurements.

By tracking cloud motions in the upper atmosphere, ground-based obser-
vations and the successive Voyager, Cassini, New Horizons, and Juno mis-
sions have provided detailed information on the wind patterns of Jupiter
and Saturn, including temporal variability and vertical profiling (e.g., Porco
et al., 2003; Sánchez-Lavega et al., 2004; Vasavada and Showman, 2005;
García-Melendo et al., 2010). Each of these planets exhibits a system of
axisymmetric (i.e. axial symmetry relative to the rotational spin axis) jet-
streams, confined in latitude and almost invariant in azimuth. These jets are
latitudinally-varying winds that alternate between prograde (eastward) and
retrograde (westward) flow with respect to the apparent mean planetary ro-
tation. Latitudinal jet profiles exhibit qualitative similarities on both Jupiter
and Saturn: there exists a broad, fast prograde equatorial jet and numerous
narrower, weaker prograde jets at higher latitudes. Despite these similari-
ties, the zonal jets are considerably broader, faster, and fewer in number on
Saturn than on Jupiter (Showman et al., 2017).

Apart from the dominant axisymmetric flow, hundreds of non-axisymmetric
eddies develop (convective storms, waves, vortices, and turbulence). Jupiter's
Great Red Spot, an anticyclonic vortex (Fletcher et al., 2010; Simon et al.,
2014), and Saturn's polar hexagon, a wavenumber-6 disturbance superim-
posed on the northern circumpolar jet (Baines et al., 2009; Sánchez-Lavega
et al., 2014; Rostami et al., 2017) constitute the largest and the most long-
lived examples. In recent decades, much effort has been made by theoreti-
cians, experimentalists, and numericists to understand the relevant dynami-
cal mechanisms responsible for the axisymmetric jets, the non-axisymmetric
eddies, and the interactions between them (Salyk et al., 2006; Del Genio and
Barbara, 2012; Galperin et al., 2014).

How do axisymmetric jets emerge in gas giants? It is well accepted that in three-dimensional (3D) systems, turbulent kinetic energy is typically transferred from larger to smaller scales through a direct energy cascade (the Kolmogorov cascade). However, this process can be reversed under certain conditions favoring coherent flow. In other words, non-axisymmetric turbulent flows (e.g. caused by flow baroclinicity, Salmon, 1980) can feed the large-scale jets. This is the case in both quasi-two-dimensional (2D) domains, where the flow is geometrically confined in depth to a shallow layer, or under rapid rotation in 3D domains, where the flow becomes invariant along the rotational spin axis and is dynamically confined to quasi-2D geostrophy. Two types of models have been proposed for jet formation in planetary atmospheres: shallow models where jet forcing is confined to the atmospheric weather layer (Cho and Polvani, 1996; Lian and Showman, 2008; Schneider and Liu, 2009; Liu and Schneider, 2010; Young et al., 2019; Spiga et al., 2020) and deep models where jet forcing extends into the convective molecular envelope of the planet (Christensen, 2002; Heimpel et al., 2005; Gastine et al., 2014; Cabanes et al., 2017). In both cases, flow bi-dimensionalization transfers kinetic energy from small to large scales via the so-called inverse energy cascade (Kraichnan, 1967). It is only recently that this theoretical framework has been generalized to 3D rotating turbulence by Sukoriansky and Galperin (2016). Using analytical tools, they showed the transition from a 3D direct energy cascade to rotation-dominated (quasi-2D) turbulence leading to an inverse energy cascade. This duality between 2D and 3D turbulence as well as geometrical and dynamical flow confinement, together with the associated transition from direct to inverse energetic cascade, is at the core of fundamental geophysical fluid dynamics (Smith et al., 1996; Pouquet and Marino, 2013; Deusebio et al., 2014).

Spherical geometry also carries an additional background potential vor-
ticity gradient $\beta$, which preferentially channels kinetic energy in the zonal
direction, promoting large-scale zonal motions (Vallis and Maltrud, 1993;
Cho and Polvani, 1996). The planetary vorticity $\beta = (2\Omega/a) \cos \varphi$ directly
results from the variation of the Coriolis parameter with latitude $\phi$, the rotation rate $\Omega$, and planetary radius $a$. This coupled action of the turbulent inverse cascade and of the $\beta$-effect in planetary atmospheric envelopes links large-scale features to smaller scale dynamics through the continuous process of scale-to-scale interactions. These specific planetary conditions define the regime of zonostrophic turbulence, coined by Sukoriansky et al. (2002).

In the framework of 2D-turbulence, Sukoriansky et al. (2002) generalized the double cascade scenario of Kolmogorov-Batchelor-Kraichnan (KBK) (i.e. the inverse cascade of energy and the direct cascade of enstrophy) to turbulence on the surface of a rotating sphere. They showed that 2D-turbulent zonostrophic regimes display universal kinetic energy spectra with a steep $-5$ slope in their axisymmetric component, while in the non-axisymmetric component the classical KBK scaling prevails with a $-5/3$ slope. This spectral anisotropy, as well as the presence of an inverse energy cascade, has been detected in Jupiter’s tropospheric winds using Cassini near-IR imaging (Choi and Showman, 2011; Galperin et al., 2014; Young and Read, 2017).

Similarly, this zonostrophic signature has been identified in 2D-turbulent numerical simulations both with a $\beta$-plane approximation (Chekhlov et al., 1996) and on the surface of a rotating sphere (Huang et al., 2001; Galperin et al., 2006). These models reproduce the basic dynamical features of generic gas giant atmospheres without reproducing specific planets, yet they lack the physical realism of real planetary systems as they explicitly assume 2D-turbulent flows. The emergence of Global Climate Models (GCMs) for gas giants (Dowling et al., 2006; Schneider and Liu, 2009; Liu and Schneider, 2010; Lian and Showman, 2010; Young et al., 2019; Spiga et al., 2020), where a thin planetary shell is simulated by solving for 3D primitive equations of motion with increasingly realistic physical parametrizations (e.g., radiative transfer, Guerlet et al., 2014; Spiga et al., 2020), now allows for in-depth spectral analysis of simulated giant planets in order to diagnose the relevant dynamical mechanisms in 3D-turbulent zonostrophic flow.

To step in this direction, the present paper (Part III of this series) uses high-resolution multi-annual 3D numerical simulations obtained by the Saturn DYNAMICO GCM described by Spiga et al. (2020) (Part II of this series). This GCM is designed to explore Saturn’s tropospheric and stratospheric dynamics with a new icosahedral dynamical core (DYNAMICO, Dubos et al., 2015) and realistic radiative transfer (Guerlet et al., 2014, Part I of this series). Following the approach that Augier and Lindborg (2013) applied to terrestrial Global Climate Models, we compute kinetic energy spectra as
well as spectral energy and enstrophy fluxes from instantaneous wind maps. This is done in two stages. First, we compute the spectral analysis of the Saturn reference simulation (SRS in the text below) described in Part II by Spiga et al. (2020), where the forcing gives rise to banded jets resulting from the growth of baroclinic eddies at small scales. Second, we run a set of idealized simulations, using only the DYNAMICO dynamical core, with a well-defined Taylor-Green forcing but with radiative transfer being switched off, which suppresses flow baroclinicity caused by equator-to-pole gradients. In these idealized simulations we vary the planetary rotation rate (as done by, e.g., Chemke and Kaspi, 2015) to bracket the properties of Saturn's zonostrophic regime, and to diagnose the relevant dynamical processes in a 3D-turbulent atmospheric layer without the effects of intermittency and spatial anisotropy induced by baroclinic forcing.

In all these simulations we identify the double cascade predicted by KBK for 2D-turbulence, as well as zonostrophic spectral anisotropy. In the SRS, semi-annual and annual energy cycles are identified. Our analysis reveals that the dynamics exhibit a transition across the Rhines scale $L_R$ (Rhines, 1975), showing a large-scale dynamical regime ($L > L_R$) dominated by long-lived energetic modes, and a small-scale dynamical regime ($L < L_R$) that is highly sensitive to the seasonal cycle. In turn, the idealized simulations emphasize the peculiar dynamics that develop within the zonostrophic inertial range, where axisymmetric kinetic energy substantially dominates over non-axisymmetric kinetic energy. As the rotation rate increases, the dynamical confinement (i.e. quasi-2D geostrophy) as well as the inverse energy cascade become more efficient, leading towards a steeper energy spectra to which the classical picture of 2D-turbulence is not applicable.

In Section 2 we remind the reader about the key characteristics of the SRS (Saturn reference simulation) and we present the numerical setup used to run our idealized simulations. The theoretical framework of zonostrophic turbulence is introduced in Section 3 as well as the typical scaling characteristics of 3D rotating turbulence. The statistical tools are detailed in Section 4. Statistical analysis is first applied to the SRS in Section 5, then to our idealized setup in Section 6. Conclusions and perspectives are in Section 7.

## 2. Numerical setup

The first part of our analysis uses the exact same SRS as is described and analyzed in detail in Spiga et al. (2020). Here we remind the reader about the
main characteristics of this simulation. In the present manuscript, latitude is denoted as $\phi$, longitude as $\lambda$, radial/vertical as $z$ and $g = 10.4 \text{ m s}^{-2}$ is the gravitational acceleration on Saturn. We use Saturn's radius $a = 58232000 \text{ m}$ and rotation rate $\Omega = 1.65121 \times 10^{-4} \text{ s}^{-1}$.

A Global Climate Model (GCM) consists of coupling a generic hydrodynamical solver (usually named the dynamical core) with physical parametrizations specific to a planet, aimed at emulating the processes unresolved by the dynamical core (radiative transfer, small-scale mixing, cloud formation, etc. ...). Developed at Laboratoire de Météorologie Dynamique, DYNAMICO is a new dynamical core that solves the primitive equations of motion with a hybrid mass-based vertical coordinate. Equations are solved in the limit of hydrostatic dynamics with the shallow-atmosphere approximation (i.e. reduced Coriolis parameter). Spatial discretizations are formulated following an energy-conserving three-dimensional Hamiltonian approach (Dubos et al., 2015). The DYNAMICO dynamical core uses a quasi-uniform icosahedral mapping of the planetary sphere to ensure excellent conservation and scalability properties in a massively parallel configuration (Dubos et al., 2015) with the projection from the unstructured icosahedral grid to a cylindrical latitude/longitude grid performed at runtime (see Spiga et al., 2020, for details and references). This approach enables us to carry out computationally-intensive simulations combining high horizontal resolutions (here $1/2^\circ$ in latitude and longitude) with a long integration time of more than ten Saturn years. Time integration is done by an explicit Runge-Kutta scheme.

Subgrid-scale dissipation in the horizontal is solved by adding a hyperdiffusion term to the vorticity, divergence and temperature equations:

$$ \left[ \frac{d\mathcal{F}}{dt} \right]_{\text{dissip}} = \frac{(-1)^{q+1}}{\tau} \nabla^{2q} \mathcal{F} \quad (1) $$

where $\mathcal{F}$ can be divergence, vorticity, or potential temperature, $q$ is the order of dissipation, and $\tau$ the damping timescale at the smallest horizontal scale resolved by the model. The selection of appropriate values for $q$ and $\tau$ and their impact on the global dynamics remains empirical, as discussed in the appendices of Spiga et al. (2020). For our idealized simulations described in section 6, we choose the same values as what the ones used by Spiga et al. (2020) for the SRS, i.e. $q = 2$ and $\tau = 10000$. Setting these hyperdiffusion parameters is a trade-off between the stability of the model integration and the degree to which small-scale perturbations are resolved.
In the vertical dimension (i.e. along radial coordinates, which correspond to the vertical in our shallow-atmosphere simulations), both the SRS in Spiga et al. (2020), and the idealized GCM simulations presented in section 6, solve the equations at 32 levels, from 3 bar to 1 mbar, from the upper troposphere to the middle stratosphere. A stress-free boundary condition is imposed at the uppermost layer, instead of an absorbing sponge layer that may result in spurious effects due to non-conservation of angular momentum (Shaw and Shepherd, 2007). At the lowest layer (i.e. at the 3-bar level), the bounding condition is stress-free and a bottom drag $u_b = -\alpha u$ is imposed with a frictional rate $\alpha = 1/\tau_R$ similar to Schneider and Liu (2009) and Liu and Schneider (2010). This large-scale frictional dissipation emulates to zeroth order, at the bottom of our model, the magneto-hydrodynamic drag (MHD) acting on jets extending at the bottom of the atmospheric layer (Liu et al., 2008) several thousand kilometers below the visible cloud layer (Galanti et al., 2019). Improving on this admittedly simplistic bottom boundary condition is an entire research goal on its own and is left for future investigations. In our model, we set a frictional region with $\tau_R = 90$ Earth days, with a non-frictional region near the equator $-33^{\circ} < \varphi < 33^{\circ}$ with $\tau_R = \infty$, as in Schneider and Liu (2009); Liu and Schneider (2010). In the vertical, subgrid-scale dissipation is handled by a combination of a Mellor-Yamada diffusion scheme for small-scale turbulence, and a dry convective adjustment scheme for organized turbulence (e.g. large-scale convective plumes).

The physical parametrizations used in our GCM have been specifically developed to simulate the physical properties of Saturn's troposphere and stratosphere, in particular a seasonal radiative transfer model optimized for Saturn, complete with internal heat flux and Saturn's ring shadowing, which has been validated using available measurements, as detailed by Guerlet et al. (2014) and Spiga et al. (2020). This is the key setting on which the SRS and the idealized simulations presented in section 6 differ. We note that the current version of the model does not parametrize moist convection, which occurs in giant planets' atmospheres (Gierasch et al., 2000) and carries substantial energy into the jet system (Sayanagi and Showman, 2007). The model also inhibits large vertical gradients of potential temperature (i.e. along a planetary radius) by linearizing deep convection around a mean vertical temperature profile (the convective adjustment scheme). These physical limits, together with the shallow-atmosphere approximation, are possibly responsible for some discrepancies between direct planetary observations and DYNAMICO simulated wind fields in a Saturn-like configuration like the
SRS (see also section 4.4 in Spiga et al., 2020). Future implementations of our Saturn DYNAMICO model will address these issues. Aware of those limitations, the present work mainly focuses on DYNAMICO's dynamical core and its capability to simulate the zonostrophic inertial range relevant for the formation of jets.

### 3. Scaling and kinetic energy distribution in 3D zonostrophic turbulence

Global Climate Models are powerful tools for modelling the complexity of planetary atmospheres in a realistic fashion. The primitive equations of motion are solved in a 3D-spherical domain, for a divergent, stratified flow in which energy is provided by the seasonal cycle of incoming sunlight (all planets) and internal heat flux (giant planets). Quasi-steady jets arise in GCM simulations after several simulated Saturn years as a result of the non-linear dynamics resolved by the model. In order to diagnose these mechanisms it is relevant to analyze the kinetic energy distribution in Fourier space. Spectral analysis has the advantage that each dynamical step leading to the formation of equilibrated jets can be interpreted in terms of typical length scales or, equivalently, wavenumbers.

When spherical harmonic functions are invoked to compute energy spectra, the basis functions $Y_n^m(\varphi, \lambda)$ are eigenfunctions of the horizontal Laplacian operator $\nabla_h^2 Y_n^m = -n(n+1)Y_n^m/a^2$, and depend on non-dimensional total and zonal indices $n$ and $m$ respectively (Boer, 1983; Boer and Shep herd, 1983). The planetary spin axis defines the zonal index $m = 0$ to be the axisymmetric mode that characterizes the zonal jets. Indices $m \neq 0$ de- fine the residual non-axisymmetric part of the flow (eddies). At any typical length scale in latitude $L$, one can attribute a typical total index $n$ given by the relation $n = 2\pi a/L$, since $\varphi \in [0, \pi]$. We represent in Fig. 1 the theoret- ical distribution of kinetic energy as a function of total index $n$ in the same fashion as presented in Galperin et al. (2010). This diagnostic aims to give a global picture of the zonostrophic regime with a scaling suitable for our Saturn reference simulation. This scaling is now described.

#### 3.1. Typical scales and the 2D horizontal anisotropy

The absorption of sunlight in a stably stratified, radiative, motionless at- mosphere leads to the development of an unstable meridional temperature gradient, which drives small-scale baroclinic eddies. These eddies result from
<figure><img src="image_2.png" /><figcaption>Figure 1: Theoretical energy distribution as a function of the total index *n* presented in a similar fashion to Galperin et al. (2010) and Galperin et al. (2019, Chapter 13). The scaling corresponds to our Saturn reference simulation. Notations are described in extenso in the text of sections 3.1 and 3.2.</figcaption></figure>
the conversion of potential energy, stored in density fronts, into kinetic energy at typical horizontal (i.e. along the spherical surface) length scales close to Rossby radii of deformation, $L_D = N_{BV}H/d\pi f$, where $d$ are positive integers defining baroclinic modes, $N_{BV} = \sqrt{-g/\rho}\nabla\rho$ is the Brunt-Väisälä frequency, and $\nabla\rho$ is the density change over depth $H$. Thus, in SRS, kinetic energy is presumably created at the typical baroclinic index $n_D \ge 2\pi a/L_D$, corresponding to the shadowed area in Fig. 1. In this figure, we also locate the smallest resolved scale corresponding to $1/2^\circ$ in latitude; this indicates that only the first Rossby deformation radius is resolved in our simulations.

It is commonly argued that rapid planetary rotation “barotropizes” baroclinic eddies, leading to an inverse kinetic energy cascade towards larger scales (Charney, 1971; Rhines, 1975; Salmon, 1980). Such a picture of energy source and energy transfer has been observed in terrestrial atmospheric and oceanic reanalysis data (Chemke and Kaspi, 2016; Chemke et al., 2016) as well as in idealised GCM simulations (Chemke and Kaspi, 2015). In the specific configuration of shallow-rotating atmospheric layer, dynamical confinement (i.e. barotropization) and geometrical confinement (i.e. weak depth extension) act together to suppress enstrophy production at large scales and to reinforce the inverse cascade at the expense of the direct one. The coexistence of geometrical and dynamical flow confinement is a topic of recent investigation insofar restricted to idealized numerical studies in a triply-periodic domain (Deusebio et al., 2014). Flow confinement means horizontal anisotropy, namely the inhibition of flow evolution in the vertical, which splits the dynamics into two subranges of 3D isotropic and quasi-2D anisotropic turbulence, as summarized in Fig. 1.

Studies suggest the existence of two typical length scales, relative to both types of confinement, that bound the region of flow anisotropisation, and therefore bound the inverse cascade.

1. The first length scale is set by the aspect ratio $S = L_f/L_z$ between the typical length in the compacted direction $L_z$ and the scale of the forcing $L_f$ (Smith et al., 1996). Celani et al. (2010) report that the energy transfer rate $\epsilon$ of the inverse cascade diminishes as the aspect ratio increases, and eventually vanishes for $S = 1/2$. Even though this value is not expected to be universal, it is useful to define a typical
latitudinal index

$$n_S = 2\pi a/H \qquad (2)$$

corresponding to this threshold (here, $S = 1/2$ and $L_Z = H$). In the scale range $n < n_S$, geometrical confinement is suspected to lead to upscale energy transfer.

2. A second length scale is provided by setting the non-dimensional Rossby number $Ro(n) = \epsilon^{1/3}/\Omega(\pi a/n)^{2/3}$ equal to one (Zeman, 1994). This defines a typical index

$$n_{Ro} = \pi a (\Omega / \epsilon^{1/3})^{3/2} \qquad (3)$$

that sets a scale range $n < n_{Ro}$ in which the flow is sensitive to the Coriolis effect, and develops barotropic modes (i.e dynamical confinement) (Smith and Waleffe, 1999; Deusebio et al., 2014). In this study, we estimate the Rossby number by integrating over all total indices $n$. Also, we will refer interchangeably to the Rossby index and the Rossby number to account for the strength of the dynamical confinement.

Fig. 1 shows that both the aspect ratio and Rossby indices ($n_S$ and $n_{Ro}$) are much smaller than the smallest scales resolved in our GCM (spectral analysis section 5 contains an estimate for $\epsilon$). Consequently, as is the case for all existing GCMs for gas giants, our simulations (with 1/2° horizontal resolution) only partially reproduce the inverse cascade inertial range and do not reproduce the 3D isotropic domain at all. Note that because kinetic energy is created within the 2D horizontal anisotropic scale-range, an upscale energy cascade is expected.

## 3.2. Spectral energy budget and zonal anisotropy

The framework of zonostrophic turbulence extends Kraichnan's classical definition of the energy transfer range (Kraichnan, 1971) to 2D turbulent flows with a $\beta$-effect (Galperin et al., 2010). The inverse energy cascade is zonally anisotropized, leading to zonal and residual spectra in the following universal form (Sukoriansky et al., 2002),

$$E_Z(n) = C_Z \beta^2 n^{-5}, \quad C_Z \approx 0.5, \qquad (4a)$$

$$E_R(n) = C_R \epsilon^{2/3} n^{-5/3}, \quad C_R \approx 5-6. \qquad (4b)$$
The residual spectrum (black dashed line in Fig. 1) closely follows the classical Kolmogorov-Batchelor-Kraichnan (KBK) theory of 2-D isotropic turbulence, while the zonal spectrum (red dashed line in Fig. 1) clearly diverges from the non-axisymmetric modes with its steeper −5 slope. The intersection of the zonal and residual spectra defines the transitional index

$$n_{\beta} = (C_{Z}/C_{R})^{3/10} (\beta^{3}/\epsilon)^{1/5} \qquad (5)$$

at which the turbulent and Rossby wave frequencies are equal. This length scale compares two processes that are competing in the system: the flow capacity to transfer energy uspscale through the value of $\epsilon$ (which depends on how flow bi-dimensionalization is efficiently set) versus its ability to channell energy in the zonal direction through parameter $\beta$ (here, the parameter $\beta$ is estimated at mid-latitude $\phi = 45^{\circ}$). In simulations assuming 2D barotropic flows, Sukoriansky et al. (2008) and Galperin et al. (2010) find $(C_Z/C_R)^{3/10} \approx 0.5$, while Galperin et al. (2014) and Cabanes et al. (2017, 2018) find $\approx 0.75$ using observed Jupiter wind fields and laboratory experiments, respectively.

In classical isotropic 2D-turbulence on a non-rotating sphere, the inverse energy cascade ranges between the injection scale $n_i$ (here $\approx n_D$), at which the flow is energized, and the large-scale friction $n_{fr}$, at which energy is dissipated. In the case of a rotating sphere, the additional $\beta$-effect splits energy transfer into two subranges $n_{fr} < n < n_{\beta}$ and $n_{\beta} < n < n_i$ (Huang et al., 2001; Galperin et al., 2010). The former corresponds to the zonostrophic inertial range where the axisymmetric flow dominates and forms a system of multiple jets (Sukoriansky et al., 2002, 2008; Galperin et al., 2010). We call the latter the residual-dominated inertial range, where the $\beta$-effect is weaker and the flow is dominated by the residual spectrum. These two subranges will be of crucial importance in the following study, as they are both fully resolved by the model. In flows with a strong $\beta$-effect, it is well accepted that the Rhines total index

$$n_R = a(\beta/2U)^{1/2} \qquad (6)$$

is a good approximation of the frictional scale at which the upscale cascade is arrested (Sukoriansky et al., 2007). Rhines' index also gives an estimate of the jets' size, with $U$ the root-mean-square (RMS) of the total velocity. The extension of the zonostrophic inertial range, and correspondingly the
strength of the zonostrophic regime, can be estimated by the zonostrophy index $R_β = n_β/n_R$. It is generally accepted that zonostrophic turbulence emerges for $R_β ≥ 2.5$ (Galperin et al., 2006).

Finally, there is an ultimate length scale at which $β$ becomes sufficiently weak that the flow reaches the domain of 2D horizontal anisotropy where the zonal and residual energy spectra are undifferentiated. This transition is observed in the modal spectrum, which is defined as

$$ \mathcal{E}(n, m) = E_R(n)/(2n + 1) = (1/2)\epsilon^{2/3}n^{-8/3} \quad (7) $$

and is the green dashed line in Fig. 1. Equating the modal spectrum to Eq. (4a), we obtain an estimate for the zonal index:

$$ n_z = (C_A/C_R)^{3/10}(2n_β/Δm)^{3/7}n_β \quad (8) $$

at which the flow transits to non-zonostrophic turbulence (here $Δm = 1$). Three primary spectral ranges are then identified as the total index $n$ increases: the zonal anisotropic, the 2D horizontal anisotropic, and the 3D isotropic ranges. The existence of two successive flow anisotropisations is specific to 3D zonostrophic turbulence. Due to our resolution limit, only the zonal anisotropy range is resolved in our simulations. Note, however, that the configuration depicted in Fig. 1 is relative to our Saturn model, and might change using other planets' properties or by changing the internal model parameters.

All spectra and length scales presented in Fig. 1 are verifiable predictions of the zonostrophic regime if a sufficiently large zonostrophy index $R_β$ is achieved. Up to now, 2D-turbulent numerical studies as well as spectral analysis of Jupiter observations have reported good agreements in both spectral slope and magnitude (Del Genio et al., 2009; Galperin et al., 2010; Sukoriansky et al., 2012; Galperin et al., 2014). In the present study, we extend the statistical analysis of zonostrophic turbulence to a 3D-turbulent model. Contrary to previous 2D barotropic studies, flow confinement in our GCM is not assumed but results from the shallow-rotating conditions that enable inverse energy cascade. Hence, geometrical $n_S$ and dynamical $n_{Ro}$ indices become two additional relevant scales of the system. The index $n_S$ is unchanged in all our simulations for a fixed atmospheric layer depth $H$. Thus, varying the rotation rate defines simulations with distinct Rossby numbers, which sets the strength of the upscale energy transfer.
# 4. Statistical flow analysis

## 4.1. Kinetic energy spectra

To compute the spectral analysis of the horizontal velocity fields $\mathbf{u}$ (i.e. a wind map at a fixed vertical position) it is appropriate to decompose the vector field in terms of two scalar functions. A decomposition into rotational (non-divergent) $\mathbf{v}$ and divergent (irrotational) $\mathbf{w}$ velocity components can be performed via the Helmholtz expression $\mathbf{u} = \nabla_h \times \psi \mathbf{e}_z + \nabla_h \phi = \mathbf{v} + \mathbf{w}$, where scalar fields $\psi$ and $\phi$ are the horizontal streamfunction and velocity potential, respectively. Using this decomposition, one obtains the vorticity $\xi$ and the divergence $\delta$ fields,

$$ \xi = \mathrm{rot}_h(\mathbf{u}) = \nabla_h^2 \psi \quad (9a) $$

$$ \delta = \mathrm{div}_h(\mathbf{u}) = \nabla_h^2 \phi. \quad (9b) $$

By definition, the scalar product of two horizontal vector fields $\mathbf{a}$ and $\mathbf{b}$ can be expressed in the spectral domain as,

$$ (\mathbf{a}, \mathbf{b})_{m,n} = \frac{a^2}{n(n+1)} \mathbb{R} \{ \mathrm{rot}_h(\mathbf{a})_{m,n} \mathrm{rot}_h^*(\mathbf{b})_{m,n} + \mathrm{div}_h(\mathbf{a})_{m,n} \mathrm{div}_h^*(\mathbf{b})_{m,n} \} \quad (10) $$

where $*$ is the complex conjugate and $\mathbb{R}\{\cdot\}$ refers to the real part (see Boer, 1983; Augier and Lindborg, 2013). It is clear from expression (10) that the 2D kinetic energy spectrum (i.e. spectral energy density per modes $m, n$) is

$$ \mathcal{E}(n, m) = \frac{(\mathbf{u}, \mathbf{u})_{m,n}}{2} = \frac{a^2}{2n(n+1)} \mathbb{R} \{ \xi_{m,n} \xi_{m,n}^* + \delta_{m,n} \delta_{m,n}^* \}, \quad (11) $$

where $\xi_{m,n}$ and $\delta_{m,n}$ are the spherical harmonic coefficients of the vorticity and divergence fields. Following Sukoriansky et al. (2002), we define the 1D energy spectrum as $E(n) = \sum_{m=-n}^{n} \mathcal{E}(n, m)$, which can be decomposed into a sum of the zonal and residual components, $E(n) = E_Z(n) + E_R(n)$. The zonal spectrum corresponds to the axisymmetric energy ($m=0$):

$$ E_Z(n) = \mathcal{E}(n, 0) \quad (12) $$

while the residual spectrum is a sum over all non-axisymmetric modes ($m \neq 0$):

$$ E_R(n) = 2 \sum_{m=1}^{n} \mathcal{E}(n, m). \quad (13) $$

In the text below, we will refer to the modal spectra $\mathcal{E}(n, m)$ as it corresponds to kinetic energy for each pair of modes ($n, m$).
### 4.2. Spectral fluxes

In order to better characterize the non-linear dynamics of eddy-eddy interactions in Saturn atmosphere, we compute non-linear enstrophy and energy fluxes, which provide further insight into the enstrophy and energy redistribution among scales. Starting with the vorticity equation,

$$ \frac{\partial \xi}{\partial t} = -(\mathbf{v} \cdot \nabla)\xi - D, \qquad (14) $$

where $\mathbf{v}$ is the rotational velocity, and $D$ is the vorticity tendency due to divergent flow and other vorticity sources and sinks. In spectral space, the enstrophy equation is

$$ \frac{\partial G_n}{\partial t} = J_n + D_n^G, \qquad (15) $$

with enstrophy $G = \frac{1}{2}\xi^2$, and divergent effects, sources, and sinks contained within $D_n^G$. The enstrophy interaction term is given by the enstrophy transfer function (Burgess et al., 2013; Read et al., 2018),

$$ J_n = -\frac{1}{4} \sum_{m=-n}^{n} \left[ \xi_{m,n}^* \{\mathbf{v} \cdot \nabla \xi\}_{m,n} + \xi_{m,n} \{\mathbf{v} \cdot \nabla \xi\}_{m,n}^* \right], \qquad (16) $$

and the associated energy transfer function is

$$ I_n = \frac{a^2}{n(n+1)} J_n. \qquad (17) $$

The enstrophy and energy fluxes can now be computed:

$$ \Pi_{n+1}^Z = -\sum_{l=1}^{n} J_l \quad \text{and} \quad \Pi_{n+1}^E = -\sum_{l=1}^{n} I_l \qquad (18) $$

respectively. The integral over all wavenumbers $N$ is zero by construction and should also be equal to zero at $n = 0$ as the energy transfer functions represent conservative processes.

## 5. 1/2° Saturn reference simulation

### 5.1. Time evolution of the kinetic energy

Spiga et al. (2020) report on simulations covering 15 Saturn years using the Saturn DYNAMICO GCM. Wind fields are output every 20 Saturn days
<figure><img src="image_3.png" /></figure>

Figure 2: Time evolution of kinetic energy as a function of Saturn's simulated years: time evolution of the zonal ($E_Z$ in red), residual ($E_R$ in black), and total energy ($E_T$ in blue), in Saturn years. Energy is in $m^2 s^{-2}$. Yellow stars correspond to times $t \approx 0.002$ years (60 days), 0.009 (220 days), 0.012 (300 days), 0.015 (380 days), 0.02 (540 days), 0.04, 0.1, 0.2, and 0.6 years. The black rectangle zooms in to the twelfth year, where the yellow stars are at times $t2 = 12.5$, $t3 = 12.78$, and $t4 = 12.8$ years.

at 32 pressure levels onto 1/2° latitude-longitude grid maps. The time evolution of the total energy $E_T = E_Z + E_R$, as well as the zonal and residual energy, $E_Z = \sum_{n=0}^N E_Z(n)$ and $E_R = \sum_{n=0}^N E_R(n)$ respectively, are reported in Fig. 2, with maximum total index $N = 360$ for 1/2° horizontal resolution. The total energy strongly increases to reach a plateau after 1 Saturn simulated year. In the steady state, the zonal energy is an order of magnitude larger than the residual energy, and remains approximately constant over the 15 simulated years. The ratios of the time-average zonal and residual energy to the total energy over the two last simulated years are $\hat{E}_Z/\hat{E}_T \approx 0.92$ and $\hat{E}_R/\hat{E}_T \approx 0.08$, showing strong flow anisotropy in the zonal direction. The energy profiles fluctuate around a mean value with two dominant frequencies corresponding to annual and semi-annual periods. A typical annual period is shown in the black square of Fig. 2, which zooms into the twelfth simulated Saturn year. In what follows, we refer to these annual and semi-annual periods as seasonal energy cycles.
<figure><img src="image_4.png" /><figcaption>Figure 3: Instantaneous (latitude-longitude) velocity and vorticity maps: zonal velocity (ms⁻¹) (left), meridional velocity (ms⁻¹) (middle) and vorticity fields (s⁻¹) (right), at altitude 94.28 km, at times t1 = 0.04, t2 = 12.5, t3 = 12.78, and t4 = 12.8 years. These particular times are indicated with yellow stars in Fig. 2.</figcaption></figure>

5.2. *Horizontal velocity and vorticity fields*

Fig. 3 shows instantaneous zonal and meridional velocity fields as well as the vorticity field at time step t1 = 0.04 Saturn year and once steady state is achieved. All wind fields are shown at 94.28 km above the 3 bar level (located on vertical cross-section Fig. 4), roughly corresponding to Saturn's visible tropospheric cloud deck. After time t1 = 0.04 Saturn years, baroclinic instabilities drive zonal and meridional motions with comparable magnitude, preferentially in the equatorial region.

In the steady state, the zonal velocity component largely dominates the flow, forming multiple axisymmetric jets with mid-latitude jets that reach 60–70 m s⁻¹. The number of jets is consistent with the most recent observations
(i.e. three per hemisphere) (Porco et al., 2005; García-Melendo et al., 2010). However, our model underestimates the jet magnitude by a factor 1.5–2. The prograde equatorial jet is only weakly super-rotating, and its intensity is approximately one order of magnitude smaller than observations. This means that jet acceleration by waves and eddies, albeit significant, remains underestimated with our model’s current configuration (Spiga et al., 2020). The three lower maps in Fig. 3 correspond to a sequence of increasing residual energy in the annual cycle observed in Fig. 2. We show times $t2 = 12.5$, $t3 = 12.78$, and $t4 = 12.8$ years (see the yellow stars Fig. 2), which correspond to the minimum, intermediate, and maximum of residual energy in the zoomed in part of Fig. 2, respectively. This sequence begins at $t2 = 12.5$ years with a pattern of steady jets, minimizing the non-axisymmetric perturbations. At $t3 = 12.78$ years, an instability grows from a jet and propagates throughout the entire southern hemisphere (referred to by Spiga et al. (2020) as “eddy bursts” following the early work by Panetta (1993)). The system finally reaches a maximum of residual energy near $t4 = 12.8$ years, before it returns to a steady jet pattern with low residual energy. The same sequence occurs in the northern and southern hemispheres within the annual cycle. These energy cycles are due to the intermittent forcing that relates to the annual and seasonal variations of sunlight absorption in the atmosphere.

### 5.3. Meridional cross-section of the velocity field

Figure 4 shows meridional cross-sections through the zonal and meridional velocities at time $t1 = 0.04$ and $t4 = 12.5$ Saturn years. At the earlier time, baroclinic instabilities force zonal and meridional motions in the topmost layers, while the bottom velocity is close to zero. A strong velocity gradient prevails in the vertical and flow barotropization (i.e. vertical invariance) is only perceptible in the equatorial region. Nonetheless, the presence of a velocity gradient during the early simulated days provides evidence for vertical anisotropy in the forcing. Once the steady state is achieved at the later time shown, the flow is dynamically constrained to have a weak vertical gradient at all latitudes (i.e. it becomes strongly barotropic). The jets extend throughout the domain depth and form a quasi-2D dynamical system in which the barotropic component of the flow accumulates most of the kinetic energy in Saturn’s global circulation. As suggested by Charney (1971), flow barotropization can be related to an upscale energy cascade into the vertical flow component. Energy accumulation into the barotropic mode reaches its maximum in the troposphere and lower stratosphere (depicted as a dashed
black square in Fig. 4) where the barotropic energy exceeds its baroclinic
counterpart by a factor 7 while it decreases to a factor 3 in the upper strato-
sphere. This might relate to the continuous energetic forcing by baroclinic
instabilities that occurs in the topmost stratospheric region. We also note
that barotropization appears to be mainly active at high- and mid-latitude
consistently with what have been observed by Chemke and Kaspi (2015) for
the Earth's atmosphere. Note that due to the shallow-atmosphere approxi-
mation, where the latitudinal projection of the rotation spin axis is neglected,
these barotropic columns align with the spherical radius. By implementing
all components of the Coriolis force in the model equations, these columns
would likely tilt in latitude. Such flow barotropization is expected to extend
down to depths of thousands of kilometres beneath the tropospheric cloud
level, probably into the region of magnetic dissipation at a depth of ~ 3,000
km for Jupiter (Kaspi et al., 2018) and ~ 9,000 km for Saturn (Guillot et al.,
2018; Galanti et al., 2019). Here, our Saturn reference simulation limits us to
a thin atmosphere of ~ 360 km in the vertical, which remains a key intrinsic
limitation of shallow GCMs.

5.4. *Spectral analysis of the wind field*

At different vertical levels we computed the horizontal velocity spherical
harmonic decomposition, following Eq. (11). Figure 5 summarizes the kinetic
energy spectra in the steady state for our Saturn reference simulation, show-
ing the zonal spectrum E<sub>Z</sub>(n), the residual spectrum E<sub>R</sub>(n), and the modal
spectrum ℰ(n, m) (for modes m = 1 to 10 only). All spectra are averaged
in time over the two last simulated years, and over several vertical levels
sampling both the troposphere and the lower stratosphere, i.e. where energy
strongly accumulates into the barotropic mode (see black dashed square in
Fig. 4).

The zonal spectrum follows the -5 slope and magnitude predicted by
Eq. (4a). At large scales, the zonal spectrum exceeds the residual spectrum
by about two orders of magnitude, providing evidence for the strong zonal
anisotropy inherent to zonostrophic turbulence. Although m = 1, 2, and 3
are more energetic and bumpy than other modes, the modal spectra preserve
a nearly n<sup>-8/3</sup> slope over a wide range of n. Likewise, the global trend in
the residual spectrum over the full range of n fits the n<sup>-5/3</sup> slope expected
from the KBK theory of 2D turbulence, and makes it possible to estimate the
inverse energy cascade rate ϵ. The bump in the residual and modal spectra
<figure><img src="image_5.png" /><figcaption>Figure 4: Instantaneous (altitude-latitude) meridional cross-section of velocity maps: zonal (left) and meridional (right) velocities (m s⁻¹) at longitude λ = 0°, times t2 = 0.04 (top) and t4 = 12.5 (bottom) Saturn years, as a function of latitude and altitude (in km, with 0 km at the 3-bar level, which is the bottom of the model). In the bottom left panel, the solid line corresponds to the level at which horizontal maps are shown. The black dashed square shows the vertical depth over which the statistical analyses are averaged.</figcaption></figure>

at moderate *n*, which makes these spectra differ at medium scales from the
*n*⁻⁵/³ and *n*⁻⁸/³ slopes respectively, will be discussed below. Using Eq. (4b)
and assuming *C*<sub>*R*</sub> = 6, we estimate *ϵ* ≈ 2 × 10⁻⁶ m² s⁻³. Considering that the
total kinetic energy growth during model spinup will be approximately linear,
i.e. *E*<sub>*T*</sub>(*t*) = *ϵ*t, and that the spinup period lasts about one Saturn simulated
year or 10⁹ s (Fig. 2), the equilibrium total kinetic energy of about 2000 m²
s⁻² in Fig. 2 corresponds very well to the total energy injected during the
first year of the simulation. With *β* = 4 × 10⁻¹² m⁻¹ s⁻¹ and *C*<sub>*Z*</sub> = 0.5, we
compute the transitional index *n*<sub>β</sub> ≃ 58, the zonal index *n*<sub>*Z*</sub> ≃ 449, and the
Rossby index *n*<sub>*Ro*</sub> ≃ 2.7 × 10⁵ and non-dimensional number *Ro* = 4 × 10⁻³ for
our Saturn reference simulation. We estimate the RMS velocity *U* ≃ 55 m s⁻¹
using the total energy in the steady state *E*̂<sub>*T*</sub>, leading to a value *n*<sub>*R*</sub> ≃ 11 for
<figure><img src="image_6.png" /><figcaption>Figure 5: Spatial kinetic energy spectra as a function of latitudinal index: vertical grey lines are the Rhines scale $n_R$ (left) and the transitional scale $n_β$ (right). The grey area corresponds to the Rossby deformation radii $n_D$, showing typical scales of the baroclinic instabilities. Normalized zonal and residual energy spectra $E_Z(n)/(\epsilon a)^{2/3}$ (red) and $E_R(n)/(\epsilon a)^{2/3}$ (black) as well as modal spectra for $m = 1$ to 10, $\mathcal{E}(n, m)/(\epsilon a)^{2/3}$ (green). Dashed lines represent theoretical predictions from expressions: Eqs (4a) in red, (4b) in black and Eq (7) in green.</figcaption></figure>

the Rhines index. In Fig. 5, the Rhines and transitional indices border the zonostrophic inertial range characterized by its zonostrophy index $R_β = 5.3$, which significantly exceeds the zonostrophic threshold $R_β \ge 2.5$ stated in Section 3.2. These statistics are summarized in Table 1. The estimated value of $n_Z$ exceeds the largest index $N = 360$ resolved in our 1/2° Saturn reference simulation, and so does not appear in Fig. 5.

By approximating the zonal and residual spectra with their predictions, Eqs. (4a) and (4b), assuming straight horizontal lines for $n \le n_{fr}$, and by integrating over indices $n$, one finds analytical expressions for $E_Z$ and $E_R$ (Sukoriansky et al., 2002; Galperin et al., 2014). This leads, in the framework
<table><caption>Table 1: Summary of diagnostic quantities that facilitates comparison between our Saturn reference simulation (SRS) and previous studies of Jupiter's tropospheric dynamics using cloud tracking from Cassini images.</caption><thead><tr><th></th><th>ε (m<sup>2</sup>s<sup>-3</sup>)</th><th>n<sub>Z</sub></th><th>n<sub>β</sub></th><th>n<sub>R</sub></th><th>R<sub>β</sub></th></tr></thead><tbody><tr><td>Saturn reference simulation</td><td>2 × 10<sup>-6</sup></td><td>449</td><td>58</td><td>11</td><td>5.3</td></tr><tr><td>Choi and Showman (2011)</td><td>5 × 10<sup>-6</sup></td><td>550</td><td>67</td><td>10.9</td><td>6.1</td></tr><tr><td>Galperin et al. (2014)</td><td>10<sup>-5</sup></td><td>430</td><td>57</td><td>10.9</td><td>5.2</td></tr></tbody></table>

of zonostrophic turbulence, to the definition of the anisotropy index

$$ \gamma = E_Z/E_R \approx (3/4)(10C_Z)^{-5/6} R_{\beta}^{10/3}. \qquad (19) $$

The associated ratios $E_Z/E_T = \gamma/(1+\gamma)$ and $\hat{E}_R/\hat{E}_T = 1/(1+\gamma)$ give $E_Z/E_T \approx 0.98$ and $E_R/E_T \approx 0.02$ with $R_\beta = 5.3$. This slightly overestimates the time-averaged ratios obtained in Section 5.1 (i.e. $E_Z/E_T \approx 0.92$ and $E_R/E_T \approx 0.8$) from spectral analysis of the simulated flow. Both estimates nevertheless point towards significant flow anisotropy.

The statistical analysis of our Saturn simulation is consistent with all basic properties of the theoretical scaling introduced in Section 3 and summarised in Fig. 1. The non-dimensional scaling used in this study allows for direct comparison with other planetary atmospheres in which zonostrophic turbulence develops. Statistical analyses of Jupiter's atmospheric circulation at cloud level have been previously described by Choi and Showman (2011) and Galperin et al. (2014) using Jupiter images from the Cassini mission. Galperin et al. (2014) calculated non-dimensional scalings similar to those used in the present study based on both their own wind field data set and Choi and Showman (2011)'s. The resulting non-dimensional quantities are reported in Table 1 and show that the Jovian troposphere and our Saturn numerical simulation share analogous dynamics with similar zonostrophy indices and spectral anisotropy. However, specific features of Saturn's dynamics appear when looking at the details of the energy spectra. In our simulation, the modal spectra $\mathcal{E}(n, m)$ show that modes $m=2$ and $3$ make a dominant contribution to the residual spectrum, suggesting strongly energetic waves that propagate in the zonal direction. Spiga et al. (2020) identified indeed wavenumber-2 and wavenumber-3 Rossby waves and Rossby-gravity waves, with possible relation to wave patterns observed by either Voyager (Achter-
berg and Flasar, 1996) or Cassini (Guerlet et al., 2018). However, as of now, a complete statistical analysis using Saturn direct observations is yet to be performed on the most recent images. Consequently, prominent modes and spectral energy bumps in our Saturn reference simulation remain to be confirmed by observations.

Our simulation captures the global statistical picture of the gas giants and specifically Saturn’s dynamics in some respects, despite the *ad hoc* conditions enforced in our GCM to compensate for unknown frictional effects and unresolved small scale dissipation. Hypeviscosity is responsible for the sudden energy drops where all spectra significantly deviate from predictions. It qualitatively defines a typical range of indices, $n_{\delta} \gtrsim 170$, for which hypeviscosity effects are predominant. Frictional effects are less evident but a likely crucial to set the large-scale dynamical regime in which a wealth of turbulent waves and eddies develop. Our model also suffers intrinsic limitations such as limited vertical resolution and the absence of moist convection, which are both likely to impact wave propagation and global atmospheric circulation.

### 5.5. *Spectral fluxes and time evolution*

We present in Fig. 6 transient dynamics showing instantaneous kinetic energy spectra as well as energy and enstrophy spectral fluxes. Time evolution of the spectral quantities are shown in Fig. 6(a-c) from an atmosphere at rest to the steady state, and in (d-f) we show the annual oscillation in energy diagnostics, which can be interpreted as the standard deviation of the mean state. These instantaneous spectra correspond to the time steps marked with yellow stars in Fig. 2. In what follows, the terms “steady state” and “mean state” are used synonymously to refer to time-averaged spectral quantities using the averaging procedure described in Sect. 5.4 (and Fig. 5). This section provides more details on non-linear interactions and general energetic pathways in Saturn’s atmospheric circulation throughout our numerical simulation. More specifically, this approach allows us to identify the scale at which kinetic energy is created in the system, how it transits, and where it ends up.

Fig. 6(a) shows the time evolution of the residual kinetic energy spectra from 60 simulated days to half a year before it reaches the steady state. Time steps are located in time on the time evolution profiles of kinetic energy in
Fig. 2. At the first time step, the spectrum is weakly energetic and approximately flat over all scales. When the flow is energized by solar absorption in the stratosphere, the spectra show that kinetic energy is predominantly increased at small scales. As mentioned above, the most relevant length scales for energy creation are the baroclinic Rossby deformation radii. However, hyperviscosity presumably damps kinetic energy in the scale range that coincides with the typical length scales of baroclinic eddies. This damping is likely to impede the formation of a sharp energy peak near $n_D$ and certainly reduce the global energy budget at all scales. Such a condition at small scales suggests that our simulation is eddy-permitting but not fully eddy-resolving. Under confinement effects, i.e. the coupled action of geometrical and dynamical confinement, the energy spectra develop a $n^{-5/3}$ slope with a bumpy front that propagates in time toward large scales, up to $n \simeq 2$. This upscale energy transfer is not arrested at the Rhines scale, and extends over a large spectral domain, i.e. $2 \lesssim n \lesssim 170$, for which the $n^{-5/3}$ slope fits the residual spectrum very well in the steady state.

Figure 6(b,c) shows snapshots and the steady state of the spectral fluxes at the same time steps. Positive spectral flux corresponds to energy transfer from small to large indices $n$, which is a downscale energy cascade, and conversely negative values correspond to an upscale energy cascade. The general trend in the steady state is clear: energy flux is negative from $n_\delta = 170$ down to $n = 2$, showing a wide inertial range with upscale energy transfer. This is consistent with the $n^{-5/3}$ spectral signature of residual kinetic energy shown in Fig. 6(a). At scales $n_\delta \gtrsim 170$, the positive energy flux suggests a downscale energy cascade over a range that perfectly matches the spectral domain where hyperviscosity directly affects the kinetic energy spectra. The index $n_\delta = 170$ at which energy flux changes sign appears to be robust at all times during spin-up, and occurs at a scale slightly larger than the first Rossby deformation radius with typical index $n_D = 225$. hyperviscosity possibly distorts energy creation, shifting to larger scale the onset of the upscale cascade. When the flow satisfies the conditions for quasi-geostrophy, theoretical studies (e.g. Charney, 1971; Salmon, 1980) also suggest the existence of an downscale enstrophy cascade. This trend seems to be consistent with our Saturn reference simulation: in Fig. 6(c), the enstrophy flux is predominantly positive and increases strongly with $n$ up to the resolution limit. At scales larger than the transitional index $n_\beta$, the enstrophy flux becomes weakly negative, showing a tendency for an upscale enstrophy cascade. Such a feature
<figure><img src="image_7.png" /><figcaption>Figure 6: Time evolution of spatial kinetic energy spectra as well as energy and enstrophy fluxes as a function of latitudinal index: vertical grey lines are the Rhines scale $n_R$ (left) and the transitional scale $n_β$ (right). The grey area corresponds to the Rossby deformation radii $n_D$, showing typical scales of the baroclinic instabilities. Normalized (top) residual kinetic energy spectra $E_R(n)/(\epsilon a)^{2/3}$, (middle) spectral fluxes of energy $\Pi^E/\epsilon$ and (bottom) spectral fluxes of enstrophy $\Pi^Z/(\epsilon^{-1}a^2)$. In all panels, spectral quantities in the steady state are averaged in time and shown as solid black lines. (Left) Transients of the instantaneous spectral quantities from the start of the simulation (colour gradient) to the steady state (black line). The colour gradient from light to dark lines corresponds to times $t \simeq 0.002$ (60 days), $0.009$ (220 days), $0.012$ (300 days), $0.015$ (380 days), $0.02$ (540 days), $0.04$, $0.1$, $0.2$, and $0.6$ years. (Right) Instantaneous spectral quantities near the minimum of the energy cycle at $t2 = 12.5$ years (coloured dashed lines) and the maximum of the energy cycle at $t4 = 12.8$ years (coloured solid lines) of the annual sequence shown in Fig. 2.</figcaption></figure>
has also been observed for Jupiter by Young and Read (2017), in which the upscale enstrophy transfer induces a weak peak in enstrophy near the Rhines scale.

In Fig. 6(d-f), we bracket the mean state with two instantaneous plots of spectral quantities close to the minimum (dashed lines) and the maximum (thin solid lines) of residual energy in the annual cycle described in Sect. 5.1 and presented in Fig. 2. Spectral instantaneous diagnostics are computed at times $t2 = 12.5$ and $t4 = 12.8$ years. The spectra in Fig. 6(d) reveal interesting energy exchanges once the steady state is achieved: the Rhines scale splits the dynamics between long-lived modes, which control the large scales for indices $n < n_R$, and transitional modes, which develop and break up within the zonostrophic inertial range (i.e. $n_R < n < n_β$) during the annual energy cycle.

On the one hand, in the large scale range $n < n_R$, the long-lived structures are dominated by two highly energetic zonal modes, $m = 2$ and $m = 3$, whose amplitude is insensitive to the annual variations. Thus, we refer to those as long-lived modes: non-stationary $m = 2$ and $m = 3$ inertial waves which are time-independent when looking at spatial spectra $E_R$. This is illustrated in Fig. 6(d), where instantaneous spectra are similar to the spectrum of the mean state in the range $n < n_R$. This is also observed in physical space when looking at meridional velocity maps (Fig. 3), for which the $m = 3$ mode is persistent over the whole energy sequence shown (the $m = 2$ mode was shown in Spiga et al., 2020, Fig. 13).

On the other hand, in the zonostrophic range, the residual spectrum oscillates around the mean state showing a strongly energetic phase with a nearly $n^{-3}$ slope, and a phase of low energy matching the KBK $n^{-5/3}$ slope. (Chertkov et al., 2007). The $n^{-3}$ spectral signature and the associated bumpy residual spectrum might be related to the instability of the zonal jets, which are barotropically unstable according to the Rayleigh-Kuo criterion (see Ingersoll and Pollard (1982) for the general perspective, and Spiga et al. (2020) for the SRS specifically). To support this, the modal spectra in Fig. 7, show that the strongly energetic phase, is related to an enrichment of a large variety of zonal modes that result from jet destabilization due to barotropic instability. This intermittent enrichment by zonal modes of the energy spectra is dominantly observed at scales that are close to but smaller than the Rhines scale, making the spectral slope at $n > n_R$ much steeper. We note that the zonal spectrum remains unchanged over time, and so are the jets.
Similarly, snapshots of the spectral energy flux over the annual cycle (Fig. 6e) show that the extended domain of the upscale cascade in Fig. 6(b) results from the averaging procedure over time, and does not account for instantaneous dynamics. Instantaneous diagnostics of the energy flux show strong perturbations and sign changes, particularly in the range $n < n_β$, which suggests a superimposition of direct and inverse energy cascades. Again, when jets become barotropically unstable they break into turbulence and feed a direct cascade of energy that appears to be of a comparable strength to the inverse cascade in the low energetic phase. Nonetheless, when the residual energy is at a maximum, the upscale cascade dominates except for a small scale range, i.e. $n \lesssim n_R$, where energy converges both from large and small scales as energy flux changes sign at the largest scales. Such a convergence of kinetic energy at the jet scale has been reported by Young and Read (2017), who computed spectral fluxes from snapshots of Jupiter's dynamics (covering 3-4 planetary rotation periods only). We also note that, within the zonostrophic inertial range, local perturbations of the energy flux (with a serrated shape) are associated with similar perturbations of the enstrophy flux, which departs from the classical picture of 2D turbulence (see Fig. 6 e-f). This erratic enstrophy production at scales larger than the transitional scale $n_β$, only visible in instantaneous fluxes, is suggestive of vorticity filamentation that develops in the zonostrophic inertial range of 3D-turbulent atmospheric layers.

In this section we have explored the dynamics of a 3D-turbulent atmospheric layer that reproduces the statistical properties of zonostrophic turbulence. Under geometric and dynamical confinements, time-averaged spectral quantities are in good agreement with the energy-enstrophy double cascade scenario. This double cascade found in our 3D simulations may be reminiscent of the double cascade described in 2D flows by Kraichnan (1967) theory. The statistical analysis of our Saturn reference simulation resembles in several respects the spectral quantities reported from Jupiter images (Young and Read, 2017) and other GCMs (Schneider and Liu, 2009; Read et al., 2018). However, a novelty of our statistical analysis is that by looking at the instantaneous spectral quantities we reveal peculiar relevance of the Rhines and transitional scales in the 3D-turbulent zonostrophic regime. The Rhines scale seems to act as an energetic barrier that splits the dynamics between long-lived modes, i.e. dominantly $m = 2$ and $m = 3$, and intermittent
<figure><img src="image_8.png" /><figcaption>Figure 7: Instantaneous modal kinetic energy spectra as a function of latitudinal index: normalized zonal energy spectra $E_Z(n)/(\epsilon a)^{2/3}$ (red) and modal spectra for $m = 1$ to 10, $\mathcal{E}(n, m)/(\epsilon a)^{2/3}$ (green). Left panel corresponds to the minimum of energy at $t2 = 12.5$ years and right panel to the maximum of energy at $t4 = 12.8$ years, in the annual energy cycle shown in Fig. 2. Dashed lines represent theoretical predictions of the modal spectra from Eq (7) and Rhines index $n_r$ (vertical dotted lines).</figcaption></figure>

modes, which strongly depend on the seasonal energy cycles. Instantaneous quantities also reveal the superimposition of the inverse and direct cascades of both energy and enstrophy, in the extended inertial range $n < n_β$. These complex dynamics are likely related to the formation of inertial modes, which develop and undergo vorticity filamentation in 3D-turbulent atmospheric layers driven by an intermittent forcing that reflects the seasonal energy cycles. In the next section we explore the role played by these seasonal energy cycles in the global flow statistics by suppressing radiative transfer computations in the GCM and implementing instead a continuous idealized Taylor-Green forcing that is time-independent.

# 6. Idealized simulations

In order to explore the regime of Saturn macro-turbulence and other open questions about the nature of the atmospheric circulation of giant planets, we run a set of four idealized numerical simulations with varying rotation rate and a well-controlled energy forcing.

We use an idealized version of our Global Climate Model based on DY-NAMICO in which atmospheric flow in a Saturn-like atmosphere is driven by a Taylor-Green forcing (Alexakis, 2015) and flow baroclinicity is suppressed (this is obtained by switching off the physical parameterizations described in
Guerlet et al., 2014, notably radiative transfer). Kinetic energy is injected every 32 dynamical time steps at all vertical levels at a prescribed horizontal wavenumber $k_i$. The forcing is given by

$$ \mathbf{f} = \begin{pmatrix} A \sin(\varphi) \cos(k_i \lambda) \sin(k_i \varphi) \cdot \mathbf{e}_{\varphi} \\ A \cos(\varphi) \sin(k_i \lambda) \sin(k_i \varphi) \cdot \mathbf{e}_{\lambda} \\ 0 \cdot \mathbf{e}_{z} \end{pmatrix} \qquad (20) $$

This approach avoids both forcing intermittency (i.e. semi-annual and annual cycles, as well as eddy bursts) and spatial anisotropy. In addition, the energy injection rate is known analytically from Eq. (10),

$$ \epsilon_i = \sum_{n=0}^{N} \sum_{m=-n}^{n} (\mathbf{u}, \mathbf{f})_{m,n}. \qquad (21) $$

The time-averaged injection rate, defined as $\hat{\epsilon}_i = \langle \epsilon_i \rangle_T$, is set to emulate the energy transfer rate $\epsilon$ of our Saturn reference simulation. The forcing magnitude $A = 0.22 \text{ m s}^{-1}$ leads to a mean energy injection rate that ranges between $5.5 \times 10^{-6} \le \hat{\epsilon}_i \le 5.6 \times 10^{-6}$ for all simulations. The injection wavenumber $k_f = 56$ project on index $n_i = 112$ when the basis functions of spherical harmonics are used. It is chosen to be both small-scale and outside the range where model hyperviscosity affects the flow statistics (i.e. $n_i < n_\delta$). One relatively unconstrained aspect of the dynamics is the strength, and latitudinal extent, of the artificial bottom drag due to the lack of detailed knowledge about exactly where and how the MHD drag acts on the atmospheric circulation. Here, we choose the drag formulation to be the same as for our Saturn reference simulation (Spiga et al., 2020) to ensure that differences among the simulations are not caused by differences in this poorly constrained parameter. Investigating the effects of bottom friction that govern flow statistics at large scales is beyond the scope of the present paper. Then rotation rate is the only varying parameter, which we set to $4\Omega$, $\Omega$, $\Omega/2$ and $\Omega/4$, and each simulation reaches equilibrium over a timescale of one Saturn year (a steady state plateau in energy is reached by 0.2 year, not shown).

## 6.1. Horizontal velocity and vorticity fields

Fig. 8 shows zonal, meridional and vorticity fields at the last remapped time step and at the same vertical level as the horizontal maps of the Saturn
<figure><img src="image_9.png" /><figcaption>Figure 8: Instantaneous (latitude-longitude) velocity and vorticity maps: flow fields are shown for four idealised simulations with varying planetary rotation rate. (Left) zonal velocity (m s⁻¹), (middle) meridional velocity (m s⁻¹), and (right) vorticity (s⁻¹) at trophospheric depth 94.28 km above 3 bar.</figcaption></figure>

reference simulation in Fig. 3 (i.e. 94.28 km above 3 bar). The simulations are ordered by decreasing rotation rate or, equivalently, by increasing Rossby non-dimensional number. In all four idealized simulations, the combination of low Rossby number with a strong β-effect leads to velocity and vorticity patterns dominated by multiple eddy-driven zonal jets. The latitudinal jet size decreases with increasing rotation rate, in keeping with Rhines’ theory. Increasing rotation is also accompanied by the formation of intense non-axisymmetric modes, which compete with the $m = 0$ dominant mode. This is particularly evident at $\Omega$ and $4\Omega$, in which the meridional velocity and vorticity fields are considerably strengthened.
<figure><img src="image_10.png" /><figcaption>Figure 9: Spatial kinetic energy spectra as well as energy and entrophy fluxes as a function of latitudinal index: spectral quantities are reported for four idealized simulations at various rotation rates. The energy injection index $n_i$ is shown with black vertical arrows. Vertical dashed lines show the Rhines ($n_R$) and transitional ($n_β$) indices from left to right in each panel. Vertical arrows correspond to the injection index $n_i = 112$. (a) Time-averaged normalized kinetic energy distribution showing zonal energy $E_Z(n)/(\epsilon_i a)^{2/3}$ (red), residuals $E_R(n)/(\epsilon_i a)^{2/3}$ (black), and modal spectra up to $m = 10$, $\mathcal{E}(n, m)/(\epsilon_i a)^{2/3}$ (green), as well as their theoretical predictions (Eqs. 4a, 4b, and 7). (b) Time-averaged (black) and instantaneous (yellow) energy fluxes $\Pi^E/\epsilon_i$. (c) Time-averaged (black) and instantaneous (blue) entrophy fluxes $\Pi^Z/(\epsilon_i^{-1}a^2)$.</figcaption></figure>
<table><caption>Table 2: Summary of energy diagnostic quantities, for comparison between our Saturn reference simulation (SRS) and the idealized simulations.</caption><thead><tr><th></th><th>n<sub>Z</sub></th><th>n<sub>β</sub></th><th>n<sub>R</sub></th><th>R<sub>β</sub></th><th>ε (m<sup>2</sup>s<sup>-3</sup>)</th><th>n<sub>Ro</sub></th><th>Ro</th><th>E<sub>rot</sub>/E<sub>div</sub></th></tr></thead><tbody><tr><td>SRS</td><td>449</td><td>58</td><td>11</td><td>5.3</td><td>2 × 10<sup>-6</sup></td><td>2.7 × 10<sup>5</sup></td><td>4 × 10<sup>-3</sup></td><td>200</td></tr><tr><td>4Ω</td><td>1238</td><td>118</td><td>31</td><td>3.9</td><td>6 × 10<sup>-6</sup></td><td>2 × 10<sup>6</sup></td><td>5 × 10<sup>-4</sup></td><td>100</td></tr><tr><td>Ω</td><td>773</td><td>85</td><td>13</td><td>6.3</td><td>3 × 10<sup>-7</sup></td><td>7 × 10<sup>5</sup></td><td>3 × 10<sup>-3</sup></td><td>1000</td></tr><tr><td>Ω/2</td><td>427</td><td>56</td><td>9.5</td><td>5.9</td><td>3 × 10<sup>-7</sup></td><td>2 × 10<sup>5</sup></td><td>5 × 10<sup>-3</sup></td><td>1000</td></tr><tr><td>Ω/4</td><td>235</td><td>37</td><td>7.3</td><td>5.1</td><td>3 × 10<sup>-7</sup></td><td>8 × 10<sup>4</sup></td><td>9 × 10<sup>-3</sup></td><td>200</td></tr></tbody></table>

## 6.2. Energy spectra of the wind fields

By computing spherical harmonic decompositions of the wind fields, we show in Fig. 9 the time-averaged residual, zonal and modal spectra for each rotation rate. Zonostrophic turbulence is a useful framework for interpreting the statistical energy distributions from these simulations. In all simulations (except $4\Omega$ that slightly departs from predictions) the spectra display a self-similar form with a well delimited zonostrophic inertial range in which zonal energy dominates residuals with a magnitude proportional to $\beta^2$ (Eq. 4a). The Rhines index locates quite effectively the most energetic scale, above which energy systematically decreases. Table 2 lists the indices and non-dimensional numbers estimated from our spectral analysis for each simulation.

The simulation at rotation $\Omega/4$ has similar statistical properties to numerical studies that model 2D-turbulent flow on the surface of a rotating sphere (Sukoriansky et al., 2002, 2007; Galperin et al., 2010). We retrieve the two distinct spectral domains shown in Fig. 1: the zonostrophic inertial range for $n_R < n < n_\beta$, and the residual-dominated range for $n_\beta < n < n_i$. The residual spectrum follows the KBK $n^{-5/3}$ slope, whose magnitude at the injection scale is estimated to be $\epsilon \simeq 3 \times 10^{-7} \text{ m}^2 \text{s}^{-3}$. No energy is injected at the jet scale, ensuring a smooth zonal spectrum and a peaked residual spectrum at injection. The “zigzag” features of the zonal spectrum at small indices $n$ are due to the hemispheric symmetry of the zonal structure, which projects onto the even indices. The modal spectrum shows equally distributed modes apart from the axisymmetric mode $m = 0$, which dominates the dynamics with a magnitude equal to $C_Z\beta^2$ (see Eq. 4a). Finally, all modes intersect at $n_z$ and the zonostrophy index, the Rossby index and the associated Rossby number are respectively $R_\beta = 5.1$, $n_{Ro} = 8 \times 10^4$, and $Ro = 9 \times 10^{-3}$.

Considering simulation $\Omega/4$ as a benchmark that reproduces the spectral
features of zonostrophic turbulence, we now explore higher rotation rates to
frame our Saturn reference simulation within an idealized 3D-turbulent set-
up. At $\Omega/2$, the zonostrophy index and Rossby index increase to $R_\beta = 5.9$
and $n_{Ro} = 2 \times 10^5$ ($Ro = 9 \times 10^{-3}$) respectively. These non-dimensional
quantities are closest among the idealized simulations to those obtained for
the Saturn reference simulation (see Table 2). Indeed, the lower rotation rate
$\Omega/2$ is balanced by the lowest kinetic energy injection rate $\epsilon$, leading to a
similar Rossby number and strength of the zonostrophic regime. Similarly
to SRS, the residual spectrum features a $n^{-5/3}$ slope that exists between
the injection scale and the transitional scale (i.e. in the residual-dominated
range $n_\beta < n < n_i$) and a $n^{-3}$ slope within the zonostrophic inertial range
that can be related to the growth of non-axisymmetric modes in the same
fashion as for SRS Fig. 6. Also, the steepening to a $n^{-3}$ slope in the residual
spectrum might stem from the large value of the friction index ($n_{fr} \sim n_R$
in the zonostrophic regime, see Sukoriansky et al., 2007) leading to erosion
and eventual destruction of KBK statistical laws at large scale, as shown
in Sukoriansky et al. (1999). If such a $n^{-3}$ spectral slope is not retrieved in
direct observations of Saturn's flow ($n^{-3}$ slope being a common feature of the
SRS and its analogue ideal simulation at $\Omega/2$), then it indicates that bottom
friction might have to be implemented in a more elaborate way in our GCM.

Moving to rotation rate $\Omega$, the zonostrophic and Rossby indices rise to
$R_\beta = 6.3$ and $n_{Ro} = 7 \times 10^5$ ($Ro = 3 \times 10^{-3}$) respectively, and the dynamical
confinement by the jets is reinforced. Again, the residual energy spectrum is
substantially steepened and the modal spectrum now features energetic zonal
modes that pile up within the zonostrophic inertial range (predominantly
with $m=1$ and $m=4$). These modes, easily discernible in the velocity maps
(Fig. 8), break the flow's zonal axisymmetry, leading to a dynamical regime
dominated by strong zonal waves. Neither the classical KBK energy cascade
nor the zonostrophic regime of Sukoriansky et al. (2002) are fully applicable
to this flow regime. By increasing the rotation rate again to $4\Omega$ we suppress
the residual-dominated inertial range ($n_\beta < n < n_i$), the transitional index
$n_\beta$ nearly coincides with the injection scale and the resulting spectra are
curved. The lack of scale separation forces energy accumulation within a
restricted spectral range bordered by the Rhines and injection scales. The
associated low zonostrophy index $R_\beta = 3.9$ suggests that the dynamics at $4\Omega$
are dominated by frictional effects that arrest the inverse cascade at relatively
small scale compare to other simulations. Consequently, all modes have a
spectrum steeper than $n^{-5}$, and the regime of zonostrophic turbulence is therefore not applicable to this simulation.

The tendency of the residual spectrum to develop a maximum at small total wavenumber, with a steeper slope at the adjacent smaller scales ($\approx n^{-3}$), is shared by all simulations with rotation rates larger than $\Omega/4$, including the SRS. It has been suggested that the bumpy residual spectrum ($\approx n^{-3}$ slope) can be related to the breaking of large scale zonal jets that are barotropically unstable and produce non-axisymmetric zonal modes. Thus, the presence of these intense zonal modes, inducing a $\approx n^{-3}$ slope, is a consequence of the dynamical regime rather than a product of the intermittent-anisotropic forcing specific to the Saturn reference simulation. One possible argument is that non-linear wave interactions play an important role in 3D turbulent models (see, e.g., Yarom and Sharon, 2014) and induce the formation of strong zonal modes that steepen the spectral signature. This was indeed not observed in previous 2D models where depth invariance was assumed and not produced spontaneously by the dynamics. Another possible explanation is that bottom friction distorts the statistics of the inverse cascade at large scale (Sukoriansky et al., 1999), but testing this hypothesis is beyond the scope of this work. In the large-scale range $n < n_R$, the dominant long-lived modes $m=2$ and $m=3$ that are observed in the Saturn reference simulation are absent in all idealized cases. As a consequence, they are presumably related to the baroclinic conditions specific to the SRS.

Finally, we note that rotation rate only weakly impacts the $n^{-5}$ spectral signature of the zonal spectrum for all simulations (except at $4\Omega$, where the zonostrophic regime is not applicable). This is also true for residual spectra at scales smaller than the injection scale, where the spectral slope is invariably close to $n^{-3}$. This is consistent with the associated downscale enstrophy cascade at indices $n > n_i$, but it could instead, and more logically, result from the hypeviscosity damping that prevails at small scales.

### 6.3. *Spectral fluxes*

Figure 9(b,c) shows time-averaged spectral fluxes as well as instantaneous spectral fluxes randomly chosen from within the steady state period. In our benchmark simulation $\Omega/4$, we retrieve the double cascade scenario of 2D-turbulent flows, characterized by a downscale cascade of enstrophy at small scales ($n \gtrsim n_i$) and an upscale energy cascade forming an extended negative plateau at large scale ($n \lesssim n_i$). Looking at the instantaneous spectral quantities (colour gradient), we find strong deviations from the mean
fluxes. If the enstrophy flux is only weakly varying around its time average, the energy flux displays strong variations, specifically intensified beyond the transitional index, i.e. for $n \lesssim n_β$. Again, the zonostrophic and residual-dominated subranges feature different dynamics. The residual-dominated range (i.e. $n_β \lesssim n \lesssim n_i$) is characterized by a stably negative energy flux that is the counterpart of the entrenched $n^{-5/3}$ slope in the residual spectrum.

On the other hand, the zonostrophic inertial range is marked by a weak energy bump in the residual spectrum, and the associated energy fluxes strongly oscillate around the mean state, sometimes becoming positive. This implies intermittent periods of upscale and downscale energy transfer, which occurs only beyond the transitional index $n_β$. Strongly varying in space and time, these direct and inverse energy cascades are not visible when looking only at averaged spectra. Moreover, as they are present in both our idealized simulations and in the SRS, the coexistence of positive and negative energy transfers at large scale does not appear to be related to the nature of the energy forcing.

Higher rotation rates show the same global tendency except that the lack of scale separation considerably reduces the residual-dominated range ($n_β \lesssim n \lesssim n_i$) within which the $n^{-5/3}$ slope develops. Due to the increasing efficiency of dynamical confinement in latitude at larger Rossby index (and hence lower Rossby number), the upscale energy flux increases by an order of magnitude between $Ω/4$ and $Ω$. Note that energy fluxes are all normalized by the energy injection rate $\Pi^E/\hat{\epsilon}_i$ so the flux shows the proportion of energy transferred with respect to the energy that was injected. Deviations from the mean fluxes also increase with increasing rotation rate, and show a succession of strengthened upscale and downscale energy transfers. The mean energy flux remains upscale, however, as it is necessary to constantly sustain the large-scale jets. Enstrophy is also strongly impacted by the intense activity of non-axisymmetric modes and, at each downscale energy peak, corresponds to a direct enstrophy cascade. It is also interesting that the $n^{-3}$ slope in the energetic spectra has positive enstrophy flux (on average) within the zonostrophic inertial range. This spectral signature ($n^{-3}$) of a downscale enstrophy cascade can be related to energy accumulation and enstrophy production at large scales. In summary, for all simulations, the amplitude of the residual spectrum at large scales ($n_R < n < n_β$) increases with rotation rate and also with the magnitudes of both the direct enstrophy cascade and the inverse energy cascade, featuring a bumpy spectrum with a nearly $n^{-3}$ slope.
7. Conclusion

We have presented a phenomenological study of dynamical energy and vorticity transfers between different horizontal scales in a 3D-turbulent Global Climate Model. This constitutes a coherent and advanced statistical analysis of the planetary flow implemented by Spiga et al. (2020) using the new dynamical core DYNAMICO, with the unprecedented combination of fine horizontal and vertical resolution, and the incorporation of atmospheric physical processes such as radiative transfer (Guerlet et al., 2014), emulating Saturn's troposphere and stratosphere.

We found spectral evidence that an inverse energy cascade occurs under rapid rotation and shallow-depth confinement in a 3D-turbulent atmospheric layer. We recover several statistical properties observed in more classical regimes of forced rotating turbulence (specifically the inverse energy cascade due to flow barotropization, Pouquet and Marino (2013); Deusebio et al. (2014); Yarom and Sharon (2014); and the dual cascade in forced rotating turbulence was analyzed by Sukoriansky and Galperin (2016) analytically) and we show some statistics specific to Saturn's atmospheric dynamics. The growth of intense baroclinic eddies at small planetary scales and the coupled action of low Rossby and a strong β-effect induces a zonal spectral anisotropy that implies the presence of a zonostrophic turbulent regime (Sukoriansky et al., 2002). A set of zonal jets emerges as the dominant flow structure, and the flow's statistical properties follow the double cascade scenario reminiscent of the 2D turbulence paradigm (Kraichnan, 1967). Scaling analysis reveals a turbulent regime in our Saturn reference simulation that is similar to that found in Jupiter's atmosphere (Galperin et al., 2014; Young and Read, 2017). Energy accumulation into the barotropic mode is also a feature shared by our Saturn simulation and Jupiter as shown in recent observations from Juno mission (Kaspi et al., 2018) and in a simulation of a Jupiter-like atmosphere by Liu and Schneider (2015). Non-axisymmetric modes with $m = 2$ and $m = 3$, which are specific to Saturn's atmospheric dynamics, are also retrieved in our spectral analysis (Spiga et al., 2020). However, although GCM simulations have advanced a lot since the pioneering shallow models of Cho and Polvani (1996), it is clear that shallow models still lack a naturally arising prograde equatorial jet, whereas they arise generically in deep convection models (Heimpel et al., 2005). GCMs are missing a fundamental geometrical constraint that leads to the dominant scaling in the equatorial region. However, at higher latitudes, only models that are restricted to a thin
shell (inner-to-outer radius ratio > 0.9) appear capable of reproducing the multiple, mid-latitude jets, which are too few and too broad in simulations using deep convective models (Vasavada and Showman, 2005).

By running idealized simulations in a simplified GCM, with a well-controlled energy forcing and for various rotation rates, we reproduce flows with a varying degree of agreement with the zonostrophic regime, framing the dynamics of our Saturn reference simulation. A set of four idealized simulations allowed us to clearly identify three inertial ranges with distinctive dynamics and specific spectral signatures:

* A residual-dominated range ($n_β < n < n_i$) characterized by a KBK $n^{-5/3}$ slope in the residual kinetic energy spectrum. The energy flux is stably negative and the entrophy flux is nearly zero.

* A zonostrophic inertial range ($n_R < n < n_β$), marked by a steeper $≈ n^{-3}$ slope with a clear maximum in the residual energy spectrum, as well as a succession of erratic upscale and downscale energy and entrophy transfers that dominate the mean dynamics. It has been suggested that the steepening of the residual energy spectrum results from the pile-up of non-axisymmetric modes at the jet scale that are produced by jet instability (likely barotropic instability). Such an energy accumulation at large scales is in turn likely responsible for the succession of positive and negative spectral fluxes at large scales. It has also been mentioned that bottom friction possibly affects the dynamics at large scale. It might, for example, be partially responsible for the energetic bump ($n^{-3}$) in the residual spectra. If such a feature is not retrieved from direct observations of Saturn's flow, then GCMs might improve their predictions by implementing a more elaborate bottom drag formulation (Sukoriansky et al., 1999).

* A large scale range, for indices $n < n_R$, dominated by friction and marked by an energy decrease. Energy transfers are nonetheless efficient up to this large scale range.

These three spectral ranges are observed in our Saturn reference simulation, assuming that kinetic energy is injected near the Rossby deformation radius. Only the dynamics within the large scale range differ from the idealized simulations, where it is marked by the presence of the two long-lived
modes $m = 2$ and $m = 3$. The main differences between the Saturn reference simulation and our idealized simulations are its intermittent-anisotropic forcing resulting from the growth of stratospheric baroclinic instabilities over seasonal cycles. As a consequence, while the non-dimensional quantities computed from the SRS are similar to those from our $\Omega/2$ idealized simulation, the ratio between rotational and divergent energy is significantly different, with $E_{rot}/E_{div} \approx 100$ for the SRS compared with $\approx 1000$ for the idealized case (see Table 2). The periodic energy injection in the topmost atmospheric layers requires a continuous adjustment of the dynamics that manifest as a strengthened role of the divergent velocity component, leading to the specific statistical features of our Saturn reference simulation.

In this study, we used only the rotational (non-divergent) part of the horizontal velocity field when calculating the spectral fluxes. This spectral flux based on the 2D vorticity equation is primarily responsible for the upscale cascade while the divergent part contribution increases in the downscale cascade range (Augier and Lindborg, 2013). In our configuration the flow is highly barotropized, the upscale cascade covers a wide spectral domain, and the energy in the rotational part of the flow is much larger than in the divergent part, by a factor $\approx 100$ to $1000$ depending on the simulation. In addition, the downscale cascade is presumably poorly resolved due to hyperviscosity effect at small scales. Nevertheless, the divergent part of the spectral flux should be examined more closely in the future and might reveal new statistical features in our Saturn reference simulation. Furthermore, the implementation of atmospheric moist convection as an additional kinetic energy source will significantly impact the dynamics, with its own spectral signature. We also note that the purpose of our study was to disclose the dynamics of Saturn reference simulation. Then, a set of idealized simulations was designed to emphasize the role and the impact on flow statistics of the different parameters previously chosen for SRS (e.g. bottom friction is the same in all simulations). However, we could have reversed our approach and explore a parameter space of dimensionless numbers to identify the most appropriate set of parameters that reproduces Saturn's atmosphere. This has been done by Wang et al. (2018) to frame the regime of terrestrial atmospheric circulation by varying the planetary Rossby number and the frictional Taylor Number. They also compute spectral analysis of the global energy budget in Read et al. (2018), including both kinetic energy and the available potential energy. Such an approach will certainly be the purpose
of a further study with DYNAMICO in order to refine the set of parameters
for our Saturn reference simulation.

A last (but not least) limitation of our approach is the horizontal and
vertical resolution. Baroclinic eddies can develop but are distorted by hy-
peviscosity. Increasing the horizontal resolution from eddy-permitting to
eddy-resolving would likely improve the representation of the global energy
budget. Ultimately, the strength of the upscale energy transfer is decisive in
setting the large scale non-axisymmetric dynamics, marked by the presence
of zonal modes. Dynamical confinement in latitude being fixed by Saturn’s
rotation rate, the strength of the upscale cascade might vary with vertical
resolution and/or depth extension through geometric confinement. However,
the exploration of various layer depths is beyond the scope of the present
study.

Acknowledgments

The authors acknowledge exceptional computing support from Grand
Équipement National de Calcul Intensif (GENCI) and Centre Informatique
National de l'Enseignement Supérieur (CINES). All the simulations presented
in this paper were carried out on the *Occigen* cluster hosted at CINES. This
work was granted access to the High-Performance Computing (HPC) re-
sources of CINES under the allocations A001-0107548, A003-0107548, A004-
0110391 made by GENCI. The authors acknowledge funding from Agence
Nationale de la Recherche (ANR), project HEAT ANR-14-CE23-0010 and
project EMERGIANT ANR-17-CE31-0007. This project has received fund-
ing from the European Union's Horizon 2020 research and innovation pro-
gramme under the Marie Skłodowska-Curie grant agreement N° 797012.
Fruitful discussions with Sandrine Guerlet, Ehouarn Millour, Thomas Du-
bos, Frédéric Hourdin and Alexandre Boissinot from our team helped refine
some discussions in the paper.

**Appendix A. Open access statement**

Datasets related to this article can be found at 10.5281/zenodo.3638105,
an open-source online data repository hosted at Zenodo.org, [dataset] Ca-
banes et al. (2020).
References

Achterberg, R. K., Flasar, F. M., 1996. Planetary-scale thermal waves in Sat-
urn's upper troposphere. Icarus 119, 350-369. doi:10.1006/icar.1996.0024.

Alexakis, A., 2015. Rotating Taylor-Green flow. Journal of Fluid Mechanics 769, 46-78. doi:10.1017/jfm.2015.82.

Augier, P., Lindborg, E., 2013. A new formulation of the spectral energy budget of the atmosphere, with application to two high-resolution General Circulation Models. Journal of the Atmospheric Sciences 70, 2293-2308. doi:10.1175/JAS-D-12-0281.1.

Baines, K. H., Momary, T. W., Fletcher, L. N., Showman, A. P., Roos-Serote, M., Brown, R. H., Buratti, B. J., Clark, R. N., Nicholson, P. D., 2009. Saturn's north polar cyclone and hexagon at depth revealed by Cassini/VIMS. Planetary and Space Science 57, 1671-1681. doi:10.1016/j.pss.2009.06.026.

Boer, G. J., 1983. Homogeneous and isotropic turbulence on the sphere. Journal of the Atmospheric Sciences 40, 154-163. doi:10.1175/1520-0469(1983)040<0154:HAITOT>2.0.CO;2.

Boer, G. J., Shepherd, T. G., 1983. Large-scale two-dimensional turbulence in the atmosphere. Journal of the Atmospheric Sciences 40, 164-184. doi:10.1175/1520-0469(1983)040<0164:LSTDTI>2.0.CO;2.

Burgess, B. H., Erler, A. R., Shepherd, T. G., 2013. The troposphere-to-stratosphere transition in kinetic energy spectra and nonlinear spectral fluxes as seen in ECMWF analyses. Journal of the Atmospheric Sciences 70, 669-687. doi:10.1175/JAS-D-12-0129.1.

Cabanes, S., Aurnou, J., Favier, B., Le Bars, M., 2017. A laboratory model for deep-seated jets on the gas giants. Nature Physics 13, 387-390. doi:10.1038/nphys4001.

Cabanes, S., Favier, B., Le Bars, M., 2018. Some statistical properties of three-dimensional zonostrophic turbulence. Geophysical & Astrophysical Fluid Dynamics 112, 207-221. doi:10.1080/03091929.2018.1467413.

Celani, A., Musacchio, S., Vincenzi, D., 2010. Turbulence in more than two and less than three dimensions. Physical Review Letters 104, 184506. doi:10.1103/PhysRevLett.104.184506.
Charney, J. G., 1971. Geostrophic turbulence. Journal of the Atmospheric Sciences 28, 1087-1095. doi:10.1175/1520-0469(1971)028<1087:GT>2.0.CO;2.

Chekhlov, A., Orszag, S. A., Sukoriansky, S., Galperin, B., Staroselsky, I., 1996. The effect of small-scale forcing on large-scale structures in two-dimensional flows. Physica D 98, 321-334. doi:10.1016/0167-2789(96)00102-9.

Chemke, R., Dror, T., Kaspi, Y., 2016. Barotropic kinetic energy and enstrophy transfers in the atmosphere. Geophysical Research Letters 43, 7725-7734. doi:10.1002/2016GL070350.

Chemke, R., Kaspi, Y., 2015. The latitudinal dependence of atmospheric jet scales and macroturbulent energy cascades. Journal of the Atmospheric Sciences 72, 3891-3907. doi:10.1175/JAS-D-15-0007.1.

Chemke, R., Kaspi, Y., 2015. Poleward migration of eddy-driven jets. Journal of Advances in Modeling Earth Systems 7, 1457-1471. doi:10.1002/2015MS000481.

Chemke, R., Kaspi, Y., 2016. The latitudinal dependence of the oceanic barotropic eddy kinetic energy and macroturbulence energy transport. Geophysical Research Letters 43, 2723-2731. doi:10.1002/2016GL067847.

Chertkov, M., Connaughton, C., Kolokolov, I., Lebedev, V., 2007. Dynamics of energy condensation in two-dimensional turbulence. Physical Review Letters 99, 084501. doi:10.1103/PhysRevLett.99.084501.

Cho, J. Y.-K., Polvani, L. M., 1996. The morphogenesis of bands and zonal winds in the atmospheres on the giant outer planets. Science 273, 335-337. doi:10.1126/science.273.5273.335.

Choi, D. S., Showman, A. P., 2011. Power spectral analysis of Jupiter's clouds and kinetic energy from Cassini. Icarus 216, 597-609. doi:10.1016/j.icarus.2011.10.001.

Christensen, U. R., 2002. Zonal flow driven by strongly supercritical convection in rotating spherical shells. Journal of Fluid Mechanics 470, 115-133. doi:10.1017/S0022112002002008.
[dataset] Cabanes, S., Spiga, A., Young, R. M. B., 2020. Jump - data collection - part i: Jets from a global climate model. Zenodo. doi:10.5281/zenodo.3638105.

Del Genio, A. D., Achterberg, R. K., Baines, K. H., Flasar, F. M., Read, P. L., Sánchez-Lavega, A., Showman, A. P., 2009. Saturn atmospheric structure and dynamics. In: Saturn from Cassini-Huygens. Springer, pp. 113–159. doi:10.1007/978-1-4020-9217-6.

Del Genio, A. D., Barbara, J. M., 2012. Constraints on Saturn's tropospheric general circulation from Cassini ISS images. Icarus 219, 689–700. doi:10.1016/j.icarus.2012.03.035.

Deusebio, E., Boffetta, G., Lindborg, E., Musacchio, S., 2014. Dimensional transition in rotating turbulence. Physical Review E 90, 023005. doi:10.1103/PhysRevE.90.023005.

Dowling, T. E., Bradley, M. E., Colón, E., Kramer, J., LeBeau, R. P., Lee, G. C. H., Mattox, T. I., Morales-Juberías, R., Palotai, C. J., Parimi, V. K., Showman, A. P., 2006. The EPIC atmospheric model with an isentropic/terrain-following hybrid vertical coordinate. Icarus 182, 259–273. doi:10.1016/j.icarus.2006.01.003.

Dubos, T., Dubey, S., Tort, M., Mittal, R., Meurdesoif, Y., Hourdin, F., 2015. DYNAMICO-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility. Geoscientific Model Development 8, 3131–3150. doi:10.5194/gmd-8-3131-2015.

Fletcher, L. N., Orton, G. S., Mousis, O., Yanamandra-Fisher, P., Parrish, P. D., Irwin, P. G. J., Fisher, B. M., Vanzi, L., Fujiyoshi, T., Fuse, T., Simon-Miller, A. A., Edkins, E., Hayward, T. L., De Buizer, J., 2010. Thermal structure and composition of Jupiter's Great Red Spot from high-resolution thermal imaging. Icarus 208, 306–328. doi:10.1016/j.icarus.2010.01.005.

Galanti, E., Kaspi, Y., Miguel, Y., Guillot, T., Durante, D., Racioppa, P., Iess, L., 2019. Saturn's Deep Atmospheric Flows Revealed by the Cassini Grand Finale Gravity Measurements. Geophysical Research Letters 46, 616–624. doi:10.1029/2018GL078087.
Galperin, B., Read, P. L., 2019. Zonal Jets: Phenomenology, Genesis, and Physics. Cambridge University Press. doi:10.1017/9781107358225.

Galperin, B., Sukoriansky, S., Dikovskaya, N., 2010. Geophysical flows with anisotropic turbulence and dispersive waves: flows with a $\beta$-effect. Ocean Dynamics 60, 427-441. doi:10.1007/s10236-010-0278-2.

Galperin, B., Sukoriansky, S., Dikovskaya, N., Read, P. L., Yamazaki, Y. H., Wordsworth, R., 2006. Anisotropic turbulence and zonal jets in rotating flows with a $\beta$-effect. Nonlinear Processes in Geophysics 13, 83-98. doi:10.5194/npg-13-83-2006.

Galperin, B., Sukoriansky, S., Young, R. M. B., Chemke, R., Kaspi, Y., Read, Peter, L., Dikovskaya, N., 2019. Barotropic and zonostrophic turbulence. In: Galperin, B., Read, P. L. (Eds.), Zonal Jets: Phenomenology, Genesis, and Physics. Cambridge University Press, Ch. 13, pp. 220-237.

Galperin, B., Young, R. M. B., Sukoriansky, S., Dikovskaya, N., Read, P. L., Lancaster, A. J., Armstrong, D., 2014. Cassini observations reveal a regime of zonostrophic macroturbulence on Jupiter. Icarus 229, 295-320. doi:10.1016/j.icarus.2013.08.030.

García-Melendo, E., Sánchez-Lavega, A., Legarreta, J., Perez-Hoyos, S., Hueso, R., 2010. A strong high altitude narrow jet detected at Saturn's equator. Geophysical Research Letters 37, L22204. doi:10.1029/2010GL045434.

Gastine, T., Heimpel, M., Wicht, J., 2014. Zonal flow scaling in rapidly-rotating compressible convection. Physics of the Earth and Planetary Interiors 232, 36-50. doi:10.1016/j.pepi.2014.03.011.

Gierasch, P. J., Ingersoll, A. P., Banfield, D., Ewald, S. P., Helfenstein, P., Simon-Miller, A., Vasavada, A., Breneman, H. H., Senske, D. A., Galileo Imaging Team, 2000. Observation of moist convection in Jupiter's atmosphere. Nature 403, 628-630. doi:10.1038/35001017.

Guerlet, S., Fouchet, T., Spiga, A., Flasar, F. M., Fletcher, L. N., Hesman, B. E., Gorius, N., 2018. Equatorial Oscillation and Planetary Wave Activity in Saturn's Stratosphere Through the Cassini Epoch. Journal of Geophysical Research: Planets 123, 246-261. doi:10.1002/2017JE005419.
Guerlet, S., Spiga, A., Sylvestre, M., Indurain, M., Fouchet, T., Leconte, J., Millour, E., Wordsworth, R., Capderou, M., Bézard, B., Forget, F., 2014. Global climate modeling of Saturn's atmosphere. Part I: Evaluation of the radiative transfer model. *Icarus* 238, 110–124. doi:10.1016/j.icarus.2014.05.010.

Guillot, T., Miguel, Y., Militzer, B., Hubbard, W. B., Kaspi, Y., Galanti, E., Cao, H., Helled, R., Wahl, S. M., Iess, L., Folkner, W. M., Stevenson, D. J., Lunine, J. I., Reese, D. R., Biekman, A., Parisi, M., Durante, D., Connerney, J. E. P., Levin, S. M., Bolton, S. J., 2018. A suppression of differential rotation in jupiter's deep interior. *Nature* 555, 227–230. doi:10.1038/nature25775.

Heimpel, M., Aurnou, J., Wicht, J., 2005. Simulation of equatorial and high-latitude jets on jupiter in a deep convection model. *Nature* 438, 193–196. doi:10.1038/nature04208.

Huang, H.-P., Galperin, B., Sukoriansky, S., 2001. Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere. *Physics of Fluids* 13, 225–240. doi:10.1063/1.1327594.

Ingersoll, A. P., 1990. Atmospheric dynamics of the outer planets. *Science* 248, 308–315. doi:10.1126/science.248.4953.308.

Ingersoll, A. P., Pollard, D., 1982. Motion in the interiors and atmospheres of Jupiter and Saturn: Scale analysis, anelastic equations, barotropic stability criterion. *Icarus* 52, 62–80. doi:10.1016/0019-1035(82)90169-5.

Kaspi, Y., Galanti, E., Hubbard, W. B., Stevenson, D. J., Bolton, S. J., Iess, L., Guillot, T., Bloxham, J., Connerney, J. E. P., Cao, H., Durante, D., Folkner, W. M., Helled, R., Ingersoll, A. P., Levin, S. M., Lunine, J. I., Miguel, Y., Militzer, B., Parisi, M., Wahl, S. M., 2018. Jupiter's atmospheric jet streams extend thousands of kilometres deep. *Nature* 555, 223–226. doi:10.1038/nature25793.

Kraichnan, R. H., 1967. Inertial Ranges in Two-Dimensional Turbulence. *Physics of Fluids* 10, 1417–1423. doi:10.1063/1.1762301.

Kraichnan, R. H., 1971. Inertial-range transfer in two- and three-dimensional turbulence. *Journal of Fluid Mechanics* 47, 525–535. doi:10.1017/S0022112071001216.
Lian, Y., Showman, A. P., 2008. Deep jets on gas-giant planets. Icarus 194, 597-615. doi:10.1016/j.icarus.2007.10.014.

Lian, Y., Showman, A. P., 2010. Generation of equatorial jets by large-scale latent heating on the giant planets. Icarus 207, 373-393. doi:10.1016/j.icarus.2009.10.006.

Liu, J., Goldreich, P. M., Stevenson, D. J., 2008. Constraints on deep-seated zonal winds inside Jupiter and Saturn. Icarus 196, 653-664. doi:10.1016/j.icarus.2007.11.036.

Liu, J., Schneider, T., 2010. Mechanisms of jet formation on the giant planets. Journal of the Atmospheric Sciences 67, 3652-3672. doi:10.1175/2010JAS3492.1.

Liu, J., Schneider, T., 2015. Scaling of off-equatorial jets in giant planet atmospheres. Journal of the Atmospheric Sciences 72 (1), 389-408.

Panetta, R. L., 1993. Zonal jets in wide baroclinically unstable regions: Persistence and scale selection. Journal of the Atmospheric Sciences 50, 2073-2106. doi:10.1175/1520-0469(1993)050<2073:ZJIWBU>2.0.CO;2.

Porco, C. C., Baker, E., Barbara, J., Beurle, K., Brahic, A., Burns, J. A., Charnoz, S., Cooper, N., Dawson, D. D., Del Genio, A. D., Denk, T., Dones, L., Dyudina, U., Evans, M. W., Giese, B., Grazier, K., Helfenstein, P., Ingersoll, A. P., Jacobson, R. A., Johnson, T. V., McEwen, A., Murray, C. D., Neukum, G., Owen, W. M., Perry, J., Roatsch, T., Spitale, J., Squyres, S., Thomas, P., Tiscareno, M., Turtle, E., Vasavada, A. R., Veverka, J., Wagner, R., West, R., 2005. Cassini imaging science: Initial results on Saturn's rings and small satellites. Science 307, 1226-1236. doi:10.1126/science.1108056.

Porco, C. C., West, R. A., McEwen, A., Del Genio, A. D., Ingersoll, A. P., Thomas, P., Squyres, S., Dones, L., Murray, C. D., Johnson, T. V., Burns, J. A., Brahic, A., Neukum, G., Veverka, J., Barbara, J. M., Denk, T., Evans, M., Ferrier, J. J., Geissler, P., Helfenstein, P., Roatsch, T., Throop, H., Tiscareno, M., Vasavada, A. R., 2003. Cassini imaging of Jupiter's atmosphere, satellites, and rings. Science 299, 1541-1547. doi:10.1126/science.1079462.
Pouquet, A., Marino, R., 2013. Geophysical turbulence and the duality of the energy flow across scales. Physical Review Letters 111, 234501. doi:10.1103/PhysRevLett.111.234501.

Read, P. L., Tabataba-Vakili, F., Wang, Y., Augier, P., Lindborg, E., Valeanu, A., Young, R. M., 2018. Comparative terrestrial atmospheric circulation regimes in simplified global circulation models. part ii: Energy budgets and spectral transfers. Quarterly Journal of the Royal Meteorological Society 144 (717), 2558–2576.

Rhines, P. B., 1975. Waves and turbulence on a beta-plane. Journal of Fluid Mechanics 69, 417–443. doi:10.1017/S0022112075001504.

Rostami, M., Zeitlin, V., Spiga, A., 2017. On the dynamical nature of Saturn’s North Polar hexagon. Icarus 297, 59–70. doi:10.1016/j.icarus.2017.06.006.

Salmon, R., 1980. Baroclinic instability and geostrophic turbulence. Geophysical & Astrophysical Fluid Dynamics 15, 167–211. doi:10.1080/03091928008241178.

Salyk, C., Ingersoll, A. P., Lorre, J., Vasavada, A., Del Genio, A. D., 2006. Interaction between eddies and mean flow in Jupiter’s atmosphere: Analysis of Cassini imaging data. Icarus 185, 430–442. doi:10.1016/j.icarus.2006.08.007.

Sánchez-Lavega, A., del Río-Gaztelurrutia, T., Hueso, R., Pérez-Hoyos, S., García-Melendo, E., Antuñano, A., Mendikoa, I., Rojas, J. F., Lillo, J., Barrado-Navascués, D., Gomez-Forrellad, J. M., Go, C., Peach, D., Barry, T., Milika, D. P., Nicholas, P., Wesley, A., 2014. The long-term steady motion of Saturn’s hexagon and the stability of its enclosed jet stream under seasonal changes. Geophysical Research Letters 41, 1425–1431. doi:10.1002/2013GL059078.

Sánchez-Lavega, A., Hueso, R., Pérez-Hoyos, S., García-Melendo, E., Rojas, J. F., 2004. Observations and models of the general circulation of Jupiter and Saturn. In: Lecture Notes and Essays in Astrophysics. Vol. 1. CreateSpace, pp. 63–86.
Sayanagi, K. M., Showman, A. P., 2007. Effects of a large convective storm on Saturn's equatorial jet. Icarus 187, 520-539. doi:10.1016/j.icarus.2006.10.020.

Schneider, T., Liu, J., 2009. Formation of jets and equatorial superrotation on Jupiter. Journal of the Atmospheric Sciences 66, 579-601. doi:10.1175/2008JAS2798.1.

Shaw, T. A., Shepherd, T. G., 2007. Angular momentum conservation and gravity wave drag parameterization: Implications for climate models. Journal of the Atmospheric Sciences 64, 190-203. doi:10.1175/JAS3823.1.

Showman, A. P., Ingersoll, A. P., Achterberg, R., Kaspi, Y., 2017. The global atmospheric circulation of Saturn. In: Saturn in the 21st Century. Cambridge University Press, pp. 295-336. doi:10.1017/9781316227220.011.

Simon, A. A., Wong, M. H., Rogers, J. H., Orton, G. S., de Pater, I., Asay-Davis, X., Carlson, R. W., Marcus, P. S., 2014. Dramatic Change in Jupiter's Great Red Spot from Spacecraft Observations. The Astrophysical Journal Letters 797, L31. doi:10.1088/2041-8205/797/2/L31.

Smith, L. M., Chasnov, J. R., Waleffe, F., 1996. Crossover from two-to three-dimensional turbulence. Physical Review Letters 77, 2467-2470. doi:10.1103/PhysRevLett.77.2467.

Smith, L. M., Waleffe, F., 1999. Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence. Physics of Fluids 11, 1608-1622. doi:10.1063/1.870022.

Spiga, A., Guerlet, S., Millour, E., Indurain, M., Meurdesoif, Y., Cabanes, S., Dubos, T., Leconte, J., Boissinot, A., Lebonnois, S., Sylvestre, M., Fouchet, T., 2020. Global climate modeling of Saturn's atmosphere. Part II: Multi-annual high-resolution dynamical simulations. Icarus 335, 113-137. doi:10.1016/j.icarus.2019.07.011.

Sukoriansky, S., Dikovskaya, N., Galperin, B., 2007. On the arrest of inverse energy cascade and the Rhines scale. Journal of the Atmospheric Sciences 64, 3312-3327. doi:10.1175/JAS4013.1.
Sukoriansky, S., Dikovskaya, N., Galperin, B., 2008. Nonlinear waves in zonostrophic turbulence. Physical Review Letters 101, 178501. doi:10.1103/PhysRevLett.101.178501.

Sukoriansky, S., Dikovskaya, N., Grimshaw, R., Galperin, B., 2012. Rossby waves and zonons in zonostrophic turbulence. In: AIP Conference Proceedings. Vol. 1439. AIP, pp. 111–122. doi:10.1063/1.3701355.

Sukoriansky, S., Galperin, B., 2016. QNSE theory of turbulence anisotropy and onset of the inverse energy cascade by solid body rotation. Journal of Fluid Mechanics 805, 384–421. doi:10.1017/jfm.2016.568.

Sukoriansky, S., Galperin, B., Chekhlov, A., 1999. Large scale drag representation in simulations of two-dimensional turbulence. Physics of Fluids 11, 3043–3053. doi:10.1063/1.870163.

Sukoriansky, S., Galperin, B., Dikovskaya, N., 2002. Universal spectrum of two-dimensional turbulence on a rotating sphere and some basic features of atmospheric circulation on giant planets. Physical Review Letters 89, 124501. doi:10.1103/PhysRevLett.89.124501.

Vallis, G. K., Maltrud, M. E., 1993. Generation of mean flows and jets on a beta plane and over topography. Journal of Physical Oceanography 23, 1346–1362. doi:10.1175/1520-0485(1993)023<1346:GOMFAJ>2.0.CO;2.

Vasavada, A. R., Showman, A. P., 2005. Jovian atmospheric dynamics: An update after Galileo and Cassini. Reports on Progress in Physics 68, 1935–1996. doi:10.1088/0034-4885/68/8/R06.

Wang, Y., Read, P. L., Tabataba-Vakili, F., Young, R. M., 2018. Comparative terrestrial atmospheric circulation regimes in simplified global circulation models. part i: From cyclostrophic super-rotation to geostrophic turbulence. Quarterly Journal of the Royal Meteorological Society 144 (717), 2537–2557.

Yarom, E., Sharon, E., 2014. Experimental observation of steady inertial wave turbulence in deep rotating flows. Nature Physics 10, 510–514. doi:10.1038/nphys2984.
Young, R. M. B., Read, P. L., 2017. Forward and inverse kinetic energy cascades in Jupiter's turbulent weather layer. Nature Physics 13, 1135–1140. doi:10.1038/nphys4227.

Young, R. M. B., Read, P. L., Wang, Y., 2019. Simulating Jupiter's weather layer. Part I: Jet spin-up in a dry atmosphere. Icarus 326, 225–252. doi:10.1016/j.icarus.2018.12.005.

Zeman, O., 1994. A note on the spectra and decay of rotating homogeneous turbulence. Physics of Fluids 6, 3221–3223. doi:10.1063/1.868053.