RESEARCH ARTICLE

The Global Shape, Gravity Field, and Libration of Enceladus

Key Points:

R. S. Park1 , N. Mastrodemos1, R. A. Jacobson1, A. Berne2 , A. T. Vaughan1,
D. J. Hemingway3 , E. J. Leonard1 , J. C. Castillo-Rogez1, C. S. Cockell1,4 , J. T. Keane1
A. S. Konopliv1 , F. Nimmo5 , J. E. Riedel1, M. Simons2 , and S. Vance1

10.1029/2023JE008054• A full quadrupole gravity field with J3and the forced libration amplitude of
0.091° ± 0.009° are recovered• A 500-m resolution global topographymodel was computed, with some local
regions having 25−100 m resolution• The results suggest that Enceladus has
a 27–33 km mean ice shell thickness,a 21–26 km ocean thickness,
and a mean core density range of
2,270–2,330 kg/m 3
Supporting Information:

Supporting Information may be found in
the online version of this article.
Correspondence to:

R. S. Park,
Ryan.S.Park@jpl.nasa.gov
Citation:

Park, R. S., Mastrodemos, N., Jacobson,
R. A., Berne, A., Vaughan, A. T.,
Hemingway, D. J., et al. (2024). The
global shape, gravity field, and librationof Enceladus. Journal of GeophysicalResearch: Planets, 129, e2023JE008054.
https://doi.org/10.1029/2023JE008054
Received 14 AUG 2023
Accepted 20 DEC 2023
Author Contributions:Conceptualization: R. S. Park, J. T.
Keane, A. S. Konopliv, J. E. Riedel
Formal analysis: R. S. Park, N.
Mastrodemos, R. A. Jacobson, A. Berne,
A. T. Vaughan, D. J. Hemingway, E. J.
Leonard, J. T. Keane, S. Vance
Funding acquisition: R. S. Park
Investigation: R. S. Park, J. C.
Castillo-Rogez
Methodology: R. S. Park, N.
Mastrodemos, R. A. Jacobson, A. Berne,
D. J. Hemingway, E. J. Leonard, J. C.
Castillo-Rogez, J. T. Keane, F. Nimmo,
S. Vance© 2024 Jet Propulsion Laboratory,
California Institute of Technology and
The Authors. Government sponsorship
acknowledged.
This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial License,
which permits use, distribution andreproduction in any medium, provided theoriginal work is properly cited and is not
used for commercial purposes.

PARK ET AL.

,

1
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, 2Seismological Laboratory, Divisionof Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA, 3Institute for Geophysics,
University of Texas at Austin, Austin, TX, USA, 4UK Centre for Astrobiology, School of Physics and Astronomy, Universityof Edinburgh, Edinburgh, UK, 5Department of Earth and Planetary Sciences, University of California Santa Cruz, Santa
Cruz, CA, USA

Abstract In order to improve our understanding of the interior structure of Saturn's small moon Enceladus,
we reanalyze radiometric tracking and onboard imaging data acquired by the Cassini spacecraft duringclose encounters with the moon. We compute the global shape, gravity field, and rotational parameters ofEnceladus in a reference frame consistent with the International Astronomical Union's definition, where
the center of the Salih crater is located at −5° East longitude. We recover a quadrupole gravity field with J3
and a forced libration amplitude of 0.091° ± 0.009° (3-σ). We also compute a global shape model using astereo-photoclinometry technique with a global resolution of 500 m, although some local maps have higher
resolutions ranging from 25 to 100 m. While our overall results are generally consistent with previous studies,
we infer a thicker 27–33 km mean ice shell, a thinner 21–26 km mean ocean thickness, and a mean core density
range of 2,270–2,330 kg/m 3.
Plain Language Summary Geodetic data, such as shape, gravity, and rotation, provide important
constraints for probing a planetary body's interior structure. We analyze radiometric tracking and onboard
imaging data acquired during close encounters of Enceladus by the Cassini spacecraft to compute geodetic
products including topographic and gravitational fields in a common reference frame. The recovered Enceladus
topography has a global resolution of 500 m, with some local regions having 25–100 m resolution. Our study
suggests that Enceladus has a 27–33 km mean ice shell thickness, a 21–26 km mean ocean thickness, and amean core density range of 2,270–2,330 kg/m 3.
1. Introduction
Enceladus is a small, 500-km-diameter moon orbiting Saturn approximately every 33 hr. It is a geologically active
body with material jetting from crustal fractures over the body's highly tectonized South Polar Terrain (e.g.,
Crow-Willard & Pappalardo, 2015; Yin & Pappalardo, 2015). Erupted material collectively forms a broad plumethat is visible over the moon's South Pole (Hansen et al., 2006; Ingersoll et al., 2020; Porco et al., 2006; Spencer
et al., 2006). Between 2004 and 2017, the Cassini spacecraft made several close encounters with Enceladus,
resulting in various types of measurements including Deep Space Network (DSN) radiometric tracking data andimaging data. The DSN data are the primary observations used to determine the gravity field of Enceladus (Iesset al., 2014), whereas the imaging data are the primary observations used for determining shape and rotationalparameters (Bland et al., 2018; Nimmo et al., 2011; Schenk & McKinnon, 2024; Tajeddine et al., 2017; Thomas
et al., 2016).
Geodetic data, such as shape, gravity, and rotation, provide important constraints for probing a planetary body's
interior structure (Durante et al., 2019; Ermakov et al., 2021; Iess et al., 2014; Park et al., 2014, 2016; Park,
Konopliv, et al., 2020; Park, Riedel, et al., 2020). Results from Cassini measurements led to the inference that
Enceladus has a differentiated interior with a rocky core and a global subsurface liquid water ocean overlaid with
an icy shell (Beuthe et al., 2016; Hemingway & Mittal, 2019; Van Hoolst et al., 2016). Enceladus's jets may erupt
material that is sourced from the ocean (Postberg et al., 2009) although the provenance of this material throughthe crust is not well known.
Here, we analyze Cassini radiometric tracking and onboard imaging data acquired during close encountersof Enceladus using the latest reconstructed ephemerides of Enceladus and the Cassini spacecraft (Jacobson
1 of 27

Project Administration: R. S. Park
Resources: R. S. Park
Software: R. S. Park, N. Mastrodemos,
R. A. Jacobson, A. Berne, A. T. Vaughan,
D. J. Hemingway
Supervision: R. S. Park
Validation: R. S. Park, N. Mastrodemos,
R. A. Jacobson, A. Berne, A. T. Vaughan,
D. J. Hemingway, E. J. Leonard, J. C.
Castillo-Rogez, F. Nimmo, S. Vance
Visualization: R. S. Park
Writing – original draft: R. S. Park,
N. Mastrodemos, R. A. Jacobson, A.
Berne, A. T. Vaughan, D. J. Hemingway,
E. J. Leonard, J. C. Castillo-Rogez, J. T.
Keane, A. S. Konopliv, F. Nimmo, J. E.
Riedel, S. Vance
Writing – review & editing: R. S. Park

10.1029/2023JE008054

et al., 2022). We compute the gravity field, shape, and rotational parameters of Enceladus in a reference frame
consistent with the International Astronomical Union's definition, where the center of the Salih crater is located at−5° East longitude (Archinal et al., 2018). In this way, both the shape and gravity field share the same reference
frame, allowing self-consistent geodesy products for exploring the interior of Enceladus. We present an updated
quadrupole gravity field with J3 and a forced libration amplitude of 0.091° ± 0.009°, including a complete rotational model of Enceladus. The interpretation of our results is generally mostly consistent with previous studies
but suggests a thicker 27–33 km mean ice shell, a thinner 21–26 km mean ocean thickness, and a mean coredensity range of 2,270–2,330 kg/m 3.
We also present an updated global shape model using a stereo-photoclinometry technique at a global resolutionof 500 m, with some local regions having resolutions as high as 25–100 m. All shape and orientation products areavailable through the Navigation and Ancillary Information Facility (NAIF) Planetary Data System (PDS) node
for interior studies and planning for future mission to Enceladus.

2. Stereophotoclinometry
We construct the shape of Enceladus using stereo-photoclinometry (SPC), a methodology that combines aspectsof photogrammetry and photoclinometry. SPC has been used over the last 20 years for both proximity navigationand shape reconstruction (Gaskell et al., 2023; Mastrodemos et al., 2012; Park et al., 2019). Global shape models
derived using SPC are built piecemeal from small digital terrain map units referred to as landmark maps. Each
such map is a square 2-dimensional grid. The location of the center of the map is determined by a control point(or landmark vector) from the body center to the map center.
Each map grid point has a local height defined relative to the control point, local slopes defined along the map
grid coordinate directions, and a relative albedo value. We combine the central control point, map grid heights,
and a local map coordinate system to compute a vector from the body center to each grid height. Control points
are not chosen based on surface features. Instead, we use a tiling scheme of overlapping maps based on the desired
surface density of control points and the degree to which neighboring maps are chosen to overlap. Typical tilingschemes have adjacent maps overlap by 25%–50%, such that each map shares part of its surface with at least 3
to 4 adjacent maps. For each map pixel, we estimate the slopes, albedo and local height from the imaging data
as follows:
1. A map grid is specified at a location on the surface, typically based on latitude and longitude coordinateswith a set number of pixels and a map scale. The map control point and therefore the overall map height are
initialized based on the pre-existing shape. Map heights and albedo are also given initial values from existing
shape (e.g., an a priori shape model, or a pre-existing lower resolution map). Slope values at each map point
are constructed by differencing the a priori map heights.
2. Based on a priori camera position and pointing all images that overlap some portions of the map are identified. These candidate images are filtered for certain geometric conditions such as incidence angle, emission
angle, the ratio of image pixel resolution to map pixel scale, and the ratio of the image area to the map area
overlap. Typical incidence angles are in the range of 10°–80°, but this range can be extended if there are notenough images. Typical emission angles are in the range of 0°–60°. Typical image resolution to map the pixel
scale ratio is in the range of 1/3–3. However, sometimes higher resolution images may be added to lower
resolution maps to improve their pointing solution. In the case of Enceladus, images with a factor of up to 20
times higher resolution have been included for certain maps. Also, if there is a sufficient number and diversityof images at a pixel scale comparable to the map scale, images with a pixel scale larger than 2 times of themap scale may be discarded. In case of manual processing, candidate images may be further filtered based on
visual appearance such as artifacts, image smear, or a large part of the image overlaying the map appearing
too dark.
3. The selected images are orthorectified on the a priori map surface model, rescaled and resampled to the map
pixel scale. From this point on, we work with these orthorectified image templates and not the original images.
4. A model brightness profile is constructed at each map pixel based on the a priori image geometry, a priori
map pixel slopes, albedo, heights, and a model phase function. The differences between the model brightness profile and the image template brightness profile are minimized in a weighted least squares sense tosolve for the new values of the slopes and albedo. Due to position and pointing errors image templates are
initially misaligned relative to each other. The derived brightness profile is used to align image templates via
PARK ET AL.

2 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 1. A partial template of 10 Cassini images projected onto a map (top row) with the reconstructed map rendered at the viewing geometry of the corresponding
image (bottom row). Notice the difference in the appearance of many surface features between the first six and the last four images due to the differences in emission
and phase angles.

cross-correlation and a new map brightness profile is solved. This process may be iterated many times until
all images are aligned to sub-pixel precision relative to each other and the landmark map. In the initial stagesof the shape reconstruction, images that correlate poorly may have to be aligned manually. Such examples
may include images of lower phase or of high emission compared to the rest. We note that the albedo values
estimated at each map pixel are not in absolute radiometric units but are relative to the mean value. We estimate albedo information so that we can separate the brightness variation due to the local topography frombrightness due to intrinsic albedo variations.
5. The slope solution is integrated into local heights using seed points from existing topography, such as an apriori shape model, existing overlapping maps and image limbs.
6. After the local heights are constructed, further constraints are added, such as the location of the map on thelimbs of images and the vectors between control points of overlapping maps.
Steps 4–6 above may have to be iterated a few times. This entire process is repeated for each map until a sufficiently large area of the surface has been modeled. A minimum of 3 images are needed to estimate the map slopes
and albedo. There is no maximum number of images allowed, but in practice redundant images may be removed
for computational considerations.
The choice of map scale is determined from the resolution range of the imaging data. Often, a single map scale
is not sufficient to model the range in the image resolution, and layers of maps at differing scales are constructedto cover the body in increasing resolution steps, starting with coarser map scales. Typically, the highest map
resolution constructed over a particular area is comparable to the highest image resolution, for areas that have aminimum of 3 such images of comparable resolution.
Once many maps have been constructed, global parameters are estimated. These parameters are the map control
points, the camera pointing angles, and the camera-satellite vector, subject to the a priori ephemeris and pointing
uncertainties. Other constraints on this solution include the location of maps in the image limbs and the overlapping maps' control points relative position. Global parameter estimation is followed by inspection of the posteriori
control point uncertainties and the residuals between control points and images, limb locations and overlapping
maps. At this stage, outliers are inspected and manually corrected.
This formulation allows for a self-consistent closed-loop process with the estimation of spacecraft, satellite ephemerides and key physical parameters; the map control points and their image location are the optical observables
that are used in the orbit determination process combining optical and tracking data.
2.1. Cassini Image Data Set
The image set we use was acquired by the Cassini spacecraft between January 2005 and November 2016. The
surface coverage is nearly complete but highly non-uniform with some areas having redundant coverage in both
number of images and spatial resolution. Some areas have marginal coverage, requiring the extension of theallowable emission angle to the 65°–75° range to have the required minimum of 3 images. Incidence and emission
angles are often highly suboptimal and many areas were imaged under very few phase angles, often only two. The
latter makes the early stages of processing manually intensive with different images showing different surface
features in the same physical area. We show an example in Figure 1. The first row in Figure 1 shows 10 of the
Cassini images used in a particular map, after rescaling and projection on the reconstructed map. The second row
PARK ET AL.

3 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

shows the map rendered at the image geometry. The first six images are all at 46° phase angle and ∼52° emission
angle, whereas the last 4 are at 73° phase angle and ∼15° emission angle. Many surface features that are present
in one of these two groups of images are barely visible or appear absent from the other group and vice versa. For
these images, the initial steps of cross-correlation and image matching have been challenging, often requiring
manual image matching via a graphical user interface. Similar image correlation problems were also mentioned
by other authors (e.g., Bland et al., 2018).
All clear filter and monochromatic narrow-angle camera (NAC) and wide-angle camera (WAC) images with asurface pixel scale between 3 km and 14 m were initially considered as candidates. A number of these images
were later discarded due to considerations such as poor correlation, low resolution where higher resolution
images were available, large emission angles, very high phase angles (i.e., larger than 120°), or a large number ofartifacts or cosmic rays. The final number of images that contributed to the shape model was 1090 NAC and 21
WAC images, ranging in pixel scale from 14 m to 2.2 km.
2.2. Initial Conditions and Set Up
In addition to the images, key ancillary data included the a priori pose of both cameras. Inertial pointing was
based on the Cassini reconstructed attitude kernels and the camera-spacecraft relative orientation available atthe NAIF website (https://naif.jpl.nasa.gov/pub/naif/CASSINI). The Enceladus-spacecraft relative position wasderived from the SAT441 Saturnian system ephemerides (Jacobson et al., 2022), DE440 planetary ephemerides(Park et al., 2021), and the associated reconstructed spacecraft trajectory from the SAT441 solution. The rotational frame was improved throughout the shape development process, and the updated versions were incorporated to rotate between inertial and body-fixed frames and reconstruct the topography model. The camera models
used were those developed by the Cassini Optical Navigation team (Gillam et al., 2007). A smooth triaxial ellipsoid with axes (a = 256.6, b = 251.4, c = 248.3) km was used for the a priori shape.
The simultaneous estimation of correlated quantities, such as control point locations, camera pointing directions
and camera-Enceladus position, requires careful conditioning to avoid systematic aliasing between the pointing
and the position correction. A priori uncertainties for the camera's inertial pointing were determined by projecting
images on the surface and measuring the scatter in the projected image locations relative to the surface features.
Based on that analysis, a 1-σ a priori pointing uncertainty of 10 −4 radians per axis was set, equivalent to ∼16 NAC
and 1.6 WAC pixels. The a priori camera position error was derived from the combined spacecraft and Enceladus
ephemerides expected uncertainties and was refined by parametric analysis that evaluated the global post-fit
residuals between topography map control points, camera pointing and position. A 1-σ a priori position error percomponent of 1 km was found to give a better goodness-of-fit while remaining consistent with the ephemeris
error. Therefore, for most images, the pointing error is the largest contributor to the input error.
We note that the projection error of an image on the surface has a different signature depending on whether it
is due to camera pointing error or camera position error. Accounting for both pointing and position errors in
the global parameter estimation is essential to correctly project and align images with each other and with the
topography maps. We scale the initial control point uncertainties in both height and in horizontal directions along
latitude and longitude to the map pixel scale by a factor of 20, that is, for the 1-km scale maps, 20 km height relative to the a priori ellipsoid and approximately 4.5° along latitude and longitude at the equator. For the apparent
height excursions of Enceladus, this uncertainty scaling is a very loose constraint that will allow the topographyheight to be driven by the imaging data rather than being artificially constrained to remain very close to the ellipsoid heights. In addition to this loose constraint, the control point solution is also subject to the location of maps
in the image limbs and, more importantly, to the constraints of the connecting vectors between control points ofoverlapping maps. In essence, overlapping maps enforce a gradual and self-consistent height adjustment over theentire modeled surface, thereby preventing outlier heights from appearing at any point in the process.
SPC requires that the reflectance model and phase function be valid and appropriate. Even though a wide rangeof such phase functions may produce good visual results, depending on the surface properties of the body, someare more suitable than others. We adopted the previously published Minnaert photometric model, which was
originally derived from the Enceladus Cassini images (see Equation 1 and Table 3 in Bland et al. (2018)). To
understand the sensitivity of the shape model to the choice of the reflectance model, we constructed three models
each with different reflectance functions and compared these to the Minnaert model. The reflectance functions
PARK ET AL.

4 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

used were the Lunar-Lambert (McEwen, 1991), Lommel-Seeliger with no phase term and the Ceres parameterized Akimov disk model (Schroder et al., 2017), as an example of a reflectance function of an icy body. These
three models when converged resulted in posteriori residuals that were larger than those found using the Minnaert
model, by factors in the range of 1%–15%. Thus, we adopted the Minnaert photometric function as our baselinemodel.
2.3. Shape Model Development
SPC is an iterative process with a few steps repeated many times until final model convergence. We use thefollowing general procedure in the construction of the Enceladus model:
1. Tile the entire surface with overlapping maps of 99-by-99 pixel format at a 1-km-pixel scale. Converge thesemaps according to the topography process described above. Initial seed heights for these maps were derived
from the a priori triaxial ellipsoid.
2. Estimate control point vectors (i.e., landmarks), camera pointing, and camera position.
3. Inspect residuals and treatment of outliers. This step includes a variety of corrective actions: manually adjusting poorly correlated images, deleting images with large numbers of outliers from selected maps, adding more
maps in areas with few overlapping maps, adding more maps in areas with large surface inhomogeneities orwith large residuals in the adjacent map overlaps, or deleting images whose large residuals could not be reconciled otherwise. The guiding principle for deciding whether an image with a large residual should be deleted
was the evaluation of redundancy versus uniqueness of the observation geometry from the SPC observational
merit evaluation: images with large residuals in a particular map for which there was redundancy in terms ofimages with similar resolution, emission, incidence, camera, and solar azimuth angles in the same map weredeleted. Images with large residuals that offer a unique observational geometry in that particular map wereretained. In general, images with residuals up to a value of 3-σ of the total root-mean-square (RMS) of theresiduals were used. For several maps, the addition of a few images at a much higher resolution than the map
scale, where available, improved the map's spatial resolution but at the expense of larger residuals for those
few high-resolution images. Our overall strategy therefore was to maximize the diversity of the observational
geometry at the expense of somewhat higher posteriori residuals.
4. Rebuild the model. Because the map's slopes depend on the accuracy of the camera geometry, after updating
the camera position and pointing, all maps were reconstructed.
Steps 2–4 above were repeated many times until convergence was reached. Metrics for convergence are theimage-control point residuals, the residuals between overlapping maps, the posteriori control point covariance, as
well as visual inspection of the maps to ensure that the map spatial resolution captures most image details. Oncea complete model is converged at the 1 km map scale, the whole process is repeated with maps at 500 m scale
globally, which represents the final global model resolution. Initial seed heights for the 500 m maps are derived
from the shape model constructed from the 1 km maps. The surface coverage and distribution of the 500 m mapson the surface are shown in Figure 2 in grayscale. Lighter-colored areas indicate regions where a higher densityof maps was necessary to improve the model.
In addition to the global model, images at resolutions higher than ∼500 m were acquired at several distinct
locations that allowed regional and local topography to be constructed at higher resolutions. These areas weremodeled to the highest map spatial resolution level allowed by the image resolution. One such contiguous areaof particular interest is the South Pole, which we modeled at 200 m resolution from 90°S to latitudes varying
from 70° to 35°S depending on the longitude. Interior to that area, a smaller area was modeled at 100 m from
90°S to the 80°S to 70°S range, again depending on the longitude. Furthermore, 13 very high resolution images,
from 59 to 14 m pixel scale, partially overlap in the Baghdad Tiger Stripes region that allowed an area centered
at approximately (79°S, 335°E) to be constructed at 50 m scale and another smaller area centered approximately
at (79°S, 325°E) to be constructed at 25 m map scale. Most of these images have high emission angle and highphase angle with parts of the images being too dark to be used. The combination of high phase and large verticalfeatures create self-shadowing, causing the quality of topographic reconstruction in these areas to be quite variable from point to point. In addition, only a few images, typically 3–5, contribute to individual topography maps,
making the topography reconstruction challenging. Figures 3a–3c show renderings of the topography at 100 m,
50 m, and 25 m scales, respectively. Even at these adverse imaging conditions, a high level of detail can be seenin the topography. For example, Figures 4a and 4b display the topography in that area rendered at the viewing andillumination geometry of the respective images, as shown in Figures 4c and 4d.
PARK ET AL.

5 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 2. The Enceladus surface coverage with 500 m maps in gray scale in the cylindrical projection.

We also model the area bounded by 132°−142°E and 1°–10°S at a high resolution of 35 m map scale. That area
has been imaged by ∼10 partially overlapping images in a pixel scale range of 28–70 m. In contrast to the Baghdad region, most of these images are of low emission and low phase angle. The topography of this area is shownin Figure 5a, rendered at the geometry of the Cassini image N1489051187, one of the high-resolution images
covering that area, which is shown in Figure 5b for comparison. The number of maps (and therefore control
points) produced for each map scale is shown in Table 1, together with the average number of images contributing
to the topography per map scale. The map scale should be considered as a good approximation of the model's
spatial resolution at that location. A summary of the image statistics is shown in Table 2. As can be seen in thelast column, there are a sufficient number of control points per image to ensure a robust estimation of the camera
pointing angles and position.
The statistics for the global topography model from the 1 km and 500 m resolution maps in the posteriori
covariance and posteriori residuals between control points and images are summarized in Table 3. For theposteriori covariance, we report the root-sum-square of the diagonal elements of the control point vectors as
the height uncertainty at that location—the most objective representation of the absolute height error in thatlocation relative to the center of the satellite. The total number of observations is the product of the images with
the number of constraints per image, that is, control points, limb points and map overlaps. Based on the above
results, we expect that the best estimate for the average height accuracy of the global model is in the 250–300 m
range (1-σ).

Figure 3. (a) The south pole topography at 100 m scale. This large scale 150 km by 150 km map has been synthesized fromall individual 100 m maps in the tiger Stripes region. (b) A map of a section of the Baghdad region, modeled at 50 m scale.
The irregular structure of the map is due to the surface overlap of images of spatial scale in the 14–59 m range. (c) Map of asubregion of the Baghdad region shown in panel (b), constructed at 25 m scale.

PARK ET AL.

6 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 4. (a) Rendered topography in the Tiger Stripes region at 25–50 m scale shown at the viewing geometry of the
Cassini image N1637462854. The average emission angle for maps contributing to the topography is 65°. (b) Rendered
topography in the Tiger Stripes region at 25–50 m scale shown at the viewing geometry of the Cassini image N1637462964.
The decrease in spatial resolution on the low and right sides of the image is due to the lack of 25 m topography in theseareas. The average emission angle for the maps contributing to the topography is 60°. (c) High-resolution Cassini image
N1637462854 at 19 m scale at an enhanced brightness level from its original version for comparison against the topography
map. (d) High-resolution Cassini image N1637462964 at 15 m scale at an enhanced brightness level from its original version
for comparison against the topography map.

Figure 5. (a) A rendered topography map in the (1°–10°S), (218°–228°W) region at 35 m scale shown from the viewing
geometry of the Cassini image N1489051187. (b) Actual Cassini nadir image N1489051187 at 34 m scale.

PARK ET AL.

7 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

Table 1
Number of Individual Topography Maps and Associated Control Points, per
Map Scale
Number of maps (control points)

Average numberof images per map

1,000

335

82

500

1,105

70

200

923

100

739

50

494

35

258

25

438

Map scale (m)10.1029/2023JE008054

2.4. SPC Merit Study
The suitability of an imaging data set for a stereophotoclinometry shapemodel is a multivariate problem that is difficult to quantify. We have derived
an empirical method for evaluating the strength of a given data set for SPC
processing that results in a scalar figure of merit for a specified region on thesurface.

In general, for SPC processing, we want a variety of illumination conditions
and a variety of viewing angles in the imagery but defining “variety” is diffi11
cult. The incident, emission, and phase angles for a set of images that havevisibility of a surface location are scalar quantities that do not fully describe4
the solar geometry. For example, the incidence angle could be 45° for several
4
pictures, which would seem deficient, but if they were a mix of morning4
and afternoon illumination angles, then there was “variety” in these images.
The phase angle could be very low in a set of candidate images if the sun is“at your back” in all of them, even though the incident angles are diverse.
Adding solar azimuth and spacecraft azimuth to the metric helps to ensure a variety of illumination and viewing
angles, but again the azimuth angles alone are not sufficient to guarantee a meaningful data set. We have foundthat each of these terms has something to contribute to the overall SPC merit, and so each of them is included in
the analysis.
25

The general criteria used for evaluating an SPC data set depend primarily on the imaging geometry (Jaumann
et al., 2022). The specific imaging angle criteria that we have implemented for assessing the SPC dataset aresummarized in Table 4; however, the discrete binning method has significant drawbacks. Instead, we propose acoverage approach to evaluating each of these criteria. For each metric, we establish a scorable range and thenoverlay the computed angle from each image as a measurement that covers some part of that range. By carefully
defining the width of the coverage that each measurement provides to each metric and combining multiple measurements with a maximum function, we can determine the coverage of each angle range without rewarding or
penalizing redundant data.
For most metric types, we use a sinusoid with a configurable width to impart coverage on the scorable arc.
Figure 6 is an illustration demonstrating how we apply coverage to the scorable range of an incidence angle metricin a notional case. Two images that have visibility of the surface patch in question have two different complementary incident angles. They both apply their coverage curves to the scorable range, and the maximum value is kept
in each location, resulting in full coverage for regions near the measurements and partial coverage as the sinusoid
tails off. This algorithm has the important feature of not rewarding duplicate measurements: a third measurementthat occurs very near one of the existing measurements could widen the coverage marginally but would not havea significant impact on the final score.
The final scalar result for each angle metric comes after all the measurements have been applied, and we divide
the summed score by the best possible score over the allowable range. We have calibrated the width of the curves
such that the resulting scores match our intuition and the stated requirements from Table 4. For example, the solar
azimuth score approaches 1.0 when there are measurements equally spaced in each quadrant of the 360° scorable
range. The score is approximately 0.5 when only two quadrants are covered.
Importantly, the “coverage” approach to scoring (instead of the discrete
Table 2“buckets” such as those described in Jaumann et al. (2022)) is ambivalent toImage Scale Distribution and Statistics of Control Points
the clock angle of the measurements. For example, two images with camera
Average number of
azimuth angles of 89° and 91° would get credit for two buckets in the binned
Image pixel scale range
Number of images
control points per image
approach, while our algorithm considers those measurements as nearly thescale >1 km
184
34
same and would be scored significantly lower. Similarly, two measurementsthat came at 1° and 89° would fall in the same bucket in the binning system
500 m < scale <1 km
296
173
but would be strongly rewarded as diverse by our coverage approach. In addi200 m < scale <500 m
401
153
tion to the sinusoid, we can also apply a square wave to the coverage arc.
100 m < scale <200 m
136
113This is used only by the nadir measurement, which we intend to be binary.
50 m < scale <100 m
69
69
If there is a single image in the scorable range, the scalar score for the nadir
scale <50 m
25
136
metric will be 1.
PARK ET AL.

8 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

There is a correlation between the angle measurements used as metrics inthe analysis. These correlations are implemented in the algorithm as basic
filtering. For example, it would be overly optimistic to score all measures of
RMS of residuals Total number
incident angle in each image set if we knew that many of those images had(m)of observations
emission angles so high that they would be disqualified from use. For all
439
41,676
metrics except nadir, we pre-filter the measurements by the following ranges:
300
97,295incidence must be between 10° and 80°, emission must be between 0° and
338
138,97180°, and phase must be between 10° and 90°. For the nadir criteria (which is ametric of emission angle), we filter incidence between 20° and 70° and phase293
204,503
below 90°. These filter bounds are wider than any of the scorable ranges
noted in Table 4 above, which is deliberate. Even though we do not score animage with an emission angle of 65° toward the emission metric, we would allow this picture to count toward theincidence metric if that angle was favorable. This is because, in practice, we find images to be usable outside of
the scorable range from Table 4, even if they do not contribute significantly to that metric.

Table 3
Statistics for the Global Topography Model
Map scaleControl point height 1-σ
uncertainty (m)

1 km

286

500 m

297

Combined295

All maps

252

Figure 7 is a graphical representation of the scoring system for a real SPC landmark on the surface of Enceladus.
This landmark is located at approximately (3.3°S, 1.3°E) and includes 58 images that contribute to the topography at 500 m per post. Each box in the figure represents a different angle metric and shows how that metric
was scored in this case. The titles indicate the scorable range, measurements marked with magenta crosses areoutside of the scorable range and contribute zero to the metric, while measurements in green crosses contributeto the score at that location according to the coverage curve assigned to that metric. Inside the scorable range, thecolor of the arc indicates the quality of the coverage in that area due to the measurements. Green regions are wellcovered, but this coverage fades off to yellow and then red as the sinusoid tails off, as described above.
We can examine some of the subplots of Figure 7 to see some intuitive behavior from this analysis. The patch of
surface we are examining is characterized by very high illumination angles, such that most of them fall outside
the scorable range (magenta crosses above 70°). Thus, despite having 58 images, only a handful contributed tocovering the desired range of incident angles, resulting in a score of 0.54. The emission angle, on the other hand,is thoroughly covered except for the very top of the range, resulting in a score of 0.96. The nadir image is anemission angle metric, but the incident angle filter leaves us with very few valid measurements and none below
10°, resulting in a score of zero. We can see that the spacecraft azimuth has very good coverage around the entire
circle, resulting in a score of 0.86, while the solar azimuth leaves two quadrants nearly uncovered and gets a scoreof 0.63. The result is an unweighted average of all six scores equal to 0.65, from which we conclude that this data
set is only moderately well-suited for SPC. It entirely lacks a nadir image and has marginal incidence and solar
azimuth angle coverage.
The next step is to apply this metric to the entire surface (see Figure 8). As always, there are several caveats to thisprocess. We could use every possible image available from the image server, but this ignores the fact that manyof these images are corrupt, badly exposed, or have erroneous a priori pose information. Instead, for the plot in
Figure 8, we have chosen to use a curated list of images that were successfully used in the process of creating theshape model. Thus, this analysis answers the question: given the set of images found to suitable for processing
Table 4
Qualitative Criteria for Evaluating an SPC Imaging Data Set
Minimum case

Optimal case

Criteria

Range (°)

Count

Bins (°)

Count

Bins (°)

Incidence

10–70

2

10–40, 40–70

5

10–20, 20–30, 30–40, 40–50, 50–60

Emission

20–60

2

20–40, 40–60

4

20–30, 30–40, 40–50, 50–60

Nadir (emission)

0–10

1

0–10

1

0–10

Phase

10–80

2

10–80

2

10–80

Spacecraft azimuth

0–360

3

0–120, 120–240, 240–360

4

0–90, 90–180, 180–270, 270–360

Solar azimuth

0–360

2

2 images, 150+ degrees apart4

0–90, 90–180, 180–270, 270–360

Note. Also see Jaumann et al. (2022) for similar criteria for stereophotoclinometry.

PARK ET AL.

9 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 6. Illustration of applying coverage from multiple images to a hypothetical incident angle metric. In the left figure,“1” represents the maximum-allowed score. Also, “Meas” illustrates the qualitative contribution of an individual picturethat can observe this patch of surface. Using the incidence angle as an example, each picture that can view this surface patch
would contribute its own incidence angle to the overall incidence angle metric.

at all, how effective should this data be for producing an SPC model? A total of 916 images contribute to thisanalysis, filtered by their suitability for modeling as just mentioned and also their mean resolution, in the rangeof 150 m to 1.5 km per pixel as is reasonable for use in mapping at 500 m.
With scores generally >0.5, most of the data set used to produce 500 m maps at Enceladus is well-suited for SPC(Figure 8). The dark, splotchy green regions come from the binary nadir metric—a large majority of the surface
does not have nadir image coverage at the required resolution.
2.5. Properties of the New Shape and Topography
The SPC shape model is computed relative to the center of mass. In this coordinate system, the center of figure islocated at (0.3895, 0.7018, 0.3187) km, with an uncertainty of 1 km. Table 5 shows the parameters for the bestfit ellipsoid and best-fit spheroid for the new SPC-based Enceladus shape model that could be used for interior

Figure 7. Sample coverage metrics for a mid-latitude surface patch on Enceladus. The resulting metric scores are Incident (0.54), emission (0.96), phase (0.98), solar
azimuth (0.63), spacecraft azimuth (0.86), and nadir (0.0). The titles indicate the scorable range, measurements marked with magenta crosses are outside of the scorable
range and contribute zero to the metric, while measurements in green crosses contribute to the score at that location according to the coverage curve assigned to thatmetric. Inside the scorable range, the color of the arc indicates the quality of the coverage in that area due to the measurements. Green regions are well covered, but thiscoverage fades off to yellow and then red as the sinusoid tails off.

PARK ET AL.

10 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 8. Overall stereophotoclinometry merit scores across the surface of Enceladus were obtained using images with a surface resolution between 150 and 1,500 m.

analysis, reference shape, astrometric data reduction, etc. The recommended best-fit ellipsoid is (256.14, 251.16,
248.68) km.
Figure 9 shows the cylindrical projection of the SPC topography for ±75° latitude (Figure 9a), stereo projection of
the Northern hemisphere topography for 90°–0° latitude (Figure 9b), and stereo projection of the Southern hemisphere topography for −90°–0° latitude (Figure 9c). The topography is computed relative to the (256.14, 251.16,
248.68)-km best-fit ellipsoid without adjusting for the difference in the center-of-mass and the center-of-figure.
The topography height relative to the best-fit ellipsoid ranges between −3.6 to 3.4 km. We note that it is helpful
to use the best-fit ellipsoid as the reference for computing the topography for visualization purpose and likely
for astrometric data reduction purpose; in general, the geoid (or the hydrostatic shape) is recommended for theinterior analysis.
Two regions merit further discussion: the mountain centered at (18°S, −68°E) and the depression centered at(38°S, 98°E). As can be seen in Figure 9, the former appears to be the highest area relative to the best fit ellipsoidat ∼3.6 km maximum elevation. The latter appears to be the lowest lying region, at ∼3.4 km below the mean
ellipsoid height. Such large height differences had not been reported in previous studies (Bland et al., 2018, 2020;
Tajeddine et al., 2017). A depression very near the one we report here appears in Tajeddine et al. (2017); their S2basin is in approximately the same location but in their work it does not appear to be deeper than 0.6 km. Similar
topographic trends and features are present when the topography is computed relative to the geoid (which is close
to the hydrostatic shape, see Section 5.1 and Figure 6 of Durante et al. (2019)). Given that these larger height
excursions were not previously reported, we performed additional analyses
for these areas to ensure the validity of the results.
Table 5Best-Fit Ellipsoid and Spheroid ParametersCases

ae (km)be (km)ce (km)Mean radius (km)Best-fit ellipsoid

256.14251.16248.68251.99Best-fit spheroid

253.63–248.67251.98Note. Both the ellipsoid and spheroid estimates are provided since someastrometry-related applications might prefer one model over the other. For
the best-fit spheroid ae = be.

PARK ET AL.

We undertook three different verification steps for these regions, eachmethod starting from a state of progressively reduced a priori information
content compared to the previous model. The first round of verification was
to manually inspect and reconstruct all maps in those areas starting with zeroa priori slope and albedo values in each map. The second round included thefirst one but also started with a smooth triaxial ellipsoid as the a priori topography, that is, all prior topography information was erased in those areas.
The third round started from a clean slate; image position and pointing were

11 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 9. Cylindrical projection of the SPC topography for ±75° latitude (a), stereo projection of the Northern hemisphere topography for 90°–0° latitude (b), andstereo projection of the Southern hemisphere topography for −90°–0° latitude (c). The horizontal map resolution is 500 m and the topography is computed relative tothe (256.14, 251.16, 248.68) km best-fit ellipsoid without adjusting for the difference in the center of mass and the center of figure. The topography height ranges from
−3.6 to 3.4 km. In the cylindrical projection, the vertical lines represent the longitude lines for every 30° increment and the middle vertical line represents the 180°E
longitude line.

reset to their a priori values, all topography maps were deleted, and the topography process was performed forthese specific areas anew starting with a smooth ellipsoid as an a priori surface. In all these tests, we were able torecreate both of these features. We have also verified that these areas exhibit similar features for the topographymodels created under the different reflectance models. We therefore believe that these features are real to within
the overall height uncertainty of the model. These topographic features are not evident in individual images but
only in their topographic reconstruction. For example, Figure 10 shows a series of 8 images projected onto a mapthat spans the central area of the lower elevation region at (38°S, 98°E). It is not obvious from the image template(top row) that there is a general lower elevation region there. The rendered height data (bottom row) includes
the depression in the topography, but again this depression is not visually apparent. A hint of such a depression
is seen in Figure 11a, which shows a rendering of the topography in that area as seen from a 60° emission angle
from the viewing geometry of image N1500059482. A synthesized map, showing the highest elevation point at(18°S, −68°E), is shown in Figure 11b. This map that spans a 128 by 128 km square was constructed from all
PARK ET AL.

12 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 10. The top row shows a sample of image templates projected onto a map in the low elevation area centered at (38°S, 98°E). The bottom row displays arendering of the reconstructed map, viewed in the same camera and sun geometry as the corresponding image in the top row.

1 km and 500 m individual maps in that area and is shown from an arbitrary vantage point at an emission angleof 45°.
Figure 12 shows the cylindrical projection of the surface slope for ±75° latitude (Figure 12a), stereo projection of
the Northern hemisphere slope for 90°–0° latitude (Figure 12b), and stereo projection of the Southern hemisphereslope for −90°–0° latitude (Figure 12c). The slope is defined as the angle between the surface normal vector andthe (256.14, 251.16, 248.68) km best-fit ellipsoid surface normal vector. Figure 13 shows a histogram of thesurface slope shown in Figure 12. The area-weighted average slope is ∼4°.
Figure 12 shows that craters exhibit large topographic gradients and are principally scattered across the Northern
Plains region as well as Sub- and Anti-Saturnian Hemispheres. Several craters in the North Plains also exhibit
central mounds (e.g., at 60°N, 150°E) suggesting significant viscous topographic relaxation (and high regional
heat flow, see Bland et al., 2012; Passey, 1983) over time. Figure 12 also highlights several curvilinear ridges andtroughs (i.e., tectonic boundaries) surrounding leading hemisphere, trailing hemisphere, and south polar terrains(Crow-Willard & Pappalardo, 2015); ∼3 km wide L-shaped ridges within the trailing hemisphere terrain (i.e.,
Corfu and Ebony Dorsa centered at 10°N, 75°E) (Crow-Willard & Pappalardo, 2015) and near the Anti-Saturnian
point (i.e., 15°N, −170°E); and three distinct “chasma” extending from the south polar tectonic boundary into thetrailing hemisphere terrain. The south polar terrain is relatively smooth, which is consistent with recent resurfacing and a young (∼1 Myr) overall age.
2.6. Comparison With Previous Enceladus Shape Models
The global topography of Enceladus has been reported in previous studies, and the methods used generally fall
into two categories. One method uses limb profiles (Nimmo et al., 2011; Thomas et al., 2016) or a combinationof limb profiles and control points, selected manually or automatically (Tajeddine et al., 2017). The other isby constructing dense photogrammetric & stereo-photogrammetric control networks (Bland et al., 2018, 2020;

Figure 11. (a) Rendered topography at the low elevation area centered at (38°S, 98°E). The topographic map is at 0.5 km
scale, viewed at 60° emission angle at a geometry similar to that of image N1500059482. (b) Rendered topography of a
128 km by 128 km map at 0.5 km scale, centered at (18°S, −68°E) showing the centrally elevated region.

PARK ET AL.

13 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 12. (a) Cylindrical projection of the slope for ±75° latitude, (b) stereo projection of the Northern hemisphere slope for 90°–0° latitude, and (c) stereo projectionof the Southern hemisphere slope for −90°–0° latitude. The slope is defined as the angle between the surface normal vector and the (256.14, 251.16, 248.68) km best-fit
ellipsoid surface normal vector.

Schenk & McKinnon, 2024). One common characteristic of these photogrammetric-based models is their reliance on similar methodologies, often using identical software packages for point matching and bundle adjustment. Our study presents the first global Enceladus topography model produced using an SPC technique using an
independent software package. Moreover, photo-clinometric techniques in general provide higher spatial resolution than SPG techniques. Our final shape model is based on 500 m scale maps with a total of approximately 55
million distinct elevation points, which is almost a factor of 2 more than the 30 million elevation points presented
in Schenk and McKinnon (2024). In addition to the global shape model, we present the highest possible regional
topography that is allowed by the imaging data. The result is high resolution regional and local topography mostly
in the south pole and the Tiger Stripes region at 200, 100, 50, and 25 m scale that have not been produced before.
Comparison of the initial conditions and bundle adjustment process shows significant differences between our
approach and previous studies. For example, Bland et al. (2020) does not solve for camera position. Their a priori
camera pointing constraint is 0.5°, which if assumed to be a 1-σ value is much larger than ours. On the other
hand, their control point constraint of 1 km in latitude, longitude, and radius is much smaller than our assumption,
which is 20 and 10 km, respectively, for the 1 and 0.5 km topography maps. Similarly, Bland et al. (2018) does not
PARK ET AL.

14 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 13. An area-weighted histogram of the surface slope is shown in Figure 12. The slope is defined as the angle between
the surface normal vector and the (256.14, 251.16, 248.68) km best-fit ellipsoid surface normal vector. The area-weighted
average slope is ∼4°.

solve for camera position and assume a 2° camera a priori pointing error and only a 150 m control point a priori
radius error. Schenk and McKinnon (2024) discuss a 0.2 and 0.3 km 2-σ a priori control point radii relative to thea priori ellipsoid but there's no mention of the camera position and pointing constraints. The biggest difference
between the bundle adjustment of previous studies and ours appears to be their constraint of the control points,
which in our opinion is very restrictive and can explain some of the larger height departures we see in a few
locations compare to previous studies.
A detailed comparison of global topography results between different models is not trivial. A qualitative examination of our Figure 9a, with the results of Tajeddine et al. (2017) Figure 1a, Bland et al. (2020) Figure 6,
and Schenk and McKinnon (2024) Figure 2c, show both similarities and differences. Heights are relative tosimilar triaxial ellipsoids, but quantitative inferences are hard to make given the difference in scales; Tajeddine
et al. (2017) has a scale ranging from −1 km to +1 km, Schenk and McKinnon (2024) −1.5 km to +1.5 km andours from −3.5 km to +3.5 km. Examples of similarities include the slight south pole depression of 300–500 m
that has been observed in all past studies, the depression centered at (30°N, 30°E), and the large depression
bounded by −150° to −70°E and 30°S to 30°N (see Figure 9a). However, in that last region, our model displays
higher short wavelength height fluctuations than in previous studies. The largest height deviations from paststudies appear in the deepest depression centered at (38°S, 98°E) and the mountain centered at (18°S, −68°E)
discussed above. Without a more detailed examination of all shape models, our current assumption for theheight differences among these models is that they are mostly due to a combination of different a priori height
constraints, different topography methodologies and possibly different degrees of resolution in the elevation
points produced.

3. Rotational Model
To orient Enceladus, we developed expressions for the rotational elements based on the SAT441 satellite ephemerides (Jacobson, 2022). The Enceladus pole is aligned with the orbit normal direction and the rotation period
matches the mean orbital period. Periodic terms in the prime meridian expression occur because it is measured
from the node of the satellite's equator on the ICRF reference plane. Where appropriate, we added terms toaccount for long period mean motion librations. We also included a term in the prime meridian expression toaccount for the observed forced libration (Nadezhdina et al., 2016; Thomas et al., 2016) and ignore the potential
biasing impact of deformation over the orbital timescale (Berne et al., 2023; Van Hoolst et al., 2016).
The inertial orientation of Enceladus in the International Celestial Reference Frame is modeled𝐴𝐴with (𝛼𝛼𝛼 𝛼𝛼𝛼 𝛼𝛼 ),where𝐴𝐴𝐴𝐴 represents the spin-pole declination,
𝐴𝐴 and 𝐴𝐴 represents the𝐴𝐴
𝐴𝐴 represents the spin-pole right ascension,rotation around the spin-pole axis. Specifically, the time series
𝐴𝐴 of (𝛼𝛼𝛼 𝛼𝛼𝛼 𝛼𝛼 ) are defined as:
PARK ET AL.

15 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

]
𝑞𝑞
∑
(
)
̇
𝛼𝛼(𝑡𝑡) = 𝛼𝛼0 + 𝛼𝛼̇ 0 𝑡𝑡 +𝐴𝐴𝑗𝑗 sin Ω𝑗𝑗 𝑡𝑡 + Ω𝑗𝑗
(1)
[

𝑗𝑗=1

]
𝑞𝑞
∑
(
)
̇
̇
𝛿𝛿(𝑡𝑡) = 𝛿𝛿0 + 𝛿𝛿0 𝑡𝑡 +𝐵𝐵𝑗𝑗 cos Ω𝑗𝑗 𝑡𝑡 + Ω𝑗𝑗
(2)
[

𝑗𝑗=1

]
[ 𝑞𝑞
∑
(
)
̇
̇
𝑊𝑊 (𝑡𝑡) = 𝑊𝑊0 + 𝑊𝑊0 𝑑𝑑 +𝐶𝐶𝑗𝑗 sin Ω𝑗𝑗 𝑡𝑡 + Ω𝑗𝑗
(3)
𝑗𝑗=1

(
)
terms, 𝛼𝛼̇ 0 , 𝛿𝛿̇ 0 , 𝑊𝑊̇ 0 are rate terms, q is the number
where (α0,δ0,W0) are constant𝐴𝐴of terms in the nutation series, (Aj,Bj,Cj) are nutation amplitudes,𝐴𝐴Ω̇ 𝑗𝑗 is the nutation frequency, Ωj is the nutation phase, t is the time in seconds since J2000,
and d is the time in days since J2000 (Archinal et al., 2018). J2000 is defined as
01-JAN-2000 12:00:00. The rate 𝐴𝐴terms 𝐴𝐴𝐴𝐴
𝐴 0 and 𝛿𝛿̇ 0 represent spin-pole's precession and the terms in the summation represent the nutation series.

The coefficients in Equations 1–3 were obtained from a least squares fit of
the SAT441 ephemeris over the 20-year time frame 2000–2020. The constantangle, W0, was chosen such that it is consistent with the International Astronomical Union's definition of Enceladus's prime meridian, which is defined by the
Figure 14. The center of the Salih crater, located at −5° East longitude in thecenter of a small crater called Salih located at −5° longitude (see Figure 14).
planetocentric coordinate system, defines the prime meridian of Enceladus.
The frequency and phase of the forced libration amplitude were set based on
the average period of Enceladus over the duration of the Cassini mission. The
amplitude of the forced libration was estimated by fitting the landmark observations, that is, landmark body-fixed
vectors and landmark image sample, line locations, and iterated with the SPC topography development until the solution converged. This process also placed the center of the Salih crater at −5° East longitude. The estimated parameters
were the forced libration amplitude, camera pointing, camera position, and landmark vectors. The landmarks usedwere those of the 1 km and 500 m maps because they are uniformly distributed on the surface in both longitude andlatitude. The baseline solution included 1,440 landmarks and images from 2005 to 2016, with a total of 95,023 landmark observations. The updated rotational model was used to rebuild the Enceladus topography. The final rotationalmodel is self-consistent with the topography solution in that it minimizes the posteriori residuals of the topographyreconstruction. Various parametric studies that included reduced data sets, such as using landmarks spanning different
latitude bands or using a subset of the images, also confirmed the baseline solution.

Table 6
Estimated Forced Libration Amplitudes Based on Different Latitudinal
Cutoff for Landmarks and a Subset of Data

Table 6 shows estimated forced libration amplitude based on different
latitudinal cutoff for landmarks and a subset of data. The result shows arobust recovery of the forced libration amplitude with the resulting valueof −0.091° ± 0.009° (3-σ). Table 6 also shows two previous forced libration
amplitude estimates from analyzing the Cassini imaging data (Nadezhdina
et al., 2016; Thomas et al., 2016).

Forcedlibration
amplitude (°)

Libration amplitudeuncertainty (°) (3-σ)

All data, all landmarks−0.091

0.009

All data, landmarks in ±60° latitudes−0.0930.012

All data, landmarks in ±30° latitudes−0.0920.021+ 0◦ .000686 sin (S02) − 0◦ .000472 sin (S03)

All data before 2013, all landmarks−0.0920.015− 0◦ .000897 sin (S04) + 0◦ .002970 sin (S05)Thomas et al. (2016)−0.12

0.014 (2-σ)

Nadezhdina et al. (2016)−0.1550.014 (σ unspecified)+ 0◦ .001127 sin (S06) + 0◦ .000519 sin (S07)(4)+ 0◦ .000228 sin (S08) + 0◦ .036804 sin (S09)Cases

Note. The result shows a robust recovery of the forced libration amplitudewith the resulting value of −0.091° ± 0.009° (3-σ). The dataset used for the
fourth-row case (i.e., “all data before 2013, all landmarks”) is similar to thedataset considered in Thomas et al. (2016) and Nadezhdina et al. (2016).

PARK ET AL.

The complete expressions for the inertial orientation angles (α,δ,W), for
Enceladus are
𝛼(𝑡) = 40◦ .592915 − 0◦ .0902111773𝑇 + 0◦ .026616 sin (S01)− 0◦ .001107 sin (S10) + 0◦ .073107 sin (S11)− 0◦ .000167 sin (S12) − 0◦ .000376 sin (S15)+ 0◦ .000248 sin (S16) − 0◦ .000137 sin (S17)

16 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

𝛿𝛿(𝑡𝑡) = 83◦ .534180 − 0◦ .0071054901𝑇𝑇 + 0◦ .004398 cos(S01) − 0◦ .000264 cos(S02)−0◦ .000185 cos(S0 − 0◦ .000093 cos(S04) − 0◦ .000068 cos(S05)−0◦ .000236 cos(S06) − 0◦ .000028 cos(S08) + 0◦ .004141 cos(S09)(5)−0◦ .000124 cos(S10) + 0◦ .008229 cos(S11) + 0◦ .000007 cos(S12)−0◦ .000039 cos(S15) + 0◦ .000026 cos(S16) − 0◦ .000016 cos(S17)
𝑊 (𝑡) = 8◦ .325383 + 262.7318870466𝑑 − 0◦ .026447 sin(S01) − 0◦ .000682 sin(S02)+ 0◦ .000469 sin(S03) − 0◦ .005118 sin(S04) + 0◦ .036955 sin(S05)− 0◦ .013111 sin(S06) + 0◦ .014206 sin(S07) − 0◦ .006687 sin(S08)− 0◦ .036404 sin(S09) + 0◦ .001082 sin(S10) − 0◦ .072604 sin(S11)(6)− 0◦ .266358 sin(S12) − 0◦ .188429 sin(S13) − 0◦ .004710 sin(S14)+ 0◦ .000337 sin(S15) − 0◦ .000183 sin(S16) − 0◦ .001724 sin(S17)− 0◦ .091295 sin(S18)with the corresponding arguments for nutation sinusoidal terms:
S01 = 335◦ .844470 + 51◦ .7682239𝑇
S02 = 355◦ .351814 + 101◦ .6467750𝑇
S03 = 9◦ .369346 + 1004◦ .8728024𝑇
S04 = 129◦ .755966 + 1223◦ .2050690𝑇
S05 = 219◦ .755966 + 1223◦ .2050690𝑇
S06 = 159◦ .835559 + 2445◦ .2902118𝑇
S07 = 249◦ .835559 + 2445◦ .2902118𝑇
S08 = 117◦ .392885 + 3667◦ .0200695𝑇



S09 = 280◦ .169482 + 7226◦ .3782354𝑇
S10 = 6◦ .997174 + 36506◦ .5422127𝑇
S11 = 196◦ .673251 + 15227◦ .2035409𝑇
S12 = 253◦ .848856 + 3258◦ .6617087𝑇
S13 = 136◦ .859155 + 9266◦ .8742489𝑇
S14 = 144◦ .630256 + 12292◦ .3910895𝑇
S15 = 9◦ .821866 + 16090◦ .5831593𝑇
S16 = 226◦ .334387 + 17383◦ .5986496𝑇
S17 = 93◦ .360491 + 18531◦ .0794323𝑇
S18 = 10◦ .9818392 + 9583937◦ .8056363𝑇

where T is centuries from J2000 and d is days from J2000 (Archinal et al., 2018).

4. Gravity Field
Among all flybys of Enceladus in the Cassini tour, the estimate of the gravity field uses only three, the E09, E12,and E19 Enceladus encounter data; only these encounters had continuous tracking throughout the close flyby
period. Figure 15 shows the ground track of these encounters for ±2 hr of the closest approach and Table 7 shows
the corresponding geometrical information.
The external gravitational potential of Enceladus is modeled using a spherical harmonic expansion:∞𝑛𝑛 ()
𝐺𝐺𝐺𝐺 ∑ ∑ 𝑅𝑅 𝑛𝑛
𝑈𝑈 (𝑟𝑟𝑟 𝑟𝑟𝑟 𝑟𝑟) =𝑃𝑃𝑛𝑛𝑛𝑛 (sin 𝜙𝜙)[𝐶𝐶𝑛𝑛𝑛𝑛 cos(𝑚𝑚𝑚𝑚) + 𝑆𝑆𝑛𝑛𝑛𝑛 sin(𝑚𝑚𝑚𝑚)],(7)𝑟𝑟 𝑛𝑛=0 𝑚𝑚=0 𝑟𝑟

PARK ET AL.

17 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 15. Cassini ground-tracks at Enceladus for encounters E09, E12, and E19. Solid circles represent the closest approach
points.

Here, G is the universal gravitational constant, M is the mass of Enceladus, R is the reference radius of Enceladus(256.6 km), n is the degree, m is the order, Pnm are the associated Legendre functions, and Cnm and Snm are theun-normalized gravity spherical harmonic coefficients (the corresponding un-normalized zonal harmonics are
Jn = −Cnm).
At each body-fixed point defined by the latitude (ϕ), longitude (λ), and radius (r), the gravitational accelerationof an external point source (e.g., Cassini spacecraft) is given by the gradient of this potential. Since the spacecraft
dynamics are influenced by Enceladus's gravity field, which is tied to its body-fixed coordinate frame, accurately
tracking the spacecraft motion allows the recovery of gravity field and rotational parameters.
Our study is an extension of the SAT441 study, where the SAT441 results were obtained within the context of the
global Saturnian satellite ephemeris development and Cassini tour reconstruction (Jacobson, 2022). The SAT441
study included the determination of the gravity fields of Saturn, Titan, and several icy satellites, and Saturn'spole orientation and precession. The E9, E12, and E19 orbits provided the data for the Enceladus gravity field
estimate. However, the data arcs covered the full Cassini orbits, and the data sets included range and spacecraft
imaging. Moreover, the SAT441 work estimated satellite ephemeris parameters as well as the spacecraft maneuvers and attitude control and acceleration due to Cassini's radioisotope thermoelectric generator. The Ka-band
Doppler was processed with a 30 s sample rate. Because E9 and E19 were near solar opposition, the Doppler datanoise was dominated by fluctuations in the troposphere and the encounter passes were weighted at 110% of the
RMS of the residuals. For E12, it was assumed that solar plasma was the source of the Doppler noise and the data
weights were derived from the data whitening algorithm.
For this work, we re-processed the E9, E12, and E19 arcs from SAT441 with a revised Enceladus orientation
model from Section 3. We added four more Doppler points to the E9 arc and three to the E19 arc and slightly
adjusted the times of the impulse burn modeled to represent the drag force due to the plume. We estimated only
Table 7
Geometric Information for Enceladus Gravity Encounters
Encounter

Time (TDB)Dist. (km)Lat. (°)

Lon. (°)

Alt.a (km)Vel. (km/s)

E09

28-APR-2010 00:11:23.031

349.1

−89.3

−158.1

97.1

6.5

E12

30-NOV-2010 11:55:05.596

294.9

61.8

−58.8

42.9

6.3

E19

02-MAY-2012 09:32:34.982

321.7

−71.9

63.0

69.7

7.5

Relative to the 251.99-km best-fit mean radius (see Table 5).

a

PARK ET AL.

18 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Table 8
Un-Normalized Dimensionless Enceladus Gravitational Harmonics (×10 6), See Equation 7
DegreeOrder

Iess et al. (2014) a

Case 1 (n = 2)

n

m

Cnm

Snm

Cnm

2

0

−5334.0 ± 34.321

9.0 ± 11.4

221520.9 ± 15.3

22.2 ± 7.3

3

0112.1 ± 22.3

–

Snm

–

−5520.55 ± 36.98

39.1 ± 22.0

−17.93 ± 13.65

9.81 ± 9.12

1517.68 ± 14.68−246.66 ± 9.31

DV (E09) mm/s

0.247 ± 0.053

0.39 ± 0.7

DV (E19) mm/s

0.256 ± 0.049

2.93 ± 0.3

Case 2b (n = 2 + J3)

Case 3 (n = 2 + J3 + J4)

Degree

Order

n

m

Cnm

2

0

−5477.45 ± 36.99

2

1

7.86 ± 13.66

2

21517.90 ± 14.70

30177.82 ± 33.42

4

0

–

Snm

Cnm

Snm

–

−5582.06 ± 51.95–

7.60 ± 9.19

−20.92 ± 15.35−10.33 ± 10.93−275.31 ± 10.87

1526.92 ± 14.95−281.63 ± 11.49–

165.59 ± 34.27––

88.71 ± 47.64–

DV (E09) mm/s

0.48 ± 0.7

0.44 ± 0.7

DV (E19) mm/s

2.89 ± 0.3

2.84 ± 0.3

Gravity harmonics from shape assuming
constant-density interior

Topography harmonics from shapec

Degree

Order

n

m

Cnm

Snm

Cnm

Snm

2

0

−7590.80–−12892.33–

21

189.85467.24323.84816.822

2

1925.82−373.393218.79−622.843

0910.13–2284.46–40

283.03–566.89–Note. The reference radius used was 256.6 km and the quoted uncertainties are the formal 1-σ uncertainties.
a
Iess et al. (2014) used 254.2 km as the reference radius. The values shown in this table are scaled using the 256.6 km reference radius. bCase 2 is our suggested solution. cThe
∞ 𝑛𝑛∑∑topography harmonics from shape can be applied to the following equation for computing the radial distance:
𝑃𝑃𝑛𝑛𝑛𝑛 (sin 𝜙𝜙)[𝐶𝐶𝑛𝑛𝑛𝑛 cos(𝑚𝑚𝑚𝑚) + 𝑆𝑆𝑛𝑛𝑛𝑛 sin(𝑚𝑚𝑚𝑚)].
𝐴𝐴
𝐴𝐴(𝜆𝜆𝜆 𝜆𝜆) = 𝑅𝑅
𝑛𝑛=0 𝑚𝑚=0

the Enceladus gravity field and the spacecraft parameters (states, maneuvers, etc.) for each arc. The SAT441
satellite ephemeris was held fixed as was the Enceladus GM.
Table 8 shows the estimated gravity field results (i.e., Cases 1, 2, and 3) along with the solution obtained by the
Cassini Radio Science Subsystem (RSS) (Iess et al., 2014). Table 8 also includes an impulsive ∆V at closest
approach to Enceladus for the E09 and E19 encounters to account for drag caused by flying through the plume.
Case 1 estimates the full quadrupole field (J2, C21, S21, C22, and S22). Case 2 estimates the full quadrupole field
and J3. Case 3 estimates the full quadrupole field, J3, and J4.
The full quadrupole field (i.e., Case 1) provides an adequate fit to the observations with a few small outlying
points near the closest approach, primarily due to the drag force affecting the spacecraft dynamics when going
through the plume-rich region. Including J3 completely flattens the residuals as shown in Figure 16, and this is
our recommended gravity field for Enceladus. Adding J4 makes an insignificant contribution to the overall result.
Figure 17 is the repeat of Case 2 but without estimating an impulsive burn for E09 and E19, showing an obvious
discontinuity in the residuals at the time of the plume crossing.

PARK ET AL.

19 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 16. Two-way (black) and three-way (green) Doppler residuals of Case 2 (with the full quadrupole field and J3): (a)E09 encounter, (b) E12 encounter, and (c) E19 encounter.

5. Enceladus Interior Modeling5.1. Topography Relative to the Geoid
By combining our shape and gravity models, we can illustrate the topography of Enceladus relative to the geoid.
The geoid is an equipotential surface defined by the gravitational field plus the effects of tides and rotation, andis given by

PARK ET AL.

20 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 17. Doppler residuals of Case 2 (with the full quadrupole field, J3, and J4) without an impulsive burn model torepresent the drag perturbation due to the plume: (a) E09 encounter, and (b) E19 encounter. There is an obvious discontinuity
in the residuals at the time of the plume crossing.
∞𝑛𝑛∑∑
𝑁𝑁(𝜆𝜆𝜆 𝜆𝜆) = 𝑅𝑅
𝑃𝑃𝑛𝑛𝑛𝑛 (sin 𝜙𝜙)(𝐶𝐶𝑛𝑛𝑛𝑛 cos 𝑚𝑚𝑚𝑚 + 𝑆𝑆𝑛𝑛𝑛𝑛 sin 𝑚𝑚𝑚𝑚)
𝑛𝑛=0 𝑚𝑚=0(8)𝜔𝜔2 𝑅𝑅2
5𝜔𝜔2 𝑅𝑅2
𝑃𝑃20 (sin 𝜙𝜙) +𝑃𝑃22 (sin 𝜙𝜙) cos 2𝜆𝜆
−
6𝑔𝑔
4𝑔𝑔

where the last two terms in Equation 8 are required to account for the combined tidal and rotational potentials andwhere ω and g are Enceladus's rotation rate and surface gravity, respectively.
The topography relative to the geoid depicted in Figure 18 suggests that Enceladus's overall shape issignificantly non-hydrostatic, with long-wavelength topography on the order of ±1,500 m. The range in
non-hydrostatic topography (i.e., roughness) is also about 3 times greater (i.e., −3–2 km range) than thatseen in other shape models (e.g., Nimmo et al., 2011; Tajeddine et al., 2017, ±1 km range). Most of the additional range in topographic variation arises from the presence of a short wavelength structure in the Southern
Hemisphere.
The ∼3 km high mountain and a ∼3 km deep basin centered at (−68°E, 18°S) and (98°E, 38°S) shown in
Figure 18 mirror similar features visible in maps of topography relative to the best-fit ellipsoid (see Figure 9a).
Both features are probably supported elastically within the ice shell to some extent (full Airy isostasy of the basinwould require a mean ice shell thickness of ∼40 km which is inconsistent with observational constraints on theinterior structure, see Section 5.3). Moreover, though the nearly antipodal position of either surface feature may
suggest an impact-driven formation mechanism, the ∼100 km diameter half-dome shape of the basin is difficult

PARK ET AL.

21 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 18. Cylindrical projection showing the topography of Enceladus expanded to degree and order 20, with respect to its geoid.

to reconcile with predictions of a relatively flat crater center if an ocean is present ∼30 km below the surface(Bray et al., 2014). We note that the basin is surrounded by a likely compressional portion of the South Polar
Terrain Boundary to the South and “chasma” (i.e., extensional fractures) to the East and West, indicating a possible tectonic origin for the feature (Yin & Pappalardo, 2015).
5.2. Geology
Correlations between the topography and geology of Enceladus can provide insight into the formation of different
terrains. As mentioned above, the prominent basins centered at (18°S, −68°E) and (38°S, 98°E) are surroundedby terrains interpreted to be compressional on one side (e.g., Yin & Pappalardo, 2015) and extensional on another(e.g., Crow-Willard & Pappalardo, 2015) Further analysis of the terrains surrounding these basins may provide
further support (or not) for a tectonic origin. Here we provide some geologic context and potential correlations
for the topographic rises observed in Figure 18.
Regions with higher crater density (Kinczyk et al., 2017; Patterson et al., 2018) appear to be associated withmoderate and broad topographic highs centered around longitudes 160° and 340°E (Figure 18). Particularly, theregion of highest crater density around (45°N, 90°E) (Kinczyk et al., 2017) is correlated with one of the highest
regions of topography. The origin of this apparent association is unclear. While there are tectonic structures
within the cratered terrains, they have been primarily interpreted as either strike-slip (Martin, 2016) or extensional (Martin et al., 2017) in origin and are therefore not typically associated with broad topographic highs. It is
possible that there is a significant tectonic history shrouded by the emplacement and presence of impact craters
and plume deposition (Patterson et al., 2018), or that these regions remain relatively topographically high due tothe lack of major tectonic features, more akin to continental lithosphere on Earth.
Another, more prominent, broad topographic rise, centered around longitude 90°E appears strongly correlatedwith the dorsa, tectonic features on the Trailing Hemisphere (Spencer et al., 2009) interpret the dorsa as cryovolcanic in origin whereas (Patthoff et al., 2022) interpret the dorsa as large-scale contractional features, which
could contribute to the high-standing nature of the regional topography. For example, as the dorsa were formed
through contractional processes, the ice shell would locally thicken, potentially creating a broad topographic
high. If this is the case, it would imply that the topographic signal from the formation of the dorsa has not yet
completely relaxed away despite inferences of past episodes of high heat flux in the surrounding region (Bland
et al., 2007; Giese et al., 2008).
The prominent topographic rise centered around (18°S, −68°E) falls just to the north-east of the enigmatic reliccratered blocks (Crow-Willard & Pappalardo, 2015; Leonard et al., 2021). While some features in this area areassociated with compression, there is also a series of relatively young fractures that run directly though this
PARK ET AL.

22 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

region (Leonard et al., 2021). While the prominent mound appears constrained by the relic cratered blocks to thesouth and two prominent troughs to the East and West, there is no apparent bounding feature to the north andapparent lack of change in the surface geology within the feature that would indicate the presence of the topographic rise. The origin of this mound remains unclear, including whether it could be compensated by the ice
shell or requires active formation.
5.3. Interior Analysis
Knowledge of the interior of Enceladus is important for understanding its origin and evolution, the nature of its
present activity, and its prospects for habitability. Estimating the thickness and laterally variable structure of theice shell is particularly important for understanding the energy budget and the nature of active eruptions in theSouth Polar Terrain.
Several groups have developed interior models for Enceladus based on previous estimates of its shape (Nimmo
et al., 2011; Tajeddine et al., 2017) and gravity field (Iess et al., 2014). Some of these models are basedon shape and gravity alone (Beuthe et al., 2016; Iess et al., 2014; McKinnon, 2015) while others additionally incorporate the estimated libration amplitude (Cadek et al., 2016; Hemingway & Mittal, 2019; Thomas
et al., 2016; Van Hoolst et al., 2016). All these models agree that Enceladus consists of a ∼2,300–2,500 kg/m 3
rocky interior surrounded by a ∼50–70 km thick hydrosphere that is partitioned into a global subsurface
liquid water ocean and an icy shell with long-wavelength topography that is largely compensated (Hemingway
et al., 2018).
Such inversion models are inherently non-unique, however, and estimates of the precise thickness and lateralvariability of the ice shell's structure are sensitive not only to the measurements but also to model assumptions. Model assumptions include the range of layer densities that are considered, the way the reference hydrostatic equilibrium figures are calculated, and details about how the compensation mechanism is modeled. For
example, several recent studies have explored a variety of distinct models of Airy-type isostatic compensation. Static isostasy models vary in whether they are aiming to: equate masses in cones of equal solid angle(Iess et al., 2014; McKinnon, 2015); equate pressures along equipotential surfaces at depth (Hemingway &Matsuyama, 2017; Hemingway & Mittal, 2019); equate forces acting at the top and bottom of the lithosphere(Čadek et al., 2017, 2019); or minimize deviatoric stresses within the lithosphere (Beuthe et al., 2016). Dynamic
models, which consider viscous or viscoelastic responses to applied loads have also been explored recently by
several authors (Beuthe, 2021a, 2021b; Čadek et al., 2017, 2019; Ermakov et al., 2017). Which of these models is
the most appropriate choice for studying the interior of Enceladus remains an open question but, for reviews, see(Beuthe, 2021a, 2021b; Beuthe et al., 2016; Cadek et al., 2019; Hemingway & Matsuyama, 2017; Hemingway
et al., 2018).
However, the goal of our analysis here is not to explore the effects of altering various model assumptions. Instead,this study focuses on how the conclusions about the interior structure of Enceladus are affected by our new
shape, gravity, and libration estimates. To that end, we carry out a direct comparison primarily against one
previous study (Hemingway & Mittal, 2019) to highlight how the new measurements shift the results when allother aspects of the modeling are held constant. While we focus mainly on this direct comparison, we will alsocomment briefly on the effects of altering our isostasy model as an example to show the distinct effects of alteringmodeling assumptions versus updating the input data.
For our interior structure analysis, we follow the methodology of Hemingway and Mittal (2019), which incorporates the (Tricarico, 2014) approach for computing hydrostatic equilibrium figures. This methodology involves
adopting a given shape model, imposing an isostasy condition, calculating the expected gravity field and libration
amplitude over a wide range of possible interior structures, and then identifying the parts of that parameter spacethat best match the observed gravity field and libration amplitude while accounting for uncertainties. Hemingway
and Mittal (2019) showed that the ice shell and ocean densities are not well constrained by shape, gravity, andlibration measurements alone and that conclusions about the ice shell and ocean thicknesses are not strongly
affected by their assumed densities. They also argued that the observations are better accommodated by an Airyrather than a Pratt isostasy compensation model. Accordingly, we do not attempt to explore the parameter space
as comprehensively as Hemingway and Mittal (2019) and instead adopt the assumption of Airy isostasy and fixed
ice shell and ocean densities of 925 and 1,020 kg/m 3, respectively. This reduces our parameter space exploration
PARK ET AL.

23 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Figure 19. Shift in the preferred interior structure due to adopting our new shape, gravity, and libration estimates. (a and e) Parts of the parameter space satisfying
the observed J2 (blue), C22 (red), and J3 (gold) gravity coefficients and the observed libration amplitude (purple). Dark and light shading correspond to 1-σ and 2-σuncertainties, respectively. (b and f) Misfit between model and observations given by Hemingway and Mittal (2019, Equation 21). (c, d, g, and h) Confidence contours(dark = 68%; intermediate = 95%; pale = 99.7%) showing model parameters that best satisfy the observations. Panels (a–d) are based on the previously publishedshape, gravity, and libration estimates (Iess et al., 2014; Tajeddine et al., 2017; Thomas et al., 2016) and correspond to panels (e–h) of Figure 4 by Hemingway and
Mittal (2019). Panels (e–h) follow the same approach but incorporate our new shape, gravity, and libration estimates.

to two dimensions, as illustrated in Figure 19 and facilitates direct comparison with Figure 4 by Hemingway and
Mittal (2019).
Although our results broadly agree with previous work, our new estimates of shape, gravity, and libration
amplitude lead us to prefer a larger ice shell thickness, smaller ocean thicknesses, and a slightly larger andlower density rocky core. Our preferred values are shown in Table 9 alongside those arising from the same
analysis but using the previously published shape, gravity, and libration data (Hemingway & Mittal, 2019;
Iess et al., 2014; Tajeddine et al., 2017; Thomas et al., 2016). Overall, our new results suggest a thicker ice
shell (mean thickness 27–33 km) and a thinner ocean (mean thickness 21–26 km), with a mean core density
in the range of 2,270–2,330 kg/m 3. The shift is due to a combination of factors, but the largest influencecomes from the fact that our new gravity model is more non-hydrostatic than that of Iess et al. (2014),
whereas our new shape model is closer to hydrostatic than that of Tajeddine et al. (2017). These differencesin the recovered gravity and shape suggest a larger compensation depth, and therefore a greater ice shell
Table 9Preferred Interior Structure Properties

PARK ET AL.

Shape, gravity, and libration estimates

Tajeddine et al. (2017), Iess et al. (2014), and Thomas
et al. (2016)

Parameter

Preferred

1−σ range

Preferred

1−σ range

Ice shell thickness (km)

21

19–24

30

27–33

Ocean thickness (km)

37

35–39

23

21–26

This work

Rocky core radius (km)194

192–195

198

197–200

Rocky core density (kg/m 3)

2,370

2,340–2,410

2,290

2,270–2,330

Moment of inertia factor

0.335

0.333–0.336

0.338

0.336–0.339

24 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

thickness. Our smaller libration amplitude likewise favors a larger ice shell thickness because the libration
amplitude is inversely proportional to the ice shell's moment of inertia, and thus its thickness. The effects
on the rocky core are more subtle, but our results suggest a slightly larger and lower density rocky core andtherefore a slightly larger moment of inertia.
Consistent with our main point of comparison (Hemingway & Mittal, 2019), the above analysis adopts the “equalpressures” conception of Airy isostasy (Hemingway & Matsuyama, 2017). As an additional check, we also carriedout the same analysis while instead employing the “equal weights” isostasy model (Cadek et al., 2019). The result
is similar in that adopting the new shape, gravity, and libration data shifts the preferred internal structure towarda larger shell thickness and smaller ocean thickness (Figure S1). A similar shift toward a larger preferred shell
thickness is not surprising because it is largely driven by the smaller libration amplitude, which is unrelated tothe isostatic compensation model. Similarly, as discussed in the previous paragraph, the shift toward a less hydrostatic gravitational field coupled with a more hydrostatic shape naturally leads to a larger compensation depth,
regardless of which isostasy model one adopts. We note, however, that under the “equal weights” isostasy model,
the combination of gravity and libration constraints do not converge quite as well. For example, the J3 constraint
suggests a shell thickness that is somewhat larger than the value suggested by the libration constraint (Figure S1).
Of course, as discussed earlier, all these results are subject to a variety of modeling assumptions, so future work
may be able to resolve such discrepancies. In any case, regardless of these sorts of model assumptions, comparedwith previously published estimates for shape, gravity, and librations, our new estimates suggest a larger shell
thickness.
While some of our results appear to differ only subtly from those of Hemingway and Mittal (2019), the increasein our preferred ice shell thickness from 21 to 30 km is significant. One of the more important implications is that
this larger shell thickness would suggest a lower heat flux. Whereas, assuming a conductive ice shell (Hemingway& Mittal, 2019) estimated a total conductive heat loss of 25–40 GW, our new shape, gravity, and libration estimates lead instead to a total heat loss of 18–28 GW.

6. Conclusions
We presented the results from an analysis of Cassini's radiometric tracking and onboard imaging data duringclose Enceladus encounters. The global shape, gravity field, and libration angle of Enceladus were computed
in a common frame consistent with the International Astronomical Union's definition, where the center of the
Salih crater is located at −5°E longitude. We recovered a quadrupole gravity field with J3 and a forced libration amplitude of 0.091° ± 0.009° (3σ). This libration estimate is lower than the previously estimated values
from Thomas et al. (2016) and Nadezhdina et al. (2016). We believe this difference is likely due to the use ofdifferent Enceladus and spacecraft ephemerides, the span of data arcs, and the techniques for extracting thelibration amplitude. A global topography model was also computed using a stereo-photoclinometry techniquewith a global resolution of 500 m, with some local regions having resolutions ranging 25–100 m. All products
are available through the PDS for interior studies and planning for future missions to Enceladus. Our overall
results suggest a thicker 27–33 km mean ice shell, a thinner ocean thickness of 21–26 km, and a mean coredensity range of 2,270–2,330 kg/m 3.

Data Availability Statement
The data that support the findings of this study are available from the PDS website. The imaging data used for
the SPC shape and libration are available in the Cassini Imaging Science Subsystem archive (Knowles, 2018)via https://doi.org/10.17189/1504135. The radiometric data used for the gravity field are available in the Cassini
RSS archive (Asmar, 2021) via https://doi.org/10.17189/1522471. The derived SPC shape models of Enceladus
in Digital Shape Kernel format and associated files are archived through the NAIF PDS node (Elson et al., 2005)
via https://doi.org/10.17189/1520114; see https://naif.jpl.nasa.gov/pub/naif/pds/data/co-s_j_e_v-spice-6-v1.0/
cosp_1000/data/. The derived shape models are available in the “dsk” folder, and the associated rotational model
kernel is available in the “pck” folder. Alternatively, the derived SPC shape models can also be found at http://
naif.jpl.nasa.gov/pub/naif/CASSINI/misc/enceladus_park2023/.

PARK ET AL.

25 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

Acknowledgments
This research was carried out at the
Jet Propulsion Laboratory, California
Institute of Technology, under a contract
with the National Aeronautics and Space
Administration. The authors would liketo thank two anonymous reviewers fortheir thoughtful reviews. Copyright 2023
California Institute of Technology. U.S.
Government sponsorship acknowledged.

PARK ET AL.

10.1029/2023JE008054

References
Archinal, B. A., Acton, C. H., A'Hearn, M. F., Conrad, A., Consolmagno, G. J., Duxbury, T., et al. (2018). Report of the IAU Working Group on
Cartographic Coordinates and Rotational Elements: 2015. Celestial Mechanics and Dynamical Astronomy, 130(3), 22. https://doi.org/10.1007/
s10569-017-9805-5
Asmar, S. W. (2021). PDS3 Cassini Radio Science User’s Guide. NASA Planetary Data System. https://doi.org/10.17189/1522471
Berne, A., Simons, M., Keane, J. T., & Park, R. S. (2023). Inferring the mean thickness of the outer ice shell of Enceladus from diurnal crustal
deformation. Journal of Geophysical Research-Planets, 128(6), e2022JE007712. https://doi.org/10.1029/2022je007712
Beuthe, M. (2021a). Isostasy with Love – I: Elastic equilibrium. Geophysical Journal International, 225(3), 2157–2193. https://doi.org/10.1093/
gji/ggab073
Beuthe, M. (2021b). Isostasy with Love: II Airy compensation arising from viscoelastic relaxation. Geophysical Journal International, 227(1),693–716. https://doi.org/10.1093/gji/ggab241
Beuthe, M., Rivoldini, A., & Trinh, A. (2016). Enceladus's and Dione's floating ice shells supported by minimum stress isostasy. GeophysicalResearch Letters, 43(19), 10088–10096. https://doi.org/10.1002/2016gl070650
Bland, M. T., Becker, T. L., Edmundson, K. L., Roatsch, T., Archinal, B. A., Takir, D., et al. (2018). A new Enceladus global control network,
image mosaic, and updated pointing kernels from Cassini's 13-year mission. Earth and Space Science, 5(10), 604–621. https://doi.
org/10.1029/2018ea000399
Bland, M. T., Beyer, R. A., & Showman, A. P. (2007). Unstable extension of Enceladus' lithosphere. Icarus, 192(1), 92–105. https://doi.
org/10.1016/j.icarus.2007.06.011
Bland, M. T., Singer, K. N., McKinnon, W. B., & Schenk, P. M. (2012). Enceladus' extreme heat flux as revealed by its relaxed craters. Geophysical Research Letters, 39(17), L17204. https://doi.org/10.1029/2012gl052736
Bland, M. T., Weller, L. A., Mayer, D. P., & Archinal, B. A. (2020). A Global Shape Model for Saturn's Moon Enceladus from a Dense
Photogrammetric Control Network. In ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, V-3-2020(pp. 579–586). Retrieved from https://isprs-annals.copernicus.org/articles/V-3-2020/579/2020/
Bray, V. J., Collins, G. S., Morgan, J. V., Melosh, H. J., & Schenk, P. M. (2014). Hydrocode simulation of Ganymede and Europa cratering trends- How thick is Europa's crust? Icarus, 231, 394–406. https://doi.org/10.1016/j.icarus.2013.12.009
Čadek, O., Běhounková, M., Tobie, G., & Choblet, G. (2017). Viscoelastic relaxation of Enceladus’s ice shell. Icarus, 291, 31–35. https://doi.
org/10.1016/j.icarus.2017.03.011
Cadek, O., Soucek, O., & Behounková, M. (2019). Is Airy isostasy applicable to icy moons? Geophysical Research Letters, 46(24), 14299–14306.
https://doi.org/10.1029/2019gl085903
Čadek, O., Souček, O., Běhounková, M., Choblet, G., Tobie, G., & Hron, J. (2019). Long-term stability of Enceladus’ uneven ice shell. Icarus,
319, 476–484. https://doi.org/10.1016/j.icarus.2018.10.003
Cadek, O., Tobie, G., Van Hoolst, T., Masse, M., Choblet, G., Lefevre, A., et al. (2016). Enceladus's internal ocean and ice shell constrained from
Cassini gravity, shape, and libration data. Geophysical Research Letters, 43(11), 5653–5660. https://doi.org/10.1002/2016gl068634
Crow-Willard, E. N., & Pappalardo, R. T. (2015). Structural mapping of Enceladus and implications for formation of tectonized regions. Journalof Geophysical Research-Planets, 120(5), 928–950. https://doi.org/10.1002/2015je004818
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
Elson, L. S., Acton, C. H., Conner, D. L., & Semenov, B. V. (2005). Cassini SPICE Kernels V1.0, CO-S/J/E/V-SPICE-6-V1.0. NASA Planetary
Data System. https://doi.org/10.17189/1520114
Ermakov, A. I., Fu, R. R., Castillo-Rogez, J. C., Raymond, C. A., Park, R. S., Preusker, F., et al. (2017). Constraints on Ceres' internal structureand evolution from its shape and gravity measured by the Dawn spacecraft. Journal of Geophysical Research: Planets, 122(11), 2267–2293.
https://doi.org/10.1002/2017JE005302
Ermakov, A. I., Park, R. S., Roa, J., Castillo-Rogez, J. C., Keane, J. T., Nimmo, F., et al. (2021). A recipe for the geophysical exploration of
Enceladus. Planetary Science Journal, 2(4), 157. https://doi.org/10.3847/psj/ac06d2
Gaskell, R. W., Barnouin, O. S., Daly, M. G., Palmer, E. E., Weirich, J. R., Ernst, C. M., et al. (2023). Stereophotoclinometry on the OSIRIS-REx
Mission: Mathematics and methods. Planetary Science Journal, 4(4), 63. https://doi.org/10.3847/psj/acc4b9
Giese, B., Wagner, R., Hussmann, H., Neukum, G., Perry, J., Helfenstein, P., & Thomas, P. C. (2008). Enceladus: An estimate of heat flux and lithospheric thickness from flexurally supported topography. Geophysical Research Letters, 35(24), L24204. https://doi.org/10.1029/2008gl036149
Gillam, S. D., Owen, W. M., Vaughan, A. T., Wang, T.-C. M., Costello, J. D., Jacobson, R. A., et al. (2007). Optical navigation for the Cassini/
Huygens mission. In Paper presented at the AIAA/AAS Astrodynamics Specialist Conference, Mackinac Island, MI.
Hansen, C. J., Esposito, L., Stewart, A. I. F., Colwell, J., Hendrix, A., Pryor, W., et al. (2006). Enceladus' water vapor plume. Science, 311(5766),
1422–1425. https://doi.org/10.1126/science.1121254
Hemingway, D. J., Iess, L., Tadjeddine, R., & Tobie, G. (2018). The interior of Enceladus. In Enceladus and the Icy Moons of Saturn. The
University of Arizona Press.
Hemingway, D. J., & Matsuyama, I. (2017). Isostatic equilibrium in spherical coordinates and implications for crustal thickness on the Moon,
Mars, Enceladus, and elsewhere. Geophysical Research Letters, 44(15), 7695–7705. https://doi.org/10.1002/2017gl073334
Hemingway, D. J., & Mittal, T. (2019). Enceladus's ice shell structure as a window on internal heat production. Icarus, 332, 111–131. https://doi.
org/10.1016/j.icarus.2019.03.011
Iess, L., Stevenson, D. J., Parisi, M., Hemingway, D., Jacobson, R. A., Lunine, J. I., et al. (2014). The gravity field and interior structure of Enceladus. Science, 344(6179), 78–80. https://doi.org/10.1126/science.1250551
Ingersoll, A. P., Ewald, S. P., & Trumbo, S. K. (2020). Time variability of the Enceladus plumes: Orbital periods, decadal periods, and aperiodic
change. Icarus, 344, 113345. https://doi.org/10.1016/j.icarus.2019.06.006
Jacobson, R. A. (2022). The orbits of the Main Saturnian satellites, the Saturnian system gravity field, and the orientation of Saturn's pole. The
Astronomical Journal, 164(5), 199. https://doi.org/10.3847/1538-3881/ac90c9
Jacobson, R. A., Brozovic, M., Mastrodemos, N., Riedel, J. E., & Sheppard, S. S. (2022). Ephemerides of the Irregular Saturnian Satellites from
Earth-based Astrometry and Cassini Imaging. The Astronomical Journal, 164(6), 240. https://doi.org/10.3847/1538-3881/ac98c7
Jaumann, R., Bell, J. F., Polanskey, C. A., Raymond, C. A., Aspaugh, E., Bercovici, D., et al. (2022). The Psyche topography and geomorphology
investigation. Space Science Reviews, 218(2), 7. https://doi.org/10.1007/s11214-022-00874-7
Kinczyk, M. J., Patterson, G. W., Perkins, R. P., Collins, G. C., Borrelli, M., Becker, T. L., et al. (2017). Evaluation of impact crater distributions
for geological terrains on Enceladus. Retrieved from https://ui.adsabs.harvard.edu/abs/2017LPI….48.2926K
Knowles, B. (2018). Cassini Imaging Science Subsystem (ISS) Data User’s Guide. https://doi.org/10.17189/1504135

26 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE008054

Leonard, E. J., Yin, A., & Pappalardo, R. T. (2021). Forming relic cratered blocks: Left-lateral shear on Enceladus inferred from ice-shell deformation in the leading hemisphere. Journal of Geophysical Research-Planets, 126(2), e2020JE006499. https://doi.org/10.1029/2020je006499
Martin, E. S. (2016). The distribution and characterization of strike-slip faults on Enceladus. Geophysical Research Letters, 43(6), 2456–2464.
https://doi.org/10.1002/2016gl067805
Martin, E. S., Kattenhorn, S. A., Collins, G. C., Michaud, R. L., Pappalardo, R. T., & Wyrick, D. Y. (2017). Pit chains on Enceladus signal the
recent tectonic dissection of the ancient cratered terrains. Icarus, 294, 209–217. https://doi.org/10.1016/j.icarus.2017.03.014
Mastrodemos, N., Rush, B., Vaughan, A., & Owen, O. (2012). Optical navigation for the Dawn mission at Vesta. In Paper presented at the 23rd
International Symposium on Space Flight Dynamics, Pasadena, CA.
McEwen, A. S. (1991). Photometric functions for photoclinometry and other applications. Icarus, 92(2), 298–311. https://doi.
org/10.1016/0019-1035(91)90053-v
McKinnon, W. B. (2015). Effect of Enceladus's rapid synchronous spin on interpretation of Cassini gravity. Geophysical Research Letters, 42(7),2137–2143. https://doi.org/10.1002/2015GL063384
Nadezhdina, I. E., Zubarev, A. E., Brusnikin, E. S., & Oberst, J. (2016). A libration model for Enceladus based on geodetic control point network
analysis. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 41, 459–462. https://doi.
org/10.5194/isprs-archives-xli-b4-459-2016
Nimmo, F., Bills, B. G., & Thomas, P. C. (2011). Geophysical implications of the long-wavelength topography of the Saturnian satellites. Journalof Geophysical Research, 116(E11), E11001. https://doi.org/10.1029/2011je003835
Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. (2021). The JPL Planetary and Lunar Ephemerides DE440 and DE441. The Astronomical Journal, 161(3), 105. https://doi.org/10.3847/1538-3881/abd414
Park, R. S., Konopliv, A. S., Asmar, S. W., Bills, B. G., Gaskell, R. W., Raymond, C. A., et al. (2014). Gravity field expansion in ellipsoidal
harmonic and polyhedral internal representations applied to Vesta. Icarus, 240, 118–132. https://doi.org/10.1016/j.icarus.2013.12.005
Park, R. S., Konopliv, A. S., Bills, B. G., Rambaux, N., Castillo-Rogez, J. C., Raymond, C. A., et al. (2016). A partially differentiated interior for(1) Ceres deduced from its gravity field and shape. Nature, 537(7621), 515–517. https://doi.org/10.1038/nature18955
Park, R. S., Konopliv, A. S., Ermakov, A. I., Castillo-Rogez, J. C., Fu, R. R., Hughson, K. H. G., et al. (2020). Evidence of non-uniform crust of
Ceres from Dawn's high-resolution gravity data. Nature Astronomy, 4(8), 748–755. https://doi.org/10.1038/s41550-020-1019-1
Park, R. S., Riedel, J. E., Ermakov, A. I., Roa, J., Castillo-Rogez, J., Davies, A. G., et al. (2020). Advanced Pointing Imaging Camera (APIC) for
planetary science and mission opportunities. Planetary and Space Science, 194, 105095. https://doi.org/10.1016/j.pss.2020.105095
Park, R. S., Vaughan, A. T., Konopliv, A. S., Ermakov, A. I., Mastrodemos, N., Castillo-Rogez, J. C., et al. (2019). High-resolution shape modelof Ceres from stereophotoclinometry using Dawn Imaging Data. Icarus, 319, 812–827. https://doi.org/10.1016/j.icarus.2018.10.024
Passey, Q. R. (1983). Viscosity of the lithosphere of Enceladus. Icarus, 53(1), 105–120. https://doi.org/10.1016/0019-1035(83)90024-6
Patterson, G. W., Kattenhorn, S. A., Helfenstein, P., Collins, G. C., & Pappalardo, R. T. (2018). The geology of Enceladus. In Enceladus and theIcy Moons of Saturn. University of Arizona Press.
Patthoff, D. A., Pappalardo, R. T., Golombek, M., Chilton, H., Crow-Willard, E., & Thomas, P. C. (2022). Thrust faulting as the origin of dorsa
in the trailing hemisphere of Enceladus. Icarus, 375, 114815. https://doi.org/10.1016/j.icarus.2021.114815
Porco, C. C., Helfenstein, P., Thomas, P. C., Ingersoll, A. P., Wisdom, J., West, R., et al. (2006). Cassini observes the active South Pole of Enceladus. Science, 311(5766), 1393–1401. https://doi.org/10.1126/science.1123013
Postberg, F., Kempf, S., Schmidt, J., Brilliantov, N., Beinsen, A., Abel, B., et al. (2009). Sodium salts in E-ring ice grains from an ocean below
the surface of Enceladus. Nature, 459(7250), 1098–1101. https://doi.org/10.1038/nature08046
Schenk, P. M., & McKinnon, W. B. (2024). New global topography of Enceladus: Hypsometry, basins, spherical harmonics, shell thickness, andtrue polar wander revisited. Icarus, 408, 115827. https://doi.org/10.1016/j.icarus.2023.115827
Schroder, S. E., Mottola, S., Carsenty, U., Ciarniello, M., Jaumann, R., Li, J. Y., et al. (2017). Resolved spectrophotometric properties of the Ceres
surface from Dawn Framing Camera images. Icarus, 288, 201–225. https://doi.org/10.1016/j.icarus.2017.01.026
Spencer, J. R., Barr, A. C., Esposito, L. W., Helfenstein, P., Ingersoll, A. P., Jaumann, R., et al. (2009). Enceladus: An Active Cryovolcanic
Satellite. In M. K. Dougherty, L. W. Esposito, & S. M. Krimigis (Eds.), Saturn from Cassini-Huygens (pp. 683–724). Springer Netherlands.
Spencer, J. R., Pearl, J. C., Segura, M., Flasar, F. M., Mamoutkine, A., Romani, P., et al. (2006). Cassini encounters Enceladus: Background andthe discovery of a south polar hot spot. Science, 311(5766), 1401–1405. https://doi.org/10.1126/science.1121661
Tajeddine, R., Soderlund, K. M., Thomas, P. C., Helfenstein, P., Hedman, M. M., Burns, J. A., & Schenk, P. M. (2017). True polar wander ofEnceladus from topographic data. Icarus, 295, 46–60. https://doi.org/10.1016/j.icarus.2017.04.019
Thomas, P. C., Tajeddine, R., Tiscareno, M. S., Burns, J. A., Joseph, J., Loredo, T. J., et al. (2016). Enceladus's measured physical libration
requires a global subsurface ocean. Icarus, 264, 37–47. https://doi.org/10.1016/j.icarus.2015.08.037
Tricarico, P. (2014). Multi-layer hydrostatic equilibrium of planets and synchronous moons: Theory and application to Ceres and to solar system
moons. The Astrophysical Journal, 782(2), 99. https://doi.org/10.1088/0004-637x/782/2/99
Van Hoolst, T., Baland, R. M., & Trinh, A. (2016). The diurnal libration and interior structure of Enceladus. Icarus, 277, 311–318. https://doi.
org/10.1016/j.icarus.2016.05.025
Yin, A., & Pappalardo, R. T. (2015). Gravitational spreading, bookshelf faulting, and tectonic evolution of the South Polar Terrain of Saturn's
moon Enceladus. Icarus, 260, 409–439. https://doi.org/10.1016/j.icarus.2015.07.017

PARK ET AL.

27 of 27

21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054, Wiley Online Library on [15/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets