The Effects of a Stably Stratified Region With Radially Varying Electrical Conductivity on the Formation of Zonal Winds on Gas Planets

The outer areas of Jupiter and Saturn have multiple zonal winds, reaching the high latitudes, that penetrate deep into the planets' interiors, as suggested by gravity measurements. These characteristics are replicable in numerical simulations by including both a shallow stably stratified layer, below a convecting envelope, and increasing electrical conductivity. A dipolar magnetic field, assumed to be generated by a dynamo below our model, is imposed. We find that the winds' depth into the stratified layer depends on the local product of the squared magnetic field strength and electrical conductivity. The key for the drop‐off of the zonal winds is a meridional circulation which perturbs the density structure in the stable layer. In the stable region its dynamics is governed by a balance between Coriolis and electromagnetic forces. Our models suggest that a stable layer extending into weakly conducting regions could account for the observed deep zonal wind structures.


Introduction
Zonal winds are alternately westwards/eastwards flows and feature across all four outer planets in our solar system.Those observed on the gas giants, Jupiter and Saturn, share some key characteristics.The dominating equatorial prograde flow on Jupiter (Saturn) spans roughly 30 • (60 • ) with an amplitude of around 100 m/s (400 m/s).This is flanked by a pair of slightly weaker retrograde jets and multiple jets reaching the highlatitude regions (Tollefson et al., 2017;García-Melendo et al., 2011).While these winds are weaker, they are still significantly stronger in amplitude than non-zonal flows.Jupiter's northern hemisphere also features an unusual prograde jet, as strong as the equatorial jet, at around 21 • latitude, introducing a strong equatorial antisymmetry into the dynamics.
Surface measurements, first from Voyager 1 and 2 (Ingersoll et al., 1981) then from Cassini (Salyk et al., 2006), have shown a strong correlation between the eddy momentum flux (or Reynolds stresses) and the zonal wind speed as a function of latitude.This confirms current theories that Reynolds stresses, which are statistical correlations of the components of the flow at small and intermediate scales, drive the zonal winds.
The extent of the winds into the jovian interior has recently been constrained using the gravity moment measurements from Juno, yielding a depth between 2, 500−3, 000km; around 96% of the planet's radius (Kaspi et al., 2018;Dietrich et al., 2021;Galanti et al., 2021).The same investigation has also been carried out for Saturn, using the Cassini measurements, suggesting the winds extend to 8, 000−9, 000km depth, around 85% of the planetary radius (Galanti et al., 2019).This is consistent with simulation-based studies where it has been found that the location of the flanking retrograde jets is usually coincident with the 'tangent cylinder', here loosely defined as the cylinder aligned with the axis of rotation with a radius corresponding to the depth at which jet quenching takes place.This has been found in numerical models studying both magnetic effects as a potential braking mechanism for the winds, with increasing electrical conductivity at depth (eg.Duarte et al. (2013)) or transition into a stably stratified region (Wulff et al., 2022).Therefore, based on these basic geometric observations of the dynamics we would expect the winds to penetrate deeper on Saturn, with its much wider equatorial jet.
A strong prograde jet flanked by two retrograde jets in the equatorial region, outside the tangent cylinder, were already reproduced in hydrodynamic simulations (Christensen, 2002;Heimpel et al., 2005;Gastine et al., 2014).However, simulations with rigid lower boundary conditions did not exhibit any zonal winds inside the tangent cylinder.Models with stress-free inner boundaries featured some high-latitude jets but failed to provide any insights into the winds' damping mechanism in the interior.
In both planets the increasing electrical conductivity at depth (e.g.French et al. (2012)), plays a crucial role in the zonal winds' downward propagation from the surface.
It has been speculated that deeply penetrating zonal winds may cause the observed secular variation (Moore et al., 2019).However, Bloxham et al. (2022) argue that a slight correction of Jupiter's rotation rate provides a better explanation, in combination with deeper flows in the dynamo region.Furthermore, considering reasonable limits for the total ohmic dissipation suggests that the winds may not penetrate into the highly conducting region of Jupiter (Liu et al., 2008;Wicht et al., 2019;Cao & Stevenson, 2017).
It was originally surmised that Lorentz forces, acting where the deep zonal flows reach the conducting region, were responsible for the braking of the winds.However, simulationbased studies such as Dietrich and Jones (2018) found that these Maxwell stresses at depth eradicate all large scale zonal flow above the conducting region, leading to zonal wind profiles with the strong flows confined to near the equator.Christensen et al. (2020) suggested that a combination of a stably stratified layer (SSL) and the magnetic effects at depth are responsible for the breaking of the zonal flows on Jupiter.They suggest that the winds decrease in the stable layer in accord with a thermal wind balance.The required density perturbation is caused by a meridional circulation which is affected by electromagnetic forces.Duer et al. (2021) present observational evidence for the existence of meridional flow associated with the winds.Gastine and Wicht (2021) conducted a global dynamo simulation with a strong radial variation of conductivity, which was successful in producing winds formed and being maintained above the highly electrically conducting region.Recently, Moore et al. (2022) also showed that dynamo simulations of Jupiter including a SSL at 90 − 95% radius produced dynamos with a dominant axial dipole component and a similar degree of complexity as the measured Jovian magnetic field.
In the context of Saturn a stable layer, shallower than the region of metallic conductivity, could help to explain both the formation of its high-latitude zonal winds and how they are quenched at depth, and its magnetic field.This is remarkably axisymmetric (Dougherty et al., 2018) and a stable layer at the top of its semi-conducting region would provide a skin-effect, reducing the smaller-scale field components (suggested by Stevenson (1979) and studied by Christensen and Wicht (2008); Stanley and Mohammadi (2008); Stanley (2010)).Furthermore, the difference in amplitude of its axial dipole field compared to the higher degree m = 0 components (Cao et al., 2020) indicates that there may be both a deeper dynamo region generating the strong dipole field, located between a dilute core and the helium rain layer, as well as a shallower layer adding the weaker latitudinally banded perturbations, operating between the helium rain region and a shallower, thin, stable layer.
However, for both planets the main uncertainty in the hypothesis is the origin, location, depth and strength of such a relatively shallow stable layer.A helium rain layer (Stevenson & Salpeter, 1977), providing a potential source of compositional stratifica-tion, is predicted to lie deeper than the extent of the zonal winds.In Jupiter, although there are some uncertainties concerning the H/He phase diagram, this would be below 86% radius based on ab initio EoS calculations of high-pressure experiments (Hubbard & Militzer, 2016;Lorenzen et al., 2011;Brygoo et al., 2021).In Saturn helium immiscibility may occur at around 65% radius, e.g.Morales et al. (2013).In both planets, however, there is not only a large uncertainty with regards to the depth of a helium rain layer but also no good estimate for its vertical extent.For the case of Jupiter the shallower regions, above where a helium rain layer is thought to reside, are potentially also more complex, based on the accurate gravity measurements from Juno, which suggests the existence of a shallow stably stratified region (Debras & Chabrier, 2019;Nettelmann et al., 2021;Debras et al., 2021), providing a potential link with the stable region associated with a quenching of the zonal winds.
In Wulff et al. (2022) we used purely hydrodynamic convection models to investigate the relationship between the degree of stratification of such a layer and the penetration of the winds, formed in the overlying convecting envelope, into the stable region below.We found that when the degree of stratification is strong, zonal flows form all the way to the higher latitudes, as is observed on both gas giants, even when imposing a no-slip boundary condition at the bottom of the stable layer.Furthermore, when encountering the SSL, the winds are quenched and geostrophy (i.e.their invariance with respect to the axis of rotation) is broken.However, the decay of the jet amplitude in this hydrodynamic study was still too gradual with depth to fit secular variation data.Furthermore, we expect that at sufficient depth the electrical conductivity will be large enough for magnetic effects to play a role.Therefore, it is crucial to investigate how this will influence both the damping of the jets in the SSL as well as their strength and latitudinal distribution in the overlying convective region.In our study we also test the concept of Christensen et al. (2020).In their simplified models the zonal flow was driven by an imposed ad-hoc force.In our models the zonal winds are driven self-consistently by the convective eddies, which implies that a potential feedback of the winds on the eddy dynamics is also accounted for.

Methods
We simulate thermal convection in a spherical shell rotating with angular velocity Ω • êz .The ratio of inner boundary radius, r i , to outer radius, r o , is 0.7.Only the upper part of the shell above 0.83r o is convectively unstable, whereas the lower part is stably stratified (described in detail in Section 2.3).We assume an exponentially varying electrical conductivity rising from a negligible value at r o to a moderate value at r i (see Section 2.4).We impose an axisymmetric dipolar magnetic field aligned with the rotation axis through a boundary condition at r i , which represents a field generated by a dynamo operating below r i .For our systematic study we use the Boussinesq approximation (i.e.incompressible flow), although we also perform additional simulations with the anelastic approximation (where a radially varying background density is prescribed).
The Boussinesq simulations are cheaper computationally and allow a wider parameter study.In this study, we keep all hydrodynamic parameters as well as the degree of stability in the SSL at fixed values, but we vary the magnetic field strength and the profile of the electrical conductivity.The anelastic simulations are carried out for a subset of these parameters in order to confirm that the trends observed also hold in the compressible models.

MHD Equations
As our primary analysis focuses on simulations that use the Boussinesq approximation, we give the governing magneto-hydrodynamic (MHD) equations here in their incompressible form (see Wulff et al. (2022) for the hydrodynamic equations under the anelastic approximation).The key features we incorporate are the radially varying mag-netic diffusivity λ(r) and dT c /dr, the imposed stratification profile, where T c is the background temperature.As we use a constant gravity, g, the equations then simplify to: where u is the velocity field, B is the magnetic field, and p is pressure.Temperature fluctuations ϑ are defined with respect to the hydrostatic reference state.We adopt a di- Boussinesq study involving a stable layer implemented in a similar way).The non-dimensionalised velocity is equivalent to a Reynolds number Re = ud/ν.The magnetic field is given in units of √ ρ o µλ i Ω, where µ is the magnetic permeability and λ is the magnetic diffusivity which we prescribe as an analytical radial profile.
The dimensionless control parameters that appear in the equations above are the Ekman number (E), Rayleigh number (Ra), Prandtl number (P r) and magnetic Prandtl number (P m).They are defined as where κ is the thermal diffusivity and α is the thermal expansivity.The magnetic Prandtl number P m, based on a reference value of the magnetic diffusivity, is kept at 0.5.However, the magnetic diffusivity at the lower boundary λ i is varied.

Hydrodynamic Control Parameters
We perform our simulations at a (nominal) Rayleigh number Ra = 6 × 10 8 , Ekman number E = 10 −5 and Prandtl number P r = 0.5.This yields a convective Rossby number of: so the Coriolis force dominates over inertia.
Some additional simulations are carried out under the anelastic approximation (see Wulff et al. (2022) for the governing equations), with polytropic index 2 and dissipation number 1, yielding a mild density stratification of ρ i /ρ o = 4. From Jones et al. (2009) and Gastine and Wicht (2012), for example, we know that the critical Rayleigh number increases with increasing density stratification.From the latter study we estimate that the increase is roughly two-fold, compared to our Boussinesq models.Therefore, to compare simulations with a similar degree of supercriticality, we double Ra for the anelastic simulations.
The values given above are based on the full shell width d.Table 1 also lists both non-dimensional numbers, re-scaled to the thickness of the convective region (the outer ∼ 57%).We also give the Rayleigh number based on the temperature (entropy for the anelastic cases) difference across the convective region alone, calculated from the horizontally averaged temperature (entropy) drop across the convecting region: where r c is the bottom of the convective region.Sim. and extreme cases B3.3 (green) and B4.0 (red).

Stably Stratified Layer
The region r > r c is fully convective, whereas at r < r s the full degree of stability has been reached, with a transition region at r s < r < r c .This is implemented by prescribing an analytic background entropy gradient profile defined, using auxiliary variable χ = (r − r c )/(r c − r s ), by: This is plotted in Figure 1a).In this study we keep r c = 2.77 = 0.831r o , r s = 2.68 = 0.804r o and A s = 200 (A s = 100 for the anelastic cases).Neutral stability is reached around 0.830r o .The ratio of the Brunt-Väisälä frequency, N , to the rotation rate, is: which is equal to 4.9, at r ≤ r s .This quantifies the effect of the restoring buoyancy force relative to the rotational forces and corresponds to a degree of stratification around the middle of the range studied in Wulff et al. (2022).This parameter is kept the same for the anelastic cases.

Magnetic Parameters
We vary the magnetic diffusivity λ, or the electrical conductivity σ = 1/λ, in this study but keep all other diffusivities (ν and κ) constant.We prescribe the magnetic diffusivity to be: For the profiles where λ would exceed 10 7 we cap it at this value to avoid numerical problems.The electrical conductivity scale height is This simple exponential profile gives the convenience of having a constant scale height throughout the shell.
In our reference model λ i = 1 and d λ = 1/ ln(10 8 ) ≈ 0.054.To investigate and distinguish the effects of a different local value of electrical conductivity and a different scale-height, we vary both d λ and λ i in a systematic parameter study (see Table 2).The electrical conductivity profiles of the extremes of the study, B3.3 and B4.0, are shown in Figure 1b).
An axial dipole field (poloidal ℓ = 1, m = 0 component) with amplitude B dip at the poles is imposed as a boundary condition at r i (negative at the North pole).The other poloidal components and the toroidal field are matched to a field in the inner core, obtained by solving the induction equation in the inner core for a constant value λ i of the diffusivity.At the outer boundary, r o , the magnetic field is matched to a potential field in the exterior.In this study we systematically vary the strength of the applied dipole, B dip .

Numerical Methods
All simulations in this study have been computed using the MHD code MagIC (available at https://github.com/magic-sph/magic).We use both the original Boussinesq version (see Wicht, 2002) and that which uses the anelastic approximation (Jones et al., 2011).
The governing equations given in section 2 are solved, with stress-free mechanical boundary conditions at both r i and r o and fixed entropy at the outer boundary and fixed entropy flux (downward in our models) at the inner boundary.This is done by expanding both velocity (or ρu in the anelastic cases) and magnetic fields into poloidal and toroidal potentials.For further details see Christensen and Wicht (2015).The potentials are expanded in Chebychev polynomials in the radial direction and spherical harmonics up to a degree ℓ max in the angular direction.
We use 145 radial grid-points for all simulations in the study.We use a non-linear mapping function (Tilgner, 1999) to concentrate the grid-points around the transition from convecting to sub-adiabatic.This ensures both the boundary between the two layers as well as the shell boundary regions are well-resolved, as illustrated in Figure 1.See the Appendix for details on the mapping For the reference case, labelled B1.1, we carried out one simulation without any imposed azimuthal symmetry, using azimuthal resolution n ϕ = 1280 and without hyperdiffusivity.For the other cases we introduced a four-fold azimuthal symmetry, reduced the number of grid-points to n ϕ = 864 and applied hyper-diffusion, where the diffusion parameters (thermal and viscous) are multiplied by the factor for ℓ ≥ ℓ hd , where ℓ hd = 250, D = 4 and β = 2.We verified that in the reference case the zonal winds formed and other features vital for our analysis did not change with imposed symmetry and hyper-diffusion.
All analysis was then based on the final stage of the simulations, which were integrated for 800,000 time-steps after they were fully equilibrated which is around 0.2τ ν (∼ 20, 000 rotations).  3 Results In our study we vary the strength of the imposed dipole field, B dip , the electrical conductivity at the inner boundary, σ i , and the conductivity scale height, d σ .The parameters are summarised in Table 2.We explore the surface zonal wind profiles, their extension into the interior and the mechanisms by which they are quenched.

Zonal Wind Distribution
The snapshot of our reference case B1.1, in Figure 2, shows that these simulations reproduce one of the key features found in the measurements of the zonal flows of the two gas giants: a set of alternating zonal jets reaching up the high latitudes.The equatorial prograde jet and its flanking retrograde jets dominate, but slightly weaker flows also persist up to the poles.These extend geostrophically, i.e. invariant with respect to z which is parallel to the rotation axis, throughout the convective region.We show the time-averaged, axisymmetric zonal flow for only one hemisphere of the hydrodynamic comparison case H in Figure 3a.Plotted on top of this is the surface profile as a function of the cylindrical coordinate s = r sin θ, i.e. the distance from the axis of rotation.
We observe that in case H, without either the additional magnetic forces or a mechanical rigid boundary condition which can act as a proxy for some force that brakes the jets, the jets are much wider and their amplitude (in this case that of the only retrograde jet present inside the tangent cylinder) only decreases slightly when reaching the stable layer.We note that a similar purely hydrodynamic simulation with a stress-free flow boundary shown in Figure 3d of Wulff et al. (2022) also shows a zonal flow pattern unlike that of Jupiter or Saturn, with a few strong jets inside the tangent cylinder (TC) that decay only weakly towards the inner boundary.The differences to the present case can be attributed to the anelastic approximation and a larger degree of stability in Wulff et al. (2022).
However, Figures 3b and c demonstrate that under the influence of finite conductivity and a large-scale magnetic field, the zonal flows develop a multiple jet structure.2) with the same colour-scale as Figure 2, with range ±6000.
On top of these are plotted the respective surface wind profiles as a function of s for the hemisphere shown.The thin vertical lines indicate the locations of the tangent cylinders associated with the bottom of the convective region, TC. d) shows the average zonal flow velocity inside the TC (defined by eq.12), as a function of the local Elsasser number evaluated at 0.8ro.See table 2 for the symbols for each case.
Furthermore, the jets are quenched effectively in the stable layer.The two cases shown, B4.0 and B3.3, are the extremes in the study.In B4.0 the conductivity starts out rather small at the inner boundary and drops rapidly with radius.This case shows strong zonal winds in the tangent cylinder reaching the polar region.In B3.3 the conductivity at r i starts out rather large and drops more weakly with radius.Here, significant jets are still found at mid-latitude, but they fade out at the high latitudes.The vertical extent of the convective region, i.e. the depth of the stable layer, is not altered in the study so the TC is in the same location and the equatorial prograde jet has the same width, with the peaks of the flanking retrograde jets being located on the TC.
The relation between the jet widths and jet amplitudes was confirmed to obey Rhines scaling well, when taking the convective region as the shell thickness (following the methodology detailed in Gastine et al. (2014)).This predicts that narrower jets are also weaker.
In order to quantify the strength of the axisymmetric zonal flow inside the tangent cylin-der (TC), we define the average surface zonal flow amplitude in this region as: where ⟨u ϕ (r o , θ)⟩ 2 is the time-averaged, axisymmetric, rms surface zonal flow and θ c = sin −1 (r c /r o ), i.e. the colatitude associated with the location of the TC at the surface.
This definition broadly captures both the extent and strength of the zonal flow and facilitates a comparison between all cases.
We observe that simulations with a stronger imposed dipole field strength, B dip , and those with higher electrical conductivity, σ, have weaker winds inside the TC.We use a local Elsasser number Λ(r) = B dip (r) 2 σ(r)/ρΩ as a proxy for the local strength of the Lorentz force, relative to the Coriolis force, although in this study we only explicitly test the dependency on B dip and σ.This expresses not only the radial variation of the electrical conductivity but also the r −3 dependence of the dipole field strength (in our definition we use the axial dipole field amplitude at the poles).
In Figure 3d we therefore plot U surf , as a function of the Elsasser number evaluated at 0.8r o , i.e. in the upper part of the stable layer, just below r s .The extremes of our parameter sweep are Λ(0.8ro ) = 1.30• 10 −5 in case B4.0, up to 1.80 • 10 −2 in case B3.3 (see Table 2).
The plot suggests that if magnetic forces remain insignificant near the SSL boundary, strong zonal winds can develop and be maintained in the overlying convecting region and are independent of the magnetic effects coming into play deeper in the stable region.However, when magnetic effects become more pronounced in the upper part of the stable layer, the zonal flow inside the TC becomes somewhat more diminished, in particular at high latitudes.Within our parameter sweep this is not a dramatic effect.
As our focus is on models that have strong jets inside the TC so we do not go beyond case B3.3.We would expect these to disappear if the semiconducting region begins at even shallower depths and Λ(0.8r o ) is increased by even just one more order of magnitude.We track the jet amplitude as a function of radius in the SSL.This is illustrated in Figure 5a) where we show the locations of the maxima/minima of the jets inside the TC for simulation B1.1.This tracking is vital as the locations of the peak velocity is no longer z-invariant in the SSL in contrast to the convective region, as can be seen in fig-

Flow Amplitude Versus Depth
ure 3b and c.
Figure 5b shows u e ϕ along the centres of these jets (we use superscript e to denote the extrema of u ϕ as a function of latitude).This also highlights the strong equatorial symmetry of these particular simulations where the northern/southern hemisphere jet pairs have almost identical velocity profiles.Finally, Figure 5c shows the same velocity profiles, with each one normalised by the respective jet velocity at r c .This plot clearly   illustrates that the relative decay with depth is rather similar for all jets, independent of their location inside the TC.
We average the radial profiles of all jets inside the TC, normalised by their velocity at the bottom of the convecting region, for each case to quantify the zonal wind decay.Figure 6b compares averaged profiles for the end-member simulations of each of the Boussinesq sets while Figure 6a shows the respective Elsasser number profiles in the SSL.
Profiles with identical Λ(r) (sets B1. and B2.) nearly perfectly overlap which highlights that the Elsasser number is the crucial parameter here; doubling the axial dipole field strength has exactly the same effect as quadrupling the electrical conductivity.When considering the other profiles shown we clearly see that in simulations with the lowest Elsasser numbers the decay of the zonal wind is very gradual.This is illustrated in fig- ure 6c where we plot the ratio of the jet amplitude at 0.8r o and the amplitude at r c , again averaging over all 8 jets to obtain one value per simulation.Therefore, magnetic effects are crucial in reducing the penetration distance of zonal winds into the SSL.

Magnetic Field Induction
Figure 4b shows the horizontally averaged induced magnetic field components for the reference case.The induced axisymmetric toroidal field is almost as strong as the dipole field at the lower boundary, for this case, but drops off rapidly with radius.The induced axisymmetric radial and latitudinal fields, ∆B r and ∆B θ , i.e. the perturbations of the imposed poloidal field, are almost three orders of magnitude smaller but do not drop off in amplitude as sharply over the SSL.
We investigate the difference in the induction when changing the imposed axial dipole field strength in sin θ .From this we see that the assumption that a stronger dipole will lead to a stronger induced field holds only when the simulations with different B dip have similar distributions of u ϕ .However, in Figure 6b we see that the zonal wind is quenched very effectively in case B1.3, so u ϕ is almost zero at mid-depth of the stable layer, while this only happens near the bottom in case B1.0.
Therefore, the Ω-effect is stronger for case B1.0 than for B1.3 in the lower third of the stable layer and the amplitude of the induced field, B ϕ , exceeds that of B1.3 significantly, as can be seen in the bottom panel of Figure 7.In fact, for case B1.0 it also exceeds the amplitude of the imposed dipole field at the lower boundary which for this case is B r (r i , θ = 0) = −0.25.Thus, the model, B1.0, which has the most interaction between the zonal flow and the electrically conducting region (i.e. over the greatest depth range) actually has the weakest induced field strength near the top of its conducting region.
The morphology of the induced B ϕ field can also be better understood by comparing the contributions of the two terms that make up the Ω-effect; the radial and the latitudinal shear of the zonal wind.This is shown on the right in Figure 7, where the rms amplitude of the two terms has been averaged horizontally to produce radial profiles.For both cases, B1.0 and B1.3, the radial shear is the more dominant term throughout the shell.Therefore, the decay of the jets with depth produces a stronger gradient than the transition between oppositely flowing jets.This also leads to the induced azimuthal field being strongest almost exactly on the zonal wind peaks.

Zonal Wind Truncation Mechanism
A thermo-magnetic wind equation can be derived by taking the ϕ-component of the curl of the Navier-Stokes equation, then averaging over azimuth and assuming steady state: where ω = ∇ × u is the vorticity and j = ∇ × B. The first three terms are from the advection term, the fourth and fifth terms are from the Coriolis force and buoyancy, respectively.The last two terms are from the Lorentz force and the viscous force.We find that in our simulations, the equation can be reduced to: as all other terms were found to be negligible.This is shown in Figure 8b and c where we plot the vertical gradient of the zonal wind (first term of eq.14) and the latitudinal temperature gradient (second term of eq.14).The two are in nearly perfect balance; the magnetic term of the thermo-magnetic wind equation is negligibly small.As in the hydrodynamic simulations in Wulff et al. (2022), the decrease of the zonal wind in the stable layer is controlled by a thermal wind balance.The associated density perturbation is caused by a meridional flow.As Lorentz forces play a critical role for the penetration of the winds into the stable layer, this should happen via their influence on the meridional circulation.To elucidate this, we consider (as in Wulff et al., 2022) the time-averaged (denoted by ⟨⟩) axisymmetric, azimuthal component of the Navier-Stokes equation, given by: This includes the: 'advective' force F Ad , Coriolis force F C and viscous force F ν .The forces associated with the Reynolds stresses and the Maxwell stresses are F R , and F M a and F M na respectively, where the former is the contribution from the large-scale (axisymmetric) magnetic field components and the latter is from the small-scale (non-axisymmetric) field components.We find that the advective force remains negligibly small and therefore omit it in Figure 9 where we show the zonal force balance.Furthermore, the Maxwell stresses arising from the correlation of the small-scale magnetic field components, F M na , also remain very small, even at depth and are thus also not shown in Figure 9.This is because in our study the stable layer suppresses small-scale flows so effectively.The conductivity distribution implies that the Lorentz forces only start acting in the SSL in these simulations, where only very weak non-axisymmetric induced magnetic field components contribute.
The first panels in Figure 9 show the Coriolis force, which is directly proportional to the s-component of the meridional flow.This meridional flow is driven in the convecting (not shown) and transition region, where the associated Coriolis force is balanced mainly by the Reynolds stresses.The Reynolds stress force is enhanced in the transition region, by the same mechanism as in the purely hydrodynamic study (Wulff et al., 2022), where radial flows and also all small-scale motion is quenched (see Figure 4).Therefore there is a sharp drop-off in the correlation of the convective flows just below r c , leading to large derivatives with respect to s and z (see eq. 15 for the definition of F R ).The large Reynolds stress force is primarily balanced by the Coriolis force F C of an enhanced meridional circulation.Inside the SSL there is a good match of F C and the force associated with the Maxwell stresses, F M a .This is the essential difference to the hydrodynamic models where only viscosity can balance F C in this region.In the MHD case, the meridional flow remains significant in the SSL, so entropy perturbations are induced (see Figure 8) and the zonal flow can be quenched more effectively in the SSL.This is broadly in agreement with the mechanism proposed by Christensen et al. (2020), where the winds were driven by an ad-hoc force rather than self-consistently by Reynolds stresses.We observe that the viscous force also plays a significant role in the SSL, as the zonal flow 12), as a function of the local Elsasser number evaluated at 0.8ro.See table 2 for the symbols for each case.
velocity is decreasing rapidly.At the much lower Ekman numbers that apply to the gas planets, viscosity is expected to play no significant role.

Anelastic Simulations
For seven models with different field strengths and conductivity profiles we replaced the Boussinesq approximation by the anelastic approximation in order to test its impact on the results (sets A1. and A2. in Table 2).Qualitatively, the zonal flows formed in these simulations are very similar to their Boussinesq counterparts, with the strongest jets being the prograde equatorial jet and its flanking retrograde jets, complemented by another four weaker jets inside the tangent cylinder (see figures 10a and b).In these simulations we also observed some time-variability in the zonal flow structure, similar to that discussed in Wulff et al. (2022), which we do not explore further within this work.We also test the relationship of the local Elsasser number and the zonal wind penetration distance for these anelastic models.We use the same analytical technique de- Table 2 for the symbols for each case.
scribed in section 3.2 and track the jet amplitudes in the stable layer.This is shown in Figure 11.As Figure 11a shows, the two sets A1. and A2. have the same 4 radially varying Elsasser number profiles (with A1.1 forming part of both sets).However, in one set the axial dipole field strength, B dip , is varied while in the other the electrical conductivity profile is varied (the conductivity scale height remains the same).Figure 11b shows that the models with the same Λ(r) have near-to identical zonal flow decay in the stable layer, with the zonal winds in models with a stronger imposed dipole strength or a greater electrical conductivity (A1.3 and A2.3 respectively) being quenched most effectively.Although the variation of the background density with radius is rather weak in our models, this suggests that our observations from the Boussinesq models also hold when there is a variable background density.Furthermore, Figure 11c, where we plot the ratio of jet amplitude at 0.8r o to that at r c , shows that the 1/ρ(r) dependency of the local Elsasser number leads to a more gradual damping of u ϕ in these models compared to cases B1. and B2., their Boussinesq equivalents.Figure 11c has the same axes as Figure 6c to highlight that these models fit on the same trend line.This is possible for this analysis as the relative decay is evaluated, while the absolute jet amplitudes are difficult to compare with the Boussinesq models.

Discussion and Conclusions
We find that the amplitude and latitudinal extent of zonal flow in the convective region, depends directly on the amplitude of the magnetic forces near the top of the underlying stable region.If these are negligible, due to both a weak dipole field strength and very weak conductivity, the zonal flow at the surface develops a structure and amplitude independent of the magnetic effects acting deep in the stable region below.If Lorentz forces are non-negligible at the bottom of the convective region, they will impact the jets formed above, in particular diminishing those inside the tangent cylinder (see Figures 3d) and 10c)).
The penetration distance of zonal flows into the stable layer is dependent on the product σB 2 at depth.For a fixed profile of σB 2 , it can be expected from Wulff et al. (2022) that the degree of stratification, N/Ω, also influences the damping of the zonal winds in the stable layer, as well as other parameters.Christensen et al. (2020) suggest that for a fixed σ and B, and in the limit of negligible viscosity, the combination (N/Ω) 2 E −1 κ Figure mensionless formulation where the reference length scale is the shell thickness d = r o − r i , where i denotes the inner boundary values and o denotes outer boundary.Time is given in units of the viscous diffusion time τ ν = d 2 /ν, where ν is the fluid viscosity.The temperature scale is normalised by the value of the gradient of the background temperature at the outer boundary |dT c /dr| o , multiplied by d (see Gastine et al. (2020) for a

H
Figure 1.a) dTc/dr profile (black) described in Section 2.3.The grey shaded region indicates rs < r < rc while the dark grey region is fully stratified.The radial grid-point separation is shown in orange (right y-axis).b) electrical conductivity σ = λ −1 for reference case B1.1 (blue)

Figure 2 .
Figure 2. A snapshot of the azimuthal flow, u ϕ , for the reference case B1.1.Both plots use the same colour-scale with a dynamic range of ±6000, where red (blue) indicates prograde (retrograde) flow.a) View onto the surface of the spherical shell from the North Pole.b) Front view of the surface flow on the left and a cut down to the bottom of the convecting layer on the right.

Figure 3 .
Figure 3. Time-averaged axisymmetric zonal flow for simulations H in panel a), B4.0 in panel b) and B3.3 in panel c) (see Table2) with the same colour-scale as Figure2, with range ±6000.

Figure 4a )
Figure 4a) shows the horizontally averaged rms velocity components for the reference simulation as a function of radius, where solid lines show the axisymmetric components (labelled with an overbar) and dashed lines the non-axisymmetric components (indicated by a prime).In the convective region, the convective flow amplitude (u ′ ϕ , u ′ θ and u ′ r ) is almost an order of magnitude weaker than the rms zonal wind amplitude (the jet peaks themselves are even stronger).Upon reaching the SSL, radial flow components are quenched most effectively and amplitudes drop by almost two orders of magnitude.At least part of the remaining radial motion seen in figure 4 may represent wave motion (gravity waves, inertial waves) and no overturning motion.Right at the SSL boundary both the latitudinal component of the meridional flow, u θ , and the horizontal components of the convective flow, u ′ ϕ and u ′ θ , increase very slightly which may be attributed to the deflection of the radial flows.However, further into the SSL all other flow components are damped.We analyse this in more detail for the zonal flow.

Figure 4 .
Figure 4. Radial profiles of time-and horizontally-averaged a) velocity and b) magnetic field strength (given in Λ) for the reference case B1.1.The dashed lines show the average nonaxisymmetric flow (field strength) and the solid lines are the axisymmetric parts, where colours indicate the three components.For B θ and Br we subtract the dipole component, of which the average amplitude is shown by the black dotted line.The dark grey (grey) shading indicates the (transition into the) SSL.

Figure 5 .
Figure 5. Illustration of the jet tracking method described, applied to the reference case B1.1.a) shows the zonal flow pattern in the SSL (r < rc).The black lines show the locations of the zonal flow extrema (denoted by superscript e) and the grey lines indicate lines of constant s. b) shows the peak amplitudes of these 8 jets as a function of radius and in c) we normalise this by the jet flow velocity at r = rc.

Figure 6
Figure 6.a) Elsasser number as a function of depth in the SSL (r < rc) for the end-members of each of the Boussinesq sets.b) shows the normalised jet amplitude profiles as shown in Figure 5c where a single profile is obtained for each case by averaging over all 8 jets.c) ratio of zonal flow amplitude at 0.8ro and the jet flow velocity at r = rc, obtained from the averaged profiles shown in b) and the remaining models omitted on this plot.See Table2for the symbols for each case.

Figure 7 .
On the left the top three plots are latitudinal profiles of u ϕ , ∆B r = B r −B dip r and B ϕ at 0.8r o .These are all as we may expect, with B1.3 having the strongest induced magnetic fields.This is due to it having the largest imposed dipole amplitude, eight times stronger than case B1.0 with the same conductivity profile, leading to a larger Ω-effect.The Ω-effect describes the induction of axisymmetric toroidal field by the shearing of the axisymmetric poloidal field by differential rotation and has the two components B r ∂ r u ϕ r and B θ sin θ r ∂ θ u ϕ

Figure 7 .
Figure 7. Left: Latitudinal profiles for cases B1.0 and B1.3.The top panels show u ϕ and ∆Br = Br − B dip r on 0.8ro.The lower panels are B ϕ at 0.8ro and ri.Right: Horizontally averaged, rms amplitude of the two terms in the Ω-effect for the same cases, where dashed(solid) corresponds to B1.0(B1.3)and dark(light) blue corresponds to the radial(latitudinal) shear term.

Figure 8 .
Figure 8.The transition region and SSL of the reference case B1.1, shown for the northern hemisphere.All terms are zonally and temporally averaged.a) u ϕ , the zonal flow, b) the vertical gradient of the zonal flow and c) the latitudinal entropy fluctuation.The thin horizontal line indicates rs, i.e. the bottom of the transition region into the SSL.

Figure 9 .
Figure 9.The zonally and temporally averaged azimuthal component of the force balance, shown for cases a) B1.0 and a) B1.3, in the transition region and SSL for one hemisphere.The panels show the Coriolis term, F C ; Reynolds stresses, F R; viscosity, F ν ; and Lorentz forces, F M a.The horizontal black line indicates rs, which is also roughly the depth at which Lorentz forces become significant and F R becomes negligible.Solid (dashed) grey contours indicate positive (negative) u ϕ .

Figure 10
Figure 10.a) and b) show the time-averaged axisymmetric zonal flow for simulations A1.0 and A2.3 (see Table 2) with the same colour-scale as Figure 2, with range ±6000.On top of these are plotted the respective surface profiles as a function of s for the hemisphere shown.The thin vertical lines indicate the locations of the tangent cylinders associated with the bottom of the convective region, TC. c) shows the average zonal flow velocity inside the TC (defined by eq. b) we first note that in case B1.3 where the dipole strength is increased (Figure9b)), the Lorentz forces already begin to act in the SSL transition region.This illustrates how they are able to impact the structure of the zonal winds in the top part of the stable region.Deeper in the stable region similar meridional circulation cells develop for both models, to balance the Lorentz forces.However, they are shifted upwards in case B1.3 relative to B1.0.In model B1.0 the winds only reach nearzero amplitude near the inner shell boundary and the transition from equator-ward (poleward) flow in the high-(mid-) latitude region to oppositely flowing meridional circulation occurs close to this boundary.In case B1.3 the winds are already quenched at around 0.74r o which is where the circulation cells are centred in this model.

Figure 10
Figure 10 is the counterpart to Figure 3, showing the dependence of the rms zonal flow amplitude at the surface, inside the TC, on the local Elsasser number.While Λ covers a smaller range than the Boussinesq study, the same trend is clearly seen: stronger winds develop in models where magnetic effects, characterised by B 2 dip σ, become significant only deeper into the stable layer.

Figure 11
Figure 11.a) Elsasser number as a function of depth in the SSL (r < rc) for all 7 anelastic models.b) shows the normalised jet amplitude profiles as shown in figure 5c where a single profile is obtained for each case by averaging over all 8 jets.c) shows ratio of zonal flow amplitude at 0.8ro and the jet flow velocity at r = rc, obtained from the averaged profiles shown in b).See Figure 2.

Table 1 .
Nominal Ekman and Rayleigh numbers are based on the full shell thickness and the surface entropy flux.Their re-scaled values, Ec and Rac, are based on the thickness of the convective region dc = ro − rc ≈ 0.57.The re-scaled Ekman number for the stable region Es is also given, based on ds = rc − ri ≈ 0.43.Ra∆ is the Rayleigh number defined by Eq. 8.As is the value of dTc/dr (d S/dr for anelastic cases) in the stable region.

Table 2 .
Simulations carried out with critical varied parameters given.The reference case is in bold.