Quenching of zonal winds in Jupiter’s interior

Significance Strong east-west directed jets dominate the circulation on Jupiter. Gravity and magnetic data show that they extend deep into the planet, but are quenched above the point where the rising electrical conductivity reaches the value of sea water. The lack of solid boundaries in gas planets requires a special mechanism for braking the jets. A simple model including a layer that is stable against radial motion explains all observational constraints when the jets are quenched near 2,000 km depth. We identify a single parameter, describing the influence of gravitational, rotational, and electromagnetic forces, which controls the results. These agree well with results from a full 3D convection model. Our model could also be applied to Saturn or to gaseous exoplanets.

Gravity and magnetic field data obtained by the Juno mission show that Jupiter's strong zonal winds extend a few thousand kilometers into the interior, but are quenched above the level where the electrical conductivity becomes significant.Here, we extend a simple linearized model [Christensen et al., Astrophys.J. 890, 61 (2020)] that explains the braking of the jets by the combination of stable stratification and electromagnetic effects.We show that in the inviscid limit, the process is essentially governed by a single parameter, which we call the MAC-number (for the forces acting on the flow-Magnetic, Archimedian, and Coriolis).The predictions for the drop-off of the zonal winds agree well with results from 3D-convection models.We run calculations that take the full range of density and electrical conductivity variations in the top 5,600 km of Jupiter into account.In order to satisfy constraints on the power driving the jets and on their effect on Jupiter's magnetic field, the top of the stable layer and the region where the jet velocity drops sharply must be near 2,000 km depth.The dissipation associated with quenching of the jets increases toward the poles, which can partly explain why the jets near ±20 • are faster than those at higher latitude.

Jupiter | interior dynamics | magnetohydrodynamics
Tracking of Jupiter's clouds has revealed that its banded surface pattern is associated with multiple vigorous jet streams, alternating with latitude between eastward and westward flow.Their velocity exceeds 100 m/s at low latitudes and is a few tens of m/s at higher latitudes.Measurements of Jupiter's gravity field by NASA's Juno mission suggest that the winds extend downward by some thousand kilometers (1).This interpretation did not go unchallenged.The nonuniqueness in the inversion of gravity anomalies allows, in principle, for an origin associated with dynamical processes in the interior that are unrelated to the surface winds (2).However, more recent analyses of the gravity signal strengthened the case for their relation to a deep continuation of the surface winds (3).They also support the assumption of the flow aligning with the planet's rotation axis rather than continuing in radial direction (4), as is expected on dynamical grounds.
The electrical conductivity in Jupiter is negligible near the surface but rises sharply with depth in the top several thousand kilometers.If high wind velocity persisted down to where the electrical conductivity becomes significant, the secular variation of Jupiter's magnetic field would be larger than observed (5) and the magnetic field morphology should be distorted by the wind shear (6).The observed secular variation has been modeled as advection by zonal flow with ≲1 cm/s (7), assumed to be unrelated to the surface winds.Strong effects associated with the surface winds would occur when a local magnetic Reynolds number (Rm) associated with their deep extension gets larger than order one at some depth (6).For a widely used conductivity model (8) Rm would exceed one near 2,500 km if the wind velocity had not started to drop above this depth.The combined constraints from the gravity and magnetic field data suggest that the jet velocity remains rather constant down to about 2,000 km but drops sharply below (9).
The jets are driven by an eddy momentum flux from small-scale convection to largescale zonal flow (10).Numerical models of convection in a rotating spherical shell (11,12) produce strong zonal winds when the boundaries are stress-free.Friction at a rigid lower boundary, which may serve as a proxy for magnetic forces impeding the flow, inhibits winds at higher latitudes, inside the so-called tangent cylinder (a cylinder aligned with the polar axis, touching the inner boundary at the equator) (13).Models which combine magnetic field generation in a deep conducting region with convective flow in a poorly conducting outer shell (14,15) indeed show strong zonal surface winds at low latitudes, but none at high latitudes.
An additional ingredient seems to be required to allow for strong winds inside the tangent cylinder.Stable stratification in the right depth range is a promising candidate.Demixing of hydrogen and helium in the gas planets is a possible cause for stable

Significance
Strong east-west directed jets dominate the circulation on Jupiter.Gravity and magnetic data show that they extend deep into the planet, but are quenched above the point where the rising electrical conductivity reaches the value of sea water.The lack of solid boundaries in gas planets requires a special mechanism for braking the jets.A simple model including a layer that is stable against radial motion explains all observational constraints when the jets are quenched near 2,000 km depth.We identify a single parameter, describing the influence of gravitational, rotational, and electromagnetic forces, which controls the results.These agree well with results from a full 3D convection model.Our model could also be applied to Saturn or to gaseous exoplanets.stratification (16), but it would only occur at depths where the jets have vanished.Stratification at shallower depth has been proposed in order to simultaneously match the gravity moments and the atmospheric element abundance (17)(18)(19), but the actual depth range and its physical origin remain uncertain.The existence of a stable radiative zone, originally suggested by Guillot and Gautier (20), has recently been supported by the finding of low abundances of alkali metals (a major source of opacity) in Jupiter's deep atmosphere (21).According to ref. 20, the radiative zone would exist in the temperature range 1,200 to 3,000 K, corresponding to about 500 to 2,000 km depth, which may be too shallow for playing a role in the truncation of the winds.We base our analysis on the hypothesis of the existence of stable stratification in the depth range where the wind velocity drops off, even though independent evidence for it remains scanty.
The dynamical role of a stable layer, in combination with electromagnetic forces, for the quenching of winds has been demonstrated in simple axisymmetric models with an imposed dipolar field where the winds are driven by an ad hoc force (22).A weak meridional circulation that is associated with the east-west winds (23), similar to Ferrel cells in the Earth's atmosphere, perturbs the density stratification in the stable layer, turning it baroclinic.In accord with a thermal wind balance the zonal velocity decreases with depth.These models also showed that the influence of electromagnetic forces, acting on the meridional flow, makes the drop-off of the zonal winds very sharp near the top of the stable region.
Some recent 3D models that cover the full range from the deep dynamo region to a nearly insulating outer convection zone included a stably stratified layer sandwiched between the two (24,25).They confirmed that such a layer allows for strong surface winds at higher latitudes which are quenched inside the stable region.These simulations are computationally very demanding.Simpler 3D models that involve only an upper convective region and a stable layer, with a magnetic field of strength B imposed, are more amenable for a systematic study, with the focus on the location and sharpness of the drop-off of the winds.Purely hydrodynamic simulations (26) showed that the penetration distance decreases with the degree of stability.Simulations with exponentially varying conductivity demonstrated that, for a fixed degree of stability, the combination B2 controls how steeply the zonal wind velocity drops in the stable layer (27).
All 3D convection simulations suffer from the fact that they are overly viscous and some control parameters differ vastly from their planetary values.In ref. 22, we developed linearized 2D Cartesian box models of jet flow in a conducting stable layer.The simplifications allowed us to reach Jupiter-like parameters.Here, we build on this concept and extend it to the inviscid limit.We isolate a single parameter combination that controls the decay of zonal winds in a stably stratified conducting region.While the original model was restricted to polar latitudes, where gravity, rotation, and dipolar magnetic field nearly align, we now extend it to arbitrary angles between the different vectors.We also run the 3D simulations of (27) at lower viscosity in order to corroborate some of the scaling laws based on the simple model.Finally, we apply the Cartesian models to Jupiter, covering the full range of density and conductivity variations from the 1-bar level to 5,600 km depth.

Linearized Cartesian Models in the Inviscid Limit
Basic Concept and Equations.Building on the concept introduced in ref. 22, we elaborate on a simple Cartesian model for the damping of zonal flow in a stably stratified region including the action of magnetohydrodynamic forces.Here, we first simplify it further by dropping viscosity and other terms in the equations that turned out to be insignificant.This reduces the problem to one that is governed by a single nondimensional control parameter.
We start out from the anelastic equations for magnetohydrodynamics (MHD) flow in a rotating system including buoyancy forces, the electromagnetic induction equation, and an advectiondiffusion equation for codensity (e.g., ref. 28).
Here, u is the velocity vector, B magnetic field, ˜ background density, J = −1 o ∇ × B electrical current density, P an effective pressure, and C the codensity (28), i.e., the (negative) density anomaly of thermal or compositional origin normalized by ˜ .Ω is the rotation frequency, F viscous force (neglected below), g gravity, thermo-compositional diffusivity, and o magnetic permeability.All are assumed to be constant, except for and ˜ , which may vary with depth.
A torus in the outer layers of the planet is approximated by a rectangular Cartesian box (Fig. 1), where z is antiparallel to gravity and x corresponds to latitude.No variation in the perpendicular y-direction, representing longitude, is assumed.The background codensity is set to zero (neutral stability) in the upper part of the box at z > 0, whereas a constant gradient in codensity is assumed in the stable lower region z < 0. The Brunt-Väisälä frequency N = (g dC/dz) 1/2 characterizes the degree of stability.We impose a uniform magnetic field B in the x-z-plane.We solve for the codensity perturbation c, assumed to be small compared to the background codensity C. We drop all nonlinear terms including momentum advection and assume a steady state.In addition to these simplifications already made in ref. 22, we drop the viscous term in Eq. 1 and neglect perturbations of the poloidal magnetic field, which is equivalent to setting J y = 0.Both had been found in ref. 22 to be negligible for planetary Fig. 1.Schematic view for the setup of the linearized models.A Cartesian rectangle is cut out from a cross-section of the outer spherical shell region of a gas planet.Its lower part (shaded) is stably stratified.The z-coordinate is antiparallel to gravity g, x points north and y points west.In general, the orientation of the rotation vector Ω and the magnetic field vector B are at oblique angles, and , respectively, relative to gravity.All three vectors are assumed constant within the rectangle.
parameter values.Because of the linearity, we can expand all variables in Fourier series in the x-direction.For each horizontal wavenumber k we obtain a set of coupled ordinary differential equations.
While the cartoon in Fig. 1 shows a more general setting that we treat later, we start out by considering the special case where g, B, and are parallel, like in ref. 22. Furthermore, we begin with the incompressible (Boussinesq) limit, i.e., ˜ = o .We assume an exponential variation of electrical conductivity, (z) = o exp(−z/d ), with o the value at the upper boundary of the stable region at z = 0, and d the conductivity scale height.This, together with the other approximations and by scaling the variables in a suitable way, allows us to reduce the system of equations to a very compact set.We use the following dependent variables, U = u y for the zonal flow, the stream function for the meridional circulation with u = ∇ × e y , b = B y for the induced toroidal magnetic field, and the codensity perturbation c.Nondimensional variables, indicated by a tilde, are obtained by scaling length with d and by converting the dependent variables according to the following scheme (which differs from the scaling used in ref. 22): We expand U , , and b in terms of sine functions, e.g.U = U (z) sin(kx), and c in cosine functions.For a given wavenumber k we obtain the following nondimensional equations, where we drop the tilde for simplicity and use shorthand d z and d zz for the first and second derivative, respectively: Eq. 7 derives from the y-component of the Navier-Stokes Eq. 1 and shows that in regions where the electrical conductivity is sufficient the meridional circulation is controlled by electromagnetic forces resulting from the induced toroidal field (assuming that here no other significant force is acting on it).Note that in our model, we do not consider the region where the zonal flows are driven, which could be predominantly in the top atmospheric layers of the planet (29).Here, the Coriolis force of the meridional circulation balances the driving Reynolds stresses (or eddy momentum flux convergence).Eq. 7 holds in this form only where Reynolds stresses are absent.The electromagnetic force in Eq. 7 can be considered to be the antagonist to the Reynolds stresses; its role is to truncate the meridional flow at depth.It has been suggested that the breaking of inertia-gravity waves that are generated in the weather layer can serve as the antagonist ( 30), but it is not clear why the downward propagating waves should break.Eq. 8 describes a thermal wind balance and derives from the y-component of the curl of Eq. 1. Eq. 9 represents the balance of advection of codensity by the meridional flow and diffusion.Eq. 10 describes the balance between the induction of toroidal magnetic field, by the shearing of the imposed poloidal field due to gradients in the zonal flow (Ω-effect), and magnetic diffusion.
For a fixed wavenumber k the equations contain only a single control parameter which we call the MAC-number to express that the system is influenced by Magnetic, Archimedian (=buoyancy), and Coriolis forces: We had found in 3D numerical simulations (27) that magnetic effects influence the decay of zonal winds in the stable layer by the parameter combination B 2 .This concurs with this combination appearing in the MAC-number.In ref. 22, we found that, for fixed and B, a parameter Ra Ω ∝ N 2 / controls the decay of the zonal flow when viscosity becomes insignificant.With the MACnumber we now combine all physical properties that influence the quenching of the zonal flow into a single parameter.
We assume that at z > 0 turbulent convection eradicates significant codensity differences, and we set c(z > 0) = 0.In ref. 22, we specified an ad hoc force in the y-direction, acting at z > 0, for driving the zonal flow.Here, we simply set the wind velocity by a boundary condition U (z = 0) = 1.U =1 also holds for z > 0, where c vanishes, according to Eq. 8.For other applied boundary conditions and the technique for solving Eqs.

7-10, see Materials and Methods.
Results for the Simple Case.We calculated results for MACnumbers from 0.05 to 500,000.Values of order one are typical for 3D numerical simulations whereas in Jupiter MAC falls into the range 10 2 to 10 5 when the upper boundary of a stable layer is between 1,500 km and 3,000 km depth and N /Ω is of order one.We fix the horizontal wavenumber at k = /20, which is a representative value for Jupiter.Fig. 2 shows depth profiles of the solution for MAC = 500.The wind velocity (Fig. 2A) drops strongly within a few conductivity scale heights below the stable layer boundary.Within a limited depth, interval there is a weak reverse flow.Deeper than six scale heights below the boundary the velocity is less than 10 −4 .The depth range over which the velocity drops depends on the MAC-number.The thin blue line in Fig. 2A represents the case for MAC = 0.1.The overshoot of the zero-line is more pronounced for the higher MAC-number, but otherwise, the general shape of the drop-off is similar.Also the meridional flow and the codensity perturbation drop to very small values at a few scale heights below the stable layer boundary (Fig. 2B).However, the toroidal magnetic field, which is induced by the velocity shear in the upper part of the A B stable layer, stays constant at greater depth in the region of rapidly increasing conductivity.Because for constant b(z) the LHS of Eq. 10 is constant, U (z) must decrease exponentially here at the same rate as the conductivity increases.To a limited degree, the toroidal field diffuses upward where it causes the meridional (vertical) flow to drop somewhat before it reaches the stable layer boundary.
To characterize the length scale of quenching of the zonal flow in the stable layer we determined the range over which U has dropped by a factor ten and call this the decay scale d 0.1 .It decreases from about six conductivity scale heights at MAC = 0.05 to 0.3 scale heights at MAC = 500,000 (Fig. 3A), although not by a simple power-law.A good empirical fit, with a mean relative error of 2.4% for MAC ≥ 1 is obtained by Changing the horizontal wavenumber k has only a moderate influence on d 0.1 (triangles in Fig. 3).The power needed for driving the jets must stay within reasonable limits.The total power balances the dissipation which is purely ohmic in our model, plus the work done by the meridional flow against buoyancy forces (see SI Appendix for their calculation).Fig. 3B shows the variation of driving power per unit surface area, scaled via as function of the MAC-number.The nondimensional power drops with MAC, but varies only little at MAC-numbers larger than 1,000, leveling off toward a constant value near 0.5.A good fit with 2.5% mean relative error is given by Comparison with 3D Convection Simulations.We test the prediction that the quenching of zonal flows in a stable layer is governed by the parameter combination in the MAC-number with 3-D convection models for a rotating spherical shell.We briefly summarize the basic concept described in ref. 27.We consider a shell with ratio of inner to outer radius r i /r o = 0.7.
The region above 0.83r o is convectively unstable, whereas at r < 0. In ref. 27, we kept the Ekman number at 10 −5 and the degree of stability in the lower layer was set such that N /Ω = 4.9.Here, we add two more cases.We lower the Ekman number to 2×10 −6 to weaken the influence of viscosity and to increase MAC by a factor of five.The same value of the MAC-number is obtained by setting N /Ω = 10.9,E = 10 −5 .We keep the Prandtl numbers at Pr = Pm = 0.5.The strength of the imposed field is set such that the local Elsasser number i B 2 (r i )/( i Ω), evaluated at the inner boundary and at the poles, is 0.25.The Rayleigh number is chosen such that in each of the simulations the convective Rossby number Ro c = E(Ra/Pr) 1/2 is kept at a value of 0.346.Given that we model strongly supercritical convection in the rapidly rotating regime, this ensures that the vigor of convection is comparable in the different cases.
Fig. 4 illustrates the flow pattern for the case with E = 10 −5 .At the surface, it is dominated by bands of zonal wind that extend parallel to the rotation axis throughout the outer convection shell.Inside the tangent cylinder, associated with the boundary of the stable layer, the wind velocity drops rapidly with depth in the lower stable region.
To compare the results of the 3D convection models with the predictions of our simplified theory we calculate the decay scale d 0.1 as follows.As in ref. 27, we average the azimuthal velocities in the center of all jets inside the tangent cylinder, normalized with the velocity at the surface, to obtain a mean radial profile ū .The relevant value of the MAC-number expressed in terms of the parameters of the 3D-simulations is obtained as where r s is the effective boundary of the stable layer.In the simple model this boundary is sharp, whereas in the 3D-models it is spread over a transition region.We choose r s = 0.82r o , where the stable codensity gradient has reached 36% of its saturated value.The values of the MAC-number in the 3Dmodels are of order one.From the averaged profiles ū we read the decay scale d 0.1 as the distance in radius over which the velocity has dropped to 10% of its value at r s .The results, normalized with d , are shown by the stars for the cases with N /Ω = 4.9 and by the circle for N /Ω = 10.9 in Fig. 3A.They agree well with the predictions from the simple inviscid and linearized Cartesian models.Although in the 3D models viscosity still plays a significant role, it hardly affects the decay scale of the winds, whereas the additional source of dissipation reduces their amplitude (SI Appendix).

Models Adapted to Jupiter
In order to make the linearized Cartesian model directly applicable to gas planets we extend it in two respects.We now study the compressible case instead of the Boussinesq limit, and we allow for arbitrary functions for the reference density ˜ (z) and electrical conductivity (z).Furthermore, we allow for an angle between rotation vector and the z-direction (the direction of gravity) and angle between magnetic field and z, to account for conditions at different latitudes (Fig. 1).The resulting equations are given in Materials and Methods.
We model the outer 5,600 km of Jupiter from the 1-bar level down to 0.92r J (r J ≈ 70,000 km is Jupiter's mean radius), and study different radii r s for the upper boundary of the stable layer.We employ analytical approximations for the electrical conductivity data in refs.8 and 22 and for the density model of (31) (SI Appendix).For the characteristic latitudinal wavelength of the zonal jets inside the tangent cylinder, we pick 10,000 km.The angle between magnetic field and the vertical at colatitude is calculated for a polar dipole via cot = −2 tan .The field strength is obtained as B = g 0 1 (r J /r s ) 3 (4 cos 2 +sin 2 ) 1/2 , with the Gauss coefficient g 0 1 = 0.41 mT.Jupiter's rotation frequency is Ω = 1.75×10 −4 s −1 .We fix the degree of stability to N /Ω = 1 and set = 1.5 × 10 −3 m 2 /s (22).For the velocity amplitude of the jets inside the tangent cylinder we take U o = 25 m/s, although closer to the equator it is larger.As the jet velocity is set to a nondimensional value of one in the calculations, the scaling factors in Eqs. 5 and 6 must be multiplied with U o d / to convert the results into physical units.The power varies quadratically with U o , hence its scaling factor (Eq. 13) is multiplied by (U o d / ) 2 .Note that in the modified scaling factors the (effective) diffusivity drops out.Fig. 5 illustrates some properties for the case with = 45 • and r s = 0.975r J in the vicinity of the stable layer boundary.The MAC-number calculated with local properties ( , B, d , and ˜ ) at r s is 147.Above r s (or z = 0) the wind velocity and also the meridional flow are aligned with the direction of rotation.U drops sharply within two conductivity scale heights below r s .The induced toroidal magnetic is aligned with the direction of the imposed poloidal field at greater depth.In terms of absolute numbers, the meridional flow is weaker than the zonal flow by at least five orders of magnitude and the toroidal field is weaker than the poloidal field by one order of magnitude.
For jets at mid-latitude ( = 45 • ) several properties of interest have been calculated for different depths of the stable layer boundary in the range from 1,540 km to 2,800 km.The electrical conductivity o rises from 10 the jet velocity drops to 10% of its surface value decreases from 400 km to 100 km (Fig. 6A).The power per unit area for driving a jet with 25 m/s velocity amplitude increases steeply with the depth of the stable layer boundary (Fig. 6B), as was found before (22).For a boundary shallower than 2,000 km, it is less than 1 W/m2 , whereas for a depth greater than 2,400 km it exceeds the mean interior heat flow from Jupiter of 7.5 W/m 2 (33).The broken lines in Fig. 6 A and B show the values calculated from Eqs. 12 and 14.Even though they are derived from a simpler Boussinesq model with B g, they fit the results of the Jupiterlike models well.Hence the MAC-number controls the essential physics of the quenching of winds in a weakly conducting stable layer also in the more realistic case.It suggests that Eqs. 12 and 14 can serve as a guide to assess the influence of changes in any of the physical parameters that play a role in the process.
The poloidal magnetic field will be modified by the zonal flow when the magnetic Reynolds number calculated with local values, Rm loc = o d U , reaches or exceeds order one at some depth.Fig. 7 shows how Rm loc varies with depth for two cases.It peaks right below the stable layer boundary and remains at roughly one tenth of the peak value throughout the stable layer.Rm loc rises strongly with increasing depth of the stable layer boundary in pace with the associated increase of electrical conductivity (Fig. 6C).Its peak value exceeds one when the boundary is deeper than 2,400 km.
We studied the dependence of the solution on latitude in the range from = 0 • to 70 • , the latter value being close to the edge of the tangent cylinder, but still within.The decay scale d 0.1 is almost insensitive to a change of latitude.As a consequence, also the peak value of the local magnetic Reynolds number is unaffected.A strong effect is found for the power requirement of the zonal flow, which becomes smaller at lower latitudes (Fig. 8).The power P varies almost as cos 2 , as shown by the broken lines.The amplitude of the meridional flow and of the induced toroidal field b also vary with colatitude, approximately as cos .Their covariation with P is expected, since the meridional flow amplitude is controlled by the driving Reynolds stresses, and and b are directly linked by Eqs. 7 and 16.
The strong dependence of P on latitude can partly be traced back to the variation of the dipole field strength with , because the scaling factor for power (Eq.13) contains B 2 .However, from the pole to = 70 • the power is found to drop by a factor of eight, whereas B 2 varies only by a factor of three.Therefore, the  change of orientation between gravity, rotation, and magnetic field vectors also has an influence on the dissipation.Simple (Boussinesq) test cases in which and are varied independently show that the angle between gravity and the magnetic field is key.Setting = 0 • , = 90 • results in a reduction of power by a factor of eight compared to = = 0 • .Increasing while keeping = 0 • has only a minor effect.
The vertical (meridional) velocity is only a small fraction of a mm/s in the stable layer.The meridional mass flux, or the product ˜ , stays constant with z according to Eq. 16 at shallow depth where b is essentially zero.Because of the sharp drop of density with radius in the outer regions of the planet, the velocity increases toward the surface.The constancy of mass flux holds as long as the driving force for the zonal winds, i.e., the Reynolds stress, remains negligible.If we assume that the winds are mostly driven by vigorous moist convection near the top of the atmosphere, down to a pressure level of 6 bars (57 km depth) (29), we can calculate the vertical velocity at this depth assuming constant ˜ below.For a stable layer boundary at 2,100 km, it is 0.5 mm/s.This gets larger (smaller) by an order of magnitude for 2,550 km (1,650 km) boundary depth.These estimates for the vertical velocity are an upper limit if Reynolds stresses in deeper parts of the convection zone contribute substantially to the driving of the jets.

Discussion and Conclusions
The MAC-number proves to be the primary control parameter for the quenching of zonal winds and the power required to drive them.We derived its influence for a simple setting, i.e., using linearized equations in the Boussinesq approximation with aligned gravity, rotation axis, and magnetic field.A comparison of the predictions of this simple model with the results of 3-D simulations of convection-driven jets, and with linearized simulations for a density-stratified fluid with different angles between the relevant vectors, shows that the MAC-number is indeed the essential parameter also in less idealized cases.We can use it to assess the consequences of variations in properties that are not well known.Double-diffusive convection in the stable layer, probably in the form of stacked sublayers, is needed to transport the internal heat of the planet without excessive temperature gradients (34).This can be parameterized by an effective diffusivity , for which we adopted the value estimated in ref. 22.However, this value is merely an educated guess.An increase or decrease of by an order of magnitude alters MAC by a factor of 10.In the relevant range for the MAC-number this implies a change of the wind decay scale and of the driving power by a factor of 1.7 at most.The same holds for variations in the stabilizing density gradient.Therefore the result is not strongly sensitive to the adopted values of these poorly constrained parameters.However, for a weakly stable stratification with N /Ω < 1 the shear flow in the stable layer may be subject to Kelvin-Helmholtz instabilities (SI Appendix).
The characteristic drop-off length scale for the winds depends moderately on the depth of the stable layer boundary.It lies in the range of 100 to 300 km (Fig. 6A), i.e., it is rather sharp and substantially smaller than the depth of the overlying region where the wind velocity is constant.In the best-fitting solutions to the observed gravity anomalies the drop-off scale is typically larger than this (e.g., extended data figure 2 in ref. 4), although when additional magnetic constraints are included it becomes sharper.The discontinuous onset of stability in our models favors a rather sharp drop-off.In the case of a gradual increase of stability with depth the decay of the jet velocity can be expected to be more smeared out.
The power required for driving the wind and the peak value of the magnetic Reynolds number associated with the winds depend strongly on the location of the stable layer.This allows us to put constraints on the depth extent of the winds in addition to those based on the gravity signal.For the dissipation (power) the internal heat flow from Jupiter is frequently considered to be an upper limit.On thermodynamic grounds, the total dissipation could be up to 3.6 times larger than the internal heat flow (32).This comprises also the dissipation associated with turbulent convection and with the dynamo process at greater depth.Therefore it appears unlikely that the dissipation linked to the zonal winds substantially exceeds the internal heat flow.The latter would be the case if the stable layer boundary was deeper than 2,500 km (Fig. 6B).This limit could be shifted somewhat, but not more than a few hundred kilometers, if the effective diffusivity or the degree of stability differed by an order of magnitude from the adopted values.Rather strong distortions of the nonaxisymmetric components of Jupiter's magnetic field by the winds must be expected if the local magnetic Reynolds number of the wind significantly exceeded one at some depth (6).Such distortions are not observed.Taking Rm loc = 3 as an upper limit requires the stable layer boundary to be no deeper than 2,600 km, which is almost identical to the constraint based on the driving power.
A shallow boundary is not in conflict with constraints from the dissipation or magnetic field properties.Explaining the observed gravity signal with a truncation depth at 1,500 km would require wind velocities that are twice larger than those observed at the cloud level (3).A similar increase of velocity from the surface to 50 km depth has been found by the Galileo probe at low latitudes (35).Should this be the case also at higher latitudes, a stable layer reaching up to 1,500 km would be possible.The dissipation associated with the winds would then be only a small fraction of the internal heat flow.
The flux of energy Q from small-scale eddies into the mean flow, per unit mass of the fluid involved, has been calculated from the observed eddy momentum flux u x u y at the cloud level (10,36), where u denotes the horizontal components of the eddy flow and the angular bracket stands for the spatial average.Here, we adopt the range of Q = (7−12) × 10 −5 W/kg (10).A similar value of has been derived with a diagnostic based on potential vorticity (37).To arrive at the power driving the zonal flow an estimate of the mass involved in the eddy motion is needed.If u x u y stays at the surface value, from a reference pressure p o down to a pressure level p, and drops to zero below, the total power is given by P = Q(p − p o )/g, where g = 25 m/s 2 is Jupiter's gravity.Setting p o = 1 bar, the pressure at the bottom of the region in which the jets are driven would be at 1.3 bar, 4 bar, and 31 bar for a power of 0.1, 1, and 10 W/m 2 , respectively, corresponding to depths of 7 km, 42 km, and 150 km.The first value is substantially smaller than the thickness of the weather layer.It seems unlikely that deeper parts of the zone of moist convection do not contribute to the eddy momentum flux, which makes a very shallow stable layer boundary, requiring little power, less likely.The middle value fits roughly the zone of moist convection.For the highest value, the zone of intense driving of the mean flow had to extend deeper than the weather layer.In addition to the other arguments this disfavors a deep stable layer boundary, provided that Reynolds stresses associated with convection below the weather layer are not important.With this reservation, the most probable depth to the stable region and the point where the zonal flows are quenched is therefore around 2,000 km.
All previous conclusions regarding the depth of the zonal winds depend on the applicability of the conductivity profile that we use (8).A more general conclusion would be that the stable layer boundary is likely located at the depth where the electrical conductivity reaches approximately 0.01 to 0.03 S/m (typical conductivity of freshwater).For example, with the conductivity model of ref. 38, the preferred depth would shift from 2,000 km to 2,400 km.
A deep meridional flow is an essential part of our model.The presence of meridional flow in Jupiter, down to a level of at least 240 bars and organized in a set of Ferrel-like cells associated with the zonal winds, has been inferred from anomalies in the ammonia concentration observed by the microwave instrument on Juno (23).These anomalies have been attributed to vertical advection of ammonia concentration gradients.Our models with a stable layer boundary near 2,000 km predict a radial velocity of the meridional flow at the 6-bar level of 0.5 mm/s.An interesting question is whether such slow flow is sufficient to cause the observed ammonia anomalies.
The significant decrease with colatitude of the dissipation associated with a jet of given velocity (Fig. 8) can be part of the reason why the low-latitude jets inside the tangent cylinder are faster than those at mid-to high-latitudes.This assumes that there is little transfer of energy between different latitudinal bands, which seems likely because the winds are driven by a local eddy momentum flux and extend along geostrophic cylinders into the region where their energy is dissipated.However, reduced dissipation is unlikely to be the sole reason for faster jets at low latitudes.The analysis from cloud-tracking shows a good correlation of the eddy momentum flux with the velocity of the jets (10).Furthermore, in the concept of zonostrophic turbulence (e.g., ref. 39) the rate of transfer depends on the parameter ∝ 2Ω sin , which decreases with latitude.Still, the latitudinal variation of dissipation may contribute to the vigor of the jets near ±20 • latitude.
The existence of a stable layer at the required depth in Jupiter is hypothetical.However, we are not aware of a suitable mechanism without the stable layer for quenching the zonal winds at a depth that is consistent with the gravity and magnetic field data, and that is compatible with reasonable limits on the driving power.Direct quenching of the zonal flows by electromagnetic forces in the absence of stable stratification was found to entirely suppress jets inside the tangent cylinder (14,15,22).In order to consolidate or falsify our model, independent evidence on the existence, depth, and degree of stratification of such layer would be highly desirable.Such evidence might come from observations of normal modes of Jupiter (40).The framework that we developed here can also be applied to other objects with similar zonal surface winds, such as Saturn and possibly gaseous exoplanets.

Materials and Methods
Anelastic Model for Arbitrary Angles Been z,  and B. Using complex notation for all variables and allowing for arbitrary functions (z) and ˜ (z) Eqs.7-10 can be generalized, where now the meridional velocity relates to the stream function by ˜ u = ∇ × ( ˜ e y ): Note that the scale height d is negative as defined here.The MAC-number (Eq.11) is calculated with the reference values of B, ˜ , , and d taken at the upper boundary of the stable layer.
Boundary Conditions.For the Jupiter-like model we map the region 0.92r J < r < r J onto the nondimensional range z l < z < z u with z = 0 at the stable layer boundary.In the simple Boussinesq models, we picked |z l | and z u large enough so that their choice did not affect the solution.We deal with a sixth-order boundary value problem and must specify six boundary conditions.These are U(0) = 1, c(0) = 0, U(z l ) = d z (z l ) = c(z l ) = 0, and b(z u ) = 0.
Solving the Linearized Equations.The boundary value problem is stiff, in particular for large values of the MAC-number and/or steep variation of conductivity and density covering many orders of magnitude.It requires special measures for a stable solution.We employ quad-accuracy for floating point numbers and a multiple shooting method (41) based on fourth-order Runge-Kutta integration with up to ten internal matching points in the stable layer.42).This study used the adapted version ( 43) from ( 27).
83r o the fluid is stably stratified.Full stability is reached below a transition region at 0.804r o .The conductivity varies exponentially with radius by a factor 10 8 .An axial dipolar magnetic field is imposed by a boundary condition at r i .We assume fourfold symmetry in longitude, i.e., we simulate only one quarter of the full sphere, which has little influence on the character of the zonal jets and saves computing time.The relevant MHD-equations in the Boussinesq form (27) contain four basic control parameters, a Rayleigh number Ra, Ekman number E, Prandtl number Pr and magnetic Prandtl number Pm (see Materials and Methods for definitions).

Fig. 5 .
Fig. 5. Solution for a Jupiter-like case with the stable layer boundary at 0.975r J for = 45 • , = −26.6 • .Isocontours of wind velocity U (Left), meridional mass flux ˜ (Middle), toroidal magnetic field b (Right).Note that the current associated with the toroidal field follows the contour lines of b.The vertical coordinate is exaggerated threefold compared to the horizontal one.The z-coordinate is in multiples of the conductivity scale height at the boundary, 158 km for 0.975r J .

Fig. 6 .
Fig.5illustrates some properties for the case with = 45 • and r s = 0.975r J in the vicinity of the stable layer boundary.The MAC-number calculated with local properties ( , B, d , and ˜ ) at r s is 147.Above r s (or z = 0) the wind velocity and also the meridional flow are aligned with the direction of rotation.U drops sharply within two conductivity scale heights below r s .The induced toroidal magnetic is aligned with the direction of the imposed poloidal field at greater depth.In terms of absolute numbers, the meridional flow is weaker than the zonal flow by at least five orders of magnitude and the toroidal field is weaker than the poloidal field by one order of magnitude.For jets at mid-latitude ( = 45 • ) several properties of interest have been calculated for different depths of the stable layer boundary in the range from 1,540 km to 2,800 km.The electrical conductivity o rises from 10 −3 S/m to 16.7 S/m between 1,540 km and 2,800 km depth; the corresponding scale heights d are 145 km and 280 km, respectively.The distance d 0.1 over which

Fig. 7 .
Fig. 7. Local magnetic Reynolds number as function of depth for two Jupiterlike cases at mid-latitude with surface velocity amplitude U = 25 m/s and stable layer boundary at 2,660 km (full line) and 2,240 km depth (dotted).

Fig. 8 .
Fig. 8. Power as function of colatitude for two Jupiter-like cases with surface velocity amplitude U = 25 m/s and stable layer boundary at 1,960 km (circles) and 2,450 km depth (squares).The broken lines represent a variation with cos 2 , anchored at the value for = 0. Thin horizontal lines as in Fig. 6B.

2 [
cos d z + sin ik]( ˜ ) = (cos d z + ik sin )b,[16] (cosd z + sin ik)U = MAC ikc, [17] (d zz + d −1 d z − k 2 )c = ik , [18] (d zz − d −1 d z − k 2 )b = − (cos d z + sin ik)U.[19] 3D Simulations.The Rayleigh number, Ekman number, Prandtl number, and magnetic Prandtl number are defined as where is thermal expansivity, d C/dr the imposed background codensity gradient in the unstable layer, d the shell thickness, and o the conductivity at r i .Simulations were conducted with the code MAGIC (https://magic-sph.github.io/).See ref. 27 for details of the simulations.All 3D simulations were carried out using the magnetohydrodynamic code MagIC which is open source and available in Zenodo ( Data, Materials, and Software Availability.