Daphné Lemasquerier

B Favier

M Le Bars

Zonal jets at the laboratory scale: hysteresis and Rossby waves resonance

Keywords: quasi-geostrophic flows, waves in rotating fluids, bifurcation

The dynamics, structure and stability of zonal jets in planetary flows are still poorly understood, especially in terms of coupling with the small-scale turbulent flow. Here, we use an experimental approach to address the questions of zonal jets formation and long-term evolution. A strong and uniform topographic β-effect is obtained inside a water-filled rotating tank thanks to the paraboloidal fluid free upper surface combined with a specifically designed bottom plate. A small-scale turbulent forcing is performed by circulating water through the base of the tank. Time-resolving PIV measurements reveal the self-organization of the flow into multiple zonal jets with strong instantaneous signature. We identify a subcritical bifurcation between two regimes of jets depending on the forcing intensity. In the first regime, the jets are steady, weak in amplitude, and directly forced by the local Reynolds stresses due to our forcing. In the second one, we observe highly energetic and dynamic jets of width larger than the forcing scale. An analytical modeling based on the quasi-geostrophic approximation reveals that this subcritical bifurcation results from the resonance between the directly forced Rossby waves and the background zonal flow.

Introduction

A recurrent feature of planetary fluid envelopes is the presence of east-west flows of alternating direction, so-called zonal jets. Zonal flows are particularly striking in the gas giants' atmospheres such as on Jupiter, where the zonation of clouds of ammonia and water ice reveals the presence of several jet streams [START_REF] Ingersoll | Dynamics of Jupiter's atmosphere. Jupiter: The Planet[END_REF][START_REF] Vasavada | Jovian atmospheric dynamics: An update after Galileo and Cassini[END_REF]. On these gas giants, it has been suggested that at the top of the clouds, the jets may contain more than 90% of the total kinetic energy (Galperin et al. 2014b) and penetrate deep into the planet's interior [START_REF] Kaspi | Jupiter's atmospheric jet streams extend thousands of kilometres deep[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]. Apart from their strength, jets on gas giants are also puzzling by their stability since the pattern has barely varied over decades [START_REF] Porco | Cassini imaging of Jupiter's atmosphere, satellites, and rings[END_REF][START_REF] Tollefson | Changes in Jupiter's Zonal Wind Profile preceding and during the Juno mission[END_REF]. On Earth, at least one zonal jet lies in each hemisphere of the atmosphere [START_REF] Schneider | The general circulation of the atmosphere[END_REF]. Perhaps surprisingly, the observation of zonal flows on the gas giants significantly predates that of the zonal flows in the Earth's oceans. This might be explained by the fact that oceanic jets only appear after a careful time averaging [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF][START_REF] Maximenko | Stationary mesoscale jet-like features in the ocean[END_REF][START_REF] Ivanov | System of quasi-zonal jets off California revealed from satellite altimetry[END_REF]. Despite their latent nature, these jets seem to penetrate deep into the ocean (e.g. Cravatte et al. 2012).

There is not yet a commonly accepted mechanism to explain the formation of zonal † Email address for correspondence: lemasquerier.pro@protonmail.com jets in planetary flows. The only consensus is that the β-effect, arising from the variation of the Coriolis effect with latitude [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF], is responsible for the anisotropisation of the turbulent flow. In his seminal paper, [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF] predicted that the β-effect would alter the inverse energy cascade expected in geostrophic turbulence, and redirect energy towards zonal modes at low wavenumbers. This work was however mainly heuristic, and since then, the dynamical process of jet formation has been the subject of intensive study. In a recent book, Galperin & Read (2019) provide a survey of the latest theoretical, numerical and experimental advancements focusing on zonal jets dynamics and their interactions with turbulence, waves, and vortices. As summarized by [START_REF] Bakas | Emergence of Large Scale Structure in Barotropic β-Plane Turbulence[END_REF], several processes can lead to zonal flows formation, such as anisotropic turbulent cascades [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]Galperin et al. 2019), modulational instability [START_REF] Lorenz | Barotropic Instability of Rossby Wave Motion[END_REF][START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF][START_REF] Manfroi | Slow Evolution of Zonal Jets on the Beta Plane[END_REF][START_REF] Berloff | A mechanism of formation of multiple zonal jets in the oceans[END_REF][START_REF] Connaughton | Modulational instability of Rossby and drift waves and generation of zonal jets[END_REF]) and mixing of potential vorticity [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF][START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF][START_REF] Read | Zonal Jet Flows in the Laboratory: An Introduction[END_REF]. Zonal flows also emerge as statistical equilibria from complex turbulent flows (Galperin & Read 2019, part VI and references therein). It is not clear yet which mechanism(s) is (are) the most relevant for planetary applications, and for which planetary flow (terrestrial ocean and atmosphere, gas giant atmospheres). For instance, and as pointed out by [START_REF] Bakas | Emergence of Large Scale Structure in Barotropic β-Plane Turbulence[END_REF], the inverse energy cascade from a small-scale forcing and its anisotropisation by the β-effect implies spectrally local interactions which are not observed in the Earth atmosphere, or at least not at low latitudes where nonlocal eddy-mean flow interactions are expected to prevail [START_REF] Chemke | The Effect of Eddy-Eddy Interactions on Jet Formation and Macroturbulent Scales[END_REF]. Non-local energy transfers towards the mean flow have also been demonstrated in the heated rotating annulus experiments [START_REF] Wordsworth | Turbulence, waves, and jets in a differentially heated rotating annulus experiment[END_REF]. On Jupiter on the contrary, the large-scale circulation seems indeed to be powered by a well defined inverse cascade emanating from the scale of baroclinic instabilities at ∼2000 km [START_REF] Young | Forward and inverse kinetic energy cascades in Jupiter's turbulent weather layer[END_REF]. Then, as a second example, robust zonal jets can form thanks to eddies even when the mixing is not sufficient to turn the initial potential vorticity profile into a staircase profile [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]. Finally, the relevance of statistical theories [START_REF] Bouchet | Statistical mechanics of two-dimensional and geophysical flows[END_REF][START_REF] Read | Zonal Jet Flows in the Laboratory: An Introduction[END_REF], where both the forcing and the dissipation are vanishing, remains to be addressed for planetary flows.

In the present study, we wish to better understand zonal jets formation thanks to an experimental setup which allows for the self-organization of the flow into a dominant and instantaneous zonal flow made of multiple jets. In the past, numerous studies focused on the characteristics of directly forced zonal flows, either through an imposed zonal acceleration [START_REF] Hide | Detached shear layers in a rotating fluid[END_REF][START_REF] Niino | An Experimental and Theoretical Study of Barotropic Instability[END_REF][START_REF] Nezlin | Experimental Configurations[END_REF][START_REF] Früh | Experiments on a barotropic rotating shear layer. Part 1. Instability and steady vortices[END_REF][START_REF] Aguiar | A laboratory model of Saturn's North Polar Hexagon[END_REF] or a radial one which converts into a zonal acceleration following the action of the Coriolis force [START_REF] Hide | On source-sink flows in a rotating fluid[END_REF][START_REF] Sommeria | Laboratory model of a planetary eastward jet[END_REF][START_REF] Solomon | Shear flow instabilities and Rossby waves in barotropic flow in a rotating annulus[END_REF]. Such a situation is relevant for certain terrestrial circulations such as the oceanic zonal currents forced by the wind, or the subtropical jet driven by the poleward motion in the Hadley cell [START_REF] Read | Zonal Jet Flows in the Laboratory: An Introduction[END_REF]. Here, in the context of self-organized large scale jets, we are interested in the formation of jets through the indirect effect of the Reynolds stresses, as a result of systematic correlations in the small-scale turbulent flow. Reproducing zonal jets without directly forcing them is experimentally challenging, namely because of the large boundary dissipation and small β-effect typically obtained in laboratory setups. Generating significant zonal motions in an indirectly forced and dissipative rotating flow requires:

-a process by which eddying turbulent motions are constantly generated; -a significant β-effect coupled with a small Ekman number E = ν ΩH 2 for the boundary dissipation to be as small as possible (ν being the kinematic viscosity, Ω the rotation rate and H the typical fluid height).

The first point can be achieved thanks to natural instabilities such as barotropic [START_REF] Condie | A convective model for the zonal jets in the atmospheres of Jupiter and Saturn[END_REF][START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF][START_REF] Read | An experimental study of multiple zonal jet formation in rotating, thermally driven convective flows on a topographic beta-plane[END_REF] or baroclinic thermal convection, in the differentially heated rotating annulus configuration [START_REF] Hide | Sloping convection in a rotating fluid[END_REF][START_REF] Bastin | A laboratory study of baroclinic waves and turbulence in an internally heated rotating fluid annulus with sloping endwalls[END_REF][START_REF] Bastin | Experiments on the structure of baroclinic waves and zonal jets in an internally heated, rotating, cylinder of fluid[END_REF][START_REF] Wordsworth | Turbulence, waves, and jets in a differentially heated rotating annulus experiment[END_REF][START_REF] Smith | Multiple Zonal Jets in a Differentially Heated Rotating Annulus* ,+[END_REF]. Note that [START_REF] Read | Jupiter's and Saturn's convectively driven banded jets in the laboratory[END_REF][START_REF] Read | Dynamics of Convectively Driven Banded Jets in the Laboratory[END_REF] alternatively used a specific convective forcing by spraying dense (salt) water at the free surface of a fresh water layer, and in a similar fashion, [START_REF] Afanasyev | On the origin of jets in the ocean[END_REF]; [START_REF] Slavin | Multiple zonal jets on the polar beta plane[END_REF]; [START_REF] Zhang | Beta-plane turbulence: Experiments with altimetry[END_REF]; [START_REF] Matulka | Zonal jets in equilibrating baroclinic instability on the polar beta-plane: Experiments with altimetry[END_REF] performed a localized forcing involving sources of buoyancy. Another method, which has the advantage of allowing a close control of the location, scale and intensity of the forcing, consists in applying a mechanical forcing, provided that, again, it does not directly force the mean flow. In that purpose, [START_REF] Whitehead | Mean flow generated by circulation on a β-plane: An analogy with the moving flame experiment[END_REF] used a vertically oscillating disk, [START_REF] Afanasyev | Quasi-two-dimensional turbulence on the polar betaplane: Laboratory experiments[END_REF]; [START_REF] Espa | Zonal jets and cyclone-anticyclone asymmetry in decaying rotating turbulence: Laboratory experiments and numerical simulations[END_REF]; [START_REF] Zhang | Beta-plane turbulence: Experiments with altimetry[END_REF]; Galperin et al. (2014a) employed an eloctromagnetic forcing, and several studies perfomed an eddy-forcing using sinks and sources of fluid [START_REF] Verdiere | Mean flow generation by topographic Rossby waves[END_REF][START_REF] Aubert | Observations of zonal flow created by potential vorticity mixing in a rotating fluid[END_REF][START_REF] Nitto | Simulating zonation in geophysical flows by laboratory experiments[END_REF][START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF][START_REF] Burin | Turbulence and jet-driven zonal flows: Secondary circulation in rotating fluids due to asymmetric forcing[END_REF].

Regarding the second point, since spatial modulation of the Coriolis parameter is difficult to set-up experimentally, the β-effect is usually achieved topographically, i.e. through the variation of the fluid height. Two principal approaches have been tested: using a sloping bottom, associated or not with a top-lid, or using the natural paraboloidal shape adopted by any fluid with a free surface in solid-body rotation. But as explained below, a β-effect alone is not sufficient for the development of large scale zonal flows and should be accompanied by the smallest possible friction. In the context of eddy-driven jets in forced-dissipative experiments, the zonostrophy index R β has been introduced to distinguish friction-dominated regimes (small R β , Earth's ocean and atmosphere) and zonostrophic regimes, i.e. regimes of strong jets (R β > 2.5, Jupiter and Saturn) [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]Galperin et al. 2019, table 13.1). This index is basically the ratio between the largest scale of the dynamics set by the large scale drag, and the scale at which the eddies start to feel the β-effect. To favour the emergence of strong jets, experimentalists should seek large R β , i.e. strong β-effect, strong flows (but still dominated by rotation), and small viscous dissipation thanks to fast rotation and/or large containers [START_REF] Read | Zonal Jet Flows in the Laboratory: An Introduction[END_REF]. Previous experimental studies [START_REF] Nitto | Simulating zonation in geophysical flows by laboratory experiments[END_REF][START_REF] Smith | Multiple Zonal Jets in a Differentially Heated Rotating Annulus* ,+[END_REF][START_REF] Zhang | Beta-plane turbulence: Experiments with altimetry[END_REF][START_REF] Read | An experimental study of multiple zonal jet formation in rotating, thermally driven convective flows on a topographic beta-plane[END_REF] lied in the range R β ∈ [0.73, 1.46] and the observed flows were not in the zonostrophic regime, but recently, [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] were able to reach R β ≈ 3.7 thanks to the fast rotation (75 RPM) of a 1m-diameter tank, thus getting closer to the regime observed on gas giants. The present study follows the work of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF]: we built a very close but significantly improved experimental setup designed to make more quantitative measurements as well as to study more precisely the dependence of the obtained jets on the forcing amplitude. We also wish to focus on their long-term evolution, if any.

The layout of the paper is as follows. In §2, we present the experimental setup. In §3, we describe the main experimental results: we observe a subcritical bifurcation between two regimes of zonal flows depending on the forcing intensity. In §4, we develop a theoretical bidimensional model based on the quasi-geostrophic approximation to explain the experimental results. In §5, we point towards experimental improvements and future work and discuss the implications at the planetary scale . A cylindrical tank of 1m in diameter and 1.6m in height, filled with 600 liters of water, is fixed on a table rotating at 75 RPM. The fluid free-surface takes a paraboloidal shape due to the centrifugally-induced pressure. The bottom plate is designed to achieve a uniform topographic β-effect. A small-scale turbulent forcing is performed by circulating water through 128 holes at the base of the tank. The forcing pattern is sketched on the right: each ring C1-C6 is controlled by an independent pump. Time-resolving PIV measurements are performed on a horizontal plane using a side green laser and a top-view camera.

Experimental methods

The experimental setup is an improved version of the setup of [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF]. Three main modifications were made compared to the initial setup. First, the vast majority of theories and numerical simulations is performed in the context of the socalled β-plane, where the Coriolis parameter is assumed to vary linearly in one direction, with β its (constant) derivative. For that reason, we designed the present setup to have a uniform topographic β effect over the whole tank, rather than a strongly varying one due to the paraboloidal free-surface. Second, in the present experiment, we are able to control the forcing amplitude with radius and we decreased the forcing scale by a factor two. Finally, and most importantly, the tank is transparent which allows for time-resolving particle image velocimetry (PIV) measurements.

Following [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF], our experiment consists of a rapidly rotating cylindrical tank filled with water with a free upper surface and a topographic β-effect induced by the parabolic increase of the fluid height with radius due to the centrifugally induced pressure. The tank, made of Plexiglas, has an external diameter of 1m, is 1 cm thick and 1.6 m high. It is covered with a top-lid also made of Plexiglas to bring the underlying air in solid body rotation, thus reducing as much as possible perturbations of the free-surface. The experimental setup is sketched in figure 1. The topographic β-effect is a source or sink of vorticity following radial motions, and is a consequence of the local conservation of angular momentum in a rapidly rotating fluid. Here, the β parameter can be written as

β = - f h dh dρ , (2.1)
where ρ is the cylindrical radius, h(ρ) is the total fluid height and f = 2Ω is the Coriolis parameter with Ω the rotation rate. Appendix A provides details about the origin of this expression. Equation (2.1) shows that for the topographic β-effect to be uniform over the whole domain, the fluid height should vary exponentially with radius. To achieve this, we choose to compensate the unalterable paraboloidal shape of the free surface using a non-flat bottom plate placed inside of the tank (figure 1 and 2). The total fluid height h above the bottom plate is the difference between the free-surface altitude h p and the bottom topography h b . In solid body rotation at a rate Ω, the water free-surface height as a function of the cylindrical radius ρ is

h p (ρ) = h min + Ω 2 2g ρ 2 = h 0 + Ω 2 2g ρ 2 - R 2 2 , (2.2)
where g is the gravitational acceleration, R is the tank radius, h min the minimum fluid height in rotation and h 0 the fluid height at rest. We want the fluid height h to have an exponential increase with ρ such that the β parameter (equation 2.1) is constant, that is

h(ρ) = h min exp - β 2Ω ρ . (2.
3)

The topography of the bottom of the tank is thus designed such that h b = h p -h. In addition, we optimized the choice of the physical parameters for h b to be the less steep possible in order to minimize the cost of production. Two additional constraints are given by the maximum rotation rate of the turntable (90 RPM) and its maximum load (1500 kg). This process led us to choose h min = 0.20 m;

(2.4) h 0 = 0.58 m;

(2.5)

Ω = 75 RPM ≈ 7.85 rad • s -1 ; (2.6) β ≈ -50.1 m -1 • s -1 . (2.7)
With these parameters, the bottom plate has the shape of a curly bracket (figure 2(c)) with a maximum height difference of 5.36 cm and a mean absolute slope of 22%. The effective fluid height is minimum at the centre, h min = 0.2 m, and increases up to h max = 0.96 m (figure 2(b)). The total volume of water, including the water located below the bottom plate is of about 600 litres. Finally, the chosen rotation rate leads to an Ekman number

E = ν Ωh 2 0 ≈ 3.78 × 10 -7 , (2.8)
where ν is the kinematic viscosity of water (ν = 10 -6 m 2 • s -1 ). Note that we assume ν to be constant, but the experiments discussed in the following were performed at room temperatures varying between 20 and 27 • C, leading to ν ∈ [0.8539, 1.0034]×10 -6 m 2 • s -1 and E ∈ [3.23, 3.80] × 10 -7 . We force small-scale fluid motions using an hydraulic system located at the base of the tank (figures 2(c,d )). This system is inspired from previous setups designed to study turbulent [START_REF] Bellani | Homogeneity and isotropy in a laboratory turbulent flow[END_REF][START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF] and zonal flows [START_REF] Verdiere | Mean flow generation by topographic Rossby waves[END_REF][START_REF] Aubert | Observations of zonal flow created by potential vorticity mixing in a rotating fluid[END_REF][START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF][START_REF] Burin | Turbulence and jet-driven zonal flows: Secondary circulation in rotating fluids due to asymmetric forcing[END_REF]. The curved bottom plate is drilled with 128 holes (64 inlets and 64 outlets) of 4 mm diameter. The forcing pattern is arranged on a polar lattice with 6 rings C 1-6 located at radii R i ∈ {0.067, 0.140, 0.214, 0.287, 0.361, 0.434} m as represented in figure 1. Each ring counts respectively 6, 12, 18, 24, 30 and 38 holes, half of them being inlets (sucking water from the tank and generating cyclones) and the other half outlets (generating anticyclones) as represented in figure 1. The holes are uniformly distributed along each ring, leading to a minimum separating distance of 7.0 cm (ring C 1 ) and a maximum separating distance of 7.6 cm (ring C 5 ). Note that there is also a spatial phase shift between each consecutive rings in order to minimize the variance in the distance between two neighbouring inlets or outlets (figure 1). All the holes of a given ring are connected to a submersible pump (TCS Micropump, M510S-V) via a network of flexible tubes (figure 2). Six submersible pumps are thus located beneath the bottom plate, and circulate water through the six rings. The resulting circulation induces no net mass flux, since the water is directly sucked from the working fluid and released in it. At this point, it is important to stress that the system was designed to minimize the direct forcing of the zonal mean zonal flow and that only the eddy momentum fluxes should be responsible for its eastward or westward acceleration. Finally, each ring is controlled by one pump independently of the others which allows us to control the forcing intensity with the radius. The pumps are controlled remotely by linking them to their drivers (TCS EQi Controllers) through the base of the tank. The drivers are controlled by a Raspberry Pi connected to a local network. We can chose the power of a given pump to be stationary, or to fluctuate randomly within a prescribed power range every 3 seconds.

To measure velocity fields, time-resolving particle image velocimetry (PIV) measurements are performed on a horizontal plane. A green laser beam (Laser Quantuum 532nm CW Laser 2 Watts) associated with a Powell lens is used to create a horizontal laser plane located 11 cm above the edge of the bottom plate (9 cm below the center of the paraboloid). The water is seeded with fluorescent red polyethylene particles of density 0.995 and 40-47 micrometers in diameter (Cospheric, UVPMS-BR-0.995). Their motion is tracked using a top-view camera (Dantec HiSense Zyla) placed above the tank (figure 1 and2). A 28 mm lens is mounted on the camera (Zeiss Distagon T* 2/28). 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. We reduced the sensor region of interest to fit the tank borders, leading to 1900 × 1900 pixels images. Optical distortion induced by the paraboloidal free-surface is corrected on DynamicStudio using a preliminary calibration performed by imaging a plate with a precise dot pattern. An illustrative movie of the particles motion during an experiment is available as supplementary movie 1. 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 figure 3). As represented in figure 1, all the devices (acquisition computer, camera, synchronizer, laser, pumps power supply, drivers and Raspberry) are attached to the rotating frame. The rotary table operates thanks to an air cushion, allowing us to reach high rotation rates even with a large non-equilibrated load (∼ 1000 kg).

A typical experimental run is as follows. We gradually increase the rotary table rotation rate from rest up to 75 RPM (∼30 min). We then wait for the water to be in solid-body rotation which takes approximately 45 min, i.e. ∼ 13 τ E where τ E = Ω -1 E -1/2 is the Ekman spin-up timescale. Note that inertial oscillations are observed even after spin-up, due to the tank's misalignment and the slight non-sphericity of its cross-section. These oscillations generate typical radial root-mean-square velocities of ∼ 4×10 -4 m•s -1 . These small amplitude and large-scale oscillations do not significantly perturb the small-scale forced geostrophic motions. We then turn on the forcing of the 6 rings simultaneously, potentially with different powers, in a stationary or random state. For a typical run, we record images for 60 minutes, corresponding to 4500 t R where t R = 2π/Ω = 0.8 s is the rotation period. We record images with framerates between 10 and 30 frames per second to resolve the fluid motions which are typically between 0.1 and 10 cm/s depending on the forcing amplitude.

The forcing was calibrated in the rotating system by measuring the root-mean square (RMS) velocity induced on the horizontal PIV plane by the forcing. This measurement was realized for each ring separately and several pump powers just after the forcing was turned on, i.e. before the jets develop. We then performed a linear fit of the induced RMS velocity as a function of power to obtain a calibration law for each pump. Details about the forcing calibration are given in appendix B. In the following, we denote the forcing amplitude U f , which corresponds to the mean of the RMS velocities of the six pumps deduced from our calibration.

Experimental results

For all the experiments performed in our setup, we observe instantaneous zonal flows independently of the number of forcing rings turned on, their power, and their state (stationary or random). However, depending on the forcing amplitude, we observe two different regimes of zonal flows described in the next sections. The results are presented on the horizontal laser plane using the polar coordinates represented in figure 1, with (u ρ , u φ ) the radial and azimuthal velocities and ζ = (∇ × u) • e z = (∂ ρ (ρu φ ) -∂ φ u ρ )/ρ the vertical component of the vorticity.

Regime I: Low-amplitude, locally forced jets

At low, stationary forcing amplitude, we observe the fast development of 5 prograde jets -in the same direction as the tank's rotation, u φ > 0 (figure 1) -and 6 retrograde ones. We will refer to this regime as Regime I.

To describe this regime, we chose a typical experiment where the pumps power are respectively P i = {7, 10, 20, 30, 45, 90}% of their nominal power, corresponding to a forcing amplitude U f = 2.4 × 10 -3 m • s -1 (see appendix B). Figure 4(a) represents the temporal evolution (Hovmller diagram) of the instantaneous azimuthal mean of the azimuthal component of the velocity u φ φ (ρ, t) -called zonal flow in the following of the paper, whereas mean flow refers to the time-averaged velocity field. The jets develop almost instantaneously (∼ 10 t R ), and reach their saturating amplitude in about 100 t R (another example is shown in figure 6 for t < 360 t R ). Supplementary movie 2 illustrates the development of regime I.

The velocity and vorticity fields obtained after saturation are represented in figures 3(a) and 5(a-c). The retrograde jets are uniform and quasi-axisymmetric, whereas the progade jets are associated with clear non-axisymmetric perturbations. Consistently with the direction of the zonal flow, the anticyclones -negative vorticity ζ -are located on the outer radius flank of the prograde jets, whereas cyclones are located on their inner radius side. In addition, the prograde jets are thinner than the retrograde ones. These Regime observations highlight the asymmetry between prograde and retrograde jets, generically observed in this type of systems (e.g. [START_REF] Scott | The structure of zonal jets in shallow water turbulence on the sphere[END_REF][START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF].

u RMS (mm/s) Ro = u RMS f R Ro ℓ = u RMS f ℓ Ro ζ = ζmax f Re = u RMS R ν Re ℓ = u RMS ℓ ν I 1.59 1.78×10 -4 1.20 × 10 -3 1.91×10 -
The saturated zonal flow profile is plotted in figure 4(b) along with its time-average, and figure 4(c) shows the zonal mean of the potential vorticity (ζ + f )/h φ . In the absence of dissipation, we expect the material conservation of potential vorticity (PV) [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF], section 4.5). In this limit, zonal flows formation can be viewed as a process of mixing of the initial potential vorticity profile f /h(ρ). [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF] showed that this profile should be turned into a staircase where the prograde jets correspond to steep gradients, and the retrograde jets correspond to weak gradients, i.e. zones of strong mixing. Here, despite the visible segregation of vorticity (figure 5(c)), the initial vorticity profile is almost not perturbed, showing that zonal jets can exist instantaneously even without this process of potential vorticity mixing. Said differently, this regime is characterized by a local Rossby number Ro ζ = ζ/f ≪ 1 (see table 1), hence the initial PV profile is not expected to be strongly modified. The instantaneous root mean square (RMS) velocity defined as

u RMS = 1 N N i=1 |u i | 2 1/2 , (3.1)
where N is the number of PIV velocity vectors, is provided in table 1 along with the global and local Rossby and Reynolds numbers of the flow. Here, u RMS is computed from the velocity fields of figure 3. It is considered "instantaneous" in opposition to the same quantity computed after a very long time average u RMS which will be used later in the paper. Table 1 shows that the flow is barely turbulent in regime I (the local Reynolds number being approximately 100), and highly constrained by rotation given the very small Rossby numbers. Finally, the zonal flow contains 23 ± 5% of the total kinetic energy.

In this regime, each prograde jet stands right above a forcing ring (see the dashed lines in figures 4 and 5). The only exception is the inner ring (ring C1), which is geometrically constrained due to its small radius and significantly perturbed by the peak at the center of the bottom plate (figure 2(c)). Despite this anomaly, the 5 other forcing rings are clearly associated with a prograde jet. This leads us to hypothesize that the prograde jets are forced locally by prograde momentum convergence towards the forcing radii. The local Reynolds stresses generated by our forcing are then balanced by viscous effects. This mechanism of zonal flow formation is reminiscent of the pioneering experiments of [START_REF] Whitehead | Mean flow generated by circulation on a β-plane: An analogy with the moving flame experiment[END_REF] and De [START_REF] Verdiere | Mean flow generation by topographic Rossby waves[END_REF]. [START_REF] Whitehead | Mean flow generated by circulation on a β-plane: An analogy with the moving flame experiment[END_REF] demonstrated that the generation of a train of Rossby waves in a rotating tank with paraboloidal free surface induces a prograde flow at the radius of the forcing, with two weak retrograde flows on both sides of the forcing. De [START_REF] Verdiere | Mean flow generation by topographic Rossby waves[END_REF] did the same observation with a forcing consisting in a ring of sink and sources able to be azimuthally translated. Corresponding theoretical studies of [START_REF] Thompson | A Prograde Jet Driven by Rossby Waves[END_REF] and [START_REF] Mcewan | Mean flows driven by weak eddies in rotating systems[END_REF] then followed and accounted for the mechanism of momentum convergence due to eddy-forcing. It is now believed to be the primary mechanism of westerlies formation in the midlatitude atmosphere (Vallis 2006, chapter 12). This mechanism will be further explored in §4.

Finally, let us mention that the relaxation dynamics of this regime is consistent with the observation of De [START_REF] Verdiere | Mean flow generation by topographic Rossby waves[END_REF]: when the forcing is stopped, the fluctuating velocities are dissipated more rapidly than the mean flow.

Regime II: High-amplitude large jets

At high, stationary forcing amplitude, regime I develops as a transient before the system reaches a statistically stationary state with stronger and broader zonal jets, hereafter called Regime II.

To describe this regime, we chose a typical experiment where the pumps power are repectively P i = {26, 33, 60, 80, 100, 100}% of their nominal power, corresponding to a forcing amplitude

U f = 4.0×10 -3 m • s -1 (see appendix B). The Hovmller diagram of this experiment is represented in figure 4(d ).
The steady jets of the first regime reorganize into 3 prograde and 3 retrograde jets. Note that in other experiments, the saturated flow can count 4 prograde jets instead of 3, as will be discussed in §3.3. This spontaneous transition from regime I to regime II is slow, and the statistically steady state is obtained after a transient of about 800 t R . Furthermore, it involves zonal flows merging events visible in the Hovmller plots of figures 4(d ) and 6. The reorganization of the jets during this transition also shows that they become more independent of the forcing pattern: in the final steady state of regime II, the jets have a typical width which is twice that of the jets in regime I, and their radial position can be shifted compared to the position of the forcing rings. A retrograde flow is even observed above some forcing rings, for instance above C1 and C3. Thus, in this regime, the system self-organizes at a global scale, and the idea of a direct local forcing is not relevant anymore. Supplementary movie 3 illustrates the development of regime II, and movie 1 shows the particles motion when the system is in steady state.

The velocity field obtained after saturation is represented in figure 3(b), and the corresponding maps of velocity and vorticity are plotted in figure 5(d-f ). The prograde jets are still meandering between cyclones on their right and anticyclones on their left, but these vortices are now large-scale ones. As can be seen in figure 3(a), the vortices forced above the inlets and outlets have a typical diameter of ∼3 cm in regime I, whereas in regime II (panel (b)), we observe fewer vortices, with a typical diameter of ∼8 cm. The instantaneous RMS velocity (table 1) is about 10 times higher than in the experiment described for regime I. The global and local Rossby numbers are still very small, i.e. the flow is still highly constrained by rotation, but the Reynolds number is multiplied by 10 hence the flow can now be considered fully turbulent. The fraction of kinetic energy contained in the zonal flow in this experiment reaches 58±8%. Figure 4(f ) shows that the PV mixing is increased in this second regime and consistently with [START_REF] Dritschel | Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers[END_REF], the prograde jets correspond to steepening of the PV profile. But again, the small vorticity of our experiment does not allow an efficient mixing process, though the jets are strong and contains most of the kinetic energy.

Nature of the transition: a first-order subcritical bifurcation

In this section, we investigate the nature of the transition between the two previously described experimental regimes.

Figure 6 shows a Hovmller diagram representing the evolution of the zonal flow profile u φ φ (ρ, t) during a single experiment as well as the corresponding evolution of the total, zonal, and fluctuating kinetic energy defined respectively as

K = 1 N N i=1 |u i | 2 , (3.2) K z = 1 N N i=1 u φ 2 φ,i , (3.3)
K f = K -K z , (3.4)
where N is the number of PIV velocity vectors. The experiment plotted in figure 6 is initialized with a stationary forcing (U f = 3.3 × 10 -3 m • s -1 ), leading to a steady state in regime I. After 360 t R , this forcing is perturbed at finite amplitude by turning the third ring into a random state. Here, it consists in changing its power every 3 seconds to random values uniformly distributed in a range centered around ±20% of its initial power. After such a perturbation, figure 6(a) shows that the system bifurcates towards the second regime through merging events and increasing zonal flow amplitude. Note that without this perturbation, the system would be locked in regime I, as shown by a separate experiment performed with the exact same forcing, at least up to t = 1875 t R .

During the transition, the fraction of kinetic energy contained in the zonal flow increases from 21±7% to 48±9% (figure 6(b)). This second value is significantly lower than the one mentioned previously for regime II since the forcing of this experiment is weaker.

After the transition, the system remains attached to this new steady state even when the forcing is set back to its initial stationary state at time t = 1600 t R . These observations demonstrate that two stable states coexist for this particular forcing and suggest that the transition between the two regimes is of subcritical nature.

To further investigate this bi-stability, we perform series of experiments where we increase or decrease the forcing step by step. We wait significantly between each step for the system to relax towards a new steady state (typically 20 minutes, i.e. 1500 t R ). We then measure the corresponding mean flow amplitude defined as the RMS velocity computed on a time-averaged velocity field:

u RMS = 1 N N i=1 |u i | 2 1/2 , (3.5)
where • denotes the time average over the whole duration of the record once the flow has reached the statistically steady state. Typically, the time average is performed over 200 to 1000 t R depending on the duration of the record. For low values of the forcing amplitude, regime I is observed with the 6 prograde jets structure and u RMS ∼ 2.5 × 10 -3 m • s -1 . As the forcing amplitude increases, the jets structure does not change but their amplitude increases smoothly. When the forcing further increases, a sharp transition occurs around U f ≈ 3.32 × 10 -3 m • s -1 corresponding to a bifurcation from regime I to regime II: both the jets size and amplitude increase abruptly (u RMS ∼ 7.5×10 -3 m • s -1 ). Once in regime II, the amplitude of the jets continues to increase with the forcing amplitude. When the forcing amplitude is gradually decreased, the bifurcation from regime II to regime I is again abrupt, but obtained at a lower forcing U f ≈ 3.11×10 -3 m • s -1 . These hysteresis experiments confirm that the two regimes coexist in a given forcing range U f ∈ [3.11, 3.32] × 10 -3 m • s -1 . The particular forcing of the experiment represented in figure 6 the bistable range in which the first regime is metastable. In §4, we propose a model to explain this hysteresis phenomenon. Finally, we note a significant variability in the mean flow amplitude in regime II. The grey points in figure 7(a), located at U f = 4 × 10 -3 m • s -1 , correspond to nine experiments where we apply the exact same forcing (P i = {26, 33, 60, 80, 100, 100}%), starting from solid-body rotation. Despite the similarity of the forcing, the flow may evolve towards different statistically steady states where the mean flow amplitude and scale are roughly the same, but the position of the jets differ. Based on 15 realizations, we have identified 3 different steady states with permutations between the location of the prograde jets, as represented in figure 7(b). The last point of the yellow curve in figure 7(a) is in configuration 2. It is located below the others probably because it had not reached its steady state when the measurements were performed (500 t R after the forcing change in contrast to 1500 t R for the grey points). Note that we have never observed spontaneous transitions between these three states, even during day-long experiments (38 000 t R ). The origin of this multi-stability is beyond the scope of the present study, but will be the focus of future investigations. In particular, its link with the multi-stability observed recently in numerical simulations of stochastically forced barotropic turbulence (Bouchet et al. 2019a,b) should be addressed.

(U f = 3.3 × 10 -3 m • s -1 ) belongs to
It is of interest to compute the transition rates between the two regimes. On figure 8, we plot the evolution of the total kinetic energy for transitions from regime I to regime II and vice versa in order to compute the corresponding timescales. Transitions in the direction II→I (decreasing power) are accompanied by an exponential decay of the total kinetic energy where K 0 is the kinetic energy in the initial steady state, and K ∞ the kinetic energy reached in the final steady state after the transition. We plot in figure 8(b) the time evolution of the normalized kinetic energy for three transitions with different initial and/or final steady states. Despite these differences, it is clear that the three transitions have the same characteristic time τ II→I ≈ 150 t R . The picture is different for transitions in the direction I→II (figure 8(c)). The kinetic energy increases in a non-trivial way before saturating in a new steady state. We compare this evolution for three transitions starting from the same initial state in the first regime (red dot I in figure 8(a)), but evolving towards three different steady states in regime II (II 1,2,3 in figure 8(a)). Note that I→II 1 corresponds to a finite amplitude perturbation of a steady state in regime I inside of the bistable range. We plot in figure 8(c) the kinetic energy normalized the same way as for II→I transitions. This time, the curves do not collapse. We observe that the closer (in terms of forcing amplitude) the second steady state, the longer the transition. The transition I→II 3 is for instance about 3 times faster than the transition I→II 1 . These differences highlight an asymmetry in the transition mechanism depending on its direction. We propose the following interpretation. In the case II→I, the transition resembles a classical relaxation following a linear dissipation process. If the Reynolds stresses sustaining the strong jets abruptly decrease when the forcing decreases, the linear friction dominates the zonal flow evolution equation ∂ t u φ φ = -α u φ φ (see equations (4.1) and (A 10)), where α is a linear friction, and the zonal flow decreases exponentially. We then expect the timescale of this transition to be of the order of the Ekman friction timescale τ E = α -1 = Ω -1 E -1/2 . In our case, τ E ≈ 206 rotation periods, which is consistent with τ II→I ≈ 150 t R determined previously. On the contrary, in the case I→II, we expect a nonlinear mechanism leading to a non-trivial increase of the zonal flow amplitude. Contrary to the linear friction, this mechanism depends on the forcing amplitude -as may be intuited from the theoretical model developed in §4. The higher the forcing, the faster we expect the transition to occur.

K = K 0 + (K ∞ -K 0 )e -t/

Theoretical model for the transition: a Rossby-wave resonance

The goal of this section is to derive a simple model to explain the transition observed in the experiment, and the associated bi-stability. To do so, we use the classical quasigeostrophic (QG) approximation to reduce the experiment to its 2D β-plane analog (Vallis 2006, p. 67). Because of the fast background rotation, or equivalently the small Rossby number of the system, the geostrophic balance dominates the experimental flow. As a consequence, 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. Nevertheless, the weakness of these effects allows their incorporation into quasi-twodimensional physical models. We derive the conventional QG model corresponding to our experimental setup in appendix A for completeness, "conventional" meaning that we retain only the linear contributions from these 3D effects. Note that in addition to this QG approximation, we make the rigid lid approximation and neglect the temporal fluctuations of the free surface, i.e. we do not take into account gravitational effects at the interface.

We use the cylindrical coordinates (ρ,φ,z) with z oriented downward, and (e ρ , e φ , e z ) the corresponding unit vectors (figure 1). We consider the flow of an incompressible fluid of constant kinematic viscosity ν and density ρ f , rotating around the vertical axis at a constant rate Ω = Ω e z , with Ω > 0 (the turntable rotates in the clockwise direction). The fluid is enclosed inside a cylinder of radius R, and the total fluid height is h(ρ). We denote the velocity field as u = (u ρ , u φ , u z ) eρ,e φ ,ez , and the vertical component of the

vorticity is ζ = (∇ × u) • e z = (∂ ρ (ρu φ ) -∂ φ u ρ )/ρ
β = - f h dh dρ , (4.2) α = E 1/2 f 2 , (4.3)
(see appendix A). In this 2D framework, we decompose the velocity into a zonal mean flow plus fluctuations using the standard Reynolds decomposition:

u φ = u φ φ (ρ, t) + u ′ φ (ρ, φ, t) = U + u ′ φ , (4.4) u ρ = u ρ φ (ρ, t) + u ′ ρ (ρ, φ, t) = u ′ ρ , (4.5) ζ = ζ φ (ρ, t) + ζ ′ (ρ, φ, t) = 1 ρ ∂(ρU ) ∂ρ + ζ ′ . (4.6)
where

• φ = 1 2π 2π 0
• dφ is the zonal mean. Here, we neglect the O E 1/2 mean radial velocity associated with the Ekman pumping, consistently with the choice of keeping only the linear Ekman friction term (see for example the discussion in [START_REF] Sansón | Nonlinear Ekman effects in rotating barotropic flows[END_REF]. The zonal average of the azimuthal component of the Navier-Stokes equation (equation (A 2)) leads to the zonal mean zonal flow evolution equation

∂U ∂t = -u ′ ρ ∂u ′ φ ∂ρ + u ′ ρ u ′ θ ρ φ R(ρ,t) + D(ρ, t), (4.7)
where D contains both the frictional and bulk dissipation of the zonal flow. Using the zero-divergence of the horizontal velocity, the source term R can be expressed as the divergence of the Reynolds stresses, or equivalently as an average vorticity flux

R(ρ, t) = - 1 ρ 2 ∂ ρ 2 u ′ ρ u ′ φ φ ∂ρ = -u ′ ρ ζ ′ φ . (4.8)
Hence, the zonal flow equation ∂U ∂t = R(ρ, t; U ) + D(ρ, t; U ), (4.9)

shows that in the absence of direct forcing, the zonal flow requires a source term which is provided through the Reynolds stresses divergence, alternatively called the eddy momentum flux. To explain the generation of the zonal flow in our experiment, this momentum flux R needs to be modeled. The Reynolds stresses are likely to be influenced by the zonal flow U that they generate through a feedback mechanism. Determining whether the feedback of the meanflow on the source term R is positive or negative would allow us to investigate the possibility of bi-stability. This is the goal of the present section. We follow the same approach as in [START_REF] Herbert | Atmospheric Bistability and Abrupt Transitions to Superrotation: Wave-Jet Resonance and Hadley Cell Feedbacks[END_REF] which focuses on transition to super-rotation based on the mechanism described by [START_REF] Charney | Multiple Flow Equilibria in the Atmosphere and Blocking[END_REF] in the framework of topographically forced zonal flows in the midlatitude atmosphere (see also [START_REF] Pedlosky | Resonant Topographic Waves in Barotropic and Baroclinic Flows[END_REF][START_REF] Held | Stationary and quasi-stationary eddies in the extratropical troposphere: Theory[END_REF][START_REF] Weeks | Transitions Between Blocked and Zonal Flows in a Rotating Annulus with Topography[END_REF][START_REF] Tian | Experimental and numerical studies of an eastward jet over topography[END_REF]. We determine the Reynolds stresses divergence R by computing the linear response to a stationary forcing on a β-plane with a background zonal flow. We show that the resulting Reynolds stresses exhibit a resonant amplification leading to a possible bi-stability. We finally compare this mechanism with the experimental observations.

Linear model for the Reynolds stresses

In this section, we determine the Reynolds stresses divergence R and its sensitivity to the zonal flow. To do so, we compute the linear response to a stationary forcing on the β-plane, in the presence of a background zonal flow U . Besides, we adopt a local approach by assuming a length scale separation between the wavelength of the forcing and the spatial variations of the zonal flow. We also assume homogeneity by considering an infinite fluid domain in both directions. This approach allows us to -forget about the geometrical effects inherent to the cylindrical geometry and work in equivalent 2D local cartesian coordinates (x, y) (see figure 1); -assume that the background flow U is constant in (x, y), which is only true locally, inside of a single jet.

For the basis (e x , e y , e z ) to be direct, with e z downward and e x zonal, in the same direction as e φ , e y has to be oriented towards the axis of rotation (figure 1). We denote u = (u, v) e x ,ey the 2D cartesian velocity components, and ζ = ∂ Note that in this cartesian framework β is now positive (d y h < 0). We have added an arbitrary stationary forcing q(x, y) representing a vorticity source. We linearize this equation around a uniform background zonal flow U = U e x by setting u = U + u ′ with |u ′ | ≪ |U |, and keeping only first order terms:

∂ζ ′ ∂t + U ∂ζ ′ ∂x + βv ′ + αζ ′ = ν∇ 2 ζ ′ + q(x, y). (4.12)
We drop the primes in the following and define the streamfunction

ψ (u = -∂ y ψ, v = ∂ x ψ and ζ = ∇ 2 ψ) such that ∂ ∂t ∇ 2 ψ + U ∂ ∂x ∇ 2 ψ + β ∂ψ ∂x + α∇ 2 ψ -ν∇ 2 ∇ 2 ψ = q(x, y). (4.13)
We perform a spatial Fourier transform of this equation in both x and y leading to

∂ ψ ∂t + [ i ω(k, l) + ω E (k, l) ] ψ = - q(k, l) k 2 + l 2 , (4.14)
where ψ and q are the Fourier coefficients associated with ψ and q, and k = (k, l) e x ,ey the wave vector. We denote ω the Rossby waves frequency, Doppler-shifted by the advection by the zonal flow

ω = kU - kβ k 2 + l 2 ,
(4.15) (Vallis 2006, p.230) and ω E the damping rate due to the viscous dissipation in the bulk and the bottom friction

ω E = α + ν(k 2 + l 2 ). (4.16)
Note that there is no gravity effects (or deformation radius) in the Rossby waves dispersion relation because we make the rigid lid approximation (appendix A). This is justified in the present work since the short wave limit is valid for the forced waves. The solution to equation (4.14) with the initial condition ψ(k, l, t = 0) = 0 is

ψ(k, l, t) = -q(k, l) (k 2 + l 2 )(iω(k, l) + ω E )
1 -e -(iω+ωE)t .

(4.17)

The inverse Fourier transform F -1 of ψ can be computed numerically to retrieve the physical streamfunction ψ at a time t for a given forcing. Similarly, the vorticity and velocity components can be computed using

           ζ(x, y, t) = F -1 -(k 2 + l 2 ) ψ(k, l, t) , u(x, y, t) = F -1 -il ψ(k, l, t) , v(x, y, t) = F -1 ik ψ(k, l, t) .
(4.18)

The Reynolds stresses term R(y, t; U ) is then easily computed as

R(y, t; U ) = - ∂ ∂y uv x = vζ x . (4.19)
The sign difference in the second equality compared to expression (4.8) in cylindrical coordinates comes from e y pointing inward whereas e ρ is pointing outward.

Comparison of the linear model with experimental results

To confirm that the reduced QG approximation is an appropriate model for the experiment, we compare the linear solution with the very beginning of experiments where only one forcing ring is turned on. Note that a good agreement is expected since in our experiments, Ro ζ = ζ/f ≪ 1 (table 1), which is the main assumption of the QG approximation. To carry this comparison, we first set all the model parameters to the experimental ones, that is

         β = f h dh dy ∼ 50 m -1 • s -1 , α = 1 2 f E 1/2 ∼ 5.6 × 10 -3 s -1 , U = 0 m • s -1 . (4.20)
When the pump is activated, the fluid is at rest in the rotating frame (solid body rotation). In addition, we design the forcing term to mimic the experimental forcing on the chosen ring: for the third ring, it corresponds to a line of 18 vortices (9 cyclones and 9 anticyclones) regularly spaced onto a perimeter of 2πR 3 ≈ 0.688 m. We chose to represent each vortex by a Gaussian one, leading to

q(x, y) = q m 18 i=1 (-1) i exp - (x -x i ) 2 + (y -y i ) 2 r 2 v , (4.21)
with q m the forcing amplitude, r v the radius of the vortices and (x i , y i ) the centre location of each vortex. The vortices' radius is set to 2/5 times the spacing between two vortices based on the experimental measurements. To estimate the forcing amplitude q m that we should use to better represent the experimental regime, we measure the vorticity linear growth rate above the forcing injection points and adjust q m so that the growth rate obtained with the linear model is comparable to the experimental one. This method leads us to use q m = 0.5 s -2 . Figure 9 shows a comparison between an experiment and the linear model 2 seconds after the third forcing ring was turned on at its maximum power. The experimental flow and the linear solution are qualitatively and quantitatively very close. They both exhibit a westward stretching of each vortex, westward meaning in the retrograde direction compared to the background rotation (decreasing φ in the experiment, decreasing x in the cartesian model). The dispersive nature of the Rossby waves emitted by the vortex are responsible for this chevron pattern pointing eastward, as explained by [START_REF] Firing | The Behavior of a Barotropic Eddy on a β-Plane[END_REF]; [START_REF] Flierl | The Application of Linear Quasigeostrophic Dynamics to Gulf Stream Rings[END_REF]; [START_REF] Chan | Analytical and Numerical Studies of the Beta-Effect in Tropical Cyclone Motion. Part I: Zero Mean Flow[END_REF] in the case of isolated vortices. In particular, the long Rossby waves which propagate westward faster than short waves are responsible for the westward stretching. This response can also be understood as the deformation of each vortex into a so-called β-plume. [START_REF] Stommel | Is the South Pacific helium-3 plume dynamically active?[END_REF] first described β-plumes when trying to understand the circulation induced by rising water from hydrothermal vents in the Pacific, and his theory was further developed by [START_REF] Davey | Flows Produced by Discrete Sources of Buoyancy[END_REF] who considered the evolution of buoyancy sources on a βplane. In both cases, the convergence or divergence of fluid around the perturbation generates cyclonic or anticyclonic motions which are subsequently elongated westward due to the emission of Rossby waves (see the review and experiments by Afanasyev & Ivanov (2019)). As a consequence, an east-west asymmetry in the radial component develops, which is clear in figure 9(b,f ): the advection of the background potential vorticity leads to a weakening of the flow on the west (left) side of each vortex, and a strengthening on their east (right) side, for both cyclones and anticyclones. More interesting for us is the fact that this asymmetry leads to a prograde momentum convergence towards the region of generation of the vortices as demonstrated by figure 9(d,h). Figure 10 compares the Reynolds stresses divergence profile in the experiment and in the linear solution. There is indeed an eastward acceleration of the zonal flow in the forcing region, located between two westward acceleration regions. This mechanism thus explains the experimental regime I, i.e. the formation of prograde jets flanked by two retrograde jets above each forcing ring. Vallis (2006, chapter 12) provides an overview of this mechanism in the framework of barotropically forced surface westerlies in the atmosphere. Our experiment then shows that this mechanism is robust since the generated prograde jets persist at later times, even in the non-linearly saturated regime.

Resonance of the Reynolds stresses and associated feedback

The previous section shows that our experiment can be successfully described by a simple 2D QG model incorporating only the β-effect and the bottom friction, at least in the linear regime. We can thus use this model to investigate the feedback that the zonal flow can have on the Reynolds stresses divergence R and study the possibility of bi-stability. This is the goal of the present section.

For simplicity, we now forget about the specific geometry of the experimental forcing, and use a generic forcing consisting of a doubly-periodic array of vortices with a wavelength comparable to the experimental one:

q(x, y) = q m cos(k f x) cos(l f y), (4.22)
with q m = 0.5 s -2 the forcing amplitude and k f = l f the forcing wavenumber. In the following, we work with k f = 63 rad • m -1 corresponding to a typical forcing wavelength of 10 cm. We show in appendix C that the forcing scale has only a small influence on the physical mechanism presented. The linear response obtained with this forcing is the superposition of the response to each forcing line, with a prograde acceleration above each forcing horizontal line, and retrograde accelerations in between.

To study the feedback of the zonal flow, we solve for the stationary linear response to this forcing (equation (4.17) with t → +∞) for various amplitudes of the background zonal flow U . We recall the local nature of our analysis: the background zonal flow U should be seen as the -uniform -zonal flow at the core of a jet. For each solution, we extract the stream function amplitude | ψ| and the Reynolds stresses divergence at y = 0, and represent them in figure 11 as a function of the zonal flow U . Given the amplitude of the stream function

| ψ| = q m k 2 f + l 2 f ω(k f , l f ) 2 + ω E (k f , l f ) 2 -1/2 , (4.23)
we expect a resonance of the linear response when ω = 0, or, in other words, when the directly forced Rossby waves are stationary, i.e.,

U = β k 2 f + l 2 f = -c (4.24)
where c is the directly forced Rossby waves phase speed in the absence of zonal flow (Vallis 2006, p.542). The resonant amplification of the response amplitude is visible in figure 11(a). Then, the amplitude of the Reynolds stresses is expected to be proportional to the squared amplitude of the streamfunction (equations (4.18) and (4.19)), leading to

|R|(U ) ∝ | ψ| 2 ∝ 1 ω 2 + ω 2 E ∝ 1 1 + U c 2 + γ 2 , (4.25)
where γ is a nondimensional parameter characterizing the Rossby waves damping

γ 2 = ω E k f c 2 . (4.26)
For the problem to be analytically tractable, we chose to model the resonant curve R(U ) plotted in figure 11(b) with a parametrized Lorentzian. We thus forget about the spatial structure of the momentum flux convergence, and set

R(U ) = R m 1 γ 2 + 1 + U c 2 . (4.27)
Doing so, we focus on the amplitude of the response rather than on the details due to our particular choice of forcing. The important physical effects of the zonal flow, the β-effect and the friction, are contained in the Lorentzian and influence the position and flatness of the resonant peak. Figure 11(b) shows that the amplitude of the Reynolds stresses is indeed largely dominated by this resonant amplification. Hence, we do not loose any important feature by modeling R with equation (4.27) which has the advantage of making the problem analytically tractable.

It is now clear that the momentum flux can lead to abrupt transitions: on the left side of the resonant peak, any increase of a zonal flow U leads to an increased prograde momentum convergence, and an increased acceleration of this zonal flow. We have thus identified a potential positive feedback mechanism, provided that it is not cancelled by the negative feedback of the viscous dissipation. Note that since the Rossby waves propagate in the westward direction, this resonance can only occur in an eastward jet (U > 0).

Stationary solutions and linear stability

The linear QG model showed the resonant amplification of the wave-induced Reynolds stresses when the zonal flow is such that the directly forced Rossby waves are stationary. Our goal is now to verify whether this feedback of the zonal flow can explain the transition and bi-stability observed in our experiment.

Consistently with our local approach, we consider a minimal model where the zonal flow U (t) is assumed to be only time-dependent, with no spatial modulation. Such a uniform zonal flow is sustained by the Reynolds stresses divergence R and dissipated by the linear friction due to the Ekman pumping:

∂U ∂t = R(U ) -αU, (4.28)
with R(U ) the Lorentzian given by equation (4.27). The stationary solutions of this equation are the roots of the third-order polynomial

P (U ) = U 3 + 2c U 2 + (1 + γ 2 )c 2 U - R m c 2 α . (4.29)
Depending on the sign of the discriminant of P , one or three stationary solutions can exist, as represented in figure 12. We denote U 1 , U 2 and U 3 those three solutions such that U 1 < U 2 < U 3 . For three stationary solutions to exist, i.e. bi-stability to be possible, a necessary but not sufficient condition is

γ 2 < 1 3 . (4.30)
When this condition is satisfied, the sufficient condition for three solutions to exist is that

R m ∈ [R 1 , R 2 ] = - 2 27 αc 9γ 2 + 1 - √ Γ , 9γ 2 + 1 + √ Γ , (4.31) (a) (b) (c) (b) (a) (c) (d)
Figure 12. Visualization of the stationary solutions (U1, U2, U3) of the zonal flow evolution equation (4.28). Except the forcing amplitude, all the parameters are fixed and equal to the experimental parameters (α = 5.6

× 10 -3 s -1 , k f = l f = 63 rad • s -1 , c = 6.3 × 10 -3 m • s -1 , γ 2 = 1.99×10 -4
). (a) Illustrative case of a small forcing amplitude Rm < R1 and a small friction or a too high friction. Note that U1 is very small, but not zero. (b) Illustrative case of a forcing in the range Rm ∈ [R1, R2] and a small friction. (c) Illustrative case of a high forcing amplitude Rm > R2 and a small friction. (d ) Amplitude of the three stationary solution as a function of the forcing amplitude Rm. The black dots represent the unstable solution whereas the white ones are stable. The shaded area is the bistable range, bounded by R1 and R2 (equation 4.31).

with Γ = 1 -3γ 2 3 . Physically, the first condition (4.30) means that bi-stability can never exist if the Rossby waves are too strongly damped. The second condition shows that even when the friction is not too high, three stationary solutions exist only for a given range of the forcing amplitude R m . As represented in figure 12, if the forcing is too high, only the super-resonant solution U 3 exists. Conversely, if the forcing is too weak, only the low-amplitude solution U 1 can exist.

To investigate the linear stability of the stationary solutions, we go back to the zonal flow evolution equation (4.28). We linearize the nonlinear operator R around the stationary state U s and denote U ′ (t) = U (t) -U s the perturbed zonal flow, to obtain the perturbations evolution equation

∂U ′ ∂t = dR dU Us U ′ -αU ′ . (4.32)
We seek U ′ under the form U ′ = U ′ 0 e σt where σ is the growth rate of the perturbation. Substituting into equation (4.32) leads to

σ = dR dU Us -α. (4.33)
The stationary solution U s is unstable if and only if the growth rate σ > 0. This condition is always verified for the second stationary solution U 2 (the sub-resonant one), whereas U 1 and U 3 are stable stationary solutions. We plot in figure 12(d ) the amplitude of the stationary solutions obtained for varying forcing amplitudes R m , all the other parameters being fixed and equal to the experimental parameters (α = 5.6

× 10 -3 s -1 , k f = l f = 63 rad • s -1 , c = 6.3 × 10 -3 m • s -1 , γ 2 = 1.99 × 10 -4
). The three stationary solutions branches are visible, and coexist only in a given range of forcing amplitude bounded by R 1 (purple line) and R 2 (blue line) given by equation (4.31). This figure also demonstrates the bi-stability and possibility of an hysteresis: for an experiment with increasing forcing, the transition U 1 → U 3 happens for R m = R 2 . At this forcing, there is a saddle-node bifurcation (S 2 ) through which U 1 loses its stability. But if the forcing amplitude is then decreased, the observed solution will remain U 3 until the forcing reaches R 1 < R 2 where there is another saddle-node bifurcation (S 1 ), and we go back to the lower branch solution U 1 . Our model predicts that close to the transition, the amplitude of the zonal flow on the lower branch is of U ≈ 2 mm/s, and U ≈ -c = 6.3 mm/s for the upper branch. Finally, this model allows us to define a typical velocity expected at the transition, to compare with the forcing RMS velocity U f used to characterize the experimental hysteresis, on figure 7. From equation (4.31), the Reynolds stresses at the transition are typically of R m,t ∼ α|c|(1 + 9γ 2 ). A typical transition velocity can be obtained supposing that the forcing is balanced by friction U t ∼ R m,t /α, leading

t ∼ |c|(1 + 9γ 2 ). (4.34)
If we use U f = U f /U t as a non-dimensional forcing, then the transition should always occur at U * f of order unity. The dependence on the β-effect, the friction, and the forcing scale are incorporated in c and γ. Note that the width of the bistable zone will however vary depending on γ.

Comparison of the experimental to theoretical transition

In this section, we report additional experimental observations that support the mechanism of Rossby waves resonance.

First, the amplitudes of the zonal flow expected in the two steady states U 1 and U 3 are very close to the experimental ones, where U 1 represents regime I and U 3 regime II, as can be seen in figure 7. Another way to see it is by comparing the zonal flow amplitude with the phase speed of the -non-advected -directly forced Rossby waves: if it is higher than -c, the system is in the super-resonant steady state, whereas if it is lower than -c, the system is sub-resonant. In the experimental setup, with k f ≈ 63 rad • m -1 and β ≈ 50 m -1 • s -1 , the phase speed of the directly forced Rossby waves is c ≈ -6.3 mm • s -1 , without advection. First, figure 13 demonstrates that the forcing indeed excites Rossby waves. The radial component of the velocity at a given radius exhibits patterns that move in the retrograde direction (i.e. decreasing φ), with the same wavelength as the forcing and at the doppler-shifted speed u φ φ + c. Then, given this typical phase speed c, figure 7 shows that in regime I the zonal flow is sub-resonant (u RMS < -c), and super-resonant in regime II (u RMS > -c), which is consistent with our model. Furthermore, the model predicts that when the forcing is decreased and U 3 (regime II) loses its stability (S 2 in figure 12(d )), the zonal flow is quasi-resonant ( u φ φ ≈ -c ) which is again compatible with the measured velocity at the transition II→I in figure 7. Note that if this equilibrium close to resonance can exist (U 3 ), our analysis only gives an explanation for its origin, ultimately, the non-linearities are responsible for locking the system into a near-resonant state. Note finally that the bi-stability range and zonal flow amplitudes are slightly varying with the range of possible forcing wavenumbers for our experiment (figure 16(b) in appendix C). For a closer match between the flow amplitudes the measured ones, we would choose k f = 57 rad • m -1 , which is reasonable given the uncertainty in the relevant forcing scale in our setup. Second, we have not yet discussed the inherent asymmetry between the prograde and the retrograde jets in our model. Since the Rossby waves propagate in the retrograde direction, their only way to become stationary and resonate with the fixed forcing is within a prograde jet. This implies that the transition only occurs in the prograde jets which is again consistent with our experimental observations. Indeed, the prograde jets govern the dynamics, by increasing in amplitude and merging, whereas the retrograde flow seems to adjust passively to the transition (see e.g. figure 6). If we suppose that the forcing imparts no net angular momentum to the fluid, then the integral of the angular momentum per unit mass, ρ(u φ + Ωρ), over the domain should be constant. If an eastward acceleration is produced by the transition, then this convergence of prograde angular momentum should be balanced by negative angular momentum elsewhere (see the Reynolds stresses profile in figure 10). The retrograde flows in our experiment seem to arise as such, which is consistent with the fact that the retrograde flow is smooth whereas the prograde one strongly interacts with vortices (figure 3).

Third, our model implies that the faster the Rossby waves, the stronger the forcing will have to be to reach the transition (increasing U t ), and conversely for slower Rossby waves.

We performed experiments at 80 RPM and 60 RPM (not shown) instead of the 75 RPM for which the bottom plate was designed. In these experiments, β is no longer uniform, but slightly varying with radius. It is higher at 80 RPM due to the increased curvature of the free-surface (β 80 ∈ [57, 73] m -1 • s -1 with a mean at 65.5 m -1 • s -1 ) and weaker at 60 RPM (β 60 ∈ [20, 30] m -1 • s -1 with a mean at 22.8 m -1 • s -1 ). As a consequence, the forced Rossby waves are respectively faster and slower (|c| = 8.3 and 2.9 mm • s -1 ). Consistently with our model, we observed that for a similar forcing (U f = 4.0 mm • s -1 ), regime II is obtained at 75 RPM whereas regime I is observed at 80 RPM. Conversely, for a forcing at which regime I is observed at 75 RPM (U f = 2.9 mm • s -1 ), regime II is obtained at 60 RPM. The transition thus occurs at larger forcing amplitudes for an increased β-effect. For a more quantitative view of the sensitivity to the β-effect, we refer the reader to appendix C. Finally, we wish to discuss the local aspect of our model relatively to the experiment. Our model only explains the local feedback mechanism inside of a given prograde jet. How the global system responds is a different question. Based on the Rhines scale L R ∼ (u RMS /β) 1/2 [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF], where we recall that u RMS takes into account all the components of the flow (equation (3.1)), we do expect an increase of the jets width during the transition, since the RMS eddy velocity increases. However, our model does not explain why the jets merge during the transition. In addition, this model does not rule out the possibility of the coexistence of regime I and regime II flows side by side. For instance, the most external forcing ring in our experiment (above C6) has a greater pressure loss because of the number of hoses (38) and it is possible that on this ring, the forcing is never super-resonant. This suggests that regional stable equilibria, where both regimes can be locally sustained in distinct regions of space, may exist. Finally, we observe that the merger events during the transition are associated with a radial shift of the jets leading to an uncorrelation between the jets position and the forcing rings. This further suggest that a local approach will not be sufficient to explain the saturation in regime II. Instead, it may be relevant to adopt a global approach based, for instance, on the turbulent properties and energy transfer of anisotropic turbulence on a β-plane. This type of approach has led to the development of the theory of zonostrophic turbulence (Galperin et al. 2019, and references therein). Due to its fast rotation and owing to the Taylor-Proudman theorem, the flow is quasi two-dimensional and may bear an inverse turbulent energy cascade. Because of the β-effect and associated Rossby waves, the energy transfer becomes anisotropic and redirected towards zonal currents. Ultimately, the large scale drag halts the expansion 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]. Determining whether such a theory is valid to explain the non-linear saturation in regime II is beyond the scope of the present work, but will be investigated in a separate study. Note also that the spectral analysis which can be found in [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] for the previous version of the experiment and corresponding DNS supports the idea that the second regime is close to zonostrophic turbulence. To conclude, we wish to underline that [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] probably only observed regime II in their experimental setup because their forcing was of larger amplitude than in the present study, thus probably always super-resonant.

Conclusions and discussion

Experimental conclusions and future work

We have described an experimental setup capable of generating robust zonal jets even in the presence of boundary dissipation. In this setup, we observed a subcritical transition between two different steady states with instantaneous zonal flows. In the first regime, obtained for a weak forcing and a moderate local Reynolds number, the jets are steadily forced by prograde momentum convergence towards the eddy-forcing regions, through the indirect action of Reynolds stresses. In the second regime, obtained for a strong forcing and larger Reynolds numbers, the jets merge into higher amplitude zonal flows at a larger scale. While the two regimes are obtained at different Reynolds numbers, they both correspond to low Rossby number QG dynamics. The two regimes coexist in a small forcing range, leading to bi-stability, and we are able to follow the corresponding hysteresis cycle. The transition is found to be due to the resonance occurring when the forced Rossby waves become stationary because of their advection by the zonal flow. Note that, in the present work, we explain the bistability with the linear resonance mechanism originally explored by [START_REF] Charney | Multiple Flow Equilibria in the Atmosphere and Blocking[END_REF], which predicts two stable states with different waves and zonal flow amplitudes. The bending of the same resonance due to weakly non-linear effects [START_REF] Benzi | The statistical properties of general atmospheric circulation: Observational evidence and a minimal theory of bimodality[END_REF][START_REF] Malguzzi | Nonlinear Amplification of Stationary Rossby Waves Near Resonance[END_REF][START_REF] Malguzzi | Nonlinear Amplification of Stationary Rossby Waves near Resonance[END_REF] is also a potential candidate to account for the observed bistability. However, in this framework of nonlinear resonance, the two regimes are expected to have similar zonal flow velocities, which is in contradiction with our observations.

In laboratory experiments, bifurcations involving multiple zonal flows steady states have been observed only a few times. [START_REF] Weeks | Transitions Between Blocked and Zonal Flows in a Rotating Annulus with Topography[END_REF]; [START_REF] Tian | Experimental and numerical studies of an eastward jet over topography[END_REF] observed bi-stability in the context of mid-latitude atmospheric jets, following the same resonance but with topography. However, in these experiments, the zonal flow is directly forced by pumping fluid in at a larger radius than where it is pumped out. Bifurcations over indirectly forced zonal flows were observed by [START_REF] Semin | Nonlinear saturation of the large scale flow in a laboratory model of the quasibiennial oscillation[END_REF] in their experimental model of the quasi-biennal oscillation, but in that case, the flow is laminar, and the low forcing amplitude state is a state with no mean flow. In the present study, we describe bi-stability between two steady states sustaining indirectly forced and multiple zonal jets. Let us mention here that bi-stability has also been observed numerically in the context of rotating thermal convection where zonal flows emerge due to the sphericity of the domain. At intermediate Ekman numbers, the saturation of the convective instability can lead to either a weak branch or a strong branch, both supporting zonal flows but which are much more vigorous on the strong branch [START_REF] Guervilly | Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model[END_REF][START_REF] Kaplan | Subcritical Thermal Convection of Liquid Metals in a Rapidly Rotating Sphere[END_REF]. The question whether this bistability may be linked to a somewhat similar resonance involving thermal Rossby waves is however far from trivial. For instance, there is no obvious reason for the convective eddies that force the zonal flow to be stationary, and to which extent they are decoupled from the zonal flow remains to be determined.

The results presented in this paper raise several questions which will be investigated in future work. First, we wish to study the specificity of the observed transition to our forcing pattern. Specifically, we want to investigate if this mechanism can hold when we break the coherence of the vortices generated by our forcing, for instance by adding elbow connectors to the inlets and outlets. We also plan to replace the pumps with more powerful ones to reach more extreme regimes: the present values of the Rossby number (table 1) offer us the possibility to force stronger flows while remaining in regimes strongly constrained by rotation.

Then, explaining the transition's origin does not explain the non-linear saturation in the second regime, and its evolution as a global system. Understanding the final equilibrated state in this regime, which is closer to the planetary ones, thus remains a challenge. For instance, we mentioned that regime II is in fact multistable, with at least three different jets configurations identified based on 15 realizations (figure 7(b)). The origin of this multistability remains to be elucidated. In particular, it would be of interest to compare it with the multistability observed in numerical simulations of stochastically forced barotropic turbulent jets, where nucleation and coalescence are associated with transitions between steady states with a different number of jets (Bouchet et al. 2019b;[START_REF] Bouchet | Zonal Flows as Statistical Equilibria[END_REF].

Future work is also needed to explain the long-term dynamics (or stability) of the zonal flows. Such a task requires to understand the complex interaction between the small-scale transient turbulence and the slowly varying zonal flow, which is difficult given the time-scale separation between the two. In this regard, laboratory experiments have a major role to play since they allow for measurements at high-resolution over long times. In our case, the long-term experiments (38 000 t R ) that we performed show no radial migration of the jets, except during merger events. The question whether this stability is due to our axisymmetric forcing, or to the uniform β-effect remains to be addressed. In this regard, it may be interesting to go back to a setup close to that of Cabanes et al. (2017) and explore the precise effect of a non-uniform β-effect by keeping the exact same forcing with a flat bottom. In the literature, drifting, merging and nucleation of jets have been described in various numerical models, e.g. in the framework of rotating thermal convection [START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF] and stochastically forced barotropic jets on the β-plane (Bouchet et al. 2019b). However, such long-term dynamics is not that common in natural systems. For instance, jets on the gas giants are remarkably steady [START_REF] Tollefson | Changes in Jupiter's Zonal Wind Profile preceding and during the Juno mission[END_REF]), even if one jovian jet seems to have broken apart into vortices in 1939-1940 [START_REF] Rogers | The Giant Planet Jupiter[END_REF][START_REF] Youssef | The dynamics of jovian white ovals from formation to merger[END_REF]. Similarly, in laboratory experiments, meridional jet migration has only been described by [START_REF] Smith | Multiple Zonal Jets in a Differentially Heated Rotating Annulus* ,+[END_REF], and the authors underline that it may be due to a long-term thermal equilibration. Despite the apparent stability of the jets in our experiments, we would like to mention that the prograde jets in regime II sometimes have an long-term fluctuating behaviour with the repetition of cycles where the jets destabilizes into vortices before recovering their initial state, as illustrated by figure 14 and supplementary movie 4. More precisely, we observed the growth of zonal perturbations of the zonal flow (figure 14(b)) followed by vortices "surfing" along the jet in the prograde direction (figure 14(c-e)). These zonal packets of vortices may correspond to envelope Rossby solitary waves, and are also reminiscent of non-linear waves called "zonons" which have been described within barotropic jets in numerical simulations [START_REF] Sukoriansky | Nonlinear Waves in Zonostrophic Turbulence[END_REF][START_REF] Galperin | Geophysical flows with anisotropic turbulence and dispersive waves: Flows with a β-effect[END_REF][START_REF] Sukoriansky | Rossby waves and zonons in zonostrophic turbulence[END_REF][START_REF] Bakas | Emergence of Large Scale Structure in Barotropic β-Plane Turbulence[END_REF]. This long-term dynamics will be the focus of future experiments.

Finally, we wish to place our experiment in the framework of zonostrophic turbulence previously mentioned (Galperin et al. 2019, and references therein). We have already underlined in the introduction that the zonostrophy index in our experiment is larger than previous experimental studies thanks to the fast rotation, and hence we are closer to the gas giants' regime. [START_REF] Cabanes | A laboratory model for deep-seated jets on the gas giants[END_REF] provide details about this index estimation, which can also be found in Appendix D for the present experiments. We have not yet discussed a second criterion which is that the forcing should act at a scale smaller than the scale at which the eddies start to feel the β-effect for a significant Kolmogorov-Kraichnan inertial range to exist and an isotropic inverse cascade to develop. In environmental flows (oceans and planetary atmospheres), typically, the forcing acts at a scale smaller than the scale of turbulence anisotropisation L β by a factor 2 to 3 (Galperin et al. 2019, table 13.1). In our experiment, we can estimate that the forcing scale is in fact about twice L β . The forcing is thus directly influenced by the β-effect, which is quite clear on the fluid response at the earliest times (figure 9). It is thus probable that we prevent an isotropic inverse energy cascade -an anisotropic cascade is probably present since the jets scale remains larger than the forcing scale. We might thus stand in the case where the process of jet formation can be considered independently from the two-dimensional inverse energy cascade. That being said, as demonstrated by [START_REF] Scott | Zonal Jet Formation by Potential Vorticity Mixing at Large and Small Scales[END_REF] in the framework of potential vorticity mixing, the late-time resulting jets profile is the same whether the forcing is performed at large scale or at small scale, the important parameter being the zonostrophy index R β . Here, R β ∈ [2.25, 2.84] (see appendix D for this estimation) which shows that our experimental regimes are not friction-dominated and explains that the observed jets are highly energetic and instantaneous. If nevertheless one wishes to reach the regime of a well-developed turbulent cascade, one way to increase the scale separation is to have stronger flows. However, increasing U by a factor 10 would only increase L β by a factor 10 2/5 ≈ 2.5 (equation (D 3)). The best option in order to study jets sustained by an inverse cascade would be to decrease even more the forcing scale, which is then a further experimental challenge.

Modeling and relevance for planetary systems

Finding an explicit expression for the Reynolds stresses is the basis of the outof-equilibrium statistical theories aiming at explaining zonal jets formation from an homogeneous turbulent flow (see e.g. [START_REF] Bakas | Emergence of Large Scale Structure in Barotropic β-Plane Turbulence[END_REF]. Indeed, it yields to a closed, deterministic system for the zonal flow dynamics. Here, a very simple framework based on the QG approximation is sufficient to explain our experimental observations but more sophisticated models where, for instance, the spatial modulations of the zonal flow are taken into account are necessary to help understand observations at a global scale.

In the present study, the key mechanism is the local resonant amplification of the Reynolds stresses by the zonal flow. It has been applied previously in two main geophysical frameworks: mid-latitude atmospheric jets and equatorial super-rotation, which is a state with a strong prograde jet at the equator. For the Earth's atmosphere, the Rossby waves resonance has been successfully employed to explain abrupt transitions of the jet stream between blocked and zonal flows [START_REF] Charney | Multiple Flow Equilibria in the Atmosphere and Blocking[END_REF], and it is now considered as a valuable candidate to explain extreme weather events in the past 20 years [START_REF] Petoukhov | Quasiresonant amplification of planetary waves and recent Northern Hemisphere weather extremes[END_REF][START_REF] Coumou | Quasi-resonant circulation regimes and hemispheric synchronization of extreme weather in boreal summer[END_REF]. Then, in the case of equatorial jets, abrupt transitions to super-rotation and bi-stability have been observed in global climate models and numerical simulations. The same wave-jet resonance feedback as for mid-latitude jets, which arises in response to a stationary equatorial heating, has been recently considered as a robust mechanism for this transition [START_REF] Arnold | Abrupt Transition to Strong Superrotation Driven by Equatorial Wave Resonance in an Idealized GCM[END_REF][START_REF] Herbert | Atmospheric Bistability and Abrupt Transitions to Superrotation: Wave-Jet Resonance and Hadley Cell Feedbacks[END_REF]. The fact that in those two frameworks, the resonance successfully explains observations at a global scale is encouraging. However, to the best of our knowledge, such a mechanism has never been studied for its potentiality to generate strong zonal jets in the broader context of an eddy-forcing, neither for its applicability to extra-equatorial jets on the gas giants or in the Earth's oceans.

First, we can briefly compare the zonal flow amplitudes relatively to the Rossby waves intrinsic phase speed (equation (4.24)). If β can be quite easily estimated for planetary systems, the relevant wavelength for the "forced" Rossby waves is not trivial. In addition, the gravity effects can be of importance for planetary applications, especially if the considered jets are shallow, and the Rossby radius of deformation should be reintroduced in the phase speed equation (4.24). For the mid-latitude jet stream, as previously mentioned, bi-stability has been observed. The zonal flow is either subresonant or super-resonant with topographically forced waves phase speed of typically 16 m/s in [START_REF] Charney | Comparison of a Barotropic Blocking Theory with Observation[END_REF]. In the Earth's oceans, if we assume typical phase speeds of 0.25 m/s (Vallis 2006, p.233), then the zonal flow could be sub-resonant since zonal jets in the ocean have typical speeds of a few centimetres per second [START_REF] Cornillon | Oceans[END_REF]. For the gas giants, we use the values reported in Galperin et al. (2019, table 13.1). The β-effect is β ≈ 3 × 10 -12 m -1 s -1 , and the deformation radius is of about 2,000 km [START_REF] Vasavada | Jovian atmospheric dynamics: An update after Galileo and Cassini[END_REF]. Using a large forcing wavelength based on the transitional scale, k f = l f = 2πL -1 β with L β ≈ 6, 000 km (smaller wavelengths would lead to lower phase speed), we obtain a phase speed of only 1.2 m/s, meaning that all the jets would be super-resonant since the typical jets velocities are of about 50 m/s [START_REF] Sánchez-Lavega | Gas Giants[END_REF]. Let us stress out that for the present mechanism to hold in planetary systems, there is the need for a partial decoupling between the forcing source and the jets, such that the forced waves advected by the zonal flow can become resonant with the forcing. For the Earth jet stream, this decoupling comes from the fact that the topography exciting the waves is fixed, just like in our experiment. For Jupiter, the forcing origin is not clear. It can take place in the weather layer due to moist convection or band-to-band horizontal contrasts in heating, but it can also arise from the deep molecular convective interior of the gas giant [START_REF] Vasavada | Jovian atmospheric dynamics: An update after Galileo and Cassini[END_REF]. At which speed these structures propagate relatively to the zonal flow is certainly not clear. All these questions require dedicated studies, but the important point is that even with an unsteady forcing, propagating azimuthally, the resonance mechanism should still hold provided that the forcing is not passively advected by the zonal flow it generates.

In addition to this simple velocities comparison, the important parameter of the model is γ 2 = (α/(k f c))

2 . This parameter compares the Rossby waves period 1/(k f c) to the friction timescale 1/α, and we have shown that it should be small (equation (4.30)) for the super-resonant solution or bi-stability to exist. The question whether such a mechanism is expected for extra-tropical jets in planetary flows is beyond the scope of the present study and would require an extensive systematic study. Besides, one should properly define the bounds of the physical parameters for the model to still be self-consistent. We recall for instance that the model is based on the linear response to a stationary forcing. Finally, determining the relevant dissipation parameter for the Rossby waves is not trivial either, and it cannot be reduced to a simple Ekman friction like in our experiment. That being said, for completeness, we illustrate in appendix C the sensitivity of the bi-stability to the model parameters (α, β, k f ). Let us briefly mention the case of the friction α, shown in figure 16(a). Interestingly, the bistable range is shifted towards lower values of the forcing amplitude as α → 0, meaning that for infinitely small friction, the super-resonant solution U 3 (regime II) would be obtained even at a very small forcing amplitude, while the sub-resonant solution U 1 (regime I) would never be observed. Regime II is thus expected in most planetary applications. scale at which the inertial term equates the β term:

L R ∝ u rms β 1/2 , ( D 1)
where u rms usually implies the root mean square eddy velocity. Under this scale, the advective term dominates and above it the β-term is dominant. Originally defined in [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF], it has often been associated with the width of zonal jets. A second scale can be defined in the context of two-dimensional turbulence by equating the eddy turnover time to the Rossby wave period

L β ∝ ǫ β 3 1/5 , (D 2)
where ǫ is the rate of upscale energy transfer within the turbulent energy cascade. This scale characterizes the threshold of turbulence anisotropisation under the action of the βeffect. [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] demonstrated that the strength of the jets and the quality of their delineation can be classified depending on the zonostrophy index R β defined as the ratio

R β = L R L β = β 1/10 (u rms ) 1/2 ǫ -1/5 ≈ βu rms Ω 2 1/10 E -1/10 (D 3)
where ǫ is estimated from the rate of energy loss due to dissipation ǫ ≈ u 2 /τ E with τ E = Ω -1 E -1/2 the Ekman spin-down timescale. 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 Rossby waves (L β ) is well separated from the scale of the final jets (L R ), i.e. for large zonostrophy index, R β > 2.5 [START_REF] Galperin | Geophysical flows with anisotropic turbulence and dispersive waves: Flows with a β-effect[END_REF]. With the experimental parameters and the typical values of u rms provided in table 1 for the two experimental regimes, we obtain R β,I ≈ 2.26 and R β,II ≈ 2.84.

Figure 1 .

 1 Figure1. Schematic of the experimental setup. A cylindrical tank of 1m in diameter and 1.6m in height, filled with 600 liters of water, is fixed on a table rotating at 75 RPM. The fluid free-surface takes a paraboloidal shape due to the centrifugally-induced pressure. The bottom plate is designed to achieve a uniform topographic β-effect. A small-scale turbulent forcing is performed by circulating water through 128 holes at the base of the tank. The forcing pattern is sketched on the right: each ring C1-C6 is controlled by an independent pump. Time-resolving PIV measurements are performed on a horizontal plane using a side green laser and a top-view camera.

Figure 2 .

 2 Figure 2. (a) Experimental setup at rest. The tank is mounted on a rotating table operating using an air cushion. It is filled with ∼600 liters of tap water so that the fluid height at rest above the bottom plate is of ∼58 cm. It is closed by a Plexiglas top-lid. (b) Experimental setup in solid body rotation at 75 RPM, with the side green laser turned on. The difference in fluid height between the center of the tank and its border is of ∼76 cm and the fluid height at the center is hmin = 20 cm. (c) View of the bottom plate through which the forcing is performed. The plate has the shape of a curly bracket for the fluid height to increase exponentially with the radius, see equation (2.3). It is drilled with 128 holes corresponding to 64 inlets and 64 outlets connected to 6 submersible pumps. (d ) View of the 6 pumps and 128 hoses placed beneath the bottom plate.

Figure 3 .Figure 4 .

 34 Figure 3. Instantaneous velocity fields computed from PIV measurements in the statistically stationary state reached in the two experimental regimes. The shaded areas in the top-right and bottom-left corners are shadow areas due to the laser refraction: no measurements are performed in these areas. The colors represent the vertical component of the vorticity ζ. Note that their is a factor ten between the color scales in the two panels. (a) Regime I. Velocity field obtained at time t =24 min = 1800 tR and averaged over 1s. (b) Regime II. Velocity field obtained at time t = 19 min = 1425 tR and averaged over 1s.

Figure 5 .

 5 Figure 5. Instantaneous maps in regime I (a-c) and II (d-f ). The black circle marks the tank boundary, and the shaded areas are the shadows where no measurements can be performed. The dashed curves in the top-left quadrant of each subplot represent the forcing rings location. (a-c) Regime I at time t = 24 min = 1800 tR and averaged over 1s. (a) Azimuthal component of the velocity u φ . (b) Radial component of the velocity uρ. (c) Vertical component of the vorticity ζ. (d-f ) Regime II at time t = 19 min = 1425 tR and averaged over 1s. (d ) u φ . (e) uρ. (f ) ζ. Note the different color scales for the two regimes.

  Figure 7(a) represents the mean flow amplitude (equation (3.5)) as a function of the forcing amplitude U f as defined in appendix B. Typical maps of the time-averaged zonal velocity are represented in panel (b).

Figure 6 .

 6 Figure 6. Experiment illustrating the bi-stability between regimes I and II: the forcing is initially of low-amplitude and stationary, so that we start in regime I (Pi = {14, 20, 46, 72, 100, 100}%, U f = 3.3 × 10 -3 m • s -1 ). At t = 360 tR, a finite perturbation is created by varying the third forcing ring power randomly around its initial value, leading to the transition to regime II (Pi = {14, 20, 46 ± 20, 72, 100, 100}%). At t = 1600 tR, the forcing is set back to its initial state but the flow remains in regime II. At time t = 2250 tR a second finite amplitude perturbation is performed (Pi = {14, 20, 46, 72 ± 20, 100, 100}%) (a). Hovmller plot: zonal flow profile as a function of time. (b) Total (K), zonal (Kz) and fluctuating (K f ) kinetic energy as a function of time.

Figure 7 .

 7 Figure 7. (a) Experimental hysteresis loop. The time-averaged RMS velocity of the flow u RMS (equation (3.5)) is plotted as a function of the forcing amplitude U f (see appendix B). The top and right axes correspond to associated non-dimensional quantities. On the right-axis, we use the Reynolds number based on the mean flow RMS velocity, Re = u RMS ℓ/ν (ℓ=7.3 cm is the distance between two forcing rings). For the top-axis, we use the typical velocity expected at the transition Ut (see section 4.4 and equation (4.34)). For each curve, the forcing is either increased (reddish) or decreased (bluish) step by step. The different colours are different experiments. The shaded area is the bistable zone. The grey points in regime II correspond to experiments initialized with the exact same forcing and for which saturation leads to three possible jets configurations represented in panel b. The last point of the yellow curve is in configuration 2, but may not have reach its stationary state. (b) Time-averaged zonal velocity maps in regime I and II. In regime II, multiple steady states can be obtained for the same forcing. The configurations 1, 2, 3 correspond to the points labeled accordingly in panel (a).

Figure 8 .

 8 Figure 8. Evolution of the total kinetic energy K during transitions from regime I → II (red curves) and II → I (blue curves). The normalized total kinetic energy is plotted as a function of time in units of rotation period. (a) Qualitative location of the transitions on the hysteresis loop. (b) Transitions from II→I. The time is initialized at the moment when the forcing was changed from a super-critical one to a sub-critical one. The normalized kinetic energy decays on a timescale τ obtained from an exponential fit (lines). (c) Transitions from I→II, starting from the same initial state.

  . Under the QG approximation, the experimental flow can be described by the classical 2D barotropic vorticity equation on the βtopographic β parameter resulting from the radial variations of the fluid height and α the linear Ekman friction parameter:

  x v -∂ y u the associated vorticity. The β-plane barotropic vorticity equation (4

Figure 9 .

 9 Figure 9. Comparison between the experimental flow and the linear QG model. (a-d ) Experimental flow measured 2 seconds after the third forcing was turned on. The experimental data, originally obtained in a cylindrical geometry, is remapped for the sake of comparison with our cartesian model. R3 is the radius of the third forcing ring. Only 8 vortices are visualized, but the third ring contains 18 vortices. (e-h) Solution of the linear model (equation (4.17)) at time t = 2 s. The model parameters (α, β, qm, k f and U ) are estimated from the experimental parameters. (a,e) Zonal velocity perturbation. (b,f ) Radial velocity perturbation (note that uρ should be compared with -v). (c,g) Vorticity. (d,h) Reynolds stresses divergence (zonal flow acceleration).

Figure 10 .

 10 Figure 10. Comparison of the Reynolds stresses divergence between the experiment and the linear model solution. (a) Experimental Reynolds stresses divergence (equation (4.8)). (b) Reynolds stresses divergence computed from the linear model (equation (4.19)). Both the experiment and the linear model show the generation of a prograde jet at the forcing location, flanked with two retrograde jets.

Figure 11 .

 11 Figure 11. (a) Streamfunction amplitude as a function of the background zonal flow U . Each symbol is a solution of the linear model with a different background zonal flow. The black line is the corresponding analytical amplitude in the case of a doubly-periodic forcing (equation (4.23)). (b) Reynolds stresses at y = 0 as a function of the background zonal flow U . Each symbol is a solution of the linear model with a different background zonal flow. The black line represents the best fit of the Lorentzian given by equation (4.27) where the amplitude Rm is a free parameter. In both panels, the vertical dashed black line shows the location of U = -c where c is the forced Rossby wave phase speed (equation (4.24)).

Figure 13 .

 13 Figure 13. Rossby waves excited by the forcing. The third forcing ring (R3 = 0.214 m) is turned on at its maximum power from time t = 0 to 75 tR. (a) Time evolution of the radial component of the velocity at a fixed radius ρ = 0.17 m. (b) Same data band-passed filtered between 0.06 and 0.12 Hz. The slope of the dashed line is the non Doppler shifted phase speed of the Rossby waves excited by the forcing (equation 4.24) c ≈ -6.3 mm s -1 . The slope of the full line is the Doppler shifted phase speed u φ φ +c. It increases in absolute value since the zonal flow increases when the forcing is turned on, and at the chosen radius the flow is retrograde, as represented on the top panel.

Figure 14 .

 14 Figure 14. Illustration of the jet's instability for an experiment in regime II (U f = 3.3 mm • s -1 ) (a) t = 2412 to 2436 tR. (b) t = 2486 to 2511 tR. (c) t = 2568 to 2580 tR. (d ) t = 2580 to 2593 tR. (e) t = 2593 to 2605 tR. This sequence is also available as supplementary movie 4.

Figure 15 .

 15 Figure15. Calibration of the experimental forcing. Left to right and top to bottom: rings C1 to C6. The total RMS velocity inside of a region of interest is plotted for each pump separately, and several fractions of the pump maximum power (dots). The lines are the result of a linear fit between the induced velocity and the pump power.

Figure 16 .

 16 Figure16. Sensitivity of the width of the bistable zone (left) and the zonal flow amplitude at the transition (right) to the model parameters. R1 and R2 are respectively the lower and upper limits of the bistable zone, in terms of forcing amplitude (see figure12). U1,max and U3,min are the zonal flow amplitude at the saddle-node bifurcation S2 and S1 respectively (figure12).(a) Varying friction coefficient α. The vertical dashed line shows the experimental friction. (b) Varying forcing wavenumber k f . The horizontal dashed lines show the mean flow amplitude at the transitions measured experimentally. (c) Varying β parameter. The vertical dashed line shows the experimental β.

Table 1 .

 1 Typical instantaneous RMS velocity (equation (3.1)), Rossby and Reynolds global and local numbers. R is the tank's inner radius R = 0.49 m, and f = 2Ω = 15.7 rad/s. For the local Rossby and Reynolds numbers, we use the distance between two forcing rings as a length scale, ℓ = 7.3 cm. Note that these values correspond to two typical experiments, but may vary in each regime depending on the forcing amplitude.

	2	671	100

cm/s 3 cm/s (a) (b)

Acknowledgments

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 and ingenuity during the conception and building of the experiment. We also thank J.-J. Lasserre for his help in setting up and calibrating the PIV system.

Declaration of interests

The authors report no conflict of interest.

Supplementary movies

Supplementary movies 1-4 are available at Appendix A. Quasi-geostrophic approximation of the experiment

We use the cylindrical coordinates (ρ,φ,z) with z oriented downward and (e ρ , e φ , e z ) the associated unit vectors (figure 1). We consider the flow of an incompressible fluid of constant kinematic viscosity ν and density ρ f , rotating around the vertical axis at a constant rate Ω = Ω e z . In our setup, Ω > 0 since the turntable rotates in the clockwise direction. We denote the velocity field u = (u ρ , u φ , u z ) eρ,e φ ,ez . The fluid is enclosed inside a cylinder of radius R. The lower boundary is a rigid plate located at z = 0 and the upper boundary is a free surface defined by z = -h(ρ). Note that here we assume that our experiment, which have a parabolic free-surface a curved bottom, can be modelled 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 (A 8)). For a bottom which is almost flat, we expect these effects to be of small amplitude, but 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.

Because of fast background rotation, or equivalently the small Rossby number of the system, the geostrophic balance dominates the experimental flow. As a consequence, the flow is quasi two-dimensional, but the curvature of the free-surface as well as the friction over the bottom (Ekman pumping) induce three-dimensional effects. Nevertheless, the weakness of these effects allows their incorporation into quasi-two-dimensional physical models, the so-called "quasi-geostrophic" models. In this section, we derive the conventional quasi-geostrophic model corresponding to our experimental setup, "conventional" meaning that we retain only the linear terms of these 3D effects. This reduced model allows us to:

-demonstrate that the free-surface curvature leads to a β-effect analogous to a linear variation of the Coriolis parameter with radius (i.e. to a β-plane);

-express the linear friction due to the Ekman pumping.

We start from the continuity and Navier-Stokes equations in the rotating frame and assume that, because the system is dominated by the geostrophic balance, the horizontal velocity field is independent of height (∂ z u ρ = ∂ z u φ = 0):

where

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 w. Indeed, as previously explained, w 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

where

The last term of the left hand side of equation (A 4), the vortex stretching term, involves the horizontal divergence of the flow which can be estimated from equation (A 3) after integration from z = -h(ρ) to z = 0 (z oriented downward) to unveil the Ekman pumping through the vertical velocity:

The vertical velocity at the free surface u z z=-h is given by the kinematic condition

since h is axisymmetric and we neglect any temporal variations of the fluid height (rigid lid approximation). The vertical velocity at the bottom u z z=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:

where δ = 2ν/f is the thickness of the Ekman layer and E = 2ν f h 2 0 is the Ekman number, h 0 being the mean fluid height. The horizontal divergence (A 6) is then

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. Substitution of the horizontal divergence (A 9) in the vorticity equation (A 4) yields

As stated before, we stand in the limit where the local Rossby number of the flow Ro = ζ/f is small, thus ζ ≪ f . Retaining only the linear part of the β-effect and Ekman pumping, we retrieve the classical 2D barotropic vorticity equation in the βplane approximation:

with β the topographic β parameter resulting from the free-surface radial variations and α the linear Ekman friction parameter:

Appendix B. Experimental forcing calibration

The forcing was calibrated in situ on the horizontal laser plane used for PIV measurements, while the system is in solid-body rotation. For each pump C i , we turn it on at a given fraction of its maximum power. We measure the corresponding velocity field, and define a region of interest (ROI) around the chosen ring, limited by two circles (rings C i-1 and C i+1 ). We measure the total RMS velocity on this ROI, 1 to 3 seconds after the forcing was turned on, i.e. when the forced vortices have reached their maximum vorticity but before the zonal jets fully develop. This measurement was realized for each ring separately and several pump powers. The corresponding data are represented in figure 15. We then performed a linear fit of the induced RMS velocity as a function of power to obtain a calibration law for each pump. In the main text, the forcing amplitude U f corresponds to the mean of the six RMS velocity deduced from our calibration, knowing the power fraction for each pump.

Appendix C. Sensitivity to the model parameters

Appendix D. Zonostrophy index of our experiments

In the context of eddy-driven jets in forced-dissipative experiments, two important length scales can be introduced. First, the Rhines scale L R [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF], which is the