# Energetics govern ocean circulation on icy ocean worlds.

Malte F. Jansen  Department of the Geophysical Sciences, The University of
Chicago, Chicago, IL  Wanying Kang  Department of Earth, Atmospheric and
Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA
Edwin Kite  Department of the Geophysical Sciences, The University of Chicago,
Chicago, IL Malte F. Jansen, Department of the Geophysical Sciences, The
University of Chicago, Chicago, IL

###### Abstract

Globally ice-covered oceans have been found on multiple moons in the solar
system and may also have been a feature of Earth’s past. However, relatively
little is understood about the dynamics of these ice-covered oceans, which
affect not only the physical environment but also any potential life and its
detectability. A number of studies have simulated the circulation of icy-world
oceans, but have come to widely different conclusions. To better understand
and narrow down these diverging results, we discuss energetic constraints for
the circulation on ice-covered oceans, focusing in particular on Snowball
Earth, Europa, and Enceladus. Energy input that can drive ocean circulation on
ice-covered bodies can be associated with heat and salt fluxes at the
boundaries as well as ocean tides and librations. We show that heating from
the solid core balanced by heat loss through the ice sheet can drive an ocean
circulation, but the resulting flows would be relatively weak and strongly
affected by rotation. Salt fluxes associated with freezing and melting at the
ice sheet boundary are unlikely to energetically drive a circulation, although
they can shape the large-scale circulation when combined with turbulent
mixing. Ocean tides and librations may provide an energy source for such
turbulence, but their strength remains highly uncertain for the icy moons,
which poses a major obstacle to predicting the ocean dynamics of icy worlds
and remains as an important topic for future research.

##  1 Introduction

Globally ice-covered oceans have been found on multiple moons in the solar
system (Carr et al., 1998; Kivelson et al., 2000; Thomas et al., 2016; Nimmo &
Pappalardo, 2016) and spark our curiosity in part due to their potential to
provide hospitable environments for life (Des Marais et al., 2008; Waite et
al., 2017; Postberg et al., 2018; Hendrix et al., 2019). Earth’s oceans may
also have been covered by a global ice sheet during the so-called “Snowball
Earth” events, and indeed eukaryotic life not only appears to have survived
through these episodes, but may have evolved significantly during them (e.g.
Hoffman et al., 2017). However, relatively little is known about these oceans
beyond their existence, and due to our inability to directly observe them at
present, we heavily rely on models to decipher their mysteries.

Although their potential biology may represent the holy grail for research on
the icy moon oceans, it is natural to start with the somewhat more tractable
problem of inferring the physical and chemical environment. In this study we
specifically focus on the ocean circulation and mixing processes, which
control the transport of heat and chemical tracers, including those that may
affect life and our ability to observe its signatures [e.g. via material
ejected in plumes].

Ocean circulation on icy moons can broadly speaking be driven by heat and salt
fluxes, tidal forcing, or magnetic forces (e.g. Soderlund et al., 2020, and
references therein). We here focus primarily on ”buoyancy driven” flow, i.e.
flows associated with temperature and salinity gradients, although we also
include a discussion of the role of tides and librations in driving vertical
mixing, which in turn affects the buoyancy field and associated flow (e.g.
Wunsch & Ferrari, 2004). Following most of the previous work on buoyancy-
driven flows, we will neglect magnetic forces, although they may be
significant on Jupiter’s moons (Gissinger & Petitdemange, 2019).

A number of studies have simulated the buoyancy-driven dynamics of ice covered
oceans both in the context of Snowball Earth (e.g. Ashkenazy et al., 2013;
Ashkenazy & Tziperman, 2016; Jansen, 2016) and icy moons (e.g. Soderlund et
al., 2014; Soderlund, 2019; Ashkenazy & Tziperman, 2021; Kang et al., 2020;
Zeng & Jansen, 2021), and have come to widely different conclusions, in
particular with regards to the characteristic current speeds in these oceans.
The Snowball Earth simulations of Ashkenazy et al. (2013), Ashkenazy &
Tziperman (2016) and Jansen (2016) broadly consistently show small Rossby-
number turbulent flows dominated by eddies and jets with characteristic
velocities on the order of $1\,$cm/s. For Europa, Soderlund et al. (2014) find
moderate Rossby-number convective turbulence, with characteristic velocity
scales on the order of111Results are reported in non-dimensionalized units as
$|U/(2\Omega D)|\sim 0.5$. Re-dimensionalizing with the rotation rate
$\Omega\sim 2\times 10^{-5}\,$s-1 and ocean depth $D\sim 1\times 10^{5}\,$m
gives $|U|\sim 2\,$m/s $1\,$m/s, whereas Ashkenazy & Tziperman (2016) instead
find largely geostrophic (i.e. low Rossby-number) turbulence and jets with
characteristic velocities on the order of $1\,$cm/s. In a parameter regime
deemed applicable to Enceladus, Soderlund (2019) finds characteristic flow
velocities on the order of222Results are reported in non-dimensionalized units
as $|U/(2\Omega D)|\sim 0.1$. Re-dimensionalizing with the rotation rate
$\Omega\sim 5\times 10^{-5}\,$s-1 and ocean depth $D\sim 1\times 10^{4}\,$m
gives $|U|\sim 0.1\,$m/s. $0.1\,$m/s, while Kang et al. (2020) and Zeng &
Jansen (2021) find velocities on the order of $0.1\,$mm/s. The flow dynamics
and associated kinetic energy levels hence vary widely across these studies,
with variations across different studies being larger than variations between
different oceans.

To shed some light on these discrepancies and to establish what insights can
be gained from first principles (i.e. without running numerical simulations
whose results are sensitive to many parameters and implicit but often unstated
assumptions), we here consider energetic constraints for the circulation of a
globally ice-covered ocean. For simplicity we limit ourselves to an ocean in a
statistical equilibrium state (as also assumed in all studies discussed in the
previous paragraph). Specifically, we assume that the net global ocean warming
or cooling is small compared to the heat fluxes through the lower and upper
ocean boundaries, and similarly that the net global mean freezing or melting
rate is small compared to the regional rates. Non-equilibrium effects could be
important (e.g. Hussmann & Spohn, 2004; Nakajima et al., 2019), but would
vastly widen the range of possible solutions. Our philosophy is that the
better-constrained equilibrium problem should serve as a null-hypothesis,
which will be rejected if and only if evidence contradicts the assumptions or
predictions of equilibrium ocean dynamics.

##  2 Energetics of the seawater Boussinesq Equations.

The weak compressibility of water allows us to employ the seawater Boussinesq
approximation, with which the dyamical equations reduce to (e.g. Young, 2010)

$\displaystyle\frac{D\mathbf{v}}{Dt}+2\mathbf{\Omega}\times\mathbf{v}+\nabla
p$ $\displaystyle=$ $\displaystyle b\hat{k}+\mathbf{\mathcal{F}}$ (1)
$\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=$ $\displaystyle 0$ (2)

where $\mathbf{v}$ is the velocity, $\mathbf{\Omega}$ is the planetary
rotation, $p$ is the pressure anomaly relative to a hydrostatic reference
state with constant density, $\rho_{0}$, $b=g(\rho_{0}-\rho)/\rho$ is
buoyancy333Notice that Young (2010) defines $b$ using a constant reference
gravity $g_{0}$ and writes the first term on the R.H.S. of Eq. (1) as b$\nabla
Z$, where Z is $\phi/g_{0}$, with $\phi$ the geopotential
($g\equiv|\nabla\phi|$), such that $\nabla Z=g/g_{0}\hat{k}$. We here absorb
the factor $g/g_{0}$ into the definition of buoyancy, assuming that g is
itself only a function of z., $\hat{k}=g^{-1}\nabla{\phi}$ is the normal
vector in the direction of gravitational acceleration, and
$\mathcal{F}=\mathcal{T}-\mathcal{D}$ is an acceleration due to tidal and
frictional forces, respectively (where we define the frictional force with a
negative sign for illustrative purposes, as it dominantly acts to decelerate
the flow).

Multiplying Eq. (1) by $\mathbf{v}$, and using $\nabla\cdot\mathbf{v}=0$
yields an equation for the kinetic energy as

$\frac{D|\mathbf{v}|^{2}}{Dt}+\nabla\cdot(\mathbf{v}p)=wb+\mathbf{v}\cdot\mathbf{\mathcal{T}}-\mathbf{v}\cdot\mathbf{\mathcal{D}}$
(3)

where $w$ is the velocity normal to the geopotential surfaces (in practice
usually well approximated by the radial velocity).

Integrating over the entire ocean volume and assuming an equilibrium state and
no-normal-flow boundary conditions444Notice that the no-normal flow boundary
condition here amounts to neglecting any possible kinetic energy injected by
jets emanating from the sea floor or ice shell (e.g. Kite & Rubin, 2016). The
effect of geothermal vents or freezing and melting at the ice-ocean interface
on buoyant plumes, however, is included via the heat and salt flux boundary
conditions., we find a balance between the conversion from potential to
kinetic energy, $wb$, KE generation by tidal forcing,
$\mathbf{v}\cdot\mathbf{\mathcal{T}}$, and frictional dissipation,
$\mathbf{v}\cdot\mathbf{\mathcal{D}}$:

$\int
wbdV+\int\mathbf{v}\cdot\mathbf{\mathcal{T}}dV=\int\mathbf{v}\cdot\mathbf{\mathcal{D}}dV$
(4)

We hence have two potential sources of kinetic energy: 1) the vertical
buoyancy flux, which converts potential to kinetic energy and is directly
related to the vertical heat and salt flux, and 2) tidal forcing. This paper
will focus primarily on buoyancy-driven circulations, that is we assume a
balance between the first term on the LHS of Eq. (4) and the dissipation on
the RHS. The potential role of tides in modulating the buoyancy-driven
circulation will also be discussed.

Importantly, Eq. (4) highlights that buoyancy forcing can drive a circulation
if and only if it is distributed such that the global mean buoyancy flux
required to balance the forcing is directed upwards.

###  2.1 Buoyancy forcing

The upward buoyancy flux is related to the upward flux of heat or
compositional buoyancy. In this manuscript we assume a water ocean with
dissolved salts, such that compositional buoyancy effects are encapsulated by
the salinity. In general we thus allow buoyancy to be some function of
potential temperature, salinity, and geopotential height (e.g. Young, 2010),
i.e.

$b=\tilde{b}(\Theta,S,z)\,,$ (5)

where $\Theta$ is potential temperature555Notice that temperature is conserved
following adiabatic motion in a Boussinesq fluid and hence temperature and
potential temperature are formally identical. However, $\theta$ should be
interpreted as potential temperature when comparing to observed fluids, $S$ is
salinity, and $z$ is depth relative to some reference geopotential height
level. Temperature and salinity evolve according to

$\displaystyle\frac{D\Theta}{Dt}$ $\displaystyle=$
$\displaystyle\kappa_{T}\nabla^{2}\Theta$ (6) $\displaystyle\frac{DS}{Dt}$
$\displaystyle=$ $\displaystyle\kappa_{S}\nabla^{2}S$ (7)

with the boundary conditions

$\displaystyle-\kappa_{T}\partial_{z}\Theta(z_{bot})=\frac{1}{\rho
c_{p}}\mathcal{Q}_{bot}\,,$
$\displaystyle-\kappa_{T}\partial_{z}\Theta(z_{top})$
$\displaystyle=\frac{1}{\rho c_{p}}\mathcal{Q}_{top}\,,$ (8)
$\displaystyle-\kappa_{S}\partial_{z}S(z_{bot})=\mathcal{S}_{bot}\,,$
$\displaystyle-\kappa_{S}\partial_{z}S(z_{top})$
$\displaystyle=\mathcal{S}_{top}\,,$ (9)

where $\kappa_{T}$ and $\kappa_{S}$ are the molecular diffusivities for heat
and salt, respectively, and $\mathcal{Q}_{bot/top}$ and
$\mathcal{S}_{bot/top}$ denote the heat and salt fluxes through the bottom and
top boundaries, respectively. For simplicity we here do not include
significant heat sources in the interior, although such sources could be
added.

If the equation of state is nonlinear, the conservation equations for $\Theta$
and $S$ cannot readily be translated into a conservation equation for
buoyancy. The most general way to account for nonlinearities in the equation
of state, which is discussed in Appendix A is to introduce the dynamic
enthalpy, which is related to the potential energy. A significantly simpler
(and more intuitive) result can be obtained if we assume that horizontal
variations in the thermal and haline expansion coefficients are small, such
that we can approximate

$b\approx\langle{b}\rangle+\langle{\alpha}\rangle(\Theta-\langle{\Theta}\rangle)+\langle{\beta}\rangle(S-\langle{S}\rangle)$
(10)

where $\alpha\equiv g^{-1}\partial_{\theta}b$ and $\beta\equiv
g^{-1}\partial_{S}b$ are the thermal and haline expansion coefficients and
$\langle{\cdot}\rangle$ denotes a horizontal average at any given depth. We
can then directly relate the total upward buoyancy flux (which represents the
source of kinetic energy in a buoyancy-driven flow) to the upward heat and
salt fluxes:

$\mathcal{B}\equiv\int wbdV\approx\int\left[g\langle{\alpha}\rangle\iint
w\Theta dA\right]dz+\int\left[g\langle{\beta}\rangle\iint
wSdA\right]dz\equiv\mathcal{B}^{\Theta}+\mathcal{B}^{S}$ (11)

where we used that $\iint wdA=0$ due to volume conservation, and we defined
$\mathcal{B}^{\Theta}$ and $\mathcal{B}^{S}$ as the total net upward buoyancy
fluxes at each level associated with heat and salt fluxes, respectively.

The total upward heat and salt fluxes can be related to the sources and sinks
of heat and salt via their conservation equations (see also sketches in Figs.
1 and 2 ). Specifically, the area-integrated vertical advective fluxes of heat
and salt are constrained by the conservation of heat and salt above or below a
respective level. In a statistically steady state, where the mean rate of
change of $\Theta$ and $S$ at any given depth is approximately zero, Eq. (6)
can be integrated over the volume above or below some depth $z$ to yield:

$\displaystyle\int(w\Theta-\kappa_{T}\partial_{z}\Theta)dA$ $\displaystyle=$
$\displaystyle\frac{1}{\rho
c_{p}}\left[\int\mathcal{Q}_{bot}\mathcal{H}(z-z_{bot})dA-\int\mathcal{Q}_{top}\mathcal{H}(z-z_{top})\
dA\right]$ (12) $\displaystyle=$ $\displaystyle\frac{1}{\rho
c_{p}}\left[\int\mathcal{Q}_{top}\mathcal{H}(z_{top}-z)dA-\int\mathcal{Q}_{bot}\mathcal{H}(z_{bot}-z)\
dA\right]$ $\displaystyle\equiv$ $\displaystyle\frac{1}{\rho
c_{p}}\mathcal{F^{Q}}$

where $\mathcal{H}$ is the Heaviside function and we defined $\mathcal{F^{Q}}$
as the the vertical heat and salt flux across any depth $z$ needed to balance
the net heat and salt fluxes through the boundaries above and below the
respective level. Similarly, Eq. ( (7)) can be integrated to obtain

$\displaystyle\int(wS-\kappa_{S}\partial_{z}S)dA$ $\displaystyle=$
$\displaystyle\int\mathcal{S}_{bot}\mathcal{H}(z-z_{bot})dA-\int\mathcal{S}_{top}\mathcal{H}(z-z_{top})\
dA$ (13) $\displaystyle=$
$\displaystyle\int\mathcal{S}_{top}\mathcal{H}(z_{top}-z)dA-\int\mathcal{S}_{bot}\mathcal{H}(z_{bot}-z)\
dA$ $\displaystyle\equiv$ $\displaystyle\mathcal{F^{S}}$

Notice that interior sources of heat or salt could be included by modifying
the definition of $\mathcal{F^{Q}}$ and $\mathcal{F^{S}}$, which generally
need to balance any net sources or sinks above or below the respective level.

Figure 1: Sketch of the vertical heat flux through the ocean, which provides
the energy source for thermally driven flows. In this manuscript we assume
that the ocean and ice sheet are in equilibrium, and that heat sources in the
interior of the ocean are negligible, such that the total heat flux from the
core ($\int Q_{bot}dA$) is equal to the net heat flux from the ocean to the
ice sheet ($\int Q_{top}dA$). Notice that the vertical heat flux in the ocean
($\mathcal{F^{Q}}$) can also be constrained from the potentially more
observable heat flux emanating from the planetary body’s surface ($\int
Q_{surf}dA$), which is given by the sum of the core heating ($\int Q_{bot}dA$)
and the tidal energy dissipation in the ice layer ($\int Q_{ice}dA$).
Specifically, the maximum vertical heat flux in the ocean is limited to
$\mathcal{F^{Q}}\leq\int Q_{surf}dA$ . Figure 2: Sketch of the vertical salt
flux through the ocean. If melting occurs under thinner ice and freezing
occurs under thicker ice, as sketched here, the salt flux required in the
ocean to balance brine rejection from freezing and melting ($\mathcal{F^{S}}$)
is upward, which implies a downward buoyancy flux that converts kinetic to
potential energy. In this scenario salt fluxes therefore do not energetically
drive a circulation. Instead, kinetic energy from an alternative source is
required to maintain a circulation that can flux salt upwards.

These definitions allow us to express the KE generation associated with heat
and salt flux forcing as:

$\mathcal{B}^{\Theta}\approx\int
g\langle{\alpha}\rangle\frac{\mathcal{F^{Q}}}{\rho c_{p}}dz+\int
g\alpha\kappa_{T}\partial_{z}\theta dV,\quad$ (14)

and

$\mathcal{B}^{S}=\int g\langle{\beta}\rangle\mathcal{F^{S}}dz+\int
g\beta\kappa_{S}\partial_{z}SdV,$ (15)

The first term on the RHS of equations (14) and (15) can be interpreted as the
potential energy source associated with heat and salt forcing at the lower and
upper boundaries, while the second term is the potential energy source or sink
due to diffusion (which provides a source of mechanical energy whenever the
stratification is statically stable, and a sink when it is unstable).

If vertical variations in $g\alpha$ and $g\beta$ are also small, thermal
forcing at the boundaries provides a source of energy that can drive a
circulation if, and only if,

$\alpha\int\mathcal{F}^{Q}dz>0\,,$ (16)

i.e. for $\alpha>0$ heating, on average, needs to occur at a greater depth
than cooling, while for $\alpha<0$ cooling needs to occur at greater depth
than heating. Similarly salt fluxes can drive a circulation only if

$\beta\int\mathcal{F}^{S}dV>0\,.$ (17)

Since generally $\beta<0$, this requires that salt needs to be removed at
greater depth and added at shallower depth. Large vertical variations in
$g\alpha$ and $g\beta$ can be accounted for using the generalized relation in
Eqs. (14) and(15), which allow us to compute the mechanical energy input
associated with any given heat and salt flux boundary conditions more
generally. When horizontal variations in the thermal and haline expansion
coefficient are non-negligible, the amount of energy that can be converted to
kinetic energy depends on the specifics of the circulation, as elaborated in
Appendix A.

###  2.2 Energy dissipation and flow properties

In equilibrium the total sources and sinks of mechanical energy have to be in
balance. The energy sources thereby provide a constraint on the energy
dissipation, which in turn provides some constraint on the flow. Relating the
energy dissipation rate to the kinetic energy of the flow itself is not
straightforward, as it depends on the characteristics of the flow field, but
we can make some progress by considering specific flow regimes and estimating
the range of parameters and scales over which the flow regimes are expected to
hold. In the following we first discuss the relationship between the kinetic
energy and the dissipation rate for a turbulent flow that is largely
unaffected by rotation, as well as the conditions under which the assumption
that rotation is negligible breaks down. We then derive an alternative
relationship between the kinetic energy and the dissipation rate in the
opposite limit of so-called geostrophic turbulence, where rotation is a
leading-order effect and flows are largely geostrophically balanced (i.e. the
leading order momentum balance is between the pressure gradient and Coriolis
acceleration).

####  2.2.1 Isotrophic turbulence and the role of rotation

We first assume a fully turbulent flow unconstrained by the influence of
rotation (i.e. high Reynolds and Rossby numbers) and we assume that, in this
limit, the dissipation rate is independent of the value of molecular
viscosity666Notice that it is unclear in how far this holds true in the
vicinity of boundaries. We here assume that most dissipation occurs in the
interior in this fully turbulent regime, but we will return to the issue of
boundary-layer dissipation below.. Specifically, we start by considering
Kolmogorov’s theory for the kinetic energy spectrum in the inertial cascade of
isotropic turbulence, which suggests that

$E(k)=\mathcal{K}\epsilon^{2/3}k^{-5/3}$ (18)

where $E(k)$ is the kinetic energy spectrum as a function of wavenumber, $k$,
$\epsilon$ is the turbulent spectral KE flux, which in turn is equal to the
dissipation rate, and $\mathcal{K}\approx 1.5$ is the Kolmogorov constant
(e.g. Vallis, 2006). Integrating over the energy inertial range yields

$E_{t}\approx 2k_{0}^{-2/3}\epsilon^{2/3}$ (19)

where $k_{0}$ is the ”injection scale” (where the KE spectrum flattens out),
and $E_{t}$ is the total turbulent kinetic energy in the inertial range.

We can use Eq. (19) to estimate the kinetic energy of turbulent flows up to
the largest scales of an isotropic turbulent energy inertial range. Isotropy
may be broken by the geometry (e.g. the vicinity to a boundary), or by the
effects of stratification or rotation. We will here consider in articular the
importance of rotation, which can break isotropy and fundamentally change the
nature of the turbulent flow throughout the watercolumn.

The role of rotation can be characterized by the Rossby number, which, using
Eq. (19) and $L\equiv 2\pi/k_{0}$, can be estimated as

$Ro_{t}=\frac{\sqrt{E_{t}}}{\Omega L}\approx\frac{\epsilon^{1/3}}{\Omega
L^{2/3}}\,.$ (20)

Equation (20) suggests a maximum length scale for turbulent flows unaffected
by rotation (c.f. Fernando et al., 1991; Jones & Marshall, 1993; Maxworthy &
Narimousa, 1994; Bire et al., 2022):

$L_{rot}\approx\frac{\epsilon^{1/2}}{\Omega^{3/2}}\,.$ (21)

Once rotation becomes important, Maxworthy & Narimousa (1994) and others find
that the convective Rossby number of a rotating plume scales as

$Ro_{rc}=\frac{\mathcal{\epsilon}^{1/2}}{\Omega^{3/2}H}=Ro_{t}^{3/2}\,,$ (22)

where $H$ is the depth of the convecting fluid (and we used Eq. (20) with
$L=H$ in the second equality). Eq. (22) is likely to be a better predictor of
the convective flow Rossby number when ocean-depth convection is strongly
affected by rotation (i.e. $L_{rot}<H$). However, both give the same
prediction for the length scale at which $Ro\approx 1$ and hence rotation
becomes important, which indeed also follows directly from dimensional
analysis if we postulate that this scale shall depend only on $\epsilon$ and
$\Omega$.

####  2.2.2 Boundary layer dissipation in geostrophic dynamics

When rotation becomes of dominant importance it is likely that much of the
energy becomes trapped in geostrophically balanced vorticies and large-scale
mean flows, which result from the up-scale kinetic energy cascade associated
with quasi-balanced turbulent motions (e.g. Vallis, 2006). The lack of a
forward energy cascade means that dissipation is likely to be limited mostly
to turbulent boundary layers near the seafloor and the ice-ocean interface.

We can estimate the energy dissipation per unit area in a turbulent boundary
layer as (e.g Jansen, 2016, and references therein):

$\int\epsilon dz=c_{D}|U_{g}|^{3}$ (23)

where the integral on the L.H.S. is over the depth of the turbulent boundary
layer, $c_{D}$ is the turbulent drag coefficient, and $U_{g}$ is the
characteristic near-surface geostrophic velocity777Notice that Eq. (23) also
follows from the Kolmogorov scaling in Eq. (19) if we take $k_{0}$ to scale as
the inverse of the bottom boundary layer depth, in which case $c_{D}\approx
E_{t}^{3/2}/|U|^{3}$, i.e. it is related to the ratio of the boundary layer
turbulent kinetic energy to the kinetic energy of the large-scale geostrophic
flow.. The value of $c_{D}$ depends on the surface roughness, with
$c_{D}\approx 0.0025$ an empirical average value that is commonly used for the
drag coefficient for Earth’s seafloor (e.g. Egbert et al., 2004; Sen et al.,
2008). The skin drag coefficient under smooth ice can be estimated to be
around $c_{D}\approx 0.002$, although rough morphology can significantly
increase this value (e.g. Brenner et al., 2021). Lacking information about the
ice roughness, and noting the order-of-magnitude nature of our estimates, we
will here use $c_{D}\approx 0.0025$ at both the sea floor and under the ice
sheet at the top. In addition to surface roughness, the drag coefficient may
depend on other flow properties, but experience from Ocean and Atmospheric
modelling has shown that drag coefficients of this order produce reasonable
results for a wide range of flows (e.g. Smagorinsky et al., 1965; Egbert et
al., 2004; Sen et al., 2008; Chen et al., 2018; Adcroft et al., 2019), which
instills confidence that a similar coefficient can be used to obtain useful
order-of-magnitude estimates for icy moon oceans.

Averaging the energy dissipation in the top and bottom boundary layers over
the depth of the whole water column gives an average energy dissipation rate
per unit volume

$\epsilon\approx\frac{2^{5/2}c_{D}}{H}E_{g}^{3/2}\,,$ (24)

where $H$ is the depth of the ocean and $E_{g}=1/2|U_{g}|^{2}$ is the
characteristic KE of the geostrophic flow. We can solve Eq. (24) to get a
crude estimate for the characteristic KE in a flow that is dominated by large-
scale balanced dynamics, where KE dissipation occurs primarily in turbulent
boundary layers:

$E_{g}\approx\frac{H^{2/3}}{2^{5/3}c_{D}^{2/3}}\epsilon^{2/3}\,.$ (25)

In general, balanced dynamics may co-exist with intermediate and high Rossby
number convective turbulence and/or tidal waves, in which case interior energy
dissipation may be important and boundary layer dissipation may be enhanced
due to the tidal flows. In this case, Eq. (25) is expected to overestimate the
KE. It is also important to note that numerical simulations of planetary
circulation typically use very high artificial viscosities for numerical
stability, which can be a large, albeit probably unphysical, source of
dissipation for balanced KE (e.g. Jansen, 2016).

##  3 Scaling laws for the kinetic energy of thermally-driven flows on
Snowball Earth, Enceladus and Europa

We can derive estimates for the kinetic energy of buoyancy-driven flows by
assuming a statistically steady state where the source of KE equals the sink.
If the primary source of energy is given by thermal forcing, Eq. (4) together
with Eq. (11) and (14) provide a scaling relation for the mean kinetic energy
dissipation per unit volume:

$\boxed{\epsilon\equiv\frac{1}{V}\int\mathbf{v}\cdot\mathcal{D}\sim\frac{\mathcal{B}^{\theta}}{V}\approx\frac{g\alpha\mathcal{F}^{Q}H}{\rho
c_{p}V}\approx\frac{Q\alpha g}{\rho c_{p}}\,,}$ (26)

where $Q$ is the average heat flux through the sea floor per unit area, and we
for now neglect variations in the thermal expansion coefficient, $\alpha$, as
well as radial variations in the gravity and the surface area throughout the
depth of the ocean. (The latter assumption leads to O(10%) errors for Europa
and Enceladus, which is not of concern here. The former may lead to larger
errors if the ocean is relatively fresh, and we will return to this issue
later.) We also ignored the energetic effect of molecular diffusion (the
second term in Eq. (14)), which is likely to be small888Notice that his
formally amounts to assuming a large Nusselt number - i.e. we assume that the
advective heat flux is large compared to the diffusive heat flux..

Eq. (26) provides a key constraint for the energetics of thermally-driven
flows. It predicts that the energy dissipation rate per unit volume increases
linearly with the bottom buoyancy flux, which in turn is given by the heat
flux multiplied by the thermal expansivity and the gravity. The results are
shown in Fig. 3 for the parameter range occupied by the icy moons and Snowball
Earth. In the following, we use this result together with the results from
sections 2.1 and 2.2 to obtain order-of-magnitude estimates for thermally-
driven ocean flows on Snowball Earth, Europa, and Enceladus.

Figure 3: Energy dissipation rate per unit volume for thermally-driven flows
predicted by the scaling in Eq. (26) as a function of the bottom heat flux
($Q$) times $\alpha/(\rho c_{p})$ and the gravitational acceleration ($g$).
The white stars mark the estimates for Snowball Earth (Ea), Europa (Eu) and
Enceladus (En), assuming the parameters given in table 1. The colorbar is
logarithmic with the contour interval a factor of $\sqrt{10}$.

To estimate the kinetic energy dissipation rate for thermally-driven flows in
the ocean’s of Snowball Earth, Europa and Enceladus, we assume an average
heating rate at the ocean floor of $Q\approx 0.1\,$W/m2 for Earth (Jansen,
2016) and $Q\approx 0.03\,$W/m2 for Europa and Enceladus (Ruiz, 2005;
Hemingway et al., 2018; Choblet et al., 2017), $g\approx 10\,$m/s2 for Earth,
$g\approx 1\,$m/s2 for Europa and $g\approx 0.1\,$m/s2 for Enceladus. The
thermal expansion coefficient is uncertain, in particular for Enceladus, as it
depends on pressure and ocean salinity. Indeed for salinities below about 20
g/kg the thermal expansion coefficient at the ice-ocean interface on Enceladus
is expected to be negative Zeng & Jansen (2021); Kang et al. (2021). For now
we will assume $\alpha\approx 2\times 10^{-4}\,$K-1 for Earth and Europa and
$\alpha\approx 1\times 10^{-4}\,$K-1 for Enceladus, although we note that this
is likely to be an upper-end estimate, and we will return to a low salinity-
scenario for Enceladus below. (Assumed parameters are also summarized in table
1.) Using further that $\rho c_{p}\approx 4\times 10^{6}\,$Jm${}^{-3}\,$K-1 we
find

$\displaystyle\epsilon$ $\displaystyle\approx$ $\displaystyle 5\times
10^{-11}\,\mbox{m${}^{2}$s${}^{-3}$}\quad\mbox{for Snowball Earth}$ (27)
$\displaystyle\epsilon$ $\displaystyle\approx$ $\displaystyle 2\times
10^{-12}\,\mbox{m${}^{2}$s${}^{-3}$}\quad\mbox{for Europa}$ (28)
$\displaystyle\epsilon$ $\displaystyle\approx$ $\displaystyle 8\times
10^{-14}\,\mbox{m${}^{2}$s${}^{-3}$}\quad\mbox{for Enceladus}\,.$ (29)

The energy dissipation rate is thus largest for Snowball Earth and smallest
for Enceladus, with differences mostly driven by the differences in the
gravity and further amplified by differences in the estimated vertical heat
flux and thermal expansivity. For comparison, the kinetic energy dissipation
rate in Earth’s deep ocean today is around $2\,$TW (Wunsch & Ferrari, 2004),
which divided by the mass of the ocean amounts to about $1.5\times
10^{-9}\,$m2s-3. This dissipation is balanced by energy input primarily from
winds and tides, which provide significantly more energy to the present-day
ocean circulation than geothermal heating.

###  3.1 The role of rotation

Given the energy dissipation rate, $\epsilon$, and the rotation rate of the
planetary body, $\Omega$, we can use Eq. (21) to estimate a critical scale
$L_{rot}$ above which we may expect turbulent motions to become strongly
affected by rotation (Fig. 4). Using the energy dissipation rates estimated
above and $\Omega\approx 8\times 10^{-5}\,$s-1 for Snowball Earth, $5\times
10^{-5}\,$s-1 for Enceladus, and $2\times 10^{-5}\,$s-1 for Europa, we find
the scale at which convection is expected to become strongly affected by
rotation from Eq. (21) to be $L_{rot}\approx 10\,$m for Snowball Earth and
Europa and $L_{rot}\approx 1\,$m for Enceladus. Uncertainties in many of our
assumed parameters are substantial (and not easily quantified), such that
these results should only be viewed as indicative of the expected order of
magnitude. Importantly, however, it is clear that $L_{rot}$ is much smaller
than the ocean depth in all three cases and we therefore expect deep
convection and any other potential large-scale flow to be strongly affected by
rotation. It is also worth noting that this length scale is significantly
smaller than the smallest length scales that can be resolved in global-scale
numerical simulations of these oceans, indicating that “large-eddy-
simulations” that can explicitly resolve part of the isotropic turbulent
inertial range are not feasible with realistic parameters (c.f. Bire et al.
(2022)).

Figure 4: Estimated length scale above which turbulence is expected to become
strongly affected by rotation, as predicted by the scaling in Eq. (21), as a
function of the turbulent dissipation rate ($\epsilon$) and the rotation rate
($\Omega$). The white stars mark the estimates for Snowball Earth (Ea), Europa
(Eu) and Enceladus (En), assuming the parameters given in table 1. The
colorbar is logarithmic with the contour interval a factor of $\sqrt{10}$.

#### Comparison to existing simulation results

The dominant importance of rotation in the oceans of Snowball Earth, Europa,
and Enceladus is qualitatively consistent with the Snowball Earth simulations
of Ashkenazy et al. (2013) and Jansen (2016), the Europa simulations of
Ashkenazy & Tziperman (2021) and the Enceladus simulations of Kang et al.
(2020) and Zeng & Jansen (2021). However, the result appears to be at odds
with simulation results for Europa’s ocean by Soderlund et al. (2014), which
show strong turbulent convection that is only relatively weakly affected by
rotation. The apparrent contradiction can readily be understood by noting that
the vertical heat flux in the simulations of Soderlund et al. (2014) (which
use fixed temperature boundary conditions) amounts to about $10^{5}\,$Wm-2—in
stark contrast to the heat flux of $0.03\,$Wm-2 assumed here. Increasing the
vertical heat flux on Europa to $10^{5}\,$Wm-2 would yield $L_{rot}\approx
25\,$km, which in turn would make significant turbulent convection that is
only relatively weakly affected by rotation plausible. Indeed, using Eq. (20)
with $L\approx H\approx 100\,$km (appropriate for Europa’s ocean) and an
energy input consistent with a heat flux of $10^{5}\,$Wm-2 (which gives
$\epsilon\approx 5\times 10^{-6}\,$m2s-3) we can estimate a Rossby number for
full-depth turbulent convection of $Ro_{t}\approx 0.4$, which is broadly
consistent with the Rossby numbers of convective motions found in the
simulations of Soderlund et al. (2014). However, using the more realistic
estimate of $\epsilon\approx 2\times 10^{-12}\,$m2s-3, Eq. ( 20) gives
$Ro_{t}\approx 3\times 10^{-3}\,$, again confirming that rotation is expected
to be very important. Using the Rossby number for rotating convection in Eq.
(22) further reduces our estimated convective Rossby number to $Ro_{rc}\approx
1\times 10^{-4}$.

###  3.2 Quasi-balanced flows

At scales much larger than $L_{rot}$, we expect the motion to be strongly
affected by rotation, leading to quasi-balanced flows that do not undergo a
forward energy cascade. Energy dissipation then may be limited mostly to
turbulent boundary layers. In this case, Eq (25) may provide a useful estimate
for the total KE and thus typical flow speeds as a function of the energy
dissipation rate, $\epsilon$, and the depth of the ocean, $H$ (Fig. 5).
Assuming that $H=3\,$km, $100\,$km, and $30\,$km for Snowball Earth, Europa
and Enceladus, respectively, we find

$\displaystyle E_{g}$ $\displaystyle\approx$ $\displaystyle 5\times
10^{-4}\,\mbox{m${}^{2}$s${}^{-2}$}\quad\mbox{for Snowball Earth}$ (30)
$\displaystyle E_{g}$ $\displaystyle\approx$ $\displaystyle 5\times
10^{-4}\,\mbox{m${}^{2}$s${}^{-2}$}\quad\mbox{for Europa}$ (31) $\displaystyle
E_{g}$ $\displaystyle\approx$ $\displaystyle 3\times
10^{-5}\,\mbox{m${}^{2}$s${}^{-2}$}\quad\mbox{for Enceladus.}$ (32)

These results suggest potential balanced flow velocities of up to a few cm/s
for Snowball Earth and Europa and a few mm/s for Enceladus. These estimates
should likely be viewed as an upper bound as we here assumed that all energy
input goes into balanced motions and we ignored any additional potential
routes to dissipation (e.g. via energy transfer to internal waves) that would
reduce the expected kinetic energy level.

Figure 5: Estimated geostrophic flow speeds ($U_{g}=\sqrt{2E_{g}}$) as a
function of the energy dissipation rate ($\epsilon$) and ocean depth ($H$) as
predicted by the scaling in Eq. (25) assuming that most dissipation happens in
turbulent boundary layers near the sea floor and ice-ocean interface with a
drag coefficient $C_{d}=0.025$. The white stars mark the estimates for
Snowball Earth (Ea), Europa (Eu) and Enceladus (En), assuming the parameters
given in table 1. The colorbar is logarithmic with the contour interval a
factor of $10^{1/4}$.

#### Comparison to existing simulation results

Flow velocities of a couple of cm/s for Snowball Earth are broadly consistent
with the GCM simulations of Ashkenazy et al. (2013) and Jansen (2016), both of
which reveal flow patterns dominated by geostrophic turbulence and mean flows.
A flow velocity of a couple of cm/s for Europa is also broadly consistent with
the simulation results of Ashkenazy & Tziperman (2021).

The Enceladus Ocean simulations with sea-floor heating of Kang et al. (2020)
and the high-salinity simulations of Zeng & Jansen (2021) have flow velocities
on the order of $0.1\,$mm/s, which is at least an order of magnitude weaker
than our estimate for the maximum KE. Some discrepancy may be expected due to
the parameterized boundary-layer drag. Instead of the quadratic boundary layer
drag assumed in Eq. (24), the simulations of Kang et al. (2020) and Zeng &
Jansen (2021) apply a relatively strong linear drag paramaterization. With a
linear drag parameterization the average dissipation rate due to drag at the
boundaries is

$\epsilon\approx 2n_{b}\frac{\gamma}{H}E_{g}$ (33)

where $\gamma$ is the linear drag coefficient, and $n_{b}=1$ in the setup of
Zeng & Jansen (2021) where frictional drag with $\gamma=1\times 10^{-4}\,$m/s
is applied only at the lower boundary, while $n_{b}=2$ in the setup of Kang et
al. (2020) where frictional drag with $\gamma=2\times 10^{-3}\,$m/s is applied
at both the upper and lower boundaries. Replacing Eq. (24) by Eq. (33) in our
estimate of the maximum geostrophic kinetic energy for Enceladus yields
$E_{g}\approx 10^{-5}\,$m2s-2 with the frictional drag of Zeng & Jansen (2021)
and $E_{g}\approx 3\times 10^{-7}\,$m2s-2 with the frictional drag of Kang et
al. (2020). This amounts to characteristic velocities of around $4\,$mm/s and
somewhat less than $1\,$mm/s respectively, which is closer to but still
significantly larger than the velocities observed in the simulations. Strong
dissipation of kinetic energy associated with the (numerically necessary)
large interior viscosity, as well as reduction of the conversion of potential
to kinetic energy due the parameterized vertical eddy diffusivity (i.e. the
second term on the R.H.S. of Eq. 14), are likely to contribute to reduced flow
speeds in the simulations of Kang et al. (2020) and Zeng & Jansen (2021).
Notice also that the choice of a relatively strong linear drag is practical in
numerical simulations (as it allows the model to equilibrate faster and
improves numerical stability), but is not justified on a physical basis. We
therefore conclude that the simulations of Kang et al. (2020) and Zeng &
Jansen (2021) likely underestimate the strength of the buoyancy-driven
circulation.

|  $Q$ [Wm-2] |  $g$ [ms-2] |  $\alpha$ [K-1] |  $H$ [m] |  $\Omega$ [s-1] |  $\epsilon$ [m2s-3] |  $L_{rot}$[m] | $Ro_{rc}$ |  $E_{g}$[m2s-2]
---|---|---|---|---|---|---|---|---|---
Snowball Earth | 0.1 | 10 |  2$\times 10^{-4}$ |  3$\times 10^{3}$ |  8$\times 10^{-5}$ |  5$\times 10^{-11}$ | 10 |  3$\times 10^{-3}$ |  5$\times 10^{-4}$
Europa | 0.03 | 1 |  2$\times 10^{-4}$ |  1$\times 10^{5}$ |  2$\times 10^{-5}$ |  1$\times 10^{-12}$ | 10 |  1$\times 10^{-4}$ |  5$\times 10^{-4}$
Enceladus | 0.03 | 0.1 |  1$\times 10^{-4}$ |  3$\times 10^{4}$ |  5$\times 10^{-5}$ |  8$\times 10^{-14}$ | 1 |  3$\times 10^{-5}$ |  3$\times 10^{-5}$
Table 1: Overview of assumed parameters for the sea floor heating rate, $Q$,
gravitational acceleration, $g$, thermal expansion coefficient, $\alpha$,
ocean depth , $H$, and Coriolis parameter, $f$ and resulting estimates for
mechanical energy input (per unit mass) due to thermal forcing,
$\epsilon=\mathcal{B}/H$, the critical length scale where rotation is expected
to become important, $L_{rot}$, the convective Rossby number for rotating
convection, $Ro_{rc}$, and an (upper-end) estimate of characteristic KE energy
of the large-scale geostrophic flow, $E_{g}$, for Snowball Earth, Europa, and
Enceladus.

###  3.3 Thermally-driven flows in a low-salinity, low-pressure ocean

If salinity and pressure are relatively low, the thermal expansion coefficient
near the freezing point is negative. This scenario was first suggested for
Europa by Melosh et al. (2004), although due to the modest pressures that are
required for a negative thermal expansion coefficient, it appears most likely
to be relevant for Enceladus, and was considered in the low-salinity Enceladus
ocean simulations of Zeng & Jansen (2021) and Kang et al. (2021). To estimate
the energetics of such an ocean we assume a stably stratified layer below the
ice sheet and above some depth $z_{strat}$, as found in the simulations of
Zeng & Jansen (2021). In the stratified layer, the temperature decreases
upwards from the “critical temperature”, $\Theta_{c}$, where
$\alpha(\Theta_{c},z_{strat})=0$ to the freezing point, $\Theta_{f}$, at the
bottom of the ice sheet, $z_{ice}$ (where $\alpha(\Theta_{f},z_{ice})<0$). For
illustrative purposes, we assume that the sea floor and ice sheet are flat
(i.e. both boundaries follow a geopotential height surface), such that
$\mathcal{F^{Q}}$ is vertically constant and equal to the total bottom heat
flux. Assuming again that horizontal variations in the thermal expansion
coefficient are small, Eq. (14) gives the KE generation as

$\mathcal{B}^{\Theta}=\int
g\langle{\alpha}\rangle\left[\frac{\mathcal{F^{Q}}}{\rho
c_{p}}+A\kappa_{T}\partial_{z}\langle{\theta}\rangle\right]dz$ (34)

In the stratified layer, $\langle{\alpha}\rangle<0$, such that the first term
in the integral in (34) amounts to an energy sink. One plausible solution is
then that the vertical heat flux through the stratified layer is balanced by
molecular diffusion (i.e. the second term in the integral in 34), which
amounts to a conversion from internal to potential energy. In this case the
stratified layer would only be a few tens of meters thick (Zeng & Jansen,
2021), and most of the ocean would be convective with $\Theta>\Theta_{c}$ and
hence $\alpha>0$. The dynamics in the convective layer would follow scaling
laws analogous to those presented above, although the very small thermal
expansion coefficient at temperatures only marginally warmer than $\Theta_{c}$
implies that the energy input and hence typical flow speeds would be even
smaller than estimated above.

It is theoretically possible that kinetic energy is transferred from the
convective layer into the stratified layer, which could allow for a turbulent
circulation throughout the depth of the ocean as long as $\int
g\langle{\alpha}\rangle dz>0$ (i.e. the net mechanical energy input is
positive), which in turn is possible if the stratified layer is sufficiently
shallow. More plausibly, tides or librations may provide a source of kinetic
energy that can balance the negative $\mathcal{B}^{\Theta}$ associated with
turbulent vertical heat transport in the stratified layer.

The simulations of Zeng & Jansen (2021) use a relatively large “eddy
diffusivity”, which is much larger than the molecular diffusivity. The
vertical heat flux through the stratified layer is balanced by this “eddy
diffusion”, such that from the perspective of the model’s energetics, the two
terms in the integral in (34) balance in the stratified layer. The simulations
of Kang et al. (2021) use an even larger eddy diffusivity, such that the
entire depth of the ocean becomes stratified with the upward heat flux and
associated downward buoyancy flux accomplished by “eddy diffusion”. However,
in both studies the model’s diffusivities, which are much larger than the
molecular value, need to be interpreted as representing turbulent mixing. If
these turbulent mixing rates are real, the associated vertical heat flux
represents turbulent advection, which amounts to a negative
$\mathcal{B}^{\Theta}$ \- that is a conversion of turbulent kinetic energy to
large-scale potential energy. An additional implied energy source is therefore
needed to generate the turbulent kinetic energy. A possible source are tides
and/or librations, to which we will return below.

##  4 Salinity Forcing

A number of studies have recently argued for the importance of salt forcing in
driving a circulation on Europa and Enceladus (Zhu et al., 2017; Kang et al.,
2020; Ashkenazy & Tziperman, 2021; Lobo et al., 2021). Since $\beta$, as
defined in Eq. (10), is always negative, salinity forcing provides a source of
mechanical energy if, and only if, salt is added (e.g. by freezing) at a
shallower depth and removed (or freshwater is added by melting) at a greater
depth (Eqs. 13 and 15). At steady state, however, freezing and melting need to
balance the convergence of ice flow, which tends to flow from regions of thick
ice to regions of thin ice. We therefore expect melting to occur where the ice
is thin, and hence the ice-ocean interface is shallow, while freezing occurs
where the ice is thick and hence the ice-ocean interface is deep (c.f. Fig.
2). In this case the required upward salt flux reduces the kinetic energy of
the ocean and hence cannot, energetically speaking, drive a circulation
(although it can still affect the circulation by shaping the horizontal
buoyancy gradients).

As for virtually all numerical simulations, the models of Zhu et al. (2017),
Kang et al. (2020), Ashkenazy & Tziperman (2021) and Lobo et al. (2021) employ
vertical diffusivities that are orders of magnitude larger than the molecular
value. In this case vertical diffusion acting against a statically stable
stratification can drive a substantial circulation, via the second term in Eq.
(15). Indeed, salinity forcing will often tend to increase the stratification
at least in some parts of the ocean, which will lead to an increase in the
energy input by the prescribed vertical diffusion. However, as in the low-
salinity Enceladus simulations of Zeng & Jansen (2021) discussed above, the
large vertical diffusivities need to be interpreted as representing mixing by
unresolved small-scale turbulence. In reality, this turbulent mixing requires
a source of kinetic energy, and the mechanistic basis for this energy source
must be established in order for the results of the simulations to be sound.
In Earth’s ocean today the source of energy for small-scale turbulent kinetic
energy ultimately comes from winds and tides (e.g. Wunsch & Ferrari, 2004).
Surface wind stress is absent in globally ice covered oceans, which leaves
tides as the most obvious source of turbulent kinetic energy in an ocean where
buoyancy forcing is primarily due to salt fluxes at the ice-ocean interface.

##  5 Tides and Turbulent Mixing

###  5.1 Theory

Tidal forcing generates kinetic energy in the ocean (via the second term on
the LHS of Eq. 4). The tidal perturbation potential is typically approximately
vertically constant throughout the depth of the ocean, thus generating
barotropic (i.e. vertically constant) tides. In Earth’s ocean these barotropic
tidal waves are largely linear until they encounter shallow shelfs, where the
tidal amplitude increases and most of the tidal energy dissipation is assumed
to occur (Wunsch & Ferrari, 2004). The icy moons do not have shelf seas
similar to Earth’s ocean, although tidal dissipation may nevertheless be
significantly enhanced in specific regions (e.g. Hay & Matsuyama, 2019).
Linear tides in the deep ocean do not contribute significantly to the
transport of heat, salt and other properties in the ocean, such that global
models of Earth’s ocean circulation can reproduce realistic large-scale
circulations and tracer distributions without explicitly accounting for tides.
However, it is possible that tidal waves on some icy moons become sufficiently
non-linear to lead to substantial rectified mean flows (e.g. Huthnance, 1981;
Brink, 2011). Since the current manuscript is focussed on buoyancy-driven
circulation, we will not further consider the possible role of rectified tidal
flows here, but note that this remains a potentially important topic for
further research.

Despite weak direct interactions with the large-scale flow, ocean tides can
have an important effect on the large-scale circulation in the presence of a
statically stable stratification (e.g. driven by freshwater fluxes at the ice-
ocean interface). Barotropic tides interacting with a rough sea-floor
topography in the presence of a statically stable stratification can transfer
their energy into baroclinic tides, i.e. internal waves that propagate both
horizontally and vertically and can be associated with substantial vertical
shears. These shears in turn can lead to instabilities and eventually the
generation of turbulence within the water column. This pathway is believed to
be one of the main drivers of 3D turbulence in Earth’s deep ocean (Vic et al.,
2019). Although most of the the turbulent kinetic energy is dissipated into
heat, some fraction, usually denoted by the ”mixing efficiency”, $\Gamma$, is
converted to potential energy via vertical mixing of the stratified water
column, i.e.:

$\mathcal{B}_{t}=\int(\overline{w^{\prime}b^{\prime}})_{t}dV=-\int\Gamma\epsilon_{t}dV$
(35)

where $-\mathcal{B}_{t}$ is the potential energy generation due to tidally-
driven turbulent mixing, $(\overline{w^{\prime}b^{\prime}})_{t}$ is the
vertical buoyancy flux associated with tidally-driven turbulence,
$\epsilon_{t}$ is the tidally generated turbulent kinetic energy dissipation
rate per unit mass, and $\Gamma\lesssim 0.2$ is the mixing efficiency (e.g.
Peltier & Caulfield, 2003).

The potential energy generated by turbulent mixing in a stratified ocean can
then drive a large-scale circulation that converts the potential energy back
to kinetic energy via a positive conversion $\overline{w}\overline{b}$ (where
the overbars denote the large-scale circulation). If turbulent mixing is the
dominant source of potential energy (i.e. when buoyancy gain and buoyancy loss
from external forcing occur at approximately the same depth) the conversion
$\overline{w}\overline{b}$ is, on average, approximately equal and opposite to
the downward buoyancy flux associated with tidally-driven turbulence:

$\mathcal{B}\equiv\int
wbdV=\int\overline{w}\overline{b}dV+\int(\overline{w^{\prime}b^{\prime}})_{t}dV\equiv\overline{\mathcal{B}}+\mathcal{B}_{t}\approx
0\,.$ (36)

where $\overline{\mathcal{B}}$ is the kinetic energy generated (and ultimately
dissipated) by the mean flow. The energy cycle can then be summarized as
follows. Tides generate turbulent kinetic energy, a fraction of which is
converted into large-scale potential energy via downward turbulent buoyancy
flux (with the remainder being dissipated into heat). This potential energy
can then be converted back to kinetic energy (and ultimately be dissipated) by
the large-scale circulation. This mechanism requires a statically stable
stratification such that turbulent mixing transports buoyancy downwards. In
the absence of a statically stable stratification, we still expect tidally-
driven turbulence to contribute to the mixing of properties, but we do not
expect it to contribute to driving a large-scale circulation.

In addition to tides, librations and injection of fluid through fissures in
the ice or solid core may provide sources of kinetic energy. Librations are
expected to play a similar role to tides, and as for tides, their effect on
the large-scale circulation is expected to depend on their ability to generate
turbulence in the ocean interior and on the presence of a statically stable
stratification. Injection of water trough fissures in the boundaries (which is
not formally included in the global kinetic energy-budget in Eq. (4), where we
assumed no-normal-flow boundary conditions) is likely to drive mechanical
turbulence primarily in the vicinity of the boundaries, but strong jets
emanating from the ice shell, as proposed by Kite & Rubin (2016) may play an
important role in the dynamics of the upper ocean999For steady planar
turbulent jets in an unstratified fluid, the peak mean flow speed has been
found to decay with distance from the orifice as $U_{max}\approx
2.5\sqrt{d/x}$, where $d$ is the width of the slot from which the jet emanates
and $x$ is the distance (e.g. Gutmark & Wygnanski, 1976). Turbulent flow
speeds are about $20-25\%$ of this value. Kite & Rubin (2016) suggest that
Enceladus’ ”tiger stripes” are associated with O($1\,$m) wide slots in which
tidal forces generate O($1\,$m/s) jets. Assuming a neutral stratification, the
results for planar turbulent jets would then suggest peak mean flows on the
order of $10\,$cm/s and associated turbulent velocities on the order of a few
cm/s to persist to a depth of a few hundred meters into the ocean. The
relatively weak depth dependence ($U_{max}\propto x^{-1/2}$) moreover
indicates that weaker but significant flows may be penetrating much deeper,
although results for steady jets are likely to become less applicable to the
oscillatory jets proposed by Kite & Rubin (2016) at greater depth.
Stratification may further limit the penetration depth of turbulent jets.. A
full investigation of the effects of such jets is beyond the scope of this
study, but remains as an interesting subject for future research.

###  5.2 Application to Snowball Earth, Europa, and Enceladus

We begin by estimating the potential role of tidally-driven turbulence in the
Snowball Earth ocean. Although little is known specifically about tidal energy
dissipation during Snowball Earth periods, average tidal dissipation since the
Precambrian was likely somewhat smaller but of similar order of magnitude as
the present-day value (Williams, 2000; Green et al., 2017). The direct effect
of the Snowball Earth ice sheet on tides was likely small (Wunsch, 2016). For
scaling purposes we therefore here assume a tidal energy input on the order of
10${}^{12}\,$W, which is comparable to the present day value of around
3.5$\times 10^{12}\,$W (Wunsch & Ferrari, 2004). Assuming a mean ocean depth
of about $2\,$km, an ocean area of 3.5$\times 10^{14}\,$m2 and density
$\rho_{0}\approx 1000\,$kg m-3, the average tidal energy dissipation per unit
mass is around $\epsilon_{t}\approx 10^{-9}\,$m${}^{2}\,$s-3, which is a
little over an order of magnitude larger than the estimated energy input due
to thermal forcing (cf. Table 1). If a significant fraction of the dissipated
tidal energy is dissipated in the stratified ocean interior and is relatively
efficiently converted to potential energy (i.e. $\Gamma\gtrsim$ 10%) it thus
may play a significant role in driving a large-scale ocean circulation.
Specifically, one may envision a scenario where freezing and melting at the
ice-ocean interface combined with tidally-driven turbulent mixing, which
pushes the light fresh water into the ocean interior, can drive a
significantly stronger large-scale ocean circulation then the thermal forcing
alone. This scenario appears to be relevant for the interpretation of
simulation results by Ashkenazy et al. (2013) and Jansen (2016), although the
turbulent vertical mixing in these simulations is parameterized with
prescribed vertical diffusivities (i.e.
$(\overline{w^{\prime}b^{\prime}})_{t}\rightarrow-\kappa\partial_{z}\overline{b}$)
and is thus not constrained by tidal energy dissipation. The present analysis
nevertheless suggests that the vertical mixing may be justifiable on energetic
grounds by the expected tidally-driven turbulence, although it would require
that a significant fraction of the tidal energy is dissipated in the
stratified ocean interior (as opposed to being highly localized in boundary
layers where the stratification is vanishingly small) and is relatively
efficiently converted to potential energy.

Tidal energy dissipation on Europa has been estimated from modeling to be on
the order of 109-10${}^{11}\,$W (Chen et al., 2014; Matsuyama et al., 2018;
Hay & Matsuyama, 2019). Assuming, as before, an ocean depth on the order of
10${}^{5}\,$m, an area of around 3$\times 10^{13}\,$m2, and a density
$\rho_{0}\approx 1000\,$kg m-3, we estimate a tidal energy dissipation per
unit volume of around $\epsilon_{t}\approx 3\times 10^{-13}-3\times
10^{-11}\,$m${}^{2}\,$s-3, which is between about an order of magnitude
smaller and about an order of magnitude larger than the energy input by
buoyancy forcing (cf. Table 1). Depending on the estimate for tidal energy
dissipation, and the fraction of the dissipated tidal energy that is converted
to potential energy, tidally-driven turbulent mixing thus may or may not play
a signifcant role in driving a large-scale ocean circulation on Europa.
Specifically, freezing and melting at the ice-ocean interface combined with
tidally-driven turbulent mixing may drive a significantly stronger large-scale
ocean circulation then the thermal forcing alone, if, and only if, ocean tidal
dissipation is near the upper end of these estimates and a large fraction of
that energy contributes to vertical mixing against a stable stratification
(i.e. $\Gamma\gtrsim$ 10%). A significant salt-driven circulation was found in
the simulations of Ashkenazy & Tziperman (2021), although turbulent vertical
mixing is parameterized with a prescribed Earth-like vertical diffusivity in
the simulations and is not constrained by tidal energy dissipation. Better
estimates of tidal energy dissipation in Europa’s ocean are needed to
determine whether a significant turbulent vertical diffusivity can be
justified for stratified regions of Europa’s ocean.

For Enceladus’ Ocean, energy dissipation associated with the eccentricity and
obliquity tides has been estimated to be relatively small, with a total
dissipation between 10-10${}^{4}\,$W suggested by Matsuyama et al. (2018), and
Hay & Matsuyama (2019). Tyler (2020) has argued that much larger dissipation
rates are possible, accounting for virtually all of the observed
O(10${}^{10}\,$W) heat flux, although the calculations generally predict the
total dissipation in the ice-ocean system, and it is likely that most of the
dissipation in the suggested scenarios would indeed occur in the ice sheet
(c.f. Beuthe, 2016). Assuming a most likely range of 10-10${}^{4}\,$W for
ocean tidal dissipation on Enceladus, an average ocean depth of around
3$\times 10^{4}\,$m and an area of around $5\times 10^{11}\,$m2, we find an
average tidal energy dissipation of around $\epsilon_{t}\approx 7\times
10^{-19}-7\times 10^{-16}\,$m${}^{2}\,$s-3 , which is at least two orders of
magnitude smaller than the estimated mechanical energy input via thermal
forcing (cf. Table 1). Mixing driven by the eccentricity and obliquity tides
thus is not likely to provide a significant source of mechanical energy for
the circulation on Enceladus’ ocean.

However, Enceladus ocean may be significantly affected by librations of the
ice shell, which is decoupled from the solid interior (Thomas et al., 2016).
It is not clear whether librations lead to relatively modest dissipation
confined to the ocean-ice boundary layer, thus not contributing significantly
to interior mixing, or whether the excitation of internal waves or elliptical
instability (an instability that can break up elliptical streamlines) can lead
to relatively strong mixing throughout the depth of the ocean (Lemasquerier et
al., 2017; Wilson & Kerswell, 2018; Rekier et al., 2019; Soderlund et al.,
2020). At the extreme end, Wilson & Kerswell (2018) argue that it is possible
that all of the O(10GW) of heating on Enceladus may be a result of librational
dissipation, which, divided by the mass of Enceladus’ ocean would give
$\epsilon_{t}\sim 10^{-9}$m${}^{2}\,$s-3—many orders of magnitude larger than
the thermal energy input. In this case it is likely that libration-driven
turbulence would play the dominant role in ocean mixing. Narrowing down the
large uncertainty in the energy dissipation rate associated with libration-
driven turbulence will be key to constraining the turbulent diffusivity, which
is an essential parameter in numerical simulations of Enceladus’ ocean (Kang
et al., 2021; Zeng & Jansen, 2021).

##  6 Conclusions

Consideration of the sources and dissipation rates of kinetic energy provides
useful constraints for the circulation of ice-covered oceans. In general,
kinetic energy can be generated by tides or by conversion from potential
energy, which in turn can be generated by heat and salt flux forcing.

A commonly assumed forcing consists of heating from the solid core balanced by
heat loss through the ice sheet, which acts as a source of potential energy
and can drive an ocean circulation, as long as the thermal expansivity of the
water is positive. In the oceans of Snowball Earth, Europa, and Enceladus, the
associated energy input, however, is orders of magnitude smaller than the wind
energy input that dominantly drives the circulation in Earth’s present-day
ocean. We predict that the resulting thermally-driven flows have flow speeds
of at most a few cm/s and will be strongly affected by rotation (i.e small
Rossby numbers). Numerical simulations of thermally-driven flows on icy moons,
which typically use artificially large viscous dissipation and sometimes
artificially large thermal forcing may misrepresent both energy sources and
sinks by multiple orders of magnitude, which can lead to widely different and
unrealistic levels of kinetic energy.

Salt fluxes associated with freezing and melting at the ice sheet boundary
only provides a source of potential energy if freezing occurs under thinner
ice than melting, which is unlikely in equilibrium where freezing and melting
needs to be in balance with ice flow. In the more likely scenario where
melting on average occurs at a shallower depth than freezing, the salt flux
forcing acts as a sink of potential energy. Turbulent vertical mixing, which
can push the lighter fresh water into the interior, is then needed to drive a
significant circulation. Such turbulent mixing, however, requires an
alternative energy source.

Ocean tides and librations may provide a key energy source for ocean
turbulence. Current estimates indicate that tidally-driven vertical mixing is
likely to be important in a Snowball Earth ocean, and could possibly play a
significant role on Europa, while librations may provide a key source of
turbulent kinetic energy on Enceladus. However, the magnitude and spatial
distribution of turbulence generated by tides and librations remains highly
uncertain, which represents a major hurdle to better constrain the
circulations of icy-world oceans. An improved understanding of ocean tides and
librations on icy moons thus remains as an important topic for future
research.

##  Appendix A A generalized expression for energy input from buoyancy forcing

To quantify the role of heat and salinity forcing on the vertical buoyancy
flux, while fully accounting for nonlinearities in the equation of state, it
is useful to introduce the dynamic enthalpy (Young, 2010):

$h^{\ddagger}(\Theta,S,z)=\int_{z}^{0}\tilde{b}(\Theta,S,z^{\prime})dz^{\prime}$
(A1)

The dynamic enthalpy, $h^{\ddagger}$, is closely related to the potential
energy, and indeed reduces to the well known expression for the potential
energy in the limit of an equation of state with no explicit depth dependence
(i.e. $b=\tilde{b}(\Theta,S)$) where $h^{\ddagger}\rightarrow bz$.

The evolution equation for $h^{\ddagger}$ is

$\frac{Dh^{\ddagger}}{Dt}=-wb+\dot{\Theta}\partial_{\Theta}h^{\ddagger}+\dot{S}\partial_{S}h^{\ddagger}$
(A2)

where $\dot{\Theta}\equiv D\theta/Dt$ and $\dot{S}\equiv DS/Dt$ are given by
Eqs. (6) and (7).

Integrating globally and assuming the global enthalpy budget to be in
equilibrium, we can relate the globally integrated vertical advective buoyancy
flux, $\mathcal{B}$, which provides a source of KE, to the thermal and
salinitiy forcing, $\dot{\Theta}$ and $\dot{S}$:

$\mathcal{B}\equiv\int
wbdV=\int\left(\dot{\Theta}\partial_{\Theta}h^{\ddagger}+\dot{S}\partial_{S}h^{\ddagger}\right)dV$
(A3)

Defining

$\widehat{g\alpha}(\Theta,S,z)\equiv\frac{1}{-z}\partial_{\Theta}h^{\ddagger}=\frac{1}{-z}\int_{z}^{0}g\alpha(\Theta,S,z^{\prime})dz^{\prime}$
(A4)

and

$\widehat{g\beta}(\Theta,S,z)\equiv\frac{1}{-z}\partial_{S}h^{\ddagger}=\frac{1}{-z}\int_{z}^{0}g\beta(\Theta,S,z^{\prime})dz^{\prime}$
(A5)

we can write Eq. (A3) as

$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle-\int\left(\dot{\Theta}\widehat{g\alpha}z+\dot{S}\widehat{g\beta}z\right)dV$
(A6) $\displaystyle=$
$\displaystyle-\int\kappa\left(\nabla^{2}\Theta\,\widehat{g\alpha}z+\nabla^{2}S\,\widehat{g\beta}z\right)dV$
$\displaystyle=$
$\displaystyle-\int\kappa\left(\nabla\cdot(\nabla\Theta\,\widehat{g\alpha}z)-\nabla\Theta\cdot\nabla(\widehat{g\alpha}z)+\nabla\cdot(\nabla
S\,\widehat{g\beta}z)-\nabla S\cdot\nabla(\widehat{g\beta}z)\right)dV$
$\displaystyle=$ $\displaystyle\int\frac{\mathcal{Q}_{top}}{\rho
c_{p}}\widehat{g\alpha}z_{top}dA-\int\frac{\mathcal{Q}_{bot}}{\rho
c_{p}}\widehat{g\alpha}z_{bot}dA$
$\displaystyle+\int\mathcal{S}_{top}\widehat{g\beta}z_{top}dA-\int\mathcal{S}_{bot}\widehat{g\beta}z_{bot}dA$
$\displaystyle+\int
g\kappa\left(\alpha\partial_{z}\Theta+\beta\partial_{z}S\right)dV$
$\displaystyle+\int
z\kappa\left(\nabla_{h}(\widehat{g\alpha})\cdot\nabla_{h}\Theta+\nabla_{h}(\widehat{g\beta})\cdot\nabla_{h}S\right)dV$

Eq. (A6) allows us to relate the globally integrated advective upward buoyancy
flux (and hence the associated KE generation) to the heat and salt fluxes
through the boundaries (the first four terms on the RHS of Eq. A6), reduced by
the diffusive buoyancy flux (the fifth term on the RHS of Eq. A6). The last
term in Eq. (A6) captures the effect of buoyancy sources or sinks that arise
from horizontal diffusive mixing of $\Theta$ and $S$ due to the nonlinearity
in the equation of state (the cabbeling effect).

The potential role of the cabbeling term in Eq. (A6) can be illustrated by
considering again the example of section 3.3, which may resemble the
conditions on Enceladus, assuming a relatively low salinity for the latter.
That is we assume a stratified layer with negative thermal expansion
coefficient above a convective deep layer with weakly positive thermal
expansion coefficient. For illustrative purposes, let us assume the upper and
lower boundaries to be horizontal, bottom and surface heat fluxes to be
spatially uniform, and we ignore radial variations in the surface area, such
that $Q_{top}=Q_{bot}\equiv Q$, as well as vertical variations in gravity.
Since variations in the thermal expansion coefficient are primarily due to
differences in the temperature, we moreover neglect the direct depth
dependence of $\tilde{b}$ and assume that variations in salinity are small,
such that $\widehat{g\alpha}=g\alpha(\theta)$. In this case, Eq. (A6)
simplifies to

$\displaystyle\mathcal{B}$ $\displaystyle=$ $\displaystyle\frac{gQ}{\rho
c_{p}}\int\left(\alpha(\Theta_{top})z_{top}-\alpha(\Theta_{bot})z_{bot}\right)dA$
$\displaystyle+\int g\kappa\alpha(\Theta)\partial_{z}\Theta dV+\int
gz\kappa\partial_{\Theta}{\tilde{\alpha}}|\nabla_{h}\Theta|^{2}dV\,.$

where we used that
$\nabla_{h}\alpha\cdot\nabla_{h}\Theta=\partial_{\Theta}{\alpha}|\nabla_{h}\Theta|^{2}$.
The last term can now be associated with the diffusive destruction of
temperature variance, which leads to a buoyancy source if
$\partial_{\Theta}{\alpha}>0$ (amounting to a positive curvature in
$b(\theta)$).

Notice that the magnitude (and potentially the sign) of the first and last
term in Eq. (A) depend of the choice of reference level. E.g. setting
$z_{bot}=0$ and $z_{top}=H$ where $H$ is the depth of the ocean (which we
assume to be constant with flat upper and lower boundaries) we find

$\displaystyle\mathcal{B}$ $\displaystyle=$ $\displaystyle\frac{gQH}{\rho
c_{p}}\int\alpha(\Theta_{top})dA$ (A7) $\displaystyle+\int
g\kappa\alpha(\Theta)\partial_{z}\Theta
dV+\int_{0}^{H}gz\kappa\partial_{\Theta}{\alpha}|\nabla\Theta|^{2}dV\,.$

Defining the reference level instead such that $z_{top}=0$ and $z_{bot}=-H$ we
get

$\displaystyle\mathcal{B}$ $\displaystyle=$ $\displaystyle\frac{gQH}{\rho
c_{p}}\int\alpha(\Theta_{bot})dA$ (A8) $\displaystyle+\int
g\kappa\alpha(\Theta)\partial_{z}\Theta
dV+\int_{-H}^{0}gz\kappa\partial_{\Theta}{\alpha}|\nabla\Theta|^{2}dV\,.$

The first term differs dramatically in the two equations, as
$\alpha(\Theta_{top})<0$ while $\alpha(\Theta_{bot})>0$. This difference is
compensated by the last term, which also changes sign, as $z>0$ in (A7) and
$z<0$ in (A8). The last term in Eq. (A6) hence can be of leading order
importance if variations in the thermal (or haline) expansion coefficient are
large. Since this term cannot be predicted based on knowledge of the boundary
conditions alone, this most general formulation is hence less immediately
useful as a predictor of flow energetics. However, if horizontal temperature
variations are relatively small, Eq. (14) remains a useful approximation and
allows us to estimate the kinetic energy input without explicit consideration
of diffusive variance destruction.

## References

  * Adcroft et al. (2019) Adcroft, A., Anderson, W., Blanton, C., et al. 2019, JAMES
  * Ashkenazy et al. (2013) Ashkenazy, Y., Gildor, H., Losch, M., et al. 2013, Nature, 495, 90
  * Ashkenazy & Tziperman (2016) Ashkenazy, Y., & Tziperman, E. 2016, Journal of Climate, 29, 869
  * Ashkenazy & Tziperman (2021) —. 2021, Nature communications, 12, 1
  * Beuthe (2016) Beuthe, M. 2016, Icarus, 280, 278
  * Bire et al. (2022) Bire, S., Kang, W., Ramadhan, A., Campin, J.-M., & Marshall, J. 2022, Journal of Geophysical Research: Planets, e2021JE007025
  * Brenner et al. (2021) Brenner, S., Rainville, L., Thomson, J., Cole, S., & Lee, C. 2021, Journal of Geophysical Research: Oceans, 126, e2020JC016977
  * Brink (2011) Brink, K. H. 2011, Journal of Marine Research, 69, 483
  * Carr et al. (1998) Carr, M. H., Belton, M. J., Chapman, C. R., et al. 1998, Nature, 391, 363
  * Chen et al. (2014) Chen, E., Nimmo, F., & Glatzmaier, G. 2014, Icarus, 229, 11
  * Chen et al. (2018) Chen, Y., Zhang, F., Green, B. W., & Yu, X. 2018, Monthly Weather Review, 146, 287
  * Choblet et al. (2017) Choblet, G., Tobie, G., Sotin, C., et al. 2017, Nature Astronomy, 1, 841
  * Des Marais et al. (2008) Des Marais, D. J., Nuth III, J. A., Allamandola, L. J., et al. 2008, Astrobiology, 8, 715
  * Egbert et al. (2004) Egbert, G. D., Ray, R. D., & Bills, B. G. 2004, Journal of Geophysical Research: Oceans, 109
  * Fernando et al. (1991) Fernando, H. J., Chen, R.-R., & Boyer, D. L. 1991, Journal of Fluid Mechanics, 228, 513
  * Gissinger & Petitdemange (2019) Gissinger, C., & Petitdemange, L. 2019, Nature Astronomy, 3, 401
  * Green et al. (2017) Green, J., Huber, M., Waltham, D., Buzan, J., & Wells, M. 2017, Earth and Planetary Science Letters, 461, 46
  * Gutmark & Wygnanski (1976) Gutmark, E., & Wygnanski, I. 1976, Journal of fluid Mechanics, 73, 465
  * Hay & Matsuyama (2019) Hay, H. C., & Matsuyama, I. 2019, Icarus, 319, 68
  * Hemingway et al. (2018) Hemingway, D., Iess, L., Tajeddine, R., & Tobie, G. 2018, Enceladus and the icy moons of Saturn, 57
  * Hendrix et al. (2019) Hendrix, A. R., Hurford, T. A., Barge, L. M., et al. 2019, Astrobiology, 19, 1
  * Hoffman et al. (2017) Hoffman, P. F., Abbot, D. S., Ashkenazy, Y., et al. 2017, Science Advances, 3, e1600983
  * Hussmann & Spohn (2004) Hussmann, H., & Spohn, T. 2004, Icarus, 171, 391
  * Huthnance (1981) Huthnance, J. 1981, Journal of Fluid Mechanics, 102, 367
  * Jansen (2016) Jansen, M. F. 2016, Journal of Physical Oceanography, 46, 1917
  * Jones & Marshall (1993) Jones, H., & Marshall, J. 1993, Journal of Physical Oceanography, 23, 1009
  * Kang et al. (2020) Kang, W., Bire, S., Campin, J.-M., et al. 2020, arXiv preprint arXiv:2008.03764
  * Kang et al. (2021) Kang, W., Mittal, T., Bire, S., Campin, J.-M., & Marshall, J. 2021, arXiv preprint arXiv:2104.07008
  * Kite & Rubin (2016) Kite, E. S., & Rubin, A. M. 2016, Proceedings of the National Academy of Sciences, 113, 3972
  * Kivelson et al. (2000) Kivelson, M. G., Khurana, K. K., Russell, C. T., et al. 2000, Science, 289, 1340
  * Lemasquerier et al. (2017) Lemasquerier, D., Grannan, A., Vidal, J., et al. 2017, Journal of Geophysical Research: Planets, 122, 1926
  * Lobo et al. (2021) Lobo, A. H., Thompson, A. F., Vance, S. D., & Tharimena, S. 2021, Nature Geoscience, 14, 185
  * Matsuyama et al. (2018) Matsuyama, I., Beuthe, M., Hay, H. C., Nimmo, F., & Kamata, S. 2018, Icarus, 312, 208
  * Maxworthy & Narimousa (1994) Maxworthy, T., & Narimousa, S. 1994, J. Phys. Oceanogr., 24, 865
  * Melosh et al. (2004) Melosh, H., Ekholm, A., Showman, A., & Lorenz, R. 2004, Icarus, 168, 498
  * Nakajima et al. (2019) Nakajima, A., Ida, S., Kimura, J., & Brasser, R. 2019, Icarus, 317, 570
  * Nimmo & Pappalardo (2016) Nimmo, F., & Pappalardo, R. 2016, Journal of Geophysical Research: Planets, 121, 1378
  * Peltier & Caulfield (2003) Peltier, W., & Caulfield, C. 2003, Annual review of fluid mechanics, 35, 135
  * Postberg et al. (2018) Postberg, F., Khawaja, N., Abel, B., et al. 2018, Nature, 558, 564
  * Rekier et al. (2019) Rekier, J., Trinh, A., Triana, S. A., & Dehant, V. 2019, Journal of Geophysical Research: Planets, 124, 2198
  * Ruiz (2005) Ruiz, J. 2005, Icarus, 177, 438
  * Sen et al. (2008) Sen, A., Scott, R. B., & Arbic, B. K. 2008, Geophysical Research Letters, 35
  * Smagorinsky et al. (1965) Smagorinsky, J., Manabe, S., & Holloway, J. L. 1965, Monthly weather review, 93, 727
  * Soderlund et al. (2014) Soderlund, K., Schmidt, B., Wicht, J., & Blankenship, D. 2014, Nature Geoscience, 7, 16
  * Soderlund (2019) Soderlund, K. M. 2019, Geophysical Research Letters, 46, 8700
  * Soderlund et al. (2020) Soderlund, K. M., Kalousová, K., Buffo, J. J., et al. 2020, Space Science Reviews, 216, 1
  * Thomas et al. (2016) Thomas, P., Tajeddine, R., Tiscareno, M., et al. 2016, Icarus, 264, 37
  * Tyler (2020) Tyler, R. H. 2020, Icarus, 113821
  * Vallis (2006) Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge, U.K.: Cambridge University Press), 745
  * Vic et al. (2019) Vic, C., Garabato, A. C. N., Green, J. M., et al. 2019, Nature communications, 10, 1
  * Waite et al. (2017) Waite, J. H., Glein, C. R., Perryman, R. S., et al. 2017, Science, 356, 155
  * Williams (2000) Williams, G. E. 2000, Reviews of Geophysics, 38, 37
  * Wilson & Kerswell (2018) Wilson, A., & Kerswell, R. R. 2018, Earth and Planetary Science Letters, 500, 41
  * Wunsch (2016) Wunsch, C. 2016, Icarus, 274, 122
  * Wunsch & Ferrari (2004) Wunsch, C., & Ferrari, R. 2004, Annu. Rev. Fluid Mech., 36, 281
  * Young (2010) Young, W. R. 2010, Journal of Physical Oceanography, 40, 394
  * Zeng & Jansen (2021) Zeng, Y., & Jansen, M. F. 2021, Planet. Sci. J., 2
  * Zhu et al. (2017) Zhu, P., Manucharyan, G. E., Thompson, A. F., Goodman, J. C., & Vance, S. D. 2017, Geophysical Research Letters, 44, 5969