of 27
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 plume
that 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 and
imaging data. The DSN data are the primary observations used to determine the gravity field of Enceladus (Iess
et al.,
2014
), whereas the imaging data are the primary observations used for determining shape and rotational
parameters (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 through
the crust is not well known.
Here, we analyze Cassini radiometric tracking and onboard imaging data acquired during close encounters
of Enceladus using the latest reconstructed ephemerides of Enceladus and the Cassini spacecraft (Jacobson
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 during
close encounters with the moon. We compute the global shape, gravity field, 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. We recover a quadrupole gravity field with
J
3
and a forced libration amplitude of 0.091° ± 0.009° (3-
σ
). We also compute a global shape model using a
stereo-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 a
mean core density range of 2,270–2,330 kg/m
3
.
PARK ET AL.
© 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 and
reproduction in any medium, provided the
original work is properly cited and is not
used for commercial purposes.
The Global Shape, Gravity Field, and Libration of Enceladus
R. S. Park
1
, N. Mastrodemos
1
, R. A. Jacobson
1
, A. Berne
2
, A. T. Vaughan
1
,
D. J. Hemingway
3
, E. J. Leonard
1
, J. C. Castillo-Rogez
1
, C. S. Cockell
1,4
, J. T. Keane
1
,
A. S. Konopliv
1
, F. Nimmo
5
, J. E. Riedel
1
, M. Simons
2
, and S. Vance
1
1
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA,
2
Seismological Laboratory, Division
of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA,
3
Institute for Geophysics,
University of Texas at Austin, Austin, TX, USA,
4
UK Centre for Astrobiology, School of Physics and Astronomy, University
of Edinburgh, Edinburgh, UK,
5
Department of Earth and Planetary Sciences, University of California Santa Cruz, Santa
Cruz, CA, USA
Key Points:
A full quadrupole gravity field with
J
3
and the forced libration amplitude of
0.091° ± 0.009° are recovered
A 500-m resolution global topography
model 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 libration
of Enceladus.
Journal of Geophysical
Research: 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
10.1029/2023JE008054
RESEARCH ARTICLE
1 of 27
Journal of Geophysical Research: Planets
PARK ET AL.
10.1029/2023JE008054
2 of 27
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
J
3
and a forced libration amplitude of 0.091° ± 0.009°, including a complete rota
-
tional 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 core
density 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 resolution
of 500 m, with some local regions having resolutions as high as 25–100 m. All shape and orientation products are
available 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 aspects
of photogrammetry and photoclinometry. SPC has been used over the last 20 years for both proximity navigation
and 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 tiling
schemes 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 coordinates
with 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 identi
-
fied. 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 not
enough 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 diversity
of images at a pixel scale comparable to the map scale, images with a pixel scale larger than 2 times of the
map 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 bright
-
ness profile and the image template brightness profile are minimized in a weighted least squares sense to
solve 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
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
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
3 of 27
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 stages
of 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 esti
-
mate albedo information so that we can separate the brightness variation due to the local topography from
brightness due to intrinsic albedo variations.
5.
The slope solution is integrated into local heights using seed points from existing topography, such as an a
priori 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 the
limbs 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 suffi
-
ciently 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 constructed
to 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 a
minimum 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 overlap
-
ping 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 ephe
-
merides 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 the
allowable 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
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.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
4 of 27
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 a
surface 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 of
artifacts 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 at
the NAIF website (
https://naif.jpl.nasa.gov/pub/naif/CASSINI
). The Enceladus-spacecraft relative position was
derived 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 rota
-
tional frame was improved throughout the shape development process, and the updated versions were incorpo
-
rated 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 ellip
-
soid 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 per
component 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 rela
-
tive 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 topography
height to be driven by the imaging data rather than being artificially constrained to remain very close to the ellip
-
soid 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 of
overlapping maps. In essence, overlapping maps enforce a gradual and self-consistent height adjustment over the
entire 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 range
of such phase functions may produce good visual results, depending on the surface properties of the body, some
are 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
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
5 of 27
used were the Lunar-Lambert (McEwen,
1991
), Lommel-Seeliger with no phase term and the Ceres parameter
-
ized 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 baseline
model.
2.3.
Shape Model Development
SPC is an iterative process with a few steps repeated many times until final model convergence. We use the
following 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 these
maps 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 adjust
-
ing 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 or
with large residuals in the adjacent map overlaps, or deleting images whose large residuals could not be recon
-
ciled 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 of
images with similar resolution, emission, incidence, camera, and solar azimuth angles in the same map were
deleted. Images with large residuals that offer a unique observational geometry in that particular map were
retained. In general, images with residuals up to a value of 3-
σ
of the total root-mean-square (RMS) of the
residuals 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 the
image-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. Once
a 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 maps
on the surface are shown in Figure
2
in grayscale. Lighter-colored areas indicate regions where a higher density
of 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 were
modeled to the highest map spatial resolution level allowed by the image resolution. One such contiguous area
of 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 high
phase angle with parts of the images being too dark to be used. The combination of high phase and large vertical
features create self-shadowing, causing the quality of topographic reconstruction in these areas to be quite varia
-
ble 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 seen
in the topography. For example, Figures
4a
and
4b
display the topography in that area rendered at the viewing and
illumination geometry of the respective images, as shown in Figures
4c
and
4d
.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
6 of 27
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 Bagh
-
dad region, most of these images are of low emission and low phase angle. The topography of this area is shown
in 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 the
last 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 the
posteriori 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 that
location 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 from
all 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 a
subregion of the Baghdad region shown in panel (b), constructed at 25 m scale.
Figure 2.
The Enceladus surface coverage with 500 m maps in gray scale in the cylindrical projection.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
7 of 27
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 these
areas. 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.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
8 of 27
2.4.
SPC Merit Study
The suitability of an imaging data set for a stereophotoclinometry shape
model 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 the
surface.
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 diffi
-
cult. The incident, emission, and phase angles for a set of images that have
visibility of a surface location are scalar quantities that do not fully describe
the solar geometry. For example, the incidence angle could be 45° for several
pictures, which would seem deficient, but if they were a mix of morning
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 found
that each of these terms has something to contribute to the overall SPC merit, and so each of them is included in
the analysis.
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 are
summarized in Table
4
; however, the discrete binning method has significant drawbacks. Instead, we propose a
coverage
approach to evaluating each of these criteria. For each metric, we establish a scorable range and then
overlay 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 meas
-
urements 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 metric
in a notional case. Two images that have visibility of the surface patch in question have two different complemen
-
tary 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 measurement
that occurs very near one of the existing measurements could widen the coverage marginally but would not have
a 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
“buckets” such as those described in Jaumann et al. (
2022
)) is ambivalent to
the clock angle of the measurements. For example, two images with camera
azimuth angles of 89° and 91° would get credit for two buckets in the binned
approach, while our algorithm considers those measurements
as nearly the
same and would be scored significantly lower. Similarly, two measurements
that came at 1° and 89° would fall in the same bucket in the binning system
but would be strongly rewarded as diverse by our coverage approach. In addi
-
tion to the sinusoid, we can also apply a square wave to the coverage arc.
This is used only by the nadir measurement, which we intend to be binary.
If there is a single image in the scorable range, the scalar score for the nadir
metric will be 1.
Map scale (m)
Number of maps (control points)
Average number
of images per map
1,000
335
82
500
1,105
70
200
923
25
100
739
11
50
494
4
35
258
4
25
438
4
Table 1
Number of Individual Topography Maps and Associated Control Points, per
Map Scale
Image pixel scale range
Number of images
Average number of
control points per image
scale >1 km
184
34
500 m < scale <1 km
296
173
200 m < scale <500 m
401
153
100 m < scale <200 m
136
113
50 m < scale <100 m
69
69
scale <50 m
25
136
Table 2
Image Scale Distribution and Statistics of Control Points
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
9 of 27
There is a correlation between the angle measurements used as metrics in
the analysis. These correlations are implemented in the algorithm as basic
filtering. For example, it would be overly optimistic to score all measures of
incident angle in each image set if we knew that many of those images had
emission angles so high that they would be disqualified from use. For all
metrics except nadir, we pre-filter the measurements by the following ranges:
incidence must be between 10° and 80°, emission must be between 0° and
80°, and phase must be between 10° and 90°. For the nadir criteria (which is a
metric of emission angle), we filter
incidence
between 20° and 70° and phase
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 an
image with an emission angle of 65° toward the emission metric, we would allow this picture to count toward the
incidence 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.
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 topog
-
raphy 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 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 that metric. 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 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 to
covering 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 an
emission 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 score
of 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 this
process. We could use every possible image available from the image server, but this ignores the fact that many
of 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 the
shape model. Thus, this analysis answers the question: given the set of images found to suitable for processing
Map scale
Control point height 1-
σ
uncertainty (m)
RMS of residuals
(m)
Total number
of observations
1 km
286
439
41,676
500 m
297
300
97,295
Combined
295
338
138,971
All maps
252
293
204,503
Table 3
Statistics for the Global Topography Model
Criteria
Range (°)
Minimum case
Optimal case
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 apart
4
0–90, 90–180, 180–270, 270–360
Note.
Also see Jaumann et al. (
2022
) for similar criteria for stereophotoclinometry.
Table 4
Qualitative Criteria for Evaluating an SPC Imaging Data Set
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
10 of 27
at all, how effective should this data be for producing an SPC model? A total of 916 images contribute to this
analysis, filtered by their suitability for modeling as just mentioned and also their mean resolution, in the range
of 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 is
located at (0.3895, 0.7018, 0.3187) km, with an uncertainty of 1 km. Table
5
shows the parameters for the best-
fit ellipsoid and best-fit spheroid for the new SPC-based Enceladus shape model that could be used for interior
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 picture
that 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.
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 that
metric. 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 this
coverage fades off to yellow and then red as the sinusoid tails off.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
11 of 27
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 hemi
-
sphere 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 the
interior 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 ellipsoid
at ∼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 S2
basin 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.
We undertook three different verification steps for these regions, each
method 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 zero
a priori slope and albedo values in each map. The second round included the
first one but also started with a smooth triaxial ellipsoid as the a priori topog
-
raphy, that is, all prior topography information was erased in those areas.
The third round started from a clean slate; image position and pointing were
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.
Cases
a
e
(km)
b
e
(km)
c
e
(km)
Mean radius (km)
Best-fit ellipsoid
256.14
251.16
248.68
251.99
Best-fit spheroid
253.63
248.67
251.98
Note.
Both the ellipsoid and spheroid estimates are provided since some
astrometry-related applications might prefer one model over the other. For
the best-fit spheroid
a
e
=
b
e
.
Table 5
Best-Fit Ellipsoid and Spheroid Parameters
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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
PARK ET AL.
10.1029/2023JE008054
12 of 27
reset to their a priori values, all topography maps were deleted, and the topography process was performed for
these specific areas anew starting with a smooth ellipsoid as an a priori surface. In all these tests, we were able to
recreate both of these features. We have also verified that these areas exhibit similar features for the topography
models 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 map
that 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
Figure 9.
Cylindrical projection of the SPC topography for ±75° latitude (a), stereo projection of the Northern hemisphere topography for 90°–0° latitude (b), and
stereo projection of the Southern hemisphere topography for −90°–0° latitude (c). The horizontal map resolution is 500 m and 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 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.
21699100, 2024, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE008054 by California Inst of Technology, Wiley Online Library on [20/06/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