11institutetext: LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne
Université, Université de Paris, 5 Place Jules Janssen, 92195 Meudon, France.
11email<EMAIL_ADDRESS>22institutetext: Laboratoire de
Planétologie de Nantes, Nantes, France 33institutetext: Laboratoire de
Météorologie Dynamique, Paris, France 44institutetext: Johns Hopkins Applied
Physics Lab, Laurel, MD, USA

# Gravitational atmospheric tides as a probe of Titan’s interior: Application
to Dragonfly

B. Charnay 11 G. Tobie 22 S. Lebonnois 33 R. D. Lorenz 44

###### Abstract

Context. Saturn’s massive gravity is expected to causes a tide in Titan’s
atmosphere, producing a surface pressure variation through the orbit of Titan
and tidal winds in the troposphere. The future Dragonfly mission could analyse
this exotic meteorological phenomenon.

Aims. We analyse the effect of Saturn’s tides on Titan’s atmosphere and
interior to determine how pressure measurements by Dragonfly could constrain
Titan’s interior.

Methods. We model atmospheric tides with analytical calculations and with a 3D
Global Climate Model (the IPSL-Titan GCM), including the tidal response of the
interior.

Results. We predict that the Love numbers of Titan’s interior should verify 1
+ $\Re(k_{2}-h_{2}$) $\sim$ 0.02-0.1 and $\Im(k_{2}-h_{2}$) $<$ 0.04. The
deformation of Titan’s interior should therefore strongly weaken gravitational
atmospheric tides, yielding a residual surface pressure amplitude of only
$\sim$ 5 Pa, with a phase shift of 5-20 hours. Tidal winds are very weak, of
the order of 3$\times$10-4 m s-1 in the lower troposphere. Finally,
constraints from Dragonfly data may permit the real and the imaginary parts of
$k_{2}-h_{2}$ to be estimated with a precision of $\pm$0.01-0.03.

Conclusions. Measurements of pressure variations by Dragonfly over the whole
mission could give valuable constraints on the thickness of Titan’s ice shell,
and via geophysical models, its heat flux and the density of Titan’s internal
ocean.

###### Key Words.:

Planets and satellites: individual: Titan - Planets and satellites:
atmospheres - Planets and satellites: interiors

##  1 Introduction

Just as Saturn’s massive gravity causes tides both in Titan’s interior as well
as its surface seas, it causes a tide in the atmosphere (Lorenz, 1992; Tokano
& Neubauer, 2002; Strobel, 2006). Tokano & Neubauer (2002) analysed these
gravitational atmospheric tides with a 3D Global Climate Model (GCM). They
found that the tidal potential produced by Saturn would generate a surface
pressure variation of $\sim$1.3 hPa through the orbit of Titan and tidal winds
in the troposphere, with a mean wind speed of $\sim$0.3-0.4 m/s at 300 m.
Walterscheid & Schubert (2006) and Strobel (2006) studied these gravitational
tidal waves with an analytical model. They found that tidal winds should
increase with altitude until saturation in the upper atmosphere, where they
would deposit energy. Walterscheid & Schubert (2006) suggested that the
vertical transport by gravitational tides could produce the haze layers in
Titan’s upper atmosphere. However, these atmospheric studies did not take into
account the deformation of Titan’s interior. Titan likely possesses an
internal water-rich ocean, as suggested from the elevated Love number
$k_{2}\sim$0.62 measured by the spacecraft Cassini (Iess et al., 2012; Durante
et al., 2019). Such a large induced gravitational potential, implies a large
deformation of the interior which should strongly impact atmospheric tides.
Measuring them could provide additional constraints on Titan’s internal
structure.

Dragonfly, a relocatable lander for Titan (Lorenz et al., 2018), is presently
under development as NASA 4th New Frontiers mission, with a view to launch in
2027 and arrival in 2034. Its objectives include assessing Titan’s
habitability, and thus probing Titan’s internal structure. In this context,
this paper indicates a new means by which Dragonfly may address this topic
through measurements of pressure variations. In Section 2, we express the
tidal potential exerted by Saturn on Titan’s atmosphere, including the
deformation of the interior. We reanalyse gravitational atmospheric tides with
analytical calculations and with a 3D GCM in Section 3. Then, we discuss the
possibility to constraint Titan’s interior with Dragonfly in Section 4. We
finish with a summary and conclusions in Section 5.

##  2 Computation of Saturn’s tidal potential for Titan’s atmosphere

###  2.1 Expression of Saturn’s tidal potential

Titan orbits around Saturn with a synchronous rotation rate. Its large orbital
eccentricity ($e=0.0292$) leads to time-dependent tides. Tides are due both to
the time variation of the distance between Titan and Saturn (radial tide) and
to the time variation of Saturn’s angular position from the subsaturnian point
(librational tide). The tidal potential produced by Saturn can be expressed as
(Sagan & Dermott, 1982; Tokano & Neubauer, 2002):

$\begin{split}&V=-\frac{GM_{\rm S}}{a}\left(\frac{R_{\rm
T}}{a}\right)^{2}\left(\frac{3}{2}\cos^{2}(\lambda)\cos^{2}(\phi)-1\right)-3e\frac{GM_{\rm
S}}{a}\left(\frac{R_{\rm T}}{a}\right)^{2}\times\\\
&\left[\frac{1}{2}(3\cos^{2}(\lambda)\cos^{2}(\phi)-1)\cos(\Omega
t)+\cos^{2}(\lambda)\sin(2\phi)\sin(\Omega t)\right]\end{split}$ (1)

where $\phi$ is the longitude from the sub-saturnian point, $\lambda$ is the
latitude, $e$ = 0.0292 is Titan’s eccentricity around Saturn, $G$ =
6.673$\times$10-11 N m2 kg-2 is the universal gravitational constant, $M_{\rm
S}$ = 5.685$\times$1026 kg is the mass of Saturn, $R_{\rm T}$ = 2575 km is
Titan’s radius, $a$ = 1.222$\times$109 m is Titan’s semi major axis, $t$ is
the time (measured from periapsis, when Titan is at the nearest point of
Saturn), $\Omega$=4.56$\times$10-6 rad s-1 is Titan’s orbital angular
velocity.

The first term in (1) is time independent. This implies a permanent
modification of Titan’s geoid and interior with no impact on the atmospheric
dynamics. We focus on the time dependent terms of (1), which are related to
Titan’s eccentricity. The dynamic tidal potential can be expressed as:

$\begin{split}&V_{\rm
dyn}=\frac{-V_{0}}{2}(3\cos^{2}(\lambda)\cos^{2}(\phi)-1)\cos(\Omega t)\\\
&+V_{0}\cos^{2}(\lambda)\sin(2\phi)\sin(\Omega t)\end{split}$ (2)

with

$V_{0}=-\frac{GM_{S}}{a}\left(\frac{R_{T}}{a}\right)^{2}3\ e$ (3)

where $V_{0}$=-12.08 m2 s-2. We can express this potential as a superposition
of prograde (eastward), retrograde (westward) and stationary waves:

$\begin{split}V_{\rm
dyn}=-V_{0}\left[\frac{7}{8}\cos^{2}(\lambda)\cos(2\phi-\Omega
t)-\frac{1}{8}\cos^{2}(\lambda)\cos(2\phi+\Omega t)\right]\\\
+V_{0}\left(\frac{3}{4}\cos^{2}(\lambda)-\frac{1}{2}\right)\cos(\Omega
t)\end{split}$ (4)

or with normalized associated Legendre polynomials $P_{r,s}(\mu)$ and
$\mu=\sin\lambda$:

$\begin{split}&V_{\rm
dyn}=-\frac{7}{6}\sqrt{\frac{3}{5}}V_{0}P_{2,2}(\mu)\cos(2\phi-\Omega t)\\\
&+\frac{1}{6}\sqrt{\frac{3}{5}}V_{0}P_{2,2}(\mu)\cos(2\phi+\Omega
t)-\frac{V_{0}}{\sqrt{10}}P_{2,0}(\mu)\cos(\Omega t)\end{split}$ (5)

The first term is a diurnal eastward wave (wavenumber=2), the second term is a
diurnal westward wave (wavenumber=2) and the last term is a stationary wave.
The eastward wave dominates at low latitudes and the stationary wave dominates
at high latitudes. We can also express the dynamic tidal potential (2) as:

$V_{\rm dyn}=V_{1}\cos(\psi-\Omega t)$ (6)

with
$V_{1}=-V_{0}\sqrt{\frac{1}{4}\left(3\cos^{2}(\lambda)\cos^{2}(\phi)-1\right)^{2}+\cos^{4}(\lambda)\sin^{2}(2\phi)}$
and
$\psi=\arctan\left[\frac{2\cos^{2}(\lambda)\sin(2\phi)}{3\cos^{2}(\lambda)\cos^{2}(\phi)-1}\right]$.
This expression directly gives the local wave amplitude as shown in Fig. 1.
The tidal wave amplitude is maximal for $\lambda$ = 0∘ N and $\phi$ = $\pm$
32∘ (modulo 180∘). It is null for $\lambda$ = $\pm$ 60∘ N and $\phi$ = 0∘ /
180∘. It is equal to $V_{0}$ at the subsaturnian and anti-saturnian points.

Figure 1: Relative amplitude ($V_{1}/V_{0}$) and phase ($\psi$) of the dynamic
tidal potential as a function of longitude and latitude. The black cross
indicates Dragonfly’s landing site.

Tyler (2008) pointed out the presence of obliquity tidal forces on a moon with
a non-zero obliquity. In such a case, the obliquity tidal potential is given
by:

$V_{\rm
obliquity}=\frac{3}{2}\Omega^{2}R_{T}^{2}\theta_{0}\cos(\lambda)\sin(\lambda)\left(\cos(\phi-\Omega
t)+\cos(\phi+\Omega t)\right)$ (7)

where $\theta_{0}$ is the obliquity in radians. $\theta_{0}$=0.0052 rad for
Titan (0.3∘). Using Kepler’s third law ($\Omega^{2}=\frac{GM_{S}}{a^{3}}$),
the obliquity tidal potential can be expressed as:

$V_{\rm
obliquity}=-V_{0}\frac{\theta_{0}}{2e}\cos(\lambda)\sin(\lambda)\left(\cos(\phi-\Omega
t)+\cos(\phi+\Omega t)\right)$ (8)

The ratio of the amplitude of the obliquity tide by the amplitude of the
eccentricity tide is approximately $\frac{\theta_{0}}{2e}tan(\lambda)$. The
obliquity tide dominates at latitudes higher than 80∘ on Titan. It is thus
necessary to take into account the obliquity tide to study Titan’s sea tides
(see Tokano (2010a); Tokano et al. (2014); Vincent et al. (2016, 2018)). In
contrast, the eccentricity tide dominates at low latitudes. Its amplitude
should be $\sim$100 times that of the obliquity tide at Dragonfly’s landing
site (7∘ N 199∘W). We therefore neglect the effect of the obliquity tide
throughout this study, which focuses on Titan’s equatorial regions.

###  2.2 Influence of the deformation of Titan’s interior

The moment of inertia inferred from Cassini gravity measurements
($C/MR^{2}\simeq 0.34$, Iess et al. (2010)) indicates that Titan’s interior is
differentiated into a outer hydrosphere and an inner rocky core (Castillo-
Rogez & Lunine, 2010; Sotin et al., 2021). The relatively high value of
$C/MR^{2}$ suggest a low density and large rocky core ($>2050-2100$ km),
potentially containing a significant fraction of organics (Néri et al., 2020).
The determination of time variations of Titan’s gravity field due to Saturn’s
tides by Cassini (Iess et al., 2012; Durante et al., 2019) provided a clear
evidence that Titan possess a liquid layer underneath its icy surface, most
likely a water ocean (Mitri et al., 2014). The high value of the Love number,
$k_{2}$ = 0.616 $\pm$ 0.067, which quantify the amplitude of induced-
gravitational potential, suggested that the ocean is dense and salty (Mitri et
al., 2014). However, as the uncertainties on $k_{2}$ are still large, the
ocean depth and density still remains unconstrained. Due to the high pressure
reached at the bottom of the hydrosphere, high-pressure ice polymorphs (phase
V and VI) are expected to form (e.g. Journaux et al., 2020). The thickness of
this high-pressure layer is conditioned by the thickness of the outer ice
shell and the composition of the ocean. For a thin outer layer ($<50-70$ km),
it is possible that no high-pressure ice layer exists and that the ocean is in
direct contact with the rocky core, which would strongly enhance the
astrobiological potential of Titan’s internal ocean. Constraining the outer
ice shell thickness would therefore provide crucial constraints on the
habitability of Titan.

The deformation of Titan’s interior in response to Saturn’s tides results in
an change of the gravitational potential induced by time-varying internal mass
redistribution, proportional to the gravitational Love number $k_{2}$
mentioned above, together with a time-varying deflection of the surface,
proportional to the displacement Love number $h_{2}$. During high tides,
$k_{2}$ results in an increase of tidal potential exerted on Titan’s
atmosphere and hence induced surface pressure, while $h_{2}$ results in a
reduction due to surface elevation. If the interior behaves as perfectly
elastic body, the internal response would be in phase with the tidal potential
produced by Saturn (Eq.1). In reality, Titan’s interior is not perfectly
elastic, part of the mechanical energy is dissipated in the interior and
surface by various friction processes (Sohl et al., 1995; Tobie et al., 2005;
Tokano et al., 2014; Lorenz et al., 2014). This results in an out-of-phase
term which amplitude depends on the mechanical properties of the interior. The
in-phase term is represented by the real part of the Love number, whereas the
out-of-phase term corresponds to the imaginary part. As a consequence, the
tidal potential exerted on Titan’s atmosphere can be expressed as:

$V_{\rm atm}=V_{1}\left[(1+\Re(k_{2}-h_{2}))\cos(\Omega
t-\psi)-\Im(k_{2}-h_{2})\sin(\Omega t-\psi)\right]$ (9)

where $k_{2}$ and $h_{2}$ are the 2nd-degree complex Love number. $\Re()$ and
$\Im()$ are the real part of these complex numbers . The local amplitude of
the tidal potential is thus reduced by a factor
$\gamma_{2}=\sqrt{(1+\Re(k_{2}-h_{2}))^{2}+\Im(k_{2}-h_{2})^{2}}$ (Sohl et
al., 1995; Tokano et al., 2014; Lorenz et al., 2014) and has a phase shift
$\psi_{2}=-\arctan\left(\frac{\Im(k_{2}-h_{2})}{1+\Re(k_{2}-h_{2})}\right)$.

Pre-Cassini interior models predicted that $\gamma_{2}$ should be small, lower
than 0.2 and probably lower than 0.1 (Sohl et al., 2003, e.g.). Here we
compute the complex Love numbers $k_{2}$ and $h_{2}$ following the approach
used in Mitri et al. (2014) assuming a viscoelastic compressible interior
(Tobie et al., 2005) and considering interior structure models reproducing the
mean density (1881 kg/m$3$) and the reduced moment of inertia (0.341) of
Titan. We consider interior structures consisting of four main layers from
center to surface: a rocky core, a high-pressure (HP) ice V–VI layer, a liquid
water ocean, an ice I layer. The tidal response is computed using the
formulation of the spheroidal deformation developed by Takeuchi & Saito
(1972), initially derived for the elastic case, extended to the viscoelastic
case by solving it in the frequency domain and by defining complex shear and
bulk modulii, equivalent to the elastic modulii used in the elastic problem
(see Tobie et al. (2005) for more details). For the liquid ocean layer, the
simplified formulation of Saito (1974) is employed assuming a quasi-static and
non-dissipative fluid material. For computing the viscoelastic deformation of
the solid layers, we assume a compressible Andrade model, defined from bulk
modulus, $K$, shear modulus, $S$, viscosity, $\eta$ and two empirical
parameters $\alpha$ and $\beta$ describing the transient response of the
viscoelastic medium (Castillo-Rogez et al., 2011). Following Castillo-Rogez et
al. (2011), we assume the relationship $\beta\sim
S^{-(1-\alpha)}\eta^{-\alpha}$ and choose a reference value of 0.25 for
$\alpha$. The rheological parameters are assumed constant in each internal
layer, except in the outer ice layer where the viscosity, $\eta$, as a
function of depth is computed from a given temperature profile, which is
either conductive or convective In the conductive case, a temperature profile
varying linearly from the surface to the bottom of the ice shell is
considered. In the convective case, the ice shell is separated in a conductive
lid and a convective isothermal sublayer. The thickness of the conductive lid,
$b_{lid}$, is determined from an imposed surface heat flux and the viscosity
of the convective layer is determined from Nusselt–Rayleigh scaling laws
(Dumoulin et al., 1999) in order to reproduce a convective heat flux equal to
the imposed surface heat flux, $\phi_{S}$, varying between 15 and 30 mW.m-2
(see appendix A in Lefevre et al. (2014) for more details). In the conductive
part of the ice shell, the viscosity is computed from the linear temperature
profile using an Arrhenius law:

$\eta(z)=\eta_{T_{b}}\exp\left(\frac{E_{a}}{R(T(z)-T_{b})}\right)$ (10)

where $T_{b}$ is the temperature at the bottom of the conductive shell,
$E_{a}$ is the activation energy ($E_{a}=50$ kJ.mol-1), $R$ is the gas
constant.

Fig. 2 shows the computed values of $1+\Re(k_{2}-h_{2})$, $\Im(k_{2}-h_{2})$
and $\Re(k_{2})$, as a function of the ice shell thickness, the convective
heat flux and the ocean density. This figure suggests that accurate
measurements of $1+\Re(k_{2}-h_{2})$ and $\Im(k_{2}-h_{2})$ could constrain
the ice shell thickness, the thermal state of the ice shell (conductive vs.
convective) and potentially give some estimate of the internal heat flux. For
a conductive ice shell, $\Re(\gamma_{2})=1+\Re(k_{2}-h_{2})$ is directly
proportional to the ice shell thickness. For a convective ice shell, the
relationship between $\Re(\gamma_{2})$ and the ice shell thickness is less
straightforward. However, by assessing both $\Re(\gamma_{2})$ and
$\Im(\gamma_{2})$, the ice shell thickness and the vigor of convection can be
determined unambiguously. Once the ice shell thickness constrained from
$\gamma_{2}$ measurements, the measurements of $\Re(k_{2})$ can be used to
constrain the density of the ocean and hence its solute content.

Figure 2: Titan’s values of $1+\Re(k_{2}-h_{2})$, $\Im(k_{2}-h_{2})$ and
$\Re(k_{2})$, as a function of the ice shell thickness, the convective flux
(given in mW/m2) and the internal ocean density. The average ocean density was
fixed at 1200 kg/m3 for the top and the middle panel. The blue area in the
middle and bottom panels shows the values compatible with Cassini data at
1-sigma ($\Re(k_{2})$ = 0.549-0.683).

As demonstrated in the next section, the gradient of surface pressure
perfectly compensate the gradient of tidal potential, with almost no phase
shift. The surface pressure variation is thus equal to $-\rho_{0}V_{atm}$,
where $\rho_{0}$ is the density of Titan’s atmosphere at the surface. The
deformation of the interior would therefore reduce the amplitude of the tidal
pressure variations by $\gamma_{2}$, passing from $\sim$64 Pa at the sub-
saturnian point, to likely less than 6.4 Pa for $\gamma_{2}<0.1$. Fig. 3 shows
the local amplitude of the tidal pressure variations ($\Delta
P_{surf}=\gamma_{2}\rho_{0}V_{1}$) and phase compared to the sub-saturnian
point ($\psi+\psi_{2}$), without interior deformation (left) and with interior
deformation (right), using $1+\Re(k_{2}-h_{2})$ = 0.08 and $\Im(k_{2}-h_{2})$
= 0.01. These values correspond to an ice shell thickness of 80 km and a
convective heat flux of 15 mW/m2. We use these values as reference in the rest
of the paper. In this case, $\gamma_{2}$ = 0.081, $\psi_{2}$ = -0.124 rad
(corresponding to a time shift or -7.6 hours) and the amplitude of pressure
variation at the sub-saturnian point is equal to $\Delta P_{surf}$ = 5.2 Pa.

Figure 3: Amplitude of the tidal pressure variations (top panel) and phase
(bottom panel) as a function of longitude and latitude, without interior
deformation (left) and with interior deformation (right). The latter is
computed with 1+$\Re(k_{2}-h_{2}$)=0.08 and $\Im(k_{2}-h_{2}$)=0.01. The black
cross indicates Dragonfly’s landing site.

##  3 Computation of atmospheric tides

###  3.1 Analytical solution

We derived the analytical solution of gravitational atmospheric tides caused
by Saturn following the tidal theory from Chapman & Lindzen (1970) developed
for the solar thermal tides and the lunar gravitational tides on Earth. This
formalism was also used by Strobel (2006) for studying the development of
gravitational tidal waves in Titan’s upper atmosphere. In Titan’s troposphere
where winds are weak (geostrophic regime), gravitational atmospheric tides can
be described by the following equations (the linearised primitive equations
with no advective term) adding the gravitational tidal potential:

$\boxed{\begin{array}[]{rcl}\frac{\partial u}{\partial
t}-2\Omega\sin(\lambda)v=-\frac{1}{R_{\rm
T}\cos(\lambda)}\frac{\partial\left(\Phi+V_{\rm atm}\right)}{\partial\phi}\\\
\frac{\partial v}{\partial t}+2\Omega\sin\lambda u=-\frac{1}{R_{\rm
T}}\frac{\partial\left(\Phi+V_{\rm atm}\right)}{\partial\lambda}\\\
\frac{1}{R_{\rm T}\cos\lambda}\left(\frac{\partial
u}{\partial\phi}+\frac{\partial
v\cos\lambda}{\partial\lambda}\right)+\frac{1}{\rho_{0}}\frac{\partial\rho_{0}w}{\partial
z}=0\\\ \frac{\partial^{2}\Phi}{\partial z\partial t}+N^{2}w=0\end{array}}$
(11)

where $u$, $v$, and $w$ are the zonal, meridional and vertical winds,
$\rho_{0}$ is the air density, $\Phi$ is the geopotential, $p$ is the
pressure.
$N=\sqrt{\frac{g}{\theta}\frac{d\theta}{dz}}=\sqrt{\frac{g}{T}(\Gamma_{d}-\Gamma)}$
is the Brunt-Väisälä frequency ($\theta$ is the potential temperature,
$\Gamma$ and $\Gamma_{d}$ are the temperature gradient and the adiabatic
temperature gradient of the atmosphere). $N\sim$ 0.001 s-1 in the first 2 km
of Titan’s troposphere (Charnay & Lebonnois, 2012). The atmosphere is forced
by a tidal potential periodic in longitude and time. We search complex
solutions of the form:

$(u,v,w,\Phi,V_{\rm atm})=Re[(\hat{u},\hat{v},\hat{w},\hat{\Phi},\hat{V}_{\rm
atm})e^{i(s\phi+2\Omega\sigma t)}]$ (12)

with $\sigma$=1 and $s$=-1/2 for the eastward mode, $\sigma$=1 and $s$=1/2 for
the westward mode, and $\sigma$=1 and $s$=0 for the stationary mode. Solutions
satisfy the equations:

$\boxed{\begin{array}[]{rcl}\hat{u}=\frac{-\sigma}{2a\Omega}\mathcal{S}^{\sigma,s}_{u}(\hat{\Phi}+\hat{V}_{atm})\\\
\hat{v}=\frac{i\sigma}{2a\Omega}\mathcal{S}^{\sigma,s}_{v}(\hat{\Phi}+\hat{V}_{atm})\\\
\mathcal{L}^{\sigma,s}\hat{\Phi}-i\frac{2R_{\rm
T}^{2}\Omega}{\sigma\rho_{0}}\frac{\partial(\rho_{0}\hat{w})}{\partial
z}=-\mathcal{L}^{\sigma,s}\hat{V}_{atm}\\\
2i\Omega\sigma\frac{\partial\hat{\Phi}}{\partial
z}+N^{2}\hat{w}=0\end{array}}$ (13)

where $\mathcal{L}^{\sigma,s}$ is the Laplace tidal operator given by:

$\mathcal{L}^{\sigma,s}\hat{\Phi}=\partial_{\mu}\left(\frac{1-\mu^{2}}{\sigma^{2}-\mu^{2}}\partial_{\mu}\hat{\Phi}\right)-\frac{1}{\sigma^{2}-\mu^{2}}\left(\frac{s(\sigma^{2}+\mu^{2})}{\sigma(\sigma^{2}-\mu^{2})}+\frac{s^{2}}{1-\mu^{2}}\right)\hat{\Phi}$
(14)

$\mathcal{S}^{\sigma,s}_{u}$ and $\mathcal{S}^{\sigma,s}_{u}$ are the
horizontal wind operators given by:

$\begin{split}&\mathcal{S}^{\sigma,s}_{u}=\left[\frac{s}{\sigma^{2}-\mu^{2}}-\frac{\mu(1-\mu^{2})}{\sigma(\sigma^{2}-\mu^{2})}\frac{\partial}{\partial\mu}\right]\\\
&\mathcal{S}^{\sigma,s}_{v}=\left[\frac{s}{\sigma(\sigma^{2}-\mu^{2})}-\frac{(1-\mu^{2})}{(\sigma^{2}-\mu^{2})}\frac{\partial}{\partial\mu}\right]\end{split}$
(15)

with $\mu=\sin\lambda$. Eigenvectors of the Laplace tidal operator are the
Hough functions $\Theta_{n}^{\sigma,s}$ which verify:

$\mathcal{L}^{\sigma,s}\Theta_{n}^{\sigma,s}+\gamma_{n}^{\sigma,s}\Theta_{n}^{\sigma,s}=0$
(16)

The eigenvalues $\gamma_{n}^{\sigma,s}$ are called the Lamb parameters, also
defined with equivalent heights $h_{n}^{\sigma,s}$ with
$\gamma_{n}^{\sigma,s}=4\Omega^{2}R_{\rm T}^{2}/(gh_{n}^{\sigma,s})$

If we decompose the tidal potential and solutions in terms of Hough functions:
$\hat{V}_{\rm
atm}=\sum\limits_{n=1}^{\infty}V^{\sigma,s}_{n}\Theta_{n}^{\sigma,s}$ and
$(\hat{u},\hat{v},\hat{w},\hat{\Phi})=\sum\limits_{n=1}^{\infty}\left(\hat{u}_{n},\hat{v}_{n},\hat{w}_{n},\hat{\Phi}_{n}\right)\Theta_{n}^{\sigma,s}e^{z/2H}$,
the problem is reduced to the resolution of an equation on the vertical
structure:

$\frac{d^{2}\hat{w}_{n}}{dz^{2}}+\left(\frac{N^{2}}{gh_{n}^{\sigma,s}}-\frac{1}{4H^{2}}\right)\hat{w}_{n}=0$
(17)

with as boundary condition at the surface $\frac{D\Phi}{Dt}=0$ which can be
expressed as:

$\frac{d\hat{w}_{n}}{dz}+\left(\frac{1}{h_{n}^{\sigma,s}}-\frac{1}{2H}\right)\hat{w}_{n}=\frac{2i\sigma\Omega
V^{\sigma,s}_{n}}{gh_{n}^{\sigma,s}}$ (18)

To compute analytically the pressure and wind variations, we first computed
the Hough functions and Lamb parameters for given values of $\sigma$ and $s$,
using the algorithm from Wang et al. (2016). Then, we expanded the tidal
potential in Hough functions. We solved the vertical structure equation for
each mode n, assuming an upward wave propagation for
$(\frac{1}{h_{n}^{\sigma,s}}-\frac{1}{2H})\geq 0$, and an evanescent wave for
$(\frac{1}{h_{n}^{\sigma,s}}-\frac{1}{2H})<0$. Since $h_{n}^{\sigma,s}<<H$ for
all modes, the wave propagation only depends on the sign of
$h_{n}^{\sigma,s}$. Then we solved the surface boundary condition and we
determined $\hat{u}_{n}$, $\hat{v}_{n}$, $\hat{w}_{n}$ and $\hat{\phi}_{n}$ in
complex form. Finally we summed results over n to get the pressure and wind
fields. Fig. 4 shows the amplitude of the surface pressure variation for the
eastward and the westward tidal waves. The surface pressure perfectly follows
the tidal forcing for both cases. Fig. 5 shows maps of tidal surface pressure
and winds at a given time for the eastward and westward modes. The maximal
tidal wind speed is 3$\times 10^{-4}$ m/s, which is three orders of magnitude
lower than the typical surface wind speed predicted by 3D GCMs (e.g. Lebonnois
et al. (2012); Charnay & Lebonnois (2012)).

Figure 4: Amplitude of the atmospheric surface pressure for the eastward (top)
and westward (left) tidal wave compared to the tidal forcing. Analytical
calculation including interior deformation with $1+\Re(k_{2}-h_{2}$)=0.08 and
$\Im(k_{2}-h_{2}$)=0.01.

Figure 5: Tidal surface pressure, and tidal zonal, meridional and vertical
surface winds as a function of longitude and latitude at t=0 for the eastward
wave (left) and the westward wave (right). Analytical calculation including
interior deformation with $1+\Re(k_{2}-h_{2}$)=0.08 and
$\Im(k_{2}-h_{2}$)=0.01.

###  3.2 Simulations with the Titan IPSL GCM

We performed 3D simulation with the Titan IPSL GCM (Lebonnois et al., 2012)
including the tidal potential caused by Saturn. We used a horizontal
resolution of 32$\times$48 (corresponding to resolutions of 11.25∘ longitude
by 3.75∘ latitude) and 55 vertical levels. Simulations were run for 10 Titan
days at the vernal equinox with thermal tides and with/without gravitational
tides. Dragonfly will operate during the southern summer, but we did not
notice influence of the season on the amplitude of gravitational tides. We
used the tidal potential without or with interior deformation (with
$1+\Re(k_{2}-h_{2})$ = 0.08 and $\Im(k_{2}-h_{2})$ = 0.01). We analysed
simulations from day 5 to day 10, removing the first 5 days during which the
atmosphere readjusts to the presence/absence of tides.

Simulated pressure and horizontal wind variations confirm the analytical
calculations from the previous subsection. The pressure field follows the
tidal potential, and tidal winds are very weak in the simulations. Fig. 6
shows the time-evolution of the surface pressure at the sub-saturnian point
(longitude 0∘ and latitude 0∘) with thermal tides (diurnal solar cycle) and
with/without gravitational tides. For the case without interior deformation
(left panel), the gravitational tides dominate the surface pressure variations
which follow well the tidal potential. The standard deviation of the surface
pressure without gravitational tides or substracting the theoretical tidal
pressure is around 4 Pa. This is less than the amplitude of the tidal pressure
variations for the case with interior deformation (i.e. $\Delta$P = 5.2 Pa).
Although small, the tidal signal is clearly identified in the simulation with
interior deformation (right panel).

We found that the gravitational tides have a negligible effect on the surface
temperature (variations of $\sim 10^{-3}$ K correlated with surface pressure
for the case with deformation of the interior), which is dominated by the
diurnal solar cycle (peak-to-peak variations of $\sim$0.5 K). Fig. 7 shows the
surface pressure variations from thermal tides for the case without
gravitational tides. It was computed by making the direct substraction of a
simulation with thermal tides minus without thermal tides. It reveals that the
surface pressure variations from thermal tides are dominated by the diurnal
mode, as also found by Tokano (2010b). It has an amplitude of 1.7 Pa, smaller
than gravitational tides and with a phase shift (the maximum occurs at $\sim$6
pm solar local time). In contrast, the pressure variation from thermal tides
on Earth are dominated by the semi-diurnal mode (Chapman & Lindzen, 1970). It
is forced in particular by the solar energy absorption by stratospheric ozone
and has two pressure maximums at $\sim$10 a.m. and $\sim$10 p.m. solar local
time.

In our simulations, a major source of variability in the lower troposphere
comes from baroclinic waves, which develop at mid-latitudes with a period of
2-5 Titan days and a pressure variation of $\sim$20 Pa (Lebonnois et al.,
2012). However, their impact is limited in the equatorial region (pressure
variation of $\sim$1 Pa) and they can be distinguished from gravitational
tides by their long period. Mitchell et al. (2011) also found that
equatorially trapped Kelvin waves with a period of $\sim$0.5 Titan day would
occur. They are likely responsible for the large chevron-shaped methane storm
observed at the equator during the equinoctial season (Turtle et al., 2011).
We did not identify these modes in our simulations. The variability caused by
planetary waves in the equatorial region may thus be higher than predicted by
our model, in particular close to the equinox when equatorial clouds and
storms form.

Fig. 8 shows the zonal tidal wind component at 35 m above the surface for the
case with interior deformation. It was computed without thermal tides to
eliminate their effect and by making the direct substraction of a simulation
with gravitational tides minus without gravitational tides. From day 5 to day
10, the wind field is sufficiently close between the two simulations
(with/without gravitational tides) that the difference of surface wind
corresponds quite well to the tidal component. The simulated tidal wind
component is extremely weak, with an amplitude of $\sim 3\times 10^{-4}$ m/s
at the sub-saturnian point. This is much lower than typical surface winds
which are around 0.5 m/s in the IPSL GCM (Lebonnois et al., 2012; Charnay &
Lebonnois, 2012). This tidal wind speed is consistent with our analytical
calculations (see Section 2.2). Friedson et al. (2009) performed 3D
simulations of Titan’s atmosphere including gravitational tides with the
Titan-CAM model. They found an amplitude of tidal winds of around 0.05 m/s
compared to 0.5 m/s in the simulations from Tokano & Neubauer (2002). Our
simulated tidal winds are respectively two orders of magnitudes and three
orders of magnitudes lower than those in the simulations from Friedson et al.
(2009) and Tokano & Neubauer (2002). The deformation of the interior is
responsible for a reduction of tidal wind speed by one order of magnitude,
meaning that large discrepancies remain between the tidal responses of GCMs
for the same forcing. These discrepancies could be related to the mean
circulation or to the thermal structure. In particular, the tidal wind
amplitude is proportional to $N^{2}$ (the squared Brunt-Väisälä frequency) and
thus very sensitive to the thermal structure. The lower troposphere may be
more adiabatic in the IPSL GCM than in the models by Tokano & Neubauer (2002)
and Friedson et al. (2009), leading to weaker tidal winds. In this regard, the
Titan IPSL GCM reproduces well the thermal structure and the zonal wind
measured by the Huygens probe in Titan’s troposphere, especially in the
planetary boundary layer (Lebonnois et al., 2012; Charnay & Lebonnois, 2012),
giving confidence in our predictions. According to our analytical and
numerical calculations, tidal winds should not be detectable in Titan’s lower
troposphere.

Figure 6: Pressure variations at longitude 0∘ and latitude 0∘N computed with
the Titan IPSL GCM. The top (bottom) panel shows the effects of gravitational
tides without (with) interior deformation. The case with interior deformation
assumed $1+\Re(k_{2}-h_{2}$)=0.08 and $\Im(k_{2}-h_{2}$)=0.01. Thermal tides
were included in all simulations. The blue line is the top panel shows the
pressure variations without gravitational tides for comparison. The analytical
tidal pressure variation is shown with a black dashed line. Figure 7: Pressure
variations caused by thermal tides at longitude 0∘ and latitude 0∘N computed
with the Titan IPSL GCM. The analytical pressure variation from gravitational
tides with interior deformation is shown with a black dashed line. Figure 8:
Component of zonal tidal winds at 35 m above the surface computed at latitude
0∘N, and longitude 0∘W (red) or 90∘W (black). The lines corresponds to the
difference between simulations with tides and without tides. Thermal tides
were not included in these simulations.

###  3.3 Impact of tides on cloud formation

The tidal pressure variation and the associated adiabatic heating/warming
induces a change of relative humidity, which can impact cloud formation and
precipitation. A statistical analysis of 15 years of data from ERA-Interim
shows that the lunar semidiurnal tide induces a variation of relative humidity
and precipitation with an amplitude of 0.039$\%$ and 0.4$\%$ respectively on
Earth (Kohyama & Wallace, 2016). One may wonder if the stronger saturnian
gravitational tides could affect cloud formation on Titan, with a preferential
formation where the amplitude of tides are maximal (i.e. close to the
subsaturnian and anti-saturnian point). This mechanism has been suggested to
explain the two maximums of cloud occurrence from Cassini close to longitudes
0∘E and 180∘E (Rodriguez et al., 2009).

The passage of the tidal wave is associated with an adiabatic warming/cooling
:

$dT=-\frac{R^{*}T}{c_{p}P}dP$ (19)

with $R^{*}$ the specific gas constant, $c_{p}$ the specific heat capacity of
Titan’s atmosphere, T the temperature and P the pressure. Above the surface,
the amplitude of temperature variations caused by tides is $\sim 1\times
10^{-3}$ K. The relative humidity of methane clouds defined as
$RH=\frac{P_{vap}}{P_{sat}}$ evolves as:

$\frac{dRH}{RH}=\frac{dP_{vap}}{P_{vap}}-\frac{dP_{sat}}{P_{sat}}$ (20)

If we assume no condensation, $\frac{dP_{vap}}{P_{vap}}=\frac{dP}{P}$. The
Clausius-Clayperon relation gives
$\frac{dP_{sat}}{dT}=\frac{L_{CH_{4}}P_{sat}}{R_{CH_{4}}T^{2}}$, where
$L_{CH_{4}}$ is the mass-latent heat of methane and $R_{CH_{4}}$ is the
specific gas constant of methane. These relations lead to:

$\frac{dRH}{RH}=\frac{dP}{P}\left(1-\frac{L_{CH_{4}}R^{*}}{C_{p}R_{CH_{4}}T}\right)$
(21)

For Titan’s tropospheric conditions,
$\frac{L_{CH_{4}}R^{*}}{C_{p}R_{CH_{4}}T}\sim 3.1$. The relative humidity is
maximal when the tidal pressure is minimal. The maximal amplitude of relative
humidity variation is around 0.01$\%$. The effect is therefore very small,
even four times smaller than the lunar tides on the Earth. We conclude that
the influence of gravitational tides on Titan’s clouds should not be
detectable in Cassini data. Longitudinal variations of the cloud occurrence
rate should rather be due to variations of surface properties (topography,
surface humidity, albedo or thermal inertia).

##  4 Implications for the Dragonfly mission

###  4.1 Dragonfly’s ability to measure atmospheric tides

The Dragonfly mission is discussed in Turtle et al. (2018), Barnes (2021) and
Lorenz et al. (2018). In addition to the pressure/tide determination approach
in the present paper, Dragonfly may also constrain the thickness of Titan’s
ice shell via the Schumann resonance (Lorenz & Le Gall, 2020), and with
seismic methods (Lorenz et al., 2019). Note that Dragonfly (now planned to
launch in 2027) is still in an early design phase, and some of the following
details may evolve. The Dragonfly initial landing site (Lorenz et al., 2021)
is just south of the Selk crater, at 7∘ N 199∘W. This is fortunately rather
close to the anti-Saturnian point, where the pressure tide amplitude is nearly
at its maximum value (see Fig. 1).

The pressure measurement on Dragonfly’s DraGMet (Dragonfly Geophysics and
Meteorology) package was specified in order to measure the boundary layer
profile during flights up to 3-4 km. The absolute accuracy specification is
therefore only about 2 mbar (200 Pa). However, the sensitivity to changes will
be considerably better, and the telemetered resolution may be as good as 1 Pa.
The precision that can be practically achieved (and thus the ability to detect
small changes) will, as on other missions, be contingent on sensor noise
performance and will likely be period-dependent and lander-activity-dependent:
tests in Phase A with representative sensors suggest that the ability to
detect variations with an amplitude of 5 Pa on a range of periods should be
possible with suitable filtering. As discussed in Chapman & Lindzen (1970),
periodic changes can be detected in long time series data by virtue of
averaging over many cycles. The Dragonfly nominal mission is 3 years, or about
70 Tsols (Titan solar days). Crudely, diurnal signals may be measurable with a
precision of $\sim$700.5 $\sim$8$\times$ better than that of individual
measurements, offering good prospects that an atmospheric tidal signal can be
measured with a precision that can discriminate interior structure models. For
comparison, the pressure sensor on the Huygens probe had an accuracy lower
than 100 Pa and a resolution lower 1 Pa, quite similar to Dragonfly’s sensors.
The Huygens probe measured surface pressure on Titan for around 30 minutes
(Harri et al., 2006). The standard deviation of those pressure readings was
$\sim$1 Pa, indicating that any variation in background pressure must have
been slower than 2 Pa/hour. Unfortunately the measurement was not long enough
to constrain the tides. The maximal pressure variation (for the ideal case
with no deformation of the interior) due to tides at the Huygens site is
$\frac{dP}{dt}=-\Omega\rho_{0}V_{0}\sim$ 1 Pa/hour. However, the precision of
Huygen’s measurements holds promise for Dragonfly’s ability to detect tides.

It is probable that the largest confounding factor in retrieving a tidal
pressure signal is that lander activities will be correlated with the time of
(solar) day. Notably, rotorcraft flight and communication with Earth, both of
which will entail appreciable energy dissipation in the lander, and thus some
change in internal temperatures including the pressure sensor, will occur
during daytime. Thus any analysis of the sort described here must be vigilant
towards imperfections in the temperature compensation of the pressure
measurement (a known challenge on other missions - e.g. Taylor et al. (2010))
which may lead to a spurious signal with a period of one Tsol (i.e. the same
period as the atmospheric tide), and it is likely that the magnitude of such
effects will not be reliably known until Titan arrival. Over three years of
landed operations, however (a tenth of a Saturnian year), the phasing of the
gravitational tidal peak, fixed with respect to Titan’s orbit around Saturn
and thus in inertial space, will change with respect to local solar time (due
to the difference in length of the solar and sidereal day), which will
mitigate this problem. It should also be recognized that while it is intended
that DraGMet measure pressure data many times per Titan day (perhaps Earth-
hourly, with some bursts of higher-rate data), the record will have step-like
jumps in pressure owing to the different elevations of various landing sites.
Thus some manual editing and adjustment preprocessing of the time series may
be necessary to most precisely extract a tidal signal.

###  4.2 Constraints on Titan’s interior

As discussed in Section 2, the degree-2 Love numbers ($k_{2}$ and $h_{2}$) and
the surface pressure variations are related together by:

$\begin{split}&P_{\rm surf}=-\rho_{0}V_{1}(1+\Re(k_{2}-h_{2}))\cos(\Omega
t-\psi)\\\ &+\rho_{0}V_{1}\Im(k_{2}-h_{2})\sin(\Omega t-\psi)\end{split}$ (22)

Assuming that Dragonfly’s pressure measurements have a gaussian noise with an
uncertainty equal to $\sigma_{P}$, the uncertainty on the tidal pressure
amplitude would be equal to $\sqrt{\frac{2}{N_{obs}}}\sigma_{P}$ and the
uncertainty on $1+\Re(k_{2}-h_{2})$ and $\Im(k_{2}-h_{2})$ would be equal to
$\sqrt{\frac{2}{N_{obs}}}\frac{\sigma_{P}}{\rho_{0}V_{1}}$, where $N_{obs}$ is
the number of measurements (Alegria & Serra, 2006). Table 1 shows the
uncertainty for these parameters after 1 or 70 Titan days, assuming a
precision of either 50 Pa (optimistic case) or 200 Pa (pessimistic case) for
each individual measurement, done every hour (so 383 measurements per Tsol).
The pessimistic case corresponds to the absolute accuracy specification. We
have chosen a precision four times better for the optimistic case, able to
detect a variation of 4 Pa amplitude after one Tsol. This is fairly consistent
with initial tests of Dragonfly’s pressure sensors. Based on these two
scenarios, we expect that Dragonfly could detect tidal pressure variations and
could measure the real and imaginary parts of $(k_{2}-h_{2})$ with a precision
of $\pm$ 0.01-0.03 over the whole mission. This precision is comparable to the
precision on $k_{2}$ from Cassini flybys and will be complementary to the
$k_{2}$ estimate. According to Fig. 2, such values would allow us to derive
the ice shell thickness with a precision of around $\pm$ 15 km and to estimate
the heat flux with a precision of around $\pm$ 5 mW/m2, using interior models.
Unfortunately, the determination of the ocean density would be limited by the
precision on $k_{2}$. Some values of density could however be excluded. For
instance 1+ $\Re(k_{2}-h_{2})$ $\geq$ 0.08 would imply a density higher than
1250 kg.m-3.

Case | Optimistic | Pessimistic
---|---|---
Number of Tsols | 1 | 70 | 1 | 70
Amplitude (Pa) |  $\pm$4 |  $\pm$0.4 |  $\pm$14 |  $\pm$1.7
$\Re(h_{2}-k_{2}$) |  $\pm$0.06 |  $\pm$0.007 |  $\pm$0.2 |  $\pm$0.03
$\Im(h_{2}-k_{2}$) |  $\pm$0.06 |  $\pm$0.007 |  $\pm$0.2 |  $\pm$0.03
Table 1: Expected precision of Dragonfly measurements for the amplitude of
tidal pressure variations, for the real and the imaginary parts of
$k_{2}-h_{2}$. Values are computed for an optimistic case (precision of 50 Pa
for each individual measurement) and a pessimistic case (precision of 200 Pa
for each individual measurement). They are given for 1 Titan day or 70 Titan
days (i.e. the whole mission), assuming measurements every hour (383
measurements per Titan day).

##  5 Summary and conclusions

In this article, we reanalysed Titan’s gravitational atmospheric tides,
including the effects of the deformation of the interior. We showed that the
tidal response of the interior, characterized by low values of
$1+\Re(k_{2}-h_{2})$ and $\Im(k_{2}-h_{2})$, should strongly decrease the
tidal potential affecting the atmosphere. Using analytical calculations and 3D
GCM simulations, we showed that the atmosphere should quickly respond to the
tidal potential, with almost no phase shift and extremely weak tidal winds. We
predict that the amplitude of the pressure variation should be $\sim$5 Pa.
This low value mitigates the possible effect of gravitational tides for
transporting energy in the upper atmosphere, as suggested by Strobel (2006).
In addition, the impact of gravitational tides on cloud formation should be
totally negligible. Yet the tidal pressure variations could be higher than the
variations caused by thermal tides or planetary waves in the equatorial
region. The methane cycle (i.e. tropospheric cloud formation and
precipitation), which was not included in our 3D simulations, could be an
additional source of variability.

Finally, we predict that Dragonfly can certainly detect the pressure signal
that would be present if tidal deformation of the interior did not occur. It
may detect even the small residual tidal pressure variations taking
interior/crustal deformation into account. The real and imaginary part of
$(k_{2}-h_{2})$ could then be inferred with a precision of $\pm$ 0.01-0.03
over the whole mission, constraining the thickness of the ice shell and (via
models) the internal heat flux with a precision of $\sim$ $\pm$ 15 km and
$\sim$ $\pm$ 5 mW/m2 respectively. These measurements are insensitive to the
ocean density, but can restrict the range of possible ocean density in
combination with current $k_{2}$ estimate. A future orbiter around Titan could
strongly improve the estimate of $k_{2}$. Combined with Dragonfly
measurements, it could then constrain the density of the internal ocean, with
implications on its composition, evolution and habitability.

###### Acknowledgements.

This work was granted access to the HPC resources of MesoPSL financed by the
Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01)
of the programme Investissements d’Avenir supervised by the Agence Nationale
pour la Recherche. This work was supported by the Programme National de
Planétologie (PNP) of CNRS/INSU, co-funded by CNES. GT acknowledges the
support from the ANR COLOSSe project. RL acknowledges the support of the
Dragonfly project, via NASA contract NNN06AA01C to the Johns Hopkins Applied
Physics Laboratory. We are grateful to Henrik Kahanpää for a useful discussion
about Huygens’ pressure measurements, and we thank the referee, Tetsuya
Tokano, for his thorough review.

## References

  * Alegria & Serra (2006) Alegria, F. C., & Serra, A. C. 2006, in 2006 IEEE Instrumentation and Measurement Technology Conference Proceedings, 1643–1647
  * Barnes (2021) Barnes, J. e. a. 2021, The Planetary Science Journal, 2, 18
  * Castillo-Rogez et al. (2011) Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, Journal of Geophysical Research: Planets, 116
  * Castillo-Rogez & Lunine (2010) Castillo-Rogez, J. C., & Lunine, J. I. 2010, Geophysical Research Letters, 37
  * Chapman & Lindzen (1970) Chapman, S., & Lindzen, R. 1970, Atmospheric Tides (D. Reidel Publ. Co., Dordrecht, Holland), 200
  * Charnay & Lebonnois (2012) Charnay, B., & Lebonnois, S. 2012, Nature Geoscience, 5, 106
  * Dumoulin et al. (1999) Dumoulin, C., Doin, M.-P., & Fleitout, L. 1999, Journal of Geophysical Research: Solid Earth, 104, 12759
  * Durante et al. (2019) Durante, D., Hemingway, D. J., Racioppa, P., Iess, L., & Stevenson, D. J. 2019, Icarus, 326, 123
  * Friedson et al. (2009) Friedson, A. J., West, R. A., Wilson, E. H., Oyafuso, F., & Orton, G. S. 2009, Planet. Space Sci., 57, 1931
  * Harri et al. (2006) Harri, A.-M., Mäkinen, T., Lehto, A., Kahanpää, H., & Siili, T. 2006, Planet. Space Sci., 54, 1117
  * Iess et al. (2010) Iess, L., Rappaport, N. J., Jacobson, R. A., et al. 2010, science, 327, 1367
  * Iess et al. (2012) Iess, L., Jacobson, R. A., Ducci, M., et al. 2012, Science, 337, 457
  * Journaux et al. (2020) Journaux, B., Kalousová, K., Sotin, C., et al. 2020, Space Science Reviews, 216, 1
  * Kohyama & Wallace (2016) Kohyama, T., & Wallace, J. M. 2016, Geochim. Res. Lett., 43, 918
  * Lebonnois et al. (2012) Lebonnois, S., Burgalat, J., Rannou, P., & Charnay, B. 2012, Icarus, 218, 707
  * Lefevre et al. (2014) Lefevre, A., Tobie, G., Choblet, G., & Čadek, O. 2014, Icarus, 237, 16
  * Lorenz et al. (2019) Lorenz, R., Panning, M., Stahler, S., et al. 2019, in Lunar and Planetary Science Conference, Vol. 50, 2173
  * Lorenz (1992) Lorenz, R. D. 1992, in Proceedings of the Symposium on Titan, Toulouse, September 1991, 119–123
  * Lorenz & Le Gall (2020) Lorenz, R. D., & Le Gall, A. 2020, Icarus, 351, 113942
  * Lorenz et al. (2014) Lorenz, R. D., Kirk, R. L., Hayes, A. G., et al. 2014, Icarus, 237, 9
  * Lorenz et al. (2018) Lorenz, R. D., Turtle, E. P., Barnes, J. W., et al. 2018, Johns Hopkins APL Technical Digest
  * Lorenz et al. (2021) Lorenz, R. D., MacKenzie, S. M., Neish, C. D., et al. 2021, The Planetary Science Journal, 2, 24
  * Mitchell et al. (2011) Mitchell, J. L., Ádámkovics, M., Caballero, R., & Turtle, E. P. 2011, Nature Geoscience, 4, 589
  * Mitri et al. (2014) Mitri, G., Meriggiola, R., Hayes, A., et al. 2014, Icarus, 236, 169
  * Néri et al. (2020) Néri, A., Guyot, F., Reynard, B., & Sotin, C. 2020, Earth and Planetary Science Letters, 530, 115920
  * Rodriguez et al. (2009) Rodriguez, S., Le Mouélic, S., Rannou, P., et al. 2009, Nature, 459, 678
  * Sagan & Dermott (1982) Sagan, C., & Dermott, S. F. 1982, Nature, 300, 731
  * Saito (1974) Saito, M. 1974, Journal of Physics of the Earth, 22, 123
  * Sohl et al. (2003) Sohl, F., Hussmann, H., Schwentker, B., Spohn, T., & Lorenz, R. D. 2003, Journal of Geophysical Research (Planets), 108, 5130
  * Sohl et al. (1995) Sohl, F., Sears, W. D., & Lorenz, R. D. 1995, Icarus, 115, 278
  * Sotin et al. (2021) Sotin, C., Kalousová, K., & Tobie, G. 2021, Annual Review of Earth and Planetary Sciences, 49, 579
  * Strobel (2006) Strobel, D. F. 2006, Icarus, 182, 251
  * Takeuchi & Saito (1972) Takeuchi, H., & Saito, M. 1972, in Methods in Computational Physics: Advances in Research and Applications, Vol. 11, Seismology: Surface Waves and Earth Oscillations, ed. B. A. BOLT (Elsevier), 217–295
  * Taylor et al. (2010) Taylor, P. A., Kahanpää, H., Weng, W., et al. 2010, Journal of Geophysical Research (Planets), 115, E00E15
  * Tobie et al. (2005) Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534
  * Tokano (2010a) Tokano, T. 2010a, Ocean Dynamics, 60, 803
  * Tokano (2010b) —. 2010b, Planetary and Space Science, 58, 814
  * Tokano et al. (2014) Tokano, T., Lorenz, R. D., & Van Hoolst, T. 2014, Icarus, 242, 188
  * Tokano & Neubauer (2002) Tokano, T., & Neubauer, F. M. 2002, Icarus, 158, 499
  * Turtle et al. (2011) Turtle, E. P., Del Genio, A. D., Barbara, J. M., et al. 2011, Geochim. Res. Lett., 38, L03203
  * Turtle et al. (2018) Turtle, E. P., Barnes, J., Trainer, M., et al. 2018, in Lunar and Planetary Science Conference, Vol. 49
  * Tyler (2008) Tyler, R. H. 2008, Nature, 456, 770
  * Vincent et al. (2018) Vincent, D., Karatekin, Ö., Lambrechts, J., et al. 2018, Icarus, 310, 105
  * Vincent et al. (2016) Vincent, D., Karatekin, Ö., Vallaeys, V., et al. 2016, Ocean Dynamics, 66, 461
  * Walterscheid & Schubert (2006) Walterscheid, R. L., & Schubert, G. 2006, Icarus, 183, 471
  * Wang et al. (2016) Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geoscientific Model Development, 9, 1477