Different Ice-shell Geometries on Europa and Enceladus due to Their Different Sizes: Impacts of Ocean Heat Transport

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 (Carr et al. 1998; Kivelson et al. 2000; Hand & Chyba 2007; Postberg et al. 2009; Thomas et al. 2016). 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, measurements and detection have been 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 so that we can more easily take measurements. As illustrated by Kang et al. (2021) and Kang & Jansen (2022), 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 Figure 1(d)). 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 (Kang & Jansen 2022).

Zoom In Zoom Out Reset image size

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). While 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 (Hurford et al. 2007; Hedman et al. 2013; Nimmo et al. 2014; Ingersoll et al. 2020) and ejecting samples of the ocean to outer space (Hansen et al. 2006; Porco et al. 2006; Howett et al. 2011; Spencer 2013). In order to sustain the strong ice topography, the ocean heat transport (OHT) that flattens the ice shell through the 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 of 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, Kang & Jansen (2022) 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 Kang & Jansen (2022) are eddy dynamics. We know from Earth's ocean that oceans 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 to Earth's 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 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: The curl of the wind stress forces meridional overturning flows, maintaining an isopycnal slope so that eddies can grow on it (McWilliams et al. 1978; Marshall & Radko 2003; Plumb & Ferrari 2005). On the contrary, the ocean on an icy satellite is sandwiched between an ice shell and a silicate core and thus experiences no wind stress. Density gradients created by the ocean–ice heat and salinity exchange, instead, drive the overturning circulation there. Because 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 (Kang & Jansen 2022). 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.

The system considered by this study is sketched in Figure 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):

$$\begin{eqnarray}&&H(\phi )={H}_{0}+{\rm{\Delta }}{{HP}}_{2}(\sin \phi ),\end{eqnarray} \tag{ 1 }$$

where H₀ = 20 km is the mean ice thickness, P₂ is the second-order Legendre polynomials, and ϕ denotes latitude. Here, the ice-thickness variation is assumed to follow the P₂ 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 Figure 1(b). One would expect the ice shell to be thinner over the poles as can be seen in Figure 1(b) because the tidal dissipation produced in the ice is stronger there (Beuthe 2018, 2019; Kang & 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₂, as long as the ice thickness is relatively simple (poleward thinning or poleward thickening) and is symmetric about the equator.

With ice flowing from thick-ice regions to thin-ice regions, ¹ this equator-to-pole ice-thickness gradient will not 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 ${{ \mathcal H }}_{\mathrm{cond}}$ is around 40 mW m⁻² on global average and is faster over regions where the ice is thin and over the poles where the ice surface temperature is low, as shown by the green curve in Figure 1(c). To balance the heat loss, the ice shell and the silicate core need to produce heat; however, the core–shell heat partition is poorly understood (Choblet et al. 2017; Beuthe 2019). 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 ${{ \mathcal H }}_{\mathrm{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 Figure 1(c) is ${{ \mathcal H }}_{\mathrm{ice}}$ with the default ice-thickness profile (Equation (1)). Besides, the equatorial freezing and polar melting required to maintain the prescribed ice thickness give rise to almost negligible latent heat release ${{ \mathcal H }}_{\mathrm{latent}}$, shown by a gray curve in Figure 1(c). At the bottom of the ice, there may also be heat exchange with the ocean (denoted by ${{ \mathcal H }}_{\mathrm{ocn}}$), manifested by ocean circulation and heat transport. In an equilibrium state, the ice-shell heat budget needs to be in balance, which means the sum of all the heating terms should vanish.

From the 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 Equation (A2) in Appendix A). Also, assuming the ice shell is in mass equilibrium, equatorial freezing and polar melting are required to prevent the ice shell from being flattened by the pressure-driven ice flow (see the dashed gray curve in Figure 1(b)). The freezing/melting will then induce salinity exchange with the subsurface ocean. Under these forcings, water over the poles becomes 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 turn affects 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 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 and test them against numerical simulations. Because we already know how the other heating terms (conductive heat loss, ice dissipation, and latent heating) depend 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 are discussed in the conclusion.

More details about the ocean circulation model, ice-flow model, and tidal dissipation model can be found in Appendix A.

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, temperature-induced density variation dominates, and the warm polar water would sink as sketched by the blue arrow in Figure 1(d) because freshwater contracts upon warming (anomalous expansion), while 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 Figure 1(d).

When considering icy satellites larger than Enceladus (most of the icy satellites of interest are), the following changes are expected:

• 1.  
  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 Figure 1(e) is 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 does not 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, ² approximately the size of Europa. Therefore, with a large-enough planetary size, downwelling will only occur over the equator regardless of the ocean salinity.
• 2.  
  The temperature difference under the ice shell between the equator and the pole ΔT will increase, as shown by the crosses in Figure 1(f). 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,

  $$\begin{eqnarray}&&{\rm{\Delta }}T={\rm{\Delta }}{T}_{f}={b}_{0}{\rm{\Delta }}P={b}_{0}{\rho }_{i}{g}_{0}{\rm{\Delta }}H,\end{eqnarray} \tag{ 2 }$$

  where ΔH is the prescribed ice-thickness difference between the equator and the poles, g₀ is the moon's surface gravity, b₀ = −7.61 × 10⁻⁴ K/dbar, and ρ_i is the ice density. With ΔH fixed, ΔT ∝ g ∝ a.
• 3.  
  Salinity forcing will weaken. The blue and brown dots in Figure 1(f) show the salinity flux (mean salinity S₀ times the freezing rate q) difference between the equator and the poles for an ocean with a 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 Appendix A, ice flow behaves like diffusion, and the flow divergence/convergence is proportional to ΔP divided by the distance square a², so

  $$\begin{eqnarray}&&q\propto {\rm{\Delta }}P/{a}^{2}\propto {a}^{-1}.\end{eqnarray} \tag{ 3 }$$
• 4.  
  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₀ ∝ 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 the melting-point ice viscosity is 10¹⁴ 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 than or comparable to that of Enceladus. In Kang & Jansen (2022), 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,

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{ocn}}\propto {a}^{3}{\rm{\Delta }}{H}^{3/2}{f}^{-1},\end{eqnarray} \tag{ 4 }$$

when the circulation depth is limited by vertical diffusion, or

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{ocn}}\propto {a}^{3}{\rm{\Delta }}{H}^{2}{f}^{-2},\end{eqnarray} \tag{ 5 }$$

when the circulation reaches the seafloor.

Given that salinity forcing tends to be dominated by temperature forcing unless the size of the icy moon is comparable to or smaller than Enceladus (see the previous section), only the temperature-induced density anomalies are considered here, following Kang & Jansen (2022). 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 ${{ \mathcal F }}_{e}$ can be represented by an equivalent diffusive process,

$$\begin{eqnarray}&&{{ \mathcal F }}_{e}={C}_{p}{\rho }_{0}{\kappa }_{e}{\partial }_{y}T\cdot (2\pi {ad})\sim 2\pi {C}_{p}{\rho }_{0}{\kappa }_{e}{\rm{\Delta }}{Td},\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&{{\rm{\Psi }}}^{\dagger }\sim (2\pi a){\rho }_{0}{\kappa }_{e}\displaystyle \frac{{T}_{y}}{{T}_{z}}\sim 2\pi {\rho }_{0}{\kappa }_{e}d({\rm{\Delta }}T{/{\rm{\Delta }}}_{v}T),\end{eqnarray} \tag{ 7 }$$

where Δ_v T is the vertical temperature contrast across the depth d.

According to Held & Larichev (1996), κ_e can be estimated by

$$\begin{eqnarray}&&{\kappa }_{e}={{kVL}}_{\beta },\end{eqnarray} \tag{ 8 }$$

where k = 0.25 is a constant, V is the rms eddy velocity, and L_β denotes the wavelength of the energy-containing eddies, which follows the Rhines scale:

$$\begin{eqnarray}&&{L}_{\beta }\sim \sqrt{V/\beta }.\end{eqnarray} \tag{ 9 }$$

This energy-containing wavelength is typically larger than the deformation radius L_d , the scale at which baroclinic instability happens:

$$\begin{eqnarray}&&{L}_{d}\sim \displaystyle \frac{{Nd}}{f}\sim \displaystyle \frac{\sqrt{\alpha g{{\rm{\Delta }}}_{v}{Td}}}{f},\end{eqnarray} \tag{ 10 }$$

f denotes the Coriolis coefficient, and $N=\sqrt{g{\rho }_{z}/{\rho }_{0}}$ denotes the Brunt–Väisälä frequency. Here, the vertical density gradient ρ_z is estimated by ρ₀ αΔ_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 ξ,

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\beta }}{{L}_{d}}\sim \displaystyle \frac{V}{U}\sim \xi \equiv \displaystyle \frac{{\rho }_{y}f/\beta }{{\rho }_{z}d}\sim \displaystyle \frac{{\rm{\Delta }}T}{{{\rm{\Delta }}}_{v}T},\end{eqnarray} \tag{ 11 }$$

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,

$$\begin{eqnarray}&&U\sim \displaystyle \frac{\alpha g{\rm{\Delta }}{Td}}{{fa}}.\end{eqnarray} \tag{ 12 }$$

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²) should be equal to the downward heat transport by residual circulation C_p Ψ^†ΔT (Jansen & Ferrari 2013) and that yields

$$\begin{eqnarray}&&{\kappa }_{v}{a}^{2}={\kappa }_{e}{d}^{2}{\xi }^{2}.\end{eqnarray} \tag{ 13 }$$

κ_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, and L_d and L_β are interchangeable ³ , so κ_e can be written as

$$\begin{eqnarray}&&{\kappa }_{e}={{kUL}}_{d}=\displaystyle \frac{k\beta }{{f}^{3}}{\left(\alpha g{\rm{\Delta }}{Td}\right)}^{3/2}\sim \displaystyle \frac{k}{{{af}}^{2}}{\left(\alpha g{\rm{\Delta }}{Td}\right)}^{3/2}.\end{eqnarray} \tag{ 14 }$$

The depth of density variations d is constrained by the vertical diffusivity κ_v through Equation (13), leading to

$$\begin{eqnarray}&&d\sim {\left(\displaystyle \frac{{\kappa }_{v}{a}^{3}{f}^{2}}{k}\right)}^{2/7}\cdot {\left(\alpha g{\rm{\Delta }}T\right)}^{-3/7}\propto {f}^{4/7}{\rm{\Delta }}{H}^{-3/7},\end{eqnarray} \tag{ 15 }$$

which indicates that the depth of the 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 to be smaller than D, then we can substitute Equation (15) and Equation (14) into Equation (6), which gives

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal F }}_{{\kappa }_{v}}\sim 2\pi \rho {C}_{p}{\rm{\Delta }}T{\left(\displaystyle \frac{k}{{{af}}^{2}}\right)}^{2/7}\\ \times \ {\left(\alpha g{\rm{\Delta }}T\right)}^{3/7}{\left({\kappa }_{v}{a}^{2}\right)}^{5/7}\propto {a}^{3}{\rm{\Delta }}{H}^{10/7}{f}^{-4/7}.\end{array}\end{eqnarray} \tag{ 16 }$$

To obtain the proportionality, the fact that ΔT ∝ ΔP ∝gΔH ∝ aΔH, g ∝ a, and β ∝ f/a is used. The dependence of ${{ \mathcal F }}_{{\kappa }_{v}}$ on a and ΔH are very similar to that of the κ_v -limited overturning circulation in the zonally symmetric case given by Kang & Jansen (2022) (see Equation (4)).

Scaling laws can also be obtained for L_d and L_β , which should characterize the size of eddies and the width of jets,

$$\begin{eqnarray}&&{L}_{\beta ,{\kappa }_{v}}={L}_{d,{\kappa }_{v}}={\left(\alpha g{\rm{\Delta }}T\right)}^{2/7}{\left({\kappa }_{v}{a}^{3}/k\right)}^{1/7}{f}^{-5/7}\propto a{\rm{\Delta }}{H}^{2/7}{f}^{-5/7},\end{eqnarray} \tag{ 17 }$$

and for the jet speed U

$$\begin{eqnarray}&&\begin{array}{l}{U}_{{\kappa }_{v}}=\beta {L}_{\beta ,{\kappa }_{v}}^{2}\sim {\left(\alpha g{\rm{\Delta }}T\right)}^{4/7}\\ \times \ {\left({\kappa }_{v}/k\right)}^{2/7}{f}^{-3/7}{a}^{-1/7}\propto a{\rm{\Delta }}{H}^{4/7}{f}^{-3/7}.\end{array}\end{eqnarray} \tag{ 18 }$$

In the above equations, the subscript κ_v indicates that this is the κ_v -limit scaling. Noticeably, the above scaling laws predict the eddy size and jet width to grow linearly with the 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

$$\begin{eqnarray}&&{\kappa }_{e}={{kVL}}_{\beta }={kU}\sqrt{U/\beta }{\xi }^{3/2}=\displaystyle \frac{k\beta }{{f}^{3}}{\left(\alpha g{\rm{\Delta }}{TD}\right)}^{3/2}{\xi }^{3/2}.\end{eqnarray} \tag{ 19 }$$

From the above equation and the vertical heat balance (Equation (13)), ξ can be solved:

$$\begin{eqnarray}&&\xi ={\left({\kappa }_{v}{a}^{2}\displaystyle \frac{{f}^{3}}{k\beta }\right)}^{2/7}{\left(\alpha g{\rm{\Delta }}T\right)}^{-3/7}{D}^{-1}.\end{eqnarray} \tag{ 20 }$$

Substituting Equation (19) and Equation (20) into Equation (6) yields the eddy heat transport:

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal F }}_{D}\sim 2\pi \rho {C}_{p}D{\rm{\Delta }}T{\left(\displaystyle \frac{k\beta }{{f}^{3}}\right)}^{4/7}\\ \times \ {\left(\alpha g{\rm{\Delta }}T\right)}^{6/7}{\left({\kappa }_{v}{a}^{2}\right)}^{3/7}\propto {a}^{3}{\rm{\Delta }}{H}^{13/7}{f}^{-8/7}.\end{array}\end{eqnarray} \tag{ 21 }$$

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 to the bottom and the ice shell. Again, this scaling is comparable to the 2D scaling obtained by Kang & Jansen (2022) (see Equation (5)).

By assumption, the penetration depth of the T/S anomalies is the entire ocean depth D, i.e.,

$$\begin{eqnarray}&&d=D.\end{eqnarray} \tag{ 22 }$$

Substituting Equation (22) and Equation (20) into the definitions of the deformation radius and Rhines scale (Equation (9), Equation (11), and Equation (10)), the following scalings are obtained:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{d,D} & = & {\left(\alpha g{\rm{\Delta }}T\right)}^{5/7}{\left({\kappa }_{v}{a}^{2}/(k\beta )\right)}^{-1/7}{{Df}}^{-10/7}\\ & \propto & {a}^{9/7}{\rm{\Delta }}{H}^{5/7}{f}^{-11/7},\end{array}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{L}_{\beta ,D}={\left(\alpha g{\rm{\Delta }}T\right)}^{2/7}{\left({\kappa }_{v}{a}^{2}/(k\beta )\right)}^{1/7}{f}^{-4/7}\propto a{\rm{\Delta }}{H}^{2/7}{f}^{-5/7},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{U}_{D}=(\alpha g{\rm{\Delta }}T){{Df}}^{-2}\beta \propto a{\rm{\Delta }}{{Hf}}^{-1}.\end{eqnarray} \tag{ 25 }$$

Just as in the κ_v -limit scaling, the dominant eddy wavenumber and jet number should be insensitive to the radius of the moon because the power of the a factors for L_d,D and L_β,D are both equal to 1. Also, slower rotation and a 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 Equation (16) and Equation (21), i.e.,

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{ocn}}\sim \min \left\{{{ \mathcal F }}_{{\kappa }_{v}},{{ \mathcal F }}_{D}\right\},\end{eqnarray} \tag{ 26 }$$

because both ocean depth and vertical diffusivity constrain how much heat can be transported by the ocean.

Shown in Figures 2(a)–(c) is the predicted ${{ \mathcal F }}_{\mathrm{ocn}}$ as a function of the icy moon's size, rotation rate, and equator-to-pole thickness difference in highly saturated colors. Bluish 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 Kang & Jansen (2022) 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 Figure 2 using dots and diamond markers, respectively. They match the prediction (Equation (26)) up to a factor of 2. The matching seems particularly good in panel (a), because the range of ${{ \mathcal F }}_{\mathrm{ocn}}$ there is much wider so that the offset appears smaller in comparison.

Zoom In Zoom Out Reset image size

Table 1. Model Parameters Used in the Ocean General Circulation Model and the Conceptual Model

Symbol	Name	Definition/Value
Physical Constants
Lf	Fusion energy of ice	334,000 J kg−1
Cp	Heat capacity of water	4000 J kg−1 K−1
Tf (S, P)	Freezing point	Equation (A2)
ρi	Density of ice	917 kg m−3
ρw	Density of the ocean	"MDJWF" scheme McDougall et al. (2003)
α	Thermal expansion coeff.	− ∂(ρ/ρ0)/∂T
β	Saline contraction coeff.	∂(ρ/ρ0)/∂S
κ0	Conductivity coeff. of ice	651 W m−1
pα	Ice-dissipation amplification factor	−1.5
ηm	Ice viscosity at freezing point	1014 Ps·s
Default Model Setup
a	Radius	2500 km
g0	Surface gravity	Equation (A1)
δ	Obliquity	31
H0	Global-mean ice thickness	20 km
H2	Equator-to-pole ice-thickness variation	3 km
H	Ice-shell thickness	Equation (A3)
D	Global-mean ocean depth	56 km
Ω	Rotation rate	2.05 × 10−5 s−1
$\bar{{T}_{s}}$	Mean surface temperature	110 K
S0	Mean ocean salinity	60 psu
P0	Reference pressure	ρi g0 H0
T0	Reference temperature	Tf (S0, P0)
νh	Horizontal viscosity	0.08 m2 s−1 (3D), (a/150 km) m2 s−1 (2D)
νv	Vertical viscosity	0.03 m2 s−1 (3D), 1 m2 s−1 (3D)
νsmag	Smagorinsky viscosity (3D only)	3
${\tilde{\nu }}_{h},{\tilde{\nu }}_{v}$	Biharmonic hyperviscosity (2D only)	$\tfrac{a}{150\,\mathrm{km}}\times$ 108 m4 s−1
κh , κv	Horizontal/vertical diffusivity	0.001 m2 s−1
(γT , γS , γM )	Water–ice exchange coeff. for T, S, & momentum	(10−5, 10−5, 10−4) m s−1
${{ \mathcal H }}_{\mathrm{cond}}$	Conductive heat loss through ice	Equation (A4)
${{ \mathcal H }}_{\mathrm{ice}}$	Tidal heating produced in the ice	Equation (A9)

Download table as:  ASCIITypeset image

The OHT is governed by the ocean's dynamic and thermodynamic states. Presented in Figure 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 days) are 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 Figure 3(a). Meanwhile, in order to sustain the ice geometry against the flattening due to ice flow (Equation (A6)), equatorial ice needs to freeze and polar ice shell needs to melt. This in turn drives the meridional salinity gradient seen in Figure 3(b). High ocean salinity and high pressure suppress the anomalous expansion behavior (contract upon warming) near the freezing point, which typically happens on a small icy moon with a fresh ocean. As a result, temperature and salinity anomalies both contribute to the high density in low latitudes, driving sinking motions there (see Figure 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 to 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 where friction can prevent radial flows from inducing strong zonal jets that can by no means be in thermal-wind balance with the weak density gradients. In an equilibrium state, the upper part of the ocean flows westward and the lower part of the ocean flows eastward (Figure 3(d)), in thermal-wind balance with the density distribution (Figure 3(c)).

Zoom In Zoom Out Reset image size

The major differences between the 3D and 2D configurations lie in the zonal flow field (Figure 3(d)): There are jets formed in 3D due to the Reynold stress associated with the baroclinic eddies, whose structure is presented in Figures 3(f) and (g). Due to the strong rotation effect, the eddy motions tend to be aligned with the direction of the rotation axis, forming convective "rolls" along the equatorial plane and wave-like structures in higher latitudes, consistent with recent icy-moon studies by Ashkenazy & Tziperman (2021), Soderlund et al. (2014), Bire et al. (2022), and Kang et al. (2020). These eddies facilitate stronger equatorward heat transport than the overturning circulation in 2D configuration, as shown by Figure 3(h). It should be noted that some previous works have jets and Taylor columns even under 2D configuration (Ashkenazy & Tziperman 2021). This difference arises from a different model configuration: In their work, the ocean is heated from below and the 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 moon's sizes, rotation rates, and ice-thickness variations, the dynamics remain qualitatively the same, as shown by Figures A1–A6 in Appendix A. On a smaller icy moon, the meridional temperature gradient weakens due to its weaker gravity (compare Figure A1(a) and Figure A2(a) with Figure 3(a)). That in turn weakens the circulation (see panel (e1)) and eddies (see panels (f) and (g)), leading to a much weaker heat transport (see panel (h)) as suggested by Equation (16) and Equation (21). According to Equation (17) and Equation (24), the number of jets remains more or less unchanged—this can be seen from the panel (d2) of Figure A1, Figure A2, and Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Fast rotation suppresses the circulation, eddies, and thereby OHT, as suggested by both 2D scalings (Equation (4) and Equation (5)) and 3D scalings (Equation (16) and Equation (21)). This trend is reflected by sensitivity tests shown in Figure A3 and Figure A4, 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 the rotation rate varies, the width of jets and eddies should change accordingly. Strong rotation leads to smaller deformation radius L_d = (NH/f) (governing the eddy size, Equation (10)) and smaller Rhines scale ${L}_{\beta }=\sqrt{U/\beta }$ (governing the jet width, Equation (9)), consistent with panels (d), (f), and (g) of Figure A3, Figure 3, and Figure A4. The dependence of 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 in Figure A5 and Figure A6.

As 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 Equation (18) and Equation (25). The averaged spacing between jets is identified by "numpy.find_peaks" and is compared against Equation (17) and Equation (24). Finally, the prominent eddy size is computed by averaging different wavelengths by the corresponding power spectrum of the U field between 40N/S and 60N/S, and the results are compared against Equation (17) and Equation (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 (Kang et al. 2020; Bire et al. 2022). As shown in Figure A7, the theory also captures the eddy and jet characteristics reasonably well, indicating that the agreement between the predicted and simulated ${{ \mathcal F }}_{\mathrm{ocn}}$ may not be a coincidence.

6.1. Scaling Theory

The dependence of ${{ \mathcal F }}_{\mathrm{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 Kang & Jansen (2022) but using the 2D scaling for ${{ \mathcal F }}_{\mathrm{ocn}}$. Here, I repeat the process for 3D scalings.

First, I need to convert ${{ \mathcal F }}_{\mathrm{ocn}}$ to the 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², the heat flux per area received by the equatorial ice shell and the heat flux leaving the polar ice shell equal to

$$\begin{eqnarray}{{ \mathcal H }}_{\mathrm{ocn}}={{ \mathcal F }}_{\mathrm{ocn}}/(\pi {a}^{{2})^\sim} \left\{\begin{array}{ll}2\rho {C}_{p}{\rm{\Delta }}T{\left(\frac{k}{{a}^{3}{f}^{2}}\right)}^{2/7}{\left(\alpha g{\rm{\Delta }}T\right)}^{3/7}{\kappa }_{v}^{5/7}\propto a{\rm{\Delta }}{H}^{10/7}{f}^{-4/7} & \mathrm{for}\,\mathrm{the}\ {\kappa }_{v}\,\mathrm{limit}\\ 2\rho {C}_{p}D{\rm{\Delta }}T{\left(\frac{k\beta }{{a}^{2}{f}^{3}}\right)}^{4/7}{\left(\alpha g{\rm{\Delta }}T\right)}^{6/7}{\kappa }_{v}^{3/7}\propto a{\rm{\Delta }}{H}^{13/7}{f}^{-8/7} & \mathrm{for}\,\mathrm{the}\ D\,\mathrm{limit}\end{array}\right..\end{eqnarray} \tag{ 27 }$$

In an equilibrium state, the ice shell needs to be in heat balance everywhere, which means

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{ice}}+{{ \mathcal H }}_{\mathrm{latent}}+{{ \mathcal H }}_{\mathrm{ocn}}={{ \mathcal H }}_{\mathrm{cond}}.\end{eqnarray} \tag{ 28 }$$

${{ \mathcal H }}_{\mathrm{ice}}$ denotes the ice dissipation, ${{ \mathcal H }}_{\mathrm{latent}}$ denotes the latent heat release, and ${{ \mathcal H }}_{\mathrm{cond}}$ denotes the conductive heat loss. The latent heat release ${{ \mathcal H }}_{\mathrm{latent}}={\rho }_{i}{L}_{f}q$ tends to be small compared to the other terms, as can be seen from Figure 1(c), except when the moon size is very small and the ice-thickness variation is very large. For simplicity, I drop ${{ \mathcal H }}_{\mathrm{latent}}$. The remaining terms can all be written as a function of the ice topography H. ${{ \mathcal H }}_{\mathrm{ice}}$ amplifies over regions with thinner ice due to the rheology feedback (Beuthe 2019; Kang & Flierl 2020) following ${{ \mathcal H }}_{\mathrm{ice}}(\phi )\,={{ \mathcal H }}_{\mathrm{ice}0}(\phi )\times {\left({H}_{0}/H(\phi )\right)}^{2}$, where H₀ is the mean ice thickness and ${{ \mathcal H }}_{\mathrm{ice}0}(\phi )$ is the ice-dissipation rate in a flat ice shell as a function of latitude ϕ. The conductive heat loss ${{ \mathcal H }}_{\mathrm{cond}}$ is inversely proportional to the local ice thickness ${{ \mathcal H }}_{\mathrm{cond}}(\phi )\,={{ \mathcal H }}_{\mathrm{cond}0}\times ({H}_{0}/H(\phi ))$. Because 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 ${{ \mathcal H }}_{\mathrm{ice}0}(\mathrm{pole})=1.25\overline{{{ \mathcal H }}_{\mathrm{ice}0}}$ for the polar box and ${{ \mathcal H }}_{\mathrm{ice}0}(\mathrm{eq})\,=0.75\overline{{{ \mathcal H }}_{\mathrm{ice}0}}$ for the equatorial box, where $\overline{(\cdot )}$ denotes the global mean. To guarantee the global heat budget balance, $\overline{{{ \mathcal H }}_{\mathrm{ice}0}}$ should be equal to ${{ \mathcal H }}_{\mathrm{cond}0}\mathop{=}\limits^{{\rm{\Delta }}}{ \mathcal H }$. Now, consider the equator-to-pole difference of the heat budget terms for an ice shell that has a mean thickness of H₀ and an equator-to-pole thickness difference of ΔH (the equatorial ice shell is thicker). Keeping the first-order terms, the heat budget simplifies to

$$\begin{eqnarray}&&\begin{array}{l}2{{ \mathcal H }}_{\mathrm{ocn}}={\rm{\Delta }}{{ \mathcal H }}_{\mathrm{ice}}-{\rm{\Delta }}{{ \mathcal H }}_{\mathrm{cond}}\\ \approx 1.25{ \mathcal H }{\left(1-\displaystyle \frac{{\rm{\Delta }}H}{2{H}_{0}}\right)}^{-2}\\ \ \ -\,0.75{ \mathcal H }{\left(1+\displaystyle \frac{{\rm{\Delta }}H}{2{H}_{0}}\right)}^{-2}-\left[{ \mathcal H }{\left(1-\displaystyle \frac{{\rm{\Delta }}H}{2{H}_{0}}\right)}^{-1}\right)\\ \left.-{ \mathcal H }{\left(1+\displaystyle \frac{{\rm{\Delta }}H}{2{H}_{0}}\right)}^{-1}\right].\end{array}\end{eqnarray} \tag{ 29 }$$

From Equations (29) and (27), ΔH can be solved numerically. The results are presented in Figure 4. The entire solution falls into the κ_v -limit regime, as denoted by the red shading (blue shading marks the D-limit regime).

If the ice-thickness variation is small (ΔH/H₀ ≪ 1), further simplification can be made to Equation (29). Taylor-expanding Equation (29) around small ΔH/H₀ and dropping high-order terms give

$$\begin{eqnarray}&&2{{ \mathcal H }}_{\mathrm{ocn}}\approx \displaystyle \frac{1}{2}{ \mathcal H }.\end{eqnarray} \tag{ 30 }$$

From Equations (30) and (27), ΔH can be solved analytically:

$$\begin{eqnarray}{\rm{\Delta }}H\approx \left\{\begin{array}{ll}{\chi }_{\kappa }{\kappa }_{v}^{-1/2}{\alpha }^{-3/10}{f}^{2/5}{a}^{-7/10}, & \mathrm{for}\,\mathrm{the}\ {\kappa }_{v}\,\mathrm{limit}\\ {\chi }_{D}{\kappa }_{v}^{-3/13}{\alpha }^{-6/13}{f}^{8/13}{a}^{-7/13}, & \mathrm{for}\,\mathrm{the}\ D\,\mathrm{limit}\end{array}\right.\end{eqnarray} \tag{ 31 }$$

where ${\chi }_{\kappa }={\left(\tfrac{{ \mathcal H }}{8\rho {C}_{p}}\right)}^{7/10}{\left(4\pi G{\rho }_{b}/3\right)}^{-13/10}{k}^{-1/5}{\left(b{\rho }_{i}\right)}^{-1}$ and ${\chi }_{D}\,={\left(\tfrac{{ \mathcal H }}{8\rho {C}_{p}D}\right)}^{7/13}{\left(4\pi G{\rho }_{b}/3\right)}^{-19/13}{k}^{-4/13}{\left(b{\rho }_{i}\right)}^{-1}$ are constants. To obtain the above solution, I have substituted g with 4π G ρ_bulk a/3, where ρ_bulk = 2500 kg m⁻³ is the bulk density.

From Equation (31), it can be seen that, when ΔH ≪ H₀, the equilibrium ice-shell-thickness variation ΔH should decrease with the moon's size and increase 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 ${\rm{\Delta }}H\propto {\left(f/a\right)}^{2/3}$ or f/a^1/2 for the κ_v limit and D limit, respectively, as shown by Kang & Jansen (2022), except that the sensitivity to rotation rate is slightly lower here.

The asymptotic scalings (shown by blue and red dashed lines in Figure 4) provide a useful approximation to the full solution of Equation (29) for a relatively small to moderate ΔH. For larger ΔH, the sensitivity of ΔH on a increases, and eventually, runaway 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 Equation (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.

Compared to Enceladus, Europa is six times larger in size, is rotating three times slower, and its ocean is likely saltier (Hand & Chyba 2007; Zolotov 2007; Zolotov & Postberg 2014; McKay et al. 2018; Kang et al. 2021) and deeper (Hand & Chyba 2007)—all of these differences suggest that Europa's ice shell may be much flatter than Enceladus'. Assuming ∣α∣ ∼ 10⁻⁵–10⁻⁴/K (Figure 1(e)), κ_v = 10⁻³ m² s⁻¹, and γ = 10⁻⁴ m s⁻¹, Equation (29) yields ΔH = 8 km for Enceladus, ⁴ in line with the strong ice topography in observations (Iess et al. 2014; Beuthe et al. 2016; Tajeddine et al. 2017; Čadek et al. 2019; Hemingway & Mittal 2019). In contrast, using Europa parameters, ΔH is estimated to be only 0.7 km! This roughly matches the constraint based on limb profile measurements (Nimmo et al. 2007). Finally, it should be noted that ΔH is sensitive to the vertical diffusivity/viscosity and is only valid under our assumptions (shell heating, temperature-dominant density variation, Maxwell ice rheology, etc.). Better understanding of the dissipative processes in the ocean driven by tides and libration motions (Rekier et al. 2019) is required to better constrain the ice-shell geometry.

6.2. 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 ${{ \mathcal H }}_{\mathrm{ice}}$ (given by Equation (A9) in Appendix A), the down-gradient ice flow ${ \mathcal Q }$ (given by Equation (A7) in Appendix A), the heat loss to space by conduction ${{ \mathcal H }}_{\mathrm{cond}}$ (given by Equation (A4) in the appendix), and the heat transmitted upward by the ocean ${{ \mathcal H }}_{\mathrm{ocn}}$, and evolves the ice thickness H over time. The total thickness tendency can be symbolically expressed as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dH}}{{dt}} & = & \displaystyle \frac{{{ \mathcal H }}_{\mathrm{cond}}(H)-{{ \mathcal H }}_{\mathrm{ice}}(H)-{{ \mathcal H }}_{\mathrm{ocn}}(H)}{{L}_{f}{\rho }_{i}}\\ & & +\displaystyle \frac{1}{a\sin \phi }{\partial }_{\phi }\left(\sin \phi { \mathcal Q }(H)\right),\end{array}\end{eqnarray} \tag{ 32 }$$

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 ${{ \mathcal H }}_{\mathrm{ocn}}$ is a new component that does not exist in Kang & Flierl (2020). Inspired by the conceptual model, ${{ \mathcal H }}_{\mathrm{ocn}}$ is parameterized as

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal H }}_{\mathrm{ocn}}=\mathrm{MIN}\left\{2\rho {C}_{p}{\rm{\Delta }}T{\left(\displaystyle \frac{k}{{a}^{3}{f}^{2}}\right)}^{2/7}{\left(\alpha g{\rm{\Delta }}T\right)}^{3/7}{\kappa }_{v}^{5/7},\right.\\ \left.2\rho {C}_{p}D{\rm{\Delta }}T{\left(\displaystyle \frac{k\beta }{{a}^{2}{f}^{3}}\right)}^{4/7}{\left(\alpha g{\rm{\Delta }}T\right)}^{6/7}{\kappa }_{v}^{3/7}\right\}\end{array}\end{eqnarray} \tag{ 33 }$$

where $H^{\prime}$ is the deviation from the prescribed global-mean ice thickness H₀ = 20 km and $\mathrm{MIN}\{\}$ selects whichever parameterization yields a smaller global standard deviation. A factor of 2 is multiplied to H' because the analysis before uses the 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 OHT, a range of ice viscosities η_m , ⁵ ∣α∣ , and κ_v are considered for Enceladus and Europa. The equilibrium ice-shell geometries are shown in Figure 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 OHT, and significant ice-thickness variations develop on both Enceladus and Europa. When the ice viscosity is not too low (η_m > 10⁻¹³ Pa·s), ice is almost completely melted over one or both poles due to the ice-rheology feedback. However, with OHT, the equilibrium ice geometry is largely flattened especially for Europa given its large size and slower rotation rate. The smoothing effect of the 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 the ice-thickness variation on Europa is likely less than 1 km. When Enceladus' parameters are used instead, the OHT'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 of understanding 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 in geometry and the nonlinear behavior of eddies at high amplitude. Among the 28 scenarios considered for Enceladus, 5 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 the presence of ocean heat redistribution.

The two scientific questions addressed here are (1) how the efficiency of the 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's 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 (Carr et al. 1998; Kivelson et al. 2000; Hand & Chyba 2007; Postberg et al. 2009; Thomas et al. 2016). 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 OHT 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 (Nimmo et al. 2007; Iess et al. 2014; Beuthe et al. 2016; Tajeddine et al. 2017; Čadek et al. 2019; Hemingway & Mittal 2019).

To drive the point home, the equilibrium ice geometry for Enceladus and Europa is 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 2 km. 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. 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 the 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 decreases with gravity, ⁶ whereas the temperature difference induced by the ice topography will increase with satellite size; even for an icy moon as small as Enceladus, the vertical temperature gradient induced by a 40 mW m⁻² 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 OHT 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 constraint to the poorly constrained icy worlds.

This work is carried out in the Department of Earth, Atmospheric and Planetary Science (EAPS) in MIT. W.K. acknowledges support as a Lorenz-Houghton Fellow by endowed funds in EAPS and helpful comments from Prof. Malte Jansen

Software: MITgcm (MITgcm-group collaboration 2010)

Our simulations are carried out using the Massachusetts Institute of Technology OGCM (MITgcm Marshall et al. 1997; MITgcm-group collaboration 2010) configured for application to icy moons.

The model integrates the nonhydrostatic 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 collaboration (2010) and Marshall et al. (1997). Such terms are typically ignored when simulating Earth's ocean because the ratio between the fluid depth and the horizontal scale is small. Instead, when the moon size is on the order of hundreds of kilometers like Enceladus, the aspect ratio is on the order of 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:

$$\begin{eqnarray}&&\begin{array}{l}g(z)\\ =\displaystyle \frac{4\pi G\left[{\rho }_{\mathrm{core}}{\left(a-{D}_{0}-{H}_{0}\right)}^{3}+{\rho }_{\mathrm{out}}({\left(a-z\right)}^{3}-{\left(a-{D}_{0}-{H}_{0}\right)}^{3})\right]}{3{\left(a-z\right)}^{2}}.\end{array}\end{eqnarray} \tag{ A1 }$$

In the above equation, G = 6.67 × 10⁻¹¹ N m⁻² kg⁻² is the gravitational constant, ${\rho }_{\mathrm{core}}=2500\,\mathrm{kg}$ m⁻³ is the assumed core density, and ρ_out = 1000 kg m⁻³ is the density of the ocean/ice layer. D₀ and H₀ are the thickness of ocean and ice on a global average, respectively. Because 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° (8.7 km). Only 30 layers (each 2 km) are used to keep the computational cost manageable. After equilibrium is reached, I interpolate the pickup files to generate initial conditions for the corresponding 3D simulation, which has a default horizontal resolution of 025×025 and 70 unevenly distributed vertical layers, whose thicknesses increase from 500 m to 2 km from top to bottom. Because the changing rotation rate has a significant impact on the size of the jets and eddies, I had to adjust the grid width along the east–west direction by a factor of 1.4 in those cases to better capture the dynamics. The model parameters are summarized in Table 1.

A.1. Diffusivity and Viscosity

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 determining density using temperature, salinity, and pressure. As demonstrated by Figure 1(e) 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:

$$\begin{eqnarray}&&{T}_{f}(S,P)={c}_{0}+{b}_{0}P+{a}_{0}S,\end{eqnarray} \tag{ A2 }$$

where a₀ = −0.0575 K/psu, b₀ = −7.61 × 10⁻⁴ K/dbar, and c₀ = 0.0901 °C. The pressure P can be calculated using hydrostatic balance P = ρ_i gH (ρ_i = 917 kg m⁻³ 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 hydrostasy (i.e., ice is floating freely on the water). The ice thickness is set to be

$$\begin{eqnarray}&&H(\phi )={H}_{0}-{H}_{2}{P}_{2}(\sin \phi ),\end{eqnarray} \tag{ A3 }$$

where H₀ is the mean ice thickness, P₂ is the second-order Legendre polynomial, and H₂ is the amplitude of the ice-thickness variation. ϕ denotes latitude. The thickness profile is shown by a solid curve in Figure 1(b) of the main text. Partial cells are 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 are taken care of by a modified version of 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 ${{ \mathcal H }}_{\mathrm{cond}}$ is represented using a 1D vertical heat conduction model,

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{cond}}=\displaystyle \frac{{\kappa }_{0}}{H}\mathrm{ln}\left(\displaystyle \frac{{T}_{f}}{{T}_{s}}\right),\end{eqnarray} \tag{ A4 }$$

where H is the thickness of ice (solid curve in Figure 1(b) 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 (Equation (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 ${{ \mathcal H }}_{\mathrm{cond}}$ = 50 mW m⁻², broadly consistent with observations (Tajeddine et al. 2017).

Ice–ocean fluxes: The interaction between the ocean and ice is simulated using MITgcm's "shelf-ice" package (Holland & Jenkins 1999; Losch 2008) 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/freshwater 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 (Equation (A2)):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dS}}_{\mathrm{ocn}-\mathrm{top}}}{{dt}} & = & \displaystyle \frac{{{qS}}_{\mathrm{ocn}-\mathrm{top}}}{\delta z}\\ \displaystyle \frac{{{dT}}_{\mathrm{ocn}-\mathrm{top}}}{{dt}} & = & \displaystyle \frac{1}{\delta z}({\gamma }_{T}-q)({T}_{f}-{T}_{\mathrm{ocn}-\mathrm{top}}).\end{array}\end{eqnarray} \tag{ A5 }$$

Here, ${S}_{\mathrm{ocn}-\mathrm{top}}$ and ${T}_{\mathrm{ocn}-\mathrm{top}}$ denote the upper-boundary salinity and temperature, γ_T = 10⁻⁵ m s⁻¹ is the water–ice exchange coefficient 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 that 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⁻³ 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 very 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 ${ \mathcal Q }$,

$$\begin{eqnarray}&&{ \mathcal Q }(\phi )={{ \mathcal Q }}_{0}{H}^{3}({\partial }_{\phi }H/a),\end{eqnarray} \tag{ A6 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal Q }}_{0}=\displaystyle \frac{2({\rho }_{0}-{\rho }_{i})g}{{\eta }_{m}({\rho }_{0}/{\rho }_{i}){\mathrm{log}}^{3}\left({T}_{f}/{T}_{s}\right)}{\int }_{{T}_{s}}^{{T}_{f}}{\int }_{{T}_{s}}^{T(z)}\\ \quad \times \ \exp \left[-\displaystyle \frac{{E}_{a}}{{R}_{g}{T}_{f}}\left(\displaystyle \frac{{T}_{f}}{T^{\prime} }-1\right)\right]\mathrm{log}(T^{\prime} )\,\displaystyle \frac{{dT}^{\prime} }{T^{\prime} }\,\displaystyle \frac{{dT}}{T}.\end{array}\end{eqnarray} \tag{ A7 }$$

Here, ϕ denotes latitude, a and g are the radius and surface gravity of the moon, T_s and T_f are the temperatures at the ice surface and the water–ice interface (equal to the local freezing point, Equation (A2)), and ρ_i = 917 kg m⁻³ and ρ₀ 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¹³–10¹⁶ Pa·s Tobie et al. 2003), and here η_m is set to 10¹⁴ Pa·s.

In steady state, the freezing rate q must equal the divergence of the ice transport thus:

$$\begin{eqnarray}&&q=-\displaystyle \frac{1}{a\cos \phi }\displaystyle \frac{\partial }{\partial \phi }(Q\cos \phi ).\end{eqnarray} \tag{ A8 }$$

As shown by the dashed curve in Figure 1(b) 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 lead to changes in local salinity and thereby a buoyancy flux.

A.5. Model of Tidal Dissipation in the Ice Shell

An 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 ${{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{mem}}$ due to the extension/compression and tangential shearing of the ice membrane, a mixed mode ${{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{mix}}$ due to vertical shifting, and a bending mode ${{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{bend}}$ 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 ${{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{flat},{\rm{x}}}$ (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¹³ Pa·s.

Ice-thickness variations are accounted for by multiplying the membrane-mode dissipation ${{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{flat},\mathrm{mem}}$ by a factor that depends on ice thickness. The membrane mode is the only mode that is amplified in thin-ice regions (see Beuthe 2019). This results in the expression

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{ice}}={(H/{H}_{0})}^{{p}_{\alpha }}{{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{flat},\mathrm{mem}}+{{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{flat},\mathrm{mix}}+{{ \mathcal H }}_{\mathrm{ice}}^{\mathrm{flat},\mathrm{bend}},\end{eqnarray} \tag{ A9 }$$

where H is the prescribed thickness of the ice shell as a function of latitude and H₀ is the global mean of H. Because 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 Figure 1(c) of the main text.

Here I provide a brief overview of 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, the ice-shell thickness H is changed over time by the melting induced by the tidal heating ${{ \mathcal H }}_{\mathrm{ice}}$ (given by Equation (A9)), the down-gradient ice flow ${ \mathcal Q }$ (given by Equation (A6)), the heat loss to space by conduction ${{ \mathcal H }}_{\mathrm{cond}}$ (given by Equation (A4)), the crack-induced cooling ${{ \mathcal H }}_{\mathrm{crack}}$ in places where the ice is sufficiently thin, and heat transmitted upward by the ocean ${{ \mathcal H }}_{\mathrm{ocn}}$. The ice-thickness tendency can be symbolically expressed as follows:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dH}}{{dt}} & = & \displaystyle \frac{{{ \mathcal H }}_{\mathrm{cond}}(H)-{{ \mathcal H }}_{\mathrm{ice}}(H)-{{ \mathcal H }}_{\mathrm{ocn}}}{{L}_{f}{\rho }_{i}}\\ & & +\displaystyle \frac{1}{a\sin \phi }{\partial }_{\phi }\left(\sin \phi { \mathcal Q }(H)\right),\end{array}\end{eqnarray} \tag{ B1 }$$

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. ${{ \mathcal H }}_{\mathrm{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 Equation (A9)). The tendency for the ice-thickness variation to increase due to the rheology feedback will be balanced by the rapid heat loss through thin ice (Equation (A4)) and the transport by ice flow (Equation (A6)). 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 ${{ \mathcal H }}_{\mathrm{ice}}$ is scaled to exactly balance the instantaneous conductive heat loss ${{ \mathcal H }}_{\mathrm{cond}}$. By so doing, the rheology feedback and thus the ice-thickness variation is maximized. Throughout the integration, the global-mean ice thickness is fixed at H₀ = 20 km.

Table 2. Parameters for Enceladus and Europa

Symbol	Name	Definition/Value
pα	Ice-dissipation amplification factor	−1.5
Parameters for Enceladus
a	Radius	252 km
g0	Surface gravity	0.113 m s−2
δ	Obliquity	27°
H0	Global-mean ice thickness	20 km Hemingway & Mittal (2019)
H2	Initial equator-to-pole ice-thickness variation	3 km
H1	Initial hemispherical asymmetry	1 km
$\bar{{T}_{s}}$	Mean surface temperature	59 K
Parameters for Europa
a	Radius	1561 km
g0	Surface gravity	1.315 m s−2
δ	Obliquity	31
H0	Global-mean ice thickness	20 km
H2	Initial equator-to-pole ice-thickness variation	3 km
H1	Initial hemispherical asymmetry	−1 km
$\bar{{T}_{s}}$	Mean surface temperature	110 K

Download table as:  ASCIITypeset image

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 Equation (19) in the main text.

The initial condition is set as follows

$$\begin{eqnarray}&&H(0)={H}_{0}-{H}_{2}{P}_{2}(\sin \phi )-{H}_{1}{P}_{1}(\sin \phi ),\end{eqnarray} \tag{ B2 }$$

where H₀ = 20 km, H₂ = 3 km, and H₁ = 1 km. P₁ and P₂ are the first and second order of Legendre polynomials.

Figures A1–A6 present 2D and 3D ocean dynamics and OHT for experiments with different satellite radius, rotation periods and different equator-to-pole ice thickness gradients. Figure A7 compares the simulated jet speeds, jet spacings and eddy sizes with the theoretical prediction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size