# Resonant stratification in Titan’s global ocean

Benjamin Idini Department of Earth and Planetary Sciences, University of
California Santa Cruz
1156 High St
Santa Cruz, CA 95064, USA  Francis Nimmo Department of Earth and Planetary
Sciences, University of California Santa Cruz
1156 High St
Santa Cruz, CA 95064, USA

###### Abstract

Titan’s ice shell floats on top of a global ocean revealed by the large tidal
Love number $k_{2}=0.616\pm 0.067$ registered by Cassini. The Cassini
observation exceeds the predicted $k_{2}$ by one order of magnitude in the
absence of an ocean, and is 3-$\sigma$ away from the predicted $k_{2}$ if the
ocean is pure water resting on top of a rigid ocean floor. Previous studies
demonstrate that an ocean heavily enriched in salts ( $S\gtrsim 200$ g/kg) can
explain the 3-$\sigma$ signal in $k_{2}$. Here we revisit previous
interpretations of Titan’s large $k_{2}$ using simple physical arguments and
propose a new interpretation based on the dynamic tidal response of a stably
stratified ocean in resonance with eccentricity tides raised by Saturn. Our
models include inertial effects from a full consideration of the Coriolis
force and the radial stratification of the ocean, typically neglected or
approximated elsewhere. The stratification of the ocean emerges from a
salinity profile where salt concentration linearly increases with depth. We
find multiple salinity profiles that lead to the $k_{2}$ required by Cassini.
In contrast with previous interpretations that neglect stratification,
resonant stratification reduces the bulk salinity required by observations by
an order of magnitude, reaching a salinity for Titan’s ocean that is
compatible with that of Earth’s oceans and close to Enceladus’ plumes.
Consequently, no special process is required to enrich Titan’s ocean to a high
salinity as previously suggested.

††journal: PSJ††software: Mathematica (Wolfram, 1999), Matplotlib (Hunter,
2007)

##  1 Introduction

Recent decades of space exploration have revealed a solar system populated
with internally heated icy worlds where large reservoirs of liquid water
accumulate in subsurface global oceans (Nimmo & Pappalardo, 2016). These
worlds signal a possibility for life beyond Earth in a location that is
accessible to future in-situ space exploration. A global ocean plays a
fundamental role in determining the potential habitability of these icy worlds
because water is required for life as we know it. Beyond detection, however,
these global oceans remain poorly understood. Ocean thickness is typically
known within broad limits ranging in the tens of percent of the icy world
radius (Sohl et al., 2003; Grindrod et al., 2008), preventing an accurate
assessment of the satellite’s liquid water inventory and thermal history. On
Earth, ocean dynamics modulate the distribution of nutrients and energy
sources required by life (Uchida et al., 2020, e.g.), but on icy worlds the
type of convection or lack thereof remains poorly known (Jansen et al., 2023).
Here we argue in favor of the stratification of Titan’s ocean based on a new
interpretation of Cassini gravity measurements where internal gravity waves in
Titan’s ocean become resonantly excited by tides raised by Saturn (see also
Luan (2019)). We hereafter refer to this proposed scenario as resonant
stratification.

Titan is the second largest solar system icy satellite and the best
characterized from the perspective of gravity measurements (i.e., moment of
inertia, $J_{2}$, $C_{22}$, and $k_{2}$ (Durante et al., 2019)), offering us a
unique opportunity to reveal the hidden interior structure and dynamics of the
global ocean within. The Cassini spacecraft unambiguously signaled the
existence of a global ocean from the observed tidal response registered in the
Love number $k_{2}=0.616\pm 0.067$ (Durante et al., 2019). The observation is
an order of magnitude larger than the predicted $k_{2}\approx 0.03$ when the
ocean is absent (Rappaport et al., 2008). Ignoring the ice shell, a global
ocean of pure liquid water produces a $k_{2}\approx 0.468$ independent of
ocean thickness if the high-pressure ice and silicates beneath the ocean
behave approximately rigidly and the total mass of the satellite is conserved
(Section 2.1.2). The presence of an overlying elastic ice shell further
restricts the motion of the ocean beneath. Estimates of Titan’s ice shell
thickness yield $d\sim 100$ km (Sohl et al., 2003; Nimmo & Bills, 2010; Luan,
2019); an ice shell this thick reduces tides down to $k_{2}\approx 0.42$
(Section 2.1.3), which is roughly $3\sigma$ away from the Cassini observation.
A thinner ice shell provides a $k_{2}$ closer to the observation, but then the
heat conducted across the ice shell exceeds the interior heat production
expected from radiogenic and tidal heating (Sohl et al., 2003; Luan, 2019),
leading to thickening of the ice shell by freezing over time. A thinner shell
is also more difficult to reconcile with the observed topography (Nimmo &
Bills, 2010). The resonant stratification presented here can self-consistently
explain the large $k_{2}$ observed by Cassini by introducing a positive
dynamical fractional correction to the non-resonant hydrostatic $k_{2}\approx
0.42$.

Figure 1: Proposed explanations to Titan’s large $k_{2}$ as observed by
Cassini. The schematic of Titan’s interior model has been extracted from de
Kleer et al. (2019).

Resonant stratification enhances $k_{2}$ by dynamic amplification of the
vertical displacement of Titan’s surface (Fig. 1). Ocean waves can produce
significant dynamic surface displacement when resonantly excited by Saturn’s
eccentricity tides, namely when Titan’s orbital frequency
($\omega_{s}=4.561\times 10^{-6}$ s-1) is close to a match with a normal mode
emerging from an aggregation of Titan’s ocean waves. Amplification results
from the constructive addition of ocean waves over many cycles of resonant
excitation, until balanced by dissipation. The effect produces dynamic gravity
that can be registered by the tracking system of a nearby passing spacecraft
(i.e., Cassini). A vertical gradient in ocean salinity promotes ocean
stratification and the emergence of internal gravity waves. These waves are
restored by buoyancy forces and organize in a spectrum of low-frequency normal
modes that asymptotically approach zero frequency from an upper bound
frequency $\omega^{2}\lesssim N^{2}$, where $N^{2}$ is the Brunt-Vaisala
frequency that typically describes the strength of ocean stratification
(Section 2.2 and Appendix A). A typical value $N^{2}\gtrsim 10^{-8}$
s${}^{-2}\gtrsim\omega_{s}^{2}$ suggests a spectrum of internal gravity waves
capable of resonance with eccentricity tides, and roughly corresponds to a
modest vertical salinity gradient where salt concentration increases by
${\gtrsim}+1$ g/kg every ${\sim}100$ km of ocean thickness (the salinity of
Earth’s oceans is roughly $35$ g/kg). Crucially, our mechanism works for a
range of mean ocean salinity values, including Earth’s and that inferred for
Enceladus (Postberg et al., 2009).

When the ocean is fully convective (i.e., unstratified) and homogeneously
mixed, ocean waves require an unphysically thin ocean ($H<1$ km) to resonate
with eccentricity tides (Matsuyama et al., 2018). An ocean this thin is not
compatible with the thermal history of Titan as expected from radiogenic and
tidal heating (Tobie et al., 2006; Grindrod et al., 2008). Only the high-
frequency tides from moon-to-moon interactions can resonantly excite a fully
mixed ocean of realistic thickness ($H\sim 100$ km)(Hay et al., 2020), but
those signals are relatively small and thus beyond Cassini’s detection
threshold.

Alternative scenarios can enhance $k_{2}$ to the value required by Cassini
(Fig. 1), but require either high salinity of the ocean or low rigidity of the
solid Titan. Previous calculations (Rappaport et al., 2008; Waite et al.,
2017) show that a heavy ocean produces the required extra gravity when its
density is increased by a high concentration of dissolved salts ( $S\sim 200$
g/kg). The required concentration of salts is an order of magnitude higher
than that in Earth’s oceans, Enceladus’ plumes, and that predicted from water-
rock interactions on Titan’s ocean floor (Postberg et al., 2009; Leitner &
Lunine, 2019). On the other hand, a low rigidity or low viscosity ocean floor
can allow large vertical displacements of the ocean floor that produce an
extra gravity signal that constructively adds to the gravity from tides at the
surface (Durante et al., 2019). In practice, reasonable estimates of the
rigidity and viscosity for silicates and high-pressure ice at tidal timescales
produce a negligible ocean floor displacement. For example, the contribution
to $k_{2}$ from the elastic icy ocean floor displacement can be estimated to
be negligible following
$k_{2,\textnormal{ice}}/k_{2}\times(\rho_{\textnormal{hp-ice}}-\rho)/\rho\sim
2\%$, where $k_{2,\textnormal{ice}}/k_{2}$ is the ratio between the Love
number when the ocean is excluded/included, respectively,
$\rho_{\textnormal{hp-ice}}$ is the density of high-pressure ice, and $\rho$
the ocean density. A viscous icy ocean floor can introduce a non-negligible
component to $k_{2}$ if the viscosity of high-pressure ice is lower than
${\sim}10^{15}$ Pa s (Rappaport et al., 2008). The viscosity of high-pressure
ice is typically comparable to or greater than this value when the temperature
of the ice is ${\sim}10\%$ lower than the melting temperature (Durham et al.,
1998), a reasonable assumption at the ocean floor.

Since resonantly excited internal gravity waves can produce the required
dynamic gravity, the main challenge for resonant stratification is how to
establish a long-lived tidal resonance in an icy satellite. The possibility of
resonant tidal excitation is suggested by the tendency of Titan’s ocean to
freeze due to secular cooling, which imposes a continuous evolution on the
frequency of ocean waves by a change in the stratification profile of the
ocean. At some point in Titan’s freezing history, the normal frequency of
ocean waves crosses the orbital frequency. After first being encountered, a
tidal resonance can be sustained over long geological timescales provided that
a stable fixed point is reached between orbital and interior evolution,
analogously to ideas previously applied to tidal resonant locking in Jupiter
and Saturn (Fuller et al., 2016; Luan et al., 2018; Lainey et al., 2020; Idini
& Stevenson, 2022). The onset of resonant ocean waves enhances Titan’s tidal
heating in the ocean (Tyler, 2011, 2014; Rovira-Navarro et al., 2019), with
the additional heating slowing or halting the freezing of the ocean.

##  2 Basic Equations

The tidal Love number $k_{\ell m}$ represents a normalization of the
gravitational tidal response by the gravitational excitation,

$k_{\ell m}=\left(\frac{\phi^{\prime}_{\ell m}}{\phi^{e}_{\ell
m}}\right)_{r=R}\textnormal{,}$ (1)

where the eccentric gravitational excitation $\phi^{e}_{\ell m}$ is derived in
Appendix B and the tidal response $\phi^{\prime}_{\ell m}$ due to eccentricity
tides is calculated numerically and analytically in this section in the
context of various models.

###  2.1 Equilibrium tide

####  2.1.1 Hydrostatic tides in a homogeneous fluid body

We first revisit the classical problem of tides in a homogeneous
incompressible body that satisfies hydrostatic equilibrium using basic
principles. The linearized equation for the conservation of momentum is

$\frac{p^{\prime}}{\rho}=\tilde{\phi}^{\prime}\textnormal{,}$ (2)

The tidal response of the icy satellite produces adiabatic perturbations
represented in the gravitational potential $\phi^{\prime}$ and pressure
$p^{\prime}$. In this expression, the potential of the gravitational pull
$\phi^{e}$ and the tidal gravitational potential $\phi^{\prime}$ are combined
into $\tilde{\phi}^{\prime}=\phi^{e}+\phi^{\prime}$ for analytical simplicity.

In an incompressible body, the density $\rho$ remains constant regardless of
forcing and the only change in the gravity field comes about from the radial
displacement $\xi_{r}$ of the surface. This diplacement is typically small
compared to the radius $R$. Thus, we can calculate the gravity field from the
integration over the volume of a thin shell of fluid with thickness $\xi_{r}$,

$\phi^{\prime}_{\ell
m}=\frac{4\pi\mathcal{G}\rho}{(2\ell+1)}\frac{R^{\ell+2}}{r^{\ell+1}}\xi_{r\ell
m}=\frac{3}{(2\ell+1)}\left(\frac{R}{r}\right)^{\ell+1}g\xi_{r\ell
m}\textnormal{,}$ (3)

where $g$ is the surface gravity acceleration calculated in hydrostatic
equilibrium $g=4\pi\mathcal{G}\rho R/3$, and $\mathcal{G}$ the gravitational
constant.

The boundary condition on the surface indicates that the fluid is free from
external pressure (i.e., $\delta p=0$). This statement translates into turning
the tidal displacement into the sole source for the pressure perturbation,
formally expressed as

$\frac{p^{\prime}}{\rho}=-\frac{\xi}{\rho}\frac{\partial p}{\partial
r}=g\xi_{r}\textnormal{.}$ (4)

Next, we can go back to equation (2) and obtain an expression for the tidal
forcing,

$\phi^{e}_{\ell m}=\frac{p^{\prime}_{\ell m}}{\rho}-\phi^{\prime}_{\ell
m}=\left(1-\frac{3}{2\ell+1}\left(\frac{R}{r}\right)^{\ell+1}\right)g\xi_{r\ell
m}\textnormal{.}$ (5)

We evaluate the tidal response and forcing at the surface ($r=R$) to obtain
the tidal Love number from the ratio between them,

$k_{\ell m}=\frac{3}{2(\ell-1)}\textnormal{,}$ (6)

which agrees with the classical result $k_{2}=1.5$ of the fluid Love number of
a uniform density body in hydrostatic equilibrium (Munk et al., 1977, e.g.).

####  2.1.2 A global ocean with a rigid ocean floor

We now consider a body with a rigid core of density $\rho_{c}$ and radius
$R_{c}$ overlaid by an ocean of density $\rho$ by introducing vertical
differentiation in the density profile. The main change compared to the
homogeneous case is a change in the expression for the surface gravity
acceleration $g=4\pi\mathcal{G}\bar{\rho}R/3$, where $\bar{\rho}$ is the
body’s mean density. The tidal gravity field still emerges solely from the
radial displacement $\xi_{r\ell m}$ of a thin region near the surface,

$\phi^{\prime}_{\ell
m}=\frac{3}{(2\ell+1)}\left(\frac{\rho}{\bar{\rho}}\right)\left(\frac{R}{r}\right)^{\ell+1}g\xi_{r\ell
m}\textnormal{.}$ (7)

Following the same steps as before, the resulting tidal Love number is
independent of core properties,

$k_{\ell m}=\frac{3}{\left(2\ell+1\right)\bar{\rho}/\rho-3}\textnormal{.}$ (8)

This equation yields $k_{2}=0.468$ when using Titan’s mean density
$\bar{\rho}\approx 1.88$ g cm-3 and $\rho=1$ g cm-3. Equation (8) reproduces
numerical results (within 5%) obtained previously for Europa when the icy
shell is ignored (Moore & Schubert, 2000), namely $k_{2}\approx 0.249$ when
using $\bar{\rho}\approx 3.01$ g cm-3. Analogous versions of Equation (8) can
be found in the literature when following an alternative derivation (Dermott,
1979; Murray & Dermott, 1999; Beuthe, 2015, e.g.).

####  2.1.3 Ice shell thickness

An elastic ice shell constrains the radial displacement of the ocean that
produces the tidal response registered in $k_{2}$. The weight of the ocean
trying to reach an equipotential surface is balanced by the resistance of the
ice shell to deformation (Kamata et al., 2015, e.g.). Here we provide a simple
argument to evaluate the role that the ice shell thickness plays in reducing
the amplitude of the tidal response.

Consider a global ocean of density $\rho$ that is trying to reach the
equipotential surface defined by the radial displacement $\xi_{req}$. The
pressure that the weight of the ocean exerts on the base of the ice shell can
be expressed by

$p\sim\rho g(\xi_{req}-\xi_{r})\textnormal{,}$ (9)

where $g$ is the gravity acceleration and $\xi_{r}$ is the radial tidal
displacement. The shear stress $\tau$ on the ice shell can be expressed as

$\tau\sim\frac{\mu\xi_{r}}{R}\textnormal{,}$ (10)

where $\mu$ is the shear modulus of the ice shell.

Next, we consider the equilibrium of forces over a meridional section dividing
the satellite into two hemispheres. The ocean pressure integrated over the
projected area of the ice shell base must be balanced by the shear stress
integrated over the section of the ice shell, namely

$\pi R^{2}\rho g(\xi_{req}-\xi_{r})\sim 2\pi d\mu\xi_{r}\textnormal{,}$ (11)

where $d$ is the ice shell thickness and $R$ is the icy satellite radius. We
can rearrange the previous expression to discover the fractional change in
$k_{2}$ due to a rigid ice shell,

$\Delta k_{2}\sim\frac{\xi_{r}}{\xi_{req}}-1\sim-\frac{2d\mu}{R^{2}\rho
g}\textnormal{,}$ (12)

which represents a competition between elastic and gravitational energy (see
also equation (11) in Goldreich & Mitchell (2010)). We may use Titan’s radius
$R\approx 2575$ km, surface gravity acceleration $g\approx 135$ cm s-2, ocean
density $\rho\approx 1$ g cm-3, and shear modulus of ice I and methane
clathrates $\mu\sim 4$ GPa, and obtain

$\Delta k_{2}\sim-10\textnormal{\%}\left(\frac{d}{100\textnormal{
km}}\right)\textnormal{.}$ (13)

An ice shell thickness $d\sim 100$ km balances radiogenic heating with heat
conduction (Appendix D) and produces $k_{2}\approx 0.42$ (the hydrostatic
response without an ice shell is $k_{2}=0.468$). We use this $k_{2}$ value as
the reference hydrostatic tidal response when computing the fractional change
$\Delta k_{2}$ due to various effects, which is compatible with previous
estimates for Titan models without salinity (Rappaport et al., 2008). The
result is valid for a shell made of ice or methane clathrates given that the
elastic modulus is similar in both cases.

###  2.2 Dynamical tides

We now solve the problem of a rotating ocean world with a stratified ocean
subjected to the full action of the Coriolis effect.

$\partial_{t}\bm{v}+2{\bf\Omega}\times\bm{v}=\nabla\tilde{\phi}-\frac{\nabla
p}{\rho}-\gamma\bm{v}\textnormal{,}$ (14)
$\nabla\cdot{\bm{v}}=0\textnormal{,}$ (15)

where $\bm{v}$ is the tidal flow, $\bm{\Omega}$ is the icy satellite spin
vector, $\rho$ is the ocean density, $p$ is pressure, and $\gamma$ is the
linear damping coefficient that represents dissipation in the ocean. In this
paper, we consider at all times that Titan is in synchronous rotation (i.e.
$\Omega=\omega_{s}$).

The resulting linearized equation of conservation of momentum is

$\partial_{t}\bm{v}+2{\bf\Omega}\times\bm{v}=\nabla\tilde{\phi}^{\prime}-\frac{\nabla
p^{\prime}}{\rho}+\frac{\rho^{\prime}}{\rho^{2}}\nabla
p-\gamma\bm{v}\textnormal{.}$ (16)

The timescale of tidal motion of ${\sim}15$ days is orders of magnitude
shorter than the timescale required to transport heat by diffusion or
convection. Tides are consequently adiabatic and the associated Lagrangian
tidal perturbations in density and pressure satisfy

$\frac{\delta p}{p}=\Gamma\frac{\delta\rho}{\rho}\textnormal{,}$ (17)

where $\Gamma$ is the adiabatic index. We can relate the Lagrangian
perturbations to the Eulerian perturbations in equation (16) using the
definitions

$\delta p=p^{\prime}+{\bm{\xi}}\cdot\nabla p\textnormal{,}$ (18)
$\delta\rho=\rho^{\prime}+{\bm{\xi}}\cdot\nabla\rho\textnormal{,}$ (19)

where ${\bm{\xi}}$ is the tidal displacement. Moving forward, we put together
the latter equations to arrive to

$\frac{\rho^{\prime}}{\rho}=\frac{p^{\prime}}{\Gamma
p}+\bm{\xi}\cdot\frac{\bm{N}^{2}}{g}\textnormal{,}$ (20)

where $N^{2}$ is the Brunt-Vaisala (BV) frequency defined in general by

$\bm{N}^{2}=g\left(\frac{\nabla p}{\Gamma
p}-\frac{\nabla\rho}{\rho}\right)\textnormal{.}$ (21)

The BV frequency represents the salinity stratification of the ocean. Here the
only relevant component in $N^{2}$ is the radial direction, given that we
consider the background state to be spherically symmetric. When $\hat{\bf
r}N^{2}>0$, the ocean is vertically stratified. A stratified ocean parcel will
develop static stability once displaced out of its equilibrium position.

According to equation (20), an Eulerian change in ocean density may emerge
from either the adiabatic response of the ocean fluid (first term in the
right-hand side) or the buoyancy of the stratified ocean parcel (second term
in the right-hand side). In the incompressible approximation used here,
Eulerian perturbations to density uniquely emerge from the buoyancy of the
stratified ocean parcel. This results from the adiabatic index
$\Gamma\equiv(\partial\log p/\partial\log\rho)_{S}$ tending to infinity: the
expected changes in pressure keep density unperturbed. As a consequence, the
BV frequency reduces to

$\bm{N}^{2}=-g\frac{\partial_{r}\rho}{\rho}\hat{\bf r}\textnormal{,}$ (22)

where the radial variation in density comes from the addition of a small
amount of salts organized in a vertical salinity gradient. As we can see from
equation (22), stratification permits local density perturbations in the ocean
regardless of its incompressibility.

Equation (22) shows a direct relationship between stratification in $N^{2}$
and the salinity gradient in $\partial_{r}\rho$. When we consider constant
stratification throughout the ocean, we can write
$\partial_{r}\rho=\Delta\rho/H$, where $\Delta\rho$ is the change in density
between top and bottom of the ocean. Our models have constant density, thus
the $\Delta\rho$ represents a virtual change in density due to the addition of
salts. The concentration of added salts $S$ is described in g of salts per kg
of water (g/kg).

Next, we rewrite the linearized momentum and continuity equations as

$\omega\tilde{\omega}{\bm{\xi}}+i2\omega{\bf\Omega}\times{\bm{\xi}}-\left(\bm{N}^{2}\cdot{\bm{\xi}}\right)\hat{\bf
r}=-\nabla\psi\textnormal{,}$ (23)
$\psi=\tilde{\phi}^{\prime}-\frac{p^{\prime}}{\rho}\textnormal{,}$ (24)
$\nabla\cdot{\bm{\xi}}=0\textnormal{,}$ (25)

where $\tilde{\omega}=\omega+i\gamma$ is a complex frequency that accounts for
tidal dissipation, $\omega=\omega_{s}$ is the tidal frequency of eccentricity
tides, and the tidal displacement $\bm{\xi}$ is periodic with a time
dependency $\propto e^{-i\omega t}$. . This set of equations is traditionally
known as the Boussinesq approximation.

A typical method of solution of equation (23) involves applying the curl to
reach the vorticity equation (Rieutord, 1987)

$\nabla\times\left(2\bm{\Omega}\times\bm{\xi}\right)+i\tilde{\omega}\nabla\times\bm{\xi}+\frac{1}{i\omega}\nabla\times\left(\left(\bm{N}^{2}\cdot{\bm{\xi}}\right)\hat{\bf
r}\right)=0\textnormal{.}$ (26)

A rigid core implies that waves produce no radial displacement near the bottom
of the ocean. We set the rigid boundary condition at the ocean bottom to

$\bm{\xi}\cdot\hat{\bf r}=0\textnormal{.}$ (27)

The free boundary at the ocean top prescribes a zero Lagrangian perturbation
of pressure (e.g., Goodman & Lackner (2009)),

$\delta p=p^{\prime}+\bm{\xi}\cdot\nabla p=0\textnormal{,}$ (28)

which leads to

$\left(g-\frac{4\pi\mathcal{G}\rho R}{2\ell+1}\right)\bm{\xi}\cdot\hat{\bf
r}+\psi=\phi^{e}_{\ell m}\textnormal{.}$ (29)

The tidal displacement field $\bm{\xi}$ is the only quantity to be determined,
which is forced by the ocean top boundary condition on the radial component of
the displacement. Previous work (Rovira-Navarro et al., 2019, e.g.) typically
sets $\bm{\xi}\cdot\hat{\bf r}$ at the surface to be equal to the equilibrium
tide displacement (i.e., no-slip boundary conditions) and/or neglects the self
gravitation of the tide (e.g., Cowling approximation). Given that we are
interested in precise tidal gravity estimates, we have retained the self
gravitation of the equilibrium tide and allowed the ocean top to displace
dynamically beyond the equilibrium tide, . No-slip boundary conditions
explicitly set the dynamical part of $\bm{\xi}\cdot\hat{\bf r}$ to zero at the
surface, removing by definition the dynamical enhancement of $k_{2}$ that we
calculate here.

The resulting set of equations is an infinite system of equations coupled in
degree by the Coriolis effect (Rieutord, 1987; Rieutord & Valdettaro, 1997;
Lockitch & Friedman, 1999; Rieutord et al., 2001; Ogilvie & Lin, 2004;
Ogilvie, 2009; Rovira-Navarro et al., 2019; Idini & Stevenson, 2021, e.g.). We
solve this system by the traditional method of projection onto vectorial
spherical harmonics (Appendix C), followed by a pseudo-spectral discretization
on the radial functions based on the analytically-tractable Chebyshev
polynomials and Gauss–Lobatto collocation points (Boyd, 2001, e.g.). We
truncate the infinite set of equations at degree $\ell_{max}=100$ and use
$N_{max}=100$ Chebyshev polynomials to represent each radial function at each
degree $\ell\leq\ell_{max}$. The Love number can then be obtained from the
tidal displacement evaluated at the surface (equation (3)), which is the only
source of tidal gravity in a homogeneous ocean.

###  2.3 Tidal heating rate in the ocean

The heating rate in the ocean can be fully determined by volume integration

$dw=-\rho\gamma\bm{v}^{*}\cdot d\bm{l}\textnormal{,}$ (30)
$\dot{E}=\rho\gamma\omega^{2}\int_{V}\bm{\xi}^{*}\cdot\bm{\xi}dV\textnormal{.}$
(31)

The linear damping coefficient $\gamma$ is unconstrained in icy worlds with
global oceans. The range of possible estimates spans from $\gamma\sim
10^{-11}$ s-1 on Enceladus to the better constrained $\gamma\sim 10^{-5}$ s-1
on Earth (Matsuyama et al., 2018). The projection of equation (31) onto
vectorial spherical harmonics is shown in Appendix C.

##  3 Numerical results

###  3.1 The predicted tidal response $k_{2}$ as a function of ocean structure

We use perturbation theory and numerical methods to calculate the dynamic
gravity produced by resonant stratification in a stratified and rotating ocean
(Section 2). We concentrate on fractional changes $\Delta k_{2}$ to the non-
resonant hydrostatic $k_{2}\approx 0.42$ obtained in a pure water ocean with a
$d\sim 100$ km ice shell resting on top. Our method of solution considers the
Coriolis force in full and avoids the thin shell approximation typically used
in tidal studies of icy satellites (Beuthe, 2016; Matsuyama et al., 2018;
Rovira-Navarro et al., 2023). The additional effort of relaxing the thin-shell
approximation allows us to study the dynamic gravity of internal gravity waves
that result from the mixing of rotational and stratification effects. We
simplify the structure of stratification by assuming a constant $N^{2}$
throughout the ocean (equation (22)), which translates into a linear increase
in salt concentration with depth starting from zero at the ocean surface. More
complicated salinity distributions are in principle possible and lead to
additional uncertainty in the inferred ocean salinity structure.

Figure 2: Titan’s $k_{2}$ enhancement (fractional correction $\Delta k_{2}$)
from dynamic gravity as a function of (a) ocean thickness $H$ and (b) the
Brunt-Vaisala frequency $N^{2}$. Peaks represent resonances with ocean normal
modes. The ocean top salinity is zero and increases linearly with depth. The
ocean bottom salinity ranges $0.8-1.5$ g/kg ($H=100$ km) and $2.3-4.4$ g/kg
($H=300$ km) for the range of $N^{2}$ shown in (b). Frictional dissipation in
(b) is $\gamma=10^{-9}~{}\rm s^{-1}$. The reference hydrostatic Love number is
$k_{2}=0.42$, a pure water ocean with an elastic ice shell with thickness
$d\sim 100$ km. The tidal frequency is equal to the rotational frequency and
the orbital frequency, $\omega=\Omega=\omega_{s}=2\pi/T_{orb}$, where
$T_{orb}=15.945$ days is Titan’s orbital period.

Our models show that resonant stratification can explain the $k_{2}$
enhancement observed by Cassini. We observe an enhancement $\Delta k_{2}$
beyond $+15\%$ when overtones composed of internal gravity waves are
resonantly excited by eccentricity tides (Fig. 2). This brings $k_{2}$ from
$3\sigma$ to $2\sigma$ away from the mean value of the observation (Fig. 2).
In this case, we have used the conservative linear damping $\gamma=10^{-9}$
s-1, but our models predict a resonant $\Delta k_{2}$ beyond $+45\%$ when
using a still realistic $\gamma=3.3\times 10^{-10}$ s-1 (Fig. 2), an
enhancement that puts $k_{2}$ at the mean value of the observation at the
saturation point of the resonance. Resonances occur at various $H$ and
$N^{2}$, preventing us from identifying a unique ocean thickness and
stratification profile based solely on $k_{2}$. In resonant stratification,
the salt concentration near the ocean floor can be as low as $<5$ g/kg in the
simplified model we use here (Fig. 2; equation (22)). The mean salinity can
then be less than that for the oceans of Earth or Enceladus (Postberg et al.,
2009), while still producing a resonant response.

Figure 3: Meridional cross section of the radial displacement on Titan’s ocean
due to tides for selected internal gravity wave normal modes of increasing
radial order shown in Fig. 2. Frictional dissipation is $\gamma=10^{-9}~{}\rm
s^{-1}$.

As a general rule, internal gravity waves become resonant when the nodes in
the radial displacement of the ocean wave perfectly fit the thickness of the
stratified ocean cavity (Fig. 3). We can achieve this resonant fit by
adjusting $N^{2}$ and changing the radial wavelength of internal gravity
waves, or by adjusting the thickness of the stratified ocean $H$. The number
of radial nodes is directly proportional to overtone order and inversely
proportional to mode frequency, with lower frequency internal gravity waves
having more radial nodes (Unno et al., 1979). Increasing the strength of
stratification $N^{2}$ shifts the spectrum of g-modes toward high frequency
and allows higher order gravity modes to become resonant with the fixed
orbital frequency (Fig. 4).

Figure 4: Heating rate produced by tidally excited internal gravity waves in
the stratified ocean, as a function of (a) ocean thickness $H$ and (b) Brunt-
Vaisala frequency $N^{2}$. The dashed line indicates Titan’s radiogenic
heating $\dot{E}_{int}\sim 3\times 10^{11}$ W (see Apprendix D). The
frictional damping is $\gamma=10^{-9}$ s-1. Overtones of g-modes are labeled
with a subindex representing the mode radial order $n$. The smaller peaks
without a label represent dissipation from resonant inertial wave
modes/attractors not discussed here (see Rovira-Navarro et al. (2019), e.g.).

###  3.2 Surface radial displacement and interior nonlinear wave breaking

The dynamical enhancement $\Delta k_{2}$ emerges from a tens-of-meter
enhancement on the radial displacement of the ocean surface.

$h_{\ell m}=g\left(\frac{\xi_{r\ell m}}{\phi^{e}_{\ell
m}}\right)_{r=R}\textnormal{,}$ (32)

resulting in the relationship

$h_{\ell
m}=\left(\frac{2\ell+1}{3}\right)\left(\frac{\bar{\rho}}{\rho}\right)k_{\ell
m}\textnormal{.}$ (33)

,

$|\xi_{r}|\approx\frac{7eh_{2}U_{22}}{g}=7eh_{2}\left(\frac{3\pi}{10}\right)^{1/2}\left(\frac{MR^{4}}{m_{s}a^{3}}\right)\approx
26\textnormal{ m,}$ (34)

where $M$ is Saturn’s mass and $m_{s}$ is Titan’s mass.

Below the ocean surface, resonant internal gravity waves produce negligible
gravity but attain a radial displacement $|\xi_{r}|\sim 1$ km in the case of
$\gamma=10^{-9}$ s-1 (Fig. 3). A lower dissipation $\gamma$ produces even
larger resonant amplitudes of tidal motion interior to the ocean (Fig. 5),
with $\gamma=10^{-10}$ s-1 reaching $|\xi_{r}|\sim 10$ km. This large radial
displacement is accompanied by a horizontal displacement that is typically
$\xi_{\perp}/\xi_{r}\sim nR/H\sim 50$ (equation (C7)) and allows the flow to
preserve continuity in an incompressible fluid. Despite the large tidal
displacement, our models of resonant tidal excitation remain far from
experiencing nonlinear wave breaking when $\gamma\gtrsim 10^{-10}$ s-1, as
stipulated in $|\xi k|\lesssim 1$, the typical criterion to avoid nonlinear
wave breaking, where $k$ is the wavenumber. For the radial displacement, this
criterion translates to $|\xi_{r}|\lesssim H/n$. In the case of the horizontal
displacement, the criterion stipulates $\xi_{\perp}\lesssim\pi R/m$, where $m$
is the azimuthal order. All g-modes shown in Figs. 2 and 3 avoid nonlinear
wave breaking by at least one order of magnitude.

###  3.3 Heating rate at saturation point

Our models also show that resonant stratification produces enough heat to
compete with the radiogenic heating generated by decaying isotopes inside
solid Titan (Fig. 4). This result is key to allow Titan to reach a stable
fixed point that may prolong the crossing of a resonance as the ocean freezes.
When in thermal steady state, the heat transported across the ice shell must
balance all interior heat sources. If heat is transported by conduction, we
can write $d\sim(3\times 10^{13})/\dot{E}_{int}$ km (Appendix D), where
$\dot{E}_{int}$ is the interior heating rate in W and $d$ is the ice shell
thickness. When the ocean freezes until balance with radiogenic heating,
$\dot{E}_{int}\sim 3\times 10^{11}$ W (Appendix D) and the ice shell thickness
grows to $d\sim 100$ km. In this scenario, the ocean plays a negligible role
in heating the interior. Resonant stratification changes this picture by
introducing an additional heating source when the ocean is still rapidly
freezing from secular cooling. In resonant stratification with
$\gamma=10^{-9}$ s-1, for example, we get a total $\dot{E}_{int}\sim 6\times
10^{11}$ W at saturation point (Fig. 4) and the ice shell thickness becomes
$d\sim 50$ km while in steady state with internal heat from radiogenic heating
and ocean tidal heating combined.

The model above is a simple conductive model, but it provides a plausible
argument in favor of catching resonant stratification via ocean freezing.
Thermal steady state cannot be attained when the ice shell thickness is very
thin near at the onset of ice shell formation ($d\lesssim 50$ km for
$\gamma=10^{-9}$ s-1); the resonance saturates at peak heating rate (Fig. 4)
before producing enough heat to balance conduction across the thin ice shell.
Ice shell growth by secular cooling stops at $d\sim 100$ km, hence resonant
stratification via ocean freezing must catch a resonance before then. This
last requirement is not significantly changed when a high abundance of ammonia
in the ocean is considered. Ammonia can lower the freezing temperature of the
ocean to the eutectic at $T\sim 180$ K if present in concentration ${\sim}15$
wt.$\%$ (Lunine & Stevenson, 1987; Grasset & Sotin, 1996), reducing the $d$
required for thermal steady state in half independently of whether resonant
stratification is in place or not (see equation (D2)).

A lower $\gamma=10^{-10}$ s-1 increases the heating rate at saturation point
to an order of magnitude above Titan’s radiogenic heating rate, and allows
$k_{2}$ to grow higher (Fig. 5). A larger $\gamma$ leads to the damping of
resonances, reducing the amplitude of the resonant $\Delta k_{2}$ and the
resonant $\dot{E}$ (Fig. 5). The hydrostatic flow is not much affected by
$\gamma$, thus the heating rate is typically increased by a larger $\gamma$
when far from resonances. Our results are relevant in the range
$\gamma=10^{-9}$ – $10^{-10}$ s-1; a $\gamma$ larger than this range damps
$\Delta k_{2}$ below the signal registered by Cassini, whereas a $\gamma$
lower than the range results in nonlinear breaking of internal gravity waves,
an effect not included in our models.

Figure 5:  Resonant dynamical gravity signal and heat production as a function
of dissipation for g-mode $g_{6}$ ($\ell=m=2$ and $n=6$). Ocean thickness is
$H=300$ km. The dashed lines are the same as in Fig. 2, and Fig. 4,
respectively for (a) and (b).

##  4 Discussion

###  4.1 Approach to a long-lived resonance via ocean freezing

In addition to thermal steady state, ocean freezing requires an approach
toward the resonance from the branch that produces a positive $\Delta k_{2}$
(Fig. 2). Approaching the resonance from the negative $\Delta k_{2}$ branch
would depart from the gravity enhancement required by Cassini. A simple
freezing model of a fully stratified ocean fails to satisfy this criterion. In
this freezing model, the ocean freezes leaving salts in the liquid phase and
redistributing them into a linear salinity profile that conserves the total
mass of salts $M_{s}$. Following equation (22), the ocean stratification in
this case follows

$N^{2}\sim-\frac{gM_{s}}{2\pi\rho R^{2}H^{2}}\textnormal{.}$ (35)

Despite its attractive simplicity, a fully stratified ocean freezes following
a trajectory that never converges toward a g-mode resonance, independently of
the $M_{s}$ assumed (Fig. 6).

In an alternative freezing model, a diffusive layer of thickness $\delta$ is
sandwiched between two layers of constant density, where the top layer has
roughly no salinity and the bottom layer is high salinity. Internal gravity
waves are excited in the stratified diffusive layer instead of the entire
ocean. In this alternative case, the ocean freezes without changing the
diffusive layer thickness $\delta$, but increasing the density contrast
$\delta\rho$ and consequently increasing $N^{2}$ (equation (22)), according to

$N^{2}\sim-\frac{gM_{s}}{4\pi\rho R^{2}h_{b}\delta}\textnormal{,}$ (36)

where $h_{b}$ is the thickness of the bottom layer with high salinity. The
freezing of an ocean with a diffusive layer (i.e., reducing $h_{b}$ in
equation (36) at constant $\delta$) suggest the possibility of crossing the
$\Delta k_{2}$ resonance from the positive branch (Fig. 6). Further
theoretical development is required to determine the effects of a diffusive
layer in the $\Delta k_{2}$ results reported here.

Figure 6: Model parameters that produce tidally excited resonances with
internal gravity waves (solid lines; see equation (A2)). The positive branch
of the $\Delta k_{2}$ resonance extends toward the bottom-left of each solid
line. The red crosses represent individual resonances identified in our
numerical simulations and shown in Fig. 3. The dashed lines represent
simplified freezing trajectories for different models of ocean stratification.
The fully stratified ocean model assumes $M_{s}\approx 8.3\times 10^{21}$ g,
which is equivalent to a $\delta\rho\sim 0.001$ g cm-3 over an ocean thickness
$H\sim 200$ km.

A stable fixed point is then established once resonant stratification succeeds
at reaching thermal steady state by halting ocean freezing with the additional
heating rate provided by the resonance (Tyler, 2011, 2014). The secular
expansion of Titan’s orbit pushes the orbital frequency away from resonance,
reducing the heating rate and promoting further ice thickening until ocean
waves tune again with the new orbital frequency. This can happen because ice
shell thickening by secular cooling can keep up with the orbital evolution.
For instance, the thermal adjustment timescale for a 100 km thick shell is
$d^{2}/\pi^{2}\kappa\approx 30$ Myr; over this timescale the orbital frequency
will have changed by 0.5% (Section 4.2). The alternative to a stable fixed
point is that Titan has by chance encountered a resonance that will last until
further orbital evolution pushes the orbital frequency out of resonance.

###  4.2 Resonant stratification via orbital evolution

In the absence of a stable fixed point, resonant stratification can still be
established by orbital evolution. Titan orbits Saturn at a slower rate than
Saturn’s rotation, leading to outward migration from tidal torques that arise
after tidal dissipation occurs inside Saturn (Lainey et al., 2020). This
outward migration imposes a slow drift in the excitation frequency of
eccentricity tides. At some point during orbital migration, the frequency of
eccentricity tides can match the frequency of a g-mode overtone trapped in the
stratified ocean, setting resonant stratification (Fig. 7). This mechanism
operates without requiring any specific freezing history for Titan’s ocean.
The resulting resonance, however, is short lived. Continuous orbital migration
will further break the resonance after pushing Titan through the resonance
width $\delta a\sim 0.02R_{s}$, where $R_{s}$ is Saturn’s radius (Fig. 7) and
we have assumed the reasonable $\gamma=10^{-9}$ s-1 (i.e., the resonance width
depends on dissipation, as shown in Fig. 5). At Titan’s current migration rate
$\dot{a}/a\sim 10^{-10}$ yr-1 (Lainey et al., 2020), the time it would take
Titan to cover the $g_{6}$-mode resonance width would be $\delta a/\dot{a}\sim
10$ Myr.

Figure 7:  Same as Fig. 2, but as a function of orbital semi-major axis $a$.
The shaded vertical line represents Titan’s current semi-major axis
$a=21R_{s}$, where $R_{s}$ is Saturn’s radius. The ocean is fixed to $H=300$
km and $N^{2}=1.4454\times 10^{-8}$ s-2 (see also Fig. 3). The tidal frequency
is equal to the orbital frequency and $\omega_{0}$ is Titan’s current orbital
frequency $\omega_{s}$. The eccentricity is fixed to Titan’s current $e$ and
we assume synchronous rotation at all times (i.e., $\Omega=\omega_{s}$).

The probability that Titan is simply passing over an ocean resonance is slim
as a result of this fast orbital migration. The astrometric observation of
Titan’s migration (Lainey et al., 2020) indicates that Titan’s orbital period
has increased ${\sim}3$ times by orbital expansion over the lifetime of the
solar system $\lambda\sim 4.5\times 10^{9}$ yr, or

$\frac{\Delta P_{s}}{P_{s}}=\frac{3\Delta
a}{2a}=\frac{3\lambda\dot{a}}{2a}\sim\frac{2}{3}\textnormal{,}$ (37)

where $P_{s}$ is Titan’s current orbital period, and $\Delta P_{s}$ and
$\Delta a$ represent orbital parameter changes on the timescale $\lambda$. The
period spacing of g-modes (Appendix A) allows us to estimate the number of
g-modes that Titan crosses over the timescale $\lambda$,

$\\#_{g}\sim 1+\frac{\Delta P_{s}}{\Delta P_{g}}\sim
1+\frac{3\lambda\dot{a}P_{s}}{2a\Delta P_{g}}\textnormal{,}$ (38)

where $\Delta P_{g}$ is the g-mode spacing. Cassini registered Titan’s
enhanced $k_{2}$ at no particular time within Titan’s interior evolution. The
probability that Cassini observed a g-mode resonance motivated uniquely by
orbital evolution (i.e., no stable fixed-point) is

$\mathcal{P}(\textnormal{g-mode})\sim\\#_{g}\left(\frac{\delta a}{\Delta
a}\right)\sim\frac{\delta
a}{a}\left(\frac{a}{\lambda\dot{a}}+\frac{3\pi}{\omega_{s}\Delta
P_{g}}\right)\textnormal{.}$ (39)

We obtain $\\#_{g}\sim 5$ and $\mathcal{P}(\textnormal{g-mode})\sim 1\%$ when
using the reasonable ocean parameters in Fig. 7. Equation (39) is only valid
for $\\#_{g}\gtrsim 2$. The low $\mathcal{P}(\textnormal{g-mode})$ indicates
that resonant stratification is unlikely to be a result of pure orbital
evolution (i.e., no ocean freezing involved), yet not impossible.

###  4.3 A mildly salty ocean versus a heavy ocean

When compared to previous studies, resonant stratification only requires a
mild concentration of solute dissolved in the ocean to explain Titan’s $k_{2}$
(Fig. 8). With $\gamma=10^{-9}$ s-1, resonant stratification yields a $k_{2}$
value 1$\sigma$ closer to the measured central value starting from the
hydrostatic $k_{2}=0.42$. When $\gamma=3.3\times 10^{-10}$ s-1, the
enhancement over the hydrostatic $k_{2}$ is 3$\sigma$, crossing the mean value
of the Cassini $k_{2}$ observation at low salinity ($S<10$ g/kg). In the
absence of resonant stratification, a heavy convective ocean requires a salt
concentration that is on average $S\sim 100-200$ g/kg higher to produce a
similar effect on $k_{2}$, depending on the exact $\gamma$ (Fig. 8). Previous
studies have argued in favor of a heavy convective ocean that holds $S\sim
200$ g/kg in salts to obtain a $1\sigma$ agreement with $k_{2}$ observations.
However, water-rock interactions at the bottom of Titan’s ocean are expected
to produce a limited $S\lesssim 10$ g/kg of salts (Leitner & Lunine, 2019).

Figure 8: Titan’s $k_{2}$ fractional enhancement as a function of ocean bulk
salinity. For an ice shell thickness $d\sim 100$ km and a given bulk salinity,
resonant stratification (RS; solid red) can occur at various combinations of
ocean thickness and B-.V. frequency, producing $\Delta k_{2}\approx+15\%$
($\gamma=10^{-9}$ s-1) and $\Delta k_{2}\approx+45\%$ ($\gamma=3.3\times
10^{-10}$ s-1) (see also Fig. 3). The hydrostatic $k_{2}$ of a fully
homogeneous ocean is shown in a dashed blue line (ice shell thickness $d\sim
100$ km) and a solid black line (no ice shell). The background shade of yellow
represents the fractional enhancement required to explain Cassini’s
observation at various confidence intervals (Durante et al., 2019). The
salinity in Enceladus’ plume is from Postberg et al. (2009) and Titan previous
studies are Rappaport et al. (2008); Iess et al. (2012); Mitri et al. (2014);
Durante et al. (2019).

Instead of water-rock interactions, the heavy convective ocean relies on a
special event to attain its relatively high salt concentration. Ammonium
sulfate ((NH4)2SO4) could form in the ocean from reactions between water-
ammonia and brine leaching upwards from a core experiencing hydration (Fortes
et al., 2007), contributing to $S\sim 200$ g/kg of dissolved solute.
Unfortunately, this scenario leads to predicted surficial expressions on Titan
that failed to be observed during the Cassini mission (Leitner & Lunine,
2019). Alternatively, magnesium sulfates can be incorporated into a heavy
ocean that is thermodynamically consistent (Vance & Brown, 2013) via a late-
delivery of salt-rich carbonaceous chondrites (Hogenboom et al., 1995).
However, it is not clear whether this delivery mechanism can provide the
required salt concentration.

Water-rock interactions on Titan’s ocean floor (Leitner & Lunine, 2019)
produce enough salts ($S\sim 10$ g/kg) to enable resonant stratification. At
$S\sim 10$ g/kg salinity, the predicted $k_{2}$ is in a $2\sigma$ agreement
with the Cassini observation when $\gamma=10^{-9}$ s-1 and within $1\sigma$
agreement when $\gamma=3.3\times 10^{-10}$ s-1 (Fig. 8). Both $\gamma$ values
are realistic and prevent nonlinear wave breaking. This bulk concentration of
salts ($S\sim 10$ g/kg) is compatible with the average salinity of Earth’s
oceans and the salinity inferred for Enceladus’ ocean from direct sampling of
E-ring particles provided by Enceladus’ plumes (Postberg et al., 2009). This
aspect favors resonant stratification over a heavy ocean because water-rock
interactions are better understood than the special event required to explain
a high concentration of salts, namely early hydration during the interior
differentiation of ices and silicates (Fortes et al., 2007) or a late delivery
of carbonaceous chondrites with a special composition (Hogenboom et al.,
1995).

###  4.4 Ocean stability to overturning convection

The heat flux at the ocean floor provided by interior heating threatens the
stability of a weakly stratified ocean. A thermal gradient can counter the
stratification produced by a chemical gradient and lead to unstable
overturning convection. Convection introduces mixing that can further erase
the chemical gradient over time when no salinity forcing is introduced. From a
simple balance of the density profile including thermal effects and a salinity
gradient (i.e., the Ledoux instability criterion), we require

$\Delta T\lesssim\frac{HN^{2}}{\alpha g}\sim 10\textnormal{ K}$ (40)

across an ocean thickness to preserve the stability of the mild stratification
$N^{2}\sim 1\times 10^{-8}$ s-2 discussed before, . This temperature gradient
is equivalent to a heat transfer ${\lesssim}2\times 10^{9}$ W across the ocean
by thermal conduction, two orders of magnitude lower than the heating rate
expected from radiogenic heating (Appendix D).

One might then presume that the heat flux from radiogenic heating should
destroy any mild stable stratification. However, this need not necessarily be
the case. For example, plumes from hydrothermal vents or volcanoes on the
surface of the silicate core may pierce through the stratified ocean in a
Rayleigh-Taylor instability (Collins & Goodman, 2007) and allow the radiogenic
heat from the solid interior to escape outward without triggering overturning
convection at the scale of the entire ocean. In this hypothetical scenario,
heat passes across the ocean in small lengthscale plumes that do not disturb
the large lengthscale structure of stratification required by resonant
stratification.

Double-diffusive convection (Radko, 2013) constitutes another hypothetical
scenario that may permit maintenance of the radiogenic heat flux without
erasing the compositional gradient. This regime is typically observed when the
temperature profile is steeper than the convective temperature profile and
less steep than the Ledoux instability criterion (Leconte & Chabrier, 2012,
e.g.). The convective temperature profile prescribes a $\Delta T\sim\alpha
HgT/c_{p}\sim 5.2$ K (Turcotte & Schubert, 2002) over an ocean with the same
properties used in equation (40), where $c_{p}\sim 4$ J g-1 K-1 is the heat
capacity. Assuming thermal equilibrium between conduction through the ocean
thickness and the radiogenic heat flux (Appendix D), we require

$N^{2}\gtrsim\frac{\alpha g\dot{E}_{int}}{4\pi\rho c_{p}R^{2}_{p}\kappa}\sim
1.6\times 10^{-6}\left(\frac{1.5\times 10^{-3}\textnormal{
cm}^{2}\textnormal{s}^{-1}}{\kappa}\right)\textnormal{ s}^{-2}$ (41)

to maintain Ledoux stability (equation (40)) of the ocean salinity gradient,
where $\kappa$ is the thermal diffusivity of the ocean. In regimes typical of
gas giant planets, numerical simulations suggest that double-diffusive
convection can increase the efficiency of heat transport by up to a factor of
${\sim}50$ compared to heat conduction (Rosenblum et al., 2011; Mirouh et al.,
2012). This enhancement can be thought of as being roughly equivalent to an
increase in $\kappa$ that allows a salinity gradient with $N^{2}\sim 10^{-8}$
s-2 to remain Ledoux stable (equation (41)). Numerical simulations confirm a
strong enhancement of heat transport by double-diffusive convection in the
regime of icy satellites (Wong et al., 2022), where, contrary to the gas giant
planets, the kinematic viscosity $\nu$ is typically larger than the thermal
diffusivity $\kappa$. Double-diffusive convection is typically accompanied by
the evolution of the compositional gradient (Radko, 2013; Wong et al., 2022),
thus further studies are required to better understand the timescales imposed
on the evolution of salinity profiles in icy satellites.

###  4.5 Predictions and future tests

Water-rock interactions can occur at the ocean bottom of other icy satellites,
thus an enhancement of $k_{2}$ by resonant stratification is also possible on
Ganymede and Europa. NASA’s Europa Clipper mission (Phillips & Pappalardo,
2014; Howell & Pappalardo, 2020) will measure Europa’s $k_{2}$ with an
expected accuracy of $2\%$ (Mazarico et al., 2023), providing us with a new
opportunity to study the interior of an icy satellite. If resonant
stratification is a common mode of operation for icy satellites, we should
also observe an important enhancement in Europa’s $k_{2}$ given that the ice
shell elasticity only provides small resistance to vertical tidal
displacements. ESA’s JUICE mission (Grasset et al., 2013) will measure
Ganymede’s $k_{2}$ with even greater precision due to the orbital design of
the mission and the improved K-band antenna onboard the spacecraft. Ganymede’s
tidal response will be measured at the excitation frequency of the various
moon-to-moon tides present in the Galilean moon system (De Marchi et al.,
2022), in addition to the conventional eccentricity tide raised by Jupiter.
This tidal response spectrum will provide a unique opportunity to measure the
potential stratification of an ocean regardless of whether resonant
stratification is in place or not, in addition to providing further
constraints to ocean thickness from sampling high-frequency moon-to-moon
tides. One quantity to look for in the tidal response spectrum is the $g$-mode
spacing, which is both sensitive to the degree of stratification and the
thickness of the stratified cavity (Appendix A).

##  5 Conclusions

We calculated Titan’s tidal response to eccentricity tides using a new
theoretical framework that includes the dynamical effects of tidally excited
waves trapped in the ocean. Our results present a new interpretation of
Cassini’s observation of Titan’s Love number $k_{2}=0.616\pm 0.067$, which is
$3-\sigma$ away from the predicted $k_{2}$ in an ocean of pure water resting
on top a rigid ocean floor. If Titan’s ocean is stably stratified, its
measured tidal response can be fully reproduced using plausible dissipation
factors ($\gamma$) without requiring a salinity greater than those of Earth or
Enceladus’s oceans (Fig. 8). This enhanced response requires the ocean to be
set in resonance with the period of the current tidal excitation, namely
Titan’s orbital period. In one possible scenario, this resonance is
encountered as the ocean progressively freezes and develops a deep salty layer
(Fig. 6); this situation yields a long-term stable thermal equilibrium with
conduction across the ice shell.

Studies on the extent to which stably stratified oceans can be maintained
against convective mixing would form a valuable theoretical addition to the
current work. Processes similar to those hypothesized to be operating at Titan
could be in play at Ganymede or Europa, and may be tested with future
spacecraft missions Europa Clipper and JUICE. The seismometer expected on the
Dragonfly mission (Barnes et al., 2021), or a future Titan orbiter (Sotin et
al., 2017, e.g.), might similarly be able to look for evidence of a
resonantly-excited ocean.

We thank David Stevenson for insightful early discussions leading to this
manuscript. We are thankful for the constructive comments of Mikael Beuthe and
one anonymous referee.

##  Appendix A The spectrum of internal gravity waves in a stratified ocean

In the absence of rotation, the spectrum of normal modes associated with
internal gravity waves (e.g., g-modes) follows (Unno et al., 1979)

${}_{n}\omega_{\ell}\simeq\frac{\sqrt{\ell(\ell+1)}}{n\pi}\int_{0}^{R}\frac{N}{r}dr\textnormal{,}$
(A1)

where $n$ is the overtone radial order, $\ell$ is degree, and $N$ the Brunt-
Vaisala frequency. For Titan, we can estimate the $N$ that produces a
resonance with eccentricity tides of tidal frequency $\omega=\omega_{s}$ using
our simplified linear model of stratification (i.e., constant $N$ across a
stratified cavity of thickness $H$),

$N=-0.424\left(\frac{n\pi}{\ln(1-H/R)\sqrt{\ell(\ell+1)}}\right)\omega_{s}\textnormal{.}$
(A2)

We obtain the numerical factor in front of the previous expression from a fit
of equation (A1) to individual resonances identified in our numerical
simulations (Fig. 6). This factor represents the shift caused on the g-mode
frequency due to the effects of Titan’s rotation with spin rate
$\Omega=\omega_{s}$.

The g-mode frequency spacing constitutes a diagnostic quantity typically used
to characterize stratified cavities inside planets and stars (Aerts et al.,
2010; Mankovich & Fuller, 2021, e.g.). Following our simple model of constant
$N$ across a stratified cavity of thickness $H$, we obtain the mode spacing

$\Delta
P_{g}\sim-0.424\left(\frac{2\pi^{2}}{N\ln(1-H/R)\sqrt{\ell(\ell+1)}}\right)\textnormal{.}$
(A3)

This expression can be used to characterize the stratification of oceans in
icy satellites when a multi-frequency $k_{2}$ is available to observation, as
currently expected from JUICE measurements of moon-to-moon tides on Ganymede
(De Marchi et al., 2022). Future icy satellite seismology could provide an
alternative observation of mode spacing from the recording of free
oscillations on the satellite’s surface (Marusiak et al., 2021), assuming that
the normal modes are excited beyond the detection threshold.

##  Appendix B The tidal forcing of eccentricity tides

The gravitational potential $\phi^{T}$ experienced by an observer at
$\mathbf{r}$ from a concentrated mass located at $\mathbf{r}^{\prime}$, is
inversely proportional to the distance between them

$\phi^{T}\propto\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}=\sum_{\ell=0}^{\infty}\frac{r^{\ell}}{r^{\prime\ell+1}}\mathcal{P}_{\ell}(\cos\alpha)\textnormal{,}$
(B1)

where $\alpha$ is the angle between the two position vectors, $\ell$ is
degree, and $\mathcal{P_{\ell}}$ are the $m=0$ case of the associated Legendre
polynomials

$\mathcal{P}_{\ell}^{m}(\mu)=\frac{(-1)^{m}}{2^{\ell}\ell!}(1-\mu^{2})^{m/2}\frac{d^{\ell+m}}{d\mu^{\ell+m}}(\mu^{2}-1)^{\ell}\textnormal{,}$
(B2)

where $m$ is azimuthal order. When the concentrated mass $M$ is that of a
planet orbiting at a semi-major axis $a$, the gravitational excitation assumes
the form

$\phi^{T}=\frac{\mathcal{G}M}{a}\sum_{\ell=2}^{\infty}\left(\frac{r}{a}\right)^{\ell}\mathcal{P}_{\ell}(\cos\alpha)\textnormal{.}$
(B3)

The two lowest degree harmonics, $\ell=0$ and $\ell=1$, are discarded as they
do not disturb the shape of the icy satellite.

The addition theorem allows us to express the gravitational excitation in
spherical coordinates, following

$\cos\alpha=\cos\theta\cos\theta_{p}+\sin\varphi\sin\varphi_{p}\cos(\varphi-\varphi_{p})\textnormal{,}$
(B4)

where $r\theta\varphi$ are spherical polar coordinates in a corotating frame
fixed to the icy satellite and the subscript $p$ denotes the position of the
planet. The gravitational excitation is now

$\phi^{T}=\frac{\mathcal{G}M}{a}\sum_{\ell=2}^{\infty}\frac{4\pi}{2\ell+1}\left(\frac{r}{a}\right)^{\ell}\sum_{m=-\ell}^{\ell}Y_{\ell}^{m}(\theta,\varphi)Y_{\ell}^{m*}(\theta_{p},\varphi_{p})\textnormal{,}$
(B5)

where we use the conventional definition of spherical harmonics

$Y_{\ell}^{m}(\theta,\varphi)=\left(\frac{(2\ell+1)}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}\right)^{1/2}\mathcal{P}_{\ell}^{m}(\cos\theta)e^{im\varphi}\textnormal{.}$
(B6)

From fundamental identities of spherical harmonics, it can be shown that

$Y_{\ell}^{m*}(\theta_{p},\varphi_{p})=\left(\frac{(2\ell+1)}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}\right)^{1/2}P_{\ell}^{m}(\cos\theta_{p})e^{-im\varphi_{p}}\textnormal{.}$
(B7)

The tidal forcing potential for $\ell m$ then becomes

$\phi^{T}_{\ell
m}=\frac{\mathcal{G}M}{a}\left(\frac{r}{a}\right)^{\ell}\left(\frac{4\pi}{2\ell+1}\frac{(\ell-m)!}{(\ell+m)!}\right)^{1/2}\mathcal{P}_{\ell}^{m}(\cos\theta_{p})e^{-im\varphi_{p}}Y_{\ell}^{m}(\theta,\varphi)\textnormal{.}$
(B8)

In the standard case of a coplanar circular orbit, we have
$\cos\theta_{p}=\varphi_{p}=0$, leading to

$\phi^{T}_{\ell m}=U_{\ell
m}\left(\frac{r}{R}\right)^{\ell}Y_{\ell}^{m}(\theta,\varphi)e^{-i\omega
t}\textnormal{,}$ (B9)

with the normalization constant

$U_{\ell
m}=\left(\frac{\mathcal{G}M}{a}\right)\left(\frac{R}{a}\right)^{\ell}\left(\frac{4\pi(\ell-m)!}{(2\ell+1)(\ell+m)!}\right)^{1/2}\mathcal{P}_{l}^{m}(0)\textnormal{.}$
(B10)

The effect of eccentricity imposes a change in the semimajor axis and a
libration in the position of the planet as seen from the corotating frame on
the icy satellite,

$a=a_{0}(1-e\cos\omega t)\textnormal{,}$ (B11) $\varphi_{p}=2e\sin\omega
t\textnormal{.}$ (B12)

Our strategy is now to expand equation (B8) to first order in $e$. The
semimajor axis dependency expands as

$\frac{1}{a^{\ell+1}}\approx\frac{1}{a_{0}^{\ell+1}}\left(1+(\ell+1)e\cos\omega
t\right)\textnormal{,}$ (B13)

and the libration of the planet expands as

$e^{-im\varphi_{p}}\approx 1-2ime\sin\omega t\textnormal{.}$ (B14)

Using the standard trigonometric definitions

$\sin z=\frac{e^{iz}-e^{-iz}}{2i}\textnormal{,}$ (B15) $\cos
z=\frac{e^{iz}+e^{-iz}}{2}\textnormal{,}$ (B16)

we obtain

$\frac{1}{a^{\ell+1}}\approx\frac{1}{a_{0}^{\ell+1}}\left(1+\frac{e(\ell+1)}{2}(e^{i\omega
t}+e^{-i\omega t})\right)\textnormal{,}$ (B17) $e^{-im\varphi_{p}}\approx
1-em(e^{i\omega t}-e^{-i\omega t})\textnormal{.}$ (B18)

The resulting tidal excitation potential in a planar orbit is

$\phi^{T}_{\ell m}\approx U_{\ell
m}\left(\frac{r}{R}\right)^{\ell}\left(1+e\left(\frac{(\ell+1)}{2}-m\right)e^{i\omega
t}+e\left(\frac{(\ell+1)}{2}+m\right)e^{-i\omega
t}\right)Y_{\ell}^{m}(\theta,\varphi)\textnormal{.}$ (B19)

The circular tide becomes static (i.e., no time dependence) in a synchronous
corotating frame where the spin of the icy world matches the orbital frequency
of the planet. Eccentricity tides propagate at the diurnal frequency in both
west and east directions for a given $m$. We observe a perfect superposition
between the east $m>0$ tide and the west $m<0$ tides. As a result, the
contributions are typically added to obtain

$\phi^{T}_{\ell m}\approx U_{\ell
m}\left(\frac{r}{R}\right)^{\ell}Y_{\ell}^{m}(\theta,\varphi)\left(1+e\left(\frac{(\ell+1)}{2}-m\right)e^{i\omega
t}+e\left(\frac{(\ell+1)}{2}+m\right)e^{-i\omega t}\right)\textnormal{.}$
(B20)

We make use of the fact that the direction of tides can be flipped by either
changing the sign in $\varphi$ or the tidal frequency $\omega_{s}$, reducing
the eccentricity tidal potential to

$\phi^{e}_{\ell m}\approx e\left(\ell+1+2m\right)U_{\ell
m}\left(\frac{r}{R}\right)^{\ell}Y_{\ell}^{m}(\theta,\varphi)e^{-i\omega
t}\textnormal{.}$ (B21)

In this paper, all quantities mimic the time dependence of the tidal forcing.
East tides correspond to $m>0$ and west tides to $m<0$. East tides propagate
in the direction of rotation, whereas west tides are counter-rotation. Notice
that the amplitude of eccentricity tides is $7e$ fold compared to the
amplitude of the static tides in the case $\ell=m=2$.

##  Appendix C Projection of dynamical tides onto vectorial spherical
harmonics

We project our equations onto vectorial spherical harmonics following the
standard decomposition

$\bm{\xi}=y_{1}{\bf Y}+y_{2}{\bf\Psi}+y_{3}{\bf\Phi}\textnormal{,}$ (C1)

where ${\bf Y}{\bf\Psi}{\bf\Phi}$ constitute an orthonormal base for the
projection of vectorial fields in spherical polar coordinates. VSH relate to
scalar spherical harmonics (equation (B6)) as

${\bf Y}=Y\hat{\bf r}\textnormal{,}$ (C2) ${\bf\Psi}=r\nabla Y\textnormal{,}$
(C3) ${\bf\Phi}=r\nabla\times{\bf Y}\textnormal{,}$ (C4)

where spherical harmonics satisfy $r^{2}\nabla^{2}Y=-\ell(\ell+1)Y$. We
further project the gravity response and pressure-gravity potentials,
respectively, as

$\tilde{\phi}^{\prime}=\tilde{\phi}^{\prime}(r)Y\textnormal{,}$ (C5)
$\psi=\psi(r)Y\textnormal{.}$ (C6)

The projection of the continuity equation (equation (25)) leads to

$y_{2}^{\ell}=\frac{\partial_{r}\left(r^{2}y_{1}^{\ell}\right)}{\ell(\ell+1)r}\textnormal{,}$
(C7)

where the upper index on $y(r)$ denotes degree instead of exponent. This
equation sets a requirement for the radial and spheroidal displacement fields.
The toroidal displacement field is automatically continuous given that it
constitutes a rotor on displacement rather than a relocation of fluid.

The projection in equation (C1) is ideal for solving the vorticity equation
(equation (26)) given the following set of convenient properties (see also
Rieutord (1987)),

$\nabla\times\left(y{\bf Y}\right)=\frac{y}{r}{\bf\Phi}\textnormal{,}$ (C8)
$\nabla\times\left(y{\bf\Psi}\right)=-\frac{\partial_{r}(ry)}{r}{\bf\Phi}\textnormal{,}$
(C9) $\nabla\times\left(y{\bf\Phi}\right)=\frac{\ell(\ell+1)}{r}y{\bf
Y}+\frac{\partial_{r}(ry)}{r}{\bf\Psi}\textnormal{.}$ (C10)

The curl of the tidal displacement in the vorticity equation can be directly
obtained from these properties,

$\nabla\times\bm{\xi}=\frac{\ell(\ell+1)}{r}y_{3}{\bf
Y}+\frac{\partial_{r}(ry_{3})}{r}{\bf\Psi}+\left(\frac{y_{1}-\partial_{r}(ry_{2})}{r}\right){\bf\Phi}\textnormal{.}$
(C11)

The stratification term in the vorticity equation can also be directly
obtained,

$\nabla\times\left(\left(\bm{N}^{2}\cdot\bm{\xi}\right)\hat{\bf
r}\right)=\frac{N^{2}y_{1}}{r}{\bf\Phi}\textnormal{.}$ (C12)

The Coriolis term in the vorticity equation requires extra effort due to the
coupling it imposes on spherical harmonics of different degree (i.e., terms
$\sin\theta\partial_{\theta}Y_{\ell}^{m}$, $\cos\theta Y_{\ell}^{m}$, and
combinations of both). The partial derivatives and trigonometric functions
applied to the spherical harmonics indicate changes in spherical harmonic
quantum numbers described in the following recursive relations (Lockitch &
Friedman, 1999),

$\sin\theta\partial_{\theta}Y_{\ell}^{m}=\ell
q_{\ell+1}Y_{\ell+1}^{m}-(\ell+1)q_{\ell}Y_{\ell-1}^{m}\textnormal{,}$ (C13)
$\cos\theta
Y_{\ell}^{m}=q_{\ell+1}Y_{\ell+1}^{m}+q_{\ell}Y_{\ell-1}^{m}\textnormal{,}$
(C14)

where

$q_{\ell}=\left(\frac{\ell^{2}-m^{2}}{4\ell^{2}-1}\right)^{1/2}\textnormal{.}$
(C15)

We can derive properties from the recursive coupling relations in equations
(C13) and (C14) that permit the projection of the Coriolis term onto radial
functions,

$\bm{\Omega}\times\left(y_{1}{\bf
Y}\right)=-\left(\frac{im\Omega}{\ell(\ell+1)}\right)y_{1}{\bf\Psi}+\frac{\Omega}{\ell(\ell+1)}\left((\ell+1)q_{\ell}y_{1}^{\ell-1}-\ell
q_{\ell+1}y_{1}^{\ell+1}\right){\bf\Phi}\textnormal{,}$ (C16)
$\bm{\Omega}\times\left(y_{2}{\bf\Psi}\right)=-im\Omega
y_{2}{\bm{Y}}-\frac{im\Omega}{\ell(\ell+1)}y_{2}{\bf\Psi}-\frac{\Omega}{\ell(\ell+1)}\left((\ell^{2}-1)q_{\ell}y_{2}^{\ell-1}+\ell(\ell+2)q_{\ell+1}y_{2}^{\ell+1}\right){\bf\Phi}\textnormal{,}$
(C17)
$\displaystyle\bm{\Omega}\times\left(y_{3}{\bf\Phi}\right)=\Omega\left((\ell-1)q_{\ell}y_{3}^{\ell-1}-(\ell+2)q_{\ell+1}y_{3}^{\ell+1}\right){\bf
Y}-\frac{im\Omega}{\ell(\ell+1)}y_{3}{\bf\Phi}$
$\displaystyle+\frac{\Omega}{\ell(\ell+1)}\left((\ell^{2}-1)q_{\ell}y_{3}^{\ell-1}+\ell(\ell+2)q_{\ell+1}y_{3}^{\ell+1}\right){\bf\Psi}\textnormal{.}$
(C18)

The vorticity of the Coriolis term results in a rather long expression. For
brevity, here we only show the ${\bf Y}{\bf\Phi}$ components that result from
combining the vorticity and continuity equations, respectively,

$\frac{q_{\ell}}{\ell^{2}}\partial_{r}\left(\frac{y_{1}^{\ell-1}}{r^{\ell-2}}\right)r^{\ell-1}+\frac{q_{\ell+1}}{(\ell+1)^{2}}\frac{\partial_{r}\left(r^{\ell+3}y_{1}^{\ell+1}\right)}{r^{\ell+2}}+\left(\frac{im}{\ell(\ell+1)}-\frac{i\tilde{\omega}}{2\Omega}\right)y_{3}=0\textnormal{,}$
(C19)
$\frac{iN^{2}}{2\Omega\omega}y_{1}+\left(\frac{i\tilde{\omega}}{2\Omega}-\frac{im}{\ell(\ell+1)}\right)\left(\frac{\partial_{r,r}(r^{2}y_{1})}{\ell(\ell+1)}-y_{1}\right)+\frac{(\ell-1)q_{\ell}}{\ell}r^{\ell}\partial_{r}\left(\frac{y_{3}^{\ell-1}}{r^{\ell-1}}\right)+\frac{(\ell+2)q_{\ell+1}}{(\ell+1)}\frac{\partial_{r}(r^{\ell+2}y_{3}^{\ell+1})}{r^{\ell+1}}=0\textnormal{.}$
(C20)

The first equation tells you that in the absence of rotation, the toroidal
displacement is zero. We can complete the set of equations required to solve
the problem by combining the momentum (equation (23)) and continuity equations
to produce the following ${\bf Y}$ component

$\frac{\partial_{r}\psi}{2\omega\Omega}+\left(\frac{\omega\tilde{\omega}-N^{2}}{2\omega\Omega}\right)y_{1}+\frac{m}{\ell(\ell+1)}\frac{\partial_{r}(r^{2}y_{1})}{r}+i(\ell-1)q_{\ell}y_{3}^{\ell-1}-i(\ell+2)q_{\ell+1}y_{3}^{\ell+1}=0\textnormal{.}$
(C21)

The projection of the ocean heating rate (equation (31)) onto spherical
harmonics follows a much simpler path. Using the base in equation (C1), we
obtain

$\dot{E}=\rho\gamma\omega^{2}\int_{R_{H}}^{R}\left(y_{1}^{2}+\ell(\ell+1)(y_{2}^{2}+y_{3}^{2})\right)r^{2}dr\textnormal{,}$
(C22)

where $y_{1}$, $y_{2}$, and $y_{3}$ are radial functions, the superscript
denotes exponent, and we have used the definitions of orthonormality of
vectorial spherical harmonics

$\int_{\hat{\Omega}}{\bf Y}\cdot{\bf Y}d\hat{\Omega}=1\textnormal{,}$ (C23)
$\int_{\hat{\Omega}}{\bf\Psi}\cdot{\bf\Psi}d\hat{\Omega}=\int_{\hat{\Omega}}{\bf\Phi}\cdot{\bf\Phi}d\hat{\Omega}=\ell(\ell+1)\textnormal{,}$
(C24) $\int_{\hat{\Omega}}{\bf
Y}\cdot{\bf\Psi}d\hat{\Omega}=\int_{\hat{\Omega}}{\bf
Y}\cdot{\bf\Phi}d\hat{\Omega}=\int_{\hat{\Omega}}{\bf\Psi}\cdot{\bf\Phi}d\hat{\Omega}=0\textnormal{,}$
(C25)

where $\hat{\Omega}$ is the solid angle and $d\hat{\Omega}=\sin\theta d\varphi
d\theta$.

##  Appendix D Balance between radiogenic heating and ice shell conduction

In the absence of ocean resonances, Titan’s interior heating rate should be
balanced by heat transfer across the ice shell. If we consider that Titan’s
ice shell operates at a temperature low enough that the most likely heat
transfer mechanism at work is conduction, the heat transferred across a thin,
spherically symmetric ice shell can be approximated by

$\dot{E}\sim\frac{4\pi R^{2}k\Delta T}{d}\textnormal{,}$ (D1)

where $R$ is the mean radius of the ice shell, $k$ is the thermal conductivity
of ice, $\Delta T$ is the temperature difference across the ice shell, and $d$
is the ice shell thickness.

The main source of Titan’s internal heat is radiogenic heating. The bulk
composition of Titan splits into rocks and ices roughly at a mass ratio 1:1
(Stevenson, 1992). We can estimate the radiogenic heating rate from
terrestrial rock samples if we assume a homogeneous distribution of radiogenic
elements along heliocentric distance. The abundances of Th-U-K isotopes and
their radiogenic decays prescribe a heating rate $H\approx 7.4\times 10^{-12}$
W/kg on Earth’s mantle (Turcotte & Schubert, 2002). The resulting radiogenic
heat production is $\dot{E}_{int}\sim 0.5Hm_{s}\sim 5\times 10^{11}$ W.

Titan’s core composition most likely departs from the composition of Earth’s
mantle due to the presence of undifferentiated iron, in which case the
previous estimate is an upper bound. An alternative to the previous estimate
comes from using the radiogenic heat production of CV chondrites $H\approx
4.5\times 10^{-12}$ W/kg (Grasset & Sotin, 1996; Spohn & Schubert, 2003),
which results in the smaller $\dot{E}_{int}\sim 3\times 10^{11}$ W.

The balance between internal radiogenic heating $\dot{E}_{int}$ and conduction
across the ice shell determines the ice shell thermal equilibrium thickness.
Titan’s icy surface is at $T_{s}\sim 90$ K. In a pure water ocean, we can use
$k\sim 2$ W m-1 K-1, Titan radius $R\approx 2575$ km, and $\Delta T\sim 180$ K
to obtain (Luan, 2019)

$d\sim\frac{4\pi R^{2}k\Delta T}{\dot{E}_{int}}\sim\left(\frac{3\times
10^{13}}{\dot{E}_{int}}\right)\textnormal{ km.}$ (D2)

Titan most likely contains a considerable amount of ammonia dissolved in its
global ocean (Lunine & Stevenson, 1987; Stevenson, 1992; Grasset & Sotin,
1996). When dissolved in water, ammonia behaves like an anti-freeze, reducing
the temperature of the eutectic to $T\sim 180K$ in an ocean with 15$\%$wt.
ammonia (Grasset & Sotin, 1996). In this scenario, the temperature gradient
across the ice shell diminishes to $\Delta T\sim 90$ K and the equilibrium ice
shell thickness can be reduced in half compared to the pure water ocean
(equation (D2)).

## References

  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer Science & Business Media)
  * Ashkenazy et al. (2018) Ashkenazy, Y., Sayag, R., & Tziperman, E. 2018, Nature Astronomy, 2, 43
  * Barnes et al. (2021) Barnes, J. W., Turtle, E. P., Trainer, M. G., et al. 2021, The Planetary Science Journal, 2, 130
  * Beuthe (2015) Beuthe, M. 2015, Icarus, 248, 109
  * Beuthe (2016) —. 2016, Icarus, 280, 278
  * Boyd (2001) Boyd, J. P. 2001, Chebyshev and Fourier spectral methods (Courier Corporation)
  * Chen et al. (2014) Chen, E., Nimmo, F., & Glatzmaier, G. 2014, Icarus, 229, 11
  * Collins & Goodman (2007) Collins, G. C., & Goodman, J. C. 2007, Icarus, 189, 72
  * de Kleer et al. (2019) de Kleer, K., McEwen, A. S., Park, R. S., et al. 2019
  * De Marchi et al. (2022) De Marchi, F., Cappuccio, P., Mitri, G., & Iess, L. 2022, Icarus, 386, 115150
  * Dermott (1979) Dermott, S. F. 1979, Icarus, 37, 575
  * Durante et al. (2019) Durante, D., Hemingway, D., Racioppa, P., Iess, L., & Stevenson, D. 2019, Icarus, 326, 123
  * Durham et al. (1998) Durham, W., Kirby, S., & Stern, L. 1998, in Solar System Ices: Based on Reviews Presented at the International Symposium “Solar System Ices” held in Toulouse, France, on March 27–30, 1995, Springer, 63–78
  * Fortes et al. (2007) Fortes, A. D., Grindrod, P. M., Trickett, S., & Vočadlo, L. 2007, Icarus, 188, 139
  * Fuller et al. (2016) Fuller, J., Luan, J., & Quataert, E. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3867
  * Goldreich & Mitchell (2010) Goldreich, P. M., & Mitchell, J. L. 2010, Icarus, 209, 631
  * Goodman & Lackner (2009) Goodman, J., & Lackner, C. 2009, The Astrophysical Journal, 696, 2054
  * Grasset & Sotin (1996) Grasset, O., & Sotin, C. 1996, Icarus, 123, 101
  * Grasset et al. (2013) Grasset, O., Dougherty, M., Coustenis, A., et al. 2013, Planetary and Space Science, 78, 1
  * Grindrod et al. (2008) Grindrod, P. M., Fortes, A. D., Nimmo, F., et al. 2008, Icarus, 197, 137
  * Hay et al. (2020) Hay, H. C., Trinh, A., & Matsuyama, I. 2020, Geophysical Research Letters, 47, e2020GL088317
  * Hogenboom et al. (1995) Hogenboom, D., Kargel, J., Ganasan, J., & Lee, L. 1995, Icarus, 115, 258
  * Howell & Pappalardo (2020) Howell, S. M., & Pappalardo, R. T. 2020, Nature Communications, 11, 1311
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Idini & Stevenson (2021) Idini, B., & Stevenson, D. J. 2021, The Planetary Science Journal, 2, 69
  * Idini & Stevenson (2022) —. 2022, The Planetary Science Journal, 3, 89
  * Iess et al. (2012) Iess, L., Jacobson, R. A., Ducci, M., et al. 2012, Science, 337, 457
  * Jansen et al. (2023) Jansen, M. F., Kang, W., Kite, E. S., & Zeng, Y. 2023, The Planetary Science Journal, 4, 117
  * Kamata et al. (2015) Kamata, S., Matsuyama, I., & Nimmo, F. 2015, Journal of Geophysical Research: Planets, 120, 1528
  * Kang (2023) Kang, W. 2023, Monthly Notices of the Royal Astronomical Society, 525, 5251
  * Kang et al. (2022) Kang, W., Mittal, T., Bire, S., Campin, J.-M., & Marshall, J. 2022, Science advances, 8, eabm4665
  * Lainey et al. (2020) Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, Nature Astronomy, 1
  * Leconte & Chabrier (2012) Leconte, J., & Chabrier, G. 2012, Astronomy & Astrophysics, 540, A20
  * Leitner & Lunine (2019) Leitner, M. A., & Lunine, J. I. 2019, Icarus, 333, 61
  * Lobo et al. (2021) Lobo, A. H., Thompson, A. F., Vance, S. D., & Tharimena, S. 2021, Nature Geoscience, 14, 185
  * Lockitch & Friedman (1999) Lockitch, K. H., & Friedman, J. L. 1999, The Astrophysical Journal, 521, 764
  * Luan (2019) Luan, J. 2019, arXiv preprint arXiv:1905.03802
  * Luan et al. (2018) Luan, J., Fuller, J., & Quataert, E. 2018, Monthly Notices of the Royal Astronomical Society, 473, 5002
  * Lunine & Stevenson (1987) Lunine, J. I., & Stevenson, D. J. 1987, Icarus, 70, 61
  * Mankovich & Fuller (2021) Mankovich, C. R., & Fuller, J. 2021, Nature Astronomy, 5, 1103
  * Marusiak et al. (2021) Marusiak, A. G., Vance, S., Panning, M. P., et al. 2021, The Planetary Science Journal, 2, 150
  * Matsuyama et al. (2018) Matsuyama, I., Beuthe, M., Hay, H. C., Nimmo, F., & Kamata, S. 2018, Icarus, 312, 208
  * Mazarico et al. (2023) Mazarico, E., Buccino, D., Castillo-Rogez, J., et al. 2023, Space Science Reviews, 219, 30
  * Mirouh et al. (2012) Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S. 2012, The Astrophysical Journal, 750, 61
  * Mitri et al. (2014) Mitri, G., Meriggiola, R., Hayes, A., et al. 2014, Icarus, 236, 169
  * Moore & Schubert (2000) Moore, W. B., & Schubert, G. 2000, Icarus, 147, 317
  * Munk et al. (1977) Munk, W. H., MacDonald, G. J., & Anderson, D. L. 1977, The rotation of the Earth: a geophysical discussion, American Association of Physics Teachers
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Nimmo & Bills (2010) Nimmo, F., & Bills, B. 2010, Icarus, 208, 896
  * Nimmo & Pappalardo (2016) Nimmo, F., & Pappalardo, R. T. 2016, Journal of Geophysical Research: Planets, 121, 1378
  * Ogilvie (2009) Ogilvie, G. I. 2009, Monthly Notices of the Royal Astronomical Society, 396, 794
  * Ogilvie & Lin (2004) Ogilvie, G. I., & Lin, D. 2004, The Astrophysical Journal, 610, 477
  * Phillips & Pappalardo (2014) Phillips, C. B., & Pappalardo, R. T. 2014, Eos, Transactions American Geophysical Union, 95, 165
  * Postberg et al. (2009) Postberg, F., Kempf, S., Schmidt, J., et al. 2009, Nature, 459, 1098
  * Radko (2013) Radko, T. 2013, Double-diffusive convection (Cambridge University Press)
  * Rappaport et al. (2008) Rappaport, N. J., Iess, L., Wahr, J., et al. 2008, Icarus, 194, 711
  * Rieutord (1987) Rieutord, M. 1987, Geophysical & Astrophysical Fluid Dynamics, 39, 163
  * Rieutord et al. (2001) Rieutord, M., Georgeot, B., & Valdettaro, L. 2001, Journal of Fluid Mechanics, 435, 103
  * Rieutord & Valdettaro (1997) Rieutord, M., & Valdettaro, L. 1997, Journal of Fluid Mechanics, 341, 77
  * Rosenblum et al. (2011) Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, The Astrophysical Journal, 731, 66
  * Rovira-Navarro et al. (2023) Rovira-Navarro, M., Matsuyama, I., & Hay, H. C. C. 2023, The Planetary Science Journal, 4, 23
  * Rovira-Navarro et al. (2019) Rovira-Navarro, M., Rieutord, M., Gerkema, T., et al. 2019, Icarus, 321, 126
  * Sohl et al. (2003) Sohl, F., Hussmann, H., Schwentker, B., Spohn, T., & Lorenz, R. 2003, Journal of Geophysical Research: Planets, 108
  * Sotin et al. (2017) Sotin, C., Hayes, A., Malaska, M., et al. 2017, in European Planetary Science Congress, EPSC2017–157
  * Spohn & Schubert (2003) Spohn, T., & Schubert, G. 2003, Icarus, 161, 456
  * Stevenson (1992) Stevenson, D. J. 1992
  * Tobie et al. (2006) Tobie, G., Lunine, J. I., & Sotin, C. 2006, Nature, 440, 61
  * Turcotte & Schubert (2002) Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge university press)
  * Tyler (2011) Tyler, R. 2011, Icarus, 211, 770
  * Tyler (2014) —. 2014, Icarus, 243, 358
  * Uchida et al. (2020) Uchida, T., Balwada, D., P. Abernathey, R., et al. 2020, Nature Communications, 11, 1125
  * Unno et al. (1979) Unno, W., Osaki, Y., Ando, H., & Shibahashi, H. 1979, Tokyo: University of Tokyo Press
  * Vance & Brown (2013) Vance, S., & Brown, J. M. 2013, Geochimica et Cosmochimica Acta, 110, 176
  * Waite et al. (2017) Waite, J. H., Glein, C. R., Perryman, R. S., et al. 2017, Science, 356, 155
  * Wolfram (1999) Wolfram, S. 1999, Assembly Automation
  * Wong et al. (2022) Wong, T., Hansen, U., Wiesehöfer, T., & McKinnon, W. B. 2022, Journal of Geophysical Research: Planets, 127, e2022JE007316