UNIVERSITY OF CALIFORNIASANTA CRUZ
LUNAR MAGNETISM, SPACE WEATHERING,
AND ICY SATELLITE INTERIORS
A dissertation submitted in partial satisfactionof the requirements for the degree of
DOCTOR OF PHILOSOPHY
in
EARTH SCIENCES
with an emphasis in PLANETARY SCIENCE
byDouglas Hemingway
September 2015
The Dissertation of Douglas Hemingway
is approved:

Professor Francis Nimmo, chair

Professor Ian Garrick-Bethell

Professor Gary Glatzmaier

Professor Jasper Halekas
Tyrus Miller
Vice Provost and Dean of Graduate Studies

Copyright © by
Douglas Hemingway
2015

Table of Contents
List of Figuresvi

Abstractviii

Acknowledgementsx

1 Introduction

1

2 Magnetic Field Direction and Lunar Swirl Morphology
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Predicted Influence of Magnetic Field Direction . . . .
2.2 Data Processing . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Lunar Prospector Magnetometer Data . . . . . . . . .
2.2.2 Clementine Reflectance Mosaics . . . . . . . . . . . . .
2.3 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Airy . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Reiner Gamma . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Source Modeling . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Lunar Prospector Data Fitting . . . . . . . . . . . . .
2.4.1.1 Airy Single Dipole Model . . . . . . . . . . .
2.4.1.2 Reiner Gamma Single Dipole Model . . . . .
2.4.1.3 Genetic Search Algorithm . . . . . . . . . . .
2.4.1.4 Airy Dipole Grid Model . . . . . . . . . . . .
2.4.1.5 Reiner Gamma Dipole Grid Model . . . . . .
2.4.2 Albedo Pattern Matching . . . . . . . . . . . . . . . .
2.4.2.1 Airy Model . . . . . . . . . . . . . . . . . . .
2.4.2.2 Reiner Gamma Model . . . . . . . . . . . . .
2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Magnetization . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Magnetizing Field . . . . . . . . . . . . . . . . . . . .
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
7
7
9
12
12
17
17
18
20
22
23
24
25
25
26
29
31
34
36
37
41
41
46
47

iii

3 Latitudinal Variation in Spectral Properties of the Lunar Mariaand Implications for Space Weathering49
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
3.2 Data Sources . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
3.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
3.3.1 Color Variation at Lunar Swirls . . . . . . . . . . . . .
53
3.3.2 Latitudinal Color Variation . . . . . . . . . . . . . . .
58
3.3.2.1 Eﬀect of Composition . . . . . . . . . . . . .
59
3.3.2.2 Isolation of the Maria . . . . . . . . . . . . .
63
3.3.2.3 Comparison with Swirl Trends . . . . . . . .
66
3.3.2.4 Mare Frigoris . . . . . . . . . . . . . . . . . .
67
3.3.2.5 Highlands Contamination . . . . . . . . . . .
71
3.3.2.6 Phase Angle Biases and LOLA Reflectance .
75
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
76
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
4 A Rigid and Weathered Ice Shell on Titan
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Gravity Field . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Topography . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Observed Admittance . . . . . . . . . . . . . . . . . . .
4.3 Model and Results . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Admittance and Flexure Model . . . . . . . . . . . . .
4.3.2 Interpretation . . . . . . . . . . . . . . . . . . . . . . .
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Ice Shell Structure and Rigidity . . . . . . . . . . . . .
4.4.1.1 Temperature Profile . . . . . . . . . . . . . .
4.4.1.2 Heat Flux . . . . . . . . . . . . . . . . . . . .
4.4.1.3 Origin of Long-Wavelength Topography . . .
4.4.1.4 Lateral Flow in the Ice Shell . . . . . . . . . .
4.4.1.5 Eﬀect on Tidal k2 . . . . . . . . . . . . . . .
4.4.2 Erosion and Deposition . . . . . . . . . . . . . . . . . .
4.4.3 Degree-2 Admittance and Fluid Love Number . . . . .
4.4.4 Degree-4 Predictions . . . . . . . . . . . . . . . . . . .
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84
85
87
87
89
93
97
97
106
112
112
113
114
114
115
117
118
118
122
123

5 Internal Structure of Enceladus
126
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.2 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.2.1 Gravity Field . . . . . . . . . . . . . . . . . . . . . . . 129
5.2.2 Topography . . . . . . . . . . . . . . . . . . . . . . . . 130
iv

5.3

5.4

5.5

5.2.3 Deviation from Hydrostatic Equilibrium . . . . . . . .
Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Separation of Hydrostatic Component . . . . . . . . . .
5.3.2 Model Admittance . . . . . . . . . . . . . . . . . . . .
5.3.3 Elastic Flexure . . . . . . . . . . . . . . . . . . . . . .
5.3.4 Moment of Inertia . . . . . . . . . . . . . . . . . . . .
Results and Discussion . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Admittance and Moment of Inertia . . . . . . . . . . .
5.4.2 Isostatic Compensation . . . . . . . . . . . . . . . . . .
5.4.3 Elastic Support . . . . . . . . . . . . . . . . . . . . . .
5.4.4 Internal Structure . . . . . . . . . . . . . . . . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A Spherical Harmonics

131
133
133
135
135
137
137
137
140
141
142
142
146

B Magnetic Potential and Source Models
150
B.1 Magnetic Scalar Potential and Vector Magnetic Fields . . . . . 150
B.2 Source Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
B.2.1 Single Dipole Source . . . . . . . . . . . . . . . . . . . 152
B.2.2 Linear Source . . . . . . . . . . . . . . . . . . . . . . . 154
C Gravitational Potential and Figures of Equilibrium
156
C.1 Gravitational Potential . . . . . . . . . . . . . . . . . . . . . . 156
C.2 Tidal and Rotational Deforming Potentials . . . . . . . . . . . 157
C.2.1 Tidal Deforming Potential . . . . . . . . . . . . . . . . 157
C.2.2 Rotational Deforming Potential . . . . . . . . . . . . . 160
C.2.3 Combined Deforming Potential . . . . . . . . . . . . . 160
C.3 Equilibrium Figures and Potentials . . . . . . . . . . . . . . . 161
C.3.1 Rotational Deformation . . . . . . . . . . . . . . . . . 162
C.3.1.1 Hydrostatic Shape . . . . . . . . . . . . . . . 162
C.3.1.2 Hydrostatic Gravity . . . . . . . . . . . . . . 163
C.3.1.3 Geoid . . . . . . . . . . . . . . . . . . . . . . 164
C.3.1.4 Reference Ellipsoid . . . . . . . . . . . . . . . 165
C.3.2 Tidal/Rotational Deformation . . . . . . . . . . . . . . 166
C.3.2.1 Hydrostatic Shape . . . . . . . . . . . . . . . 166
C.3.2.2 Hydrostatic Gravity . . . . . . . . . . . . . . 168
C.3.2.3 Geoid . . . . . . . . . . . . . . . . . . . . . . 169
C.3.2.4 Reference Ellipsoid . . . . . . . . . . . . . . . 170
Bibliography

172v

List of Figures
2.1 Magnetic Fields Due to Vertical vs. Horizontal Dipoles . . . .
2.2 Magnetic Field Maps at Airy and Reiner Gamma . . . . . . .
2.3 Eastward Magnetic Field at Airy . . . . . . . . . . . . . . . .
2.4 Albedo to Field Direction Correlation at Reiner Gamma . . .
2.5 Locations of LP-MAG data points at Reiner Gamma . . . . .
2.6 Dipole Grid Model Evolution for Airy and Reiner Gamma . .
2.7 Grid Model Error Evolution . . . . . . . . . . . . . . . . . . .
2.8 Final Grid Model for Airy . . . . . . . . . . . . . . . . . . . .
2.9 Lunar Prospector Data versus Grid Model for Airy . . . . . .
2.10 Final Grid Model for Reiner Gamma . . . . . . . . . . . . . .
2.11 Lunar Prospector Data versus Grid Model for Reiner Gamma
2.12 Reiner Gamma Source Model from Nicholas et al. (2007) . . .
2.13 Albedo-Derived Model for Airy . . . . . . . . . . . . . . . . .
2.14 LP-MAG versus Albedo-Derived Model for Airy . . . . . . . .
2.15 Magnetic Fields Due to Two Horizontal Dipoles . . . . . . . .
2.16 Albedo-Derived Model for Reiner Gamma . . . . . . . . . . .
2.17 LP-MAG versus Albedo-Derived Model for Reiner Gamma . .
2.18 Rock Magnetization Implied by Grid Models . . . . . . . . . .

11
19
20
22
23
30
31
32
33
34
35
36
38
39
42
44
45
46

3.1 Spectral Characteristics at Reiner Gamma . . . . . . . . . . .
3.2 Spectral Characteristics at Mare Ingenii . . . . . . . . . . . .
3.3 Spectral Characteristics at Mare Marginis . . . . . . . . . . .
3.4 Color variation as a function of latitude . . . . . . . . . . . . .
3.5 Latitudinal Color Variation for 12 Compositional Bins . . . . .
3.6 Latitudinal Trends by Iron Content . . . . . . . . . . . . . . .
3.7 Topographic Roughness Histogram . . . . . . . . . . . . . . .
3.8 Mare Region Mask . . . . . . . . . . . . . . . . . . . . . . . .
3.9 Clementine and LOLA Reflectance Profiles . . . . . . . . . . .
3.10 Latitudinal Color Variation Across the Maria . . . . . . . . . .
3.11 Latitudinal Trend Slopes versus Swirl and Impact Trend Slopes
3.12 Swirl-like Color Variation Across the Maria . . . . . . . . . . .
3.13 Composition and Color of Mare Frigoris . . . . . . . . . . . .

54
55
56
60
62
63
65
65
66
67
68
68
70

vi

3.14 Latitudinal Trends With and Without Mare Frigoris . . . . . .
3.15 Color Variation of Maria >50 km from Highlands . . . . . . .
3.16 Phase Angle Biases in Clementine Data . . . . . . . . . . . . .
3.17 Hypothesized Optical Evolution Process . . . . . . . . . . . .

70
73
77
80

4.1 Radar-derived Elevation Data for Titan . . . . . . . . . . . . .
4.2 Degree-3 Topography Coeﬃcients . . . . . . . . . . . . . . . .
4.3 Titan’s degree-3 gravity and topography . . . . . . . . . . . .
4.4 Degree-3 Admittance Estimates . . . . . . . . . . . . . . . . .
4.5 Ice Shell Loading and Flexure . . . . . . . . . . . . . . . . . .
4.6 Admittance as a Function of Top Load . . . . . . . . . . . . .
4.7 Model predictions of admittance and erosion . . . . . . . . . .
4.8 Degree-3 Topography and Gravity Maps . . . . . . . . . . . .
4.9 Implied Erosion as a Function of Admittance . . . . . . . . . .
4.10 Titan’s Fluid Love Number . . . . . . . . . . . . . . . . . . .
4.11 Degree-2 Gravity Maps . . . . . . . . . . . . . . . . . . . . . .
4.12 Degree-4 Admittance Prediction . . . . . . . . . . . . . . . . .
4.13 Degree-4 Topography and Gravity Maps . . . . . . . . . . . .

91
91
94
95
100
107
109
110
112
120
121
123
124

5.1
5.2
5.3
5.4
5.5

130
131
138
143
144

Enceladus’ Gravity Field . . . . . . . . . . . . . . . . . . . . .
Enceladus’ Topography . . . . . . . . . . . . . . . . . . . . . .
Converging on Admittance and Moment of Inertia . . . . . . .
Admittance and Elastic Thickness . . . . . . . . . . . . . . . .
Moment of Inertia, Mantle Thickness, and Core Density . . . .

vii

Abstract
Lunar Magnetism, Space Weathering,and Icy Satellite Interiors
byDouglas Hemingway
An enduring mystery since Apollo is that, in spite of the Moon’s lack of a
global magnetic field, the surface is nevertheless dotted with regional magnetic
fields strong enough to be detected from orbit. Did the Moon once have an
intrinsic global field that magnetized parts of the crust but has since decayed
away? This is a question of fundamental importance to understanding the formation and evolution of solid planetary bodies, and yet it remains unanswered
due in part to limitations in our knowledge of these crustal magnetic anomalies.
Adding to the puzzle, many of these magnetic anomalies are accompanied byenigmatic optical features, known as swirls, which may hold the key to understanding "space weathering"—a process by which airless bodies change color
over time due to exposure to solar wind and micrometeoroids. Here we show
both that swirl morphology provides information about the structure of theunderlying magnetic sources, and that the color of the lunar surface varies
systematically with latitude in a way that allows us to distinguish between
the eﬀects of solar wind ion and micrometeoroid bombardment, addressing adecades-old problem in remote sensing, and aiding in the interpretation of thespectra of airless bodies throughout the solar system.
The remarkable diversity of the outer solar system’s satellites provides important clues about the formation and evolution of the solar system. Many ofviii

the satellites have surprisingly young surfaces, owing in some cases to on-goinggeologic activity. Moreover, the existence of subsurface oceans within some ofthe satellites raises the intriguing possibility of extant habitable environmentsin the outer solar system. Determining the properties of their ice shells andthe structures of their deep interiors places fundamental constraints on how
the icy satellites formed and evolved, and on what governs their behavior today. Using gravity and topography data from Cassini, we develop analyticalmodels showing that Titan’s ice shell may be very rigid and therefore unlikely
to be geologically active. In contrast, we also model the internal structureof the tiny, but highly geologically active moon Enceladus, and find that a
subsurface liquid ocean is likely.

ix

Acknowledgements
The text of this dissertation includes modified reprints of previously published material. Chapter 2 is a slightly modified reprint of Hemingway andGarrick-Bethell (2012); the co-author (Garrick-Bethell) supervised the research. Chapter 3 is a slightly modified version of Hemingway et al. (2015);
the first co-author (Garrick-Bethell) supervised the research; the second coauthor (Kreslavsky) supported the research and carried out the determinationof phase angles associated with the Clementine mosaics; the text of this paper is used with the permission of all co-authors. Chapter 4 is adapted from
Hemingway et al. (2013); the first co-author (Nimmo) supervised the research;
the second (Zebker) and third (Iess) co-authors, respectively, supplied the topography and gravity models used in the analysis; the text is used with thepermission of all co-authors. Chapter 5 describes research I carried out insupport of the multi-author publication Iess et al. (2014); only my own contributions, under the supervision of Professor Nimmo, are included in thisdissertation. Material contributed by other authors of Iess et al. (2014) is notincluded here.
I am extremely grateful for the unhesitating support of my academic advisors, Ian Garrick-Bethell and Francis Nimmo, for their advice and guidance,
and for the time they shared so generously with me over the past few years. Ixam also grateful to Erik Asphaug, who has always been very supportive andhelpful with general advice, especially during my first two years of grad school.
For their time and helpful advice, I am also thankful to Gary Glatzmaier andJasper Halekas, who have kindly served on my reading committee for the last
three years. Discussions and technical support from Lon Hood, Bob Lin, andJasper Halekas were very important during my first major research work, detailed in Chapter 2. The work detailed in Chapter 3 benefitted enormously
from numerous discussions with Carle Pieters and Misha Kreslavsky. An informative conversation with Bob Pappalardo helped me first find my way into thefield of icy satellites, and ultimately the work described in Chapters 4 and 5.
Technical discussions with Howard Zebker, Bill McKinnon and Dave Stevenson were also valuable, and led to some of the ideas discussed in Chapters 4
and 5.
I am grateful to have had an amazing support network of friends and colleagues in the Department of Earth & Planetary Sciences. Of course Jennifer
Fish deserves special mention, but really the excellent environment is madeby all of the staﬀ, faculty, and grad students, who genuinely take an interestin each other’s success. I appreciate everyone’s time, and the interest people
have taken in my work, especially my lab mates Naor Movshovitz (with whom
I have shared countless hours of valuable technical discussion time, both at
the white board and at the picnic table), Viranga Perera, Erinna Chen, Tina
Dwyer, Delia Santiago, and Carver Thomason. Finally, I am thankful to myfamily for their support, and especially to Lynne Harden, who has also beena critical sounding board so very often over the past four and a half years.

xi

Chapter 1
Introduction
Gravity and magnetic fields are among the most fundamental and powerfultools for probing the interiors of planetary bodies (Connerney, 2007; Wieczorek , 2007). An internally generated geomagnetic field is a clear indicationof sustained movement of an electrically conducting fluid within the interior.
Even if no such field is present, remanent magnetization in crustal rocks can
put constraints on the material properties of the crust and may indicate that a
geomagnetic field once existed, but has since decayed away, as may be the casefor the Moon (Garrick-Bethell et al., 2009). In addition to constraining theproperties of the crust and providing clues about the deeper interior, crustal
magnetic fields may also aﬀect the evolution of the surface itself by influencing
space weathering—the gradual optical evolution of exposed surfaces—throughmagnetic deflection of solar wind (Hood and Schubert, 1980). Gravity is fundamentally about mass distribution and is therefore a crucial tool for interior
structure modeling. With a combination of gravity and topography data, we
can measure the thickness and elastic properties of a body’s lithosphere andplace constraints on the deeper interior.
1This work divides broadly into two parts. In the first part (Chapters 2
and 3), we use optical and magnetic field data to study “lunar swirls”, a classof albedo anomalies associated with strongly magnetized parts of the crust.
Swirls oﬀer insights into both the geometry of lunar crustal magnetic sources
and the role played by solar wind in space weathering. In the second part
(Chapters 4 and 5), we develop analytical models and use them in combination
with topography and gravity field data to construct interior structure models
for two of Saturn’s icy moons, Titan and Enceladus, helping to constrain thelikely level of geologic activity on the former and confirming the presence of a
subsurface liquid ocean on the latter.
Although the Moon does not possess a global magnetic field today (Nesset al., 1967), parts of its crust are nevertheless magnetized strongly enough to
be detected from orbit (Dyal et al., 1970; Hood et al., 2001). The origin of thisremanent magnetization is unknown, but it could be the result of crustal rocks
having cooled in the presence of a dynamo field that is now extinct (GarrickBethell et al., 2009) or transient episodes of intense magnetization associatedwith basin-forming impact events (Hood and Artemieva, 2008; Mitchell et al.,2008). Curiously, many of the crustal magnetic anomalies are co-located with
enigmatic collections of morphologically complex bright markings known as“swirls” (Hood and Schubert, 1980). The mechanism for the formation of swirls
is not known, but their close association with magnetic anomalies suggests thatthey may be useful for studying the characteristics and origins of lunar crustal
magnetism. In Chapter 2, we use Lunar Prospector magnetometer data andClementine reflectance mosaics to test the hypothesis that swirls are formed
where crustal magnetic anomalies, acting as mini magnetospheres, shield portions of the surface from the darkening eﬀects of solar wind ion bombardment,2thereby leaving patches that appear bright compared with their surroundings(Hood and Schubert, 1980). We find evidence from two swirls with dissimilarmagnetization orientations, Airy and Reiner Gamma, that the bright partsof swirls correspond to dominantly horizontal magnetic fields and that the
intra-swirl dark lanes correspond to vertical fields at the surface—just what
we would predict if the swirls are in fact formed due to magnetic deflection ofsolar wind. We further show that source models constrained by the morphology of the swirls produce fields that are consistent with the Lunar Prospector
magnetometer observations. We conclude that solar wind deflection is thelikely mechanism for the formation of swirls and suggest that swirl morphology provides a previously unrecognized clue to the geometry of the underlying
magnetic source bodies (Hemingway and Garrick-Bethell , 2012).
Aside from its connection to lunar magnetic anomalies, space weathering
is an important phenomenon of its own as it complicates the interpretation ofspectroscopic observations of airless bodies everywhere (Hapke, 2001). Solarwind and micrometeoroids are thought to be the dominant agents of space
weathering, but their relative contributions are not yet well understood. In
Chapter 3, building from the findings of Garrick-Bethell et al. (2011), we characterize the unique spectral signature of lunar swirls and compare it against
the spectra of the lunar surface globally. Using Clementine mosaics, we finda previously unrecognized systematic latitudinal variation in the near-infrared
spectral properties of the lunar surface, and show that the characteristics ofthis latitudinal trend match those observed at swirls. We propose that reduced
solar wind flux, which should occur both at swirls and toward higher latitudes,
is the common mechanism behind these color variations. This model helps us
quantify the distinct eﬀects of solar wind and micrometeoroid weathering and3

could aid in interpreting the spectra of airless bodies throughout the solar
system (Hemingway et al., 2015).
Beginning with Chapter 4, we turn to the problem of modeling the interior
structures of icy satellites, specifically using gravity and topography data fromthe Cassini mission. Cassini has made numerous flybys of several of Saturn’s
moons, most notably, the giant Titan and the small but geologically hyperactive Enceladus—the subjects of Chapters 4 and 5, respectively. A handful ofthe flybys have been dedicated to determination of the moons’ low order gravity fields via radio tracking during close approaches (Iess et al., 2010, 2012,2014), resulting in good constraints up to spherical harmonic degree 3 in thecase of Titan. For Enceladus, only the degree-2 terms and the degree-3 zonal
term have been determined. The shape of Enceladus has been determinedthrough compilation of dozens of limb profiles obtained from the Imaging Science Subsystem (Nimmo et al., 2011). This technique has not been applicable
to Titan, however, due to the thick and opaque atmosphere. Instead, Titan’s
topography has been determined by compiling elevation data obtained from
Cassini ’s RADAR system (Stiles et al., 2009; Zebker et al., 2009, 2012). In
both of these chapters, we employ an admittance analysis which characterizes
the degree to which the topography is compensated, and at what depth, providing an estimate of the thickness of the icy crusts and their elastic properties.
In the case of Titan (Chapter 4), the gravity and topography show a strong
and unexpected inverse correlation at spherical harmonic degree 3 (long wavelengths that are not aﬀected by tidal and rotational distortion), suggestingthat the high standing topography is associated with negative mass anomalies
at depth. We develop analytical models that allow for such a condition, but
only if the ice shell is substantially rigid (and therefore unlikely to be geologi4

cally active), if it has been loaded primarily from below, and if the surface has
experienced several hundred meters of erosion and deposition, with sedimentbeing redistributed globally (Hemingway et al., 2013).
For Enceladus (Chapter 5), the gravity data are so limited that we areforced to take on the task of analyzing the degree-2 signal, which is complicated by the eﬀects of tidal and rotational distortion. Making the assumptionthat the compensation mechanism behaves isotropically, we are able to selfconsistently separate the degree-2 gravity and topography signals into theirhydrostatic (i.e., tidal/rotational) and non-hydrostatic components in a waythat allows us to estimate the degree of compensation, the depth at which it occurs, and the moment of inertia. We conclude that Enceladus is substantiallydiﬀerentiated, and that the topography is highly compensated, suggesting that
it is supported isostatically, consistent with the presence of a global subsurface
liquid layer (Iess et al., 2014).
Three appendices are included at the end of this document, and are as
much for my own reference as anything else. In support of Chapters 2, 4, and5, Appendix A describes the conventions and notation I use when working withspherical harmonics, along with some discussion of the various normalization
schemes that are in use. Appendix B details supporting equations for Chapter 2, including conversion from magnetic scalar potential to magnetic flux
density and modeling of magnetic fields generated by individual dipoles andlinear sources. Finally, in support of Chapters 4 and 5, Appendix C provides
a summary of many of the most relevant relations concerning gravitational
potential and equilibrium figures arising from tidal and rotational distortion,
including some derivations that I find useful to have written out in my own
preferred notation.
5

Chapter 2
Magnetic Field Direction and
Lunar Swirl Morphology
This chapter is a modified reprint of Hemingway, D and I. Garrick-Bethell(2012), Magnetic field direction and lunar swirl morphology: Insights from
Airy and Reiner Gamma, J. Geophys. Res. 117, E10012.

Abstract
Many of the Moon’s crustal magnetic anomalies are accompanied by highalbedo features known as swirls. A leading hypothesis suggests that swirls
are formed where crustal magnetic anomalies, acting as mini magnetospheres,shield portions of the surface from the darkening eﬀects of solar wind ion bombardment, thereby leaving patches that appear bright compared with their surroundings. If this hypothesis is correct, then magnetic field direction should
influence swirl morphology. Using Lunar Prospector magnetometer data andClementine reflectance mosaics, we find evidence that bright regions corre6

spond with dominantly horizontal magnetic fields at Reiner Gamma and thatvertical magnetic fields are associated with the intra-swirl dark lane at Airy.
We use a genetic search algorithm to model the distributions of magnetic source
material at both anomalies and we show that source models constrained by
the observed albedo pattern (i.e., strongly horizontal surface fields in brightareas, vertical surface fields in dark lanes) produce fields that are consistent
with the Lunar Prospector magnetometer measurements. These findings support the solar wind deflection hypothesis and may help to explain not only
the general form of swirls, but also the finer aspects of their morphology. Our
source models may also be used to make quantitative predictions of the nearsurface magnetic field, which must ultimately be tested with very low altitude
spacecraft measurements. If our predictions are correct, our models could
have implications for the structure of the underlying magnetic material andthe nature of the magnetizing field.

2.1

Introduction

2.1.1

Background

Although the Moon does not now possess a global magnetic field (Ness et al.,1967), remanent crustal magnetization has been identified on the surface andseveral stable regional magnetic fields have been detected from orbit (Dyalet al., 1970; Coleman et al., 1972; Lin, 1979; Lin et al., 1998; Hood et al.,2001). Curiously, many, but not all of these crustal magnetic anomalies are accompanied by sinuous patterns of anomalously high surface reflectance known
as swirls (Hood et al., 1979; Hood and Williams, 1989; Richmond et al., 2005;

7

Blewett et al., 2011).
Compared with their surroundings, swirls are optically immature, exhibit
spectrally distinct space weathering trends (Garrick-Bethell et al., 2011), andare depleted in hydroxyl molecules (Kramer et al., 2011). Swirls appear to
overprint local topography, having no detectable topographic or textural expression of their own (Neish et al., 2011). So far not identified anywhere else in
the solar system, swirls are unique natural laboratories where space weathering and crustal magnetism intersect. As such, their study could help address
important questions in lunar science including the Moon’s dynamo history(Garrick-Bethell et al., 2009; Dwyer et al., 2011; Le Bars et al., 2011), therelative influences of solar wind and micrometeoroid bombardment on space
weathering (Hapke, 2001; Vernazza et al., 2009), and the production and distribution of water over the lunar surface (Pieters et al., 2009; Kramer et al.,2011).
Electrostatic migration of dust has been proposed as a possible mechanism
for swirl formation (Garrick-Bethell et al., 2011), as have meteoroid or comet
impacts (Schultz and Srnka, 1980; Pinet et al., 2000; Starukhina and Shkuratov , 2004). Another model suggests that crustal magnetic fields act as mini
magnetospheres, deflecting the solar wind and protecting portions of the surface from the optical maturation and darkening eﬀects of proton bombardment(Hood and Schubert, 1980; Hood and Williams, 1989). This chapter aims to
make and test predictions based on this solar wind deflection model for swirl
formation.

82.1.2

Predicted Influence of Magnetic Field Direction

The solar wind deflection model suggests that solar wind ions are magneticallydeflected due to the Lorentz force. Neglecting the induced electric field, thecross product of particle velocity and magnetic field in the Lorentz Law means
the magnetic deflection force is maximized when particle velocity is perpendicular to the magnetic field and zero when it is parallel. Hence, if the highalbedo of swirls is the result of inhibited space weathering due to magnetic
deflection of solar wind ions, magnetic field direction should influence swirl
morphology. While the solar wind incidence angle varies, the ion flux at thesurface and thus any darkening eﬀects will be greatest when the ion trajectories are vertical. We therefore predict that portions of the crust that areshielded by dominantly horizontal magnetic fields should receive maximum
protection from the solar wind while portions of the crust associated withvertically oriented magnetic fields should experience protection only at highsolar wind incidence angles, when darkening eﬀects would be minimal anyway.
This suggests that the magnetic fields directly over the bright swirls should bedominantly horizontal and that away from swirls and in the intra-swirl darklanes, the fields may be either closer to vertical or too weak to oﬀer the surface
any protection from solar wind darkening.
Because magnetic field strength decreases rapidly with distance from thesource, swirl morphology should depend mainly on the very low altitude structure of the magnetic field. The structure of the field at higher altitudes, where
field strength is weaker, has less influence on solar wind deflection. Unfortunately, spacecraft observations are typically limited to higher altitudes andtherefore do not directly capture the structure of the near-surface magnetic

9

field. It is therefore important to understand the way in which the field patterns change with observation altitude, including the role that is played by the
direction of magnetization. Figure 2.1 illustrates the magnetic field due to asingle dipole that is either vertically oriented (panels a and c) or horizontallyoriented (panels b and d). The upper panels (a and b) illustrate magnetic
field lines as seen in vertical cross section while the lower panels (c and d)show profiles of the horizontal component of the magnetic field as observed
along the dashed lines shown in the upper panels (in both cases, the model
dipole is placed at the origin and arbitrarily assigned a magnetic moment of1012 Am2 ). Directly over the vertically oriented dipole, field lines are verticalat any altitude (Figure 2.1a). Peaks in the horizontal field strength appear toeither side and are separated by a distance equal to the observation altitude(Figure 2.1c). Directly over the horizontally oriented dipole, field lines arehorizontal at any altitude (Figure 2.1b). As the observer moves along an axis
parallel to the dipole direction, horizontal field strength decreases, reachingpzero (i.e., vertical field lines) at a horizontal distance of 1/ 2 times the observation altitude, before temporarily increasing again slightly (Figure 2.1d).
Appendix B gives a quantitative treatment of the way in which magnetic fields
vary with position relative to the source.
Our study uses Lunar Prospector magnetometer data collected at altitudesof ⇠ 18 km and higher. While we cannot measure the magnetic field structure
below these altitudes, we can distinguish between vertically and horizontallyoriented magnetizations (we assume that strong crustal magnetic anomalies
are approximately dipolar). This means we can select anomalies exhibitingorientations that allow us to test our predictions. Specifically, we can lookfor intra-swirl dark lanes at an anomaly with approximately vertical magneti10

ab

Vertical Dipole Field

40

35

Altitude above source (km)Altitude above source (km)

35
30
25
20
15
10
5
0
−30

−20

−10

0

10

Horizontal position (km)

20

15
10

−20

d

−10

0

10

Horizontal position (km)20

30

Altitude = 25 km

4

Altitude = 25 km
−20

−10

0

10

80

20

2
30

Altitude = 10 km

60
40
20
−20

−10

0

5 km

10

600

20

30

Horizontal field strength (nT)

Horizontal field strength (nT)

20

6

2

Altitude = 5 km

0
−30

−20

−10

0

10

100

30

Altitude = 10 km

0
−30

−20

−10

0

10

800

20

30

Altitude = 5 km

400

200

20

50

600

400

0
−30

25

0
−30

30

4

0
−30

30

5

c
0
−30

Horizontal Dipole Field

40

5/√2 km

200
−20

−10

0

10

Horizontal position (km)

20

30

0
−30

−20

−10

0

10

Horizontal position (km)20

30

Figure 2.1: Upper panels: magnetic field lines due to a single dipole located at theorigin and oriented either vertically (a), or horizontally (b), in the plane of the page;
Lower panels: profiles of the horizontal component of the magnetic field shown at
the various altitudes represented by dashed lines in the upper panels. In both cases,
the dipolar source has a magnetic moment of 1012 Am2 .

11zation because, in that case, field lines should appear vertical at any altitude
over the dark lane (Figure 2.1, panels a and c). Likewise, we can look for a correlation between high albedo and strongly horizontal fields at anomalies withapproximately horizontal magnetization because field lines at such anomalies
should appear horizontal over the brightest areas regardless of observation altitude (Figure 2.1, panels b and d). We find that the Airy anomaly is similar
to the case illustrated in the left panels of Figure 2.1 while the Reiner Gamma
anomaly is similar to the case illustrated in the right panels of Figure 2.1.

2.2

Data Processing

2.2.1

Lunar Prospector Magnetometer Data

Our study uses Lunar Prospector 3-axis magnetometer measurements (level1B data) obtained from NASA’s Planetary Data System (ppi.pds.nasa.gov).
Two distinct approaches are available for producing maps of crustal magnetic
fields from orbital measurements. One involves direct mapping using data selected from a sequence of orbits that cross the region of interest (Hood et al.,
1981, 2001), resulting in maps of the field at the spacecraft altitude, while
the other involves developing a global spherical harmonic model of the scalarpotential field (see Appendix B) and using it to produce maps at arbitrary altitudes (Purucker and Nicholas, 2010). The spherical harmonic models are well
suited to global mapping and have the advantages of automatically accounting
for varying spacecraft altitudes and guaranteeing that the mapped field is apotential field. To date, the best available spherical harmonic models extend
to degree 170, corresponding to a horizontal wavelength of ⇠ 64 km (Purucker
12and Nicholas, 2010). However, because we are interested in small-scale swirl
morphology, we instead employ a direct mapping technique similar to thatof Hood et al. (2001), allowing us to map features at wavelengths comparable with the spacecraft altitude (as low as ⇠ 18 km). As discussed by Hood(2011), the direct mapping approach is suitable for regional mapping and has
several advantages including that it allows us to perform model fitting directly
to the minimally processed magnetometer measurements rather than to a derived field model. While we focus here on direct mapping, we did compare
our maps against maps we derived from the Purucker and Nicholas (2010)
model coeﬃcients. We found that in the latter, crustal anomaly fields appear
morphologically similar to the same features in our maps, but are typicallybroader in horizontal extent and exhibit lower field strength.
Each of our maps incorporates data from a sequence of consecutive orbitsthat cross the anomaly of interest. The polar orbiting spacecraft obtained
global coverage twice each lunation, with consecutive orbits spaced ⇠ 1 degree
in longitude (⇠ 30 km at the equator) and with along-track measurements
recorded at 5-second intervals, corresponding to ⇠ 8 km spacing in latitude.
In order to capture the undisturbed signal of the crustal magnetic anomalies,
avoiding times when the field is distorted by the impinging solar wind (Kurata
et al., 2005; Purucker , 2008), we use only those magnetometer measurements
collected when the Moon was protected from solar wind while passing throughthe Earth’s magnetotail (excluding times when the field was disturbed by theplasma sheet) or when the Moon was in the solar wind but the spacecraft wasin the lunar wake (on the dark side of the Moon and away from the terminatorby at least 20°). We examined data from the lowest (< 50 km) altitude phaseof the mission, between February and July 1999 (the final six months of the13Day(s) of
1999
35,36,37
49,50
63,64
76,77,78
90,91
104,105
117,118
131,132
145,146
158,159
172,173
186,187
199,200
131,132
91,11849,50,77,104

Lunar Prospector Data for Airy
Solar wind Minimum
Mean
Number of
condition solar zenith
altitudeMeasureangle (°)
(km)
ments
wind
56
27.5
734
wake
131
18.7
772
wind
29
28.3
759
wake
150
18.5
615
tail
9
28.3
610
wake
151
18
678
tail
26
28.3
776
wake
134
34.4
782
wind
52
38.5
838
wake
111
34.3
627
wind
77
38.6
824
wind
88
35.9
660
wind
102
37
785
wake
136
34.4
48
tail
15
28.4
82
wake
134
18.5
130

Orbit
segments
15
15
14
12
11
13
14
14
15
12
15
12
14
4
7
12

Table 2.1: Lunar Prospector data used for Airy. The upper part of the table lists the
13 orbit sequences (each consisting of segments of between 11 and 15 consecutive
orbits) we examined in a 15° × 15° window centered on the anomaly. The solar
wind condition column indicates whether the spacecraft was exposed directly to thesolar wind (wind), or protected from it by being either in the lunar wake (wake)or the Earth’s magnetotail (tail). The last three rows of the table show the 260
measurements retained after discarding data collected outside of wake and tail times,
after combining like-altitude orbit sequences, and after cropping to a 3.25° × 4° study
area.

mission). The upper portions of Tables 2.1 and 2.2 list the 13 orbit sequences
we examined for each of the Airy and Reiner Gamma anomalies, respectively.
In each case, of the 13 orbit sequences we examined, 5 took place in the lunar
wake and 2 in the Earth’s magnetotail.
We examined orbit segments spanning 15° of latitude approximately centered on each anomaly. Having eliminated measurements taken outside ofwake and tail times, we assume that any remaining external fields are steady14Day(s) of
1999
40,41
54,55
67,68
81,82
95,96
108,109
122,123
136,137
149,150
163,164
177,178
190,191
204,205150
136,163,164
54,81,82,109

Lunar Prospector Data for Reiner Gamma
Solar wind Minimum
Mean
Number of
condition solar zenith
altitudeMeasureangle (°)
(km)
ments
wind
51
27.7
820
wake
139
19
769
wind
24
28.5
832
wake
162
18.6
795
wind
3
28.6
834
wake
158
18.3
634
wind
29
28.5
833
wake
135
34.8
703
tail
56
39.3
817
wind + wake
109
34.5
798
wind + tail
82
39.4
842
wind
83
36.3
838
wind
107
37.7
769
tail
57
39.7
48
wake
110
34.3
79
wake
141
18.3
118

Orbit
segments
15
14
15
15
16
12
15
13
15
15
15
15
14
4
7
10

Table 2.2: Lunar Prospector data used for Reiner Gamma. The upper part of thetable lists the 13 orbit sequences (each consisting of segments of between 12 and 16
consecutive orbits) we examined in a 15° × 15° window centered on the anomaly. The
solar wind condition column indicates whether the spacecraft was exposed directly
to the solar wind (wind), or protected from it by being either in the lunar wake
(wake) or the Earth’s magnetotail (tail). The last three rows of the table show the
245 measurements retained after discarding data collected outside of wake and tailtimes, after combining like-altitude orbit sequences, and after cropping to a 3.25° ×
4° study area.

15

over these 15° of latitude and that fields from isolated crustal sources span nomore than a few degrees. After subtracting the mean background field, we
assume the remaining fields are due to crustal sources. We then convert the
spacecraft position into selenographic spherical coordinates, taking the Lunar
Orbiter Laser Altimeter reference radius of 1737.4 km (Smith et al., 2010) as
zero altitude, and transform the magnetometer measurements into local east,
north, and radial components. We next combine data from consecutive orbits
to produce magnetic field maps by fitting regular square meshes to the datausing Delaunay triangulation. In order to capture the structure of the signal,
the grid cells must be no larger than half the spacecraft altitude, which canbe as low as ⇠ 18 km. We therefore use 0.25° × 0.25° (about 7.6 km × 7.6 km
at the equator) as the grid spacing when fitting to the Lunar Prospector magnetometer (LP MAG) data. This grid resolution is also finer than the spacing
between the observations (⇠ 8 km in latitude, up to ⇠ 30 km in longitude)
meaning that no observed signal variations are lost in the gridding process.
Using still finer resolution has no eﬀect on the resulting linearly interpolated
surface and is undesired as it increases computation time unnecessarily. Conversely, using a larger (more coarse) grid spacing leads to undesired smoothingin latitude, and as cell size approaches 1°, undesired smoothing occurs in longitude as well. Despite our eﬀorts to remove external fields and avoid times withsignificant field distortion, portions of some of the orbit sequences still appear
to be contaminated by transient signals and are therefore discarded. In orderto obtain maps with improved spatial coverage, we combine measurementsfrom diﬀerent orbit sequences if the observation altitudes are suﬃciently similar (diﬀering by less than ⇠ 1 km). Since the spacecraft altitude varies only
slightly (< 1 km) over the scale of our maps (> 120 km), we do not perform16

upward or downward continuation of the signal to some fixed altitude. Instead, our maps represent the magnetic field at the slowly varying spacecraft
altitude. Finally, we define study areas within the boundaries of the processed
data: the Airy study area is illustrated in the left panels of Figure 2.2 while
the Reiner Gamma study area is illustrated in the right panels of Figure 2.2.
As summarized in the lower portions of Table 2.1 and Table 2.2, this resultsin Airy study area maps at altitudes of ⇠ 18 km, ⇠ 28 km, and ⇠ 34 km, using
a total of 260 LP MAG measurements, and Reiner Gamma study area maps
at altitudes of ⇠ 18 km, ⇠ 34 km, and ⇠ 40 km, using a total of 245 LP MAG
measurements (only the lowest altitude maps are presented here).

2.2.2

Clementine Reflectance Mosaics

Our study uses Version 2 Clementine 750 nm reflectance mosaics producedby the USGS (United States Geological Survey) map-a-planet service (mapaplanet.org). Compared with the Version 1 (V1) maps, the Version 2 (V2)maps use a newer geodetic control network (Archinal et al., 2006), refining
the horizontal registration by as much as 10 km in some locations. We have
manually adjusted the V2 reflectance values to match those of the V1 maps,
which have been more carefully controlled to represent true reflectance.

2.3

Observations

We focus our discussion here on two specific examples: Airy and Reiner
Gamma. The Airy anomaly exhibits a magnetic field resembling the caseillustrated in Figure 2.1a while the Reiner Gamma anomaly field resembles
the case illustrated in Figure 2.1b. The distinct orientations of these two17

anomalies reveal diﬀerent aspects of the relationship between magnetic field
direction and swirl morphology.

2.3.1

Airy

First described by Blewett et al. (2007), the Airy swirl is found near Airycrater in the lunar nearside highlands. At the Lunar Prospector spacecraft altitude, the magnetic field lines tend to point inward toward the middle of theanomaly, becoming increasingly vertical and downward pointing at the center(Figure 2.2a), consistent with a source magnetization that is pointed mainly
downward. This inferred magnetization direction is based on examination of
the vector components and is supported by the models described in section 2.4.
Here, contour maps of the total magnetic field strength illustrate only thatthe albedo anomaly is approximately centered on the magnetic anomaly (Figure 2.2b). However, maps of the horizontal field alone (Figure 2.2c) reveal
structure that is more closely related to albedo morphology. For example, thedark lane through the center of the albedo anomaly forms an approximate
plane of symmetry in the observed horizontal magnetic field map. Even morestriking is the alignment between the dark lane in the center of the anomaly
and the line representing zero east-west magnetic field strength (dashed whiteline in Figure 2.3). The brightest parts of the swirl are organized into tworoughly parallel lobes on either side of the dark lane. The center-to-center
horizontal distance between the lobes is approximately 8-10 km. Peaks in thehorizontal field strength are also organized into two lobes on either side of thedark lane but with a peak-to-peak separation of roughly 30 km (Figure 2.2c).
As illustrated in Figure 2.1c, a closer spacing of horizontal field strength peaks

18

Airy
−16.5

9

a

8.5

Reiner Gammad

Latitude (degrees N)

Latitude (degrees N)

−17

−17.5

−18

−18.5

8

7.5

7

6.5

−19

10km
−19.5

−16.5

1.5

6

2

2.5

3

3.5

Longitude (degrees E)

4

4.5

5

10km
299.5

9

b

8.5

300

300.5

301

301.5

Longitude (degrees E)

302

302.5

303

302

302.5

303

302

302.5

303

e

Latitude (degrees N)

Latitude (degrees N)

−17

−17.5

−18

−18.5

−19

−19.5

−16.5

10
2nT 4 6 8
1.5

2

2.5

22
20
18
16
12 14
3

8

7.5

44

7

36
6.5

20

6

3.5

Longitude (degrees E)

4

4.5

4nT
299.5

5
9

c

8.5

300

28

12
300.5

301

301.5

Longitude (degrees E)

fLatitude (degrees N)

Latitude (degrees N)

−17

−17.5

−18

3
−18.5

−19

−19.5

3nT
1.5

2

6

12

6

12

9

9

9

6

8

7.5

24
16

7

8nT

6.5

3

32 40

6
2.5

3

3.5

Longitude (degrees E)

4

4.5

299.5

5

300

300.5

301

301.5

Longitude (degrees E)

Figure 2.2: Magnetic field maps derived from Lunar Prospector magnetometer dataover Clementine 750 nm reflectance maps at the Airy (left panels) and Reiner Gamma(right panels) anomalies. Panels (a) and (d) show the direction of field lines at the
spacecraft altitude (⇠ 18 km in both cases), with arrow lengths showing relative
horizontal field strength. Panels (b) and (e) show contours of the total magnetic
field strength and panels (c) and (f) show contours of the horizontal component
alone. The Airy maps shown here are derived from LP MAG data collected on days
49-50, 77, and 104 of 1999 while the Reiner Gamma maps are derived from datacollected on days 54, 81-82, and 109 of 1999.

19

−17

Latitude (degrees N)

0

+2

−16.5

-2

+4

-4

+6

-6
-8
-10

+8

−17.5

+10nT

-12nT
−18

−18.5

−19

−19.5

1.5

2

2.5

3

3.5

Longitude (degrees E)

4

4.5

5

Figure 2.3: East-west component of magnetic field as measured by Lunar Prospector
over Clementine 750 nm reflectance at Airy. Contour line colors and labels indicate
magnitude of eastward field strength, with warm colors indicating east-pointing fields
and cool colors indicating west-pointing fields. The dashed white line indicates zero
field strength in the east-west direction. The map is derived from LP MAG datacollected on days 49-50, 77, and 104 of 1999.

is expected at lower altitudes, potentially allowing for better alignment withthe swirl’s bright lobes.

2.3.2

Reiner Gamma

Reiner Gamma, the type example for lunar swirls, is found in Oceanus Procellarum near the western limb of the lunar nearside. At the Lunar Prospector
spacecraft altitude, the magnetic field lines are south pointing (Figure 2.2d),
consistent with a source magnetization that is mainly horizontal and northpointing. This inferred magnetization direction is based on examination of
the vector components, is supported by the models described in section 2.4,
and is in agreement with Kurata et al. (2005). Contour maps of total magnetic
field strength, seen here (Figure 2.2e) and elsewhere (e.g. Hood et al. (2001)),
illustrate only that the albedo anomaly is approximately co-located with themagnetic anomaly; there is no clear relationship between total field strength
20and the morphology of the albedo anomaly. However, contour maps of thehorizontal component alone (Figure 2.2f) show that regions of high horizontal field strength correspond well with the bright swirl. Based on SELENE
magnetometer data, Shibuya et al. (2010) have reported similar findings forReiner Gamma and other anomalies. Patches of high albedo tend not to occur
where fields are weak, or where fields are strong but lack a large horizontal
component. This observation is consistent with the hypothesis that darkening
due to solar wind ion bombardment will be minimized where magnetic fields
are strongly horizontal.
We tested the degree to which magnetic field direction is related to reflectance by measuring the correlation between reflectance and the angle thefield makes with the vertical. Figure 2.4 plots reflectance versus magnetic field
direction using data points gathered from all six retained sets of LP MAG measurements. Since the field strength varies with altitude, values were normalized
to the maximum total field strength for each map before being combined into
the set of 245 data points used to generate the scatter plot. Figure 2.5 illustrates the locations of the data points and corresponding reflectance values
used to make the scatter plot. The reflectance map in Figure 2.5 has beensmoothed by a 0.25° × 0.25° window moving average filter in order to reduce
its eﬀective resolution to be comparable to the LP MAG data resolution. This
removes the high frequency component of the albedo signal that could not possibly be captured in the lower frequency magnetic field data. The data points
in Figure 2.4 are color-coded according to normalized total field strength. Coolcolors represent weak fields and warm colors represent strong fields. A blue
regression line is fit to the weakest third of the data points and illustrates that,as expected for weak fields, reflectance values are low regardless of magnetic
21

1

0.14

0.9

0.135

0.8

0.13

0.7

Normalized Total Field Strength

750nm Reflectance

0.145

0.125

0.6

0.12

0.5

0.115

0.4

0.11

0.3

0.105

0.2

0.1

0.1

0.095
10

20

30

40

50

60

Angle from vertical (degrees)

70

80

0

Figure 2.4: Correlation between magnetic field direction and reflectance at Reiner
Gamma. ‘Angle from vertical’ is 90° where field lines are horizontal. Data points arecolor-coded by total field strength with cool colors representing weak fields and warm
colors representing strong fields. Blue and red regression lines are fit to the weakest
and strongest thirds of the data points, respectively, demonstrating that magnetic
field direction becomes increasingly important with increasing field strength.

field direction. A red regression line is fit to the strongest third of the data
points and illustrates that, for stronger fields, reflectance values are low where
field lines are vertical and high where field lines are horizontal. Dashed lines
are used to show 95% confidence intervals on the estimated slopes, indicating that the observed diﬀerences in slope are statistically significant. These
trends are precisely what we would expect if the albedo anomalies owe theirbrightness to magnetic deflection of solar wind.

2.3.3

Discussion

The examples of Airy and Reiner Gamma illustrate distinct aspects of thesolar wind deflection phenomenon: the Airy case shows that dark lanes may
be associated with vertical magnetic fields (Figures 2.2a–c and Figure 2.3),
while the Reiner Gamma example demonstrates that bright areas may be
associated with strongly horizontal fields (Figures 2.2d–f and Figure 2.4). In
22

9
0.145
0.14
0.135

8

0.13
7.5

750nm Reflectance

Latitude (degrees N)

8.5

0.125
0.12

7

0.115
6.5

0.11
0.105

6
299.5

300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

0.1

Figure 2.5: Locations of the 245 LP MAG data points (green circles) used to generate
Figure 2.4. Background is a Clementine 750 nm reflectance mosaic with resolution
reduced to be comparable with the LP MAG data resolution.

principle, however, both of these eﬀects should be present at both locations.
That is, we should expect to see alignment between the bright lobes of Airyand peaks in horizontal magnetic field strength and we should expect to see
vertical field lines over the dark lanes of Reiner Gamma. But as discussed in
section 2.1.2, we expect the near-surface field patterns to dominate solar wind
deflection and we do not yet have observations at suﬃciently low altitudes
to map the near-surface field directly. However, we can use source models todetermine whether the near surface field pattern we predict (horizontal overbright areas and vertical over dark lanes) is consistent with the observational
constraints we do have.

2.4

Source Modeling

Using a combination of techniques, including a genetic search algorithm, wedeveloped subsurface magnetization source models to gain insight into the possible distribution of source material at the Airy and Reiner Gamma anomalies.
We use these models to support our interpretation of the observations discussed
23

above and to allow us to further test the plausibility of our hypothesis—that
areas of high albedo should coincide with dominantly horizontal near-surface
magnetic fields.

2.4.1

Lunar Prospector Data Fitting

As a first step, we use the Lunar Prospector magnetometer (LP MAG) measurements to determine the most probable characteristics of the magnetic
source (the approximate distribution of magnetic material as well as the magnetic moment magnitude and orientation). The results are spatially coarse as
the data are insensitive to variations on scales smaller than the observation
altitudes. The results are also inherently non-unique as trades can be madebetween depth and magnetic moment magnitude, for example. Previous studies have modeled the source of the central part of the Reiner Gamma anomaly
as a single dipole (Kurata et al., 2005) or as a grid of dipoles (Nicholas et al.,2007) using inversion techniques (Von Frese et al., 1981; Purucker et al., 1996;
Dyment and Arkani-Hamed , 1998; Parker , 2003) to obtain models that best
fit the LP MAG data in a least-squares sense. Here, we employ an alternative
approach involving iterative forward modeling to identify the characteristics ofthe best-fitting solutions. Our approach demonstrates that a range of diﬀerent solutions can deliver similarly good results. Throughout our modeling, we
make the simplifying assumption that the source material is coherently magnetized in a single direction; much larger remanent magnetizations would be
required to produce the observed fields if the source materials were not unidirectionally magnetized, and the inferred values are already large (as discussed
in section 2.5.1).

24

2.4.1.1

Airy Single Dipole Model

We begin with a single dipole model to find the magnetic moment directionthat best fits, in a least squares sense, the 260 LP MAG measurements of
the Airy study area. We define the ‘eﬀective error’ at each data point as the
maximum of the three vector component residuals. The least squares solution
thus minimizes residuals in all three vector components simultaneously. We
vary the burial depth from 0 to 20 km in 1-km increments and the magnetic
moment magnitude from 0 to 8 ⇥ 1012 Am2 in increments of 5 ⇥ 1011 Am2 . We
vary the inclination (angle positive downward from the horizontal) from -90° to+90° and the declination (angle positive clockwise from north) through 360°,
each with 1° resolution. We find the measurements are well accommodated by adipole buried between 10 and 20 km below the surface, with magnetic moment
magnitude between 4⇥1012 Am2 and 7⇥1012 Am2 , inclination between 79° and
81° and declination between -14° and +29° (i.e., pointed steeply downward andslightly to the north). Since this method delivers equally good solutions overa wide range of depths and magnetic moment magnitudes, we do not attempt
to constrain depth and magnitude at this stage. In any case, because we areassuming a single dipole, the model fit overestimates the true source depth.
2.4.1.2

Reiner Gamma Single Dipole Model

For Reiner Gamma, we again begin with a single dipole oriented so that it bestfits, in a least squares sense, the 245 LP MAG measurements at the Reiner
Gamma study area. Using the ‘eﬀective’ error metric described for the Airymodel above, we search for the least squares best fit by varying the burial
depth from 0 to 20 km in 1-km increments, the magnetic moment magnitude

25

from 0 to 2 ⇥ 1013 Am2 in increments of 1012 Am2 , the inclination from -90° to+90° and the declination through 360°, with 1° resolution for both inclinationand declination. We find that the measurements are well accommodated by adipole that is buried between 8 and 15 km below the surface with magnetic
moment between 9 ⇥ 1012 Am2 and 16 ⇥ 1012 Am2 , inclination between -6°
and +5° and declination between -18° and -3° (i.e., lying approximately in theplane of the lunar surface and pointing slightly west of north). For comparison,
Kurata et al. (2005) modeled this part of the anomaly as a single dipole buried
11.1 km below the surface with magnetic moment 11.3 ⇥ 1012 Am2 , inclination+1.3° and declination -11°. Once again, depth and magnitude are not wellconstrained by the data and are overestimated by this single dipole model.
2.4.1.3

Genetic Search Algorithm

Next, for each of the two anomalies, we expand the parameter space by replacing the single dipole with a grid of dipoles, separated by 0.25° in both
latitude and longitude (a distance comparable with half the minimum spacecraft altitude), resulting in a total of 208 dipoles. We allow each of the 208
dipoles’ magnitudes to vary independently. We allow the burial depth of thegrid to vary as a whole but burial depth does not vary between dipoles and
we do not allow the dipoles to be above the lunar surface. Likewise, we allow
the magnetization direction to vary as a whole but direction does not vary
between dipoles. This opens a parameter space that is too large to explore
completely, even with coarse resolution. Instead, we employ a heuristic searchtechnique known as a “genetic algorithm” that iteratively adjusts each of the
211 independent parameters (208 independent dipole moment magnitudes plus
depth of the grid, inclination and declination), gradually progressing toward26

improving least-squares solutions.
Genetic algorithms have a diverse range of applications (Holland , 1992;
Goldberg, 1989) including large-scale optimization problems such as ours: tofind a source model (defined by 211 parameters) that produces minimal error
between predicted and observed magnetic fields. Genetic algorithms employ
concepts borrowed from gene-centered biological evolution (Dawkins, 1976) inorder to iteratively (i.e., over many generations) progress towards solutions
with increasing degrees of ‘fitness’. In our case, this ‘fitness’ is a measure of
how well the magnetic field produced by the source model matches the Lunar
Prospector magnetometer (LP MAG) observations.
Our algorithm begins by generating a population of individual source models with randomly distributed characteristics. The characteristics of each individual source model are defined by its ‘genes’, with each individual having
all 211 distinct genes (208 dipole moment magnitudes plus depth of the grid,
inclination and declination). At each iteration, or generation, the individual
members of the population are evaluated according to how well they predictthe magnetic field at each of the LP MAG observation points. We computethe error between the prediction and the observation using the ‘eﬀective error’
metric described above (maximum of the three vector component residuals).
We then rank the individuals from lowest to highest sum of squared eﬀective
errors. The individuals with the lowest errors are then selected as the ‘parents’ of the next generation. The next iteration begins by generating a new
population to replace the previous generation. Each individual in the newgeneration is formed by setting each of its 211 genes equal to the corresponding gene from one of its parents in a process known as crossover (analogous
to chromosomal crossover). For example, an individual’s gene for inclination
27

will match the inclination gene of one of its parents (chosen at random fromamong the parents). Because only the best-performing individuals contribute
their genes to the next generation, individuals of each new generation tend
to inherit the best characteristics of the previous generation. We then apply
random mutations (i.e., randomly adjust several isolated genes). This step ensures variation in the population and allows for the possibility of introducing
advantageous genes not possessed by the parents. Finally, members of the newgeneration are evaluated and their best performers are selected as the parentsof the next generation. With each generation, genes that result in large model
errors tend to be discarded while genes that result in smaller model errors tend
to be retained. Inheritance and selection lead the gene pool to be increasingly
rich in genes that form good source models while crossover and mutations ensure variation, allowing for innovations that can lead oﬀspring to outperformtheir parents. The result is that the population progresses gradually towards
optimality in terms of minimum model error.
While this explanation captures the essence of the algorithm, our implementation includes additional details such as gene mutation rate (the probability that any given gene will mutate) and how “random” mutations are
distributed (i.e., it is not useful to have sudden changes in burial depth on theorder of 1000 km, for example, so mutations must be limited to reasonableadjustments according to some distribution). We experimented with various
population sizes and numbers of parents selected at each generation (there
must be at least two parents but more than two is also allowed). For computational eﬃciency, a small population size is preferred. The algorithm wasfound to operate eﬀectively with a population of just 10 individuals in which
the top-performing 3 contribute their genes to the next generation. We allow
28

mutations to occur with a probability of 0.1 and when they do occur, thegene is adjusted from its current value according to a normal distribution withsome pre-defined standard deviation (which depends on whether the gene hasunits of kilometers, degrees or Am2 ). We also apply smaller mutations with ahigher probability in order to ensure a measurable degree of variation in the
population. The initial population is generated based on a seed model; eachindividual in the initial population is formed by applying normally distributed
adjustments to each of the genes in the seed model. While the choice of seedmodel and parameter settings influences the eﬃciency of the algorithm (i.e.,
how quickly it reaches a good solution), the end results tend to be very similar
across a wide range of parameters and initial conditions. The algorithm wasallowed to iterate through 600 generations (Figure 2.6). Based on the trends
in error evolution (Figure 2.7), we would not expect additional iterations toyield substantial improvements.
2.4.1.4

Airy Dipole Grid Model

For the Airy dipole grid model, we initiate the genetic search algorithm with adownward pointing magnetization in accord with the best fitting single dipole(80° inclination, 8° declination) and an arbitrary burial depth of 10 km. Again,
the genetic search algorithm attempts to find a solution that minimizes thesum of squares over the ‘eﬀective’ error at each of the 260 LP MAG data
points. The dipoles in the final model grid have inclination 79°, declination
20° and a total magnetic moment of 4.7 ⇥ 1012 Am2 distributed as shown in
Figure 2.8. The final grid has a burial depth of ⇠ 650 meters, but again, depth
is not well constrained by the LP MAG data as it can be traded against magnetic moment magnitude and the lateral extent of the source; an equally good29

Figure 2.6: Snapshots at generations 1, 100, 300, and 600 during the evolution of thedipole grid source models for Airy (left panels) and Reiner Gamma (right panels).
Each panel shows the distribution of magnetic moments from the best-performing
model of the specified generation.

30

Airy Model Error

9000

4

8000

3

6000

2.5

5000

Error

Error

Reiner Gamma Model Error

3.5

7000

4000

2

1.5

3000

1

2000

0.5

1000
0

4

x 10

0

100

200

300

Iteration

400

500

600

0

0

100

200

300

Iteration

400

500

600

Figure 2.7: Evolution of dipole grid source model error over 600 generations for Airy(left panel) and Reiner Gamma (right panel). The error shown is the sum of squared‘eﬀective’ errors (as described in section 2.4.1.1) over all LP MAG data points forthe best-performing model in each generation.

fit could be obtained with a source that is deeper, stronger and more horizontally concentrated. Although the model is spatially coarse, the distributionof magnetic moments appears to be consistent with a source structure that iscentered on the anomaly and elongated in a north-south direction. Figure 2.9
shows that the resulting model field compares very well with the observed
field. The magnetization direction is close to that of the hypothetical field
illustrated in Figure 2.1a, indicating that the intra-swirl dark lane may wellcoincide with vertical field lines at the surface.
2.4.1.5

Reiner Gamma Dipole Grid ModelAs with the Airy model, we replace the single dipole model with a grid of
208 dipoles and use the genetic algorithm to find a solution that best fits the
245 LP MAG data points. The final model dipole grid is magnetized withinclination +4° and declination -12° at a depth of 1.6 km (but again, depth is
not well constrained by the LP MAG data), and has total magnetic moment

31

Figure 2.8: The Airy source model obtained via the genetic search algorithm described in section 2.4.1.3. In the left panel, each square represents a single dipole
covering 0.25° × 0.25° (roughly 5.5 ⇥ 107 m2 ) with the color indicating each dipole’s
total magnetic moment (typical values in the center are ⇠ 1.25 ⇥ 1011 Am2 per
dipole). The right panel shows the same information as contours over Clementine
albedo, suggesting approximate agreement between the source structure’s longitudinal axis and that of the swirl’s dark lane.

⇠ 1.4 ⇥ 1013 Am2 distributed as shown in Figure 2.10. The source model isspatially coarse but is consistent with a source structure that is elongatedin the east-west direction and is most intense in the brightest part of thealbedo anomaly. Figure 2.11 shows that the resulting model field compares
well with the observed field. For comparison, Nicholas et al. (2007) modeled
the Reiner Gamma source as a grid of dipoles separated by 0.1° and placedat the surface, coincident with the albedo anomaly. Those authors assumed anorthward magnetization and solved for the magnitude and sign at each of thedipoles, obtaining the distribution of magnetization illustrated in Figure 2.12.
In spite of the very diﬀerent methods employed, the two results resemble one
another, predicting strong magnetization in the brightest parts of the anomaly.

32

4

10

−17

5

−17.5

0

−18

−5

−18.5
−19

−10

5

−17
−17.5

0

−18
−18.5

−5

−19
−19.5

2

3

4

Longitude (degrees E)

5

−16.5

−5

−17.5
−18

−10

−18.5

−15

−19
−19.5

−20
2

3

4

Longitude (degrees E)

5

12

−17

10

−17.5

8

−18

6

−18.5

4

−19
−19.5

2
2

3

4

Longitude (degrees E)

5

2

3

4

5

−16.5

10

−17

5

−17.5

0

−18

−5

−18.5
−19
−19.5

−10
2

3

4

5

−16.5
5

−17
−17.5

0

−18
−18.5

−5

−19
−19.5

2

3

4

5

−16.5

0

−17

−5

−17.5
−18

−10

−18.5

−15

−19
−19.5

−20
2

3

4

Longitude (degrees E)
14

−16.5

−19.5

Longitude (degrees E)
0

−17

5

−19

Longitude (degrees E)−16.5

Radial Field Strength (nT)

Latitude (degrees N)

5

Latitude (degrees N)

4

North Field Strength (nT)3

Latitude (degrees N)

2

10

−18.5

Longitude (degrees E)

Latitude (degrees N)

Latitude (degrees N)−16.5

Longitude (degrees E)

Latitude (degrees N)

5

Latitude (degrees N)

3

East Field Strength (nT)2Longitude (degrees E)

Total Field Strength (nT)

−19

−18

East Field Strength (nT)

5

15

−17.5

North Field Strength (nT)

10

−18.5

20

−17

Radial Field Strength (nT)

−18

−16.5

5

14

−16.5

12

−17

10

−17.5

8

−18

6

−18.5

4

−19
−19.5

Horizontal Field Strength (nT)

15

−17.5

Latitude (degrees N)

20

−17

−19.5

Total Field Strength (nT)

−16.5

−19.5

Latitude (degrees N)

Model

Horizontal Field Strength (nT)

Latitude (degrees N)

LP MAG

2
2

3

4

Longitude (degrees E)

5

Figure 2.9: Comparison of the observed field at ⇠ 18 km over Airy (left panels;
derived from LP MAG data collected on days 49-50, 77, and 104 of 1999) and themodel field (right panels; obtained by fitting a grid of dipoles to the LP MAG datausing the genetic search algorithm).

33

Figure 2.10: The Reiner Gamma source model obtained via the genetic search algorithm described in section 2.4.1.3. In the left panel, each square represents a single
dipole covering 0.25° × 0.25° (roughly 5.7 ⇥ 107 m2 ) with the color indicating each
dipole’s total magnetic moment (typical values in the center are ⇠ 3.4⇥1011 Am2 per
dipole). The right panel shows the same information as contours over Clementine
albedo.

2.4.2

Albedo Pattern Matching

The results of the genetic search algorithm suggest that the distribution of
the underlying source material coincides roughly with the shape of the albedo
anomaly: a north-south distribution for Airy and an east-west distribution forReiner Gamma. Here, we refine our source models by applying the constraintthat the near-surface field be structured according to our hypothesis: strongly
horizontal over the brightest parts of swirls and vertical in the intra-swirl darklanes. This allows for greatly improved spatial resolution since Clementinereflectance mosaics are available at 256 pixels/degree whereas LP MAG data
are limited to ⇠ 1 4 pixels/degree. If our hypothesis is correct, source models
constrained by the albedo pattern should produce fields that are consistent
with the LP MAG observations made at higher altitudes. Below, we showthat even simple source models are suﬃcient to accomplish this, suggestingthat our models are highly plausible.

34

LP MAG

Model

0

7

−5

6.5

−10

6
302

Longitude (degrees E)

Latitude (degrees N)

9

−10

7.5
7

−20

6.5

−30

6
300

301

302

Longitude (degrees E)

303

Latitude (degrees N)

9
8

10

7.5

0

7

−10

6.5

−20
302

Longitude (degrees E)

303

9
8

30

7.5

20

7
6.5

10

6
300

301

302

Longitude (degrees E)

303

Total Field Strength (nT)

−5

6.5

−10

6
301

302

303
10

8.5

0

8

−10

7.5
7

−20

6.5

−30

6
301

302

303

9

−40

30

8.5

20

8

10

7.5

0

7

−10

6.5

−20
−30

6
300

301

302

Longitude (degrees E)

40

8.5

0

7

300

−30

6

5

7.5

Longitude (degrees E)

20

301

10

8

9

−40

30

8.5

300

15

Longitude (degrees E)

0

8

303

9

300
10

8.5

302

8.5

303

North Field Strength (nT)301

Latitude (degrees N)

5

East Field Strength (nT)10

Radial Field Strength (nT)

Latitude (degrees N)

15

7.5

301

Longitude (degrees E)

9300

Latitude (degrees N)

300

8.5
8

10

6

303

Latitude (degrees N)

302

Latitude (degrees N)

301

Longitude (degrees E)

303

9Latitude (degrees N)

300

20

7
6.5

East Field Strength (nT)

10

6

30

7.5

North Field Strength (nT)

20

7
6.5

8

Radial Field Strength (nT)

30

7.5

40

8.5

Horizontal Field Strength (nT)89Latitude (degrees N)

Total Field Strength (nT)

40

8.5

Horizontal Field Strength (nT)

Latitude (degrees N)

9

40

8.5
8

30

7.5

20

7
6.5

10

6
300

301

302

Longitude (degrees E)

303

Figure 2.11: Comparison of the observed field at ⇠ 18 km over Reiner Gamma (leftpanels; derived from LP MAG data collected on days 54, 81-82, and 109 of 1999)
and the model field (right panels; obtained by fitting a grid of dipoles to the LP
MAG data using the genetic search algorithm).

35

9

Latitude (degrees N)

0.14

Clementine V1 750nm Reflectance

0.145

8.5

0.135

8

0.13
7.5

0.125
0.12

7

0.115
6.5

0.11

0.105

6
299.5

300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

0.1

Figure 2.12: Distribution of source magnetization at Reiner Gamma as estimatedby Nicholas et al. (2007) assuming the source layer is 40 km thick. The contour
interval is 0.1 A/m with warm colors positive, cool colors negative and the dashedwhite contour representing zero magnetization. The contours are superimposed over
Clementine V1 750 nm reflectance. Note: our study instead uses Clementine V2
mosaics, which diﬀer in horizontal registration by ⇠ 7 km in the vicinity of Reiner
Gamma.

2.4.2.1

Airy Model

Given the downward-pointing magnetization at Airy (section 2.4.1.4) and theassociated field pattern (Figure 2.1c), the morphology of the Airy swirl suggests a roughly linear source structure following the dark lane. In order tobe consistent with our hypothesis, the vertically magnetized source material
cannot substantially underlie the bright lobes since this would mean verticalfield lines in bright areas. We therefore replace the grid of dipoles described in
section 2.4.1.4 with an array of 32 dipoles arranged as a two-segmented linear
feature following the dark lane (Figure 2.13, left panel). In cross-section, thefield pattern from such a line source resembles the pattern from a single dipole(illustrated in panels a and c of Figure 2.1). However, for a linear source,
the peaks in the horizontal field profile are shifted outward slightly such thatthe peak-to-peak separation is 2/√3 times the altitude above the source rather
than being equal to the altitude above the source (details are given in Ap36

pendix B). If we assume the two bright lobes at Airy (separated by ⇠ 9 km)
are aligned with peaks in the surface horizontal magnetic field profile, then thedepth of this line source must be ⇠ 7.8 km. Because we assume an infinitely
narrow source, this represents the maximum possible source depth consistent
with the hypothesis that the bright lobes align with peaks in the horizontal
magnetic field profile. The LP MAG data are too spatially coarse to constrain
the distribution of magnetic moments within the source structure so we assign
equal strength to each of the 32 dipoles in the array. With the depth fixed
and the dipole moment magnitudes set to be equal, the parameter space isreduced to three dimensions (total magnitude, inclination, and declination),allowing for a complete search (varying magnitude in increments of 1011 Am2and inclination and declination in 1° increments) to find the solution that bestfits, in a least squares sense, the LP MAG data. The resulting best-fit solution
has the dipoles pointed steeply downward (inclination 75°, declination -20°)
with total magnetic moment 3.8 ⇥ 1012 Am2 (Table 2.3). This source model
produces a magnetic field pattern at the surface that mimics the morphologyof the albedo anomaly, with the field being strongly horizontal in high albedo
areas and vertical in the intra-swirl dark lane (Figure 2.13). Remarkably, thissame simple source model produces a field pattern at the spacecraft altitude
that substantially resembles the LP MAG observations (Figure 2.14).
2.4.2.2

Reiner Gamma Model

The Airy model (2.4.2.1) suggests that the dark lane is aligned with a cusp
between two peaks in the horizontal magnetic field profile. Based on the Airy
example, we may view the two dark lanes at Reiner Gamma as suggestiveof cusps in the near-surface horizontal magnetic field profile (see Harnett and37

750nm Reflectance

0.225

−17.5

0.22

−18

0.215
−18.5

0.21

−19

−19.5

180
−17

0.23

Latitude (degrees N)

Latitude (degrees N)

−17

2

2.5

3

3.5

Longitude (degrees E)

4

4.5

5

160
140

−17.5

120
100

−18

80
−18.5

60
40

−19

0.205

1.5

200

−16.5

0.235

−19.5

0.2

Horizontal Field Strength (nT)

−16.5

20
1.5

2

2.5

3

3.5

Longitude (degrees E)

4

4.5

5

Figure 2.13: Left panel: Airy study area showing Clementine albedo map with red
circles indicating the locations of the source model’s 32 dipoles; Right panel: resulting
horizontal magnetic field strength predicted at the surface. White arrows indicate
where field lines become vertical in the intra-swirl dark lane. Figure 2.14 compares
the model field predicted at the spacecraft altitude with the LP MAG observations.

Albedo Pattern-Based Model for Airy
Latitude (°N)

Longitude (°E)

Latitude (°N)

Longitude (°E)

-17.4892.826-17.906

3.322-17.5122.862-17.951

3.330-17.5342.899-17.9963.338-17.557

2.935-18.0413.347-17.579

2.972-18.0873.355-17.602

3.008-18.132

3.363-17.6243.044-18.177

3.372-17.647

3.081-18.2233.380-17.669

3.117-18.268

3.389-17.692

3.153-18.3133.397-17.7143.190-18.3593.405-17.7373.226-18.404

3.414-17.7593.263-18.4493.422-17.782

3.299-18.4943.430-17.815

3.305-18.5403.439-17.860

3.313-18.585

3.447

Table 2.3: Final dipole array model for Airy as described in section 2.4.2.1 (Figure 2.13). The model consists of 32 dipoles at the indicated latitudes and longitudes,
all buried 7.8 km below the surface. The total magnetic moment is 3.8 ⇥ 1012 Am2 .
All dipoles in the array have inclination 75°, measured positive downward from thehorizontal and declination -20°, measured positive clockwise from north.

38

4

10

−17

5

−17.5

0

−18

−5

−18.5
−19

−10

5

−17
−17.5

0

−18
−18.5

−5

−19
−19.5

2

3

4

Longitude (degrees E)

5

−16.5

−5

−17.5
−18

−10

−18.5

−15

−19
−19.5

−20
2

3

4

Longitude (degrees E)

5

12

−17

10

−17.5

8

−18

6

−18.5

4

−19
−19.5

2
2

3

4

Longitude (degrees E)

5

2

3

4

5

−16.5

10

−17

5

−17.5

0

−18

−5

−18.5
−19
−19.5

−10
2

3

4

5

−16.5
5

−17
−17.5

0

−18
−18.5

−5

−19
−19.5

2

3

4

5

−16.5

0

−17

−5

−17.5
−18

−10

−18.5

−15

−19
−19.5

−20
2

3

4

Longitude (degrees E)
14

−16.5

−19.5

Longitude (degrees E)
0

−17

5

−19

Longitude (degrees E)−16.5

Radial Field Strength (nT)

Latitude (degrees N)

5

Latitude (degrees N)

4

North Field Strength (nT)3

Latitude (degrees N)

2

10

−18.5

Longitude (degrees E)

Latitude (degrees N)

Latitude (degrees N)−16.5

Longitude (degrees E)

Latitude (degrees N)

5

Latitude (degrees N)

3

East Field Strength (nT)2Longitude (degrees E)

Total Field Strength (nT)

−19

−18

East Field Strength (nT)

5

15

−17.5

North Field Strength (nT)

10

−18.5

20

−17

Radial Field Strength (nT)

−18

−16.5

5

14

−16.5

12

−17

10

−17.5

8

−18

6

−18.5

4

−19
−19.5

Horizontal Field Strength (nT)

15

−17.5

Latitude (degrees N)

20

−17

−19.5

Total Field Strength (nT)

−16.5

−19.5

Latitude (degrees N)

Model

Horizontal Field Strength (nT)

Latitude (degrees N)

LP MAG

2
2

3

4

Longitude (degrees E)

5

Figure 2.14: Comparison of the observed field at ⇠ 18 km over Airy (left panels;
derived from LP MAG data collected on days 49-50, 77, and 104 of 1999) and themodel field (right panels; obtained as described in 2.4.2.1). Figure 2.13 shows thehorizontal component of the model field predicted at the surface.

39

Winglee (2003) for a related discussion). As illustrated in Figure 2.1d, ahorizontally magnetized source will produce strongly horizontal fields directly
above the source while producing cusps (where the horizontal field strength
drops to zero) on either side of the source. This suggests that the horizontallymagnetized source material underlies the brightest parts of the Reiner Gamma
albedo feature and that the dark lanes may be aligned with the cusps. The
appearance of the two dark lanes and the relatively bright material between
them may be explained by a superposition of two sources (see Figure 2.15 andcompare with Figure 2.1).
As with the case illustrated in Figure 2.1 panels b and d, the field illustrated in Figure 2.15 is strongly horizontal at any altitude over the twohorizontally magnetized sources. However, in this case there is an additional
region of elevated horizontal field strength where the side lobes interfere constructively between the two sources. Taking this pattern as a cue, we replace
the grid of dipoles described in section 2.4.1.5 with an array of 55 dipoles arranged as curvilinear structures beneath the three brightest parts of the swirl(Figure 2.16, left panel). In cross-section, the field pattern from the two approximately linear sources adjacent to the dark lanes resembles the patternfrom the two-dipole case illustrated in Figure 2.15, but the cusps in the horizontal magnetic field profile are shifted outward slightly such that they aredisplaced laterally (in this case to the north and south) from the center of thesources by a distance equal to the altitude above the source rather than 1/√2
times the altitude above the source (details are given in Appendix B). If weassume the dark lanes at Reiner Gamma (which are displaced ⇠ 5 km from thecenters of the bright lobes) are aligned with the cusps in the surface horizontal
magnetic field profile, then the depth of the source must be ⇠ 5 km. Because40we assume infinitely narrow source structures, this represents the maximum
possible source depth consistent with the hypothesis that the dark lanes alignwith cusps in the horizontal magnetic field profile. The LP MAG data are too
spatially coarse to constrain the distribution of magnetic moments within thesource structure so we assign equal strength to each of the 55 dipoles in thearray. As with the Airy case, we perform a complete search to find the solutionthat best fits, in a least squares sense, the LP MAG data. The resulting bestfit solution has the dipoles pointed with inclination +2° and declination -8°(i.e., pointed nearly horizontally and slightly west of north) and total magnetic
moment 1.0 ⇥ 1013 Am2 (Table 2.4). This source model produces a magnetic
field pattern at the surface that mimics the morphology of the albedo anomaly,
with the field being strongly horizontal in high albedo areas and vertical in
the intra-swirl dark lanes (Figure 2.16). Remarkably, this same source model
produces a field pattern at the spacecraft altitude that agrees with the Lunar
Prospector observations (Figure 2.17).

2.5

Discussion

2.5.1

Magnetization

The models described in section 2.4.1 suggest that the source material is concentrated under the central parts of the albedo anomalies. Based on the source
material distributions illustrated in Figures 2.8 and 2.10, we can compute theimplied magnetizations at Airy and Reiner Gamma for various assumed layer
thicknesses. Figure 2.18 illustrates that even when the magnetized layer isassumed to be 10 km thick, typical magnetizations at the Reiner Gamma

41

40

Altitude above source (km)

35
30
25
20
15
10
5
0
−30

−20

−10

0

10

Horizontal position (km)20

30

Altitude = 20 km

10Horizontal field strength (nT)

5
0
−30

−20

−10

0

10

20

30

Altitude = 5 km

600
400
200
0
−30

−20

−10

0

10

20

30

Altitude = 2 km

10000
5000
0
−30

−20

−10

0

10

Horizontal position (km)20

30

Figure 2.15: Upper panel: magnetic field lines due to a pair of dipoles separated
by 15 km and oriented horizontally in the plane of the page; Lower panel: profilesof the horizontal component of the magnetic field shown at the various altitudes
represented by dashed lines in the upper panel. Each of the dipoles has a magnetic
moment of 1012 Am2 .

42

Albedo Pattern-Based Model for Reiner Gamma
Latitude (°N)

Longitude (°E)

Latitude (°N)

Longitude (°E)

7.649

301.460

7.073

301.1437.673

301.4187.056

301.079

7.693

301.3717.043

301.0147.710

301.3207.035

300.9497.723

301.2647.031

300.884

7.731

301.205

7.031

300.8207.735

301.1447.035

300.7597.735

301.0807.043

300.700

7.731

301.0157.056

300.6447.723

300.9507.073

300.5937.710

300.8857.093

300.5467.693

300.8217.117

300.504

7.673

300.7597.542

301.7747.649

300.700

7.538

301.8197.622

300.6457.543

301.8787.592

300.5937.555

301.9437.559

300.546

7.573

302.001

7.525

300.505

7.593

302.0437.488

300.4697.612

302.0627.451

300.4397.628

302.0567.315

301.5257.636

302.0247.278

301.495

7.637

301.9737.241

301.4597.629

301.9107.207

301.4187.615

301.847

7.174

301.3717.596

301.794

7.144

301.3197.576

301.759

7.117

301.2647.557

301.7497.093

301.205

Table 2.4: Final dipole array model for Reiner Gamma as described in section 2.4.2.2(Figure 2.16). The model consists of 55 dipoles at the indicated latitudes andlongitudes, all buried 5 km below the surface. The total magnetic moment is
1.0 ⇥ 1013 Am2 . All dipoles in the array have inclination 2°, measured positive
downward from the horizontal and declination -8°, measured positive clockwise from
north.

43

9

9

0.14

7.5

750nm Reflectance

0.13

0.125
0.12

7

0.115
6.5

Latitude (degrees N)

900

0.135

8

300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

800
700

7.5

600
500

7

400
300
200

0.105
299.5

8

6.5

0.11

6

1000

8.5

Horizontal Field Strength (nT)

8.5

Latitude (degrees N)

1100

0.145

100

6

0.1

299.5

300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

Figure 2.16: Left panel: Reiner Gamma study area showing Clementine albedo map
with red circles indicating the locations of the source model’s 55 dipoles; Rightpanel: resulting horizontal magnetic field strength predicted at the surface. White
arrows indicate where field lines become vertical in one of the intra-swirl dark lanes.
Figure 2.17 compares the model field predicted at the spacecraft altitude with the
LP MAG observations.

anomaly are on the order of 1 A/m. If the magnetized layer is only 1 km
thick, the implied magnetization approaches 10 A/m at Reiner Gamma. For
comparison, Nicholas et al. (2007) predict a minimum magnetization of 1 A/m
for a layer 1 km thick at Reiner Gamma and Wieczorek et al. (2012) calculate⇠ 2 A/m for the same layer thickness assuming an anomaly that produces a
10 nT field at 30 km altitude. Terrestrial mid-ocean ridge basalts can havemagnetizations of ⇠ 1

6 A/m (Pariso and Johnson, 1991) and the magnetic

lineations in the Martian southern highlands require rock magnetization of⇠ 20 A/m (Connerney, 1999), suggesting that our lower bounds for magnetization at the Airy and Reiner Gamma anomalies are reasonable. However, if
we suppose the source structures are further horizontally concentrated, as thealbedo-pattern-constrained models of section 2.4.2 suggest (Figures 2.13 and2.16), the magnetizations would have to be even greater.

44

LP MAG

Model

0

7

−5

6.5

−10

6
302

Longitude (degrees E)

Latitude (degrees N)

9

−10

7.5
7

−20

6.5

−30

6
300

301

302

Longitude (degrees E)

303

Latitude (degrees N)

9
8

10

7.5

0

7

−10

6.5

−20
302

Longitude (degrees E)

303

9
8

30

7.5

20

7
6.5

10

6
300

301

302

Longitude (degrees E)

303

Total Field Strength (nT)

−5

6.5

−10

6
301

302

303
10

8.5

0

8

−10

7.5
7

−20

6.5

−30

6
301

302

303

9

−40

30

8.5

20

8

10

7.5

0

7

−10

6.5

−20
−30

6
300

301

302

Longitude (degrees E)

40

8.5

0

7

300

−30

6

5

7.5

Longitude (degrees E)

20

301

10

8

9

−40

30

8.5

300

15

Longitude (degrees E)

0

8

303

9

300
10

8.5

302

8.5

303

North Field Strength (nT)301

Latitude (degrees N)

5

East Field Strength (nT)10

Radial Field Strength (nT)

Latitude (degrees N)

15

7.5

301

Longitude (degrees E)

9300

Latitude (degrees N)

300

8.5
8

10

6

303

Latitude (degrees N)

302

Latitude (degrees N)

301

Longitude (degrees E)

303

9Latitude (degrees N)

300

20

7
6.5

East Field Strength (nT)

10

6

30

7.5

North Field Strength (nT)

20

7
6.5

8

Radial Field Strength (nT)

30

7.5

40

8.5

Horizontal Field Strength (nT)89Latitude (degrees N)

Total Field Strength (nT)

40

8.5

Horizontal Field Strength (nT)

Latitude (degrees N)

9

40

8.5
8

30

7.5

20

7
6.5

10

6
300

301

302

Longitude (degrees E)

303

Figure 2.17: Comparison of the observed field at ⇠ 18 km over Reiner Gamma (leftpanels; derived from LP MAG data collected on days 54, 81-82, and 109 of 1999) andthe model field (right panels; obtained as described in 2.4.2.2). Figure 2.16 shows
the horizontal component of the model field predicted at the surface.

45

3

10

Rein

2

Magnetization (A/m)

10

er G

Airy

amm

a

1

10

0

10

−1

10

−1

10

0

10

Layer Thickness (km)1

10Figure 2.18: Rock magnetization implied by source models illustrated in Figures 2.8
and 2.10 versus assumed layer thickness. Solid lines are based on the maximum
individual dipole moments found in those models (⇠ 2.1 ⇥ 1011 Am2 for Airy, ⇠6.5 ⇥ 1011 Am2 for Reiner Gamma) and the dashed lines are based on characteristic
values of 1.25 ⇥ 1011 Am2 for Airy and 3.4 ⇥ 1011 Am2 for Reiner Gamma.

2.5.2

Magnetizing Field

The strong magnetizations at the Airy and Reiner Gamma anomalies could
be the result of the source material having cooled in a long lasting global
magnetic field, perhaps generated by a core dynamo (Garrick-Bethell et al.,2009; Dwyer et al., 2011; Le Bars et al., 2011; Shea et al., 2012). As we haveshown here, however, the Airy anomaly (located at approximately 17°S, 3°E)
is magnetized with a steep downward inclination while the Reiner Gamma
anomaly (located at approximately 7°N, 59°W) is magnetized with almost
zero inclination and points approximately toward the north. If these two
anomalies acquired their magnetizations by cooling in a dipolar dynamo field,
they could not have formed contemporaneously. Instead, the two anomalies
may have formed during diﬀerent global field orientation epochs. Alternatively,
the Moon’s dynamo field may have had substantial higher order components
46

(i.e., beyond dipolar).
Another possibility that avoids the diﬃculties associated with inconsistentmagnetization directions is that, rather than being the thermal remanent magnetization signatures of an extinct core dynamo, the Moon’s crustal magnetic
anomalies may instead be the result of shock remanent magnetization that occurs following basin-forming impact events when magnetohydrodynamic shockwaves converge near the basin antipode producing strong transient fields andpotentially magnetizing iron-rich ejecta materials (Hood et al., 2001; Halekas
et al., 2001; Richmond et al., 2005; Hood and Artemieva, 2008).

2.6

Conclusions

Our examination of swirls at Airy and Reiner Gamma, two magnetic anomalies with dissimilar orientation, suggests that magnetic field direction and swirl
morphology are related in the way we predict based on the solar wind deflection
hypothesis: the Reiner Gamma case delivers evidence that swirls are brightest where magnetic field lines are dominantly horizontal and the Airy casedemonstrates a connection between dark lanes and vertically oriented fieldlines. These findings support the solar wind deflection model for swirl formation, implying that diﬀerential solar wind darkening is largely responsible forcreating the albedo anomalies. Although our source models do not represent
unique solutions, they agree with observational constraints while plausibly accounting for the alternating bright and dark bands at both Airy and Reiner
Gamma. Our model results suggest that swirl morphology could potentiallybe used to infer small-scale structure in the near-surface magnetic field as wellas the layout and burial depth of the magnetic source material. For Airy and47

Reiner Gamma, we infer maximum source burial depths of ⇠ 8 km and ⇠ 5 km,
respectively, and horizontally concentrated sources with strong magnetizations(⇠ 10 A/m or greater for a layer 1 km thick). Examination of additional swirls
may help to further our understanding of how magnetic field direction relates
to swirl morphology, but ultimately, very near-surface magnetic field and solar
wind flux measurements (i.e., from altitudes of hundreds of meters or less) will
be required to confirm our predictions.

Acknowledgements
This research was performed as part of the World Class University Projectat Kyung Hee University and sponsored by the Korean Ministry of Education, Science and Technology. Additional support was provided by the NASA
Ames Research Center/University of California Santa Cruz University Aﬃliated Research Center. We would like to thank Lon Hood, Robert Lin, Jasper
Halekas, Greg Delory, Robert Lillis, and Michael Purucker for their comments
and suggestions and Joseph Nicholas for the data used to generate Figure 2.12.

48

Chapter 3
Latitudinal Variation in Spectral
Properties of the Lunar Maria and
Implications for Space Weathering
This chapter is a modified version of Hemingway, D. J., I. Garrick-Bethell, and
M. A. Kreslavsky (2015), Latitudinal Variation in Spectral Properties of the
Lunar Maria and Implications for Space Weathering, Icarus (in press).

Abstract
Space weathering alters the optical properties of exposed surfaces over time,complicating the interpretation of spectroscopic observations of airless bodies
like asteroids, Mercury, and the Moon. Solar wind and micrometeoroids arelikely the dominant agents of space weathering, but their relative contributions
are not yet well understood. Based primarily on Clementine mosaics, we reporta previously unrecognized systematic latitudinal variation in the near-infrared49

spectral properties of the lunar maria and show that the characteristics of thislatitudinal trend match those observed at ’lunar swirls’, where magnetic fields
alter local solar wind flux without aﬀecting the flux of micrometeoroids. We
show that the observed latitudinal color variations are not artifacts of phaseangle eﬀects and cannot be accounted for by compositional variation alone.
We propose that reduced solar wind flux, which should occur both at swirls
and toward higher latitudes, is the common mechanism behind these color
variations. This model helps us quantify the distinct eﬀects of solar wind andmicrometeoroid weathering and could aid in interpreting the spectra of airlessbodies throughout the solar system.

3.1

Introduction’Space weathering’ refers to the processes by which the optical properties ofairless bodies change due to exposure to solar wind and micrometeoroid impacts. However, the diﬃculties of reproducing space-weathering conditions in
the laboratory, or returning weathered samples to Earth, make it challenging to determine precisely how space weathering operates (Pieters et al., 2000;
Hapke, 2001; Vernazza et al., 2009; Pieters et al., 2012; Domingue et al., 2014).
Remote sensing measurements, studies of lunar samples, and laboratory experiments have established that solar wind ion and micrometeoroid bombardment
weaken spectral absorption features and cause the lunar surface to darken andredden (increase in spectral continuum slope in the visible and near-infrared)
with time. These changes appear to be due to some combination of the formation of impact glasses and agglutinates (Adams and McCord , 1971), theregolith’s disintegration into increasingly finer soils (Pieters et al., 1993), and50

the accumulation of nanophase iron (Hapke, 2001; Sasaki et al., 2001; Nobleet al., 2007). Larger impacts also expose fresh material, which then graduallymatures until the reflectance spectrum reaches a steady state, which we call’equilibrium color’ for simplicity.
The equilibrium color varies considerably across the lunar surface, due
primarily to diﬀerences in mineralogy. This is most obvious in the dichotomy
between the bright, anorthositic highlands and the darker basaltic maria. However, as we will argue, the presence of ’lunar swirls’ suggests that equilibriumcolor may also be influenced by the flux of weathering agents, rather than justtheir total accumulation (see sections 3.3.1 and 3.4). If this is the case, thenequilibrium color may also vary with latitude. Both solar wind and micrometeoroids originate primarily from within the ecliptic plane, which is inclined
from the Moon’s equator by just 1.5°. Hence maximum flux of these weathering agents occurs near the equator, with flux decreasing as incidence angle
increases towards the poles.
This paper’s central observation is that, when we examine imagery fromacross the lunar surface, we find that the equilibrium color does vary systematically with latitude. In section 3.3.2, we show that this latitudinal color trend
persists across a range of distinct compositions and that it is not an artifactof phase angle biases in the Clementine mosaics. Interestingly, the spectral
properties of the latitudinal color trend match the characteristic color variation found at lunar swirls. In section 3.3.1, we quantify the characteristics ofthe swirl-related color variation and, in section 3.3.2, we show that it is statistically equivalent to the observed latitudinal color trends, with a transition
toward higher latitudes being attended by the same color change that occurs
towards brighter parts of swirls. Finally, in section 3.4, we argue that the best51

candidate for a common mechanism behind these color variations is alteredsolar wind flux. We present a qualitative empirical model illustrating how thishypothesis comports with the observations and we discuss the possible implications with respect to the interpretation of spectral data, particularly at highlatitudes.

3.2

Data Sources

In this study, we primarily use mosaics based on imagery from the 750 nm
and 950 nm channels of the Clementine UVVIS (ultraviolet-visible) experiment (Nozette et al., 1994; Eliason et al., 1999), available from the USGS
(www.mapaplanet.org). As a point of comparison, we also examine 1064 nm
reflectance from the Lunar Orbiter Laser Altimeter (LOLA) experiment onboard the Lunar Reconnaissance Orbiter (LRO) (Lucey et al., 2014). In discussing the observed trends in the Moon’s spectral properties, we may use theword ’color’ in a general sense to refer to combinations of albedo and the ratio
between 950 nm and 750 nm reflectance (e.g., as a proxy for continuum slope).
Parts of our analysis require isolating portions of the lunar surface according to composition and/or topographic roughness. For composition, we
use results from the Lunar Prospector Gamma Ray Spectrometer (Lawrence
et al., 2002; Prettyman et al., 2006), specifically in order to identify FeO andTiO2 content in the regolith. The topographic roughness metric we use is theinterquartile range of the along-profile second derivative of elevation, at 1.8-km
baseline (Kreslavsky et al., 2013), derived from Lunar Orbiter Laser Altimeter(LOLA) data. The latter is used to distinguish between the smooth maria andthe rougher highlands.
52

3.3

Analysis

Before discussing the observed latitudinal color variation, we revisit the characteristic color signature observed at swirls, developing a new parameterizationthat will allow for a quantitative comparison between swirls and the newly
observed latitudinal trends.

3.3.1

Color Variation at Lunar Swirls

Lunar swirls are enigmatic collections of sinuous bright markings, often interposed with narrow dark lanes, that are co-located with many of the Moon’s
crustal magnetic anomalies (Figure 3.1A). The bright parts of swirls superficially resemble optically immature surfaces such as fresh impact craters (Luceyet al., 2000a; Wilcox et al., 2005; Blewett et al., 2011). However, it has beenshown (Garrick-Bethell et al., 2011) that swirls exhibit spectral trends that are
distinct from those associated with impact-related brightening (Lucey et al.,2000a). The two trends can be distinguished from one another, using Clementine UVVIS (ultraviolet-visible) mosaics (Nozette et al., 1994; Eliason et al.,1999), by plotting 750 nm reflectance against the 950 nm to 750 nm reflectanceratio (Garrick-Bethell et al., 2011; Blewett et al., 2011), the former representing
albedo and the latter serving as a proxy for both the near-infrared continuum
slope and the 1 mm absorption feature found in iron-bearing silicate minerals.
Both the swirl- and impact-related color variations involve changes in bothalbedo and the 950 nm/750 nm band ratio, but the impact-related variation isaccompanied by a proportionally greater change in the 950 nm/750 nm band
ratio (Figure 3.1B), as originally reported by Garrick-Bethell et al. (2011).
In order to establish a quantitative basis for comparison with the latitu53

0.15

A

0.12

7

0.11

301

301.5

302

Longitude (degrees E)

302.5

303

9

1

0.85
0.8
0.75

0.09

0.7
0.09

0.76

9

2

300.5

pe m

em

6
300

0.1

10km

0.1

0.11

0.12 0.13 0.14 0.15
750nm Reflectance

0.16

8.5

8

0.74

8

7.5

0.73

7

0.72

6.5

0.71

0.18

0.295
0.29
0.285
0.28

Beta

Alpha

0.75

Latitude (degrees N)

8.5

0.17
0.3

D

C

Latitude (degrees N)

swirl
equilibrium
color

lop

6.5

0.9

s
d,

7.5

`

rl tr
end
, slo

n
tre

0.13

ing

swi

ct

8

0.95

B

pa
im

Clementine V1 750nm Reflectance

0.14

g_
sin
rea
inc

Latitude (degrees N)

8.5

backgroundequilbirium color
incr
eas

1

950nm−to−750nm Reflectance Ratio

9

7.5

0.275
0.27

7

0.265
0.26

6.5

0.255
6
300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

0.7

6
300

300.5

301

301.5

302

Longitude (degrees E)

302.5

303

0.25

Figure 3.1: (A) 750 nm reflectance (albedo) at the Reiner Gamma swirl. (B) Spectral
characteristics (albedo versus 950 nm/750 nm band ratio) of pixels sampled fromthe five rectangular boxes, with corresponding colors, in (A). The red point cloudcorresponds to a representative background region whereas the dark blue point cloudcorresponds to the brightest part of the swirl. Each of the point clouds displays asteep trend that follows the progression from fresh impact craters (lower right part ofeach point cloud) to darker background soils (upper left part of each point cloud); thedashed black line is a best fit through the red point cloud and, averaging over three
separate mare swirl areas, has a slope of m2 = 5.7±0.5. The black circles represent
the steady state (equilibrium) colors reached with maximum optical maturity in eachrectangle (since the majority of pixels are likely near saturation maturity, we takethe median point in each cluster to represent the equilibrium color). The solid blackline is a best fit through the centroids (black circles) of each cluster and has a slopeof m1 = 1.6±0.2, averaged over the three separate mare swirl areas. (C) and (D)are maps of the ↵ and parameters, defined according to equations (3.1) and (3.2),
respectively.

54

−32.5

0.2

A

0.17

−34.5

0.16

−35

0.15

−35.5

0.14

−36

0.13

−36.5

0.12

−37.5
162

Clementine 750nm Reflectance

0.18

−34

15km
163

164

165

166

167

−32.5

C

0.95
0.9
0.85
0.8

0.1

0.7

0.84

−32.5

0.1

0.12

0.14
0.16
750nm Reflectance

0.18

0.2
0.35

D

0.83

−33

−33.5

0.82

−33.5

0.33

0.81
0.8

Alpha

−34
−34.5
−35

0.79

−35.5

0.78

−36

0.77

−36.5
−37
−37.5
162

163

164

165

Longitude (degrees E)

166

167

0.34

−34

0.32

−34.5

0.31

Beta

−33

1

0.75

0.11

Longitude (degrees E)

B

1.05
950nm−to−750nm Reflectance Ratio

−33.5

−37

Latitude (degrees N)

1.1

0.19

Latitude (degrees N)

Latitude (degrees N)

−33

−35

0.3

−35.5

0.29

−36

0.28

0.76

−36.5

0.27

0.75

−37

0.26

0.74

−37.5
162

163

164

165

Longitude (degrees E)

166

167

0.25

Figure 3.2: (A) 750 nm reflectance at the Ingenii swirl. (B) Spectral characteristics(albedo versus 950 nm/750 nm band ratio) of pixels sampled from the five rectangularboxes, with corresponding colors, in (A). (C) and (D) are maps of the ↵ andparameters, defined according to equations (3.1) and (3.2), respectively.

55

16

0.17

1.1

A
0.16

Clementine 750nm Reflectance

15

0.14
0.13

14.5

0.12
14

0.11

10km
82.5

83

83.5

Longitude (degrees E)

84

84.5

16

C

0.95
0.9
0.85
0.8
0.75

0.1

0.7

0.82

16

0.1

0.12

0.14
0.16
750nm Reflectance

0.18

0.325

15.5

0.8
0.79

Alpha

15

0.78
14.5

0.77
0.76

14

0.32
0.315

15

0.31
0.305

14.5

0.3
0.295

14

0.29

0.75
13.5
82

82.5

83

83.5

Longitude (degrees E)

84

84.5

0.2
0.33

D

0.81

15.5

Latitude (degrees N)

1

Beta

13.5
82

950nm−to−750nm Reflectance Ratio

0.15

Latitude (degrees N)

Latitude (degrees N)

15.5

B

1.05

0.285
13.5
82

0.74

82.5

83

83.5

Longitude (degrees E)

84

84.5

Figure 3.3: (A) 750 nm reflectance at the Marginis swirl. (B) Spectral characteristics(albedo versus 950 nm/750 nm band ratio) of pixels sampled from the five rectangularboxes, with corresponding colors, in (A). (C) and (D) are maps of the ↵ andparameters, defined according to equations (3.1) and (3.2), respectively.

56

dinal trends we discuss in section 3.3.2, we parameterize the color variations
that are characteristic of swirls, averaging over three diﬀerent mare swirl areas: Reiner Gamma in western Oceanus Procellarum, Mare Ingenii on thefarside, and Mare Marginis on the eastern limb. In each case, in the albedo
versus band ratio diagrams, we found similar steep trends associated withthe transition between impact craters and background soils, and shallower
trends associated with the transition between dark and bright parts of swirls(Figure 3.1B, Figure 3.2B, Figure 3.3B), in accord with Garrick-Bethell et al.
(2011). The distinct color variations associated with impacts and swirls allow
us to define parameters that clearly separate the two trends (Figure 3.1C and
D). The impact-related progression from bright craters to the more maturebackground soils can be characterized by an impact maturity parameter↵ = R750

✓

R950
R750

◆

/m1

(3.1)

where R750 and R950 are the Clementine 750 nm and 950 nm reflectances,respectively, and wherem1 =1.6±0.2is the slope of the swirl-related trends (±1 ), averaged from the three separate
mare swirl areas. The swirl-related trend slope is used in equation (3.1) so that
impact maturity (↵) is not aﬀected by swirl-related color variations. Equation (3.1) resembles previously developed optical maturity parameters (Luceyet al., 2000a; Wilcox et al., 2005) except that here, the goal is explicitly to isolate the impact-related color variation from that associated with swirls, and so
the constants are diﬀerent. Similarly, we can represent the swirl-related color

57variation, which we regard as distinct from optical maturity, as= R750

✓

R950
R750

◆

/m2

(3.2)

where
m2 =

5.7±0.5is the typical slope of the impact-related trends (±1 ). The impact-related
trend slope is used in equation (3.2) so thatis not aﬀected by impact-relatedcolor variations.
Although the ↵ andvalues vary according to local composition, the slopesof the impact- and swirl-related trends do not vary significantly across diﬀerent mare regions. The values given here for m1 and m2 are therefore largely
composition independent, at least within the maria. The ↵ parameter is designed to have constant values along the swirl-related trends such that swirl
features do not influence the value of ↵ and so maps generated for the ↵ parameter show impact features but not swirl features (Figure 3.1C). Conversely,
theparameter is designed to have constant values along the impact-related

trends such that maps of theparameter highlight swirl features while muting

impact features (Figure 3.1D).

3.3.2

Latitudinal Color Variation

When we examine imagery from across the lunar surface, we find that thereflectance spectra vary systematically with latitude. The eﬀect is not obvious
when we examine the Moon as a whole, likely because the spectra are sostrongly aﬀected by composition, which varies considerably across the surface.

58However, when we account for variations in composition, the latitudinal trends
emerge. As we will show, the latitudinal trends are especially pronounced
within the maria, and may help to account for the unusually high albedo ofMare Frigoris—the highest latitude mare region.
Color anomalies have been identified previously in polar regions (Yokota
et al., 2011; Zuber et al., 2012) and latitudinal variation in space weathering has
been considered previously (Yokota et al., 2011; Hendrix et al., 2012), however,this is the first observation of a broad systematic latitudinal color trend thatis visible across the lunar maria. Moreover, the spectral characteristics ofthis trend match those found at lunar swirls (section 3.3.1), suggesting a link
between the two phenomena that may be helpful for understanding space
weathering (see section 3.4).
In addition to controlling for composition, below we examine the possible
contributing eﬀects of highland contamination in the maria and phase angle
biases in the Clementine mosaics; we show that the observed latitudinal trends
cannot be artifacts of such eﬀects.
3.3.2.1

Eﬀect of Composition

Because surface color is strongly aﬀected by composition, we separate thedata into bins according to TiO2 and FeO content (as measured by the Lunar Prospector Gamma Ray Spectrometer (Lawrence et al., 2002; Prettyman
et al., 2006)), and examine the spectral characteristics separately for eachcompositional bin. Figure 3.4A shows, in the same parameter space we used
to characterize swirls, that pixels sampled from higher latitudes tend to havehigher albedo (750 nm reflectance) and lower 950 nm to 750 nm band ratios,
than those sampled from lower latitudes.
59

B

1.15

0.95
0.9
0.85

1

0.8

0.95

0.75

Cosine Latitude

1.1
1.05

0.9

0.7

0.85

0.65

0.8

0.6

0.75
0.7

15−17 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

0.12
0.14
0.16
750nm Reflectance1.1

950nm−to−750nm Reflectance Ratio

950nm−to−750nm Reflectance Ratio

1.15

1

A

0.2

0.22

1

0.95

0.9

0.85

0.8

0.55

0.18

1.05

0.5

0.75
0.08

0.1

0.12

0.14

0.16

0.18

0.2

750nm Reflectance

0.22

0.24

0.26

0.28

Figure 3.4: Color variation as a function of latitude. (A) Spectral characteristics ofpixels sampled from regions with between 15 and 17 wt% FeO and between 1 and 2
wt% TiO2 . Points are color-coded according to cosine latitude, illustrating that the
optical characteristics of the surface vary systematically with latitude. Equatorial
regions (red points) tend to have lower albedo and higher 950 nm/750 nm reflectanceratios than regions farther from the equator (blue points). The black trend line
represents a least squares best fit through the data points. The solid black circle
corresponds to the color predicted by the best-fit line at the equator and the whitecircle represents the color predicted at latitudes of ±70°. (B) Compilation of the 59
latitudinal trends we examined (grey lines), one for each compositional bin, comparedwith the trends discussed in section 3.3.1: the solid red, green and blue lines represent
the swirl-related color trends observed at Reiner Gamma, Ingenii, and Marginis,
respectively; the dashed lines represent the impact-related color trends observed ateach of those swirls.

60

While Figure 3.4A shows just one illustrative example, we examined all
combinations of FeO and TiO2 content using bins of 2 wt% width in FeO
and 1 wt% width in TiO2 , set at half bin width increments (see Figure 3.5
for additional examples). Because the latitudinal trend can only emerge if
pixels are sampled from a wide range of latitudes, we discard compositional
bins whose pixels do not reach latitudes of at least ±50° or do not span atleast 50 total degrees of latitude. To ensure that slopes are estimated fromrobust samples only, we also discard bins whose pixels cover less than 0.25%
of the lunar surface. Finally, to avoid artifacts in the Clementine mosaics
that are related to severe illumination conditions at high latitudes, we also
discard data within 20° of the poles. After applying these selection criteria, 59
compositional bins were retained for analysis. Figure 3.4B shows a summaryof all 59 latitudinal trend lines (grey lines) along with the swirl- and impactrelated trends discussed in section 3.3.1.
The eﬀect of iron content is apparent when the trend lines are color-codedby weight % FeO, as determined by the Lunar Prospector Gamma Ray Spectrometer (Lawrence et al., 2002) (Figure 3.6). Lower iron regions (pale orange
lines) tend to plot farther up and to the right in Figure 3.6, in accord with
Wilcox et al. (2005), whereas high iron regions (dark lines) plot farther down
and to the left. This means that both ↵ andtend to be larger for low iron

regions (in section 3.3.2.5, we discuss the implications of this observation withrespect to highland contamination).
Although the latitudinal trend slopes are largely similar across diﬀerent
compositional bins, they become less consistent when iron content is very
low, such as in the highlands. This is not surprising given that several studies
suggest that the production of nanophase iron plays an important role in space
61

0.8
0.75

0.12

0.14
0.16
750nm Reflectance

0.18

0.2

0.7

0.22

1.1

950nm−to−750nm Reflectance Ratio

1.15

1.1
1.05
1
0.95
0.9
0.85
0.8

0.7

950nm−to−750nm Reflectance Ratio

0.1

0.8

1.15

0.75

14−16 wt% FeO
1.5−2.5 wt% TiO2
0.08

0.1

0.12

0.18

0.2

0.8

1.1

1

0.9
0.85
0.8

0.7

13−15 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

0.12

0.18

0.2

1.15
1.1

1
0.95
0.9
0.85
0.8
0.75
0.7

13−15 wt% FeO
0.5−1.5 wt% TiO2
0.08

0.1

0.12

0.18

0.5
±60

0.2

0.22

0.08

0.1

0.12

0.14
0.16
750nm Reflectance

0.18

0.2

0.08

0.1

0.12

0.14
0.16
750nm Reflectance

0.18

0.2

0.22

0.14
0.16
750nm Reflectance

0.18

0.2

0.22

0.14
0.16
750nm Reflectance

0.18

0.2

0.22

0.18

0.2

0.22

1.1

1
0.95
0.9
0.85
0.8

18−20 wt% FeO
1.5−2.5 wt% TiO2
0.08

0.1

0.12

1.15

15−17 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

1.1
1.05
1
0.95
0.9
0.85
0.8
0.75

0.12

0.14
0.16
750nm Reflectance

0.18

0.2

0.7

0.22

17−19 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

0.12

1.15

1

0.9
0.85

0.7

18−20 wt% FeO
2.0−3.0 wt% TiO2

1.05

0.7

0.22

0.95

0.8

0.7

0.75

1.05

0.75
0.14
0.16
750nm Reflectance

16−18 wt% FeO
1.5−2.5 wt% TiO2

0.9

1.1

0.8

1.15

0.85
0.8

0.9
0.85

0.22

1

1.15

1.05

0.2

0.95

0.7

0.22

0.18

1.05

0.75
0.14
0.16
750nm Reflectance

0.14
0.16
750nm Reflectance

0.9

1.15

0.95

0.12

0.85

1.1
1.05

0.1

1

0.7

0.22

0.08

1
0.95

0.75

0.95

0.75
0.14
0.16
750nm Reflectance

16−18 wt% FeO
2.0−3.0 wt% TiO2

1.05

1.15

0.75

950nm−to−750nm Reflectance Ratio

0.08

0.9
0.85

0.75

950nm−to−750nm Reflectance Ratio

950nm−to−750nm Reflectance Ratio

0.7

14−16 wt% FeO
2.0−3.0 wt% TiO2

1.1
1.05

950nm−to−750nm Reflectance Ratio

0.9
0.85

1
0.95

950nm−to−750nm Reflectance Ratio

1
0.95

1.05

950nm−to−750nm Reflectance Ratio

1.05

1.15
950nm−to−750nm Reflectance Ratio

1.1

950nm−to−750nm Reflectance Ratio

1.15

1.1

950nm−to−750nm Reflectance Ratio

950nm−to−750nm Reflectance Ratio

1.15

15−17 wt% FeO
0.5−1.5 wt% TiO2
0.08

0.6

0.1

1
0.95
0.9
0.85
0.8
0.75

0.12

0.7
±50

1.1
1.05

0.14
0.16
750nm Reflectance

Cosine Latitude±40

Latitude (°)

0.18

0.2

0.22

0.8

0.9
±30

0.7

17−19 wt% FeO
0.5−1.5 wt% TiO2
0.08

0.1

0.12

0.14
0.16
750nm Reflectance1
0

Figure 3.5: Latitudinal color variation across 12 distinct compositional bins. Equatorial regions (red points) tend to have lower albedo and higher 950 nm/750 nm
reflectance ratios than regions farther from the equator (blue points). Each of theblack trend lines represents the best-fit color variation between the equator (solid
black circle) and latitudes of ±70° (white circle). Although only 12 are shown here,
we examined 59 compositional bins in total.

62

25

1.15

20

1.05

1

15

FeO wt%950nm−to−750nm Reflectance Ratio

1.1

0.95

10

0.9

0.85

5

0.8

0.75
0.08

0.1

0.12

0.14

0.16

0.18

0.2

750nm Reflectance

0.22

0.24

0.26

0.28

0

Figure 3.6: Compilation of latitudinal trend lines for all 59 compositional bins, colorcoded by iron content (Lawrence et al., 2002). Dark trend lines correspond to highiron bins while pale trend lines correspond to low iron bins. For comparison, theswirl- and impact-related trends are also shown: the solid red, green and blue lines
represent the swirl-related color trends observed at Reiner Gamma, Ingenii, andMarginis, respectively; the dashed lines represent the impact-related color trends
observed at each of those swirls.

weathering (Hapke, 2001; Sasaki et al., 2001; Noble et al., 2007). A paucity ofiron in the highlands may prevent a clear appearance of the type of latitudinal
trends we observe in the maria. In order to obtain as clear a signal as possible,
we therefore focus our subsequent analysis on the maria.
3.3.2.2

Isolation of the Maria

Because the latitudinal trends are most consistent within the maria, and because we are interested in comparing these trends with the spectral characteristics we analyzed for three mare-based swirls (section 3.3.1), we restrict theremainder of our analysis to the lunar maria.
As previously discussed, we use Lunar Prospector Gamma Ray Spectrometer (GRS) data to determine both iron and titanium content in the soils we
examine. Unfortunately, in comparison to the Clementine data, the relatively

63

coarse resolution of the GRS data (the eﬀective footprint is roughly 45 km wide
(Lawrence et al., 2002)) means that we cannot use GRS-based iron estimates
to isolate the maria—choosing a conservatively high threshold for iron content
would result in the exclusion of significant portions of the maria, while choosing a low threshold would result in the inclusion of highland pixels. Inclusionof low-iron highland pixels is undesired as it could artificially skew the spectral
characteristics towards higher 750 nm reflectance and higher 950 nm/750 nm
reflectance ratios (see Figure 3.6 and section 3.3.2.5).
We cannot isolate the maria using iron estimates obtained from Clementinebased spectral techniques (e.g., (Lucey et al., 2000b; Wilcox et al., 2005)) because these techniques rely on the same parameter space we use in our own
spectral analysis, namely 750 nm reflectance versus the 950 nm/750 nm reflectance ratio. Attempting to isolate the maria according to such spectralbased iron estimates would necessarily bias our results, artificially altering theslope of the best-fit latitudinal trends shown in Figure 3.4.
To distinguish mare from highland terrain, we instead adopt an independent metric that characterizes topographic roughness at 1.8-km baseline(Kreslavsky et al., 2013). A histogram of the topographic roughness (Figure 3.7) shows a bimodal distribution, reflecting the dichotomy of kilometerscale roughness between smoother maria and rougher highlands, with mariapeaking at a roughness value of 2.3 m/km2 and highland terrain peaking at aroughness value of approximately 19 m/km2 . We identify mare terrain as areas exhibiting roughness values below 5 m/km2 . After applying a de-speckling
filter, designed to remove isolated features less than approximately one degree
in diameter, we obtain a mask that we use to exclude non-mare pixels from oursubsequent analysis. Figure 3.8 illustrates the boundaries of our mare mask.
64

8

4

x 10

7

Number of Pixels

6
5
4
3
2
1
0
0

10

20
30
LOLA 1.8−km Roughness (m/km2)

40

50Figure 3.7: Histogram of topographic roughness with red line indicating the threshold below which we consider pixels to belong to the lunar maria; pixels above thethreshold are discarded in our subsequent analysis.

0.35

60

0.3

30

0.25

0

0.2

−30

0.15

−60
−90
−180

Clementine 750nm Reflectance

Latitude (degrees N)

90

0.1
−135

−90

−45

0

45

Longitude (degrees E)

90

135

180

Figure 3.8: Clementine 750 nm reflectance with red lines indicting the boundaries ofthe mare regions as identified by our topographic roughness criterion (regions withless than 5 m/km2 roughness at 1.8-km baseline).

65

90

A

60

60

30

30

0
−30

Latitude (degrees N)

Latitude (degrees N)

90

B

0
−30

−60

−60

−90

−90
0.1
0.15
0.2
Clementine 750nm Reflectance

0.15
0.2
0.25
LOLA 1064nm Reflectance

Figure 3.9: (A) Clementine 750 nm reflectance profile across the lunar maria (asidentified by surface roughness, see Figure 3.8). (B) Profile of LOLA 1064 nm
reflectance for the same region. In both profiles, reflectance tends to be lowest near
the equator and higher towards higher northern and southern latitudes. The blacklines and grey bands represent the mean and standard deviation at each latitude.

After excluding data from the lunar highlands, a clear latitudinal trend
emerges in the Clementine 750 nm reflectance profile (Figure 3.9A): maresurfaces are darkest near the equator and become increasingly bright towards
higher latitudes, with the profile being approximately symmetric about theequator. This is a straightforward visualization of the latitudinal trend andallows for simple comparison with diﬀerent datasets (e.g., reflectance from theLunar Orbiter Laser Altimeter; see section 3.3.2.6 and Figure 3.9B), but inorder to facilitate comparison with the swirl-related spectral trends, we now
return to plots of 750 nm reflectance versus the 950 nm/750 nm reflectanceratio.
3.3.2.3

Comparison with Swirl Trends

After excluding the highlands, 28 of the original 59 compositional bins survive
the selection criteria we described in section 3.3.2.1. Figure 3.10 illustrates not
only that the latitudinal color trends persist across a range of diﬀerent compo66

1.15

increasing `

latitudin

1.05

al trends, slo

1

pe m

3

0.95

inc

0.9

swir

as
re

nds

ren
tt
ds

0.8

0.75
0.08

l tre

c
pa

_

0.85

im

ing

950nm−to−750nm Reflectance Ratio

1.1

0.09

0.1

0.11

0.12

0.13

0.14

750nm Reflectance

0.15

0.16

0.17

0.18

Figure 3.10: Summary of latitudinal trend lines (grey lines) for the 28 compositional
bins that correspond to the maria. The solid red, green and blue lines represent
the swirl-related color trends observed at Reiner Gamma, Ingenii, and Marginis,
respectively. The dashed lines represent the impact-related color trends observed ateach of those swirls. The latitudinal trends have slopes of 1.5 ± 0.3, whereas the
swirl-related trends have slopes of 1.6 ± 0.2.

sitions, but also that these trends are remarkably similar to those associatedwith lunar swirls. That is, low latitude regions tend to have lower albedo andhigher 950 nm/750 nm band ratios when compared with high latitude regions,and these changes occur in the ratio m3 = –1.5 ± 0.3, a slope that is indistinguishable from the swirl-related slope, m1 = –1.6 ± 0.2, at the one-sigma level(Figure 3.11). Figure 3.12 further demonstrates this result by showing thata transition to higher latitudes is attended by the same type of color changethat occurs towards brighter parts of swirls (Figure 3.1D)—that is, increasing.
3.3.2.4

Mare Frigoris

Mare Frigoris, spanning much of the near side at roughly 60°N, is the highest
latitude and visibly brightest mare region, meaning that it makes a significant
contribution to the latitudinal trend we observe (Figure 3.12). The high
67

2

2

A

0

0

−1

−1

Spectral Trend Slope

Spectral Trend Slope

B

1

1

−2
−3
−4
−5

−2
−3
−4
−5

−6

−6

−7

−7

−8
11

12

13

14

15

16

FeO wt%

17

18

19

−8

20

0

0.5

1

1.5

2

TiO2 wt%

2.5

3

Figure 3.11: Latitudinal trend slopes (crosshairs) compared with swirl trend (pale redband) and impact trend (pale green band) slopes (with ±1 variability illustrated
by the width of each band). All 28 mare-based slope estimates (the same trends
illustrated in Figure 3.10) are shown versus the FeO bins (A) and the TiO2 bins (B).
In each case, the crosshairs show the width of the FeO or TiO2 bin and the ±1
uncertainty on the estimated slope.

90

Mare Frigoris

0.31

Mare
Crisium30

0.29

Oceanus
Procellarum

Beta

0

0.3

0.28

−30

0.27

−60

−90
−180

−135

−90

−45

0

45

Longitude (degrees E)

90

135

180

B

60

Latitude (degrees N)

60

Latitude (degrees N)

90

0.32

A

30
0
−30

0.26

−60

0.25

−90

0.25

0.3
Beta

0.35

Figure 3.12: (A) Map of , computed according to equation (3.2), over the lunar
maria (as identified by surface roughness, see section 3.3.2.2). (B) Latitudinal profileof across the lunar maria. The black line and grey band represent the mean andstandard deviation at each latitude, averaged over all longitudes.

68

(high albedo and low 950 nm/750 nm band ratio) of Mare Frigoris could be
attributed in part to its low FeO and TiO2 content (Lucchitta, 1978; Pieters,
1978; Staid and Pieters, 2000). However, the FeO and TiO2 abundances typical of Mare Frigoris (14±2 wt% FeO; 0.9±0.4 wt% TiO2 , see crosshairs in
Figure 3.13A) are also found in lower latitude mare regions, which, by comparison, appear darker and have higher 950 nm/750 nm band ratios (andtherefore lower ). The albedo and band ratio of Mare Frigoris diﬀer significantly from those of low-latitude maria of equal composition (Figure 3.13B,
see also lower left panel in Figure 3.5).
Given that Mare Frigoris is apparently compositionally unexceptional, its
spectral characteristics cannot be attributed merely to low FeO and TiO2 content. While it remains possible that they are the result of anomalous surface
contamination (see section 3.3.2.5), our results suggest that the high albedo
and low 950 nm/750 nm band ratio of Mare Frigoris are better understood aspart of the broad latitudinal trend we observe across the lunar maria. Although
Mare Frigoris accounts for an important component of the observed latitudinal
trend, the analysis is not significantly aﬀected by its removal (Figure 3.14).
After manually removing those mare pixels belonging to Mare Frigoris, fiveof the original 28 compositional bins must be excluded from analysis due totheir reduced area and latitudinal span, leaving a total of 23 bins. After Mare
Frigoris is removed, fewer high latitude pixels are included. Nevertheless, thelatitudinal trends remain, and the estimated latitudinal trend slope is unaffected (m3 = –1.5 ± 0.3 with or without Mare Frigoris). Hence, anomalous
spectral characteristics of Mare Frigoris alone cannot account for the latitudinal color trend, but the latitudinal color trend can account for the spectral
characteristics observed at Mare Frigoris.
69

6

Clementine 950nm−to−750nm Reflectance Ratio

A

5

TiO2 wt%

4

3

2

1

0
0

5

10

15
FeO wt%

20

25

1.15

B

1.1
1.05
1
0.95
0.9
0.85
0.8
0.75
0.7

30

0.08

0.1

0.12
0.14
0.16
0.18
Clementine 750nm Reflectance

0.2

0.22

Figure 3.13: (A) FeO and TiO2 content (based on Lunar Prospector Gamma Ray
Spectrometer data (Lawrence et al., 2002; Prettyman et al., 2006)) of Mare Frigoris(blue points) compared to other mare regions (red points). The black crosshairs
illustrate the mean and standard deviation FeO and TiO2 abundances for Mare
Frigoris (14±2 wt% FeO; 0.9±0.4 wt% TiO2 ). (B) Pixels from Mare Frigoris (bluepoints) have higher albedo and lower band ratio than low latitude mare regions (redpoints, restricted to within ±15° of the equator) of equal composition. The meanand standard deviation of each cluster is marked with crosshairs.

1.15

1.15

A

1.05

1

0.95

0.9

0.85

0.8

0.75
0.08

B

1.1

950nm−to−750nm Reflectance Ratio

950nm−to−750nm Reflectance Ratio

1.1

1.05

1

0.95

0.9

0.85

0.8

0.09

0.1

0.11

0.12

0.13

0.14

750nm Reflectance

0.15

0.16

0.17

0.18

0.75
0.08

0.09

0.1

0.11

0.12

0.13

0.14

750nm Reflectance

0.15

0.16

0.17

0.18

Figure 3.14: Comparison of latitudinal trends with and without pixels from Mare
Frigoris. (A) Summary of all 28 mare-based latitudinal trend lines, including datafrom Mare Frigoris. (B) Summary of the 23 mare-based latitudinal trend lines thatremain after Mare Frigoris pixels are excluded. The mean slope diﬀers by less than
0.04 between the two cases (much less than 1 = 0.3). In both cases, the swirl- andimpact-related trends are shown for comparison as solid and dashed lines, respectively.

703.3.2.5

Highlands Contamination

Diﬀusive contamination near mare/highland boundaries and contaminationfrom major crater rays could be partly responsible for the latitudinal trends
we observe if, by chance, such contamination occurs preferentially at higher
latitudes, and if it produces the same relative changes in albedo and the
950 nm/750 nm band ratio. Here, we examine these factors and conclude
that their contributions to the observed latitudinal trend are not significant.
Impact-induced diﬀusive mixing of highland and mare materials occurs
near the mare/highland boundaries, and may be significant within 10–20 kmof the boundary (Mustard et al., 1998). The edges of the maria may therefore
be relatively bright compared with other mare regions of similar underlying
composition. If the morphology of high latitude mare regions is such that theirperimeter length to surface area ratios tend to be larger (e.g., Mare Frigoris isnarrow compared with Mare Imbrium), the amount of highland contamination
may be relatively greater at high latitudes, potentially contributing to thelatitudinal color trends we observe. If true, this would imply that the color
variation we observe for any given compositional bin (e.g., Figure 3.4A) is, tosome extent, a function of each pixel’s distance from the nearest highlands,
with pixels close to the highlands being brighter than pixels that are farther
from the highlands, where contamination is weaker.
We test this hypothesis in two ways. First, we restrict our latitudinalcolor trend analysis to mare regions that are at least 50 km from the nearest
highlands. Mustard et al. (1998) show that diﬀusive mare-highland mixing issignificant only within <10 km of the mare-highland boundary and they findno evidence for mixing beyond 40 km from the boundary. We use the great-71

circle distance (in kilometers) from the center of each mare pixel to the centerof its nearest highland neighbor pixel, where mare/highland surfaces are identified according to our topographic roughness metric (see section 3.3.2.2). Afterexcluding mare regions within 50 km of the nearest highlands, and again applying the selection criteria described in section 3.3.2.1, only 13 compositional
bins remain. The surface area coverage is considerably more limited than in
our full-mare analysis (Figure 3.10), but the resulting latitudinal trend slopes(-1.2 ±0.3) are not significantly diﬀerent (see Figure 3.15A and B and comparewith Figure 3.4). Note that the 50-km threshold we use here is also larger thanthe Lunar Prospector Gamma Ray Spectrometer footprint (⇠ 45 km), helping
to improve the compositional accuracy of the bins.
As a further test, we once again plot 750 nm reflectance versus the 950
nm/750 nm reflectance ratio, except that in this case, we color-code points
not by cosine latitude, but instead by each pixel’s distance from the highlands, from 50 km up to 100 km (Figure 3.15C). No significant color trend isobvious in Figure 3.15C and a least squares best fit line through the data (blackline) predicts only a slight tendency for pixels 50 km from the mare/highland
boundary (white circle) to be brighter than pixels 100 km from the boundary(black circle). This is not to say that mare pixels are not typically brighternear mare/highland boundaries, only that (within a given compositional bin)
this eﬀect is relatively small compared to the eﬀect of varying latitudes (Figure 3.15A). While Figure 3.15C shows just one illustrative example, we repeated this analysis for each of the 13 compositional bins described above.
The resulting highland-distance trends do not resemble the latitudinal trends
and, in some cases, even oppose them (Figure 3.15D). Across the 13 bins, thehighland-distance color trends are often poorly determined (large uncertain72

1.1

0.9

1.05

0.85

1

0.8

0.95

0.75

0.9

0.7

0.85

0.65

0.75
0.7

0.6

15−17 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

0.12

0.14

0.16

750nm Reflectance

0.18

0.2

0.22

950nm−to−750nm Reflectance Ratio

0.5

1.15

95

1.1

90

1.05

85

1

80

0.95

75

0.9

70

0.85

65

0.8

60

0.75
0.7

15−17 wt% FeO
1.0−2.0 wt% TiO2
0.08

0.1

0.12

0.14

0.16

0.18

0.2

0.22

0.95

0.9

0.85

0.09

0.1

0.11

0.1

0.11

0.12

0.13

0.14

0.15

0.16

0.17

0.18

0.12

0.13

0.14

0.15

0.16

0.17

0.18

750nm Reflectance

1.15

D
1.1

1.05

1

0.95

0.9

0.85

0.8

55

750nm Reflectance

1

0.75
0.08

100

C

1.05

0.8

0.55

Distance from Highlands (km)

0.8

B
1.1

Cosine Latitude

0.95

950nm−to−750nm Reflectance Ratio

1.15

950nm−to−750nm Reflectance Ratio

950nm−to−750nm Reflectance Ratio

1.15

1

A

50

0.75
0.08

0.09

750nm Reflectance

Figure 3.15: Color variation of mare surfaces that are at least 50 km from the nearest
highlands (compare with Figure 3.4, which characterizes the entire lunar surface).
(A) Spectral characteristics of pixels sampled from regions with between 15 and17 wt% FeO and between 1 and 2 wt% TiO2 . Points are color-coded according
to cosine latitude, as in Figure 3.4. (B) Compilation of latitudinal trends acrossthe 13 compositional bins compared to the swirl- and impact-related trends (red,green, and blue lines). (C) Spectral characteristics of the same pixels as in (A),except that points are color-coded according to distance from the highlands: blue
points correspond to locations that are 50 km from the mare/highland boundary
whereas red points correspond to locations that are 100 km or more from the nearest
highlands (pixels within 50 km of the nearest highlands are excluded). The black
trend line represents a least squares best fit through the data points with the end
points indicating the color predicted by the best-fit line at locations 50 km (whitecircle) and 100 km (black circle) from the mare/highland boundary. (D) Compilationof the 13 highland-distance trend lines (grey lines) compared to the swirl- and impactrelated trends (red, green, and blue lines).

73

ties), inconsistent (slopes are highly variable, with standard deviation ±3.9),
and generally weak (most of the trend lines are relatively short). For example,
the average change in 750 nm reflectance between 50 km and 100 km away fromthe mare/highland boundary is 0.001 ± 0.004, whereas the average change in
750 nm reflectance between the equator and latitudes of ±70° is 0.05 ± 0.01
(ranges are 1 ). Although the eﬀects of highland contamination cannot be
said to make no contribution whatsoever, we conclude that they cannot be
responsible for the stronger and more consistent latitudinal trends we observe(Figures 3.4, 3.10, 3.14, 3.15A and 3.15B).
Ray systems from major young impacts could also contribute to the highalbedo of Mare Frigoris and other high latitude maria, however, ray systems
similarly cross the low latitude mare regions and should have similar eﬀects
there. For example, rays from Copernicus, Aristarcus, and Kepler all extend
across low latitude mare regions. The young crater Anaxagoras (⇠ 50 km indiameter) is the main contributor to contamination of Mare Frigoris but lies
more than 300 km from its edge, at the closest point. Significant brightening
associated with the similarly young, but larger, Copernicus crater (⇠ 100 km
in diameter) is limited to within ⇠ 300 km, suggesting that Anaxagoras contamination is unlikely to dominate the spectra of Mare Frigoris. Moreover,although contamination from crater ray systems may aﬀect the spectra of certain parts of the maria, there is no reason to expect such contamination tovary systematically with latitude in such a way as to create the latitudinal
trends we report here. Werner and Medvedev (2010) found no concentrationof rayed craters at high latitudes and, due to the present impact rate inhomogeneity (Le Feuvre and Wieczorek , 2011), we would actually expect lesscontamination from young crater ejecta at high latitudes. Crater ray contam74

ination could be substantially responsible for the trends we observe only ifsuch contamination favors high latitudes by chance. Although we cannot ruleout this possibility, we regard it as unlikely, especially given the symmetry of
the latitudinal trends (Figure 3.9) and the consistency of their slopes across a
range of diﬀerent compositions (Figure 3.10, Figure 3.14, Figure 3.15B).
Finally, even if highland contamination were to have a significant eﬀect
on the observed spectral characteristics, the eﬀect would more likely involvean increase (rather than a decrease) in the 950 nm/750 nm reflectance ratio.
As illustrated in Figure 3.6, low-iron regions, such as the highlands, tend toplot farther up and to the right compared with high-iron regions, such as themaria, meaning that spectral trends due to highland contamination shouldhave a positive slope, unlike the negative latitudinal trend slopes we observe
(m1 =

1.5).

3.3.2.6

Phase Angle Biases and LOLA Reflectance

The Clementine spacecraft imaged the lunar surface over a range of phase angles (sun-surface-spacecraft angle), from typically less than 25° at the equator,
to as high as ⇠ 100 ° at high latitudes. The Clementine mosaics are photometrically normalized and calibrated such that they provide an estimate ofreflectance reduced to standard illumination conditions (Eliason et al., 1999),
but the imperfect calibration and imperfect knowledge of photometric behavior, leave small but systematic errors in the mosaics. In particular, the 950 nm
reflectance mosaics exhibit artificially low reflectance values when based on imagery collected at relatively high phase angles (imperfections in the 750 nm
reflectance mosaics are insignificant by comparison). This can be seen bycomparing adjacent parts of the mosaic that are based on imagery collected75

at diﬀering phase angles (Figure 3.16). From equation (3.2), and bearing in
mind that m2 is negative, it follows that this has the eﬀect of making relatively higher phase angle areas artificially low in(and therefore appearing

less like the bright parts and more like the dark parts of swirls). Because
the high latitude Clementine data correspond to systematically higher phaseangles,therefore tends to be underestimated at high latitudes. This is the

opposite of the trend we observe, indicating that the reported latitudinal eﬀect
on color must be even stronger than it appears in Figure 3.4, Figure 3.10 and
Figure 3.12.
The conclusion that phase angle biases are not the source of the observed
latitudinal color variation is further supported by examination of an independent dataset obtained from the Lunar Orbiter Laser Altimeter (LOLA). LOLA
was recently used to measure lunar surface reflectance at 1064 nm (Lucey et al.,2014). Because LOLA uses its own laser light source to illuminate the surface, reflectance estimates are obtained at a constant phase angle of nearlyzero. Figure 3.9B shows LOLA 1064 nm reflectance across the lunar maria,
confirming that, independent of phase angle, higher latitude maria appear to
exhibit higher reflectance, as seen in the Clementine data.

3.4

Discussion

The latitudinal color variation we observe is unlike the color trends associated with impacts but statistically equivalent to those observed at swirls (Figure 3.10, Figure 3.11), suggesting a common mechanism. It has been proposedthat the swirl-related color trend could be the result of magnetic and/or electric field-related alteration of the regolith microstructure (Pieters et al., 2014),
76

14

61

30

65

15
10

10

9

19

20

14

21

22

23

Longitude (degrees E)

24

25

26

63
62

59

61
58

60
59

57

58
57

56

56
55
−12

5

−11

−10

61

−9

−8

−7

Longitude (degrees E)

−6

−5

−4

0.17

0.096

0.165

Clementine 950nm Reflectance

60

0.094

12

0.092

11

0.09

0.088

10

Latitude (degrees N)

Latitude (degrees N)

13

0.086

9

0.16

0.155

59

0.15
58

0.145
0.14

57

0.135
0.13

56

0.125

0.084
8
18

19

20

21

22

23

Longitude (degrees E)

24

25

55
−12

26

14

55

Clementine 950nm Reflectance11

Latitude (degrees N)

Latitude (degrees N)

20

8
18

60

25

12

Clementine Phase Angle (degrees)

Clementine Phase Angle (degrees)

64
13

−11

−10

−9

−8

−7

Longitude (degrees E)

−6

−5

−4

61

0.275

0.12

0.34
0.335

60

11

0.265

10

0.33
0.325

59

0.32

Beta

12

Beta

Latitude (degrees N)

0.27

Latitude (degrees N)

13

58

0.315
0.31

57

0.305

0.26

8
18

0.3

56

9

0.295

19

20

21

22

23

Longitude (degrees E)

24

25

26

55
−12

0.255

−11

−10

−9

−8

−7

Longitude (degrees E)

−6

−5

−4

0.29

Figure 3.16: Phase angle eﬀects on 950 nm reflectance mosaics and derived . Leftpanels show a low-latitude region in Mare Tranquillitatis. Right panels show a highlatitude region in Mare Frigoris. Top panels show the phase angle corresponding to
the imagery collected across the regions. The middle panels show Clementine 950 nm
reflectance, illustrating a bias toward lower reflectance values for imagery collected
at relatively higher phase angles. The bottom panels show the derived parameter, illustrating that the phase angle bias means that is underestimated for higher
phase angle imagery.

77electrostatic sorting of fine dust (Garrick-Bethell et al., 2011), or similar mixing of compositionally distinct materials (Blewett et al., 2011). However, it is
not clear how such mechanisms could also produce the latitudinal color variation we report here. Instead, we propose that reduced solar wind flux isthe most likely common mechanism behind the latitudinal and swirl-relatedcolor trends. Solar wind flux is reduced both with increasing latitude (due tothe increasing incidence angle), and very likely at swirls, where strong magnetic fields may be partially shielding the surface from the impinging solar
wind (Hood and Schubert, 1980; Hood and Williams, 1989; Hemingway andGarrick-Bethell , 2012). Micrometeoroid flux, in contrast, varies with latitude
but should be unaﬀected by magnetic fields (Richmond et al., 2003) and therefore cannot account for the color variations observed at swirls.
The observed latitudinal and swirl-related variations in spectral characteristics (i.e., those described by theparameter) suggest a dependence on

average flux rather than on the total accumulation of solar wind ions. A single
grain saturates after only ⇠ 100, 000 years (Hapke, 2001) such that even a surface shielded by as much as 90% would saturate after only ⇠ 1 million years.
Consistent with this constraint, swirls often exhibit diﬀuse morphologies, with
gradual transitions from the bright parts to darker background soils. In theseregions, mean solar wind flux could vary smoothly between the bright partsof swirls and the darker background soils. Hence, if swirls are the result of
locally altered solar wind weathering, their color must be controlled not by
total accumulation of solar wind ions, but instead by some flux-dependentequilibrium—for instance between regolith gardening and solar-wind-induced
alteration of exposed grain surfaces.
These observations lead us to propose a new model of how space weath78

ering operates for solar wind and micrometeoroids. Over length scales of afew kilometers, the surface exhibits considerable variability in ↵, owing to theappearance of fresh craters (low ↵) among the darker background soils (high↵), whereas variability inis low (Figure 3.1B, Figure 3.2B, Figure 3.3B).

Our interpretation is that the equilibrium

is reached so rapidly that we do

not observe variability in its value over ⇠ kilometer scales, whereas the evolution of ↵ occurs gradually, as impact craters transition into mature soils overlonger timescales. We propose a model in which fresh impact craters begin
with high albedo and low 950 nm/750 nm band ratios (the hypothetical "hyperfresh" point illustrated in Figure 3.17, which should depend only on local
mineralogy), then rapidly evolve to an equilibrium value of(controlled by so-lar wind flux), and finally follow the steeper progression towards larger valuesof ↵, as they gradually mature toward the local equilibrium color. Both therapid weathering to a solar wind flux-dependentand the gradual weathering

toward the local saturation level of ↵ involve decreases in albedo and increases
in the 950 nm/750 nm band ratio. However, the initial rapid weathering process (change in ) has a proportionally greater eﬀect (by a factor of ⇠ 3.6 )on albedo (Figure 3.17). In this model, the trends are parallel, as the soil
matures toward distinct solar wind flux-dependent equilibrium colors, rather
than converging to a common point, contrary to expectations for a simple
mixing process (see also (Staid and Pieters, 2000; Wilcox et al., 2005)).
The diﬀerences in optical eﬀects may reflect diﬀerences in the way solar wind and micrometeoroid bombardment aﬀect lunar soils. For example,
the darkening may be due primarily to the accumulation of nanophase reduced iron (Hapke, 2001; Noble et al., 2007) generated both by solar windand micrometeorites, whereas the increased 950 nm/750 nm band ratio, which
79

)

ing

s
rea

inc

l
dua
gra eringhs
at
we crater
of
oils
s
o
int ow `)
l
(at

l
dua
gra eringhtsa
we cratersofsoil )ot
in igh `ht(a_

950nm-to-750nm Reflectance Ratio

increasing `
lled by solar windflux(controrapid decreashypothetical“hyperfresh”
surface(transition from
fresh craters tomature soils)

e in `

750nm Reflectance

Figure 3.17: Inferred optical evolution process. For a given composition, all impact
craters large enough to excavate fresh material should begin with the same characteristic color (the hypothetical "hyperfresh" point at right). In this model, following
the impact event, the freshly exposed material experiences rapid weathering until itreaches an equilibrium value of , controlled by solar wind flux. This rapid color
change occurs primarily in albedo, but is accompanied by a small change in the
950 nm/750 nm band ratio. Subsequently, the crater material matures gradually(increases in ↵) until it reaches the local saturation maturity level of the surrounding, well-developed soils. This gradual maturation involves additional darkening,but is also accompanied by a proportionally greater band ratio increase than during
the initial mode of rapid weathering.

80

reflects both suppression of spectral features and increased continuum slope,
may also be influenced by impact vitrification (Adams and McCord , 1971) andthe soil’s impact-induced disintegration into increasingly finer grains (Pieters
et al., 1993), neither of which depend on solar wind.
In spite of the uncertainties regarding the precise mechanism for the observed latitudinal trends, our results help to quantify the eﬀects of latitudedependent space weathering, which may need to be accounted for when interpreting spectral measurements at diﬀerent latitudes (Zuber et al., 2012;
Cohen et al., 2014). The production and retention of hydroxyl (OH) groups,
for instance, reportedly varies with latitude (Clark , 2009; Pieters et al., 2009;
Sunshine et al., 2009; McCord et al., 2011; Hendrix et al., 2012) as well as
at swirls (Kramer et al., 2011). The latitudinal variation in space weathering
eﬀects we report here may influence the interpretation of the spectral observations behind such findings.
Our results also suggest that it may be possible to use variations in surface
color to quantify the reduction in solar wind flux at swirls. For example,
the magnitude of color variation at the Reiner Gamma swirl is approximately
equivalent to the color variation observed between the equator and 60° latitude,
suggesting that this magnitude of color variation corresponds to a ⇠ 50%
reduction in solar wind flux. This type of analysis, in conjunction with othermethods, may in turn help in estimating the strength of surface fields at swirls
and in making predictions that can eventually be tested with near surface
magnetic field and solar wind flux measurements.

81

3.5

Conclusions

Our analysis reveals a systematic latitudinal variation in the near-infrared
spectral properties of the lunar surface, and in particular, across the maria.
Specifically, low latitude mare regions tend to be darker and have higher
950 nm/750 nm reflectance ratios than high latitude mare regions, such as thenoticeably bright Mare Frigoris. This latitudinal trend persists across a rangeof distinct compositions, confirming that it is not an artifact of the fact that
regions with the highest iron and titanium content happen to be concentrated
at low latitudes. Our analysis also shows that the trends are not significantly
aﬀected by contamination from the highlands, nor can they be artifacts ofphase angle biases, as confirmed by our comparison with the constant phaseangle LOLA reflectance data. Furthermore, the spectral characteristics of thelatitudinal color trend are statistically equivalent to those observed at lunar
swirls; higher latitude regions appear more like the bright parts of swirls. We
propose that reduced solar wind flux, which should occur both at swirls andtoward higher latitudes, is the common mechanism behind the observed color
variations. We suggest a process by which freshly exposed materials initiallyexperience rapid changes in color until reaching a solar wind flux-dependentequilibrium, followed by a more gradual period of optical maturation driven
mainly by micrometeoroid impacts. If correct, this model could help quantify
the distinct eﬀects of solar wind versus micrometeoroid weathering and may
be important for the interpretation of spectral observations made at diﬀerent
latitudes on the Moon and other airless bodies throughout the solar system.

82

Acknowledgements
This work was supported by the BK21 plus program through the National
Research Foundation (NRF), funded by the Ministry of Education of Korea.
Partial support was provided by the UCSC/NASA Ames University AﬃliatedResearch Center, Aligned Research Program; the Alfred P. Sloan Foundation;
and NASA grant NNX13AJ51G. We thank Carle Pieters for important discussions and advice, Paul Lucey for supplying the LOLA reflectance data, andthe two reviewers for their helpful comments and suggestions.

83

Chapter 4
A Rigid and Weathered Ice Shellon TitanThis chapter is adapted from Hemingway, D., F. Nimmo, H. Zebker, and L. Iess(2013), A rigid and weathered ice shell on Titan, Nature, 500 (7464), 550–552,
doi:10.1038/nature12400.

Abstract
Several lines of evidence suggest that Saturn’s largest moon, Titan, has a
global subsurface ocean beneath an outer ice shell 50 to 200 kilometers thick(Béghin et al., 2010; Bills and Nimmo, 2011; Iess et al., 2012; Tobie et al.,2006). If convection is occurring (Mitri and Showman, 2008; Tobie et al.,2005), the rigid portion of the shell is expected to be thin; similarly, a weak,
isostatically-compensated shell has been proposed to explain the observed topography (Choukroun and Sotin, 2012; Nimmo and Bills, 2010). Here we
report a strong inverse correlation between gravity (Iess et al., 2012) and to84

pography (Zebker et al., 2012) at long wavelengths that are not dominated bytides and rotation. We argue that negative gravity anomalies (mass deficits)
produced by crustal thickening at the base of the ice shell overwhelm positive
gravity anomalies (mass excesses) produced by the small surface topography,
giving rise to this inverse correlation. We show that this situation requires a
substantially rigid ice shell with an elastic thickness exceeding 40 kilometers,
and hundreds of meters of surface erosion and deposition, consistent with recent estimates from local features (Moore et al., 2013; Neish et al., 2013). Our
results are therefore not compatible with a geologically active, low-rigidity ice
shell. After extrapolating to wavelengths that are controlled by tides and rotation, we suggest that Titan’s moment of inertia may be even higher (that is,
Titan may be even less centrally-condensed) than is currently thought (Iesset al., 2010).

4.1

Introduction

Combined studies of gravity and topography can yield useful information about
the near-surface interior structure of solid planetary bodies (Wieczorek , 2007).
In particular, the ratio between the gravity and topography at a particular
wavelength, known as the spectral admittance, is generally a function of the
thickness of the lithosphere and the degree to which the topography is compensated by subsurface density anomalies. For example, small admittances arise
when topography is compensated isostatically such that the additional gravity
associated with surface topography is oﬀset by the reduced gravity associatedwith the corresponding low-density isostatic roots. Conversely, large admittances usually indicate incomplete isostatic compensation. Large admittances85are typical at short wavelengths, for which elastic flexure of the lithosphere isan eﬀective means of supporting the topography. The wavelength at which theprimary support mechanism transitions from flexural to isostatic is generally
a function of the eﬀective elastic layer thickness, which, in turn, tells us about
the thermal structure of the lithosphere.
In the case of Titan, our limited knowledge of the gravity field restricts ouranalysis to the longest of wavelengths (mainly degrees 2 and 3). The degree2 gravity field is well determined, but is dominated by the eﬀects of tidal
and rotational distortion. The gravitational eﬀects of isostatically supported
topography at degree 2 are small by comparison and, in any case, diﬃcult
to distinguish from the tidal and rotational eﬀects. In this analysis, we will
therefore focus on the degree-3 gravity and topography. The degree-3 signal
is useful because, although it is small compared to those at degrees 2 and4, it is not directly aﬀected by tides, rotation, or tidal heating (Nimmo andBills, 2010). Though our analysis will focus primarily on degree 3, we will also
discuss the implications of extrapolating our results to degree 2.
Titan is thought to consist of a large, rocky inner core, surrounded bylayers of H2 O in various phases, and covered by a dense nitrogen atmosphere.
The outer-most layer of the H2 O mantle is thought to consist primarily of ice-I
and estimates for its thickness range from 50 to 200 km (Béghin et al., 2010;
Bills and Nimmo, 2011; Iess et al., 2012; Tobie et al., 2006). The surface isgeomorphologically varied, with seas of liquid hydrocarbons, dendritic river
networks, and large fields of sand dunes. Several surface features have been
described as icy volcanoes, suggesting ongoing geologic activity in the ice shell.
However, evidence for such “cryovolcanism” is not entirely conclusive and there
is considerable disagreement about the present level of geologic activity on
86Titan. Obtaining constraints on the loading history and elastic thickness of
the ice shell will help in understanding the evolution of Titan’s surface andthe likely level of ongoing geologic activity (Moore and Pappalardo, 2011).
Section 4.2 describes the observed gravity and topography, as well as theimplied admittance. In section 4.3, we describe our model and use it to interpret the observed admittance, leading to estimates for the ice shell’s elastic thickness and the magnitude of surface erosion that has taken place. In
section 4.4, we discuss the implications of our results before concluding insection 4.5.

4.2

Observations

4.2.1

Gravity Field

The description of Titan’s gravitational field is given by Iess et al. (2012) as
g
g
, Slm
such that the
non-normalized, dimensionless potential coeﬃcients Clm

gravitational potential at the reference radius (R) can be written

U (✓, ) =1l

GM X X gg(C cos m + Slm
sin m ) Plm (cos ✓)R l=0 m=0 lm

where Plm (cos ✓) are the non-normalized associated Legendre functions, ✓ iscolatitude,gg
and Slm
are non-normalized spherical haris longitude and Clm

monic coeﬃcients of degree l and order m. In terms of the non-normalized
coeﬃcients, the radial component of the gravitational acceleration at the ref-87

erence radius (R) can be written

gr (✓, ) =1l

GM X Xgg(l + 1) (Clmcos m + Slm
sin m ) Plm (cos ✓)R2 l=0 m=0In our admittance analysis (section 4.2.3), we will use fully normalized,gg, S̄lm, such that
dimensionalized gravitational acceleration coeﬃcients, C̄lmgr (✓, ) =

l
1 X
X

g
g
C̄lm
cos m + S̄lm
sin mP̄lm (cos ✓)l=0 m=0

where P̄lm (cos ✓) are the fully normalized associated Legendre functions (seeAppendix A). The fully normalized and dimensionalized gravity coeﬃcients,gg, S̄lm, are related to the non-normalized, dimensionless potential coeﬃC̄lm
g
g
, Slm
by
cients Clm

C̄, S̄

g
=
lm

✓

(2

(l m)!
0m )(2l + 1)
(l + m)!

◆ 12(l + 1)GM{C, S} glm
2
R

(4.1)

where the square root term does the normalization, the (l + 1) term arises
from the diﬀerentiation associated with converting from potential to gravity,
term generates dimensional coeﬃcients (which we will express in terms
the GMR2
of mGal = 10 5 ms 2 ), and where

0m is the Kronecker delta.

Three solutions for Titan’s low-order gravity field have been obtained via
radio tracking of the Cassini spacecraft (Iess et al., 2010, 2012), with the results being well constrained up to spherical harmonic degree 3. In the firsttwo (SOL1a and SOL1b), results from six Cassini gravity flybys were analyzed separately and then combined into multi-arc solutions. Whereas SOL1a
attempts to model the gravity field only up to degree 3, SOL1b attempts to88

model the field up to degree 4, but only as a means of verifying the robustness of the degree-3 solution. It is found that the degree-3 solution diﬀers
only modestly between SOL1a and SOL1b. For the last solution (SOL2), a
global model extending to l = 3 was derived using Pioneer and Voyager dataand satellite ephemerides in addition to Cassini observations. In spite of thediﬀerent approaches, the SOL2 field closely resembles the SOL1a field. The
SOL1b solution matches less closely, likely because it also attempts to include
the degree-4 component of the field. Although all three solutions give a consistent estimate of the periodic (tidal) k2 Love number (with a central valueof 0.6), the static part of the degree-4 field is currently not well constrained
due to the limited nature of the observations. Our calculations will be based
primarily on the SOL1a gravity field.

4.2.2

Topography

Titan’s topography is likewise described by Zebker et al. (2012) in terms of
h
h
, Slm
and can be written as
non-normalized spherical harmonic coeﬃcients Clm

H(✓, ) =

l
1 X
X

h
h
Clm
cos m + Slm
sin m

Plm (cos ✓)l=0 m=0

where, again, Plm (cos ✓) are the non-normalized associated Legendre functions,✓ is colatitude,hhand Slm
are non-normalized sphericalis longitude and Clm

harmonic coeﬃcients of degree l and order m. To convert from non-normalizedhh, S̄lm, we write
coeﬃcients to fully normalized coeﬃcients, C̄lm

C̄, S̄

h
lm =

✓

(2

(l m)!
0m )(2l + 1)
(l + m)!

89◆ 12{C, S}hlm(4.2)

and the topography is described (in meters) byH(✓, ) =

1 X
l
X

h
h
C̄lm
cos m + S̄lm
sin mP̄lm (cos ✓)l=0 m=0

where P̄lm (cos ✓) are the fully normalized associated Legendre functions (seeAppendix A).
The topography models are derived from a combination of radar altimetry
and analysis of the overlapping regions of radar images, in a technique known as
SARtopo (Zebker et al., 2009; Stiles et al., 2009; Zebker et al., 2012), and wereused to derive shape solutions up to l = m = 11. Several distinct solutions were
produced, depending on where the harmonic expansion was truncated. We
denote these solutions Deg4-exp, Deg5-exp, ... Deg11-exp, where the number
indicates the highest degree and order used to fit the observations.
Due to the large gaps in Cassini radar coverage (Figure 4.1), topographymodels with power beyond degree 6 are not adequately constrained unless
an a priori restriction is applied (e.g., minimize rms deviation from best-fit
sphere). Even when a priori constraints are applied, the coeﬃcients tend to
be less stable when the model’s expansion limit exceeds degree 6 (Figure 4.2).
For our purposes, we prefer to use the highest resolution data available
without requiring a priori constraints in the model fits. We therefore primarily
use the Deg6-exp model (Zebker et al., 2012) in our analysis, the coeﬃcientsfor which are given in Table 4.1.

90

600

60

400

30

200

0

0

−30

−200

−60

−400

−90

Elevation (m)Latitude (degrees N)

90

−600
0

90

180

270

360

Longitude (degrees E)

Figure 4.1: Cassini radar-derived elevation data for Titan. Elevation is given relative
to the 2575-km reference sphere.
40

Coefficient value (m)

20

0

−20

−40

−60

C30
C31
S31
C32
S32
C33
S33

−80
3

4

5

6
7
Model expansion limit

8

9

10

Figure 4.2: Degree-3 normalized topography model (Zebker et al., 2012) coeﬃcients(with 1- error bars) as a function of the maximum spherical harmonic degree allowed
when fitting the data.

91

Table 4.1: Fully normalized Deg6-exp topography model (Zebker et al., 2012) coefficients (in meters).
Term
C00
C10
C11
S11
C20
C21
S21
C22
S22
C30
C31
S31
C32
S32
C33
S33
C40
C41
S41
C42
S42
C43
S43
C44
S44

Estimate
2574750.0
0.0
23.1
18.5
-169.5
-17.8
31.8
120.8
20.1
-4.9
2.8
-21.3
-8.8
-2.9
0.0
14.3
-39.3
22.1
75.9
4.5
-26.8
16.7
-50.2
0.0
-47.3

1 error
6.0
5.8
5.2
5.8
4.5
5.4
7.0
4.6
4.6
3.8
3.7
5.6
5.9
5.9
7.2
7.2
3.7
3.2
4.2
4.5
4.5
0.0
0.0
0.0
0.0

Term
C50
C51
S51
C52
S52
C53
S53
C54
S54
C55
S55
C60
C61
S61
C62
S62
C63
S63
C64
S64
C65
S65
C66
S66

92

Estimate
28.6
-26.9
-14.0
37.1
-12.4
0.0
0.0
0.0
0.0
0.0
0.0
-2.5
8.9
-61.0
-24.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

1 error
3.3
2.3
4.7
6.2
6.2
0.0
0.0
0.0
0.0
0.0
0.0
3.1
2.5
3.8
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

4.2.3Observed Admittance

When we compare the degree-3 gravity and topography signals, a striking
negative correlation is apparent (Figure 4.3). We quantify this relationship by
employing an admittance analysis (McKenzie, 1994; Wieczorek , 2007), whichmeasures the gravity-to-topography ratio in a way that is not biased by uncorrelated signals in the gravity data. The degree-l admittance, Z (l), is givenbyZ(l) =

Dhg (l)Dhh (l)(4.3)

where Dij (l) represents the cross-power spectrum between fields i and j andis given by
Dij (l) =

l
X

j
j
i
i
Clm
Clm
+ Slm
Slm

(4.4)

m=0

The correlation between the two signals is given by
Dhg (l)(l) = pDhh (l)Dgg (l)(4.5)

If some fraction of the gravity signal is not correlated with the surface topography, then the coherence ( 2 ) will be less than one. However, the crucial
advantage of equation (4.3) is that any such gravity noise does not aﬀect theestimated admittance, Z. For the case of Titan, contributions to gravity fromdeeper interfaces (such as the silicate interior) are likely to be important, while
contributions to the surface topography from these processes are likely negligible. An approach like that embodied in equation (4.3), which is unaﬀectedby noisy gravity, is essential for interpreting the limited observations available
at Titan.
For each of the three gravity models, we computed admittances based on
93

(a)

(b)
Figure 4.3: Titan’s degree-3 gravity and topography. (a) Gravity field derived fromthe potential coeﬃcients of the SOL1a gravity field representation of Iess et al. (2012)(multi-arc analysis, 3 ⇥ 3 gravity field). (b) Spherically referenced topography basedon degree-6 harmonic expansion (Zebker et al., 2012). The two signals display astrong negative correlation (with correlation coeﬃcient = 0.61) and give rise toan admittance of 32 mGal/km, based on a Monte Carlo analysis accounting for theuncertainty in the two signals. Maps are shown in Mollweide projection, centeredon the anti-Saturnian point (180° longitude).

94Topography

Degree−3 Admittance (mGal/km)
−40

−20

0

−24±9 (−0.51)−13±10 (−0.30)

Deg4−exp

Gravity

+20

SOL1a
SOL1b

−22±10 (−0.43)SOL2−19±8 (−0.45)

SOL1a−8±9 (−0.21)

Deg5−exp−15±8 (−0.33)−32±11 (−0.61)−16±12 (−0.34)

Deg6−exp−27±12 (−0.49)

SOL1b
SOL2
SOL1a
SOL1b
SOL2

Figure 4.4: Admittance estimates for nine sets of gravity and topography data.
Three distinct gravity models (Iess et al., 2012) and three distinct topography models(Zebker et al., 2012) (truncated at harmonic degrees 4, 5 and 6) were used to produce
a total of nine separate admittance estimates. Each estimate is based on a Monte
Carlo analysis in which the admittance was computed for each of 100,000 distinct setsof gravity and topography coeﬃcients, distributed according to the 1 uncertainties
in the model coeﬃcients. For each Monte Carlo analysis, diamonds are plottedat the mean admittance estimate, and 1 error bars are shown to represent the
distribution. To the right of the error bars are the mean ± standard deviationof the admittance estimates in milligals per kilometer (mGal/km), followed by thecorrelation coeﬃcients ( ) in parentheses.

95three separate topography models (Figure 4.4). The admittance and correlation estimates in Figure 4.4 were obtained through a Monte Carlo analysis.
For each gravity/topography combination, the results were obtained from a
distribution of 100,000 distinct admittance and correlation estimates, each ofwhich was based on gravity and topography coeﬃcients that were generated
randomly and distributed according to the 1- uncertainties in the model coeﬃcients. Individual correlations and admittance estimates were computed
according to equations (4.3) and (4.5). In spite of the variability between
the estimates, the degree-3 admittance appears to be substantially negative,
and the Monte Carlo analysis shows that this result is robust to the model
uncertainties (Figure 4.4).
The admittance estimated based on the Monte Carlo analysis will have aslightly smaller magnitude than the admittance estimated directly from thecoeﬃcients (i.e., when uncertainties are ignored). This is because, as long as
there is uncertainty in the topography coeﬃcients, the mean of the distributionof Dhh (see equation 4.4) will always be greater than the value of Dhh obtained
directly from the estimated topography coeﬃcients (because if x is normally
distributed, then E(x2 ) > [E(x)]2 ). For example, if uncertainties are ignored,
the admittance computed directly from the SOL1a gravity (Iess et al., 2012)
coeﬃcients and the Deg6-exp topography (Zebker et al., 2012) coeﬃcients is
39 mGal/km, whereas when uncertainties are accounted for using a Monte
Carlo analysis, the mean estimated admittance is

32 mGal/km. We adopt

the latter value because it is more conservative—more negative admittances
would require higher magnitudes of erosion and/or larger elastic thicknesses(see section 4.3.2).

96

4.3

Model and Results

The negative correlation we observe between the gravity and topography signals (section 4.2.3) is surprising—negative admittances are rare. They canoccur in convecting systems with strong viscosity stratification (Richards andHager , 1984), but it is unclear why Titan’s ice shell would have such layering,and the shell is suﬃciently thin that—as for Enceladus (Roberts and Nimmo,2008)—convective features should have much shorter horizontal length scales
than degree 3. Below we construct a model that allows for negative admittances under certain conditions.

4.3.1

Admittance and Flexure Model

Here, we will assume that the topography is supported by some combinationof shell thickness variations (Airy isostasy, e.g., Nimmo and Bills (2010)) andelastic flexure. An admittance model requires that we describe how the gravity(g) is aﬀected by topography (h).
Gravity Anomalies
The gravity anomaly at degree l due to a thin surface layer of amplitude hl
and density ⇢c is given by (Jeﬀreys, 1976)glt =(l + 1)
4⇡Ghl ⇢c
(2l + 1)

97

(4.6)

and similarly, the gravity anomaly due to a thin layer (a "root") of thickness
rl and density contrast⇢ = ⇢mglb =⇢c at the base of the shell, is given by✓(l + 1)
4⇡Grl ⇢ 1
(2l + 1)

d
R

◆l+2

(4.7)

where the mean thickness of the shell is d, the radius of the body is R, and ⇢mis the density of the material underlying the shell (i.e., the subsurface ocean).
In the short-wavelength limit (l1), equation (4.6) reduces to the usual

flat-plate formula, as required.
When the net gravity anomaly and the surface topography (hl ) are known,
the theoretical admittance is given by

Zl =

glt +
hl

glb

(4.8)

In the remainder of this development, we drop the subscripts from both hand r and take it as understood that these parameters correspond to a specificwavelength.
In practice, it will be diﬃcult to observe r and therefore to computeglb

according to equation (4.7). Instead, we would like to find an expression for r
in terms of h, which can be more readily observed. This is generally possible
because, for a finite elastic thickness, there will be a balance between the
overburden pressure of positive surface topography (⇢c gh), the buoyancy of
the root ( ⇢gr), and the restoring forces due to flexure.

98Lithospheric Deflection, Cartesian case
If the deflection of the initial ice shell is w, in a Cartesian system, and if geoid
undulations are neglected for the moment, this pressure balance can be writtenas
Dr4 w =

⇢gr

⇢c gh

(4.9)

Here, D represents flexural rigidity and is given by

D=

ET 3
12(1 ⌫ 2 )(4.10)

where T is the eﬀective elastic thickness, E is Young’s modulus and ⌫ is
Poisson’s ratio. We treat w as positive upward bending, h as positive upward
relief above the reference ice shell surface, and r as positive downward relief
from the base of the ice shell. The relationship between r and h depends onthe elastic properties of the shell and the thickness of loads applied at the top(dt ) and bottom (db ) of the shell (Figure 4.5); dt is the thickness of material
added at the surface (a negative value would indicate erosion), and db is the
thickness of material added at the base of the ice shell (a positive value would
indicate basal freezing, a negative value, basal melting). Our model represents
the equilibrium state achieved after the lithosphere has finished deflecting in
response to the applied load(s). The model also assumes that the ice shellproperties do not change over time.
Our formulation is similar to those of Turcotte et al. (1981) and McGovern
et al. (2002), however, our sign convention diﬀers slightly and, for simplicity,
we assume that material added at the top or bottom of the shell is also of
density ⇢c . Our formulation also diﬀers in that we handle top and bottom loads

99

simultaneously with w being the total deflection resulting from the combined
eﬀects of top and bottom loading.
h

w

dt

d

w

r

db

Figure 4.5: Illustration of the influence of top loading (pale pink material, dt ) andbottom loading (pale blue material, db ) on ice shell flexure (w), surface relief (h) androot thickness (r).

From Figure 4.5, we haveh = w + dt

(4.11)w(4.12)

r = dbThen, assuming the loads are periodic and in-phase, we can solve (4.9) forw, obtainingw=⇢db ⇢c dt⇢m + µ(4.13)

Here, we have introduced a parameter, µ, which has units of density andwill serve as a shorthand for the flexural rigidity at a particular wavelengthand gravity. In a Cartesian system,µ(k) =ET 3 k 4
12(1 ⌫ 2 )g100

(4.14)

where k is a wavenumber. The advantage of using this shorthand will become
clear when we move from a Cartesian to a spherical system.
It is useful to define a compensation function, C(⇢), that expresses thedegree of compensation under flexural support compared with the case of pureisostasy. This can be defined as the ratio of the deflection, w, according to(4.13) to the zero-rigidity deflection, w0 , obtained from (4.13) when µ = 0.
That is, C = w/w0 , orC(⇢) =

1
1 + µ⇢(4.15)

When the elastic thickness, T , is zero, C = 1 (fully compensated). The
parameter ⇢ is the density contrast that is resisting the flexure (i.e., related
to buoyancy, overburden pressure, or both). In the isostatic limit,⇢c h, whereas in the top loading case,⇢r =⇢r = ⇢c hC( ⇢), while in the bottom

loading case, C(⇢c ) ⇢r = ⇢c h, as we will see. The theoretical value of C = 0
corresponds to the zero compensation case which occurs when the ice shellis infinitely rigid (i.e., as µ ! 1). In this case, deflection (w) becomes zero(equation 4.13) and so, from (4.11) and (4.12), h = dt and r = db . In this
scenario, h and r are independent of one another and so both dt and db mustbe specified in order to predict admittance. However, as long as C > 0, there
will be some finite deflection and it will be possible to obtain r as a functionof h.
If C > 0 and both h and dt are specified, then from (4.11), (4.12), (4.13)
and (4.15), it can be shown that"#dt⇢c h 1 dht+
r=⇢ C(⇢c )h101(4.16)

We have factored out ⇢c h/ ⇢ in order to facilitate direct comparison with
the isostatic case and because it will be convenient to do so when calculating
admittance using (4.6), (4.7) and (4.8). In the case where no material has been
added to the surface (i.e., dt = 0, so that loading is purely from the bottom),this expression reduces toC(⇢c ) ⇢r = ⇢c h

(4.17)

Similarly, if C > 0 and only h and db are specified, it can be shown that
"
db
⇢c h
h
r=
⇢ C(⇢c )#db+ 1 C( ⇢)h(4.18)In the case where no material has been added to the base (i.e., db = 0, sothat loading is purely from the top), this expression reduces to⇢r = ⇢c hC( ⇢)(4.19)

Lithospheric Deflection, spherical case
The foregoing gives correct values for r in the Cartesian case, which is anappropriate approximation for short wavelength loads. However, in order tointerpret admittance at very long wavelengths, we must consider the sphericalcase. Assuming the icy crust behaves as a thin elastic shell of radius R,
equation (4.9) becomesDr6 w + 4Dr4 w + ET R2 r2 w + 2ET R2 w
= R4 (r2 + 1⌫) ( ⇢gr

102⇢c gh + ⇢m ghg )(4.20)

where D is as in (4.10). Here, we have adopted a modified version of theapproaches of Turcotte et al. (1981) and McGovern et al. (2002), which themselves follow the earlier derivation of Kraus (1967). The final term in (4.20)
accounts for the elevation or depression of the geoid (hg , which we treat as
positive upward) that occurs with loading of the ice shell. Here, we adopt an
approximation in order to obtain hg , namely, we assume that 1

d l+2
⇡ 1,
R

and that the geoid and gravitational acceleration do not change with depth
in the shell (as has been pointed out by Belleguic et al. (2005), this was animplicit assumption of Turcotte et al. (1981)).
Having obtained hg , and using (4.11) and (4.12), we rewrite (4.20) as
2

✓

4 1

3⇢m
(2l + 1)⇢◆ 1

0 T 2 r6 +4r4
1
3
(
)
2
ET @ R2 12(1 ⌫ 2 ) + r + 2 A+ ⇢m 5 w =R2 g
r2 + 1 ⌫⇢db⇢c d twhere ⇢ is the mean density of the body. If w is expressed in spherical harmonics, we can replace r2 withl(l + 1) (e.g., Turcotte et al. (1981), equations 17

and 18) and solve for w, recovering equation (4.13),w=⇢db ⇢c dt⇢m + µbut with the flexural rigidity parameter now beingµ(l) =✓

13⇢m(2l + 1)⇢◆ 1

0 T 2 l3 (l+1)3 4l2 (l+1)2[]+ l(l + 1)
ET @R2 12(1 ⌫ 2 )R2 gl (l + 1) (1 ⌫)2

1

A (4.21)

For a spherical system, it is also necessary to account for the ratio of
surface areas at the top and bottom of the shell (since the buoyancy of theroot depends on its volume, not its thickness). This eﬀect complicates the103derivations but it can be shown that, if dt and h are specified, then the surface
area correction leads to"#✓
dt⇢c h 1 dht+
1r=⇢ C(⇢c )
h

d
R

◆ 2

(4.22)

If db is specified rather than dt , we obtain
"
db
⇢c h
h
r=
⇢ C(⇢c )

#
db
+1 C
h

✓

⇢ 1

d
R

◆2 ! ✓

d
R

1

◆ 2

(4.23)As required, both of these expressions reduce to the Cartesian equivalents
as R ! 1. This correction also partially relaxes the simplifying assumption
we made to obtain the geoid height so that we now assume 1
d l+2
⇡ 1.
R

than 1d l⇡ 1 rather
R

Having obtained an expression for r as a function of h (which is possible
as long as C > 0), we can now substitute (4.22) into (4.7) and combine with
(4.6) and (4.8) to get an expression for admittance that depends on h and dt ,but not r"(l + 1)
4⇡G⇢c 1
Z(l) =(2l + 1)

1 dht
dt
+
C(⇢c )
h

!✓

1

d
R

◆l #

(4.24)

Equation (4.24) implicitly accounts for the root thickness (r), bottom load
thickness (db ), and mantle density (⇢m ), such that these terms do not appear
in the final expression. Model admittance is sensitive to mantle density only
insofar as the ratio ⇢m /⇢ influences the geoid, the eﬀect of which is capturedin C(⇢c ) via (4.21) and (4.15).
If db is specified rather than dt , we instead substitute (4.23) into (4.7) and104

obtain"(l + 1)
4⇡G⇢c 1
Z(l) =(2l + 1)db
hC(⇢c )

!
db
+1 C
h

✓
⇢ 1

d
R

◆2 ! ✓◆l #d
1
R(4.25)

When C = 1 (pure isostasy), both of these expressions reduce to"(l + 1)
Z(l) =
4⇡G⇢c 1(2l + 1)

✓

d
R

1

◆l #

(4.26)

d/R)l must always be lessand admittance will always be positive since (1than 1.

Finally, in the limit of an infinitely rigid shell (C = 0), for which we mustspecify both dt and db , it can be shown that the admittance is"(l + 1)
4⇡G⇢c 1
Z(l) =(2l + 1)db ⇢d t ⇢c✓1dR◆l+2 #(4.27)

In general (i.e., when 0 < C < 1), admittance depends on topography (h),the amount of top and/or bottom loading (dt and/or db ), as well as the elastic
thickness, T , and the mean shell thickness, d. For a given wavelength, meanshell thickness (d), elastic thickness (T ), and a fixed, positive h, admittanceis a positive linear function of dt (Figure 4.6) and crosses zero when dt =
✓
⌘ 1⇣
⌘◆
⇣
l
1
d
1
1 . When both top and bottom loading
1 R
h 1
C(⇢c )
have taken place, and if h is known, admittance may be either positive or
negative, and is uniquely defined if either dt or db is specified.
Based on equation (4.25), it can be shown that, when loading occurs purely
from the top (i.e., db = 0), admittance is independent of surface relief, h, and105

is given by"(l + 1)
4⇡G⇢c 1
Z(l) =(2l + 1)

✓

C

⇢ 1

d
R

◆2 ! ✓

which is necessarily positive, since both C and (1

1

d
R

◆l #

(4.28)

d/R)l are always less than

1. Based on (4.24), when loading occurs purely from the bottom (i.e., dt = 0),
admittance is, again, independent of surface relief, h, and is given by

Z(l) ="

1
C(⇢c )(l + 1)
4⇡G⇢c 1(2l + 1)

which is positive as long as C(⇢c ) > (1

✓

1

d
R

◆l #

(4.29)

d/R)l . For degree 3, this is always

true for the parameters given in Table 4.2.
Hence, from Figure 4.6, we see that degree-3 admittance may be negative
only if dt /h is negative (i.e., when erosion has occurred at topographic highs)
or if substantially diﬀerent parameter values are adopted. To obtain a negative
degree-3 admittance with pure bottom loading would require an increase of⇠ 30% in the ratio of Young’s modulus to the ice shell density (E/⇢c ). It isalso possible to obtain a negative admittance without erosion for suﬃciently
large elastic thicknesses and suﬃciently short wavelengths (e.g., for l = 6 whenT > 350 km, or for l = 9 when T > 200 km).

4.3.2

Interpretation

Having developed our admittance model, and having seen that negative top
loading (that is, erosion at topographic highs) is required in order to give riseto negative admittance at degree 3, we can now use equation (4.24) to see howadmittance varies as a function of top load (dt ) and the elastic layer thickness
106

100

80

Admittance (mGal/km)

60

40

20

0

−20

−40

−60

−80

−100
−1500

−1000

−500

0

500

1000

1500

Top load, dt (m)Figure 4.6: Admittance as a function of top load (dt , where negative dt indicateserosion) when h = 66 m, T = d = 200 km and assuming the properties listed inTable 4.2.

Table 4.2: Parameter values assumed for admittance calculations

Parameter
Poisson’s ratio for ice
Young’s modulus for ice
Crustal (ice shell) density
Mantle (subsurface ocean) density
Titan’s mean density
Titan’s radius
Acceleration due to gravity at the surface

107

Symbol
⌫
E
⇢c
⇢m
⇢
R
g

Assumed Value
0.25
9 GPa
920 kg/m3
1000 kg/m3
1880 kg/m3
2575 km
1.35 m/s2

(T ). This allows us to constrain the relationship between elastic thickness andthe amount of erosion/deposition that must have taken place. Figure 4.7a
illustrates the admittance corresponding to Figure 4.3 and shows that morethan 100 m of surface erosion is required even for very large (⇠ 400 km) elastic
thicknesses. If the elastic thickness is less than 40 km, more than ⇠ 1 km of
erosion is required. Figure 4.7b plots combinations of d, T and dt that satisfy
the observed admittance. For an ice-shell thickness of 200 km, for example,
more than 200 m of erosion are required.
Our results suggest that the negative admittance we observe at degree-3
is the result of negative gravity anomalies from large roots dominating overthe positive gravity anomalies from the associated topography. We tested
this scenario by computing, everywhere over the surface of Titan, the gravity
anomaly implied by the observed topography and then comparing the resultwith the observations (Figure 4.8c). The gravity anomaly is obtained by multiplying the observed degree-3 topography (Figure 4.8a) by equation (4.24),
assuming T = d = 200 km and a degree-3 erosion amplitude of 293 m (i.e,293 m of erosion at the topographic peaks and 293 m of deposition in the valleys). This is the amount of erosion required to produce

39 mGal/km, theadmittance obtained directly (i.e., neglecting uncertainties) from the SOL1a
gravity (Iess et al., 2012) and the Deg6-exp topography (Zebker et al., 2012)(see section 4.2.3). Figure 4.8b shows the resulting gravity anomaly, computed
everywhere over the surface. The gravity field predicted through this procedure resembles the observed field (compare panels (b) and (c) in Figure 4.8).
Assuming the mantle density given in Table 4.2, the implied root thickness
amplitude is ⇠ 1.4 km.
As illustrated in Figure 4.6, admittance is approximately a direct linear
108

100

400

80
60

300
40
250

20

200

0
−20

150

−40
100
−60
50

−80

0
−1500

(a)

Admittance (mGal/km)

Elastic Shell Thickness, T=d (km)

350

−1000

−500

0

500

−100

Top load, dt (m)
0

500

00
−4

00
−6

00
−8

−160
−14 0
−1200
0
−10 0
00

600

−200
0
20

−400

400
−600
300

−800
−1000

Top load (m)Shell Thickness, d (km)

−

200
−1200
100
−1400
0

(b)

0

50

100

150

200

250

300

350

400

−1600

Elastic Thickness, T (km)Figure 4.7: Model predictions of admittance and erosion. (a) Degree-3 admittancepredicted by our model for various combinations of elastic thickness T and top load
dt , here assuming that T equals the shell thickness d. Larger values of T correspondto smaller values of C (that is, less compensation). Negative top load indicates erosion at topographic highs and deposition at topographic lows. The solid black line
indicates the admittance corresponding to Figure 4.3 ( 32 mGal/km); the dashedlines indicate the 1 uncertainty on that estimate. Whereas the rheology of ice
implies that T < 0.5d (section 4.4.1), assuming that T ⇡ d leads to conservative
estimates of erosion. (b) Top load amplitude, dt , required to produce the observed
degree-3 admittance for various combinations of total shell thickness d and elastic
thickness T . Dashed white lines highlight the likely range of Titan’s ice shell thickness, 50–200 km. Both (a) and (b) were generated assuming a degree-3 topographic
amplitude h of 66 m.

109

(a)

(b)

(c)
Figure 4.8: Degree-3 topography and gravity maps: (a) Deg6-exp topography (Zebker et al., 2012); (b) gravity computed as described in section 4.3.2 (assuming
T = d = 200 km and 293 m of erosion); (c) SOL1a gravity (Iess et al., 2012).

110

function of top load (dt ). Conversely, the top load required to produce a given
admittance can be obtained by solving equation (4.24) for dt :dt = h 1✓

1
C(⇢c )

1◆ 1 "✓

1

Z(l) (2l + 1)
4⇡G⇢c (l + 1)

◆✓
1

d
R

◆ l

1

#!

(4.30)

Based on (4.30), Figure 4.9 shows how the implied erosion (negative dt )varies with the estimated admittance given various combinations of shell thickness (d) and elastic thickness (T ). The dashed black line corresponds toZ(3) =

32 mGal/km, the admittance estimate obtained from the Monte

Carlo analysis (i.e., accounting for uncertainties) based on the SOL1a gravity(Iess et al., 2012) field and the Deg6-exp topography (Zebker et al., 2012) solution. This is also the admittance assumed in generating Figure 4.7b. The±11 mGal/km uncertainty in that admittance estimate (Figure 4.4) translates
to ±81 m uncertainty in the erosion estimate when T = d = 200 km. Diﬀerent
combinations of T and d lead to slightly diﬀerent uncertainties, but roughly±30% is typical.
Although, as we will argue in section 4.4.1, the ice shell is not likely to be
entirely elastic, adopting T = d leads to more conservative estimates of themagnitude of erosion. For example, the magnitude of erosion required to giverise to Z(3) =32 mGal/km is ⇠ 241 m when we assume T = d = 200 km

and ⇠ 577 m when we assume T = 100 km, d = 200 km.

111

500

0
T=200km

Top load, dt (m)

−500

−1000

d=200km

0km
d=10
0km
km
T=10
=200
d
m
0k
T=10

0km

10

−1500

0km

d=

5

T=
−2000

d=

m

0k

20

km

50

T=

−2500

−3000
−70

−60

−50

−40

−30

−20

−10

0

10

Admittance (mGal/km)

Figure 4.9: Magnitude of implied erosion as a function of admittance given several
combinations of mean shell thickness (d) and elastic thickness (T ). Negative top
loads correspond to surface erosion. The dashed black line indicates the mean admittance estimate obtained from the Monte Carlo analysis based on the SOL1a gravity(Iess et al., 2012) and the Deg6-exp topography (Zebker et al., 2012) ( 32 mGal/km).

4.4

Discussion

4.4.1

Ice Shell Structure and Rigidity

Throughout our analysis, we assume the ice shell to be in an equilibrium state
where the various forces (flexure within the elastic part of the shell, weight ofthe overlying topography, and buoyancy of the root) are in balance. This isreasonable because the vertical response time of the shell should be fast compared with the loading timescale. The survival of large-impact basin rims andother topographic features (Neish et al., 2013; Moore et al., 2013, 2014) rulesout erosion amplitudes greater than ⇠ 1 km, implying a shell elastic thickness
greater than 40 km (Figure 4.7). This high rigidity could be the result of acold (ammonia-rich) subsurface ocean (Moore and Pappalardo, 2011; Grasset112

et al., 2000) (section 4.4.1.1), a clathrate-rich shell (Tobie et al., 2006) or low
heat flux from the interior (section 4.4.1.2). The addition of a rigid shell wouldslightly reduce the k2 tidal Love number, but not enough to conflict with themeasured value (Iess et al., 2012) (section 4.4.1.5). The implied elastic thickness rules out a vigorously convecting shell, limits the potential for widespread
cryovolcanism (Moore and Pappalardo, 2011; Lopes et al., 2007), and permits
the survival of lateral shell thickness variations (section 4.4.1.4). To generate
the observed topography, our model requires shell thickness variations at degree 3; possible sources include tidal heating in a laterally heterogeneous shell(Běhounková et al., 2012), and redistribution of material via non-Newtonian
lateral flow (Nimmo, 2004) (section 4.4.1.3). A rigid conductive shell is also
only weakly dissipative, potentially helping to explain Titan’s high present-day
orbital eccentricity (Tobie et al., 2006; Nimmo and Bills, 2010).
4.4.1.1

Temperature Profile

Roughly speaking, ice will undergo a transition from elastic to viscous behavior
at temperatures in the range 160

180 K, depending on the exact strain rate

and grain size assumed (Nimmo et al., 2002). For a conductive ice shell withbasal temperature Tb = 270 K, the elastic thickness (T ) will then be 38-50%of the total shell thickness (d), while if Tb = 210 K, then T will be 58-75% of
d. While the transition from elastic to viscous behavior will occur over somefinite region, that region will be thin because of the very strong variation inviscosity with depth. Hence, a two-layer model is a good approximation.
In order to obtain conservative estimates for the elastic thickness and surface erosion, we assumed T = d in parts of our analysis (e.g., Figure 4.7a).
Such a scenario would require that Tb ⇡ 180 K, which may not be likely, but
113is at least possible in the limit of an ammonia-water eutectic ocean, for which
the freezing point is ⇠ 176 K.
4.4.1.2

Heat Flux

Assuming a heat flux of F ⇡ 4 mW/m2 through Titan’s (conductive) ice shell(Nimmo and Bills, 2010), and allowing thermal conductivity to vary withtemperature (Klinger , 1980), we can estimate the shell’s elastic thickness (T )according to
567
ln
T ⇡
F

✓

Tz
Ts

◆

(4.31)

where Ts is the surface temperature and Tz is the temperature at which theshell transitions from elastic to viscous behavior. The resulting estimatedelastic thickness is T ⇡ 82

98 km. Despite the highly approximate nature ofthis analysis, it yields an elastic thickness that is consistent with our estimates(Figure 4.7b). Larger elastic thicknesses are possible, but would correspondto a lower heat flux.
4.4.1.3

Origin of Long-Wavelength Topography

Like any tidally locked satellite, Titan’s shape is primarily determined bytidal and rotational distortion. Centrifugal forces from the rotation lead toan equatorial bulge (flattening at the poles) and the permanent tide leads toelongation along the Titan-Saturn axis. The degree-2 and degree-4 topography
will also be aﬀected by variations in ice shell thickness that arise due to tidal
heating (Nimmo and Bills, 2010). As discussed below (section 4.4.1.4), nonNewtonian flow within the lower part of the ice shell could cause degree-3
shell thickness variations to develop from a pattern that is initially confined to114

degrees 2 and 4. This may, in part, explain the source of the observed degree-3
topography. Since ongoing tidal heating will support the maintenance of shell
thickness variations at degrees 2 and 4, those variations could persist even as
lower crustal flow continues to generate shell thickness variations at degree 3.
Heterogeneities in the ice shell could also contribute to a departure from thepurely degree-2 and -4 pattern predicted from tidal heating and could thus beresponsible for part of the shell thickness variations, and therefore topography,
at degree 3.
4.4.1.4

Lateral Flow in the Ice Shell

Shell thickness variations lead to flow in the lowermost, low-viscosity part ofthe shell, which will tend to smooth out any such variations. For a Newtonian
fluid, the timescale (⌧ ) for removal of variations is given by (Nimmo, 2004)⌧=⌘bg ⇢ 3k2(4.32)

where ⌘b is the viscosity at the base of the shell, g is the acceleration due togravity,⇢ is the density contrast between the shell and the fluid underneath,is the eﬀective channel thickness in which flow occurs, and k is the wavenumber(k = l/R, where l is spherical harmonic degree and R is the planetary radius).
Assuming a linear temperature gradient and a thermal conductivity (Klinger ,1980) which goes as c/T , where c = 567 W/m, the eﬀective channel thickness( ) is given by (Nimmo, 2004)=

Rg Tb dQ ln (Tb /Ts )(4.33)

where Rg is the gas constant, d is the shell thickness, Q is the activation energy
115and Tb and Ts are the basal and surface temperatures, respectively. Finally,
the viscosity ⌘b is given byQ⌘b = ⌘ref exp
Rg

✓

1
Tb

1
Tref

◆

(4.34)

where the viscosity of ice is ⌘ref at a reference temperature Tref .
Table 4.3 gives ⌧ as a function of Tb for spherical harmonic degree 3. Here,
we have assumed ⌘ref = 1014 Pa s at Tref = 273 K, Q = 60 kJ/mol, Ts = 90 K,
d = 100 km, g = 1.35 m/s2 ,⇢ = 80 kg/m3 and R = 2575 km. For therange of Tb values explored,= 3.4 km. Rheological parameters (Goldsby and

Kohlstedt, 2001) are subject to some uncertainty; nonetheless, the results ofTable 4.3 serve to illustrate the main conclusion, which is that flow is slow if
the ocean is suﬃciently cold (Tb . 220 K). A temperature of 220 K corresponds
H2 O system (Kargel , 1998).

to 25 wt% ammonia in a simple NH3

Table 4.3: Viscosity and timescale for removal of degree-3 ice shell thickness variations as a function of temperature.

Tb (K) ⌘b (Pa s)
273
1.0 ⇥ 1014250
1.1 ⇥ 1015
230
1.4 ⇥ 1016
210
2.8 ⇥ 1017

⌧ (Myr)
0.55
6.34
77.6
1491

In practice, the rheology of ice may be non-Newtonian, in which case ourflow timescales will be underestimates (Nimmo, 2004), permitting larger valuesof Tb . An important consequence of non-Newtonian flow is that mode-coupling
occurs: an initially degree-2 or degree-4 pattern (e.g., due to tidal heating) willdevelop a degree-3 component as flow proceeds, thus potentially explaining theobserved degree-3 signal (see 4.4.1.3).

116

Finally, we note that if tidal heating is indeed occurring, a balance may
develop wherein shell thickness variations are being generated by tidal heating
just as quickly as lateral flow is removing those variations. Such an equilibrium
situation could be stable even if the relaxation timescales are relatively short,again permitting higher values of Tb .
4.4.1.5

Eﬀect on Tidal k2

The tidal Love number, k2 , of Titan has been measured (Iess et al., 2012),
with a 2- lower bound on k2 of 0.413-0.439. It is therefore important to checkthat the kind of rigid elastic lid that we are proposing does not contradict
the observations. To do so, we constructed a highly simplified model for theinterior of Titan (Table 4.4).
Table 4.4: Simple model of Titan’s interior, used to determine the eﬀect of a rigid
shell on the tidal Love number, k2 .

Layer
Solid Core
High-Pressure Ice
Ocean
Rigid Outer Shell

Outer
Radius
(km)
2110
2275
2575 d2575

Rigidity(GPa)

Viscosity(Pa s)

Density(kg/m3 )

3
3
0
3

1021
1021
0
1021

2600
1000
1000
920

This model is not meant to be realistic, but suﬃces to demonstrate ourresults. The density structure approximately satisfies the nominal momentof inertia constraint, while the low rigidity in the inner layers is designed
to reproduce the observed k2 . We followed Moore and Schubert (2003) incalculating the model k2 values. For rigid shell thicknesses of d =5, 100 and
200 km we obtained k2 values of 0.568, 0.519 and 0.413, respectively. Hence,adding a rigid shell of thickness 100 km or 200 km reduces k2 by 9% or 27%,117

respectively—not enough to contradict the observed k2 values (Iess et al.,2012).

4.4.2

Erosion and Deposition

Our lower bound on the extent of erosion and deposition (200 m) over global
length scales is compatible with estimated local amounts based on impact
crater degradation (Neish et al., 2013; Moore et al., 2013, 2014), but somewhat
larger than estimates from fluvial incision (Black et al., 2012). The impliedvertical erosion/deposition rate is of the order of a meter per million years,
comparable to the lower end of aeolian deposition rates measured on Earth(Patterson et al., 1999); erosion may be occurring by physical comminution,
dissolution or sublimation. Because Titan’s topography is high at the equator,
we predict maximum erosion occurring at the equator and sediment transport,
via fluvial or aeolian processes, predominantly towards the poles. An alternative possibility is transport in the vapor phase, if the mobile material has asublimation temperature close to that of the surface temperature.

4.4.3

Degree-2 Admittance and Fluid Love Number

Degree-2 admittance analysis of Titan’s ice shell is complicated by the factthat the body is tidally and rotationally distorted. Tidal/rotational distortion
dominates the degree-2 gravity signal and also makes a large contribution to
the degree-2 topography. If we assume that basal freezing, uplift and erosion
processes act similarly at degrees 2 and 3, then ice shell thickness variations
should be responsible for a portion of the observed degree-2 gravity and topography signals.
118

Separating the degree-2 gravity and topography signals into their hydrostatic (i.e., tidal/rotational) and non-hydrostatic (i.e., due to ice shell thickness
variations) parts is necessarily an iterative process. We begin by estimating
Titan’s fluid Love number, h2f , from the observed degree-2 gravity signal (J2 )according to
h2f = 1 +6 g
J2
5 R! 2(4.35)

where R is Titan’s mean radius, g is the mean surface gravity, and ! is theangular frequency of rotation. Based on the observed gravity field (Iess et al.,2012), we obtain h2f ⇡ 2.0 (moment of inertia factor ⇠ 0.34). This allows us
to predict the expected hydrostatic topography, hT according toR2 ! 2 1
3 cos2hT = h2f
g
2

+1

1

cos2 ✓

5
6(4.36)

which we then subtract from the observed topography to get the non-hydrostatic
topography, hshell (i.e., that which is due to variations in ice shell thickness).
We then multiply hshell by equation (4.24) to estimate the gravity signal due toan ice shell with anomalously deep roots (the large root size will be forced implicitly by our choices of C and dt , both of which will be a function of mean shell
thickness, d and elastic thickness, T ). Having obtained this gravity anomaly
due to ice shell thickness variations ( gshell ), we conclude that the portion of
the gravity signal that is due to tidal distortion is

gtidal =gtotalgshell . Fi-

nally, we use the newly obtained tidal gravity field to get an updated estimate
for h2f , again using (4.35). After 3-4 iterations, our estimate for h2f converges to the fourth decimal place, allowing us to separate, in a self-consistentway, the tidal/rotational and ice shell thickness contributions to the degree-2

119

gravity signal. The final estimate for h2f depends on the assumed mean shell
thickness (d) and elastic thickness (T ), as illustrated in Figure 4.10. If weassume that T = 100 km and d = 200 km, then we obtain h2f t 2.15 (momentof inertia factor ⇠ 0.36).
2.2

11

2

3

2.1

2.

2.1

2.16
2.14

400

2.12
14

2.1

2.

300

15

2.

17

2.

6

18

2.

9
2.1

2.08

1
2.

200

2.06
654321987654321987654321987654321987654321987654321987
2212111.1.1.100.11098.079.68057.40635.0241.3020.918.7968.574.6352413654321908721654321987654321987654321987654321987

Fluid Love number, h2f

09

2.
1

2.18
2.

2.0
7

500

Shell Thickness, d (km)

2.
08

600

2.04

100

2.02

50

100

150

200

250

300

350

400

2

Elastic Thickness, T (km)Figure 4.10: Estimate for Titan’s fluid Love number, h2f , as a function of mean shell
thickness (d) and elastic thickness (T ).

Using this fluid Love number, we can estimate the non-hydrostatic portionsof the degree-2 topography and gravity signals. Figure 4.11 shows how theobserved degree-2 gravity field (a) compares with the predicted field (b), where
the predicted field is the sum of the estimated hydrostatic gravity (c), basedon h2f = 2.15, and the gravity anomaly expected from the estimated ice shell
thickness variations (d), assuming an erosion amplitude of 577 m (obtained
from Figure 4.7b assuming T = 100 km and d = 200 km). The amplitude of theestimated non-hydrostatic gravity is ⇠ 2 mGal while the estimated hydrostatic
gravity amplitude is ⇠ 21 mGal. This is one measure of Titan’s departure fromhydrostatic equilibrium.
120

(a)

(b)

(c)

(d)
Figure 4.11: Degree-2 gravity maps centered on the sub-Saturnian point (180° longitude): (a) SOL1a gravity (Iess et al., 2012); (b) total predicted gravity signal; (c)gravity signal caused by tidal/rotational distortion assuming h2f = 2.15; (d) gravity
signal caused by ice shell thickness variations assuming T = 100 km and d = 200 km,
and therefore 577 m of erosion.

1214.4.4Degree-4 Predictions

Assuming once again that the observed topography is the result of uplift due
to basal freezing and that the surface has experienced a similar magnitude oferosion at degrees 3 and 4, we can predict the admittance at degree 4. We
first obtain the magnitude of degree-3 erosion (dt ) over a range of values for T
and d from Figure 4.7b. Using this same value of dt for the degree-4 erosion
amplitude, we then use equation (4.24) to estimate the admittance over thesame range of values for T and d. Figure 4.12 illustrates that the degree-4
admittance should be negative if the elastic thickness, T , accounts for most of
the total shell thickness, d.
If, for example, T = d = 200 km, and we assume the same magnitude oferosion at degrees 3 and 4 (in this case, 241 m), then based on the observed
topography (Zebker et al., 2012) (Figure 4.13a), we obtain a degree-4 admittance of

5.1 mGal/km and we can compute the implied degree-4 gravity

anomaly everywhere over the surface (Figure 4.13b). Although the amplitudes are similar, our result is spatially unlike the reported degree-4 gravity
field (Iess et al., 2012) (Figure 4.13c). However, as noted in section 4.2.1, thedegree-4 gravity field is not currently regarded as reliable; future gravity flybys
are expected to improve the determination of the degree-4 field by a factor of
two, providing a better test of our prediction. Note that if we assume insteadthat T = 100 km and d = 200 km (the corresponding erosion amplitude being 577 m), the greater compensation leads to a muted gravity signal and an
admittance that approaches zero. Hence a weak degree-4 gravity signal, or adegree-4 admittance that is positive or only weakly negative, may be an indication that the elastic layer accounts for a smaller portion of the total shell122thickness (Figure 4.12).
20

600
16

14

15
12
10

2

8

400

10

0

6
4

−6

−2

5

−4

300

−−420

0

200

Degree−4 Admittance (mGal/km)

Shell Thickness, d (km)

500

−5

100
−6

50

100

150

200

250

300

350

400

−10

Elastic Thickness, T (km)Figure 4.12: Degree-4 admittance predicted for a range of elastic thicknesses (T )and total shell thicknesses (d).

4.5

Conclusions

Titan’s degree-3 gravity and topography exhibit a strong inverse correlation
and, consequently, the degree-3 admittance is substantially negative, evenwhen model uncertainties are taken into account. Although negative admittances are not common, our model predicts negative admittances at degree 3
under a particular set of conditions: when the ice shell is substantially rigid(with an elastic thickness greater than ⇠ 40 km), when it has been loaded
primarily from below (e.g., due to basal freezing), and when some few hundred meters of erosion has taken place at topographic highs. The conclusionthat the ice shell is rigid rules out vigorous convection and limits the potential
for widespread cryovolcanism. The substantial erosion required by our model
suggests a tendency for sediment to be transported from the high-standing
123

(a)

(b)

(c)
Figure 4.13: Degree-4 topography and gravity maps centered on the sub-Saturnian
point (180° longitude): (a) Topography (Zebker et al., 2012); (b) Predicted gravity
signal caused by ice shell thickness variations assuming T = d = 200 km and therefore
241 m of erosion; (c) SOL1b gravity (Iess et al., 2012) (not currently considered
reliable at degree 4).

124topography near the equator to the lower elevations found toward the poles.
If similar erosion and deposition processes are taking place at degree 2,then the present-day degree-2 topography is smaller than what should be expected from tidal and rotational distortion alone. This would imply that theobserved degree-2 gravity signal is an underestimate of the hydrostatic gravity
and that, consequently, Titan’s fluid Love number has been underestimated.
This would indicate that Titan has a higher moment of inertia (and is thus
even less centrally condensed) than previously thought, reinforcing the needto understand how so large a body could have accreted without undergoing
more complete diﬀerentiation. Our analysis suggests that, if the elastic thickness exceeds half the total shell thickness, the degree-4 admittance will also be
negative. Future Cassini gravity fly-bys will improve the determination of thedegree-4 gravity field by a factor of two and will therefore provide a partial
test of this prediction. If a non-negative degree-4 admittance is observed, that
may be an indication that a smaller portion (less than half) of the ice shell iselastic.

Acknowledgements
We thank the Cassini radar science team, M. Manga, D. Stevenson, R. Pappalardo, W. McKinnon, and two anonymous reviewers for their suggestions.
Portions of this work were supported by NASA grants NNX13AG02G and
NNX11AK44G.

125

Chapter 5Internal Structure of Enceladus
This chapter details analysis I carried out, along with Francis Nimmo, in
support of Iess, L. et al. (2014), The Gravity Field and Interior Structure of
Enceladus, Science, 344 (6179), 78–80, doi:10.1126/science.1250551.

Abstract
With water-rich jets actively erupting from its South Polar Terrain (Spitaleand Porco, 2007), Saturn’s tiny moon Enceladus is one of the most compelling
exploration targets in the solar system. A subsurface liquid sea or ocean is often presumed to be the source of the jets. However, without better knowledgeof Enceladus’ gravity field, it has been diﬃcult to constrain interior structure models suﬃciently to confirm the presence of an internal ocean. Here we
combine the newly determined gravity field (Iess et al., 2014), obtained from
Doppler tracking of the Cassini spacecraft, with a shape model determinedfrom limb profile analyses (Nimmo et al., 2011), to constrain the interior structure of Enceladus, confirming that a subsurface ocean is likely. We employ an126admittance analysis that self-consistently separates the hydrostatic and nonhydrostatic components of the gravity and topography signals in a way thatallows us to determine both the moment of inertia for Enceladus and the depth
at which the topography may be compensated isostatically. The observations
are most easily accommodated by a diﬀerentiated Enceladus, with moment ofinertia ⇠ 0.335M R2 , surrounded by an H2 O mantle consisting of an outer ice
shell, 20 40 km thick, and a subsurface liquid ocean with an average thickness> 15 km.

5.1

Introduction

Enceladus, embedded in Saturn’s E-ring, is a small (⇠ 250 km radius), but
distinctly geologically active world (Porco et al., 2006). In particular, its South
Polar Terrain (SPT) is cross-cut by four major fissures, known as the “Tiger
Stripes”, along which more than 100 distinct geysers have been identified.
The geysers spray water, rich in salts and organic compounds, out from thesmall moon, resulting in a bright, “snow” covered surface, and a persistentvapor plume that extends outward some 500 km (Porco et al., 2006; Spitaleand Porco, 2007). Shear heating (Nimmo et al., 2007) has been proposed asa mechanism for the vapor production and a subsurface liquid sea or ocean,although not directly observed, is considered to be likely (Spencer and Nimmo,2013).
Reliable internal structure models have not previously been available tohelp confirm or refute the existence of a subsurface ocean. Although good
shape models have been available for some time (Nimmo et al., 2011), determination of the gravity field, which requires dedicated flybys and modeling127

of unknown non-gravitational disturbances (e.g., from plume drag), has only
recently been completed (Iess et al., 2014). Below, we combine the shape models with the newly determined gravity field in order to construct an interior
structure model for Enceladus.
In Chapter 4, we used the spectral admittance, Z (l), to quantify the relationship between Titan’s gravity and topography, and used the result to draw
conclusions about the thickness and elastic behavior of the ice shell. Our analysis was focused mainly on the degree-3 admittance since analysis of degree-2
admittance is complicated by the eﬀects of tidal and rotational distortion. In
the case of Enceladus, the gravity data are even more limited: the degree2 gravity field is well constrained, but determination of the degree-3 field islimited to the zonal term (C30 ) alone. Here, we instead take on the task of decomposing the degree-2 gravity and topography signals into their hydrostaticand non-hydrostatic components. The hydrostatic components of the degree-2
signal are controlled by tidal and rotational distortion. The remainder of thedegree-2 gravity and topography signals are due to additional mass and shapeanomalies that may arise, for example, from an irregular core shape (Thomas
et al., 2007; McKinnon, 2013), or from lateral variations in the thickness ofthe outer ice shell (Nimmo and Bills, 2010).
In section 5.2, we briefly describe the observed gravity and topography
and the degree to which the degree-2 terms deviate from the expectation fora relaxed, hydrostatic body. In section 5.3, we review the relevant theoretical
tools for modeling the hydrostatic gravity and topography, for computing admittances under various conditions of isostasy and flexure, and for determining
the moment of inertia. In section 5.4, we discuss the results before concluding
in section 5.5.
128

5.2

Observations

5.2.1

Gravity Field

The gravitational field of Enceladus is described by Iess et al. (2014) as nongg
, Slmsuch that the gravinormalized, dimensionless potential coeﬃcients Clm

tational potential at the reference radius (R) can be written

U (✓, ) =1l

GM X X g
g(Clm cos m + Slm
sin m ) Plm (cos ✓)R l=0 m=0

where Plm (cos ✓) are the non-normalized associated Legendre functions (seeAppendix A), ✓ is colatitude,gg
and Slm
are non-normalized
is longitude and Clm

spherical harmonic coeﬃcients of degree l and order m.
Enceladus’ low-order gravity field was determined (Iess et al., 2014) via radio tracking of the Cassini spacecraft during three separate flybys in which thespacecraft passed within 100 km of the moon, twice over the southern hemisphere (E9, E19) and once over the northern hemisphere (E12). Doppler datafrom the three flybys were analyzed separately and combined into a multi-arc
solution (Figure 5.1). Constraining the gravity field was especially challenging in the case of Enceladus due to its small size (⇠ 250 km radius) but also
due to non-gravitational eﬀects including the spacecraft’s interaction with thevapor plume, which causes appreciable drag in the ⇠ 20 seconds around closest approach during passes over the southern hemisphere. Accounting for thisgeﬀect, Iess et al. (2014) were able to constrain the major degree-2 terms (C20gg
and C22) to 1% uncertainty along with the degree-3 zonal term (C30
) to 20%

uncertainty.

129

90

20
30

10

0

0
−10

−30

Geoid (meters)

Latitude (degrees N)

30
60

−20
−60
−30
−90
0

45

90

135

180

225

270

315

360

Longitude (degrees E)

Figure 5.1: Enceladus’ geoid with respect to the reference ellipsoid (see Appendix C).
Also shown are the Cassini ground tracks from the three gravity flybys: E9 (orange),
E12 (blue), and E19 (green).

5.2.2

Topography

Titan’s topography is described by Nimmo et al. (2011) in terms of fully norhh
, S̄lm, such that the topography at
malized spherical harmonic coeﬃcients C̄lm

any point on the surface is given byH(✓, ) =

l
1 X
X

h
h
cos m + S̄lm
sin m
C̄lm

P̄lm (cos ✓)l=0 m=0

where P̄lm (cos ✓) are the fully normalized associated Legendre functions (seeAppendix A). To facilitate direct comparison with the non-normalized gravhh, Slmfrom theity coeﬃcients, we obtain the non-normalized coeﬃcients Clm

normalized coeﬃcients with

{C, S}hlm =s(20m ) (2l + 1)

(l m)!
C̄, S̄
(l + m)!

hlm

so that the topography is described byH(✓, ) =

l
1 X
X

h
h
cos m + Slm
sin m
Clm

l=0 m=0130

Plm (cos ✓)

(5.1)

1000

60
500

30
0

0

−30

−500

Topography (m)Latitude (degrees N)

90

−60
−90
0

−1000
45

90

135

180

225

270

315

360

Longitude (degrees E)

Figure 5.2: Enceladus’ topography for harmonic degrees 3 through 8.

where, again, Plm (cos ✓) are the non-normalized associated Legendre functions.
Enceladus’ topography was determined through analysis of limb-profiles(Nimmo et al., 2011). 43 separate limb profiles were measured from images ofEnceladus, providing 43 arcs of topographic data, which were then combined toconstrain a global shape model up to spherical harmonic degree 8 (Figure 5.2).
hh, about 2% for C22, and
The shape model’s uncertainty is less than 5% for C20h.
about 1% for C30

5.2.3

Deviation from Hydrostatic Equilibrium

For tidally-locked (synchronous) satellites, the degree-2 shape and gravity field
are dominated by rotational flattening and elongation due to the permanenttide (see Appendix C). If the body is in hydrostatic equilibrium, its shape andgravity field conform to the expectation for a body whose outer layer is aninviscid fluid, in which case, only the C20 and C22 terms are significant, andhave the approximate ratio
C20
10⇡ 3.33
=
C22
3

131

While this approximation is excellent in most cases, it is an overestimate
in the case of very fast-rotating bodies, such as Enceladus (Tricarico, 2014).
For Enceladus, a ratio of ⇠ 3.25 may be more appropriate (Mckinnon, 2015).
We will proceed here assuming the typical hydrostatic ratio of 10/3 but, for ananalysis that builds on ours while taking this eﬀect into account, see Mckinnon(2015).
For Enceladus, the gravity coeﬃcients have the ratiogC20⇡ 3.51±0.06
g
C22

and the topography coeﬃcients have the ratiohC20⇡ 4.20±0.28hC22

Whereas the topography deviates substantially from the hydrostatic expectation, the gravity field deviates only modestly. This is an indication thatthe excess non-hydrostatic topography may be largely compensated by internal density anomalies, as would be expected for topography that is supported
isostatically. Perhaps the most conspicuous feature of the gravity field is thenegative gravity anomaly over the South Polar Terrain (Figure 5.1). Such agravity anomaly is expected given the large topographic depression in that
region (Figure 5.2), but the magnitude of the gravity anomaly in relation to
the magnitude of the topographic depression will be diagnostic of the degree
to which the topography is compensated isostatically.

132

5.3

Theory5.3.1

Separation of Hydrostatic Component

We begin by presuming that the gravity and topography coeﬃcients can be
decomposed into their hydrostatic and non-hydrostatic parts as follows
g
g,hyd
g,nh
= Clm
+ Clm
Clm

(5.2)

h,hyd
h,nh
h
= Clm
+ Clm
Clm

(5.3)

g,hyd
h,hyd
where Clmand Clm
are the hydrostatic components: the parts of thesignals that result from tidal and rotational distortion. The hydrostatic partsof the gravity and topography signals are a function of the body’s moment ofinertia. If we make an assumption about the moment of inertia factor (↵), we
can obtain the fluid Love number, h2f , via the Radau-Darwin approximation(Murray and Dermott, 1999)
2
I
=
↵=
2
MR
3

1

2
5

s

5
h2f

1

!

(5.4)

The fluid Love number describes the degree to which a body deforms in
response to the tidal and rotational deforming potentials. Given the fluidLove number and using the present-day rotation rate, we can estimate thehydrostatic shape and gravity field (Appendix C). In terms of non-normalized
spherical harmonic coeﬃcients, these are
h,hyd
=
C20

5 R2 ! 2
h2f
6 g

133

(5.5)

h,hydC22=

1 R2 ! 2
h2f
4 g

(5.6)

The change in shape, and therefore mass distribution, also produces acorresponding change in the gravitational potential such that
g,hyd
=
C20

5 R! 2
(h2f
6 g

1)

(5.7)

g,hyd
=
C22

1 R! 2
(h2f
4 g

1)

(5.8)

Of course, we don’t yet know the moment of inertia (↵), so any such calculation is only a guess. But if we make the key assumption that the nonhydrostatic gravity-to-topography ratios should be equal (i.e., that the admittance is isotropic), then the correct moment of inertia is the one that yields

Z20 = Z22
where the component-wise admittance is defined asg,nhClm
Zlm = h,nh
Clm

(5.9)

More specifically,g+ 56 R!g (h2f
C20
2

Z20 =

Z22 =

h
C20
+6

5 R2 ! 2
g

1)

h2f

g
C22

1 R! 2(h2f 1)4 g
1 R2 ! 2
h
C22
h2f
4 g

The task now becomes finding a value of h2f for which Z20 = Z22 .

134(5.10)(5.11)

5.3.2Model Admittance

We already derived the model admittance in Chapter 4 (section 4.3.1). However, since we are now working with dimensionless gravitational potential coeﬃcients rather than acceleration, admittance here has units of m 1 , and wemust alter the equations slightly. Specifically, the admittances we use here
will be those of the previous chapter (and Hemingway et al. (2013)) divided= G 43 ⇡R¯⇢ (l + 1).
by GMR(l+1)
2Much of our analysis here will make use of the simple case in which we
assume that surface loads are compensated completely, for which the degree-l
admittance is simply (c.f. equation (4.26) and section 5.3.3)3⇢c1
Z (l) =
(2l + 1) R¯
⇢

✓

1

d
R

◆l !

(5.12)

where R and ⇢¯ are, respectively, the radius and bulk density of Enceladus, ⇢cis the density of the surface load, which we take to be the same as the densityof the icy crust, and d is the compensation depth.
In the special case where we assume a surface load that is entirely uncompensated (i.e., supported rigidly, see section 5.3.3), thenZ (l) =

5.3.3

3⇢c
(2l + 1) R¯⇢(5.13)

Elastic Flexure

More generally, we may need to account for the eﬀects of elastic flexure (bending and membrane stresses) in the lithosphere. If the lithosphere is at least
partially elastic and is deflected downward by the addition of a surface load,135

then the degree-l admittance is given by (c.f. equation (4.28))
3⇢c
Z (l) =
1
(2l + 1) R¯
⇢

✓

Ct 1

d
R

◆l !

(5.14)

where Ct is the degree of compensation assuming top-loading (where Ct = 1
means complete isostatic compensation and Ct = 0 means zero compensation, or perfectly rigid support). Similarly, if loading occurs purely from thebottom (i.e., under-plating) then the degree-l admittance is given by (c.f.
equation (4.29))
1
Cb

3⇢c
1
Z (l) =
(2l + 1) R¯
⇢

✓

1

d
R

◆l !

(5.15)

where Cb is the degree of compensation assuming pure bottom-loading. The
compensation terms Ct and Cb become important when considering the eﬀects
of elastic flexure (see Chapter 4, section 4.3.1 for the derivations) and aredefined as
Cb =1
1 + ⇢µc(5.16)

and
Ct =
where

⇢ = ⇢m1
1 + µ⇢ 1

dR(5.17)

2⇢c is the density contrast at the water/ice interface (⇢m isthe density of the subsurface ocean “mantle”), where µ describes the flexural
rigidity at a particular wavelength and gravity, and is defined in equation (4.21)
asµ(l) =✓

13⇢m(2l + 1)⇢◆ 1

0 T 2 l3 (l+1)3 4l2 (l+1)2[]+ l(l + 1)
ET @R2 12(1 ⌫ 2 )R2 gl (l + 1) (1 ⌫)
136

2

1

A (5.18)

5.3.4

Moment of Inertia

For a simple, spherically symmetric, two-layer body, the moment of inertia isgiven byI=8⇡ ⇢o R5 + (⇢i15d)5⇢o ) (R(5.19)

where the subscripts (i) and (o) refer to the inner and outer layers, respectively,
and where d is the thickness of the mantle so that the core radius is Rcore =Rd. We will often discuss the normalized moment of inertia (also called themoment of inertia factor), ↵ = MIR2 , which can be written2↵=⇢o + (⇢i5¯⇢✓⇢o ) 1dR◆5 !

(5.20)

where ⇢¯ is the bulk density of Enceladus.

5.4

Results and Discussion

5.4.1

Admittance and Moment of Inertia

Considering a range of possible moments of inertia, we use equations (5.4),(5.10), and (5.11) to compute the corresponding Z20 and Z22 admittances andfind that the two independent estimates converge at ⇠ 3.4 ⇥ 10 7 m 1 whenthe moment of inertia factor is ⇠ 0.335 (Figure 5.3).
g), we
Because we also have the zonal part of the degree-3 gravity field (C30

can compute Z30 and use it as an independent check on our admittance estimate (note that the degree-3 terms have no hydrostatic component so Z30 is obtained directly from equation (5.9) and does not depend on moment of inertia).
Because they reflect behavior at diﬀerent wavelengths, however, the degree-3
137

−7

18

x 10

16
14

Uncompensated (top−loading)

Admittance (/m)

12
10
8
6
4

Compensated (d=40km)

2

Compensated (d=20km)

0
−2
0.3

0.31

0.32

0.33
0.34
0.35
Moment of Inertia Factor

0.36

0.37

Figure 5.3: Region of 1 agreement (dark blue) and 2 agreement (pale blue) between the admittances estimated separately from the (2,0) (blue line), (2,2) (redline), and (3,0) (green line) terms. The (3,0) term has been multiplied by 0.99
to obtain the equivalent degree-2 value (see text, section 5.4.1). The solid line at⇠ 14 ⇥ 10 7 m 1 represents the admittance expected for uncompensated (rigidly
supported) topography. The dashed horizontal lines near 2 and 4 ⇥ 10 7 m 1 represent the admittances expected for fully compensated topography with compensation
depths (ice shell thicknesses) of 20 and 40 km, respectively.

138

Table 5.1: Parameter values assumed for admittance calculations

Parameter
SymbolPoisson’s ratio for ice⌫
Young’s modulus for ice
E
Crustal (ice shell) density⇢cMantle (subsurface ocean) density⇢mEnceladus’ mean density⇢
Enceladus’ radiusR
Enceladus’ surface gravitygAssumed Value
0.25
9 GPa
920 kg/m3
1000 kg/m3
1609 kg/m3
252.1 km0.114 m/s2

admittance (Z30 ) cannot be compared directly with the degree-2 admittances(Z20 and Z22 ). The relevant comparison is between the degree-2 estimates (Z20
and Z22 ) and the degree-3 estimate, scaled by the ratio Z (2) /Z (3), which we
obtain by evaluating equation (5.12) at both l = 2 and l = 3. Doing this calculation requires assuming a compensation model. For simplicity, we assume
complete isostatic compensation (C = 1) and then we solve equation (5.12) for
d. Given that (ignoring uncertainties for the moment) Z30 = 3.003⇥10 7 m 1 ,and assuming the parameters given in Table 5.1, we obtain d ⇡ 29 km. The
⇣
⌘
Z(2)
⇡
0.99
so
that
Z
⇡ 2.973 ⇥ 10 7 m 1 . This scaled
result is that Z(2)
30
Z(3)
Z(3)
degree-3 admittance (green line in Figure 5.3) plots very close to the intersec-

tion of the Z20 and Z22 lines (blue and red lines, respectively). This result thus
provides an a posteriori check that the assumption of isotropic admittance isappropriate. Table 5.2 lists the gravity and topography coeﬃcients separated
into their hydrostatic and non-hydrostatic parts, along with the resulting admittances.
Taking uncertainties into account, agreement between the three estimatesoccurs at the 1

level for normalized moments of inertia (↵ = MIR2 ) in therange 0.333–0.338 and admittances in the range 2.4

3.6 ⇥ 10 7 m 1 (Fig-ure 5.3). Agreement at the 2 level occurs in the range ↵ = 0.330
139

0.340 and

Table 5.2: Gravity and topography coeﬃcients separated into their hydrostatic andnon-hydrostatic components, along with the resulting admittance estimate. A moment of inertia factor of ↵ = 0.335 is assumed in computing the hydrostatic gravity
and topography.
CoeﬃcientObservedg
C20

5.435 ⇥ 10

h
C20

3

3, 846 m

g
C22

1.550 ⇥ 10h
C220.115 ⇥ 10hC30

5.000 ⇥ 10

Non-Hydrostatic

3

2, 569 m
3

917 m

g
C30

Hydrostatic

1.500 ⇥ 10

3

771 m

0.050 ⇥ 10
0.115 ⇥ 10384 m

Z20 = 3.411 ⇥ 10

7m 1

Z22 = 3.421 ⇥ 107m 1

Z30 = 3.003 ⇥ 10

7m 1

3

146 m

3

Admittance3

1, 277 m

384 m

Z (2) = 1.8

0.436 ⇥ 103

4.2 ⇥ 10 7 m 1 . The regions of agreement (blue shaded areas in

Figure 5.3) were obtained by performing a Monte Carlo analysis in which eachadmittance value is calculated 100,000 times with the gravity and topography
coeﬃcients being distributed randomly according to their formal uncertainties.
The regions in which there is overlap between the resulting 1 ranges of thethree separate admittance estimates (based on the Z20 , Z22 , and Z30 estimates)
is considered to be the region of agreement at the 1 level (dark blue shaded
area in Figure 5.3). The same procedure was carried out to determine theregion of agreement at the 2 level (pale blue shaded area in Figure 5.3).

5.4.2

Isostatic Compensation

Admittance is a function of the degree and depth of compensation (e.g., equation (5.14)), with large admittances indicating little compensation or verylarge compensation depths, and small admittances indicating substantial (oreven complete) compensation and/or shallow compensation depths. The observed degree-2 admittance, Z (2) ⇡ 3.4 ⇥ 10 7 m 1 , is much smaller than
the admittance given by equation (5.13) for fully uncompensated (that is,
rigidly supported) topography, ⇠ 14 ⇥ 10 7 m 1 , suggesting that the topog140

raphy is highly compensated (Figure 5.3). If the topography is compensatedcompletely (C = 1), then we can use equation (5.12) to estimate the depth ofcompensation. The observed admittance compares well with the expectation
for complete compensation at depths of roughly 20

40 km (dashed lines in

Figure 5.3).

5.4.3

Elastic Support

If there is some measure of elastic strength in the lithosphere, then our compensation depth estimates must be altered. If we assume the ice shell has been
loaded from the top, we can use equation (5.14) to see how admittance varies
with compensation depth (d) and the elastic layer thickness (T ). Enceladus isso small that, even a very thin elastic layer provides significant flexural support
and leads to large admittances (Figure 5.4a). To produce the admittance weobserve, for the case of top loading, the elastic layer must be less than ⇠ 25 min thickness. If the topography is supported in part by elastic flexure, thecompensation depth must be even shallower than the 20

40 km we estimatedfor the fully isostatic case.
Conversely, if the topography arises due to buoyant loads emplaced at thebase of the ice shell (bottom loading), larger elastic thicknesses require deeper
compensation depths in order to be compatible with the observed admittances(Figure 5.4b). This is because bottom loading in the presence of an elastic
lithosphere leads to reduced buoyant uplift, leaving the low-density isostaticroots to influence the gravity signal, reducing the admittance. If the ice shell
has appreciable elastic strength, and was loaded primarily from below, thecompensation depth must be greater than the 20

141

40 km we estimated forthe fully isostatic case. Not all points in Figure 5.4b are possible, however,
because the maximum mantle thickness corresponding to the low end of theadmittance range is less than about 60 km (see Figure 5.5a). Consequently,
we can rule out elastic thicknesses larger than ⇠ 300 m.

5.4.4

Internal Structure

The moment of inertia implied by the convergence illustrated in Figure 5.3,
⇠ 0.335M R2 , indicates a substantial degree of diﬀerentiation. We approximate
Enceladus as a simple two-layer body: a silicate core surrounded by an H2 O
mantle, with density somewhere between that of ice and liquid water. Given
the observed bulk density constraint, ⇢¯ = 1609 kg/m3 , we use equation (5.20)
to obtain a mantle thickness of roughly 50
density around 2250

65 km (Figure 5.5a) and a core2500 kg/m3 (Figure 5.5b). The low core density is anindication that the (presumably silicate) core may be substantially hydrated.

5.5

Conclusions

In spite of the limited nature of the gravity observations, we have been able todraw some useful conclusions. The combination of the large non-hydrostatic
topography and the modest non-hydrostatic gravity suggests substantial compensation, consistent with lateral thickness variations in an ice shell that is
supported isostatically, i.e., “floating” in a subsurface liquid layer. We obtaina self-consistent interior structure model, with an average compensation depth
of 20

40 km, when we assume the moment of inertia is ⇠ 0.335M R2 . Given

the known bulk density (1609 kg/m3 ), our results are consistent with an Enceladus consisting of a hydrated silicate core of radius 185
142

200 km and density

100
90

70
60
50
40
30
20

7
−0 07
9e e−
53 634
4.1 3.5

Shell Thickness, d (km)

80

10
0
0

(a)

20

40

60

80

100

800

1000

Elastic Thickness, T (m)10090

Shell Thickness, d (km)

80
70
60
50
40

7
−0
9e −07
3
15 4e
4. 63
7
5
−0 7
3.
6e −0
2
e
38 2
2. 792
1.

30
20
10

(b)

0
0

200

400

600

Elastic Thickness, T (m)Figure 5.4: Admittance as a function of compensation depth (d) and elastic layer
thickness (T ) assuming the topography arises from: (a) surfaces loads that deflect
the lithosphere downward; or (b) under-plating of the ice shell that deflects thelithosphere upward. The dark blue bands correspond to the admittance estimates
within the 1 region of agreement while the pale blue bands correspond to the 2region. Because compensation takes place at the base of the ice shell, compensation
depth and shell thickness are equivalent.

143

120

rhom = 920 kg/m3
3

rhom = 1000 kg/m

Mantle Thickness (km)

100

80

60

40

20

(a)

0
0.3

0.32

0.34
0.36
Moment of Inertia Factor

4500

0.38

0.4

rhom = 920 kg/m3
rhom = 1000 kg/m3

Core Density (kg/m3)

4000

3500

3000

2500

2000

(b)

0.3

0.32

0.34
0.36
Moment of Inertia Factor

0.38

0.4

Figure 5.5: Mantle Thickness (a) and Core Density (b) as a function of Moment ofInertia for two diﬀerent assumed mantle densities. The dark blue bands illustrate
the moment of inertia factor within the 1 region of agreement while the pale blue
bands correspond to the 2 region (refer to Figure 5.3).

144

2250

2500 kg/m3 , a subsurface ocean layer with mean thickness 15

and an outer ice shell 20

25 km,

40 km thick.

Acknowledgements
We thank the Cassini radio science team and Luciano Iess for inviting us towork with the newly acquired Enceladus gravity data.

145

Appendix A
Spherical Harmonics
Because several diﬀerent conventions are in widespread use, I make my usage
and notations explicit here to avoid ambiguity. Much of this Appendix follows
Wieczorek (2007), although my notation diﬀers slightly.
Any square-integrable function defined over the surface of a sphere, may be
represented by a linear combination of spherical harmonic functions as follows

F (✓, ) =l
1 X
X(Clm cos m + Slm sin m ) Plm (cos ✓)(A.1)

l=0 m=0

where ✓ is colatitude,

is longitude, Plm (cos ✓) are the associated Legendre

functions of degree l and order m, and where Clm and Slm are the corresponding
expansion coeﬃcients. Equation (A.1) is sometimes written more compactlyasF (✓, ) =1 X
lX

flm Ylm (✓, )(A.2)

l=0 m= lwhere flm are the expansion coeﬃcients, and where Ylm (✓, ) are the real146

spherical harmonic functions, defined byYlm (✓, ) =

8
>
>
<P

lm (cos ✓) cos m>>:Plm (cos ✓) sin |m|

m

0

(A.3)

m<0

I will mainly use the former convention, but the latter will be especiallyconvenient for the derivations in Appendix C.
To ensure accurate computations, even across a range of spherical harmonic
degrees, it is often necessary to normalize the spherical harmonic functions.
When working with normalized spherical harmonics, I write

F (✓, ) =l
1 X
XC̄lm cos m + S̄lm sin mP̄lm (cos ✓)(A.4)

l=0 m=0

where P̄lm (cos ✓) are the normalized associated Legendre functions and whereC̄lm and S̄lm are the corresponding normalized expansion coeﬃcients. The
more compact form is written

F (✓, ) =l
1 XXf¯lm Ȳlm (✓, )(A.5)

l=0 m= lwhere f¯lm are the normalized expansion coeﬃcients, and where Ȳlm (✓, ) are
the normalized real spherical harmonic functions, defined byȲlm (✓, ) =8>><P̄lm (cos ✓) cos m>>:P̄lm (cos ✓) sin |m|

m

0

(A.6)

m<0

Several diﬀerent normalizations are in use, but I will discuss only two here.

147

The first is the Schmidt semi-normalization, widely used in geomagnetismP̄lm (µ) =
where

0m

s

(2

m0 )

(l m)!
Plm (µ)(l + m)!

(A.7)

is the Kronecker delta. The second is what I may call thegeodesy 4⇡-normalization, or full normalization, which I use when working
with topography and gravity fields

P̄lm (µ) =s(2m0 ) (2l + 1)

(l m)!
Plm (µ)(l + m)!

(A.8)

The unnormalized associated Legendre functions are given by

Plm (µ) = 1

µ2

m
2

dm
Pl (µ)dµm(A.9)

where Pl (µ) is the degree-l Legendre polynomial, and is given by

Pl (µ) =1 dlµ2
2l l! dµl

1l

(A.10)

Note that caution is required when using the Matlab function legendre()since its conventions diﬀer from the above, and so appropriate corrections
must be made (carefully). For instance, the Matlab definitions also include theCondon-Shortly phase factor, ( 1)m . When requesting normalized associatedLegendre functions from Matlab, this is not an issue because the phase factor
appears both in the definition of Plm and in the normalization, such that they
cancel. For the Schmidt normalization, the normalized associated Legendre
functions obtained from Matlab are the same as in equation (A.7), but for the148

geodesy 4⇡ normalization, the Matlab result givesP̄lm (µ) Matlab =

s✓

1
l+
2

◆

(l m)!
Plm (µ)(l + m)!

(A.11)

Hence, we need the following adjustment when using the full (or geodesy
4⇡) normalizationP̄lm (µ) =

p
2 (2

m0 )

P̄lm (µ) Matlab(A.12)

When using unnormalized associated Legendre functions from Matlab, thefollowing correction is required to account for the Condon-Shortly phase factor
Plm (µ) = ( 1)m (Plm (µ))Matlab(A.13)

Finally, if the function F (✓, ) is known over ✓ and , the fully normalized
expansion coeﬃcients can be obtained via
1f¯lm =
4⇡ˆ ˆF (✓, ) Ȳlm (✓, ) sin ✓d✓d✓

149

(A.14)

Appendix B
Magnetic Potential and Source
ModelsB.1

Magnetic Scalar Potential and Vector Magnetic Fields

Magnetic fields originating within a planetary body, and measured externally,
can be specified by a magnetic scalar potential, V , which satisfies Laplace’s
equationr2 V = 0
and can be represented as a linear combination of spherical harmonic functionsas

V (r, ✓, ) = a

l ⇣ ⌘
1 X
X
a l+1
l=0 m=0r(Clm cos m + Slm sin m ) Plm (cos ✓)150

where r is the radius at which the potential is to be expressed, ✓ is colatitude,
is longitude, Plm (cos ✓) are the associated Legendre functions of degree l andorder m, Clm and Slm are the corresponding expansion coeﬃcients representing
the magnetic scalar potential at a reference radius a.
The vector magnetic field (also called magnetic flux density or magnetic
induction, with units of Teslas) outside the reference radius, a, is given by thenegative gradient of the scalar magnetic potential

B=

rV

In spherical coordinates, the magnetic field vector components are therefore
given by1 X
l⇣ a ⌘l+2
X@V=(l + 1)(Clm cos m + Slm sin m ) Plm (cos ✓)@rr
l=0 m=0

Br =

B✓ =

B =

@V
=
r@✓

@V
r sin ✓@

=

l ⇣ ⌘
1 X
X
a l+2
l=0 m=0r(Clm cos m + Slm sin m )l1
1 X X ⇣ a ⌘l+2
m
(Clm sin m
sin ✓ l=0 m=0r@Plm (cos ✓)@✓Slm cos m ) Plm (cos ✓)Note that these expressions diﬀer slightly from those of Connerney (2007)
due, I believe, to typographical errors in the latter.

B.2

Source Models

Here, we calculate the locations of the maxima and minima in the horizontal magnetic field profiles due to both vertically and horizontally magnetized151sources, as discussed in Chapter 2. We first consider the single dipole sources
illustrated in Figure 2.1 and then repeat the calculations for the line sources
discussed in section 2.4.

B.2.1

Single Dipole Source

If we take the horizontal and vertical positions in the top panels of Figure 2.1
to be the x and z coordinates, respectively, and if the magnetic dipoles areeach located at the origin with magnetic moment m (with units of Am2 ), thenN), at a position,
the resulting magnetic field, B (with units of Teslas or Am

r = xî + y ĵ + z k̂, is given byB=µ03(m · r)r
4⇡

mr2

1
r5

(B.1)

Equation (B.1) is adapted from Blakely (1995), pp. 75. Here, we employ
SI units with µ0 being the magnetic permeability of free space (which has avalue of 4⇡ ⇥ 10 7 AN2 ). For points in the x-z plane, the horizontal componentof the magnetic field (Bh ) is identical to the magnitude of the x-component(Bx ) since, by symmetry, magnetic field lines cannot cross the x-z plane (i.e.,cannot have a y-component).
Considering the case of the vertically oriented dipole (m =the horizontal component of the magnetic field becomesµ0
Bh = |Bx | =
4⇡

3mzx
5

(x2 + z 2 ) 2

!



T

0 0 m

),

(B.2)

Equation (B.2) demonstrates that field lines directly over a vertically oriented dipole are vertical at any altitude because Bh = 0 when x = 0 (this is152also evident in Figure 2.1a). Diﬀerentiating equation (B.2) with respect to x
and setting the result to zero yields the positions of the peaks in the horizontal
field strength profile1
xBh max = ± z
2

(B.3)

This means that the separation between peaks (2x) is equal to the altitude
above the source (z), which is evident in Figure 2.1c.
Considering the case of the horizontally oriented dipole (m =



T

m 0 0

and still being restricted to points in the x-z plane, the vertical (z) componentof the magnetic field becomesµ0
Bv = |Bz | =
4⇡

3mxz
5

(x2 + z 2 ) 2

!

(B.4)

Equation (B.4) demonstrates that field lines directly over a horizontallyoriented dipole are horizontal at any altitude because Bz = 0 when x = 0(Figure 2.1b). Again, starting from equation (B.1), and considering only points
in the x-z plane, the horizontal component of the magnetic field due to ahorizontally oriented dipole that is aligned with the x-axis becomesµ0
Bh = |Bx | =4⇡m (2x2z2)
5(x2 + z 2 ) 2

!

(B.5)

Setting equation (B.5) to zero shows that the horizontal field strength drops
to zero when1xBh =0 = ± p z
2

(B.6)

Hence, the cusps in the horizontal magnetic field profile are laterally disp
placed from the source by a distance of 1/ 2 times the altitude above the

153

),

source (this is evident in Figure 2.1d).

B.2.2

Linear Source

Instead of a single dipole, consider a source that is a linear structure such asa long, thin, uniformly magnetized cylinder. If this linear source is infinitely
long and coincident with the y-axis, and the magnetization per unit length ism0 , then the magnetic field is given byB=µ0
4(m0 · r)r
4⇡

2m0 r2

1
r4

(B.7)

Equation (B.7) is adapted from Blakely (1995), pp. 96. Again, considering
only points in the x-z plane, the horizontal component of the magnetic field
due to the vertically magnetized line source becomesµ0
Bh = |Bx | =4⇡✓4m0 zx(x2 + z 2 )2

◆

(B.8)

As with the dipole case, field lines directly over a vertically magnetizedline source are vertical at any altitude. Diﬀerentiating equation (B.8) withrespect to x and setting the result to zero yields the positions of the peaks in
the horizontal field profile1
xBh max = ± p z
3

(B.9)

Considering the case of a horizontally magnetized line source, the vertical(z) component of the magnetic field becomesµ0
Bv = |Bz | =
4⇡

✓

154

4m0 xz(x2 + z 2 )2◆

(B.10)

Hence, field lines directly over a horizontally magnetized linear source arehorizontal at any altitude. Again, starting from equation (B.7), and considering only points in the x-z plane, the horizontal component of the magnetic field
due to a linear source coincident with the y-axis and magnetization parallel tothe x-axis becomesµ0
Bh = |Bx | =4⇡✓2m0 (x2 z 2 )(x2 + z 2 )2◆

(B.11)

Setting equation (B.11) to zero shows that the horizontal field strength
drops to zero whenxBh =0 = ±z

(B.12)

Hence, the cusps in the horizontal magnetic field profile are laterally displaced from the line source by a distance equal to the altitude above thesource.

155

Appendix C
Gravitational Potential andFigures of Equilibrium
C.1

Gravitational Potential

If U is the work per unit mass done by an object’s gravitational field, then theacceleration due to gravity is given by

rU

g=

(C.1)

where U is the gravitational potential, or Newtonian potential, and is givenby
U (r) =G

⇢ (r0 )
dV 0
0|
|r
r
V

ˆ

(C.2)

where the integral is performed over the object’s entire volume V , the density
at position r0 is ⇢ (r0 ), and where G is the universal gravitational constant.
Everywhere outside the object, the gravitational potential satisfies Laplace’s

156

equation (Kaula, 2000)(C.3)

r2 U (r) = 0

and the potential can be expressed as a linear combination of spherical harmonic functions as

U (r, ✓, ) =1l

GM X Xr l=0 m= l

✓

R0
r

◆l

g
Clm
Ylm (✓, )(C.4)

where r is the radius at which the potential is to be expressed, ✓ is colatitude,gare the degree-l and order-m spherical harmonic exis longitude, and Clm

pansion coeﬃcients representing the dimensionless gravitational potential at areference radius R0 (see Appendix A).
Following equation (C.1), to get the radial part of the gravitational acceleration, we diﬀerentiate with respect to r, obtaining
1l

GM X Xgr (r, ✓, ) = 2
r l=0 m= l

✓

R0
r

◆l

g
(l + 1) Clm
Ylm (✓, )(C.5)

where we have followed the usual convention of making gravity positive downward.

C.2

Tidal and Rotational Deforming Potentials

C.2.1

Tidal Deforming Potential

A satellite in synchronous orbit around its parent body will experience permanent tidal distortion due to the variation in potential caused by the nearby
tide-raising body. Following from Murray and Dermott (1999), the potential,157

V , at a point, P , is
V =

GM
d

where M is the mass of the tide-raising body and d is the distance between
P and the center of the tide-raising body. If the center-to-center distance
between the two bodies is a, the body’s radius is R, and P makes an angle
with the centerline between the two bodies, then (by the cosine rule)d 2 = a2 + R 2

2aR coswhich can be written as✓ ◆2
Rd=a 1+
a

!1/2
✓ ◆
R
cos
2
a

Expanding this, we obtain
V tid =

"
✓ ◆
R
GM
1+
cos
a
a

✓ ◆2
R
1
+
3 cos2a2

1 + ...

#(C.6)

The first term in (C.6) is a constant and therefore does not contribute totidal distortion. The second term in (C.6) relates to circular motion of thesatellite around its parent body and therefore does not contribute to tidal
distortion. The tide-raising part of the potential begins with the third term.
Neglecting higher order terms, the tide-raising part of the potential is
V tid =

GM R2 13 cos2
a3 2

1 =

GM R2
P2 (cos )a3

If the tide-raising body is much more massive than the synchronous satellite

158

(M

m), and the angular frequency is !, then ! 2 = GM/a3 . Hence
V tid =

! 2 R2 P2 (cos )Assuming the tidal axis is perpendicular to the spin axis and converting tocoordinates of colatitude, ✓, and longitude, , we have cos= cos sin ✓, sothat
V tid =

(C.7)

! 2 R2 P2 (cos sin ✓)This can be re-written in terms of the (non-normalized) spherical harmonic
functions, Y20 and Y22 . First, we look at P2 (cos sin ✓)P2 (cos sin ✓) =1
3 cos2 sin2 ✓2

1which, after some manipulation, can be written

P2 (cos sin ✓) =

3
cos 2
4

But because Y20 = P2 (cos ✓) = 12 (3 cos2 ✓3 (113 cos2 ✓4cos2 ✓

111) and Y22 = P22 (cos ✓) cos (2 ) =cos2 ✓) cos (2 ), we have1P2 (cos sin ✓) = Y22
4

Hence
V

tid

2

=! R

2



1
Y20
2

159

1
Y20
2

1
Y22
4

(C.8)

C.2.2

Rotational Deforming Potential

For a body rotating with angular velocity !, the centrifugal acceleration is
acf = ! 2 x, where x = R sin ✓ is a vector pointed outward from, and perpendicular to, the axis of rotation and reaching the surface at colatitude ✓. Sinceacf =rV cf , the centrifugal potential is thus
V cf =

Because P2 (cos ✓) = 12 (3 cos2 ✓1 2 2 2
! R sin ✓21), we can write this as1
V cf = ! 2 R2 [P2 (cos ✓)
3

1]

(C.9)

Since the last term in (C.9) is a constant, it does not contribute to rotational distortion. The relevant part of the rotational potential is thus
1
V rot = ! 2 R2 P2 (cos ✓)3(C.10)

This can be re-written simply in terms of Y20 = P2 (cos ✓) as1
V rot = ! 2 R2 Y20
3

C.2.3

(C.11)

Combined Deforming Potential

A body subjected to both tidal and rotation distortion will experience a total
disturbing potential V tot = V rot + V tid . Combining (C.7) and (C.10), we have

V

tot

2

=! R

2



1
P2 (cos ✓)3160

P2 (cos sin ✓)Or, in terms of spherical harmonic functions, we combine (C.8) and (C.11)
to get
V

tot

2

=! R

2



1
1
Y20 + Y20
3
2

1
Y22
4

where the first term is the result of rotation and the second and third terms
are due to the permanent tide. This simplifies to

V

C.3

tot

2

=! R

2



5
Y20
6

1
Y22
4

(C.12)

Equilibrium Figures and Potentials

A body subjected to a deforming potential, V , will deform to some degree,depending on the body’s rigidity and internal mass distribution, attributes
reflected in the Love number h. In determining the equilibrium shape of arelaxed hydrostatic body, we are concerned with the fluid Love number hf ,and the deformed shape is given by (Monk and MacDonald , 1960)
H hyd (✓, ) = R

hfVg(C.13)

The change in shape will result in a further change in the gravitational
potential, kf V , which defines the Love number kf , so that the non-central
parts of the body’s hydrostatic gravity are given by
U hyd (✓, ) = kf V

(C.14)

where, for a fluid body, hf = 1 + kf .
For a perfectly fluid body, we can relate the fluid Love number, h2f , to161the body’s normalized moment of inertia using the Radau-Darwin equation(Murray and Dermott, 1999), which can be written
2
I
=
M R2
3

2
5

1

s

C.3.1

Rotational Deformation

C.3.1.1

Hydrostatic Shape

5
h2f

1

!

(C.15)

A planet with no permanent tidal elongation will have an equilibrium shapedetermined by rotational distortion. Here we treat the disturbing potential
as a degree-2 harmonic and so we are concerned with the degree-2 fluid Love
number h2f . Substituting equation (C.11) into equation (C.13), we obtainHhyd(✓, ) = R1
! 2 R2h2f
Y20
g
3(C.16)

If the hydrostatic topography is represented in spherical harmonics as
H hyd (✓, ) =

1 X
l
X

h,hyd
Clm
Ylm (✓, )l=0 m= lthen the non-zero coeﬃcients for the hydrostatic shape are
h,hyd
=R
C00

h,hyd
=
C20

1 R2 ! 2h2f
3 g(C.17)(C.18)

Note that the ratio between the centrifugal and gravitational accelerations162

at the surface of the body is often written

q=

R! 2
g

(C.19)

so we can also write
1
Rqh2f
3

h,hyd
=
C20

The above hydrostatic shape (predicted based on rotational distortion, accounting for the fluid Love number h2f ) is often defined by the semi-axes of anellipsoid (a = b > c). This is done simply by evaluating the above hydrostatic
shape at three locations: at the spin pole, when ✓ = 0 (gives the magnitudeof the c-axis), and at two points on the equator: when ✓ = ⇡2 ,and when ✓ = ⇡2 ,= ⇡2 (b-axis),= 0 (a-axis). The resulting semi-axes are✓1
a = b = R 1 + qh2f
6
✓

1
qh2f
3

c=R 1

◆

◆

From which we haveacR

C.3.1.2

1= qh2f2

Hydrostatic Gravity

Substituting equation (C.11) into equation (C.14), we obtain

U

hyd

2

2

(✓, ) = ! R k2f



1
Y20
3

(C.20)

If the hydrostatic gravitational potential is represented (at the reference163

radius) in spherical harmonics as
U hyd (✓, ) =1l

GM X X g,hydC
Ylm (✓, )R l=0 m= l lm

g,hyd
are the dimensionless gravitational potential coeﬃcients, then,where Clm

, the only non-zero coeﬃcient for thenoting that the surface gravity is g = GMR2hydrostatic gravitational potential is
g,hyd
=
C20

1 !2R
k2f
3 g

(C.21)

2

Again using q = R!g , and noting that the negative (2,0) coeﬃcient is also
called J2 , so we can writeJ2 =C.3.1.3

1
g,hyd
C20
= qk2f
3

Geoid

The geoid is the equipotential surface corresponding to the average potential ata reference radius: Uref = GM/R. The shape of the geoid is a function of theactual mass distribution within the body (captured in the potential coeﬃcients,g, and including rotational deformation) plus the eﬀects of rotation on theClmequipotential surface.
For a body with no permanent tide, the geoid is given by

H

geoid

(✓, ) = R + R

l
1 X
X

l=2 m= l

g
Clm
Ylm

1 R2 ! 2Y20
3 g(C.22)

Hence, the spherical harmonic coeﬃcients for the geoid are equal to the164potential coeﬃcients multiplied by R, with the exception of the (0,0) and(2,0) coeﬃcients, which are
h,geoid=R
C00◆✓
1 R! 2h,geoidg
C20= R C20
3 g

C.3.1.4

Reference Ellipsoid

The C20 term, which arises primarily from rotational distortion, typically dominates the geoid. To highlight those parts of the geoid that are caused bymass anomalies rather than rotational distortion, the geoid is sometimes shownwith respect to a reference ellipsoid. In general, an ellipsoid is triaxial withsemi-axes a > b > c, but for a body subjected only to rotational distortion,a = b > c (see section C.3.2.4 for the reference ellipsoid corresponding toa synchronous satellite, which experiences both rotational and tidal distortion). The reference ellipsoid is sometimes expressed as a degree-2 shape thatmatches the geoid in its C20 term alone. That is

H

ref

✓

(✓, ) = R 1 +

✓

g
C20

1 !2R
3 g

◆

Y20

◆

(C.23)

which is evaluated at three locations to get the semi-axes of the “equivalent”
ellipsoid
a=H

ref

⇣⇡
2

⌘

✓

,0 = R 1

1 g
1 !2R
C20 +
2
6 g

◆

✓
◆
⇣⇡ ⇡ ⌘
1 g
1 !2R
,
=R 1
C +
b=H
2 2
2 20 6 g
◆
✓
1 !2R
g
ref
c = H (0, 0) = R 1 + C20
3 g
ref

165

The ellipsoid’s surface is given byabc
H ref (✓, ) = q(bc sin ✓ cos )2 + (ca sin ✓ sin )2 + (ab cos ✓)2(C.24)

Strictly speaking, an ellipsoid has power at all even harmonics, meaningthat the reference ellipsoid given by equation (C.24) is not identical to thepurely degree-2 shape described by equation (C.23), except at the pole and at
the equator. Because we are usually working with spherical harmonics, equation (C.23) is easier to use and so the geoid with respect to the “reference ellipsoid” is often obtained by subtracting equation (C.23) from equation (C.22).
That isH

geoid-ref(✓, ) = R

l
1 X
X

g
Clm
Ylm

g
RC20
Y20

(C.25)

l=2 m= lHence, the spherical harmonic coeﬃcients for the geoid with respect tothe “reference ellipsoid” are equal to the dimensionless gravitational potential
coeﬃcients multiplied by R, with the exception of the (2,0) coeﬃcient, whichis zero.

C.3.2

Tidal/Rotational Deformation

C.3.2.1

Hydrostatic ShapeA synchronous, or tidally-locked, satellite will have an equilibrium shape determined by both tidal and rotational distortion. Since we are considering onlya degree-2 disturbing potential, we are concerned with the degree-2 fluid Love
number h2f . Substituting equation (C.12) into equation (C.13), we obtain

H

hyd


1
! 2 R2h2f
Y22(✓, ) = R +
g
4
166

5
Y20
6

(C.26)

If the hydrostatic topography is represented in spherical harmonics asHhyd(✓, ) =

1 X
l
X

h,hyd
Clm
Ylm (✓, )l=0 m= lthen the non-zero coeﬃcients for the hydrostatic shape are
h,hyd
=R
C00

(C.27)

h,hyd
=
C20

5 R2 ! 2h2f
6 g(C.28)

h,hyd
C22
=

1 R2 ! 2h2f
4 g(C.29)The ratio between the degree-2 coeﬃcients is thus
h,hyd
C20

h,hyd
C22

=

10
3

Note also that, because Y20 = P2 (cos ✓) = 12 (3 cos2 ✓P22 (cos ✓) cos (2 ) = 3 (11) and Y22 =2

cos2 ✓) cos (2 ), and using q = R!g , equation (C.26)

is sometimes writtenHhyd(✓, ) = R + Rqh2f



1
3 cos22+1

1

cos2 ✓

5
6

The above hydrostatic shape (predicted based on tidal and rotational distortion, accounting for the fluid Love number h2f ) is often defined by thesemi-axes of a triaxial ellipsoid (a, b, c). This is done simply by evaluating theabove hydrostatic shape at three locations: at the spin pole, when ✓ = 0 (gives
the magnitude of the c-axis), and at two points on the equator: when ✓ = ⇡2 ,167

= ⇡2 (b-axis), and when ✓ = ⇡2 ,= 0 (a-axis). The resulting semi-axes are7
a = R 1 + qh2f
6

◆

✓
b=R 1

1
qh2f
3

◆

c=R 1

5
qh2f
6

✓

✓

◆

From which we haveacR

= 2qh2fandab
C.3.2.2c=4cHydrostatic Gravity

Substituting equation (C.12) into equation (C.14), we obtain

U

hyd

2

2

(✓, ) = ! R k2f



5
Y20
6

1
Y22
4

(C.30)

If the hydrostatic gravitational potential is represented (at the reference
radius) in spherical harmonics asUhyd(✓, ) =

1l

GM X X g,hydC
Ylm (✓, )R l=0 m= l lm

g,hyd
are the dimensionless gravitational potential coeﬃcients, then,where Clm

, the non-zero coeﬃcients for thenoting that the surface gravity is g = GM
R2

168

hydrostatic gravitational potential are
g,hyd
=
C20

5 !2R
k2f
6 g

(C.31)

g,hyd
=
C22

1 !2R
k2f
4 g

(C.32)

2

Note that, using the shorthand q = R!g , we can also write
J2 =

5
g,hyd
C20
= qk2f
6

The ratio between the degree-2 coeﬃcients is
g,hyd
C20
=
g,hyd
C22

C.3.2.3

10
3

Geoid

The geoid is the equipotential surface corresponding to the average potential
at a reference radius: Uref = GM/R. The shape of the geoid is a functionof the actual mass distribution within the body (captured in the potential
coeﬃcients, and including tidal and rotational deformation) plus the eﬀects oftides and rotation on the equipotential surface.
The synchronous satellite geoid is given by
H geoid (✓, ) = R + R

l
1 X
X

g
Clm
Ylm

l=2 m= l

5 R2 ! 2
1 R2 ! 2
Y20 +
Y22
6 g
4 g

(C.33)

Hence, the spherical harmonic coeﬃcients for the geoid are equal to thedimensionless gravitational potential coeﬃcients multiplied by R, with the169

exception of the (0,0), (2,0) and (2,2) coeﬃcients, which are
h,geoid=R
C00

h,geoid=R
C20h,geoid=R
C22

C.3.2.4

✓

g
C20

✓

5 R! 2
6 g

1 R!
g
C22
+
4 g

2

◆

◆

Reference Ellipsoid

The C20 and C22 terms, which arise primarily from tidal and rotational distortion, typically dominate the geoid. To highlight those parts of the geoidthat are caused by mass anomalies rather than tidal/rotational distortion, thegeoid is sometimes shown with respect to a reference ellipsoid. The reference
ellipsoid is a triaxial ellipsoid with semi-axes a > b > c. It is sometimes expressed as a degree-2 shape that matches the geoid in its C20 and C22 terms
only. That is

H

ref

✓

✓
g
(✓, ) = R 1 + C20

5 !2R
6 g

◆

Y20 +

✓

1 !2R
g
C22
+
4 g

◆

Y22

◆

(C.34)

which is evaluated at three locations to get the semi-axes of the “equivalent”
ellipsoid
a=H
b=H

ref

ref

⇣⇡
2

⌘

✓

,0 = R 1✓⇣⇡ ⇡ ⌘
,
=R 1
2 2

c=H

1 g
7 !2R
g
C20 + 3C22
+
2
6 g

ref

1 g
C
2 20
✓

g
(0, 0) = R 1 + C20

170

◆

1 !2R
3 g
◆
5 !2R
6 g

g
3C22

◆

The ellipsoid’s surface is given byabc
H ref (✓, ) = q(bc sin ✓ cos )2 + (ca sin ✓ sin )2 + (ab cos ✓)2(C.35)

Strictly speaking, an ellipsoid has power at all even harmonics, meaningthat the reference ellipsoid given by equation (C.35) is not identical to thepurely degree-2 shape described by equation (C.34), except along the three
Cartesian axes (where a, b, c are defined). Because we are usually working
with spherical harmonics, equation (C.34) is easier to use and so the geoidwith respect to the “reference ellipsoid” is often obtained by subtracting equation (C.34) from equation (C.33). That isH

geoid-ref(✓, ) = R

1 X
l
X

g
Clm
Ylm

g
RC20
Y20

g
RC22
Y22

(C.36)

l=2 m= lHence, the spherical harmonic coeﬃcients for the geoid with respect to the“reference ellipsoid” are equal to the dimensionless gravitational potential coeﬃcients multiplied by R, with the exception of the (2,0) and (2,2) coeﬃcients,
which are zero.

171

Bibliography
Adams, J. B., and T. B. McCord (1971), Alteration of Lunar Optical Properties: Age and Composition Eﬀects, Science, 171 (3971), 567–571.
Archinal, B. B. A., M. R. Rosiek, R. L. Kirk, and B. L. Redding (2006), The
Unified Lunar Control Network 2005, Tech. rep., USGS, Flagstaﬀ, AZ.
Béghin, C., C. Sotin, and M. Hamelin (2010), Titan’s native ocean revealed
beneath some 45km of ice by a Schumann-like resonance, Comptes Rendus
Geosci., 342 (6), 425–433, doi:10.1016/j.crte.2010.03.003.
Belleguic, V., P. Lognonné, and M. Wieczorek (2005), Constraints on theMartian lithosphere from gravity and topography data, J. Geophys. Res.,
110 (E11), E11,005, doi:10.1029/2005JE002437.
Bills, B. G., and F. Nimmo (2011), Rotational dynamics and internal structureof Titan, Icarus, 214 (1), 351–355, doi:10.1016/j.icarus.2011.04.028.
Black, B. A., J. T. Perron, D. M. Burr, and S. A. Drummond (2012), Estimating erosional exhumation on Titan from drainage network morphology,
J. Geophys. Res., 117 (E8), 1–17, doi:10.1029/2012JE004085.
Blakely, R. J. (1995), Potential Theory in Gravity and Magnetic Applications,Cambridge University Press.
Blewett, D. T., B. R. Hawke, N. C. Richmond, and C. G. Hughes (2007),
A magnetic anomaly associated with an albedo feature near Airy crater inthe lunar nearside highlands, Geophys. Res. Lett., 34 (24), 1–6, doi:10.1029/
2007GL031670.
Blewett, D. T., E. I. Coman, B. R. Hawke, J. J. Gillis-Davis, M. E. Purucker, and C. G. Hughes (2011), Lunar swirls: Examining crustal magnetic
anomalies and space weathering trends, J. Geophys. Res., 116 (E02002), doi:
10.1029/2010JE003656.
Běhounková, M., G. Tobie, G. Choblet, and O. Čadek (2012), Tidally-induced
melting events as the origin of south-pole activity on Enceladus, Icarus,
219 (2), 655–664, doi:10.1016/j.icarus.2012.03.024.
172

Choukroun, M., and C. Sotin (2012), Is Titan’s shape caused by its meteorology and carbon cycle?, Geophys. Res. Lett., 39 (4), 1–5, doi:10.1029/
2011GL050747.
Clark, R. N. (2009), Detection of adsorbed water and hydroxyl on the Moon.,
Science, 326 (5952), 562–4, doi:10.1126/science.1178105.
Cohen, B. A., P. O. Hayne, D. A. Paige, and B. T. Greenhagen (2014), Lunar
Flashlight: Mapping Lunar Surface Volatiles Using a Cubesat, in Annu.
Meet. Lunar Explor. Anal. Gr., vol. 35812, p. 3031, doi:10.1126/science.
1187726.
Coleman, P. J., B. R. Lichtenstein, C. T. Russell, L. R. Sharp, and G. Schubert(1972), Magnetic fields near the moon, in Third Lunar Sci. Conf., pp. 2271–
2286.
Connerney, J. E. (1999), Magnetic Lineations in the Ancient Crust of Mars,
Science, 284 (5415), 794–798, doi:10.1126/science.284.5415.794.
Connerney, J. E. P. (2007), Planetary Magnetism, in Treatise Geophys., chap.
10.07, pp. 243–280.
Dawkins, R. (1976), The Selfish Gene, Oxford Univ. Press, Oxford, U.K.
Domingue, D. L., C. R. Chapman, R. M. Killen, T. H. Zurbuchen, J. A.
Gilbert, M. Sarantos, M. Benna, J. A. Slavin, D. Schriver, P. M. Trávníček,T. M. Orlando, A. L. Sprague, D. T. Blewett, J. J. Gillis-Davis, W. C.
Feldman, D. J. Lawrence, G. C. Ho, D. S. Ebel, L. R. Nittler, F. Vilas, C. M.
Pieters, S. C. Solomon, C. L. Johnson, R. M. Winslow, J. Helbert, P. N.
Peplowski, S. Z. Weider, N. Mouawad, N. R. Izenberg, and W. E. McClintock(2014), Mercury’s Weather-Beaten Surface: Understanding Mercury in theContext of Lunar and Asteroidal Space Weathering Studies, Space Sci. Rev.,
181 (1-4), 121–214, doi:10.1007/s11214-014-0039-5.
Dwyer, C. A., D. J. Stevenson, and F. Nimmo (2011), A long-lived lunar
dynamo driven by continuous mechanical stirring, Nature, 479 (7372), 212–
214, doi:10.1038/nature10564.
Dyal, P., C. W. Parkin, and C. P. Sonett (1970), Apollo 12 Magnetometer:
Measurement of a Steady Magnetic Field on the Surface of the Moon, Science, 169 (3947), 762–764.
Dyment, J., and J. Arkani-Hamed (1998), Equivalent source magnetic dipoles
revisited, Geophys. Res. Lett., 25 (11), 2003–2006, doi:10.1029/98GL51331.

173

Eliason, E., C. Isbell, E. Lee, T. Becker, L. Gaddis, A. McEwen, and M. Robinson (1999), The Clementine UVVIS Global Lunar Mosaic, Tech. rep., U.S.
Geological Survey / Planetary Data System.
Garrick-Bethell, I., B. P. Weiss, D. L. Shuster, and J. Buz (2009), Early lunar
magnetism., Science, 323 (5912), 356–9, doi:10.1126/science.1166804.
Garrick-Bethell, I., J. W. Head, and C. M. Pieters (2011), Spectral properties,
magnetic fields, and dust transport at lunar swirls, Icarus, 212 (2), 480–492,
doi:10.1016/j.icarus.2010.11.036.
Goldberg, D. E. (1989), Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading, Mass.
Goldsby, D. L., and D. L. Kohlstedt (2001), Superplastic deformation of ice:
Experimental observations, J. Geophys. Res., 106 (B6), 11,017–11,030, doi:
10.1029/2000JB900336.
Grasset, O., C. Sotin, and F. Deschamps (2000), On the internal structureand dynamics of Titan, Planet. Space Sci., 48 (7-8), 617–636, doi:10.1016/
S0032-0633(00)00039-8.
Halekas, J. S., D. L. Mitchell, R. P. Lin, S. Frey, L. L. Hood, M. H. Acuña, andA. B. Binder (2001), Mapping of crustal magnetic anomalies on the lunar
near side by the Lunar Prospector electron reflectometer, J. Geophys. Res.,106 (E11), 841–852.
Hapke, B. (2001), Space weathering from Mercury to the asteroid belt, J.
Geophys. Res., 106 (E5), 10,039–10,073.
Harnett, E. M., and R. M. Winglee (2003), 2.5-D fluid simulations of thesolar wind interacting with multiple dipoles on the surface of the Moon, J.
Geophys. Res., 108 (A2), doi:10.1029/2002JA009617.
Hemingway, D., and I. Garrick-Bethell (2012), Magnetic field direction andlunar swirl morphology: Insights from Airy and Reiner Gamma, J. Geophys.
Res., 117 (E10012), doi:10.1029/2012JE004165.
Hemingway, D., F. Nimmo, H. Zebker, and L. Iess (2013), A rigid andweathered ice shell on Titan, Nature, 500 (7464), 550–552, doi:10.1038/
nature12400.
Hemingway, D. J., I. Garrick-Bethell, and M. A. Kreslavsky (2015), Latitudinal
Variation in Spectral Properties of the Lunar Maria and Implications forSpace Weathering, Icarus, doi:10.1016/j.icarus.2015.08.004.
174Hendrix, A. R., K. D. Retherford, G. Randall Gladstone, D. M. Hurley,
P. D. Feldman, A. F. Egan, D. E. Kaufmann, P. F. Miles, J. W. Parker,
D. Horvath, P. M. Rojas, M. H. Versteeg, M. W. Davis, T. K. Greathouse,
J. Mukherjee, A. J. Steﬄ, W. R. Pryor, and S. A. Stern (2012), The lunar
far-UV albedo: Indicator of hydration and weathering, J. Geophys. Res. E
Planets, 117 (12), 1–8, doi:10.1029/2012JE004252.
Holland, J. (1992), Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence, MIT Press, Cambridge, Mass.
Hood, L., and G. Schubert (1980), Lunar Magnetic Anomalies and Surface
Optical Properties, Science, 208, 49–51, doi:10.1126/science.307.5713.1191.
Hood, L., and C. Williams (1989), The lunar swirls-Distribution and possible
origins, in Lunar Planet. Sci. Conf. Proc., vol. 19, pp. 99–113.
Hood, L., C. Russell, and P. Coleman (1981), Contour maps of lunar remanent
magnetic fields, J. Geophys. Res., 86, 1055–1069.
Hood, L. L. (2011), Central magnetic anomalies of Nectarian-aged lunar impact basins: Probable evidence for an early core dynamo, Icarus, 211 (2),1109–1128, doi:10.1016/j.icarus.2010.08.012.
Hood, L. L., and N. A. Artemieva (2008), Antipodal eﬀects of lunar basinforming impacts: Initial 3D simulations and comparisons with observations,
Icarus, 193 (2), 485–502, doi:10.1016/j.icarus.2007.08.023.
Hood, L. L., P. Coleman, and D. Wilhelms (1979), The Moon: Sources of the
Crustal Magnetic Anomalies, Science, 204 (4388), 53–57.
Hood, L. L., A. Zakharian, J. Halekas, D. Mitchell, R. Lin, M. Acuña, and
A. Binder (2001), Initial mapping and interpretation of lunar crustal magnetic anomalies using Lunar Prospector magnetometer data, J. Geophys.
Res., 106 (E11), 27,825–27,839, doi:10.1029/2000JE001366.
Iess, L., N. J. Rappaport, R. A. Jacobson, P. Racioppa, D. J. Stevenson,
P. Tortora, J. W. Armstrong, and S. W. Asmar (2010), Gravity field, shape,and moment of inertia of Titan., Science, 327 (5971), 1367–9, doi:10.1126/
science.1182583.
Iess, L., R. A. Jacobson, M. Ducci, D. J. Stevenson, J. I. Lunine, J. W. Armstrong, S. W. Asmar, P. Racioppa, N. J. Rappaport, and P. Tortora (2012),
The Tides of Titan, Science, 337, 457–459, doi:10.1126/science.1219631.

175

Iess, L., D. J. Stevenson, M. Parisi, D. Hemingway, R. A. Jacobson, J. I.
Lunine, F. Nimmo, J. W. Armstrong, S. W. Asmar, M. Ducci, and P. Tortora (2014), The Gravity Field and Interior Structure of Enceladus, Science,
344 (6179), 78–80, doi:10.1126/science.1250551.
Jeﬀreys, H. (1976), The Earth: its origin, history, and physical constitution,
6th ed., 535–566 pp., Cambridge University Press, New York.
Kargel, J. S. (1998), Physical Chemistry of Ices in the Outer Solar System, inSol. Syst. Ices, edited by B. Schmitt, C. de Bergh, and M. Festou, pp. 3–32,
Kluwer Academic Publishers, Dordrecht.
Kaula, W. M. (2000), Theory of Satellite Geodesy: Applications of Satellites
to Geodesy, Dover.
Klinger, J. (1980), Influence of a Phase Transition of Ice on the Heat and Mass
Balance of Comets, Science, 209 (4453), 271–272.
Kramer, G. Y., S. Besse, D. Dhingra, J. Nettles, R. Klima, I. Garrick-Bethell,
R. N. Clark, J.-P. Combe, J. W. Head, L. a. Taylor, C. M. Pieters, J. Boardman, and T. B. McCord (2011), M3 spectral analysis of lunar swirls and thelink between optical maturation and surface hydroxyl formation at magnetic
anomalies, J. Geophys. Res., 116 (E00G18), doi:10.1029/2010JE003729.
Kraus, H. (1967), Thin Elastic Shells, Wiley, New York.
Kreslavsky, M. A., J. W. Head, G. A. Neumann, M. A. Rosenburg, O. Aharonson, D. E. Smith, and M. T. Zuber (2013), Lunar topographic roughness
maps from Lunar Orbiter Laser Altimeter (LOLA) data: Scale dependence
and correlation with geologic features and units, Icarus, 226 (1), 52–66, doi:
10.1016/j.icarus.2013.04.027.
Kurata, M., H. Tsunakawa, Y. Saito, H. Shibuya, M. Matsushima, and
H. Shimizu (2005), Mini-magnetosphere over the Reiner Gamma magnetic anomaly region on the Moon, Geophys. Res. Lett., 32 (24), 5–8, doi:
10.1029/2005GL024097.
Lawrence, D. J., W. C. Feldman, R. C. Elphic, R. C. Little, T. H. Prettyman, S. Maurice, P. G. Lucey, and A. B. Binder (2002), Iron abundances on the lunar surface as measured by the Lunar Prospector gammaray and neutron spectrometers, J. Geophys. Res., 107 (E12), 5130, doi:
10.1029/2001JE001530.
Le Bars, M., M. A. Wieczorek, O. Karatekin, D. Cébron, and M. Laneuville(2011), An impact-driven dynamo for the early Moon, Nature, 479 (7372),
215–218, doi:10.1038/nature10565.
176

Le Feuvre, M., and M. A. Wieczorek (2011), Nonuniform cratering of the Moon
and a revised crater chronology of the inner Solar System, Icarus, 214 (1),
1–20, doi:10.1016/j.icarus.2011.03.010.
Lin, R. P. (1979), Constraints on the origins of lunar magnetism from electron reflection measurements of surface magnetic fields, Phys. Earth Planet.
Inter., 20, 271–280.
Lin, R. P., D. L. Mitchell, D. W. Curtis, K. A. Anderson, C. W. Carlson,
J. McFadden, M. H. Acuña, L. L. Hood, and A. Binder (1998), Lunar Surface
Magnetic Fields and Their Interaction with the Solar Wind: Results from
Lunar Prospector, Science, 281 (5382), 1480–1484, doi:10.1126/science.281.
5382.1480.
Lopes, R., K. Mitchell, E. Stofan, J. Lunine, R. Lorenz, F. Paganelli, R. Kirk,
C. Wood, S. Wall, L. Robshaw, a.D. Fortes, C. Neish, J. Radebaugh, E. Reffet, S. Ostro, C. Elachi, M. Allison, Y. Anderson, R. Boehmer, G. Boubin,
P. Callahan, P. Encrenaz, E. Flamini, G. Francescetti, Y. Gim, G. Hamilton, S. Hensley, M. Janssen, W. Johnson, K. Kelleher, D. Muhleman, G. Ori,
R. Orosei, G. Picardi, F. Posa, L. Roth, R. Seu, S. Shaﬀer, L. Soderblom,
B. Stiles, S. Vetrella, R. West, L. Wye, and H. Zebker (2007), Cryovolcanicfeatures on Titan’s surface as revealed by the Cassini Titan Radar Mapper,
Icarus, 186 (2), 395–412, doi:10.1016/j.icarus.2006.09.006.
Lucchitta, B. K. (1978), Geologic Map of the North Side of the Moon.
Lucey, P. G., D. T. Blewett, G. J. Taylor, and B. R. Hawke (2000a), Imagingof lunar surface maturity, J. Geophys. Res., 105 (E8), 20,377–20,386.
Lucey, P. G., D. T. Blewett, and B. L. Jolliﬀ (2000b), Lunar iron and titanium
abundance algorithms based on final processing of Clementine ultravioletvisible images, J. Geophys. Res., 105 (E8), 20,297–20,305.
Lucey, P. G., G. A. Neumann, M. A. Riner, E. Mazarico, D. E. Smith, M. T.
Zuber, D. A. Paige, D. B. J. Bussey, J. T. S. Cahill, A. McGovern, P. Isaacson, L. M. Corley, M. H. Torrence, H. J. Melosh, J. W. Head, and E. Song(2014), The global albedo of the Moon at 1064 nm, J. Geophys. Res. Planets,
119, 1665–1679, doi:10.1002/2013JE004592.Received.
McCord, T. B., L. A. Taylor, J.-P. Combe, G. Kramer, C. M. Pieters, J. M.
Sunshine, and R. N. Clark (2011), Sources and physical processes responsible
for OH/H2O in the lunar soil as revealed by the Moon Mineralogy Mapper(M3), J. Geophys. Res., 116, E00G05, doi:10.1029/2010JE003711.

177

McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A.
Wieczorek, R. J. Phillips, G. A. Neumann, O. Aharonson, and J. W. Head
(2002), Localized gravity/topography admittance and correlation spectraon Mars: Implications for regional and global evolution, J. Geophys. Res.,
107 (E12), 1–25, doi:10.1029/2002JE001854.
McKenzie, D. (1994), The Relationship between Topography and Gravity on
Earth and Venus, Icarus, 112, 55–88.
McKinnon, W. B. (2013), The shape of Enceladus as explained by an irregularcore: Implications for gravity, libration, and survival of its subsurface ocean,
J. Geophys. Res. Planets, 118, doi:10.1002/jgre.20122.
Mckinnon, W. B. (2015), Eﬀect of Enceladus’s rapid synchronous spin on
interpretation of Cassini gravity, Geophys. Res. Lett., 42, doi:10.1002/
2015GL063384.
Mitchell, D., J. Halekas, R. Lin, S. Frey, L. Hood, M. Acuña, and a. Binder(2008), Global mapping of lunar crustal magnetic fields by Lunar Prospector,
Icarus, 194 (2), 401–409, doi:10.1016/j.icarus.2007.10.027.
Mitri, G., and A. P. Showman (2008), A model for the temperature-dependenceof tidal dissipation in convective plumes on icy satellites: Implications for
Europa and Enceladus, Icarus, 195 (2), 758–764, doi:10.1016/j.icarus.2008.
01.010.
Monk, W., and G. MacDonald (1960), The Rotation of the Earth, Cambridge
University Press.
Moore, J. M., and R. T. Pappalardo (2011), Titan: An exogenic world?, Icarus,
212 (2), 790–806, doi:10.1016/j.icarus.2011.01.019.
Moore, J. M., A. D. Howard, and P. M. Schenk (2013), Bedrock denudationon Titan: estimates of vertical extent and lateral debris dispersion, in 44th
Lunar Planet. Sci. Conf., p. 1763.
Moore, J. M., A. D. Howard, and A. M. Morgan (2014), The landscape of
Titan as witness to its climate evolution, J. Geophys. Res. Planets, 119,2060–2077, doi:10.1002/2014JE004608.
Moore, W. B., and G. Schubert (2003), The tidal response of Ganymede andCallisto with and without liquid water oceans, Icarus, 166 (1), 223–226, doi:
10.1016/j.icarus.2003.07.001.
Murray, C. D., and S. Dermott (1999), Solar System Dynamics, Cambridge
University Press.
178

Mustard, J. F., L. Li, and G. He (1998), Nonlinear spectral mixture modelingof lunar multispectral data: Implications for lateral transport, J. Geophys.
Res., 103 (E8), 19,419–19,425.
Neish, C., D. Blewett, D. Bussey, S. Lawrence, M. Mechtley, and B. Thomson (2011), The surficial nature of lunar swirls as revealed by the Mini-RF
instrument, Icarus, 215 (1), 186–196, doi:10.1016/j.icarus.2011.06.037.
Neish, C., R. Kirk, R. Lorenz, V. Bray, P. Schenk, B. Stiles, E. Turtle,
K. Mitchell, and A. G. Hayes (2013), Crater topography on Titan: Implications for landscape evolution, Icarus, 223 (1), 82–90, doi:10.1016/j.icarus.
2012.11.030.
Ness, N. F., K. W. Behannon, C. S. Scearce, and S. C. Cantarano (1967),
Early Results from the Magnetic Field Experiment on Lunar Explorer 35,
J. Geophys. Res., 72 (23), 5769–5778.
Nicholas, J. B., M. E. Purucker, and T. J. Sabaka (2007), Age spot or youthfulmarking: Origin of Reiner Gamma, Geophys. Res. Lett., 34 (2), 1–4, doi:
10.1029/2006GL027794.
Nimmo, F. (2004), Non-Newtonian topographic relaxation on Europa, Icarus,
168 (1), 205–208, doi:10.1016/j.icarus.2003.11.022.
Nimmo, F., and B. Bills (2010), Shell thickness variations and the longwavelength topography of Titan, Icarus, 208 (2), 896–904, doi:10.1016/j.
icarus.2010.02.020.
Nimmo, F., R. T. Pappalardo, and B. Giese (2002), Eﬀective elastic thickness
and heat flux estimates on Ganymede, Geophys. Res. Lett., 29 (7), 1158.
Nimmo, F., J. R. Spencer, R. T. Pappalardo, and M. E. Mullen (2007), Shear
heating as the origin of the plumes and heat flux on Enceladus., Nature,447 (7142), 289–91, doi:10.1038/nature05783.
Nimmo, F., B. G. Bills, and P. C. Thomas (2011), Geophysical implicationsof the long-wavelength topography of the Saturnian satellites, J. Geophys.
Res., 116 (E11), E11,001, doi:10.1029/2011JE003835.
Noble, S. K., C. M. Pieters, and L. P. Keller (2007), An experimental approachto understanding the optical eﬀects of space weathering, Icarus, 192 (2), 629–
642, doi:10.1016/j.icarus.2007.07.021.
Nozette, S., P. Rustan, L. P. Pleasance, D. M. Horan, P. Regeon, E. M. Shoemaker, P. D. Spudis, C. H. Acton, D. N. Baker, J. E. Blamont, B. J. Buratti,
M. P. Corson, M. E. Davies, T. C. Duxbury, E. M. Eliason, B. M. Jakosky,
179

J. F. Kordas, I. T. Lewis, C. L. Lichtenberg, P. G. Lucey, E. Malaret,
M. A. Massie, J. H. Resnick, C. J. Rollins, H. S. Park, A. S. Mcewen, R. E.
Priest, C. M. Pieters, R. A. Reisse, M. S. Robinson, R. A. Simpson, D. E.
Smith, T. C. Sorenson, R. W. Vorder Breugge, and M. T. Zuber (1994), The
Clementine Mission to the Moon: Scientific Overview, Science, 266 (5192),
1835–1839.
Pariso, J. E., and H. P. Johnson (1991), Alteration Processes at Deep Sea
Drilling Project/Ocean Drilling Program Hole 504B at the Costa Rica Rift:
Implications for Magnetization of Oceanic Crust, J. Geophys. Res., 96 (B7),11,703–11,722.
Parker, R. L. (2003), Ideal bodies for Mars magnetics, J. Geophys. Res.,108 (E1), 5006, doi:10.1029/2001JE001760.
Patterson, D. B., K. A. Farley, and M. D. Norman (1999), He-4 as a tracer
of continental dust: A 1.9 million year record of aeolian flux to the west
equatorial Pacific Ocean, Geochim. Cosmochim. Acta, 63 (5), 615–625.
Pieters, C. M. (1978), Mare basalt types on the front side of the moon: A
summary of spectral reflectance data, in 9th Lunar Planet. Sci. Conf., pp.
2825–2849.
Pieters, C. M., E. M. Fischer, O. Rode, and A. Basu (1993), Optical Eﬀects
of Space Weathering: The Role of the Finest Fraction, J. Geophys. Res.,
98 (E11), 20,817–20,824.
Pieters, C. M., L. A. Taylor, S. K. Noble, L. P. Keller, B. Hapke, R. V. Morris,
C. C. Allen, D. S. McKay, and S. Wentworth (2000), Space weathering on
airless bodies: Resolving a mystery with lunar samples, Meteorit. Planet.
Sci., 35, 1101–1107.
Pieters, C. M., J. N. Goswami, R. N. Clark, M. Annadurai, J. Boardman,
B. Buratti, J.-P. Combe, M. D. Dyar, R. Green, J. W. Head, C. Hibbitts,
M. Hicks, P. Isaacson, R. Klima, G. Kramer, S. Kumar, E. Livo, S. Lundeen, E. Malaret, T. McCord, J. Mustard, J. Nettles, N. Petro, C. Runyon, M. Staid, J. Sunshine, L. A. Taylor, S. Tompkins, and P. Varanasi(2009), Character and spatial distribution of OH/H2O on the surface ofthe Moon seen by M3 on Chandrayaan-1, Science, 326 (5952), 568–72, doi:
10.1126/science.1178658.
Pieters, C. M., E. Ammannito, D. T. Blewett, B. W. Denevi, M. C. De Sanctis,
M. J. Gaﬀey, L. Le Corre, J.-Y. Li, S. Marchi, T. B. McCord, L. A. McFadden, D. W. Mittlefehldt, A. Nathues, E. Palmer, V. Reddy, C. A. Raymond,180

and C. T. Russell (2012), Distinctive space weathering on Vesta from regolith
mixing processes., Nature, 491 (7422), 79–82, doi:10.1038/nature11534.
Pieters, C. M., D. P. Moriarty, and Ian Garrick-Bethell (2014), Atypical regolith processes hold the key to enigmatic lunar swirls, in 45th Lunar Planet.
Sci. Conf.
Pinet, P. C., V. V. Shevchenko, S. D. Chevrel, Y. Daydou, and C. Rosemberg (2000), Local and regional lunar regolith characteristics at Reiner
Gamma Formation: Optical and spectroscopic properties from Clementine
and Earth-based data, J. Geophys. Res., 105 (E4), 9457–9475.
Porco, C. C., P. Helfenstein, P. C. Thomas, A. P. Ingersoll, J. Wisdom,
R. West, G. Neukum, T. Denk, R. Wagner, T. Roatsch, S. Kieﬀer, E. Turtle,
A. McEwen, T. V. Johnson, J. Rathbun, J. Veverka, D. Wilson, J. Perry,J. Spitale, A. Brahic, J. A. Burns, A. D. Delgenio, L. Dones, C. D. Murray,
and S. Squyres (2006), Cassini observes the active south pole of Enceladus.,
Science, 311 (5766), 1393–1401, doi:10.1126/science.1123013.
Prettyman, T. H., J. J. Hagerty, R. C. Elphic, W. C. Feldman, D. J. Lawrence,
G. W. McKinney, and D. T. Vaniman (2006), Elemental composition of
the lunar surface: Analysis of gamma ray spectroscopy data from Lunar
Prospector, J. Geophys. Res., 111 (E12007), doi:10.1029/2005JE002656.
Purucker, M. (2008), A global model of the internal magnetic field of the Moonbased on Lunar Prospector magnetometer observations, Icarus, 197 (1), 19–23, doi:10.1016/j.icarus.2008.03.016.
Purucker, M. E., and J. B. Nicholas (2010), Global spherical harmonic models of the internal magnetic field of the Moon based on sequential andcoestimation approaches, J. Geophys. Res., 115 (E12), 1–14, doi:10.1029/
2010JE003650.
Purucker, M. E., T. J. Sabaka, and R. A. Langel (1996), Conjugate gradientanalysis: A new tool for studying satellite magnetic data sets, Geophys. Res.
Lett., 23 (5), 507–510.
Richards, M. A., and B. H. Hager (1984), Geoid anomalies in a dynamic Earth,
J. Geophys. Res., 89 (B7), 5987, doi:10.1029/JB089iB07p05987.
Richmond, N. C., L. L. Hood, J. Halekas, D. Mitchell, R. Lin, M. Acuña, andA. B. Binder (2003), Correlation of a strong lunar magnetic anomaly with ahigh-albedo region of the Descartes mountains, Geophys. Res. Lett., 30 (7),1395, doi:10.1029/2003GL016938.

181

Richmond, N. C., L. Hood, D. Mitchell, R. Lin, M. Acuña, and A. Binder(2005), Correlations between magnetic anomalies and surface geology antipodal to lunar impact basins, J. Geophys. Res., 110 (E5), 1–11, doi:
10.1029/2005JE002405.
Roberts, J. H., and F. Nimmo (2008), Tidal heating and the long-term stabilityof a subsurface ocean on Enceladus, Icarus, 194 (2), 675–689, doi:10.1016/j.
icarus.2007.11.010.
Sasaki, S., K. Nakamura, Y. Hamabe, E. Kurahashi, and T. Hiroi (2001),
Production of iron nanoparticles by laser irradiation in a simulation of lunarlike space weathering., Nature, 410 (6828), 555–7, doi:10.1038/35069013.
Schultz, P. H., and L. J. Srnka (1980), Cometary Collisions on the Moon andMercury, in Lunar Planet. Sci. Conf. 11th, pp. 1009–1011.
Shea, E. K., B. P. Weiss, W. S. Cassata, D. L. Shuster, S. M. Tikoo, J. Gattacceca, T. L. Grove, and M. D. Fuller (2012), A long-lived lunar core dynamo,
Science, 335 (6067), 453–6, doi:10.1126/science.1215359.
Shibuya, H., H. Tsunakawa, F. Takahashi, H. Shimizu, M. Matsushima, andK. M.-l. Team (2010), Near surface magnetic field mapping over the swirls
in the SPA region using Kaguya LMAG data, in Eur. Planet. Sci. Congr.,
vol. 5.
Smith, D. E., M. T. Zuber, G. A. Neumann, F. G. Lemoine, E. Mazarico, M. H.
Torrence, J. F. McGarry, D. D. Rowlands, J. W. Head, T. H. Duxbury,
O. Aharonson, P. G. Lucey, M. S. Robinson, O. S. Barnouin, J. F. Cavanaugh, X. Sun, P. Liiva, D.-d. Mao, J. C. Smith, and A. E. Bartels (2010),
Initial observations from the Lunar Orbiter Laser Altimeter (LOLA), Geophys. Res. Lett., 37 (18), 1–6, doi:10.1029/2010GL043751.
Spencer, J. R., and F. Nimmo (2013), Enceladus: An Active Ice World in
the Saturn System, Annu. Rev. Earth Planet. Sci., 41 (1), 693–717, doi:
10.1146/annurev-earth-050212-124025.
Spitale, J. N., and C. C. Porco (2007), Association of the jets of Enceladus
with the warmest regions on its south-polar fractures., Nature, 449 (7163),
695–7, doi:10.1038/nature06217.
Staid, M. I., and C. M. Pieters (2000), Integrated Spectral Analysis of Mare
Soils and Craters: Applications to Eastern Nearside Basalts, Icarus, 145,
122–139, doi:10.1006/icar.1999.6319.

182

Starukhina, L. V., and Y. G. Shkuratov (2004), Swirls on the Moon andMercury: meteoroid swarm encounters as a formation mechanism, Icarus,
167 (1), 136–147, doi:10.1016/j.icarus.2003.08.022.
Stiles, B. W., S. Hensley, Y. Gim, D. M. Bates, R. L. Kirk, A. Hayes,
J. Radebaugh, R. D. Lorenz, K. L. Mitchell, P. S. Callahan, H. Zebker,
W. T. Johnson, S. D. Wall, J. I. Lunine, C. A. Wood, M. Janssen, F. Pelletier, R. D. West, and C. Veeramacheneni (2009), Determining Titan
surface topography from Cassini SAR data, Icarus, 202 (2), 584–598, doi:
10.1016/j.icarus.2009.03.032.
Sunshine, J. M., T. L. Farnham, L. M. Feaga, O. Groussin, F. Merlin, R. E.
Milliken, and M. F. A’Hearn (2009), Temporal and spatial variability of lunar
hydration as observed by the Deep Impact spacecraft., Science, 326 (5952),
565–8, doi:10.1126/science.1179788.
Thomas, P., J. Burns, P. Helfenstein, S. Squyres, J. Veverka, C. Porco,
E. Turtle, A. Mcewen, T. Denk, and B. Giese (2007), Shapes of the saturnian icy satellites and their significance, Icarus, 190 (2), 573–584, doi:
10.1016/j.icarus.2007.03.012.
Tobie, G., O. Grasset, J. I. Lunine, A. Mocquet, and C. Sotin (2005), Titan’s
internal structure inferred from a coupled thermal-orbital model, Icarus,
175 (2), 496–502, doi:10.1016/j.icarus.2004.12.007.
Tobie, G., J. I. Lunine, and C. Sotin (2006), Episodic outgassing as the origin
of atmospheric methane on Titan., Nature, 440 (7080), 61–4, doi:10.1038/
nature04497.
Tricarico, P. (2014), Multi-Layer Hydrostatic Equilibrium of Planets and Synchronous Moons: Theory and Application To Ceres and To Solar System
Moons, Astrophys. J., 782 (2), 99, doi:10.1088/0004-637X/782/2/99.
Turcotte, D. L., R. J. Willemann, W. F. Haxby, and J. Norberry (1981), Roleof Membrane Stresses in the Support of Planetary Topography, J. Geophys.
Res., 86 (1), 3951–3959.
Vernazza, P., R. P. Binzel, A. Rossi, M. Fulchignoni, and M. Birlan (2009),
Solar wind as the origin of rapid reddening of asteroid surfaces, Nature,458 (7241), 993–5, doi:10.1038/nature07956.
Von Frese, R. R. B., W. J. Hinze, and L. W. Braile (1981), Spherical Earth
Gravity and Magnetic Anomaly Analysis by Equivalent Point Source Inversion, Earth Planet. Sci. Lett., 53, 69–83.

183

Werner, S. C., and S. Medvedev (2010), The Lunar rayed-crater population Characteristics of the spatial distribution and ray retention, Earth Planet.
Sci. Lett., 295 (1-2), 147–158, doi:10.1016/j.epsl.2010.03.036.
Wieczorek, M. A. (2007), Gravity and Topography of the Terrestrial Planets, in Treatise Geophys., vol. 10, chap. 10.05, pp. 165–206, doi:10.1016/
B978-044452748-6/00156-5.
Wieczorek, M. A., B. P. Weiss, and S. T. Stewart (2012), An Impactor Origin
for Lunar Magnetic Anomalies, Science, 335 (6073), 1212–1215, doi:10.1126/
science.1214773.
Wilcox, B. B., P. G. Lucey, and J. J. Gillis (2005), Mapping iron in the lunar
mare: An improved approach, J. Geophys. Res., 110 (E11001), doi:10.1029/
2005JE002512.
Yokota, Y., T. Matsunaga, M. Ohtake, J. Haruyama, R. Nakamura, S. Yamamoto, Y. Ogawa, T. Morota, C. Honda, K. Saiki, K. Nagasawa, K. Kitazato, S. Sasaki, A. Iwasaki, H. Demura, N. Hirata, T. Hiroi, R. Honda,
Y. Iijima, and H. Mizutani (2011), Lunar photometric properties at wavelengths 0.5-1.6 micron acquired by SELENE Spectral Profiler and their dependency on local albedo and latitudinal zones, Icarus, 215 (2), 639–660,
doi:10.1016/j.icarus.2011.07.028.
Zebker, H., B. Stiles, S. Hensley, R. Lorenz, R. L. Kirk, and J. Lunine (2009),
Size and shape of Saturn’s moon Titan., Science, 324 (5929), 921–3, doi:
10.1126/science.1168905.
Zebker, H. A., L. Iess, S. D. Wall, R. D. Lorenz, J. I. Lunine, and B. W. Stiles(2012), Titan’s Figure Fatter, Flatter Than Its Gravity Field, in AGU Fall
Meet. 2012, pp. P23F–01.
Zuber, M. T., J. W. Head, D. E. Smith, G. A. Neumann, E. Mazarico, M. H.
Torrence, O. Aharonson, A. R. Tye, C. I. Fassett, M. A. Rosenburg, and H. J.
Melosh (2012), Constraints on the volatile distribution within Shackleton
crater at the lunar south pole., Nature, 486 (7403), 378–81, doi:10.1038/
nature11216.

184