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

# A laboratory model for deep-seated jets on the gas giants

Simon Cabanes, Jonathan M. Aurnou, Benjamin Favier, Michael Le Bars

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

## To cite this version:

Simon Cabanes, Jonathan M. Aurnou, Benjamin Favier, Michael Le Bars. A laboratory model for deep-seated jets on the gas giants. Nature Physics, 2017, 13, 10.1038/nphys4001. hal-01478957

HAL Id: hal-01478957

https://hal.science/hal-01478957

Submitted on 30 Oct 2018

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_3.png" /></figure>
# A laboratory model for deep-seated jets on the gas giants

Simon Cabanes<fnref n="1" />*, Jonathan Aurnou<fnref n="1" />,², Benjamin Favier<fnref n="1" /> and Michael Le Bars<fnref n="1" />

The strong east-west jet flows on the gas giants, Jupiter<fnref n="1" /> and Saturn², have persisted for hundreds of years. Yet, experimental studies cannot reach the planetary regime and similarly strong and quasi-steady jets have been reproduced in numerical models only under simplifying assumptions and limitations. Two models have been proposed: a shallow model where jets are confined to the weather layer and a deep model where the jets extend into the planetary molecular envelope. Here we show that turbulent laboratory flows naturally generate multiple, alternating jets in a rapidly rotating cylindrical container. The observed properties of gas giants' jets are only now reproduced in a laboratory experiment emulating the deep model. Our findings demonstrate that long-lived jets can persist at high latitudes even under conditions including viscous dissipation and friction and bear relevance to the shallow versus deep models debate in the context of the ongoing Juno mission³.

Azimuthally directed (that is, zonal, east-west) jet flows are one of the dominant characteristics in the surficial cloud features observed on the gas giants, Jupiter and Saturn. An essential question of planetary dynamics and structure is whether these jet motions exist only within the shallow troposphere or extend through the molecular envelope that exists above the deeper dynamo region³. Determining the depth of these atmospheric jets is one of the prime directives of the NASA (National Aeronautics and Space Administration) Juno mission, which entered into low-altitude Jovian orbit in August 2016⁴. Despite the long-lived scientific interest in these flows, dominant multiple jets have been problematic in fully three-dimensional (3D) numerical models of convection. In particular, multiple banded flows are not found in the most recent, high-resolution models that couple the molecular envelope to the deeper dynamo region. In these models, magnetic dissipation damps the higher-latitude deep jets out of existence⁵⁻⁷. Similarly, dissipation has also proved overly important in laboratory experiments carried out to date. Laboratory approaches were analysed in the framework of the shallow-layer model and strong viscous damping by the container boundaries only allows for the formation of weak zonal jets with tenuous instantaneous signatures⁸⁻¹⁴. Thus, it has yet to be demonstrated, as proposed for the gas giant planets¹⁵, that deep zonally dominant jet flows can exist in the presence of boundary dissipation.

We have developed a new laboratory experimental device that is capable of generating strong zonal jets despite viscous friction on the boundaries (Fig. 1a). The working fluid is water, contained in a 1.37-m-high by 1-m-diameter cylindrical tank. The depth of the fluid layer is $h_o = 50$ cm at rest, and the tank's rotation rate is $\Omega = 7.85 \text{ rad s}^{-1}$ (75 revolutions per minute). Once equilibrated at $\Omega$, the water's free surface takes the shape of a paraboloid, with the fluid layer depth ranging from $h_{\min} \approx 20$ cm on the axis of

rotation to $h_{\max} \approx 90$ cm at the tank's outer radius. This rotating surface shape is analogous to the large-scale curvature of a deep spherical planetary fluid layer¹⁶,¹⁷. In addition, the rotation provides strong Coriolis forces, as exist in planetary settings. Once solid body rotation is reached, a submersible pump situated at the base of the tank is turned on, and small-scale turbulence is injected at the base of the fluid layer. The pump continuously circulates water through a lattice of 32 outlets (4 mm diameter) and 32 inlets (2 mm diameter) arranged on a flat base plate without any axisymmetric features (see the injection pattern in Supplementary Fig. 1a). Typical root-mean-square (r.m.s.) fluctuating velocities are in the range $u_{r.m.s.} \approx 1-5 \text{ cm s}^{-1}$. This small-scale turbulence is analogous to the convective turbulence that exists in deep planetary interiors¹⁸,¹⁹ and constitutes an appropriate source of energy to supply deep jets²⁰. Inspired by several preceding experimental attempts¹⁰⁻¹⁴,²¹,²² this laboratory device bridges the gap to the planetary 'zonostrophic' regime and suggests the possibility of deep-seated jets on the gas planets.

In natural settings, zonal jets develop from rapidly rotating turbulence in the presence of strong boundary curvature²³. In particular, the Coriolis force must dominate the fluid's inertia (that is, low Rossby number $Ro$), which, in turn, must dominate the viscous dissipative effects (large Reynolds number $Re$ and small Ekman number $E$) (see Supplementary Table 1). Such flows, because they are essentially two-dimensional (2D), transfer energy upscale via a so-called inverse cascade²⁴,²⁵ driven by Reynolds stresses¹⁴. In addition, the spherical geometry of planetary fluid layers is ideal for the formation of zonal jet flows: it induces a zonal anisotropization of the inverse energy cascade following the so-called $\beta$-effect. Energy is consequently channelled into zonal jets with characteristic latitudinal width approximated by the Rhines scale $L_{Rh} = \pi\sqrt{2U_{r.m.s.}^T/\beta}$, where $U_{r.m.s.}^T$ is the total r.m.s. velocity used as the characteristic system scale flow velocity²³,²⁵. The spherical curvature of the fluid shell, expressed via $\beta$, acts to halt the 2D energy cascade, leading to the formation of Rhines scale zonal jets.

All these physical ingredients can be reproduced following two distinct scenarios, depending on the nature of the anisotropic $\beta$ parameter and on the jets dynamical aspect ratio $\eta = h/L_{Rh}$ on which jets develop (with $h$ and $L_{Rh}$ being the vertical and horizontal scales of jets). In the first scenario, jets are confined into a shallow troposphere where shallow turbulence is essentially 2D because of depth confinement ($\eta \ll 1$). The zonal anisotropization of the inverse cascade is due to the variation of the Coriolis force with latitude $\theta$ (the planetary $\beta$-effect: $\beta = 2\Omega \cos\theta/r$, with $r$ being the planetary radius) and the relevant axis of rotation is the local normal to the planetary surface. In the second scenario, jets are deep-seated and extend into a 3D convective deep layer where rotation dynamically constrains large-scale motions to be invariant along the

<footnote n="1">¹CNRS, Aix Marseille University, Centrale Marseille, IRPHE, 13013 Marseille, France. ²Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, California 90095, USA. *e-mail: cabanes.simon@gmail.com</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1** | Experimental and numerical systems. **a**, Schematic of the experimental device. The tank is 1.37 m high by $R=0.5$ m in radius and is spun at up to 75 revolutions per minute ($\Omega=7.85 \text{ rad s}^{-1}$) with $\simeq 400 \text{ l of water such that } R_0=7 \times 10^{-3}, R_e=3 \times 10^4 (E=2.3 \times 10^{-7})$ and $\eta=5.15$, with a mean zonostrophy index $R_\beta \simeq 3.7$ (see Methods for more details). The free surface of the fluid has a depth $h(r)$ ranging from $h_{\text{min}}=20$ cm to $h_{\text{max}}=90$ cm with a depth at rest $h_0=50$ cm. The non-axisymmetric injection pattern is shown in Supplementary Fig. 1a. **b**, Direct numerical simulation of the experiment. Sidecut snapshot at $t \simeq 1,000$ rotations showing the instantaneous non-dimensional azimuthal velocity field at $R_0=10^{-3}, R_e=10^3 (E=10^{-6})$, corresponding to a numerical mean zonostrophy index $R_\beta \simeq 3.4$.</figcaption></figure>

rotation axis ($\eta \gtrsim 1$). In such a flow, dynamical confinement leads to bidimensionality, the resulting inverse cascade is dependent on the rotation rate and 3D wave motions may play an important role in the dynamics²². In deep spherical fluid layers the zonal anisotropy arises due to cylindrically radial variations in the depth of the fluid layer²³ (the topographic $\beta$-effect: $\beta=2\Omega h^{-1} dh/dr$) and the relevant axis of rotation is the planetary spin axis. The distinction between shallow and deep-seated jets has never been clear in previous laboratory devices with moderate aspect ratio ($0.57 \le \eta \le 1$). Mid-latitude Jovian jet aspect ratios are $\eta \simeq 5$ and $\eta \simeq 0.004$ for the deep and shallow scenarios, respectively. In our laboratory device we emulate a gas planet deep-layer scenario with $\eta=5.15$ (see Supplementary Table 1 for comparative non-dimensional analysis).

Atmospheric zonal flows are detected in all four Solar System gas planets (see for example, Fig. 2a for Jupiter), and, recently, zonal jet structures have also been observed in Earth's oceans²⁶. However, gas giant jets are known to be incredibly steady over time whereas oceanic flows are strongly fluctuating and zonal features are revealed only following careful time-averaging of the flow fields. These differing regimes of jet flow behaviours can be characterized in terms of the zonostrophy index²⁷ $R_\beta$ (see Methods for details on its calculation). This index is the ratio of the Rhines scale and the transitional scale $L_\beta \simeq (\epsilon/\beta^3)^{1/5}$ above which the $\beta$-effect affects the energy cascade, with $\epsilon$ being the inverse energy transfer rate through the system²³. Simplistically, $R_\beta = L_{Rh}/L_\beta$ represents the strength of the zonal flows relative to the strength of the underlying mixing. If mixing effects are strong (large $L_\beta$), the zonal jets tend to be undulatory and unsteady in time, or fail to form altogether²⁸. As $L_\beta$ is decreased, at a given value of $L_{Rh}$, quasi-steady zonal jets are predicted to become the prevalent fluid motion in the system²³,²⁸. Such flows, argued to exist at $R_\beta \gtrsim 2.5$, are in the so-called zonostrophic regime²⁷. Deep Jovian mid-latitude jets exceed this threshold²⁹ with a typical value $R_\beta \simeq 5$, while Earth's oceans

are limited to $R_\beta \simeq 1.5$. In previous laboratory experiments, strong boundary dissipation and relatively low $\beta$ led to large $L_\beta$, non-zonostrophic flows (see Supplementary Table 1 for comparative non-dimensional analysis).

In our experimental device, we create for the first time in the laboratory, the appropriate conditions for generating gas planet-like deep-seated zonal flows (that is, laboratory non-dimensional parameters $R_0 = 7 \times 10^{-3}$, $R_e = 3 \times 10^4$, $E = R_0/R = 2.3 \times 10^{-7}$, mean aspect ratio $\eta = 5.15$ and large values of the mean zonostrophy index $R_\beta \simeq 3.7$ that exist in deep zonostrophic regime; see Methods and Supplementary Table 1). Experimental measurements of the surface velocities are acquired by recording the tracks of small floating particles from a top lid camera onboard the rotating frame and analysing their paths using a Lagrangian tracking method (see Methods). To validate our laboratory surface velocity measurements, we have carried out high-resolution direct numerical simulations of the experiment. Using the same upper fluid surface shape as in the 75 r.p.m. experiments, but at slightly higher Ekman number, $E = 10^{-6}$, the simulated flow is forced at the base of the tank via an array of vortices that approximate the small-scale basal pumping in the experiment, starting from the fluid initially at rest in the rotating frame (see Methods and Supplementary Fig. 1b). Figure 1b shows an instantaneous sidecut rendering of the azimuthal velocity field after 1,000 simulated rotation times. The large-scale surface zonal flows are deep-seated, validating our laboratory surface velocimetry and confirming that the flow is dynamically quasi-2D due to the dynamical effect of rotation ($R_0 \ll 1$). Further, it reveals that locally injected turbulence can drive large-scale zonal jets, even in the presence of boundary dissipation (all boundaries except the upper paraboloid being no-slip, see Supplementary Fig. 2c,d).

Our laboratory experiment produces six zonostrophic jets with strong, instantaneous signatures, similar to those on gas planets (Fig. 2). In Fig. 2b, we show experimental particle tracks acquired from 4,900 to 5,350 rotations. The tracks are coloured on the basis of their instantaneous zonal velocity, red (blue) being prograde (retrograde). The zonally elongated tracks show that the jets exist not only in time average, but instantaneously at all times (see also Supplementary Movie 1). This instantaneous signature is novel in a rapidly rotating turbulent jets experiment and has not yet been attained in numerical models of deep jet generation with boundary dissipation. Compiling all the velocity data in Fig. 2b yields the time-averaged zonal velocity field map shown in Fig. 2c, revealing strong zonation. We calculate the surface average of the zonostrophy index in our experiment (that is, with an estimate of $\epsilon$ using the total horizontal kinetic energy, see equation (6) in Methods), to be $R_\beta = 3.7$, confirming the relevance of our experiment to the gas giants zonostrophic regime reached for $R_\beta \ge 2.5$.

Figure 3a shows the time- and azimuthal-averaged zonal velocity profile. The solid red line denotes the spatiotemporal-averaged jet flows as a function of the cylindrical radius, whereas the r.m.s. values of the velocity fluctuations are demarcated by the dashed blue line. Using kinetic energy spectrum analysis in our deep model we retrieve in Fig. 3b the theoretical²⁵ anisotropic zonal flow spectrum $E_Z(k) = C_Z \beta^2 k^{-5}$ and the 2D residual (non-zonal) energy spectrum from velocity fluctuations $E_R(k) = C_R \epsilon^{2/3} k^{-5/3}$ (where $k$ is the 2D horizontal wavenumber and $C_R=5$, see Methods). Considering our radial evolution of $\beta$ we obtain that the zonal constant $C_Z$, which best fits our zonal spectrum, ranges from 1.7 to 3.7. This constitutes the first experimental measurement that confirms the estimate from direct observations of Jupiter's clouds that gives $C_Z \simeq 2$ (ref. 29). However, since both of these studies are based on a constant $\beta$ approximation, one cannot discriminate between shallow and deep $\beta$-effects for Jupiter using this spectral analysis. Similarly to the Jupiter observations, we report that energy is maximum at jets scale (that is, at the Rhines scale $k_{Rh} = 2\pi/L_{Rh}$) where zonal kinetic energy strongly dominates over the residual fluctuating flow. Transition
<figure><img src="image_5.png" /></figure>

**Figure 2** | Polar view and zonal velocities. **a,** Southern polar view of Jupiter's atmosphere from the Cassini spacecraft mission (image credit: NASA/JPL/Space Science Institute). **b,** Experimental polar view of 4,900 Lagrangian particle tracks, collected between 4,900 and 5,350 rotations. Tracks are coloured on the basis of their instantaneous zonal azimuthal velocity direction, with red (blue) being prograde (retrograde). The radial velocity components of the tracks are reported in Supplementary Fig. 2b. **c,** Experimental polar view of the time-averaged zonal velocity field computed from 4,900 Lagrangian tracks. Acquisition and time-averaging procedure are detailed in the Methods.

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

**Figure 3** | Velocity profiles and kinetic energy spectra. **a,** Radial dependence of the mean zonal velocity $\langle U_\psi \rangle / \Omega R$ (thick red line) and the r.m.s. velocity averaged in azimuth $\bar{u}_{r,\text{m.s.}} / \Omega R$ (dashed blue line) computed from the particle tracks shown in Fig. 2c. **b,** Experimental (red solid line) and numeric (blue solid line) zonal kinetic energy spectra $E_Z(k)$ as a function of total wavenumber $k = \sqrt{k_r^2 + k_\psi^2}$, with $k_r$ and $k_\psi$ the radial and azimuthal wavenumbers. The green solid line is the 2D spectrum of the horizontal residual flow $E_R(k)$ obtained from the numerical simulation. In all cases, the spectra are in good agreement with their theoretical predictions, reported as dashed and dotted black lines, using the mean value of $\beta = 73 \text{ m}^{-1}\text{s}^{-1}$ and $\epsilon = 6.3 \times 10^{-7} \text{ m}^2\text{s}^{-3}$ based on the Ekman layer dissipation. Black vertical lines report the Rhines scale $k_{\text{Rh}} = \sqrt{\beta/2U_{r,\text{m.s.}}^\text{T}}$ and the transitional wavenumber $k_\beta = 0.73(\beta^3/\epsilon)^{1/5}$, which, for a mean $\beta$ value, leads to $R_\beta = 3.7$ in the experiment. The vertical grey area shows an estimate of the scale range of kinetic energy injection, with a lower limit corresponding to numerical forcing and an upper limit corresponding to the estimated experimental injection scale (that is, half of the lattice grid size, 5 cm). Details of the calculation method are given in the Methods.

to a strongly anisotropic flow where energy becomes preferentially channelled into the zonal direction is well observed at the predicted transitional scale $k_\beta = 1/L_\beta$. We find that more than 50% of the total surface kinetic energy is concentrated in the zonal velocity component, and this value goes up to 70% in the sharpest prograde jet (see Methods for more details).

Focusing on the first 1,000 rotations of another 75 r.p.m. experiment, Fig. 4 shows the evolving width of the jets over time, with merging events occurring typically within the first 500 to 600 rotations. After this, strongly energetic, long-lived zonal jets reach a quasi-steady state (see Supplementary Fig. 3a for complementary results). The steadiness of the multiple jets system is an important feature shared by gas planets as well as our laboratory and numerical models. Even in the more dissipative no-slip numerical simulation (where $E = 10^{-6}$), dissipation does not prevent the formation of intense jets (see Supplementary Fig. 3b). Comparison between the experimental and numerical cases, including the simulation

with stress-free boundary conditions (see Supplementary Fig. 3c), highlights that the boundary dissipation is a key ingredient in determining the equilibrated scale and strength of the deep-seated jets, as well as the typical time for reaching a steady-state regime.

Our experiments demonstrate for the first time that robust, quasi-steady zonostrophic, deep ($\eta \ge 1$) jets are realizable in a relatively simple laboratory-scale device that brings together the physical ingredients of planetary zonal flow systems: dynamically constraining rapid rotation, highly turbulent flows, and large fluid layer depth variations. Zonal jet formation then occurs even in the presence of bottom boundary friction, in sharp contrast to the results of overly viscous 3D simulations in which the mid- to high-latitude jets are damped out$^{5,30}$. We further claim that dissipation is an essential ingredient for jet equilibration$^{15}$. On the basis of our findings, we hypothesize that deep planetary jets will also form in the presence of a magnetic dissipation region, as exists within the gas planets. Indeed, a magneto-dissipative layer localized at the
<figure><img src="image_7.png" /></figure>

**Figure 4** | Temporal evolution of the zonal velocity. Azimuthal and time-windowed average of laboratory zonal velocities, $\langle U_\psi \rangle / \Omega R$, acquired within the first 1,000 rotations. At time $t=0$, the fluid is at rest in the rotating frame and the pump is primed. The time average is computed over a time window of 50 rotations sliding by steps of 3 rotations. The calculation method is given in the Methods.

transition between conductive and non-conductive regions will not be likely to impact the ability for geophysical turbulence to generate deep zonostrophic zonal jets extending across the gas giants' molecular envelopes. Ultimately, our results suggest that the forthcoming Juno and Cassini measurements will support the existence of deep-seated zonal jets on the gas giants.

## Methods

Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper.

Received 2 May 2016; accepted 28 November 2016; published online 23 January 2017

## References

1. Bagenal, F., Dowling, T. E. & McKinnon, W. B. *Jupiter: The Planet, Satellites and Magnetosphere* Vol. 1 (Cambridge Univ. Press, 2006).

2. Del Genio, A. D. et al. *Saturn from Cassini-Huygens* 113–159 (Springer, 2009).

3. Ingersoll, A. P. Atmospheric dynamics of the outer planets. *Science* **248**, 308–315 (1990).

4. Bolton, S. The Juno mission. *Proc. Int. Astron. Union* **6**, 92–100 (2010).

5. Heimpel, M. & Gómez Pérez, N. On the relationship between zonal jets and dynamo action in giant planets. *Geophys. Res. Lett.* **38**, L14201 (2011).

6. Gastine, T., Heimpel, M. & Wicht, J. Zonal flow scaling in rapidly-rotating compressible convection. *Phys. Earth Planet. Inter.* **232**, 36–50 (2014).

7. Jones, C. A dynamo model of Jupiter's magnetic field. *Icarus* **241**, 148–159 (2014).

8. Colin De Verdiere, A. Mean flow generation by topographic Rossby waves. *J. Fluid Mech.* **94**, 39–64 (1979).

9. Condie, S. A. & Rhines, P. B. A convective model for the zonal jets in the atmospheres of Jupiter and Saturn. *Nature* **367**, 711–713 (1994).

10. Afanasyev, Y. D. & Wells, J. Quasi-two-dimensional turbulence on the polar beta-plane: laboratory experiments. *Geophys. Astrophys. Fluid Dyn.* **99**, 1–17 (2005).

11. Di Nitto, G., Espa, S. & Cenedese, A. Simulating zonation in geophysical flows by laboratory experiments. *Phys. Fluids* **25**, 086602 (2013).

12. Smith, C. A., Speer, K. G. & Griffiths, R. W. Multiple zonal jets in a differentially heated rotating annulus. *J. Phys. Oceanogr.* **44**, 2273–2291 (2014).

13. Zhang, Y. & Afanasyev, Y. Beta-plane turbulence: experiments with altimetry. *Phys. Fluids* **26**, 026602 (2014).

14. Read, P. et al. An experimental study of multiple zonal jet formation in rotating, thermally driven convective flows on a topographic beta-plane. *Phys. Fluids* **27**, 085111 (2015).

15. Liu, J. & Schneider, T. Mechanisms of jet formation on the giant planets. *J. Atmos. Sci.* **67**, 3652–3672 (2010).

16. Schopp, R. & Colin De Verdiere, A. Taylor columns between concentric spheres. *Geophys. Astrophys. Fluid Dyn.* **86**, 43–73 (1997).

17. Heimpel, M., Aurnou, J. & Wicht, J. Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model. *Nature* **438**, 193–196 (2005).

18. Ingersoll, A. P. & Pollard, D. Motion in the interiors and atmospheres of Jupiter and Saturn: scale analysis, anelastic equations, barotropic stability criterion. *Icarus* **52**, 62–80 (1982).

19. Aurnou, J. et al. Rotating convective turbulence in Earth and planetary cores. *Phys. Earth Planet. Inter.* **246**, 52–71 (2015).

20. Showman, A. P., Gierasch, P. J. & Lian, Y. Deep zonal winds can result from shallow driving in a giant-planet atmosphere. *Icarus* **182**, 513–526 (2006).

21. Solomon, T., Holloway, W. & Swinney, H. L. Shear flow instabilities and Rossby waves in barotropic flow in a rotating annulus. *Phys. Fluids A* **5**, 1971–1982 (1993).

22. Yarom, E. & Sharon, E. Experimental observation of steady inertial wave turbulence in deep rotating flows. *Nat. Phys.* **10**, 510–514 (2014).

23. Vallis, G. K. & Maltrud, M. E. Generation of mean flows and jets on a beta plane and over topography. *J. Phys. Oceanogr.* **23**, 1346–1362 (1993).

24. Kraichnan, R. H. Inertial ranges in two dimensional turbulence. *Phys. Fluids* **10**, 1417–1423 (1967).

25. Rhines, P. B. Waves and turbulence on a beta-plane. *J. Fluid Mech.* **69**, 417–443 (1975).

26. Maximenko, N. A., Bang, B. & Sasaki, H. Observational evidence of alternating zonal jets in the world ocean. *Geophys. Res. Lett.* **32**, L12607 (2005).

27. Sukoriansky, S., Dikovskaya, N. & Galperin, B. On the arrest of inverse energy cascade and the Rhines scale. *J. Atmos. Sci.* **64**, 3312–3327 (2007).

28. Scott, R. K. & Dritschel, D. G. The structure of zonal jets in geostrophic turbulence. *J. Fluid Mech.* **711**, 576–598 (2012).

29. Galperin, B. et al. Cassini observations reveal a regime of zonostrophic macroturbulence on Jupiter. *Icarus* **229**, 295–320 (2014).

30. Jones, C. A. & Kuzanyan, K. M. Compressible convection in the deep atmospheres of giant planets. *Icarus* **204**, 227–238 (2009).

## Acknowledgements

This work was supported by the French Agence Nationale pour la Recherche (Grant No. ANR-13-JS05-0004-01) and by the Aix-Marseille University Foundation (A*MIDEX project No. ANR-11-IDEX-0001-02). J.A.'s sabbatical in Marseille was supported by S. A. Kremen, by the invited professor programmes of Ecole Centrale Marseille and Aix-Marseille University, and by the Labex MEC (grant ANR-11-LABX-0092). We also acknowledge support from IDRIS for computational time on Turing (Project No. 100508) and from the HPC resources of Aix-Marseille University (project Equip@Meso No. ANR-10- EQPX-29-01 of the programme Investissements d'Avenir).

## Author contributions

This work is based on an initial design concept by J.A. All authors contributed to the final design and the building of the experimental set-up. S.C. ran the experiment and carried out data processing. B.F. performed the numerical simulations. All authors contributed to data analysis and to writing the letter, with a major contribution from S.C.

## Competing financial interests

The authors declare no competing financial interests.
## Methods

**Experimental device.** Our device is based on several previous set-ups$^{10–14,21,22}$. It consists of a 1-m-diameter by 1.37-m-high cylindrical tank of polythene that is attached to the rotating table by a rigid aluminium super-structure. A 32-cm-high bottom table, made of PVC, encloses our basal injection system. We drilled 64 holes of, respectively, 2 mm and 4 mm diameter connected to 32 injection and 32 suction tubes distributed following a non-axisymmetric square grid (see Supplementary Fig. 1a). This small-scale turbulence injection system is used to mimic convective turbulence in deep planetary interiors, in the manner of previous experimental devices$^{8,22,31,32}$. These tubes are connected to an immersed AQUASON-506-M-AL pump giving an operating flow of $2 < Q < 5 \text{ m}^3/\text{h}$. Consequently, the injected flow at the bottom of the domain is mainly three-dimensional (3D) and weakly affected by rotation (injection velocity $1 < \bar{u} < 3 \text{ m s}^{-1}$ and Rossby number $0.1 < Ro < 0.4$). An acrylic top lid carries a GoPro HERO3+ camera that records the surface flow at 59 frames per second. A lighting apparatus, including a powerful light-emitting diode lamp and a diffuser, sits on top of the aluminium super-structure. The whole device is attached to an ABRT 7686 rotary table using robust high-precision axial roller bearings with an air bearing in the horizontal axis for smooth, precise rotary motion. The table is 1.4 m in diameter with maximum loading of 1,000 kg and a typical operating speed of 1–90 r.p.m.

**Experimental data acquisition.** Assuming the flow to be quasi-geostrophic (that is, nearly invariant along the rotation axis, see Fig. 1b), we perform free-surface visualization. Lagrangian surface velocities are measured using a particles tracking code$^{33}$. We follow the surface motions of 90 black floating particles taken to be smaller than the dynamical flow scales (3 mm diameter spherical particles) and to be at the same density as the carrier fluid (for better optical contrast, the water is tinted white via a titanium dioxide suspension).

After camera calibration and correction of parabolic free-surface deformation, we obtain zonal and radial flow velocity components along Lagrangian tracks. We record the floaters' motions over 5,000 rotations. A quasi-steady state is achieved within 600 rotations (see Fig. 4 and Supplementary Fig. 3). We define a Cartesian grid of $142 \times 142$ square elements of 7 mm each. Temporal mean zonal velocity $\langle U_\phi \rangle$ and azimuthal-averaged zonal velocity $\langle U_\psi \rangle$ are obtained in each grid element from time-averaging over 450 rotations once steady state is achieved (see Figs 2c and 3a). Figure 4 and Supplementary Fig. 3 represent azimuthal and time-windowed averages (the time average is computed over 50 rotations and sliding by steps of 3 turns) from the initial time at which we activate the pump. Velocity fluctuations reported in Fig. 3a are obtained by computing the r.m.s. velocity $u_{\phi \text{, m.s.}}^q = \sqrt{1/N \sum_{i=1}^N [(U_\phi)_{li} - U_\phi^q]^2}$, where N is the total number of tracks in the i'th grid elements and $U_{\psi \text{, n.i.}}^q$ is the zonal velocity of the n'th track. By using the same procedure for the radial component of the flow, we obtain the radial mean velocity $\langle U_\psi \rangle$ and radial fluctuations $u_{\psi \text{, m.s.}}^r$. We report that the radial mean flow is nearly zero and radial fluctuations are of the same order as zonal fluctuations $u_{\psi \text{, m.s.}}^q \approx u_{\psi \text{, m.s.}}^r$. The total kinetic energy corresponds to the sum $E_T = \langle U_\psi \rangle^2 + (u_{\psi \text{, m.s.}}^q)^2 + (u_{\psi \text{, m.s.}}^r)^2 + (u_{\psi \text{, m.s.}}^r)^2$ and $E_{\text{fluct.}} = (u_{\psi \text{, m.s.}}^q)^2 + (u_{\psi \text{, m.s.}}^r)^2$ is the energy of the fluctuating velocity. Within the jets, at the maximum of each zonal velocity peak, the kinetic energy of the fluctuations reaches 43% of the total kinetic energy ($0.3 < E_{\text{fluct.}}/E_T < 0.43$), while the zonal kinetic energy exceeds 57% of the total ($0.57 < E_{\psi}/E_T < 0.7$).

**Numerical method.** The Navier–Stokes equations are solved in their weak variational form$^{34}$ for an incompressible fluid in three dimensions with the spectral element code Nek5000 (http://nek5000.mcs.anl.gov). The spectral element method combines the geometric flexibility of finite-element methods with the accuracy of spectral methods. It is particularly well adapted to our problem involving turbulent flows in complex geometries. Nek5000 has, for example, already been used in the context of rotating flows in ellipsoidal containers$^{35}$. The computational domain is decomposed into non-overlapping hexahedral elements. Within each element, the velocity is decomposed onto Lagrange polynomial interpolants of order N, based on the tensor-product of Gauss–Lobatto–Legendre quadrature points. For all the simulations discussed in this paper, numerical convergence was checked by fixing the number of elements and increasing the degree N of the polynomial decomposition. The temporal discretization is based on a semi-implicit formulation in which the nonlinear and rotation terms are treated explicitly and all remaining linear terms are treated implicitly. Note that our solution is dealiased following the 3/2 rule. The code is efficiently parallelized using MPI (Message Passing Interface) and we use up to 2,048 processors for the highest resolution considered in this study. The entire cylindrical domain is discretized using 53,760 elements and a spectral order of N = 11. The element distribution is denser close to the external and bottom boundaries. This, in addition to the Gauss–Lobatto–Legendre point distribution close to the element boundaries, ensures an appropriate resolution of the Ekman boundary layers with approximately six grid points to describe them in all cases.

The vertical extent of the cylindrical domain is stretched so that the upper boundary matches to the analytical paraboloid shape that depends on the rotation

rate and on the gravitational acceleration. This surface is assumed to be stress-free and impenetrable. To mimic the experimental forcing, a small-scale steady volumic force is defined in Cartesian coordinates by

$$F = \begin{pmatrix} AS(r) \sin(\pi k_f x) \cos(\pi k_f y) \sin(\pi k_f z) \\ AS(r) \cos(\pi k_f x) \sin(\pi k_f y) \sin(\pi k_f z) \\ 0 \end{pmatrix} \quad (1)$$

if $0 \le z \le 1/k_f$ and $F=0$ elsewhere. The function $S(r)$ is given by

$$S(r) = \frac{1}{2} \left[ \tanh(60(r - 0.15)) - \tanh(60(r - 0.92)) \right] \quad (2)$$

A is the ratio of the force amplitude to the centrifugal acceleration and is related to the Rossby number (see Supplementary Fig. 1b). This forcing imposes a 3D small-scale turbulent flow at the base of the numerical domain, mimicking the turbulence in the experiments. The numerical simulations presented correspond to $A=0.002$ and $k_f=12$. Note that the typical jets observed in our simulations have a corresponding radial wavenumber $k=4-6$ so that we reach a reasonable scale separation between the energy injection mechanism and the large-scale zonal flow. We have checked that the results do not depend on $k_f$ although it is difficult to modify it while keeping a large Reynolds number and a low Rossby number. Thus, we considered only $8 < k_f < 16$.

**Non-dimensional analysis and zonostrophic regime.** Part of the experimental challenge was to drive the whole apparatus, containing 400 l of water, at a spin rate of $\Omega=75$ r.p.m. for the purpose of reducing as low as possible the Ekman and Rossby numbers of a developed turbulent flow. The Rossby number is the ratio of inertial to Coriolis forces defined as $Ro = U_{\psi \text{, m.s.}}^q / 2\Omega h_0$, where we use $U_{\psi \text{, m.s.}}^q$ the total r.m.s. velocity as the characteristic system scale velocity and $h_0$ the fluid layer depth at rest. To simultaneously maintain a large Reynolds number, defined as $Re = U_{\psi \text{, m.s.}}^q h_0 / v$ with $v$ being the kinematic viscosity, and a low Rossby number, one has to rapidly rotate the system for Coriolis to dominate over inertial forces. The ratio between Rossby and Reynolds numbers leads to the Ekman number $E = v / 2\Omega h_0^2$. Our rapidly rotating deep-layer experiment strongly reduces this Ekman number to values that are almost unachievable for 3D direct numerical simulations. The Ekman number appears to be crucial to investigate jets saturation mechanisms. Another non-dimensional parameter is the dynamical aspect ratio of the jets that we define as $\eta = h_0 / L_{Rh}$, with $L_{Rh}$ being the typical jets width. This parameter compares vertical ($h_0$) and horizontal ($L_{Rh}$) characteristic length scales, setting whether or not the jets develop in a 'shallow' or 'deep' domain. In the limit of $\eta \to 0$, the system tends to be 2D as the horizontal scale is much larger than the jet depth. Two-dimensional flows prevent direct cascade of energy and kinetic energy is consequently transferred to large-scale structures$^{24}$. Thus, in shallow-layer turbulence, the dynamic confines to a quasi-2D system at a finite value of $\eta \ll 1$, where the upscale cascade of energy is favoured. On the contrary, in deep-layer models, the vertical length scale can exceed the horizontal length scale, leading to 3D turbulent structures. Due to strong background rotation (that is, low Rossby number), however, large-scale motions are nearly invariant in depth and become essentially 2D (ref. 23). The quasi-2D components of the flow field tend to drive energy upscale via a so-called inverse cascade while 3D motions persist at smaller scales.

We report comparative values of the four non-dimensional quantities estimated from recent experiments$^{11–14}$, as well as for Jupiter and Earth's oceans in Supplementary Table 1. This table shows that our experimental approach gets closer to Jupiter's parameter range (low Ro, E and large Re) than any previous experiments. Moreover, we designed our experimental device to promote deep-seated jets and exclude the possibility of a shallow-like scenario. Our set-up is characterized by a large dynamical aspect ratio ($\eta=5.15$), which clearly disqualifies the possibility of 2D turbulence due to depth confinement. Quasi-geostrophic structures are then required to cascade energy upscale and, consequently, jet formation sensitively depends on the Rossby number. The geostrophic organization of the flow is well supported by our numerical simulation (Fig. 1). This distinction between shallow and deep-seated jets is less clear in all other previous experiments with moderate aspect ratio ($0.57 \le \eta \le 1$) and in which the role of geostrophy is not well defined. The experimental device designed by Smith et al.$^{12}$ has a typical aspect ratio of $\eta=4.2$ leading to deep-seated jets. However, the associated zonostrophic index is much lower and large-scale baroclinic vortices dominate the dynamics.

In the presence of background rotation, the free surface of a fluid takes the shape of a paraboloid. In our cylindrical domain, radial variation of the fluid depth is defined as

$$h(r) = h_{\text{min}} + \left( \frac{\Omega^2}{2g} \right) r^2 \quad (3)$$

where r is the cylindrical radius and g is the gravitational acceleration. The centre of the tank, of minimum depth $h_{\text{min}} = h_0 - \Omega^2 R^2 / 4g$, represents the planetary pole, while the edge of the domain corresponds to lower latitudes. Note that the paraboloidal surface shape is the first-order Taylor expansion of the spherical
curvature near the polar region. Under rapid rotation, the classical turbulent cascade, which usually dissipates kinetic energy from large to small scales, reverses direction to feed large-scale coherent structures. Similarly to planetary spherical configuration, depth variation h(r) induces a zonal anisotropization of the reversed energy cascade (the so-called zonostrophic turbulence²³) following the topographic β-effect quantified by

$$ \beta = \frac{2\Omega}{h} \left| \frac{dh}{dr} \right| \qquad (4) $$

where topographic paraboloidal β is zero at the pole and increases with radius.

The presence of this β-effect selects a transitional wavenumber at which the inverse cascade transitions to Rossby-wave-dominated dynamics:

$$ k_{\beta} \approx \left( \frac{\beta^3}{\epsilon} \right)^{1/5} \qquad (5) $$

where ϵ is the upscale energy transfer rate. Following Vallis and Maltrud²³, the anisotropy wavenumber k_β defines the scale at which the eddy turnover timescale τ_e = ϵ⁻¹/³k⁻²/³, with k being the eddy wavenumber, becomes comparable with the inverse of the Rossby wave frequency τ_R₀ = k²/βk_ψ, with k_ψ being the azimuthal wavenumber. k_β is an estimation of the smallest scale at which anisotropic Rossby wave propagation, due to boundary curvature together with rotation, affects the turbulent inverse cascade. The other main scale of β-turbulence comes out of the typical advection timescale τ_a = (Uk)⁻¹, where U ≈ U_T^T_{r.m.s.} is the characteristic system-scale flow velocity that advects turbulent eddies of scale k, and from which we recover the Rhines scale²⁵ k_Rh = √β/2U.

In the context of planetary inverse cascades, dimensional analysis can be reduced to the single zonostrophy index defined as the ratio between transitional wavenumber and Rhines wavenumber, R_β = k_β/k_Rh (ref. 27). This index can be seen as an estimate of the inverse turbulent cascade efficiency to grow zonal structures on a characteristic timescale τ_e compared with the turbulent structures turnover time τ_a. Large β-effect coupled with a strong Coriolis force leads to large index value of R_β > 2.5, corresponding to τ_e < τ_a. This threshold distinguishes between the zonostrophic and partially zonostrophic regimes.

It is not easy to access the value of ϵ in laboratory experiments. It can be estimated as ϵ ≈ (U_T^T_{r.m.s.})²/2ε_T based on the total horizontal kinetic energy and the Ekman spinup time¹⁴ τ_E = h₀/√v²Ω. Using U_T^T_{r.m.s.} = 0.035 m s⁻¹ and U_T^T_{r.m.s.} = 0.015 m s⁻¹ for respectively our experiment and the numerical simulation, we obtain ϵ = 3.4 × 10⁻⁶ m² s⁻¹ and ϵ = 6.3 × 10⁻⁷ m² s⁻¹. We end up with zonostrophic index values in the experiment 4.9 < R_β < 5.2, using experimental energy transfer rate and with values of 4.5 < R_β < 4.8 in the numerical simulation. This zonostrophic index is estimated in the range of β from 60 m⁻¹ s⁻¹ to 89 m⁻¹ s⁻¹ in the cylindrical radius range 0.2 < r/R < 0.9 in which jets form. One can rewrite the zonostrophy index in the form

$$ R_{\beta} = 2^{7/10} (\beta U_{r.m.s.}^{T} \tau_{E}^{2})^{1/10} \qquad (6) $$

An overall view of the zonostrophic regime is given by spectral analysis where the total energy spectrum E(k) can be decomposed into zonal and residual (or non-zonal) components E(k) = E_z(k) + E_R(k) with the total wavenumber in the horizontal plane k = √k_x² + k_ψ², where k_x and k_ψ are, respectively, the radial and azimuthal wavenumbers. The zonal spectrum can also be written as the axisymmetric component of the total energy spectrum E_z(k) = E_z(k_x, k_ψ = 0). Following Sukoriansky et al.³⁶, in the zonostrophic regime, the zonal and residual spectra are

$$ E_z(k) = C_z \beta^2 k^{-5} \qquad E_R(k) = C_R \epsilon^{2/3} k^{-5/3} \qquad (7) $$

where numerical simulations³⁷ give C_z ≃ 0.5 and C_R ≃ 5 − 6. More recently, Galperin et al.²⁹ re-evaluated the zonal constant to C_z = 2 from spectral analysis of high-resolution images of Jupiter's clouds. Using expressions (7) one can estimate the wavenumber at which zonal and residual spectra intersect, recovering the transitional wavenumber mentioned above,

$$ k_{\beta} = \left( \frac{C_z}{C_R} \right)^{3/10} \left( \frac{\beta^3}{\epsilon} \right)^{1/5} \qquad (8) $$

In their simulations, Galperin et al.³⁸ estimated that (C_z/C_R)³/10 ≃ 0.5, which might lead to a more accurate estimation of k_β than from dimensional analysis equation (5). By considering the recent estimate of C_z = 2 based on Jupiter's zonal

spectrum, one can re-calculate the prefactor (C_z/C_R)³/10 = 0.73 in equation (8). This notably reduces the value of the transitional wavenumber and consequently the corresponding zonostrophic index. Similarly, we use our range of β along the cylindrical radius and we report an estimate of the zonal constant that best fits our zonal spectrum. We obtain 1.7 < C_z < 3.7 with our experimental data and 0.42 < C_z < 0.9 with our numerical data using equation (7). Interestingly, the zonal constant ranges evaluated from our experiment and our numerical simulation include respectively the values computed from direct observations of Jupiter and previous numerical simulations reported²⁹,³⁸. It is important to have in mind that Galperin et al.²⁹ reduce the zonal constant to a single value while β also varies with Jupiter's latitude. Using the new expression (8) for k_β, one can rewrite the expression of the zonostrophy index in the form²⁹,³⁸:

$$ R_{\beta} = 2^{1/2} \left( \frac{C_z}{C_R} \right)^{3/10} \left( \frac{\beta (U_{r.m.s.}^T)^5}{\epsilon^2} \right)^{1/10} \qquad (9) $$

The zonostrophic index values in the main text are based on C_z = 2, which is in good agreement with our zonal kinetic energy spectrum and our mean β value.

We compute the zonal spectra E_z defined by equation (7) from numerical and experimental data using 1D Fourier transforms along the radius of the mean zonal velocity (U_ψ) (that is, temporal and azimuthal average of the zonal velocity). Being limited to Lagrangian tracers to access surface velocity in the experiment, we compute the residual spectrum from numerical simulation alone, using a horizontal square grid at a fixed depth h = 10 cm. The residual spectrum E_R defined by equation (7) is obtained from 2D Fourier transform of the non-zonal flow component defined as velocity deviations from the mean flow.

Figure 3b shows good agreement of the zonal and residual spectra with the theoretical predictions of equation (7). Residual and zonal energy spectra intersect at the transitional wavenumber predicted in equation (8) and kinetic energy cascades upscale to accumulate in the zonal component at Rhines scale²⁵. Injection scales presented in Fig. 3b range from an upper wavenumber defined as the half lattice grid size in the experiment and a lower wavenumber defined by numerical forcing. Using our estimation of the inverse energy cascade rate ϵ = 6.3 × 10⁻⁷ m² s⁻¹ and C_k = 5, we obtain a good prediction of our kinetic energy residual spectrum, which confirms our zonostrophic turbulence approach. Ultimately, with (C_z/C_R)³/10 ≃ 0.73 and 60 m⁻¹ s⁻¹ < β < 89 m⁻¹ s⁻¹, we compute the zonostrophy index from equation (9). We obtain 3.6 < R_β < 3.8 in our laboratory experiment and 3.3 < R_β < 3.5 in the numerical simulation, which confirms that we are well into the zonostrophic regime. Mean values of R_β = 3.7 and R_β = 3.4, using the mean value β = 73 m⁻¹ s⁻¹, are reported in the main text. Values of R_β for planetary objects and other experimental attempts are detailed in Supplementary Table 1. Expression (9) has been used to compare all experiments and Jupiter's observations, reducing the zonostrophic index of previous experiments by a factor of 0.73.

**Data availability.** The data that support the plots within this paper and other findings of this study are available from the corresponding author on request.

## References

31. Sommeria, J., Meyers, S. D. & Swinney, H. L. Laboratory model of a planetary eastward jet. Nature 337, 58–61 (1989).

32. Aubert, J., Jung, S. & Swinney, H. L. Observations of zonal flow created by potential vorticity mixing in a rotating fluid. Geophys. Res. Lett. 29, 1876 (2002).

33. Kelley, D. H. & Ouellette, N. T. Using particle tracking to measure flow instabilities in an undergraduate laboratory experiment. Am. J. Phys. 79, 267–273 (2011).

34. Deville, M. O., Fischer, P. F. & Mund, E. H. High-Order Methods for Incompressible Fluid Flow (Cambridge Univ. Press, 2002).

35. Favier, B., Grannan, A. M., Le Bars, M. & Aurnou, J. M. Generation and maintenance of bulk turbulence by libration-driven elliptical instability. Phys. Fluids 27, 066601 (2015).

36. Sukoriansky, S., Galperin, B. & Dikovskaya, N. Universal spectrum of two-dimensional turbulence on a rotating sphere and some basic features of atmospheric circulation on giant planets. Phys. Rev. Lett. 89, 124501 (2002).

37. Huang, H.-P., Galperin, B. & Sukoriansky, S. Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere. Phys. Fluids 13, 225–240 (2001).

38. Galperin, B., Sukoriansky, S. & Dikovskaya, N. Zonostrophic turbulence. Phys. Scripta 2008, 014034 (2008).