Issue 
A&A
Volume 670, February 2023



Article Number  A15  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244278  
Published online  31 January 2023 
Direct driving of simulated planetary jets by upscale energy transfer
^{1}
MaxPlanckInstitut für Sonnensystemforschung,
JustusvonLiebigWeg 3,
37077
Göttingen, Germany
email: boening@mps.mpg.de
^{2}
GeorgAugustUniversität Göttingen,
FriedrichHundPlatz 1,
37077
Göttingen, Germany
Received:
15
June
2022
Accepted:
5
December
2022
Context. The precise mechanism that forms jets and largescale vortices on the giant planets is unknown. An inverse cascade has been suggested by several studies. Alternatively, energy may be directly injected by smallscale convection.
Aims. Our aim is to clarify whether an inverse cascade feeds zonal jets and largescale eddies in a system of rapidly rotating, deep, geostrophic sphericalshell convection.
Methods. We analyze the nonlinear scaletoscale transfer of kinetic energy in such simulations as a function of the azimuthal wave number, m.
Results. We find that the main driving of the jets is associated with upscale transfer directly from the small convective scales to the jets. This transfer is very nonlocal in spectral space, bypassing largescale structures. The jet formation is thus not driven by an inverse cascade. Instead, it is due to a direct driving by Reynolds stresses, statistical correlations of velocity components of the smallscale convective flows. Initial correlations are caused by the effect of uniform background rotation and shell geometry on the flows and provide a seed for the jets. While the jet growth suppresses convection, it increases the correlation of the convective flows, which further amplifies the jet growth until it is balanced by viscous dissipation. To a much smaller extent, energy is transferred upscale to largescale vortices directly from the convective scales, mostly outside the tangent cylinder. There, largescale vortices are not driven by an inverse cascade either. Inside the tangent cylinder, the transfer to largescale vortices is even weaker, but more local in spectral space, leaving open the possibility of an inverse cascade as a driver of largescale vortices. In addition, largescale vortices receive kinetic energy from the jets via forward transfer. We therefore suggest a jet instability as an alternative formation mechanism of largescale vortices. Finally, we find that the jet kinetic energy scales approximatively as ℓ^{−5}, the same as for the socalled zonostrophic regime.
Key words: planets and satellites: interiors / hydrodynamics / instabilities / turbulence
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model.
Open access funding provided by Max Planck Society.
1 Introduction
The flow at the surfaces of the giant planets Jupiter, Saturn, Uranus, and Neptune is dominated by strong zonal jets. In addition, Jupiter’s atmosphere features prominent eddies over a larger range of sizes, the great red spot being the largest. Recent observations of the polar regions by the Juno mission also revealed that Jupiter’s poles are surrounded by larger cyclonic eddies (Adriani et al. 2018; Gavriel & Kaspi 2021). Saturn’s polar hexagon may also be explained by eddies shaping a zonal wind jet (Yadav & Bloxham 2020).
Until the Juno mission, it was an open question as to whether the jets are features of a thin weather layer, as on Earth, or if they reach deep into the planet. The recent gravity measurements of the Juno and Cassini missions show that the winds reach depths of about 3000 km on Jupiter (Guillot et al. 2018; Iess et al. 2018; Kaspi et al. 2018) and about 8000 km on Saturn (Galanti et al. 2019; Iess et al. 2019). The jets are therefore deep and not restricted to a weather layer, but the questions remain as to which mechanism prevents them from penetrating even deeper and at precisely what depth they decay (e.g., Kong et al. 2018; Wicht et al. 2020; Dietrich et al. 2021, and references therein). Recently, it has been suggested that a stably stratified layer could be responsible (Christensen et al. 2020; Wicht & Gastine 2020).
Although deeply penetrating, the jets could still be driven by motions in a shallow weather layer (e.g., Showman 2007; Schneider & Liu 2009; Lian & Showman 2010; Read et al. 2020; Cabanes et al. 2020). Indeed, the dynamics of rotating systems predicts that the jets organize on geostrophic cylinders with minimal variation along the rotation axis, independent of where along the cylinders the flow is driven.
Traditionally, there were two competing theories regarding how giant planetary jets are driven. Both attempt to explain how kinetic energy is transported from a small convective driving scale to the jet scale. The first theory invokes a kinetic energy cascade and the second a direct driving of jets by smallscale convection (see Fig. 1). Recently, some evidence has accumulated in favor of the direct driving idea (Galperin & Read 2019), but there is still no final conclusion on the precise driving mechanism of the jets.
Cascades are known to transport energy from the driving scale to the scales where the energy is dissipated. In a forward cascade, the energy is transferred downscale to smaller scales (Kolmogorov 1941), which is more common for threedimensional dynamics (see, e.g., Alexakis & Biferale 2018, for a review on cascades). An inverse cascade transfers energy upscale to larger scales and is more common for twodimensional dynamics (Kraichnan 1967). Responsible for any energy exchange between different scales is the nonlinear inertial (or advective) term in the Navier–Stokes equation. The range of scales over which this transport to the dissipative scale dominates is therefore called the inertial range. A cascade has two characteristic properties (e.g., Frisch 1995, p. 104). First, the dynamics should be scaleinvariant in the inertial range. This includes a constant scaletoscale flux of energy. Second, the word “cascade” implies a continuous transfer between scales, which should be local in wavenumber space, meaning that several intermediate steps occur between the driving and dissipation scales and that eddyeddy interactions play a major role Richardson (1922, p. 66).
An inverse cascade could therefore explain the transport of energy from the driving scale to the large jet scale (see Fig. 1a). This relates to the idea of eddy merging: starting at the driving scale, eddies would continue to merge to form larger and larger features. One idea is that the merging would stop at the socalled Rhines length scale, where turbulent inertial dynamics give way to Rossby waves (Rhines 1975). In a somewhat similar scenario called the zonostrophic regime (e.g., Galperin et al. 2006; Sukoriansky et al. 2007), a twodimensional inverse cascade operates at intermediate length scales. These studies assume a beta plane approximation or forced turbulence on a rotating spherical surface. Since this inverse cascade is isotropic in spectral space, it does not favor the creation of zonal structures. At scales larger than the PelinovskyVallisMaltrud scale (Pelinovsky 1978; Vallis & Maltrud 1993; see also Sukoriansky et al. 2007), the inverse cascade becomes anisotropic; due to the presence of the beta effect, kinetic energy is preferentially transferred to zonal wave numbers, which leads to the driving of jets. Again, the inverse cascade stops at the Rhines scale in this scenario. Numerical experiments have shown that the anisotropy caused by the variation in Coriolis force with latitude results in a spectral anisotropization of the inverse cascade, which helps explain why flows at the Rhines scale preferably organize in axisymmetric jets rather than in large eddies (see Vasavada & Showman 2005; Sukoriansky et al. 2007, for overviews). The Rhines scale indeed successfully describes the jet width for Jupiter, for Saturn, for laboratory experiments, and for deepshell numerical simulations (e.g., Heimpel & Aurnou 2007; Gastine et al. 2014). This inverse cascade picture of jet formation is traditionally rather connected to twodimensional or shallow general circulation models (e.g., Rhines 1975; Williams 1978; Vallis & Maltrud 1993; Showman 2007; Cabanes et al. 2020). Young & Read (2017) measured the upscale transfer of kinetic energy in observations of flows in Jupiter’s atmosphere obtained from cloud tracking by the Cassini spacecraft and find considerable evidence for the eddyeddy transfer being due to an inverse cascade.
The eddies on Jupiter seem to be shallower, yet also penetrate deeper than the weather layer (e.g., Parisi et al. 2021). Eddy merging and thus a cascade provide a possible driving scenario for these features as well. An inverse cascade to largescale vortices could be expected from threedimensional simulations of rapidly rotating RayleighBénard convection (e.g., Rubio et al. 2014; Favier et al. 2014; Kunnen et al. 2016; Guervilly & Hughes 2017; Novi et al. 2019; Maffei et al. 2021), where the largescale flows have a strong quasitwodimensional component, similar to the sphericalshell case. These authors, however, employ a slightly less strict definition of a cascade by not including the necessity of local interactions. In agreement with these simulations, Siegelman et al. (2022) find an upscale transfer to large scales in the Juno data of Jupiter’s polar vortices. The locality of the transfer has, however, not been constrained, and it is an open question as to whether Jupiter’s polar vortices are driven by an inverse cascade.
An alternative explanation to the cascade model is a direct driving of the jets by Reynolds stresses (see Fig. 1b). Reynolds stresses describe a statistical correlation of smallscale flow components over the jet scale. This picture is usually connected to deepshell simulations of rotating convection (e.g., Aurnou & Olson 2001; Christensen 2001; Heimpel et al. 2005; Kaspi et al. 2009; Gastine & Wicht 2012; Yadav & Bloxham 2020). It is consistent with cloud tracking observations of the surface flows of Jupiter (e.g., Salyk et al. 2006; Ingersoll et al. 2021; Duer et al. 2021), which show that the observed Reynolds stresses are consistent with a direct flux of angular momentum from the eddies into the jets. The corresponding upscale energy flux from the eddies to the jets is about three times stronger than the eddyeddy upscale flux (Young & Read 2017) and thus supports the idea of primarily a direct driving. To achieve a correlation that drives the jets, azimuthal and meridional flows should be sufficiently correlated. The latitudinal scale of this correlation should reflect the jet width. A straightforward way to achieve this is a tilt of geostrophic convective features. At mild parameters, these features may still assume the form of convective columns, while at more extreme parameters, the columns may only exist in a statistical sense. The shear provided by the jets is a very efficient way to cause the tilt. A small initial flow can thus lead to runaway growth until viscous effects provide a balance. These ideas were introduced by Busse under the name mean flow instability (see, e.g., Busse 2002, for an overview). Eddymean flow interactions play the main role in this picture, which is supported by quasilinear simulations (e.g., Tobias et al. 2011; Srinivasan & Young 2012; Marston et al. 2016).
What has been missing is an analysis of the scaletoscale transfer of kinetic energy, which would allow us to determine whether an inverse cascade is operating in deepshell simulations of rapidly rotating convection. We provide such an analysis here. Our main aim is therefore to answer the question of whether in these simulations energy is transferred into the jets in steps involving intermediate and largescale eddies or whether the transfer occurs directly from smallscale convection. Such an analysis would indicate how a deep driving of the jets could be constrained observationally. Since buoyancy does not drive jets, an upscale transfer of kinetic energy is in any case expected from the presence of jets and has been measured in simulations of convection in a spherical shell (Reshetnyak 2013), leaving both a direct driving and a driving by an inverse cascade as possibilities.
We restrict ourselves to studying the transfer between different azimuthal wave numbers, m. Our analysis addresses the question of how kinetic energy, which is driven by buoyancy at typical convective wave numbers m > 0, finds its way into the axisymmetric jets, which have wave number m = 0. In particular, we are interested in the question of which wave numbers drive the jets and whether there is a cascade involving more than one step from the convective driving scale to m = 0. We leave a detailed discussion of the jet width and the Rhines scale (Rhines 1975; Gastine et al. 2014) to future studies, for which an analysis of energy transfer between harmonic degrees, ℓ, is more appropriate (e.g., Young & Read 2017).
Fig. 1 Schematic illustration of two possible mechanisms for jet driving. Top row (a): the inverse cascade idea, where eddies merge until the jet scale is reached. The first step in this process illustrates the required anisotropy necessary to prefer either cyclonic or anticyclonic eddies, here illustrated by a sorting of both eddy types. Bottom row (b): the direct driving via Reynolds stresses, where an initial tilt of convective features (second column) yields weak zonal flows that then amplify the tilt until a viscous balance is reached. 
2 Numerical simulations of deep sphericalshell convection
2.1 Governing equations and numerical implementation
We analyzed the results of direct numerical simulations of convectiondriven flow in a rotating spherical shell with inner radius r_{i} and outer radius r_{o} using the pseudospectral magnetohydrodynamic code MagIC (Wicht 2002a; Gastine & Wicht 2012; Lago et al. 2021, version 5.10). We solved the hydrodynamic equations in the anelastic approximation in a nondimensional form. The momentum equation in a corotating frame is (1)
the continuity equation is (2)
and the entropy equation is (3)(4) (5)
and the viscous dissipation into heat is (6)
Here u is the velocity field, p the modified pressure (including the centrifugal force) and S the entropy. Further, e_{r} and e_{z} are the unit vector in the radial direction and along the axis of rotation, respectively. The quantities , and are purely radial profiles of gravity, density and temperature characterizing the hydrostatic and adiabatic background state on which the fluctuations evolve according to Eqs. (1)–(3). We use a spherical shell with a radius ratio of r_{i}/r_{o} = 0.7 and assume a polytropic ideal gas with a polytropic index of n = 2. The contribution of the shell mass to the total mass can be neglected and thus the gravity profile is . In most of this work, we use spherical polar coordinates (r, θ, φ), while at some locations, we also employ cylindrical coordinates (s, z, φ). Equations (1)–(3) were nondimensionalized by using the shell thickness d = r_{o} − r_{i} as a length scale and the viscous timescale τ_{visc} = d^{2}/v as a timescale. The nondimensional velocity u is therefore expressed in terms of a Reynolds number.
We fix the nondimensional input parameters, such as Ekman (E), Rayleigh (Ra), Prandtl (Pr), and dissipation (Di) numbers to (7) (8) (9) (10)
where Ω is the rotation rate, dS / dr_{o} is the entropy gradient at the outer boundary (the same as at the inner boundary), v is the viscosity, and κ is the thermal diffusivity. In addition, α_{o}, Θ_{o}, and g_{o} are the reference values of thermal expansion coefficient, temperature, and gravity at the outer boundary. For this choice of parameters, the total density contrast across the shell is ρ_{i}/ρ_{o} ≈ 6 and the convection is quite vigorous, yet still rotationally constrained (RaE^{2}/Pr < 1). Similarly for both boundaries, we apply freeslip conditions for the flow and fixed entropy flux conditions, which are relatively standard for simulating planetary jet formation (e.g., Heimpel et al. 2016).
2.2 Rapidly rotating convection with moderately turbulent geostrophic flows
We analyzed a simulation that has been integrated long enough to reach a statistically steady state. In particular, we made sure that the zonal flows were well established and reached a quasiconstant amplitude and structure.
We obtained the final simulation of the statistically steady state from a spinup phase from a preliminary run with eightfold azimuthal symmetry. To guarantee that a statistically steady state was reached, we simulated the preliminary run for over three viscous timescales with a final radial resolution of 217 Chebyshev grid points. We then lifted the azimuthal symmetry by running the simulation for an additional 0.4 viscous times in the entire azimuthal domain to be sure that a statistically steady state was reached. The final simulation used 1024 grid points in longitude and 512 grid points in latitude, resolving spherical harmonic degrees until ℓ ≤ 341. After this equilibrated state has been reached, we continued the simulation for another 0.6 viscous timescales and used 416 realizations from this interval for the subsequent statistical analysis of the steadystate energy transfer.
We show example snapshots from the statistically steady state of this run in Fig. 2. The overall pattern of the flow is generally as expected for rapidly rotating convection in a spherical shell (e.g., Christensen 2001). Due to the Taylor–Proudman theorem, the axisymmetric zonal flows are predominantly geostrophic, and therefore invariant along the z axis, as well as relatively timeinvariant. The nonaxisymmetric, turbulent smallscale convection shows a complex spatial and time dependence with elongated features in the zdirection.
The jet profile is shown in Fig. 2a. It consists of an equatorial prograde jet, flanked by a midlatitude retrograde jet starting roughly at the tangent cylinder (TC; the cylindrical radius of the shell at the equatorial plane), which is again flanked by a weakly prograde polar jet. The jet profile is largely invariant along the z axis, as expected from the Taylor–Proudman theorem. Small deviations due to thermal winds exist. The simulated jet profile is simpler than the actual profile on Jupiter or Saturn (e.g., Vasavada & Showman 2005; Genio et al. 2009). Given that the simulations have to be run for multiple viscous timescales until a statistically steady state is reached, we restrict this study to one case that is somewhat simpler in nature but nevertheless includes the essential physics.
Figure 3 shows kinetic energy spectra, which we averaged over a time interval from the statistically steady state. For the spectra, the flow was divided into toroidal and poloidal components (e.g., Christensen & Wicht 2015). Toroidal flows are flows with radial vorticity and poloidal flows are flows with horizontal divergence including flows driven by convection. The spectra were computed as a function of harmonic degree ℓ and azimuthal wave number m according to Christensen & Wicht (2015). The toroidal kinetic energy for m = 0 is due to the jets, while the poloidal kinetic energy for m = 0 is associated with meridional flows. The zigzag shape of the toroidal kinetic energy in Fig. 3a is due to the hemispheric symmetry of the zonal winds, because the m = 0 component dominates E_{kin}(ℓ) for small ℓ. The peak in the poloidal kinetic energy reflects the convective driving scale, which depends in a complex way on the system parameters and the geometry.
In the following, we distinguish different characteristic azimuthal scales with different properties in the fluid motions based on Figs. 2 and 3 as well as based on our results (Sect. 4). The flows consist of a strong zonal flow and of a weaker meridional flow (both at the jet scale, m = 0), of largescale eddies (0 < m ≲ 20), of smallscale convection (20 ≲ m ≲ 100), and of flows at the classical dissipation scales (m ≳ 100). Toroidal zonal flows are dominant at large scales, while poloidal and toroidal flows have a similar amplitude at the convective and classical dissipation scales (see Fig. 3).
The simulated flows and in particular the zonal winds are predominantly geostrophic, that is, variations in the directions of the rotation axis are minimal. Geostrophy is enforced by a primary force balance between Coriolis force and pressure gradient. This is clearly evident from Fig. 4, which shows the square root of the shellaveraged squared forces in spectral space (rms force spectra). The Coriolis force dominates at all scales, as it is typical for sphericalshell convection at small Ekman numbers.
Recently, the idea of a zonostrophic regime has been discussed in the context of jet formation in twodimensional forced turbulence on a βplane or rotating spherical surface (e.g., Galperin et al. 2019), in shallow global climate models (Cabanes et al. 2020), and in an experimental study of forced rotating turbulence (Lemasquerier et al. 2023). The zonostrophic regime is characterized by strong jets, which develop between the (large) Rhines scale and the (sufficiently smaller) PelinovskyVallisMaltrud scale (Galperin et al. 2006; Sukoriansky et al. 2007). In the corresponding spectral region, the kinetic energy spectrum is dominated by zonal flows and shows a characteristic ℓ^{−5} scaling. Interestingly, an approximate ℓ^{−5} scaling can be seen in our solution for ℓ < 20 (see Fig. 3a). The scale ℓ = 20 marks the scale where buoyancy starts to dominate inertia in the rms force spectra (see Fig. 4), as expected for thermal Rossby waves (e.g., Roberts 1968; Busse 1970). This is consistent with the interpretation that the ℓ^{−5} range is dominated by Rossby waves (Sukoriansky et al. 2007).
We therefore estimated the Rhines scale, ℓ_{Rh}, and the Pelinovsky–Vallis–Maltrud scale, ℓ_{β}, and find (11) (12)
where we used β = Ω/R = E^{−1}/r_{o} (see also Galperin et al. 2019) and we estimated ℓ_{β} as the crossover scale, where an average Rossbywave period τ_{RW} = 2πℓ/Ω equals the turnover timescale τ_{t} = 2π(U_{ℓ}ℓ/R)^{−1}, where (similar to Rhines 1975) and M is the total mass of the shell. We therefore find a zonostrophy index of ℓ_{β}/ℓ_{Rh} = 6.4, which implies that the ℓ^{−5} spectral range is sufficiently large for the zonostrophic regime (e.g., Sukoriansky et al. 2007).
We however also find differences to the zonostrophic regime. Firstly, the estimated ℓ_{β} ≈ 30 is very close to the forcing scale, which is around ℓ_{F} ≈ 40. Strictly speaking, there has to be a minimum scale separation of at least a factor of two between the forcing scale and ℓ_{β} to permit an isotropic inverse cascade to develop (e.g., Galperin et al. 2006, 2019; Sukoriansky et al. 2007). From zonostrophic theory, one might therefore expect that the inverse cascade or the characteristic ℓ^{−5} scaling do not show clearly in our case. On the other hand, it is possible that the forcing scale does not matter much for the development of the ℓ^{−5} scaling, as was also speculated by Lemasquerier et al. (2023) but which is in contradiction to a statement by Galperin et al. (2006).
Further research is necessary to firmly conclude on the ℓ^{−5} power law scaling of jets in deepshell convection, in particular an analysis of spectral scaling and spectral transfer in harmonic degree, ideally including a suite of simulations spanning a larger range of parameters and with a number of jets inside the TC. At small scales, it is possible to fit an ℓ^{−5/3} slope to our spectra, but this slope is not unique due to the curved shape of the spectrum.
Fig. 2 Typical simulation snapshots in the statistically steady state. Panel (a) shows the jet profile, . Panels (b) and (c) show the fluctuations around the zonal mean of the zonal and meridional winds, and , at the central meridian. Panel (d) shows at the top boundary. Panels (b)–(d) are saturated at half of the maximum value in the domain shown. 
Fig. 3 Kinetic energy spectra timeaveraged over an interval that is part of the statistically steady state of our simulation as a function of harmonic degree (a) and azimuthal order (b). Both panels show separate spectra for the toroidal and poloidal flow components. See Sect. 2.2 for a brief discussion of the possible power laws shown in panel (a). 
Fig. 4 Square root of the shellaveraged squared forces in spectral space as a function of harmonic degree ℓ (also known as rms force spectra). The Coriolis force is dominating at all scales as is typical for rapidly rotating sphericalshell convection. 
3 Kinetic energy transfer between azimuthal wave numbers in a spherical shell
3.1 Scalebyscale balance of kinetic energy
We here derive the expressions for the balance and flux of kinetic energy at azimuthal wave numbers m using the Cartesian formulation of Alexakis & Biferale (2018) as a guideline. To do this, we perform a Fourier transform in φ of all variables using the normalization of Lago et al. (2021). The kinetic energy is (13)
where V is the full spherical shell volume and (14)
is the kinetic energy at azimuthal wave number m. Here, u^{(m)} is the flow in physical space filtered in such a way as to retain only the mth azimuthal Fourier component. The equations in this section are valid both for decomposing into spherical and cylindrical coordinates. By taking the time derivative of the kinetic energy (Eqs. (13) and (14)) and inserting Eq. (1), we obtain the kinetic energy balance (15)
where the individual contributing terms are due to nonlinear advection (T), work done by pressure (F_{P}), by buoyancy (F_{Buoy}), by the Coriolis force (F_{C}), and by viscous dissipation (D), (16) (17) (18) (19) (20)
We here use 〈a, b〉 for the vector product a · b and all expressions are functions of time. The Coriolis force merely transfers energy between the s and φ flow components and therefore does no net work. Consequently, the integrand of F_{C}(m) is zero at every location. The work done by the pressure force is zero at any m because we integrate over the entire shell.
In addition to the global energy balance, we also consider the distribution of the local contributions to Eq. (15) in the meridional plane (r, θ). To this end, we replace the volume integrals in Eqs. (16) to (20) by a simple integration over the φcoordinate for each (r, θ) and denote the result with the same letter for simplicity. More formally, for any of the quantities Q = E_{kin}, T, F_{Buoy}, we write (21)
3.2 Kinetic energy transfer from scale to scale
Only the nonlinear advection term T(m) is able to transport kinetic energy from one spectral scale to another. We evaluate this scaletoscale transfer of kinetic energy by T(m, m′), which is a scalebyscale version of T(m) and given by (22)
The relationship between T(m) and T(m, m′) is given by (23)
The energy transfer function, T(m, m′), describes how nonlinear advection transports kinetic energy from scale m′ to scale m. Positive T(m, m′) means a transfer of energy from wave number m′ to m, which is also known as forward transfer if it is to smaller scales (m > m′) and upscale transfer if it is to larger scales (m < m′). We note that different choices for the sign of T appear across the literature (e.g., Rubio et al. 2014; Favier et al. 2014; Alexakis & Biferale 2018). Due to energy conservation, we have T(m, m′) = −T(m′, m).
3.3 Energy transfer to the jets
The m = 0 component of the kinetic energy represents mainly the zonal winds, since the meridional circulation, the second axisymmetric flow contribution, is orders of magnitude smaller. Since the axisymmetric zonal winds are a toroidal flow, they do not experience a driving by buoyancy. Therefore, advection mainly feeds energy into the jets from the eddies (m > 0). Using simple orthogonality relations and a partial integration, one can show that (24)
where we have used the summation convention over i and the overbar denotes an azimuthal average. As a consequence, the total energy transfer to the m = 0 flow is (25)
is the Reynolds stress from the entire flow field. Here, the quantity (27)
is the Reynolds stress, which is responsible for the energy transfer from a fluid scale m′ > 0 to the jet scale m = 0 (see Eq. (24)). The energy transfer to the jets is entirely due to Reynolds stresses. The value T(0, 0) describes the generation and redistribution of zonal flow by meridional circulation. This is very small locally and averages to zero over the whole convective shell.
Fig. 5 Scalebyscale balance (a) and scaletoscale flux (b) of kinetic energy in the statistically steady state. The values T(0) ≈ D(0) ≈ 7 × 10^{9} were excluded from the plotting range in panel (a) for better visibility. 
3.4 Scaletoscale flux of kinetic energy
In addition to the scalebyscale balance of kinetic energy, it is interesting to consider the flux of kinetic energy between different ranges of scales, which is mediated by the advection term T(m). This can be obtained by a cumulative sum over m, (28) (29) (30)
where we used T (m, m′) = −T(m′, m) in Eq. (29). The scaletoscale flux Π(m) is the kinetic energy received by scales from scales by advection and thus characterizes the energy flux across scale m. Positive Π(m) means transfer across m from larger to smaller scales, which is also known as forward flux, and negative Π(m) is transfer across m to larger scales, which is also known as upscale (or inverse) flux.
In total, the nonlinear advection term therefore acts neither as a source nor as a sink of kinetic energy because lim_{m→∞} Π(m) = 0. The total kinetic energy balance for the entire shell reads (31)
and buoyancy is the global energy source and dissipation the global energy sink.
When the simulation reaches a statistically steady state, the temporal mean of the kinetic energy does not change any more to a significant extent. From here on, we consider temporal averages of the above terms over sufficient time such that their values are close to their longterm means. The time derivative of the timeaveraged kinetic energy is thus very small and can be neglected, ∂_{t}E_{kin}(m) = 0, and the balances (15) and (28) become (32) (33)
Equations (32) and (33) show that the advection term balances the source and sink terms of the kinetic energy at each m. It transports kinetic energy away from the buoyantly forced scales to the dissipation scales.
4 Results for the statistically steady state
4.1 Scaletoscale flux: Upscale flux at large scales, forward flux at small scales
As a first result, we show the scalebyscale spectral balance of kinetic energy in Fig. 5a. Buoyancy work drives the kinetic energy at all scales, in agreement with Cartesian simulations of Rayleigh–Bénard convection (e.g., Verma 2019, Fig. 3). The dominant contribution (79%) to the convective driving comes from the small convective scales. The amplitude of the convective driving drops steeply for larger wave numbers but retains a sizable value for the smallest wave numbers. At all scales except m = 0 and m ≳ 200, buoyancy is the strongest driver of the motions. The small buoyant driving at m = 0 injects energy into the meridional flow.
From the small convective scales, nonlinear advection transfers kinetic energy away to other scales (T is negative). It is transferred to the jets and to largescale eddies, as well as to the classical dissipation scales, where T is positive. As a consequence, there is a net upscale flux of kinetic energy at larger scales (Π < 0 at m < 60). The scaletoscale flux is relatively constant across the largest scales (m < 20), so that the constant scaletoscale flux criterion for an inverse cascade is satisfied (e.g., Young & Read 2017; Alexakis & Biferale 2018).
At smaller scales, there is a net forward flux of kinetic energy to ever smaller scales with Π positive and nearly constant for 80 ≲ m ≲ 150. Such an approximately constant forward flux might similarly be associated with the inertial range of a forward cascade.
Viscous dissipation extracts energy from all scales as expected. The small convective scales experience the strongest dissipation (46%), while the classical dissipation scales give rise to only a fraction of the dissipation (18%). Dissipation at the jet scale contributes 26% to the total dissipation. Since we use freeslip boundary conditions, there is no Ekman boundary layer and the dissipation in general takes place in the bulk of the domain. There is some tendency for increased dissipation toward the outer boundary of the shell, likely because the convective motions become stronger in that region. The strong role of dissipation at the convective and large scales (m < 100) is most likely due to the use of a turbulent viscosity. The viscosity of actual planetary interiors is such that the Ekman number is much smaller than in our study, which results in viscosity acting only on the smallest scales. In an actual planet, nonlinear advection to these unresolved scales is likely to play the role that the turbulent viscosity plays in our simulations.
Fig. 6 Scaletoscale transfer of kinetic energy, T (m, m′), in the statistically steady state. Note that the figure has logarithmic axes and two colorbars; the left one indicates values for m = 0 or m′ = 0 and the right one for m, m′ > 0. A positive value of T(m, m′) shows that kinetic energy is transferred from wave number m′ to wave number m. We use the regions indicated by gray lines and labeled with letters (a) to (e) to refer to specific regions of the plot in the text and show the different spatial contributions to T(m, m′) for those regions in Fig. 7. 
4.2 Results for the scaletoscale transfer of kinetic energy
4.2.1 Upscale transfer of kinetic energy into jets occurs directly from most buoyant scales via Reynolds stresses
We show results for the energy transfer function in Fig. 6. We highlight different features in the bidimensional wavenumber space (m, m′) that are indicative of different modes of energy transfer by gray lines or polygons and discuss their meaning below. Overall, we find the largest values for the energy transfer function for m = 0 and m′ = 0. In regions that involve m = 0 or m′ = 0 the amplitude of the transfer is roughly 30 times larger than elsewhere. From the large amplitude in the region labeled (a), it is clear that this is energy transfer into the jets, rather than into the much weaker meridional flow. This transfer is caused by statistical correlations of flow components (Reynolds stresses; see Sect. 3.3). The energy transfer by Reynolds stresses arises predominantly from the small convective scales that experience the strongest convective driving (region (a) in Fig. 6). This transfer is extremely nonlocal in wavenumber space. It takes place directly from the convectively forced scales to the jets and bypasses a transfer into scales of intermediate size, from which the energy could have been transferred to the jets in a later step. We thus find that the jets are in the statistically steady state directly driven by the small convective scales rather than by an inverse cascade. We further underpin this conclusion in Sect. 5 by analyzing the evolution of the simulation until the statistically steady state. The steadystate transfer of kinetic energy from the convective scales to the jet resembles to some extent the energy transfer in simulations of rapidly rotating RayleighBénard convection (e.g., Rubio et al. 2014; Favier et al. 2014; Kunnen et al. 2016). These authors, however, seem to associate this upscale transfer with an inverse cascade despite the nonlocality of the transfers in their simulations.
Figure 7 shows the spatially resolved contribution to the energy transfer function. Here, we have summed the local contribution T(m, m′, r, θ) over all (m, m′) that are part of a spectral region labeled in Fig. 6. From Fig. 7a, we take that the injection of energy into the jets is rather localized at the centers of the prograde and retrograde jet flow.
4.2.2 Nonlocal upscale transfer to large vortices
For the transfer of energy to largescale eddies, we observe two different types of behavior. Firstly, we find a region of upscale energy transfer from the scales with the largest buoyancy driving to comparatively larger scales (regions (b) and (b′) in Fig. 6). This way, kinetic energy is transferred to larger scales approaching the jet scale. In this region, the more nonlocal transfers (region b) are stronger than the transfer in the more local region (b′) by a factor of over 50, when summed over the entire region. The total transfer from region (b) is about seven times smaller than the transfer from region (a).
Region (b) has a small negative contribution to the upscale flux shown in Fig. 5b. The transfers are however quite nonlocal and can be understood as follows. The dominant convective scale of m′ ≈ 35 can be expected to have a strong nonlinear interaction with other scales of similar size, say m″ ≈ 20–50. The resulting third coupling partner then necessarily has the scale m = m′ ± m″. As rotational effects induce upscale transfers in addition to the usual forward transfers, this results in transfers to m = m′ − m″ ≈ 0–15. This feature is therefore most likely a side effect of the transfer from the convective scales to the jet scale not always exactly matching m = 0. The correlations of convective motions are not always perfect and may thus instead of feeding m = 0 also feed m = 1, 2, … to a lesser extent. This possibility is supported by the maximum of the upscale transfer for each m being always near the strongest convective scale, m′ ≈ 35. This scenario seems more plausible than interpreting region (b) as a signature of an inverse cascade because it explains the nonlocality of the transfers. We therefore conclude that the upscale transfers in region (b) are likely not associated with an inverse cascade.
In region (b′), we find local positive transfers that could be a signature of an inverse cascade. In this region, however, the transfers are significantly weaker than in the neighboring region (b) and contribute barely to the upscale flux. We therefore do not interpret this feature as a signature of a global inverse cascade, which operates in the entire domain. The transfers in region (b′) are predominantly caused by flows inside the TC. We therefore show the energy transfer function integrated only over the region inside the TC in Fig. 8. We find that the transfers in region (b′) inside the TC are closer to the diagonal and therefore more local. Our results therefore indicate the possibility that largescale vortices are driven by an inverse cascade inside the TC. The jets, however, are clearly also driven directly by the small convective scales also inside the TC (see region a in Fig. 8).
Fig. 7 Spatial contribution, T_{i}(r, θ), to the energy transfer function, T(m, m′), summed over different values of m and m’ according to the masks (i) = (a),(b),(c),(d),(e) shown in Fig. 6. The panels display results for the direct driving of the jets by the convective scales via upscale energy transfers and Reynolds stresses (a), the upscale transfer to large eddies (b), the feeding of largescale vortices from the jet energy (c), the smallscale forward cascade (d), and the local forward transfer likely associated with a forward enstrophy transfer (e). All panels are saturated at three times their maximum value, except for panel (a), which is saturated at 18 times its maximum value. 
4.2.3 Transfer from the jets to largescale vortices: Jet instability as an alternative formation mechanism
In addition to the upscale transfer, we find that largescale vortices receive kinetic energy from the jets (see region (c) in Fig. 6 and panel (c) in Fig. 7). This forward transfer from the jets is quite confined in m compared to the rather broad upscale transfer region (b).
The driving of largescale vortices by a transfer from the jets is further supported by a comparison to the toroidal energy spectra. We find that there is a small peak in toroidal kinetic energy at those scales and centered around m = 5 (see Fig. 3b). The spatially resolved contribution to the kinetic energy, E_{kin}(m, r, θ) further supports this conclusion. Figure 9 shows that for m = 3, …, 8, there is a noticeable excess power near the location of our retrograde jet at s ≈ 0.5.
These features also appear in snapshots of our simulation (see the region of latitudes between ±50 and ±65 degrees in both hemispheres in Fig. 10). At these latitudes, structures in m = 3, …, 8 are visible in both (shown) and (not shown) that appear phaseshifted so as to form vortices.
Over the entire shell, the transfer of kinetic energy from the jets to these largescale vortices is smaller than the work done by buoyancy in our simulations (see a comparison between F_{Buoy}(m) and T(m) at 3 ≤ m ≤ 9 in Fig. 5). In the region of the lower part of the shell at the retrograde jet, where the m = 3, …, 8 wave numbers are strongest, however, the energy transfer due to T(m, 0) is larger than the driving by buoyancy. This is evident from a comparison of Figs. 7c and 11a.
Figure 11b shows the spatial dependence of the transfer to the largescale vortices from all other scales m′ > 0, (34)
We find that the driving of the largescale vortices by other scales except the jet scale is weaker than the driving from the jets in the deeper regions, but stronger in the nearsurface and equatorial regions.
We conclude that the deeper largescale vortices are predominantly driven by a mechanism that transfers energy from the jets to these scales. The origin of this transfer of kinetic energy from the jets to the largescale vortices is likely an instability of the jets. We checked the zonal flow instability criterion from Wicht (2002b, which is the same as for barotropic instability), which however was not fulfilled and this instability can be ruled out.
Fig. 8 Energy transfer function, T(m, m′), in the statistically steady state, integrated over the region inside the TC (s < r_{i}). 
Fig. 9 Spatial dependence of the contribution to the kinetic energy in the statistically steady state, E_{kin}, for a sum over azimuthal wave numbers m = 3,…, 8 (see label (c) in Fig. 6). See Eq. (21) for a definition of E_{kin}(m, r, θ). 
4.2.4 Forward cascade to small scales
Toward the classical dissipation scales (m, m′ > 60), we find in agreement with the energy flux in Fig. 5 that kinetic energy is transferred to ever smaller scales (region (d) in Fig. 6). This transfer occurs in almost all of the spherical shell and in a stronger way toward the surface, where convective motions experience the strongest driving (see Figs. 7d and 11a).
The logarithmic axes in Fig. 6 give the impression that the width of the region of forward transfer shrinks with m. However, this is not the case and the width of this region is indicative of the convective driving. The forward transfer is not strictly local but slightly nonlocal, meaning not only to directly neighboring scales but mostly to scales in a close neighborhood. Such a somewhat nonlocal transfer is expected for a driving that is broad in wavenumber space like our driving by buoyancy (see F_{buoy} in Fig. 5 and Kuczaj et al. 2006). For turbulent convection, similar results have been observed in several Cartesian studies (e.g., Mishra & Verma 2010; Moll et al. 2011; Verma et al. 2017; Alexakis & Biferale 2018) in both rapidly rotating and nonrotating systems (e.g., Favier et al. 2014; Verma et al. 2017). The nonlocality of the transfers (i.e., the distance to the diagonal) can be explained by the necessity of the two wave numbers m and m′ to have a coupling partner, for which only m − m′ and m + m′ are an option. A closer inspection of Fig. 6 shows that the main coupling partners for m are around m′ ≈ m ± 35, which is the mean distance of region (d) from the diagonal. This can be explained by the main convective driving and main convective power being at the same scales m − m′ = m″ ≈ 35 (see the poloidal spectrum in Fig. 3b and F_{Buoy} in Fig. 5). The forward transfer region observable in our simulation thus resembles a relatively local forward cascade (e.g., Alexakis & Biferale 2018).
4.2.5 Local forward transfer
Furthermore, there is a very narrow region of forward transfer of kinetic energy close to the diagonal, which at first sight looks like the sign of a forward cascade (region (e) in Fig. 6). Since this pattern starts at m ~ 10 and since it strongly involves the interaction with the m = 1 component, we suppose that this feature is related to the dynamics of the depthaveraged geostrophic component of the flow (see Rubio et al. 2014). This is supported by panel (e) in Fig. 7, which shows that the interaction is due to a very largescale component with a strong alignment in z. Similar results appear in Maltrud & Vallis (1993, Figs. 7a,c) in forced twodimensional turbulence and in Rubio et al. (2014, bottom left panel in Fig. 5) in rapidly rotating RayleighBénard convection. These authors attributed similar features to the signature of a forward enstrophy transfer in the energy transfer function. In twodimensional flows, the domainintegrated enstrophy is an invariant in addition to the kinetic energy and also shows forward cascading behavior (e.g., Kraichnan 1967; Alexakis & Biferale 2018). The forward enstrophy cascade shows up in the kinetic energy transfer as a local forward transfer (see also Verma et al. 2005), but without contributing to the kinetic energy flux (see discussion in Maltrud & Vallis 1993 and Kraichnan 1967). In agreement with this, we find that this region does not have a significant contribution to the scaletoscale flux. Forward enstrophy transfer likely results from stretching of vortices by the background largescale shear and is associated with the upscale flux of energy to the largest scales (e.g., compare to Maltrud & Vallis 1993; Chen et al. 2006; Xiao et al. 2009). In line with this, we find this transfer in our simulation to be strongest in the region of strongest shear in the zonal flow (see panel (e) in Fig. 7). These arguments support an attribution of this feature to forward enstrophy transfer. A separate analysis of the spectral transfer of enstrophy is necessary to firmly conclude this, which is left to future studies.
Fig. 10 Snapshots of the statistically steady state at the shell center (r/r_{o} = 0.85). The fluctuations around the meridional flow, , clearly show structures of m ≈ 5 at higher latitudes. 
Fig. 11 Spatial dependence of the driving of kinetic energy of largescale eddies by buoyancy, F_{Buoy,c} (a), and kinetic energy transfer to the same eddies from other eddies, T_{c,m′>0} (b), both in the statistically steady state. For both panels, we summed T and F_{Buoy} over m as given by label (c) in Fig. 6 and summed T over m′ > 0. 
5 Results for the spinup phase
Once zonal flows have been established, the direct driving via Reynolds stresses dominates the feeding of energy into the jets. However, the mechanism that initializes the runaway growth could have a different origin. Busse (2002) speculates that the boundary curvature establishes a small initial tilt because Rossby waves in the equatorial region tend to drift faster in prograde direction than at depth. This could explain the simple structure typically found outside the TC in simulations of convection in a deep rotating shell. In these simulations, a prograde equatorial jet dominates a flanking retrograde jet. However, for the multiple jets developing in many simulations inside the TC, a more complex explanation is required. Busse originally envisioned nested layers of convective columns, but this has never been confirmed in any numerical simulation.
The Rhines scale more or less correctly predicts the jet width inside but not outside the TC (Gastine et al. 2014). This could indicate that the jet formation mechanisms are different inside and outside the TC. Inside the TC, an initial cascade could provide the jet structure while outside TC, it is not present. Once the jets are strong enough, however, the direct driving via Reynolds stress would dominate due to the dominating effect of the zonal flows on the dynamics inside and outside the TC. A similar behavior has for example been observed in forced twodimensional turbulence on a sphere (Nozawa & Yoden 1997a,b; Huang & Robinson 1998).
In order to explore this possibility, we run a spinup simulation that starts from an initial state with developed convective flows but without jets and evolves for about 0.21 viscous timescales. The initial condition is established somewhat artificially by forcing the axisymmetric contributions of flow, entropy and pressure to remain zero for a period of 0.15 viscous timescales. The exception are the spherically symmetric contribution to entropy and pressure, which describes the diffusive background state. The manipulated simulation was long enough to establish a new statistically steady state.
Figure 12 gives an overview how the total zonal flow energy (a), the surface zonal flows, the poloidal energy (b), the volumeaveraged Reynolds stress (c) and correlation C_{R} of convective flows (d), and the surface zonal jets (e) develop during the spin up. The volumeaveraged Reynolds stress (Christensen 2002; Gastine & Wicht 2012), (35)
where V is the shell volume, relies on the correlation of nonaxisymmetric flow contributions in s and the azimuthal direction, (36)
The manipulated initial state is characterized by particularly strong convective motions (b), which are suppressed during the spinup by the growing zonal flows.
Even without zonal flows, the manipulated simulation shows a correlation C_{R} and thus some Reynolds stress. The pattern of this Reynolds stress is illustrated in Fig. 13a. Correlation and stress are mostly positive (red) in the region outside the TC. This is consistent with the theory outlined by Busse (2002) that Rossby waves would travel faster in prograde direction closer to the equator and therefore cause a consistent tilt with a positive correlation C_{R}. Inside the TC, the Reynolds stress is much weaker and the correlation is mostly negative. The negative sign is also consistent with the Rossby wave theory, which would travel retrograde inside the TC and faster closer to the TC. We thus speculate that this Rossby wave theory largely explains the initial correlation when no zonal winds are present.
Figure 12 demonstrates that jets first appear outside the TC where the Reynolds stress is strongest. Over the spinup process, the correlation increases (d) along with the zonal winds (a and e). The Reynolds stress amplitude (c), however, decreases, since the convective flow (b) becomes weaker.
Figure 13b shows that the Reynolds stress distribution after about ten turnover times (or one percent of a viscous dissipation time) remains largely unchanged. However, at the end of the spinup process (panel c), the pattern close to and inside the TC has clearly changed. The direction of the tilt of convective feature and thus of the correlation now changes roughly where the retrograde jets have their peak. The initial correlation and tilt structure is thus not the only factor that determines the end result. The simple Rossby wave picture may explain the initial correlation but fails to explain the evolution, at least inside the TC.
Figure 14 shows the energy transfer to the jet scale, T(m = 0, m′), for the first time interval directly after the start of the spinup. This interval has a short duration of less than a convective turnover time, which is at the order of 10^{−3} viscous timescales. During this period, the flows thus largely stem from the initial condition. We find that almost all azimuthal scales contribute to the initial growth of the jets, with a major contribution from the small convective scales (here around m ≈ 25; slightly larger scales compared to the statistically steady state). This shows that from the start of the spinup, the jets are not driven by an inverse cascade, but rather by direct feeding from convective structures. During the spinup, the pattern of the transfer function evolves with the merging and growth of the jets to its steadystate version as the Reynolds stresses evolve (see the plot for the steadystate energy transfer in Fig. 14). We also computed the transfer function T(m, m′) for m, m′ > 0 during spinup (not shown). In the early phase of the spinup, the results are too noisy to permit a conclusion because they have been averaged over a too short period in time (intervals I and II). At the later stages, when patterns appear in the transfer function, they are qualitatively rather similar to the steady state (Fig. 6). Some temporal evolution takes place; the transfer from the jets to largescale vortices (label c in Fig. 6) appears for example only at the latest stage.
Fig. 12 Temporal development of the jets during spinup. We show the kinetic energy of zonal flows (a) and poloidal flows (b, mostly convection), shellaveraged rms Reynolds stresses (c and d), and the surface jet profile (e). See Eqs. (35) and (36) for definitions of the spatial averages of the Reynolds stresses and correlation coefficients. This run was performed using a convective initial condition without jets. The duration of a convective turnover time is about 10^{−3} τ_{visc}. We show results for the spatially resolved Reynolds stresses and the energy transfer to the jets during the intervals marked as gray bands in Figs. 13 and 14. 
Fig. 13 Reynolds stresses R_{s,φ} at different times during the spinup: at the initial condition without jets and uniform background rotation (a), after about 100 turnover times and in the growth phase of the jets (b, same as interval III in panel (e) in Fig. 12) and in the statistically steady state with fully developed jets (c). 
Fig. 14 Kinetic energy transfer to the jets, T(0, m′), for the spinup run. The transfer functions were computed for the first interval directly after the start of the spinup (interval I). The time window is shown as light gray band in Fig. 12. For comparison, we also show results for the statistically steady state from Fig. 6. 
6 Conclusions
In this study, we have determined the spectral transfer of kinetic energy in rapidly rotating sphericalshell convection for moderately turbulent flows. We ran the simulations until a statistically steady state was safely reached and studied the resulting statistical properties of the fluid motions to understand whether an inverse cascade may power jets and largescale vortices.
We find rich dynamics that depend strongly on the azimuthal wave number, that is, on the azimuthal length scale. The dominant energy transfer is upscale and goes directly from the small convective scales into the jets. It is due to statistical correlations of the convective flows (Reynolds stresses). This transfer is very nonlocal in wavenumber space and thus not due to an inverse cascade, although it drives a relatively constant upscale scaletoscale flux at large scales.
The energy transfer to the largescale vortices is about an order of magnitude smaller than the transfer to the jets. Largescale vortices receive energy through different mechanisms. The main driver is buoyancy, but there is also upscale transfer from smallscale convection and forward transfer from the jets. Outside the TC, the upscale transfer is again directly from the small convective scales and thus likely not due to an inverse cascade, although it does contribute to an upscale scaletoscale flux. Inside the TC, the upscale transfer is more local, and we therefore cannot rule out an inverse cascade as a driver of largescale vortices. The transfers inside the TC are, however, quite weak. The forward transfer from the jets to largescale vortices mostly happens in deeper regions and goes to specific wave numbers around m = 5. This transfer takes place predominantly in a region at depth close to the jet center of the retrograde jet, where kinetic energy shows a peak at the same wave numbers around m = 5. We therefore propose an instability of the jets as an alternative mechanism for the driving of largescale vortices.
We also find two qualitatively different regions of relatively local forward transfer. As expected, kinetic energy is transferred from the small convective scales to the classical dissipation scales in a forward cascade. In addition, we find a region of very local forward transfer that does not have a significant contribution to the scaletoscale energy flux. This is likely the signature of local forward enstrophy transfer, a question that is left to future studies.
Analyzing the spinup and growth of the jets more closely, we show that the growing jets suppress convection to some extent but result in an increased correlation in the convective flows, which is likely due to the growing shear. The increasing correlation compensates for the suppression of the convective flows and helps maintain strong Reynolds stresses, which drive the jets. The initial number of jets is continuously reduced by the merging of the jets. This continues until a balance with (turbulent) viscosity is reached.
We used a manipulated initial condition with developed convective flows but where the jet growth has been suppressed for a long simulation time. In this initial state, correlations of the smallscale convection remain due to the effect of uniform rotation and shell geometry. They drive the Reynolds stresses that kick off the jet growth. However, flow correlation and Reynolds stresses change over the spinup process. While the initial correlation seems consistent with the Rossby wave picture of Busse (2002), the evolution and final jet structure obviously require an additional ingredient. Whether a cascade still plays a role remains unclear, at least inside the TC. There is no clear support for a cascade from our transfer analysis.
Previous analyses of simulations with multiple jets indicate that the Rhines scale explains the jet width inside but not outside the TC. We speculate that the strong correlation caused by the consistent tilt of convective features outside the TC dictates the length scale. Consequently, a direct driving of the jets always dominates. Inside the TC, additional factors come into play, but the applicability of the Rhines scale does not necessarily imply the action of a cascade. The results from the spinup may be used in future studies to improve rotatingturbulence models of stellar differential rotation and planetary jets, for example using the socalled A effect (e.g., Rüdiger 1989; Kitchatinov 2013; Barekat et al. 2021).
We find that the classical picture of a stepbystep transfer of kinetic energy from one neighboring scale to another, as in a true inverse cascade, is not a driver of jet formation in our simulation. Instead, eddymean flow interactions are the main drivers of the jets, with eddyeddy interactions playing only a minor role, which supports the results of quasilinear simulations (Tobias et al. 2011; Srinivasan & Young 2012; Marston et al. 2016). Our results are consistent with the dominance of eddymean flow interactions over eddyeddy interactions in observations (Young & Read 2017). This result may be underpinned by a comment by Frisch (1995, p. 251), who stated that an inverse cascade is typically associated with unbounded domains, as the generation of coherent vortical structures at the finite box size seems to ruin the classical picture of the inverse cascade. Accordingly, simulations of inverse cascades in a box need to suppress the growth of the boxscale mode via a strong damping of the largest scales (e.g., Chen et al. 2006; Xiao et al. 2009) to actually obtain an inverse cascade. This is also in agreement with simulations of rapidly rotating Rayleigh–Bénard convection (e.g., Rubio et al. 2014; Favier et al. 2014), where energy accumulates at two boxfilling vortices that are in the statistically steady state and directly driven by upscale transfers from the injection scale. These authors seem to associate their results with an inverse cascade despite their finding of nonlocal upscale transfers to the largest scale.
Here, we analyzed the transfer between azimuthal wave numbers, m, in order to be able to infer, in particular, the transfer to the jets and the dependence on convective scale as well as on latitude. In future studies, we may also analyze the spectral transfer in harmonic degree, ℓ, as considered by Young & Read (2017). This may provide answers to the questions of how the jets merge in the simulations and what sets the resulting number of jets and their width (e.g., Rhines 1975; Gastine et al. 2014). Additionally, such a study may help explain the scaling of the kinetic energy spectrum of the jets, for which our results are similar to the ℓ^{−5} scaling in zonostrophic turbulence (see, e.g., Rhines 1975; Sukoriansky et al. 2007; Galperin et al. 2019). Broadly speaking, our simulation may well be compared to the zonostrophic regime, because of the ℓ^{−5} scaling and because the relevant spectral range is sufficiently large. On the other hand, we also find significant differences compared to the zonostrophic regime, most importantly the absence of an isotropic inverse cascade. This might be because the β effect is already a relevant effect at the forcing scale. In future studies, it would therefore be interesting to study the spectral properties for simulations that span a larger range of parameters, including varying degrees of stratification that would lead to smaller forcing scales. Finally, we note that our simulations do not include a largescale drag, unlike simulations of the zonostrophic regime (see, e.g., Sukoriansky et al. 2007; Galperin et al. 2019); instead, the viscous force is due to viscous friction in the bulk of the domain. It is therefore possible that the ℓ^{−5} scaling applies to a wider class of flows than previously thought.
Our simulation of a rather thick spherical shell produces a main equatorial jet but only two flanking jets in each hemisphere. It would be interesting to analyze more Jupiterlike simulations in a thinner shell that yield a richer jet structure in the future. This would also help clarify the possible differences inside and outside the TC and help us understand what determines the jet width.
The shear exerted by fully developed jets always leads to a significant tilt of convective features. It is therefore no surprise that the direct driving of the jets dominates. This will very likely also be the case for a more complex jet system.
Finally, an analysis of observational data would be very interesting, such as from the Cassini or Juno missions (e.g., Galperin et al. 2014; Young & Read 2017; Siegelman et al. 2022). A comparison between observations, deep spherical shell simulations, and shallow general circulation models (e.g., Schneider & Liu 2009; Cabanes et al. 2020) may help resolve the question of whether the jet driving is deep or shallow.
Acknowledgements
The authors thank the referee, Peter Read, for very helpful and insightful comments which greatly improved the manuscript. V.B. thanks P. Read, R. Yadav, JinHan Xie, K. Julien, M. Rheinhardt, A. Barekat, M. Käpylä, and J. Warnecke for discussions. V.B. thanks R. Cameron for pointing him to the topic of kinetic energy transfer. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the Priority Program SPP 1992 “The Diversity of Exoplanets”. This work used NumPy (Oliphant 2006; van der Walt et al. 2011), matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), and SHTns (Schaeffer 2013).
References
 Adriani, A., Mura, A., Orton, G., et al. 2018, Nature, 555, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Alexakis, A., & Biferale, L. 2018, Phys. Rep., 767769, 1 [Google Scholar]
 Aurnou, J. M., & Olson, P. L. 2001, Geophys. Res. Lett., 28, 2557 [NASA ADS] [CrossRef] [Google Scholar]
 Barekat, A., Käpylä, M. J., Käpylä, P. J., Gilson, E. P., & Ji, H. 2021, A&A, 655, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Busse, F. H. 1970, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 2002, Phys. Fluids, 14, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 Cabanes, S., Spiga, A., & Young, R. M. 2020, Icarus, 345, 113705 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, S., Ecke, R. E., Eyink, G. L., et al. 2006, Phys. Rev. Lett., 96, 084502 [Google Scholar]
 Christensen, U. R. 2001, Geophys. Res. Lett., 28, 2553 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U. R. 2002, J. Fluid Mech., 470, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, U. R., & Wicht, J. 2015, in Treatise on Geophysics, ed. G. Schubert (Amsterdam: Elsevier), 245 [CrossRef] [Google Scholar]
 Christensen, U. R., Wicht, J., & Dietrich, W. 2020, ApJ, 890, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Dietrich, W., Wulff, P., Wicht, J., & Christensen, U. R. 2021, MNRAS, 505, 3177 [NASA ADS] [CrossRef] [Google Scholar]
 Duer, K., Gavriel, N., Galanti, E., et al. 2021, Geophys. Res. Lett., 48, e95651 [NASA ADS] [CrossRef] [Google Scholar]
 Favier, B., Silvers, L. J., & Proctor, M. R. E. 2014, Phys. Fluids, 26, 096605 [NASA ADS] [CrossRef] [Google Scholar]
 Frisch, U. 1995, Turbulence: The Legacy of A. N. Kolmogorov (Cambridge: Cambridge University Press) [Google Scholar]
 Galanti, E., Kaspi, Y., Miguel, Y., et al. 2019, Geophys. Res. Lett., 46, 616 [Google Scholar]
 Galperin, B., & Read, P. L., 2019, Zonal Jets: Phenomenology, Genesis, and Physics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Galperin, B., Sukoriansky, S., Dikovskaya, N., et al. 2006, Nonlinear Process Geophys., 13, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Galperin, B., Young, R. M., Sukoriansky, S., et al. 2014, Icarus, 229, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Galperin, B., Sukoriansky, S., Young, R. M. B., et al. 2019, in Zonal Jets: Phenomenology, Genesis, and Physics, eds. B. Galperin, & P. L. Read (Cambridge: Cambridge University Press), 220 [CrossRef] [Google Scholar]
 Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [Google Scholar]
 Gastine, T., Heimpel, M., & Wicht, J. 2014, Phys. Earth Planet. Interiors, 232, 36 [Google Scholar]
 Gavriel, N., & Kaspi, Y. 2021, Nat. Geosci., 14, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Genio, A. D. D., Achterberg, R. K., Baines, K. H., et al. 2009, in Saturn from CassiniHuygens, eds. M. K. Dougherty, L. W. Esposito, & S. M. Krimigis (Dordrecht and Heidelberg: Springer), 113 [CrossRef] [Google Scholar]
 Guervilly, C., & Hughes, D. W. 2017, Phys. Rev. Fluids, 2, 113503 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, 227 [Google Scholar]
 Heimpel, M., & Aurnou, J. 2007, Icarus, 187, 540 [NASA ADS] [CrossRef] [Google Scholar]
 Heimpel, M., Aurnou, J., & Wicht, J. 2005, Nature, 438, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Heimpel, M., Gastine, T., & Wicht, J. 2016, Nat. Geosci., 9, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, H.P., & Robinson, W. A. 1998, J. Atmos. Sci., 55, 611 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364 [Google Scholar]
 Ingersoll, A. P., Atreya, S., Bolton, S. J., et al. 2021, Geophys. Res. Lett., 48, e2021GL095756 [CrossRef] [Google Scholar]
 Kaspi, Y., Flierl, G. R., & Showman, A. P. 2009, Icarus, 202, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, Y., Galanti, E., Hubbard, W. B., et al. 2018, Nature, 555, 223 [Google Scholar]
 Kitchatinov, L. L. 2013, in IAU Symposium, 294, Solar and Astrophysical Dynamos and Magnetic Activity, eds. A. G. Kosovichev, E. de Gouveia Dal Pino, & Y. Yan, 399 [Google Scholar]
 Kolmogorov, A. 1941, Akad. Nauk SSSR Dokl., 30, 301 [Google Scholar]
 Kong, D., Zhang, K., Schubert, G., & Anderson, J. D. 2018, Proc. Natl. Aca. Sci. U.S.A., 115, 8499 [NASA ADS] [CrossRef] [Google Scholar]
 Kraichnan, R. H. 1967, Phys. Fluids, 10, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 Kuczaj, A. K., Geurts, B. J., & McComb, W. D. 2006, Phys. Rev., 74, 016306 [NASA ADS] [Google Scholar]
 Kunnen, R. P. J., OstillaMónico, R., van der Poel, E. P., Verzicco, R., & Lohse, D. 2016, J. Fluid Mech., 799, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Lago, R., Gastine, T., Dannert, T., Rampp, M., & Wicht, J. 2021, Geosci. Model Dev., 14, 7477 [NASA ADS] [CrossRef] [Google Scholar]
 Lemasquerier, D., Favier, B., & Le Bars, M. 2023, Icarus, 390, 115292 [NASA ADS] [CrossRef] [Google Scholar]
 Lian, Y., & Showman, A. P. 2010, Icarus, 207, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Maffei, S., Krouss, M. J., Julien, K., & Calkins, M. A. 2021, J. Fluid Mech., 913, A18 [NASA ADS] [CrossRef] [Google Scholar]
 Maltrud, M. E., & Vallis, G. K. 1993, Phys. Fluids A: Fluid Dyn., 5, 1760 [NASA ADS] [CrossRef] [Google Scholar]
 Marston, J. B., Chini, G. P., & Tobias, S. M. 2016, Phys. Rev. Lett., 116, 214501 [NASA ADS] [CrossRef] [Google Scholar]
 Mishra, P. K., & Verma, M. K. 2010, Phys. Rev. E, 81, 056316 [NASA ADS] [CrossRef] [Google Scholar]
 Moll, R., Graham, J. P., Pratt, J., et al. 2011, ApJ, 736, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Novi, L., von Hardenberg, J., Hughes, D. W., Provenzale, A., & Spiegel, E. A. 2019, Phys. Rev. E, 99, 053116 [NASA ADS] [CrossRef] [Google Scholar]
 Nozawa, T., & Yoden, S. 1997a, Phys. Fluids, 9, 2081 [NASA ADS] [CrossRef] [Google Scholar]
 Nozawa, T., & Yoden, S. 1997b, Phys. Fluids, 9, 3834 [NASA ADS] [CrossRef] [Google Scholar]
 Oliphant, T. E. 2006, A Guide to NumPy, 1 (Trelgol Publishing USA) [Google Scholar]
 Parisi, M., Kaspi, Y., Galanti, E., et al. 2021, Science, 374, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Pelinovsky, E. N. 1978, Okeanologiya, 18, 192 [Google Scholar]
 Read, P. L., Young, R. M. B., & Kennedy, D. 2020, Geosci. Lett., 7, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Reshetnyak, M. Y. 2013, Izvestiya, Phys. Solid Earth, 49, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Rhines, P. B. 1975, J. Fluid Mech., 69, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, L. F. 1922, Weather Prediction by Numerical Process (Cambridge University Press) [Google Scholar]
 Roberts, P. H. 1968, Philos. Trans. Roy. Soc. A, 263, 93 [Google Scholar]
 Rubio, A. M., Julien, K., Knobloch, E., & Weiss, J. B. 2014, Phys. Rev. Lett., 112, 144501 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1989 Differential rotation and stellar convection. Sun and the solar stars (Berlin: Akademie Verlag) [CrossRef] [Google Scholar]
 Salyk, C., Ingersoll, A. P., Lorre, J., Vasavada, A., & Del Genio, A. D. 2006, Icarus, 185, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Schaeffer, N. 2013, Geochem. Geophys. Geosyst., 14, 751 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, T., & Liu, J. 2009, J. Atmos. Sci., 66, 579 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P. 2007, J. Atmos. Sci., 64, 3132 [NASA ADS] [CrossRef] [Google Scholar]
 Siegelman, L., Klein, P., Ingersoll, A. P., et al. 2022, Nat. Phys., 18, 357 [CrossRef] [Google Scholar]
 Srinivasan, K., & Young, W. R. 2012, J. Atmos. Sci., 69, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Sukoriansky, S., Dikovskaya, N., & Galperin, B. 2007, J. Atmos. Sci., 64, 3312 [NASA ADS] [CrossRef] [Google Scholar]
 Tobias, S. M., Dagon, K., & Marston, J. B. 2011, ApJ, 727, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Vallis, G. K., & Maltrud, M. E. 1993, J. Phys. Oceanogr., 23, 1346 [NASA ADS] [CrossRef] [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Vasavada, A. R., & Showman, A. P. 2005, Rep. Progr. Phys., 68, 1935 [CrossRef] [Google Scholar]
 Verma, M. K. 2019, Physica Scripta, 94, 064003 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, M. K., Ayyer, A., Debliquy, O., Kumar, S., & Chandra, A. V. 2005, Pramana, 65, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, M. K., Kumar, A., & Pandey, A. 2017, New J. Phys., 19, 025012 [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 [CrossRef] [Google Scholar]
 Wicht, J. 2002a, Phys. Earth Planet. Interiors, 132, 281 [CrossRef] [Google Scholar]
 Wicht, J. 2002b, Icarus, 155, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Wicht, J., & Gastine, T. 2020, Nat. Commun., 11, 2886 [NASA ADS] [CrossRef] [Google Scholar]
 Wicht, J., Dietrich, W., Wulff, P., & Christensen, U. R. 2020, MNRAS, 492, 3364 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, G. P. 1978, J. Atmos. Sci., 35, 1399 [NASA ADS] [CrossRef] [Google Scholar]
 Xiao, Z., Wan, M., Chen, S., & Eyink, G. L. 2009, J. Fluid Mech., 619, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., & Bloxham, J. 2020, Proc. Natl. Acad. Sci. U.S.A., 117, 13991 [NASA ADS] [CrossRef] [Google Scholar]
 Young, R. M. B., & Read, P. L. 2017, Nat. Phys., 13, 1135 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Schematic illustration of two possible mechanisms for jet driving. Top row (a): the inverse cascade idea, where eddies merge until the jet scale is reached. The first step in this process illustrates the required anisotropy necessary to prefer either cyclonic or anticyclonic eddies, here illustrated by a sorting of both eddy types. Bottom row (b): the direct driving via Reynolds stresses, where an initial tilt of convective features (second column) yields weak zonal flows that then amplify the tilt until a viscous balance is reached. 

In the text 
Fig. 2 Typical simulation snapshots in the statistically steady state. Panel (a) shows the jet profile, . Panels (b) and (c) show the fluctuations around the zonal mean of the zonal and meridional winds, and , at the central meridian. Panel (d) shows at the top boundary. Panels (b)–(d) are saturated at half of the maximum value in the domain shown. 

In the text 
Fig. 3 Kinetic energy spectra timeaveraged over an interval that is part of the statistically steady state of our simulation as a function of harmonic degree (a) and azimuthal order (b). Both panels show separate spectra for the toroidal and poloidal flow components. See Sect. 2.2 for a brief discussion of the possible power laws shown in panel (a). 

In the text 
Fig. 4 Square root of the shellaveraged squared forces in spectral space as a function of harmonic degree ℓ (also known as rms force spectra). The Coriolis force is dominating at all scales as is typical for rapidly rotating sphericalshell convection. 

In the text 
Fig. 5 Scalebyscale balance (a) and scaletoscale flux (b) of kinetic energy in the statistically steady state. The values T(0) ≈ D(0) ≈ 7 × 10^{9} were excluded from the plotting range in panel (a) for better visibility. 

In the text 
Fig. 6 Scaletoscale transfer of kinetic energy, T (m, m′), in the statistically steady state. Note that the figure has logarithmic axes and two colorbars; the left one indicates values for m = 0 or m′ = 0 and the right one for m, m′ > 0. A positive value of T(m, m′) shows that kinetic energy is transferred from wave number m′ to wave number m. We use the regions indicated by gray lines and labeled with letters (a) to (e) to refer to specific regions of the plot in the text and show the different spatial contributions to T(m, m′) for those regions in Fig. 7. 

In the text 
Fig. 7 Spatial contribution, T_{i}(r, θ), to the energy transfer function, T(m, m′), summed over different values of m and m’ according to the masks (i) = (a),(b),(c),(d),(e) shown in Fig. 6. The panels display results for the direct driving of the jets by the convective scales via upscale energy transfers and Reynolds stresses (a), the upscale transfer to large eddies (b), the feeding of largescale vortices from the jet energy (c), the smallscale forward cascade (d), and the local forward transfer likely associated with a forward enstrophy transfer (e). All panels are saturated at three times their maximum value, except for panel (a), which is saturated at 18 times its maximum value. 

In the text 
Fig. 8 Energy transfer function, T(m, m′), in the statistically steady state, integrated over the region inside the TC (s < r_{i}). 

In the text 
Fig. 9 Spatial dependence of the contribution to the kinetic energy in the statistically steady state, E_{kin}, for a sum over azimuthal wave numbers m = 3,…, 8 (see label (c) in Fig. 6). See Eq. (21) for a definition of E_{kin}(m, r, θ). 

In the text 
Fig. 10 Snapshots of the statistically steady state at the shell center (r/r_{o} = 0.85). The fluctuations around the meridional flow, , clearly show structures of m ≈ 5 at higher latitudes. 

In the text 
Fig. 11 Spatial dependence of the driving of kinetic energy of largescale eddies by buoyancy, F_{Buoy,c} (a), and kinetic energy transfer to the same eddies from other eddies, T_{c,m′>0} (b), both in the statistically steady state. For both panels, we summed T and F_{Buoy} over m as given by label (c) in Fig. 6 and summed T over m′ > 0. 

In the text 
Fig. 12 Temporal development of the jets during spinup. We show the kinetic energy of zonal flows (a) and poloidal flows (b, mostly convection), shellaveraged rms Reynolds stresses (c and d), and the surface jet profile (e). See Eqs. (35) and (36) for definitions of the spatial averages of the Reynolds stresses and correlation coefficients. This run was performed using a convective initial condition without jets. The duration of a convective turnover time is about 10^{−3} τ_{visc}. We show results for the spatially resolved Reynolds stresses and the energy transfer to the jets during the intervals marked as gray bands in Figs. 13 and 14. 

In the text 
Fig. 13 Reynolds stresses R_{s,φ} at different times during the spinup: at the initial condition without jets and uniform background rotation (a), after about 100 turnover times and in the growth phase of the jets (b, same as interval III in panel (e) in Fig. 12) and in the statistically steady state with fully developed jets (c). 

In the text 
Fig. 14 Kinetic energy transfer to the jets, T(0, m′), for the spinup run. The transfer functions were computed for the first interval directly after the start of the spinup (interval I). The time window is shown as light gray band in Fig. 12. For comparison, we also show results for the statistically steady state from Fig. 6. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.