Resonant Stratification in Titan’s Global Ocean

Titan’s ice shell floats on top of a global ocean, as revealed by the large tidal Love number k 2 = 0.616 ± 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σ 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 (salinity S ≳ 200 g kg−1) can explain the 3σ 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 the 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.


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 ± 0.067 (Durante et al. 2019).The observation is an order of magnitude larger than the predicted k 2 ≈ 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 ≈ 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 ∼ 100 km (Sohl et al. 2003;Nimmo & Bills 2010;Luan 2019); an ice shell this thick reduces tides down to k 2 ≈ 0.42 (Section 2.1.3),which is roughly 3σ 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 ≈ 0.42.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 (ω s = 4.561 × 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 ω 2 ≲ 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 ≳ 10 −8 s −2 ≳ ω 2 s 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 ≳ + 1 g/kg every ∼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 ∼ 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 (salinity S ∼ 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,ice /k 2 × (ρ hp-ice − ρ)/ρ ∼ 2%, where k 2,ice /k 2 is the ratio between the Love number when the ocean is excluded/included, respectively, ρ hp-ice is the density of high-pressure ice, and ρ the ocean density.The previous estimate assumes that the tidal response of the ocean floor can be crudely approximated by the tidal response of an oceanless Titan model after introducing a correction for the relatively higher surface density of high-pressure ice.A viscous icy ocean floor can introduce a non-negligible component to k 2 if the viscosity of high-pressure ice is lower than ∼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 ∼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(Tyler , 2014;;Rovira-Navarro et al. 2019), with the additional heating slowing or halting the freezing of the ocean.

BASIC EQUATIONS
The tidal Love number k ℓm represents a normalization of the gravitational tidal response by the gravitational excitation, where the eccentric gravitational excitation ϕ e ℓm is derived in Appendix B and the tidal response ϕ ′ ℓm due to eccentricity tides is calculated numerically and analytically in this section in the context of various models.We concentrate in the ℓ = m = 2 Love number k 2 , which matches the spherical harmonic of the Cassini observation.We do not discuss obliquity tides in this manuscript because they excite m = 1 spherical harmonics that do not contribute to the Cassini k 2 observation.

Equilibrium tide
Here we reproduce previously known results of the tidal response of icy satellites with oceans using simple principles.Our objective is to illuminate the effects contributing to the amplitude of the equilibrium tide and to produce an estimate for Titan within a few percent of accuracy.We start by deriving the hydrostatic k 2 in a homogeneous fluid body.This derivation shows how all the tidal gravity in k 2 emerges from a thin layer of displaced fluid near the surface.However, Titan strongly departs from this estimate because its density profile is not homogeneous and its tidal response is not fully hydrostatic.Our second derivation introduces the effects of a nonhomogeneous density profile and an inner core that responds rigidly to diurnal tides (i.e., tidal frequency ω = ω s ) rather than hydrostatically.Finally, our last derivation shows how an elastic ice shell covering the ocean reduces the tidal response and provides a rough estimate that agrees with more sophisticated modeling.After considering all these effects acting jointly, the result is what we consider the diurnal equilibrium tide of Titan k 2 ≈ 0.42; a tide restored purely by static forces where inertial forces are neglected (i.e., no dynamics).

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 The tidal response of the icy satellite produces adiabatic perturbations represented in the gravitational potential ϕ ′ and pressure p ′ .In this expression, the potential of the gravitational pull ϕ e and the tidal gravitational potential ϕ ′ are combined into φ′ = ϕ e + ϕ ′ for analytical simplicity.
In an incompressible body, the density ρ remains constant regardless of forcing and the only change in the gravity field comes about from the radial displacement ξ 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 ξ r , where g is the surface gravity acceleration calculated in hydrostatic equilibrium g = 4πGρR/3, and G the gravitational constant.
The boundary condition on the surface indicates that the fluid is free from external pressure (i.e., δp = 0).This statement translates into turning the tidal displacement into the sole source for the pressure perturbation, formally expressed as Next, we can go back to equation ( 2) and obtain an expression for the tidal forcing, (5) We evaluate the tidal response and forcing at the surface (r = R) to obtain the tidal Love number from the ratio between them, 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.).

A global ocean with a rigid ocean floor
We now consider a body with a rigid core of density ρ c and radius R c overlaid by an ocean of density ρ 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πG ρR/3, where ρ is the body's mean density.The tidal gravity field still emerges solely from the radial displacement ξ rℓm of a thin region near the surface, Following the same steps as before, the resulting tidal Love number is independent of core properties, This equation yields k 2 = 0.468 when using Titan's mean density ρ ≈ 1.88 g cm −3 and ρ = 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 ≈ 0.249 when using ρ ≈ 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.).

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 ρ that is trying to reach the equipotential surface defined by the radial displacement ξ req .The pressure that the weight of the ocean exerts on the base of the ice shell can be expressed by where g is the gravity acceleration and ξ r is the radial tidal displacement.The shear stress τ on the ice shell can be expressed as where µ 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 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, which represents a competition between elastic and gravitational energy (see also equation ( 11) in Goldreich & Mitchell (2010)).An alternative derivation of equation ( 12), including all the correct numerical factors, can be found in Beuthe (2015).We may use Titan's radius R ≈ 2575 km, surface gravity acceleration g ≈ 135 cm s −2 , ocean density ρ ≈ 1 g cm −3 , and shear modulus of ice I and methane clathrates µ ∼ 4 GPa, and obtain An ice shell thickness d ∼ 100 km balances radiogenic heating with heat conduction (Appendix D) and produces k 2 ≈ 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 ∆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.Instead of being purely elastic, the ice shell covering the ocean can be viscous near the melting temperature.Viscous ice can in principle flow at certain timescale and reduce the resistance that the ice shell imposes on the tidally excited ocean motion.In practice, this effect is negligible at the tidal timescale ∼16 days unless large portions of the d ∼ 100 km ice shell thickness is in convection and thus have relatively low viscosity (∼10 14 Pa s).The compensated long-wavelength topography of Titan suggests that the ice shell is unlikely to be convecting (Nimmo & Bills 2010).

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.The equation of conservation of momentum and continuity are where v is the 3D incompressible tidal flow, Ω is the icy satellite spin vector, ρ is the ocean density, p is pressure, and γ 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.Ω = ω s ).
We derive the linearized form of the conservation of momentum equation by introducing linear perturbations of the form φ ≈ ϕ 0 + φ′ , where quantities with the zero subindex indicate the background state and primed quantities indicate perturbations due to tides.In the following, we drop the zero subindex for convenience, thus all nonprimed quantities indicate the background state.The resulting linearized equation of conservation of momentum is The timescale of tidal motion of ∼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 where Γ is the adiabatic index.We can relate the Lagrangian perturbations to the Eulerian perturbations in equation ( 17) using the definitions where ξ is the tidal displacement.Moving forward, we put together the latter equations to arrive to where N 2 is the Brunt-Vaisala (BV) frequency defined in general by 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 rN 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 ( 21), 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 Γ ≡ (∂ log p/∂ log ρ) S tending to infinity: the expected changes in pressure keep density unperturbed.As a consequence, the BV frequency reduces to 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 ( 23), stratification permits local density perturbations in the ocean regardless of its incompressibility.Equation ( 23) shows a direct relationship between stratification in N 2 and the salinity gradient in ∂ r ρ.When we consider constant stratification throughout the ocean, we can write ∂ r ρ = ∆ρ/H, where ∆ρ is the change in density between top and bottom of the ocean.Our models have constant density, thus the ∆ρ 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 where ω = ω + iγ is a complex frequency that accounts for tidal dissipation, ω = ω s is the tidal frequency of eccentricity tides, and the tidal displacement ξ is periodic with a time dependency ∝ e −iωt .We have used v = ∂ t ξ = −iωξ.This set of equations is traditionally known as the Boussinesq approximation.
A typical method of solution of equation ( 24) involves applying the curl to reach the vorticity equation (Rieutord 1987) 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 The free boundary at the ocean top prescribes a zero Lagrangian perturbation of pressure (e.g., Goodman & Lackner (2009)), which leads to The tidal displacement field ξ 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 ξ•r at the surface to be equal to the equilibrium tide radial 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, relaxing the no-slip assumption.No-slip boundary conditions explicitly set the dynamical part of ξ • 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 ℓ max = 100 and use N max = 100 Chebyshev polynomials to represent each radial function at each degree ℓ ≤ ℓ 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.

Tidal heating rate in the ocean
The heating rate in the ocean can be fully determined by volume integration of the work done by the dissipative force in the equation of motion (equation ( 14)).This work per unit volume follows where dl is an infinitesimal line element along a fluid parcel motion.We can use dl = vdt and integrate ẇ across the ocean volume to obtain the time averaged ocean heating rate (Chen et al. 2014;Rovira-Navarro et al. 2023, e.g.) The linear damping coefficient γ is unconstrained in icy worlds with global oceans.The range of possible estimates spans from γ ∼ 10 −11 s −1 on Enceladus to the better constrained γ ∼ 10 −5 s −1 on Earth (Matsuyama et al. 2018).The projection of equation ( 32) onto vectorial spherical harmonics is shown in Appendix C.

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 ∆k 2 to the non-resonant hydrostatic k 2 ≈ 0.42 obtained in a pure water ocean with a d ∼ 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 ( 23)), 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.Our models show that resonant stratification can explain the k 2 enhancement observed by Cassini.We observe an enhancement ∆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σ to 2σ away from the mean value of the observation (Fig. 2).In this case, we have used the conservative linear damping γ = 10 −9 s −1 , but our models predict a resonant ∆k 2 beyond +45% when using a still realistic γ = 3.3 × 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 ( 23)).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.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).

Surface radial displacement and interior nonlinear wave breaking
The dynamical enhancement ∆k 2 emerges from a tens-of-meter enhancement on the radial displacement of the ocean surface.In our models, this fractional enhancement to ocean surface radial displacement ∆h 2 equals ∆k 2 for any given combination of ocean parameters.We derive this result from combining equations ( 1 resulting in the relationship A direct inspection of equation ( 34) leads to ∆h 2 = ∆k 2 after substitution of the h ℓm and k ℓm that include a base value and a fractional correction.The Love number h ℓm describes the radial displacement of the ocean surface as a function of the equilibrium tide.A h ℓm = 1 implies that the surface follows the shape of the equilibrium tide when the self gravity of the tide is ignored.In a homogeneous fluid body, the average density is ρ = ρ and equation (34) recovers the classical result h 2 = 5/2.In the case of simple Europa interior models, equation (34) reproduces previous numerical results for h 2 when the core/mantle are approximately rigid (Moore & Schubert 2000).In general, equation ( 34) is valid for simple models where the tidal gravity originates entirely from the radial displacement of the surface, as is approximately the case for icy satellites with global oceans overlying roughly rigid rocky interiors.When we consider that Titan's equilibrium tide produces a surface displacement |ξ r | ≈ 26 m, the dynamical enhancement from resonant stratification shown in Fig. 3 leads to tides with surface displacement below |ξ r | ⪅ 40 m in the most dramatic case.We can calculate Titan's equilibrium tide Love number h ℓm using equations ( 8) and (33), Titan responds to Saturn's gravity with an equilibrium tide h 2 ≈ 1.47, which is reduced by −10% when a d ∼ 100 km ice shell is included (Section 2.1.3).We can obtain the radial displacement of Titan's equilibrium tide from equations ( 33), (B13), and (B24), 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 |ξ r | ∼ 1 km in the case of γ = 10 −9 s −1 (Fig. 3).A lower dissipation γ produces even larger resonant amplitudes of tidal motion interior to the ocean (Fig. 5), with γ = 10 −10 s −1 reaching |ξ r | ∼ 10 km.This large radial displacement is accompanied by a horizontal displacement that is typically ξ ⊥ /ξ r ∼ nR/H ∼ 50 (equation ( C31)) 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 γ ≳ 10 −10 s −1 , as stipulated in |ξk| ≲ 1, the typical criterion to avoid nonlinear wave breaking, where k is the wavenumber.For the radial displacement, this criterion translates to |ξ r | ≲ H/n.In the case of the horizontal displacement, the criterion stipulates ξ ⊥ ≲ π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.

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 ∼ (3 × 10 13 )/ Ėint km (Appendix D), where Ėint is the interior heating rate in W and d is the ice shell thickness.When the ocean freezes until balance with radiogenic heating, Ėint ∼ 3 × 10 11 W (Appendix D) and the ice shell thickness grows to d ∼ 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 γ = 10 −9 s −1 , for example, we get a total Ėint ∼ 6 × 10 11 W at saturation point (Fig. 4) and the ice shell thickness becomes d ∼ 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 ≲ 50 km for γ = 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 ∼ 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 ∼ 180 K if present in concentration ∼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 (D51)).
A lower γ = 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 γ leads to the damping of resonances, reducing the amplitude of the resonant ∆k 2 and the resonant Ė (Fig. 5).The hydrostatic flow is not much affected by γ, thus the heating rate is typically increased by a larger γ when far from resonances.Our results are relevant in the range γ = 10 −9 -10 −10 s −1 ; a γ larger than this range damps ∆k 2 below the signal registered by Cassini, whereas a γ lower than the range results in nonlinear breaking of internal gravity waves, an effect not included in our models.In addition to thermal steady state, ocean freezing requires an approach toward the resonance from the branch that produces a positive ∆k 2 (Fig. 2).Approaching the resonance from the negative ∆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 ( 23), the ocean stratification in this case follows 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 δ 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 δ, but increasing the density contrast δρ and consequently increasing N 2 (equation ( 23)), according to 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 ( 38) at constant δ) suggest the possibility of crossing the ∆k 2 resonance from the positive branch (Fig. 6).Further theoretical development is required to determine the effects of a diffusive layer in the ∆k 2 results reported here.
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(Tyler , 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 /π 2 κ ≈ 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.
The two ocean freezing scenarios described above assume a radially isotropic distribution of salts that linearly increases with depth over certain radial distance, either the ocean thickness H or the diffusive layer thickness δ.More complicated distributions of salts are in principle possible.The distribution of salts can show lateral variations in the presence of nonuniform ice shell thickness due to alternating regions of melting and freezing below the ice shell (Ashkenazy et al. 2018;Lobo et al. 2021;Kang et al. 2022;Kang 2023).Future investigations are required to better understand the impact of various distributions of salts in the dynamics of tidally excited internal gravity waves and the ocean freezing path that leads to a stable resonance.

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 δa ∼ 0.02R s , where R s is Saturn's radius (Fig. 7) and we have assumed the reasonable γ = 10 −9 s −1 (i.e., the resonance width depends on dissipation, as shown in Fig. 5).At Titan's current migration rate ȧ/a ∼ 10 −10 yr −1 (Lainey et al. 2020), the time it would take Titan to cover the g 6 -mode resonance width would be δa/ ȧ ∼ 10 Myr.
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 ∼3 times by orbital expansion over the lifetime of the solar system λ ∼ 4.5 × 10 9 yr, or Idini & Nimmo f r e e z i n g o c e a n : f u l l y s t r a t i fi e d freezing ocean: diffusive layer where P s is Titan's current orbital period, and ∆P s and ∆a represent orbital parameter changes on the timescale λ.The period spacing of g-modes (Appendix A) allows us to estimate the number of g-modes that Titan crosses over the timescale λ, where ∆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 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 × 10 −8 s −2 (see also Fig. 3).The tidal frequency is equal to the orbital frequency and ω 0 is Titan's current orbital frequency ω s .The eccentricity is fixed to Titan's current e and we assume synchronous rotation at all times (i.e., Ω = ω s ).
uniquely by orbital evolution (i.e., no stable fixed-point) is We obtain # g ∼ 5 and P(g-mode) ∼ 1% when using the reasonable ocean parameters in Fig. 7. Equation ( 41) is only valid for # g ≳ 2. The low P(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.

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 γ = 10 −9 s −1 , resonant stratification yields a k 2 value 1σ closer to the measured central value starting from the hydrostatic k 2 = 0.42.When γ = 3.3 × 10 −10 s −1 , the enhancement over the hydrostatic k 2 is 3σ, 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 ∼ 100 − 200 g/kg higher to produce a similar effect on k 2 , depending on the exact γ (Fig. 8).Previous studies have argued in favor of a heavy convective ocean that holds S ∼ 200 g/kg in salts to obtain a 1σ agreement with k 2 observations.However, water-rock interactions at the bottom of Titan's ocean are expected to produce a limited S ≲ 10 g/kg of salts (Leitner & Lunine 2019).
Instead of water-rock interactions, the heavy convective ocean relies on a special event to attain its relatively high salt concentration.Ammonium sulfate ((NH 4 ) 2 SO 4 ) 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 ∼ 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 ∼ 10 g/kg) to enable resonant stratification.At S ∼ 10 g/kg salinity, the predicted k 2 is in a 2σ agreement with the Cassini observation when γ = 10 −9 s −1 and within 1σ agreement when γ = 3.3 × 10 −10 s −1 (Fig. 8).Both γ values are realistic and prevent nonlinear wave breaking.This bulk concentration of salts (S ∼ 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).

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 across an ocean thickness H ∼ 300 km to preserve the stability of the mild stratification N 2 ∼ 1 × 10 −8 s −2 discussed before, where α ∼ 2 × 10 −4 K −1 is the thermal expansivity and g ≈ 135 cm s −2 is the surface gravity.This temperature gradient is equivalent to a heat transfer ≲2 × 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 ∆T ∼ αHgT /c p ∼ 5.2 K (Turcotte & Schubert 2002) over an ocean with the same properties used in equation ( 42), where c p ∼ 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 to maintain Ledoux stability (equation ( 42)) of the ocean salinity gradient, where κ 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 ∼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 κ that allows a salinity gradient with N 2 ∼ 10 −8 s −2 to remain Ledoux stable (equation ( 43)).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 ν is typically larger than the thermal diffusivity κ.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.

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).

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 ± 0.067, which is 3 − σ 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 (γ) 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.
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 ) 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.

B. THE TIDAL FORCING OF ECCENTRICITY TIDES
The gravitational potential ϕ T experienced by an observer at r from a concentrated mass located at r ′ , is inversely proportional to the distance between them where α is the angle between the two position vectors, ℓ is degree, and P ℓ are the m = 0 case of the associated Legendre polynomials 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 The two lowest degree harmonics, ℓ = 0 and ℓ = 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 α = cos θ cos θ p + sin φ sin φ p cos(φ − φ p ), (B7) where rθφ 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 ) where we use the conventional definition of spherical harmonics ) From fundamental identities of spherical harmonics, it can be shown that The tidal forcing potential for ℓm then becomes In the standard case of a coplanar circular orbit, we have cos θ p = φ p = 0, leading to with the normalization constant 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, The resulting tidal excitation potential in a planar orbit is 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 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.

Ψ. (C42)
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 ≈ 7.4 × 10 −12 W/kg on Earth's mantle (Turcotte & Schubert 2002).The resulting radiogenic heat production is Ėint ∼ 0.5Hm s ∼ 5 × 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 ≈ 4.5 × 10 −12 W/kg (Grasset & Sotin 1996;Spohn & Schubert 2003), which results in the smaller Ėint ∼ 3 × 10 11 W.
The balance between internal radiogenic heating Ėint and conduction across the ice shell determines the ice shell thermal equilibrium thickness.Titan's icy surface is at T s ∼ 90 K.In a pure water ocean, we can use k ∼ 2 W m −1 K −1 , Titan radius R ≈ 2575 km, and ∆T ∼ 180 K to obtain (Luan 2019) d ∼ 4πR 2 k∆T Ėint ∼ 3 × 10 13 Ėint km. (D51) 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 ∼ 180K in an ocean with 15%wt.ammonia (Grasset & Sotin 1996).In this scenario, the temperature gradient across the ice shell diminishes to ∆T ∼ 90 K and the equilibrium ice shell thickness can be reduced in half compared to the pure water ocean (equation (D51)).

Figure 1 .
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).

Figure 2 .
Figure 2. Titan's k 2 enhancement (fractional correction ∆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 γ = 10 −9 s −1 .The reference hydrostatic Love number is k 2 = 0.42, a pure water ocean with an elastic ice shell with thickness d ∼ 100 km.The tidal frequency is equal to the rotational frequency and the orbital frequency, ω = Ω = ω s = 2π/T orb , where T orb = 15.945days is Titan's orbital period.

Figure 3 .
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 γ = 10 −9 s −1 .

Figure 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 Ėint ∼ 3 × 10 11 W (see Apprendix D).The frictional damping is γ = 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.).

Figure 5 .
Figure 5. Resonant dynamical gravity signal and heat production as a function of dissipation for g-mode g 6 (ℓ = 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).
Approach to a long-lived resonance via ocean freezing

Figure 6 .
Figure 6.Model parameters that produce tidally excited resonances with internal gravity waves (solid lines; see equation (A2)).The positive branch of the ∆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 ≈ 8.3 × 10 21 g, which is equivalent to a δρ ∼ 0.001 g cm −3 over an ocean thickness H ∼ 200 km.
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 × 10 −8 s −2 (see also Fig.3).The tidal frequency is equal to the orbital frequency and ω 0 is Titan's current orbital frequency ω s .The eccentricity is fixed to Titan's current e and we assume synchronous rotation at all times (i.e., Ω = ω s ).

a
= a 0 (1 − e cos ωt), (B14) φ p = 2e sin ωt.(B15) Our strategy is now to expand equation (B11) to first order in e.The semimajor axis dependency expands of the planet expands as e −imφp ≈ 1 − 2ime sin ωt.≈ 1 − em(e iωt − e −iωt ).(B21) of the fact that the direction of tides can be flipped by either changing the sign in φ or the tidal frequency ω s , reducing the eccentricity tidal potential toϕ e ℓm ≈ e (ℓ + 1 + 2m) U ℓm r R ℓ Y m ℓ (θ, φ)e −iωt .(B24)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 ℓ = m = 2.C.PROJECTION OF DYNAMICAL TIDES ONTO VECTORIAL SPHERICAL HARMONICSWe project our equations onto vectorial spherical harmonics following the standard decompositionξ = y 1 Y + y 2 Ψ + y 3 Φ,(C25)where YΨΦ constitute an orthonormal base for the projection of vectorial fields in spherical polar coordinates.VSH relate to scalar spherical harmonics (equation (B9)) asY = Y r, (C26) Ψ = r∇Y , (C27) Φ = r∇ × Y,(C28)where spherical harmonics satisfy r 2 ∇ 2 Y = −ℓ(ℓ + 1)Y .We further project the gravity response and pressure-gravity potentials, respectively,