How does ice shell geometry shape ocean dynamics on icy moons?

TL;DR

This study reveals that the geometry of an icy moon’s ice shell induces a poleward under-ice temperature gradient, which translates into a buoyancy contrast that drives an eddy-dominated, baroclinically unstable ocean circulation. By developing a scaling framework that couples isopycnal slope and eddy transport, the authors predict meridional heat flux and stratification across a broad parameter space, including surface/topographic effects and vertical diffusion, and validate these scalings against high-resolution simulations. Applying the framework to Enceladus, Europa, and Titan, they quantify how topography and salinity govern ocean heat transport, stratification, tidal heating, and the potential for observable ice-shell rotation changes. The results provide a pathway to constrain subsurface-ocean properties from ice thickness observations and future mission data, linking ocean dynamics to tidal dissipation and habitability considerations.

Abstract

A poleward-thinning ice shell can drive circulation in the subsurface oceans of icy moons by imposing a meridional temperature gradient--colder at the equator than the pole--through the freezing point suppression due to pressure. This temperature gradient sets a buoyancy gradient, whose sign depends on the thermal expansion coefficient determined by ocean salinity. Together with vertical mixing, this buoyancy forcing shapes key oceanic features, including zonal currents in thermal wind balance, baroclinic instability of those currents, meridional heat transport by eddies, and vertical stratification. We use high-resolution numerical simulations to explore how variations in ice shell thickness affect these processes. Our simulations span a wide range of topographic slopes, pole-to-equator temperature differences, and vertical mixing strengths, for both fresh and salty oceans. We find that baroclinic eddies dominate large-scale circulation and meridional heat transport, consistent with studies assuming a flat ice-ocean interface. However, sloped topography introduces new effects: when lighter water overlies denser water along the slope, circulation weakens as a stratified layer thickens beneath the poles. Conversely, when denser water lies beneath the poles, circulation strengthens as topography increases the available potential energy. We develop a scaling framework that predicts heat transport and stratification across all simulations. Applying this framework to Enceladus, Europa, and Titan, we infer ocean heat fluxes, stratification, and tidal energy dissipation and showing large-scale circulation constrains tidal heating and links future observations of ice thickness and rotation to subsurface ocean dynamics.

Figures (7)

• Figure 1: Setup of our numerical simulations. Panel (A) illustrates an idealized icy moon with a poleward-thinning ice shell. This ice shell thickness gradient creates a meridional temperature gradient at the water-ice interface ("Cold" and "Warm" in Panel A) due to freezing point suppression by high pressure. The resulting temperature contrast $\Delta T$ leads to a buoyancy contrast $\Delta b = \alpha g \Delta T$, where $\alpha$ is the thermal expansion coefficient and $g$ is the gravitational acceleration. The sign of $\Delta b$ depends on the sign of $\alpha$. This buoyancy gradient drives ocean circulation, which, under background rotation (indicated by arrows in Panel A), gives rise to baroclinic eddies -- referred to as "ocean weather systems" in zhang2024ocean. Panel (B) depicts the numerical setup. Simulations are conducted in a Cartesian coordinate system. The poleward-thinning ice shell is represented as a tilted upper boundary with an elevation difference $\Delta H$ (see Eq. \ref{['eq:z-top']}). A buoyancy contrast, consistent with the imposed temperature difference, is prescribed as a boundary condition at the ocean top (Eq. \ref{['eq:b-top']}) and diffuses into the ocean interior via vertical diffusion ($-\kappa_v \partial_z b$ in Panel B). The Coriolis force is implemented in Cartesian coordinates by conserving the angle between the rotation vector $\mathbf{f}$ (green arrows in Panel B) and gravity, following bire2022exploring, zhang2024ocean. Parameters for each simulation are listed in Table \ref{['tab:simulation']}.
• Figure 2: Ocean circulation and buoyancy transport in simulations h4$^+$ and h4$^-$ (see Table \ref{['tab:simulation']} for parameters). Panel A1 corresponds to h4$^+$; Panel A2 to h4$^-$. In Panels A1 and A2, solid contours show buoyancy (contour interval = $|\Delta b|/10$), with buoyant water overlying denser water. Green and purple arrows denote eddy buoyancy flux $(\overline{v'b'}, \overline{w'b'})$ and diffusive flux $(-\kappa_h \partial_y \overline{b}, -\kappa_v \partial_z \overline{b})$, normalized by the maximum $\overline{v'b'}$ (listed top right). The equivalent overturning circulation, diagnosed from Eq. \ref{['eq:diagnose-psi-star']}, is shown by $\psi^\star$ shading (normalized by its maximum; values listed). Negative $\psi^\star$ (blue) indicates counterclockwise overturning; positive (red) indicates clockwise overturning. Panel B shows vertically integrated buoyancy transport: solid lines for h1$^+$--h4$^+$ and c1$^+$ (varying $\Delta H$; legend), dashed lines for h1$^-$--h4$^-$ and c1$^-$.
• Figure 3: Zonal-mean buoyancy and buoyancy flux patterns from our numerical simulations. Each panel corresponds to one simulation (see Table \ref{['tab:simulation']} for parameters) and shows the same variables as Fig. \ref{['fig:phenomenology']}(A) The middle column shows the control experiment for $\Delta b>0$ (top 3 panels) and $\Delta b<0$ (bottom 3 panels). For convenience the panels presenting the control are repeated to ease comparison as parameters are changed moving left and right.
• Figure 4: Tests of scaling theory against numerical simulations. Panel (A) shows a scaling for the isopycnal slope (Eq. \ref{['eq:scaling-kappa-over-psi']}). Panel (B) and Panel (C) show the scaling for the eddy-driven overturning $\psi^\star$ for positive and negative $\Delta b$, respectively (Eq. \ref{['eq:scaling-psi']}). Panel (D) shows the scaling for the vertical integrated meridional buoyancy flux (proportional to heat flux). Here, $s$ represents the slope of the median isopycnal, $s_m$; $s_\mathrm{top}\equiv \Delta H / ((\pi/2)a)$ is the slope of the top topography; $\psi^\star$ represents the mean overturning diagnosed using Eq. \ref{['eq:diagnose-psi-star']}; $\mathrm{sgn}(\Delta b)$ is the sign of the prescribed equator-to-pole buoyancy difference $\Delta b$; $\mathcal{F}_b$ is the maximum of the vertical integrated meridional buoyancy flux. In Panel (B) and (C), $\psi^\star$ is normalized using the buoyancy difference $\Delta b$, the moon radius $a$, and the $\beta$ factor ($\beta \equiv 2\Omega / a$). The line and the formula at the upper left corner of each panel shows the best fit; the red dashed line and blue dotted lines represent the fit for positive and negative $\Delta$ cases, respectively. The marker and color of each simulation is the same as Fig. \ref{['fig:zonal-mean']}. The parameters used in these simulations are set out in Table \ref{['tab:simulation']}.
• Figure 5: Meridional buoyancy fluxes in our simulations compared against predictions from our scaling framework (Equations \ref{['eq:scaling-s-positive']}, \ref{['eq:scaling-s-negative']}, \ref{['eq:scaling-psi']}, and \ref{['eq:scaling-f']}). Panel (A) plots the predicted vertically integrated meridional buoyancy transport $\mathcal{F}_b$ against the measured value (defined as the maximum along all latitudes) from each of our simulations. The solid line and the formula on the left upper corner show the best fit, and is almost the 1:1 line. Panels (B,C,D) show how $\Delta H$, $\kappa_v$, and $|\Delta b|$ influences $\mathcal{F}_b$: the dashed and dotted lines represent the predictions of our scalings for positive and negative $\Delta b$, respectively. The parameters used in these simulations are tabulated in Table \ref{['tab:simulation']}.