Different ice shell geometries on Europa and Enceladus due to their different sizes: impacts of ocean heat transport

On icy worlds, the ice shell and subsurface ocean form a coupled system -- heat and salinity flux from the ice shell induced by the ice thickness gradient drives circulation in the ocean, and in turn, the heat transport by ocean circulation shapes the ice shell. Therefore, understanding the dependence of the efficiency of ocean heat transport (OHT) on orbital parameters may allow us to predict the ice shell geometry before direct observation is possible, providing useful information for mission design. Inspired by previous works on baroclinic eddies, I first derive scaling laws for the OHT on icy moons, driven by ice topography, and then verify them against high resolution 3D numerical simulations. Using the scaling laws, I am then able to make predictions for the equilibrium ice thickness variation knowing that the ice shell should be close to heat balance. Ice shell on small icy moons (e.g., Enceladus) may develop strong thickness variations between the equator and pole driven by the polar-amplified tidal dissipation in the ice, to the contrary, ice shell on large icy moons (e.g., Europa, Ganymede, Callisto etc.) tends to be flat due to the smoothing effects of the efficient OHT. These predictions are manifested by the different ice evolution pathways simulated for Enceladus and Europa, considering the ice freezing/melting induced by ice dissipation, conductive heat loss and OHT as well as the mass redistribution by ice flow.


INTRODUCTION
Many of the icy satellites in the outer solar system are likely to contain a subsurface ocean underneath their ice shell due to tidal dissipation (Scharf 2006), which may lead to a suitable environment for life to thrive. Enceladus (a satellite of Saturn) and Europa (a satellite of Jupiter), in particular, have been confirmed to have a global subsurface ocean by data brought back by the Galileo and Cassini missions (Postberg et al. 2009;Thomas et al. 2016;Carr et al. 1998;Kivelson et al. 2000;Hand & Chyba 2007). As two of the most enigmatic targets to search for extraterrestrial life (Des Marais et al. 2008;Hendrix et al. 2019), Enceladus and Europa are to be further explored in the future (e.g., Europa Clipper and JUICE ). Thus far, the measurements and detection have be carried out above the ice shell for the most part, and are likely to remain that way in the near future. Therefore, it becomes crucial to infer the subsurface condition using the information collected above the surface and to put better constraints on the ice shell geometry so we know where to send our landing or drill missions.
These motivate us to study the interaction between the subsurface ocean and the ice shell that we can more easily take measurement of. As illustrated by Kang et al. (2021) and , the interaction happens in a mutual way: the variation of ice thickness on Enceladus can drive ocean circulation by inducing salinity flux through freezing/melting and by changing the local freezing point; in turn, ocean circulation can converge heat to regions covered by relatively thick ice, flattening the ice shell (sketched in Fig. 1d). Such interaction makes it possible to infer the subsurface ocean properties from the information about the ice shell geometry or to infer the ice shell geometry based on our understanding of ocean dynamics .
Thus far, based on the surface measurements (libration, shape, gravity etc.) done by Cassini, Enceladus' ice shell has been revealed to be around 20 km thick on global average and to become significantly thinner toward the poles (Iess et al. 2014;Beuthe et al. 2016;Tajeddine et al. 2017;Čadek et al. 2019;Hemingway & Mittal 2019). Whilst the equatorial ice shell is 30 km thick, the ice shell thickness over the south pole is only 6 km (Hemingway & Mittal 2019); all geysers gather around this area, opening up periodically under tidal stress (Hedman et al. 2013;Hurford et al. 2007;Nimmo et al. 2014;Ingersoll et al. 2020) and ejecting samples of the ocean to outer space (Porco et al. 2006;Hansen et al. 2006;Howett et al. 2011;Spencer et al. 2013). In order to sustain the strong ice topography, the ocean heat transport (OHT) that flattens the ice shell through ice pump mechanism (Lewis & Perkin 1986) cannot be arbitrarily strong, which in turn puts constraints on the ocean salinity and the partition of heat production between the silicate core and the ice shell (Kang et al. 2021).
Europa's ice shell geometry is not as well constrained, but evidence has been found in favor a relatively thin (<15 km Hand & Chyba 2007) and flat (Nimmo et al. 2007) ice shell, in line with the evidence that Europa geysers are not as concentrated (Roth et al. 2014;Jia et al. 2018;Arnold et al. 2019;Huybrighs et al. 2020). As a separate line of evidence,  found that the OHT-induced per-area heat flux scales with the satellite's radius to the power of 0.5 or 1 (depending on whether the magnitude of vertical diffusivity is sufficient to communicate the entire ocean column from the ice to the seafloor), which also supports a rather flat ice shell on Europa given its large size.
However, what has been missing in the framework proposed by Kang et al. (2021) and  are eddy dynamics. We know from earth ocean that ocean can be filled by eddies of various scales and that they play an important role in transporting heat, momentum and tracers (Thompson 2008;Volkov et al. 2008; Thompson et al. 2014). Scaling laws that govern the eddy diffusivity, tracer transports and equilibrium stratification have been found under forcings that are relevant for earth ocean or atmosphere (Held & Larichev 1996;Karsten et al. 2002;Jansen & Ferrari 2013). In the context of icy satellites, recent works by Ashkenazy & Tziperman (2021) and Ashkenazy & Tziperman (2016) demonstrate that eddies also exist in the subsurface ocean of snowball Earth and Europa, and a set of sensitivity tests run under 3D configuration in Kang et al. (2021) also shows strong baroclinic eddies. Finding scaling laws for the eddy heat transport in icy moon oceans and making predictions about the equilibrium ice shell geometry are the main goals of this work.
The works on the eddy transports across the the antarctic circumpolar current (ACC) on earth provide useful insights for the OHT scaling on icy moons. However, a few differences between icy moon oceans and the terrestrial ocean should be noted. The overturning circulation near the ACC is largely driven by surface wind stress especially: curl of wind stress forces meridional overturning flows, maintaining isopycnal slope so that eddies can grow on it (McWilliams et al. 1978;Marshall & Radko 2003;Plumb & Ferrari 2005). To the contrary, the ocean on an icy satellite is sandwiched between an ice shell and a silicate core, and thus experience no wind stress. Density gradients created by the ocean-ice heat and salinity exchange, instead, drive the overturning circulation there. Since both the top and bottom boundaries are frictious on icy satellites, they can provide the drag required to balance the Coriolis acceleration associated with the overturning flows along the boundaries . The equilibrium temperature and salinity profiles adjusted only by the zonally-symmetric overturning circulation may undergo baroclinic instability (Charney 1947), which further boosts the OHT.
2. THE COUPLED OCEAN-ICE SYSTEM.
The system considered by this study is sketched in Fig. 1 -a 56-km deep ocean covered by an ice shell that is about 20 km thick. A poleward-thinning ice shell is assumed, given the fact that the tidal dissipation in the ice shell amplifies over the poles (Beuthe 2019).
where H 0 = 20 km is the mean ice thickness, P 2 is the 2nd order Legendre polynomials, and φ denotes latitude. Here, the ice thickness variation is assumed to follow the P 2 profile for simplification, and unless otherwise mentioned, the equator-to-pole thickness difference ∆H is set to 3 km. This default ice thickness profile is shown by the solid curve in Fig. 1b. One would expect ice shell to be thinner over the poles as can be seen in Fig. 1b, because the tidal dissipation produced in the ice is stronger there (Beuthe 2018(Beuthe , 2019Kang & Flierl 2020). In a situation where the ice shell is thinner over the equator (∆H < 0), the results here can still apply after reversing the sign of the circulation and heat transport. Also, we expect the qualitative results obtained in this work to hold when the ice thickness variation follows a profile other than P 2 , as long as the ice thickness is relatively simple (poleward-thinning or poleward-thickening) and is symmetric about the equator. . Panel (a) sketches the primary sources of heat and heat fluxes, which include: heating due to tidal dissipation in the ice Hice, the heat flux from the ocean to the ice Hocn and the conductive heat loss to space H cond . Ocean heat transport is shown by the horizontal arrow. Panel (b) shows the default ice shell thickness profile considered here a black solid curve, which is thinner over the poles because ice dissipation amplifies going poleward (Beuthe 2019). The gray dashed curve shows the freezing (positive) and melting rate (negative) required to maintain a steady state based on an upside-down shallow ice flow model (see appendix for details). In this calculation, the default 2500 km radius is considered. Panel (c) shows the profiles of Hice, H cond and H latent given the information in panel (b). Panel (d) sketches the key physical processes in an ocean covered by an ice shell with varying thickness (see main text for description). Panel (e) shows how thermal expansion coefficient under the ice shell varies with the satellite's size (gravity), at 10 psu (blue) and 60 psu (brown) ocean salinities. Panel (f) shows the salinity forcing (equatorial minus polar salinity flux, dots) and the temperature forcing (the freezing point difference under the equatorial and polar ice shell, crosses) as a function of the moon's radius.
With ice flowing from thick ice regions to thin ice regions 1 , this equator-to-pole ice thickness gradient won't last, unless freezing/melting can constantly enhance the poleward-thinning ice topography. The ice freezing/melting, in turn, is governed by the ice shell heat budget. Driven by the hundreds of degree temperature difference between the water-ice interface and the ice surface, the ice constantly loses heat in form of heat conduction. On Enceladus and Europa, where the ice shell is likely around 20 km thick, the heat loss rate H cond is around 40 mW/m 2 on global average, and is faster over regions where ice is thin and over the poles where the ice surface temperature is low, as shown by the green curve in Fig. 1c. To balance the heat loss, the ice shell and the silicate core needs to produce heat, however, the core-shell heat partition is poorly understood (Beuthe 2019;Choblet et al. 2017). In this work, I will focus on the shell-heating scenario and discuss the potential impacts of bottom heating only toward the end. The ice heat production due to tidal flexing (denoted by H ice ) amplifies over the polar regions even if the ice is completely flat (Beuthe 2018), and this polar-amplifying pattern is further enhanced by the poleward thinning ice geometry through the ice rheology feedback, which concentrates heat production toward regions with thinner ice shell and hence weaker mechanical strength (Beuthe 2019). Shown by a red curve in Fig. 1c is H ice with the default ice thickness profile (Eq. 1). Besides, the equatorial freezing and polar melting required to maintain the prescribed ice thickness gives rise to almost negligible latent heat release H latent , shown by a gray curve in Fig. 1c. At the bottom of the ice, there may also be heat exchange with the ocean (denoted by H ocn ), manifested by ocean circulation and heat transport. In an equilibrium state, the ice shell heat budget needs to be in balance, which states as the vanishing of the sum of all the heating terms.
From ocean's perspective, however, it is not the heat budget of the ice shell that matters, but the heat and salinity fluxes from the ice. These fluxes can be derived, as long as the ice geometry is given. In direct contact with ice, the ocean temperature at the water-ice interface will be relaxed toward the local freezing point, which is lower under a thick ice shell because of the high pressure (see Eq. A2 in the appendix). Also, assuming ice shell is in mass equilibrium, equatorial freezing and polar melting is required to prevent the ice shell from being flattened by the pressure-driven ice flow (see the dashed gray curve in Fig. 1b). The freezing/melting will then induce salinity exchange with the subsurface ocean. Under these forcings, water over the poles become warmer and fresher than the water at low latitudes. The resultant density variations drive ocean circulation and eddies, transporting heat down-gradient from the poles to the equator, which in turns affect the heat budget of the ice shell, leading to changes in the ice geometry. Therefore, the critical task for predicting the ice shell geometry is to determine the dependency of the ocean heat transport (OHT) on the ice shell geometry via their heat and salinity fluxes under various satellite parameters (such as ocean salinity, gravity, size etc.).
In section 4 and section 5, I will derive scaling laws for the dependency of OHT on the ice shell geometry under various satellite parameters, test them against numerical simulations. Since we already know how the other heating terms (conductive heat loss, ice dissipation and latent heating) depends on ∆H, this generic formula for OHT as a function of ∆H and the satellite parameters (such as ocean salinity, gravity, size etc.) allows one to solve for the equilibrium ∆H for a specific icy satellite under the condition that the ice shell heat terms should balance one-another. This is done in section 6. It should be noted that, throughout this work, the heating produced in the silicate core is assumed to be zero and all heat is assumed to be generated in the ice shell. The impacts of core heating is discussed in the conclusion.
More details about the ocean circulation model, ice flow model, and the tidal dissipation model can be found in the appendix.

WHY SIZE MATTERS?
Using parameters relevant for Enceladus, Kang et al. (2021) show that the circulation that arises from the surface heat and freshwater forcing can go either direction depending on the ocean salinity: in the low-salinity limit, temperatureinduced density variation dominates, and the warm polar water would sink as sketched by the blue arrow in Fig. 1d because fresh water contracts upon warming (anomalous expansion); whilst in the high-salinity limit, the anomalous expansion is suppressed, and both salinity-and temperature-induced density gradients contribute to downwelling at low-latitudes, as sketched by the orange arrow in Fig. 1d.
When considering icy satellites larger than Enceladus (most of the icy satellites of interest are), the following changes are expected: • The thermal expansion coefficient will become more positive, and eventually anomalous expansion will be completely suppressed even if the ocean is relatively fresh. Shown in Fig. 1e are the dependence of thermal expansion coefficient under the ice shell on the icy moon's radius, a, for two different salinities, 10 psu and 60 psu. Anomalous expansion doesn't occur in an ocean with 60 psu salinity regardless of the size of the satellite, yet it does occur in a 10-psu ocean, but only when a < 1500 km 2 , approximately the size of Europa. Therefore, with large enough planetary size, downwelling will only occur over the equator regardless of the ocean salinity.
• The temperature difference under the ice shell between the equator and the pole ∆T will increase, as shown by the cross markers in Fig. 1f. The ocean temperature adjacent to the ice shell should be close to the local freezing point (denoted by T f ), and T f decreases with pressure linearly.
where ∆H is the prescribe ice thickness difference between the equator and the poles, g 0 is the moon's surface gravity, b 0 = −7.61 × 10 −4 K/dbar and ρ i is the ice density. With ∆H fixed, ∆T ∝ g ∝ a.
• Salinity forcing will weaken. The blue and brown dots in Fig. 1f show the salinity flux (mean salinity S 0 times the freezing rate q) difference between the equator and the poles for an ocean with mean salinity of 10 psu and 60 psu, respectively. The freezing rate q is set to balance the divergence of ice flow. As derived in the appendix A4, ice flow behaves like diffusion, and the flow divergence/convergence is proportional to ∆P divided by the distance square a 2 , so q ∝ ∆P/a 2 ∝ a −1 . (3) • The same density gradient will drive a stronger ocean circulation and heat transport as a result of the stronger gravity -fixing the bulk density of the satellite, surface gravity g 0 ∝ a. The stronger heat transport will then flatten the ice shell more efficiently.
Given the above reasoning, one can see that, on larger icy moons, the OHT is likely to be 1) more efficient in flattening the ice shell and 2) dominantly controlled by temperature variations. Assuming the ice thickness varies by 30% meridionally and a melting-point ice viscosity of 10 14 Pa·s, a back-of-the-envelope calculation will show that the salinity-induced density anomaly is only comparable to the temperature-induced one if the moon's size is smaller or comparable to that of Enceladus. In , the authors have demonstrated these points in a zonally symmetric framework, ignoring eddy transports. Their results show that the meridional OHT should scale with the moon's radius a, the equator-to-pole ice thickness difference ∆H and the Coriolis coefficient f as follows, when the circulation depth is limited by vertical diffusion, or when the circulation reaches the seafloor.

SCALING LAWS FOR OCEAN HEAT TRANSPORT (OHT) CONSIDERING EDDIES.
Given that salinity forcing tends to be dominated by temperature forcing unless the size of the icy moon is comparable or smaller than Enceladus (see the previous section), only the temperature-induced density anomalies is considered here, following . For smaller icy satellites, the salinity-driven circulation may add to or cancel out the temperature-driven circulation depending on the ocean salinity, and as a result, the OHT could be off by one order of magnitude (Kang et al. 2021).
The eddy heat transport F e can be represented by an equivalent diffusive process, where κ e is the equivalent diffusivity of baroclinic eddies and d denotes the depth that the density anomaly penetrates downward from the surface. The corresponding residual circulation can be written as where ∆ v T is the vertical temperature contrast across the depth d.
According to Held & Larichev (1996), κ e can be estimated by where k = 0.25 is a constant, V is the rms eddy velocity and L β denotes wavelength of the energy-containing eddies, which follows the Rhines scale This energy-containing wavelength is typically larger than the deformation radius L d , the scale at which baroclinic instability happens: f denotes Coriolis coefficient, and N = gρ z /ρ 0 denotes the Brunt-Vasala frequency. Here, the vertical density gradient ρ z is estimated by ρ 0 α∆ v T divided by d.
According to the scaling proposed by Held & Larichev (1996), the ratio L β /L d and V /U (U denotes the zonal jet speed) is governed by the supercriticality ξ, where β ∼ f /a is the meridional gradient of f . In the above equation, d denotes the depth, to which circulation and dynamics can penetrate. The thermal wind balance connects U with d, A second constraint can be provided by the balance of vertical heat transport. In an equilibrium state, the horizontally integrated vertical diffusion of heat ρC p κ v (∆ v T /d)(2πa 2 ) should be equal to the downward heat transport by residual circulation C p Ψ † ∆T (Jansen & Ferrari 2013) and that yields κv-limit. -I first consider the κ v -limited situation, where, because of the low vertical diffusivity, the densest isentrope initiated from the equator bends over and reaches the poles without intersecting with the bottom. In this scenario, d < D, ∆ v T = ∆T , and ξ ∼ 1. As a result, V and U , L d and L β are interchangeable 3 , so κ e can be written as The depth of density variations d is constrained by vertical diffusivity κ v through Eq. (13), leading to which indicates that the depth of density variation is insensitive to the moon size, but increases with the rotation rate and decreases with ice thickness gradients. To make sure that the above solution is consistent with the assumptions, one needs to check whether d is indeed smaller than the ocean depth D. If this is not satisfied, the system should be in the D-limit instead. If d indeed turns out smaller than D, then we can substitute Eq. (15) and Eq. (14) into Eq. (6), which gives To obtain the proportionality, the facts that ∆T ∝ ∆P ∝ g∆H ∝ a∆H, g ∝ a and β ∝ f /a are used. The dependence of F κv on a, ∆H are very similar to that of the κ v -limited overturning circulation in the zonally symmetric case given by  (see Eq. 4). Scaling laws can also be obtained for L d and L β , which should characterize the size of eddies and the width of jets, and for the jet speed U In the above equations, subscript κ v indicates that this is the κ v -limit scaling. Noticeably, the above scaling laws predict eddy size and jet width to grow linearly with moon's radius a, meaning the dominant wavenumber and the number of jets will be insensitive to the moon's radius, but the jets will be stronger (U ∝ a). Besides, slower rotation and stronger ice thickness gradient will make the eddies and jet stronger and larger in size.
D-limit. -For those scenarios, where isentropes intersect with the seafloor, d = D, ∆ v T < ∆T , and ξ > 1 (supercritical). The equivalent eddy diffusivity should then be written as From the above equation and the vertical heat balance (Eq. 13), ξ can be solved Substituting the Eq. (19) and Eq. (20) into Eq. (6) yields the eddy heat transport The above scaling exhibits a stronger dependence on ∆H and f than the ξ = 1 scaling, and it becomes dependent on the ocean depth D because the vertical diffusion now is strong enough to communicate the bottom and the ice shell. Again, this scaling is comparable to the 2D scaling obtained by Kang & Jansen (2022) (see Eq. 5).
By assumption, the penetration depth of the T/S anomalies is the entire ocean depth D, i.e., Substituting Eq.22 and Eq.20 into the definitions of deformation radius and Rhines scale (Eq. 9, Eq. (11), Eq. 10), the following scalings are obtained: Just as in the κ v -limit scaling, the dominant eddy wavenumber and jet number should be insensitive to the radius of the moon since the power of the a factors for L d,D and L β,D are both equal to one. Also, slower rotation and stronger ice thickness gradient will enhance the jet speed and increase the sizes of jets and eddies.
Combining the two scalings together, the final OHT should be equal to the lower value between the two scalings given by Eq. (16) and Eq. (21) because both ocean depth and vertical diffusivity constrain how much heat can be transported by the ocean.

EXAMINE THE THEORY USING NUMERICAL SIMULATIONS.
Shown in Fig. 2(a-c) is the predicted F ocn as a function of the icy moon's size, rotation rate and equator-to-pole thickness difference in highly saturated colors. Blueish colors denote fresh ocean and reddish colors denote salty ocean. Default parameters can be found in Table 1. For comparison, the corresponding 2D scalings given by  are shown in lighter colors. To verify these scaling laws, I ran 2D and 3D numerical simulations varying the three parameters (a, f and ∆H) and diagnose the OHT from the model output. The 2D and 3D results are marked on Fig. 2 using dots and diamond markers, respectively. They match the prediction (Eq. 26) up to a factor of 2. The  (a) to (c) bottom shows the dependence on the moon's size a, the rotation rate reflected by the Coriolis coefficient f , and on the equator-to-pole ice thickness difference ∆H. The lines in highly saturated colors present the 3D scaling given by Eq. (16) and Eq. (21), and lines in lighter colors present the 2D scaling given by . Scattered on top are the diagnosed OHT from 3D numerical experiments (diamond markers) and 2D numerical experiments (dots). Different colors are used to differentiate different ocean salinities: from blueish color to reddish color, salinity increases. Default parameters used in the scaling and the numerical experiments can be found in Table 1.
matching seems particularly good in panel (a), because the range of F ocn there is much wider so that the offset appears smaller in comparison. The OHT is governed by the ocean's dynamic and thermodynamic states. Presented in Fig.3 are the model solutions obtained under 2D and 3D configurations. In these simulations, moon's radius a = 2500 km, ocean salinity S = 60 psu, equator-to-pole ice thickness difference ∆H = 3 km, and Europa's rotation period (3.5 day) is used. Due to the pressure gradient induced by the poleward thinning ice geometry, the freezing point under the ice is higher over the poles compared to the equator, driving the poleward warming pattern seen in Fig.3-a. Meanwhile, in order to sustain the ice geometry against the flattening due to ice flow (Eq.A7), equatorial ice needs to freeze and polar ice shell needs to melt. This in turn drives a meridional salinity gradient seen in Fig.3-b. High ocean salinity and high pressure suppress the anomalous expansion behavior (contract upon warming) near the freezing point, which typically happens on small icy moon with fresh ocean. As a result, temperature and salinity anomalies both contribute to the high density in low latitudes, driving sinking motions there (see Fig.3-e). The circulation is forced to be mostly aligned with the direction of rotation in the interior. This is because any motions moving closer or away from the rotating axis will lead to eastward/westward acceleration by virtue of angular momentum conservation, and the resultant zonal flows will be too strong to be in thermal-wind balance with the weak density variation. Flows across the direction of the rotation axis concentrate near the two rough boundaries at the top and bottom, because only in appearance friction, radial flows can be sustained without inducing imbalance in the zonal momentum budget. In an equilibrium state, the upper part of the ocean flows westward and the lower part of the ocean flows eastward (Fig.3-d), in thermal wind balance with the density distribution ( Fig.3-c).
The major differences between 3D and 2D configurations lie in the zonal flow field (Fig.3-d): there are jets formed in 3D due to the Renold stress associated with the baroclinic eddies, whose structure is presented in Fig.3-f,g. Due to the strong rotation effect, the eddy motions tend to be aligned with the direction of rotation axis, forming convective "rolls" along the equatorial plane and wave-like structures in higher latitudes, consistent with recent icy moon studies Ashkenazy & Tziperman (2021); Soderlund et al. (2014); Bire et al. (2022); . These eddies facilitate stronger equatorward heat transport than the overturning circulation in 2D configuration, as shown by Fig.3h. It should be noted that some previous works have jets and Taylor columns even under 2D configuration (Ashkenazy & Tziperman 2021). This difference arises from the different model configuration: in their work, the ocean is heated from below and no-slip boundary condition is applied to the top and the bottom; whereas in this work, the bottom heating is assumed to equal to zero and the boundary drag is parameterized by a linear drag toward zero.
At different the moon's sizes, rotation rates and the ice thickness variations, the dynamics remain qualitatively the same, as shown by Fig.6-Fig.11 in the appendix. On a smaller icy moon, the meridional temperature gradient weakens due to its weaker gravity (compare Fig.6a and Fig.7a with Fig.3a). That in turn weakens the circulation (see panel e1) and eddies (see panel f,g), leading to a much weaker heat transport (see panel h) as suggested by Eq.16 and Eq.21. lat a = 2500 km, 2D Figure 3. Ocean circulation and thermodynamic state under 2D and 3D configurations. The left column (panels a1-e1) shows temperature, salinity, density, zonal flow and meridional streamfunction from a zonally symmetric 2D simulation. The second column (panels a2-e2) shows the same thing for the 3D simulation. Panel (f) and panel (g) present the temperature and zonal flow anomalies from the zonal mean in a plane view. Panel (h) presents the vertically and zonally integrated meridional ocean heat transport diagnosed from the 2D (dashed) and 3D model (solid). In this defualt setup: moon's radius a = 2500 km, ocean salinity S = 60 psu, equator-to-pole ice thickness difference ∆H = 3 km, and Europa's rotation period (3.5 day).
According to Eq. (17) and Eq. (24), the number of jets remain more of less unchanged -this can be seen from the panel-d2 of Fig.6, Fig.7 and Fig.3. Fast rotation suppresses the circulation, eddies and thereby OHT, as suggested by both 2D scalings (Eq.4 and Eq.5) and 3D scalings (Eq.16 and Eq.21). This trend is reflected by sensitivity tests shown in Fig.8 and Fig.9, where a shorter rotation period of 1.37 days (Enceladus' rotation period) and a longer rotation period of 10 days are employed, respectively. Also, when rotation rate varies, the width of jets and eddies should change accordingly. Strong rotation leads to smaller deformation radius L d = (N H/f ) (governing the eddy size, Eq.10) and smaller Rhines scale L β = U/β (governing the jet width, Eq.9), consistent with panel (d,f,g) of Fig.8, Fig.3 and Fig.9. The dependence of Kang OHT on the ice thickness variation ∆H is quite intuitive. When the ice is flatter, the temperature/salinity variations, eddy amplitude, eddy size and heat transport all decrease, as shown by Fig.10 and Fig.11.
As a further verification of the theory, the jet speed is diagnosed from the simulations by subtracting the averaged U from the maximum zonal mean zonal flow, and is compared against the predictions given by Eq. (18) and Eq. (25). The averaged spacing between jets is identified by 'numpy.find_peaks', and is compared against Eq. (17) and Eq. (24). Finally, the prominent eddy size is computed by averaging different wavelength by its corresponding power spectrum of U field between 40N/S and 60N/S, and the results are compared against Eq. (17) and Eq. (24), because the eddies' energy-containing scale follows the Rhines scale (Held & Larichev 1996). All these diagnostics are carried out in high-latitudes inside the tangent cylinder -a cylinder whose sides are parallel to the moon's axis of rotation and are tangential (hence the name) to the ocean's floor at the equator, because the dynamics outside the tangent cylinder are very different Bire et al. 2022). As shown in Fig.12, the theory also catch the eddy and jet characteristics reasonably well, indicating that the agreement between the predicted and simulated F ocn may not be a coincidence.
6. THE EQUILIBRIUM EQUATOR-TO-POLE ICE THICKNESS GRADIENT.

Scaling Theory
The dependence of F ocn on orbital parameters and ∆H can be used to predict the equilibrium ice thickness variation using the fact that the ice shell heat budget should be closed. This analysis has been done by  but using the 2D scaling for F ocn . Here, I repeat the process for 3D scalings.
First, I need to convert F ocn to heat flux anomaly received by the ice shell. Assuming that the heat transported from the polar regions to the equatorial regions by the ocean is evenly distributed over a half hemisphere with a surface area of πa 2 , the heat flux per area received by the equatorial ice shell and the heat flux leaving the polar ice shell equal to (αg∆T ) 6/7 κ 3/7 v ∝ a∆H 13/7 f −8/7 for D-limit (27) In an equilibrium state, the ice shell needs to be in a heat balance everywhere, which means H ice denotes the ice dissipation, H latent denotes the latent heat release and H cond denotes the conductive heat loss. The latent heat release H latent = ρ i L f q tends to be small compared to the other terms, as can be seen from Fig. 1c is the ice dissipation rate in a flat ice shell as a function of latitude φ. The conductive heat loss H cond is inversely proportional to the local ice thickness H cond (φ) = H cond0 × (H 0 /H(φ)). Since, the ice dissipation over the polar regions is roughly twice as strong as that over the equator in absence of ice topography (Beuthe 2018), I choose H ice0 (pole) = 1.25H ice0 for the polar box and H ice0 (eq) = 0.75H ice0 for the equatorial box, where (·) denotes global mean. To guarantee global heat budget balance, H ice0 should be equal to H cond0 ∆ = H. Now, consider the equator-to-pole difference of the heat budget terms for an ice shell that has a mean thickness of H 0 and an equator-to-pole thickness difference of ∆H (equatorial ice shell is thicker). Keeping the first order terms, the heat budget simplifies to (29) From Eqs. (29) and (27), ∆H can be solved numerically. The results are presented in Fig. 4. The entire solution falls into the κ v -limit regime, as denoted by the red shading (blue shading marks the D-limit regime).

----limit
κ v ---D-limit Figure 4. Predicted equator-to-pole ice thickness difference in equilibrium solved from Eqs. (29) and (27). Red shading marks the κv-limit regime and blue shading marks the D-limit regime. Highly saturated colors represent 3D scalings and lighter colors represent 2D scalings obtained in . When ∆H > H0/2 (H0 is the mean ice thickness), the polar ice shell thickness approaches zero, and the small ∆H assumption required for Eq. (30) no longer holds. Those scenarios are considered to be run-away poleward-thinning and mask those with gray shadings.
If the ice thickness variation is small (∆H/H 0 1), further simplification can be made to Eq. (29). Taylor expanding Eq. (29) around small ∆H/H 0 and dropping high order terms gives From Eqs. (30) and (27), ∆H can be solved analytically: ∆H ≈ χ κ κ −1/2 v α −3/10 f 2/5 a −7/10 , for κ v -limit constants. To obtain the above solution, I have substituted g with 4πGρ bulk a/3, where ρ bulk = 2500 kg/m 3 is the bulk density. From Eq. (31), it can be seen that, when ∆H H 0 , the equilibrium ice shell thickness variation ∆H should decrease with the moon's size and increases with the rotation frequency, following ∆H ∝ f 2/5 a −7/10 or f 8/13 a −7/13 depending on the dynamic regime of the ocean. The results here again qualitatively agree with the 2D scaling results, which follow ∆H ∝ (f /a) 2/3 or f /a 1/2 for κ v -limit and D-limit, respectively, as shown by , except that the sensitivity to rotation rate is slightly lower here.
The asymptotic scalings (shown by blue and red dashed lines in Fig. 4) provide a useful approximation to the full solution of Eq. (29) for relatively small to moderate ∆H. For larger ∆H, the sensitivity of ∆H on a increases, and eventually, a run-away poleward-thinning happens when a ≈ 200 km (masked by dark gray shading), due to the strengthening of the ice-rheology feedback. Besides the dependence on a and f , it can be seen from Eq. (31) that a higher ocean salinity (leading to larger α and stronger salinity-driven circulation), a stronger turbulent diffusivity, and (provided sufficient turbulent mixing) a deeper ocean can also reduce ∆H.

Numerical Results for Enceladus and Europa
To demonstrate the potential impacts of the size of the icy satellite on its equilibrium ice shell geometry, I integrate an ice evolution model forward using Enceladus' and Europa's parameters, respectively. The model is modified based on Kang & Flierl (2020). It calculates the melting induced by the tidal heating H ice (given by Eq.15 in the appendix), the down-gradient ice flow Q (given by Eq.13 in the appendix), the heat loss to space by conduction H cond (given by Eq.4 in the appendix), and heat transmitted upward by the ocean H ocn , and evolves the ice thickness H over time. The total thickness tendency can be symbolically expressed as where L f and ρ i are the latent heat of freezing and density for ice, a is the moon's radius and φ denotes latitude. When the ice shell is thinner than 3 km, I assume that the ice shell will crack open under the tidal stress, and the resultant geysers will carry away large amounts of heat, preventing further melting. In the evolution model, I overwrite the thickness tendency with zero when H < 3 km, to implicitly represent this additional heat loss. The ocean-ice heat exchange H ocn is a new component that doesn't exist in Kang & Flierl (2020). Inspired by the conceptual model, H ocn is parameterized as where H is the deviation from the prescribed global mean ice thickness H 0 = 20 km and MIN{} selects whichever parameterization yields a smaller global standard deviation. A factor of 2 is multiplied to H' because the analysis before all use equator to pole thickness difference, which is twice as large as the equator/pole's deviation from the mean. Also, salinity-driven heat transport is assumed to be roughly comparable with the temperature-driven one, and that accounts for the other factor of 2 multiplied to the above formula.
To account for the uncertainties associated with the ice shell rheology and the efficiency of ocean heat transport, a range of ice viscosities η m 5 , |α| and κ v are considered for Enceladus and Europa. The equilibrium ice shell geometries are shown in Fig. 5. For both moons, the equilibrium geometry also varies with the ice and ocean properties. The equilibrium ice shell tends to be flatter with smaller η m and higher |α| and κ v .
The bottom row assumes no ocean heat transport, and significant ice thickness variations develop on both Enceladus and Europa. When the ice viscosity is not too low (η m > 10 −13 Pa·s), ice is almost completely melt through over one or both poles due to the ice-rheology feedback. However, with ocean heat transport, the equilibrium ice geometry is largely flattened especially for Europa given its large size and slower rotation rate. The smoothing effect of ocean is particularly clearly shown in the right row, where ice flow is ineffective. If Europa's ocean is saltier than 50 psu as suggested by the strong magnetic induction signal (Hand & Chyba 2007), α should be closer to the upper bound, leading us to the conjecture that ice thickness variation on Europa is likely less than 1 km. When Enceladus' parameters are used instead, the ocean heat transport's impact on ice shell geometry is more limited. The equilibrium ice shell geometry obtained under the influence of OHT is similar to those without, unless the vertical diffusivity is very high. This sensitivity highlights the importance to understand the dissipative processes in the ocean driven by tides and libration motions (Rekier et al. 2019). It should also be noted that the OHT scaling obtained in this work may not apply to scenarios with very strong ice thickness variation due to the changes of geometry and the nonlinear behavior of eddies at high amplitude. Among the 28 scenarios considered for Enceladus, five develop the hemispheric asymmetry seen in observation (Iess et al. 2014;Hemingway & Mittal 2019), suggesting that the symmetry-breaking mechanism proposed by Kang & Flierl (2020) could work in presence of ocean heat redistribution.
The two scientific questions addressed here are 1) how the efficiency of the ocean heat transport (OHT) forced by the ice thickness variations varies with the icy moon's orbital parameters and 2) how the OHT in turn affects the equilibrium ice geometry. To do so, I derive scaling laws for the OHT on icy moons, inspired by previous theoretical work on the baroclinic eddies in the context of earth ocean or atmosphere (Held & Larichev 1996;Karsten et al. 2002;Jansen & Ferrari 2013). These scaling laws are then verified by 3D general circulation simulations run for various planetary radii, rotation rates and associated ice thickness variations. It is found that heat convergence toward the thick-ice regions is more efficient on icy satellites with greater sizes and slower rotation rates. Therefore, those icy moons' ice shells are expected to be flatter.
Enceladus and Europa are two icy satellites in the solar system that are known to contain a global subsurface ocean (Postberg et al. 2009;Thomas et al. 2016;Carr et al. 1998;Kivelson et al. 2000;Hand & Chyba 2007). Despite their similar global-mean ice thickness and per-area heat production rate, Europa's ice shell is likely to undergo a drastically different evolution path from that of Enceladus. Due to its larger size and slower rotation rate, the ocean heat transport on Europa is likely to be much more efficient and thus the equilibrium ice thickness variation is predicted to be lower than 1 km, in line with the so-far available observations (Iess et al. 2014;Beuthe et al. 2016;Tajeddine et al. 2017;Čadek et al. 2019;Hemingway & Mittal 2019;Nimmo et al. 2007).
To drive the point home, the equilibrium ice geometry for Enceladus and Europa are solved by integrating an ice evolution model, where OHT is parameterized based on the scaling laws. All Europa scenarios with OHT parameterization form a rather flat ice shell with thickness variation below 2km. Most equilibrium ice geometries obtained using Enceladus parameters, to the contrary, exhibit strong thickness variations, unless the ice sheet is very mobile (low ice viscosity) or the assumed vertical diffusivity and the thermal expansion coefficient are both high. Some Enceladus scenarios even form the significant hemispheric asymmetry seen in observations (Iess et al. 2014;Beuthe et al. Kang 2016;Tajeddine et al. 2017;Čadek et al. 2019;Hemingway & Mittal 2019), indicating that the symmetry breaking mechanism proposed by Kang & Flierl (2020) can work in presence of OHT.
It should be noted that other factors, such as ice viscosity, vertical mixing in the ocean and thermal expansion coefficient (determined by ocean salinity), also have significant impacts on the equilibrium ice geometry for Enceladus. These factors are thus far poorly constrained. More work along this line will improve the prediction of equilibrium ice shell geometry. Also, in this work, all heat is assumed to be produced in the ice shell. With heat produced in the silicate core, the ocean's stratification will change: a fresh ocean on a small moon with negative α will become more stratified, whereas a salty ocean on a large moon with positive α will become less stratified and even globally convective. In appearance of convection (salty/high pressure), the heating delivered to the ice shell by convective Taylor plumes is not going to be evenly distributed if the heating released from the seafloor is (Soderlund et al. 2014;Soderlund 2019;Bire et al. 2022). This may induce topography in the ice shell as well and needs to be investigated in the future. However, we expect the effect of bottom heating to play a less important role as the satellite size increases, because the vertical temperature gradient induced by bottom heating decrease with gravity 6 , whereas the temperature difference induced by ice topography will increase with satellite size; and even for an icy moon as small as Enceladus, the vertical temperature gradient induced by a 40 mW/m 2 bottom heating is likely one order of magnitude smaller than that induced by the observed ice thickness variation (Kang et al. 2021).
Despite these uncertainties, the qualitative result that ocean heat transport is more efficient at limiting ice-shell thickness variations on large satellites is likely to be robust. By connecting the equilibrium ice shell geometry with the icy moon's orbital parameters and the ocean properties, this work may bring a bit more constraints to the poorlyconstrained icy worlds.
This work is carried out in the Department of Earth, Atmospheric and Planetary Science (EAPS) in MIT. WK acknowledges support as a Lorenz-Houghton Fellow by endowed funds in EAPS and helpful comments from Prof. Malte Jansen. Our simulations are carried out using the Massachusetts Institute of Technology OGCM (MITgcm MITgcm-group 2010; Marshall et al. 1997) configured for application to icy moons.
The model integrates the non-hydrostatic primitive equations for an incompressible fluid in height coordinates, including a full treatment of the Coriolis force in a deep fluid, as described in MITgcm-group (2010); Marshall et al. (1997). Such terms are typically neglected when simulating Earth's ocean because the ratio between the fluid depth and horizontal scale is small. Instead, when the moon size is order hundreds of kilometers like Enceladus, the aspect ratio is order 0.1 and so not negligibly small. The size of each grid cell shrinks with depth due to spherical geometry and is accounted for by switching on the "deepAtmosphere" option of MITgcm. Also, the gravity will vary with depth as well. This is accounted for using the following profile of gravity.
In the above equation, G = 6.67 × 10 −11 N/m 2 /kg 2 is the gravitational constant, ρ core = 2500 kg/m 3 is the assumed core density and ρ out = 1000 kg/m 3 is the density of the ocean/ice layer. D 0 and H 0 is the thickness of ocean and ice on global average. Since it takes several tens of thousands of years for our solutions to reach equilibrium, all of our experiments are first run under a zonally symmetric 2D configuration with a moderate resolution of 2 degree (8.7 km). Only 30 layers (each 2 km) are used to keep the computational cost manageable. After equilibrium is reached, I interpolate the pick up files to generate initial conditions for the corresponding 3D simulation, which has a default horizontal resolution of 0.25×0.25 degree and 70 unevenlly distributed vertical layers, whose thicknesses increase from 500 m to 2 km from top to bottom. Since changing rotation rate has significant impacts on the size of the jets and eddies, I had to adjust the grid width along east-west direction by a factor of 1.4 in those cases to better capture the dynamics.

A.1. Diffusivity and Viscosity.
Vertical diffusivity affects the energetics of the ocean (Young 2010;Jansen et al. 2022). To account for the mixing of heat and salinity by unresolved turbulence, in our calculations, I set the explicit vertical diffusivity to 0.001 m 2 /s in both 2D and 3D simulations, following Kang et al. (2021). This is roughly 4 orders of magnitude greater than molecular diffusivity, but broadly consistent with dissipation rates suggested by Rekier et al. (2019) for Enceladus, according to the scaling suggested by Wunsch & Ferrari (2004). In all experiments, horizontal diffusivity is set to be equal to the vertical diffusivity regardless of the resolution and the size of the icy moon.
Viscosity is necessary to keep the model stable. For the coarse resolution 2D simulations, the horizontal viscosity are set to a 150 km m 2 /s (a is the radius of the moon). Additionally, a bi-harmonic hyperviscosity of a 150 km × 10 8 m 4 /s is employed to further damp numerical noise induced by our use of stair-like ice topography. The same viscosities are used in . For the high resolution 3D simulations, the explicit horizontal and vertical viscosity is set to much smaller values (0.08 m 2 /s and 0.03 m 2 /s respectively), and I use the widely-applied Smagorinsky viscosity scheme (Smagorinsky 1963). Unlike the fixed viscosity scheme, Smagorinsky scheme determines the viscosity based on the resolved dynamics, and as a result, numerical noise will be damped while dynamics can be kept to a larger extent. The Smagorinsky viscosity constant is set to 3 by default. As mentioned before, resolution in the x-direction is increased (decreased) by a factor of 1.4 under a higher (lower) rotation rate. In those experiments, horizontal viscosity is adjusted proportional to the x-grid width.
In the coarse resolution model, convection cannot be resolved, so parameterization is needed. Following Kang et al. (2021), I set the diffusivity to a much larger value in convectively unstable regions, to represent the vertical mixing associated with convective overturns. This convective diffusivity κ conv is set to increase from 1 m 2 /s to 30 m 2 /s as gravity and convection strengthens with the satellite radius. Similar approaches are widely used to parameterize convection in coarse resolution ocean models (see, e.g. Klinger et al. (1996)) and belong to a family of convective adjustment schemes. Our results turn out to be insensitive to κ conv , as long as the convective timescale D 2 /κ conv < 1 yr is much shorter than the advective time scale M half /Ψ ≈ 10 2 -10 3 yrs (M half is half of the total mass of the ocean and Ψ is the maximum meridional streamfunction in kg/s).

A.2. Equation of state and the freezing point of water
To make the dynamics as realistic as possible, the "MDJWF" equation of state (EOS McDougall et al. 2003) is adopted when it comes to determine density using temperature, salinity and pressure. As demonstrated by Fig. 1e of the main text, the thermal expansion coefficient α at the freezing point is negative at the ice-ocean interface when the moon size is small (low pressure) and when the ocean is fresh. This anomalous expansion can suppress the convection driven by bottom heating (Zeng & Jansen 2021;Kang et al. 2022) and can alter the direction of ocean circulation (Kang et al. 2021).
The freezing point of water T f is assumed to depend on local pressure P and salinity S as follows, where a 0 = −0.0575 K/psu, b 0 = −7.61 × 10 −4 K/dbar and c 0 = 0.0901 degC. The pressure P can be calculated using hydrostatic balance P = ρ i gH (ρ i = 917 kg/m 3 is the density of the ice and H is the ice thickness).

A.3. Boundary conditions
The ocean is encased by an ice shell with meridionally-varying thickness, assuming hydrostacy (i.e., ice is floating freely on the water). The ice thickness is set to be where H 0 is the mean ice thickness, P 2 is the 2nd order Legendre polynomial, and H 2 is the amplitude of the ice thickness variation. φ denotes latitude. The thickness profile is shown by a solid curve in Fig.1b of the main text.
Partial cells is switched on to better represent the ice topography: water is allowed to occupy a fraction of the height of a whole cell with an increment of 10%. Interactions between the ice shell and the ocean is taken care of by a modified version of the MITgcm's "shelfice" module (Losch 2008). The ocean is forced by heat and salinity fluxes from the ice shell at the top. Diffusion of heat through the ice Heat loss to space by heat conduction through the ice H cond is represented using a 1D vertical heat conduction model, where H is the thickness of ice (solid curve in Fig.1b of the main text), the surface temperature is T s and the ice temperature at the water-ice interface is the local freezing point T f (Eq. A2). The surface temperature T s is set to the radiative equilibrium temperature, which can be computed given the incoming solar radiation and obliquity (δ = 3 • ) and assuming an albedo of 0.81. Typical heat losses averaged over the globe are H cond = 50 mW/m 2 , broadly consistent with observations (Tajeddine et al. 2017).
Ice-ocean fluxes The interaction between ocean and ice is simulated using MITgcm's "shelf-ice" package (Losch 2008;Holland & Jenkins 1999) with some modifications.
At the water-ice interface, we consider the response of the ocean to a prescribed ice freezing rate while ignoring the possible response of the ice to the water-ice heat/salinity exchange. The freezing/melting induces a salinity/fresh water flux into the ocean (we assume the ice salinity to be zero); meanwhile, the ocean temperature at the upper boundary is relaxed to the local freezing point T f determined by the local salinity and pressure (Eq. A2).
Here, S ocn−top and T ocn−top denote the upper boundary salinity and temperature, γ T = 10 −5 m/s are the water-ice exchange coefficients for temperature and salinity, δz = 2 km is the thickness of the water-ice "boundary layer" and q is the freezing rate in m/s (note that q is orders of magnitude smaller than γ T ). The "boundary layer" option is switched on to avoid possible numerical instabilities induced by an ocean layer which is too thin. In addition to the above conditions on temperature and salinity, the tangential velocity is relaxed back to zero at a rate of γ M = 10 −3 m/s at the upper and lower boundaries.

A.4. Ice flow model
The prescribed freezing rate q is computed using the divergence of the ice flow, assuming the ice sheet geometry is in equilibrium. Here, an upside-down land ice sheet model is used following Ashkenazy et al. (2018). The ice flows down its thickness gradient, driven by the pressure gradient induced by the spatial variation of the ice top surface, somewhat like a second order diffusive process. At the top, the speed of the ice flow is negligible because the upper part of the shell is so cold and hence rigid; at the bottom, the vertical shear of the ice flow speed vanishes, as required by the assumption of zero tangential stress there. This is the opposite to that assumed in the land ice sheet model. In rough outline, I calculate the ice flow using the expression below obtained through repeated vertical integration of the force balance equation (the primary balance is between the vertical flow shear and the pressure gradient force), using the aforementioned boundary conditions to arrive at the following formula for ice transport Q, where Here, φ denotes latitude, a and g are the radius and surface gravity of the moon, T s and T f are the temperature at the ice surface and the water-ice interface (equal to local freezing point, Eq. A2), and ρ i = 917 kg/m 3 and ρ 0 are the ice density and the reference water density. E a = 59.4 kJ/mol is the activation energy for diffusion creep, R g = 8.31 J/K/mol is the gas constant and η m is the ice viscosity at the freezing point. The latter has considerable uncertainty (10 13 -10 16 Pa·s Tobie et al. 2003), and here η m is set to 10 14 Pa·s. In steady state, the freezing rate q must equal the divergence of the ice transport thus: As shown by the dashed curve in Fig.1b of the main text, ice melts in high latitudes and forms in low latitudes at a rate of a few kilometers every million years. A more detailed description of the ice flow model can be found in Kang & Flierl (2020) and Ashkenazy et al. (2018). Freezing and melting leads to changes in local salinity and thereby a buoyancy flux.
A.5. Model of tidal dissipation in the ice shell Icy moon's ice shell is periodically deformed by tidal forcing and the resulting strains in the ice sheet produce heat. I follow Beuthe (2019) to calculate the ice dissipation rate. Instead of repeating the whole derivation here, I only briefly summarize the procedure and present the final result. Unless otherwise stated, parameters are the same as assumed in Kang & Flierl (2020).
Tidal dissipation consists of three components (Beuthe 2019): a membrane mode H mem ice due to the extension/compression and tangential shearing of the ice membrane, a mixed mode H mix ice due to vertical shifting, and a bending mode H bend ice induced by the vertical variation of compression/stretching. Following Beuthe (2019), I first assume the ice sheet to be completely flat. By solving the force balance equation, I obtain the auxiliary stress function F , which represents the horizontal displacements, and the vertical displacement w. The dissipation rate H flat,x ice (where x = {mem, mix, bend} ) can then be written as a quadratic form of F and w. In the calculation, the ice properties are derived assuming a globally-uniform surface temperature of 60K and a melting viscosity of 5 × 10 13 Pa·s.
Ice thickness variations are accounted for by multiplying the membrane mode dissipation H flat,mem ice , by a factor that depends on ice thickness. The membrane mode is the only mode which is amplified in thin ice regions (see Beuthe (2019)). This results in the expression: where H is the prescribed thickness of the ice shell as a function of latitude and H 0 is the global mean of H. Since thin ice regions deform more easily and produce more heat, p α is negative. Because more heat is produced in the ice shell, the overall ice temperature rises, which, in turn, further increases the mobility of the ice and leads to more heat production (the rheology feedback). The tidal heating profile corresponding to p α = −1.5 is the red solid curve plotted in Fig.1c of the main text.
B. IDEALIZED ICE EVOLUTION MODEL.
Here I provide a brief overview for the idealized model used to evolve the ice shell of Enceladus and Europa. Interested readers are referred to (Kang & Flierl 2020) and its supplementary material for more detail.
In this model, ice shell thickness H is changed over time by the melting induced by the tidal heating H ice (given by Eq. A9), the down-gradient ice flow Q (given by Eq. A7), the heat loss to space by conduction H cond (given by Eq. A4), the crack-induced cooling H crack in places where ice is sufficiently thin, and heat transmitted upward by the ocean H ocn . The ice thickness tendency can be symbolically expressed as following, where L f and ρ i are the latent heat of freezing and density for ice, a is the moon's radius and φ denotes latitude. Physical constants and parameters for Enceladus and Europa are stated in Table.2. H ice is polar-amplified, and as a result, the polar ice shell tends to be thinner, which in turn increases the heat production over the pole (see Eq. A9). The tendency for ice thickness variation to increase due to the rheology feedback will be balanced by the rapid heat loss through thin ice (Eq. A4) and the transport by ice flow (Eq. A7). An additional heat sink is activated only when the ice thickness is less than H crack = 3 km to prevent further melting, and crudely represents the effect of cracks and geysers that carry the extra heat away. For all time, the global tidal dissipation H ice is scaled to exactly balance the Kang instantaneous conductive heat loss H cond . By so doing, the rheology feedback and thus the ice thickness variation are maximized. Throughout the integration, the global mean ice thickness is fixed at H 0 = 20 km. The ocean-ice heat exchange is a new process I introduced. Inspired by the conceptual model, the heat flux coming from the ocean is parameterized by Eq. 19 in the main text.
The initial condition is set as follows where H 0 = 20 km, H 2 = 3 km and H 1 = 1 km. P 1 and P 2 are the first and second order of Legendre Polynomials. initial equator-to-pole ice thickness variation 3 km H1 initial hemispherical asymmetry -1 km Ts mean surface temperature 110K    Figure 9. Same as Fig. 3 except rotation period is set to 10 day instead of 3.5 day.  From blueish color to reddish color denotes increasing ocean salinity.