Icarus 345 (2020) 113713

Contents lists available at ScienceDirect

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

The gravity field and interior structure of Dione
Marco Zannoni a, *, Douglas Hemingway b, c, Luis Gomez Casajus a, Paolo Tortora aaDipartimento di Ingegneria Industriale, Universit�
a di Bologna, Forlì, Italy
Department of Earth & Planetary Science, University of California Berkeley, Berkeley, CA, USAc
Department of Terrestrial Magnetism, Carnegie Institution for Science, Washington, DC, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Interiors
Orbit determination
Satellites, composition
Saturn, satellites

During its mission in the Saturn system, Cassini performed five close flybys of Dione. During three of them, radio
tracking data were collected during the closest approach, allowing estimation of the full degree-2 gravity field by
precise spacecraft orbit determination.
The gravity field of Dione is dominated by J2 and C22, for which our best estimates are J2 � 106 ¼ 1496 � 11
and C22 � 106 ¼ 364.8 � 1.8 (unnormalized coefficients, 1-σ uncertainty). Their ratio is J2/C22 ¼ 4.102 � 0.044,
showing a significative departure (about 17-σ) from the theoretical value of 10/3, predicted for a relaxed body inslow, synchronous rotation around a planet. Therefore, it is not possible to retrieve the moment of inertia directlyfrom the measured gravitational field.
The interior structure of Dione is investigated by a combined analysis of its gravity and topography, which
exhibits an even larger deviation from hydrostatic equilibrium, suggesting some degree of compensation. The
gravity of Dione is far from the expectation for an undifferentiated hydrostatic body, so we built a series of threelayer models, and considered both Airy and Pratt compensation mechanisms. The interpretation is non-unique,but Dione’s excess topography may suggest some degree of Airy-type isostasy, meaning that the outer ice shell isunderlain by a higher density, lower viscosity layer, such as a subsurface liquid water ocean. The data permit abroad range of possibilities, but the best fitting models tend towards large shell thicknesses and small ocean
thicknesses.

1. Introduction

era. From ground optical and spectroscopic observations, the surface
was determined to be composed of almost pure ice. The mass of Dione
was estimated by exploiting the mutual orbital perturbations with
Enceladus, using ground-based astrometric observations (e.g. Kozai,
1976; Harper and Taylor, 1993; Vienne and Duriez, 1995; Dourneau andBaratchart, 1999).
The knowledge about Dione advanced significantly with the Saturn
flybys of the twin Voyager probes, in 1980 and 1981. Voyager 1 pro­vided an almost global image coverage of the surface, revealing evi­
dence of past geological activity, including resurfacing events and,possibly, cryovolcanic features (Smith et al., 1981; Plescia, 1983). The
shape reconstructed from Voyager images, along with the mass given by
Earth-based astrometric measurements, provided a first measure of thebulk density of Dione, at about 1440 kg/m3 (Campbell and Anderson,1989). This relatively high value was compatible with a mixture of about
55% by mass water ice (density 1000 kg/m3) and 45% rock-metal(density 3000 kg/m3). Subsequently, Jacobson (1995) obtained an
improved estimation of the mass of Dione by measuring the orbital

With a mean radius of 561 km, Dione is the fourth-largest moon of
Saturn. It was discovered in 1684 by the Italian astronomer Giovanni
Domenico Cassini, during observations made at the Paris Observatory.
With a semi-major axis of approximately 6.26 Saturn radii (RS)(377,400 km), Dione is in a 1:2 mean-motion resonance with the smaller
moon Enceladus (252 km radius, 3.95 RS (237,900 km) semi-majoraxis). This resonance causes an orbital libration with a period of about
11 yr, and a circulation with a period of about 3.8 yr (Murray and
Dermott, 1999). The resonance maintains a non-zero orbital eccentricity
of both Enceladus and Dione, about 0.0047 and 0.0022, respectively.
Dione has two co-orbital moons, Helene (17.6 km radius, discovered in1980 from ground-based observations) and Polydeuces (1.3 km radius,
discovered in 2004 from Cassini images), orbiting around the Lagrangian
points L4 and L5, respectively.
Due to the small size and the large distance from the Earth, little was
known about Dione’s internal structure before the advent of the space* Corresponding author.
E-mail address: m.zannoni@unibo.it (M. Zannoni).

https://doi.org/10.1016/j.icarus.2020.113713
Received 19 August 2019; Received in revised form 19 February 2020; Accepted 20 February 2020
Available online 27 February 2020
0019-1035/© 2020 Elsevier Inc. All rights reserved.

M. Zannoni et al.

Icarus 345 (2020) 113713

perturbations caused on Helene using ground-based astrometric andVoyager images. In preparation for the Cassini tour, Jacobson (2004)
provided an updated estimate of the ephemerides and physical param­
eters of the Saturn system bodies, including the mass of Dione. An
extensive data set was used, comprising Earth-based astrometry, Pioneer
11 and Voyagers’ radiometric and optical data, and Cassini optical dataacquired before the arrival in the Saturn system.
During its 13-year tour of the Saturn system, Cassini performed fiveclose encounters of Dione, four of which were dedicated to the deter­
mination of its mass and gravity field, with the objective of constraining
its internal structure. During the first two flybys, referred to as D1(October 2005) and D2 (April 2010) according to the numbering scheme
used by the Cassini project, radiometric data were collected only before
and after (but not during) the closest approach (C/A), only allowing theestimation of the moon’s mass. Using Cassini data acquired during theSaturn tour up to June 2006, including the data acquired during D1,
Jacobson et al. (2006) provided an updated estimate of Saturn’s gravity
and pole orientation, and the masses of the satellites. In particular, themass estimation of Dione improved by 1 order of magnitude, with theinformation coming mainly from astrometric and spacecraft imaging ofHelene. The corresponding bulk density of Dione increased slightly to(1475.7 � 3.6) kg/m3. However, in the absence of measurements of thehigh-degree gravity harmonics, the internal structure could not be
inferred. The first flyby dedicated to the determination of Dione’s
gravity field was D3 (December 2011). The analysis of the Doppler dataacquired during the closest approach produced the first estimation of
Dione’s J2 and C22, suggesting that the moon is not compatible with thecondition of hydrostatic equilibrium (Iess et al., 2012a). However, given
the limited amount of data, the solution wasn’t fully reliable. For thisreason, during the extended mission, Cassini performed two other flybys
of Dione with tracking during the closest approach, D4 (June 2015) andD5 (August 2015), to better characterize the moon’s internal structure.
This paper presents the first estimation of Dione’s quadrupole gravity
field, obtained from the analysis of Doppler data acquired during all
Cassini flybys of the moon. The main characteristics of Dione’s Cassini
flybys are summarized in Table 1, while the corresponding groundtracks are displayed in Fig. 1. The flybys dedicated to gravity in­vestigations, D3 and D5, provide good spatial coverage for the retrievalof the quadrupole gravity field. D3 was nearly equatorial, with aninclination at the C/A of about 175� . In order to de-correlate the esti­
mation of J2 and C22, D5 was designed to be nearly polar, with aninclination at the C/A of about 96� . D3 also flew over the leading
hemisphere, whereas D5 was over the trailing hemisphere. Being sepa­
rated by only two months, D4 had an orbital geometry similar to D5.
However, while D3 was characterized by a very low altitude at the C/A,
<100 km, all the other flybys had a much higher altitude, of about 500
km, thus significantly reducing the sensitivity to the gravity field. The
noise on X-band Doppler measurements is mainly due to the solar
plasma and Earth troposphere (Asmar et al., 2005; Iess et al., 2012b).
The former is correlated with the Sun-Earth-Probe (SEP) angle, which
was larger than 50� during all encounters. The Doppler noise levelaround the C/A of the different flybys varies between a minimum of
0.021 mm/s, and a maximum of 0.036 mm/s.
This paper is organized as follows: Section 2 describes the data
analysis approach for the estimation of Dione’s gravity field, along withthe spacecraft dynamical model, and the data selection and calibration
procedure. Section 3 provides a geophysical interpretation of the results,
by means of a combined analysis of Dione’s estimated gravity andtopography. Finally, Section 4 summarizes our findings and conclusions.
2. Gravity analysis
2.1. Introduction
The gravity field of Dione was estimated by precisely reconstructing
the trajectory of Cassini during the five close flybys of the moon. The
estimation of the gravity field of Dione was based on the same procedureand techniques adopted in the previous gravity analyses of Saturn’s
moons performed by the Cassini Radio Science Team (Iess et al., 2012c;
Iess et al., 2014; Tortora et al., 2016; Durante et al., 2019).
The main difference from past gravity analyses relates to the updateof the ephemerides of Dione. Doppler data acquired around the peri­
center of a flyby are very sensitive to the relative position of the Cassini
spacecraft with respect to the moon. Outside the sphere of influence ofDione, which has a radius of about 2000 km, the data are sensitive
mainly to the relative position of the spacecraft with respect to Saturn.
During a close encounter, Cassini stays inside the sphere of influence ofthe moon for about 10 min. As a result, the orbit of Dione must be known
at a level currently not met by the JPL satellite ephemerides, and so it
must be estimated and updated as a part of the orbit determination
procedure.
In previous work by Tortora et al. (2016) and others, the orbit of themoon under study was numerically integrated for the entire time spancovered by the data, from an epoch prior to the first flyby, to after thelast flyby. This approach ensures that the satellite trajectory is dynam­
ically coherent.
However, as in the work of Durante et al. (2019) for Titan, in thiswork it was not possible to obtain a satisfactory fit of the data by esti­
mating a single, coherent, orbit of Dione. This may be an indication of an
incomplete dynamical model of the Saturn system, given the longtimespan covered by Cassini data (D1 and D5 are separated by about 10
years), the poor sampling of the orbit with time (5 encounters), and thehigh level of accuracy of the data. Possible areas of improvement of the
dynamical model are: the proper modelling of a time-variable or
longitudinally-dependent component of the gravity of Saturn (Iess et al.,2019); the ephemerides of Enceladus, whose gravitational perturbations
on Dione are relevant because of the orbital resonance; the evolution of
Saturn’s pole; the tidal interaction between Saturn and its moons. In
particular, regarding the latter, recent measurements suggest that the
tidal dissipation of Saturn is higher than predicted by standard tidal
theories, and that it is not constant between the different satellites, aspredicted by the resonance locking tidal theory (Lainey et al., 2012,2017, 2020; Fuller et al., 2016). However, the accurate modelling of the
motion of Dione inside the Saturn system was beyond the scope of thispaper. Hence, the same approach followed by Durante et al. (2019) was
adopted, estimating an updated orbit of the moon for each encounter.
This over-parametrization causes an increase of the uncertainties of thequadrupole gravity coefficients up to 40%, but it ensures an unbiased
estimate of Dione’s gravity field.

Table 1
Summary of Cassini flybys of Dione. For each encounter, the table reports the date of Cassini’s closest approach (C/A), the minimum altitude reached, the orbital
inclination, the Sun-Earth-Probe (SEP) angle, the number of Doppler points (60 s integration time), and the root mean square (RMS) of Doppler residuals at 60 s.
Flyby

Date of C/A

Altitude (km)Rel. Vel. (km/s)

Inclination (deg)SEP (deg)

Doppler points

RMS at 60s (mm/s)

D1
D2
D3
D4
D5

11 Oct. 2005
07 Apr. 201012 Dec. 201116 Jun. 201517 Aug. 2015

495
503
98
511
476

9.1
8.3
8.7
7.3
6.4

120
0.5
175
80
96

69
163
53
154
93

829
410
1309
989
1144

0.035
0.027
0.035
0.038
0.035

2

M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 1. Cassini ground track on Dione during D1, D2, D3, D4 and D5, considering a time interval of �15 min around the C/A (indicated by a circle). The ticks are
separated by 60 s. The ground tracks are represented over a visible map of Dione produced by Paul Schenk (Lunar and Planetary Institute) from Cassini ISS data
(NASA, JPL).

2.2. Dynamical model

modeled using the even zonal harmonics J2-J10 and J3, as provided by
the reference satellite ephemerides. For consistency, the corresponding
rotational model of Saturn was adopted.
In estimating Dione’s gravity coefficients (Eq. (1)), the minimum
field capable of fitting the Doppler data to the noise level was a fulldegree 2. Higher degree and order fields were also estimated, to assess
the stability of the solution. The time-variable gravity field of Dione
caused by eccentricity tides was neglected, due to the very low orbital
eccentricity and the limited data coverage. Considering only the en­counters with coverage at C/A, for a k2 tidal Love number of 0.5, theexpected variation of J2 and C22 due to eccentricity tides is about 0.1%
and 0.3%, respectively. These values are below the sensitivity of themeasurements by almost a factor 2. Nevertheless, the k2 tidal Love
number of Dione was also estimated, as a stability test.
Regarding the rotational model of Dione, we adopted a dynamicallydefined, perfectly synchronous frame, which points always to the empty
focus of the orbit (Murray and Dermott, 1999). In addition, to assess thestability of the solution, we alternatively assumed the rotational models
suggested by IAU (see Section 2.5).
The dynamical model of Cassini included also the main nongravitational accelerations: the Solar Radiation Pressure (SRP), andthe thermal recoil pressure due to the anisotropic thermal emission,
mainly caused by the three on-board Radioisotope Thermoelectric
Generators (RTG). For both the accelerations, the models adopted by the
Cassini navigation team were implemented (Di Benedetto, 2011). The
accelerations due to the albedo and infrared thermal emission of Dione
were neglected, being at least 1 order of magnitude smaller than the SRP(Di Benedetto, 2011).
Finally, during D4, two attitude maneuvers, executed with thrusters,
were performed about 12 h and 4 h before closest approach. These
maneuvers were modeled as impulsive changes in the spacecraft ve­
locity vector, starting from the values reconstructed by the Cassini
navigation team.

The adopted dynamical model included all the relevant accelerations
acting on Cassini and on Dione, mainly the relativistic gravitational ac­
celeration due to Saturn, its main satellites, the Sun, the other planets of
the Solar System, the Moon, and Pluto. The masses and the states of theplanets, the Moon, and Pluto were retrieved from the latest planetary
ephemerides produced by JPL (DE438). The masses and the states of theSaturn satellites were retrieved from the latest satellite ephemerides
produced by the JPL (SAT389). Both planetary and satellite ephemer­
ides can be retrieved from ftp://ssd.jpl.nasa.gov/pub/eph/. Different
versions of satellite ephemerides were also adopted, to test the stabilityof the solution.
Dione’s gravity field was modeled as a linear combination of
spherical harmonic functions:Uðr; θ; ϕÞ ¼

�llmax X
l �
GM X
Rref
ðClm cosmϕ þ Slm sinmϕÞPlm ðcosθÞr l¼0 m¼0 r(1)

where Clm and Slm are the unnormalized degree-l and order-m dimen­
sionless gravitational potential coefficients for a reference radius Rref ¼
560 km, Plm are the unnormalized associated Legendre functions (e.g.,
Wieczorek, 2015), lmax is the maximum degree considered in the
expansion, GM is the product of the gravitational constant and Dione’s
mass (see Table 2), and the potential is evaluated here at radial positionr, co-latitude θ, and longitude ϕ. In our discussion, we often represent
the zonal terms in the expansion as Jl ¼
Cl0. In addition, the setup
included the accelerations due to the extended gravity field of Saturn,Table 2Estimated values and 1-σ formal uncertainties of Dione’s quadrupole gravity
unnormalized coefficients. The adopted a priori values and uncertainties are also
shown. The reference radius for the spherical harmonics is 560 km. The degree-1
terms are assumed zero, so that the origin of the Dione body-fixed frame is themoon’s center of mass. The estimated value of the ratio J2/C22 and the corre­
lation coefficient between J2 and C22 are also shown.
Unit
GM
J2
C21
S21
C22
S22
J2/C22
corr J2-C22

3

km /s
�106
�106
�106
�106
�106

A priori
2

73.116
1430
0
0
365
0
3.9

�
�
�
�
�
�
�
0.0

2.3. Data selection and calibration

Solution
0.025
500
60
20
130
20
1.9

73.1118
1496
0.6
4
364.8
14.2
4.102

�
�
�
�
�
�
�
0.50

The observable used in the estimation procedure was the range-rate,
derived from the Doppler shift of a highly stable microwave carriertransmitted between Cassini and the ground antennas of NASA’s DeepSpace Network (DSN). The range observables, derived from the roundtrip light time of a modulated code, were not used, because the infor­
mation content provided is mainly related to planetary and satellite
ephemerides, while the sensitivity to gravity field is limited.
The count time of Doppler data was chosen as a trade-off between the

0.0025
11
6.8
20
1.8
1.9
0.044

3

M. Zannoni et al.

Icarus 345 (2020) 113713

sensitivity to gravity spherical harmonics and numerical considerations.
In fact, at the pericenter rp of a close flyby, the spatial scale Δl of aspherical harmonic of degree l is (Milani and Gronchi, 2010):πΔl ¼ rplthe solution. No hydrostatic equilibrium constraint between J2 and C22
was imposed. For C21, S21, and S22 we used a different strategy. From
MacCullagh’s theorem, these gravity coefficients are related to amisalignment between the adopted Dione-fixed frame and its principal
axes of inertia. Because the data are not sufficiently sensitive to therotational state of Dione, the a priori uncertainties of C21, S21, and S22
were set to a value corresponding to a rotation of about 1� . Larger values
were also used, up to a misalignment of 20� , to assess the stability of thesolution.
In addition, a set of “local” parameters, affecting only a single
encounter, was estimated. For each encounter, they include the initialstate of Cassini and Dione, a constant correction to the SRP acceleration,
constant Doppler bias for the three-way passes, and the impulsive ΔV
due to the maneuvers executed during D4. The a priori uncertainties for
Cassini’s position and velocity were 10 km and 0.1 m/s, respectively.

(2)

The time interval corresponding to the spatial scale at pericenter isobtained by dividing by the relative velocity:
Δt ¼

Δl π rp
¼
vp
l vp

(3)

The obtained time interval represents the theoretical maximum
sampling time to correctly reconstruct the gravity field of degree l.
Among all the Cassini flybys of Dione, the smallest time interval asso­
ciated with the degree 2 field is 120 s, obtained during D3. This value
represents a lower bound, because the orbit is hyperbolic. Therefore,
Doppler data were integrated over 60 s, sufficiently smaller than theminimum time interval to be sensitive to the low degree gravity field andsufficiently large to avoid numerical noise issues (Zannoni and Tortora,
2013).
During the encounters, Doppler data at X (8.4 GHz) and Ka band
(32.5 GHz) were acquired by the antennas of NASA’s DSN, phasecoherent to a common X-band (7.2 GHz) uplink. The analysis also used
tracking data from standard navigation passes, covering about 5 daysaround the closest approach. This marginally improves the uncertainties
in the gravity field estimation because of the improved reconstruction of
the orbits of Cassini and Dione during the encounter. When available, X/
Ka Doppler data were preferred over the standard X/X measurements,
because they are less affected by the dispersive sources of noise, likesolar plasma and Earth’s ionosphere. When two-way Doppler data were
not available, three-way data were also used. However, a bias on threeway data, constant per tracking pass, was estimated, accounting for a
possible offset between the clocks of different DSN complexes. The
additional path delay due to the Earth’s troposphere was corrected using
the standard GPS-based calibrations or, when available, the advanced
calibrations based on Water Vapor Radiometers (Bar-Sever et al., 2007).
The data were analyzed using JPL’s orbit determination program
MONTE (Evans et al., 2018), currently used for the operations of several
NASA deep space missions and for past radio science data analysis (e.g.
Iess et al., 2018; Iess et al., 2019). The mathematical formulation of
MONTE is described in detail in (Moyer, 1971) and (Moyer, 2000). Data
were weighted using the observed RMS of the residuals, constant foreach pass. Data acquired below 15 degrees of elevation, as viewed fromthe ground station, were discarded because of possible residual cali­
bration errors of the Earth’s troposphere and ionosphere, including dataaround the closest approach of D4 (between C/A - 1 h and C/A þ 10
min).

2.5. Results
The estimated gravity field coefficients of Dione are reported inTable 2, while Fig. 2 shows the estimated values of J2 and C22 in the C22J2 plane. Dione’s quadrupole is dominated by J2 and C22, as expected bya satellite in synchronous rotation around its planet. However, the ratio
J2/C22 is 4.102 � 0.044, about 17-σ away from the ideal hydrostatic
value of 10/3. Therefore, Dione’s gravity field is significantly nonhydrostatic, meaning that the moment of inertia cannot be inferred
directly from either J2 or C22 using the Radau-Darwin approx­
imation—the interpretation requires a more sophisticated approach (seeSection 3). The correlation coefficient between J2 and C22 is 0.5:
although this value is not ideal (due to the limited coverage of the fly­bys), the important point is that the estimation of both parameters was
not artificially constrained.
The estimated value of C21 and S21 are null within 1-σ . These values
correspond to misalignments between the assumed spin axis and the
maximum inertia axis of 0.02� � 0.17� and 0.3� � 1.5� around the yand x axes, respectively. However, the uncertainty of S21 does not
significantly improve with respect to the adopted a priori value, con­firming a relatively low sensitivity of Cassini’s tracking data to therotational state of Dione, in particular to rotations around its assumed
prime meridian. S22 is larger, being about 7-σ away from zero. This
corresponds to a misalignment between the principal axis of inertia andthe prime meridian used in the analysis (pointing to the empty focus of
the orbit of Dione around Saturn) of about 1.12� � 0.15� .
The stability of the solution has been assessed by perturbing the
adopted dynamical model, the data selection, and the estimation setup(such as the a priori covariances). Given the strategy to update thesatellite ephemerides only locally, different sets of a priori ephemerides
were also adopted. Since the nominal solution was obtained using all
available data, different combinations of encounters were also tested,such as using only one encounter or removing one encounter from thedataset. Additionally, we considered different rotational models of
Dione, in particular the ones suggested by IAU (Seidelmann et al., 2002;
Archinal et al., 2018). In all cases, the estimated values were compatible
with the reference solution within 1-σ, and the residuals were of very
similar quality.
While a quadrupole gravity field is fully sufficient to fit the data to
the noise level (Fig. 3), the neglected higher degree components of thepotential may introduce a bias in the estimation of J2 and C22. To test therobustness of the reference solution, we estimated also a gravity field ofdegree and order 4. However, given the number and the geometry of the
Cassini flybys, an unconstrained estimation of the higher degrees is notpossible. Thus, the a priori uncertainties on the normalized coefficientsof degree l were set using the Kaula rule K/l2 (Kaula, 1963). Thisempirical law can successfully describe the gravity power spectrum of
the rocky planets, the Moon, and Vesta (Ermakov et al., 2018). A good
agreement was also found for Titan, even if the gravity field is available
only up to degree 5 (Durante et al., 2019). As of today, there are no2.4. Estimation
The data analysis was carried out using a multi-arc approach, inwhich radiometric data obtained during the different encounters areanalyzed together to produce a single solution of a set of “global” pa­rameters, which do not vary among the arcs (Milani and Gronchi, 2010).
The multi-arc approach has been successfully applied to the analysis ofradio science data of several deep space missions (Iess et al., 2012c; Iesset al., 2014; Modenini and Tortora, 2014; Tortora et al., 2016; Zannoni
et al., 2018; Durante et al., 2019; Serra et al., 2019, Gomez Casajus et al.,2020). The parameters were estimated using a weighted least-squares
batch filter, which determines corrections to an a-priori dynamical
model to minimize the difference between the real and the simulated
measurements (Bierman, 1977).
The set of global parameters includes the gravitational parameter(GM) of Dione, its full degree-2 gravity field, and Cassini’s RTG accel­
eration at a reference epoch. The a priori uncertainties of Cassini’s RTG
acceleration, Dione’s GM, J2, and C22 were chosen to avoid constraining
4M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 2. Amplitude of latitudinal verses longitudinal variations in the gravitational field (J2 vs C22). Dashed line represents the expectation for a body in perfecthydrostatic equilibrium assuming various possible moments of inertia ranging from 0.3 to 0.4 (0.4 corresponds to a uniform body). The theoretical hydrostatic values
are computed following the approach of Tricarico (2014).
Fig. 3. Two-way range-rate residuals around theclosest approach of D3, D4, and D5 encounters (ver­tical line), in mm/s. The quality of the residuals isgood, which is an indication of a good orbital fit: noevident signatures are visible, the mean is approxi­
mately zero. The RMS is about 27 μm/s (D3), 21 μm/s
(D4), and 36 μm/s (D5). As a comparison, the twoway range-rate signal due to the degree 2 field isabout 20 mm/s (D3), 12 mm/s (D4), and 14 mm/s
(D5). The closest approach of D4 is not covered
because the elevation was lower than 15� , as viewed
from the ground station.

5

M. Zannoni et al.

Icarus 345 (2020) 113713

geophysical arguments or empirical evidence to justify its applicability
to the mid-sized icy moons of Saturn. However, even increasing thecoefficient K up to a very large value of 10 3, the quadrupole remains
compatible within 1-σ with the reference solution, confirming its
stability.
The tides raised by Saturn produce a time variable component of thegravity field, which can be modeled using the k2 tidal Love number.
Even if the orbital eccentricity of Dione is small and the coverage pro­vided by the Cassini encounters is limited, to test the robustness of thesolution we tried also to estimate both the real and imaginary compo­
nents of k2. The estimated gravity coefficients remained compatible withthe reference solution within 1-σ , while the estimated component of theLove number are Re(k2) ¼ 0.01 � 0.58 and Im(k2) ¼ 0.04 � 0.70,
statistically equivalent to zero, confirming that the tidal response of
Dione at the timescale of its orbital period is not observable using Cas­sini’s tracking data.
Constraining the ratio J2/C22 to the ideal hydrostatic equilibrium
value of 10/3, the residuals show a large signature at the closest
approach of D3 and D5 (Fig. 4), confirming that Cassini’s data are notcompatible with Dione being in hydrostatic equilibrium.

3. Interpretation
3.1. Basic observations
Because their interiors are weak on long timescales, bodies as large
as Dione are expected to have relaxed to hydrostatic equilibrium. Thatis, they should exhibit near spherical symmetry, with some small de­partures from symmetry arising due to centrifugal and tidal forces. The
magnitude of these asymmetries (captured by the J2 and C22 gravity
coefficients) is a function of the rotation rate, the mass and proximity ofthe parent body, and the body’s internal radial density structure andtherefore its moment of inertia. In general, one can use the RadauDarwin relation (Darwin, 1899; Murray and Dermott, 1999) tocompute the hypothetical hydrostatic J2 and C22 gravity coefficients fora range of possible moments of inertia, which can then be compared tothe measured J2 and C22 to assess the degree to which the body’s rela­
tively stiff exterior supports a departure from the hydrostatic expecta­tion. In our analysis, we employ the slightly more accurate approach of
Tricarico (2014) in computing the expected hydrostatic figure (dashedblack line in Fig. 2, with slope ~3.307).
The statistically significant departure from hydrostatic equilibrium(J2/C22 ¼ 4.102 � 0.044) makes it impossible to determine the precise
moment of inertia (and hence radial density structure) directly from themeasured gravitational field. It is, however, clear that Dione is far fromthe expectation for an undifferentiated hydrostatic body, for which theFig. 4. Two-way range-rate residuals around the closest approach of D3, D4, and D5 encounters (vertical line), in mm/s, constraining the ratio J2/C22 to the idealhydrostatic equilibrium value of 10/3. Strong signatures are visible around the closest approach of D3 and D5, not present in the unconstrained solution (Fig. 3).
6

M. Zannoni et al.

Icarus 345 (2020) 113713

gravitational potential coefficients would be J2 ¼ 2127 � 10 6 and C22¼ 649 � 10 6 (the upper right terminus of the dashed black line in Fig. 2,
where the moment of inertia factor is 0.4). The measured gravitational
potential is more consistent with some degree of differentiation and a
moment of inertia factor of roughly 0.32–0.34 (Fig. 2; more precise
values are obtained below). For a simple two-layer interior model,
consisting of a rocky core surrounded by an envelope of water ice (withdensity 930 kg/m3), this moment of inertia range corresponds to coreradii and densities of 380–430 km and 2100–2700 kg/m3, respectively(Fig. 5). Hence 56–66% of Dione’s mass (31–46% of its volume) is in the
rocky core. Assuming (non-hydrated) rock and ice densities of 3500 and
930 kg/m3, respectively, yields a rock:ice ratio very close to 1:1 by mass(or about 1:4 by volume).
To take the analysis further, we can combine the above gravitational
field with a model of the shape. The radius and shape of Dione have been
determined via analyses of limb profiles (e.g., Thomas et al., 2007;
Thomas, 2010; Nimmo et al., 2011). Incorporating the latest limb profile
observations (P. Thomas, personal communication), and repeating theanalysis of Nimmo et al. (2011), we obtain an updated model for the
degree-2 shape (Table 3, F. Nimmo, personal communication).
Fig. 6 shows our shape model along with previously publishedmodels (Thomas, 2010; Nimmo et al., 2011). Also shown in Fig. 6 are thegeoid (small purple ellipse) and the expectations for a perfectly hydro­
static Dione (dashed line), assuming various moments of inertia (i.e.,corresponding to a range of possible radial density structures). The geoid
coefficients are approximated byNðθ; ϕÞ ¼ Rref

∞ XlXðClm cosmϕ þ Slm sinmϕÞPlm ðcosθÞl¼0 m¼0

þω2 R2ref4g

Table 3The shape of Dione, based on analysis of Cassini ISS limb profiles, and repre­
sented here with unnormalized spherical harmonic expansion coefficients up todegree-2 (F. Nimmo, personal communication). In order to distinguish from thegravitational potential coefficients (Clm), we have used the abbreviation Hlm torepresent the degree-l and order-m shape coefficients (uncertainties are 1-σ).
Unit
R
H20
H22
H20/H22

Value
m
m
m

561.4
1834
374
4.9

�
�
�
�

0.4
134
13
0.4

on the leading and trailing faces when measured with respect to theobserved geoid (Fig. 7). This is unusual and the reason for it is not
obvious. However, it is worth emphasizing that the most prominent
topographic feature is the elevated region found near 60� E and justnorth of the equator (Fig. 7b), and coinciding with parts of the brightchasmata found on the trailing hemisphere. A less pronounced, but
broad topographic high also exists on the leading hemisphere. Althoughlimb profile coverage is incomplete, it is sufficient to reveal these fea­tures clearly (Fig. 8).
3.2. Isostatic compensation and interior models
The existence of considerable non-hydrostatic topography indicatesthat the exterior of Dione has been cold and rigid enough to support theassociated stresses since the formation of that topography, such that thefigure has not completely relaxed to hydrostatic equilibrium. The factthat the corresponding non-hydrostatic gravity is smaller by compari­son, however, is an indication that this non-hydrostatic topography is atleast partly compensated (e.g., isostatically). That is, the topography’s
contribution to the gravitational field is partly offset by internal mass
anomalies, likely resulting from lateral density variations and/or relief
along internal density boundaries. The relationship between the nonhydrostatic gravity and the non-hydrostatic topography is a functionof the degree and depth of this compensation, and therefore tells us
about the shallower internal structure of Dione. The challenge, however,is to isolate these non-hydrostatic signals from the total observed shapeand gravitational field, which are strongly affected by rotational andtidal deformation.
Following an approach developed for Enceladus (Hemingway et al.,2013; supplement of Iess et al., 2014), we consider Dione to be a mostlyhydrostatic body (whose shape and gravitational field are dominated byrotational and tidal deformation), superimposed with some nonhydrostatic topography (of unspecified origin). That is,5ω2 R2refP20 ðcosθÞ
6g

P22 ðcosθÞcos2ϕ(4)where ω is Dione’s rate of spin and g ¼ 0.232 m/s2 is its surface gravity.
In spite of the considerable uncertainties in the shape model, it isclear from Fig. 6 that, compared with the gravitational field (and thecorresponding geoid), the measured shape exhibits a greater departurefrom the hydrostatic expectation, with the ratio between the maindegree-2 coefficients being H20/H22 ¼ 4.9 � 0.4 (recall that the corre­
sponding ratio for the measured gravitational field is J2/C22 ¼ 4.102 �
0.044). Because the H22 component of the shape is smaller than thecorresponding term for the geoid, the figure exhibits topographic highs

hyd
obs
nh
H20
¼ H20
þ H20
hyd
obs
nh
H22
¼ H22
þ H22

(5)

For a moment of inertia factor of 0.33 (shown later to be represen­tative of the best fitting interior models), it can be shown (see methods
by Hemingway and Mittal, 2019) that the degree-2 figure and gravity
break down into their hydrostatic and non-hydrostatic parts as shown in
Table 4 and that the non-hydrostatic gravity expected if the topography
were uncompensated is Cnh
188 � 106 and C22nh � 64 � 106. The
20 �
non-hydrostatic values in Table 4 being 63–65% of the uncompensatedvalues is one measure of the degree of compensation.
We model Dione as consisting of a rocky (though not necessarily
purely rock) core surrounded by an H2O envelope. To model thecompensation of the surface topography, we consider the end member
cases of Airy and Pratt isostatic equilibrium. For Airy compensation, in
which the topography is supported by lateral variations in the thicknessof the outer ice shell, this ice shell must be underlain by a higher density,
lower viscosity material—most naturally a liquid water ocean. Hence,for our Airy models, we partition the H2O envelope into liquid and solid
phases. For Pratt compensation, the topography is compensated by

Fig. 5. Dione’s core radius and density for a range of likely moments of inertia(0.32–0.34), assuming a simple two-layer model consisting of a rocky coresurrounded by a water ice envelope with density 930 kg/m3. The total radiusand bulk density are 564.1 km and 1478 kg/m3, respectively.
7

M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 6. Amplitude of latitudinal verses longitudinalvariations in the shape (H20 vs H22) for various shapemodels, including our new shape model (gold), incomparison with our computed geoid (purple).
Dashed line represents the expectation for a body inperfect hydrostatic equilibrium assuming various
possible moments of inertia ranging from 0.3 to 0.4(0.4 corresponds to a uniform interior). The theoret­
ical hydrostatic values are computed following theapproach of Tricarico (2014). The negative value of
H20 is shown in order to maintain correspondencewith plots of J2 vs C22. (For interpretation of thereferences to color in this figure legend, the reader isreferred to the web version of this article.)

Fig. 7. Dione’s shape (Nimmo et al., 2011) with respect to the observed geoid (an equipotential surface), expanded to spherical harmonic degree 8 (a-b) andspherical harmonic degree 2 (c-d). Panels (b) and (d) show the topography as contours over a visible map of Dione produced by Paul Schenk (Lunar and Planetary
Institute) from Cassini ISS data (NASA, JPL). Compared with the geoid, Dione’s shape exhibits excess polar flattening and topographic highs on the leading andtrailing faces (see also Fig. 8). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)�
obs
conforming to the observed values Hobs
20 ; H22 . The models are param­
eterized according to the mean thicknesses and densities of the two outerlayers, yielding a four-dimensional parameter space. For each point in
the parameter space, the mean radius and density of the innermost layer
is constrained by the known total radius (561.4 km) and bulk density(1478 kg/m3). For each interior model, we use the numerical approach��hydof Tricarico (2014) to compute the hydrostatic terms Hhyd
20 ; H22 . The�nhnhremaining non-hydrostatic topography H20is then assumed to be; H22lateral density variations that persist through to some compensation
depth. For these models, we partition the H2O envelope into an upper
layer, in which there are lateral density variations, and a lower layer
with some uniform density that is slightly greater than that of the var­iable density upper layer; both layers are assumed to be in the solid
phase.
Following the approach of Hemingway and Mittal (2019), weconstruct a series of three-layer models with the exterior shape8M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 8. Limb profile elevation data (P. Thomas, personal communication) shown relative to the geoid, with warm and cool colors indicating high and low elevations,
respectively (compare with our Fig. 7 and with Figure 2d by Nimmo et al., 2011). The background image is a visible map of Dione produced by Paul Schenk (Lunarand Planetary Institute) from Cassini ISS data (NASA, JPL) and is intended to provide geographic context. Although limb profile coverage is sparse, topographic highs
noticeably coincide with the chasmata on the trailing hemisphere (around 60–80� E) and, to a lesser extent, the leading hemisphere (around 270� E), whereas there isa topographic low near the anti-meridian (around 180� ). (For interpretation of the references to color in this figure legend, the reader is referred to the online versionof this chapter.)shell thicknesses can be traded off against larger ocean thicknesses(Fig. 9c). Similar trade-offs exist between core radius and density, but
the preferred range is approximately 400 � 25 km and 2400 � 200 kg/
m3 (Fig. 9d). When a wider range of possible ice and ocean densities is
considered, the range of best fitting parameters widens only slightly(Fig. 9e–f); hence, the ice and ocean densities are not well constrained.
An analysis by Beuthe et al. (2016), based on our preliminary results(Hemingway et al., 2016), delivered a similar range of best fittinginterior models (compare their Fig. 2b with our Fig. 9c), in spite of theiruse of somewhat different methods for computing the equilibrium fig­
ures and the isostatically compensating basal topography.
In the Pratt model, not considered in previous work, the topography
is compensated instead by lateral density variations in the outermost
layer. Such density variations could be the result of differences inporosity within the outermost layer, for example. Assuming nominal
mean densities of 925 kg/m3 for the outermost layer, and 1030 kg/m3
for the intermediate layer (which we call the mantle), the best fittinginterior models correspond to a compensation depth of roughly 160 km,
though again there is a range of possibilities as compensation depth canbe traded to some degree against the thickness of the underlying mantle.
The range of likely core radii and densities is very broad. The range ofpossibilities widens when we consider a broader range of mean layer
densities. In all cases, however, the best fitting results correspond tocompensation depths >60 km. Confining the density anomalies to ashallower layer would require more pronounced density anomalies,
leading to values of J2 and C22 that would exceed the observed value (i.
e., as we move below the curves shown in Fig. 10a).

Table 4
Decomposition of principal degree-2 terms into their hydrostatic and nonhydrostatic components.
H20
H22
C20
C22

Unit

Observed

Hydrostatic

m
m
x106
x106

1834
374
1464
365

1555
469
1342
406

Non-hydrostatic
279
95
122
40

compensated isostatically. We compute the compensating basal topog­
raphy (when assuming Airy compensation) or density variations (whenassuming Pratt compensation) using the equal pressures isostasy
approach of Hemingway and Matsuyama (2017). The shape of theinnermost layer (the core) is assumed to conform to the hydrostatic
expectation. This is a reasonable assumption if Dione’s core is suffi­
ciently unconsolidated to be weak on long timescales. Departures from ahydrostatic figure for the core are possible if the material is sufficiently
strong (e.g., McKinnon, 2013; Tajeddine et al., 2014) but, as in previous
work (Hemingway and Mittal, 2019), we regard a hydrostatic core as thesimplest assumption since there is no reason to prefer non-hydrostatic
core topography of any particular orientation. Allowing for a wider
range of core shapes would widen the uncertainties in our internal
structure models but would not bias our results in any particular di­rection. Finally, we compute the resulting gravitational field (Hemi­
ngway and Mittal, 2019, Eq. (8)), taking into account the finite
amplitude effects (Wieczorek and Phillips, 1998), and compare theresult with the observed gravitational field (Table 2). We use the misfit
between the model and measured gravitational fields to construct a
probability density function across the parameter space (see section 2.6
of Hemingway and Mittal, 2019), indicating which parameter values aremost likely, given the observations. We carry out this exercise assuming
either Airy (Fig. 9) or Pratt (Fig. 10) isostasy.
In the Airy isostasy model, the topography is supported by lateral
thickness variations in an ice shell floating on a liquid water ocean.
Assuming nominal ice shell and ocean densities of 925 kg/m3 and 1020
kg/m3, respectively, the best fitting interior models tend towards large
shell thicknesses (as much as 140 km) and small ocean thicknesses.
However, the data permit a considerable range of possibilities as smaller

3.3. Discussion
Because we do not evaluate the probabilities in absolute terms, theprobability density functions in Figs. 9 and 10 cannot be compareddirectly. Instead, only the misfits (Figs. 9b and 10b) can be compareddirectly. Although both the Airy and Pratt end member scenarios permit
solutions within the measured uncertainties (Figs. 9a, 10a), the Airy
scenario yields much smaller misfits and thus accommodates the datamuch better. Moreover, it may be difficult for the Pratt mechanism to
dominate given that it requires the lateral density anomalies to persist to
depths of several tens of kilometers. If the density anomalies are due to9M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 9. The range of likely interior structures for Dione assuming Airy isostasy. (a) Curves illustrate the range of parameters that satisfy the observed J2 (blue) or C22(red) gravitational potential coefficients, with shaded bands indicating 1-σ uncertainties. (b) Misfit between model and observed gravitational potential. (c) Prob­
ability density function, with 68% (dark), and 95% (intermediate), and 99.7% (pale) confidence contours showing the parts of the parameter space that best fit theobserved gravitational potential. (d) Corresponding range of best fitting core radii and densities. Panels (a–d) assume ice and ocean densities of 925 kg/m3 and 1020
kg/m3, respectively. (e–f) Same as (c–d) except allowing for a wider range of possible ice and ocean densities (ice densities of 850–950 kg/m3; ocean densities of
1000–1100 kg/m3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)10M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 10. The range of likely interior structures for Dione assuming Pratt isostasy. (a) Curves illustrate the range of parameters that satisfy the observed J2 (blue) or
C22 (red) gravitational potential coefficients, with shaded bands indicating 1-σ uncertainties. (b) Misfit between model and observed gravitational potential. (c)Probability density function, with 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours showing the parts of the parameter space that best fit theobserved gravitational potential. (d) Corresponding range of best fitting core radii and densities. Panels (a–d) assume ice and mantle densities of 925 kg/m3 and1030 kg/m3, respectively. (e–f) Same as (c–d) except allowing for a wider range of possible densities (mean outer layer densities of 850–950 kg/m3; mantle densitiesof 1000–1300 kg/m3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)variations in porosity, they may not be able to reach such depths given
the high overburden pressures (>2 MPa per 10 km of depth, reaching~13 MPa at 60 km) and increasing temperatures with depth (~1 K/km
assuming a 120 km thick conductive shell), both of which effects would
tend to close those pores (Besserer et al., 2013). Dione’s excess topog­
raphy may therefore require some degree of Airy-type isostasy, meaningthat the outer ice shell could be underlain by a higher density, lower
viscosity layer—most straightforwardly interpreted as a subsurface
liquid water ocean (Fig. 11a).
If Dione does indeed harbor an internal liquid water ocean (Fig. 11a),then the approximately known temperature of the ice/ocean interface (i.
e., near 270 K) allows us to place a lower bound on the rate of heat loss.
In the most conservative case, where the entire ice shell behaves con­ductively, with conductivity going as T 1, the temperature structure canbe described by� �� �
� � R 1 R 1rdTbTðrÞ ¼ Ts
Ts(6)where Tb and Ts are the basal and surface temperatures, respectively, d isthe mean ice shell thickness, and R is Dione’s mean radius. The heat flux
at the surface is given by
� ��
�
c
Tb
d
F ¼ ln1(7)TsdR
where c is an empirical constant capturing the temperature dependenceof the thermal conductivity, k ¼ c/T, where we take c ¼ 651 W/m(Petrenko and Whitworth, 1999, p.43). Taking the mean surface
11

M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 11. Depictions of the interior of Dione, to scale. (a)Interior when Airy compensation is assumed. The model
matches the observed external shape and reproduces theobserved gravitational field exactly. The model comprises anouter ice shell with a mean thickness of 120 km, a subsurface
liquid water ocean with a mean thickness of 30 km, and a 411
km radius core, with a density of ~2300 kg/m3; the moment ofinertia factor is ~0.332. (b) Interior when Pratt compensation
is assumed. The best fitting Pratt model reproduces theobserved gravitational field within the uncertainties. The
model comprises a ~400 km radius core, with a density of~2400 kg/m3, and an H2O envelope with a mean thickness of~160 km; lateral density anomalies persist throughout thislayer (but see text); the moment of inertia factor is ~0.327.
Surface texture is a global visual map of Dione produced by
Paul Schenk (Lunar and Planetary Institute) from Cassini ISS
data (NASA, JPL).

temperature to be Ts ¼ 87K, and assuming an ice shell thickness of 120
km, the total conductive heat loss is approximately 19 GW (~4.8 mW/
m2). Although studies of elastic flexure (Hammond et al., 2013) andcrater relaxation (White et al., 2017) have suggested episodes of evengreater heat flows in Dione’s past, such an intense rate of heat loss isdifficult to sustain without the ocean freezing out. For Dione, assuming
an ice shell thickness of 120 km, ocean freezing occurs at a rate of ~40 m
per million years per 1 GW of net heat loss, meaning that a 30-km thickocean would freeze completely within ~40 Myrs if this 19 GW of heatloss were not balanced by internal heat production. Assuming a chon­
dritic heating rate of 4.5 � 10 12 W/kg (Spohn and Schubert, 2003) andDione’s rock mass of ~5.5 � 1020 kg (again assuming the chondritic rockmass density is 3500 kg/m3), radiogenic heat production is ~2 GW.
Given the recent finding that Saturn may be more dissipative thanpreviously thought (Lainey, 2016; Lainey et al., 2012, 2017, 2020; Fuller
et al., 2016), tidal heating is potentially large enough to balance the heatbudget, but more work is required to determine whether, how, andwhere tidal heating within Dione’s interior could be sufficient to prevent
an ocean from freezing completely. Maintaining an internal liquid water
ocean would be even more difficult if the ice shell were undergoing solidstate convection—a possibility that cannot be ruled out given that, un­less very high effective viscosities are assumed, such a thick ice shellyields large Ralyeigh numbers (e.g., Kamata, 2018).
If the ice shell is entirely conductive, then from Eq. (6), and again
assuming Ts ¼ 87K and an ice shell thickness of 120 km, we expect theupper 57–83 km to be cold enough (i.e., below 140–180 K) to behave
elastically, even on long timescales (e.g., Nimmo et al., 2002). On theother hand, our Airy isostasy model assumes effectively no elastic sup­port for the long wavelength (degree-2) topography. In principle,
including the effect of elastic support would shift our shell thickness
estimates to smaller values if the topography were generated at thesurface, or to larger values if the topography were generated due toheterogeneous freezing and melting at the ice/ocean interface (Hemi­
ngway and Mittal, 2019). Our model for complete Airy compensationand the implied present-day temperature structure are thus mutuallyconsistent only if the topographic anomalies formed when the effective
elastic layer thickness was negligible (i.e., at times of higher heat flow;
Hammond et al., 2013; White et al., 2017) or if the elastic layer is suf­
ficiently fractured to make long wavelength elastic support ineffective.
For example, following McKinnon (2013) and Beuthe et al. (2016), we
can estimate the tangential lithospheric stresses resulting from thedeflection required by the non-hydrostatic topography, assuming a
shear modulus of 3.5 GPa. For the case of bottom loading, which would
imply ~300 m of deflection to produce the degree-2 topography(Fig. 7c), we estimate stresses of ~30 MPa—far in excess of the ~1–2
MPa failure limit for intact ice (Schulson and Duval, 2009). In the case of
top loading, the lithospheric deflection required to satisfy the gravity
observations depends on the compensation depth: complete Airycompensation at a depth of 120 km produces an effect similar to that ofpartial compensation at a shallower depth. However, given that the nonhydrostatic gravity is ~0.65 of the uncompensated value, it can beshown that, even in the limit of very shallow compensation, the required
deflection is >1 km, leading to stresses of >100 MPa. Hence, regardlessof the loading history and temperature structure, elastic support is un­
likely to be effective for the observed degree-2 topography.
Since the high-standing non-hydrostatic topography is found on theleading and trailing faces, Airy isostasy would imply that these are thethickest parts of the ice shell (Fig. 12). Why the ice shell should bethickest on the leading and trailing faces, however, is not clear. Het­
erogeneous tidal dissipation within the ice shell should lead to lateralshell thickness variations (Ojakangas and Stevenson, 1989; Beuthe,
2013; Hemingway and Mittal, 2019), but the thickest parts of the shell
are expected to be at the prime- and anti-meridians, where tidal dissi­
pation is weakest, not on the leading and trailing faces. This 90-degree
longitudinal shift in the locations of the thickest parts of the shellmight be suggestive of non-synchronous rotation or an ice shell reor­
ientation event, perhaps contributing to the mismatch between thetheoretical and observed apex-antapex asymmetry in cratering (LeliwaKopystynski et al., 2012). It is also worth noting that the highest
standing topography, measured relative to the geoid, coincides with theprominent series of chasmata found on the trailing hemisphere (Fig. 7b;
Fig. 8). These features may be expected to exhibit lower densities andthus may be partly compensated in the Pratt sense. However, it is notclear whether such density anomalies could persist to the depths (>60km) required for this to significantly account for the observed
compensation.
Finally, we note that our best fitting core density of 2300–2400 kg/
m3 is very close to the inferred core density for Enceladus at 2340–2410
kg/m3 (Hemingway and Mittal, 2019) in spite of the two bodies having
somewhat different bulk densities (1478 kg/m3 for Dione, 1609 kg/m3
for Enceladus). Similarly, for Titan, whose bulk density is 1881 kg/m3,
the core density may be as low as ~2500 kg/m3, if one assumes a hy­
drostatic interior and a density of <1000 kg/m3 for the H2O envelope(Durante et al., 2019). Although we do not have a ready explanation forthe similarities in these core densities, we bring attention to thisobservation because, if it is not a mere coincidence, it could have im­plications for the formation processes of the Saturnian satellites and icy
moons in general.
4. Conclusions
We presented an estimation of the gravity field of Dione, obtained byanalyzing the Cassini Doppler tracking data acquired during five close
flybys of the moon. A full degree 2 field was sufficient to fit the data to
the noise level. The estimated values of the principal quadrupole terms,
J2 � 106 ¼ 1496 � 11 and C22 � 106 ¼ 364.8 � 1.8 (unnormalized
coefficients, 1-σ uncertainty), and their ratio, J2/C22 ¼ 4.102 � 0.044,
indicate a significant departure from the expectation for a body that has
12

M. Zannoni et al.

Icarus 345 (2020) 113713

Fig. 12. Dione’s ice shell structure assuming Airy isostasy (where the effects of topography beyond spherical harmonic degree 2 are excluded). The thickest parts of
the ice shell correspond with the regions of highest-standing non-hydrostatic topography, roughly centered on the leading and trailing faces.

relaxed to hydrostatic equilibrium. The departure from hydrostatic
equilibrium means that the moment of inertia cannot be inferred
directly from the gravitational field, but a combined analysis of gravity
and topography suggests a substantial degree of differentiation, with amoment of inertia factor of approximately 0.33. The analysis further
demonstrates that the high-standing topography is largely compensatedby some combination of lateral density anomalies and the deflection ofinternal density interfaces (i.e., Pratt and Airy isostatic support,
respectively), the latter mechanism being consistent with the presence ofan internal liquid water ocean. The Airy model better accommodates theobservational constraints, but an internal liquid water ocean would
imply rapid heat loss requiring some 17 GW of tidal heating to keep theocean from freezing.
Further insights into Dione’s interior may come from an improve­
ment of the topography or gravity models, and from the development of
more sophisticated geophysical and geochemical models. In particular,
reconstructing a coherent trajectory of Dione during the entire timespanof the Cassini mission has the potential to decrease the uncertainty on thegravity coefficients by approximately a factor of 2, though we note thatthe uncertainties in our interior models are rather dominated by un­
certainties in the shape models derived from limb profile analyses. Afterthe spectacular end of the Cassini mission in 2017, there are no currently
planned missions to study Dione or the other mid-sized icy moons ofSaturn. However, a future mission dedicated to a comprehensive char­acterization of these bodies could shed further light on their interiors,
leading to an improved understanding of the formation and evolution of
the Saturn system, and of icy moons in general.

the JPL Planetary Radar and Radio Science Group, the NASA/JPL DeepSpace Network and their operations personnel involved in the acquisi­tion of Doppler data analyzed here.
Data availability
The Doppler data and ancillary information used in this analysis arearchived in NASA’s Planetary Data System. https://pds.nasa.gov
References
Archinal, B.A., Acton, C.H., Hilton, J.L., Williams, I.P., Duxbury, T., Klioner, S.A.,
Consolmagno, G.J., Seidelmann, P.K., Tholen, D.J., A’Hearn, M.F., Kirk, R.L.,
Ping, J., Oberst, J., Meech, K., Thomas, P.C., McCarthy, D., Hestroffer, D.,
Conrad, A., 2018. Report of the IAU working group on cartographic coordinates androtational elements: 2015. Celest. Mech. Dyn. Astron. 130, 1–46. 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 observations. Radio
Sci. 40, 1–9. https://doi.org/10.1029/2004RS003101.
Bar-Sever, Y.E., Jacobs, C.S., Keihm, S., Lanyi, G.E., Naudet, C.J., Rosenberger, H.W.,
Runge, T.F., Tanner, A.B., Vigue-Rodi, Y., 2007. Atmospheric media calibration forthe deep space network. Proc. IEEE 95, 2180–2192. https://doi.org/10.1109/
JPROC.2007.905181.
Besserer, J., Nimmo, F., Roberts, J.H., Pappalardo, R.T., 2013. Convection-driven
compaction as a possible origin of Enceladus’s long wavelength topography. J.
Geophys. Res. E Planets 118, 908–915. https://doi.org/10.1002/jgre.20079.
Beuthe, M., 2013. Spatial patterns of tidal heating. Icarus 223 (1), 308–329. https://doi.
org/10.1016/j.icarus.2012.11.020.
Beuthe, Mikael, Attilio Rivoldini, and Antony Trinh. 2016. “Enceladus’ and Dione’s
Floating Ice Shells Supported by Minimum Stress Isostasy.” Geophysical Research
Letters 43. doi:https://doi.org/10.1002/2016GL070650.
Bierman, G.J., 1977. Factorization Methods for Discrete Sequential Estimation. Academic
Press, New York.
Campbell, J.K., Anderson, J.D., 1989. Gravity field of the Saturnian system from Pioneer
and Voyager tracking data. Astron. J. 97, 1485. https://doi.org/10.1086/115088.
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.
Di Benedetto, M., 2011. The non-gravitational accelerations of the Cassini spacecraft andthe nature of the Pioneer anomaly. PhD Thesis, Rome.
Dourneau, G., Baratchart, S., 1999. Masses and densities of the four inner major Saturn’ssatellites. Astron. Astrophys. 350, 680–684.
Durante, D., Hemingway, D.J., Racioppa, P., Iess, L., Stevenson, D.J., 2019. Titan’s
gravity field and interior structure after Cassini. Icarus 326, 123–132. https://doi.
org/10.1016/j.icarus.2019.03.003.
Ermakov, A.I., Park, R.S., Bills, B.G., 2018. Power laws of topography and gravity spectraof the solar system bodies. J. Geophys. Res. Planets 123, 2038–2064. https://doi.
org/10.1029/2018JE005562.
Evans, S., Taber, W., Drain, T., Smith, J., Wu, H.C., Guevara, M., Sunseri, R., Evans, J.,
2018. MONTE: the next generation of mission design and navigation software. CEAS
Sp. J. 10, 79–86. https://doi.org/10.1007/s12567-017-0171-7.

Acknowledgements
The authors are grateful to Luciano Iess, John W. Armstrong, andSami W. Asmar of the Cassini Radio Science Team, for the useful dis­cussions and suggestions regarding the procedures for Cassini data
analysis. We thank Peter Thomas for providing updated limb profiles
and Francis Nimmo for using these to generate the updated degree-2
shape model used in our analysis. M.Z., L.G.C., and P.T. acknowledge
support from the Italian Space Agency through the Agreement 2017-10H.O in the context of the NASA/ESA/ASI mission, and Caltech and theJet Propulsion Laboratory for granting the University of Bologna a li­
cense to an executable version of MONTE Project Edition S/W. D.H. wasfunded by the Miller Institute for Basic Research in Science at the Uni­
versity of California Berkeley, and by the Carnegie Institution for Sci­
ence in Washington, DC. We thank the Cassini Project, A. Anabtawi, and13

M. Zannoni et al.

Icarus 345 (2020) 113713
Lainey, V., et al., 2020. New tidal paradigm in giant planets supported by rapid orbital
expansion of Titan, Nature Astronomy (under Review).
Leliwa-Kopystynski, J., Banaszek, M., Wlodarczyk, I., 2012. Longitudinal asymmetry ofcraters’ density distributions on the icy satellites. Planetary and Space Science 60(1), 181–192.
McKinnon, W., 2013. The shape of Enceladus as explained by an irregular core:
Implications for gravity, libration, and survival of its subsurface ocean. Journal ofGeophysical Research: Planets 118. https://doi.org/10.1002/jgre.20122.
Milani, A., Gronchi, G., 2010. Theory of orbit determination. Cambridge University
Press, Cambridge. https://doi.org/10.1017/CBO9781139175371.
Modenini, D., Tortora, P., 2014. Pioneer 10 and 11 orbit determination analysis shows nodiscrepancy with Newton-Einstein laws of gravity. Phys. Rev. D 90, 022004. https://doi.org/10.1103/PhysRevD.90.022004.
Moyer, T. D., “Mathematical formulation of the Double-Precision Orbit Determination
Program (DPODP),” Technical Report 32–1527, Jet Propulsion Laboratory, Pasadena(CA), 1971.
Moyer, T.D., 2000. Formulation for Observed and Computed Values of Deep Space
Network Data Types for Navigation. Wiley, Hoboken (NJ). https://doi.org/10.1002/
0471728470.
Murray, C.D., Dermott, Stanley F., 1999. Solar system dynamics. Cambridge University
Press.
Nimmo, F., R. T. Pappalardo, and B. Giese (2002), Effective Elastic Thickness and Heat
Flux Estimates on Ganymede. Geophys. Res. Lett. 29 (7): 1158. doi: https://doi.
org/10.1029/2001GL013976.
Nimmo, F., Bills, B.G., Thomas, P.C., 2011. Geophysical implications of the longwavelength topography of the Saturnian satellites. J. Geophys. Res. 116 (E11),
E11001 https://doi.org/10.1029/2011JE003835.
Ojakangas, G.W., Stevenson, D.J., 1989. Thermal state of an ice shell on Europa. Icarus
81 (2), 220–241. https://doi.org/10.1016/0019-1035(89)90053-5.
Petrenko, V.F., Whitworth, R.W., 1999. Physics of Ice. OUP Oxford.
Plescia, J.B., 1983. The geology of Dione. Icarus 56 (2), 255–277. https://doi.org/
10.1016/0019-1035(83)90038-6.
Schulson, E.M., Duval, P., 2009. Creep and Fracture of Ice. Cambridge Univ, Press,
Cambridge, U. K.
Seidelmann, P.K., Abalakin, V.K., Bursa, M., Davies, M.E., De Bergh, C., Lieske, J.H.,
Oberst, J., Simon, J.L., Standish, E.M., Stooke, P., Thomas, P.C., 2002. Report of the
IAU/IAG working group on cartographic coordinates and rotational elements of theplanets and satellites: 2000. Celest. Mech. Dyn. Astron. 82, 83–110.
Serra, D., et al., 2019. A solution of Jupiter’s gravitational field from Juno data with the
ORBIT14 software. MNRAS 490 (1), 766–772. https://doi.org/10.1093/mnras/
stz2657.
Smith, B.A., et al., 1981. Encounter with Saturn: Voyager 1 imaging science results.
Science 212 (4491), 163–191. https://doi.org/10.1126/science.212.4491.163.
Spohn, T., Schubert, G., 2003. Oceans in the icy Galilean satellites of Jupiter? Icarus 161(2), 456–467. https://doi.org/10.1016/S0019-1035(02)00048-9.
Thomas, P., Burns, J., Helfenstein, P., Squyres, S., Veverka, J., Porco, C., Turtle, E.,
Mcewen, A., Denk, T., Giese, B., 2007. Shapes of the Saturnian icy satellites and theirsignificance. Icarus 190 (2), 573–584. https://doi.org/10.1016/j.
icarus.2007.03.012.
Tajeddine, R., Rambaux, N., Lainey, V., Charnoz, S., Richard, A., Rivoldini, A.,
Noyelles, B., 2014. Constraints on Mimas’ interior from Cassini ISS libration
measurements. Science 346 (6207), 322–324. https://doi.org/10.1126/
science.1255299.
Thomas, P.C., 2010. Sizes, shapes, and derived properties of the Saturnian satellites after
the Cassini nominal mission. Icarus 208 (1), 395–401. https://doi.org/10.1016/j.
icarus.2010.01.025.
Tortora, P., Zannoni, M., Hemingway, D.J., Nimmo, F., Jacobson, R.A., Iess, L., Parisi, M.,
2016. Rhea gravity field and interior modeling from Cassini data analysis. Icarus
264, 264–273. https://doi.org/10.1016/j.icarus.2015.09.022.
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(2), 99. https://doi.org/10.1088/0004-637X/782/2/99.
Vienne, A., Duriez, L., 1995. TASS 1.6: ephemerides of the major Saturnian satellites.
Astron. Astrophys. 297 (588–588).
White, O.L., Schenk, P.M., Bellagamba, A.W., Grimm, A.M., Dombard, A.J., Bray, V.J.,
2017. Impact crater relaxation on Dione and Tethys and relation to past heat flow.
Icarus 288, 37–52. https://doi.org/10.1016/j.icarus.2017.01.025.
Wieczorek, M., 2015. Gravity and Topography of the Terrestrial Planets. Treatise on
Geophysics 153–193. https://doi.org/10.1016/B978-0-444-53802-4.00169-X.
Wieczorek, M.A., Phillips, R.J., 1998. Potential anomalies on a sphere: applications to thethickness of the lunar crust. J. Geophys. Res. 103 (E1), 1715–1724.
Zannoni, M., Tortora, P., 2013. Numerical error in interplanetary orbit determination
software. J. Guid. Control. Dyn. 36, 1008–1018. https://doi.org/10.2514/1.59294.
Zannoni, M., Tommei, G., Modenini, D., Tortora, P., Mackenzie, R., Scoubeau, M.,
Herfort, U., Carnelli, I., 2018. Radio science investigations with the Asteroid impact
mission. Adv. Sp. Res. 62, 2273–2289. https://doi.org/10.1016/j.asr.2017.12.003.

Fuller, J., Luan, J., Quataert, E., 2016. Resonance locking as the source of rapid tidal
migration in the Jupiter and Saturn moon systems. Mon. Not. R. Astron. Soc. 458,
3867–3879. https://doi.org/10.1093/mnras/stw609.
Gomez Casajus, L., et al., 2020. Updated Europa Gravity Field and Interior Structure froma Reanalysis of Galileo Tracking Data, Icarus (under review).
Hammond P., N., Phillips B., C., Nimmo, F., Kattenhorn A., S., 2013. Flexure on Dione:
Investigating Subsurface Structure and Thermal History. Icarus 223 (1), 418–422.
https://doi.org/10.1016/j.icarus.2012.12.021. In press.
Harper, D., Taylor, D.B., 1993. The orbits of the major satellites of Saturn. Astron.
Astrophys. 268, 326–349.
Hemingway, D., Nimmo, F., Iess, L., 2013. Enceladus’ Internal Structure Inferred from
Analysis of Cassini-Derived Gravity and Topography, P53E-03 (in AGU Fall Meeting2013).
Hemingway, D.J., Matsuyama, I., 2017. Isostatic equilibrium in spherical coordinates
and implications for crustal thickness on the Moon, Mars, Enceladus, and elsewhere.
Geophys. Res. Lett. 44, 7695–7705. https://doi.org/10.1002/2017GL073334.
Hemingway, D.J., Mittal, T., 2019. Enceladus’s ice shell structure as a window oninternal heat production. Icarus 332, 111–131. https://doi.org/10.1016/j.
icarus.2019.03.011.
Hemingway, D.J., Zannoni, M., Tortora, P., Nimmo, F., Asmar, S.W., 2016. Dione’s
internal structure inferred from Cassini gravity and topography. In Lunar andPlanetary Science Conference 47, 1314.
Iess, L., Budnik, F., Colamarino, C., Corbelli, A., Di Benedetto, M., Fabbri, V., Graziani,
A., Hunt, R., James, N., Lanucara, M., Madd�e, R., Marabucci, M., Mariotti, G.,
Mercolino, M., Racioppa, P., Simone, L., Tortora, P., Westcott, M., Zannoni, M.,
2012a. ASTRA: Interdisciplinary Study on Enhancement of the End-To-End Accuracy
for Spacecraft Tracking Techniques, in: 63rd International Astronautical Congress.
pp. 1–11.
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., 2012b. The gravity field of
the Saturnian satellites Enceladus and Dione, in: American Geophysical Union, Fall
Meeting 2012.
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., 2012c. The tides of Titan.
Science 337, 457–459. https://doi.org/10.1126/science.1219631.
Iess, L., Stevenson, D.J., Parisi, M., Hemingway, D.J., Jacobson, R.A., Lunine, J.I.,
Nimmo, F., Armstrong, J.W., Asmar, S.W., Ducci, M., Tortora, P., 2014. The gravity
field and interior structure of Enceladus. Science 344, 78–80. https://doi.org/
10.1126/science.1250551.
Iess, L., Folkner, W.M., Durante, D., Parisi, M., Kaspi, Y., Galanti, E., Guillot, T.,
Hubbard, W.B., Stevenson, D.J., Anderson, J.D., Buccino, D.R., Gomez Casajus, L.,
Milani, A., Park, R.S., Racioppa, P., Serra, D., Tortora, P., Zannoni, M., Cao, H.,
Helled, R., Lunine, J.I., Miguel, Y., Militzer, B., Wahl, S., Connerney, J.E.P., Levin, S.
M., Bolton, S.J., 2018. Measurement of Jupiter’s asymmetric gravity field. Nature
555. https://doi.org/10.1038/nature25776.
Iess, L., Militzer, B., Kaspi, Y., Nicholson, P.D., Durante, D., Racioppa, P., Anabtawi, A.,
Galanti, E., Hubbard, W.B., Mariani, M.J., Tortora, P., Wahl, S., Zannoni, M., 2019.
Measurement and implications of Saturn’s gravity field and ring mass. Science 2965,
eaat2965. https://doi.org/10.1126/science.aat2965.
Jacobson, R.A., 1995. The orbits of the minor Saturnian satellites. Bull. Am. Astron. Soc.
27, 22.
Jacobson, R.A., 2004. The orbits of the major Saturnian satellites and the gravity field ofSaturn from spacecraft and earth-based observations. Astron. J. 128, 492–501.
https://doi.org/10.1086/421738.
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.W., 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.
Kamata, S., 2018. One-Dimensional Convective Thermal Evolution Calculation Using a
Modified Mixing Length Theory: Application to Saturnian Icy Satellites. Journal ofGeophysical Research: Planets 123, 93–112. https://doi.org/10.1002/
2017JE005404.
Kaula, W.M., 1963. Determination of the Earth’s gravitational field. Rev. Geophys.
https://doi.org/10.1029/RG001i004p00507.
Kozai, Y., 1976. Masses of satellites and oblateness parameters of Saturn. Publ. Astron.
Soc. Japan 28, 675–691.
Lainey, V., 2016. Quantification of tidal parameters from Solar System data. Celest.
Mech. Dyn. Astron. 126 (1), 1–12. https://doi.org/10.1007/s10569-016-9695-y.
Lainey, V., et al., 2012. Strong tidal dissipation in Saturn and constraints on Enceladus’
thermal state from astrometry. Astrophys. J. 752 (14) https://doi.org/10.1088/
0004-637X/752/1/14.
Lainey, V., Jacobson, R.A., Tajeddine, R., Cooper, N.J., Murray, C., Robert, V., Tobie, G.,
Guillot, T., Mathis, S., Remus, F., Desmars, J., Arlot, J.E., De Cuyper, J.P., Dehant, V.,
Pascu, D., Thuillot, W., Le Poncin-Lafitte, C., Zahn, J.P., 2017. New constraints on
Saturn’s interior from Cassini astrometric data. Icarus 281, 286–296. https://doi.
org/10.1016/j.icarus.2016.07.014.

14