# The Physical Oceanography of Ice-Covered Moons

Krista M Soderlund, Marc Rovira-Navarro, Michael Le Bars, Britney E Schmidt, Theo Gerkema

► To cite this version:

Krista M Soderlund, Marc Rovira-Navarro, Michael Le Bars, Britney E Schmidt, Theo Gerkema. The Physical Oceanography of Ice-Covered Moons. Annual Review of Marine Science, 2024, 16 (1), pp.25-53. 10.1146/annurev-marine-040323-101355. hal-04238916

HAL Id: hal-04238916

https://hal.science/hal-04238916v1

Submitted on 12 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
# The Physical Oceanography of Ice-Covered Moons

Krista M. Soderlund<fnref n="1" />, Marc Rovira-Navarro<fnref n="2" />, Michael Le Bars<fnref n="3" />, Britney E. Schmidt<fnref n="4" /> and Theo Gerkema<fnref n="5" />

Annu. Rev. Mar. Sci. 2023. xxx:1–31

https://doi.org/10.1146/((please add article doi))

Copyright © 2023 by the author(s).
All rights reserved

## Keywords

ocean worlds, physical oceanography, ocean dynamics, ice-ocean interactions

## Abstract

In the outer solar system, a growing number of giant planet satellites are abodes for global oceans hidden below an outer layer of ice. These planetary oceans are a natural laboratory for studying physical oceanographic processes in settings that challenge traditional assumptions made for the Earth's oceans. While some driving mechanisms are common to both systems, such as buoyancy-driven flows and tides, others such as libration, precession, and electromagnetic pumping are likely more significant for moons in orbit around a host planet. We here review these mechanisms and how they may operate across the solar system, including their implications for ice-ocean interactions. Future studies should not only continue to advance our understanding of each of these processes, but also how they may act together in concert. This interplay also has strong implications for habitability as well as testing oceanic hypotheses with future missions.

<footnote n="1">¹ Institute for Geophysics, Jackson School of Geosciences, University of Texas at Austin, USA, 78758; email: krista@ig.utexas.edu, ORCID:0000-0002-7901-3239</footnote>

<footnote n="2">² Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ, USA, 85721; email: mroviranavarro@arizona.edu, ORCID:0000-0002-9980-5065</footnote>

<footnote n="3">³ CNRS, Aix Marseille Univ, Centrale Marseille, IRPHE, Marseille, 13013, France; email: lebars@irphe.univ-mrs.fr, ORCID:0000-0002-4884-6190</footnote>

<footnote n="4">⁴ Departments of Astronomy and of Earth & Atmospheric Sciences, Cornell University, Ithaca, NY, USA, 14853; email: britneys@cornell.edu, ORCID:0000-0001-7376-8510</footnote>

<footnote n="5">⁵ NIOZ Royal Netherlands Institute for Sea Research, Department of Estuarine and Delta Systems EDS, P.O. Box 140, 4400 AC Yerseke, the Netherlands; email: theo.gerkema@nioz.nl, ORCID:0000-0002-7702-1788</footnote>
# 1. INTRODUCTION

For centuries, the subject of physical oceanography was limited to Earth's oceans and seas. The discovery of oceans buried beneath the icy crusts of moons in the outer solar system has widened the horizons beyond Earth. In the last decades, observations and process-based studies have started to unveil the complex and diverse dynamics possible in these distant oceans. Yet, many questions remain in this new era of planetary oceanography. The aim of this paper is to review our current understanding of the dynamics of ice-covered oceans and underline the main outstanding questions.

To study ocean dynamics, it is necessary to know their physical characteristics and understand their geophysical context. Data about extraterrestrial oceans, however, is scarce. For some moons, observations have provided unambiguous evidence of the existence of subsurface oceans (Europa, Ganymede, Enceladus, and Titan); for others, observations or theoretical models either suggest or allow for possible oceans (Figure 1). Our knowledge of the characteristics of such oceans is limited (at best) to constraints on ocean thickness and composition, which affects efforts to understand planetary ocean dynamics. Section 2 provides a concise overview of the characteristics of ocean-bearing moons and their interiors. An in-depth review on the topic can be found in Nimmo & Pappalardo (2016).

Beyond the physical characteristics of ice-covered oceans, we can infer their dynamics using computational and analytical models, laboratory experiments, and Earth analogs

<figure><img src="image_2.png" /><figcaption>Figure 1: Confirmed, tentative, and theoretical ocean worlds in the outer solar system. For the confirmed ocean worlds, internal structure models were adapted from Vance et al. (2018). Additional information about the moons can be found in Supplementary Tables 1 and 2.</figcaption></figure>
(both the oceans and the core). However, laboratory experiments are limited by scale, Earth analogs by sometimes stark differences between Earth's ocean and those of icy ocean worlds, and numerical ocean models by existing computational resources. For the latter, this makes it impossible to directly solve the equations governing ocean dynamics. Instead, different assumptions are employed to simplify the system and study specific processes. In Section 3, we review the equations governing ocean dynamics and the simplifying approximations often used.

Different processes can result in ocean currents. We distinguish between buoyancy- (Section 4), mechanically- (Section 5), and electromagnetically- (Section 6) driven flows. Buoyancy-driven flows are caused by the exchange of heat and salt between the ocean and the underlying mantle/core and overlying ice shell. The decay of radioactive material and tidal heating within the mantle results in a basal heat flux. Salts and other impurities may be introduced by water-rock reactions below the seafloor, while melting and/or freezing of the overlying ice shell will freshen or introduce salts locally. These processes alter the density of the ocean and drive ocean circulation. Mechanically-driven flows are the consequence of the moons' orbit and their rotation. Icy moons experience gravitational tides that deform them and stir their oceans. Tidal torques also cause periodic changes in the rotation of the moons (e.g., libration, precession) that can drive ocean currents. Finally, electromagnetic interactions with their host planet's magnetic field is an additional mechanism for ocean currents. We finish by providing perspectives from Earth (Section 7), summarizing the key points, and highlighting future issues for the field.

## 2. ICE-COVERED OCEAN WORLDS

### 2.1. The Jovian System

The *Voyager* and *Galileo* missions revealed that the Jovian moons exhibit a gradient in geological activity. Io, the rocky innermost moon, is the most volcanically active world in the Solar System and is devoid of impact craters; Europa's icy surface has few craters, exhibits global resurfacing, and is tectonized; Ganymede presents heavily cratered plains interrupted by bright, less cratered younger terrain textured by ridges and troughs; and Callisto is heavily cratered and shows little evidence of geological activity (e.g., Schenk 2010). The gradient in geological activity is linked to the moons' orbits. The three innermost Galilean moons participate in a Laplace resonance (Peale 1999). The moons mutually exert a coherent periodic gravitational perturbation to each other that forces their orbital eccentricity. The non-zero eccentricity results in tidal forces; dissipation converts part of tidal energy into internal heat — a process known as tidal heating. Io experiences the highest amount of tidal heating followed by Europa; Ganymede and Callisto are not experiencing significant tidal heating at present but Ganymede might have in the past (Downey et al. 2020, Bland et al. 2009).

The density of the Galilean moons decreases with distance away from Jupiter; the water/rock mass ratio changes from 5% for Europa to roughly 50% for Callisto (see Soderlund et al. 2020, and references therein). Geodetic measurements by the *Galileo* mission indicated that Europa has a 80 – 170 km thick H₂O layer above a rocky mantle that could either be enriched in metals or fully differentiated into a silicate envelope and a metallic core, with the second scenario favoured. Ganymede is differentiated into a core with a radius of ~500 – 1000 km, surrounded by a silicate mantle and a H₂O layer ~800 km thick. Callisto may not be fully differentiated; gravity data is consistent with both two-
layered models consisting of a metal-rock-ice core beneath a ~350 km hydrosphere and three-layered models where the metallic elements have separated from the silicates.

The small density contrast between liquid water and ice precluded deducing the state of the hydrosphere from gravity data alone. The detection of induced magnetic fields around the moons by *Galileo* demonstrated the presence of an electrically conductive layer — i.e., a salty ocean— beneath the surfaces of Europa (Kivelson et al. 2000), Ganymede (Kivelson et al. 2002), and Callisto (Khorana et al. 1998). A Europan subsurface ocean is also consistent with its surface geology (e.g., Pappalardo et al. 1999). The case of a subsurface ocean in Ganymede is further supported by *Hubble Space Telescope* aurora observations (Saur et al. 2015); for Callisto, however, the magnetic signal may also be caused by its ionosphere (Hartkorn & Saur 2017). Europa's subsurface ocean is in direct contact with the silicate shell. In contrast, the high pressures attained within Ganymede's hydrosphere results in the formation of high-pressure ices (Journaux et al. 2020).

Ocean thickness and electrical conductivity cannot be inferred independently from existing magnetic data, and the dominant salts and salinities of subsurface oceans are poorly constrained (see Soderlund et al. 2020, for a review). For Europa, the preponderance of geologic evidence suggests ice shells in the tens of kilometers (e.g., Pappalardo et al. 1999), though other estimates range between a few km (Walker & Rhoden 2022) to more than ~50 km (Vilella et al. 2020). In the absence of significant tidal heating at present, the outer shells of Ganymede and Callisto are likely thicker. The feedback between interior and orbital dynamics likely drives temporal variations of shell thickness over geological timescales (e.g., Hussmann & Spohn 2004), and potentially spatial variations as well (e.g., Ojakangas & Stevenson 1989).

## 2.2. The Saturnian System

Titan represents roughly 95% of the mass of Saturn's moons, with the remaining in a collection of smaller icy satellites. How this system architecture formed is still debated and has important implications for the moon ages (e.g., Sasaki et al. 2010, Ćuk et al. 2016). The *Cassini-Huygens* mission provided constraints on the interior structure of Titan and characterized its methane-based hydrological cycle (see Sotin et al. 2021, for a review). The density and moment of inertia indicate that Titan has a ~400 km H₂O thick layer above a core with a density of ~2500 kg m⁻³. Measurements of Titan's tidal response and obliquity indicate that it has an ocean with a density higher than 1100 and as high as 1350 kg m⁻³. The lower end can be met by pure water with 3 wt% ammonia; higher densities could be explained by the presence of ammonium sulfate, magnesium sulfate, or chlorides. Shape and gravity data indicate an ice shell thickness in the range of 40–100 km, with potential variations in thickness or density. Shell thickness variations could be caused by non-uniform heating and be as high ~12 km, while density variations could be driven by Titan's methane cycle. Thickness variations are suggestive of a cold conductive shell, consistent with the limited tidal heating required to sustain the moon's relatively high free (i.e., not forced via a resonance) orbital eccentricity. These variations may suggest that oceanic heat flux into the ice shell peaks at high latitudes with a minima near the equator (Kvorka et al. 2018).

Saturn's mid-sized moons are compositionally similar with varied geology (e.g., Castillo-Rogez et al. 2018). Mimas and Rhea show no sign of recent activity, Tethys and Dione have evidence for past resurfacing, while Enceladus features a somewhat two-faced surface, with
tectonized young regions and older cratered plains. Despite being just 250 km in radius (as opposed to 2574 km for Titan and 1563 km for Europa), Enceladus is highly active, with plumes of water venting hundreds of km above its south polar limb (Porco et al. 2006). *Cassini* showed that Enceladus possesses a porous, low density core ($\rho \sim 2300 \text{ kg m}^{-3}$) below a global ocean and an ice shell (see Hemingway et al. 2018, for a review). Shape and gravity modeling suggest that heat flux at the top of the ocean is characterized by a strong degree-2 zonal component (Čadek et al. 2019). The detection of salt grains, silicate particles, $H_2$, and $N_2$ in the material ejected by the plumes indicates that the subsurface ocean is in contact with the core and has ongoing hydrothermal activity (Hsu et al. 2015, Waite et al. 2017). Gravity and shape data are compatible with a subsurface ocean in Dione beneath a $\sim 60–140 \text{ km}$ thick shell (Zannoni et al. 2020). The librations of Mimas indicate the moon might have a subsurface ocean (e.g., Rhoden & Walker 2022), which contrasts with its inactive surface and high free eccentricity. The formation and longevity of the oceans of Saturn's mid-sized moons depend on their orbital evolution, which is still poorly understood (e.g., Neveu & Rhoden 2019).

## 2.3. Beyond Saturn

Without the benefit of dedicated missions to the Uranian and Neptunian systems, our knowledge of these moons widely derives from data returned by the *Voyager 2* mission (e.g., Smith et al. 1986, 1989). The Uranian moons have densities similar to those of Saturn's mid-sized moons. Miranda and Ariel experienced resurfacing events as recently as 100 Ma ago and heat fluxes ($\sim 10–100 \text{ mW m}^{-2}$) higher than those experienced by Europa and Enceladus at present (Beddingfield et al. 2015, Peterson et al. 2015). While none of the moons are in an orbital resonance at present, they have gone through different resonances in the past (e.g., Čuk et al. 2020). These resonances likely caused the resurfacing events and high heat fluxes recorded on Miranda and Ariel's surfaces and might have resulted in the formation of subsurface oceans. Subsurface oceans are also possible in Titania and Oberon from radiogenic heating alone (Bierson & Nimmo 2022). Ammonia has been detected in Miranda and Ariel (Bauer et al. 2002, Cartwright et al. 2020), which may help in the preservation of these oceans until present (cf. Castillo-Rogez et al. 2023).

Neptune's moon Triton differs from other icy satellites as it is a captured Kuiper Belt object in a retrograde orbit tilted with respect to the ecliptic (e.g., McKinnon 1984). Triton has a young surface with extensive resurfacing (e.g., Smith et al. 1989) and plumes that might be driven by solid-state greenhouse effect or internal heating (e.g., Hofgartner et al. 2022). Triton likely experienced strong tides during its capture that melted its interior (Ross & Schubert 1990), which could be maintained by a combination of radiogenic and tidal heating until present (e.g., Nimmo & Spencer 2015).

## 2.4. Earth From Afar

In modelling ice-covered ocean worlds, we are thus in a data-poor, possibility-rich universe. This is akin to modeling Earth's oceans without detailed knowledge of ocean salinity, structure, or bathymetry; a modeller could do little more than assume an unstratified or uniformly stratified ocean of nearly uniform depth. As terrestrials, we know how disparate these idealized assumptions are from what our real ocean looks like. It is illuminating to reflect on how this complex reality affects our ocean's dynamics, and to tentatively explore which of these factors may be at work in ice-covered ocean worlds.
<figure><img src="image_3.png" /></figure>

Figure 2: Buoyancy frequency *N* along a transect in the Pacific Ocean (for the definition of *N*, see Equation 3 below). This transect follows the meridian at longitude 179°E, except for an eastward bend in the Bering Sea. Notice that the color scale is logarithmic. Data supplied by the WOCE Hydrographic Program. Figure reproduced from Gerkema (2019).

By way of example, **Figure 2** shows a meridional transect from the Pacific Ocean. The bathymetry is uneven, exhibiting ridges and seamounts throughout, as well as a deep trench. The effect of bathymetry on ocean dynamics is wide-ranging. For example, continental slopes act as a wave guide for tides and shape the pattern of tides, including their resonance characteristics. Internal tides reflecting from a slope can result in intensification of internal tidal beams, making ridges and seamounts hotspots of internal tide breaking and shear instabilities, and hence of mixing. Abyssal mixing plays a key role in maintenance of the large-scale meridional overturning circulation (Wunsch & Ferrari 2004). Below the surface, both ice and seafloor topography might act as internal wave generators, and might trigger significant material and heat exchange. However, our ability to study these processes is hindered by our limited knowledge of such roughness.

Another conspicuous feature in **Figure 2** is the ocean's non-uniform stratification (*N*) and its horizontally layered structure. The top layer is shaped by solar radiation, and atmosphere-ocean interactions create an upper mixed layer in the upper 50 m or so (barely discernible in the figure). Beneath this thin layer, temperature undergoes a rapid change with depth in tropical regions and at mid-latitudes in summer. This layer of high *N* is called the thermocline or pycnocline, which occupies the upper few hundreds of meters. The deeper part of the ocean, beneath 1 km, always has temperatures lower than 5°C, even in the tropics. For moons of the outer planets, the structure of *N* is likely to be very different than that in **Figure 2**, and might instead resemble more exotic Earth environments in the polar regions where the stratifying component shifts from temperature to salinity (e.g., Vance & Goodman 2009, Soderlund et al. 2020, Lawrence 2022).

### 3. GOVERNING EQUATIONS AND COMMON APPROXIMATIONS

The equations governing ocean dynamics are the mass, momentum, and salinity conservation equations, the heat equation, an equation of state, Poisson's equation, the induction equation, and the no magnetic monopole requirement (Equation 1, respectively). The full equations are given in Supplementary Materials; below, we simplify using the anelastic approximation (which filters sound waves from the equations of motion) and assume a lin-
ear equation of state. Density $\rho$, pressure $p$, temperature $T$, potential temperature $\theta$, and salinity $S$ are written in terms of a known background assumed to be adiabatic, isohaline, and hydrostatic (denoted by a 0 subscript) and a perturbation (denoted by primes). The governing equations are then given by Ingersoll (2005):

$$ \nabla \cdot (\rho_0 \mathbf{u}) = 0, \qquad 1a. $$

$$ \frac{D\mathbf{u}}{Dt} + 2\mathbf{\Omega} \times \mathbf{u} = -\frac{1}{\rho_0} \nabla p' + \frac{\rho'}{\rho_0} \mathbf{g} - \nabla \phi' + \nu \nabla^2 \mathbf{u} + \frac{1}{\mu_0 \rho_0} (\nabla \times \mathbf{B}) \times \mathbf{B} - \frac{d\mathbf{\Omega}}{dt} \times \mathbf{r}, \quad 1b. $$

$$ \frac{DS'}{Dt} = Q_1, \qquad 1c. $$

$$ \frac{D\theta'}{Dt} = \frac{\theta_0 Q_2}{T_0 c_p}, \qquad 1d. $$

$$ \rho' = \rho_0(-\alpha\theta' + \beta S' + \frac{1}{\rho_0 c_s^2} p'), \qquad 1e. $$

$$ \nabla^2 \phi' = 4\pi G \rho', \qquad 1f. $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) - \nabla \times \left( \frac{1}{\mu_0 \sigma} \nabla \times \mathbf{B} \right), \qquad 1g. $$

$$ \nabla \cdot \mathbf{B} = 0. \qquad 1h. $$

Here, **u** is velocity field, **g** is gravitational acceleration, $\phi = \phi_0 + \phi'$ is gravitational potential with both static and perturbation components, $G$ is the gravitational constant, **B** is magnetic field, $\Omega$ is rotation rate of the reference frame of the moon, $D/Dt$ is material derivative, **r** is position vector, $\nu$ is kinematic viscosity, $\alpha$ is thermal expansion coefficient, $\beta$ is saline contraction coefficient, $c_p$ is specific heat capacity, $c_s$ is speed of sound, $\sigma$ is electrical conductivity, and $\mu_0$ is the vacuum magnetic permeability. $Q_1$ and $Q_2$, respectively, represent sinks and sources of salinity and heat, including irreversible (diffusive) processes (see Supplementary text). The domain, initial, and boundary conditions must also be specified; the latter couples the dynamics of the ocean to the ice shell and deep interior.

Commonly used approximations are listed in Table 1, such as the anelastic (ii) and Boussinesq (iii) approximations. Nearly all planetary oceanography studies neglect magnetic fields (i), eliminating Equations 1g-1h and the Lorentz force from (1b) (cf. Tyler 2011a, Gissinger & Petitdemange 2019, Vance et al. 2020). Additionally, in some studies, the thin-layer (vii) and spherical approximations (viii) are utilized, where both gravity and the centrifugal acceleration are assumed to be constant throughout the water column and its non-radial components are ignored (see Section 5).

Buoyancy- and mechanically- driven flows are often studied separately. Convective motions are the consequence of diabatic processes (here loosely understood as exchange of heat and salinity) and determine the stratification structure of the ocean. In contrast, to study mechanically-driven flows, diabatic processes are neglected and the stratification structure is taken as a given. It must be noted, however, that the two are closely linked.
<table><caption>Table 1: Common approximations. $c_s$ and $c_{surf}$ are the sound and surface wave speed, respectively. $\varepsilon$, $\Omega$ and $R$ are body's oblateness, rotational frequency, and radius. $H$ and $N$ are the ocean's thickness and buoyancy frequency. $\eta_b$ is the displacement of the ocean bottom due to mechanical forces, $\phi^{(SG)}$ is the gravitational potential arising from the body's self-gravitation. $L$ and $\omega$ are characteristic wavelength and frequency of the problem; for tides, they correspond to the body's radius $R$ and the tidal frequency.</caption><thead><tr><th>ID</th><th>Approximation</th><th>Typically used if</th><th>Filtered/Ignored Phenomena</th></tr></thead><tbody><tr><td>(i)</td><td>No magnetism</td><td>Electrically insulating</td><td>Magnetic induction, Lorentz force</td></tr><tr><td>(ii)</td><td>Anelastic approximation</td><td>ωL &ll; c<sub>s</sub></td><td>Sound waves</td></tr><tr><td>(iii)</td><td>Boussinesq approximation</td><td>(ii), ρ<sub>o</sub> is constant</td><td>No density variations beyond buoyancy</td></tr><tr><td>(iv)</td><td>Linear equation of state</td><td>α, β are constant</td><td>Cabbeling, thermobaricity</td></tr><tr><td>(v)</td><td>Isohaline</td><td>Q<sub>1</sub> &asymp; 0</td><td>Variations in salinity</td></tr><tr><td>(vi)</td><td>Isentropic</td><td>Q<sub>2</sub> &asymp; 0</td><td>Non-adiabatic, non-reversible processes</td></tr><tr><td>(vii)</td><td>Thin-layer approximation</td><td>H &ll; R</td><td>Radial gravity variations</td></tr><tr><td>(viii)</td><td>Spherical-body approximation</td><td>ε &ll; 1</td><td>Horizontal gravity component</td></tr><tr><td>(ix)</td><td>Mass-less approximation</td><td>φ<sup>(SG)</sup> &ll; φ<sub>0</sub></td><td>Self-gravity</td></tr><tr><td>(x)</td><td>Rigid-mantle/core</td><td>η<sub>b</sub> &ll; 1</td><td>No ocean bottom displacements</td></tr><tr><td>(xi)</td><td>Constant rotational rate</td><td>Ω = Ω<sub>0</sub></td><td>Poincaré force, spin-driven flows</td></tr><tr><td>(xii)</td><td>Linearization</td><td>Small amplitude waves</td><td>Turbulence, wave-breaking</td></tr><tr><td>(xiii)</td><td>Traditional approximation</td><td>H/L &ll; 1</td><td>Horizontal component rotation vector</td></tr><tr><td>(xiv)</td><td>Shallow-water approximation</td><td>H/L &ll; 1</td><td>Internal inertial waves</td></tr><tr><td>(xv)</td><td>Unstratified ocean</td><td>N &asymp; 0</td><td>Internal gravity waves</td></tr><tr><td>(xvi)</td><td>Non-rotating approximation</td><td>Ω/ω &ll; 1</td><td>Rossby waves</td></tr><tr><td>(xvii)</td><td>Equilibrium approximation</td><td>(xv), (xvi), ωL &ll; c<sub>surf</sub></td><td>All types of ocean waves</td></tr></tbody></table>

Mixing and dissipation due to mechanically-driven flows play a key role in the thermal (and haline) structure of the ocean. Moreover, coupling between the ocean and the ice shell and exchange between the two require understanding both simultaneously.

For buoyancy-driven flows, the criterion for convective instability requires the super adiabaticity of the fluid to exceed any compositional gradients: $\nabla_T - \nabla_{ad} > (\beta/\alpha)\nabla_S$, known as the Ledoux instability criterion. Here, $\nabla_T$, $\nabla_{ad}$, and $\nabla_S$ are the radial temperature, adiabatic temperature, and salinity gradients, respectively. In the absence of compositional gradients, this simplifies to the Schwarzschild criterion: $\nabla_T - \nabla_{ad} > 0$. For these flows generally, the perturbing potential is ignored and the ocean geometry is assumed to be fixed (ix,x), the rotational rate is assumed to be constant (xi), and $Q_2$ is assumed to consist of thermal diffusion only; other sources and sinks of heat (e.g., viscous and ohmic heating) are ignored but can still be present via the boundary conditions. $Q_1$ is kept as a generic sink/source of salinity, hence the difference in form compared to the heat equation. Furthermore, the Boussinesq approximation (iii) is often used instead of the anelastic approximation (ii). Here, the background density and potential temperature are assumed to be constant and the density perturbation to be independent of pressure perturbations (see Tritton (1998), chapter 14 appendix, for a detailed discussion). This reduces non-magnetic (i) Equation 1 to

$$ \nabla \cdot \mathbf{u} = 0, \qquad 2a. $$
$$ \frac{D\mathbf{u}}{Dt} + 2\Omega_0 \times \mathbf{u} = -\frac{1}{\rho_0}\nabla p' - \alpha T' \mathbf{g} + \beta S' \mathbf{g} + \nu \nabla^2 \mathbf{u}, \qquad 2b. $$

$$ \frac{\partial S'}{\partial t} + (\mathbf{u} \cdot \nabla)S' = Q_1, \qquad 2c. $$

$$ \frac{\partial T'}{\partial t} + (\mathbf{u} \cdot \nabla)T' = \kappa \nabla^2 T', \qquad 2d. $$

with thermal diffusivity $\kappa = k/c_p\rho$ and thermal conductivity $k$.

For mechanically-driven flows, it is normally assumed that processes are isohaline (v) and isentropic (vi). Doing so and neglecting small terms, the non-magnetic (i) Equation 1 reduces to (Gerkema & Zimmerman 2008):

$$ \nabla \cdot \mathbf{u} = 0, \qquad 3a. $$

$$ \frac{D\mathbf{u}}{Dt} + 2\Omega \times \mathbf{u} = -\frac{1}{\rho_0}\nabla p' + \frac{\rho'}{\rho_0}\mathbf{g} - \nabla\phi' + \nu\nabla^2\mathbf{u} - \frac{d\Omega}{dt}\times\mathbf{r}, \qquad 3b. $$

$$ \frac{g}{\rho_0}\frac{D\rho'}{Dt} - N^2 u_r = 0, \qquad 3c. $$

$$ \nabla^2 \phi' = 4\pi G \rho', \qquad 3d. $$

where $u_r$ is the radial velocity and $N^2 = -(g/\rho_0)(d\rho_0/dr + \rho_0 g/c_s^2)$ is the square of the buoyancy frequency.

Equation 3 still contains rich and complex dynamics. The equations are often linearized (xii) to study ocean waves, which removes different phenomena that can cause turbulence. The traditional approximation (TA) (xiii) and shallow water approximation (xiv) are also often employed. The TA consists of neglecting the terms of the Coriolis force that arise from the non-radial component of the rotation vector to separate vertical and horizontal dynamics. In the shallow-water approximation, the inertial term in the radial momentum equation is ignored, transforming it into a force balance between the pressure gradient and the gravitational and Coriolis accelerations; or when used in combination with (xiii) by simply hydrostatic balance. (xiii) and (xiv) are typically invoked when the horizontal lengthscale of ocean motions is much larger than the ocean thickness; however, their validity involves more subtle aspects (e.g., Gerkema et al. 2008). If the ocean is unstratified (xv), these assumptions lead to the classic Laplace Tidal Equations (LTE):

LTE: Laplace Tidal Equations

$$ \partial_t(\eta_s - \eta_b) + \nabla \cdot (H\bar{\mathbf{u}}) = 0, \qquad 4a. $$

$$ \partial_t \bar{\mathbf{u}} - f \bar{\mathbf{u}} \times \mathbf{e}_r = -\frac{1}{\rho_0} \nabla p' - \nabla \phi' + D(H, \bar{\mathbf{u}}), \qquad 4b. $$

$$ p' = g\rho_0\eta_s + q, \qquad 4c. $$

where $\eta_s$ and $\eta_b$ are radial displacements of the ocean surface and bottom, $\bar{\mathbf{u}}$ is average horizontal velocity, $f = 2\Omega \cos \theta_c$ is the Coriolis parameter with $\theta_c$ being co-latitud, and $q$ is pressure exerted by the ice shell. $D(H, \bar{\mathbf{u}})$ parameterizes dissipation within the system, and the form of this function is uncertain. For Earth, $D(H, \bar{\mathbf{u}})$ can be tuned to fit a large
set of observations (e.g., Green & Nycander 2013). For icy moons, we do not have the benefit of such observations. In analogy with Earth, it can be assumed that energy is dissipated at the ocean boundaries due to bottom-drag, $D = c_D |\bar{\mathbf{u}}| \bar{\mathbf{u}} / H$, with $H$ the ocean thickness and $c_D$ the bottom drag coefficient (e.g., Hay & Matsuyama 2017). However, the dissipation mechanisms at play in subsurface oceans might be different than on Earth's. Because of this, many studies have considered simpler linear dissipation terms such as $\alpha_R \bar{\mathbf{u}}$, with $\alpha_R$ being the Rayleigh coefficient (e.g., Tyler 2011b), or $\nu_{turb} \Delta \bar{\mathbf{u}}$, with $\nu_{turb}$ being the horizontal eddy viscosity (Chen et al. 2014), and varied the dissipation parameters over orders of magnitude.

If the forcing frequency ($\omega$) is much higher than the rotational frequency of the body ($\Omega$), the rotation of the moon can be neglected (xvi). Finally, if the surface gravity wave speed is much larger than the speed at which the perturbation propagates, mechanically-driven ocean dynamics can be ignored all-together and the ocean follows the equilibrium tide (xvii), given by $p'/\rho_0 + \phi' = 0$.

# 4. BUOYANCY-DRIVEN FLOWS

Motivated by regions of disrupted ice on Europa's surface, known as chaos terrains, early work on Europa considered whether hydrothermal venting on the seafloor could lead to melt-through events at the satellite's surface (Thomson & Delaney 2001, Goodman et al. 2004, Goodman & Lenferink 2012). These plumes were envisioned to traverse the ocean in narrow columns confined by the Coriolis force, a consequence of the Taylor-Proudman theorem where rapid rotation limits shearing in the direction parallel to the rotation axis. While melt-through events are not favored energetically (e.g., Lowell & DuBose 2005, Vance & Goodman 2009), the feasibility of columnar plumes extending across the ocean remains debated. Central to this debate are the degree to which rotation organizes the flow into vertical convection columns and how ocean salinity may limit their ascent, as will be discussed below.

Not limited to Europa, the question of hydrothermal vents and their oceanic manifestation is also of fundamental importance for Enceladus given the nature of ejected materials from its south pole (e.g., Hsu et al. 2015, Kang et al. 2022a, Schoenfeld et al. 2023). Liquid water fluxes from the high-pressure ice layers of Titan and Ganymede have also been hypothesized (Kalousova & Sotin 2020, Kalousová et al. 2023), making plume dynamics important across known ocean worlds. Local convective dynamics are not isolated, however, and can have global consequences through the development of zonal (east-west) flows, meridional overturning circulations, and modulation of heat and salt flux patterns at the ice-ocean interface.

## 4.1. The Role of Rotation

In the limit of rapid rotation where the leading order force balance is between the Coriolis and pressure gradient terms, quasi-geostrophic turbulence associated with vertically stiff convection will lead to an inverse energy cascade and alternating zonal jets with prograde flow at the equator and with widths outside the equatorial region set by the Rhines scale (e.g., Heimpel et al. 2005, Ashkenazy & Tziperman 2021, Bire et al. 2022, Hay et al. 2023). As the influence of rotation decreases, convective turbulence becomes increasingly three-dimensional, which changes the mechanism of zonal wind generation. This change
<figure><img src="image_4.png" /><figcaption>Figure 3: Simulated configurations of Europa's ocean from two proposed models. (a,b) Mean superadiabatic temperature normalized by the vertical contrast across the ocean (panel a) and mean zonal velocity as a function of radius and latitude (panel b) from a direct numerical simulation. (c,d) Mean temperature (panel c) and mean zonal velocity (panel d) from a large-eddy simulation. For both models, red (blue) denotes warm (cool) temperatures and prograde (retrograde) zonal flows. Velocities are given in dimensional units above the color bar and in dimensionless Rossby number (Ro) units below. The comparison between temperature plots should focus on their respective patterns rather than quantitative values, which should not be compared one-to-one between the models. Panels a and b adapted from Soderlund (2019) (CC BY 4.0); panels c and d adapted with permission from Bire et al. (2022) (CC BY-NC-ND 4.0).</figcaption></figure>

is typically associated with a relative increase in the influence of the inertial force, which enhances mixing in both the temperature and angular momentum fields. The resulting angular momentum homogenization leads to retrograde jets at large cylindrical radii and prograde jets closer to the rotation axis (e.g., Aurnou et al. 2007, Gastine et al. 2013). As demonstrated systematically by Soderlund (2019), among others, intermediate regimes also exist with multiple jets and retrograde equatorial flow. When rotation is no longer first order such that convection is largely isotropic, zonal flow generation can become insignificant (Gastine et al. 2013). As highlighted in **Figure 3**, significant debate exists in the literature about which convective regime(s) are most relevant to icy ocean worlds (Soderlund et al. 2014, Soderlund 2019, Amit et al. 2020, Ashkenazy & Tziperman 2021, Kvorka & Čadek 2022, Bire et al. 2022, D.G. Lemasquerier, C. Bierson, & K.M. Soderlund, manuscript in review). This debate stems predominantly from the inability to resolve turbulence across the relevant length and temporal scales, uncertainty in which scaling laws are most appropriate, and because bulk dynamic and heat transfer behavioral transitions do not necessarily coincide (e.g., Cheng et al. 2018, Aurnou et al. 2020, Bire et al. 2022, Hawkins et al. 2023).

An array of latitudinal heat transfer configurations is also found as a function of buoyancy forcing, rotation rate, and fluid properties, characterized by the dimensionless Rayleigh Ra, Ekman E, and Prandtl Pr numbers. Kvorka & Čadek (2022) identify three regimes of cooling behaviors in a large suite of numerical thermal convection models that span at
**Rayleigh number,** *Ra*: ratio of buoyancy to diffusion

**Ekman number,** *E*: ratio of viscous to Coriolis forces

**Prandtl number,** *Pr*: ratio of momentum to thermal diffusivities

**Rossby number,** *Ro*: ratio of inertia to Coriolis forces

**Convective Rossby number,** *Ro*<sub>c</sub>: ratio of thermal wind velocity to the free-fall velocity

**Natural Rossby number,** *Ro*<sup>*</sup>: ratio of rotational timescale to the buoyancy timescale

least one order of magnitude for each *Ra*, *E*, and *Pr* as well as both no-slip and free-slip mechanical boundary conditions, encompassing the models of Soderlund (2019) and Amit et al. (2020). “Polar cooling” where heat flux at the upper boundary peaks at high latitudes is limited to a relatively narrow band of parameter space where the modified transitional number falls in the range $Ra_G^* = RaE^{12/7}Pr^{-1} \in \langle 1, 10 \rangle$, with “equatorial cooling” outside this range. It is worth noting that these transitions exhibit some sensitivity to mechanical boundary conditions and do not correspond with the commonly used convective Rossby number (e.g., Aurnou et al. 2020), which is closely related to the natural Rossby number (e.g., Bire et al. 2022). The thermal convection simulations by Bire et al. (2022) additionally indicate that the relative thickness of the ocean is also significant for the characteristics of heat transfer.

Heat flux into the ocean may also not be spatially uniform (e.g., Běhounková et al. 2021, Choblet et al. 2017a,b). This has been investigated for global-scale variations associated with tidal heating (Lemasquerier et al. 2022, D.G. Lemasquerier, C. Bierson, & K.M. Soderlund, manuscript in review) as well as smaller spatial scales associated with convection in the underlying region (rocky core or high-pressure ice layer) (Kang et al. 2022a, Terra-Nova et al. 2023). It appears that global-scale basal heat flux anomalies can be translated efficiently to the upper boundary, although the length scale where this translation becomes negligible remains unclear.

## 4.2. The Role of Ocean Salinity

The thermal properties of water affect buoyant flows by changing its density, but ocean salinity can have an even stronger control on ocean dynamics. The thermal expansivity of water is negative at low pressures and low salinity, such that the point of maximum density can occur at a temperature warmer than the freezing point. Although the density difference is small in these limited conditions, warming of the fluid counter-intuitively increases its density and cooling decreases its density, meaning cold water rises (Townsend 1964).

If the ocean salinity is low, at pressures relevant to part of Europa's hydrosphere and the entirety of Enceladus and other small moons, this anomalous density contrast driven by thermal expansion is more significant than for deeper oceans (Vance & Goodman 2009, Lawrence 2022, J.D. Lawrence, B.E. Schmidt, P.M. Washam, J.J. Buffo, C Chivers, & S.M. Miller, manuscript in review; hereafter Lawrence et al. 2023a). From this, Melosh et al. (2004) suggested that a buoyant “stratosphere” could develop in the ocean where upwelling motions are inhibited from the ice-ocean interface. Depending on the degree of underlying turbulence, this layer may be eroded away, which depends on the competition between the erosion rate (a mechanical forcing) and the cooling rate of the stable layer that restores its buoyancy. In addition, ice formation at the interface driven by this cold ventilation could create brine rejection that would further erode the layer or drive mixing that could disrupt the stratification.

While originally proposed for Europa, this phenomena is especially pertinent to Enceladus (e.g., Soderlund 2019, Zeng & Jansen 2021, Kang et al. 2022b, Lawrence et al. 2023a). If a stable layer were to persist, it would suppress heat transfer to the ice shell and delay the transport of tracers that would need to be diffused rather than advected. The detection of nanosilicate particles in Enceladus' plume suggests that ice-ocean-rock exchange occurs and, further, their size implies a short transport timescale on the order of months to years (Hsu et al. 2015), which is supported by convective particle entrainment calculations
(Schoenfeld et al. 2023).

However, this density change due to thermal expansivity is also pressure dependent, such that the mass of the planet and radius of the ice shell affect the sign of $\alpha$. The role of pressure is important as it changes the point at which the ice shell and ocean can exchange mass through melting via freezing point depression, or refreezing through supercooling. The adiabatic lapse rate determines the change in temperature that is incurred by a parcel of water moving from one depth to another, in the absence of heat or material exchange with the surrounding environment, and scales with gravity. For ocean worlds, this means that the freezing point depression is more significant for larger planets and thicker ice, i.e. a parcel of water traversing a 1 m change in depth on Earth would experience the same change as over 8.5 m on Europa, and ~85 m on Enceladus (Lawrence 2022, Lawrence et al. 2023a). All things being equal, this means that more significant topography or thickness variations can be maintained by the ice shells of smaller ocean worlds than on Europa, Ganymede, or Titan, which is consistent with the variability of Enceladus' ice shell (e.g., Čadek et al. 2019). It also means that the regime of ice-ocean interactions may change over time as the ice shell thickness changes (Lawrence 2022, Lawrence et al. 2023a). Kang & Jansen (2022) arrive at similar conclusions using ocean circulation models. They find that large icy moons with strong gravity tend to have stronger ocean heat transport under the same ice shell topography, leading to an equilibrium ice shell geometry that is flatter on moons with larger size than those that are smaller, and vice versa.

On Earth, freezing point suppression under ice shelves creates an exchange of material through both melting and ice production. Under roughly equilibrium conditions, such as those found beneath the Ross Ice Shelf, Antarctica, where the ocean temperatures are always near the freezing point and the regional circulation is heavily dependent on thermo-haline processes, melting of deeper ice is balanced by freezing at shallower depths, via marine ice accretion (direct freezing onto the ice shelf base), or frazil and platelet ice production (free-floating ice crystals in the water column formed from relief of supercooling) (e.g., Lewis & Perkin 1986, Wolfenbarger et al. 2022). This process is called the “ice pump” and ice formation acts to erode topographic differences across the ice shelf, filling in crevasses (e.g., Lawrence et al. 2023b, Lawrence 2022), affecting the base of thinner regions of the ice shelf hundreds of kilometers away (Koch et al. 2015), and influencing sea ice formation and the sub-ice platelet ice layer (Robinson et al. 2014). The ratio $\alpha\Delta T/\beta\Delta S$ describes whether an environment produces melting or refreezing (Lawrence 2022, Lawrence et al. 2023a), and ice pumps can form when $\alpha\Delta T/\beta\Delta S < 1$ as occurs under ice shelves (Figure 4). On the other hand, for low salinity but high pressures, a narrow regime in which $\alpha\Delta T/\beta\Delta S > 1$ exists where salinity does not control buoyancy. Figure 4 shows that all but the very freshest oceans and deepest ice shell depths for Ganymede and Titan would permit an ice pump to form.

Some recent progress has been made in capturing such ice ocean interactions into planetary ocean models, but none are complete. Ashkenazy & Tziperman (2021) used large-eddy simulations of Europa's ocean and included the effects of melting and freezing of the ice shell, finding Taylor columns, alternating zonal jets with prograde flow at low latitudes, colder equatorial waters, and a nearly uniform ice shell thickness to be driven by a $\alpha\Delta T/\beta\Delta S$ ratio much less than unity. For Enceladus, in steady-state, large variations in the ice shell would be required to drive melting or refreezing. While these variations must be maintained against viscous relaxation, ice flow is also likely slow due to low gravity and cold temperatures through much of the ice. Since both thermal and compositional buoyancy
<figure><img src="image_5.png" /></figure>

Figure 4: Three ice-ocean interaction regimes, determined by the temperature of maximum density and $\alpha\Delta T/\beta\Delta S$ upon meltwater mixing. In the lacustrine regime, $\alpha$ is negative and both $\Delta T$ and $\Delta S$ upon melting contribute to positive buoyancy and ice pumping occurs. In the marine regime, $\alpha$ is positive, but $\beta$ dominates and the freshening effect upon melting causes ice pumping. The white shaded rectangle indicates observed pressure and salinity ranges for sub-ice shelf oceans in the marine regime (Ross and Amery ice shelves; RIS and AIS). At higher pressures and lower salinities where $\alpha$ is positive, however, the ocean is too fresh for the dilution from melt input to overcome the temperature change, and no ice pumping occurs (e.g., subglacial Lake Vostok beneath the East Antarctic Ice Sheet). Ranges for ice/ocean interface pressures and ocean salinities are approximated with shaded regions for Enceladus, Europa, Ganymede, and Titan. Adapted from Lawrence (2022).

would aid in the production of a freshened upper layer, the resulting density structure in the shallow ocean leads to a meridional overturning circulation with a lens of relatively fresh water below the thinnest ice at the south pole (Lobo et al. 2021). The structure of this circulation depends on the magnitude and meridional extent of stable stratification in the ocean (e.g., Kang et al. 2022b).

While the above studies have considered both compositional and thermal buoyancy anomalies, they have not taken their different diffusion rates into account (i.e., double-diffusive convection). If the temperature gradient is unstable but the salinity gradient is stabilizing, a “diffusive” regime develops that can form staircased layers. Conversely, if the temperature gradient is stable but the salinity gradient in destabilizing, “finger” plumes may develop. Wong et al. (2022) demonstrate that layering by double-diffusive convection is possible in Europa’s ocean, although this may be transient depending on the amplitudes of the buoyancy anomalies. Vertical ocean transport efficiency would be reduced while such layers are present.
# 5. MECHANICALLY-DRIVEN FLOWS

Tides raised by the gravitational pull of the planet and other moons in the system can significantly deform icy moons, producing tidal torques that change the spin of the moon. Precisely measuring the spin-state and tidal deformations of a moon allows evaluation of its internal structure. While early studies assumed oceans to be an inviscid, static layer that transfers pressure between the crust and deep interior (e.g., Van Hoolst et al. 2008), tides and changes in the spin might generate intense flows that can significantly alter the tidal and rotational response of the moon and contribute to the heat budget.

Although tides and spin changes are intimately coupled, flows driven by tides and rotation have largely been studied separately. For the former, changes in the spin state from tidal torques are ignored; for the latter, the tidal bulge is assumed to be stationary and the geometry of the ocean to remain unchanged (rigid shell and mantle/core).

## 5.1. Tides

5.1.1. **The Tidal Potential.** The tides experienced by a moon depend on its orbit and spin. The tidal potential can be expanded as a Fourier series and in terms of spherical harmonics with each component characterised by its frequency $\omega$ and spherical harmonic degree $n$ (total wavenumber) and order $m$ (zonal wavenumber). The amplitude of tidal components quickly decreases with increasing degree, hence it is common to consider only degree 2 terms. The leading term in the tidal potential is a degree 2 order 2 component and has an amplitude $\propto \Omega^2 R^2$, where $\Omega$ is the moon's rotation rate. For synchronous bodies, we expect a static tide that produces a permanent tidal bulge and does not directly drive ocean currents. The remaining terms are the result of the moon's non-zero eccentricity $e$ and obliquity $\theta_O$ (i.e., the inclination of the moon's rotational axis with respect to the orbital plane). The leading terms have a frequency $\Omega$ and amplitudes $\propto \Omega^2 R^2 e$ and $\propto \Omega^2 R^2 \sin\theta_O$, known as the *eccentricity* and *obliquity* tides, respectively. The eccentricity tide contains $m=0$ and $m=2$ components arising from the change in moon-planet distance and from the longitudinal variation of the sub-planet point throughout one orbit, respectively. Obliquity tides only have an $m=1$ component which follows from the latitudinal libration of the sub-planet point. Both eccentricity and obliquity tides can be split into westward and eastward propagating components. The next-in-order terms are often ignored given the small eccentricity and obliquity of icy moons (**Supplementary Tables 1-2**) — a complete expansion of the tidal potential for planet tides is given in Renaud et al. (2021). Other moons in the system also raise tides (Hay et al. 2020). As the moon's mass is much smaller than that of the planet, the amplitude of moon-moon tides is lower than planet tides. However, moon-moon tides have a richer frequency spectrum and introduce additional notes into the tidal response.

5.1.2. **The Laplace Tidal Equations.** Most works on icy moon ocean dynamics have relied on the LTE (Equation 4), and the ocean's response can be understood in terms of eigenmodes. The excitation of the modes depends on the spatial structure and frequency of the tidal force and a series of non-dimensional parameters (e.g., Rovira-Navarro et al. 2023). For a non-self-gravitating surface ocean overlying a rigid mantle/core, the tidal response is
**Lamb Parameter, c:**

ratio between the squares of the speeds at which the tidal perturbation and surface gravity waves propagate

controlled by the Lamb parameter: $\epsilon = \Omega^2 R^2 / gH$. For small $\epsilon$, two types of modes can be distinguished: class 1 and class 2 modes with eigenfrequencies given respectively by $\omega^2 / \Omega^2 = n(n+1) / \epsilon$ and $\omega / \Omega = 2m / (n(n+1))$. Class 1 modes are surface gravity waves modified by rotation, gravity acting as the restoring force; class 2 modes, also known as Rossby waves, arise due to the variations of the Coriolis force with latitude.

For icy moons $\epsilon \ll 1$, which implies that resonant surface gravity waves are excited if $\omega^2 / \Omega^2 \gg 1$. This is not the case for planet-tides, which have a characteristic frequency of $\Omega$. This implies that resonant gravity modes would only be excited if subsurface oceans were much thinner ($\lesssim 1$ km). On the other hand, the frequency of some of the components of moon-moon tides is much higher. Because of this, moon-moon tides can resonantly excite the ocean, leading to energy dissipation that can exceed radiogenic heating and other tidal sources (Hay et al. 2020) and have distinct observational signatures (Hay et al. 2022). Away from resonances, the ocean response is well-approximated by the equilibrium tide; the ocean adjusts almost immediately to changes in the geoid ($\eta_s = -\phi'/g$), and both characteristic flow velocities and tidal heating are small. Relaxing the mass-less and rigid mantle/core assumptions (see **Table 1**) modifies the resonant condition and increases the amplitude of the equilibrium tide, which can be now written in terms of a degree-dependent admittance $Z_n$, $\eta_s = -Z_n\phi'/g$ (Matsuyama 2014, Matsuyama et al. 2018). An ice shell does not suppress surface gravity waves. The effect of the overlying ice shell depends on the moon's effective rigidity, $\gamma = 2\mu_{\text{ice}}(1 + \nu_{\text{ice}})H_{\text{ice}}/g\rho_{\text{ice}}R^2$, with $\nu_{\text{ice}}$ the Poisson's ratio of ice (Matsuyama et al. 2018, Beuthe et al. 2016). For moons with a low effective rigidity, $\gamma \lesssim 1$ (*soft shell moons*), the overlying ice shell has a small effect on the ocean's tidal response, while the opposite applies to moons with high effective rigidity (*hard shell moons*). In terms of gravity wave resonances, the effective rigidity increases the eigenfrequencies of surface gravity waves, reducing the strength of planet-tides (**Figure 5**). The ice shell also decreases the admittance $Z_n$ so that $\gamma \to \infty$, $Z_n \to 0$.

The westward component of the obliquity tide can excite Rossby waves. For an ice-free, inviscid ocean, Rossby waves are divergence-free ($\eta_s = 0$) and their amplitude is independent of ocean thickness and proportional to the moon's obliquity (Tyler 2008). Rossby waves can drive strong ocean currents provided the moon's obliquity is high (e.g., $0.1 \text{ m s}^{-1}$ for Europa provided its obliquity is 0.1 deg), and lead to energy dissipation higher than dissipation in the moon's solid layers (**Figure 5**). Except for Titan, the obliquity of icy moons has not been measured. Obliquity can be estimated assuming that the moon is in a Cassini state (Chen et al. 2014). From these estimations, it has been suggested that Rossby waves are relevant for Europa's (Tyler 2014) and Triton's (Nimmo & Spencer 2015) subsurface oceans and in the orbital evolution of Callisto (Downey et al. 2020). However, the above conclusions apply to oceans of constant thickness. Ocean thickness variations alter the eigenfrequency of Rossby waves and hinder their excitation by obliquity tide (Rovira-Navarro et al. 2023).

**Cassini State:**

Rotational state in which (1) rotation and orbital periods are equal; (2) the rotation axis has a constant inclination with respect to the ecliptic; and (3) the three axes related to spin, orbit and precession are co-planar

**5.1.3. Beyond the Laplace Tidal Equations.** LTE use is often justified by the small value of the ratio between ocean thickness and the characteristic length-scale of the tidal force (the body's radius) and the fact that the ocean is experiencing vigorous convection that precludes the formation of a density gradient. However, the ratio between ocean thickness and radius is as high as $\sim 2 \cdot 10^{-1}$ for some icy moons, subsurface oceans could be partially stratified (see Section 4), and flow instabilities might be caused by tides.

Wave-motions beyond the LTE can be studied using a linearized version of Equation 4. The stratification of the ocean in density and angular momentum makes it possible for
<figure><img src="image_6.png" /></figure>

Figure 5: Non-dimensional tidal heating due to the eccentricity (left) and obliquity (right) tide. Dimensional tidal heating for some icy moons is indicated. In all cases, self-gravity is ignored and a rigid core assumed. Tidal heating depends on the linear-drag coefficient; here, $\alpha_R/\Omega = 10^{-3}$ is employed. Regions of intense energy dissipation (yellow regions) are associated with strong ocean currents. Adapted from Rovira-Navarro et al. (2023).

internal gravito-inertial to be excited. The properties of these waves are markedly different from those of surface gravity and Rossby waves (Maas 2005). If the ocean is unstratified in density, only internal inertial waves can be excited. Internal inertial waves propagate at a constant angle with respect to the rotational axis, $\arcsin(\omega/2\Omega)$, which for diurnal tides equals $30^\circ$. Upon reflection, this angle remains constant but the wavelength might change. Depending on the ocean's geometry, this can lead to wave-focusing; waves converge towards some particular trajectories known as wave-attractors (Rieutord et al. 2001). The shape of the attractor depends on the forcing frequency and the radius ratio between the ocean's bottom and top. For an inviscid ocean, a singularity arises along wave-attractors in which velocity grows unbounded. An additional singularity also develops at the critical latitude (Stewartson & Rickard 1969). When the non-zero viscosity of the ocean is considered, these singularities give rise to internal shear layers where substantial energy can be dissipated (Rieutord et al. 2001). Rovira-Navarro et al. (2019) and Rekier et al. (2019) showed that the flow fields of Europa and Enceladus may be dominated by rays emanating from the critical latitude that might feed an attractor. The amount of energy dissipation within the shear layers varies by orders of magnitude depending on the ocean geometry, but remains well below the amount of radiogenic and tidal heat generated in the solid layers.

If the ocean is also stratified in density, the full complexity of gravito-inertial waves comes into focus. The modes excited by the tidal force depend on the stratification structure and the forcing and rotational frequencies. Two types of modes can be distinguished: short-wavelength modes that feature prominent internal shear layers, and long wavelength modes that do not feature internal shear layers (Dintrás et al. 1999). The first set of modes are similar to those discussed for internal inertial waves. However, while the wave-attractors characteristic of internal inertial waves span the whole domain, for gravito-inertial modes, wave-attractors are curbed and only fill parts of the fluid domain. Further insight can be obtained using the shallow-water approximation, which transforms Equation 3 into a LTE-type equation with an equivalent depth that depends on the ocean's stratification profile and is smaller than the total ocean depth (e.g., LeBlond & Mysak 1978). The eigenfrequencies of baroclinic modes are lower than that of surface gravity waves, allowing internal wave resonances to be excited by planet tides (Tyler 2011a). Using a simplified two-layer model,

**Critical Latitude:**

Ocean bottom latitude where the excited internal wave ray is tangent to the bottom ocean boundary
Rovira-Navarro et al. (2023) showed that internal gravity waves are more relevant for moons with high effective rigidity, such as Enceladus. Here, internal wave resonances can drive strong ocean currents (~1 m s⁻¹) and tidal heating can be as high as the moon’s observed thermal output. Such strong currents challenge the linearization assumption. The velocity shear associated with the baroclinic mode further makes the system prone to the Kelvin-Helmholtz instability, which can lead to ocean mixing and hence result in complex feedbacks between buoyancy and mechanically-driven flows.

Note that in the discussion so far, internal waves are excited at the global scale by tidal forcing. A local excitation is also possible, due for example to the barotropic tide flowing across boundary topography. Such a mechanism is key to closing the Earth’s oceanic energy budget; it accounts for approximately 30% of the total energy dissipated by ocean tides (e.g., Garrett & Kunze 2007, Xie et al. 2023). No study relevant for subsurface oceans, involving different shell thicknesses and topographies along both top and bottom boundaries as well as various relative amplitudes of stratification and rotation, have been performed so far. Focusing on the simple spin-up of a constant density fluid, the seminal experiment of Burmann & Noir (2018), however, highlights the primordial role of inertial waves radiated from boundary roughness.

Reintroducing non-linearities in the momentum equation leads to even more complicated and poorly explored dynamics. Non-linear self-interaction of a tidally excited internal wave can drive zonal flows (Tilgner 2007, Morize et al. 2010), whose associated shear might even become large enough to generate local turbulence (Favier et al. 2014, Sauret et al. 2014). Wave focalisation around attractors or because of their interactions with an ambient mean flow might also lead to a strong increase in their amplitude and to breaking (e.g., Staquet & Sommeria 2002, for internal gravity waves). Finally, the tidal flow can trigger a so-called parametric instability involving two gravito-inertial waves, complementary in terms of wavenumbers and frequencies (e.g., Le Bars et al. 2015). Even if the tidal excitation is of small amplitude, the resonant feature of this tidal instability might result in intense, three-dimensional flows whose saturation is still debated (Le Reun et al. 2017, 2018). For now, these various nonlinear tidal processes have been studied either in local theoretical models ignoring boundaries (e.g., Barker & Lithwick 2013), or in global experimental and numerical models in idealized geometries, including a few in deep spherical shells (e.g., Aldridge et al. 1997, Lacaze et al. 2005). Their assessment in the thin spherical shell geometry relevant for icy satellite oceans remains a challenging prospect, mainly because the overwhelming viscous dissipation in accessible laboratory and computer regimes renders the fully turbulent state inaccessible. Also, while these tidal flows still operate in the presence of convection (Cébron et al. 2010, de Vries et al. 2023), the interactions of these possibly simultaneous processes deserve additional studies.

## 5.2. Spin-driven Flows

5.2.1. Libration. The tidal bulge supports a periodic torque along the elliptical orbit, tending to re-align it with the central planet. This produces a small oscillation of the moon's spin called physical longitudinal libration. For synchronously rotating moons, the libration frequency equals the orbital and the mean rotational ones; but a rich libration spectrum settles in from gravitational couplings with other neighboring bodies (e.g., Richard et al. 2014). In structural libration models, ocean flow is neglected: the fluid follows a simple solid body rotation independent of the libration of its boundaries. This is indeed the largely dominant
flow. Yet libration might generate intense flows and dissipation in the liquid ocean, hence significantly contributing to its heat budget. As described below and in more detail in Le Bars et al. (2015), there is a large variety of possible flows depending on the moon spin, orbit, mass, etc. These intense flows might be located in the ocean boundary layers, i.e. near the ice-ocean interface or the ocean-mantle interface, or in the ocean bulk. Also, these flows might reach a turbulent quasi-steady state, or involve successive cycles of growth and collapse.

Librating the solid boundary of a rotating fluid induces successive spin-up and spin-down phases through time. Let's first consider an axisymmetric geometry: coupling between solid and liquid is viscous only, and for a rapidly rotating fluid, contained into the thin Ekman boundary layer. Hence, the spin-up and spin-down of the solid boundary are only felt within this Ekman layer, where it propagates by viscosity from the no-slip boundary, while the bulk fluid outside the Ekman layer keeps rotating at its initial spin. The laminar torque associated to this differential rotation has been shown to be the dominant cause of dissipation at the linear order in Enceladus by Rekier et al. (2019), even if it remains largely insufficient to explain the endogenic heat flow measured by Cassini.

These oscillating boundary flows, however, are prone to centrifugal instability, a classical instability occurring in a homogeneous fluid owing to the dynamical effects of rotation and of streamline curvature (e.g. Drazin & Reid 2004); boundary flows can then become turbulent. This mechanism has been studied by Noir et al. (2009), who showed the emergence of longitudinal rolls, then of boundary turbulence, above some given threshold in terms of libration amplitude, defined as the product of the libration frequency and the maximum angular displacement. Enceladus and numerous moons actually fall in the turbulent regime, leading to strongly increased heating: while the exact quantification of this turbulent dissipation is still unknown, Wilson & Kerswell (2018) provide a rigorous upper bound that would fulfill *Cassini* measurements when the effects of tidal distortion and boundary roughness are included.

In addition to boundary layer flows, libration of an axisymmetric container can also resonantly excite internal waves that propagate in the system bulk. This mechanism has been observed in rotating experiments, in both a sphere (Aldridge & Toomre 1969) and in spherical shells (Koch et al. 2013, Hoff et al. 2016, Lemasquerier et al. 2017). Those later works report more intense inertial mode resonances in shells thanks to energy focusing towards wave attractors, as already described for tidal forcing. Yet, the efficiency of the libration resonant excitation has been challenged in the limit of low viscosity and small libration amplitude by Zhang et al. (2013), so the applicability of direct wave forcing in axisymmetric oceans remains open. For instance, Rekier et al. (2019) do not notice any significant resonance in their linear, laminar, axisymmetric libration models. Note, however, that turbulent Ekman layers can efficiently excite waves, even for a driving libration frequency beyond the frequency range of internal wave existence (Sauret et al. 2013).

Libration studies in spherical shells (Koch et al. 2013, Hoff et al. 2016, Lemasquerier et al. 2017) also systematically report the emergence of a prograde, geostrophic zonal jet within the Stewartson layer along the tangent cylinder. Geostrophic jet formation is a generic process due to the non-linear self interaction of any oscillating flow. In addition to internal boundary layer flows that produce the Stewartson jet, non-linear self interaction of oscillatory flows and geostrophic jet formation also take place in the outer boundary layer of librating systems (Busse 2010, Sauret et al. 2010) and in the bulk flow associated to internal waves (Calkins et al. 2010). Torques and dissipation associated to those zonal flows, and

**Ekman layer:**
boundary layer in
which the
momentum equation
experiences balance
between the viscous
force and Coriolis
acceleration
**Stewartson layer:**
detached shear
layers in a rotating
container, whose
general structure
and thickness
depend on the
Ekman number
(Stewartson 1966)

**Tangent cylinder:**
imaginary right
cylinder tangent to
the ocean's bottom
boundary at the
equator

their consequences for the celestial motions of planets and moons, remain to be investigated. Even more interesting for dissipation, these jets become unstable at sufficiently large forcing amplitudes, generating bulk filling turbulence (Hoff et al. 2016).

As for tidal flows, accounting for asymmetry of the ocean leads to even more intense and complex flows. All processes described above persist; resonant wave and zonal flow excitations are amplified by the topographic forcing (e.g., Noir et al. 2012, Rekier et al. 2019). In ellipsoidal shells, libration can also excite a so-called parametric instability, similar to the one described in the previous section for tidal forcing. The basic mechanism for instability is a resonance between two internal waves of the system and the libration base flow in the bulk induced by topographic coupling (see details in Le Bars et al. 2015). This instability exponentially grows and saturates in a complex turbulent state made of a superposition of columnar structures, wavy patterns, and overturning flows (e.g., **Figure 6**). Quantification of turbulent saturation is a forthcoming challenge to fully appreciate the role of libration in subsurface oceans. In fact, this question echoes the broader, longstanding issue of rotating turbulence (Godeferd & Moisy 2015), where two different models, the strong and weak turbulence models, interact and compete (the same is true for stratified flows, but will not be detailed here, see Brouzet et al. 2016). In the strong turbulence model, i.e. for strong forcing, the flow is quasi-geostrophic and concentrates in large columnar structures aligned with the rotation axis following an inverse energy cascade (see also Section 4.1). In the weak turbulence model, i.e. for small forcing, energy transfers through successive triadic internal wave resonances extending the parametric resonance instability described above, leading to a fully three-dimensional state. The predominance of geostrophic vs. wave turbulence in subsurface oceans remains an open question, while both lead to very different estimates in terms of typical flow scales, torques, and dissipation (Le Reun et al. 2019).

Finally, most studies of libration-driven flows have considered an isodensity fluid: studies coupling libration and buoyancy effects are required to better assess their role in subsurface oceans. A priori, all previous results should apply to an adiabatic ocean efficiently mixed by a strong, small-scale convection; they should also apply to a stably stratified ocean, replacing inertial waves by gravito-inertial ones. The other main remaining challenge stands in developing fully coupled models involving both the nonlinear ocean dynamics and the evolution of the overlying ice crust, including its realistic tidal deformation and its local growth or melting following the spatiotemporal heat flux variations associated with oceanic turbulence.

**5.2.2. Precession.** Precession corresponds to the mechanical forcing associated with rota-
tion of the spin axis of a planet or moon along the normal to the ecliptic. The amplitude
of this forcing is equal to the product of the precession rate and the sine of the obliquity.
Precession is assumed subdominant compared to tides and libration for driving present
subsurface ocean flow. But this might have been different in the past, and precession can
combine with libration and tidal instabilities to induce even stronger flows.

Fluid dynamics of precession has been historically studied because of its possible rel-
evance for planetary core flows and dynamos (e.g., Malkus 1968). It is largely similar to
the dynamics already described above for other mechanical forcings, with some additional
complexities and specificities (see the review in Le Bars et al. 2015). Precession first drives
a “Poincaré” base flow, which consists of a solid body rotation at a constant rate lagging
behind the solid boundaries motion. This Poincaré base flow can resonate with the internal
eigenmodes of the fluid, in particular with the simplest one called the spin over mode or Free
<figure><img src="image_7.png" /></figure>

Figure 6: Sketch and visualizations of the turbulent flow in a meridional plane of the UCLA libration set-up by Lemasquerier et al. (2017), for a full ellipsoid and two different ellipsoidal shells (the rigid inner core is highlighted by a yellow dotted curve). Strong vertical shear zones and superposed wavy patterns are visualised in a laser sheet by light reflection from small elongated flakes. These images are obtained for a libration frequency equal to four times the spin rate, but for different libration angles and rotation rates depending on the limits of the experimental set-up. More cases and details are given in Lemasquerier et al. (2017) (see in particular their figure A1). The purpose here is to illustrate that the turbulent state persists no matter what the ocean depth is.

Core Nutation by astronomers, leading to strong differential rotation with the boundaries (Malkus 1968, Nobili et al. 2021). Second, precession excites internal waves which organize in oblique shear layers, whose nonlinear interactions, as well as nonlinear interactions in the boundary layers, drive intense axisymmetric geostrophic jets. Also, the shear between the fluid and the boundaries due to the Poincaré flow can drive turbulence in the Ekman layers, leading to significantly increased dissipation (Cébron et al. 2019). Finally, precession can excite a parametric instability, whose exact forcing mechanism is still controversial (Nobili et al. 2021). In spherical geometries, the forcing comes from the excited oblique shear layers (Lin et al. 2015). But accounting for the polar flattening and the equatorial bulge of real planets and moons, this conical shear instability competes with two other mechanisms described by Kerswell (1993), corresponding to the periodic shear and to the elliptical distortion of the rotating streamlines of the Poincaré flow in non-spherical domains. Note that the different mechanisms are not exclusive, but may be present at the same time and interact with each other. As described for the libration case, the turbulent state of precession flows is still unknown.
## 6. ELECTROMAGNETICALLY-DRIVEN FLOWS

A novel mechanism for driving ocean circulation has been hypothesized for Europa through coupling with Jupiter's magnetosphere (Gissinger & Petitdemange 2019). Here, the tilt of the planet's dipolar magnetic field component causes the satellite to experience a periodic magnetic oscillation that drives fluid motions through a Lorentz force. These flows would be concentrated in a retrograde equatorial jet and associated low-latitude meridional overturning. Ohmic dissipation within the ocean is also expected, although it is likely small compared to radiogenic and tidal heating. The resulting electromagnetic torque applied to the ocean depends on the phase lag between the planet's field and the satellite's induced magnetic field, and the resulting flow speeds are argued to depend on the planet's magnetic field strength at the satellite, so this effect decreases with increasing distance from the planet. Electromagnetic pumping may therefore be significant at Europa, small at Ganymede, and negligible at Callisto. In the Saturnian system, this effect is likely small due to the planet's nearly axisymmetric magnetic field. The large tilts of Uranus' and Neptune's dipolar magnetic field component would lead to a significant magnetic oscillation, but this would be offset by the weaker total field strength; the role of electromagnetic pumping is thus unclear for ice giant systems and warrants further investigation.

## 7. PERSPECTIVES FROM EARTH

Observations from beneath terrestrial ice shelves suggest further lessons for planetary oceanography. In systems that are out of equilibrium, i.e. where some process perturbs the temperatures in the upper ocean from the near steady state they generally experience, progressively fresher boundary layers may be built up over time due to ice melt without reaching the freezing point. This occurs because of the difference between the diffusion rates of heat and salt. However, advection of heat laterally also becomes important in such regimes and can drive widening of ice shell topography, as opposed to deepening or erasing it. This has been observed in rapidly melting regions of Antarctic Ice shelves (e.g., Dutrieux et al. 2014, Schmidt et al. 2023) and Petermann Gletscher, Greenland (e.g., Washam et al. 2018).

Melting along ice surfaces is not uniform, but instead depends on stratification and turbulent heat flow. In fully turbulent conditions, the surface of melting ice shelves develops a characteristic golfball-like texture of “scallops” and other quasi-circular or elongate melt structures (e.g., Lawrence et al. 2023b). These melted features create drag within the boundary layer. Such ice-ocean coupling can either inhibit or enhance boundary layer-driven processes, depending on conditions.

Dynamics of plumes along ice-ocean interfaces are increasingly important for understanding Earth's cryosphere (e.g., Hewitt 2020, Straneo & Cedenese 2015), but thus far have not been studied in depth on planetary bodies. Consider a buoyant plume (say a mixture of ocean water and meltwater sourced from deep draft ice) rising along a sloped ice-ocean interface: as the parcel moves up, increased buoyancy increases its rate of ascent. At the same time, this increased flow speed also drives turbulence in the surrounding ocean, which drives mixing of the ocean water into the parcel. This feedback is key: under the right conditions, the rising plume never actually reaches the ice-ocean interface and instead breaks away from the ice lower in the water column, driving further mixing and disrupting stratification. Other feedbacks between heat transport, ice slope and shape, current speed, and stratification can strongly influence how the ice and ocean exchange material and heat,
affecting ocean dynamics as well as the ice (e.g., Schmidt et al. 2023, Lawrence et al. 2023b). Depending on conditions, ice topography such as keels and crevasses can be erased (through freezing or melting) or enhanced (by melting).

These dynamics suggest complex and important feedbacks between temperature, com-
position, pressure, ice shell thickness, and topography (on many scales) as well as circulation
and mechanical forcing between the ice and the ocean, such that many of the simplifying
assumptions made in studying planetary oceanographic regimes are unable to realistically
constrain these interactions. In particular, the mass flux between the ice and the ocean,
which is thought to be important in determining not only the interior structure of ocean
worlds but also their habitability, may be driven by these interactions. Moreover, stratifi-
cation and topography strongly affect each other, as well as ocean circulation both locally
and on broad scales, such that self-consistently capturing these phenomena are important
future advances for both polar and planetary oceanography.

We further argue that the central iron-alloy core of Earth (and other planets) also pro-
vides a useful analog for planetary oceanography as both are global reservoirs of a rotating,
low-viscosity liquid (see Landeau et al. 2022, for a review of potential processes driving
the geodynamo). Moreover, core convection is driven by both thermal and compositional
buoyancy anomalies, which can operate in tandem or in opposing directions. Mechanically-
driven flows are also at play, especially for the early geodynamo. The modification of flows
by magnetic fields is additionally a principal element. Indeed, many of the models and
laboratory experiments presented above leveraged codes, experiments, and ideas developed
by the deep Earth community. To conclude, the momentum for planetary oceanography
will continue to accelerate in the decades ahead through multi-disciplinary collaborations.

SUMMARY POINTS

1. Oceans beyond Earth are known to exist in at least four ice-covered moons in the outer solar system (Europa, Ganymede, Enceladus, and Titan), with observational hints for several others and theoretical possibilities for more still.

2. These oceans are rich dynamical systems with flows excited and modulated by buoyancy, tides, libration, precession, and electromagnetic pumping. Each driving mechanism is inherently complex, with many permutations themselves, so few studies have yet crossed these artificial boundaries.

3. Ice-ocean interactions, pertaining to the exchange of material, heat, and momentum, are a critical and as yet understudied part of the ocean world system that requires tracking the full array of oceanographic processes.

FUTURE ISSUES

1. Future missions to the outer solar system are in development to explore the habitability of Europa (*Europa Clipper*), Ganymede and the Jovian system more broadly (*JUICE*), and Titan (*Dragonfly*). Flagship-class missions to the Uranian system and Enceladus were also identified as high priorities in the *Origins, Worlds, and Life: Decadal Strategy for Planetary Science and Astrobiology 2023-2032*.

2. As these missions are being designed and launched, it is essential to develop testable
hypotheses for ocean dynamics (e.g., ice shell thickness variations, motional mag-
netic field induction, non-synchronous rotation), noting that the data returned by
these future missions will reflect all driving mechanisms in aggregate.

3. For buoyancy-driven flows, future work should unite direct numerical simulations, large-eddy simulations, and laboratory experiments with theory to understand how convective processes operate at extreme (i.e., more realistic) input parameters.

4. For tidally-driven flows, future work should explore the role of ice shell and basal topography, push the parameter regime of numerical and laboratory experiments towards more realistic regions, and explore feedbacks between ocean tidal dissipation and the interior and orbital evolution of icy worlds.

5. For spin-driven flows, future work should study turbulence in the relevant parameter regime, the effect of boundary roughness, and combine forcings together (e.g., tides in rotating and stratified flows, tides and libration, libration and convection).

6. Work aimed at understanding the interactions and feedbacks between buoyancy-, mechanically-, and electromagnetically- driven flows is necessary to understand subsurface ocean dynamics (e.g., role of ocean mixing in the ocean's overturning circulation).

7. With such critical implications for mass exchange between the ocean and ice, accu-
rately coupled ice-ocean studies is a frontier for both Earth and planetary oceanog-
raphers. Better parameterizations for boundary flows, turbulent mixing, and heat
and salt fluxes between these reservoirs are sorely needed, and realistic treatment
of melting and freezing along dynamic ice shell topography is crucial.

8. Without direct access to subsurface oceans, further exploration of Earth-analog environments and laboratory experiments at subsurface ocean conditions are key to provide required context for the physical oceanography of ice-covered moons.

JUICE: JUpiter ICy
moon Explorer

DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.

ACKNOWLEDGMENTS

KMS and BES were funded by the NASA Astrobiology program grant Oceans Across
Space and Time (Grant No. 80NSSC18K1301). MRN was supported by NASA grant
No. 80NSSC20K0570 issued through the Solar System Workings program. MLB acknowl-
edges the financial support from the French Agence Nationale de Recherche, project ANR
COLOSSE (grant No. ANR-2020-CE49-0010).

LITERATURE CITED

Aldridge K, Seyed-Mahmoud B, Henderson G, van Wijngaarden W. 1997. Elliptical instability of the earth's fluid core. *Phys. Earth Planet. Int.* 103(3-4):365-374

Aldridge KD, Toomre A. 1969. Axisymmetric inertial oscillations of a fluid in a rotating spherical container. *J. Fluid Mech.* 37(2):307-323
Amit H, Choblet G, Tobie G, Terra-Nova F, Čadek O, Bouffard M. 2020. Cooling patterns in rotating thin spherical shells—Application to Titan’s subsurface ocean. *Icarus* 338:113509

Ashkenazy Y, Tziperman E. 2021. Dynamic europa ocean shows transient taylor columns and convection driven by ice melting and salinity. *Nature Comm.* 12(1):1–12

Aurnou JM, Heimpel MH, Wicht J. 2007. The effects of vigorous mixing in a convective model of zonal flow on the Ice Giants. *Icarus* 190:110–126

Aurnou JM, Horn S, Julien K. 2020. Connections between nonrotating, slowly rotating, and rapidly rotating turbulent convection transport scalings. *Phys. Rev. Res.* 2(4):043115

Barker AJ, Lithwick Y. 2013. Non-linear evolution of the tidal elliptical instability in gaseous planets and stars. *Mon. Notices Royal Astron. Soc.* 435(4):3614–3626

Bauer JM, Roush TL, Geballe TR, Meech KJ, Owen TC, et al. 2002. The near infrared spectrum of miranda: Evidence of crystalline water ice. *Icarus* 158(1):178–190

Beddingfield C, Burr D, Emery J. 2015. Fault geometries on uranus’ satellite miranda: Implications for internal structure and heat flow. *Icarus* 247:35–52

Běhounková M, Tobie G, Choblet G, Kervazo M, Melwani Daswani M, et al. 2021. Tidally induced magmatic pulses on the oceanic floor of Jupiter’s moon Europa. *Geophys. Res. Lett.* 48(3):e2020GL090077

Beuthe M, Rivoldini A, Trinh A. 2016. Enceladus’s and Dione’s floating ice shells supported by minimum stress isostasy. *Geophys. Res. Lett.* 43(19):10,088–10,096

Bierson CJ, Nimmo F. 2022. A note on the possibility of subsurface oceans on the uranian satellites. *Icarus* 373:114776

Bire S, Kang W, Ramadhan A, Campin JM, Marshall J. 2022. Exploring ocean circulation on icy moons heated from below. *J. Geophys. Res. Planets* 127(3):e2021JE007025

Bland MT, Showman AP, Tobie G. 2009. The orbital–thermal evolution and global expansion of ganymede. *Icarus* 200(1):207–221

Brouzet C, Ermanyuk EV, Joubaud S, Sibgatullin I, Dauxois T. 2016. Energy cascade in internal-wave attractors. *Europhys. Lett.* 113(4):44001

Burmann F, Noir J. 2018. Effects of bottom topography on the spin-up in a cylinder. *Phys. Fluids* 30(10):106601

Busse F. 2010. Mean zonal flows generated by librations of a rotating spherical cavity. *J. Fluid Mech.* 650:505–512

Čadek O, Souček O, Běhounková M, Choblet G, Tobie G, Hron J. 2019. Long-term stability of Enceladus’ uneven ice shell. *Icarus* 319:476–484

Calkins MA, Noir J, Eldredge JD, Aurnou JM. 2010. Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry. *Phys. Fluids* 22(8):086602

Cartwright RJ, Beddingfield CB, Nordheim TA, Roser J, Grundy WM, et al. 2020. Evidence for ammonia-bearing species on the uranian satellite ariel supports recent geologic activity. *Astro-phys. J. Lett.* 898(1):L22

Castillo-Rogez J, Weiss B, Beddingfield C, Biersteker J, Cartwright R, et al. 2023. Compositions and interior structures of the large moons of Uranus and implications for future spacecraft observations. *J. Geophys. Res. Planets* 128(1):e2022JE007432

Castillo-Rogez JC, Hemingway D, Rhoden A, Tobie G, McKinnon WB. 2018. Origin and evolution of Saturn’s mid-sized moons. In *Enceladus and the Icy Moons of Saturn*, eds. PM Schenk, RN Clark, CJA Howett, AJ Verbiscer, JH Waite. 285

Cébron D, Laguerre R, Noir J, Schaeffer N. 2019. Precessing spherical shells: flows, dissipation, dynamo and the lunar core. *Geophys. J. Int.* 219(Supplement 1):S34–S57

Cébron D, Maubert P, Le Bars M. 2010. Tidal instability in a rotating and differentially heated ellipsoidal shell. *Geophys. J. Int.* 182(3):1311–1318

Chen E, Nimmo F, Glatzmaier G. 2014. Tidal heating in icy satellite oceans. *Icarus* 229:11–30

Cheng JS, Aurnou JM, Julien K, Kunnen RPJ. 2018. A heuristic framework for next-generation models of geostrophic convective turbulence. *Geophys. Astrophys. Fluid Dyn.* 112(4):277–300
Choblet G, Tobie G, Sotin C, Běhounková M, Čadek O, et al. 2017a. Powering prolonged hydrothermal activity inside Enceladus. *Nature Astron.* 1(12):841

Choblet G, Tobie G, Sotin C, Kalousová K, Grasset O. 2017b. Heat transport in the high-pressure ice mantle of large icy moons. *Icarus* 285:252–262

Čuk M, Dones L, Nesvorný D. 2016. Dynamical evidence for a late formation of Saturn's moons. *Astrophys. J.* 820(2):97

Čuk M, Moutamid ME, Tiscareno MS. 2020. Dynamical history of the uranian system. *Planet. Sci. J.* 1(1):22

de Vries NB, Barker AJ, Hollerbach R. 2023. The interactions of the elliptical instability and convection. *Phys. Fluids*

Dintras B, Rieutord M, Valdettaro L. 1999. Gravito-inertial waves in a rotating stratified sphere or spherical shell. *J. Fluid Mech.* 398:271–297

Downey BG, Nimmo F, Matsuyama I. 2020. Inclination damping on Callisto. *Mon. Not. R. Astron. Soc.* 499(1):40–51

Drazin PG, Reid WH. 2004. Hydrodynamic stability. Cambridge university press

Dutrieux P, Stewart C, Jenkins A, Nicholls KW, Corr HFJ, et al. 2014. Basal terraces on melting ice shelves. *Geophys. Res. Lett.* 41(15):5506–5513

Favier B, Barker A, Baruteau C, Ogilvie G. 2014. Non-linear evolution of tidally forced inertial waves in rotating fluid bodies. *Mon. Notices Royal Astron. Soc.* 439(1):845–860

Garrett C, Kunze E. 2007. Internal tide generation in the deep ocean. *Annu. Rev. Fluid Mech.* 39:57–87

Gastine T, Wicht J, Aurnou JM. 2013. Zonal flow regimes in rotating anelastic spherical shells: An application to giant planets. *Icarus* 225:156–172

Gerkema T. 2019. An introduction to tides. Cambridge University Press

Gerkema T, Zimmerman J. 2008. An introduction to internal waves. *Lecture Notes, Royal NIOZ, Texel* 207:207

Gerkema T, Zimmerman JTF, Maas LRM, van Haren H. 2008. Geophysical and astrophysical fluid dynamics beyond the traditional approximation. *Rev. Geophys* 46(2)

Gissinger C, Petitdemange L. 2019. A magnetically driven equatorial jet in Europa's ocean. *Nature Astron.* 1

Godeferd FS, Moisy F. 2015. Structure and dynamics of rotating turbulence: a review of recent experimental and numerical results. *Applied Mechanics Reviews* 67(3)

Goodman JC, Collins GC, Marshall J, Pierrehumbert RT. 2004. Hydrothermal plume dynamics on Europa: Implications for chaos formation. *J. Geophys. Res.* 109:E03008

Goodman JC, Lenferink E. 2012. Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds. *Icarus* 221:970–983

Green JAM, Nycander J. 2013. A comparison of tidal conversion parameterizations for tidal models. *J. Phys. Oceanogr.* 43(1):104–119

Hartkorn O, Saur J. 2017. Induction signals from Callisto's ionosphere and their implications on a possible subsurface ocean. *J. Geophys. Res. Space Phys.* 122(11):11,677–11,697

Hawkins EK, Cheng JS, Abbate JA, Pilegard T, Stellmach S, et al. 2023. Laboratory models of planetary core-style convective turbulence. *Fluids* 8(4):106

Hay HC, Matsuyama I. 2017. Numerically modelling tidal dissipation with bottom drag in the oceans of Titan and Enceladus. *Icarus* 281:342–356

Hay HCFC, Fenty I, Pappalardo RT, Nakayama Y. 2023. Turbulent drag at the ice-ocean interface of Europa in simulations of rotating convection: Implications for nonsynchronous rotation of the ice shell. *J. Geophys. Res. Planets* e2022JE007648

Hay HCFC, Matsuyama I, Pappalardo RT. 2022. The high-frequency tidal response of ocean worlds: Application to europa and ganymede. *J. Geophys. Res. Planets* 127(5):e2021JE007064E2021JE007064 2021JE007064

Hay HCFC, Trinh A, Matsuyama I. 2020. Powering the Galilean satellites with moon-moon tides.
Heimpel MH, Aurnou JM, Wicht J. 2005. Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model. *Nature* 483:193–196

Hemingway DJ, Iess L, Tajeddine R, Tobie G. 2018. The interior of Enceladus. In *Enceladus and the Icy Moons of Saturn*. University of Arizona Press, 57–78

Hewitt IJ. 2020. Subglacial plumes. *Annu. Rev. Fluid Mech.* 52:145–169

Hoff M, Harlander U, Egbers C. 2016. Experimental survey of linear and nonlinear inertial waves and wave instabilities in a spherical shell. *J. Fluid Mech.* 789:589–616

Hofgartner JD, Birch SP, Castillo J, Grundy WM, Hansen CJ, et al. 2022. Hypotheses for Triton's plumes: New analyses and future remote sensing tests. *Icarus* 375:114835

Hsu HW, Postberg F, Sekine Y, Shibuya T, Kempf S, et al. 2015. Ongoing hydrothermal activities within Enceladus. *Nature* 519:207–210

Hussmann H, Spohn T. 2004. Thermal-orbital evolution of Io and Europa. *Icarus* 171(2):391–410

Ingersoll AP. 2005. Boussinesq and anelastic approximations revisited: Potential energy release during thermobaric instability. *J. Phys. Oceanogr.* 35(8):1359–1369

Journaux B, Kalousová K, Sotin C, Tobie G, Vance SD, et al. 2020. High-pressure ices in large ocean worlds. *Space Sci. Rev.* 216(1):7

Kalousová K, Soderlund KM, Solomonidou A, Sotin C. 2023. Structure and evolution of Ganymede's hydrosphere. In *Ganymede*, eds. M Volwerk, MA McGrath, X Jia, T Spohn. Cambridge University Press

Kalousova K, Sotin C. 2020. The Effect of a Clathrate Layer on the Evolution and Dynamics of Titan's Crust. *LPI* (2326):2110

Kang W, Jansen M. 2022. On icy ocean worlds, size controls ice shell geometry. *Astrophys. J.* 935(2):103

Kang W, Marshall J, Mittal T, Bire S. 2022a. Ocean dynamics and tracer transport over the south pole geysers of enceladus. *Mon. Not. R. Astron. Soc.* 517(3):3485–3494

Kang W, Mittal T, Bire S, Campin JM, Marshall J. 2022b. How does salinity shape ocean circulation and ice geometry on enceladus and other icy satellites? *Science Adv.* 8(29):eabm4665

Kerswell R. 1993. The instability of precessing flow. *Geophys. Astro. Fluid Dyn.* 72(1-4):107–144

Khurana KK, Kivelson MG, Stevenson DJ, Schubert G, Russell CT, et al. 1998. Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto. *Nature* 395(6704):777–780

Kivelson M, Khurana K, Volwerk M. 2002. The permanent and inductive magnetic moments of Ganymede. *Icarus* 157(2):507–522

Kivelson MG, Khurana KK, Russell CT, Volwerk M, Walker RJ, Zimmer C. 2000. Galileo magnetometer measurements: A stronger case for a subsurface ocean at Europa. *Science* 289(5483):1340–1343

Koch I, Fitzsimons S, Samyn D, Tison JL. 2015. Marine ice recycling at the southern McMurdo Ice Shelf, Antarctica. *J. Glaciol.* 61(228):689–701

Koch S, Harlander U, Egbers C, Hollerbach R. 2013. Inertial waves in a spherical shell induced by librations of the inner sphere: experimental and numerical results. *Fluid Dynamics Research* 45(3):035504

Kvorka J, Čadek C, Tobie G, Choblet G. 2018. Does Titan's long-wavelength topography contain information about subsurface ocean dynamics? *Icarus* 310:149–164

Kvorka J, Čadek O. 2022. A numerical model of convective heat transfer in Titan's subsurface ocean. *Icarus* 376:114853

Lacaze L, Le Gal P, Le Dizes S. 2005. Elliptical instability of the flow in a rotating shell. *Phys. Earth Planet. Int.* 151(3-4):194–205

Landeau M, Fournier A, Nataf HC, Cébron D, Schaeffer N. 2022. Sustaining earth's magnetic dynamo. *Nature Reviews Earth & Environment* 3(4):255–269

Lawrence JD. 2022. Antarctic ice shelves as ocean world analogs. Ph.D. thesis, Georgia Institute of Technology
Lawrence JD, Washam PM, Stevens C, Hulbe C, Horgan HJ, et al. 2023b. Crevasse refreezing and signatures of retreat observed at Kamb Ice Stream grounding zone. *Nature Geosci.* :1–6

Le Bars M, Cébron D, Le Gal P. 2015. Flows driven by libration, precession, and tides. *Annu. Rev. Fluid Mech.* 47(1):163–193

Le Reun T, Favier B, Barker AJ, Le Bars M. 2017. Inertial wave turbulence driven by elliptical instability. *Phys. Rev. Lett.* 119(3):034502

Le Reun T, Favier B, Le Bars M. 2018. Parametric instability and wave turbulence driven by tidal excitation of internal waves. *J. Fluid Mech.* 840:498–529

Le Reun T, Favier B, Le Bars M. 2019. Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence. *J. Fluid Mech.* 879:296–326

LeBlond PH, Mysak LA. 1978. Waves in the ocean. Amsterdam: Elsevier

Lemasquerier D, Bierson CJ, Soderlund KM. 2022. Can convection in icy moons' oceans translate interior tidal heating patterns to the ice-ocean boundary?, In *Fall Meeting 2022*. AGU

Lemasquerier D, Grannan A, Vidal J, Cébron D, Favier B, et al. 2017. Libration-driven flows in ellipsoidal shells. *J. Geophys. Res. Planets* 122(9):1926–1950

Lewis EL, Perkin RG. 1986. Ice pumps and their rates. *J. Geophys. Res.* 91:11756–11762

Lin Y, Marti P, Noir J. 2015. Shear-driven parametric instability in a precessing sphere. *Phys. Fluids (1994-present)* 27(4):046601

Lobo AH, Thompson AF, Vance SD, Tharimena S. 2021. A pole-to-equator ocean overturning circulation on enceladus. *Nature Geosci.* 14(4):185–189

Lowell RP, DuBose M. 2005. Hydrothermal systems on Europa. *Geophys. Res. Lett.* 32:L05202

Maas LRM. 2005. Wave attractors: linear yet nonlinear. *Int. J. Bifurc. Chaos* 15(09):2757–2782

Malkus W. 1968. Precession of the earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the earth's dynamo. *Science* 160(3825):259–264

Matsuyama I. 2014. Tidal dissipation in the oceans of icy satellites. *Icarus* 242:11–18

Matsuyama I, Beuthe M, Hay HC, Nimmo F, Kamata S. 2018. Ocean tidal heating in icy satellites with solid shells. *Icarus* 312:208–230

McKinnon WB. 1984. On the origin of Triton and Pluto. *Nature* 311(5984):355–358

Melosh HJ, Ekholm AG, Showman AP, Lorenz RD. 2004. The temperature of Europa's subsurface water ocean. *Icarus* 168:498–502

Morize C, Le Bars M, Le Gal P, Tilgner A. 2010. Experimental determination of zonal winds driven by tides. *Phys. Rev. Lett.* 104(21):214501

National Academy of Sciences. 2022. Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032

Neveu M, Rhoden AR. 2019. Evolution of Saturn's mid-sized moons. *Nature Astronomy* 3(6):543–552

Nimmo F, Pappalardo RT. 2016. Ocean worlds in the outer solar system. *J. Geophys. Res. Planets* 121(8):1378–1399

Nimmo F, Spencer J. 2015. Powering Triton's recent geological activity by obliquity tides: Implications for Pluto geology. *Icarus* 246:2–10Special Issue: The Pluto System

Nobili C, Meunier P, Favier B, Le Bars M. 2021. Hysteresis and instabilities in a spheroid in precession near the resonance with the tilt-over mode. *J. Fluid Mech.* 909

Noir J, Cébron D, Le Bars M, Sauret A, Aurnou J. 2012. Experimental study of libration-driven zonal flows in non-axisymmetric containers. *Phys. Earth Planet. Int.* 204:1–10

Noir J, Hemmerlin F, Wicht J, Baca S, Aurnou J. 2009. An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans. *Phys. Earth Planet. Int.* 173(1-2):141–152

Ojakangas GW, Stevenson DJ. 1989. Thermal state of an ice shell on Europa. *Icarus* 81(2):220–241

Pappalardo RT, Belton MJS, Breneman HH, Carr MH, Chapman CR, et al. 1999. Does Europa have a subsurface ocean? Evaluation of the geological evidence. *J. Geophys. Res. Planets*
104(E10):24015-24055

Peale SJ. 1999. Origin and evolution of the natural satellites. *Annu. Rev. Astron. Astrophys.* 37(1):533-602

Peterson G, Nimmo F, Schenk P. 2015. Elastic thickness and heat flux estimates for the uranian satellite ariel. *Icarus* 250:116-122

Porco CC, Helfenstein P, Thomas PC, Ingersoll AP, Wisdom J, et al. 2006. Cassini observes the active South Pole of Enceladus. *Science* 311(5766):1393-1401

Rekier J, Trinh A, Triana SA, Dehant V. 2019. Internal energy dissipation in Enceladus's sub-surface ocean from tides and libration and the role of inertial waves. *J. Geophys. Res. Planets* 124(8):2198-2212

Renaud JP, Henning WG, Saxena P, Neveu M, Bagheri A, et al. 2021. Tidal Dissipation in Dual-body, Highly Eccentric, and Nonsynchronously Rotating Systems: Applications to Pluto-Charon and the Exoplanet TRAPPIST-1e. *Planet. Sci. J.* 2(1):4

Rhoden AR, Walker ME. 2022. The case for an ocean-bearing mimas from tidal heating analysis. *Icarus* 376:114872

Richard A, Rambaux N, Charnay B. 2014. Librational response of a deformed 3-layer titan perturbed by non-keplerian orbit and atmospheric couplings. *Planet Space Sci.* 93:22-34

Rieutord M, Georgeot B, Valdettaro L. 2001. Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum. *J. Fluid Mech.* 435:103-144

Robinson NJ, Williams MJM, Stevens CL, Langhorne PJ, Haskell TG. 2014. Evolution of a super-cooled ice shelf water plume with an actively growing subice platelet matrix. *J. Geophys. Res. Oceans* 119(6):3425-3446

Ross MN, Schubert G. 1990. The coupled orbital and thermal evolution of triton. *Geophys. Res. Lett.* 17(10):1749-1752

Rovira-Navarro M, Matsuyama I, Hay HCFC. 2023. Thin-shell tidal dynamics of ocean worlds. *Planet. Sci. J* 4(2):23

Rovira-Navarro M, Rieutord M, Gerkema T, Maas LR, van der Wal W, Vermeersen B. 2019. Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus? *Icarus* 321:126-140

Sasaki T, Stewart GR, Ida S. 2010. Origin of the Different Architectures of the Jovian and Saturnian Satellite Systems. *Astrophys. J.* 714(2):1052-1064

Saur J, Duling S, Roth L, Jia X, Strobel DF, et al. 2015. The search for a subsurface ocean in Ganymede with Hubble Space Telescope observations of its auroral ovals. *J. Geophys. Res. Space Phys.* 120(3):1715-1737

Sauret A, Cébron D, Le Bars M. 2013. Spontaneous generation of inertial waves from boundary turbulence in a librating sphere. *J. Fluid Mech.* 728

Sauret A, Cébron D, Morize C, Le Bars M. 2010. Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity. *J. Fluid Mech.* 662:260-268

Sauret A, Le Bars M, Le Gal P. 2014. Tide-driven shear instability in planetary liquid cores. *Geophys. Res. Lett.* 41(17):6078-6083

Schenk P. 2010. Atlas of the galilean satellites. Cambridge University Press

Schmidt BE, Washam P, Davis PED, Nicholls KW, Holland DM, et al. 2023. Heterogeneous melting near the Thwaites Glacier grounding line. *Nature* 614(7948):471-478

Schoenfeld AM, Hawkins EK, Soderlund KM, Vance SD, Leonard E, Yin A. 2023. Particle en-trainment and rotating convection in Enceladus' ocean. *Communications Earth & Environment* 4(1):28

Smith BA, Soderblom LA, Banfield D, Barnet C, Basilevsky AT, et al. 1989. Voyager 2 at Neptune: Imaging science results. *Science* 246(4936):1422-1449

Smith BA, Soderblom LA, Beebe R, Bliss D, Boyce JM, et al. 1986. Voyager 2 in the uranian system: Imaging science results. *Science* 233(4759):43-64

Soderlund KM. 2019. Ocean dynamics of outer solar system satellites. *Geophys. Res. Lett.*
46(15):8700-8710

Soderlund KM, Kalousová K, Buffo JJ, Glein CR, Goodman JC, et al. 2020. Ice-Ocean Exchange Processes in the Jovian and Saturnian Satellites. *Space Sci. Rev.* 216(5):1-57

Soderlund KM, Schmidt BE, Wicht J, Blankenship DD. 2014. Ocean-driven heating of Europa's icy shell at low latitudes. *Nature Geosci.* 7:16-19

Sotin C, Kalousová K, Tobie G. 2021. Titan's interior structure and dynamics after the cassini-huygens mission. *Annu. Rev. Earth Planet. Sci.* 49(1):579-607

Staquet C, Sommeria J. 2002. Internal gravity waves: from instabilities to turbulence. *Annu. Rev. Fluid Dyn.* 34(1):559-593

Stewartson K. 1966. On almost rigid rotations. part 2. *J. Fluid Mech.* 26(1):131-144

Stewartson K, Rickard JA. 1969. Pathological oscillations of a rotating fluid. *J. Fluid Mech.* 35(4):759-773

Straneo F, Cedenese C. 2015. The dynamics of greenland's glacial fjords and their role in climate. *Annu. Rev. Mar. Sci.* 7:89-112.

Terra-Nova F, Amit H, Choblet G, Tobie G, Bouffard M, Čadek O. 2023. The influence of heterogeneous seafloor heat flux on the cooling patterns of Ganymede's and Titan's subsurface oceans. *Icarus* 389:1152-323

Thomson RE, Delaney JR. 2001. Evidence for a weakly stratified European ocean sustained by seafloor heat flux. *J. Geophys. Res.* 106:12355-12365

Tilgner A. 2007. Zonal wind driven by inertial modes. *Phys. Rev. Lett.* 99(19):194501

Townsend AA. 1964. Natural convection in water over an ice surface. *Q. J. R. Meteorol. Soc.* 90(385):248-259

Tritton DJ. 1998. Physical fluid dynamics. Oxford: Oxford University Press

Tyler RH. 2008. Strong ocean tidal flow and heating on moons of the outer planets. *Nature* 456(7223):770-772

Tyler RH. 2011a. Magnetic remote sensing of Europa's ocean tides. *Icarus* 211:906-908

Tyler RH. 2011b. Tidal dynamical considerations constrain the state of an ocean on Enceladus. *Icarus* 211(1):770-779

Tyler RH. 2014. Comparative estimates of the heat generated by ocean tides on icy satellites in the outer Solar System. *Icarus* 243(Supplement C):358-385

Van Hoolst T, Rambaux N, Karatekin Ö, Dehant V, Rivoldini A. 2008. The librations, shape, and icy shell of europa. *Icarus* 195(1):386-399

Vance SD, Goodman JC. 2009. Oceanography of an ice-covered moon. In *Europa*, eds. RT Papalardo, WB McKinnon, KK Khurana. Tucson: University of Arizona Press, 459-482

Vance SD, Panning MP, Stahler S, Cammarano F, Bills BG, et al. 2018. Geophysical investigations of habitability in ice-covered ocean worlds. *J. Geophys. Res.* 123:180-205

Vance SD, Styczinski MJ, Bills BG, Cochrane CJ, Soderlund KM, et al. 2020. Magnetic Induction Responses of Jupiter's Ocean Moons Including Effects from Adiabatic Convection. *J. Geophys. Res. Planets* e2020JE006418

Vilella K, Choblet G, Tsao WE, Deschamps F. 2020. Tidally heated convection and the occurrence of melting in icy satellites: Application to Europa. *J. Geophys. Res. Planets* 125(3):e2019JE006248

Waite JH, Glein CR, Perryman RS, Teolis BD, Magee BA, et al. 2017. Cassini finds molecular hydrogen in the Enceladus plume: Evidence for hydrothermal processes. *Science* 356:155-159

Walker ME, Rhoden AR. 2022. Tidal heating at europa using the multifrequency analysis of tidal heating toolkit. *Planet. Sci. J.* 3(7):149

Washam P, Münchow A, Nicholls KW. 2018. A decade of ocean changes impacting the ice shelf of Petermann Gletscher, Greenland. *J. Phys. Oceanogr.* 48(10):2477-2493

Wilson A, Kerswell RR. 2018. Can libration maintain Enceladus's ocean? *Earth Planet. Sci. Lett.* 500:41-46

Wolfenbarger NS, Buffo JJ, Soderlund KM, Blankenship DD. 2022. Ice shell structure and composition of ocean worlds: Insights from accreted ice on Earth. *Astrobiology* 22(8):937-961
Wong T, Hansen U, Wiesehöfer T, McKinnon WB. 2022. Layering by double-diffusive convection in the subsurface oceans of Europa and Enceladus. *J. Geophys. Res. Planets* 127(12):e2022JE007316

Wunsch C, Ferrari R. 2004. Vertical mixing, energy, and the general circulation of the oceans. *Annu. Rev. Fluid Mech.* 36:281–314

Xie X, Liu X, Chen Z, Wang Y, Chen D, et al. 2023. Pure inertial waves radiating from low-frequency flows over large-scale topography. *Geophys. Res. Lett.* 50(4):e2022GL099889

Zannoni M, Hemingway D, Gomez Casajus L, Tortora P. 2020. The gravity field and interior structure of dione. *Icarus* 345:113713

Zeng Y, Jansen MF. 2021. Ocean circulation on Enceladus with a high-versus low-salinity ocean. *Planet. Sci. J.* 2(4):151

Zhang K, Chan KH, Liao X, Aurnou JM. 2013. The non-resonant response of fluid in a rapidly rotating sphere undergoing longitudinal libration. *J. Fluid Mech.* 720:212–235