A Global Simulation of the Dynamo, Zonal Jets, and Vortices on Saturn

The fluid dynamics in planet Saturn gives rise to alternating east-west jet streams, large cyclonic and anticyclonic vortices, and a dipole-dominant magnetic field which is highly axisymmetric about the planetary rotation axis. Modelling these features in a self-consistent manner is crucial for understanding the dynamics of Saturn's interior and atmosphere. Here we report a turbulent high-resolution dynamo simulation in a spherical shell which produces these features simultaneously for the first time. A crucial model ingredient is a long-hypothesised stably stratified layer (SSL), sandwiched between a deep metallic hydrogen layer and an outer low-conductivity molecular layer, born out of limited solubility of Helium inside metallic Hydrogen at certain depths. The model spontaneously produces polar cyclones and significant low and mid latitude jet stream activities in the molecular layer. The off-equatorial low-latitude jet streams partially penetrate into the SSL and interact with the magnetic field. This helps to axisymmetrize the magnetic field about the rotation axis and convert some of the poloidal magnetic field to toroidal field, which appears as two global magnetic energy rings surrounding the deeper dynamo region. The simulation also mimics a distinctive dip in the fifth spherical harmonic in Saturn's magnetic energy spectrum as inferred from the Cassini Grand Finale measurements. Our model highlights the role of an SSL in shaping the fluid dynamical and magnetic features of giant planets, as exemplified at Saturn.


I. INTRODUCTION
The gas giant planet Saturn, well known for its spectacularly bright icy rings, has many more defining properties. On its visible surface, the different colored gases and clouds sketch the profiles of the planet-circulating east-west winds which alternate in direction from the equator to the poles and reach speeds of several 100s of meters per second (Lavega 1982, Sanchez-Lavega et al. 2000. One of these zonal jets takes the form of a regular hexagon in the north polar region (Godfrey 1988). Embedded among these jets are large vortical storms which can rotate either in the clockwise or anticlockwise direction (Trammell et al. 2014(Trammell et al. , 2016. A polar cyclonic vortex is found at each geographic pole (Sánchez- Lavega et al. 2006). In addition, the angle between Saturn's spin axis and the axis of the dominant magnetic dipole (<0.007 • ) is almost 3 orders of magnitude smaller than that of the Earth (Cao et al. 2020, Dougherty et al. 2018). Furthermore, its non-dipolar magnetic field is also highly axisymmetric about the rotation axis, much more so than any other object in the solar system (Cao et al. 2020, Dougherty et al. 2018. These properties of Saturn are a direct consequence of the rotating fluid dynamics occurring in its interior and atmosphere, and provide a great opportunity for building and constraining theoretical models of Saturn's interior fluid dynamics. Through magnetohydrodynamic inter- * To whom correspondence should be addressed: rakesh yadav@fas.harvard.edu † hcao@epss.ucla.edu ‡ jeremy bloxham@harvard.edu actions, the electrical conductivity of a fluid inside a gas giant planet like Saturn will likely act as an important factor in demarcating regions with different flow properties (Stevenson 2003). It has been hypothesised (Liu et al. 2008) that due to the varying electrical conductivity of liquid hydrogen inside Saturn (and Jupiter) (French et al. 2012, Nellis 2000, the magnetohydrodynamic interactions largely separate Saturn's interior into an outer strong zonal flow region and an inner one with much reduced zonal flows. The detailed mechanism involved in this separation is an active area of research (Christensen et al. 2020, Gastine andWicht 2021).
A significant body of work has been dedicated to understanding how the zonal jets and vortices are created in the outer layer of Saturn. Broadly speaking, there are two interpretations. In one, which we may call a "bottom up" approach see Kulowski et al. (2021) for a recent discussion of this issue), the zonal jets are considered to be surface manifestations of deep (thousands of kilometers) concentric, axially-aligned, rotating cylinders generated through the convective, rotating turbulence. This idea was put forth by Busse (1976) and has been extensively explored in many numerical simulations studies (e.g., Aurnou and Olson 2001, Christensen 2001, Gastine et al. 2014a, Heimpel et al. 2005, Jones and Kuzanyan 2009, Kaspi et al. 2009, Yadav and Bloxham 2020. In the other approach, the zonal jets are thought to exist in a thin atmospheric layer (100s of kilometers or less), similar to the Earth's atmosphere, within which the shallow water approximation is typically employed (Cabanes et al. 2020, Cho and Polvani 1996a,b, Lian and Showman 2010, Liu and Schneider 2010, Williams 1978. Here, the driving sources are commonly assumed to be latent heat released through cloud condensation and solar irradia-tion. The gravity data collected by the Cassini Grand Finale has shown that Saturn's gravity perturbations can be explained by slightly modified surface zonal winds penetrating to nearly 9000 km . Therefore, the zonal jets are very deep indeed. However, teasing out the jet-driving mechanism remains difficult because even a shallow production mechanism confining in the atmosphere of Saturn could also lead to very deep jets Showman 2008, Showman et al. 2006).
Saturn's highly axisymmetric magnetic field is exceptional in light of the Cowling's theorem (Cowling 1933) which states that a dynamo in an incompressible fluid cannot maintain a purely axisymmetric magnetic field. Although, several of the original requirements have been relaxed in subsequent work (Hide andPalmer 1982, Kaiser 2018), we can safely assume that a dipole dominant magnetic field perfectly aligned with the spin axis cannot be generated by a turbulent dynamo itself. We can extend this line of thought to assume that even the magnetic fields like the one existing in Saturn, where the dipole tilt angle is extremely small, likely require special mechanisms for becoming highly axisymmetric. Stevenson (1979) theoretically modelled the interaction of Helium and Hydrogen at physical conditions relevant to giant planets (also see Salpeter 1973). He proposed that Helium concentration may be at or close to saturation as hydrogen transitions from molecular to metallic with increasing depths. This would lead to differentiation ('rain out') of Helium and the development of a stably stratified layer in the outer parts of the metallic hydrogen layers (Hubbard and DeWitt 1985, Salpeter 1973, Stevenson 1982a. Recent ab initio studies suggest immiscibility of Helium in hydrogen for pressures lower than 1 Mbar at Saturn-like conditions (e.g., see Lorenzen et al. 2009, Morales et al. 2013, which is around 65% of Saturn's radius. However, a good estimate of the thickness of such a stable layer remains highly uncertain (Fortney et al. 2011, Stevenson andSalpeter 1977). Using such a layer as an ingredient, a mechanism was proposed to axisymmetrize a planetary dynamo: "If a conducting fluid shell is undergoing spin-axisymmetric differential rotation and overlies the dynamo generating region of a planet then it is capable of greatly reducing the non-spin-axisymmetric components of the generated field" (Stevenson 1982b).
This idea has been tested by several subsequent models. A kinematic dynamo study showed that an SSL can help axisymmetrize the magnetic field, but not always (Love 2000). Fully non-linear dynamo simulations performed with an overlying SSL showed that the SSL was very effective at removing non-axisymemtric small-scale magnetic fields (Christensen and Wicht 2008). Another study imposed latitudinal heat-flux variations to excite axisymmetric flow inside an SSL on top of a dynamo region and showed that the flows inside the SSL can help either axisymmetrize or destabalize the magnetic field depending on their properties (Stanley 2010, Stanley andMohammadi 2008). In the most recent study (Yan and Stanley 2021) in this context, a thick SSL with an imposed latitudinally-varying heat-flux perturbation is used to promote Saturn-like dipole dominant magnetic fields with very small dipole tilt angles of about 0.066 • which is still about one order of magnitude larger than the latest observational upper bound (Cao et al. 2020).
A few studies have investigated models without an SSL (and hence also without any latitudinally-varying boundary heat-flux) in which only the extremely steep drop in the electrical conductivity in the outer layers of Saturn (and Jupiter) is considered. It is generally found that strong zonal flows are excited (in the equatorial region) in the outer layer where the electrical conductivity is lower Jones 2018, Duarte et al. 2013). The dipole dominant solutions, a requirement for modelling Saturn's dynamo, are less stable in these setups. However, we recently showed (Yadav et al. 2022) that such setups without an SSL can produce dipole-dominant solutions with exceptionally small dipole tilts of ≈ 0.0008 • if the non-linearities in the simulation (promoted by low convective supercriticalities) are small enough. The axial quadrupole in that model, however, is many orders of magnitude smaller than that observed at Saturn.
To summarize, the zonal jets on Saturn are typically modelled separately from the dynamo, while the interior dynamo is modelled either with an SSL or with an overlying low-conductivity layer. In this study we try to consolidate these approaches and present a global model where we simulate the deep dynamo, an SSL, and the molecular hydrogen layer on top. It must be noted, however, that our "global" model still lacks the meteorological physics (e.g. cloud formation, rapid density variation, and solar insolation) occurring in the outer few hundred kilometers of Saturn. The model details given below will make it clear that our main aim is to investigate the non-linear interactions between different dynamically important regions in giant planets -specifically the interaction of an SSL with other regions. We do not attempt to match the exact background state of an interior model proposed for Saturn. To date, the state of Saturn's interior is an active area of research; for recent developments see Dewberry et al. (2021), , , , Movshovitz et al. (2020).

A. MHD Equations
For our simulations, we use the same tool and methodology as used by Gastine and Wicht (2021) where they investigate the effect of a stably stratified layer on the fluid dynamics and the resulting dynamo mechanism in the context of Jupiter. We refer the reader to their study for a detailed description of the methodology employed to model a sandwiched SSL; also see Dietrich and Wicht (2018). Here we only highlight some of the essential de-tails.
To simulate some of the magnetohydrodynamical processes operating inside a giant planet, we first assume the model planet has a perfectly spherical shape and is bounded by an inner boundary at r i and an outer boundary at r o . The spherical shell rotates about a fixed spin axis (ẑ) with angular velocity Ω. The motions of the fluid inside the spherical shell are modelled using the LBR (Lantz-Braginsky-Roberts) version of the anelastic approximation (Braginsky andRoberts 1995, Lantz andFan 1999). This approximation assumes that a thermodynamic quantity, for instance, x(r, θ, φ, t), is a combination of a background state,x(r), and a perturbation, x (r, θ, φ, t). We work with a non-dimensional system of variables and use the shell thickness d = r o − r i as the length scale, and the viscous diffusion time d 2 /ν, where ν is the kinematic viscosity, as the time scale. The magnetic field, B, is normalized by √ρ o µ o λ i Ω, whereρ o is density on the outer boundary, µ o is magnetic permeability, and λ i is magnetic diffusivity on the inner boundary. Note that the square of the magnetic field magnitude in this unit is known as the Elsasser number. The entropy, s, and the fluid density, ρ, is scaled by their values at the outer boundary.
The following system of equations govern the time evolution of velocity ( u), entropy, and magnetic field: ρT ∂s ∂t + u · ∇s + u r ds dr = 1 P r ∇ · ρT ∇s + P r Di In the above equations, p is pressure, T is temperature,α is the dimensionless expansion coefficient, g is gravity vector (gravity magnitude increases linearly with radius), λ norm is the radially-varying magnetic diffusivity normalized by its value at the inner boundary. The traceless rate-of-strain tensor S in Eqn. 2 is given by where δ ij is the identity matrix. The viscous heating contribution in Eqn. 3 is given by The non-dimensionalization process introduces the following fundamental control parameters in the equations: where subscript 'o' denotes the values at the outer boundary.

B. Model Setup
Unlike the study of Gastine & Wicht (Gastine and Wicht 2021) where an interior model of Jupiter was taken as a strong guiding factor, here we only model interactions among some features which we consider important for the interior dynamics of Saturn-like planets. Therefore, our study is an exploratory one. The primary aim being to study how a thick SSL might affect the workings of the planetary dynamo and the deep atmosphere. The background profile of the thermal variables is given in Supplementary Figure 1.
We set the inner boundary of the shell at r i = 0.27r o . The convection in the fluid shell is driven by an imposed entropy contrast. Within the spherical shell, we assume that there is a stably stratified layer between 0.62r o and 0.8r o which resists large-scale overturning convection and provides a restoring force to radial fluid motions. This is achieved by imposing a constant positive radial gradient (a free parameter in the model) in the background entropy profile within the SSL and a value of -1 elsewhere in the shell. The imposed entropy radial gradients in the convective layer and in the SSL are joined smoothly with each other using hyperbolic tangent functions. We refer the reader to (Gastine and Wicht 2021) for more details of this implementation strategy (also see Gastine et al. 2020, Takehiro andLister 2001). An important parameter in this context is the ratio of the Brunt-Väisälä frequency and the rotation rate (Ω). With the control parameters chosen in the simulation, this ratio is about 1.5 in the SSL. In the study of Gastine and Wicht (2021), this ratio is about 10. The lower value utilized in our setup allows for a greater penetration of convection into the SSL (e.g., see Takehiro and Lister 2001).
The viscosity and the thermal conductivity of the fluid are assumed constant throughout the spherical shell. The electrical conductivity, however, has a strong radial dependence: it stays constant from r i to about 0.78r o and exponentially decreases by more than five orders of magnitude within a thickness of about 0.1r o ; outside of 0.9r o , the fluid is assumed to be electrically insulating.
With this model setup, there effectively exist four dynamically important regions which self-consistently interact with each other (see Fig. 1a): a convective dynamo region (from 0.27r o to 0.62r o ), an electrically conducting SSL (0.62r o to 0.8r o ), a 'semi' conducting and convective region from 0.8r o to 0.9r o , and a purely hydrodynamic and convective region (from 0.9r o to r o ).

C. Control Parameters and Boundary Conditions
The choice of the control parameters is largely driven by the available computing resources, a need for lowering the fluid viscosity as much a possible (quantified by the Ekman number) and pushing the (magnetic) Reynolds number to high values. Given these constraints, the main control parameters we could model are: Ek = 7 × 10 −7 , Ra = 2.2 × 10 11 , P r = 1, P m = 0.4, Di = 6. The mechanical boundaries at r i and r o are impenetrable and stress free. For the magnetic field, the lower boundary at r i acts as a conductor having the same electrical conductivity as the fluid in the dynamo region; the radial level 0.9r o and higher act as an electrical insulator for the magnetic field. The entropy is assumed constant on each boundary at r i and r o with a prescribed entropy contrast between the two.
Due to the complexities involved in running such a demanding low-viscosity simulation, several compromises need to made: (1) It is not possible to run such a simulation with a small seed magnetic field and let the dynamo build up the large scale field since it will take impractically large computing power. Therefor, we start the simulation from a case having a saturated dynamo state. (2) The Ekman number needs to be decreased in incremental stages to reach 7×10 −7 . In the reported simulation here, we started with a simulation at Ek = 1.3 × 10 −6 and decreased the Ekman number to 7 × 10 −7 in stages in the first ≈15% of the total time evolution. (3) The Rayleigh number was increased to 2.2 × 10 11 from 4 × 10 10 in the first ≈50% of the time evolution. (4) The P r was constant throughout the simulation but the P m was lowered from 0.6 to 0.4 in the first ≈50% of the time evolution. (5) Due to the low viscosity associated with Ek = 7 × 10 −7 , the flow exists on very small length scales which were not possible to be resolved given the computing grid we used. Therefore, a numerical technique called hyperdiffusion was used throughout the simulation where flows on length scales smaller than a set value are suppressed more than those with larger length scales. This is a commonly utilized strategy in low Ekman number simulations (Gastine and Wicht 2021, Heimpel et al. 2022, Kuang and Bloxham 1999, Soderlund et al. 2012). The utilized hyperdiffusion was such that the Ekman number gradually increases by about 100 after spherical harmonic (SH) degree ∼70 in the first half and after SH degree ∼100 in the second half of the simulation.
The simulation could only be evolved for about 0.0045 viscous diffusion time, or about 0.011 magnetic diffusion time for the parameters mentioned above. It translates to about 1000 rotations of the spherical shell. The relatively short time span of the simulation, measured by the diffusion time, is an unavoidable drawback at such extreme parameters. Therefore, the statistical steadiness of the simulation is not guaranteed. The short time evolution of the simulation still consumed about 5.5 million core hours which is a rather large amount of computing time. Despite these limitations, we believe the richness of the results obtained warrants a documentation and a discussion.

D. Simulation Code and Grid Size
The Eqns. 1 through 5, operating under the above control parameters and boundary conditions, are solved using the open source code MagIC (https://magicsph.github.io) (Gastine and Wicht 2012). It uses the pseudo-spectral approach and projects the primitive variables in the latitudinal-longitudinal directions onto the spherical harmonic functions while the radial component is projected onto the Chebyshev polynomials. It also uses the toroidal-poloidal decomposition to maintain strict divergenceless condition for the relevant quantities. MagIC utilizes the open source library SHTns (Schaeffer 2013) to performe the spherical harmonic transformations. The grid size, given as number of grid points in longitudinal, latitudinal, and radial, respectively is [2112,1056,480] for the first 90% of the simulation evolution and increases to [2560,1280,640] in the final 10%.

III. RESULTS AND DISCUSSION
In this section we present our results. We begin by looking at the properties of the surface and internal zonal flow. We then describe some properties of large vortices that are produced near the simulation surface. In the later half of the section, we present the properties of the magnetic field produced by the dynamo region. We end the section with a comparison of the magnetic field morphology with the one observed by Cassini at Saturn.

A. Surface zonal flow
We begin by first describing the flow properties arising in the simulation. In Fig. 1b we plot the instantaneous, axisymmetric zonal flow on the outer boundary of the simulation. As a reference we also plot the corresponding values obtained for Saturn from cloud tracking (Read et al. 2009, Sanchez-Lavega et al. 2000, Sánchez-Lavega et al. 2006, Vasavada et al. 2006. Both data are given in terms of the Rossby number, Ro = u/(Ωr o ), on the left axis where the length scale is the planetary radius and in m/s on the right axis. The peak strength of the simulation zonal flow in the equatorial region is larger than the Saturn's value by about 65%. The latitudinal profile of the broad, eastward equatorial zonal jet is similar to Saturn's equatorial jet. The flows in the first westward jet in the simulation are faster as well (by about a factor of two) and peak slightly closer to the equator than the similar jet on Saturn. Prior studies (Aurnou and Olson 2001, Christensen 2001, Gastine et al. 2014a, Heimpel et al. 2005, Jones and Kuzanyan 2009, Yadav and Bloxham 2020 investigating the origin of zonal flows in deep spherical shells show that the location of the first westward zonal jet is largely governed by the location of the lower boundary of the spherical shell, with deeper shells promoting westward jets at higher latitudes. In our simulation, the location of the top of the SSL acts as the lower barrier (or a pseudo 'boundary') for the outer hydrodynamic layer. Therefore, one may expect that if the SSL was located slightly deeper, then the westward jet would shift to higher latitudes and would be more in agreement with the Saturn's values; however, it should be mentioned that a deeper shell depth might also promote a stronger or weaker equatorial jet.
The two off-equatorial eastward jets peaking at around 45 • north and south are reduced in strength by a factor of about five compared to the equatorial jet, similar to the observations in Saturn's northern hemisphere (Fig. 1b). At latitudes polewards of about 45 • , longitudinally coherent zonal jets are not present in our simulations, although a longitudinal averaging does produce a tiny eastward peak at around 55 • . For latitudes within about 15 • of both rotational poles, we see some zonal jet activity with a weak, but broad and coherent eastward zonal flow in both hemispheres (discussed later; see Fig. 3). Note that, on the other hand, Saturn's cloud deck features coherent zonal jets at all latitudes.

B. Zonal flow in the deep interior
In Fig. 1c we present the variation of the zonal flow as a function of depth. Since the fluid shell is rotating along a fixed axis and has a low Ekman number, we can expect the features to be highly influenced by the planetary spin in accordance with the Proudman-Taylor theorem (Proudman 1916, Taylor 1923. Indeed, in the hydrodynamic layer (r > 0.8r o ) overlying the SSL, the zonal flow (Fig. 1b) in the mid to low latitudes largely projects downwards along spin-aligned cylinders. At higher latitudes, the zonal flow lacks such cylindrical coherence. The strong spin-axis aligned zonal jets in such conditions are thought to be maintained by the Reynolds stresses in rotating, turbulent, convective flows (e.g., Busse 1994. FIG. 2. Radial velocity (given in Reynolds number) on the equator and at various depths; the color map is saturated at ±5000 and the red (blue) shades correspond to +ve (-ve) values.

C. Comparison of the zonal flow with earlier studies
As we briefly alluded to in the introduction, zonal jet formation has been extensively investigated with nonmagnetic convection in rapidly rotating spherical shells. It has been shown that the flow boundary conditions on the spherical shell boundaries play a deciding role. For models where both the bottom and the top boundaries are stress free (or free slip), conherent zonal jets with great strength are readily excited at all latitude when the convective flow reaches sufficiently high Reynolds numbers (Aurnou and Olson 2001, Christensen 2001, Gastine et al. 2014a, Heimpel et al. 2005, Jones and Kuzanyan 2009, Yadav and Bloxham 2020, Yadav et al. 2022). If both boundaries are assumed to be no-slip (or rigid), then it is much harder to excite coherent zonal jets (e.g., see Yadav et al. 2016). If the boundary conditions are mixed, i.e. no-slip condition at the bottom and freeslip condition at the top, then only one zonal jet outside the tangent cylinder (a spin-aligned cylinder touching the inner boundary) is formed (Aurnou and Heimpel 2004, Heimpel et al. 2022, Jones and Kuzanyan 2009). The quenching of zonal flow in the presence of rigid boundaries is largely due to the enhanced friction provided by the Ekman boundary layers (e.g., see Jones 2007). In the absence of Ekman boundary layers, the zonal jets can saturate at much larger strengths, which is only limited by the much lower volume friction (Christensen 2002, Jones 2007). This also helps explain the case for mixed boundaries: any local cylindrical section of the fluid domain outside the tangent cylinder only experiences freeslip boundaries on its northern and southern ends, which allows strong zonal flows. For flows inside the tangent cylinder, the bottom no-slip boundary provides enough friction to damp coherent zonal flows.
The zonal jet features mentioned for the mixed boundary conditions have also been found in some dynamo models where electrically conducting fluid in the deeper part of the shell ('dynamo' region) and low conductivity fluid in the shallower layer ('hydro' region) are simulated simultaneously. If these dynamo simulations produce strong dipole-dominant magnetic fields, then the Lorentz force effectively renders the dynamo region into a rigidly rotating shell which provides a drag 'boundary condition' (mimicking some aspects of a no-slip boundary condition) for the fluid in the shallower low-conductivity region on top. A magnetic tangent cylinder is setup (Dietrich and Jones 2018) at the interface of the dynamo and the hydrodynamic region which determines the width of the broad, eastward equatorial jet (Duarte et al. 2013, Gastine et al. 2014b, Yadav et al. 2022. In our simulation, some of the zonal jets present in the outer convective shell manage to penetrate significantly into the stably stratified layer (Fig. 1c). Remarkably, the strong westward jet flanking the equatorial jet manages to pierce through the SSL almost cylindrically. The adjacent higher-latitude, eastward jets in the north and in the south penetrate to a lesser degree, are less cylindrical in the SSL, and join with a much weaker flow at the equator. Note that since the SSL provides only a radial restoring force, it can indeed be expected that zonal flow may penetrate into the SSL. The recent study of Gastine and Wicht (2021) shows much less penetration of zonal flows into the SSL. We attribute this difference to the less 'stiff' SSL with lower B-V frequency used in our study (see Model Setup subsection).
The convective dynamo region (r < 0.62r o ) also features zonal flows which are much weaker in strength but are maintained steadily. Such zonal flows with a broad westward drift in the mid to low latitudes and an eastward drift in the higher latitudes are commonly produced in the geodynamo simulations and are maintained by a thermal wind balance (Aubert 2005). The weak westward flow from the dynamo region penetrates into the lower half of the SSL and is contained in the mid to lower latitudes. Note that a weak westward flow is also present in the dynamo region of Gastine and Wicht (2021).

D. Radial flow features
The radial flows present in the simulation are shown in Fig. 2. The SSL featuring much weaker radial flow (with much fainter colors) is immediately evident. The radial flows within the SSL are predominantly internal gravity waves (IGWs). The flows in the inner dynamo region, showcasing spin-aligned, radially-stretched convective plumes, appear very much like those seen in low Ekman number geodynamo simulations (e.g., see Yadav et al. 2016).

E. Zonal jet truncation
Using a detailed force balance analysis, Gastine and Wicht (2021) showed that the primary mechanism that limits the strongest zonal flows to the outer low conductivity layers on top of the SSL is similar to the one outlined by the theoretical model of Christensen et al. (2020). As the zonal winds reach fluid with sufficiently high conductivity, they result in, via induction, Lorentz forces that drive a meridional circulation pattern. In the SSL, this circulation modifies the background stratification and the entropy profile, promoting a thermal wind balance which acts to truncate the zonal winds. Given the broad similarities between our model setup and the results we obtain with the one in Gastine and Wicht (2021), it is likely that the zonal flow truncation follows a similar mechanism.

F. Large cyclones and anticyclones
Before we move on to describe the properties of the dynamo generated magnetic field, we highlight the rich dynamics present in the outer layers of the simulation. Figure 3a shows the equatorial jet, and the westward and eastward jets (present in both hemispheres) at higher latitudes. The edges of these jets are cyclonic and anticylonic shear regions which promote and sustain strings of cyclones and anticyclones. We recently demonstrated such a mechanism of vortex generation in spherical shell with convectively driven flows (Yadav et al. 2022). This mechanism is known to exists on both Saturn and Jupiter (e.g., see Vasavada and Showman 2005). Most of these large scale vortices in the simulation are deep, spinaligned structures existing in the low conductivity region on top of the SSL. Although we do not have a general overview of how deep are the vortices on giant planets, a recent orbit of Juno on Jupiter was close enough to show that at least some of the vortices could be 100s of kilometers deep (Bolton et al. 2021).
Note that the current simulation lacks the preferential production of giant anticyclones as reported in our earlier study (Yadav et al. 2022). As we hypothesized in that study, large plumes from a deep dynamo region are needed to excite anticyclonic activity in the overlying hydrodynamic layer. In the simulation we report here, such a mechanism is not possible since the thick SSL on top of the dynamo region does not allow any plume to pass through it. This may help explain why, Jupiter, which may only have a thin SSL (e.g. see Debras andChabrier 2019, Wahl et al. 2017), possibly allowing some plumes to go through, has an anticyclonic great red spot and a high preference for anticyclones, while Saturn lacks such preference. The region from approximately 45 • -75 • north (south) contain no coherent jet activity and the vortices that form do not have a preferred rotation sense (see the animation of the Fig. 3 provided in the SI). In the polar region (Fig. 3b), however, a coherent cyclone is present with an accompanying serpentine eastward jet (similarly in the south). The polar cyclone is not stationary and meanders about the pole. Such polar cyclones were recently reported in convectively driven spherical shell models (Garcia et al. 2020, Heimpel et al. 2022, Yadav and Bloxham 2020. The serpentine jet visible in the simulation's polar region is somewhat reminiscent of the polygonal jet, shaped as a hexagon, on Saturn (Godfrey 1988). However, it lacks the remarkable coherence seen on Saturn and in our previous study (Yadav and Bloxham 2020) where we reported polygonal jets in a convectively driven thin spherical shell.

G. Dipole dominant magnetic field
Let us now turn to the properties of the dynamo generated magnetic field in the simulation. The radial component of the magnetic field at various depths is shown in Fig. 4a. On the surface of the simulation, the magnetic field has very smooth character and is highly dipoledominant (carrying about 95% of the magnetic energy). At a depth in the lower half of the outer molecular layer, some small scale features start appearing due to the increasing electrical conductivity of the fluid (Fig. 1a). The magnetic field is still dipole-dominant. As we reach the middle of the SSL, where the electrical conductivity has saturated to its maximum value, large deviations in morphology compared to the surface develop; here, the dipole carries only about 2% of the magnetic energy. The most noticeable deviation being the depletion of radial magnetic field inside the tangent cylinder to the inner boundary at r i . At radial depths reaching the top of the deep dynamo, where electrical conductivity is high and convection is taking place, the magnetic field is dominated by small scale structure. The magnetic field has completely changed in character and the dipole mode, which was dominant near the surface, now carries only about 0.5% of the magnetic energy. A breakdown of the poloidal magnetic energy across various length scales and at different radii is provided in the Fig. 7. Remarkably, the radial magnetic field inside the tangent cylinder has reversed in direction compared to that at the surface. Such weakening of radial magnetic field inside the tangent cylinder has been seen in some geodynamo simulations (see a recent investigation and discussion presented in . The axisymmetric toroidal and poloidal magnetic field shown in Fig. 4b and 4c provide further insights into the magnetic field morphology. The poloidal field lines demonstrate the transition from a complex field inside the dynamo region to a largely dipole-dominant one in the outer layers. It is evident that the reversal of radial magnetic field inside the tangent cylinder is a consequence of the poloidal field lines getting pulled back into the dynamo region as opposed to connecting to the other hemisphere. The counter rotating magnetic loop established inside the tangent cylinder is likely a consequence of the weak, but coherent, eastward zonal flow present inside the tangent cylinder and the associated meridional circulation ).

H. Rings of magnetic energy
In terms of the toroidal magnetic field shown in Fig. 4b  and 4c, the most remarkable features are found in the SSL. At low latitude, each hemisphere contains a strong toroidal magnetic field region, which is followed by another broader and much stronger one at higher latitudes. In Fig. 5a we present an orthographic view to elucidate where the magnetic energy is concentrated. Not surprisingly, the deep dynamo region (r < 0.62r o ) stands out as a central sphere carrying most of the magnetic energy. Moreover, two distinct magnetic energy rings hovering above the deep dynamo, one in north and another in south, are clearly visible. The field lines shown in Fig. 5b clearly demonstrate that these magnetic energy rings are composed primarily of toroidal field lines that are wrapped up in the azimuthal direction. Note that the two patches visible near the equator in Fig. 4c are significantly weaker, and, therefor, are filtered out by the transparency function adopted in Fig. 5a. Azimuthally wrapped magnetic field structures were previously reported in a thin shell dynamo simulation of a Sun-like star (Brown et al. 2010). There, such structures were not as azimuthally coherent and were part of the dynamo-wave process, unlike in our study where they form outside of the convective dynamo region and are relativley stable.
The local zonal flow shown on a meridional slice in Fig 5a demonstrates that the location of the magnetic energy rings roughly coincides with the termination of the first off-equatorial eastward jet in both hemispheres. This gives us a hint that the interaction of the jets with the electrically conducting fluid in the SSL is likely producing strong shear which generates the toroidal magnetic field through the classical Ω-effect (Krause and Rädler 1980). Such a mechanism essentially bends poloidal magnetic fields into toroidal morphology at the expense of the kinetic energy stored in the zonal flows. Following (Brown et al. 2011), we quantify the Ω-effect as: where the overbar denotes azimuthal averaging. In Fig. 5c, we plot the axisymmetric component of λ(r)Ω on a meridional plane. In the deep dynamo layer, the contours of the toroidal magnetic field have no discernible correlation with the Ω-effect color map. In the SSL, however, the color map almost perfectly aligns with the contours in the mid to low latitude regions. This strong correlation indicates that the magnetic energy rings present in the SSL are due to the Ω-effect. Therefore, the production of the toroidal magnetic field in the SSL is intimately linked to the termination of the zonal jets in the SSL. Another way to inspect the interaction of the zonal jets and magnetic field is to look at their behavior at a certain depth. In Figure 6, we plot the axisymmetric components of the zonal flow, the radial magnetic field, and the azimuthal magnetic field at the top of the SSL (0.8r o ). The profile of the radial magnetic field (Fig. 6b) and the zonal flow (Fig. 6a) at this depth motivates us to look at smaller scale fluctuation in the radial field. field, we separate the smaller scale fluctuations and compare it to the azimuthal magnetic field at the same depth (Fig. 6d). It is immediately evident that there are clear correlations among these magnetic field components and the zonal flow; for instance, the first off-equatorial peaks are located around 25 • north and south for zonal flow as well as the zonal magnetic field and the smaller scaled radial magnetic field. Similar correlations were envisioned by Cao and Stevenson (2017). Note that the peaks in the azimuthal magnetic field and the smaller scaled radial magnetic field are shifted by about a quarter of the phase. Such an interaction was modelled by Cao and Stevenson (2017) using mean-field electrodynamics in the context of jet-magnetic field interaction at Jupiter and Saturn; also see Galanti and Kaspi (2021) for application of this idea for constraining the depth of the winds in the two gas giants.

I. Comparing the magnetic field spectrum with Cassini observations
The magnetic power spectrum displayed in Fig. 8a shows that the spherical harmonic degrees from 1 through 7 carry most of the energy in their axisymmetric component. For degrees eight and higher, the behavior is reversed. The magnetic field spectrum compares favorably ( Fig. 8b) with the one inferred for Saturn based on the Cassini Grand Finale magnetic field measurements (Cao et al. 2020, Dougherty et al. 2018. A good agreement with Saturn's magnetic field up to spherical harmonic degree eight, spanning almost six orders in magnitude of the magnetic power, is remarkable. Beyond degree eight, the simulation generally contains more magnetic energy at  (Cao et al. 2020). Note that both spectra are normalized by the power in the dipole mode. different length scales compared to Cassini observations. Note that in Fig. 8b we have normalized the spectra so that the dipole power is at unity for comparison. Due to the mismatch in control parameters of the simulation and those at Saturn (for instance, much higher fluid viscosity in the simulation), it is not entirely straightforward to compare the magnetic field magnitudes.
We must note here that the spectrum of the simulation shown in Fig. 8 is simply a time shot and the similarity with Saturn's spectrum might not endure. As might be expected for a dynamical non-linear system, the spectrum may evolve into other states in which the characteristic dip at harmonic degree five goes away and a peak could appear instead; the latter is a common feature in dipole-dominant dynamo simulations (e.g., see Duarte et al. 2013). For Saturn, since we do not have a long record of its magnetic field morphology, its magnetic field spectrum might have been significantly different at different epochs.

IV. SUMMARY & DISCUSSION
The rich phenomenology of the fluid dynamics and magnetic field observed at Saturn has motivated many to propose novel driving mechanisms. To explain the exceptional level of axisymmetry in the Saturn's dynamo, Stevenson (1979) proposed the presence of a stably stratified layer (SSL) in which zonal flows may be able to reduce the non-axisymmetric components of a dynamo generated field deep within the planet. We set out to investigate how such a layer would impact the overall dynamics of the planet. Assuming a spherical shell geometry, we model a convective dynamo layer from 0.27 to 0.62, an electrically conducting and stably stratified layer form 0.62 to 0.8, and a convective, low-conductivity outermost layer from 0.8 to 1 (Fig. 1a).
The outermost layer self-consistently generates coherent, alternating zonal flows that reach the mid-latitudes ( Fig. 1b and 1c). These strong jets spawn large cyclonic or anticyclonic vortices depending on the local shear direction in the region where jets change direction (Fig. 2a). In higher latitudes, the zonal jet activity is diminished and is replaced by (anti)cyclonic vortices with no clear preference (Fig. 2b). A cyclonic polar storm exists on both rotation poles. We attribute this unprecedentedly rich dynamics to the low Ekman number (7 × 10 −7 ) and high level of turbulence (Reynolds numbers reaching 5000 or higher) achieved in the simulation. The multiple zonal jets are largely confined to the low conductivity layer above the SSL.
The deep dynamo generates a magnetic field with the octupole harmonic degree ( = 3) mode being the dominant carrier of energy within the dynamo region. However, after passing through the SSL and the low conductivity outer layers to the surface, the magnetic field morphology becomes strongly dipole dominant ( Fig. 4 and Fig. 7). Typically, the dipole mode is also tilted with a small angle of about 0.5 degrees with respect to the rotation axis.
The eastward oriented zonal jet situated around midlatitudes on the surface, projects cylindrically downwards along the direction of the rotation axis and interacts with the magnetic field in the outer half of the SSL. This leads to an efficient local Ω-effect where the poloidal magnetic field gets converted to toroidal morphology. The resultant toroidal magnetic field stands out in the magnetic energy space as two rings (one in each hemisphere) surrounding the deep dynamo region (Fig. 5).
At certain epochs in the simulation, the magnetic field spectrum on the surface compares favorably, for spherical harmonic degrees eight or less, with that observed at the present-day Saturn. On smaller length scales, the simulation produces more magnetic energy than the present-day Saturn (Fig. 8).
Although the features described above do not precisely match those observed at Saturn, the several qualitative similarities we obtain is made remarkable by the fact that they are all produced in one single simulation without any special care given to tune the results. Therefore, the idea of a thick SSL inside Saturn just above the deep dynamo appears to take us very close to many of Saturn's observed features.
We can certainly imagine a few directions in which the current model can be improved. One way to 'hide' magnetic energy from the surface is to convert the poloidal magnetic field into toroidal morphology which can not escape the outer non-conducting fluid layers. Our simulation lacks coherent zonal flow from mid to high latitudes. Since low Ekman numbers promote strong, coherent zonal flows (Cabanes et al. 2017), we can expect a similar simulation at even lower Ekman number to produce coherent zonal jets at even higher latitudes. These zonal flows might help to convert more of the poloidal energy to toroidal energy at various length scales, helping to put the simulation more in line with the Saturn's internal magnetic field for spherical harmonic degrees higher than eight.
The dipole tilt angle in the simulation is also not as small as the one on Saturn. We recently showed that dynamos with moderate Ekman numbers and low convective supercriticalities are able to produce dipole tilts even smaller than Saturn (Yadav et al. 2022). However, the model of Yadav et al. (2022) lacked many other properties of Saturn, e.g. multiple zonal jets were not present, there were no large scale vortices near the surface, and the magnetic energy in the quadrupolar components was much too low as compared to Saturn.
Based on the results and discussion above, we believe that the introduction of the SSL is an additional and crucial ingredient needed to promote more Saturn-like features. The presence of an SSL partially separates the outer hydrodynamic layer from the inner dynamo. This is likely needed for promoting multiple off-equatorial zonal jets in outer layers since earlier models without an SSL only produced one equatorial jet outside the magnetic tangent cylinder (Dietrich and Jones 2018, Duarte et al. 2013, Yadav et al. 2022. The SSL also provides a relatively quiescent region where zonal jets, impinging from the outer hydrodynamic layer, can interact with the magnetic field and selectively suppress some of the poloidal field depending on the zonal flow strength and length scale. This, we believe, appears crucial for inducing a selective dip in the 5th spherical harmonic degree component of the poloidal magnetic field. Therefore, future improvements will likely need to combine aspects of both studies. For instance, the parameter regime where dynamos with extremely small dipole tilt exist are likely broader at lower Ekman numbers (Yadav et al. 2022). There might be a parameter space where the dynamo region in the deep interior itself generates a dipole-dominant solution with very small dipole tilts, and the SSL and zonal jets in the outer layers help further reduce the non-axisymmetric features to minuscule values. With a broader parameter study in future, we hope to answer some of these questions. Acknowledgements: The work was supported by the NASA Cassini Data Analysis Program (Grant No. 80NSSC21K1128). The computing resources were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and the Research Computing, Faculty of Arts & Sciences, Harvard University. Some of the intermediate simulations were carried out using the computing support provided by the NASA Juno project. RKY thanks Dr. Ankit Barik for helping to extract the background pressure profile of the simulation.

Data
Availability: The simulation input file that can be used to reproduce the results is archived on the Harvard Dataverse: "https://doi.org/10.7910/DVN/QTHFVS".
The simulation code used is open access and is available here: "https://github.com/magic-sph/magic/". The large sized output data, which is not practical to be archived on an open access repository, is available upon request to the authors.

LEGENDS FOR THE SI ANIMATIONS
SI Animation 1: The animations show the evolution of the trajectories of massless particles driven by the local latitudinal-longitudinal flow components (see Fig. 3). The left panel shows a global overview of the simulation from a northern mid-latitude viewpoint. The right panel views the same animation from the north pole. The trajectories are visualized at radius 0.93r o . The animation, lasting about 3 seconds, shows roughly 12 rotations of the simulation in physical time. The highly non-linear and complex evolution of the zonal jets and vortices is clearly visible in the animation.
SI Animation 2: The animation (lasting for about 1m6s) shows the simulation from the viewpoint of a hypothetical orbit around it. The left panel shows the instantaneous magnetic energy as volumetric rendering, along with a meridional slice of the instantaneous local zonal flow as color map (see Fig. 5a). We have made the smaller magnitude magnetic field either completely transparent or translucent in order to highlight the regions with intense magnetic field. The two rings of azimuthal magnetic energy as well as the central dynamo region are made clearly visible in this way. The right panel shows the instantaneous stream lines of the magnetic field (see Fig. 5b). The grey spherical surface is located at radius 0.7r o . The rings of magnetic energy visible in the left panel are now visible as magnetic field lines which are mainly wrapped around the sphere.