Space Sci Rev (2013) 178:25–55
DOI 10.1007/s11214-013-9962-0
The Scientific Measurement System of the Gravity
Recovery and Interior Laboratory (GRAIL) Mission
Sami W. Asmar · Alexander S. Konopliv · Michael M. Watkins · James G. Williams ·
Ryan S. Park · Gerhard Kruizinga · Meegyeong Paik · Dah-Ning Yuan ·
Eugene Fahnestock · Dmitry Strekalov · Nate Harvey · Wenwen Lu · Daniel Kahan ·
Kamal Oudrhiri · David E. Smith · Maria T. Zuber
Received: 5 October 2012 / Accepted: 22 January 2013 / Published online: 21 February 2013
© Springer Science+Business Media Dordrecht 2013
Abstract The Gravity Recovery and Interior Laboratory (GRAIL) mission to the Moon uti-
lized an integrated scientific measurement system comprised of flight, ground, mission, and
data system elements in order to meet the end-to-end performance required to achieve its
scientific objectives. Modeling and simulation efforts were carried out early in the mission
that influenced and optimized the design, implementation, and testing of these elements.
Because the two prime scientific observables, range between the two spacecraft and range
rates between each spacecraft and ground stations, can be affected by the performance of
any element of the mission, we treated every element as part of an extended science in-
strument, a science system. All simulations and modeling took into account the design and
configuration of each element to compute the expected performance and error budgets. In
the process, scientific requirements were converted to engineering specifications that be-
came the primary drivers for development and testing. Extensive simulations demonstrated
that the scientific objectives could in most cases be met with significant margin. Errors are
grouped into dynamic or kinematic sources and the largest source of non-gravitational er-
ror comes from spacecraft thermal radiation. With all error models included, the baseline
solution shows that estimation of the lunar gravity field is robust against both dynamic and
kinematic errors and a nominal field of degree 300 or better could be achieved according to
the scaled Kaula rule for the Moon. The core signature is more sensitive to modeling errors
and can be recovered with a small margin.
Keywords Gravity · Moon · Remote sensing · Spacecraft
Acronyms and Abbreviations
AGC Automatic Gain Control
S.W. Asmar () · A.S. Konopliv · M.M. Watkins · J.G. Williams · R.S. Park · G. Kruizinga · M. Paik ·
D.-N. Yuan · E. Fahnestock · D. Strekalov · N. Harvey · W. Lu · D. Kahan · K. Oudrhiri
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
e-mail: sami.asmar@jpl.nasa.gov
D.E. Smith · M.T. Zuber
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology,
Cambridge, MA 02139-4307, USA
26 S.W. Asmar et al.
AMD Angular Momentum Desaturation
C&DH Command & Data Handler
CBE Current Best Estimate
CG Center of Gravity
CPU Central Processing Unit
DSN Deep Space Network
DOWR Dual One-way Range
ECM Eccentricity Correction Maneuver
EOP Earth Orientation Platform
ET Ephemeris Time
GPS Global Positioning System
GR-A GRAIL-A Spacecraft (Ebb)
GR-B GRAIL-B Spacecraft (Flow)
GRACE Gravity Recovery and Climate Experiment
GRAIL Gravity Recovery and Interior Laboratory
GSFC Goddard Space Flight Center
ICRF International Celestial Reference Frame
IERS International Earth Rotation and Reference Systems Service
IR Infra Red
IPU Instrument Processing Unit (GRACE mission)
JPL Jet Propulsion Laboratory
KBR Ka-Band Ranging
KBRR Ka-Band Range-Rate
LGRS Lunar Gravity Ranging System
LLR Lunar Laser Ranging
LOI Lunar Orbit Insertion
LOS Line of Sight
LP Lunar Prospector
mGal milliGal (where 1 Gal = 0.01 m s−2)
MGS Mars Global Surveyor
MIT Massachusetts Institute of Technology
MOS Mission Operations System
MIRAGE Multiple Interferometric Ranging and GPS Ensemble
MMDOM Multi-mission Distributed Object Manager
MONTE Mission-analysis, Operations, and Navigation Toolkit Environment
MPST Mission Planning and Sequence Team
MRO Mars Reconnaissance Orbiter
MWA Microwave Assembly
NASA National Aeronautics and Space Administration
ODP Orbit Determination Program
OPR Orbital Period Reduction
OSC Onboard Spacecraft Clocks
OTM Orbit Trim Maneuver
PDS Planetary Data System
PM Primary Mission
PPS Pulse Per Second
RSB Radio Science Beacon
RSR Radio Science Receiver
SCT Spacecraft Team
The GRAIL Scientific Measurement System 27
SDS Science Data System
SIS Software Interface Specification
SRIF Square Root Information Filter
SRP Solar Radiation Pressure
TAI International Atomic Time
TCM Trajectory Correction Maneuver
TDB Barycentric Dynamic Time
TDS Telemetry Delivery System
TDT Terrestrial Dynamic Time
TLC Trans-Lunar Cruise
TSF Transition to Science Formation
TSM Transition to Science Maneuver
TTS Time Transfer System
USO Ultra-stable Oscillator
UTC Universal Time Coordinated
VLBI Very Long Baseline Interferometry
1 Introduction and Heritage
The Gravity Recovery and Interior Laboratory (GRAIL) mission is comprised of two space-
craft, named Ebb and Flow, flying in precision formation around the Moon. The mission’s
purpose is to recover the lunar gravitational field in order to investigate the interior structure
of the Moon from the crust to the core. The spacecraft were launched together on Septem-
ber 10, 2011 and began science operations and data acquisition on March 1, 2012. Zuber
et al. (2013, this issue) presents an overview of the mission including scientific objectives
and measurement requirements. Klipstein et al. (2013, this issue) describes the design and
implementation of the GRAIL payload. Hoffman (2009) described GRAIL’s flight system
and Roncoli and Fujii (2010) described the mission design.
This paper illustrates how a team of scientists and engineers prepared to meet GRAIL
scientific objectives and data quality requirements through simulations and modeling of the
design and configuration of the flight and ground systems. It details dynamic and kinematic
models for estimating error sources in the form of non-gravitational forces and how these
models were applied, along with the lunar gravity model, to elaborate computer simulations
in the context of an integrated scientific measurement system. This paper also documents the
methods, tools, and results of the simulations. This work was carried out at the Jet Propulsion
Laboratory (JPL) prior to the science orbital phase and reviewed by expert peers from differ-
ent institutions; the knowledge is based on the combined experiences of the team members
with gravity observations on numerous planetary missions. This effort demonstrated that the
mission was capable of meeting the science requirements as well as paved the way to the
operational tools and procedures for the actual science data analysis.
The GRAIL concept was derived from the Gravity Recovery and Climate Experiment
(GRACE) Earth mission and utilized a modified GRACE payload called the Lunar Gravity
Ranging System (LGRS); the GRAIL and GRACE spacecraft are unrelated. For an overview
of the GRACE mission see Tapley et al. (2004a, 2004b); for a description of the GRACE
payload, see Dunn et al. (2003); and for error analysis in the GRACE system and measure-
ments, see Kim and Tapley (2002).
Despite the high heritage, there are significant differences between the GRAIL and
GRACE science payloads, listed in Table 1. GRACE is equipped with a Global Positioning
28 S.W. Asmar et al.
Table 1 Functional differences between the GRAIL and GRACE Missions
GRAIL GRACE
Target body Moon Earth
Launch vehicle Delta II, USA Rockot, Russia
Nominal prime mission duration 3 months 5 years
Orbiter mass (kg) 313 487
Launch date 9/10/11 3/17/02
Prime mission mean orbital altitude (km) 55 470
Gravity coefficients 420 120
Timing synchronization method RSB GPS
Science-quality accelerometer N Y
Adjustable mass for accelerometer at CG N Y
Center of gravity calibrations for antenna Y Y
Inter-spacecraft links Ka-/S-band Ka-/K-band
Spacecraft separation distance (km) 85–225 170–270
Attitude control Reaction wheels Magnetic torque
Thrusters gas Hydrazine Nitrogen
Star cameras per spacecraft 1 2
Science processor Single String Redundant
Star camera software host C&DHa IPUb
USOs per spacecraft 1 2
Absolute timing accuracy DSN: millisecond GPS: nanosecond
Relative timing accuracy TTS: picosecond GPS: picosecond
Communication stations DSN German stations
aC&DH is GRAIL’s Command and Data Handling Subsystem.
bIPU is GRACE’s Instrument Processing Unit.
System (GPS) receiver for timing synchronization, and accelerometers for non-gravitational
force calibrations, while GRAIL is not. Furthermore, GRACE inter-spacecraft ranging uti-
lizes two radio links at K- and Ka-bands (∼26 GHz and ∼32 GHz, respectively) in order to
calibrate the effects of charged particles in the Earth ionosphere, while GRAIL utilizes only
one Ka-band link. In lieu of GPS time synchronization, which is not available at the Moon,
GRAIL introduced two elements, a second inter-spacecraft link at S-band (∼2.3 GHz) for
a Time Transfer System (TTS), and a one-way X-band (∼8.4 GHz) link transmitted from
each spacecraft’s Radio Science Beacon (RSB) to the Deep Space Network (DSN) stations.
With these differences, the GRAIL observable time tagging and synchronization is handled
differently from the GRACE GPS-based system as will be discussed below.
Furthermore, while the GRACE observables are referenced to a geocentric frame,
GRAIL measurements are referenced to Ephemeris Time (ET) and the solar system barycen-
tric frame of Barycentric Dynamic Time (TDB). Finally, since GRAIL does not carry an
accelerometer, attention was given in the design, assembly, and testing of the spacecraft sys-
tem in order to minimize on the non-gravitational forces acting on the spacecraft, including
the solar radiation pressure, lunar albedo and spacecraft outgassing.
All radio signals in the science payload, illustrated in Fig. 1, the Ka-band inter-spacecraft
link, the S-band TTS inter-spacecraft link, and the X-band RSB link to Earth, are referenced
The GRAIL Scientific Measurement System 29
Fig. 1 The GRAIL radio links: Ka-band and S-band inter-spacecraft links, X-band one-way downlink to
ground stations, and two-way S-band links for telecommunications and navigation
on one Ultra-Stable Oscillator (USO) per spacecraft. The navigation and telecommunica-
tions telemetry and command functions are handled by a separate two-way S-band link with
the DSN. This is a spacecraft system function not linked to the science payload or the USO.
The science data quality simulations did not incorporate the utility of this telecommunica-
tions link since no science performance requirements were imposed on it, but in reality, the
science team collaborated with the navigation team to assess its usability to enhance the
science results.
2 Simulations Tool and Data Levels
Over several decades, NASA’s JPL has developed techniques, algorithms, and software tools
to conduct investigations of planetary gravitational fields and applied them to practically ev-
ery planet in the solar system and several satellites of the outer planets. JPL relies primarily
on the Orbit Determination Program (ODP) whose formulation is detailed in Moyer (2003).
The ODP has enabled precision navigation for the vast majority of deep space missions and,
due to its criticality to the success of these missions, has received rigorous development
and testing as well as continued improvements (a new tool called Mission-analysis, Opera-
tions, and Navigation Toolkit Environment, or MONTE, has replaced the ODP for mission
navigation purposes).
30 S.W. Asmar et al.
Fig. 2 A functional flowchart of the MIRAGE software tool as used in the simulation process
The GRAIL scientists at JPL use a version of the ODP called Multiple Interferometric
Ranging and GPS Ensemble (MIRAGE), which originated from a GPS version of the ODP
developed for the TOPEX mission (described in Guinn and Wolff 1993, and Leavitt and
Salama 1993) and further developed for gravity field analysis, (Fahenstock 2009). Figure 2
shows the MIRAGE flowchart process utilized for GRAIL and the various programs that
process the generalized inputs, the spacecraft path integration, computation of dynamic pa-
The GRAIL Scientific Measurement System 31
rameter partials, and the data observables. This figure documents the necessary interfaces
between the software elements and input/output files as well as the relevant computational
parameters and has been a key figure for the simulations peer-review process. There are three
subsets of programs that integrate the spacecraft motion, process the spacecraft observations
and filter or estimate the spacecraft state and related parameters using the observations.
To determine the spacecraft dynamical path, the program numerically integrates the
spacecraft Cartesian state by including all known forces acting on the spacecraft, such
as gravity, solar pressure, lunar albedo, and spacecraft thrusting. The spacecraft state and
the force model partial derivatives (e.g., gravity harmonics) that are later estimated are
integrated using the variable order Adams method described in Krogh (1973). The non-
rotating International Celestial Reference Frame (ICRF) defines the inertial coordinate sys-
tem, which is nearly equal to the Earth’s mean equator and equinox at the epoch of J2000.
The GRAIL data are categorized in 3 levels, also shown in of Zuber et al. (2013, this
issue). Level 0 is the raw data acquired by the spacecraft science payload, the LGRS, and
DSN Doppler. Level 1 is the expanded, edited and calibrated data. Level 1 processing is the
conversion from Level 0 files to Level 1 files. Level 1 processing also applies a time tag
conversion, time of flight correction, and phase center offset, as well as generates instan-
taneous range-rate and range-acceleration observables by numerical differentiation of the
biased range observables. Level 2 is the gravity field spherical harmonic expansion; level 2
processing refers to the production of Level 2 data. The simulations described herein emu-
late the generation of Levels 1 and 2 GRAIL mission data.
3 Gravity Model Representation
Gravitational fields provide a key tool for probing the interior structure of planets. The lunar
gravity, when combined with topography, leads to geophysical models that address impor-
tant phenomena such as the structure of the crust and lithosphere, the asymmetric lunar
thermal evolution, subsurface structure of impact basins and the origin of mascons, and the
temporal evolution of crustal brecciation and magmatism. Long-wavelength gravity mea-
surements can place constraints on the presence of a lunar core.
A gravitational field represents variations in the gravitational potential of a planet and
gravity anomalies at its surface. It can be mathematically represented via coefficients of a
spherical harmonic expansion whose degree and order reflect the surface resolution. A field
of degree 180, for example, represents a half-wavelength, or spatial block size, surface reso-
lution of 30 km; for degree n, the resolution is 30×180/n km. The gravitational potential in
spherical harmonic form is represented in the body-fixed reference frame with normalized
coefficients (Cnm, Snm) is represented after Heiskanen and Moritz (1967) and Kaula (1966)
as:
U = GM
r
+ GM
r
∞∑
n=1
n∑
m=0
(
Re
r
)n
P nm(sinϕlat )
[
Cnm cos(mλ) + Snm sin(mλ)
] (1)
G is the gravitational constant, M is the mass of the central body, r is the radial distance
coordinate, m is the order, P nm are the fully normalized associated Legendre polynomials,
Re is the reference radius of the body, ϕlat is the latitude, and λ is the longitude. The gravity
coefficients are normalized so that the integral of the harmonic squared equals the area of
a unit sphere, and are related to the un-normalized coefficients by Kaula (1966), where δ is
32 S.W. Asmar et al.
the Kronecker delta:
(
Cnm
Snm
)
=
[
(n − m)!(2n + 1)(2 − δ0m)
(n + m)!
]1/2(
Cnm
Snm
)
= fnm
(
Cnm
Snm
)
(2)
There exist singularities at the pole in the partials of the gravity acceleration with respect to
the spacecraft position when using the Legendre polynomials as a function of latitude. To
accommodate this, MIRAGE uses a nonsingular formulation of the gravitational potential,
including recursion relations given by Pines (1973), in calculation of the acceleration and
partials.
The gravitational potential also accounts for tides caused by a perturbing body. The
second-degree tidal potential acting on a satellite at position r relative to the central body,
with the perturbing body (e.g., Sun and Earth for GRAIL) at position rp , is:
U = k2 GMp
R
R6
r3r3p
[
3
2
(rˆ · rˆp)2 − 12
]
(3)
where k2 is the second degree potential Love number, Mp is the mass of the perturbing body
causing the tide, and R is the equatorial radius of the central body. Tides raised on the Moon
by the Sun are two orders-of-magnitude smaller than tides raised by the Earth. The acceler-
ation due to constant lunar tides is modeled using a spherical harmonics representation:
Cnm − iSnm = knm2n + 1
∑
j
GMj
GM
Rn+1M
rn+1mj
Pnm(sinϕj )e−imλj (4)
Simplifying, the non-dissipative tides contribute time-varying components to second degree
and order normalized coefficients as follows (McCarthy and Petit 2003):
J 2 = −k20
√
1
5
GMpR
3
GMr3p
[
3
2
sin2 ϕp − 12
]
C21 = k21
√
3
5
GMpR
3
GMr3p
sinϕp cosϕp cosλp
S21 = k21
√
3
5
GMpR
3
GMr3p
sinϕp cosϕp sinλp (5)
C22 = k22
√
3
20
GMpR
3
GMr3p
cos2 ϕp cos 2λp
S22 = k22
√
3
30
GMpR
3
GMr3p
cos2 ϕp sin 2λp
Here, ϕp and λp are the latitude and longitude of the perturbing body on the surface of
the central body. Separate Love numbers have been used for each order, though they are ex-
pected to be equal (k20 = k21 = k22). Degree-3 Love number solutions have been investigated
and their effect is barely detectable.
The tidal potential consists of a variable term and a constant or permanent term. De-
pending on choice of convention, the constant term may or may not be included in the
corresponding gravity coefficient. The MIRAGE-generated gravity fields do not include the
The GRAIL Scientific Measurement System 33
Fig. 3 Possible lunar core motion and the relationship between different frames of reference
permanent part of the tide. Our formulation assumes an elastic Moon and does not include
the frequency-dependent dissipation terms. The elasticity does not affect the overall simula-
tion results and was not considered in this study. The k2 estimate uncertainty from lunar laser
ranging and spacecraft tracking is between 6–8 percent. The GRAIL results will determine
the k2 Love number to better than 1 percent.
The acceleration due to the gravitational potential must be rotated from the body-fixed
principal axis frame to the inertial frame using the lunar physical libration angles included in
a planetary ephemeris database (e.g., JPL DE421) described in Williams et al. (2008). Three
Euler angles describe lunar orientation: the angle along the J2000 equator from the J2000
equinox to the intersection of the lunar equator with the J2000 equator, the angle between
the two equators, and the angle along the lunar equator from the intersection of equators to
the lunar meridian of zero longitude (Newhall and Williams 1997).
On the basis of re-analysis of Apollo seismic observations, Weber et al. (2011) proposed
that the Moon has a solid inner core surrounded by a fluid outer core. Given an oblate inner
core, a time-varying signature could result from the monthly motion of the lunar core equator
relative to the lunar body-fixed or mantle frame (Williams 2007), affecting the degree 2 and
order 1 spherical harmonics and the second-degree tidal potential changes due to the Earth
and Sun. Figure 3 illustrates the Moon’s expected core motion; a point on the core equator
moves relative to the body-fixed equator with a period of one month.
34 S.W. Asmar et al.
Due to the pole offset of the core and mantle frame, the core motion introduces a monthly
signature in the C21 and S21 gravity coefficients as follows:
C21 = α21 cos(ω˙t + ϕ), (6)
S21 = β21 cos(ω˙t + ϕ) (7)
where C21,S21 is the monthly gravitational potential oscillation due to a possible solid
inner core with an axis of rotation tilted relative to the mantle’s axis, included in all simu-
lations, ω˙ is the frequency and ϕ is the phase of this periodic signature. For the latter, we
assume a priori knowledge when estimating the amplitudes of the C21 and S21 signatures
(α21 and β21) along with the gravity field and tidal Love number. If the inner core had an
equilibrium figure for tide and spin distortion, then the ratio of amplitudes for C21 and S21
signatures would be 4. While this ratio is not assumed, it has been used to set requirements
for amplitude uncertainties. We investigated both the uncertainty of the core amplitudes and
the differences of the estimated values with the a priori values. These estimated amplitudes
plus the tidal Love numbers encapsulate the results of GRAIL’s science investigations ad-
dressing the deep interior.
4 Model Estimation and Dynamical Integration
JPL’s gravity field estimation process relies on two primary data types: a link between the
spacecraft and Earth, which is a one-way X-band link, and an inter-spacecraft link called
the Ka-band Range (KBR). The latter’s first derivative, the Ka-band Range Rate (KBRR),
precisely measures the relative movement of Ebb and Flow, which permits estimation of the
lunar gravity field. The combined measurement of two sets of ranging data, one measured
by Ebb and a second by Flow, is called the Dual One-Way Range (DOWR) measurement.
Ebb and Flow are tracked from Earth by the DSN, which produces Doppler data used to
determine the absolute position of each spacecraft:
zd = ρˆse · ρ˙se, (8)
zs = ρˆba · ρ˙ba (9)
where ρse represents the vector from spacecraft to the DSN station and ρba represents the
vector from Ebb to Flow.
The estimation of the gravity field follows the same steps as the orbit determination
process in navigation but involves many more parameters and methods that may constrain
the gravity field and other model parameters to obtain the most realistic solution. Although
the planetary gravity field solutions often require a Kaula power law constraint (Kaula 1966),
the uniform and global coverage of the KBRR data does not require a constraint in our
simulations except for solutions of high degree (i.e., degree ∼270) where a small power-
type constraint was applied.
Letting r and v be the position and velocity vectors of the spacecraft relative to the central
body, the software integrates the second order differential equations
¨r = f (r, v, q) = ∇U(r) + fpm + fin-pm + fin-obl + fsrp + falb + fatt + frel + · · · (10)
Here, f (r, v, q) is the total acceleration of the spacecraft and q are all the constant ( ˙q = 0)
model parameters to be estimated (e.g., gravity harmonic coefficients). Contributions to the
The GRAIL Scientific Measurement System 35
total acceleration include the acceleration of the spacecraft relative to the central body due
to the gravitational potential of the central body ∇U(r), the spacecraft acceleration due to
other solar system bodies treated as point masses fpm, the indirect point mass acceleration
of the central body in the solar system barycentric frame due to the other planets and natural
satellites fin-pm, the indirect oblateness acceleration of the central body (e.g., Moon) due to
another body’s oblateness (e.g., Earth) fin-obl , the acceleration of the spacecraft due to solar
radiation pressure fsrp , the acceleration due to lunar albedo falb , the acceleration due to
spacecraft gas thrusting for attitude control maneuvers (usually for de-spinning angular mo-
mentum wheels) fatt , and the pseudo-acceleration due to general relativity corrections frel .
Other accelerations also exist and may include spacecraft thermal forces, infrared radiation,
tides, and empirical, usually periodic, acceleration models. Specific acceleration models that
have been taken into account are described below.
4.1 Acceleration Due to Solar Radiation Pressure
Each spacecraft is modeled with five single-sided flat plates to model the acceleration due to
solar radiation pressure (SRP) as detailed in Fahnestock et al. (2012) and Park et al. (2012).
For each plate, the acceleration is computed as:
asrp = CSs
msr2sp
(Fnuˆn + Fruˆs), (11)
Fn = −A(2κdvd + 4κsvs cosα) cosα, (12)
Fr = −A(1 − 2κsvs) cosα. (13)
The acceleration due to SRP is on the order of 10−10 km/s2. It is separable from the effect
of gravity in the estimation process. With a ray-tracing technique to model self-shadowing
on the spacecraft bus and on-board telemetry of the power system to detect entry and exit
from lunar shadow, the SRP accelerations can be determined to a few percent level.
4.2 Acceleration Due to Spacecraft Thermal Radiation:
For a flat plate component, the acceleration due to spacecraft thermal re-radiation is:
astr = −2 × 10
−6Aσsb
3msc
εT 4uˆn. (14)
This is used to convert from any given plate’s surface temperature to its acceleration
contribution.
4.3 Acceleration Due to Lunar Albedo and Thermal Emission
The element of acceleration on a spacecraft due to lunar radiation pressure from a point P
on the surface of the Moon can be computed (from Park et al. 2012) as:
dalrp = H(Fnuˆn + Fr rˆps)cosψ
πr2ps
dAplanet . (15)
For reflected sunlight (albedo):
H = CSm cosψs
msr2ms
N∑
=0
∑
m=0
(
CAm cosmλp + SAm sinmλp
)
Pm(sinϕp), (16)
36 S.W. Asmar et al.
and for thermal emission (infrared):
H = C
4msr2ms
N∑
=0
∑
m=0
(
CEm cosmλp + SEm sinmλp
)
Pm(sinϕp). (17)
The albedo map is a constant field whereas the thermal map is a function of local lunar
time because of topographic variation; the thermal map derived using the measurements
from Lunar Reconnaissance Orbiter’s Diviner Lunar Radiometer Experiment data. For this
reason, the following simplified thermal emission model was derived for the simulation of
the total error budget:
H =
⎧
⎨
⎩
CσsbT
4
max cosψs
4msr2csL
, if ψs ≤ 89.5◦,
CσsbT
4
min
4msr2csL
, otherwise,
(18)
where Tmax = 382.86 K and Tmin = 95 K. Thermal maps were computed at the local noon-
time when the Sun is at 0◦ longitude and 0◦ latitude.
4.4 Acceleration Due to Un-modeled Forces
The acceleration due to un-modeled forces is used to represent the errors in the non-
gravitational forces from solar pressure, spacecraft thermal radiation, lunar radiation, and
spacecraft outgassing and is represented as the periodic acceleration formulation:
auf = (Pr + Cr1 cos θ + Cr2 cos 2θ + Sr1 sin θ + Sr2 sin 2θ)eˆr
+ (Pt + Ct1 cos θ + Ct2 cos 2θ + St1 sin θ + St2 sin 2θ)eˆt
+ (Pn + Cn1 cos θ + Cn2 cos 2θ + Sn1 sin θ + Sn2 sin 2θ)eˆn, (19)
where eˆr , eˆt , and eˆn represent the radial, transverse, and normal unit-vectors, respectively
and θ denotes the angle from the ascending node of the spacecraft orbit on the EME2000
plane to the spacecraft. The periodic acceleration is nominally set to zero in the initial tra-
jectory integration and is used to estimate the errors in the non-gravitational accelerations.
The terms Pi represent the constant accelerations during the time interval that the corre-
sponding periodic acceleration model is active. The terms (Ci1, Si1) and (Ci2, Si2) represent
the once-per-orbit and twice-per-orbit acceleration amplitudes, respectively.
In addition to integrating the spacecraft position and velocity, MIRAGE integrates the
variational equations to estimate the epoch state and constant parameters. Following nomen-
clature in Tapley et al. (2004a, 2004b), the nominal trajectory is given by:
X∗(t) =
⎛
⎝
r∗(t)
v∗(t)
q∗
⎞
⎠ . (20)
The first order differential equation to integrate in order to determine the nominal orbit
is:
X˙∗(t) =
⎛
⎝
v∗
f (r∗, v∗, q∗)
0
⎞
⎠ = F (X∗, t). (21)
The GRAIL Scientific Measurement System 37
The variation of the trajectory from its nominal path is x(t) = X(t) − X∗(t) and the
linearized equations:
x˙(t) = A(t)x(t) =
(
∂F (t)
∂X(t)
)∗
x(t) (22)
The integrated solution is the state transition matrix Φ(t, t0), which relates the deviation
from the nominal path at epoch t0 to the deviation from the nominal path at time t for the 6
position and velocity epoch parameters matrix (U6×6) and the p constant model parameters
(V6×6):
x(t) = Φ(t, t0)x(t0) =
[
U6×6 V6×6
0p×6 Ip×p
]
x(t0). (23)
The second order differential equations that MIRAGE integrates for each GRAIL space-
craft include the 3 position variables of Eq. (10), 18 variables representing the changes in
position and velocity due to small changes in epoch position and velocity which define the
matrix U6×6, and three equations for each dynamic parameter or constant from being esti-
mated. For a complete gravity field of degree and order n, the total number of gravity field
parameters is given by (n − 1)(n + 3), or, for example, 32,757 parameters for a 180 degree
and order field.
5 Processing and Filtering of Observations
After numerical integration, MIRAGE processes Doppler and range observations. Following
Tapley et al. (2004a), the general form of the observation equation is
Y = G(X, t) + ε, (24)
where Y is the actual observation, G(X, t) is a mathematical expression to calculate the
modeled observation value, and ε is the observation error. The DSN Doppler data is not an
instantaneous velocity measurement, but is processed in similar fashion to a range observ-
able and is given by a differenced range measurement for two-way Doppler as
G(X, t) = ((r12 + r23)e − (r12 + r23)s
)
/t + · · · (25)
where r12 is the uplink range transmitted by the ground station and received at the space-
craft, and r23 is the downlink range from the spacecraft to the earth station, with subscripts
denoting the end and start of the Doppler count interval, t . To process a Doppler obser-
vation, we must solve the light time equation in a solar system barycentric frame, i.e., find
the original transmit time at the first station and the receive time at the spacecraft using
an iterative procedure. Equation (25) requires DSN calibrations for Earth ionospheric and
tropospheric refraction (Mannucci et al. 1998), and corrections for relativistic propagation
delay due to the Sun and planets, solar plasma delays due to the solar corona of the Sun, and
any measurement biases.
The dual one-way phase measurement between Ebb and Flow can be converted to a
biased range, by an algorithm first developed by Kim (2000). Our lunar gravity recovery
process ingests instantaneous range-rate, modeled as a projection of the velocity difference
vector, r˙12, along the line-of-sight unit vector,

e12.
G(X, t) = ρ˙ = r˙12 • e12 (26)
38 S.W. Asmar et al.
Processing observables also requires the linearized form of Eq. (24). Given an observ-
able Y , we compute a nominal observable Y ∗(t) based on an input nominal orbit, and cal-
culate an observation residual y:
y = Y − Y ∗(t) (27)
Using the state transition matrix to map to the epoch time, Eq. (24) is then written as
y =
(
∂G
∂X
)
Φ(t, t0)x0 + ε = Hx0 + ε. (28)
Based on the vector of residuals y and partials matrix H , the MIRAGE filter solves for a
state X that minimizes these ε error terms.
The calculation of the nominal DSN Doppler observable and related partials in Eq. (28)
involves the precise location of the Earth station in a solar system barycentric ICRF frame
as shown in Yuan et al. (2001). The Earth-fixed coordinate system is consistent with the
International Earth Rotation and Reference Systems Service (IERS) terrestrial reference
frame labeled ITRF93 as shown in Boucher et al. (1994). The rotation of the Earth-fixed
coordinates of the DSN locations to the Earth centered inertial system requires a series of
coordinate transformations due to precession as in the IAU 1976 model described in Lieske
et al. (1977) and nutation of the mean pole as in the IAU 1980 nutation theory described in
Wahr (1981) and Seidelmann (1982) plus daily corrections to the model from the JPL Earth
Orientation Platform (EOP) product of Folkner et al. (1993), rotation of the Earth as in Aoki
et al. (1982) and Aoki and Kinoshita (1983) and UTC-UT1R corrections of the JPL EOP file,
and polar motion of the rotation axis. The JPL EOP product is derived from the Very Long
Baseline Interferometry (VLBI) and Lunar Laser Ranging (LLR) observations and includes
Earth rotation and polar motion calibrations and, in addition, nutation correction parameters
necessary to determine inertial station locations to the level of a few centimeters.
The body-fixed ITRF93 DSN station locations have been determined with VLBI mea-
surements and conventional and GPS surveying. The coordinate uncertainties are about 4 cm
for DSN stations that have participated in regular VLBI experiments, and about 10 cm for
other stations; Folkner (1996) also provides the antenna phase center offset vector for each
DSN station. These DSN station locations are consistent with the NNR-NEWVAL1 plate
motion model (Argus and Gordon 1991). The variations of DSN station coordinates caused
by solid Earth tide, ocean tide loading, and rotational deformation due to polar motion are
corrected according to the IERS standards for 1992 (McCarthy and Petit 2003).
Once the observation equations are found, MIRAGE estimates the spacecraft state and
other parameters using a weighted Square Root Information Filter (SRIF), see Lawson and
Hanson (1995). SRIF computation time dominates MIRAGE processing, and for the larger
planetary gravity fields of the Moon we run on two Beowulf Linux clusters (a 28-node
machine with 112 CPU cores and a 45-node machine with 360 CPU cores). In normal form,
the least-squares solution is given by:
xˆ = (HT WH + P −1ap
)−1
HT Wy (29)
W is the weight matrix for the observations and Pap is the a priori covariance matrix of the
parameters being estimated. In the MIRAGE SRIF filter, the solution equation is kept in the
form:
Rxˆ = z (30)
The GRAIL Scientific Measurement System 39
R is the upper triangular square-root of the information array and R and z are related to the
normal equations as:
RT R = HT WH + P −1ap , (31)
z = (RT )−1HT Wy (32)
and the covariance P of the solution (inverse of the information array) is given by:
P = R−1(R−1)T (33)
We separate observations for gravity field determination into disjoint time spans called
data arcs. Two-day-long data arcs are typical. The parameters estimated in the arc-by-arc
gravity solutions consist of arc-dependent local variables: spacecraft state, solar radiation
pressure coefficients, etc., and global variables common to all data arcs: gravity coefficients,
tide parameters, etc. Merging the global parameter portion of a sequence of data arc square
root information arrays produces a solution equivalent to solving for a single set of global
parameters plus independent arc-specific local parameters (Kaula 1966).
When solving for a large number of parameters, convergence is very sensitive to a priori
values and uncertainties. If the spacecraft initial state is poorly known and a filter tries to
solve for both the trajectory and a high-resolution gravity field at the same time, the iteration
may never converge. In order to avoid this problem, the local parameters are first estimated,
and once a solution is obtained, the global parameters are estimated.
For each spacecraft, the local parameters consist of the spacecraft initial state, the solar
radiation pressure scale factor, two constant SRP scaling terms orthogonal to the spacecraft-
to-Sun vector, fifteen periodic acceleration terms for every two hours, four inter-satellite
range-rate measurement correction terms for every two hours, and constant Earth-based
Doppler bias and drift rate. Local parameters are used to constrain non-gravitational effects
and measurement biases and are chosen based on experience. The global parameters consist
of three inter-satellite range-rate time-tag biases, degree 2 and 3 Love numbers, degree 2 and
order 1 amplitudes of periodic tidal signature, Moon’s mass (GM), and a 150 × 150 gravity
field (approximately 23,000 parameters). The time-tag biases represent the offset between
the DSN time and a KBRR time-tag derived from the spacecraft clock.
Due to the accumulation of spacecraft angular momentum, maneuvers for Angular Mo-
mentum Desaturations (AMD) take place periodically. AMD maneuvers disrupt the quiet
environment for gravity measurement and break the arc of data to be processed. Since we
expect maneuvers, and to avoid numerical noise limitations on trajectory integration, we
postulate 2-day arcs in our simulations. As described in Park et al. (2012), for each 2-day
arc, we first estimate and re-estimate local parameters for each arc until convergence. Hav-
ing converged on local parameters, we then compute SRIF arrays containing both local and
global parameters for each arc, combine, and estimate, re-compute, re-combine, re-estimate,
repeating until convergence.
6 Modeling Parameters
The input parameters to the simulations of the GRAIL mission are discussed below, grouped
in the categories of data noise, data coverage, data arcs, orbital parameters, dynamic errors,
and kinematic errors. To show the types of issues the simulation team was addressing, Ta-
ble 2 lists a summary of parameters relevant to the simulation results and our model confi-
dence in each one.
40 S.W. Asmar et al.
Table 2 Confidence level in GRAIL parameters relevant to science simulations
Parameter Assumption Note Confidence
level
Orbit initial
conditions
Orbit conditions and spacecraft
alignment are favorable
Inclination and node differences
between spacecraft match
requirements (0.02◦)
High
Instrument noise
spectrum
Spectrum includes thermal noise
and USO jitter
GRACE analysis and performance
modified for GRAIL
High
DSN data amount
and noise
Tracking coverage is sufficient and
noise characterization valid,
includes USO
8 hours per day per spacecraft.
DSN noise of 0.05 mm/s at 10-s
integration time
High
KBR data continuity No hardware resets Tested with 5-min gaps once per
day; show no impact
High
Time tag offset
between payload and
DSN time
Known to 100 ms, stable to 100
micro-seconds over 2 days
Convergence confirmed in science
simulation
Medium to
high
Temperature of
spacecraft and
payload elements
Linear dependence on beta angle Tested conservatively, small error
contribution
High
Propellant leakage Constant and small Preliminary information from
spacecraft team
Medium
Outgassing Small after cruise Data from previous spacecraft High
Lunar surface
radiation
From lunar mission experience Published models Medium
Fuel slosh Very small Use of a propellant tank diaphragm High
Solar radiation
pressure
Constant reflectivity properties per
arc and un-modeled errors <2 %
Currently investigating variability
over an arc
Medium to
high
Lunar librations Modeled with Lunar Laser
Ranging data
Known to a few milliseconds of arc High
Lunar core signature Monthly periodic Phase not known Medium
6.1 Data Noise
GRAIL simulations tools create DSN Doppler and inter-spacecraft Ka-band range rate data
and apply noise to both data types. The simulations do not include DSN range data since
this data type does not significantly improve GRAIL orbital accuracy. Since data noise levels
are non-Gaussian, the applied noise and data weights assigned during follow-up parameter
estimation are not always identical.
As discussed earlier, only the X-band Doppler link was included for simulation purposes,
not the communications and navigation two-way S-band link between the DSN and the
GRAIL spacecraft. The shorter wavelength X-band is less susceptible to ionosphere and
interplanetary plasma noise. Expressed in units of velocity, our studies assume 0.05 mm/s
DSN one-way X-band link residual noise and data weight at an integration time of 10 s,
when simulating Doppler data and when filtering simulated data. This assumption is slightly
more conservative than the typical noise level of 0.03 mm/s at 10 s experienced with the
Mars Global Surveyor X-band performance.
The non-Gaussian residual noise associated with the payload’s KBRR is added to the
simulated inter-spacecraft data as a function of frequency. The long wavelength noise for
5-s samples is 0.4 µm/s for long wavelengths and then transitions to 1.0 µm/s at the short
wavelength. However, in the filter, a constant data weight of 1 µm/s white noise is applied.
The GRAIL Scientific Measurement System 41
6.2 Data Coverage and Data Arcs
As a baseline, we simulated 8 hours of DSN daily tracking data for each spacecraft, non-
overlapping, for a total of 16 hours. This DSN coverage provides information on absolute
orbit for Ebb or Flow and improves long wavelength gravity field solutions, including the
lunar core parameters, but contributes minimally to the global and regional science require-
ments. The coverage of the KBRR data is assumed to be continuous. Obtaining 16 hours of
DSN coverage per day, every day, for one mission is considered very challenging due to the
loading on the DSN but the requirements were accepted since the prime mission duration is
relatively short, on the order of 3 months.
Since successive momentum dumps occur typically two days apart, GRAIL simulations
assume a two-day data arc length, starting from an epoch of 4 March 2012 (actual epoch var-
ied). Longer arcs are typically desirable but the momentum dumps are their natural bound-
aries.
6.3 Orbital Parameters
During the 82-day Science Phase, the Moon rotates three times underneath the GRAIL orbit.
The collection of gravity data over one complete rotation, 27.3 days, is called one mapping
cycle. Ebb and Flow are in a common near-polar, near-circular orbit with a mean altitude of
approximately 55 km during the prime mission. However, as described in Roncoli and Fujii
(2010) the periapsis altitude ranges from approximately 16 km to 51 km above a reference
lunar sphere. The Ebb-Flow separation distance is designed to slowly vary. For approxi-
mately the first half of the mission, they drift apart and their separation distance increases
from ∼85 to ∼225 km and then, with only one small orbit trim maneuver, they drift to-
wards each other and the distance decreases to ∼65 km near the end of the mission. The
shorter separation distance is optimum for data exploring the local and regional spatial fea-
tures while the segment around the maximum separation is optimum for the determination
of the global studies such as the lunar core parameters, which are the Love number and the
periodic signature of degree 2. The separation distances are designed to ensure that there is
no degradation of the Ka-band signal due to multipath off the lunar surface and, according
to Roncoli and Fujii (2010), a shorter spacecraft separation is required because of the lower
spacecraft altitude.
The spacecraft separation contains a drift in order to reduce the resonance effects corre-
sponding to the harmonic of the separation distance. Resonance effect degradation occurs
at harmonic degrees of the form N = 360/(D/30) where D is the spacecraft separation.
This corresponds to degrees 54, 108, and 162 for a 200-km separation distance. For a 50-km
separation, the resonance occurs much later, starting at degree 216.
The spacecraft inclination varies between approximately 88.4 and 89.85 degrees with a
twice-per-month periodic signature. The average inclination, approximately 89.1 degrees,
is offset from a perfectly polar orbit to improve the determination of low degree harmonic
coefficients, but kept to a minimum to reduce the gap in data coverage at the poles.
The GRAIL science orbital phase is limited in part by a solar beta angle constraint of
49 degree imposed by the capability of the electrical power system; the spacecraft cannot
generate sufficient power from the solar arrays for angles below this constraint. Figure 14
of Roncoli and Fujii (2010) illustrates the time history of the solar beta angle as well as the
relationship between beta angle and the duration of solar eclipse during the science phase;
eclipse durations are a maximum at the beginning and end and no solar eclipses when the
solar beta angle is near 90 degrees near the middle. For the simulations the Sun-angle is
42 S.W. Asmar et al.
represented after Park et al. (2012) with the inclination computed with respect to the lunar
pole vector, as:
β = 90◦ − cos−1(eˆh · rˆsm) (34)
To separate the non-gravitational signature from gravitational effects, a β angle of 90◦ is
optimal. The spacecraft enters terminator crossing at β ∼ 76◦, and for β-angles less than
this, modeling non-gravitational forces becomes more difficult, as the perturbations change
rapidly due to partial shadowing of the orbit.
6.4 Dynamic Errors
The MIRAGE filter estimates both local parameters that are dependent for each data arc
and global parameters that are common to all data arcs. The local dynamic parameters that
are estimated include three dimensionless parameters for the solar pressure model of each
spacecraft, a constant scale factor for the force along the sun-spacecraft direction with a
nominal value of 1.0 and an a priori uncertainty of 0.10 or 10 % of the solar pressure force,
and a scale factor for each off normal directions with a nominal value of 0.0 and an a priori
uncertainty of 0.02 or 2 % of the overall solar pressure force. The solar pressure is modeled
as a box-wing plate model with appropriate specular and diffuse coefficient values for each
plate. Since both spacecraft are nearly identical, the solutions for the solar pressure scale
factors are expected to be nearly equal. Several gravity solutions were also generated with
strong a priori correlations between the two spacecraft solar pressure solutions to force them
to be nearly equal. This constraint improved the core parameter uncertainties by about 30 %
but had little effect on the global and regional gravity requirements. The baseline approach,
however, is to treat the solar pressure solution of each spacecraft independently.
To account for un-modeled residual solar pressure errors, 15 periodic coefficients are
estimated for each spacecraft for each arc with a priori amplitude equal to ∼2 % of the solar
pressure force, or 3 × 10−12 km/s2. The coefficients include constant, once per revolution
and twice per revolution amplitudes for the radial, orbit normal, and along the velocity
directions.
The lunar albedo model is not part of the baseline simulations results but albedo errors
were independently investigated and found to be minimal. The albedo surface representation
is given by the 10th degree spherical harmonic expansion of Floberhagen et al. (1999).
The lunar surface thermal re-radiation is also investigated and is similar in size to albedo.
Another non-gravitational force to be considered for GRAIL is the thermal radiation force
as a result of heating on the spacecraft; this force is assumed to be small.
In the estimation of the gravity field, we assume a nominal gravity field and a truth gravity
field. For the early simulation, the lunar gravity model derived from the Lunar Prospector
mission to degree and order 150 and designated LP150Q in Konopliv et al. (2001) was used
as both the nominal and truth models. Simulations since then have also used the smaller
LP100J lunar gravity model to degree 100 as the nominal model to test convergence to the
LP150Q truth model for different modeling assumptions.
The global dynamic parameters that are estimated include a gravity field to a given de-
gree, the second-degree Love numbers, and the periodic amplitudes of C21 and S21 for core
detection. In order to reduce computation time, the globally estimated gravity field was to
degree 150 and extrapolated the results to degree 180. With current assumptions about data
quality, it is expected that the gravity field will be recovered to higher than degree and order
300, which would make the Moon the body with the highest known gravity resolution in the
universe.
The GRAIL Scientific Measurement System 43
6.5 Kinematic Errors
In addition to dynamic errors, which directly affect the spacecraft orbit because of a force,
there are kinematic errors that affect the KBRR or DSN Doppler observables directly. Most
of the kinematic corrections are estimated as local parameters that affect only that data arc.
For every GRAIL orbit period of ∼2 hours, a KBRR bias, drift, and cosine and sine once
per revolution amplitudes are estimated. For a 2-day arc this amounts to 48 parameters
being estimated. No a priori constraint is applied to these parameters. For each arc, one
DSN Doppler bias and one drift parameter are estimated to correct the USO frequency for
the X-band one-way Doppler.
The very important kinematic error of timing offsets is addressed in a separate section
below. Other kinematic errors are investigated by introducing noise or systematic trends in
the KBRR residuals. The KBRR noise spectrum includes error contributions from the USO
and spacecraft attitude jitter and the simulation account for them with a spectrum ranging
from 0.4 µm/s at the long-wave-length to 1 µm/s at the short wavelength. Errors due to the
offset of the phase center from the line connecting the center-of-mass from each spacecraft
are minimized by actively pointing the spacecraft to align the phase centers. The changes in
the temperature of the payload antenna and related hardware are modeled as a systematic
trend in the observables. The current models we have investigated are sinusoidal once per
orbit tones of 80 µm, a twice per orbit tone, and a triangular shaped twice per orbit signature.
Amplitudes of all of these depend on the Sun angle. There is also a small error due to a shift
reactive to the Ka-band system. These errors have a small effect on the core parameters and
a negligible effect on the global and regional requirements. Errors due to sloshing of the fuel
in the spacecraft fuel tank are expected to be negligible.
7 Modeling System Contributions
GRAIL utilized an integrated scientific measurement system comprised of flight, ground,
mission, and data system elements in order to meet the end-to-end performance required
for achieving the scientific objectives. Simulations leading to end-to-end error budgets were
used to optimize the design, implementation, and testing of these elements. Because the
inter-spacecraft range and range-rate observables and the range-rate between each spacecraft
and ground stations can be affected by the performance of all elements of the mission, they
were all treated as parts of an extended science instrument or a science system.
7.1 Flight System
The flight system is the spacecraft, which is comprised of the orbiters and science payload.
The orbiter is comprised of the spacecraft bus, heritage from the Experimental Satellite
System 11 (XSS-11) mission, and subsystems described in Hoffman (2009). At the heart of
the flight system is the LGRS payload. Its design and performance are detailed in Klipstein
et al. (2013, this issue). Key components of the GRAIL payload, such as the USO and RSB,
are discussed here in the context of the quality of radio links, Doppler data, timing effects,
and end-to-end error budgets.
As discussed above, to assess the ability of the flight system to meet the science require-
ments, all known sources of dynamic and kinematic error were assessed for their possible
contributions to the uncertainty of simulated solutions for the gravity field and time vary-
ing tidal and core parameters. Fahnestock et al. (2012) described how early calculations
44 S.W. Asmar et al.
of these effects revealed that the acceleration in reaction to spacecraft thermal re-radiation
was among the dominating sources of solution uncertainty. Thus detailed modeling of the
spacecraft’s geometry, surface material, optical properties and thermal state was investigated
along with approaches to the interface of such modeling with the science data processing.
First, a model of each spacecraft was constructed with plate surfaces defined for each
spacecraft to cover the external area, not accounting for protuberances such as the thrusters
and cameras. The two spacecraft are virtually identical but externally mirrored about the
body-frame X–Z plane (see Fig. 4). Thermal modeling was performed by Lockheed Martin
using the Thermal Desktop software, and for each spacecraft for the thermal loading con-
ditions of varying beta angles. This included incident sunlight, lunar albedo and infra-red
(IR) emission, and equipment power modes during the data acquisition period. Steady-state
nodal temperature output, repeatable from orbit to orbit, was averaged over the plate sur-
faces to produce averaged temperature profiles for each plate surface over one orbit. These
were used with the best available material optical properties, most measured after thermal
vacuum bake-out, to compute body frame acceleration profiles over one orbit. These were
then double interpolated over phase angle to obtain body frame acceleration at any epoch in
the science mission.
As detailed in Fahnestock et al. (2012), the sensitivities of the thermal re-radiation accel-
eration history model with respect to assumptions and inputs were examined for impact on
science. These included quantifying the difference spatial averaging of nodal temperature
output versus no spatial averaging over the solar array surfaces, the sensitivity to extract-
ing power from the solar arrays, the sensitivity to 3σ variations in the optical properties of
every type of used material, and sensitivity to additive worst-case temperature biases on all
surfaces.
A 5 K global bias bounded most changes resulting from realistic input variations. The
difference in predicted thermal re-radiation acceleration history between that case and the
nominal case was taken as indicative of the magnitude of un-modeled accelerations that
would act on the spacecraft, which was used to create an a priori error model for the constant,
and once and twice per orbit sinusoidal, periodic acceleration parameters in the simulations.
For the gravity field solution, the modeled thermal re-radiation acceleration history,
a telemetry-derived thermal re-radiation acceleration history, or a hybrid of the two was
included in the nominal dynamical model. For the first, the solar array bus open circuit volt-
age and short circuit current telemetry channels were utilized to determine the actual epoch
of transit into and out of the Moon’s shadow, and then the modeled history was shifted and
stretched in time to match these transit times. For the telemetry-derived history, given n
surfaces in the spacecraft model, each of m < n instantaneous surface temperatures was set
equal to the average of instantaneous readings of one or more temperature sensors, selected
based on geometric proximity to that surface. Each of the n − m remaining surface tem-
peratures was set equal to the closest associated one of the m surface temperatures, plus a
bias computed as the time average, over one orbit at β = 90◦ epoch, of the difference of the
pre-flight thermal modeling output for the two surfaces in question. The β = 90◦ epoch had
no eclipsing and was when the spacecraft were the most thermally and dynamically quiet,
so we chose to tailor our mapping from sensor temperatures to surface temperatures to this
time period. The hybrid history, illustrated in Fig. 5, was a combination of the modeled his-
tory for the Y and Z body axes directions, and the telemetry-derived history for body X
axis direction. The apparent best methodology to use in data processing appears to be the
hybrid thermal re-radiation acceleration history and a priori error model derived from it, as
described earlier.
The GRAIL Scientific Measurement System 45
Fig. 4 Three views illustrating the placement of surfaces (numbered in boxes) on Flow’s highest fidelity
thermo-optical model, used for thermal re-radiation acceleration calculation. Ebb is similar, and temperature
sensor locations are also noted
46 S.W. Asmar et al.
Fig. 5 Spacecraft body-frame components of hybrid thermal re-radiation acceleration history
The SRP acceleration was also included in the dynamical modeling in MIRAGE along
with lunar albedo and IR emission pressure accelerations. For computing all three, a five-
plate model of each spacecraft was obtained from the highest resolution model. This was
done through combination of surfaces, first for common material and orientation, then for
only common orientation, but always in a manner that kept the models roughly equivalent
in four quantities: power coming into the spacecraft, acceleration due to SRP, power com-
ing out of the spacecraft, and acceleration due to thermal re-radiation. All modeling of the
non-gravitational accelerations within the spacecraft dynamics was sufficiently accurate to
improve the orbit determination and data calibration as well as the gravity field solution.
7.2 Ground System
The ground system is comprised of the stations of the DSN, which was used exclusively by
GRAIL, and the infrastructure facilities at JPL that support the DSN operations and transfer
of data to users. The DSN’s primary service to GRAIL and other missions is the telecommu-
nications aspect, telemetry and commanding as well as navigation radio-metric data types.
DSN services relevant to science are those providing precision measurements of the sig-
nal carrier frequency/phase for the purpose of Doppler observables. Two DSN subsystems
Tracking and Radio Science are examined. As described earlier, GRAIL utilizes two bands
for links to the DSN, a two-way S-band for telemetry, command, and navigation and a one-
way X-band un-modulated carrier from the payload’s RSB for science. As a result, DSN
stations capable of S-band uplink and downlink as well as X-band downlink are required to
support the mission, a criterion that narrows the available stations to a few 34-m diameter
stations throughout the network but still enables the mission to get required coverage.
The Tracking Subsystem The Doppler data observable is generated in real-time at the DSN
stations. The tracking receiver is a closed-loop system that finds the carrier frequency via a
built-in algorithm and tracks it, aided by a prediction file for initial acquisition, producing
the receiver’s one-time real-time computation of the frequency and Automatic Gain Control
The GRAIL Scientific Measurement System 47
(AGC). Within its design threshold for dynamic conditions and signal-to-noise ratio, the
output of the receiver is useful with a quantifiable Doppler noise that ranges between 0.02
and 0.1 mm/s. For GRAIL simulations, the X-band, USO driven link is associated with
Doppler noise of 0.05 mm/s. If the threshold is exceeded, the receiver loses lock and data
are not recoverable.
Radio Science Receiver Designed specifically for Radio Science experiments described in
an overview by Asmar (2010), the DSN’s Radio Science Receiver (RSR) is at the heart of an
open-loop reception/recoding subsystem that preserves the raw qualities of the electromag-
netic wave propagating from the spacecraft source to the DSN stations. The digital receiver
neither locks onto the carrier signal nor makes real-time decisions about its frequency or
amplitude. Instead, it down-converts the signal in a predictions-driven heterodyne method
and records the raw complex samples into files for users’ post-pass processing. Asmar et al.
(2005) describes the RSR usage and typical performance derived from other missions.
The use of the RSR proved to be critical for enhancing the quantity and quality of GRAIL
X-band Doppler data in two ways: (1) practically double the amount of RSB data received
by the DSN by an unofficial use of the concept of multiple-spacecraft per aperture, where
the DSN station scheduled to track Ebb, for example, also views Flow and vice versa, and
(2) contribute to understanding the various timing effects, as explained at length in Sect. 8.
GRAIL funded the development of 3 portable RSRs for use throughout the DSN in support
of GRAIL data acquisition and timing synchronization.
7.3 Mission System
The mission operations system is comprised of the JPL and DSN infrastructure as well as
the Lockheed Martin operations. Since the mission design affects science data quality, the
following factors had to be very carefully considered: spacecraft altitude and spatial varia-
tion of the altitude, spacecraft separation distance, orbit inclination, ground track separation,
mission duration, number of maneuvers, time separation of maneuvers, AMD separation,
and amount of ground station coverage. Mission design, navigation, deep space stations,
and the ground data system critically contribute to the quality of science data. Specifics of
each factor were described above and additional details on the mission system can be found
in Roncoli and Fujii (2010), Hatch et al. (2010), Beerer and Havens (2012), and Zuber et al.
(2013, this issue).
7.4 Data System
The GRAIL Science Data System (SDS) is comprised of all project hardware and software
tools that contribute to the quality of the science data. Sine the SDS team is the first team to
assess the quality of the data on a daily basis, it provides immediate feedback to the Mission
System on the health of the spacecraft, payload, or ground system in case action is required
to address an anomaly.
Zuber et al. (2013, this issue) provides a functional block diagram of the SDS that shows
the data flow from all the sources to the final science users and the archives of the Planetary
Data System. For science data processing, delivery, and archiving, the SDS is organized
to provide daily support including weekends in order to handle the data volume as well as
prevent any oversight of anomalies for any extended period of time.
48 S.W. Asmar et al.
Fig. 6 Timing effects in the science data system
8 Data Timing and Synchronization
We have described how the process of determining the lunar gravity field starts with the
inter-spacecraft Ka-band phase measurements used to compute the DOWR observables that
are in turn converted to instantaneous range, range rate and range acceleration measure-
ments. Very accurate timing of the measurements is crucial to achieve the high accuracy
gravity results. The timing of the measurement has two components: absolute timing knowl-
edge allows us to assign a measurement to a position around the Moon, and the relative tim-
ing knowledge between the payload clocks needed for the proper alignment of the LGRS
phase measurements of both spacecraft at a common coordinate time. Following the formu-
lation used in Kim (2000) and Thomas (1999), our analysis shows the aggregate errors in
range to be below 1 µm and in range-rate to below 1 µm/s.
The GRAIL data are time-stamped on each spacecraft by the payload with time derived
from the payload clock, namely the USO; Ebb and Flow each carries one USO that drifts
independently from the other. The data are then passed to the spacecraft’s central computer
for packetizing prior to transmission to Earth as telemetry, and the computer puts a time-
stamp on the packets derived from the spacecraft clock, which is independent from the
USO. Finally, the data packets are received by the DSN and time-stamped at arrival, one-
way light time after transmission, with the DSN time, which is derived from yet another
independent clock. Counting two USOs, there are four independent clocks to synchronize
in the post-processing in order to prevent errors in the gravity field solution, and this is
carried out by the SDS team in the Level 1 processing stage as outlined in Fig. 6. The SDS
team estimates the necessary time tag correction by combining information from available
sources: the telemetry packets received at the DSN, the absolute frequency observed at the
RSR, and the synchronization from the onboard inter-spacecraft TTS (note that additional
observations were obtained at the DSN’ RSR after launch of the inter-spacecraft TTS signals
The GRAIL Scientific Measurement System 49
edge to deduce the absolute time on-board and correlate it to the known USO drift resulting
in micro-second level of accuracy, orders of magnitude better than time correlation packet
information).
For the purpose of the simulations, it is assumed that the a priori clock offset knowledge
is 100 milliseconds, constant over one month, and that the reconstruction of the time tag
offset is ∼20 microseconds. To achieve such accuracies, the RSB was added to the science
payload to transmit a one-way X-band un-modulated sine-wave signal generated by the USO
and recorded by the DSN’s RSR. The RSR measurement of the frequency bias is <10−5 Hz
and standard deviation = 10−3 Hz.
The DSN clocks are synchronized with the highly stable Coordinated Universal Time
(UTC) standard. GRAIL’s data processing, on the other hand, utilizes Barycentric Dynam-
ical Time (TDB). The timing analysis derives a time correlation between the LGRS clocks
and the TDB time scale, to be provided as a Level 1 ancillary data product called CLK1B.
To produce CLK1B relating LGRS and TDB, we preprocess the timing data types and run
them through a non-causal Kalman filter.
The LGRS clock on each spacecraft is driven by the USO for maximum stability, which
is 3 × 10−13 over integration times of 1 to 100 seconds, expressed in Allan Deviation. This
clock, however, does not report the absolute time but reports readout with respect to the
clock startup epoch with errors from the drift of the USO. Relying on the nearly quadratic
behavior of the onboard clock and an assessment of relativistic contributions, we believe
that this system enables determining the relative time on Ebb vs. Flow with a bias < 10−7 s
and standard deviation = 9 × 10−11 s.
The Onboard Spacecraft Clock (OSC) is derived from a crystal oscillator with inferior
stability to the USO. The onboard computer tags LGRS timing data packets with OSC time,
including the LGRS 1 Pulse Per Second (PPS) packet. Ebb and Flow transmit time corre-
lation packets to DSN stations where the arrival time is recorded in UTC, which provides
a time correlation between the OSC and UTC. The DSN uses very stable hydrogen maser
clocks and time-stamps the arrival of telemetry and tracking data in the UTC frame, which
is tied to the International Atomic Time (TAI) frame. Based on DSN monitoring reports,
the real-time timing performance of DSN time-tags is at the microsecond level and post-
processing analysis improves the performance to the 10−9 second level.
By combining the LGRS/OSC and OSC/UTC time correlation products, a time correla-
tion between LGRS time and UTC can be determined and the OSC clock drops out. Because
OSC error is under one microsecond over intervals shorter than one second, the stability
characteristics of the OSC do not limit LGRS and UTC correlation accuracy. Considering
possible unknown timing delays in packet transmission, we expect a measurement bias of
up to 100 milliseconds, and standard deviation of up to 30 milliseconds.
9 Relativistic Effects
Turyshev et al. (2013) has developed a realization of astronomical relativistic reference
frames in the solar system and its application to the GRAIL mission. A model was devel-
oped for the necessary space-time coordinate transformations for light time computations
addressed practical aspects of the implementation and all relevant relativistic coordinate
transformations needed to describe the motion of the GRAIL spacecraft and to compute ob-
servable quantities. Relativistic effects contributing to the double one-way range observable,
which is derived from one-way signal travel times between the two GRAIL spacecraft were
accounted for and a general relativistic model for this fundamental observable of GRAIL,
50 S.W. Asmar et al.
accurate to 1 µm and range-rate to 1 µm/s were also developed. The formulation justifies
the basic assumptions behind the design of the GRAIL mission and may also be used in
post-processing to further improve the results after the mission is complete.
It was recognized early during GRAIL’s development phase that due to the expected high
accuracy of ranging data, models of its observables must be formulated within the frame-
work of Einstein’s general theory of relativity in order to avoid significant model discrep-
ancy. The ultimate observable model must correctly describe all the timing events occurring
during the science operations of the mission for the links to Earth as well as the inter-
spacecraft links. The model must take into account the different times at which the events
have to be computed, involving the time of transmission of the Ka-band signal at one of the
spacecraft, say Ebb, at the reception of this signal by its twin, Flow. In addition, the model
must include a description of the process of transmitting S-band and X-band signals from
both spacecraft and reception of this signal at a DSN tracking station.
Relevant points regarding relativistic corrections at the level of accuracy required by
GRAIL include: (1) for a spacecraft around the Moon, we can model proper time treating
the Moon as a point mass; (2) JPL’s long-standing ODP models designed for proper time
of a station on Earth are already sufficiently accurate, with no changes required; (3) up to
a constant bias, computing one-way range from DOWR requires a pair of corrections from
one-way light time to instantaneous distance. It suffices to iteratively solve for light time in
terms of instantaneous distance, by re-computing transmission position bearing in mind the
elapsed light time, in the presence of the Shapiro delay.
Turyshev et al. (2013) also notes that measuring the signal frequency involves computing
three numbers: the derivative of proper time at the receiver with respect to coordinate time
of reception, the derivative of proper time at the transmitter with respect to coordinate time
of transmission, and the derivative of coordinate time of reception with respect to coordinate
time of transmission. The first number must be modified to account for the fact that the clock
at the DSN receiver attempts to synchronize with UTC time, rather than simply acting as a
TDB receiver placed on the surface of the earth. The effect of the Earth’s and the Moon’s
gravity on the third term will be below our level of error; if we did choose to include them
it would certainly suffice to use a point mass.
10 Results of Simulations
10.1 A Priori Assumptions and Kaula Constraints
The a priori uncertainties for the models used in the simulations are summarized in Table 3.
Furthermore, we assume that the KBRR data have σs = 1 µm/s uncertainty at 5-second count
time, and the DSN Doppler tracking uncertainty at 10-second count time σd = 0.05 mm/s,
obtained from previous payload and flight experience. The KBRR data give an average
accuracy of 2×10−10 km/s2, which is equivalent to 0.002 mGal in the gravity measurement.
Gravity field estimates often require a constraint by the Kaula rule (Kaula 1966). In this
study, we assume the Kaula constraint to be 2.5 × 10−4/n2. Note that the Kaula constraint
approximates the root-mean-square (RMS) of the spherical harmonics coefficients, i.e.,
RMSn =
√
σ 2n
2n + 1 , (35)
for the degree variance σ 2n =
∑n
m=0(C
2
nm + S2nm).
The GRAIL Scientific Measurement System 51
Table 3 GRAIL simulations a
priori uncertainty model Parameters A priori uncertainty
Position 1000 km in each direction
Velocity 1 m/s in each direction
Overall SRP scaling factor 10 %
Orthogonal SRP scaling factors 2 % in each direction
Doppler correction terms open
LGRS correction terms open
Periodic acceleration Variable
Lunar gravity field LP150Q solution
k2, α21, and β21 open
Table 4 Baseline core-signature estimates
k2 × 104 σk2 × 104 α21 × 1010 σα21 × 1010 β21 × 1011 σβ21 × 1011
Requirement 2.00 1.00 2.50
Baseline 0.93 0.94 0.62 0.67 −2.01 1.77
10.2 Dynamic Conditions and Core Signature Model Results
Table 4 shows the estimates of the core signature terms for the baseline case, where i rep-
resents the difference between the truth and the estimated parameter i and σi represents the
estimated 1-σ uncertainty of parameter i. The effect of different data arc lengths, which are
bounded by spacecraft maneuvers with a baseline of 2 days, is that the longer the arc the
better the results for the core signature, as expected (Park et al. 2012). The baseline assumes
Ebb and Flow are tracked by the DSN for at least 8 hours per spacecraft per day. Longer
tracking is better and the diversified coverage over different DSN complexes provides addi-
tional orbit information due to geometric parallax.
The remaining possible error contributions from un-modeled non-gravitational accelera-
tions are applied as an a priori periodic acceleration model that impacts the formal uncertain-
ties of the estimated low-degree gravity field and time-varying core signature. The minimum
periodic acceleration is chosen to be 3 × 10−13 km/s2, as shown in Fig. 7, and scaled up ac-
cordingly for the period of lower spacecraft separation; with this minimum a priori periodic
acceleration, all science requirements are met. We showed in Park et al. (2012) that varying
the a priori periodic acceleration to the level of 1 × 10−12 km/s2 still allows satisfying all
science requirements but with smaller error margin for the core parameters.
10.3 Kinematic Error Results
Table 5 shows the effect of kinematics errors on the estimated core parameters. The error
in the estimated parameters represents the effect due to kinematic errors and the formal
uncertainties are the same as in the baseline case. The temperature control case shows the
effect of the error in the LGRS measurement frequency and signal path-length due to thermal
variation. The time-tag error shows the contribution to the total error of estimating the KBRR
time tag with an initial time-tag offset of 100 milliseconds. Lastly, the attitude pointing error
shows the effect of a 3-σ single-spacecraft attitude pointing error, which translates to about
0.06 µm/s on LGRS data.
52 S.W. Asmar et al.
(A)
(B)
Fig. 7 A priori periodic acceleration with 3×10−13 km/s2 minimum acceleration constant term in panel (A)
and once-per-orbit and twice-per-orbit terms in panel (B) (Park et al. 2012)
The GRAIL Scientific Measurement System 53
Table 5 Effect of kinematics
errors on estimated core
parameters
Cases k2 × 104 α21 × 1010
Temperature control 0.05 0.01
Time-tag error 0.02 0.02
Attitude pointing error 0.20 0.10
Fig. 8 Root mean square of the baseline simulation results
10.4 Summary
The most significant result of the simulations is the RMS of the estimated gravity field
corresponding 1-σ formal uncertainties and the difference between the truth and estimated
gravity fields, shown in Fig. 8. Also shown in the same figure are the global and regional
science requirement lines generated based on the surface acceleration accuracy required for
global and regional science requirements, which are satisfied with the baseline assumptions.
The difference between the truth and estimated gravity fields is smooth and is bounded
by the formal uncertainty indicating a correct filter setup and a stable filter solution. The
colored measurement noise is well bounded by the white noise assumption that was used
in the estimation process since the recovered values are well represented by the formal
uncertainties. The linear extrapolation of the estimated uncertainties indicates that a nominal
gravity field of degree 300 or better can be determined according to the Kaula rule. The
largest source of non-gravitational error comes from spacecraft thermal radiation, which is
characterized with variable a priori error constraint model derived from orbit geometry and
expected force magnitude.
54 S.W. Asmar et al.
With all error models included, detailed and numerous simulations show that estimating
the lunar gravity field is robust against dynamic and kinematic errors and meets the high
accuracy lunar gravity requirements by at least an order of magnitude. A nominal lunar
gravity field of degree 300 or better can be achieved according to the scaled Kaula rule
for the Moon. The core signature is more sensitive to modeling errors and depends on how
accurately the spacecraft dynamics can be modeled; the requirement can be achieved with a
small margin.
Acknowledgements The GRAIL mission is supported by the NASA Discovery Program under contracts
to the Massachusetts Institute of Technology and the Jet Propulsion Laboratory. The work described in this
paper was mostly carried out at Jet Propulsion Laboratory, California Institute of Technology, under contract
with the National Aeronautics and Space Administration. The authors thank colleagues who have contributed
to this work or reviewed it, especially at JPL: Duncan McPherson, Ralph Roncoli, William Folkner, Kevin
Barltrop, Charles Dunn, William Klipstein, Randy Dodge, William Bertch, Daniel Klein, Dong Shin, Stefan
Esterhausin, Slava Turyshev, Tom Hoffman, Charles Bell, Hoppy Price, Neil Dahya, Joseph Beerer, Glen
Havens, Robert Gounley, Ruth Fragoso, Susan Kurtik, Behzad Raofi, and Dolan Highsmith. From Lockheed
Martin Space Systems Company (Denver): Stu Spath, Tim Linn, Ryan Olds, Dave Eckart, and Brad Haack,
Kevin Johnson, Carey Parish, Chris May, Rob Chambers, Kristian Waldorff, Josh Wood, Piet Kallemeyn,
Angus McMechan, Cavan Cuddy, and Steve Odiorne. From the NASA Goddard Space Flight Center: Frank
Lemoine and David Rowlands, and from the University of Texas: Byron Tapley and Srinivas Bettadpur.
References
S. Aoki, H. Kinoshita, Note on the relation between the equinox and Guinot’s non-rotating origin. Celest.
Mech. 29, 335–360 (1983)
S. Aoki, B. Guinot, G.K. Kaplan, H. Kinoshita, D. McCarthy, P.K. Seidelmann, The new definition of univer-
sal time. Astron. Astrophys. 105, 359–361 (1982)
D.F. Argus, R.G. Gordon, No-net-rotation model of current plate velocities incorporating plate motion model
NUVEL-1. Geophys. Res. Lett. 18, 2039–2042 (1991)
S.W. Asmar, Radio as a science tool. Proc. IEEE 98, 10 (2010)
S.W. Asmar, J.W. Armstrong, L. Iess, P. Tortora, Spacecraft Doppler tracking: noise budget and achievable
accuracy in precision radio science observations. Radio Sci. 40 (2005). doi:10.1029/2004RS003101
J.G. Beerer, G.G. Havens, Operation the dual-orbiter GRAIL mission to measure the Moon’s gravity, in
SpaceOps 2012 Conference, Stockholm, Sweden, June 2012
C. Boucher, Z. Altamimi, L. Duhem, Results and analysis of the ITRF93. IERS Technical Note, 18, Obser-
vatoire de Paris, 1994
C. Dunn, W. Bertiger, Y. Bar-Sever, S. Desai, B. Haines, D. Kuang, G. Franklin, I. Harris, G. Kruizinga,
T. Meehan, S. Nandi, D. Nguyen, T. Rogstad, J.B. Thomas, J. Tien, L. Romans, M. Watkins, S.C. Wu,
S. Bettadpur, J. Kim, Instrument of GRACE: GPS augments gravity measurements. GPS World 14,
16–28 (2003)
E.G. Fahnestock, Comprehensive gravity and dynamics model determination of binary asteroid systems, in
American Astronomical Society, DPS Meeting 41, #50.11 (2009)
E.G. Fahnestock, R.S. Park, D.-N. Yuan, A.S. Konopliv, Spacecraft thermal and optical modeling impacts
on estimation of the GRAIL lunar gravity field, in AIAA/AAS Astrodynamics Specialist Conference,
Minneapolis, MN, AIAA, August 13–16, 2012, pp. 2012–4428
R. Floberhagen, P. Visser, F. Weischede, Lunar albedo forces modeling and its effect on low lunar orbit and
gravity field determination. Adv. Space Res. 23, 378–733 (1999)
W.M. Folkner, DSN station locations and uncertainties. JPL TDA Progress Report, 42-128, 1-34, Jet Propul-
sion Laboratory, California Institute of Technology, Pasadena, CA, 1996
W.M. Folkner, J.A. Steppe, S.H. Oliveau, Earth orientation parameter file description and usage. Interof-
fice Memorandum 335.1-11-93 (internal document), Jet Propulsion Laboratory, California Institute of
Technology, Pasadena, CA, 1993
J. Guinn, P. Wolff, TOPEX/Poseidon operational orbit determination results using global positioning satel-
lites, in AAS/AIAA Astrodynamics Specialists Conference (1993). AAS-93-573
S. Hatch, R. Roncoli, T. Sweetser, GRAIL trajectory design: lunar orbit insertion through science, in AIAA
Guidance, Navigation, and Control Conference, Toronto, Ontario, Canada, August 2010. AIAA 2010-
8385
The GRAIL Scientific Measurement System 55
W.A. Heiskanen, H. Moritz, Physical geodesy. Bull. Géod. 86(1), 491–492 (1967)
T.L. Hoffman, GRAIL: gravity mapping the Moon, in Aerospace Conference, 2009 IEEE, Big Sky, MT, 7–14
March 2009. ISBN 978-1-4244-2622-5
W.M. Kaula, Theory of Satellite Geodesy (Blaisdell, Waltham, 1966). 124 pp.
J. Kim, Simulation study of a low-low satellite-to-satellite tracking mission. Ph.D. dissertation, Univ. of Texas
at Austin, May 2000
J. Kim, B. Tapley, Error analysis of a low-low satellite-to-satellite tracking mission. J. Guid. Control Dyn.
25(6), 1100–1106 (2002)
W.M. Klipstein, B.W. Arnold, D.G. Enzer, A.A. Ruiz, J.Y. Yien, R.T. Wang, C.E. Dunn, The lunar gravity
ranging system for the gravity recovery and interior laboratory (GRAIL) mission. Space Sci. Rev. (2013,
this issue)
A.S. Konopliv, S.W. Asmar, E. Carranza, D.N. Yuan, W.L. Sjogren, Recent gravity models as a result of the
lunar prospector mission. Icarus 150, 1–18 (2001)
F.T. Krogh, Changing stepsize in the integration of differential equations using modified divided differences.
JPL Tech. Mem. No. 312, Section 914 (internal document), Jet Propulsion Laboratory, California Insti-
tute of Technology, Pasadena, CA, 1973
C.L. Lawson, R.J. Hanson, Solving Least Squares Problems. SIAM Classics in Applied Mathematics, vol. 15
(Society for Industrial and Applied Mathematics, Philadelphia, 1995)
R. Leavitt, A. Salama, Design and implementation of software algorithms for TOPEX/POSEIDON ephemeris
representation, in AIAA/AAS Astrodynamics Specialists Conference (1993). AAS-93-724
J.H. Lieske, T. Lederle, W. Fricke, B. Morando, Expressions for the precession quantities based upon the IAU
(1976) system of astronomical constants. Astron. Astrophys. 58, 1–16 (1977)
A.J. Mannucci, B.D. Wilson, D.-N. Yuan, C.H. Ho, U.J. Lindqwister, T.F. Runge, A global mapping technique
for GPS-derived ionospheric total electron content measurements. Radio Sci. 33(3), 565–582 (1998)
D.D. McCarthy, G. Petit (eds.), IERS Conventions, IERS Technical Note, vol. 32 (2003)
T.D. Moyer, Formulation for Observed and Computed Values of Deep Space Network Data Types for Navi-
gation (Wiley, Hoboken, 2003). 576 pp.
X.X. Newhall, J.G. Williams, Estimation of the lunar physical librations. Celest. Mech. Dyn. Astron. 66,
21–30 (1997)
R.S. Park, S.W. Asmar, E.G. Fahnestock, A.S. Konopliv, W. Lu, M.M. Watkins, Gravity recovery and interior
laboratory simulations of static and temporal gravity field. J. Spacecr. Rockets 49, 390–400 (2012)
S. Pines, Uniform representation of the gravitational potential and its derivatives. AIAA J. 11, 1508–1511
(1973)
R. Roncoli, K. Fujii, Mission design overview for the gravity recovery and interior laboratory (GRAIL)
mission, in AIAA Guidance, Navigation, and Control Conference, Toronto, Ontario, Canada, August
2010. AIAA 2010-8383
P.K. Seidelmann, 1980 IAU theory of nutation: the final report of the IAU working group on nutation. Celest.
Mech. 27, 79–106 (1982)
B.D. Tapley, B. Schutz, G. Born, Statistical Orbit Determination (Elsevier, Boston, 2004a). 547 pp.
B.D. Tapley, S. Bettadpur, M. Watkins, C. Reigber, The gravity recovery and climate experiment: mission
overview and early results. Geophys. Res. Lett. 31 (2004b). doi:10.1029/2004GL019920
B.J. Thomas, An analysis of gravity-field estimation based on intersatellite dual-1-way biased ranging. JPL
Publication 98–15, May 1999
S.G. Turyshev, V.T. Toth, M.V. Sazhin, General relativistic observables of the GRAIL mission. Phys. Rev. D
87, 024020 (2013)
J.M. Wahr, The forced nutations of an elliptical, rotating, elastic, and oceanless Earth. Geophys. J. R. Astron.
Soc. 64, 705–727 (1981)
R.C. Weber, P.-Y. Lin, E.J. Garnero, Q. Williams, P. Lognonne, Seismic detection of the lunar core. Science
331, 309–313 (2011). doi:10.1126/science.1199375
J.G. Williams, A scheme for lunar inner core detection. Geophys. Res. Lett. 34, L03202 (2007).
doi:10.1029/2006GL028185
J.G. Williams, D.H. Boggs, W.M. Folkner, DE421 lunar orbit, physical librations, and surface coordinates.
JPL IOM 335-JW, DB, WF-20080314-001, March 14, 2008
D.-N. Yuan, W. Sjogren, A. Konopliv, A. Kucinskas, Gravity field of mars: a 75th degree and order model.
J. Geophys. Res. 106(E10), 23377–23401 (2001)
M.T. Zuber, D.E. Smith, D.H. Lehman, T.L. Hoffman, S.W. Asmar, M.M. Watkins, Gravity recovery and
interior laboratory (GRAIL): mapping the lunar interior from crust to core. Space Sci. Rev. (2013, this
issue). doi:10.1007/s11214-012-9952-7