domain: sdu phys

text:
D Lemasquerier
email: lemasquerier.pro@protonmail.com

Benjamin Favier
email: benjamin.favier@univ-amu.fr

M Le Bars
email: lebars@irphe.univ-mrs.fr

M Le

Zonal jets experiments in the gas giants' zonostrophic regime

Keywords:

Intense east-west winds called zonal jets are observed in the atmospheres of Jupiter and Saturn and extend in their deep interior. We present experimental results from a fully three-dimensional laboratory analog of deep gas giants zonal jets. We use a rapidly rotating deep cylindrical tank, filled with water, and forced by a small-scale hydraulic circulation at the bottom. A topographic 𝛽-effect is naturally present because of the curvature of the free surface. Instantaneous turbulent zonal jets spontaneously emerge from the small-scale forcing, equilibrate at large scale, and can contain up to 70% of the total kinetic energy of the flow once in a quasi-steady state. We show that the spectral properties of the experimental flows are consistent with the theoretical predictions in the zonostrophic turbulence regime, argued to be relevant to gas giants. This constitutes the first fully-experimental validation of the zonostrophic theory in a completely three-dimensional framework. Complementary, quasi-geostrophic (QG) simulations show that this result is not sensitive to the forcing scale. Next, we quantify the potential vorticity (PV) mixing. While PV staircasing should emerge in the asymptotic regime of the gas giants, only a moderate PV mixing occurs because of the strong forcing and dissipation, as confirmed by QG simulations at smaller Ekman number. We quantify the local PV mixing by measuring the equivalent of a Thorpe scale, and confirm that it can be used to estimate the upscale energy transfer rate of the flow, which otherwise needs to be estimated from a much more demanding spectral analysis.

Introduction

Zonal jets are east-west currents observed in a large variety of geophysical flows, from the atmospheres of gas giants [START_REF] Vasavada | Jovian atmospheric dynamics: An update after Galileo and Cassini[END_REF] to terrestrial oceans [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF] and atmospheres [START_REF] Mitchell | Terrestrial Atmospheres[END_REF]. They are also expected to develop in buried fluid layers such as the outer liquid cores of telluric planets [START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF] or the subsurface oceans of icy moons [START_REF] Soderlund | Ice-Ocean Exchange Processes in the Jovian and Saturnian Satellites[END_REF]. The vast occurrence of zonal jets underlines their generic nature, which relies on three basic physical effects: the dominance of rotation on the flow, the presence of turbulent motions forced by various processes, and a 𝛽-effect which represents the variation of the Coriolis force with latitude. Apart from their clear geophysical relevance, zonal jets are outstanding fluid dynamical features where a turbulent flow self-organizes at large-scale. The properties of the associated statistically steady-state in terms of number, width, intensity and stability of the jets are described by several theories based on the phenomenology of turbulence [START_REF] Vallis | Generation of Mean Flows and Jets on a Beta Plane and over Topography[END_REF][START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF], potential vorticity mixing [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF], wave-mean flow interaction [START_REF] Quinn | Modulational Instability in Basic Plasma and Geophysical Models[END_REF] or statistical physics [START_REF] Bouchet | Statistical mechanics of two-dimensional and geophysical flows[END_REF][START_REF] Farrell | Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection[END_REF], but a common and robust physical framework is still lacking.

In the present study, we specifically focus on the regime of turbulent zonal jets observed on the gas giants, Jupiter and Saturn, where the jets' dynamics is not complicated by topographic features or seasonal variations. Jupiter's zonal jets are famous because they are responsible for the banded appearance of the planet, through the spreading of clouds of ammonia and water ices. Jupiter's intense zonal winds have been measured by various spacecrafts using cloud tracking, from Voyager in 1979 to Juno, which is in orbit around Jupiter since 2016 (Bolton, Lunine, Stevenson, Connerney, Zonal jets experiments in the gas giants' zonostrophic regime [START_REF] Bolton | The Juno Mission[END_REF]. Direct measurements revealed the remarkable stability of these strong winds over decades [START_REF] Tollefson | Changes in Jupiter's Zonal Wind Profile preceding and during the Juno mission[END_REF], and allowed to demonstrate that they contain up to 90% of the kinetic energy of the turbulent flow, at least at the level of the cloud layer (Galperin, Young, Sukoriansky, Dikovskaya, Read, Lancaster and Armstrong, 2014b). Recently, the inversion of the gravity measurements of the Juno spacecraft showed that the Jovian jets extend down to about 3,000 kilometers beneath the clouds for Jupiter [START_REF] Kaspi | Jupiter's atmospheric jet streams extend thousands of kilometres deep[END_REF][START_REF] Guillot | A suppression of differential rotation in Jupiter's deep interior[END_REF], and 9,000 kilometers for Saturn [START_REF] Galanti | Saturn's Deep Atmospheric Flows Revealed by the Cassini Grand Finale Gravity Measurements[END_REF] based on Cassini measurements. This result is consistent with independent constraints provided by the secular variation of Jupiter's magnetic field [START_REF] Moore | Time variation of Jupiter's internal magnetic field consistent with zonal wind advection[END_REF]. Despite the wealth of new constraints provided by Juno, comprehensive idealized models still need to be developed to complement observations and better understand the jets origin, properties, and three-dimensional structure.

From the fluid dynamics point of view, gas giants exhibit very intense flows, characterized by a very large Reynolds number (𝑅𝑒), a vanishing Ekman number (𝐸) and a small Rossby number (𝑅𝑜). The Reynolds number quantifies the ratio of inertial to viscous forces, the Ekman number the ratio of viscous to Coriolis forces, and the Rossby number the ratio of inertial to Coriolis forces. Denoting 𝑈 a typical velocity scale, 𝐻 a typical flow length scale, 𝜈 the molecular viscosity of the fluid and Ω the planet's rotation rate, these non-dimensional parameters can be expressed as

𝑅𝑒 ∼ 𝑈 𝐻 𝜈 , 𝐸 ∼ 𝜈 Ω𝐻 2 , 𝑅𝑜 ∼ 𝑈 Ω𝐻 . ( 1
)
The aforementioned conditions hence mean that Jupiter's flows are both highly turbulent and rotationally constrained (see Table 1 and Fig. 1). Because rotation is dominant, the turbulent motions are quasi two-dimensional (2D) and may bear an inverse cascade of kinetic energy [START_REF] Young | Forward and inverse kinetic energy cascades in Jupiter's turbulent weather layer[END_REF][START_REF] Siegelman | Moist convection drives an upscale energy transfer at Jovian high latitudes[END_REF] feeding large-scale features (vortices and jets). We denote 𝜖 the corresponding rate of energy transfer. The zonation of Jupiter's turbulence then arises because of the so-called 𝛽-effect, representing either the latitudinal variation of the projection of the planet's rotation vector on the local vertical (shallow, atmospheric 𝛽-effect), or the variation of the fluid height when measured parallel to the rotation axis (deep, topographic 𝛽-effect). Because of the 𝛽-effect and associated Rossby waves, the turbulent energy transfer to large scales becomes anisotropic and redirected towards zonal currents. Ultimately, the large scale drag halts the development of the inverse cascade [START_REF] Sukoriansky | Universal Spectrum of Two-Dimensional Turbulence on a Rotating Sphere and Some Basic Features of Atmospheric Circulation on Giant Planets[END_REF][START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF]. A fourth nondimensional parameter, known as the zonostrophy index 𝑅 𝛽 [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF], is used as an indicator of the strength of the zonation of the turbulent flow. It is defined as the ratio between the Rhines scale [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF]) 𝐿 𝑅 ,

𝐿 𝑅 ∝ ( 2𝑈 𝛽 ) 1∕2 , ( 2
)
assumed to represent the scale of the jets, and the transitional scale [START_REF] Pelinovsky | Wave turbulence on beta-plane[END_REF][START_REF] Vallis | Generation of Mean Flows and Jets on a Beta Plane and over Topography[END_REF]) 𝐿 𝛽 ,

𝐿 𝛽 ∝ ( 𝜖 𝛽 3 ) 1∕5 ,
(3) which is the scale at which a turbulent eddy turnover time equates a Rossby wave period. The transitional scale is the smallest scale above which the turbulent flow becomes anisotropic due to the 𝛽-effect. The zonostrophy index is then

𝑅 𝛽 ∼ 𝐿 𝑅 𝐿 𝛽 ∼ 2 1∕2 ( 𝑈 5 𝛽 𝜖 2
) 1∕10 (4) [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF]. The regime of strong and rectilinear jets -so-called zonostrophic regime -is obtained when the scale at which the eddies start being deformed by the 𝛽-effect is well separated from the scale of the final jets, i.e. for a large zonostrophy index. From 2D simulations on the sphere, [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF]; [START_REF] Galperin | Geophysical flows with anisotropic turbulence and dispersive waves: Flows with a 𝛽-effect[END_REF] show that the zonostrophic regime is observed when 𝑅 * 𝛽 ≳ 2.5 whereas the regime is frictiondominated when 𝑅 * 𝛽 ≲ 1.5, as represented in Fig. 1(b). Here, 𝑅 * 𝛽 = (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 𝑅 𝛽 ≈ 0.5𝑅 𝛽 , and 𝐶 𝑍 and 𝐶 𝐾 are Non-dimensional parameters and typical values. The Ekman number, 𝐸, Reynolds number 𝑅𝑒 and zonostrophy index 𝑅 𝛽 are independent, the Rossby number 𝑅𝑜 is indicated for completeness but can be expressed as 𝑅𝑜 = 𝑅𝑒 × 𝐸. For the experiments and QG simulations, assuming dissipation in the Ekman boundary layers (equation (7)), these parameters can be evaluated either using the fluid height ℎ 0 (first column) or the upscale energy transfer rate 𝜖 (second column). 𝑢 rms is the total and time-averaged rms velocity measured once in statistically steady-state. This is relevant only for the experiments and QG simulations, not the planetary flows. The planetary parameters used to compute the non-dimensional parameters for the gas giants and Earth's atmosphere and oceans are provided in Appendix B. 𝑎 see Table 13.1 in Galperin and Read (2019) supposedly universal constants determined from spectral analysis (see section 4). Jovian mid-latitude jets have 𝑅 * 𝛽 ≈ 5 and are expected to be in the regime of zonostrophic turbulence (Galperin et al., 2014b). On the contrary, Earth's oceans are under the threshold with 𝑅 * 𝛽 ≈ 1.5 [START_REF] Galperin | Barotropic and Zonostrophic Turbulence[END_REF]. Consistently, Jovian jets are strong, instantaneous and contain most of the kinetic energy of the flow (Galperin et al., 2014b) whereas oceanic jets are weak and meandering, hence a careful time-averaging is required to reveal them [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF]. 𝛽 ) (Reynolds number, Ekman number, and zonostrophy index). For the sake of comparison, we used similar definitions for all the experimental studies (first column of Table 1). The data used to plot this figure are provided in Appendix B. We recall that 𝑅 * 𝛽 ≈ 0.5𝑅 𝛽 and takes into account the spectral prefactor (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 (see section 4). The horizontal dashed lines represent the threshold in 𝑅 * 𝛽 between the different regimes determined by [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF]; [START_REF] Galperin | Geophysical flows with anisotropic turbulence and dispersive waves: Flows with a 𝛽-effect[END_REF] .

Modeling turbulent zonal jets in regimes relevant to gas giants requires to fulfill these four dynamical constraints simultaneously (𝑅𝑒 ≫ 1, 𝐸 ≪ 1, 𝑅𝑜 ≲ 1, 𝑅 * 𝛽 ≳ 2.5). This is particularly challenging, both from the numerical and experimental point of view because of computational and technical constraints. This regime can only be approached by simplified models relying on physical approximations whose relevance needs to be systematically addressed. For instance, current 3D numerical models cannot reach the asymptotic Reynolds and Ekman numbers of gas giants flows, and when trying to do so, it becomes extremely costly to simulate the dynamics over the very long radiative and frictional timescales, necessary for the system to equilibrate (see e.g. [START_REF] Heimpel | Simulation of deep-seated zonal jets and shallow vortices in gas giant atmospheres[END_REF][START_REF] Gastine | Stable stratification promotes multiple zonal jets in a turbulent Jovian dynamo model[END_REF]. On the contrary, experiments usually allow to reach more extreme regimes than numerical models, and study the dynamics of real, fully nonlinear and developed flows in a statistically steady-state, even if they cannot reproduce the very small viscous effects of Jovian flows. Of course, experiments come with their own limitations, such as necessary physical boundaries, partial measurements, and the difficulty in incorporating magneto-hydrodynamical or compressibility effects.

In a previous paper [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], we described a new experimental setup following on from [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] where zonal jets spontaneously emerge from a mechanically-forced turbulent flow in a rapidly rotating water tank. This setup is a one layer system, with a deep water height, and hence the zonal jets are deep and barotropic. Here, we would like to stress out that a barotropic (one layer) experimental setup is still relevant to model a priori strongly baroclinic systems, such as oceans and atmospheres. Indeed, planetary, quasi-geostrophic flows can undergo a phenomenon of barotropization, and the barotropic modes can end up containing a significant fraction of the energy, even if the energy feeding these modes initially comes from baroclinic modes and instabilities [START_REF] Charney | Geostrophic Turbulence[END_REF][START_REF] Rhines | The dynamics of unsteady currents[END_REF][START_REF] Salmon | Two-layer quasi-geostrophic turbulence in a simple special case[END_REF] (see section 6.3 for further discussion). In the experiment of [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], the zonal jets emerge due to the radiation of Rossby waves and the associated transport and deposition of momentum. In addition, we observed a transition between two regimes of multiple zonal jets. We showed that the transition arises from a Rossby waves resonance due to their advection by the background zonal flow. We showed that the first regime (hereafter Regime I), where the jets are individual and locally forced, may be relevant for jets in the terrestrial oceans. The second regime (hereafter Regime II) is accompanied by a coarsening and an intensification of the zonal jets, and is relevant for the gas giants. In the present study, we focus on the saturated turbulent and statistically steady state obtained in Regime II, far from the transition, and we compare our experimental results with two theories which aim to explain zonal jets properties at a global scale.

The first theory on which we focus is the so-called theory of zonostrophic turbulence (see Galperin et al., 2019). In this framework, zonal jets are described as emerging from the anisotropisation of two-dimensional turbulence in the presence of a 𝛽-effect. However, this description may only be relevant when the four aforementioned dynamical constraints are fulfilled, and when sufficient scale separations are achieved such that turbulent inertial ranges can exist. Following [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF], our experimental setup is specifically designed to favor a large 𝑅 𝛽 , fast flows (𝑅𝑒 ≫ 1), but still dominated by rotation (𝑅𝑜 ≪ 1), and small viscous dissipation (𝐸 ≪ 1), thus getting closer to the regime observed on gas giants. As shown in Fig. 1(b), previous experimental studies lied in the range 𝑅 * 𝛽 ∈ [1,[START_REF] Chekhlov | The effect of small-scale forcing on large-scale structures in twodimensional flows[END_REF] and the observed flows were not in the zonostrophic regime. Let us stress out that, because of the definition of the zonostrophy index and the power 1/10 (see Table 1), if all the other parameters are unchanged, a rms velocity ten times larger (and hence a Reynolds ten times larger) only increases the zonostrophy index by a factor 10 1∕10 ∼ 1.26! Here, we explore more extreme regimes, including 𝑅 * 𝛽 ≳ 2.5. The experimental measurements can hence be compared with predictions in the so-called zonostrophic regime relevant to gas giants. The spectral analysis performed in [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] on a previous version of the experiment supports the idea that predictions from zonostrophic turbulence are retrieved experimentally. Kinetic energy spectra were however computed for the zonal flow only and not the fluctuations, due to the use of particle tracking velocimetry. With the new setup [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], the temporal and spatial resolution of our PIV measurements allows us to quantify turbulence statistics for both the zonal and fluctuating components of the flow, and to complete the analysis in the framework of zonostrophic turbulence.

The second theoretical framework to which we compare our experimental results is potential vorticity (PV) mixing [START_REF] Dunkerton | A barotropic model of the angular momentum-conserving potential vorticity staircase in spherical geometry[END_REF][START_REF] Mcintyre | Potential-vorticity inversion and the wave-turbulence jigsaw: Some recent clarifications[END_REF][START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF][START_REF] Scott | The generation of zonal jets by large-scale mixing[END_REF]. In our case of a barotropic fluid with varying height, the PV is 𝑞 = (𝜁 +𝑓 )∕ℎ, where 𝜁 is the relative vorticity, 𝑓 = 2Ω is the "planetary" vorticity (Ω is the rotation rate), and ℎ is the fluid height. In the framework of PV mixing, zonal jets are described as emerging from the mixing of the scalar field 𝑞, which is materially conserved in the absence of dissipation. Retrograde jets correspond to well-mixed regions of PV, whereas prograde jets correspond to steep gradients of PV, thus defining a staircase-like profile in the radial direction if multiple jets are present. The mechanism of staircase formation was described in [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF] by analogy with the mixing of a stratified fluid, the so-called "Phillipseffect" [START_REF] Phillips | Turbulence in a strongly stratified fluid-is it unstable[END_REF]. The idea is that of a feedback mechanism: stratification is reduced in the mixed regions, the restoring force at the origin of gravity waves (buoyancy) is thus weakened, which allows for further mixing. On the contrary, where interfaces are forming, steeper stratification locally develop, and gravity waves are enhanced, thus inhibiting mixing across the interface. For a PV gradient instead of a density gradient, the mechanism is analogous but invokes Rossby waves whose propagation depends on how steep the local PV gradient is. That being said, subtleties arise from the fact that PV is different from density since there is a direct relationship between the flow dynamics and PV through the vorticity.

Galperin, Hoemann, Espa and Nitto (2014a) proposed to use the analogy between PV and density staircases to quantify PV mixing, through the method of PV monotonizing. The idea is to use local and instantaneous PV (or density) profiles to quantify the turbulent overturn by measuring the equivalent of a Thorpe scale [START_REF] Thorpe | The Turbulent Ocean[END_REF][START_REF] Gargett | Determining Thorpe Scales from Ship-Lowered CTD Density Profiles[END_REF]. The initial PV (density) profile, due to background rotation and 𝛽-effect, is monotonous and decreasing with radius. Locally, turbulent eddies can cause overturns and bring fluid parcels with higher PV (higher density) outward (above) compared to a fluid parcel with smaller PV (smaller density); an unstable configuration. By sorting the non-monotonous PV profile into a monotonous one, one can define a Thorpe scale 𝐿 𝑇 , which is the rootmean-squared displacement of the fluid parcels needed to bring them back to a stable position. Galperin et al. (2014a) show experimentally that similarly to stratified flows where 𝐿 𝑇 is close to the Ozmidov scale [START_REF] Ozmidov | On the turbulent exchange in a stably stratified ocean[END_REF][START_REF] Thorpe | Turbulence and mixing in a Scottish Loch[END_REF], which is the scale at which turbulent eddies feel the stratification, 𝐿 𝑇 is commensurate with the transitional scale, 𝐿 𝛽 , at which turbulent eddies feel the 𝛽-effect. If a robust relationship exists between the Thorpe and transitional scales, then sorting PV would constitute a powerful tool to estimate 𝐿 𝛽 and indirectly, the turbulence intensity. The only requirement would then be to measure instantaneous profiles of potential vorticity of the flow (see for instance [START_REF] Cabanes | Revealing the Intensity of Turbulent Energy Transfer in Planetary Atmospheres[END_REF], for an application to Jupiter and Saturn). However, the three experiments used in (Galperin et al., 2014a) are far from the zonostrophic regime: the forcing, which is performed using magnets and a saline solution, accelerates a westward zonal flow locally but also directly (all the magnets are aligned and have the same polarity). To what extent the scaling between the Thorpe and the transitional scale holds in our experiments where the jets are strongly turbulent and self-developed from a fluctuating forcing with no azimuthal mean, thus remains to be addressed.

The present paper is organized as follows. In §2, we briefly recall the main characteristics of the experimental setup, and we present the numerical quasi-geostrophic model which is used to complement the experimental results. In §3, we describe the fully developed flows obtained in both the experiments and simulations. In §4, we compute kinetic energy spectra to compare with zonostrophic turbulence predictions and we measure a turbulent energy dissipation rate from the spectra. QG simulations are used to extend the experimental results to smaller forcing scales and hence better scale separation. In §5, we address the question of potential vorticity mixing. At a global scale, we quantify PV staircaising and extend our experimental results with simulations performed at a smaller Ekman number. At a local scale, we compute the Thorpe scale associated with PV mixing, which provides a second independent way of estimating the upscale energy transfer rate. We conclude and discuss our results in §6.

Experimental and numerical methods

Experimental set-up

Description of the setup

The experimental set-up is sketched in Fig. 2(a) and is extensively described in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF] and [START_REF] Lemasquerier | Experimental and Numerical Study of Jupiter's Dynamics: Jets, Vortices and Zonostrophic Turbulence[END_REF]. For consistency, we summarize here its principal features. It consists in a rotating cylindrical water tank, of radius 𝑅 = 0.49 m, filled with about 600L of water. The fast rotation of the tank induces a centrifugal force leading to a strongly deformed, paraboloidal free-surface. Since the fluid height varies with radius, the relative vorticity of a fluid column moving radially outward (stretching) or inward (squeezing) is respectively increased or decreased. This effect is the so-called topographic 𝛽-effect and is responsible for introducing anisotropy in the system in the form of zonation along the azimuthal direction. Here, the topographic 𝛽 parameter can be written as

𝛽 = - 𝑓 ℎ dℎ d𝜌 , ( 5
)
where 𝜌 is the cylindrical radius, ℎ(𝜌) is the total fluid height and 𝑓 = 2Ω is the Coriolis parameter with Ω the rotation rate. The vast majority of theoretical developments are performed in the 𝛽-plane approximation, where 𝛽 is spatially homogeneous -for planetary applications, this amounts to a local approach where the global variations of 𝛽 are neglected. Equation (5) shows that for the topographic 𝛽-effect to be uniform over the domain, the fluid height should vary exponentially with radius. To achieve this, we compensate the unalterable paraboloidal shape of the free surface using a curved bottom plate placed inside of the tank (Fig. 2(c)). The total fluid height ℎ above the bottom plate is the difference between the free-surface altitude ℎ 𝑝 and the bottom topography altitude ℎ 𝑏 . The bottom topography has been chosen such that, in solid body rotation at a fixed rate Ω, the total water height as a function of the cylindrical radius 𝜌 is

ℎ(𝜌) = ℎ 0 + Ω 2 2𝑔 ( 𝜌 2 - 𝑅 2 2 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ℎ 𝑝 (𝜌) -ℎ 𝑏 (𝜌) = ℎ min exp ( - 𝛽 2Ω 𝜌 ) , ( 6
)
where 𝑔 is the gravitational acceleration, 𝑅 is the tank radius, ℎ 0 the fluid height at rest, ℎ min the minimum fluid height once in rotation and 𝛽 < 0 the spatially uniform 𝛽-parameter (eq.( 5)). We worked with Ω = 75 RPM (rotation frequency of 1.25 Hz), ℎ 0 = 0.58 m, ℎ min = 0.20 m and 𝛽 = -50.1 m -1 s -1 . Note that our experimental set-up was designed to work at a single rotation rate of 75 RPM at which the 𝛽-effect is then uniform across the tank. However, we performed few experiments at rotation rates of 60 and 80 RPM to investigate the consequences of a modified 𝛽-effect. At these rotation rates the 𝛽-effect is no more uniform but nevertheless varies in a limited range (|𝛽 80 | ∈ [57, 73] m -1 s -1 with a mean at 65.5 m -1 s -1 and |𝛽 60 | ∈ [20, 30] m -1 s -1 with a mean at 22.8 m -1 s -1 ).

Once solid-body rotation is reached, the flow is forced by circulating the tank's water through 128 holes drilled into the curved bottom plate (see Fig. 2(a) and (c)). The holes are distributed on a polar lattice composed of six rings each connected to a separate pump. Half of the holes are inlets, and the other half are outlets, with as many inlets as outlets along each ring. Because of the dominant rotation, this forcing generates an array of cyclonic and anticyclonic vortices, which subsequently break into turbulence. The zonal average of the forcing is zero by construction, both along each ring and in total, meaning that we do no directly accelerate the zonal flow. Additional details on the experiment can be found in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF] and D. Lemasquerier's thesis manuscript [START_REF] Lemasquerier | Experimental and Numerical Study of Jupiter's Dynamics: Jets, Vortices and Zonostrophic Turbulence[END_REF]. Two different sets of pumps were used to reach different forcing amplitudes. They are described in Appendix B of [START_REF] Lemasquerier | Experimental and Numerical Study of Jupiter's Dynamics: Jets, Vortices and Zonostrophic Turbulence[END_REF].

To measure velocity fields, time-resolving particle image velocimetry (PIV) measurements are performed on a horizontal plane. A green horizontal laser sheet (532 nm) crosses the water layer 11 cm above the bottom plate (9 cm below the center of the paraboloid). The water is seeded with fluorescent red polyethylene particles and their motion is tracked using a top-view camera placed above the tank and embarked in solid-body rotation with the tank and the laser (Fig. 2(a)). The particles emit an orange light (607 nm) so that using a high-pass filter on the lens allows to filter out the green laser reflections on the free-surface and tank sides, leading to a better image quality and hence better PIV measurements. The images are acquired using DANTEC's software DynamicStudio. Optical distortion induced by the paraboloidal free-surface is corrected on DynamicStudio using a preliminary calibration performed by imaging a home-made calibration target with a precise dot pattern [START_REF] Lemasquerier | Experimental and Numerical Study of Jupiter's Dynamics: Jets, Vortices and Zonostrophic Turbulence[END_REF]. A movie of the particles motion during a typical experiment is available as supplementary movie 1 in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. The velocity fields are deduced from these images using the MATLAB program DPIVSoft developed by [START_REF] Meunier | Analysis and treatment of errors due to high velocity gradients in particle image velocimetry[END_REF]. We consider 32×32 pixels boxes on 1900×1900 pixels images and obtain 100×100 velocity vector fields (40% overlap between the boxes). Note that due to the refraction of the laser plane by the tank sides, there are two shadow zones where measurements are not possible (see the grey areas in Fig. 5(a)).

Note that the present setup was built from scratch but inspired from the setup used in [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF]. In addition to the time-resolved PIV measurements mentioned above, the main improvements of our setup compared to the one of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] include:

-A spatially-uniform 𝛽-effect to be able to compare with theories derived in the 𝛽-plane approximation.

-A control of the spatial distribution of the forcing. In [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF], all the injection and suction points were linked to the same pump. Here, the forcing pattern is polar such that we can control independently the six forcing rings using six different pumps and modulate the intensity of the forcing with radius.

-An increase in the number of forcing points from 64 to 128 in order to decrease the forcing scale and have a better scale separation with the transitional scale 𝐿 𝛽 .

Non-dimensional parameters

The experimental flows we consider are defined by five dimensional parameters, the rotation rate, Ω, the mean fluid height, ℎ 0 , the molecular viscosity, 𝜈, the initial gradient of potential vorticity 𝛽 and the root-mean-squared (rms) velocity, 𝑢 rms . Here, 𝑢 rms is an output parameter measured as the total rms velocity once in steady state. It is linearly related to the forcing amplitude, which is an input parameter set by choosing the power of the pumps (see Lemasquerier , 2021). With two dimensions (length and time), our experiments can be characterized by three independent nondimensional parameters, which we chose to be the Ekman number, 𝐸, the Reynolds number, 𝑅𝑒 and the zonostrophy index 𝑅 𝛽 , defined in the introduction, where we chose 𝑢 rms and ℎ 0 as the typical fluid velocity and length scales 𝑈 and 𝐻, respectively.

In our laboratory experiments (but not in planetary flows), a consequence of dominant rotation is that one can assume that dissipation mainly occurs in the Ekman boundary layer forming along the frictional bottom plate. In such case, the upscale energy transfer 𝜖 can be expressed as a function of the Ekman spin-down timescale 𝜏 𝐸 = Ω -1 𝐸 -1∕2 :

𝜖 ∼ 𝑢 2 rms 2𝜏 𝐸 = 𝑢 2 rms (𝜈Ω) 1∕2 2ℎ 0 . ( 7
)
We will see later that 𝜖 can actually be measured, and agree relatively well with this estimate (section 4). In the remaining of the paper, and as indicated in Table 2, we make the difference between the estimated energy transfer rate, 𝜖 𝐸 , the transfer rate measured from the spectra, 𝜖 𝑆 , and the transfer rate measured from local potential vorticity mixing, 𝜖 𝑇 . Let us define a frictional lengthscale 𝐿 𝐸 as 𝐿 𝐸 ∼ 2𝑢 rms 𝜏 𝐸 = 𝑢 3 rms 𝜖 -1 . The zonostrophy index can then be alternatively expressed as:

𝑅 𝛽 = 𝐿 𝑅 𝐿 𝛽 = 2 1∕2 ( 𝐿 𝐸 𝐿 𝑅 ) 1∕5 = 2 1∕2 ( 𝐿 𝐸 𝐿 𝛽 ) 1∕6 . ( 8
)
These expressions show that 𝑅 𝛽 simultaneously compares the transitional scale, the Rhines scale, and the friction scale:

• if 𝑅 𝛽 < 1 (𝐿 𝐸 < 𝐿 𝑅 < 𝐿 𝛽 ), then the 𝛽-effect is weak and the scale of the flow is defined by the frictional scale 𝐿 𝐸 which is then smaller than 𝐿 𝑅 and 𝐿 𝛽 . In other words, since the frictional scale is smaller than the transitional scale, the energy is dissipated before the flow is fully affected by the 𝛽-effect. In that case, one would expect the flow to remain nearly isotropic and show features of the classical Kolmogorov-Batchelor-Kraichnan (KBK) turbulence (see [START_REF] Boffetta | Two-Dimensional Turbulence[END_REF], and references therein).

• if 𝑅 𝛽 > 1 (𝐿 𝛽 < 𝐿 𝑅 < 𝐿 𝐸 ), then the 𝛽-effect is strong, and the turbulence is anisotropic since its equilibrium scale is larger than the transitional scale 𝐿 𝛽 . The large scales are then limited at a value less than 𝐿 𝐸 .

Table 1 provides the expression of the non-dimensional parameters of the problem (𝑅𝑒,𝐸,𝑅 𝛽 ) as a function of (𝑢 rms , ℎ 0 , Ω, 𝜈, 𝛽) or alternatively (𝑢 rms , 𝜖, Ω, 𝜈, 𝛽), where 𝜖 is approximated by equation ( 7). As previously mentioned, to model gas giants flows, both 𝑅𝑜 should be small and 𝑅𝑒 should be large, which can only be achieved with an asymptotically small 𝐸, hence the need for a large tank, rotating fast. In our experiment, these constraints are fulfilled, as indicated by the values reported in Table 1. Note that our Ekman number is much higher than that of Jupiter, and our Reynolds is much smaller, but these discrepancies compensate in the Rossby number which has the good order of magnitude. Physically, this means that we have the good ratio of inertial to Coriolis forces, but our experimental flow is overly damped by bulk viscous effects which are vanishingly small on Jupiter. We are nevertheless in the relevant regime, turbulent while still dominated by rotation (𝑅𝑒 ≫ 1, 𝑅𝑜 ≪ 1).

List of experiments

In Table 3, we report the physical and non-dimensional parameters of the experiments that are used in the following of the paper. Because of experimental constraints, it is not possible to vary independently the three non-dimensional parameters 𝐸, 𝑅𝑒 and 𝑅 𝛽 . First, the essential control parameter that we have is the forcing intensity (i.e. indirectly 𝑢 rms ), which allows us to explore significantly different Reynolds numbers. However, increasing 𝑢 rms also increases 𝑅 𝛽 by decreasing the relative importance of large-scale drag. Second, we cannot explore a large range of Ekman numbers because the rotation rate of the experiment is fixed, as mentioned above. Third, for the three experiments performed at a different rotation rate, both 𝐸 and 𝑅 𝛽 are modified since changing the rotation rate also modifies the 𝛽-effect by changing the fluid height.

Quasi-geostrophic numerical model

To complement and strengthen the experimental results, we perform numerical simulations. The goal is in particular to explore the sensitivity of the results to the forcing scale and pattern, which is difficult to do experimentally. Indeed, it would require to design and build a new bottom plate for each forcing pattern, and the numbers of inlets and outlets (128) would be difficult to increase significantly. To numerically model the experiment, 3D direct numerical simulations (DNS) are not adapted neither as they are very computationally demanding, in particular because the simulation has to resolve a turbulent flow with both very large scale structures (the jets) and thin Ekman boundary layers that are essential for the long term dynamics (see [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF]. Instead, we model the experiment using a quasi-geostrophic (QG), two-dimensional numerical model. QG models take advantage of the fast background rotation, or equivalently the small Rossby number of the system: since the geostrophic balance dominates the experimental flow, the flow is quasi two-dimensional. The curvature of the free-surface as well as the friction over the bottom (Ekman pumping) induce three-dimensional effects, but their weakness allows their incorporation into quasi-two-dimensional physical models. With a QG model, we can exactly match the experimental conditions in terms of Reynolds, Rossby and Ekman numbers at a moderate computational cost, such that we can run several simulations over time scales comparable with that of the experiments (thousands of rotation times, or equivalently tens of frictional time scales). This would have been inconceivable using 3D direct numerical simulations.

In Appendix C, we derive the QG equations used in our numerical code. The difference with the conventional QG model is that, following Sansón andVan Heijst (2000, 2002) we retain higher order, non-linear terms for the Ekman pumping and 𝛽-effect. We consider a polar domain delimited by the outer and inner radii 𝑟 𝑜 = 𝑅 and 𝑟 𝑖 = 0.05𝑅, where 𝑅 is the total radius. The inner boundary, absent in the experiment, is introduced to avoid the coordinate singularity at the center of the domain. We use the polar coordinates (𝜌,𝜙) and denote (𝐞 𝜌 , 𝐞 𝜙 ) the associated unit vectors (Fig. 2). We solve separately the axisymmetric (zonal flow) and non-axisymmetric motions (see e.g. [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF]. To this end, we perform a Reynolds decomposition of the velocity field into an azimuthal average plus some fluctuations, denoted with a prime:

⟨𝑋⟩ 𝜙 = 1 2𝜋 ∫ 2𝜋 0 𝑋𝑑𝜙, ( 9
)
𝑢 𝜙 = ⟨𝑢 𝜙 ⟩ 𝜙 + 𝑢 ′ 𝜙 = 𝑈 𝜙 (𝜌, 𝑡) + 𝑢 ′ 𝜙 (𝜌, 𝜙, 𝑡) (10) 𝑢 𝜌 = ⟨𝑢 𝜌 ⟩ 𝜙 + 𝑢 ′ 𝜌 = 𝑈 𝜌 (𝜌, 𝑡) + 𝑢 ′ 𝜌 (𝜌, 𝜙, 𝑡), ( 11
)
𝜁 = ⟨𝜁 ⟩ 𝜙 + 𝜁 ′ = 𝜕 𝜌 (𝜌𝑈 𝜙 )∕𝜌 + (𝜕 𝜌 (𝜌𝑢 ′ 𝜙 ) -𝜕 𝜙 𝑢 ′ 𝜌 )∕𝜌 (12)
In the following, all the variables denoted with a tilde are non-dimensional using 1∕𝑓 as the time-scale, and the radius of the domain 𝑅 as the length-scale. For the fluctuations (non-axisymmetric modes), we solve the 2D barotropic vorticity equation and associated Poisson equation for the streamfunction:

𝜕 ζ 𝜕 t +  ( q, ψ) - 𝐸 1∕2 𝑅 2 h ∇ ψ ⋅ ∇ q = 𝐸 𝑅 2 ∇2 ζ - 𝐸 1∕2 𝑅 2 h ζ ( ζ + 1) + F , ( 13
) ζ = - 1 h ∇2 ψ + 1 h2 ∇h ⋅ ∇ ψ + 𝐸 1∕2 𝑅 h2  ( h, ψ), ( 14
)
where q = ( ζ +1)∕ h is the potential vorticity, ψ is a modified streamfunction defined in Appendix C, 𝐸 𝑅 = 𝜈∕(Ω𝑅 2 ) = (ℎ 0 ∕𝑅) 2 𝐸 is the Ekman number based on the radius, and F is a forcing term described below.  is the non-dimensional Jacobian operator in cylindrical coordinates

 (𝑎, 𝑏) = 1 ρ ( 𝜕𝑎 𝜕 ρ 𝜕𝑏 𝜕𝜑 - 𝜕𝑏 𝜕 ρ 𝜕𝑎 𝜕𝜑 ) . ( 15
)
Equations ( 13)-( 14) are solved for all the non-axisymmetric modes on the polar domain with free-slip boundary conditions for the outer and inner boundaries. The free-slip boundary conditions are employed to avoid high computational costs due to the accumulation of strong vorticity near the boundaries, particularly for the lowest 𝐸 cases. We verified that a no-slip boundary condition on the outer boundary leads to qualitatively similar results. For the axisymmetric mode ( Ũ𝜌 , Ũ𝜙 ), we directly solve for the zonal flow evolution equation obtained by performing a zonal average of the zonal component of the Navier-Stokes equations [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF]:

𝜕 Ũ𝜙 𝜕 t + Ũ𝜌 𝜕 Ũ𝜙 𝜕 ρ - Ũ𝜙 Ũ𝜌 ρ + Ũ𝜌 = - ⟨ ũ′ 𝜌 𝜕 ũ′ 𝜙 𝜕 ρ - ũ′ 𝜙 ũ′ 𝜌 ρ ⟩ 𝜙 + 𝐸 𝑅 2 ( ∇2 Ũ𝜙 - Ũ 2 𝜙 ρ2 ) , ( 16
) Ũ𝜌 = 𝐸 1∕2 𝑅 2 h Ũ𝜙 . ( 17
)
The derivation of the previous equations as well as the numerical methods employed to solve them can be found in Appendix C. Note that the simulations are performed with the same fluid height profile as in the experiment, h( ρ) = hmin exp((|𝛽|𝑅∕𝑓 ) ρ), i.e. with the same topographic 𝛽-effect. However, we use the expression of the Ekman pumping over a flat surface, and thereby neglect the bottom boundary curvature. This is justified by the fact that the experimental parameters were carefully chosen such that the bottom topography is as small as possible in amplitude (resulting in a maximum height difference of 5.36 cm and a mean absolute slope of 22%). Furthermore, the greatest slope is located at the center of the bottom plate, which is not resolved numerically. If the curvature of the bottom plate were to be considered, the Ekman pumping term would contain a geometrical factor involving the local slope of the bottom topography [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF][START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF].

For now, we have introduced the forcing as an additional source of vorticity (term F in equation ( 13)). The goal is to reproduce the experimental forcing such that the QG numerical model can be used as a guide and complement the experimental exploration. In the experiment, because of the Coriolis effect, each inlet or outlet generates respectively a small cyclone or anticyclone right above it. This process can be modeled as a stationary source of vorticity in the form of positive or negative Gaussian sources of vorticity of radius 𝓁 𝑓 distributed on a prescribed array. We thus define 𝑁 forcing points distributed over the numerical domain, and at each point, we place a Gaussian source of vorticity such that

F ( x, ỹ) = F0 𝑁 ∑ 𝑖=1 (-1) 𝑖 exp ⎛ ⎜ ⎜ ⎝ - [ x -x𝑖 l𝑓 ] 2 - [ ỹ -ỹ𝑖 l𝑓 ] 2 ⎞ ⎟ ⎟ ⎠ , ( 18
)
where ( x, ỹ) are non-dimensional Cartesian coordinates, the pairs ( x𝑖 , ỹ𝑖 ), 𝑖 ∈ [ [1, 𝑁]] are the center of each forcing vortex, l𝑓 is their radius and F0 is the forcing amplitude. These vorticity sources are distributed over a prescribed stationary array. Both polar and Cartesian forcing arrays were tested. For the polar array, we chose the number of forcing rings arbitrarily, then, the vorticity sources are distributed such that the distance between two rings and between two adjacent vortices is approximately the same (this condition cannot be rigorously verified because of the periodicity of the domain in the azimuthal direction). There are as many positive as negative sources on each ring such that the zonally-averaged forcing term is zero by construction (no direct acceleration of the zonal flow). For the Cartesian forcing, we define the array on the square which encloses the circular domain. The number of lines of the Cartesian array is chosen a priori and the sign of the vorticity sources alternates along lines and columns. The array is then cropped to keep only sources fully within the circular domain. In any case, there is the same number of positive or negative vorticity sources such that there is no net angular momentum introduced by our forcing. In terms of amplitude, in the simulations, we work with a uniform amplitude over the domain. Experimentally, a uniform forcing requires to force more strongly at higher radii because the fluid height increases, which is why being able to control each forcing ring independently is important.

The parameters of the simulations discussed in the present paper are summarized in Table 4. We discuss one simulation with a polar forcing (QG1a) designed to match the experiment, two complementary simulations at smaller Ekman (QG1b,c) and three simulations with a Cartesian forcing and a progressively reduced distance between the vorticity sources to decrease the forcing scale compared to the experiment (QG2,3,4). In the following, all the numerical results are shown in dimensional form using the experimental parameters to allow for a direct comparison with the experimental measurements.

Table 4

Parameters of the quasi-geostrophic simulations used to complement the experimental results. 𝑁 is the number of vorticity sources, Δ 𝑓 is the average spacing between vorticity sources, 𝓁 𝑓 is their radius, 𝐹 0 is the forcing amplitude (equation (18)) and 𝐸 𝑅 the Ekman number based on the radius. For the Cartesian cases, 𝑁 is the number of sources on the square in which the disk in inscribed, and the number in parenthesis gives the exact number of sources on the disk. The experimental parameters are given as a reference, Δ 𝑓 and 𝓁 𝑓 have been estimated by measuring the vorticity above forcing injection and suction points at the very beginning of an experiment (see Figure 9 in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]).

Label Forcing pattern

𝑁 Δ𝑓 (Δ 𝑓 ) l𝑓 (𝓁 𝑓 ) F0 (𝐹 0 ) 𝐸 𝑅 Exp Polar 128 (6 rings) 1.4 × 10 -1 (7 cm) 5.6 × 10 -2 (2.8 cm) 2 × 10 -3 (0.5 𝑠 -2 ) 5.3 × 10 -7 QG1a Polar 540 (12 rings) 6.9 × 10 -2 (3.4 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 1.2 × 10 -7 QG1b Polar 540 (12 rings) 6.9 × 10 -2 (3.4 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 6 × 10 -8 QG1c Polar 540 (12 rings) 6.9 × 10 -2 (3.4 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 2 × 10 -8 QG2 Cartesian 24 2 (408) 8.3 × 10 -2 (4.2 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 1.2 × 10 -7 QG3 Cartesian 48 2 (1716) 4.2 × 10 -2 (2.1 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 6.1 × 10 -8 QG4 Cartesian 72 2 (4404) 2.8 × 10 -2 (1.4 cm) 1.1 × 10 -2 (5.4 mm) 5 × 10 -3 (1.2 𝑠 -2 ) 4.1 × 10 -8

Qualitative observations

Experimental flows

After reaching solid-body rotation, the forcing is turned on to reach uniform rms velocity fluctuations over the domain. We work with forcing amplitudes above the threshold of the transition identified in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. In this regime, zonal jets progressively emerge from the forced turbulent flow, with a transient lasting about 1000 rotation periods (800s), as seen on the space-time diagram of Fig.12(a). During this transient, the jets, which initially emerge at the scale of the forcing rings, drift and merge to self-organize at a larger scale. A statistically-steady state is then achieved where three to one prograde jets are observed depending on the forcing amplitude. Fig. 3(a,c) show two instantaneous velocity fields as measured by PIV once the statistically steady state is achieved. Both experiments are in Regime II but correspond to different forcing amplitude, leading to 𝑅𝑒 = 10.8 × 10 3 and 19.7 × 10 3 respectively (experiments J and A of Table 3). Qualitatively, our experiments performed at higher Reynolds number lead to broader and intensified zonal jets, as can be seen in Fig. 3(b,d) which shows the zonal flow profile for the same experiments. Since the prograde jets width increase with the forcing amplitude, for our most extreme experiments, only one prograde jet can fit in the experimental domain, whereas two to three prograde jets can be obtained for the smallest 𝑅𝑒. As the forcing is increased, we also notice that the prograde jets become less "wavy". Our hypothesis is that for moderate 𝑅𝑒, the flow is close to the transition threshold and thus quasi-resonant, as described in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. Farther from the transition threshold, the resonant Rossby waves are less prone to develop and interact with the jet.

In [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], we show that in a classical quasi-geostrophic model of the experiment, the zonal flow evolution equation ( 16) reduces to

𝜕𝑈 𝜙 𝜕𝑡 = - 1 𝜌 2 𝜕⟨𝜌 2 𝑢 ′ 𝜌 𝑢 ′ 𝜙 ⟩ 𝜙 𝜕𝜌 + . (19)
Here,  represents both bulk and parietal viscous dissipation. The first term on the right-hand-side is the Reynolds stresses divergence, and is responsible for accelerating the zonal flow. In [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], we show that at the very beginning of an experiment, the zonal jets emerge following convergence of prograde momentum at the location of the forcing rings because of the radiation of Rossby waves (see Figure 9(d) in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]).

Once in a fully turbulent, quasi-steady state, we loose this very regular pattern seen in the transient and the associated velocity correlations. However, the zonal jets remain eddy-driven, as illustrated on Figure 4, where we represent the time-averaged radial profile of the Reynolds stresses, for experiment J and the reference simulation QG1a. It is clear that the Reynolds stresses decrease in a prograde jet (prograde momentum convergence), and increase in retrograde flows (prograde momentum divergence).

Independence on the forcing scale

The QG numerical simulations can be used to verify that the zonal flow in Regime II is not sensitive to the detail of the forcing pattern in terms of length scale and distribution (polar versus Cartesian).

The first step was to find a numerical setup for which the obtained solution reproduces well the experimental observations. With this in mind, we performed simulations with a polar forcing pattern. The first important result is that in our QG numerical simulations, we are able to retrieve the transition between the two regimes identified experimentally in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], with a transition between locally and globally forced jets (not shown). As this transition is not the focus of the present paper, we focus only on simulations in Regime II. By slightly tuning the forcing scale and amplitude, we define a reference simulation which reproduces best the experimental flow, labelled QG1a in Table 4. The corresponding space-time diagram is shown in Appendix A, Fig. 12(c), and shows the emergence of zonal jets with properties very close to the experimental ones. Fig. 5 shows the vorticity and azimuthal velocity fields for experiment J and four QG simulations, once in the statistically steady state. Comparing panels (a) and (b) of Fig. 5 shows that the simulation QG1a exhibits three prograde jets of similar amplitude and width compared to the experiment. This reference simulation shows that for the QG simulations to best reproduce the experimental results, we should use forcing patterns with a typical scale smaller than the experimental pattern (540 inlets and outlets arranged on 12 forcing rings in the present case). This is not surprising given that the QG model imposes a purely barotropic and stationary forcing, whereas in the experiment, three-dimensional effects affect the forcing injections and suctions which are not perfect gaussian sources of vorticity aligned along the rotation axis, and which are very likely to generate smaller and more fluctuating structures.

Fig. 5(c) shows the results of simulation QG2 performed with a Cartesian forcing pattern, of approximately the same length scale as in the polar case discussed above, QG1a. The resulting statistically steady state is similar to what is obtained with a polar array. This demonstrates that the polar distribution of the forcing does not significantly influence the late-time jet profile. We note however that whereas the jets do not drift in the experiments and the polar simulation QG1a, the cartesian forcing seems to allow for a slow radial drift of the prograde jets (Fig. 12(e)). This is probably because the forcing cannot be strictly equal on both sides of each jet given its Cartesian distribution. This might also explain the drift reported in [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] where a Cartesian forcing pattern was used in the experiment. We wish to underline that in our cartesian forcing simulations, the jets drift even if the 𝛽-effect is uniform. Hence, the drift mechanism should be captured with models in the 𝛽-plane approximation (see e.g., [START_REF] Cope | The Dynamics of Geophysical and Astrophysical Turbulence[END_REF]). Second, while useful experimentally, it appears that the polar forcing pattern artificially locks the jets at a fixed radial position. Note that we performed experiments where the power of each forcing ring fluctuates randomly around a prescribed mean, and experiments where the forcing alternates between two groups of three rings, but we never observed a long-term drift. Finally, Fig. 5(d,e) show the result of simulations performed with a Cartesian pattern of smaller forcing scale (QG3 and QG4). We recall that for these simulations, we decrease the forcing scale by dividing the distance between vorticity sources by a factor 2 (see Table 4). The decrease in the forcing scale is evident on the vorticity maps which exhibit finer structures, however, qualitatively, the late-time jets amplitude and width do not vary significantly. Note that comparing the number and position of the prograde jets is not robust, given that the system is multistable: as illustrated in Fig. 6, the number and position of the jets can vary for different realisations of the exact same simulation. This multistability is also observed in the experiment. Hence, we argue that the fact that simulation QG4 loses the most central prograde jet is not significant. More generic, idealized simulations and systematics should be performed to make this result more robust. Notably, [START_REF] Scott | Zonal Jet Formation by Potential Vorticity Mixing at Large and Small Scales[END_REF] demonstrated in the framework of potential vorticity mixing that the late-time resulting zonal jets profile is the same whether the forcing is performed at large scale or at small scale, the important parameter being the zonostrophy index 𝑅 𝛽 only. This study could be extended to our geometry, but this is beyond the scope of the present work. Here, our point is that modifying the experimental setup to perform a forcing at smaller scale does not seem to be extremely valuable, since the jets dynamics appears to be unaffected by changes to the forcing scale.

Spectral analysis

Spectral properties and signature of the zonostrophic regime

We now turn to a spectral analysis of the flow observed experimentally and numerically. The goal of this section is to assess whether Regime II is consistent with the predictions in the regime of zonostrophic turbulence [START_REF] Galperin | Geophysical flows with anisotropic turbulence and dispersive waves: Flows with a 𝛽-effect[END_REF], and hence relevant to the regimes of turbulent zonal jets observed on Jupiter [START_REF] Choi | Power spectral analysis of Jupiter's clouds and kinetic energy from Cassini[END_REF]Galperin et al., 2014b). This is not evident a priori given that zonostrophic turbulence was described from simulations of two-dimensional turbulence on the sphere. Whether such regime can be achieved in a fully three-dimensional flow -yet constrained by rotation -still remains to be addressed.

Zonal and residual kinetic energy spectra

In the zonostrophic turbulence regime, the turbulent flow shows profound anisotropy, and the kinetic energy spectra computed from the axisymmetric component of the flow 𝐸 𝑧 (hereafter, zonal spectrum) and from the residual non-axisymmetric component 𝐸 𝑟 (residual spectrum) follow different scalings similar to those in 𝛽-plane turbulence [START_REF] Huang | Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere[END_REF][START_REF] Sukoriansky | Universal Spectrum of Two-Dimensional Turbulence on a Rotating Sphere and Some Basic Features of Atmospheric Circulation on Giant Planets[END_REF][START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF][START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF]. Assuming isotropy for intermediate wave number modes (scales smaller than the transitional scale but larger than the forcing scale), the residual spectrum is expected to be compatible with the inverse cascade branch of the Kolmogorov-Batchelor-Kraichnan theory of two-dimensional turbulence: where 𝑘 is the total wavenumber, and 𝐶 𝐾 ≈ 6 is the universal Kolmogorov-Kraichnan constant [START_REF] Boffetta | Two-Dimensional Turbulence[END_REF]. On the contrary, the zonal spectrum follows a steeper slope:

𝐸 theo 𝑟 (𝑘) ∼ 𝐶 𝐾 𝜖 2∕3 𝑘 -5∕3 , ( 20
)
𝐸 theo 𝑧 (𝑘 𝑟 ) ∼ 𝐶 𝑍 𝛽 2 𝑘 -5 𝑟 , ( 21
)
with 𝑘 𝑟 the wavenumber in the direction orthogonal to the zonal flow, i.e. the radial wavenumber in our case. 𝐶 𝑍 is a constant of order unity [START_REF] Chekhlov | The effect of small-scale forcing on large-scale structures in twodimensional flows[END_REF][START_REF] Huang | Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere[END_REF], whose value was shown to lie around 0.5 by numerical simulations on the sphere [START_REF] Sukoriansky | Universal Spectrum of Two-Dimensional Turbulence on a Rotating Sphere and Some Basic Features of Atmospheric Circulation on Giant Planets[END_REF], around 2 for Jupiter ( Galperin et al., 2014b) and from 1.7 to 3.7 in the experiments of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF]. The theory of zonostrophic turbulence shows that the kinetic energy spectra can provide a useful insight on the exchange of energy between scales as well as a diagnostic tool to determine to what extent a given flow is in the regime of zonostrophic turbulence. In this section, we compare the predictions from this theory with our experimental measurements, complemented by QG simulations.

Given the rotational symmetry of the experimental 2D velocity fields obtained from PIV, we perform a Bessel-Fourier decomposition relevant to the polar system of coordinates [START_REF] Wang | Fourier analysis in polar and spherical coordinates[END_REF], for which details can be found in the Appendix D. Owing to the periodicity in 𝜙, the finite domain in radius, 𝜌 ≤ 𝑅, and the zero-value boundary condition at 𝜌 = 𝑅, the relevant basis functions Ψ 𝑛𝑚 to decompose our fields are separable in polar coordinates and consist in a complex exponential for the angular part, and normalized Bessel functions for the radial part:

Ψ 𝑛𝑚 (𝜌, 𝜙) = 1 √ 2𝜋𝑁 𝑛𝑚 𝐽 𝑚 (𝑘 𝑛𝑚 𝜌)𝑒 𝑖𝑚𝜙 , ( 22
)
with

𝑁 𝑛𝑚 = 𝑅 2 2 𝐽 2 𝑚+1 (𝛼 𝑛𝑚 ). ( 23
)
𝐽 𝑚 is the Bessel function of the first kind of order 𝑚, 𝛼 𝑛𝑚 is the 𝑛 th positive zero of 𝐽 𝑚 , and the radial wavenumber 𝑘 𝑛𝑚 = 𝛼 𝑛𝑚 ∕𝑅 takes discrete values since the domain is bounded. For each Ψ 𝑛𝑚 , 𝑚 is the number of periods in the angular direction, and 𝑛 -1 corresponds to the number of zero crossings in the radial direction. The value of 𝑘 𝑛𝑚 is thus an indication of the scale of the basic patterns, similarly to the normal Fourier transform. Any function, including our velocity fields 𝒖 = (𝑢 𝜌 , 𝑢 𝜙 ) can be decomposed on this basis such that

𝒖(𝜌, 𝜙) = ∞ ∑ 𝑛=1 +∞ ∑ 𝑚=-∞ û𝑛𝑚 Ψ 𝑛𝑚 (𝜌, 𝜙), ( 24
)
where the double hat indicates the double transform. The Bessel-Fourier transform coefficients are

û𝑛𝑚 = ∫ 𝑅 0 ∫ 2𝜋 0 𝒖(𝜌, 𝜙)Ψ * 𝑛𝑚 (𝜌, 𝜙)𝜌𝑑𝜌𝑑𝜙, ( 25
)
where the star denotes the complex conjugate. For our experimental results, the PIV fields are obtained on a regular Cartesian grid. To evaluate the coefficients û𝑛𝑚 = ( û𝜌,𝑛𝑚 , û𝜙,𝑛𝑚 ), we perform a spline interpolation of the velocity field on a polar grid where the radial and azimuthal parts of the transform are separable. Since the velocity field is discrete, the coefficients (25) are computed using fast Fourier transform and discrete Hankel transform algorithms. The angular part is done using the Matlab fft function, and the radial part is performed using the Matlab algorithm provided by Guizar-Sicairos and Gutiérrez-Vega ( 2004). Using the Parseval relation which arises from the orthogonality of the basis functions, kinetic energy spectra can be computed directly in the spectral space, and are proportional to û𝑛𝑚 û * 𝑛𝑚 (see Appendix D). We denote 𝐸 𝑛𝑚 the kinetic energy contained at the wavenumber 𝑘 𝑛𝑚 . We distinguish the kinetic energy contained in the zonal mode characterized by 𝑚 = 0 and denoted 𝐸 𝑧 , and the residual kinetic energy spectra 𝐸 𝑟 , which is the sum of the contribution of all the non-zonal modes, 𝑚 ≠ 0. Fig. 7(a,b) show the zonal and residual kinetic energy spectra for experiments M and B of Table 3. The spectra are computed once the statistically steady state is achieved and time-averaged over 50 statistically-independent spectra spanning 600 rotation times. Exp. B has a Reynolds number about twice higher compared to Exp. M (𝑅𝑒 ∼15,900 versus 7310). The energy spectra provide different information depending whether we are looking at scales smaller or larger than the forcing scale. Experimentally, the forcing scale should lie around the mean distance between two injections on the bottom plate, i.e. Δ 𝑓 ≈ 7 cm (𝑘 𝑓 ∼ 90 rad m -1 ). At scales smaller than the forcing scale (shaded area, 𝑘 > 𝑘 𝑓 ), the residual spectra exhibit a slope close to 𝑘 -4 . This slope is steeper than the -3 slope expected from the direct cascade of enstrophy, but this behavior is commonly observed in both numerical simulations and experiments of 2D turbulence, and may be due either to friction effects, or to the fact that we do not use any logarithmic correction (see [START_REF] Boffetta | Two-Dimensional Turbulence[END_REF], for a review). At scales larger than the forcing scale (𝑘 < 𝑘 𝑓 ). We observe that the residual energy spectra are compatible with a -5∕3 slope consistent with the presence of an inverse cascade of energy. This slope stops as soon as 𝑘 𝛽 is reached, meaning that the isotropic inverse cascade does not continue past the anisotropisation threshold. The dark green and blue lines in Fig. 7(a,b) show the zonal energy spectra. During the transient (thin dark blue line in panel (b)), the zonal energy first peaks close to the forcing scale, meaning that the zonal flow is first locally forced. Once in the statistically steady state and for both experiments M and B, the zonal spectra peak at a larger scale meaning that there is a progressive energy transfer towards scales larger than the forcing scale. As predicted in the zonostrophic regime, the spectrum follows a steep 𝑘 -5 slope. This means that the zonal flow profile is extremely smooth with localised sharp features leading to a rapidly decaying spectral slope. We note that the 𝑘 -5 slope seems to continue below the forcing scale. Even if this was also observed in previous studies [START_REF] Zhang | Beta-plane turbulence: Experiments with altimetry[END_REF][START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF][START_REF] Cabanes | Some statistical properties of three-dimensional zonostrophic turbulence[END_REF], we argue that it is not a robust feature given that it corresponds to scales smaller than the energy injection scale. For experiments M and B, the zonal energy peaks at 𝑘 ∼ 30 rad m -1 and 24 rad m -1 respectively, corresponding to wavelengths of 21 cm and 26 cm which represent the large radial wavelength of the zonal jets. At the largest scales, the zonal 𝑚 = 0 mode contains more energy than all the residual modes. More precisely, the time-averaged total kinetic energy in the flow is of 7.94 × 10 -5 m 2 s -2 in Exp. M and the zonal mode contains 58% of it, whereas for Exp. B, the total kinetic energy is of 3.73 × 10 -4 m 2 s -2 and the zonal mode contains 68% of it. This is significantly larger than our experiments in Regime I, where the zonal flow contains at most ∼ 20% of the total kinetic energy (not shown, see [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF] for more details).

Rate of upscale energy transfer

Previously, we estimated the upscale energy transfer of the inverse cascade, 𝜖 𝐸 , by assuming that dissipation consists in Ekman friction only (equation ( 7)). The amplitude of the theoretical residual spectra can be used to actually measure 𝜖 (equation ( 20)) which is otherwise difficult to do experimentally. To do so, we take the commonly accepted value 𝐶 𝐾 ≈ 6 for the Kolmogorov-Kraichnan constant [START_REF] Boffetta | Two-Dimensional Turbulence[END_REF] and deduce 𝜖 from the fit of the -5∕3 slope on the spectra. The deduced values, denoted 𝜖 𝑆 , are reported in Table 5, and have the correct order of magnitude compared to the estimated 𝜖 𝐸 (Table 3). The precise comparison between the two, shown in Fig. 8(a), shows a small deviation from a pure proportionality relationship. For the zonal spectra, since 𝛽 is known, we report in Table 5 the values of 𝐶 𝑍 corresponding to our fits (equation ( 21)). 𝐶 𝑍 is supposedly universal in the zonostrophic regime. According to the numerical simulations reported in [START_REF] Sukoriansky | Universal Spectrum of Two-Dimensional Turbulence on a Rotating Sphere and Some Basic Features of Atmospheric Circulation on Giant Planets[END_REF]; [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF]; [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF], it should lie around 0.5. Our experimental measurements show that 𝐶 𝑍 ∼ 0.1 -0.3 for all of our experiments despite very different upscale energy transfer rate, which confirms the idea of a universal constant. Its value is smaller for our weakest experiments (L,M,N), probably because these experiments are only slightly super-critical with respect to the transition from Regime I to II described in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. We also note that our value for 𝐶 𝑍 is smaller than that of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] who found 𝐶 𝑍 ∼ 1.7 -3.7. This difference may come from the fact that the 3. The thin dark blue lines correspond to 𝐸 𝑧 computed during the transient of the experiment. For 𝐸 theo 𝑟 , we take 𝐶 𝐾 = 6 and find 𝜖 that best fits the experimental spectra. For 𝐸 theo 𝑧 , since 𝛽 is known, we find 𝐶 𝑍 that best fits the data. (c-e) Spectra derived from QG numerical simulations with decreasing forcing scale compared to the experiment.

experiments of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] are not in the 𝛽-plane framework since 𝛽 is strongly varying with radius. The 𝛽 parameter to use in equation ( 21) is thus uncertain, whereas in our case 𝛽 is uniform and uniquely defined. In addition, the forcing scale in [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] is twice larger than the present one, and close to the transitional scale as we discuss later. Consistently, [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF] report that in the case of a too small scale separation between the forcing scale and the transitional scale, the spectra exhibit non-universal behaviors.

To complement the experiments, we performed the same spectral analysis on the QG simulations. The results are reported in Table 5. The parameters measured on the reference simulation, QG1a, are quantitatively consistent with the spectral parameters measured from the experimental results, for both 𝜖 𝑆 and 𝐶 𝑍 .

Table 5

Parameters measured from the spectral analysis. We use 𝐶 𝐾 = 6 and deduce 𝜖 𝑆 from equation (20) using a fit of the residual spectra for wavenumbers between the forcing wavenumber and the peak energy wavenumber. The -5∕3 slope of the fit is imposed, and the prefactor is determined using Matlab's fit function 𝑎 with least squares fit type. The confidence interval corresponds to the 95% confidence bounds. The constant 𝐶 𝑍 is determined by a fit on the -5 slope of the zonal spectra, since 𝛽 is known (equation ( 21

Zonostrophy index

The intersection of the theoretical -5∕3 (equation ( 20)) and -5 (equation ( 21)) slopes defines a spectral transitional wavenumber 𝑘 𝑆 𝛽 , which can be expressed as

𝑘 𝑆 𝛽 = ( 𝐶 𝑍 𝐶 𝐾 ) 3∕10 ( 𝛽 3 𝜖 𝑆 ) 1∕5 ⏟⏞⏞⏟⏞⏞⏟ 𝑘 𝜖 𝛽 . ( 26
)
This scale thus defines the threshold of spectral anisotropy. In our case, since 𝐶 𝐾 = 6 and 𝐶 𝑍 ∼ 0.1 -0.3, the prefactor (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 is of about 0.29 -0.41. The spectral transitional wavenumber, 𝑘 𝑆 𝛽 is thus less than half that based on the values of 𝜖 𝑆 and 𝛽 only, denoted 𝑘 𝜖 𝛽 . Both 𝑘 𝑆 𝛽 and 𝑘 𝜖 𝛽 are reported in Table 5 for our selected experiments. We additionally report the values of the Rhines wavenumber 𝑘 𝑅 , based on the total rms velocity. The Rhines and transitional scales are also indicated on the experimental spectra, Fig. 7(a,b). Finally, the last two columns of Table 5 report the zonostrophy index with or without using the prefactor on the transitional wavenumber: 𝑅 𝑆 𝛽 = 𝑘 𝑆 𝛽 ∕𝑘 𝑅 and 𝑅 𝜖 𝛽 = 𝑘 𝜖 𝛽 ∕𝑘 𝑅 . The spectral zonostrophy index of our experiments lies between 1.15 and 2.06, whereas without the prefactor on 𝑘 𝛽 , we obtain values between 3.69 and 4.71. The QG simulations again support this result, with 𝑅 𝑆 𝛽 ∼ 1.61 and 𝑅 𝜖 𝛽 ∼ 3.55 for the reference simulation QG1a. The moderate values of 𝑅 𝑆 𝛽 deduced from our spectral analysis question the extent to which our experiments are in the zonostrophic regime. Using 2D barotropic simulations on the sphere, [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF]; [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF] found 𝐶 𝑍 ∼ 0.5 and 𝐶 𝐾 ∼ 6 (prefactor (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 ∼ 0.5), and a lower bound for the zonostrophic regime 5). Horizontal error bars account for the standard deviation when measuring the Thorpe scale.

𝑅 𝑆

𝛽 ≳ 2.5. With values below 2.5, we would then be in the transitional regime between the dissipation-dominated regime and the zonostrophic turbulence. But since we measure 𝐶 𝑍 ∼ 0.1 -0.3, the prefactor is reduced, and the threshold should be rescaled as 𝑅 𝑆 𝛽 ≳ 1.8. In addition, the 2.5 threshold in 𝑅 𝛽 is only indicative (the boundary between the regimes is not a strict proportionality relationship, and implies an offset, see Fig. 13.2 in (Galperin et al., 2019)). We also recall that the boundaries between zonostrophic and friction-dominated regimes have been obtained in the very specific case of a purely two-dimensional flow on the sphere, which significantly differs from our setup. That being said, in [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF]; [START_REF] Sukoriansky | On the Arrest of Inverse Energy Cascade and the Rhines Scale[END_REF], the authors in fact use the spectra as a diagnostic tool to determine the flow regime. In other words, the simulations considered in the zonostrophic regime are those where the predicted spectra are recovered, with a universal behavior and a 𝐶 𝑍 ∼ 0.5 constant. Consistently, we argue that one should rather look at the spectra to determine in which regime we stand, and not the absolute value of the zonostrophy index alone. The important differences between our 3D experiments and their 2D simulations on the sphere is very likely to induce significant differences in the threshold of the zonostrophic regime. The fact that we recover a -5 steep slope for the zonal spectrum, and that the zonal flow contains more energy than the fluctuations at large scales is somewhat more robust and physically meaningful that the absolute value of 𝑅 𝛽 alone. The zonostrophy index is nevertheless useful to compare the degree of zonostrophy of different experiments, as soon as it is computed the same way for each experiment. This is what was done to plot Fig. 1, where we compare various experiments on zonal jets.

Influence of the forcing scale

An important point that we want to discuss is the influence of the forcing scale. [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF] underline that the zonostrophic regime is not defined by a single inequality on 𝑅 𝛽 , but rather a series of inequalities expressed as: 30∕𝑅 ≤ 8𝑘 𝐸 ≤ 2𝑘 𝛽 ≤ 𝑘 𝑓 , in their specific 2D simulations on the sphere, where 𝑅 is the radius of the sphere, and we recall that 𝑘 𝑓 is the forcing wavenumber. The last inequality requires the forcing to act at a scale smaller than half the scale at which the eddies start to feel the 𝛽-effect. Physically, this constraint can be seen as the need for a significant Kolmogorov-Kraichnan inertial range to exist such that an isotropic inverse cascade is able to develop between 𝑘 𝑓 and 𝑘 𝛽 . Practically, when the forcing scale is too large, [START_REF] Galperin | Anisotropic turbulence and zonal jets in rotating flows with a 𝛽-effect[END_REF] report that non-universal behaviors are observed on the spectra, without giving much additional details. In environmental flows (oceans and planetary atmospheres), typically, the forcing acts at a scale smaller than the scale of turbulence anisotropisation 𝐿 𝛽 by a factor 2 to 3 (Galperin et al., 2019, table 13.1). In our experiment, and as discussed previously, we can estimate that the forcing scale is in fact about one to twice 𝐿 𝜖 𝛽 . The forcing is thus already directly influenced by the 𝛽-effect, which is quite clear on the fluid response at the earliest times of our experiments (see Fig. 9 in [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]). The differences in the forcing mechanism and forcing scale may explain the rather non-universal estimates of 𝐶 𝑍 obtained in the past from simulations on the sphere, Cassini measurements of Jupiter's zonal flows and experiments. Second, it is probable that with a rather large forcing scale, we prevent the initial isotropic inverse energy cascade -an anisotropic cascade is present since the jets scale remains larger than the forcing scale. Since we do obtain the scaling expected in the zonostrophic regime, our experiments demonstrate that the scale of the forcing is not such an important parameter for the zonal flow to develop at large scales, and that a steep zonal spectra may develop even if 𝑘 𝑓 ∼ 𝑘 𝛽 . To verify this hypothesis, we compute kinetic energy spectra from simulations QG2,3 and 4. We recall that for these simulations, we decrease the forcing scale by dividing the distance between vorticity sources by a factor 2 (see Table 4). The spectra corresponding to the steady state flows are represented in Fig. 7(c-e). It is clear that as the forcing scale is reduced the scale separation between 𝑘 𝑓 (beginning of the shaded area) and 𝑘 𝛽 (intersection between the residual and zonal spectra) is increased, and a slope consistent with a -5/3 power law appears in that range.

To have a better scale separation between the forcing scale and the transitional scale, another possibility is to decrease the rotation rate, because it would decrease the free surface curvature, decrease the topographic 𝛽-effect, and hence increase 𝐿 𝛽 . We did not explore various rotation rates in the present study because the bottom plate was designed for a rotation rate of 75 RPM only. Second, one should keep in mind that when reducing the 𝛽-effect, the width and distance between the zonal jets will increase, and consequently finite-size effects might become too important. Third, decreasing the rotation rate will increase the Rossby and Ekman numbers, and the system will become less rotationallyconstrained.

To sum up, the important conclusions of our spectral analysis are that:

1. Our experiments in Regime II are consistent with the picture of zonostrophic turbulence, both in terms of slope and amplitude of the theoretical spectra, which confirms their relevance to gas giants applications. This constitutes the first fully-experimental validation of the zonostrophic theory in a completely three-dimensional framework; 2. The absence of scale separation between 𝑘 𝑓 and 𝑘 𝛽 does not impede the establishment of a steep 𝑘 -5 slope for the zonal spectrum. Said differently, the zonal velocity profile obtained is not strongly sensitive to the forcing scale, at least in the range explored in the present study; 3. If one wishes to reach the regime of a well-developed isotropic inverse energy cascade, the goal should be to decrease the forcing scale, which is then a further experimental challenge.

Global and local potential vorticity mixing

We now turn to an analysis of the potential vorticity mixing, and use it to quantify turbulent dissipation rates independently of the spectral analysis presented in the previous section.

Staircasing

In the inviscid limit, the quasi-geostrophic model of the experimental flow reduces to the material conservation of the potential vorticity (PV) 𝑞 defined as

𝑞(𝜌, 𝜙, 𝑡) = 𝜁 (𝜌, 𝜙, 𝑡) + 𝑓 ℎ(𝜌) . ( 27
)
As mentioned in the introduction, in the framework of PV mixing, narrow prograde jets correspond to strong gradients of PV whereas large retrograde jets coincide with regions of weak gradients, leading to the establishment of a PV staircase. In Fig. 9(a,d), we plot instantaneous PV maps for the experiments J and A, which correspond to increasing Reynolds numbers. Because the fluid height ℎ increases exponentially with radius, the background PV profile, 𝑓 ∕ℎ, obtained when the flow is at rest in the rotating frame, is maximum at the center and decreases smoothly towards the edge of the tank. When zonal jets develop, the PV maps show that the PV does not decrease smoothly radially, but instead thin areas of sharp gradients develop. In Fig. 9(a,d), we represent the associated time-averaged and azimuthallyaveraged PV profiles, ⟨𝑞⟩ 𝜙 (𝜌). The steps are more readily visible on these profiles, and each of them corresponds to the presence of a prograde jet, whereas the regions in between are supporting a retrograde flow. Despite significantly increasing 𝑅𝑒, and thus increasing the degree of mixing, the "staircasing" remains moderate even for our most extreme experiments.

Performing a Eulerian zonal average is the most straightforward way of obtaining the mean PV profile as a function of radius. This procedure works well when the jets are straight with no or weak meanders. However, [START_REF] Dritschel | Jet sharpening by turbulent mixing[END_REF] suggest that if the zonal jets are meandering, the zonal average can smear out the locally sharp PV gradients, and thereby underestimate the staircasing. Instead of performing a zonal Eulerian mean, the authors rather suggest to perform an average along the lines of constant PV to represent radial profiles of PV or velocity. Since our zonal jets are sometimes strongly meandering, especially at moderate forcing amplitude, we apply their procedure to compute PV profiles. The method is the following: PV contours are computed from an instantaneous PV map, then each contour is assigned an equivalent latitude (or radius). The equivalent radius for a given PV contour  is defined as the radius of the circular contour which encloses the same area as . We can then plot the PV of each contour as a function of the equivalent radius, denoted 𝑞(𝜌 eq ). This procedure was performed using codes developed by David Dritschel (personal communication). The corresponding results are plotted as pink lines on Fig. 9(b,e). The insert on panel (b) shows that the PV gradients are only very slightly stronger in the prograde jets with this method. Our conclusions are hence unchanged: the degree of staircasing is only moderate in our experiments, even in our most extreme cases.

To provide a more quantitative estimate of the degree of staircaising, we use the method employed by [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]. For each mean PV profile ⟨𝑞⟩ 𝜙 , we define an equivalent staircase, 𝑞 𝑠 such that the area below the curve is preserved, but with constant PV between each jump (see the dashed-dotted lines in Fig. 9(b,e)). Note that the PV jumps are selected as local maxima in the PV radial derivative such that the difference between the real PV gradient and the background PV gradient is larger than an arbitrary value of 50 m -2 s -1 (we use Matlab's findpeaks function). This threshold was chosen arbitrarily such that the identified peaks correspond to dynamically relevant jets, and not to small fluctuations in the zonal velocity profile (see the stars in Figures 10(d-f)). We then compute the integral quantities, 𝐼 1 , 𝐼 2 and 𝐼 3 to quantify how far the actual profile is from the equivalent staircase and the background PV:

𝐼 1 = ∫ 𝑅 0 |⟨𝑞⟩ 𝜙 -𝑞 𝑠 | d𝜌, 𝐼 2 = ∫ 𝑅 0 |𝑞 𝑠 -𝑓 ∕ℎ| d𝜌, 𝐼 3 = ∫ 𝑅 0 |⟨𝑞⟩ 𝜙 -𝑓 ∕ℎ| d𝜌.
The degree of staircaising can be evaluated by the quantities 𝐼 12 = 1 -𝐼 1 ∕𝐼 2 and 𝐼 23 = 𝐼 3 ∕𝐼 2 , which both take the value of 0 if ⟨𝑞⟩ 𝜙 = 𝑓 ∕ℎ (no staircaising) and 1 if ⟨𝑞⟩ 𝜙 = 𝑞 𝑠 (perfect staircaising). For the experiments J and A represented on Figure 9, we obtain 𝐼 12 = 0.14 and 0.12 and 𝐼 23 = 0.15 and 0.16 respectively. Based on Figure 13 in [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF], these values would correspond to a zonostrophy index between approximately 3.5 and 4.5, which is consistent with the zonostrophy index estimated independently for our experiments. To achieve high degrees of staircaising, higher zonostrophy indices are needed, as well as vanishing forcing and vanishing dissipation. Even if strong, instantaneous jets form in our experiment, they are both strongly forced and strongly dissipated. The strong viscous friction constrains us to strongly force the flow in order to achieve sufficiently high Reynolds. It is hence not surprising that the degree of staircaising is weak in the experiments.

The QG simulations can again be employed to verify this hypothesis. Numerically, the viscous dissipation of the flow can be reduced by decreasing the Ekman number, everything else remaining unchanged (something that is not possible in the experiment). Simulations QG1b and QG1c have the exact same parameters as the reference simulation QG1a, but the Ekman number was reduced from 𝐸 𝑅 = 1.2 × 10 -7 down to 6 × 10 -8 and 2 × 10 -8 . Fig. 10 shows instantaneous potential vorticity maps and time-averaged profiles for these three simulations. The inserts show that the slopes in the mean PV profile increase as the Ekman number is decreased. For the degree of staircaising, we obtain respectively 𝐼 12 = [0.17, 0.46, 0.38] and 𝐼 23 = [0.25, 0.54, 0.48]. For the two simulations at smaller Ekman number, the staircases become more pronounced. Once again, these values are consistent with the results reported in [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]. Note that the large dispersion of both 𝐼 12 and 𝐼 23 for the same zonostrophy index (see their Fig. 13) can explain why the lowest Ekman case QG1c have a staircaising degree slightly smaller than the intermediate case QG1b.

Local mixing: the Thorpe scale

Even if strong, global staircaising does not develop, it is possible to quantify the local potential vorticity mixing in each experiment by measuring the equivalent of a mixing length, the so-called Thorpe scale. To this end, we follow the idea developed by Galperin et al. (2014a). We recall that the Thorpe scale is the rms displacement of fluid parcels when the PV profile is sorted into a monotonous profile. We apply this procedure to instantaneous PV profiles along the radius, for all the possible angles 𝜙. Examples of PV profiles before and after sorting are shown in Fig. 9(c,f), and the resulting Thorpe scale is indicated at the top on each panel. The increase in the Thorpe scale and thus in the turbulent mixing is clear when the Reynolds number of the experiment is increased: experiment J (Fig. 9(c)) shows only few deviations from a monotonous PV profile, meaning that overturns are not frequent and the PV is not efficiently mixed, whereas strong deviations are observed for experiment A (Fig. 9(f)). The Thorpe scales measured for every experiment and numerical simulation are reported in Table 5.

The advantage of the Thorpe scale is that it is very easily measured a posteriori, for a given flow, and that it is an indirect measure of turbulence intensity [START_REF] Thorpe | The Turbulent Ocean[END_REF]. For stratified turbulence, with a buoyancy frequency 𝑁, the Thorpe scale can be used to estimate the so-called Ozmidov scale, 𝐿 𝑂 [START_REF] Ozmidov | On the turbulent exchange in a stably stratified ocean[END_REF][START_REF] Thorpe | Turbulence and mixing in a Scottish Loch[END_REF] which corresponds to the scale at which the turbulent turnover time is equal to the internal waves period, i.e. the scale at which the wave effect begins to dominate turbulence:

𝐿 𝑂 = ( 𝜖 𝑁 3 ) 1∕2 . ( 28
)
For stratified turbulence, 𝐿 𝑂 is thus the largest isotropic scale of the system, and is analogous to the transitional scale, 𝐿 𝛽 , for 𝛽-plane turbulence. Hence, if 𝐿 𝑇 allows to estimate 𝐿 𝑂 in stratified turbulence, it may allow to estimate 𝐿 𝛽 in 𝛽-plane turbulence. This is of interest because 𝜖 can be estimated from 𝐿 𝛽 (equation ( 3)). Here, we wish to confront this hypothesis to our experimental measurements, where we have measured 𝜖 and thus 𝐿 𝜖 𝛽 in a completely independent way using the spectral analysis. Fig. 11 compares the Thorpe scale measured in our experiments and QG simulations with the Rhines scale 𝐿 𝑅 and the transitional scales 𝐿 𝜖 𝛽 and 𝐿 𝑆 𝛽 . These figures show that there is a very good correlation between 𝐿 𝑇 and both 𝐿 𝑅 and 𝐿 𝛽 , which holds for the whole range of 𝑅𝑒 explored:

𝐿 𝑇 ≈ 𝑎𝐿 𝜖 𝛽 , 𝐿 𝑇 ≈ 𝑏𝐿 𝑅 , 𝐿 𝑇 ≈ 𝑐𝐿 𝑆 𝛽 , ( 29
)
The fact that the correlation works equally well for those three scales arises from the fact that 𝐿 𝑅 ∼ 𝑅 𝜖 𝛽 𝐿 𝜖 𝛽 , and the zonostrophy index does not vary sufficiently for the two scales to behave independently. In our experimental set, 𝑅 3)). The comparison between the measured 𝜖 𝑆 and 𝜖 𝑇 is represented in Fig. 8(b). There is a non-negligible dispersion around the 1:1 relationship, but the overall tendency is good, suggesting that PV sorting can effectively be used to estimate turbulent energy transfers.

Conclusions and Discussion

We have presented experimental results from a laboratory setup designed to reach extreme regimes of zonal jets relevant to gas giants [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF][START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. Two regimes of turbulent zonal jets have been identified by [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF]. Here, we analyzed the statistical properties of the intense, experimental 𝛽-plane turbulence obtained in the super-resonant Regime II, for large forcing amplitudes. We showed that strong, instantaneous, turbulent zonal jets develop and contain up to 70% of the total kinetic energy of the flow in the most extreme experiments. These experiments are complemented by 2D quasi-geostrophic numerical simulations which allow us to reach the extreme experimental conditions in terms of Reynolds and Ekman numbers, and complement them by addressing various forcing and even more extreme regimes.

Zonostrophic turbulence and potential vorticity mixing

From a spectral analysis of the flow, we showed that the measured turbulent flow shares the properties of the socalled zonostrophic turbulence, relevant to the gas giants, in which the zonal flow alone contains more kinetic energy than the remaining of the flow. Along with [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF], these are the first experiments able to reach such extreme regimes (see Fig. 1). While we retrieve the predicted steep slope for the zonal spectrum, it is hard to argue that a clear inverse cascade of energy develop in the experiments. Thanks to QG simulations, we show that this is due to the relatively large forcing scale of our experiment, which is in particular too close to the transitional scale. In other words, the forced eddies are directly influenced by the 𝛽-effect, and no isotropic inverse cascade of energy can robustly develop. Yet, this absence of inverse cascade does not impede the flow from feeding larger scales in an anisotropic way (i.e. feeding zonal structures), as expected for scales larger than the transitional scale. Interestingly, the forcing scale does not seem to have a significant influence on the final zonal jets profile obtained once in the quasi-steady state. This is consistent with the conclusions of [START_REF] Scott | Zonal Jet Formation by Potential Vorticity Mixing at Large and Small Scales[END_REF], obtained in idealized simulations of potential vorticity mixing, and means that the relatively large forcing scale is not a significant limitation of our current experimental setup to study turbulent zonal jets.

At this point, we would like to raise the question of the "locality" of the zonal jets forcing. In a first paper [START_REF] Lemasquerier | Zonal jets at the laboratory scale: Hysteresis and Rossby waves resonance[END_REF], we showed that the jets developing in our experiment are initially due to wave-mean flow interactions, and emerge via the transport and deposition of momentum by Rossby waves, i.e. a streaming process. If this is particularly relevant for Regime I, not described here, the nature of the transition towards the second regime, a Rossby wave resonance, leads us to the hypothesis that wave-mean flow interactions are probably also of importance in the second regime, even in the highly turbulent cases. The acceleration of the zonal flow by a streaming process is non-local, and implies a direct transfer between waves and the mean flow. This picture, where both local cascades and nonlocal direct transfers coexist, is in line with quantifications of energy transfer in zonostrophic turbulence simulations, where both processes are at play (Galperin et al., 2019).

In terms of potential vorticity mixing, despite the strength and rectilinear shape of the zonal jets in our most extreme experiments, we showed that they are not accompanied by a strong potential vorticity mixing. This is because the zonostrophy index of our experiments is still "moderate" relatively to the cases leading to strong staircaising [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]. Once again, QG simulations show that this is due to experimental limitations, since decreasing the viscous friction allow to reach higher degrees of staircaising. This result is important by itself because it bridges the gap between experiments and asymptotic numerical models. Furthermore, it shows that it is not necessary to be in the asymptotic regime of vanishing friction and forcing to obtain dominant turbulent zonal jets. In other words, strong, instantaneous zonal jets can be obtained even if the background PV gradient is moderately mixed, i.e. even if the relative vorticity associated with the jets is small compared to the background, planetary vorticity. Note that PV mixing studies (e.g. [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF] showed that there is a relation between the jets spacing (the staircase width) and the jets intensity (the PV jump at the interfaces). We could not test this relation in our experiments, because we are confronted with spatial confinement, and because only a small number of jets develops: the measure of the jets spacing is then compromised by finite-size effects. Next, we showed that it is possible to estimate the local intensity of the potential vorticity mixing by measuring the equivalent of the Thorpe scale [START_REF] Thorpe | The Turbulent Ocean[END_REF]. We verified experimentally that the measured Thorpe scale is proportional to the transitional scale, and that this relation can be used to deduce the upscale energy transfer rate. We showed that we obtain a good estimate of the upscale energy flux measured using spectral analysis. This result underlines the consistency between the zonostrophic turbulence and potential vorticity mixing theories.

Advantages and drawbacks of different methods for estimating the upscale turbulent energy transfer rate

Our experiments allow to address how efficiently kinetic energy spectra and PV mixing can be used to quantify the upscale energy transfer rate 𝜖 of a given turbulent flow with a 𝛽-effect. In the present paper, we have used three different methods:

1. The first and simplest estimate that can be made is based on the assumption that, in the experiment, the energy is dissipated by Ekman friction at a rate 𝛼 = Ω𝐸 1∕2 , and hence 𝜖 can be estimated simply by measuring the total kinetic energy of the flow: 𝜖 𝐸 ∼ 𝑢 2 rms 𝛼∕2. 2. The second method, and perhaps the most rigorous one, consists in computing the residual kinetic energy spectra and fitting the -5∕3 slope (equation ( 20)), 𝑠, which is proportional to 𝜖 2∕3 : 𝜖 𝑆 ∼ (𝑠∕𝐶 𝐾 ) 3∕2 . 3. The third method consists in measuring a Thorpe scale by sorting local and instantaneous PV profiles. The deduced scale, assuming that it is correlated with 𝐿 𝛽 , allows to retrieve 𝜖 thanks to the relation 𝜖 𝑇 ∼ 𝛽 3 (𝐿 𝑇 𝛽 ) 5 . Method 1 gives a good order of magnitude for 𝜖, but we did not find an exactly linear relationship between 𝜖 𝐸 and 𝜖 𝑆 (Fig. 8(a)), probably because bulk viscous dissipation and side friction are not completely negligible in our experiment. In addition, the main dissipation mechanism in natural flows is far from being a simple Ekman friction, and this is rather a practical scaling for idealized experimental or numerical model than a true generic method. While this first method assumes that dissipation is due to bottom viscous boundary layers only, the two other methods do not make any hypothesis on the nature of the dissipation process. Method 2, which we consider as the reference one, is probably the most robust but requires to have high resolution spatial and temporal observations to be able to compute energy spectra. For method 3, we show that the Thorpe scale is strongly correlated with the transitional scale, suggesting that the turbulent transfer rate could be efficiently retrieved in natural flows using potentially simpler measurements of instantaneous potential vorticity profiles. The success of this method however requires a precise estimate of the scaling prefactor for the relationship between 𝐿 𝑇 and 𝐿 𝛽 . Here, we found 𝐿 𝛽 ∼ 2.1𝐿 𝑇 , but this factor may potentially vary. We show nevertheless that this scaling holds for the whole range of Reynolds numbers explored. A complete systematic study, where 𝐸 and 𝑅 𝛽 are varied independently would allow to determine the validity of this scaling depending on the flow regime. Note that the Thorpe scale is routinely used in oceanography to measure the vertical mixing of stratified turbulence. Our results support the analogy raised by Galperin et al. (2014a), i.e. that an equivalent method can be used to measure the lateral mixing due to 𝛽-plane turbulence, as done by [START_REF] Cabanes | Revealing the Intensity of Turbulent Energy Transfer in Planetary Atmospheres[END_REF] for Jupiter and Saturn.

Implication for Jupiter's cloud layer dynamics

What are the implications of our results regarding the dynamics observed in the cloud layer of Jupiter? We recall that the theory of zonostrophic turbulence was developed in a purely two-dimensional framework, thanks to numerical simulations on the sphere, i.e. indirectly assuming that Jupiter's turbulence is shallow and confined in its weather layer. It is thus not obvious a priori that this theory is still relevant in the case of Jupiter and Saturn jets, which are now strongly believed to be deep, extending as deep as 3,000 and 9,000 km respectively [START_REF] Kaspi | Jupiter's atmospheric jet streams extend thousands of kilometres deep[END_REF][START_REF] Galanti | Saturn's Deep Atmospheric Flows Revealed by the Cassini Grand Finale Gravity Measurements[END_REF][START_REF] Kaspi | Comparison of the Deep Atmospheric Dynamics of Jupiter and Saturn in Light of the Juno and Cassini Gravity Measurements[END_REF]. Quantitative analyses of the Cassini 70-days movie of Jupiter's clouds support the picture of an inverse energy cascade in Jupiter's weather layer [START_REF] Choi | Power spectral analysis of Jupiter's clouds and kinetic energy from Cassini[END_REF]Galperin et al., 2014b). Using the same data, [START_REF] Young | Forward and inverse kinetic energy cascades in Jupiter's turbulent weather layer[END_REF] showed that there is indeed an inverse transfer of energy from scales of 2,000-3,000 km, up to 20,000 km. However, they also identify a direct transfer of energy from 2,000 km towards small scales, which is not expected in classical 2D turbulence. In addition, they show that most energy transfer occurs between eddies and the zonal flow directly, i.e. non-locally. The scale of transition between direct and inverse energy transfer is close to the first Rossby radius of deformation, which supports the idea that the forcing is due to baroclinic instabilities. However, one should keep in mind that the velocity fields from which these conclusions are drawn are measured in the weather layer only. It is possible that the forcing of the jets at depth occurs via small scale convective motions, while baroclinic processes dominate in the shallow weather layer, with coupling mechanisms unexplored to date. Note that dual cascades to small scales and large scales have been reported for several physical systems, such as magnetohydrodynamics in two and three space dimensions, surface capillary waves and rotating-stratified turbulence (see [START_REF] Pouquet | Dual constant-flux energy cascades to both large scales and small scales[END_REF] for an overview, and [START_REF] Galperin | Quasinormal scale elimination theory of the anisotropic energy spectra of atmospheric and oceanic turbulence[END_REF]).

In the present experimental setup, the zonal jets are barotropic jets extending parallel to the rotation axis in a deep layer of water, which is hence representative of the deep jets configuration contrary to previous shallow water experiments. Our setup allows to model a fully three-dimensional system, without any a priori assumption on the twodimensionality of the flow. The fact that jets can emerge spontaneously from a deep, fast rotating, fully-3D turbulent flow, even in the presence of a large bottom drag (relatively to that on Jupiter), supports the deep hypothesis of zonal jets origin. Furthermore, the spectral analysis shows that, even for deep jets, the spectral properties derived in the framework of two-dimensional turbulence can be recovered. Therefore, the fact that these scalings are also measured on Jupiter [START_REF] Choi | Power spectral analysis of Jupiter's clouds and kinetic energy from Cassini[END_REF]Galperin et al., 2014b) does not allow to discriminate between shallow and deep jets. One possible line of investigation would be to continue the spectral analysis to confirm or not the presence of a turbulent inverse cascade of energy by computing structure functions. This would also allow for a direct comparison with the aforementioned turbulent statistics of Jupiter, and in particular, we could verify if we also retrieve a direct energy cascade despite the quasi-two-dimensionality of our fast-rotating experimental flows. This analysis was out of reach with our current experimental datasets because better time and space-resolved PIV measurements are required. Finally we wish to conclude by underlining that if Jupiter's zonal jets are deep and emerge from deep thermal convection, it is possible that the turbulence observed in the cloud layer is somehow "superimposed" to a deeper turbulence from which the zonal jets emerge. In such a case, comparing the residual and zonal kinetic energy spectra should be taken with caution given that we may be looking at the superposition or coexistence of flows having different physical origins.

As already mentioned in the introduction, we want to underline that exploring barotropic experimental or numerical models is still relevant for systems which exhibit a predominantly baroclinic circulation, such as oceans and atmospheres. Indeed, planetary, quasi-geostrophic flows can undergo a phenomenon of barotropization. This process refers to the fact that following baroclinic instabilities, the energy initially cascades from the scale of meridional temperature gradients down to the radius of deformation (commensurate with the most unstable wavenumber). There is subsequently a conversion to barotropic energy, because the energy tends to feed the gravest mode, which is a barotropic (vertically-invariant) mode. Finally, an inverse cascade to larger modes occur [START_REF] Rhines | The dynamics of unsteady currents[END_REF][START_REF] Salmon | Two-layer quasi-geostrophic turbulence in a simple special case[END_REF]. Said differently, in QG flows, the energy cascade to larger horizontal scales is generally accompanied by a cascade to larger vertical scales, and the barotropic modes can end up containing a significant fraction of the energy (see [START_REF] Charney | Geostrophic Turbulence[END_REF] and Chapter 12 in [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF]). For the Earth, both GCM simulations and reanalysis of atmospheric and oceanic data confirmed this tendency to strong barotropization (see Galperin et al., 2019, and references therein). In baroclinic flows, the barotropization occurs at a scale close to the internal radius of deformation, which should be small enough to leave room for an inverse cascade in the barotropic modes. In our experiment, instead of relying on baroclinic instabilities to inject energy, we directly (mechanically) force the flow at a small and well-controlled scale, and the barotropic flow grows by itself because of the quasi-geostrophy of the system.

Finally, in the present manuscript, we focus on gas giants dynamics at midlatitudes, where the large variation of the Coriolis force with radius is responsible for a strong 𝛽-effect and hence strong zonal flows. At increasing latitudes, the gradient of the Coriolis force decreases and tends towards zero at the pole. Consequently, the zonal jets weaken, and the dynamics becomes more isotropic, with a predominance of vortices rather than jets. In addition, the vortices become predominantly cyclonic, contrary to what is observed at midlatitudes. This transition in the dynamics with increasing latitude has been first suggested by [START_REF] Theiss | Equatorward Energy Cascade, Critical Latitude, and the Predominance of Cyclonic Vortices in Geostrophic Turbulence[END_REF] for Earth's oceans, and later observed in shallow-water turbulence models on the sphere (e.g. [START_REF] Cho | The emergence of jets and vortices in freely evolving, shallow-water turbulence on a sphere[END_REF]; [START_REF] Scott | Forced-Dissipative Shallow-Water Turbulence on the Sphere and the Atmospheric Circulation of the Giant Planets[END_REF]). When the Juno spacecraft arrived at Jupiter, visible and infrared observations from above the poles revealed an incredible dynamics consisting in persistent polygonal patterns of cyclones around both poles [START_REF] Adriani | Clusters of cyclones encircling Jupiter's poles[END_REF]. A natural extension of this work is hence to look at how the dynamics evolve for a progressively less dominant 𝛽-effect.

Advantages and limitations of experimental approaches of zonal jets

We wish to conclude this study by a broader discussion on the approach and methods employed in the prospect of better understanding gas giants dynamics. The regimes reached experimentally (𝐸 ∼ 3 × 10 -7 , 𝑅𝑒 ∼ 10 4 , 𝑅𝑜 ∼ 10 -3 ) are not impossible to reach in direct numerical simulations, but at very large cost. To give an order of magnitude, a single DNS of 1,000 rotation times at an Ekman number three times larger than in the experiment (𝐸 ∼ 10 -6 ) represents 13 minutes of an experiment, but would require 13 days of computation on 2048 CPU cores (650,000 CPU hours). This type of DNS was performed in [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF], but it is clear that numerical systematic studies in these regimes are inconceivable. The high cost of these DNS is inherent to 3D simulations of geostrophic turbulent flows, where both large-scale structures and small-scale turbulent eddies and inertial waves are present and need to be resolved simultaneously. On the contrary, one experimental realisation "costs" about 2 days (one for the actual experiment, and one for saving and post-processing the images through the PIV algorithm), which allows for multiple realisations and exploration of the parameter space. In addition, the dynamics of the large-scale jets is slow, and results from cumulative effects from the underlying turbulence. Studying their long-term dynamics requires to wait for very long times. Experimentally, we can easily reach several thousands of rotation times, and time-resolved particle image velocimetry allows for high-resolution records of the interactions between the turbulent flow and the slowly evolving large-scale jets.

Of course, idealized numerical models allow to circumvent these difficulties. This is the case for instance of reduced two-dimensional models, such as shallow-water models and quasi-geostrophic models including the one used in the present study. Statistical simulations and quasilinear models (e.g. [START_REF] Constantinou | Emergence and Equilibration of Jets in Beta-Plane Turbulence: Applications of Stochastic Structural Stability Theory[END_REF] where eddy-eddy interactions are neglected can also be used. Recently, rare events algorithm have also been employed to study multistability and spontaneous transitions among zonal jets [START_REF] Bouchet | Rare Event Algorithm Links Transitions in Turbulent Flows with Activated Nucleations[END_REF]. However, the assumptions underlying each of these models and their relevance for gas giants dynamics need to be systematically addressed.

Despite their clear advantages, experiments also come with their own limitations, as underlined throughout this paper:

• The spatial confinement of the experiment is problematic for studying zonal jets equilibration. For instance, finite-size effects "discretize" the evolution of jets spacing when varying any control parameter. Due to confinement, jets are not free to evolve in space, and this may for instance impede long-term drift or nucleation of jets. One could nevertheless argue that jets on the planets are also confined, but there are still about ten prograde jets on Jupiter, leading to a scale separation between the size of the domain and the jets. We should hence seek to build experiments where the scale of the tank is large compared to the scale of the jets, which itself is large compared to the scale of injection. This is challenging, because when trying to reach more extreme regimes, we increase the turbulence intensity or decrease the Ekman friction, which increases the jets scale. Imposing a fast rotation on larger scale containers is also a technical challenge.

• Experiments at high rotation rates allow to reach regimes with small friction (here 𝐸 ∼ 3 × 10 -7 ), but never asymptotically small. We stand in a strongly forced-dissipative regime, and it is difficult to bridge the gap with idealized models where both the forcing and the dissipation are vanishing, as discussed in the case of potential vorticity mixing.

• The experimental forcing is performed at small-scale, but it is difficult to reach the same scale separation as in numerical simulations and less straightforward to change the forcing properties. Exploring numerically the effect of the forcing properties (scale, spatio-temporal stationarity...) is hence required to extrapolate experimental results to more realistic planetary conditions, as proposed here using QG simulations.

• Finally, a last difficulty is that we cannot vary independently the Ekman and Reynolds numbers and the zonostrophy index of the flow. Changing the rotation rate modifies the Ekman number but also the zonostrophy index because it changes the free surface shape. One way of avoiding this is to use a sloping bottom and a rigid lid instead of a free surface for the 𝛽-effect, but it introduces a supplementary friction which can kill the zonal flow.

These limitations justify that we employed idealized numerical simulations to reproduce the experimental conditions and then explore one by one some effects artificially introduced by the experimental constraints (forcing nature, spatial confinement, high viscosity, boundary conditions, etc). Note that it is also important to consider the addition of physical effects which will hardly be incorporated in experiments, such as magneto-hydrodynamical dissipation of zonal flows or compressibility effects.

Finally, let us underline that the goal of the idealized experimental and numerical models presented here is to shed light on fundamental mechanisms governing the dynamics of the system, rather than be quantitatively predictive regarding the specific case of Jupiter. Jupiter exhibits a wealth of dynamical processes, occurring at very different temporal and spatial scales. Our understanding of Jupiter and the gas giants as global systems is still complicated by numerous sources of uncertainty and technical limitations. The continuous improvement of technical and computational capabilities allows experimental and numerical models to get closer to the planetary regimes. However, understanding Jupiter by forward modelling requires observations of sufficiently good quality and coverage to which the models outputs can be compared. The recent, accurate observations of Jupiter and Saturn from the Juno and Cassini missions constitute a very good opportunity in this regard, while introducing new and exciting challenges for planetary modellers at the same time.

B. Data from previous experiments

In Table B, we provide the dimensional and non-dimensional parameters that are used to plot Fig. 1, where the regimes reached by experiments on zonal jets are compared. Note that various forcing types are employed. The definition of the Ekman number 𝐸, Reynolds number 𝑅𝑒 and and zonostrophy index 𝑅 𝛽 are those of Table 1.

C. Quasi-geostrophic model of the experiment C.1. Derivation

We use the cylindrical coordinates (𝜌,𝜙,𝑧) with 𝑧 oriented downward and (𝐞 𝜌 , 𝐞 𝜙 , 𝐞 𝑧 ) the associated unit vectors (Fig. 2). We consider the of an incompressible fluid of constant kinematic viscosity 𝜈 and density 𝜌 𝑓 , rotating around the vertical axis at a constant rate 𝛀 = Ω 𝐞 𝑧 . In our experimental setup, Ω > 0 since the turntable rotates in the clockwise direction. We denote the velocity field 𝐮 = (𝑢 𝜌 , 𝑢 𝜙 , 𝑢 𝑧 ) 𝐞 𝜌 ,𝐞 𝜙 ,𝐞 𝐳 . The fluid is enclosed inside a cylinder of radius 𝑅. The lower boundary is a rigid plate located at 𝑧 = 0 and the upper boundary is a free surface defined by 𝑧 = -ℎ(𝜌). Note that here we assume that our experiment, which has a parabolic free-surface and a curved bottom, can be modeled with a flat bottom and an exponential free-surface. Doing so, we neglect the influence of the shape of the bottom topography on the vertical velocity (see equation ( 37)). The experimental parameters were carefully chosen such that the bottom topography is as small as possible in amplitude (resulting in a maximum height difference of 5.36 cm and a mean absolute slope of 22%). Thus, one should keep in mind that the presently derived model is only valid for relatively smooth bottom topographies for which we can use the expression of the Ekman pumping over a flat surface. If the topography was of high amplitude, then the local inclination of the bottom boundary would enter the QG model because it modifies the Ekman pumping [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF].

We start from the continuity and horizontal Navier-Stokes equations in the rotating frame:

𝜕𝑢 𝜌 𝜕𝑡 + 𝑢 𝜌 𝜕𝑢 𝜌 𝜕𝜌 + 𝑢 𝜙 𝜌 𝜕𝑢 𝜌 𝜕𝜙 - 𝑢 2 𝜙 𝜌 -𝑓 𝑢 𝜙 = - 1 𝜌 𝑓 𝜕𝑃 𝜕𝜌 + 𝜈 ( ∇ 2 𝑢 𝜌 - 𝑢 𝜌 𝜌 2 - 2 𝜌 2 𝜕𝑢 𝜙 𝜕𝜙 ) , ( 30
) 𝜕𝑢 𝜙 𝜕𝑡 + 𝑢 𝜌 𝜕𝑢 𝜙 𝜕 𝜌 + 𝑢 𝜙 𝜌 𝜕𝑢 𝜙 𝜕𝜙 + 𝑢 𝜙 𝑢 𝜌 𝜌 + 𝑓 𝑢 𝜌 = - 1 𝜌 𝑓 1 𝜌 𝜕𝑃 𝜕𝜙 + 𝜈 ( ∇ 2 𝑢 𝜙 - 𝑢 𝜙 𝜌 2 + 2 𝜌 2 𝜕𝑢 𝜌 𝜕𝜙 ) , ( 31
) 1 𝜌 𝜕(𝜌𝑢 𝜌 ) 𝜕𝜌 + 1 𝜌 𝜕𝑢 𝜙 𝜕𝜙 + 𝜕𝑢 𝑧 𝜕𝑧 = 0, (32)
where

∇ 2 ⋅ = 𝜕 2 𝜌 ⋅ +𝜕 2 𝜙 ⋅ ∕𝜌 2 + 𝜕 𝜌 ⋅ ∕𝜌. The Coriolis parameter is 𝑓 = 2Ω and 𝑃 = 𝑝 + 𝜌 𝑓 𝑔𝑧 -𝜌 𝑓 𝑓 2 𝜌 2 ∕8
is the reduced pressure incorporating the gravity and centrifugal effects. Note that if we neglect the vertical dependence of the horizontal velocity, we keep it for the vertical velocity 𝑢 𝑧 . Indeed, as previously explained, 𝑢 𝑧 is expected to strongly vary close to the top and bottom boundaries, and we want to take into account these effects on the horizontal velocity divergence.

The curl of the Navier-Stokes equation leads to the vorticity equation

𝜕𝜁 𝜕𝑡 + 𝑢 𝜌 𝜕𝜁 𝜕𝜌 + 𝑢 𝜙 𝜌 𝜕𝜁 𝜕𝜙 + (𝜁 + 𝑓 ) 𝛁 ℎ ⋅ 𝒖 = 𝜈∇ 2 𝜁 , ( 33
)
where 𝜁 = (𝛁 × 𝒖) ⋅ 𝒆 𝑧 = (𝜕 𝜌 (𝜌𝑢 𝜙 ) -𝜕 𝜙 𝑢 𝜌 )∕𝜌 is the vertical component of the vorticity and 𝛁 ℎ ⋅ 𝒖 is the horizontal divergence

𝛁 ℎ ⋅ 𝒖 = 1 𝜌 𝜕(𝜌𝑢 𝜌 ) 𝜕𝜌 + 1 𝜌 𝜕𝑢 𝜙 𝜕𝜙 . ( 34
)
The last term of the left hand side of equation ( 33), the vortex stretching term, involves the horizontal divergence of the flow which can be estimated from equation (32) after integration from 𝑧 = -ℎ(𝜌) to 𝑧 = 0 (𝑧 oriented downward) to unveil the Ekman pumping through the vertical velocity:

𝛁 ℎ ⋅ 𝒖 = - 1 ℎ(𝜌) ∫ 0 𝑧=-ℎ 𝜕𝑢 𝑧 𝜕𝑧 d𝑧 = 𝑢 𝑧 | | |𝑧=-ℎ -𝑢 𝑧 | | |𝑧=0 ℎ(𝜌) . ( 35
)
The vertical velocity at the free surface 𝑢 𝑧 | | |𝑧=-ℎ is given by the kinematic condition

𝑢 𝑧 | | |𝑧=-ℎ = - ( 𝜕ℎ 𝜕𝑡 + 𝑢 𝜌 𝜕ℎ 𝜕𝜌 + 𝑢 𝜙 𝜌 𝜕ℎ 𝜕𝜙 ) = -𝑢 𝜌 𝜕ℎ 𝜕𝜌 , ( 36
)
since ℎ is axisymmetric and we neglect any temporal variations of the fluid height (rigid lid approximation). The vertical velocity at the bottom 𝑢 𝑧 | | |𝑧=0 results from the no-slip boundary condition generating an Ekman pumping. According to linear Ekman theory, for a flat bottom and small Rossby number, the vertical velocity at the top of the boundary layer is proportional to the relative vorticity in the interior flow (see section 5.7 in [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF]:

𝑢 𝑧 | | |𝑧=0 = - 1 2 𝛿𝜁 = - 1 2 𝐸 1∕2 ℎ 0 𝜁 , ( 37
)
where 𝛿 = √ 2𝜈∕𝑓 is the thickness of the Ekman layer and 𝐸 = 2𝜈∕(𝑓 ℎ 2 0 ) is the Ekman number, ℎ 0 being the mean fluid height. The horizontal divergence ( 35) is then

∇ ℎ ⋅ 𝐮 = - 𝑢 𝜌 ℎ dℎ d𝜌 + 𝐸 1∕2 2 ℎ 0 ℎ 𝜁 . ( 38
)
The squeezing and stretching of vorticity is hence due to both the changes in the fluid depth and the vertical velocity induced by the Ekman boundary layer.

= 𝜈∇ 2 𝜁 , (39)
We switch to non-dimensional variables using 1∕𝑓 as the timescale and the radius of the tank, 𝑅, as the length-scale, and we denote the non-dimensional variables with a tilde such that

𝜁 = ζ 𝑓 , ( 40
) 𝑢 𝜌 = ũ𝜌 𝑓 𝑅 (41) 𝑢 𝜙 = ũ𝜙 𝑓 𝑅 (42) 𝜌 = ρ𝑅 (43) 𝑡 = t∕𝑓 (44) ℎ = h𝑅 (45)
We keep the tilde in the following to better identify non-dimensional variables, and to avoid confusion with the experimental variables and parameters which are always given in dimensional forms first. Equation (39) becomes

𝜕 ζ 𝜕 t + ũ𝜌 𝜕 ζ 𝜕 ρ + ũ𝜙 ρ 𝜕 ζ 𝜕𝜙 - ũ𝜌 h d h d ρ ( ζ + 1) + 𝐸 1∕2 𝑅 2 h ( ζ + 1) ζ = 𝐸 𝑅 2 ∇2 ζ + F , ( 46
)
where 𝐸 𝑅 is the Ekman number based on the radius of the tank, 𝐸 𝑅 = 2𝜈∕(𝑓 𝑅 2 ) = (ℎ 0 ∕𝑅) 2 𝐸, and we have introduced a forcing term F . To close this equation, we now need an expression for the horizontal components of the velocity. To do so, we use the definition of 𝜁 to rewrite the expression of the horizontal divergence (equation ( 38)) as a zero-divergence for a modified velocity field:

1 ρ 𝜕( ρ ũ𝜌 ) 𝜕 ρ + 1 ρ 𝜕 ũ𝜙 𝜕𝜙 = - ũ𝜌 h 𝜕 h 𝜕 ρ + 𝐸 1∕2 𝑅 2 h ζ (47) ⟹ 𝜕 𝜕 ρ ( ρ [ h ũ𝜌 -ũ𝜙 𝐸 1∕2 𝑅 2 ]) + 𝜕 𝜕𝜙 ( h ũ𝜙 + ũ𝜌 𝐸 1∕2 𝑅 2 ) = 0. (48)
This allows us to define a streamfunction ψ such that

h ũ𝜌 -ũ𝜙 𝐸 1∕2 𝑅 2 = 1 ρ 𝜕 ψ 𝜕𝜙 , ( 49
) h ũ𝜙 + ũ𝜌 𝐸 1∕2 𝑅 2 = - 𝜕 ψ 𝜕 ρ , ( 50
)
or equivalently

ũ𝜌 = 1 h 1 1 + 𝐸 𝑅 h-2 ( 1 ρ 𝜕 ψ 𝜕𝜙 - 𝐸 1∕2 𝑅 2 h 𝜕 ψ 𝜕𝜌 ) = 1 h ( 1 ρ 𝜕 ψ 𝜕𝜙 - 𝐸 1∕2 𝑅 2 h 𝜕 ψ 𝜕 ρ ) + (𝐸 𝑅 h-2 ), (
) ũ𝜙 = 1 h 1 1 + 𝐸 𝑅 h-2 ( - 𝐸 1∕2 𝑅 2 h ρ 𝜕 ψ 𝜕𝜙 - 𝜕 ψ 𝜕 ρ ) = 1 h ( - 𝐸 1∕2 𝑅 2 h ρ 𝜕 ψ 𝜕𝜙 - 𝜕 ψ 𝜕 ρ ) + (𝐸 𝑅 h-2 ), 51
where we have neglected terms of order greater or equal to (𝐸 𝑅 h-2 ), which is justified in the limit where we stand since h is of order unity, and 𝐸 𝑅 ≪ 1. Physically, this approximation means that the Ekman boundary layers are very thin compared to the fluid height. Substituting the horizontal velocities with their expressions (51) and ( 52) into the vorticity equation ( 46), we obtain in its condensed form the final vorticity equation

𝜕 ζ 𝜕 t +  ( q, ψ) - 𝐸 1∕2 𝑅 2 h ∇ ψ ⋅ ∇ q = 𝐸 𝑅 2 ∇2 𝜁 - 𝐸 1∕2 𝑅 2 h ζ ( ζ + 1) + F ( 53
)
where  is the non-dimensional Jacobian operator in cylindrical coordinates

 (𝑎, 𝑏) = 1 ρ ( 𝜕𝑎 𝜕 ρ 𝜕𝑏 𝜕𝜙 - 𝜕𝑏 𝜕 ρ 𝜕𝑎 𝜕𝜙 ) , ( 54
)
and we introduced the potential vorticity

q = ζ + 1 h . ( 55
)
Taking the curl of ( 51)-( 52), we obtain the modified Poisson equation that links the vorticity and the streamfunction, which closes the system of equations:

ζ = - 1 h ∇2 ψ + 1 h2 ∇h ⋅ ∇ ψ + 𝐸 1∕2 𝑅 h2  ( h, ψ). ( 56
)
Note that the potential vorticity q is a materially conserved quantity if we drop the forcing term and neglect viscous effects. Equation (53) can indeed be recast as

𝜕 q 𝜕 t + 1 h ρ ( 𝜕𝜓 𝜕𝜙 𝜕 q 𝜕 ρ - 𝜕𝜓 𝜕 ρ 𝜕 q 𝜕𝜙 ) = 0 (57) ⇔ 𝐷 q 𝐷 t = 0. (58)

C.2. Justification for keeping higher order, non-linear terms

In the case of QG models derived for rotating convection (e.g. [START_REF] Cardin | Chaotic thermal convection in a rapidly rotating spherical shell: Consequences for flow in the outer core[END_REF][START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF], the Ekman pumping effects are incorporated the same way as we did, except that their geometry is more complicated and the no-slip boundaries cannot be considered as flat. Apart from these geometrical factors, the main difference lies in the fact that in (33), it is common in the rotating convection community to assume that the local vorticity 𝜁 is negligible compared to the planetary vorticity 𝑓 , and retain only the linear terms for the topographic 𝛽-effect and Ekman pumping. In that case, the terms 𝐸 1∕2 𝑅 ∕2 h 𝛁 ψ ⋅ 𝛁 q and ζ 2 𝐸 1∕2 𝑅 ∕2 h are removed from equation (53). These two terms hence represent corrections due to nonlinear Ekman effects. The first term corrects the potential vorticity advection, while the second one is a correction of vortex stretching effects. Following Sansón andVan Heijst (2000, 2002), we argue that these terms should be kept in our simulations since they are at least of same order as the term 𝐸 𝑅 ∕2∇ 2 ζ which represents the bulk viscous effects.

This can be verified by introducing the local Rossby number, based on the local vorticity, 𝑅𝑜 𝜁 = 𝜁 ∕𝑓 = ζ , which is different from the global Rossby number defined in the main text, 𝑅𝑜 = 𝑢 rms ∕𝑓 𝑅. In our case, we are in a regime where both the Ekman and global Rossby numbers are small, meaning that rotation dominates respectively viscous effects and inertia. Both conditions are mandatory to legitimately assume a two-dimensionalization of the flow. However, the local Rossby number may not be small and the local vorticity 𝜁 associated with small turbulent eddies may be of same order as the rotation rate 𝑓 ∕2, which is indeed verified in our simulations. This justifies the fact that we keep supplementary non-linear terms compared to other studies. Using the (non-dimensional) Ekman spin-down time scale 𝐸 -1∕2 as the reference time scale, the different terms of equation ( 53) are indeed of order

𝑎. 𝜕 ζ 𝜕 t ∼ 𝑅𝑜 𝜁 𝐸 1∕2 ∼ 3 × 10 -5 , ( 59
) 𝑏.  ( q, ψ) ∼ (1 + 𝑅𝑜 𝜁 )𝑅𝑜 ∼ 1 × 10 -4 , ( 60
) 𝑐. 𝐸 1∕2 2 h ∇ ψ ⋅ ∇ q ∼ 𝐸 1∕2 2 (1 + 𝑅𝑜 𝜁 )(𝑅𝑜) ∼ 3 × 10 -8 , ( 61
) 𝑑. 𝐸 2 ∇2 ζ ∼ 𝐸 2 𝑅𝑜 𝜁 ∼ 1 × 10 -8 , ( 62
) 𝑒. 𝐸 1∕2 2 h ζ 2 ∼ 𝐸 1∕2 2 𝑅𝑜 2 𝜁 ∼ 3 × 10 -6 , ( 63
) 𝑓 . 𝐸 1∕2 2 h ζ ∼ 𝐸 1∕2 2 𝑅𝑜 𝜁 ∼ 3 × 10 -5 , ( 64
)
where we used 𝐸 ∼ 10 -7 , 𝑅𝑜 𝜁 ∼ 10 -1 , 𝑅𝑜 ∼ 10 -3 . The two non-linear terms discussed (c. and e.) are greater or equal to the bulk viscous effects (d.) and should not be neglected in our case.

C.3. Forcing

For now, we have introduced the forcing as an additional source of vorticity (term F in equation ( 53)). The goal is to reproduce the experimental forcing such that the QG numerical model can be used as a guide and complement the experimental exploration.

In the experiment, because of the Coriolis effect, each inlet (sucking water from the tank) or outlet generates respectively a small cyclone or anticyclone right above it. This process can be modeled as a stationary source of vorticity in the form of positive or negative Gaussian sources of vorticity of radius 𝓁 𝑓 distributed on a prescribed array. We thus define 𝑁 forcing points distributed over the numerical domain, and at each point, we place a Gaussian source of vorticity such that

F ( x, ỹ) = F0 𝑁 ∑ 𝑖=1 (-1) 𝑖 exp ⎛ ⎜ ⎜ ⎝ - [ x -x𝑖 l𝑓 ] 2 - [ ỹ -ỹ𝑖 l𝑓 ] 2 ⎞ ⎟ ⎟ ⎠ , ( 65
)
where ( x, ỹ) are non-dimensional cartesian coordinates, the pairs ( x𝑖 , ỹ𝑖 ), 𝑖 ∈ [ [1, 𝑁]] are the centre of each forcing vortex, 𝓁 𝑓 their radius and 𝐹 0 the forcing amplitude. Various forcing configurations were explored, as listed in Table 4.

C.4. Zonal flow evolution equation

We perform a Reynolds decomposition of the flow by writing the velocity field as an azimuthal (zonal) average plus some fluctuations: )

⟨𝑋⟩ 𝜙 = 1 2𝜋 ∫ 2𝜋 𝜙=0
. ( 70
)
Equation ( 70) shows that the zonal flow is driven by non-linear interactions between eddies, . This term corresponds to the divergence of horizontal Reynolds stresses, sometimes referred to as the eddy momentum flux.

C.5. Numerical methods

The numerical model used in this paper is directly adapted from a previous code used to study convection in a plane layer (Favier, Silvers and Proctor, 2014;[START_REF] Matthews | Compressible magnetoconvection in three dimensions: Planforms and nonlinear behaviour[END_REF].

The non-axisymmetric and axisymmetric motions are solved separately. For the non-axisymmetric motions (azimuthal Fourier modes |𝑚| > 0), we solve the vorticity-stream-function system (53)-( 56) with free-slip boundary conditions at both the outer and inner boundaries. The unknown ψ and ζ are decomposed into their Fourier components up to degree 𝑚 = 2048 and we use a pseudo-spectral method in the azimuthal direction with the 2∕3 rd -rule dealiasing method. In the radial direction, we use centered finite differences of fourth-order on 1024 uniformly-distributed points. The zonal mode (𝑚 = 0) is computed separately by solving equation (70). Time integration is performed with an implicit Crank-Nicolson scheme for the linear operator, applied directly in the Fourier space. An explicit third-order Adams Bashforth scheme is used for the non-linear terms. We use an adaptative time-step so that the CFL stability condition is verified, with a safety factor of 0.1. The code is parallelized using MPI, and we typically ran the simulations over 64 CPUs. For the simulations with 𝐸 𝑅 = 1.25 × 10 -7 , the typical computational time is of ∼20 hours for ∼ 3, 000 rotation times.

D. Discrete Bessel-Fourier transform

The experimental velocity fields being obtained on a discrete grid, the Bessel-Fourier transform coefficients need to be computed using numerical quadrature, and we will thus obtain a finite set of discrete coefficients { û𝑛𝑚 |𝑛 = 1, 2,⃛, 𝑁 𝜌 and 𝑚 = -𝑁 𝜙 ∕2 + 1,⃛, 0, 1, 2,⃛, 𝑁 𝜙 ∕2}. (72)

The associated discrete version of the Parseval relation is

𝑁 𝜙 ∑ 𝑗=1 |𝑓 𝑗 | 2 = 1 𝑁 2 𝜙 𝑁 𝜙 ∑ 𝑚=1 f𝑚 f * 𝑚 . ( 73
)
The kinetic energy per azimuthal wavenumber can then be computed as

∀
) 78
Note that û𝑛𝑚 is in m 3 s -1 , and thus 𝐸 𝑛𝑚 is in m 3 s -2 , which is homogeneous to a kinetic energy per wavenumber. The corresponding wavenumber is 𝑘 𝑛𝑚 = 𝛼 𝑛𝑚 ∕(2𝜋𝑅). The zonal kinetic energy spectrum 𝐸 𝑧 is the spectrum of the axisymmetric mode 𝑚 = 0 only, and the residual energy spectra 𝐸 𝑟 is the sum of the contributions of the remaining modes:

𝐸 𝑧 = 𝐸 𝑛0 , ( 79
)
𝐸 𝑟 = ∑ 𝑚,𝑚≠0 𝐸 𝑛𝑚 . ( 80
)
Note that each azimuthal mode, 𝑚, has a different corresponding wavevector, 𝑘 𝑛𝑚 because the zeros of a Bessel function, 𝛼 𝑛𝑚 , depends on its order 𝑚. To perform the summation leading to the residual spectrum, we hence use data binning. More precisely, the magnitudes of all spectra with 𝑚 ≠ 0 are combined and sorted into wavevector bins, which are defined to correspond to the zonal wavevector 𝑘 𝑛0 . Then, the amplitudes contained in each bin are added together. This process provides a single residual spectrum associated to a unique wavevector. Read (2019). The depth of the fluid layer for Jupiter and Saturn is taken from [START_REF] Galanti | Combined magnetic and gravity measurements probe the deep zonal flows of the gas giants[END_REF]. The kinematic viscosity used to compute the Ekman and Reynolds numbers is assumed to be 𝜈 ∼ 3 × 10 -7 m 2 s -1 for Jupiter and Saturn [START_REF] Gastine | Zonal flow scaling in rapidly-rotating compressible convection[END_REF], 1.5 × 10 -5 m 2 s -1 for the atmosphere, and 1 × 10 -6 m 2 s -1

for the oceans.

Ref.

Figure 1 :

 1 Figure 1: Location of previous zonal jets laboratory experiments in the non-dimensional parameter space (𝑅𝑒, 𝐸, 𝑅 * 𝛽 ) (Reynolds number, Ekman number, and zonostrophy index). For the sake of comparison, we used similar definitions for all the experimental studies (first column of Table 1). The data used to plot this figure are provided in Appendix B. We recall that 𝑅 * 𝛽 ≈ 0.5𝑅 𝛽 and takes into account the spectral prefactor (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 (see section 4). The horizontal dashed lines represent the threshold in 𝑅 * 𝛽 between the different regimes determined by Sukoriansky et al. (2007); Galperin et al. (2010) .

Figure 2 :

 2 Figure 2: Experimental set-up. (a) Schematic of the experimental setup. The forcing pattern is sketched at the bottom: each ring C1-C6 is controlled by an independent pump. (b) Picture of the set-up once in solid body rotation. (c) Picture of the bottom forcing plate designed to achieve a uniform topographic 𝛽-effect.

Figure 3 :

 3 Figure 3: (a,c) Instantaneous velocity fields measured from PIV in the statistically steady state of (a) experiment J and (c) experiment A of Table 3. The vectors magnitude scale is the same for the two plots. (b,d) Time-averaged zonally-averaged zonal flow in (b) experiment J and (d) experiment A of Table3. The shaded area represents the envelope of all the instantaneous profiles. Dashed lines represent the radii of the forcing rings.

Figure 4 :

 4 Figure 4: Comparison between the zonal flow (black line) and the Reynolds stresses or eddy momentum fluxes (colored line) in the experiment J (a) and the simulation QG1a (b). The left 𝑦-axis corresponds to the zonal flow velocity, and the right 𝑦-axis is for the Reynolds stresses. The colour of the line represents the divergence of the Reynolds stresses. When it is positive (red), it indicates a convergence of prograde momentum, and a divergence when it is negative (blue).

Figure 5 :

 5 Figure 5: Instantaneous maps of vorticity 𝜁 (left side of each panel) and azimuthal velocity 𝑢 𝜙 (right side of each panel). (a) Experiment J. (b) References simulation QG1a with a polar forcing. (c,d,e) Cartesian forcing simulations (QG2, QG3 and QG4) with a progressively reduced forcing scale.

Figure 6 :

 6 Figure 6: Maps of azimuthal velocity once in steady state, illustrating the multistability. (a-c) Experiments: three realizations with the same forcing as Exp.J. (d-f) QG simulations: three realizations with the parameters of simulation QG1a. The initial condition of the simulations is a fluid at rest, plus a random noise of non-dimensional amplitude 1 × 10 -3 added to the vorticity field. For each realization, the seed of the initial noise was changed.

Figure 7 :

 7 Figure7: Kinetic energy spectra for two different experiments and three QG numerical simulations. We separate the zonal (𝐸 𝑧 ) and residual (𝐸 𝑟 ) contributions. The dashed-dotted lines correspond to the theoretical predictions given by equations (20) and (21). The shaded area correspond to scales smaller than the forcing scale, 𝑘 > 𝑘 𝑓 . (a) Experiment M of Table3.(b) Experiment B ofTable 3. The thin dark blue lines correspond to 𝐸 𝑧 computed during the transient of the experiment. For 𝐸 theo

Figure 8 :

 8 Figure 8: Comparison between estimates of the upscale energy transfer rates. Circles: experiments. Diamonds: QG simulations. (a) Upscale energy transfer rate measured on the spectra, 𝜖 𝑆 , versus the one estimated assuming pure dissipation by Ekman friction, 𝜖 𝐸 (equation (7)). (b) 𝜖 𝑆 versus the transfer rate deduced from the Thorpe scale, 𝜖 𝑇 , assuming that 𝐿 𝑇 𝛽 = 𝐿 𝑇 ∕0.47 (Fig.11(b)) and 𝜖 𝑇 = (𝐿 𝑇 𝛽 ∕2𝜋) 5 𝛽 3 . Vertical error bars account for the uncertainty in 𝜖 𝑆 due to the uncertainty in the slope measured on the spectra (see Table5). Horizontal error bars account for the standard deviation when measuring the Thorpe scale.

Figure 9 :

 9 Figure 9: Illustration of potential vorticity mixing for experiments J (first row) and A (second row). (a,d) Instantaneous PV maps. (b,e) PV profiles. The dashed black line is the background PV profile, 𝑓 ∕ℎ, the continuous black line is the time and zonally-averaged PV profile ⟨𝑞⟩ 𝜙 (𝜌), and the purple line is the equivalent latitude PV profile 𝑞(𝜌 𝑒 𝑞). The grey lines show the variability of the PV profiles before performing the zonal average. The grey dashed-dotted line is the equivalent staircase 𝑞 𝑠 defined in the text. The inserts show the radial derivative of each PV profile. The stars show the locations of PV jumps selected to compute the equivalent staircase. (c,f) Instantaneous PV profile at a fixed angular position 𝑞(𝜌, 𝜙 0 ). Grey line: before sorting, Blue line: after sorting. The Thorpe scale deduced from the sorting process for that specific angle 𝜙 0 is indicated at the top of each panel.

Figure 10 :

 10 Figure 10: Potential vorticity maps and profiles in simulations QG1a (a,d), QG1b (b,e) and QG1c (c,f). (a-c) Instantaneous PV maps once in quasi-steady state. (d-f) PV profiles. The dashed black line is the background PV profile, 𝑓 ∕ℎ, the thick black line is the time and zonally-averaged potential vorticity profile ⟨𝑞⟩ 𝜙 (𝜌). The grey lines show PV profiles for different angles before performing the zonal and time averages. The grey dashed-dotted line is the equivalent staircase 𝑞 𝑠 defined in the text. The inserts show the radial derivative of the mean PV profile. The stars show the locations of PV jumps selected to compute the equivalent staircase.

Figure 11 :

 11 Figure 11: Comparison between the Thorpe scale 𝐿 𝑇 and (a) the Rhines scale (equation (2)), (b) the transitional scale based on the upscale energy transfer, and (c) the spectral transitional scale (equation (26)). Circles: experiments. Diamonds: QG simulations. The equation at the top of each panel is the best linear fit of the experimental data, represented by the dashed line.

  For simplicity, we denote 𝑓 𝑗𝑖 = 𝑓 (𝜌 𝑖 , 𝜙 𝑗 ) where 𝜌 𝑖 and 𝜙 𝑗 are the discrete radial and azimuthal positions respectively(𝑖 ∈ [[1, 𝑁 𝜌 ]], 𝑗 ∈ [[1, 𝑁 𝜙 ]]).Discrete Fourier transformThe angular part of the transform is performed using the Matlab fft function. For each radius, we thus compute a discrete FFT in the azimuthal direction. The direct and inverse discrete transforms read respectively f𝑚 (𝜌) =

Figure 12 :

 12 Figure 12: Space-time diagrams and zonal flow profiles for a typical experiment in regime II (a,b), and for QG simulations QG1a (c,d), QG2 (e,f), QG3 (g,h) and QG4 (i,j) (see Table4).

  Figure 12: Space-time diagrams and zonal flow profiles for a typical experiment in regime II (a,b), and for QG simulations QG1a (c,d), QG2 (e,f), QG3 (g,h) and QG4 (i,j) (see Table4).

Table 1

 1

Table 2

 2 Symbols and definitions used in the paper

	Symbol Expression 𝑘 𝑅 ( 𝛽∕2𝑢 rms ) 1∕2 𝜖 𝐸 (𝑢 2 rms ∕2)Ω𝐸 1∕2 𝜖 𝑆 fit on residual spectra 𝜖 𝑇 (𝑘 𝑇 𝛽 ) 5 𝛽 3 𝑘 𝐸 𝛽 ( 𝛽 3 ∕𝜖 𝐸 ) 1∕5 𝑘 𝜖 𝛽 ( 𝛽 3 ∕𝜖 𝑆 ) 1∕5 𝑘 𝑆 𝛽 ( 𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 ( 𝛽 3 ∕𝜖 𝑆 ) 1∕5 Transitional wavenumber estimated on the spectral slopes intersection Definition Rhines wavenumber Inverse energy transfer rate estimated a priori assuming Ekman friction only Inverse energy transfer rate measured on the spectra using Eq.(20) Inverse energy transfer rate deduced from the Thorpe scale Transitional wavenumber estimated using 𝜖 𝐸 Transitional wavenumber estimated using 𝜖 𝑆 𝑅 𝐸 𝛽 𝑘 𝐸 𝛽 ∕𝑘 𝑅 Zonostrophy index estimated using 𝑘 𝐸 𝛽 𝑅 𝜖 𝛽 𝑘 𝜖 𝛽 ∕𝑘 𝑅 Zonostrophy index based on 𝑘 𝜖 𝛽 𝑅 𝑆 𝛽 𝑘 𝑆 𝛽 ∕𝑘 𝑅 Zonostrophy index based on 𝑘 𝑆 𝛽

Table 3

 3 Parameters of the experiments and QG simulations used in the present paper. 𝑢 rms is averaged over the last 500 rotation times of each experiment (typically between 𝑡 = 2500 𝑡 𝑅 and 𝑡 = 3000 𝑡 𝑅 ). 𝜖 𝐸 is an estimate of the upscale energy transfer rate based on dissipation in the Ekman boundary layers (equation (7)). The Ekman number, 𝐸, the Reynolds number 𝑅𝑒 and the zonostrophy index 𝑅 𝛽 are defined in Table1. The superscript ⋅ 𝐸 indicates that 𝑅 𝐸 𝛽 is estimated a priori using 𝜖 𝐸 (see Table2).

	Label A B C D E F G H I J K L M N QG1a QG1b QG1c QG2 QG3 QG4	𝛽 (m -1 s -1 ) (cm s -1 ) 𝑢 rms 50.1 3.44 50.1 2.73 65.5 3.23 22.8 3.60 50.1 3.01 50.1 3.36 50.1 3.32 50.1 1.64 50.1 1.85 50.1 1.88 22.8 1.64 50.1 1.19 50.1 1.26 50.1 0.92 50.1 2.21 50.1 3.06 50.1 2.96 50.1 2.68 50.1 2.40 50.1 2.45	𝜖 𝐸 (m 2 s -3 ) 28.7 × 10 -7 18.1 × 10 -7 26.1 × 10 -7 28.0 × 10 -7 21.8 × 10 -7 27.2 × 10 -7 26.6 × 10 -7 6.46 × 10 -7 8.25 × 10 -7 8.51 × 10 -7 5.82 × 10 -7 3.40 × 10 -7 3.84 × 10 -7 2.07 × 10 -7 6.90 × 10 -7 9.00 × 10 -7 4.88 × 10 -7 10.1 × 10 -7 5.67 × 10 -7 4.92 × 10 -7	𝐸 (×10 -7 ) (×10 3 ) (×10 -3 ) 𝑅𝑒 𝑅𝑜 3.78 20.0 7.57 3.78 15.9 6.00 3.55 18.8 6.55 4.73 20.9 9.87 3.78 17.4 6.60 3.78 19.5 7.37 3.78 19.2 7.28 3.78 9.48 3.59 3.78 10.7 4.05 3.78 10.9 4.12 4.73 9.52 4.50 3.78 6.88 2.60 3.78 7.31 2.77 3.78 5.36 2.03 1.27 22.2 2.82 0.60 44.5 2.67 0.20 74.8 1.50 1.27 26.8 3.42 0.63 34.3 2.14 0.42 43.3 1.81	𝑅 𝐸 𝛽 4.98 4.86 5.04 4.71 4.91 4.96 4.96 4.62 4.68 4.69 4.36 4.48 4.50 4.37 5.32 4.93 6.59 5.42 5.76 6.01

  )). The transitional wavenumber based on the spectral slopes intersection is denoted 𝑘 𝑆 𝛽 whereas the transitional wavenumber based on 𝜖 only is denoted 𝑘 𝜖 𝛽 (equation (26)). 𝑘 𝑅 is the Rhines wavenumber (equation (2)), 𝑅 𝑆 𝛽 = 𝑘 𝑆 𝛽 ∕𝑘 𝑅 and 𝑅 𝜖 𝛽 = 𝑘 𝜖 𝛽 ∕𝑘 𝑅 . 𝐿 𝑇 is the Thorpe scale measured experimentally and discussed in section 5.

	Label A B C D E F G H I J K L M N QG1a 25.4 ∈ [20.7, 31.1] 0.25 𝜖 𝑆 𝐶 𝑍 (×10 -7 m 2 s -3 ) 70.8 ∈ [61.3, 81.7] 0.32 30.3 ∈ [26.5, 34.7] 0.27 55.6 ∈ [43.8, 70.5] 0.13 96.4 ∈ [78.9, 117] 0.86 47.0 ∈ [42.2, 52.4] 0.27 66.3 ∈ [58.2, 75.5] 0.27 49.0 ∈ [41.6, 57.7] 0.25 7.69 ∈ [6.87, 8.61] 0.13 8.73 ∈ [5.78, 13.2] 0.21 8.37 ∈ [5.66, 12.4] 0.18 5.09 ∈ [4.24, 6.13] 0.38 4.44 ∈ [4.03, 4.90] 0.11 4.78 ∈ [4.17, 5.47] 0.10 1.70 ∈ [1.12, 5.83] 0.06 QG1b 49.8 ∈ [36.1, 68.8] 0.57 QG1c 20.0 ∈ [16.2, 24.9] 0.37 QG2 84.1 ∈ [69.3, 102] 0.43 QG3 46.2 ∈ [41.8, 51.1] 0.41 QG4 42.3 ∈ [38.0, 47.1] 0.67	𝑘 𝑆 𝛽 (rad m -1 ) (rad m -1 ) (rad m -1 ) 𝑘 𝜖 𝛽 𝑘 𝑅 46.3 112.2 27.2 52.6 132.9 30.5 43.9 138.3 32.1 36.7 65.7 18.2 48.3 121.7 29.1 44.6 113.7 27.6 46.5 120.7 27.6 56.5 174.8 39.3 62.7 170.5 36.9 60.4 171.9 36.6 51.8 118.4 26.7 59.3 195.1 46.1 56.2 192.3 45.0 59.6 236.4 52.2 53.3 137.7 30.6 59.6 120.4 32.3 62.5 144.4 31.80 49.2 108.4 33.6 54.5 122.2 28.6 64.3 124.4 29.07	𝑅 𝑆 𝛽 1.72 4.16 2.61 𝑅 𝜖 𝛽 𝐿 𝑇 (cm) 1.73 4.39 2.11 1.38 4.35 2.11 2.06 3.69 3.51 1.67 4.22 2.23 1.63 4.16 2.72 1.69 4.39 2.60 1.44 4.47 1.99 1.70 4.63 1.57 1.65 4.71 1.80 1.96 4.49 2.56 1.29 4.24 1.38 1.26 4.32 1.35 1.15 4.54 0.98 1.61 3.55 1.89 1.69 3.78 2.53 2.02 3.91 1.99 1.58 4.09 2.38 2.08 4.21 2.94 2.15 4.97 2.74
	𝑎 https://mathworks.com/help/curvefit/fit.html

  𝐿 𝜖 𝛽 = (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 𝐿 𝑆 𝛽 , with (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 ≈ 0.29 -0.41. Hence, we expect 𝑐 ≈ (𝐶 𝑍 ∕𝐶 𝐾 ) 3∕10 𝑎, which is again what is observed. That being said, considering 𝐿 𝜖 𝛽 alone, this result is a demonstration that 𝐿 𝑇 can be used as a proxy of the typical length scale of the flow. If one assumes 𝐿 𝑇 ∼ 0.47𝐿 𝜖 𝛽 as determined from Fig.11(b), we can define a Thorpe estimate of the turbulent dissipation rate 𝜖 𝑇 computed from the relation between 𝐿 𝛽 and 𝜖 (equation (

𝜖 𝛽 ≈ 3.69 -4.71, hence, we should have 𝑎 ≈ 3.69 -4.71 𝑏, which is indeed what is observed. Similarly,

  Note that the zonally-averaged radial velocity, Ũ𝜌 , is not zero because of the Ekman pumping. Instead, from (47) we have Taking the zonal mean of the 𝜙-component of the momentum equation (31) gives an equation for the zonal flow evolution:

				𝑋𝑑𝜙,						(66)
	ũ𝜙 = ⟨ ũ𝜙 ⟩ 𝜙 + ũ′ 𝜙 = Ũ𝜙 ( ρ, t) + ũ′ 𝜙 ( ρ, 𝜙, t) ũ𝜌 = ⟨ ũ𝜌 ⟩ 𝜙 + ũ′ 𝜌 = Ũ𝜌 ( ρ, t) + ũ′ 𝜌 ( ρ, 𝜙, t).						(67) (68)
	Ũ𝜌 =	𝐸 𝑅 1∕2 2 h Ũ𝜙 .							(69)
	𝜕 Ũ𝜙 𝜕 t + Ũ𝜌	𝜕 Ũ𝜙 𝜕 ρ -	Ũ𝜙 Ũ𝜌 ρ + Ũ𝜌 = -⟨ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ũ′ 𝜌 𝜕 ũ′ 𝜙 𝜕 ρ -ũ′ 𝜙 ũ′ 𝜌 ρ ⟩ 𝜙	+	𝐸 𝑅 2	(	∇2	Ũ𝜙 -	Ũ 2 𝜙 ρ2
				

  𝑚, 𝐸 𝑚 (𝜌) = 1 Weyn and Menchaca, 2017, for more details).Discrete Hankel transformFor each azimuthal mode 𝑚, its Fourier transform coefficient f𝑚 (𝜌) has a radial structure onto which we perform a discrete Hankel transform of order 𝑚. We use the Matlab algorithm provided by Guizar-Sicairos and[START_REF] Guizar-Sicairos | Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields[END_REF]. Briefly, for each mode 𝑚, the discrete Hankel transform coefficients are computed following where f𝑚𝑖 = f𝑚 (𝜌 𝑖 ) is the radial structure of each Fourier mode, 𝑅 is the maximum radius, 𝑉 is the maximum radial wavenumber, which is different for each mode 𝑚 (𝑉 𝑚 = 𝛼 𝑁 𝜌 +1,𝑚 ∕(2𝜋𝑅)), and 𝑆 = 2𝜋𝑅𝑉 . Here again, 𝛼 𝑛𝑚 is the 𝑛 th zero of the Bessel function of the first kind of order 𝑚, 𝐽 𝑚 . The corresponding discrete Parseval relation is Taking into account both the Fourier and Hankel transforms, the kinetic energy for each mode (𝑚, 𝑛) can then be expressed as

				𝑁 2 𝜙	[ û𝑚	û * 𝑚 + v𝑚	v * 𝑚	]	(74)
	(see Durran, ∀ 𝑚, f𝑛𝑚 =	1 𝜋𝑉 2	𝑁 𝜌 ∑ 𝑖=1	f𝑚𝑖 𝑚+1 (𝛼 𝑚𝑖 ) 𝐽 2	𝐽 𝑚	( 𝛼 𝑛𝑚 𝛼 𝑚𝑖 𝑆	)	,	(75)
	∀ 𝑚, f𝑚𝑖 =	1 𝜋𝑅 2		𝑁 𝜌 ∑ 𝑛=1		f𝑛𝑚 𝑚+1 (𝛼 𝑚𝑛 ) 𝐽 2	𝐽 𝑚	( 𝛼 𝑛𝑚 𝛼 𝑚𝑖 𝑆	)	(76)
	∀ 𝑚,	𝑁 𝜌 ∑ 𝑖=1	| f𝑚𝑖 | 2 2𝜋 2 𝑉 2 𝐽 2 𝑚+1 (𝛼 𝑚𝑖 )	=	𝑁 𝜌 ∑ 𝑛=1	| f𝑚𝑛 | 2 2𝜋 2 𝑅 2 𝐽 2 𝑚+1 (𝛼 𝑚𝑛 )	.	(77)
	∀ (𝑛, 𝑚), 𝐸 𝑛𝑚 =	1 𝑁 2 𝜙	1 2𝜋 2 𝑅 3 𝐽 2 𝑚+1 (𝛼 𝑛𝑚 )	[	û𝑛𝑚	û * 𝑛𝑚 + v𝑛𝑚	v *

𝑛𝑚

]

. (

Table 6 :

 6 Dimensional and non-dimensional parameters of present and previous experimental studies of zonal jets. 𝑢 rms is the total root-mean-squared velocity, ℎ 0 the mean fluid thickness, 𝑅 the tank radius, Ω the rotation rate, and the 𝛽-effect arises from the paraboloidal free surface, a sloping bottom, or both depending on the study. When a range of parameters was explored, we indicate the values of the most extreme cases. For rotating annuli, 𝑅 is the gap width, not the external radius. The non-dimensional parameters definitions are provided in Table1. For the planetary flows, 𝑢 rms , 𝛽 and 𝑅 𝛽 are taken from Table13.1 in Galperin and

D. Lemasquerier, B. Favier and M. Le Bars: Preprint submitted to Elsevier

Aknowledgements

The authors acknowledge funding by the European Research Council under the European Union's Horizon 2020 research and innovation program through Grant No. 681835-FLUDYCO-ERC-2015-CoG. The authors are most grateful to E. Bertrand and W. Le Coz for their help during the conception and building of the experiment, and to J.-J. Lasserre for his help in setting up the PIV system. We also thank David Dritschel for providing numerical codes to compute equivalent latitude potential vorticity profiles. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources. This work was performed using HPC resources from GENCI-IDRIS (Grants 2021-A0100407543 and 2022-A0120407543).

Data Availability

The data used to plot Figures 1, 8 and 11 are provided in Tables 3, 5 and B.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A. Hovmoller diagrams and zonal flows for experiment and QG simulations

On Figure 12, we plot space-time diagrams and zonal flow profiles for a typical experiment and four QG numerical simulations.