Ocean Circulation on Enceladus With a High Versus Low Salinity Ocean

Previous studies that have considered the ocean circulation on Enceladus have generally assumed the salinity to be Earth-like. However, according to observations and geochemical constraints, the salinity of Enceladus' ocean is likely to be lower, and importantly, it is probably low enough to reverse the sign of thermal expansivity. We investigate the ocean circulation and stratification of Enceladus' ocean using a combination of theoretical arguments and simulations using the MITgcm. We find that, if the salinity is high, the whole ocean is unstratified, and convection dominates the entire ocean. However, if the salinity is low enough, there exists a stratified layer in the upper ocean, whose thickness depends on the magnitude of the turbulent vertical diffusivity, which remains poorly constrained. Such a layer can suppress the vertical flux of heat and tracers, thereby affecting the heat flux to the ice shell and leading to a vertical tracer mixing time scale across the stratified layer of at least hundreds of years. This time scale is inconsistent with a previous estimate of vertical ocean mixing of several years, based on the size of detected silica nanoparticles in the plumes, leading us to conclude that either the salinity of Enceladus' ocean is higher than previously suggested or the interpretation of silica nanoparticle observations has to be reconsidered.


INTRODUCTION
Strong evidence suggests that Enceladus maintains a global ocean (e.g., Postberg et al. 2011;Patthoff & Kattenhorn 2011;Thomas et al. 2016), and the possible existence of liquid water in contact with a rocky interior makes it a hot target in the search for life in the solar system. Material from Enceladus' ocean is continuously ejected into space, where some of it forms Saturn's E-ring (Schmidt et al. 2008;Kempf et al. 2008), which makes it the only extraterrestrial ocean that is known to be so accessible for sampling. Understanding the ocean on Enceladus can also assist in understanding oceans on other icy moons in and outside the solar system, and help us better predict their habitability.
The ocean on Enceladus is estimated to be about 40 km in depth on average, and covered by a global ice shell (Thomas et al. 2016). The ice shell is about 20 km deep on average, with the thickest part at the equator estimated to be more than 30 km and the thinnest part at the south pole less than 10 km (Beuthe et al. 2016;Čadek et al. 2019;Hemingway & Mittal 2019). To maintain such a global ocean under the ice shell, as well as to explain the heat loss rate of around 10 GW in the south polar region (Spencer et al. 2006;Howett et al. 2011;Spencer et al. 2013), there must be an energy source inside Enceladus.
The energy source is likely to be associated primarily with tidal dissipation. In general, tidal dissipation is expected to occur in the ice shell, in the ocean, and in the solid core. However, tidal heating in the ocean is believed to be negligible compared to that in the ice shell and inner solid core (Chen & Nimmo 2011;Tyler 2011;Beuthe 2016;Hay & Matsuyama 2017). Libration can also generate heat in the ocean and may result in a total heating of up to O(0.1 GW) (e.g., Wilson & Kerswell 2018;Rekier et al. 2019;Soderlund et al. 2020), which, however, is still smaller than the estimated tidal dissipation rate in the ice shell and the solid core. The total tidal dissipation rate in the ice shell (not including dissipation within liquid-water conduits) has recently been suggested to be on the order of 1 GW, with the maximum at the south pole where the ice shell is thinnest Beuthe 2019). Tidal dissipation in the ice shell has also been suggested as an explanation for the north-south asymmetry (Kang & Flierl 2020) and equator-to-pole variations  in the ice shell thickness, as well as for the sustained plumes associated with the south-polar tiger stripes (Kite & Rubin 2016). Continuous high-temperature hydrothermal activity suggests vigorous tidal heating in the solid core Hsu et al. 2015). The tidal dissipation rate in the core is likely to reach O(10 GW), and is believed to be strongest at the pole and weakest at the substellar and anti-substellar point at the equator (Choblet et al. 2017). Tidal energy dissipation in the solid core will generate heat, which must be transported outwards to the ice shell by the ocean circulation. Although the fraction of tidal heating coming from the solid core and the ice shell remains highly uncertain, some heating from the bottom solid core is expected, which fundamentally shapes the ocean circulation on Enceladus.
Previous studies have looked into possible scenarios for the ocean circulation on Enceladus and other icy moons, with both ocean-only models and ice-ocean coupled models (e.g., Soderlund et al. 2014;Travis & Schubert 2015;Soderlund 2019;Amit et al. 2020;Ashkenazy & Tziperman 2020;Kang et al. , 2021. One characteristic feature of the ocean circulation on Enceladus may be hydrothermal convection columns, which are expected to be aligned parallel to the rotation axis, and could extend from the sea floor to the ice-ocean interface (Goodman et al. 2004;Goodman & Lenferink 2012;Soderlund 2019;Ashkenazy & Tziperman 2020;. For another icy moon, Europa, hydrothermal convection plumes, quasi-3D turbulence, and baroclinic eddies have been suggested, based on different estimated parameter regimes, characterizing the strengths of the forcing that drives convection (Goodman et al. 2004;Goodman & Lenferink 2012;Soderlund et al. 2014), and depending on whether salinity changes from freezing and melting are taken into consideration (Ashkenazy & Tziperman 2020). Noting the substantial discrepancies among the simulations, the large uncertainty in external parameters, and the inherent computational challenges in modelling the global ocean in a realistic parameter regime, the ocean circulation regime on Enceladus remains uncertain.
All previous simulations of Enceladus' ocean assumed the salinity to be roughly similar to Earth's ocean and thus high enough that the thermal expansivity (α = −(1/ρ)(∂ρ/∂T )) is always positive. However, the ocean on Enceladus is likely to be fresher than Earth's ocean. Although some studies have suggested the salinity to be higher than 20 g kg −1 (e.g., Ingersoll & Nakajima 2016) or around 20 g kg −1 (Kang et al. 2021), based on dynamical considerations, geochemical evidence and modelling suggest that a salinity of around 20 g kg −1 , dominated by NaCl, is an upper bound (e.g., The area below the black line is where thermal expansivity of sea water is negative. The salinity at the intersection of the black and red lines is the maximum salinity for which negative α can exist under the given pressure. The calculation of the density is based on Jackett & Mcdougall (1995), and the calculation of the freezing point is based on Fofonoff & Millard Jr (1983). Postberg et al. 2009;Hsu et al. 2015;Glein et al. 2018). If the salinity is less than 20 g kg −1 , the thermal expansivity of sea water is negative near the freezing point ( Figure 1). In this case, a stably stratified layer may exist in the upper ocean on Enceladus, similar to a scenario that has been proposed for Europa by Melosh et al. (2004). To the best of our knowledge, no General Circulation Model (GCM) simulations have thus far been carried out to confirm the existence and effect of such a stably stratified layer.
Here we investigate the role of salinity on the dynamics of Enceladus' ocean, focusing in particular on the ocean's role in transporting heat and tracers from the sea floor to the ice shell. Ocean mixing can play a key role in transporting heat and constituents from the ocean-rock interface to the ice shell, thus affecting ice melting as well as properties observable in plumes. According to the detected size of silica nanoparticles and the growth rate expected from Ostwald ripening, the vertical mixing time scale has been estimated to be at most several years 3 (Hsu et al. 2015). One goal of this study is to analyze whether this vertical mixing time scale is consistent with the expected ocean circulation and transport on Enceladus. We present theoretical arguments and carry out global ocean simulations using the 3D general circulation model MITgcm. Section 2 describes theoretical predictions for the ocean circulation on Enceladus. Section 3 shows results from numerical simulations. Section 4 provides concluding remarks.

The Thermal Expansion Coefficient and Implications for Ocean Stratification
The equation of state for saltwater describes the density as a non-linear function of temperature (T ), salinity (S), and pressure (P ) (Jackett & Mcdougall 1995). For salinities and pressures found in Earth's ocean, the density increases with salinity and decreases with temperature. Assuming a heat flux from the bottom due to tidal energy dissipation in the rocky interior, convection may then be expected in the whole ocean. As a result, the ocean would likely be well-mixed with very small temperature gradients (Figure 2(a)), consistent with the simulation results of Soderlund (2019), Ashkenazy & Tziperman (2020), and .
However, for low enough temperature and salinity (and at the modest pressures expected on Enceladus), the thermal expansion coefficient, α, is negative, such that density increases with temperature. Under such circumstance, a critical temperature, T c , exists at which the density of water reaches a maximum at a given salinity and pressure (Figure 1). In this regime, heating from below will not trigger convection until the temperature is heated above T c . However,

Stratified layer
Convective layer Expected for Earth-like high salinity ocean Expected for low salinity ocean T c Figure 2. Schematic for the expected vertical structure of the ocean on Enceladus with different salinities. (a) and (b) show the vertical structure and temperature profile of high and low salinity oceans, respectively. Red lines indicate the temperature profile, and the blue dashed line in (b) indicates the critical temperature Tc, which decreases as pressure increases. The high salinity ocean is expected to be virtually unstratified, with convection throughout the ocean and a small negative vertical temperature gradient around the freezing point. In the low salinity ocean, we expect two layers: the upper stratified layer with linear vertical temperature profile, and the lower convective layer that has a weak vertical temperature gradient. The temperature at the interface between these two layers is at the critical temperature Tc, where the thermal expansion coefficient changes sign. (See Figure 3 and associated discussions for the horizontal structure of the heat flux.) the temperature at the upper boundary will be kept at the freezing point, T f , due to the surface ice shell. A stably stratified layer is thus expected to form in the upper ocean, where the thermal expansion coefficient is negative and the temperature decreases upwards from T c to T f . This is similar to the stratified layer that has been hypothesized to exist in the ocean of Europa under the low salinity assumption by Melosh et al. (2004), although we note that the lower gravity and hence lower pressures on Enceladus make the conditions more favorable for the existence of a stably stratified layer. Below the stratified layer, the temperature is above T c so that the thermal expansivity is positive. This layer is expected to be qualitatively similar to the high salinity scenario and we expect it to be characterized by convection. The temperature at the interface between these two layers should be near the critical temperature T c . A schematic for the expected vertical stratification in a low salinity ocean is shown in Figure 2(b).
Since convection cannot occur in the stratified layer, the vertical heat flux is expected to be dominated by diffusion (driven by either molecular diffusion or small-scale turbulence). We can then estimate the depth of the stratified layer using the equation of heat diffusion, Q = c p ρκ z,heat ∂T /∂z, where Q is the vertical heat flux, c p is the heat capacity, ρ is the density of the water, and κ z,heat is the vertical thermal diffusion coefficient. Given the temperature contrast across the stratified layer, ∆T = T f − T c , we can estimate the depth of the stratified layer H as If we assume the total bottom heating to be around 20 GW (Choblet et al. 2017), the mean vertical heat flux is Q ≈ 0.03 W m −2 at the top of the ocean. This value is consistent with the magnitude estimated inČadek et al.
(2019) (i.e. tens of milliwatts per squaremeter). If we further assume the salinity to be 8.5 g kg −1 for Enceladus' ocean (Glein et al. 2018), we find T f ≈ −0.6 • C, T c ≈ 1.4 • C, so that |∆T | ≈2.0 K. If we use the magnitude of the molecular thermal diffusivity, around 10 −7 m 2 s −1 , the depth of the stratified layer is H ≈ 30 m. In general, vertical mixing in the ocean can be intensified through turbulence, in which case the molecular diffusivity should be replaced by the turbulent diffusivity. If we take the magnitude of the turbulent diffusivity in Earth's ocean, around 10 −5 m 2 s −1 (Munk & Wunsch 1998), the depth of the stratified layer would be H ≈ 3 km. The question of whether and by how much vertical mixing in the ocean of Enceladus is enhanced by turbulence is hence important and is discussed in the following section.

Tidal Dissipation and Turbulent Mixing in a Stratified Ocean
Vertical mixing in a stably stratified ocean increases the potential energy of the water column as buoyancy is fluxed downward (mixing lighter water downwards and denser water upwards), and hence requires a source of energy. If we know the rate of turbulent kinetic energy dissipation in the ocean's interior, we can estimate the vertical turbulent diffusivity based on the energy required to mix a stably stratified ocean (e.g., Wunsch & Ferrari 2004;Yang et al. 2017): Here κ z is the vertical turbulent diffusivity, ε is the turbulent kinetic energy dissipation per unit volume, Γ ≈ 20% is the "mixing efficiency", i.e. the fraction of the kinetic energy dissipation that contributes to the generation of potential energy (Peltier & Caulfield 2003;Wunsch & Ferrari 2004), and N is the Brunt-Väisälä frequency. Given that N 2 = −(g/ρ θ )(∂ρ θ /∂z) ≈ αg∆T /H where ρ θ is potential density and g is gravity, we have where E is the total turbulent kinetic energy dissipation rate and A is the horizontal area, and hence ΓE/A is the energy used for vertical mixing per unit area, in units of W m −2 . On Enceladus, the turbulent kinetic energy that drives vertical mixing of the stably stratified water column can be derived from tidal disspation and libration. The tidal dissipation in the ocean on Enceladus is likely to be relatively weak, and is suggested to be 10 1 -10 4 W, based on different models (Chen et al. 2014;Matsuyama et al. 2018;Hay & Matsuyama 2019). Libration can result in significant dissipation up to O(0.1 GW) in the ocean (Wilson & Kerswell 2018). The dissipation caused by libration is likely to be concentrated primarily in the surface Ekman boundary layer (e.g., Greenspan 1968;Liao & Zhang 2008;Cébron et al. 2012;Lemasquerier et al. 2017). If mixing is confined to a thin boundary layer under the ice shell, it is expected to have little effect on the stratified layer underneath. However, some energy may be transferred into internal waves which can break and contribute to mixing in the interior (e.g., Werne & Fritts 1999;Wunsch & Ferrari 2004). Whether the libration-driven dissipation is concentrated in the surface boundary layer, or can significantly contribute to interior mixing, remains uncertain. We here treat the value of 0.1 GW as an upper-bound limit for the interior turbulent energy dissipation that contributes to mixing in the stratified layer, although we consider a much smaller value to be more likely.
We therefore estimate that the total turbulent energy dissipation rate that contributes to mixing of the stratified layer may be anywhere in the range of E ∼ 10 1 − 10 8 W, and thus the energy input to vertical mixing per unit area, ΓE/A, is around 3 × 10 −12 to 3 × 10 −5 W m −2 . Following the estimate in Section 2.1 for the conditions in the stratified layer of the low salinity ocean, we get ∆T ≈ −2.0 K and α ≈ −4 × 10 −5 K −1 , where we assumed a pressure of 20 bar and a temperature near the freezing point. We can then estimate κ z to be around 3 × 10 −10 to 3 × 10 −3 m 2 s −1 . For the upper limit (3×10 −3 m 2 s −1 ), the depth of the stratified layer would be expected to be around 800 km according to Equation (1), which is deeper than the whole ocean. This indicates that if libration can drive strong vertical turbulent mixing in the interior, the whole ocean may become stably stratified, as seen in the low salinity simulations of Kang et al. (2021) where a vertical diffusivity of 5 × 10 −3 m 2 s −1 is assumed. For the lower limit, the turbulent diffusivity is weaker than the molecular value (10 −7 m 2 s −1 for thermal diffusion and 10 −9 m 2 s −1 for tracer diffusion), indicating that turbulence would be unable to significantly enhance vertical mixing in the interior.

Conceptual Models of Ocean Heat Transport
To illustrate the role of ocean dynamics for the transport of heat from the rocky core to the ice shell, we here discuss conceptual models of ocean heat transport based on possible directions of anisotropy in mixing brought about by rotation and gravity. Anisotropy brought about by gravity tends to align radially, while anisotropy due to rotation tends to align along the rotation axis (like Taylor columns). We therefore suggest four conceptual models in which heat is transported along/normal to the direction of rotation/gravity.  In the first model "Parallel to Rotation", we assume heat is transported along the axis of rotation, representing slantwise convection (Figure 3(a)). On Enceladus, convection is likely to be parallel to the rotation axis-or more generally along surfaces of constant total angular momentum (Goodman et al. 2004;Goodman & Lenferink 2012;Soderlund 2019;. In this limit, the heat flux inside the tangent cylinder that encircles the solid core (the blue dashed lines in Figure 3(a)), can be calculated by matching the bottom heat flux to the corresponding area at the ice-ocean interface. At the intersection of the tangent cylinder and the ocean surface, the heat flux to the ice-ocean interface diverges in this model, because heat from a finite bottom area is transported to an infinitesimal surface area. Outside the tangent cylinder (at low-latitudes), there is no heat flux at the surface of the ocean since slantwise convective columns, following the axis of rotation, cannot reach the surface there (Figure 3(e)).

Parallel to Rotation
In the second model "Perpendicular to Rotation", we consider the limit case where all heat is transported in columns perpendicular to the rotation axis (Figure 3(b)). This model represents the lateral heat transport by Taylor columns, which may be dominant on icy moons, especially at low latitudes (e.g., Cardin & Olson 1994;Christensen 2002;. Opposite to the first model, there is no heat flux at high latitudes while the heat flux at low latitudes is calculated by matching the bottom heat flux to the corresponding area at the ice-ocean interface (Figure 3(f)). To the best of our knowledge, a vanishingly small heat flux to the ice shell at the equator or the pole has not been observed in any prior simulations of icy moon oceans, although many previous simulations do show differing heat fluxes at low-versus high-latitudes (e.g., Soderlund 2019; Amit et al. 2020;, with strongly differing patterns resulting from different assumptions for various parameters. Specifically, the results of Amit et al. (2020) suggest that heat flux parallel to the axis of rotation may be dominant at relatively large Rossby numbers, while heat flux perpendicular to the axis of rotation may be expected to dominate at relatively small Rossby numbers. In general, we may not expect one single model/mechanism to fully explain the heat transport process, but, depending on parameter assumptions, some combination of these models may be applicable.
In the third model, "Parallel to Gravity", we assume that heat is transported to the surface radially. This model is expected to be relevant if heat transport is dominated either by radial convection, or by 3-dimensional turbulence or molecular diffusion, which generate an approximately isotropic diffusive transport. If the depth of the ocean is relatively small compared to the horizontal scale over which temperature varies, the heat flux in isotropic diffusion is expected to be dominated by the radial component. In this limit, the pattern of the heat flux at the top of the ocean will match the bottom heat flux, reduced by a constant factor to account for the differences between the area at the bottom and the surface (Figure 3(c) & (g)).
In the fourth model "Perpendicular to Gravity", we assume that heat is well mixed horizontally (normal to the direction of gravity). If horizontal mixing is very efficient (e.g. due to the presence of geostrophic turbulence), any gradients in the deep ocean heating should be homogenized by ocean mixing and a relatively uniform surface heat flux is expected (Figure 3(d)). Under this circumstance, the surface heat flux should be equal to the average of the bottom heat flux (reduced by a constant factor) everywhere (Figure 3(h)).
The heat flux to the surface ice shell is likely to play a role in determining the ice thickness distribution. With more heat transported to the surface at the pole than at the equator (the "Parallel to Rotation" and "Parallel to Gravity" models), ice is likely to be melting at the pole but forming at the equator. This heat flux pattern is therefore qualitatively consistent with the observed ice thickness distribution, and is also qualitatively similar to the heat flux in the low Ekman number simulation of Soderlund (2019), and the relatively large Rossby number simulations of Amit et al. (2020). With strong horizontal ocean heat transport (the "Perpendicular to Gravity" model), the thickness of the ice shell may be expected to be uniform (Ashkenazy et al. 2018). This heat flux pattern is qualitatively similar to that found in the simulations of Europa's ocean in Ashkenazy & Tziperman (2020). The "Perpendicular to Rotation" and the "Perpendicular to Gravity" model cannot explain the observation that the ice shell is thicker at the pole than the equator on Enceladus (Čadek et al. 2019;Hemingway & Mittal 2019), and would require the ice thickness distribution to be shaped by tidal dissipation in the ice shell itself (e.g., Kang & Flierl 2020;).

Vertical Tracer Mixing Time Scale
In the high salinity regime (or the convective layer in the low salinity regime), geochemical tracers and small particles can be transported from the bottom to the surface of the ocean (or the top of the convective layer) through convective processes relatively effectively. The vertical velocity in a rotating plume can be estimated as w ∼ B/f , where B = gαQ/ρc p is the buoyancy flux and f ≈ 1 × 10 −4 s −1 is the Coriolis parameter (Jones & Marshall 1993). The convective mixing time scale is then estimated to be where D ≈ 40 km is estimated as the depth of the ocean. If we assume the heat flux Q ≈ 0.03 W m −2 and thermal expansivity α to be around 10 −6 − 10 −4 K −1 , we find τ conv ≈ 50 − 500 years. For an ocean with Earth-like salinity (35 g kg −1 ), the thermal expansivity is relatively large and the lower limit (tens of years) is more applicable. However, in the low salinity regime, there is no convection in the stratified layer, where vertical transport of tracers is achieved mainly by diffusion (either turbulent or molecular). The diffusive time scale in the stratified layer can be estimated by the vertical diffusion equation: where C is the tracer concentration and κ z,tracer is the vertical diffusivity of the tracer. Through scaling analysis of Equation (5), we find the diffusive mixing time scale where H is given by the depth of the stratified layer as predicted by Equation (1) (if it is smaller than the ocean depth) or the ocean depth if the stratified layer occupies the whole ocean. We choose Q ≈ 0.03 W m −2 and |∆T | ≈2.0 K for a salinity of 8.5 g kg −1 , following Section 2.1. Using the range for thermal and tracer vertical diffusivity estimated in Section 2.2, we find the maximum mixing time scale to be τ max dif f ≈ 3.5 × 10 5 years when the turbulent diffusivity is κ z,heat = κ z,tracer ≈ 1.4 × 10 −4 m 2 s −1 such that the depth of the stratified layer is the same as the depth of the ocean, and the minimum time scale to be τ min dif f ≈ 250 years when κ z,heat = κ z,tracer = 10 −7 m 2 s −1 , i.e. the tracer diffusivity is enhanced by turbulent mixing to a similar value as the molecular thermal diffusivity. Vertical mixing through the stratified layer is therefore expected to take at least hundreds, and possibly up to hundreds of thousands of years. Notice that there may also be a buoyancy-driven circulation in the stratified layer (due to inhomogeneities in the ice shell thickness and freshwater/salinity fluxes from melting/freezing of the ice shell; e.g., Ashkenazy & Tziperman 2020; Lobo et al. 2021;Kang et al. 2021), which results in advective transport of the tracer. However, any such circulation would be diffusively controlled (i.e., a balance exists between vertical advection and diffusion), such that the advective time scale remains constrained by the diffusive time scale. The presence of a stratified layer in a low salinity ocean is therefore likely to greatly increase the vertical mixing time scale compared to a fully convective high salinity ocean.

Experimental Design
We perform numerical simulations using the Massachusetts Institute of Technology General Circulation Model (MITgcm) to solve the non-hydrostatic equations for a Boussinesq fluid in a rotating spherical shell, where all sphericity terms are preserved, including all components of the Coriolis force (see Adcroft et al. 2018, and APPENDIX B). Radius, rotation rate and gravity are set to be the same as Enceladus (Table 1), and vertical variation of gravity is also taken into consideration (Figure 4(a), see APPENDIX B).
At the bottom of the ocean, we apply a fixed bottom heat flux pattern, following the zonal-mean tidal forcing pattern in Choblet et al. (2017) with a total flux of 20 GW (Figure 4(b), see APPENDIX B). The bottom heat flux boundary condition is qualitatively similar to that in  and Ashkenazy & Tziperman (2020), but differs from Soderlund (2019) who applies fixed and uniform bottom and surface temperatures. The prescribed bottom heat flux is better constrained from tidal models and observations than the bottom temperature distribution and provides a key constraint on the energetics of the circulation (Jansen 2016). We apply a linear bottom drag with a drag coefficient r b ∼ 10 −4 m s −1 (such that the velocity in the bottom layer is relaxed to zero with a rate of 10 −7 s −1 ). This value is likely to be unrealistically large, but we found a relatively large value to be needed to avoid slow spin-up of unrealistic mean flows. Our choice of a linear bottom drag is qualitatively similar to  (although they apply an even larger value of r b = 2 × 10 −3 m s −1 ), and is also practically similar to Ashkenazy & Tziperman (2020) who apply a no-slip bottom boundary condition together with vertical viscosity, but is different from Soderlund (2019) who applies stress-free boundary conditions. We believe that some bottom drag is preferable as it provides a physical constraint on the bottom zonal flows. At the top model layer, we apply a linear restoring of temperature towards the freezing point T f (−2 • C for salinity of 35 g kg −1 and −0.6 • C for salinity of 8.5 g kg −1 ) with a restoring time of 1 month (30 days). A 30 day restoring time scale is fast compared to the typical advective time scale, so that this boundary condition effectively amounts to a fixed temperature boundary condition. Our upper boundary condition can implicitly capture the heat exchange between the surface global ice shell and the ocean, but the fresh water flux and brine rejection associated with melting and freezing are not included in our simulations. Salinity changes associated with freezing and melting at the surface of the ocean can change the density of sea water, thus further affecting ocean circulation. A pole-to-equator overturning circulation driven by salinity gradients, which in turn are caused by ice formation at low-latitudes and melting at the pole, has been suggested in Lobo et al. (2021) and has also been found in coupled ice-ocean simulations using the MITgcm (Ashkenazy & Tziperman 2020;Kang et al. , 2021. Although freezing and melting at the bottom of the ice shell is likely to affect the ocean circulation, it is important to note that such forcing cannot energetically drive a circulation unless melting (which reduces the density) happens at a higher pressure (i.e. greater depth) than freezing (e.g. Wunsch & Ferrari 2004). Instead, assuming that the ice shell is in equilibrium, any melting or freezing of the ice 4 must be balanced by the viscous ice flow (e.g., Lefevre et al. 2014;Corlies et al. 2017;Ashkenazy et al. 2018;Cadek et al. 2019), which tends to advect ice from thicker regions to thinner regions. As a result, melting needs to   The salinity forcing then cannot drive a circulation but instead acts to stabilize the stratification in regions where the ice shell is relatively thin. The energy source for the global overturning circulation in Lobo et al. (2021), and the "salt-driven" circulation in the simulations of Ashkenazy & Tziperman (2020), , and Kang et al. (2021), comes from parameterized turbulent vertical mixing where the ocean is stably stratified. Whether a strong "salt-driven" circulation is possible on Enceladus thus remains an open question that directly ties in with the question of how much turbulent kinetic energy is available for vertical mixing, which is here left for future work. We simulate a 40 km deep ocean over a zonal range of 15 • with zonally periodic boundary conditions and a meridional range from 85.5 • S to 85.5 • N, with free-slip, no-normal flow conditions at the meridional boundaries. The longitudinal extent of the domain is limited to save computational resources and the poles are masked with land to ensure numerical stability. The vertical resolution is 1000 m. The horizontal resolution is 1 • in the zonal and 0.95 • in the meridional direction in most simulations, which is around 4 km × 4 km near the equatorial surface. This horizontal resolution is not sufficient to adequately resolve single convective columns. The expected minimum horizontal scale of the convective columns can be estimated using the length scale where rotation becomes important: l r ∼ B 1/2 f −3/2 ≈ 0.2 m (Jones & Marshall 1993), which is many times smaller than the horizontal resolution. This indicates that in our simulations, convection is affected by the resolution and parameterized turbulent diffusivities and viscosities, a situation that is unavoidable in global-scale simulations, where resolutions of O(0.1 m) are computationally impossible to achieve. In order to provide at least some insight into how resolution can affect the simulation results, we perform one additional simulation with a finer horizontal resolution of 0.5 • in the zonal and 0.475 • in the meridional direction.
The horizontal turbulent diffusivity applied in our models is 0.25 m 2 s −1 in our 1 • × 0.95 • simulations, chosen for numerical stability. In the high horizontal resolution simulation, this value is decreased to 0.1 m 2 s −1 following Kolmogorov scaling (see APPENDIX B). The turbulent vertical diffusivity is likely to be small in the stratified layer of the low salinity ocean, as discussed in Section 2.2. For the low salinity ocean we therefore use an anisotropic diffusion with a smaller vertical diffusivity set to κ z = 5 × 10 −5 m 2 s −1 . This value is similar to the vertical diffusivity in Earth's ocean and was chosen here to ensure numerical stability and to be able to explicitly resolve the stratified layer. The turbulent vertical diffusivity in the high salinity ocean is not constrained by the energetic argument discussed in Section 2.2, as no energy is required to mix the unstratified ocean. The horizontal and vertical scale of the grid in our simulations are of similar order (10 3 m), so that an isotropic diffusion regime is plausible. However, due to the effect of gravity and rotation, highly anisotropic turbulent diffusion may also be justified. We have therefore performed simulations with both isotropic and anisotropic diffusion for the high salinity ocean. All viscosities are set by fixing P r = 10.
We carry out 4 simulations to examine the influence of different factors on ocean circulation (Table 2). In order to examine the role of salinity, we carry out two simulations, HSaniso and LSaniso, with two different salinities: 35 g kg −1 for the high salinity case (similar to Earth's ocean) and 8.5 g kg −1 for the low salinity case (Glein et al. 2018). We perform simulations with both isotropic and anisotropic diffusion (HSiso versus HSaniso) in high salinity oceans to test the robustness of our results. The isotropic diffusion is consistent with Soderlund (2019), and the anisotropic diffusion is consistent with , as well as with our low salinity simulation. In order to examine the influence of resolution, we set up a simulation with doubled horizontal resolution (HSiso05 compared with HSiso). We only test the effect of resolution in the high salinity ocean because the low salinity ocean takes a much longer time to reach an equilibrium state, due to the diffusive adjustment of the stratified layer, which makes the simulation computationally very expensive. All simulations are integrated to a near equilibrium state in which the energy imbalance is less than 2%. The presented results are 10-year-averages for high salinity cases and 250-year-averages for the low salinity case after this equilibrium has been reached. The longer averaging time in the low salinity case was chosen to account for a much larger low-frequency variability in this simulation.
To study tracer mixing processes in the ocean, we carry out three tracer simulations initialized from the equilibrium states of HSaniso, HSiso and LSaniso. We initialize two passive tracers (i.e. with no effect on ocean density and thus dynamics) at the bottom of the ocean at 0 • and 60 • S to study the evolution of tracer concentration. The turbulent diffusivity for the tracers is set to be the same as the thermal turbulent diffusivity. Each tracer simulation is run for 1500 years, and 10-year-averages are used for analysis.

Vertical Stratification and Circulation
In our high salinity simulations, the ocean stratification is weakly convectively unstable at all depths ( Figure 5 In HSiso and HSiso05, the resolved convection at mid-and high-latitudes is, however, weak (Figure 6(f) & (i)). This result can be understood by noting that the large vertical diffusivity and viscosity in these two cases leads to a relatively low Rayleigh number (Ra = αgδT D 3 /νκ where δT is the temperature contrast between the bottom and surface of the ocean and ν is viscosity) of Ra ≈ 1 × 10 6 in HSiso and Ra ≈ 1.5 × 10 7 in HSiso05. Importantly, the Rayleigh number in these two cases is smaller than the critical value that has been suggested for rotating convection: Ra S = 8.7E −4/3 where E = ν/(2ΩD 2 ) is the Ekman number (Cheng et al. 2018), which is E ≈ 1.5 × 10 −5 in HSiso and E ≈ 6 × 10 −6 in HSiso05. The critical Rayleigh number is thus Ra S ≈ 2 × 10 7 in HSiso and Ra S ≈ 8 × 10 7 in HSiso05, which is larger than the respective Rayleigh numbers. As a result, radial heat flux is dominated by the parameterized turbulent diffusion.
Temperature gradients are very small in all high salinity simulations, with the whole ocean near the freezing point. The bottom to surface buoyancy contrast is αgδT ≈ 10 −7 m s −2 in HSaniso and αgδT ≈ 10 −8 m s −2 in HSiso and HSiso05 (Figure 6(a), (d) & (g)). For comparison, the prescribed buoyancy contrast in Soderlund (2019) is around 10 −5 m s −2 for Enceladus, which is at least two orders of magnitude larger than our result, consistent with a radial heat flux in the simulations of Soderlund (2019) that is many orders of magnitude larger than the observed surface heat loss.   Figure 6. Ocean potential temperature and flow fields in the equilibrium state of different simulations. From left to right, the three columns show time-averaged zonal-mean fields of potential temperature θ, zonal velocity U and vertical kinetic energy W 2 , respectively. Note that the colorbar is saturated in panel (j). The surface temperature in LSaniso is very low (-0.6 • C) and increases with depth linearly in the stratified layer, as shown in Figure 5(e).
The zonal flow at mid-latitudes is constrained by the temperature structure via the thermal wind relationship, in particular the equator-to-pole temperature gradient, with warmer temperatures at the poles and colder temperatures at low latitudes, which leads to a negative (westward) vertical current shear at mid-latitudes (second colomn in Figure 6). In HSaniso, there exists superrotation in the equatorial upper ocean (Figure 6(b)), while such superrotation does not exist in HSiso and HSiso05 (Figure 6(e) & (h)), suggesting that the equatorial dynamics are sensitive to poorly constrained simulation parameters. The mechanism for the superrotation is associated with upward eddy momentum fluxes, consistent with previous studies (Aurnou & Olson 2001;Kaspi 2008;Ashkenazy & Tziperman 2020, see APPENDIX C for a more detailed discussion of the mechanism for superrotation). Note that the simulation results of HSiso and HSiso05 with different resolutions are qualitatively similar but differ quantitatively (compare the 2 nd and 3 rd rows in Figure 6).
As predicted, there are two different vertical layers in the low salinity ocean (Figure 5(d)-(f)). The bottom layer is a convective layer, which is slightly negatively stratified with positive thermal expansivity and small temperature variation around the critical point T c (Figure 5(e) & Figure 6(j)). In the time-averaged equilibrium state, convection in the convective layer in the low salinity ocean (Figure 6(l)) is weaker than in the high salinity ocean (Figure 6(c)), as expected due to the smaller thermal expansivity and thus smaller buoyancy input. The upper layer is the stratified layer with stable stratification and negative thermal expansivity, so that convection is blocked at the bottom of this layer (Figure 6(l)). The vertical temperature gradient is relatively large in this layer, with a profile that is linearly increasing from the freezing point T f at the surface of the ocean to the critical temperature T c at the interface between the two layers. We can apply Equation (1) to estimate the expected depth of the stratified layer. In our simulation, the temperature variation is |∆T | ≈2.0 K, the vertical heat flux is Q ≈0.03 W m −2 , and the vertical thermal diffusivity is κ z,heat = 5×10 −5 m 2 s −1 , such that the predicted depth of the stratified layer is around 14 km. This estimated depth matches well with the simulation results ( Figure 5(d)-(f)). Notice that the specific thickness of the stratified layer is sensitive to our assumption for the vertical diffusivity, which is poorly constrained. However, the fact that the numerical simulations support the theoretical predictions of Section 2 lends support to the more general scaling arguments discussed there.

Ocean Heat Transport
The mechanisms governing ocean heat transport vary widely across our simulations. With anisotropic diffusivity and viscosity (HSaniso) the vertical heat flux in the high salinity simulation is dominated by explicitly resolved convection (Figure 7(a)). With isotropic diffusivity and viscosity (HSiso and HSiso05 ) instead, the vertical heat flux is dominated by parameterized turbulent diffusion (Figure 7(b) & (c)) due to the large vertical diffusivity.
The choice of anisotropic versus isotropic diffusivity and viscosity also affects the heat flux at the ice-ocean interface, and by comparing the ocean surface heat flux from simulation results with the conceptual models in Section 2.3 we find that different aspects appear to play a role in different configurations. In HSaniso, the peaks at around 60 • match well with the "Perpendicular to Rotation" model, indicating the important role of heat transport by the horizontal currents around Taylor columns. Heat fluxes are further enhanced in the equatorial region, associated with particularly strong "equatorial rolls", which have previously been discussed in . Notice that the strong peaks near the poles are likely to be associated with the artificial boundaries at 85 • N/S (Figure 8(a)). Although the result that vertical heat transport is dominated by convection and associated Taylor columns is likely to be robust, the specific patterns of the convection are likely to depend on model resolution and poorly constrained parameters, and the resolution remains inadequate to resolve the natural scale of rotating convection (c.f. Section 3.1). In particular, the highly anisotropic viscosity and diffusion coefficients clearly affect the dynamics significantly. Although unresolved turbulence is likely to be anisotropic, even in a high salinity scenario, both the strength and alignment of the anisotropy remain unknown.
In HSiso and HSiso05, the surface heat fluxes are very close to the "Parallel to Gravity" model ( Figure 8(b) & (c)), consistent with the result that vertical heat flux is dominated by parameterized isotropic turbulent diffusion. The heat flux at the equator is somewhat larger than predicted by the "Parallel to Gravity" model because it is enhanced by heat transport through the "equatorial rolls". In HSiso05, there are small peaks outside the tangent cylinder, which are caused by convection (c.f. Figure 6(i)). These peaks qualitatively resemble those predicted by the "Parallel to (O(10 mW m −2 )) and heat loss rates increasing with latitude, although the spatial variations of the heat fluxes are smaller in our simulations.
In the low salinity case, LSaniso, the vertical heat transport is dominated by convection in the bottom convective layer, but is dominated by diffusion in the upper stratified layer (Figure 7(d)). The surface heat flux is almost uniform, and well approximated by the "Perpendicular to Gravity" model ( Figure 8(d)). The approximately uniform surface heat flux is a result of the relatively strong horizontal mixing in the stably stratified layer where the horizontal diffusivity is much larger than the vertical diffusivity. The peak at the equator is associated with grid-scale vertical advection, which is likely to be an artifact of numerical instability.

Ocean Tracer Mixing
The three tracer simulations illustrate the mixing processes in the ocean and the time it takes for constituents to be mixed from the bottom to the surface of the ocean. In the high salinity ocean, tracers can reach the surface of the ocean within tens of years. The convective mixing time scale is within an order of magnitude but somewhat shorter than the estimate according to Equation (4) because the scaling underestimates the vertical velocities in the simulations. Horizontal mixing across the hemisphere takes around 1000 years, while global mixing across the equator remains incomplete by the end of our 1500-year simulation (Figure 9). The mixing time scales are only moderately dependent on the choice of isotropic versus anisotropic diffusion. Note that the poleward drift of the highest concentration center and the tilting pattern arise due to spherical geometry effects on the diffusive transport as well as advection by meridional currents.
In the low salinity simulation, LSaniso, the mixing time scale in the bottom convective layer is similar to that in the high salinity simulations. However, the tracers do not penetrate significantly into the stratified layer, even after 1500 years (Figure 9). The vertical mixing time scale in the stratified layer can be estimated using Equation (6). In our simulation, κ z,heat = κ z,tracer = 5×10 −5 m 2 s −1 , so that τ dif f ≈ 4 × 10 12 s ≈ 1 × 10 5 years. This is consistent with the simulation results showing little tracer penetration into the stratified layer after 1500 years. Note that this specific time scale is sensitive to our assumption for the vertical diffusivity, which is poorly constrained. However, the fact that the numerical simulations are consistent with the scaling arguments in Section 2, increases our confidence in the theoretical prediction that the vertical mixing time scale in the stratified layer would be between hundreds and hundreds of thousands of years, depending on the assumed turbulent vertical diffusivity.

CONCLUSION
We find that the salinity of Enceladus' ocean fundamentally determines the vertical stratification and circulation of the ocean, thereby affecting the heat and tracer transport from the rocky core to the surface ice shell. If salinity is high, the ocean is unstratified and is dominated by convection over the whole depth. Tracers can be transported from the bottom to the surface within tens of years. In a low salinity ocean (below a critical point of about 20 g kg −1 ), a stably stratified layer exists in the upper ocean, whose depth may be anywhere from tens of meters to the full depth of the ocean, depending on assumptions about the strength of tidally (and/or librationally) driven turbulence that can contribute to mixing of the stratified layer. The stratified layer strongly suppresses vertical mixing, leading to a vertical mixing time scale between hundreds and hundreds of thousands of years. In the presence of significant horizontal mixing, the stratified layer further leads to horizontally homogeneous heat transport to the ice shell, independent of any spatial structure in the heat flux underneath. More detailed simulation results are sensitive to model resolution, as well as the magnitudes of the assumed turbulent viscosities and diffusivities, which remain poorly constrained.
Our simulations are based on an ocean-only model, hence the effect of freezing and melting on salinity as well as variations in the ice thickness are not included. When the ice thickness varies, the temperature at the ice-ocean interface (equal to the local freezing point T f ) decreases with depth. However, this temperature variation is around 0.1 K, assuming an ice shell thickness variation of around 20 km (c.f. Figure 2 in Kang et al. 2021), which is much smaller than the difference between the freezing point T f and the critical temperature T c in a low salinity ocean (T c − T f >1 K at a salinity of 8.5 g kg −1 ). Spatial variations in the temperature contrast across the stratified layer and associated variations in layer thickness are hence expected to be relatively small. We moreover expect that horizontal variations in the layer thickness and height would be kept small by baroclinic instability and circulations that would arise in the presence of significant horizontal density gradients. A crude sketch of the temperature structure that may be expected in the presence of varying ice thickness is shown in Figure 10, although we note that this prediction is speculative and needs to be tested. Future work is required to investigate dynamics of the ocean on Enceladus when the ocean is coupled with an ice shell of spatially varying thickness. Care must be taken in any such study to investigate the sensitivity of the results to model parameters and resolution, which are always based on compromises enforced by limited computational resources.
Our results indicate a contradiction in estimates of salinity and vertical mixing time scale of Enceladus' ocean in previous studies. Studies on the geochemistry of Enceladus' ocean have generally indicated a relatively low salinity, likely less than 20 g kg −1 (e.g. Postberg et al. 2009;Hsu et al. 2015;Glein et al. 2018). In this case we expect a stably stratified layer to form and result in a vertical mixing time scale between the rocky sea floor and the ice shell of at least hundreds of years. However, the detection of silica nanoparticles in plumes has been argued to set an upper limit of several years on the mixing time scale, based on the size and growth rate of the particles (Hsu et al. 2015). Our Figure 10. Sketch of the expected temperature structure of the low salinity ocean with varying ice shell thickness. Grey dashed lines indicate the contour lines of the temperature, which increases with depth in the stratified layer. The temperature is fixed near the local freezing point T f at the ice-ocean interface and critical temperature Tc at the interface between the stratified layer and the convective layer. The temperature at the surface of the ocean is colder at the equator but warmer at the poles because the freezing point T f decreases with pressure, but this temperature contrast is small compared to the temperature contrast across the stratified layer.

Stratified layer
results suggest that these inferences are not compatible with each other, indicating that Enceladus' ocean is either saltier than previously suggested or the interpretation of silica nanoparticles needs to be reconsidered. One possible explanation is that melting of ice in the polar regions due to stronger tidal heating ) may freshen the water expelled in the plumes, thus possibly leading to an underestimate of the bulk ocean salinity by observations of plume constituents. However, more work is needed to establish a consistent picture of Enceladus' ocean salinity and mixing time scale.

ACKNOWLEDGMENTS
We are grateful to Jun Yang, Edwin S. Kite, Wanying Kang, Dorian S. Abbot, Mikael Beuthe and two anonymous reviewers for helpful discussions and comments. Y.Z. thanks the Department of Atmospheric and Oceanic Sciences and School of Physics at Peking University for the financial support during the summer exchange program at University of Chicago. This work was completed with resources provided by the University of Chicago Research Computing Center.

A. CALCULATION OF SURFACE HEAT FLUXES IN FOUR CONCEPTUAL MODELS
In this part, we will explain the calculation of surface heat fluxes in the four conceptual models shown in Figure 3. In the "Parallel to Rotation" model, we relate points at the ice-ocean interface to those on the sea floor that lie on the same line parallel to the axis of rotation (Figure 11(a)): where the subscript t indicates the top of the ocean and the subscript b indicates the bottom of the ocean, r is radius and θ is latitude. Setting θ b = 0 we can solve for the latitude of the tangent cylinder at the top of the ocean: If the bottom heating is entirely transported to the surface by the columnized convection, we moreover find where Q is the heat flux. From Equation (A1) we get dθ b /dθ t = (r t sin θ t )/(r b sin θ b ), so that Equation (A3) can be solved for the heat flux at the ocean surface, Q t (θ t ), inside the tangent cylinder. The surface heat flux is zero outside the tangent cylinder, so that it can generally be written as where θ b is determined by Equation (A1) and θ c1 t is given by Equation (A2). Notice that Equation (A4) has a singularity at the latitude where the surface intersects with the tangent cylinder.
In the "Perpendicular to Rotation" model ( Figure 11(b)), following a similar argument as in the "Parallel to Rotation" model, we have which gives a critical latitude By substituting Equation (A5) into Equation (A3) (which is also valid in this model), we have (a) (b) Figure 11. Schematic for the "Parallel to Rotation" model (a) and the "Perpendicular to Rotation" model (b). The rotation axis is the z-axis, and other symbols are explained in the text.
where θ b is determined by Equation (A5) and θ c2 t is given by Equation (A6). In the "Parallel to Gravity" model, heat fluxes at the sea floor can be matched directly to surface fluxes at the same latitude, reduced only by a factor that accounts for the radial increase in surface area: In the "Perpendicular to Gravity" model, the surface heat flux is independent of latitude and simply given by the total heat flux divided by the surface area: (A9) Note that the detailed patterns of the surface heat flux vary with the bottom heat flux patterns, except for the uniform surface flux in the model where mixing is "Perpendicular to Gravity", which is independent of the bottom heat flux pattern. However, Equations (A4), (A7), (A8) and (A9) can be used to compute the surface heat fluxes in the respective limit cases for arbitrary bottom heat flux distributions, and many of the main features of the surface heat flux distribution are robust for moderate changes in the bottom heat flux pattern.

B. ADDITIONAL SIMULATION DETAILS
Due to the deep ocean (40 km) compared to the planetary radius (252 km), our GCM simulations do not use the thin-shell approximation typically applied in Earth-like simulations (Adcroft et al. 2018). As a result, the variation of grid cell area with depths at the same latitude and longitude is taken into consideration, which allows vertical fluxes to be calculated more precisely. The vertical variation of gravity is also taken into account. We calculate the gravity as a function of depth as g(r) = g c r 2 c r 2 + 4πGρ w r 2 r rc r 2 dr , where G is the gravitational constant, r c is the radius at the surface of the solid core, i.e. the bottom of the ocean, ρ w =1000 kg m −3 is the density of water (and also ice, for simplicity) and g c = 0.127 m s −2 is the gravity at the bottom of the ocean, chosen such as to keep the gravity at the surface of the ice shell g(r s ) to be the observed value of 0.113 m s −2 . The vertical profile of the gravity is shown in Figure 4(a). Various studies have estimated the amplitude and patterns of tidal dissipation in the solid core on Enceladus (e.g. Beuthe 2013; Choblet et al. 2017;Beuthe 2018Beuthe , 2019. Since small variations in the bottom heat flux pattern are not expected to affect the main results of this study, we assume a simplified bottom heating pattern that reflects the tidal disturbing potential. The global total energy input is assumed to be 20 GW, with the maximum heat flux twice of the minimum (Choblet et al. 2017). The magnitude of the bottom heat flux may affect the convective mixing time scale (τ conv ∝ Q −1/2 ; Equation (4)), the depth of the stratified layer (H ∝ Q −1 ; (Equation (1)), and the diffusive time scale in the stratified layer (τ dif f ∝ Q −2 ; (Equation (6)). Generally, a larger bottom heat flux will result in a shorter vertical convective mixing time scale, a shallower stratified layer and hence a shorter vertical diffusive mixing time scale in that layer. However, unless the bottom heat flux differs by orders of magnitude from the value assumed here, our main conclusions are expected to be robust. The heating pattern is then suggested to be: where F total = 20 GW is the total tidal heating, C 0 = (0.5 + √ 3/8) 5/π ≈ 0.904 is a constant, Θ is co-latitude, λ is longitude, and Y 20 and Y 22 are degree-2 spherical harmonic functions. Since we are only simulating part of the zonal range of the global ocean (15 • in longitude), we cannot adequately simulate zonal variations. We therefore apply a zonally symmetric heating profile given by the zonal mean of Equation (B2), which is (Figure 4(b)): In simulations with varying horizontal resolution, the viscosity coefficient is modified such that the Kolmogorov scale is proportional to the grid scale (Vallis 2017). The Kolmogorov scale L ν is the length scale at which the viscosity becomes important, and scales as: where ε is the turbulent energy cascade rate, which should be independent of resolution, and ν is the turbulent viscosity. Setting L ν ∼ L grid , we can get the relationship:

HSiso:
High ! ! " "# (d) Figure 12. Momentum budget and equatorial eddies in HSaniso (left) and HSiso (right). Panels (a)-(c) show ten-year-averaged zonal-mean zonal acceleration due to Coriolis force, dissipation, and advection terms, respectively. Black lines in (a)-(c) are zero contour lines. Panel (d) shows snapshots of vertical velocity in the equatorial plane. Note that the equatorial convective patterns are grid-scale because the horizontal resolution is not able to resolve the intrinsic scale of the convective plumes, but within the constraint imposed by this limitation, the equatorial velocity structure in HSaniso shows a vertical tilt, which is associated with upward momentum transport.
We here choose L grid as the horizontal grid scale, which in our model is generally larger than the vertical grid scale. The horizontal grid length in case HSiso05 is half of that in case HSiso, so that the viscosity is correspondingly changed from 2.5 m 2 s −1 to 2.5 × 0.5 4/3 ≈ 1.0 m 2 s −1 . We set other parameters in case HSiso05 by keeping the diffusivity and viscosity isotropic and P r=10 (Table 2).

C. MECHANISMS FOR SUPERROTATION
There exists superrotation in the upper ocean in HSaniso but not in HSiso and HSiso05. Here we compare HSaniso and HSiso to analyze the mechanism for equatorial superrotation. The two isotropic simulations, HSiso and HSiso05, are similar (not shown).
The zonal-mean zonal momentum budget is (Vallis 2017): where the overbar indicates the zonal and temporal average; u, v, w are zonal, meridional and vertical velocity, respectively; θ is latitude, and F x is dissipation. On the right-hand-side, the terms in the first line are the Coriolis force terms, the term in the second line is the dissipation term, and the terms in the third line are the advection terms (including metric terms). In an equilibrium state, the left-hand-side should be zero. Superrotation in HSaniso, which is largest around r = 219 km (around 13 km in depth), is driven by momentum flux convergence (Figure 12(c)) and counteracted by Coriolis terms and frictional dissipation (Figure 12(a) & (b)). In HSiso, however, such a momentum flux convergence is absent. A cross section of vertical velocity along the equator (Figure 12(d)) shows that the upward momentum flux in HSaniso is driven by vertically tilted eddies (Aurnou & Olson 2001;Kaspi 2008). Although similar eddies appear to exist in HSiso, they are weaker and not tilted, and hence do not carry momentum upwards.