The effect of salinity on ocean circulation and ice-ocean interaction on Enceladus

Observational data suggest that the ice shell on Enceladus is thicker at the equator than at the pole, indicating an equator-to-pole ice flow. If the ice shell is in an equilibrium state, the mass transport of the ice flow must be balanced by the freezing and melting of the ice shell, which in turn is modulated by the ocean heat transport. Here we use a numerical ocean model to study the ice-ocean interaction and ocean circulation on Enceladus with different salinities. We find that salinity fundamentally determines the ocean stratification. A stratified layer forms in the low salinity ocean, affecting the ocean circulation and heat transport. However, in the absence of tidal heating in the ice shell, the ocean heat transport is found to always be towards lower latitudes, resulting in freezing at the poles, which cannot maintain the ice shell geometry against the equator-to-pole ice flow. The simulation results suggest that either the ice shell on Enceladus is not in an equilibrium state, or tidal dissipation in the ice shell is important in maintaining the ice shell geometry. The simulations also suggest that a positive feedback between cross-equatorial ocean heat transport and ice melting results in spontaneous symmetry breaking between the two hemispheres. This feedback may play a role in the observed interhemispheric asymmetry in the ice shell.


INTRODUCTION
Enceladus is believed to have a global ocean of about 40 km in depth under a global ice shell (e.g., Postberg et al. 2011;Patthoff & Kattenhorn 2011;Thomas et al. 2016).The energy source that maintains the global liquid ocean and a heat loss of about 10 GW at the south polar region (Spencer et al. 2006;Howett et al. 2011;Spencer et al. 2013) is believed to primarily come from tidal dissipation.In general, tidal dissipation is expected to occur in the ice shell, in the ocean, and in the solid core.Tidal heating in the ocean has been thought to be negligible compared to that in the ice shell and inner solid core (e.g., Chen & Nimmo 2011;Tyler 2011;Chen et al. 2014;Beuthe 2016;Hay & Matsuyama 2017;Matsuyama et al. 2018;Hay & Matsuyama 2019).Although recent work suggests dissipation of internal tides may generate more heat in a stratified ocean, if resonances occur (Rovira-Navarro et al. 2023), the stratification in Enceladus's ocean is unlikely to be strong enough to support a heat generation rate that is comparable to the heat loss rate (e.g., Kang et al. 2022c).Tidal dissipation in the ice shell, as well as the turbulent dissipation of water in conduits at the south pole (Kite & Rubin 2016), has been argued to potentially reach O(10 GW) (Beuthe 2019;Souček et al. 2019).However, the rheology of ice on Enceladus is poorly understood, which strongly affects the tidal dissipation rate in the ice shell.The low density of Enceladus suggests that the core is porous, which possibly supports strong tidal dissipation in the core, reaching O(10 GW) (Iess et al. 2014;Roberts 2015;Čadek et al. 2016;Choblet et al. 2017).Libration has been argued to cause dissipative heating in the ocean of up to O(0.1 GW) (Lemasquerier et al. 2017;Wilson & Kerswell 2018;Rekier et al. 2019;Soderlund et al. 2020).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.
The ice shell on Enceladus is about 20 km thick on average, and is not flat: it is thickest at the equator (over 30 km) and thinnest at the south polar region (less than 10 km).There is also an interhemispheric asymmetry in the ice shell with ice being thinner at the south pole than at the north pole (Beuthe et al. 2016;Čadek et al. 2019;Hemingway & Mittal 2019).The significant gradient in the ice shell thickness suggests an equator-to-pole ice flow.If the ice shell on Enceladus is in an equilibrium state, the equator-to-pole mass transport of the ice flow must be balanced by ice freezing at the equator and melting at the pole.Ice shell models generally predict that a freezing/melting rate of a few millimeters per year2 on Enceladus is required to balance the equator-to-pole ice flow, or equivalently, O(0.01 W m −2 ) of latent heat ( Čadek et al. 2019;Kang & Flierl 2020;Kang et al. 2022c).
Freezing and melting of the ice is determined by the heat budget.The associated latent heat (Q LH , with melting defined as positive) has to match the imbalance between the heat flux from the ocean to the ice (Q oi ) and the upward heat flux within the ice shell.Since the ice shell on Enceladus is unlikely to be convecting (Nimmo et al. 2018), we assume that the upward heat flux in the ice is purely due to conduction (Q cond ).The energy balance at the ice-ocean interface is then: The conductive heat flux at the ice-ocean interface in turn is modulated by the ice shell thickness, as well as the tidal heating in the ice shell (Q i ).Without strongly localized tidal heating in the ice shell, conductive heat loss is strongest at the pole where the ice is thinnest, which would favor ice freezing at the pole -the opposite of what is needed to maintain the ice geometry.Tidal dissipation in the ice shell could be a possible answer to the problem, since it is strongest at the pole and is amplified where the ice is thin (e.g., Beuthe 2018Beuthe , 2019;;Kang & Flierl 2020), and thus could potentially allow for melting at the pole.Alternatively, the heat flux from the ocean, Q oi , could be amplified at the pole.The global ocean heat flux to the ice comes from the tidal heating in the bottom solid core that peaks at the pole (Choblet et al. 2017), but the spatial pattern is modulated by the ocean circulation, which is the focus of this paper.
Previous studies have explored the global ocean circulation on Enceladus and other icy moons (Soderlund et al. 2014;Travis & Schubert 2015;Soderlund 2019;Amit et al. 2020;Lobo et al. 2020;Kang et al. 2020;Zeng & Jansen 2021;Ashkenazy & Tziperman 2021;Bire et al. 2022;Kvorka & Čadek 2022;Kang 2022;Kang & Jansen 2022;Kang et al. 2022a,c;Jansen et al. 2023).These studies have used a range of different methods (e.g., conceptual models, ocean-only numerical simulations, and ice-ocean coupled numerical simulations), in either 2-D or 3-D domains, with a wide range of parameters and different boundary conditions.However, there have not yet been simulations with fully coupled ice-ocean thermodynamics in a 3-D global domain.The published simulation results vary substantially in their ocean circulation and ocean-ice heat fluxes, likely as a result of the large uncertainties in the external parameters and the difficulties in modeling the global ocean in a realistic parameter regime.The real regime of ocean circulation and heat transport on Enceladus therefore remains uncertain.
One important factor determining the ocean circulation and heat transport is salinity (S).Salinity determines the equation of state of water, thus affecting ocean stratification and circulation.Most importantly, if the salinity is lower than a critical value, the thermal expansivity (α ≡ −(1/ρ)(∂ρ/∂T ) where ρ is density and T is temperature) is negative near the freezing point.As a result, in such a low salinity ocean, a stably stratified layer forms in the upper ocean with bottom heating (Melosh et al. 2004;Zeng & Jansen 2021).The threshold for the negative thermal expansivity to exist under the pressure found in Enceladus's ocean is around 20 g kg −1 (Zeng & Jansen 2021;Kang et al. 2022c).There are different methods to estimate the salinity of Enceladus's ocean, based on the observations of salt-rich ice grains (Postberg et al. 2009), the existence of silica nanoparticles (Hsu et al. 2015), using models of boiling water in ice conduits (Ingersoll & Nakajima 2016), equilibrated water-rock interaction (Zolotov 2007), ocean geochemistry (Glein et al. 2018), and ocean heat transport (Kang et al. 2022c).The estimates vary from study to study, giving a wide range of 2-40 g kg −1 .Therefore, both "high salinity" (S > 20 g kg −1 ) and "low salinity" (S < 20 g kg −1 ) scenarios should be carefully considered when studying the ocean circulation on Enceladus.
Salinity can also influence the overturning circulation driven by surface density contrasts, and modulate the ocean heat transport, which further affects the ice geometry.Kang et al. (2022c) model the ocean circulation with prescribed ice thickness and prescribe melting rates at the pole and freezing rates at the equator as required to balance the estimated ice flow induced by the ice thickness gradient.They find that the surface density contrast between the equator and the pole drives an overturning circulation in each hemisphere.The direction of the circulation depends on the mean salinity of the ocean, but the ocean heat transport is found to always be equatorward because the freezing point at the pole (with thinner ice and hence lower pressure) is higher than at the equator (with thicker ice and hence higher pressure), and the heat transport is always from the warmer pole to the colder equator.The equatorward ocean heat transport favors freezing at the pole and melting at the equator and thereby acts to reduce the equator-to-pole ice shell thickness gradient, which is known as the ice pump effect (Lewis & Perkin 1986).
The partitioning of tidal heating between the ice shell and the bottom solid core may also affect the ocean heat transport, further influencing the asymmetric ice shell geometry between the two hemispheres.Kang et al. (2022a) prescribe an ice geometry with an interhemispheric asymmetry where ice is slightly thinner in the southern hemisphere.They incorporate fully coupled ice-ocean thermodynamics in a 2-D domain with tidal dissipation in both the ice shell and the bottom solid core.They find that when the heating is dominated by tidal dissipation in the ice shell, there is stronger melting at the south pole than the north, and the overturning circulation transports heat to  the south, thus enhancing the asymmetry; instead, when the heating is primarily in the core, there is stronger freezing at the south pole than the north, and the overturning transports heat to the north, thus suppressing the asymmetry.
Here we study the effect of salinity on the ocean circulation and ice-ocean interaction on Enceladus, focusing on the ocean heat transport.We carry out 3-D global ocean simulations using the Massachusetts Institute of Technology General Circulation Model (MITgcm) (Adcroft et al. 2018), with prescribed ice geometry and fully coupled ice-ocean thermodynamics.We exclude tidal heating in the ice to focus specifically on whether ocean heat transport from the solid core can maintain the observed ice topography and the associated pattern of freezing and melting.Section 2 describes the model configuration and experimental design.Section 3 introduces the numerical simulation results, focusing on ocean stratification, circulation, and heat transport.Section 4 compares the simulation results with previous studies and provides concluding remarks.

EXPERIMENTAL DESIGN
We perform numerical simulations using the MITgcm to solve the non-hydrostatic equations for a Boussinesq fluid in a rotating spherical shell in a 3-D domain, where all sphericity terms are preserved, including all components of the Coriolis force (Adcroft et al. 2018).Radius, rotation rate and gravity are set to be the same as Enceladus, and the vertical variation of gravity is also taken into consideration (see Table 1 and Zeng & Jansen 2021).We use a non-linear equation of state (Jackett & Mcdougall 1995) so that the thermal expansivity becomes negative near the freezing point in a low salinity ocean.
We simulate an ocean with a depth of 40 km on average 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.At the surface of the ocean, we prescribe the ice geometry and simulate the ice-ocean thermodynamics.In order to include the ice pump effect, we set the ice geometry to be an idealized (zonally-invariant and interhemispherically symmetric) approximation to Enceladus's inferred ice topography (Tajeddine et al. 2017;Čadek et al. 2019), with ice thinning towards the pole.We do not allow the ice shell thickness to evolve (which would occur on a timescale that is much longer than the duration of the simulations) but compare the freezing and melting rate calculated in the simulations with the expected rate that is required to maintain the ice geometry against the ice flow.This approach is the same as used in Kang et al. (2022a) but different from Kang et al. (2022c) (where the freezing and melting rates are prescribed).At the bottom of the ocean we apply a prescribed heat flux based on Choblet et al. (2017) (Fig. 1).The detailed model description can be found in APPENDIX A.
We carry out simulations with both high salinity (35 g kg −1 , same as Earth ocean) and low salinity (8.5 g kg −1 , suggested by Glein et al. 2018) to study the effect of salinity on ocean circulation and ice-ocean interaction on Enceladus.For each salinity, we use two different horizontal turbulent diffusivities (κ h ), 0.05 m 2 s −1 and 0.5 m 2 s −1 , to test how sensitive the results are to this poorly constrained parameter, which is meant to mimic the effect of unresolved eddies and directly affects the meridional heat transport.The experimental setups are summarized in Table 2.All simulations are integrated to a quasi-equilibrium state where the total energy imbalance is less than 3% of the seafloor heating rate (see APPENDIX A.3).

SIMULATION RESULTS
Ocean stratification is highly dependent on the mean ocean salinity.The high salinity ocean is mostly unstratified or weakly stratified below the ice-ocean interface.Instead, as previously found in ocean-only simulations with a flat upper boundary and constant salinity, a stably stratified layer develops in the low salinity ocean (Fig. 2).In the stratified layer, the temperature decreases with height due to bottom heating, which, due to the negative thermal expansivity, amounts to a statically stable stratification (Zeng & Jansen 2021).The stratified layer maintains a large vertical thermal contrast due to the difference between the freezing point and the critical temperature at which the thermal expansivity changes sign (around 2 • C difference with a salinity of 8.5 g kg −1 , Fig. following the ice shell topography, indicating that ocean circulation is efficient in removing any large horizontal density gradients that would be associated with an undulating stratified layer.The temperature and salinity structure is strongly affected by the surface ice shell in both high and low salinity simulations.In the upper ocean, the ocean is in contact with the ice shell, and the temperature is near the freezing point at the ice-ocean interface, which decreases with increasing pressure and hence with increasing ice thickness.Therefore, the temperature increases poleward in the upper ocean (Fig. 2a).The salinity mostly increases poleward except for the equatorial region in the low salinity ocean (Fig. 2b).The density mostly increases poleward in all simulations, although the role that temperature and salinity play in setting this structure differs in the high versus low salinity oceans (Fig. 2c).In the high salinity ocean, the salinity and temperature gradients act in the opposite way, but the contribution from the salinity anomaly is about 2-3 times larger than that of the temperature anomaly, so that the density increases poleward.In the low salinity ocean, the thermal expansivity is negative, so that both temperature and salinity anomalies contribute to the poleward-increasing density, with the contribution from the temperature anomaly 3-5 times larger than the salinity anomaly.
The density structure controls the zonal flow via the thermal wind relationship, which can be written as (see e.g.Kaspi et al. 2009;Bire et al. 2022): where ⃗ Ω is the planetary rotation rate, u is the zonal velocity, Z is the axis parallel to rotation, r is radius, ϕ is latitude, and b ≡ −gδρ/ρ 0 is buoyancy, where g is gravity, δρ is the density anomaly and ρ 0 is a constant reference density.The thermal wind balance holds very well in all simulations (Fig. 3), except in regions where there is strong grid-scale convection (c.f.Fig. 4b & c).With density increasing poleward (i.e.−∂b/∂ϕ > 0 in the northern hemisphere; and opposite in the southern  hemisphere), at mid-and high-latitudes we therefore find ∂u/∂Z > 0, and there are eastward jets in the upper ocean with the jet cores aligned with the axis of rotation.Near the equator, sin ϕ is close to zero and therefore cos ϕ∂u/(r∂ϕ) > 0, so that the zonal velocity decreases equatorward and there is a westward current near the equator (Fig. 4a1-a3).An exception is that in the low salinity simulation with small horizontal diffusivity (LS l ), an eastward flow exists near the equatorial surface (Fig. 4a4) as a result of the locally reversed density gradient (Fig. 2b4 & c4).
In both high and low salinity oceans, there is a slantwise Hadley-like circulation that ascends from the seafloor at the equator, following the tangent cylinder (a surface of constant angular momentum), and descends at higher latitudes.In the low salinity ocean, there is another reversed circulation cell at mid-latitudes along with an eastward jet, which bears similarity to the Ferrel cell in Earth's atmosphere.We also see grid-scale slantwise convection near the surface at high latitudes (Fig. 4b).Although the radial stratification (parallel to the gravity vector) is near neutral or stable as maintained by the convective adjustment scheme in the model, the stratification along surfaces of constant angular momentum (mostly parallel to the rotational axis) can be unstable, which can cause slantwise convection.
A strong small-scale overturning cell exists near the equatorial surface in the high salinity simulation with high horizontal diffusivity (HS h , Fig. 4b1), which is associated with a positive feedback loop.There is a strong but narrow (only several degrees in latitude) melting region at and just south of the equator (Fig. 5a1).As a result, a local density minimum and sharp density gradient exists at the equatorial surface (Fig. 2c1).Following the thermal wind relationship, this density gradient is associated with a strong westward jet near the equatorial surface (Fig. 4a1), and the corresponding frictional drag is associated with an Ekman transport that drives a localized overturning circulation.The local circulation cell then provides heat to maintain strong melting, forming a positive feedback loop.Note that this overturning is associated with a slight asymmetry around the equator where the peak of the melting and hence the salinity minimum is slightly displaced towards the southern hemisphere.
An even stronger asymmetry exists in the high salinity simulation with low horizontal diffusivity (HS l ).There are overturning cells across the equator (Fig. 4b2), and the salinity profile is asymmetric between the two hemispheres as a result of asymmetric freezing and melting patterns (Fig. 2b2).The advection of heat and associated freezing and melting generates a positive feedback loop that can result in spontaneous symmetry breaking between the two hemispheres: once we get more melting in one hemisphere, the density becomes smaller, driving an overturning circulation across the equator, with sinking in the hemisphere with higher density and rising in the hemisphere with lower density.Due to the bottom heating, the temperature increases with depth at low latitudes, so that the overturning transports warmer water to the hemisphere with lower density, enhancing melting in that hemisphere.Note that in our simulation, this positive feedback seems to be confined to outside the tangent cylinder at low latitudes, although the details of the results are likely affected by poorly constrained model parameters and insufficient resolution.There is almost no asymmetry in the low salinity simulations, likely because the density gradient is dominated by the temperature gradient in the low salinity ocean.The temperature gradient is largely determined by the surface freezing point, which in turn is determined by the ice shell topography.As long as the ice shell topography remains symmetric (as assumed in this study) the temperature field also stays approximately symmetric (Fig. 2a).
The ocean circulation drives ocean heat transport and modulates the surface density structure through the freezing and melting rate, and the density structure in turn modulates the ocean circulation.The heat budget at the ice-ocean interface is Q cond + Q LH = Q oi , where the conductive heat loss (Q cond ) in the simulation is prescribed along with the ice geometry (see APPENDIX A.2).If the ocean is in an equilibrium state, the ocean heat flux to the ice shell (Q oi ) is balanced by the bottom heating (Q b ) when globally integrated.However, locally, the ocean heat transport determines the latent heat (i.e., the freezing and melting rate) by modulating the bottom heating pattern and providing an ocean heat flux to the ice shell.
In all simulations, there is strong freezing at the pole and melting at lower latitudes (low latitudes in high salinity simulations and mid-latitudes in low salinity simulations; see Fig. 5a).Locally, the dominant balance is between the latent heat flux (Q LH ) and the sensible heat flux from the ocean to the ice (Q oi ), because the spatial variations in the meridional heat flux convergence are much larger than the variations in the bottom heating and the conductive heat loss (Fig. 5b).This latent heating profile, where freezing happens at higher latitudes than melting, is a result of the ocean heat transport mostly towards lower latitudes, induced by the ice pump effect.The meridional -90 -60 -30 0 30 60 90 -0.5 temperature gradient is strongest in the upper ocean where the ocean is in contact with the ice shell, so that the temperature gradient is largely determined by the freezing point at the ice-ocean interface (Fig. 6).With polar-thinning ice geometry, the freezing point increases poleward so that the temperature is generally warmer at the pole than at the equator.As a result, any down-gradient heat flux is transporting warm water towards lower latitudes, so that melting occurs at lower latitudes than freezing.The strong freezing at the pole is moreover associated with a positive feedback: freezing at the pole leads to a density increase due to salt rejection, which destabilizes the stratification, thus promoting convection.Because the temperature decreases with depth under the ice shell at the poles (due to horizontal mixing with colder low-latitude water), the convection brings up colder water from below, which supports further freezing (see APPENDIX B.2).
The meridional heat flux is larger in the low salinity ocean than in the high salinity ocean (Fig. 5b).This result can be understood by noting that the meridional heat flux has three components: diffusive heat flux, resolved eddy heat flux, and heat flux by the mean circulation (See APPENDIX B.1).All three fluxes become larger when there is a stratified layer.The meridional heat transport by diffusion and resolved eddies both increase with the vertical integral of the meridional temperature gradient.The meridional temperature gradient near the surface is similar in both high and low salinity oceans, because the temperature here is controlled by the freezing point, and hence is determined by the ice geometry, which is the same in both simulations.The meridional temperature gradient in the unstratified convective layer is small in all simulations.However, in the low salinity ocean, the temperature gradient in the stratified layer is comparable to that near the surface, so that the vertical integral of the meridional temperature gradient (i.e. the vertical extent of the meridional temperature gradient) is larger in the low salinity ocean (Fig. 6), resulting in a higher meridional heat flux for a given eddy diffusivity.Additionally, the heat transport by the mean circulation is proportional to the strength of the overturning multiplied by the temperature contrast between the poleward and equatorward branches of the overturning circulation.Since the vertical temperature  contrast is much larger in the low salinity ocean with a stratified layer, the heat transport by the mean circulation is also expected to be larger in the low salinity ocean.
Comparing simulations with different horizontal diffusivities indicates that, although freezing at the pole is found robustly across all simulations, the detailed pattern of melting is sensitive to model parameters.In the high salinity simulations, there is persistent melting at the equator in the simulation with higher diffusivity (HS h ), while in the simulation with lower diffusivity (HS l ), the melting latitude varies temporally.In the northern hemisphere, there is an intermittent melting region outside the tangent cylinder (at low latitudes), associated with the spontaneous symmetry breaking mechanism (Fig. 7).The system switches between an asymmetric state and a more symmetric state with a period of about 1000-2000 years, but the cross-equatorial overturning never reverses in this simulation (although the sign of the circulation depends on the initial conditions -see APPENDIX C).In the southern hemisphere, there is intermittent melting in a region inside the tangent cylinder (at mid-latitudes).The southern hemisphere melting only persists for several hundreds of years and is likely associated with an instability near the tangent cylinder (Fig. 7).The strong freezing at the pole and the pattern of melting in the low salinity simulations are less sensitive to horizontal diffusivities.In both high-and low-salinity ocean simulations, when the horizontal diffusivity κ h is high, the meridional heat transport is dominated by diffusion; when κ h is low, the advective heat transport becomes more important, while the overall heat flux changes relatively little (Fig. 5b).The results indicate that increased eddy heat transport, associated with increased eddy activity (Fig. 8), tends to compensate for the decrease in the prescribed parameterized diffusivity.

DISCUSSION AND CONCLUSION
We find that ice topography and salinity fundamentally affect the ocean stratification, circulation, and heat transport on Enceladus.In both high and low salinity scenarios, the upper ocean temperature, salinity and density all increase poleward.The temperature is near the freezing point at the ice-ocean interface, which varies with pressure and is thus set by the ice shell thickness variation.The salinity pattern is strongly affected by the salinity flux through freezing and melting.The density anomalies are dominated by the salinity anomalies in the high salinity simulations and are dominated by the temperature anomalies in the low salinity simulations, which is qualitatively consistent with Kang et al. (2022c).
In a low salinity ocean, with tidal heating at the bottom, a stably stratified layer forms in the upper ocean, as previously found by Zeng & Jansen (2021) in simulations with a flat upper boundary and no thermodynamic ice-ocean interface.This global thermally-driven stratified layer is caused by bottom heating and the negative thermal expansivity near the freezing point in a low salinity ocean, and is thus different from locally stable salt-stratified layers that form by surface melting, like the "freshwater lens" at the polar surface ocean discussed in Lobo et al. (2020), and the stratification in the "shell-heating" scenario in Kang et al. (2022a,c).Note that the depth of the stratified layer here is deeper than in the more idealized simulations of Zeng & Jansen (2021), but is thinner than in the "core-heating", low ocean in Kang et al. (2022a,c), where the entire ocean is stratified.The depth of the stratified layer H is proportional to the vertical temperature contrast |∆T r |, vertical diffusivity κ r , and the inverse of the vertical diffusive heat flux (Zeng & Jansen 2021).The vertical turbulent diffusivity in Kang et al. (2022a,c) is set 100 times larger than in our simulations and in Zeng & Jansen (2021), such that the depth of the stratified layer exceeds the depth of the ocean in Kang et al. (2022a,c), leaving the entire ocean stratified.The smaller difference between our simulations and those in Zeng & Jansen (2021) can be understood by noting that all heat flux in the stratified layer is diffusive in the simulations in Zeng & Jansen (2021), while in the simulations presented in this manuscript, some vertical heat transport is accomplished by advection.Therefore, Q dif f here is somewhat smaller and H is somewhat larger than in Zeng & Jansen (2021).This small discrepancy is likely due to the much larger meridional temperature gradient in this work, which in turn arises from the freezing point variations due to the ice shell topography.
In our simulations, the convergence of the meridional heat flux is much stronger than the variation of the bottom heating and surface conductive heat loss.As a result, the ocean heat transport becomes the primary control on the pattern of ice freezing and melting, which is consistent with the coreheating scenario in Kang et al. (2022a) but is different from Zeng & Jansen (2021), which uses a flat ice shell.With ice thickness decreasing from equator to pole, the freezing point increases poleward, resulting in a much stronger meridional temperature gradient and hence a stronger meridional heat flux in this work and Kang et al. (2022a) compared to Zeng & Jansen (2021) where the ice shell is flat.
The meridional heat flux is mostly towards lower latitudes in both high-and low-salinity oceans, due to the ice pump effect.This heat transport can be understood qualitatively by noting that the freezing point increases poleward due to the thinning of the ice shell, and the ocean heat transport tends to be directed down the meridional temperature gradient.Consequently, we find freezing at the poles, which cannot maintain the ice geometry against the poleward ice flow, as also found in Kang et al. (2022a,c).The result that freezing happens at the poles indicates that either the ice shell on Enceladus is not in an equilibrium state, or that tidal heating in the ice shell is important in maintaining the ice geometry.
The stratified layer in the low salinity ocean tends to enhance the meridional heat transport by the mean overturning circulation through increasing the vertical temperature contrast, and enhance the meridional heat transport by diffusion and resolved eddies through increasing the vertical extent of the meridional temperature gradient.The first mechanism is also found in Kang et al. (2022c).In their core-heating scenario, the ocean heat transport is stronger in the low salinity ocean compared to the high salinity ocean, despite a weaker overturning.
In our high salinity simulation with small diffusivity, we find a positive feedback where enhanced melting in one hemisphere drives a cross-equatorial overturning circulation that further enhances the melting.This feedback can result in spontaneous symmetry breaking between the two hemispheres.Although the specific simulation results here cannot explain the observed ice thickness distribution, it is possible that a similar mechanism could play a role in the interhemispheric asymmetry of Enceladus's ice shell thickness.Indeed, a similar mechanism was previously identified by Kang et al. (2022a).In the core-heating scenario in Kang et al. (2022a), they however find that this ocean circulation acts to weaken the asymmetry of the ice shell, because melting occurs where the ice is thick in their simulations, which use an asymmetric ice shell topography.The combined results therefore indicate that the positive feedback may only hold as long as the asymmetry remains relatively small.With asymmetric tidal heating in the ice shell, Kang et al. (2022a) instead find enhanced melting in the hemisphere where the ice shell is thin, such that the ocean heat transport feedback can further enhance the asymmetry of the ice shell.
Besides Enceladus, previous studies have also investigated the heat transport in oceans on other icy moons (e.g., Soderlund 2019;Amit et al. 2020;Ashkenazy & Tziperman 2021;Bire et al. 2022;Kvorka & Čadek 2022).In studies with direct numerical simulation (DNS) methods, simulations have been carried out without parameterizations of eddies and convection, and without including the effect of ice topography and ice-ocean thermodynamics (e.g., Soderlund 2019;Amit et al. 2020;Bire et al. 2022;Kvorka & Čadek 2022).These studies reveal how the ocean circulation and ocean heat transport vary over a wide non-dimensional parameter regime.However, care must be taken when interpreting the real icy moon ocean regime with these simulations because the non-dimensional parameters in the numerical simulations and real icy moon oceans differ by many orders of magnitudes due to numerical constraints (Jansen et al. 2023).More directly comparable to the present work are studies using General Circulation Models (GCMs), where the effect of sub-grid-scale eddies (and in some cases convection) are parameterized (e.g., Ashkenazy & Tziperman 2021;Kang et al. 2022a).In Ashkenazy & Tziperman (2021), the authors use MITgcm to study the ocean circulation and heat transport on Europa considering the effect of salinity and ice-ocean thermodynamics with a flat ice shell.They find that the meridional ocean heat transport is sufficient to compensate for the local imbalance between bottom heating and surface heat loss, thus allowing for a uniform ice thickness.This result is qualitatively consistent with our finding that ocean heat transport will tend to remove any ice thickness variations.However, the ice geometry is flat in the simulations of Ashkenazy & Tziperman (2021) so that there is no ice pump effect, and the planetary parameters of Europa and Enceladus are also significantly different, which makes it difficult to quantitatively compare the results of Ashkenazy & Tziperman (2021) to our work.
In global GCM ocean simulations of Enceladus and other icy moons, there are many processes that are difficult, or computationally impossible, to resolve.For example, slantwise convection cannot be adequately resolved but arises in the form of grid-scale convection in our simulations.Additionally, the effects of tides and librations are not explicitly included in our model -we only consider their effect in providing energy for eddy mixing in the choice for the vertical eddy diffusivity (see Section 2.2 in Zeng & Jansen 2021).Therefore, we need parameterizations for these processes, which introduce external parameters (e.g., eddy diffusivties and viscosities) that are poorly constrained.Our results indicate that at least some of the results are sensitive to these parameterizations.The development of improved parameterizations for these unresolved processes in global icy moon ocean simulations therefore should be a high priority for future work.
We thank Wanying Kang, Edwin S. Kite, and two anonymous reviewers for helpful discussions and comments.This work was completed with resources provided by the University of Chicago Research Computing Center.In our simulations, the vertical resolution is 300-400 m near the boundaries, about 1000 m in the upper ocean, where the stratified layer is formed in a low salinity ocean, and 5000 m in the deep ocean.The varying vertical resolution is designed to save computational resources and to better resolve iceocean boundaries and the stratified layer in the low salinity ocean.The horizontal resolution is 1 • in the zonal and 0.95 • in the meridional direction, which is around 4 km × 4 km near the equatorial surface, decreasing to around 300 m × 3 km near the seafloor at the highest latitudes.
Sub-grid scale processes are parameterized by turbulent viscosities and diffusivities.However, the turbulent viscosities and diffusivities are poorly constrained in simulations for ocean circulation on Enceladus.The vertical turbulent diffusivity in a stably stratified ocean can be constrained by the energy required to mix the water column (Wunsch & Ferrari 2004).On Enceladus, this energy can come from tidal dissipation in the ocean and/or librations, based on which the vertical turbulent diffusivity in the stratified layer is estimated to be anywhere below 3 × 10 −3 m 2 s −1 , although likely much smaller, according to Zeng & Jansen (2021).Here we choose a vertical turbulent diffusivity of 5 × 10 −5 m 2 s −1 so that the estimated depth of the stratified layer is around 14 km and can be explicitly resolved.We further set the vertical viscosity to ν r = 5 × 10 −4 m 2 s −1 to keep the vertical Prandtl number P r = 10 (Soderlund 2019).To ensure numerical stability, we apply a griddependent horizontal viscosity (ν h ), increasing linearly with the square of the zonal grid scale (L 2 x ).The horizontal viscosity maximizes at about 6 m 2 s −1 near the equatorial surface and minimizes at about 0.5 m 2 s −1 near the seafloor at the highest latitudes.We test two different horizontal diffusivities, 0.05 m 2 s −1 and 0.5 m 2 s −1 .These values are motivated by Kang et al. (2020) who estimate the horizontal eddy diffusivity to be κ h = 0.8 m 2 s −1 using a scaling for geostrophically adjusted convection (Jones & Marshall 1993).
The horizontal grid scale is not sufficient to resolve convective plumes.The length scale where rotation becomes important is estimated to be l r ≈ 0.2 m, many times smaller than our horizontal resolution, which means accurately resolving convective plumes is computationally impossible (Zeng & Jansen 2021).We therefore apply a convective adjustment by increasing the vertical diffusivity from the background value κ r = 5 × 10 −5 m 2 s −1 to κ conv =1 m 2 s −1 whenever the stratification becomes unstable to parameterize the effect of convection.The value of κ conv is chosen such that the parameterized convective time scale (τ conv = D 2 o /κ conv ≈ 50 years where D o is the averaged of the ocean) is similar to the convective time scale estimated for the high salinity ocean (see Zeng & Jansen 2021).

A.2. Boundary conditions
We set the ice geometry based on Enceladus's topography and gravitational field as inferred by Čadek et al. (2019).We only apply the Y 20 component of the spherical harmonic function, which means we exclude zonal variations and hemispheric asymmetry in the ice geometry for simplicity (Fig. 1a).
At the bottom of the ocean, we apply a fixed bottom heat flux pattern as thermal boundary condition, following the zonal-mean tidal forcing pattern in Choblet et al. (2017) with a total flux of 20 GW (Fig. 1b).We are applying a smoothly varying bottom heating pattern, although there could be strong local hot spots in the core heating (e.g., Choblet et al. 2017).While this may affect the details of the simulation results, it is not likely to affect the main conclusions.As recently shown by Kang et al. (2022b), warm water from a localized strong heat source will be quickly mixed with the surroundings due to baroclinic instability before it reaches the ice shell.As a result, the small-scale structure of the heating pattern is expected to have relatively little effect on the freezing and melting patterns at the ice-ocean interface.There is no salinity source or sink at the bottom of the ocean.We apply a linear bottom drag with a drag coefficient r b ∼ 10 −4 m s −1 as dynamic boundary condition (see Zeng & Jansen (2021) for a more detailed discussion of the bottom boundary condition).The same linear drag is also applied at the surface.The surface thermal and salinity boundary conditions are determined by the heat and salinity flux due to the ice-ocean interaction.
We apply the "shelfice" package in MITgcm to calculate the ice-ocean interaction (Adcroft et al. 2018).The heat flux balance at the ice-ocean interface is between the conductive heat loss, the latent heat, and the ocean heat flux to the ice-ocean boundary.The conductive heat loss in the shelfice model is assumed proportional to the inverse of the ice shell thickness and the temperature difference between the freezing point and the ice surface temperature, the latter of which is assumed to be spatially constant in our simulations.The magnitude of the conductive heat loss is chosen such as to exactly balance the total heat flux at the sea floor, as is necessary to obtain an equilibrated ocean state.The ocean heat flux to the ice shell is calculated by a turbulent exchange flux through the ice-ocean boundary layer, where the heat flux is assumed proportional to the temperature contrast between the local ocean temperature and the temperature at the bottom of the ice shell (assumed to be the freezing point).The freezing and melting rate is calculated based on the latent heat required to balance the heat budget at the ice-ocean interface.In the shelfice model, the freezing and melting rate is converted to a salinity flux and does not affect the freshwater volume.The volume flux associated with meltwater input is negligible in our simulations: the latent heat is Q LH ∼ O(0.1 W m −2 ) in all simulations (Fig. 5a), so that the volume flux is w q = Q LH /ρL = 3 × 10 −10 m s −1 where ρ is the density of the water and L is the latent heat of fusion.This flux is negligible compared with the typical vertical velocity in the ocean (w o ∼ O(10 −6 m s −1 )).
Our model does not account for the latitudinal variation of the ice surface temperature and the dependence of the heat diffusion coefficient of ice on the local ice temperature (both parameters are assumed constant in our model).We carried out one sensitivity test using the parameters of simulation HS h but with the conductive heat loss through the ice shell computed using a latitudinallyvarying surface temperature and a temperature-dependent heat diffusion coefficient, following Kang et al. (2022c).By comparing the results of this simulation with HS h in the main text, we find that the results do not qualitatively change (Fig. 9).Therefore, we do not expect these variations to qualitatively affect our results, as long as the conductive heat loss is increasing towards the pole, as a result of the polar-thinning ice shell.
It is not known whether Enceladus's ocean and ice shell are in an equilibrium state or not.The adjustment time scale of the ice shell can be estimated as τ i ∼ D i /w q ≈ 2 × 10 6 years where D i = 20 km is the average ice thickness.For comparison, the (parameterized) convective time scale of the ocean is τ conv ≈ 50 years, the overturning time scale of the ocean is τ adv ∼ V o /Ψ ≈ 3×10 5 years where V o is the volume of the ocean and Ψ is the stream function, and the diffusive time scale of the stratified layer in the low salinity ocean is τ dif f ∼ H 2 /κ r ≈ 10 5 years, where H is the thickness of the stratified layer and κ r is the radial diffusivity.We find that the adjustment time scale of the ice shell is much longer than the adjustment time scale of the ocean.As a result, although we do not know if the ice shell on Enceladus is equilibrated, it is more likely that the ocean is near an equilibrium state.

A.3. Integration of the simulations
The ice-ocean coupled 3-D global ocean simulations are difficult to integrate to a fully equilibrated state.Therefore, we present quasi-equilibrium states where the total energy imbalance is less than 3% of the seafloor heating rate.Although there are still temperature and salinity trends in parts of the ocean, we expect that the remaining small imbalance does not qualitatively affect the general results discussed in this paper.The presented results are 500-year averages for HS h , LS h , and LS l , and 4000-year averages for HS l which has larger low-frequency variability.
We start the integration of all simulations with a uniform salinity (35 g kg −1 for high salinity and 8.5 g kg −1 for low salinity).The initial temperature for the high salinity ocean is set to be near the freezing point (about −2 • C) everywhere.For the low salinity ocean, we set the temperature in the upper ocean where the ocean is in contact with the ice shell to be near the freezing point (about −0.6 • C).At the bottom of the ocean, we set the temperature to the critical temperature where the thermal expansivity changes sign and becomes positive (about 1.4 • C).In between, we calculate the depth of the stratified layer using Eq. 1 in Zeng & Jansen (2021) by assuming that the diffusive flux is equal to the total bottom heating, and then prescribe a linear temperature profile to connect the upper and bottom ocean layers, so that the initial conditions are reasonably close to the expected equilibrium state.Small-amplitude random noise is added to generate zonal variability.
This initial estimate still differs significantly from the final solution and the adjustment is very slow in the low salinity ocean.Therefore, we apply an acceleration method for the low salinity ocean simulations.For the temperature, we take the horizontally averaged profile for the temperature tendency (∂T /∂t), and accelerate the simulation by adding this tendency multiplied by a long time period (several thousands of years) to the original temperature profile.We only apply the tendency in layers where it has the same sign as the global temperature tendency: if the global ocean is net warming, we only apply positive temperature tendencies, and vice versa.After each acceleration step is applied, we integrate the model, initialized from the new profile, for about 200 years, and then repeat the acceleration step with the new tendencies.We repeat this procedure until the global temperature tendency changes sign, at which point we continue to integrate the model to a global thermal quasi-equilibrium state (defined as described above).It takes about 20 cycles (several tens of thousands of years) and a final integration of another several thousands of years for the low salinity ocean to reach a quasi-equilibrium state.During the acceleration steps, we did not accelerate the salinity, as the amplification of salinity tendencies from transient freezing and melting does not systematically nudge the simulated solution closer to its equilibrium state.Instead we apply a single adjustment step to the deep ocean salinity once a thermal quasi-equilibrium has been reached.In an equilibrium state, extrema of salinity can only exist at the surface of the ocean where there are salinity sources and sinks, so that every isoline of salinity should start or end at the surface.We hence find the lowermost isoline that intersects with the ice-ocean boundary and set the salinity below this isoline equal to this value.After this adjustment, we run the simulation for an additional 1000 years and use the last 500 years for analysis.Notice that the equilibrated solution is expected to be independent of the spin-up procedure and hence the acceleration method (unless multiple equilibria exist).While our low salinity simulation has not fully reached equilibrium, the remaining drifts are very small, suggesting that the results are not likely to differ qualitatively from the true equilibrium.

A.4. Sensitivity to domain width
In our simulations, we only simulate a zonal domain of 15 • in longitude to save computational resources.However, the effect of horizontal eddies can be important in both ocean heat transport (see Fig. 5 and relevant discussions in the main text) and momentum transport (Ashkenazy & Tziperman 2021).Therefore, we carried out one sensitivity test with the same parameters as in HS h , but with the zonal width of the domain extended from 15 • to 360 • (Fig. 10).To accelerate the simulation, we set the initial condition to be the zonal-mean field of the quasi-equilibrium state of HS h with small random perturbations added to the temperature field.This simulation has been run for 120 years, at which point a quasi-equilibrium state is reached where the total heat imbalance is smaller than 7% of the seafloor heating rate.The results of this simulation are qualitatively similar to those in HS h , in particular showing a very similar freezing and melting pattern.where r b is the radius at the bottom of the ocean, r oi is the radius at the ice-ocean interface, v is the meridional velocity, T is the temperature, λ is longitude, λ x is the longitudinal extent of the domain, and the overbar indicates a time-average.The time-and zonally-averaged advection can be written as where the bracket denotes the zonal average and primes denote deviations from the zonal and temporal average -which we will refer to as the eddy component.Eq.B1 can then be written as If we assume that the eddy heat transport is down-gradient and is proportional to the mean gradient of temperature, we can apply an equivalent "eddy diffusivity" κ e to replace the eddy term: [v ′ T ′ ] = −κ e ∂[T ]/r∂ϕ (Vallis 2017).Then Eq.B2 becomes where on the right hand side the first term indicates heat transport by the mean overturning circulation, while the second term indicates heat transport by resolved eddies and diffusion (which may be thought of as representing the effect of unresolved eddies).

B.2. Ocean heat budget of the surface layer
To better understand how the parameterized upright convection is connected with the surface freezing and melting pattern, we compare the vertical diffusive heat flux near the ice-ocean interface (at the bottom of the model layer in contact with the ice shell) with the ocean heat flux to the ice shell Q oi (Fig. 11).The parameterized vertical convective diffusivity (κ conv = 1 m 2 s −1 ) is much larger than the background vertical diffusivity (κ r = 5 × 10 −5 m 2 s −1 ), so that the vertical diffusive heat flux primarily shows the effect of the parameterized upright convection.In all simulations, near the polar surface, Q oi approximately matches the parameterized upright convective heat flux, indicating that the convection brings up colder water from below and supports freezing near the poles.

C. MULTIPLE EQUILIBRIA IN HIGH SALINITY SIMULATIONS
We have found a positive feedback which can result in spontaneous symmetry breaking in the high salinity simulations.This mechanism, however, should not have a preference between the northern and southern hemisphere.We hence carried out one sensitivity test with the same parameters as the HS l simulation but using the "mirrored" salinity field of the quasi-equilibrium state in HS l as the initial condition.In this sensitivity test, we observe that the cross-equatorial overturning circulation has the opposite sign.Consistently, we find stronger melting outside the tangent cylinder at low latitudes in the southern hemisphere (in contrast to the northern hemisphere in HS l ), and we still see the switching between an asymmetric state and a more symmetric state (Fig. 12).The nonlinear dynamics controlling the multiple quasi-equilibria could be an interesting topic for future studies.

Figure 1 .
Figure 1.Prescribed ice topography (a), corresponding conductive heat loss (b, blue line), and bottom heating (b, red line) per unit area, as a function of latitude.The global integral of bottom heating balances the conductive heat loss.The apparent offset between the blue and red lines in panel (b) is due to the different surface areas at the bottom and ice-ocean interface.

Figure 3 .
Figure 3. Thermal wind balance in the quasi-equilibrium state of all simulations.(a) shows 2 ⃗ Ω • ∇[u], and (b) shows −∂[b]/r∂ϕ, respectively, where [•] denotes a zonal average.The first to the fourth columns are results of the high salinity simulation with high horizontal diffusivity (HS h ), high salinity simulation with low horizontal diffusivity (HS l ), low salinity simulation with high horizontal diffusivity (LS h ), and low salinity simulation with low horizontal diffusivity (LS l ), respectively.

Figure 4 .
Figure 4. Zonal-mean zonal velocity (u, a), stream function (Ψ, with clockwise circulation defined as positive, b), and zonal-mean vertical kinetic energy (w 2 , in logarithmic scale, c) in the quasi-equilibrium state of all simulations.The first to the fourth columns are the results of the high salinity simulation with high horizontal diffusivity (HS h ), high salinity simulation with low horizontal diffusivity (HS l ), low salinity simulation with high horizontal diffusivity (LS h ), and low salinity simulation with low horizontal diffusivity (LS l ), respectively.The stream function is calculated as Ψ = − λx 0 r r b vr cos ϕdλdr where λ x = 15π/180 is the longitudinal extent of the domain, v is the meridional velocity, λ is longitude, and r is radius.The black dashed lines denote the boundary of the tangent cylinder (a cylinder that is parallel to the rotational axis and is tangent to the bottom core).

Figure 5 .Figure 6 .
Figure5.Zonal-mean surface heat flux (a) and zonally-and vertically-integrated meridional ocean heat transport (b, see APPENDIX B.1 for details) in the quasi-equilibrium state of all simulations.The first to the fourth columns are results of the high salinity simulation with high horizontal diffusivity (HS h ), high salinity simulation with low horizontal diffusivity (HS l ), low salinity simulation with high horizontal diffusivity (LS h ), and low salinity simulation with low horizontal diffusivity (LS l ), respectively.In (a), the black dashed line is the bottom heating (Q b ), the blue line is the latent heating (Q LH ), the red line is the conductive heat loss (Q cond ), and the black solid line is the surface ocean-to-ice heat flux (Q oi ).In (b), the blue line is the diffusive heat flux (F Dif ), the red line is the advective heat flux (F Adv ), and the black line is the total heat flux (F Q ).

Figure 7 .
Figure 7. Temporal variability in the high salinity simulation with low horizontal diffusivity (HS l ).Zonalmean latent heating (a), zonal-mean surface temperature (b), zonal-mean surface salinity (c, a reference salinity of 34.97 g kg −1 is subtracted), zonal-mean surface density (d, a reference density of 1029.13 kg m −3 is subtracted), and stream function at r = 225.375km (e) as a function of time and latitude.The black dashed lines are the latitudes where the tangent cylinder intersects the ocean surface in (a)-(d) and the latitudes where the tangent cylinder intersects r = 225.375km in (e).The first 4000-year average is used as the quasi-equilibrium state of the simulation HS l .

Figure 8 .
Figure 8. Snapshots of meridional eddy velocity v ′ (deviations from time-and zonal-mean) in all simulations at r = 225.375km.The first to the fourth columns are results of the high salinity simulation with high horizontal diffusivity (HS h ), high salinity simulation with low horizontal diffusivity (HS l ), low salinity simulation with high horizontal diffusivity (LS h ), and low salinity simulation with low horizontal diffusivity (LS l ), respectively.Note that the magnitude of v ′ in (a) is 10 4 times smaller than in (b-d).
Parameterization of sub-grid scale processes

Figure 9 .
Figure 9. Simulation with parameters in HS h but including surface temperature variation and temperaturedependent heat diffusion coefficient of the ice when prescribing the conductive heat loss (following Kang et al. 2022c).(a): zonal-mean zonal velocity u. (b): stream function Ψ, with clockwise circulation defined as positive.(c): surface heat budget, where the black dashed line is the bottom heating (Q b ), the blue line is the latent heating (Q LH ), the red line is the conductive heat loss (Q cond ), and the black solid line is the surface ocean-to-ice heat flux (Q oi ).(d): zonal-mean temperature T (the contour interval is 0.01 • C below -2.14 • C and 0.001 • C above -2.14• C). (e): zonal-mean salinity S (a reference salinity of 34.97 g kg −1 is subtracted, and the contour interval is 10 −3 g kg −1 ).(f): zonal-mean potential density ρ σ (a reference density of 1029.13 kg m −3 is subtracted, and the contour interval is 10 −3 kg m −3 below 5 × 10 −3 kg m −3 and 3 × 10 −4 kg m −3 above 5 × 10 −3 kg m −3 ).

Figure 10 .
Figure 10.Simulation with the same parameters as in HS h but with the zonal width of the domain extended from 15 • to 360 • .(a): zonal-mean zonal velocity u. (b): stream function Ψ, with clockwise circulation defined as positive.(c): surface heat budget, where the black dashed line is the bottom heating (Q b ), the blue line is the latent heating (Q LH ), the red line is the conductive heat loss (Q cond ), and the black solid line is the surface ocean-to-ice heat flux (Q oi ).(d): zonal-mean temperature T (the contour interval is 0.01 • C below -2.14 • C and 0.001 • C above -2.14• C). (e): zonal-mean salinity S (a reference salinity of 34.97 g kg −1 is subtracted, and the contour interval is 10 −3 g kg −1 ).(f): zonally-and vertically-integrated meridional ocean heat transport.The blue line is the diffusive heat flux (F Dif ), the red line is the advective heat flux (F Adv ), and the black line is the total heat flux (F Q ).Note that the color bar range in (b) and the vertical axis range in (f) are multiplied by 24 compared to Fig. 4b1 and Fig. 5b1 to account for the larger zonal width of the domain and allow for easier comparison.

Figure 11 .
Figure 11.Vertical heat fluxes near the surface of the ocean (the uppermost ocean model layer in contact with the ice shell).The blue lines show the radial diffusive heat flux (including the parameterized convective heat flux) at the bottom of the surface ocean model layer, and the black line shows the ocean heat flux to the ice shell (Q oi ).The thinner lines are the raw data and the thicker lines are moving-averages of 9 grids to smooth grid-scale variabilities.The first to the fourth columns are results of the high salinity simulation with high horizontal diffusivity (HS h ), high salinity simulation with low horizontal diffusivity (HS l ), low salinity simulation with high horizontal diffusivity (LS h ), and low salinity simulation with low horizontal diffusivity (LS l ), respectively.

Table 1 .
Parameters for Enceladus simulations

Table 2 .
Experimental Design