Icarus 326 (2019) 123–132

Contents lists available at ScienceDirect

Icarusjournal homepage: www.elsevier.com/locate/icarus

Titan's gravity field and interior structure after Cassini
a,⁎

b,c

a

a

d

Daniele Durante , D.J. Hemingway , P. Racioppa , L. Iess , D.J. Stevenson

Ta

Department of Mechanical and Aerospace Engineering, Sapienza University of Rome, Italy
Department of Earth and Planetary Science, University of California, Berkeley, CA, USA
Department of Terrestrial Magnetism, Carnegie Institution for Science, Washington, DC, USAdDivision of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USAbc

ARTICLE INFOABSTRACT

Keywords:
Geophysics
Orbit determination
Radioscience
Titan interior

Since its arrival at Saturn in 2004, Cassini performed nine flybys devoted to the determination of Titan's gravity
field and its tidal variations. Here we present an updated gravity solution based on the final data set collectedduring the gravity-dedicated passes, before Cassini's plunge into Saturn's atmosphere. The data set includes anadditional flyby (T110, March 2015, primarily devoted to imaging Titan's north polar lakes) carried out with the
low-gain antenna. This flyby was particularly valuable because the closest approach occurred at a high latitude(75°N), over an area not previously sampled.
Previously published gravity results (Iess et al., 2012) indicated that Titan is subject to large eccentricity tides
in response to the time varying perturbing potential exerted by Saturn. The magnitude of the response quadrupole field, expressed in the tidal Love number k2, was used to infer the existence of an internal ocean. The new
gravity field determination provides an improved estimate of k2 of about 0.62, accurate to a level of a few
percent. The value is higher than the simplest models of Titan suggest and the interpretation is unclear; possibilities include a high density ocean (as high as 1300 kg/m3), a partially viscous response of the deeper region,or a dynamic contribution to the tidal response. The new solution includes higher degree and order harmoniccoefficients (up to 5) and offers an improved map of gravity anomalies. The geoid is poorly correlated with the
topography, implying strong compensation. In addition, the updated geoid and its associated uncertainty could
be used to refine the gravity-altimetry correlation analysis and for improved interpretation of radar altimetricdata.

1. Introduction

the other instruments that require pointing. Also, during gravity passes,
thruster firings, commonly used for attitude maneuvers, were absent inorder to preserve the coherency of Cassini's orbit, allowing an accurate
reconstruction of the dynamical system. In addition to the nine gravitydedicated Titan passes, during the high-altitude T110 pass the gravity
data were acquired using the omnidirectional Low-Gain Antenna (LGA)
while Cassini's Visual and Infrared Mapping Spectrometer (VIMS) tooka high-resolution regional map of Titan's north polar lakes.
The analysis of the first four gravity flybys (Iess et al., 2010) reported accurate measurements of Titan's gravity field up to degree 3.
The inclusion of two additional flybys (Iess et al., 2012) shed light on
Titan's tidal response to eccentricity tides, revealing its large deformability over timescales of the orbital period, an indication of a global
ocean at depth. Here, we update Titan's gravity field and its internal
structure by including the data collected during three additional
gravity-dedicated flybys and the T110 pass.
The main details of the passes when Cassini collected gravity data

On 15 September 2017, Cassini plunged into Saturn's atmosphereand completed its long-lasting mission. After 7 years spent in deep
space and 13 years orbiting Saturn, the spacecraft collected an incredible amount of scientific data about the Saturnian system. One of
the main objectives of the mission was the study of Titan, Saturn'slargest moon. Also, Titan was crucial for navigating Cassini around thesystem because it was exploited to perform large orbital correction
maneuvers with low propellant cost. For this reason, Cassini performed
124 flybys of Titan during the mission (labelled with a T, which stands
for Titan, followed by the number of the pass, e.g., T011), nine of whichwere entirely dedicated to the study of its gravity field, with the aim ofinferring its interior structure. During the gravity-dedicated Titan
passes, Cassini pointed its High-Gain Antenna (HGA) toward Earth,ensuring the best possible data quality. Since Cassini lacks a scan
platform, this attitude does not simultaneously allow the operation of⁎

Corresponding author.
E-mail address: daniele.durante@uniroma1.it (D. Durante).

https://doi.org/10.1016/j.icarus.2019.03.003
Received 16 November 2018; Received in revised form 28 February 2019; Accepted 6 March 2019
Available online 08 March 20190019-1035/ © 2019 Elsevier Inc. All rights reserved.

Icarus 326 (2019) 123–132

D. Durante, et al.

are reported in Table 1. Each pass has a different orbital geometry, but
the spatial coverage of Titan is not uniform (see Fig. 1a); the equatorial
region of the planet has been sampled more than the polar regions. To
determine Titan's tidal response to Saturn at its orbital frequency, themean anomaly is the key quantity (Iess et al., 2012). Cassini's samplingof Titan's mean anomalies is very broad (see Fig. 1b) and allows us tofurther constrain Titan's tidal behavior. The Sun-Earth-Probe (SEP)
angle is very important in determining the Doppler noise in the radio
link, due to the effect of solar plasma noise turbulence (Asmar et al.,2005; Iess et al., 2014). The Doppler noise budget is influenced also by
a daily variability due to the troposphere conditions at the groundstation. For these reasons, the Doppler noise varies from a minimum of
0.017 mm/s during T011 to a maximum of 0.068 mm/s during T045
(see Table 1).

arc. The initial values for Cassini's state vectors are taken from the latest
release of the navigation team for each pass, with a priori uncertainties
on position and velocity of, respectively, 2 km and 2 cm/s in all the
direction, large enough not to constrain the solution. Particular attention should be paid to Titan's state vectors and the details are discussed
below. We combine the information from all the arcs to solve for Titan's
gravity field coefficients up to degree and order 5, its gravitationalparameter, the tidal Love number k2, the density of Titan's atmosphereat a reference altitude and its scale height, and the spacecraft bodyfixed acceleration from the onboard RTGs. All these quantities have apriori uncertainties which do not constrain the estimate. The solar radiation pressure acceleration is not estimated but considered1 with anuncertainty of 5%. The uncertainty on Titan's rotation model from SAR
imaging (Meriggiola et al., 2016) is much lower than that achievablewith gravity measurements and thus the rotational parameters have notbeen estimated.
The analysis of T110 involves one additional challenge, related tothe use of Cassini's LGA. Since the main purpose of this Titan flyby was
to take a high-resolution regional map of Titan's north polar lakes with
Cassini's VIMS, the HGA was not pointed toward the Earth, and thephase-coherent radio link was established through the LGA. As thespacecraft attitude changes over time, the phase center of the LGA
moves with respect to the spacecraft's center of mass, therefore generating a large Doppler signal which must be removed (Barbaglio et al.,2012). As a result, the location of the phase center of this antenna withrespect to Cassini's center of mass must be included in the estimation
process.
An important issue in the analysis of Cassini's gravity data is theneed to update Titan's ephemeris. During a typical Titan flyby, Cassini's
positioning with respect to Titan has a formal uncertainty of at most a
few meters in the radial direction during Titan's closest approach (according to orbital geometry and available data). The analysis of the first
four passes (Iess et al., 2010) was carried out by adopting the latest
ephemeris available from the navigation team whereas the first six
gravity passes (Iess et al., 2012) were analyzed by updating Saturnian
ephemerides (in particular Titan) in order to match the observed
Doppler data with Titan's updated ephemeris and gravity field. The
updated ephemerides had to be dynamically consistent between all thepasses. The inclusion of the latest data does not allow us to obtain asatisfactory fit of the Doppler data together with an update of Titan's
ephemeris. This may be related to the need for a more complete dynamical model for the Saturnian satellites due to the long timespancovered by Cassini's observations (10 years). The inclusion of effects
such as Saturn's tidal dissipation (Fuller et al., 2016) or frequency dependent tides, Saturn's normal modes, or an accurate modelling of Saturn's pole motion may be crucial to obtain a dynamically consistentmotion of Saturn's satellites to the level of accuracy required by thedata. Thus, the approach taken hereafter is to update Titan's ephemerisfor each of Cassini's passes, without requiring the updated ephemeris tobe dynamically coherent between the passes. This approach allows us toovercome the difficulties inherent to requiring a global ephemeris, allowing for deviations in Titan's orbit. Of course, the larger number of
estimated parameters implies a penalty on the uncertainties of thegravity coefficients. We use the latest available Saturn's satellite ephemerides, SAT389 (see NASA's Navigation and Ancillary Information
Facility website), as a reference ephemeris. Titan's position has an apriori uncertainty of 5 m in the radial direction, and 50 and 500 m in
the along-track and out-of-plane directions (Boone et al., 2017). The
limited dataset cannot provide a more accurate estimate of Titan's orbit.
The required corrections are of the same order of magnitude as theassumed a priori uncertainties. Our reconstruction of Titan's orbit isvalid only in the timeframe of each arc with the intent of obtaining an

2. Data analysis
Titan's updated gravity field exploits the full data set, which includes mostly two-way Doppler data. After Cassini's Ka-band translator
failed in 2003, two- or three-way Doppler data are available only with
X-band uplink (7.2 GHz) and X or Ka-band downlink (respectively 8.4
and 32.0 GHz). The X-band uplink makes the radio link especiallysensitive to charged particle effects, which increase dramatically when
the SEP angle is small (i.e., near solar conjunctions; see Asmar et al.,2005; Iess et al., 2014). In the case where multiple data types areavailable in the same time interval, the data selection privileges twoway Doppler data over three-way data, and Ka-band over X-band in
downlink. The three-way data are about 12% of the total Doppler
points. We discarded data acquired when the elevation angle of thespacecraft, as viewed from the ground station, is below 20° to prevent
bias in the solution that may come from uncalibrated troposphere. We
weighted the data from each tracking pass according to the observed
noise RMS. Each observation arc spans about 2 to 5 days about the timeof Cassini's closest approach (C/A) to Titan, and avoids orbit correction
maneuvers and thruster firing (used to desaturate Cassini's reaction
wheels). Range data have not been included because they do not contribute substantially to the determination of Titan's gravity field; also,
since range data are more affected by systematic errors, we prefer asolution with Doppler data only. However, range data have been used
afterward to test the robustness of the solution based on Doppler data
only. While, as expected, range data do not change the estimatedgravity coefficients, they are used to update the ephemerides of Saturnin the solar system.
The Doppler data have been compressed at 60 s and analyzed with
MONTE, the JPL orbit determination code (Evans et al., 2018). Cassini's
dynamical model includes the gravitational effect in a relativistic context due to the Solar System planets including Saturn and its satellites.
For Saturn, the zonal gravity field up to degree 10 has been included, inagreement with the dynamical model used for Titan's ephemerides (seeJacobson et al., 2006 and subsequent JPL releases). Titan's gravity field
coefficients are included up to degree and order 5. The rotational elements have been adopted from Meriggiola et al. (2016), obtained fromthe analysis and georeferencing of Cassini RADAR images. We included
non-gravitational accelerations from the solar radiation pressure andthe anisotropic thermal emission produced by the three onboard
Radioisotope Thermoelectric Generators (RTGs). In addition, low-altitude flybys (in particular T022 and T068) require the inclusion of an
exponential model for the density of Titan's upper atmosphere since thedrag acceleration produced at those altitudes is large enough to perturb
Cassini's trajectory in a non-negligible way (the peak acceleration in
T022 is of the order of 3 × 10−6 m/s2).
The Doppler data from the different passes are collected in a multiarc, least-squares estimation filter, which provides the corrections to
the dynamical model's parameters and, after a few iterations, allows us
to fit the Doppler data to the noise level. In each arc, we solve for both
Cassini's and Titan's state vectors at the beginning of each observation

1
Consider parameters are those which are not estimated but whose a priori
uncertainty contributes to the covariance matrix.

124

Icarus 326 (2019) 123–132

D. Durante, et al.

Table 1
Summary of Titan flybys during which Cassini collected gravity data. For each pass, the date of Cassini's closest approach (C/A), the minimum altitude reached, thelatitude and longitude of the C/A, the Sun-Earth-Probe (SEP) angle, the number of Doppler points (both two- and three-way), and the root mean square (RMS) of
Doppler residuals at 60 s are reported.

T011
T022
T033
T045
T068
T074
T089
T099
T110
T122

Date of C/A

AltitudeLatitudeLongitude

SEP angle

Number of Doppler points

RMS of Doppler residuals (at 60 s)

27 Feb. 2006
28 Dec. 2006
29 June 2007
31 July 2008
20 May 201018 Feb. 201117 Feb. 201306 Mar. 201416 Mar. 2015
10 Aug. 2016

1812 km
1297 km1933 km1614 km1397 km
3651 km1978 km1500 km
2274 km
1698 km

0.2°S
45.4°N
8.4°N
43.5°S
48.9°S
1.0°N
21.0°N
31.1°S
74.8°N
12.4°N

255.6°E
355.9°E
63.1°E
162.7°E
241.1°E
113.4°E
203.1°E
181.0°E
263.1°E
234.4°E

150.2°
130.6°
46.4°
30.0°
120.3°
131.2°
106.0°
111.0°
108.6°
112.7°

1051
1151
1508
1288
2226
1847
1479
1860
1836
1858

0.017 mm/s
0.024 mm/s
0.028 mm/s
0.068 mm/s
0.025 mm/s
0.063 mm/s
0.033 mm/s
0.022 mm/s
0.064 mm/s
0.061 mm/s

Fig. 1. (a) Cassini's ground tracks on Titan, for 3 h around the closest approach time. The location of the C/A is indicated with circles. The equatorial region is very
well mapped, whereas the polar regions have more limited coverage. The flybys at the highest latitude are T110 for the northern hemisphere, and T068 and T099 forthe southern hemisphere. The underlying image is a mosaic of Titan's surface from Cassini's ISS team (NASA/JPL-Caltech/SSI). (b) Titan's mean anomaly at the timeof Cassini's closest approach. A good sampling is crucial for the determination of Titan's tidal response.

unbiased estimate of the gravity field.

Table 2Estimated values and 1-σ uncertainties for Titan's gravitational parameter (GM),unnormalized dimensionless gravitational potential coefficients (whereJl = − Cl0) corresponding to a reference radius of 2575 km, and the tidal Love
number k2 for the reference solution.

3. Gravity field results
Titan's gravitational potential, evaluated at radius r, colatitude θ,and longitude ϕ, can be expressed as a linear combination of sphericalharmonic functions as

GM (km3/s2)
6l

U (r , , ) =GM
Rref l(Clm cos m + Slm sin m ) Plm (cos )r l=0 m =0 r
(1)where Clm and Slm represent the unnormalized degree-l and order-mdimensionless gravitational potential coefficients at some reference
radius (which we take to be Rref = 2575 km), and where Plm(cosθ) are
the unnormalized associated Legendre functions. G is the universal
gravitational constant and M is Titan's mass. Note that all Sl0 = 0 andthat C00 ≡ 1. Because the gravitational potential is resolved about thecenter of mass, the degree-1 terms are all zero.
The gravity field coefficients estimated for Titan are reported inTable 2 along with Titan's tidal Love number k2 whereas Fig. 2 reports
the root mean square of the fully normalized gravity field coefficients
Cl, showing that the power spectrum of the gravity coefficients is larger
than the uncertainty up to degree 5. Indeed, truncating the gravity field
expansion to a smaller degree cannot provide a good fit of all Doppler
data. On the other hand, higher degrees cannot be determined with
sufficient accuracy and have been neglected. However, only a few
coefficients are determined with good accuracy, such as the degree 2hydrostatic terms and the sectorial harmonics of degree 3 and 4. Still,
an expansion up to degree and order 5 of the gravity field is an improvement with respect to the previous solution (Iess et al., 2012),

8978.1383 ± 0.0003

J2 (×10 )C21 (×106)
C22 (×106)
J3 (×106)
C31 (×106)
C32 (×106)
C33 (×106)
J4 (×106)
C41 (×106)
C42 (×106)
C43 (×106)
C44 (×106)
J5 (×106)
C51 (×106)
C52 (×106)
C53 (×106)
C54 (×106)
C55 (×106)

33.089 ± 0.609
0.513 ± 0.215
10.385 ± 0.084−0.179 ± 0.720
1.481 ± 0.254
0.183 ± 0.153−0.222 ± 0.017−1.077 ± 1.844−0.842 ± 0.299
0.183 ± 0.107−0.012 ± 0.039−0.014 ± 0.003
1.118 ± 2.022
0.361 ± 0.406−0.097 ± 0.118−0.016 ± 0.019
0.007 ± 0.004
0.000 ± 0.001

k2

0.616 ± 0.067

S21 (×106)
S22 (×106)

0.612 ± 0.359−0.064 ± 0.066

S31 (×106)
S32 (×106)
S33 (×106)

0.811 ± 0.402−0.027 ± 0.099−0.226 ± 0.019

S41 (×106)
S42 (×106)
S43 (×106)
S44 (×106)

0.191 ± 0.717
0.198 ± 0.106−0.062 ± 0.033−0.012 ± 0.004

S51 (×106)
S52 (×106)
S53 (×106)
S54 (×106)
S55 (×106)

0.267 ± 0.604
0.044 ± 0.094−0.004 ± 0.012−0.002 ± 0.004
0.000 ± 0.001

which resolved Titan's gravity field only up to degree 3 (SOL1a andSOL2) or degree 4 (SOL1b). The uncertainties on the gravity field
coefficients reported in Table 2 are generally up to a factor of 2 larger
than those reported by Iess et al. (2012), even if a broader dataset has
been used. This apparent discrepancy lies in the complexity of the125Icarus 326 (2019) 123–132

D. Durante, et al.

reduces the volume in which the relevant density anomalies arise), andthe lower strength of ice (which reduces the peak deviatoric stresses
that are assumed to be universal in the scaling argument given above).
Unfortunately, Titan is the only icy satellite whose gravity field has
been confidently resolved up to degree 5, and a comparison within thisclass of bodies is not possible at present. Note also that terrestrial planets as a group are actually much more similar to each other (in respectof a core and mantle structure) than the large icy satellites, and that
Europa, a body of much current interest, is in neither class (it has far
less ice than rock).
The direction of Titan's principal axes of inertia can be inferred fromthe degree 2 coefficients. For a synchronous rotator, the axis of least
inertia is expected to be aligned approximately with the empty focal
point of Titan's orbit around Saturn (Murray and Dermott, 1999). A
small 3-2-1 rotation can be defined to relate the body-fixed reference
frame (a frame whose z-axis aligns with the spin pole and x-axis aligns
with the prime meridian; we adopted the frame defined in Meriggiola
et al., 2016) to the principal axes reference frame. The three rotations
are related to the fully normalized gravity coefficients (e.g., Liu andChao, 1991) by:

Fig. 2. Root mean square of the fully normalized gravity coefficients,22
m (Clm + Slm ). The values and 1-σ uncertainties are reported and compared
2l + 1
to a Kaula's rule, Cl = K/l2, with K = 10−5.

S21= 2.85° ± 1.67°
3 C20 + C22
C21= 0.55° ± 0.23°
y =3 C20 C22
S22
= 0.18° ± 0.18°
z =
2 C22

Cl =

x =

adopted dynamical model: firstly, the maximum rank of our gravity
solution is larger, meaning that more coefficients are included in theestimation process; secondly, the approach selected for Titan's ephemeris is different, with our approach leading to about a 20% increase of
the uncertainties, especially in the low-degree field. The tidal Love
number k2, however, is determined with a better accuracy due to thebroader sampling in longitude of Titan's mean anomaly.
Fig. 2 plots the power law Cl = K/l2, known as the Kaula rule
(Kaula, 1963), with K = 10−5, showing good agreement with the estimated power of Titan's gravity field for degrees larger than 2 (J2 andC22 are dominated by the hydrostatic component). The Kaula rule is anempirical law able to represent the power spectrum of the gravitational
field of rocky planets, but it has never been observed for an icy satellite
like Titan. Note that a Kaula constraint has not been imposed into theestimation process. This empirical rule conveys very little information
about the physical nature of the density anomalies responsible for thefield, since a power law like this can arise in a number of physicallyplausible ways and is seen for all terrestrial bodies (including Earth)despite the substantial differences among these bodies that are known
to exist by other methods, most obviously the fact that Earth has plate
tectonics. A fit of the gravity coefficients (excluding degree 2) yields
Cl = 0.97 × 10−5/l2. The coefficient (K ≈ 10−5) might convey more
physical information than the spectral slope since it is affected by the
degree of compensation, the volume of the region responsible for thedensity anomalies, and the typical stress levels that exist within the
body. For example, if we suppose that there is a peak deviatoric stressthat scales as Δρ × g × R (where Δρ is a typical density anomaly, g is
gravitational acceleration, and R is radius or some comparable large
length scale) then we might expect Cl ∝ (R × g)−1 ∝ g−2 since g scales as
R. Note however that this scaling should not be used when going fromrocky bodies to icy bodies since the expected peak stresses are lower in
ice and the mean density is lower, so the partial success of this scaling
for terrestrial bodies does not imply it should work for icy bodies. Recently, Ermakov et al. (2018) reviewed the Kaula rule of different terrestrial bodies in the literature to achieve a more general scaling ruleamong them. Their result for the gravity spectra is Cl ∝ g−1.72. However,
if we were to apply this suggested rescaling rule to Titan, we obtain an
overestimate of Titan's gravity spectrum by a factor of 30 to 50 (Titan's
surface gravity is 0.14 times Earth's gravity). This is not surprisingconsidering the large degree of compensation (discussed further below),
the closely related likelihood of a rather thin outer shell of ice (which(2)

all compatible with zero at most at 3-σ. The direction of the axis of least
inertia is compatible at 1-σ with the expectation for a synchronous
rotator.
Due to the non-negligible eccentricity of Titan's orbit (e = 0.028),
tides raised by Saturn change with time, producing a variable quadrupole component. Titan's tidal response to Saturn can be described byits k2 tidal Love number, which is crucial to understanding the interior
structure of Titan. Cassini's flybys dedicated to gravity science occurredat different Titan mean anomalies (see Fig. 1b), probing the variable
quadrupole field at different locations, and providing a very accurate
measurement of its degree-2 tidal Love number. The best estimate ofTitan's tidal response is k2 = 0.616 ± 0.067, which is within 1-σ of the
all the three solutions given by Iess et al. (2012), but has a smaller
uncertainty due to the larger data set used for its estimation. However,
the accuracy on k2 is somehow degraded due to the local ephemeris
approach used for Titan's orbit (a covariance analysis shows a degradation of about 20% in the uncertainty with respect to a global Titan
ephemeris). However, the consistency with the previous solution confirms Titan's large deformability over timescales of its orbital period,
consistent with a global ocean beneath an outer ice shell that is sufficiently thin that it does not greatly impede the deformation that would
be expected for a hydrostatic fluid overlying the core. The implicationsof this value are discussed further in Section 4.
The reference solution does not account for the tidal dissipation,and only the real component of the Love number k2 is estimated. By
including the time lag into the tidal model, the gravity coefficients remain compatible at 1-σ with our reference solution. The estimate for thetime lag is Δt = 14,000 ± 20,000 s, which corresponds to an angularlag angle of δ = 3.7° ± 5.4° or to an imaginary part of the complex
Love number of k2im = 0.082 ± 0.118. The estimated quality factor
can be inferred from the real and imaginary part of the tidal Love
number k2 from Q = k2re/k2im. The inverse of the quality factor is 1/
Q = 0.133 ± 0.192, implying Q ≈ 7.6. The time lag is compatible withzero at 1-σ, thus we can draw no strong conclusions about the internaldissipation from Cassini gravity data.
The robustness of our reference solution has been tested, with theaim of assessing whether the formal uncertainties we quoted in Table 2
are statistically meaningful and to confirm that biases in the estimated126

Icarus 326 (2019) 123–132

D. Durante, et al.

over the other, we therefore find the moment of inertia factor thatyields, via the Radau-Darwin relation (see references above), hydroandC22hydthatminimizestaticvalues
J2hyd

hyd 2
obs
(J2obs J2hyd )2 + (C22C22
) . The resulting best fitting moment ofinertia factor is approximately 0.341, corresponding to a potential fluidLove number, k2f, of approximately 1.01, with J2hyd = 33.21 × 10−6
and C22hyd = 9.96 × 10−6 (cf. Table 2).
Note that, whereas for fast spinning bodies like Enceladus, the hydrostatic ratio becomes slightly smaller than 10/3 as higher ordercorrections become important (Tricarico, 2014), the first order approximation of 10/3 is highly accurate for Titan (to within 0.02%).
Another correction can account for Titan's orbital eccentricity(Matsuyama, 2011), giving a correction smaller than 0.3%. Thus, we
can state that for a hydrostatic Titan, the value of 10/3 is correct towithin 0.3%, ten times below the current accuracy. Our estimate of theratio is J2/C22 = 3.186 ± 0.077, compatible with the hydrostatic value
at 2-σ. Though, as we discuss further below, having the 10/3 ratio is anecessary, but not sufficient condition for hydrostatic equilibrium.
In the simplest terms, Titan may be regarded as a mainly rocky coresurrounded by an H2O mantle consisting of various solid and liquid
phases. Adopting a simple two-layer model (rocky core, H2O mantle),
taking the moment of inertia factor to be 0.341, and assuming the
average density across the H2O layers is 1000 kg/m3, the known radius(2575 km) and bulk density (1881 kg/m3) can be used to solve for thecore radius and density (e.g., see eqs. 1 and 2 in Hemingway et al.,2018). Under these assumptions, we obtain an H2O layer thickness of~400 km, leaving a large rocky core with a density around 2500 kg/m3.
For the sake of generality, we carry out this calculation over a wide
range of different moments of inertia and H2O layer densities (Fig. 5),
allowing for the possibility of, for example, the higher ocean densities
permitted by the large tidal k2 (Iess et al., 2012; Stiles, et al., 2008;
Meriggiola et al., 2016; and see below). Such higher density outerlayers would imply a significantly smaller and denser core, even when
adopting the same moment of inertia. Moreover, the possibility of a
non-hydrostatic interior (see below) permits moments of inertia smaller
than 0.341, corresponding to still smaller and denser cores (Fig. 5).
The measured J2 and C22 tell us that the hydrostatic, or fluid, k2
Love number is close to unity (see above). This is substantially larger
than the tidal k2, which tells us that some major portion of Titan does
not behave like a fluid at tidal timescales. The value for the tidal k2 for athin shell and rigid interior beneath the ocean consistent with themoment of inertia is about 0.42 (Rappaport et al., 2008; Iess et al.,2012). This assumes that the outermost shell has little effect on the
response; any effect it could have because of its elasticity would only
decrease the predicted k2. It follows that there are only two candidates
for increasing the theoretical value to the observed value: either increase the ocean density (which also decreases the core size to becompatible with the moment of inertia) or allow the core to have asignificant non-elastic (i.e., viscous) response at tidal timescales. Some
combination of these two is possible of course. The interpretation is
non-unique but the choice of a viscous core response (which is discussed in Iess et al., 2012) requires a non-Maxwell like response since a
Maxwell model that simultaneously increases the real part of k2 also
increases the imaginary part to a level inconsistent with the data. The
core need not exhibit a Maxwell rheology, so this possibility is not
excluded. As discussed by Iess et al. (2012) and considered by
Rappaport et al. (2008) and Baland et al. (2014), an ocean density of
1200–1300 kg/m3 is cosmochemically possible and yields a k2 of 0.55
to 0.6. We lack sufficient understanding of the rheology of deeper regions to assess quantitatively the possible contribution of a non-elastic
response for that region. It is also possible that the outer shell is reducing the Love number (which makes it harder to explain the result weobtained with a dense ocean alone) and a dynamic response to the tidal
perturbation has also been suggested. Most interest in dynamic tides
has centered on the possibility of additional tidal heating (Matsuyama

Fig. 3. Comparison of estimated values and 1-σ uncertainties for J2, C22, and k2
between our reference solution and other solutions: SAT337 and SAT375 update a different version of Saturn's satellite ephemerides, whereas SAT389 has
3-times larger a priori uncertainty for Titan's position and velocity; DRAG estimate a drag coefficient for each of the low-altitude passes instead of Titan'sbase density and scale height; IAU 2000 and IAU 2009 adopt the rotation
models suggested by IAU.

values are absent. Since the selection of the reference ephemeris andthe a priori covariance given to Titan are crucial in our approach, arobustness test has been performed by starting from older versions of
the ephemerides (e.g., SAT337 and SAT375) or by enlarging the a priori
covariance by a factor of 3 (see Fig. 3). All the solutions are statisticallyconsistent at 1-σ, confirming the robustness of Titan's gravity field solution with respect to a change in Titan's ephemeris. Concerning Titan's
atmosphere, we verified that the estimation of drag coefficients in lowaltitude flybys instead of Titan's base density and scale height yields astatistically compatible solution (within 1-σ, see Fig. 3). We tested also
the impact of the assumed rotational model by adopting rotation
models advised by the International Astronomical Union. The models
given in Seidelmann et al. (2002) and that proposed by Stiles, et al.
(2008) and later adopted by IAU (Archinal et al., 2011; Archinal et al.,2018) give gravity fields statistically consistent at 1-σ with our reference solution (see Fig. 3).
4. Implications for Titan's interior
Titan's gravitational field conforms closely to the expectation for abody that has relaxed to hydrostatic equilibrium, with a moment ofinertia factor in the vicinity of 0.341 (Fig. 4a; see below). The momentof inertia factor is given by I/MR2, where I is the mean moment ofinertia, M is Titan's mass, and R is Titan's mean radius. For a relaxedhydrostatic body, the departure from spherical symmetry arises mainly
from rotational and tidal deformation and the associated redistributionof mass gives rise to the large J2 and C22 terms in the gravitational
potential. For a relaxed hydrostatic body, these terms have a theoreticalratio of J2/C22 = 10/3, to first order (the dashed line in Fig. 4a). The
moment of inertia values indicated along the dashed line in Fig. 4a canbe computed via the Radau-Darwin equation (Darwin, 1899; Murray
and Dermott, 1999), which relates the moment of inertia factor to thefluid Love number k2f, which can in turn be related to J2 and C22 (e.g.,
see eqs. 16, 21–23 in Hemingway et al., 2018). If Titan were perfectlyhydrostatic, one could thus obtain k2f, and therefore the moment ofinertia factor, directly from either of the J2obs or C22obs terms, where thesuperscript ‘obs’ refers to the observed values (Table 2). Because of thedeparture from the hydrostatic expectation, however, J2obs and C22obs
correspond to slightly inconsistent moment of inertia factors. In order
to approximately satisfy both observations, with no preference for one
127Icarus 326 (2019) 123–132

D. Durante, et al.

0.360
0.358
0.356
0.354
0.352
0.350
0.348
0.346
0.344
0.342
0.340
0.338
0.336
0.334
0.332
0.330

38
37

J2 (x106)

36
35
34
33
32
31
30

b 400

Hydrostatic
Iess+ 2012
This work

9

9.5

10

10.5

11

350

-H 2,0 (meters)

a 39

11.5

C 2,2 (x106)

300

250
0.40

200

150

0.30

100
40

45

0.32

0.34

50

0.36

55

0.38
Hydrostatic
Zebker+ 2009 triaxial
Zebker+ 2012 (l<=6)Corlies+ 2017
This work geoid

60

65

70

75

80

85

H 2,2 (meters)

Fig. 4. (a) The newly measured (red) dominant quadrupole gravitational potential coefficients, J2 and C22, in comparison to the results of Iess et al. (2012) (blue),with 1-σ uncertainty ellipses, and (b) the corresponding degree-2 terms of Titan's figure, -H20 and H22, in relation to the measured geoid (small purple ellipse,
highlighted by purple arrow); three distinct shape models are shown for comparison (blue: the triaxial ellipsoidal solution of Zebker et al., 2009; red: the l ≤ 6solution of Zebker et al., 2012; and yellow: Corlies et al., 2017). In both panels, the observations are shown in relation to the expectation for hydrostatic equilibrium(dashed line, with numbers indicating normalized mean moments of inertia, computed via the Radau-Darwin relation). Whereas the shape departs substantially fromthe expectation for hydrostatic equilibrium, the gravitational field departs only slightly from the hydrostatic expectation. (For interpretation of the references to color
in this figure legend, the reader is referred to the web version of this article.)of inertia factor of 0.341. Because the measured gravitational potential
is so close to the hydrostatic expectation, this hypothetical hydrostatic
equilibrium figure is nearly identical to the geoid (small purple ellipse
in Fig. 4b). The geoid is a surface where the potential, including tidal
and rotational effects, has a constant value of GM/Rref. The geoid can beapproximated to first order as

4500
900 kg/m 31000 kg/m 3

Core Density (kg/m 3)

0.33
4000

1100 kg/m 3
1200 kg/m 31300 kg/m 3
1400 kg/m 3

0.30

3500

l

N ( , ) = Rref

0.31

0.34
3000

0.32

2500
0.36
2000
1400

1600

1800

2
5 2Rref

0.33

0.35

2000

(Clm cos m + Slm sin m ) Plm (cos )l=0 m=0

6g

0.34
0.35
0.36
0.37
0.37
0.38
0.38
0.39
0.40
0.39
0.40
2200
2400
2600

P20 (cos ) +

2R 2
ref

4g

P22 (cos ) cos 2(3)where the last two terms are required to account for the combined tidal
and rotational potentials and where ω and g are Titan's rotation rate andsurface gravity, respectively. We do not show the topography with respect to the geoid because the result is visually indistinguishable fromthe non-hydrostatic topography shown in Fig. 6d. Thus, measured withrespect to the geoid, or almost equivalently, to a hydrostatic equilibrium figure with moment of inertia factor 0.341, Titan's topography isrelatively high in low latitudes—a notable exception being the topographic low near the Xanadu province—and is lowest near the poles,
especially on the trailing hemisphere (Fig. 6d).
The non-hydrostatic topography complicates the simple picture
above because it should also contribute to the gravitational field, andyet there appears to be no significant non-hydrostatic gravity, at least atdegree 2 (Fig. 4a). This is a strong indication that the > 500 m of nonhydrostatic topography are highly compensated, meaning that thegravity anomalies associated with this topography must be offset byinternal density anomalies. This is expected in the limit of a weak lithosphere, because the topography can only be supported isostatically(i.e., by buoyancy). The presence of a liquid water ocean beneath thefloating ice shell provides a natural mechanism for Airy-type isostaticcompensation in the sense of lateral shell thickness variations (e.g.,
Nimmo and Bills, 2010) wherein the low density isostatic roots offset
the effect of the mass corresponding to the surface topography. It has
also been suggested that Pratt-type isostatic compensation could play arole in the form of hydrocarbon-cycle-related density anomalies near
Titan's polar regions (Choukron and Sotin, 2012), which would likewise
reduce the contribution to the measured gravity field made by the nonhydrostatic topography.

Core Radius (km)Fig. 5. Range of possible core radii and densities in relation to a range ofmoments of inertia (constant along dashed curves) assuming various meandensities for the H20 mantle (solid curves). A two-layer body is assumed (rocky
core, H2O mantle) and is subject to the constraints of Titan's known total radius(2575 km) and bulk density (1881 kg/m3).

et al., 2018), which is not desirable for explaining Titan, but non-resonant dynamic tides might be responsible for a significant contribution
to the Love number.
Unlike its gravitational field, Titan's figure (Zebker et al., 2009;
Corlies et al., 2017) departs substantially from the hydrostatic expectation, exhibiting considerable excess flattening. Using the latest
shape model (Corlies et al., 2017), the dominant degree-2 terms havethe ratio −H20/H22 = 5.5 ± 1.2, far larger than the hydrostatic valueof 3.333 (Fig. 4b). As with the gravitational potential, we compute thehydrostatic equilibrium figures, shown along the dashed line in Fig. 4b
for various moments of inertia, again using the Radau-Darwin relation(e.g., see eqs. 18, 21–23 in Hemingway et al., 2018). Whereas Fig. 6bshows the Corlies et al. (2017) topography with respect to a reference
sphere, with radius 2575 km, Fig. 6d shows the topography with respect
to the hydrostatic equilibrium figure expected for an assumed moment
128

-5

0

45°E

5
2.

2.5

45°E

90°E

90°E

135°E

.5

-2

135°E

180°E

0

180°E

0

225°E

225°E

-7
.5

270°E

270°E

8

315°E

15
12.5
10
7.5
5
2.5
0
-2.
5
-5

10

-5

315°E

0

2.5

360°E

360°E

90°S
0°

60°S

30°S

0°

30°N

60°N

90°S
0°
90°N

60°S

30°S

0°

30°N

60°N

90°N

0

d

250

0

45°E

0

-500

-400

45°E

-75

0

-500

90°E

-250

90°E

-250

-300

0

-2
5

-500

-50
0

135°E

180°E

180°E

Topography (meters)
-100
0
100

135°E

-200

225°E

225°E

200

0
0
25

0

-250

270°E

0

400

-250

270°E

300

-250

315°E

315°E

250

500

0

50
-2

2.5

17.5

360°E

0

360°E

Fig. 6. Titan's gravity and topography. Top panels (a–b) show full gravitational acceleration and topography measured with respect to a sphere of radius 2575 km. Panels (c–d) show only the non-hydrostatic gravity andtopography (i.e., departures from the hydrostatic expectation, where the hydrostatic component is computed via the Radau-Darwin relation) where we have assumed a moment of inertia factor of 0.341.

5

c

0

-2.5

-10

-12.5

6

0

90°S
0°

60°S

30°S

0°

30°N

60°N

90°S
0°
90°N

60°S

30°S

0°

30°N

-5

Gravitational Acceleration (mGal)
-4
-2
0
2
4

25

-6

-2.5

.5
17

60°N

-8

250

250
500

90°N

-2.
5

0

12.5

-5

129

2.5

5
2.
5 .5
7
10 2.5
1
15

2.5

5
-2.5

.
-7

5

5
2.5
0
.
-2 5

0

10
7.5

-2
.

15

2.5

.5
-12

-2
.5

0
-1

0

0
25

25

0

-5

0

0
-5
25

-10

D. Durante, et al.

Icarus 326 (2019) 123–132

Icarus 326 (2019) 123–132

D. Durante, et al.

90°N

-8

a

60°N

Gravitational Acceleration (mGal)
-4
-2
0
2
4

-6

6

8

10

-500

90°N

l=2

0.25

0.5

-0.5

-0.25

30°N

-400

60°N

30°N
0.5

Topography (meters)
-100
0
100

-200

75
100

5

12

-0.25

0
-25
-50
-75
-100
-125
-150

0.5

0

45°E

90°E

c

60°N

135°E

180°E

225°E

270°E

315°E

360°E

l=3

1

75

25

60°S

90°S
0°
90°N

100

50

0

60°S

500

0

75

30°S
0.25

-0.5

400

50

5
0.2

30°S

300

l=2

125

0°

200

-150
-125
-100
-75
-50
-25
0
25

10

0°

-300

-175

b

0

-10

-175

90°S
0°
90°N

45°E

90°E

d

60°N

135°E

180°E

225°E

270°E

-100

315°E

360°E

l=3
0
-5

50
0
10

-2

-50

30°N

0

-100
-50

0

-2
2

-2
-1
0

-5

0

-1

0

0

50

-10

0

30°S
2

30°S

10

0°

0

-3

150

-1

3
2
1

-2

0°

-150

-1

0

1

30°N

1

100

1

60°S

60°S
50

-100
0

-50

225°E

270°E

315°E
0
-0.5

e

-3
-2
-2 .5

315°E

l=5

2

45°E

90°E

135°E

180°E

225°E

0°
0

30°S

-2
5

-1
-15 25
0

-100

150

100

l=5

315°E

-75
-50

0
15

-250
10

50
75

1.5

45°E

5
-2

100

-125

90°S
360°E
0°

0

75

50

-50

60°S

-1
-0.5

270°E

360°E

75 50
25
-25 0

0
-1

25 50 75

1

315°E

100
125

-100
-5-75
0

-2
-1
-1 .5

1
1.
25

-2

-1
.5

0.5

270°E

5
17

0

-1.5

1

225°E

-125
50
-1

1.5

-2.5

0
.5
-0
-1
5

-1.

-1 0.5
0

-50
-100
-150

180°E

0
-10

1

2
0.5 1
0
-0.5-1
.5
-1 -2
5

0

250

-100

1.

-1

0
20

-75

-2

0

.5

1.5

-1

0.5

1

-2

h

135°E
-50
-75
-1-100-150 25

30°N

-1.5

1 .5
1
2

-0.5

90°E

.5
-0

5

5
0.5

1.5

60°N

45°E

0

-1

0
1 .5
1
2 .5

-1.5

1
0.5
0

1.5

-1.5

2.

-0
.

1.5

90°S
0°

90°S
360°E
0°
90°N

0
-5

0

100
12
5

270°E

30°S
-2

25
0

-2

225°E

2

-1

-2

200

0.5
1
1.5
2
2
3 .5
1

2

-20

0
-10
-150

180°E

-2 5
0

1.5

-1

-2

0

0°

60°S

60°S

100
50

-50

0

50

135°E

2.5

-1
.5

30°N

0

-2
00
-1
50

150

0.5

-0.5

0.5 11.5 2

-50

0

360°E

l=4

00
-1

-0.5

315°E

250

150
100

-2

90°E

270°E

50

1

g

60°N

30°S

225°E

-250

5

0

45°E

180°E

200
1500
100
5

0
-0.5-1
-1.5

1.5

-1
.

-1

90°S
0°
90°N

135°E

-150
00
-2

-2

-3

.5
-2

2
2.5
3

60°S

3
2

0.5

-10

0°

.5
2
1.51.5
0

30°S

1
1.5

150

-

0

.5
-11 .5- -0 0

50

1

0°

10
0

90°E

30°N

-1.5 -1-0.5

1.5

f

60°N

5
1. -1

50

0

-2

30°N

l=4

5

2.

1 1.5 2

3

60°N

-150
-100
-50
0

0
-20

0.5

45°E

200

180°E

-1
50

135°E

25
50
75
1
1200
15 5
0

90°E

50
100
150
20
0
25
0

45°E

0

90°S
360°E
0°
90°N
0.5

200

0

90°S
0°
90°N

-100

90°E

135°E

100 125

150

025

180°E

225°E

25
-50-25-75

270°E

0

315°E

360°E

Fig. 7. Non-hydrostatic gravity (left column) and topography (right column) for degrees 2–5. Whereas the signals beyond degree 2 (c–h) are independent of theassumed moment of inertia, the degree-2 signal (a–b) required subtracting the contributions from a hydrostatic body with moment of inertia factor 0.341.

Even in the limit of perfect isostatic equilibrium, however, the topography's contribution to the gravity signal should generally not be
zero because of the finite thickness of the ice shell, which attenuates thecompensating effect. Indeed, the relationship between the non-hydrostatic gravity and topography, typically quantified by the spectral admittance and correlation (e.g., Wieczorek, 2015; our Fig. 8), can often
be used to estimate the compensation depth (i.e., the ice shell thickness). In the present case, however, Titan's non-hydrostatic topography
does not contribute strongly to the non-hydrostatic gravity—the admittance and correlation spectra are compatible with zero (purple line
in Fig. 8). That is, while some non-hydrostatic gravity is evident, at least

beyond degree 2 (Fig. 7, left column), it is not spatially well correlatedwith the non-hydrostatic topography (Fig. 7, right column), precluding
direct quantification of the compensation depth via a spectral admittance analysis and therefore a meaningful estimate of the ice shell
thickness, at least on the basis of gravity and topography data alone.
This lack of spatial correlation could be an indication that the ice
shell is very thin, making isostatic compensation very effective in reducing the non-hydrostatic gravity. Alternatively, the topography could
be over-compensated. That is, the isostatic roots may be larger thanwhat would be expected in the limit of isostatic equilibrium. This could
be the case if the ice shell thickness variations arise from thinning/
130

Icarus 326 (2019) 123–132

D. Durante, et al.

a 20

b 1.5

10

1

0.5

-10

Correlation

Admittance (mGal/km)

0

-20
-30

0

-0.5

-40
Iess+ 2012/Zebker+ 2012 (l<=6)Iess+ 2012/Corlies+ 2017
This work/Zebker+ 2012 (l<=6)This work/Corlies+ 2017

-50
-60

2

3

4

-1

-1.5

5

2

3

4

5

Spherical Harmonic Degree

Spherical Harmonic Degree

Fig. 8. (a) Spectral admittance and (b) correlation between the non-hydrostatic parts of the gravitational acceleration and the topography, shown for several
combinations of shape and gravity models (see legend, inset; note that the Iess et al. (2012) nominal [SOL1a] solution was limited to degree 3), along with 1-σ error
bars. Admittance and correlation are computed from the normalized gravitational acceleration and topography coefficients following Wieczorek (2015). The hydrostatic component is computed assuming a moment of inertia factor of either 0.341 (solid lines) or 0.32 (dashed lines).

thickening taking place due to melting/freezing at the base of the ice
shell and if the ice shell is thick and stiff enough that elastic flexural
support is important—though a rigid ice shell would also reduce thetidal k2 Love number (see Hemingway et al., 2013). In fact, gravity/
topography correlation could even be negative if the ice shell is sufficiently stiff. An anti-correlation between Titan's degree-3 gravity andtopography signals was previously noted by Hemingway et al. (2013)
based on shape and gravity models available at the time (Iess et al.,2012; Zebker et al., 2012; blue lines in Fig. 8). When adopting the
updated shape and gravity models, however, the correlation at degree 3
becomes significantly weaker (Figs. 7 and 8). Likewise, the correlationsat degrees 4 and 5 are compatible with zero. Although strong compensation could help explain why the non-hydrostatic topography does
not contribute much to the non-hydrostatic gravity, the appearance of
uncorrelated non-hydrostatic gravity (Fig. 7, left column) suggests thepresence of additional internal heterogeneities that have little ornothing to do with the topographic anomalies (Fig. 7, right column),
perhaps pointing instead to anomalies within the rocky part of the interior (see also Palguta et al., 2006).
A possible exception to the lack of correlation occurs at degree 2,
where the non-hydrostatic gravity and topography could be spatiallycorrelated, depending on what is assumed for the moment of inertia
factor. The assumption of hydrostatic equilibrium amounts to choosing
the moment of inertia factor as the point on the hydrostatic line in
Fig. 4a that is nearest to the measured gravity or, equivalently, tominimizing the non-hydrostatic part of the degree-2 gravity (i.e.,
minimizing the magnitude of the signal shown in Fig. 7a). Under thisassumption, the non-hydrostatic degree-2 gravity is negligible and so
the degree-2 admittance is close to zero (Fig. 8a). Moments of inertia
larger or smaller than 0.341 are, however, possible if we relax the assumption of hydrostatic equilibrium. As mentioned earlier, the measured J2/C22 ratio being close to 10/3 is a necessary, but not sufficientcondition for hydrostatic equilibrium. Indeed, at least one study favors
a moment of inertia factor in the range 0.31–0.33 in order to satisfy
Titan's measured obliquity (Baland et al., 2014). In that scenario, thedeparture from hydrostatic equilibrium would involve more significant
non-hydrostatic gravity, which would be positively correlated with thenon-hydrostatic topography (i.e., both would involve excess power inthe l,m = 2,0 and l,m = 2,2 terms). This is illustrated by the dashedlines showing the positive admittance (Fig. 8a) and correlation (Fig. 8b)computed from the degree-2 non-hydrostatic gravity and topography

when the moment of inertia factor is assumed to be 0.32. This scenariodoes, however, require even the non-hydrostatic parts of J2 and C22 to
be in a ratio close to 10/3. Although such a ratio for the non-hydrostatic
gravity terms would be fortuitous, the possibility cannot be excluded.
5. Conclusion
The analysis of ten Cassini flybys of Titan dedicated to gravity science proved that the acquired data have been adequate for a valuable
determination of the gravity field of Saturn's largest moon. The gravity
coefficients up to degree 5 have been reported, increasing the resolutionof previous Cassini results, and providing insight into Titan's interior
structure. The new solution confirms the previous result of Titan's large
response to tidal forcing from Saturn: the confidence in the tidal Love
number k2 has been further increased, and the tidal time lag has been
determined to be close to zero at 1-σ, indicating a phase angle lag
smaller than 5.4°.
The new solution likewise confirms that, despite Titan's considerable excess flattening, its interior has relaxed to a state compatible withhydrostatic equilibrium with a moment of inertia factor close to 0.341,
though the presence of such significant non-hydrostatic topography is areminder that hydrostatic equilibrium is not guaranteed, precluding adefinitive determination of the moment of inertia. Comparison of the
non-hydrostatic parts of the gravitational potential with the topography
show that Titan's topography is highly compensated, as expected if theice shell is thin (permitting effective compensation) or stiff (permitting
overcompensation that remains effective even for a thicker ice shell).
The weak spatial correlation between the non-hydrostatic gravity andtopography beyond degree 2, however, may suggest the presence ofinternal mass anomalies at depth. A more strongly differentiated interior, with moment of inertia factor < 0.341, is permitted but would
imply less compensation at degree 2.
Although new observational constraints on Titan's gravitational
field must now await the next spacecraft mission to the Saturn system,
increasingly comprehensive models, coupled with constraints fromcomplementary observations (e.g., Titan's rotational state) and laboratory experiments, have the potential to yield further insights into Titan's
interior and evolutionary history and to resolve many of the remaining
puzzles. These include Titan's extreme excess flattening in spite of its
apparently hydrostatic gravitational field, the nature of Titan's topography at various length scales and how it is maintained in spite of the131Icarus 326 (2019) 123–132

D. Durante, et al.

basal ice being near its melting temperature (and therefore likely toundergo viscous relaxation), as well as the nature of those mass
anomalies that are uncorrelated with the observed topography.

and A. J. Verbiscer, Eds. Tucson, Arizona: University of Arizona Press, pp. 57–77.
Iess, L., Rappaport, N.J., Jacobson, R.A., Racioppa, P., Stevenson, D.J., Tortora, P.,
Armstrong, J.W., Asmar, S.W., 2010. Gravity field, shape, and moment of inertia ofTitan. Science 327, 1367–1369. https://doi.org/10.1126/science.1182583.
Iess, L., Jacobson, R.A., Ducci, M., Stevenson, D.J., Lunine, J.I., Armstrong, J.W., Asmar,
S.W., Racioppa, P., Rappaport, N.J., Tortora, P., 2012. The tides of Titan. Science
337, 457–459. https://doi.org/10.1126/science.1219631.
Iess, L., Di Benedetto, M., James, N., Mercolino, M., Simone, L., Tortora, P., 2014. Astra:interdisciplinary study on enhancement of the end-to-end accuracy of spacecraft
tracking techniques. Acta Astronautica 94, 699–707. https://doi.org/10.1016/j.
actaastro.2013.06.011.
Jacobson, R.A., Antreasian, P.G., Bordi, J.J., Criddle, K.E., Ionasescu, R., Jones, J.B.,
Mackenzie, R.A., Meek, M.C., Parcher, D., Pelletier, F.J., Owen, W.M., Roth, D.C.,
Roundhill, I.M., Stauch, J.R., 2006. The gravity field of the Saturnian system fromsatellite observations and spacecraft tracking data. Astron. J. 132, 2520–2526.
https://doi.org/10.1086/508812.
Kaula, W.M., 1963. Determination of the Earth's gravitational field. Rev. Geophys. 1 (4),507–551.
Liu, H.S., Chao, B.F., 1991. The Earth's equatorial principal axes and moments of inertia.
Geophys. J. Int. 106, 699–702. https://doi.org/10.1111/j.1365-246X.1991.
tb06341.x.
Matsuyama, I., 2011. Is Titan in Hydrostatic Equilibrium? In: EPSC-DPS2011-344-1.
Matsuyama, I., Beuthe, M., Hay, H.C.F.C., Nimmo, F., Kamata, S., 2018. Ocean tidal
heating in icy satellites with solid shells. Icarus 312, 208–230. https://doi.org/10.
1016/j.icarus.2018.04.013.
Meriggiola, R., Iess, L., Stiles, B.W., Lunine, J.I., Mitri, G., 2016. The rotational dynamics
of Titan from Cassini RADA images. Icarus 275, 183–192. https://doi.org/10.1016/j.
icarus.2016.01.019.
Murray, C.D., Dermott, S.F., 1999. Solar System dynamics. Cambridge University Press0521-57295-9 ISBN.
NASA's Navigation and Ancillary Information Facility (NAIF) n.d. Saturn's satellite
ephemerides can be found online: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/
Nimmo, F., Bills, B.G., 2010. Shell thickness variations and the long-wavelength topography of titan. Icarus 208, 896–904. https://doi.org/10.1016/j.icarus.2010.02.020.
Palguta, J., Anderson, J.D., Shubert, G., Moore, W.B., 2006. Mass anomalies on
Ganymede. Icarus 180, 428–441. https://doi.org/10.1016/j.icarus.2005.08.020.
Rappaport, N.J., Iess, L., Wahr, J., Lunine, J.I., Armstrong, J.W., Asmar, S.W., Tortora, P.,
Di Benedetto, M., Racioppa, P., 2008. Can Cassini detect a subsurface ocean in titan
from gravity measurements? Icarus 194, 711–720. https://doi.org/10.1016/j.icarus.
2007.11.024.
Seidelmann, P.K., Abalakin, V.K., Bursa, M., Davies, M.E., De Bergh, C., Lieske, J.H.,
Simon, J.L., Standish, E.M., Stooke, P., Thomas, C., 2002. Report of the IAU/IAG
Working Group on Cartographic Coordinates and Rotational Elements of the Planets
and Satellites: 2000. Celestial Mechanics Dynamical. Astronomy 82, 83–110. https://doi.org/10.1023/A:1013939327465.
Stiles, B.W., Kirk, R.L., Lorenz, R.D., Hensley, S., Lee, E., Ostro, S.J., Allison, M.D.,
Callahan, P.S., Gim, Y., Iess, L., Perci del marmot, P., Hamilton, G., Johnson, W.T.K.,
West, R.D., the Cassini RADAR team, 2008. Determining Titan's spin state from
Cassini RADAR images. Astron. J. 135, 1669–1680. https://doi.org/10.1088/00046256/135/5/1669.
Tricarico, P., 2014. Multi-layer hydrostatic equilibrium of planets and synchronous
moons: theory and application to Ceres and to solar system moons. Astrophys. J. 782,
99. https://doi.org/10.1088/0004-637X/782/2/992.
Wieczorek, M.A. (2015). “Gravity and topography of the terrestrial planets” in Treatise on
Geophysics, 2nd edition, Elsevier. doi:https://doi.org/10.1016/B978-0-444-53802-4.
00169-X
Zebker, H.A., Stiles, B., Hensley, S., Lorenz, R., Kirk, R.L., Lunine, J., 2009. Size and shapeof Saturn's moon Titan. Science 324 (5929), 921–923. https://doi.org/10.1126/
science.1168905.
Zebker, H.A., Iess, L., Wall, S.D., Lorenz, R.D., Lunine, J.I., Stiles, B.W., 2012. Titan's
figure fatter, flatter than its gravity field. AGU Fall Meeting 2012, P23F–01.

Acknowledgments
The work of D.D., P.R., and L.I. has been supported in part by the
Italian Space Agency under ASI-INAF contract 2017-10-H.O. D.J.H. was
supported in part by the Miller Institute for Basic Research in Science at
the University of California, Berkeley.
References
Archinal, B.A., A'Hearn, M.F., Bowell, E., Conrad, A., Consolmagno, G.J., Courtin, R.,
Fukushima, T., Hestroffer, D., Hilton, J.L., Krasinsky, G.A., Neumann, G., Oberst, J.,
Seidelmann, P.K., Stooke, P., Tholen, D.J., Thomas, P.C., Williams, I.P., 2011. Reportof the IAU working group on cartographic coordinates and rotational elements: 2009.
Celestial Mechanics Dynamical Astronomy 109 (2), 101–135. https://doi.org/10.
1007/s10569-010-9320-4.
Archinal, B.A., Acton, C.H., A'Hearn, M.F., Conrad, A., Consolmagno, G.J., Duxbury, T.,
Hestroffer, D., Hilton, J.L., Kirk, R.L., Klioner, S.A., McCarthy, D., Meech, K., Oberst,
J., Ping, J., Seidelmann, P.K., Tholen, D.J., Thomas, P.C., Williams, I.P., 2018. Reportof the IAU working group on cartographic coordinates and rotational elements: 2015.
Celestial Mechanics Dynamical Astronomy 130, 22. https://doi.org/10.1007/
s10569-017-9805-5.
Asmar, S.W., Armstrong, J.W., Iess, L., Tortora, P., 2005. Spacecraft Doppler tracking:noise budget and accuracy achievable in precision radio science investigations. Radio
Sci. 40. https://doi.org/10.1029/2004RS003101.
Baland, R., Tobie, G., Lefèvre, A., Van Hoolst, T., 2014. Titan's internal structure inferredfrom its gravity field, shape, and rotation state. Icarus 237, 29–41. https://doi.org/
10.1016/j.icarus.2014.04.007.
Barbaglio, F., Armstrong, J.W., Iess, L., 2012. Precise Doppler Measurements for
Navigation and Planetary Geodesy Using Low Gain Antennas: Test Results from
Cassini. (ISSFD 23rd).
Boone, D., J. Bellerose, and D. Roth (2017). Resolution of Orbit Determination Prediction
Instabilities at Titan During Cassini's Solstice Mission. ISSFD-2017–109.
Choukron, M., Sotin, C., 2012. Is Titan's shape caused by its meteorology and carbon
cycle? Geophys. Res. Lett. 39. https://doi.org/10.1029/2011GL050747.
Corlies, P., A. G. Hayes, S. P. D. Birch, R. Lorenz, B. W. Stiles, R. Kirk, V. Poggiali, H.
Zebker, and L. Iess (2017). Titan's topography and shape at the end of the Cassini
Mission. Geophys. Res. Lett., 44 (23), 11, 754–11, 761, doi:https://doi.org/10.1002/
2017GL075518
Darwin, G.H., 1899. The theory of the figure of the earth carried to the second order insmall quantities. Mon. Not. R. Astron. Soc. 60, 82–124. https://doi.org/10.1093/
mnras/60.2.82.
Ermakov, A.I., Park, R.S., Bills, B.G., 2018. Power laws of topography and gravity spectraof the solar system bodies. Journal of Geophysical Research: Planets 123, 2038–2064.
https://doi.org/10.1029/2018JE005562.
Evans, S., Taber, W., Drain, T., Smith, J., Wu, H., Guevara, M., Sunseri, R., Evans, J.,
2018. MONTE: the next generation of mission design and navigation software. CEAS
Space Journal 10, 79–86. https://doi.org/10.1007/s12567-017-0171-7.
Fuller, J., Luan, J., Quataert, E., 2016. Resonance locking as the source of rapid tidal
migration in the Jupiter and Saturn moon system. Mon. Not. R. Astron. Soc. 458,
3867–3879. https://doi.org/10.1093/mnras/stw609.
Hemingway, D.J., Nimmo, F., Zebker, H., Iess, L., 2013. A rigid and weathered ice shell on
Titan. Nature 500 (7464), 550–552. https://doi.org/10.1038/nature12400.
Hemingway, D. J., L. Iess, R. Tajeddine, and G. Tobie (2018). “The interior of Enceladus,”in Enceladus and the Icy Moons of Saturn, P. M. Schenk, R. N. Clark, C. J. A. Howett,132