2.3 Coordinated Universal Time (UTC) clock used by DSN . . . . . . . . . . 6

3 LGRS DOWR measurement errors due to timing errors 6

3.1 LGRSDOWRmeasurementmodel ...................... 6

3.2 LGRS DOWR measurement errors due to expected timing performance oftheGRAILsystem ............................. 9

4 Algorithms for GRAIL timing 13

1 Introduction

The objective of the GRAIL mission is to determine with high accuracy the lunar gravity field for scientific research. The input data for the gravity field determination process are Ka phase measurements between the two GRAIL spacecraft which are used to compute the DOWR measurement. The DOWR measurement is then converted to instantaneous range, range rate and range acceleration measurements, which serve as inputs to the gravity field estimation process. Very accurate timing of the GRAIL measurements are crucial to achieve the high accuracy measurements needed for a high quality gravity field.

In this memorandum the GRAIL timing system and its performance will be discussed in terms of errors in the LGRS measurements. For this purpose we derived the LGRS model with second order terms and included the timing effects due to general relativity. The GRAIL mission can not use the GPS system for timing like the GRACE mission [6] did, but has to rely on the combination of several measurements to achieve the desired timing accuracy. Furthermore, we will asses how the design of the GRAIL timing system and its performance will impact the errors in the LGRS measurements. Based on this assessment, algorithms and processing choices will be presented and developed to insure that the required timing precision is achieved for the GRAIL mission.

2 Description of the GRAIL timing system

Figure 1: GRAIL clocks, models and measurements used for timing

This section describes the clocks and timing systems used by the GRAIL mission. The GRAIL mission uses the following three timing systems: 1) The LGRS clock, which is driven by an USO to create a highly stable clock, 2) the Onboard Spacecraft Clock (OSC), which is used in the timing correlation between the LGRS/OSC clocks and the time correlation between the OSC/DSN clocks and 3) the DSN clock, which is synchronized with highly stable time standard UTC. In figure 1 all the clocks used for GRAIL timing are shown with information or data flows between the different clocks and measurement systems. The ultimate goal is to provide a time correlation between the LGRS clocks and the TDB time scale, which will be provided in the Level1 product CLK1B.

2.1 LGRS clock

Each spacecraft has an LGRS clock which is driven by an USO to create a very stable clock. The LGRS clock is used for timing of the LGRS Ka phase measurement and the ranging data from the TTS. The LGRS clock has no notion of absolute time. Instead, the LGRS clock reading is with respect to the clock startup epoch. The required stability of the USO is listed in Table 1. The stability is given in terms of the Allan Deviation [1] which is the frequency change divided by the reference frequency (δf/f) averaged over a specific time interval.

Table 1: Allan Deviation requirements for GRAIL USO

Integration Time Allan Deviation (seconds)

13 _× 10−¹³ 10 3 _× 10−¹³ 100 3 _× 10−¹³ 1000 6 _× 10−¹³ 57600 7 _× 10−¹¹

2.2 Onboard spacecraft clock

The OSC is run by a crystal oscillator and has inferior stability characteristics compared to the LGRS USO. The OSC is used for time tagging all spacecraft data and the arrival of the LGRS data packets (which includes the LGRS 1 Pulse Per Second (PPS) packet) by the onboard computer. By time tagging the arrival of LGRS data packets and the arrival of the LGRS 1 PPS, a time correlation can be generated between the LGRS clock and the OSC. Furthermore, the OSC is used to time stamp spacecraft time correlation packets which are transmitted to a DSN station where the arrival time is recorded in UTC and thus providing a time correlation between the OSC and UTC. By combining the LGRS/OSC and OSC/UTC time correlation products, a time correlation between the LGRS and UTC can be determined and the OSC clock drops out. Hence the stability characteristics of the OSC do not affect the LGRS and UTC time correlation because the OSC error over short intervals (< 1 second) are < 1 micro second.

2.3 UTC clock used by DSN

The DSN uses very stable clocks which are maintained by the DSN Frequency and Timing Subsystem (FTS) [7]. The DSN time stamps the arrival of telemetry and radio metric tracking data in the UTC time scale, which is tied to the International Atomic Time (TAI) time scale. Based on FTS reports the real-time timing performance is at the 10−⁶ second level and post processing analysis improves the performance to the 10−⁹ second level.

3.1 LGRS DOWR measurement model

This subsection describes the mathematical model for the LGRS measurement in the presence of time tag errors of the phase measurement, phase errors due to the instrument, general relativity and random phase errors. Two mathematical descriptions are considered in this section, both were developed for the GRACE mission [4],[8]. We re-derived the mathematical model to account for differences between the GRAIL and GRACE timing systems. Additionally, we must consider general relativity effects when bounding errors in the LGRS measurement model.

The timing of an LGRS measurement has two components. First, absolute timing knowledge allows us to properly assign a LGRS measurement to a position around the moon. Secondly, the relative timing knowledge between the onboard LGRS clocks allows the proper alignment of the LGRS phase measurements of both spacecraft at a common coordinate time. This alignment is needed to form the DOWR measurement, which is processed by the science team.

For GRACE the absolute and relative timing is determined from the onboard GPS measurements [6]. On the other hand the GRAIL absolute timing is determined by the DSN and the relative timing is determined by the TTS onboard GRAIL. The main difference between the GRAIL and GRACE timing is that for GRACE the absolute and relative timing is available continuously during an orbit, while for GRAIL the relative timing is available continuously but the absolute timing is only available when the GRAIL spacecraft is visible from the Earth and tracked by the DSN. Furthermore, even when the DSN tracks a spacecraft, we only generate an absolute timing measurement once every 10 minutes. Based on planned DSN tracking and lunar orbit geometry, we expect less than 100 absolute timing measurements per day.

Given the limited absolute timing information for the GRAIL mission, the LGRS measurement model is re-derived and formulated in terms of the relative and absolute timing error. For this discussion the formulation and notation in Kim [4] was chosen because the Kim [4] formulation was in terms of the total time tag error. The Brooks

[8] formulation used time tag errors after an estimate of the time tag correction was applied. Furthermore both Kim [4] and Brooks [8] assume that the time tag errors of the phase measurements are small in the development of the LGRS measurement model. However, due to limited absolute timing information, the time tag error of the phase measurements may not be small and a bias may exist in the relative timing as well. In the re-derivation of the LGRS measurement model a second order expansion is used for the phase approximation to assess the LGRS measurement errors due to larger phase time tag errors and a bias in the relative timing.

Finally the LGRS measurement model described in Kim [4] and Brooks [8] did not consider general relativity. In this derivation the treatment of the light time is such that an assessment of not modelling general relativity on the LGRS measurement can be made and the effect of general relativity on phase and clock drift is considered as well. Under the assumptions and conditions described the derivation of the LGRS measurement model is as follows.

The single carrier frequency phase measurement for a spacecraft at a given coordinate time t can be written as

ϕ²₁(t +Δt1)= ϕ1(t +Δt1) ₋ ϕ²(t +Δt1)+ N₁² + d²₁ + �² (1)

1

ϕ¹₂(t +Δt2)= ϕ2(t +Δt2) ₋ ϕ¹(t +Δt2)+ N₂¹ + d¹₂ + �¹ (2)

2

where

ϕ²₁(t +Δt1) = differential phase measurement at spacecraft 1

ϕ¹₂(t +Δt2) = differential phase measurement at spacecraft 2

t = coordinate time

Δt1, Δt2 = time tag error at spacecraft 1 and 2 with respect to t

ϕ1(t +Δt1) = spacecraft 1 reference phase

ϕ²(t +Δt1) = received phase transmitted by spacecraft 2

ϕ2(t +Δt2) = spacecraft 2 reference phase

ϕ¹(t +Δt2) = received phase transmitted by spacecraft 1

N₁²,N₂¹ = integer ambiguities

d²₁,d¹₂ = phase shift due to instrument, offset, multipath, etc.

�²₁,�¹₂ = random measurement noise

It should be noted that the phase delay due to the ionosphere has been omitted from the formulation described in [4] because of the lack of a lunar ionosphere.

Each phase consists of the reference phase ϕ_i which corresponds to the constant reference frequency, and the phase error δϕi(t) due to oscillator drift, frequency instability and phase change due to general relativity and hence may be written as:

ϕ1(t)= ϕ₁(t)+ δϕ1(t) (3)

ϕ2(t)= ϕ₂(t)+ δϕ2(t) (4) Additionally, the received phase ϕⁱ(t) can be identified with the transmitted phase ϕi:

ϕⁱ(t)= ϕi(t ₋ τ_jⁱ) (5) where τ_jⁱ is the light time from the i-th satellite to the j-th satellite. It should be noted that the light times of the two phase signals are different from each other since the two satellites are moving. Substituting (3), (4), (5) into (1) and (2) yields

ϕ²₁(t +Δt1)= ϕ₁(t +Δt1)+ δϕ1(t +Δt1) ₋ ϕ₂(t +Δt1 ₋ τ₁²) ₋ δϕ2(t +Δt1 ₋ τ₁²)+ N₁² + d²₁ + �²₁ (6)

ϕ¹₂(t +Δt2)= ϕ₂(t +Δt2)+ δϕ2(t +Δt2) ₋ ϕ₁(t +Δt2 ₋ τ₂¹) ₋ δϕ1(t +Δt2 ₋ τ₂¹)+ N₂¹ + d¹₂ + �¹₂ (7)

Adding (6) and (7) gives the DOWR phase measurement Θ(t) definition:

Θ(t) _≡ ϕ²₁(t +Δt1)+ ϕ₂¹(t +Δt2)

= ϕ₁(t +Δt1) ₋ ϕ₁(t +Δt2 ₋ τ₂¹)+ ϕ₂(t +Δt2) ₋ ϕ₂(t +Δt1 ₋ τ₁²)

δϕ1(t +Δt1) ₋ δϕ1(t +Δt2 ₋ τ₂¹)+ δϕ2(t +Δt2) ₋ δϕ2(t +Δt1 ₋ τ₁²)+

(N₁² + N₂¹)+(d₁² + d¹₂)+(�₁² + �¹₂) (8)

Unlike in Kim [4], we choose to approximate our phase measurements by a second order expansion about coordinate time t. Thus the second order expansion for phase terms in equation (8) about the coordinate time t may be written as

1

ϕ_i(t +Δti)= ϕ_i(t)+ ϕ˙₁(t)Δti + ϕ¨ ₁(t)Δt² _i (9)

2 1 ₂

δϕi(t +Δti) δϕi(t)+ δϕ˙i(t)Δti + δϕ¨1(t)Δt_i +Δi (10)

≈ 2

where ϕ˙₁(t) is the constant nominal frequency fi and ϕ¨₁(t) = 0 because ϕ˙₁(t) is per definition a constant. δϕ˙is the combined frequency offset δf due to the LGRS clock frequency offset and general relativity . The δϕ¨1(t) represents the frequency drift δf^˙i, which consists of a linear fit to the oscillator frequency and the frequency drift due general relativity. The last term Δi in equation (8) represents the phase error resulting from integrating the frequency residuals of the linear fit of the oscillator frequency. To assess the size Δi we used five minute spaced GRACE clock solutions. From these clock solutions, the clock drift was computed for each five minute sample as well as a one hour clock drift using only hourly samples. The hourly clock drift estimates were then differenced with the half our clock drift estimates from the five minute samples. The size of Δi was then computed by multiplying the speed of light, the half hour clock drift difference and the light time of 10−³ sec. of Based on the GRACE clocks, the size of Δi results in range errors less than 1 micro meter. Assuming the GRACE and GRAIL USOs are similar then Δi can be neglected. Thus the different phase and phase error terms in equation (8) can be written as the following second order expansions:

ϕ_i(t +Δti)= ϕ_i(t)+ fiΔti (11)

ϕ_i(t +Δtj ₋ τ_jⁱ)= ϕ_i(t)+ fi(Δtj ₋ τ_jⁱ) (12)

δϕi(t +Δti) _≈ δϕi(t)+ δfiΔti + ¹₂ δf^˙iΔt_i ² (13)

δϕi(t +Δtj ₋ τ_jⁱ) _≈ δϕi(t)+ δfi(Δtj ₋ τ_jⁱ)+ ₂¹ δf^˙i(Δtj ₋ τ_jⁱ)² (14)

Substituting (11),(12),(13),(14) into (8) and combining where possible in terms of the relative timing error (Δt1 ₋ Δt2), equation (8) can be written as:

Θ(t) _≈ (f1τ₂¹ + f2τ₁²)+(δf1τ₂¹ + δf2τ₁²)+

(f1 ₋ f2)(Δt1 ₋ Δt2)+(δf1 ₋ δf2)(Δt1 ₋ Δt2)+

(δf^˙1Δt2 ₋ δf^˙2Δt1)(Δt1 ₋ Δt2)+ ¹(δf^˙1 + δf^˙2)(Δt1 ₋ Δt2)² +

211

δf^˙1(2Δt2 ₋ τ₂¹)τ₂¹ + δf^˙2(2Δt1 ₋ τ₁²)τ₁² +

22
(N₁² + N₂¹)+(d₁² + d¹₂)+(�₁² + �¹₂) (15)

The first term of equation (15) represents the true phase measurement and all the other terms represent errors due to phase time tag errors, general relativity, phase noise, frequency offset of the oscillator and random phase noise. The first four terms are identical to the terms described in Kim [4] who used a linear expansion. All other terms are the result from the second order expansion and many terms are a function of the relative timing error (Δt1 ₋ Δt2). The second order expansion also shows a dependency on the absolute timing error terms Δt1 and Δt2. It should be noted the DOWR phase has an unknown bias due to the fact that the phase integer ambiguities (N₁² +N₂¹) are not known. Thus the resulting DOWR range measurement has an unknown bias. In section

3.2 the performance of the GRAIL timing system will be discussed and an assessment of the LGRS measurement errors

3.2 LGRS DOWR measurement errors due to expected timing performance of the GRAIL system

This subsection describes the LGRS measurement errors due to expected timing performance of the GRAIL system and the effects of general relativity for the GRAIL system in lunar orbit as described in section 4.1. The expected timing performance listed in Table 2 is used to assess the error terms in the DOWR measurement as described in equation (15). The error assessment is made in terms of the DOWR which can be defined from Θ(t) in equation (15) as

Θ(t)

Rdowr(t) _≡ c (16)

f1 + f2

Table 2: Expected timing performance of the GRAIL timing system

Timing Bias Variability (seconds) (seconds) (1σ)

absolute (Δti) i=1,2 1 _× 10−¹ 1 _× 10−⁴

relative (Δt1 ₋ Δt2)1 _× 10−⁷ 9 _× 10−¹¹

where c is the speed of light. The error terms in equation (15) will be multiplied by a factor

^{c 3} × ¹⁰⁸ = _O(10−³) (17)

f1 + f2 ≈ 6 _× 10¹¹

to assess the error terms in range.

The first term of equation (15) only represents the true phase measurement for the nominal constant frequencies f1 and f2. The second term of equation (15) represents the part of the true measurement due to the oscillator frequency offset plus drift and the frequency drift due to general relativity. When the actual measurement is made, the actual frequencies (f1 + δf1) and (f2 + δf2) are used. For GRAIL we have a direct measurement of (f1 + δf1) and (f2 + δf2) using RSR with an accuracy of 10−³ Hz at X-band frequency. Therefore we need to combine the first and second term of equation

(15) to derive the total truth phase measurement. We start by defining:

f1 + δf1 = f1 + δf^¯1 + εf1 = f^¯1 + εf1 (18)

f2 + δf2 = f2 + δf^¯2 + εf2 = f^¯2 + εf2 (19)

where f^¯1, f^¯2 are the noise free frequency measurements by the RSR and δf^¯1, δf^¯2 represent the best estimate of the actual frequency offset plus a linear drift of the oscillator and the frequency offset due to general relativity. The RSR frequency measurement errors are described by εf1 and εf2. Substitution of equations (18) and (19) into the first two terms of equation (15) plus converting to DOWR range, results in

cc c

_¯¯Θ(t)=_¯¯(f^¯1τ₂¹ + f^¯2τ₁²)+ _¯¯(εf1τ₂¹ + εf2τ₁²) (20)f1 + f2 f1 + f2 f1 + f2

Using the RSR frequency accuracy of 10−³ Hz, Ka frequency = 32 _× 10−⁹ Hz and a light time of 10−³ sec then the second term of equation (20)

c 3 _× 10⁸

_¯¯(εf1τ₂¹ + εf2τ₁²) _{≈ 32 × 109 ×} 2 _× 10−³ _× 10−³ = _O(10−⁷) meter (21)

f1 + f2

can be neglected. Now we rewrite the light times τ₁² and τ₂¹ as:

τ₂¹ = τ + δτ₂¹ (22)

τ₁² = τ + δτ₁² (23) where τ is the euclidean instantaneous light time. Substitution of equations (22) and

(23) into equation (20) yields

_{f¯1 +} ^c_f¯2 Θ(t) _{≈ f¯1 +} ^c_f¯2 (f^¯1 + f^¯2)τ + _{f¯1 +} ^c_f¯2 (f^¯1δτ ₂¹ + f^¯2δτ₁²)

_≡ cτ + ρT OF
= ρ + ρT OF (24)

where δτ_jⁱ are the differences between the actual light time from satellite i to j (including general relativity) and the euclidean instantaneous light time τ. ρ is the euclidean instantaneous inter-satellite range. The ρT OF represents the range change due to spacecraft motion and general relativity. For GRACE the ρT OF was approximated in Kim [4] with sufficient accuracy. However, most assumptions for the GRACE approximation in

[4] were invalid for GRAIL. Therefore the MIRAGE software is used for the ρT OF by computing the one way ranges between the spacecraft which takes into account satellite motion and general relativity without approximations. Before GRAIL becomes operational, we will need to verify that our MIRAGE calculation of τ_jⁱ has an accuracy of better than 0.01 ps, because if our light time calculation has error less than 0.01 pico seconds, then the DOWR error due to light time errors ετ₂¹ and ετ₁² can be written as:

(f^¯1ετ₂¹ + f^¯2ετ₁²)

c _¯¯≈O(10−⁶) meter (25)

f1 + f2

which can neglected.

Terms three and four of equation (15) depends on the relative synchronization error (Δt1 ₋ Δt2). The expected relative synchronization performance is listed in Table 2. It should be noted that the bias in the relative synchronization does not add errors to the DOWR measurement but adds an additional bias to the DOWR measurement. Therefore, we can only evaluate the variability in relative clock synchronization when bounding terms four and five of equation (15). Again we will combine terms three and four of equation (15) to rewrite these terms with the observed frequencies from the RSR measurement. Using definitions (18) and (19) plus converting to DOWR range we may write terms three and four as:

_(f¯_{1 − f}¯₂₎

c

c (Δt1 ₋ Δt2)+ (εf1 ₋ εf2)(Δt1 ₋ Δt2) (26)¯¯¯¯

f1 + f2 f1 + f2

We know for GRAIL that the baseband frequency ( f^¯1 ₋ f^¯2) _≈ 6 _× 10⁵ Hz and the Ka frequency f^¯i _≈ 32 _× 10⁹ Hz. Thus we can bound the first term of equation (26) by:

c ^(f¯1 − ^f¯2)(Δt1 ₋ Δt2) _≈ 3 _× 10^{8 6} × ¹⁰⁵ = _O(10−⁶) meter (27)

¯¯_{6 × 10}10 _× ⁹ _× ¹⁰⁻¹¹

f1 + f2

I should be noted that the ( f^¯1 ₋ f^¯2) is not constant in this context unlike the derivation in Kim [4]. Therefore any variability will have an effect due to the bias in (Δt1 ₋ Δt2). However, based on the expected clock performance we do not expect a large effect on the LGRS measurement. Furthermore, the level-2 processing will have model parameters to observe drifts in the LGRS measurement in case the LGRS clock does not have the expected performance. Finally, for term one in equation (26), the differential effect of general relativity on the frequency of both GRAIL clocks can be neglected because the frequency difference is < 3 Hz based on the maximum relative clock drift in figure 3. The second term in equation (26) may be bounded by

c 10−³

_¯¯(εf1 ₋εf2)(Δt1 ₋Δt2) _≈ 2_×3_×10⁸ _{6 × 1010 ×}9_×10−¹¹ = _O(10−¹⁴) meter (28)

f1 + f2

Before bounding term five of equation (15) we first rewrite this term as

cc

(δf^˙1Δt2 ₋ δf^˙2Δt1)(Δt1 ₋ Δt2)= (δf^˙1 ₋ δf^˙2)Δt2(Δt1 ₋ Δt2)+

f1 + f2 f1 + f2
c

δf^˙1(Δt1 ₋ Δt2)² (29)f1 + f2

The first term of equation (29) may be bounded

c 10−⁹

_{f1 + f2} (δf^˙1 ₋ δf^˙2)Δt2(Δt1 ₋ Δt2) _≈ 2 _× 3 _× 10⁸ _{6 × 1010} Δt2 _× 10−⁷

=Δt2 _×O(10−¹⁷) meter (30)

For this bound a δf^˙i < 10−⁹ was used based on performance of the GRACE clocks. Even though the bound is a function of the time tag error Δt2, the error it self is _� 10−⁶ meter for large time tag errors. The second order term of equation (29) and term six in equation (15) can both be bounded by (Δt1 ₋ Δt2)²

c ^(Δt1 − ^Δt2)2 _≈ 3 _× 10^{8 10−14} = _O(10−¹⁶) meter (31)

f1 + f2 6 _× 10¹⁰

Since all the multipliers for the (Δt1 ₋ Δt2)² terms in equation (15) are _� 1, we can neglect all these terms.

Finally terms seven and eight of equation (15) do also depend directly on the absolute time synchronization errors Δt1, Δt2 and the time of flight. Only the error bound for term seven of equation (15) is discussed because term eight is similar and can be bounded in the same way. Thus the seventh term can be written as:

1 _cδf^˙1(2Δt2 ₋ τ₂¹)τ₂¹ _{= c} δf^˙1 _{τ21Δt2 −} 1 _c δf^˙1 _(τ21)2

2 f1 + f2 f1 + f2 2 f1 + f2
10−⁹ 1 10−⁹

_≈ 3 _× 10⁸ _{6 × 1010} Δt2 _{− 2 ×} 3 _× 10⁸ _{6 × 1010 ×O}(10−⁶)

=5 _× 10−¹²Δt2 + _O(10−¹⁷) meter (32)

which shows that even for an absolute synchronization error of days the error in the DOWR measurement is less than 10−⁶ meters. For this bound a δf^˙1 < 10−⁹ was used based on performance of the GRACE clocks. Furthermore, the δf^˙1 < 10−¹² due to general relativity based on data in figure 2. As listed in Table 2, we do not expected absolute synchronization errors larger than 0.1 seconds, hence the seventh term should not provide a significant error source.

In summary with the expected GRAIL timing performance, we do not expect that the errors in the DOWR measurement are not going to exceed 10−⁶ meters.

4.1 Modelling LGRS clock rate due to general relativity

This section describes how the LGRS clock rate due to general relativity is modeled for the GRAIL mission. For this discussion it is assumed that the LGRS clock is perfect thus its drift with respect to coordinate time would be zero if it were infinitely far from any mass. According to the theory of general relativity this perfect clock in lunar orbit will have a clock rate which differs from the TDB clock rate. The change in clock rate due to general relativity is described in [5] and can be written as:

dτ U 1 v²

_dt =1 _{− c2 − 2 c}+ L (33)

2

where τ is proper time, which is the perfect clock reading in lunar orbit, t is coordinate time (TDB), U is the gravitational potential at the GRAIL spacecraft in lunar orbit, c is the speed of light, v is the Solar-System barycentric velocity of the GRAIL spacecraft and L is a constant value. The value of L is chosen such that the clock rate of TDB is equal to the average clock rate of a clock on the Earth surface at sea level with respect to TDB [5]. Thus L may be written in term of the Solar-System barycentric gravitational potential of a clock at Earth’s sea level UE and the barycentric velocity vE of that clock [5].

11 ₂

L = ₂ <UE + v_E > (34)

c2 The brackets <> denote a long term average. The reason for selecting this L is that our notion of time is established by atomic clocks on the Earth’s surface also known as TAI. Choosing another L would cause a linear trend between TDB and TAI. The DSN uses UTC for timing, which is offset by a constant with respect to TAI but this constant changes every time a new leap second is introduced.

For GRAIL equation (33) is integrated to get proper time as a function of coordinate time, using the GRAIL spacecraft ephemerides and the planetary ephemerides DE421[3] for the ten principle solar bodies (point mass). In this integration L =1.550520 _× 10−⁸ is used [5]. In Figure 2 the GRAIL clock drift rate due to general relativity is shown

Figure 2: GRAIL clock drift rate due to general relativity during the science phase

Figure 3: Differential GRAIL clock drift rate due to general relativity during the science phase

for the first 30 days of the GRAIL science phase, where the spacecraft separation starts at 85 km and ends at 225 km. This figure contains the effects of the GRAIL orbit, lunar orbit and the Earth/Moon system orbit around the Sun. The maximum drift rate is about 2 _× 10−⁹ sec/sec. This number will be used in other sections to bound LGRS and TTS measurement errors due to general relativity.

The relative clock drift rate between the two GRAIL clocks for the first 30 days of the GRAIL science phase is shown in Figure 3. From this figure it can be seen that the differential clock drift is less than _±40 _× 10−¹² sec/sec. It should be noted that the differential clock rate is shown for the orbit geometry of the first 30 days of the science phase. However, an orbit geometry to maximize the differential clock drift error puts an upper bound of _±70 _× 10−¹² sec/sec [9]. The upper bound of 70 _× 10−¹² sec/sec will be used to bound LGRS and TTS measurement errors due to differential clock drift effects of general relativity

4.2 Modelling the LGRS clock offset due to LGRS clock imperfections

This subsection describes the model used for the LGRS clock drift due to imperfections in the LGRS clock. This model does not include LGRS clock drift due to general relativity. The model described in this subsection when added to the general relativity contribution from subsection 4.1 gives the total LGRS clock drift as seen from the TDB time scale.

Unlike GRACE, for which GPS measurements allow us to track clocks continuously, we do not have continuous clock offset measurements for GRAIL. GRAIL may have only 20 LGRS clock offset measurements per day. Moreover, these measurements are biased and a have variability of 20 milli seconds (1 σ). The bias in the clock offset measurement are due to unknown timing delays in the spacecraft. The variability in the clock offset measurement is caused by unknown data packet traffic on the spacecraft bus, which causes a variable transfer time of the time correlation packet from the onboard computer to the spacecraft antenna.

As mentioned in subsections 2.1 and 2.3, we know that the LGRS clock and the clocks that realize the UTC time tags at the DSN are very stable. Therefore, a simple polynomial may be sufficient to describe the LGRS clock drift with respect to TDB, where TDB time tags are determined from UTC using general relativity by integration of equation 33 for a DSN station ( see subsection 4.1). It should also be noted that the clock offset measurements also contain a contribution due to general relativity. The general relativity contribution is computed as described in subsection 4.1 and removed from the measurements before fitting a second order polynomial to the measurements.

Using 60 days of clock offset measurements as observed by GRACE, an assessment is made of modeling a second order polynomial for the clock drift error. The GRACE clock serves as a proxy for the LGRS clock since both clocks have similar stability characteristics. The residual for the second order fit is shown in Figure 4.

From Figure 4 it can be seen that throughout the 60 days, the time tag error using a second order polynomial for the LGRS clock drift would be less than 2 micro seconds. Therefore the error in the clock offset is 2 micro seconds or less due to the the non-quadratic LGRS clock drift behavior. The use of higher order polynomials was found not to reduce the time tag errors significantly, therefore we will use a second order polynomial to model the LGRS clock drift due to imperfections in the LGRS clock. If clock drift should prove larger than we hope, we will allow our second order polynomial coefficients to random walk.

Figure 4: Time tag residuals for GRACE-2 clock after second order polynomial fit

4.3 Algorithm for inter-satellite LGRS clock offset determination

This subsection describes the algorithm and assumptions for the extraction of the inter-satellite LGRS clock offset using the GRAIL TTS inter-satellite ranging system. A complete description of the TTS ranging algorithm can be found in [2]. For this discussion it is assumed that an inter-satellite range measurement is available from both GRAIL spacecraft and that the objective is to determine the LGRS clock offset at a given coordinate time. It should be noted that the relative LGRS clock relationship derived in this section is the inverse of the relative coordinate time clock offset in equation (15) (see section 3.1). The relative LGRS clock offset derived in this section will be used as an observation equation in the determination of the LGRS clock error model described in section 4.4.5. The inversion of the relative LGRS clock offset is part of the complete clock offset model inversion described in section 4.4.6

The objective is to determine the clock offset O12(T ) of the two LGRS clocks at a given coordinate time T which may be written as

t˜S1(T ) ₋ t˜S2(T )= O12(T ) (35)

where t˜1(T ) is the satellite 1 LGRS clock reading at coordinate time T and t˜2(T ) is the satellite 2 LGRS clock reading at coordinate time T . Our first assumption is that the interpolation errors of the inter-satellite range are small for both spacecraft. Therefore

with interpolation it is possible to change the time of the data and thus change the clock reading without having to restart the range code generation. At satellite 1 we measure:

t˜¹ _S1(T1) ₋ t˜¹ _S2(T1 ₋ τ₁²)+ ν (36)

where τ₁² is the light time of the ranging signal from spacecraft 2 to spacecraft 1 and ν is random noise (assumed to be 10 nsec for GRAIL). In this measurement t˜1(T1) is known and τ² is unknown. Expanding t˜¹ (T1 ₋ τ₁²) about T1

1 S2

t˜_S¹₂(T1 ₋ τ₁²)= t˜_S¹₂(T1) ₋ t˜^˙_S¹₂(T1)τ₁² +Δ_S¹₂ + _O((τ₁²)²) (37)

where t˜^˙_S¹₂(T1) is a linear fit to the clock drift of spacecraft 1 with respect to coordinate time at T1. This clock rate is the sum of the clock rate due to general relativity and the LGRS clock error rate. Δ¹ _S2 represents the residuals about the linear fit of the clock drift of spacecraft 1. The last term in equation (37) represents the second order term caused by general relativity. The size Δ¹ _S2 is based on the same GRACE clock analysis used for equation (10) and can be bounded by 10−¹² sec over the GRAIL light time of τ₁² = 10−³ sec. Assuming the GRAIL USO is similar to the GRACE USO, we can neglect Δ¹ _S2. To bound the last term in equation (37) we determined an upper bound of _O((10−¹²) for ¨¹

t˜_S2 based on general relativity and shown in figure 2. Thus the last term of equation

(37) can be bounded by

1 12

¨

_O((τ₁²)²) _≈ t˜_S2(T1)(τ₁²) < _O(10−¹²) _×O(10−⁶)= _O(10−¹⁸) sec (38)

2

Substituting (37) into (36) we measure

11

t˜_S1(T1) ₋ t˜_S2(T1)+ t˜^˙_S¹₂(T1)τ₁² + ν = O12(T1)+ t˜^˙_S¹₂(T1)τ₁² + ν (39)

Thus with no knowledge of τ₁² we know O12(T1) to 10−³ sec since the clock drift rate is

t˜^˙1 _S2(T1)=1+ r where r = _O(10−⁶) sec/sec (40)

Parameter r in equation (40) is determined from actual GRACE clock error drift observation and found to be _O((10−⁶) assuming the GRACE and GRAIL USOs are similar. The clock drift due to general relativity is of _O((10−⁹) thus the bound for parameter r is _O((10−⁶). Since the LGRS clocks of satellite 1 and satellite 2 are synchronized to 10−³ sec level we can resample R2. Thus we have now the following range measurements R1,R2 (in terms of seconds)

R1 = t˜¹ _S1(T1) ₋ t˜¹ _S2(T1 ₋ τ₁²)+ ν1 (41)

R2 = t˜² _S2(T2) ₋ t˜² _S1(T2 ₋ τ₂¹)+ ν2 (42)

where _|t˜¹ (T1) ₋ t˜² (T2)_| < 10−³ sec and given the clock rate in equation (40) then

S1S2

T1 ₋ T2 = _O(10−³). Substituting (39) into (41) and (42) results in

R1 = O12(T1)+ t˜^˙1 _S2(T1)τ₁² + ν1 (43) ₋R2 = O12(T2) ₋ t˜^˙2 _S1(T2)τ₂¹ ₋ ν2 (44) Expanding O12(T2) about T1 results in O12(T2)= O12(T1)+ O^˙12(T1)(T1 ₋ T2)+ΔO12 (45) where term ΔO12 can be neglected using similar arguments as for Δ¹ _S2 in equation (37). Adding equations (43), (44) and substituting (45) results in ˙

R1 ₋ R2 _≈ 2O12(T1)+ O^˙12(T1)(T1 ₋ T2)+ t˜^˙1 _S2(T1)τ₁² ₋ t˜² _S1(T2)τ₂¹ + ν1 ₋ ν2 (46) To understand the error terms we rewrite the light times τ₁² and τ₂¹ as: τ₁² = τ(T )+ δτ₁² (47) τ₂¹ = τ(T )+ δτ₂¹ (48)

where τ is the euclidean instantaneous light time at a given coordinate time T . We may assume that the LGRS clock rates t˜^˙_S²₁(T2), t˜^˙_S¹₁(T1) are identical since T1 and T2 are already synchronized to 10−³ sec after one iteration and using ^¨ t˜_S¹₂ from equation (38)

we can show t˜^˙2 _S1(T2)= t˜^˙1 _S1(T1)+ _O(10−³) _×O(10−¹⁰) sec (49) Substituting (47),(48),(49) into (46) results in:

˜

R1 ₋ R2 _≈ 2O12(T1)+ O^˙12(T1)(T1 ₋ T2)+(t˜^˙_S¹₂(T1) ₋ t^˙_S¹₁(T1))τ(T1)+

t˜^˙_S¹₂(T1)δτ ₁² ₋ t˜^˙_S¹₁(T1)δτ₂¹ + ν1 ₋ ν2 (50) _≈ 2O12(T1)+ O^˙12(T1)(T1 ₋ T2)+(t˜^˙_S¹₂(T1) ₋ t^˙_S¹₁(T1))τ(T1)+

˜r2δτ₁² ₋ r1δτ₂¹ + δτ₁² ₋ δτ ₂¹ + ν1 ₋ ν2 (51)

where r1 _≡ t˜^˙_S¹₁(T1) ₋ 1 (52) r2 _≡ t˜^˙_S¹₂(T1) ₋ 1. (53)

Substituting the expected values for clock drifts (40)), clock synchronization error and light time, then the terms of (51) can be approximated as follows: O^˙12(T1)(T1 ₋ T2)= _O(10−⁶) _×O(10−³)= _O(10−⁹) sec (54) (t˜^˙1 _S2(T1) ₋ t˜^˙1 _S1(T1))τ (T1)= _O(10−⁶) _×O(10−³)= _O(10−⁹) sec (55) r2δτ ₁² ₋ r1δτ ₂¹ = _O(10−⁶) _×O(10−⁷)= _O(10−¹³) sec (56) δτ₁² ₋ δτ ₂¹ = _O(10−⁷) sec (57) ν1 ₋ ν2 = _O(10−⁸) sec (58)

The first observation that can be made is that the clock drift difference multiplied by the light time or synchronization error is determining the size of the of the error terms. It can also be seen that the clock drifts due to general relativity do not play a role because these drifts are nearly identical (difference < 70 psec/sec) and cancel out when computing the terms of equation (51).

The term in equation (54) can be made small by iterating equation (51), so it can be neglected. The iteration procedure will be described later when the algorithm of the inter-satellite clock offset extraction is discussed.

The second term (55) is a variable term because the light time τ(T1) varies over the science phase from about 0.3 to 0.7 milli seconds. If the clock drift rates are of order _O(10−⁶) then this term needs to be computed and used as a correction since the size of this term is larger than the expected variability in the TTS measurement of 9 _× 10−¹¹ sec (1σ). The relative clock rate due to general relativity is of order _O(10−¹¹) and can be neglected

The third and fourth terms (56) represent the change in light time due to satellite motion and general relativity. The size of the third term is well below the expected noise of the TTS measurement of 9 _× 10−¹¹ sec (1σ) therefore it can be neglected. The fourth term needs to be computed. Given our expected orbit determination accuracy, residual error should fall at or below _O(10−¹¹).

The fifth term (58) is the error due the random noise in the range measurements. Given the size of _O(10−⁸) sec, smoothing is required. We can use carrier phase smoothing since both phase and range measurements are available from the TTS system. The noise on the phase measurements are significant lower than the range measurements. However, the phase measurements have an unknown bias due to the phase ambiguity. With carrier phase smoothing the bias in the phase measurements is computed using the range measurements. The computed phase bias is then used to adjust all phase measurements to form new range measurements with less noise.

The procedure for carrier phase smoothing is as follows. Given are range R(t˜i) and phase Φ(t˜i) measurements (expressed in range) at the same observation time t˜i which may be written as

R(t˜i)= R^¯(t˜i)+ νi i =1,N (59) ¯

Φ(t˜i) = Φ(t˜i)+ K + εi i =1,N (60) ¯¯

where K is the phase bias, R(t˜i) is the truth range and Φ(t˜i) is the truth phase range. The range and phase noise are denoted by νi and εi. Equations (59) and (60) can now be used to solve for K under the assumption that R^¯(˜Φ(˜

ti)= ^¯ti)

N N

1 ^� 1 ^�

K = (R(t˜i) ₋ Φ(t˜i))+ (νi ₋ εi) (61)

N N

i=1 i=1

For large N the second term of equation (61) will approach zero. N can be chosen over time intervals for which continuous phase measurements are available. Finally the estimate for K is removed from the phase range measurement to form the new range measurements.

In conclusion the algorithm for extracting the inter-satellite clock offset is as follows. First perform carrier phase smoothing for both spacecraft’s range and phase measurements according to

R^˜1(t˜¹ _i ) = Φ(t˜¹ _i ) ₋ K1 i =1,N (62)

R^˜2(t˜² _i ) = Φ(t˜² _i ) ₋ K2 i =1,N (63)

where t˜_i ¹) and t˜_i ²) are the observation times reported by the LGRS clock. Next iteratively apply equation (51) to minimize the synchronization error in coordinate time (T1 ₋ T2) shown in equation(54). The iteration procedure starts by computing the first estimate

_O˜1 _t1

of the clock offset ₁₂(˜_i )

_O˜11 ₍˜12 _(˜˙_˜˙

₁₂(t˜_i )= ¹₂R1(t˜_i ) ₋ R^˜2(t˜_i )) ₋ ¹₂t ¹ _{S2 −} t_S¹₁)τ (64)

where we assume the LGRS clock drifts t˜^˙_S¹₁, t˜^˙_S¹₂ and the light time τ at observation time t˜¹ are known. In general the estimate for the clock offset O^˜1 t¹ _i ) is not equal to (t˜¹ t_i ²).

i 12^(˜i ₋ ^˜We can now use the difference to update t˜² _i by defining

=(t˜¹ t² O¹ t¹ (65)

^δ1 i ₋ ^˜i ⁾ ₋ ^˜12^(˜i ⁾

where 1 denotes the iteration number of the t˜² _i update. Using δ1 to update t˜_i ²) we can now update the calculation for the relative time offset

O^˜2 t¹ = ¹(R^˜1(t˜¹ R2(t˜^{2 1}(t˜^˙1 t˜^˙1 (66)

12^(˜i ⁾₂i ⁾ ₋ ^˜i ^{+ δ}1⁾ _{− 2}S2 ₋ S1^)τ

Thus the complete iteration procedure of M iterations for equation (64) until t˜² _i converges, may be written as

δM =(t˜¹ t² O^M−¹(t˜_i ¹) (67)

i ₋ ^˜i ⁾ ₋ ^˜12

M

_O˜M+1_(˜1 ¹₍ ˜12 ^{� 1}_(˜˙¹ _˜˙¹

₁₂ t_i )= R1(t˜_i ) ₋ R^˜2(t˜_i + δk) ₋ t_{S2 −} t_S1)τ (68)

22

k=1

It should be noted that the convergence of t˜_i ² depends on the phase range noise εi. The t2 ^�M t2

evaluation of R^˜2(˜_i + _i=1 δk) is done with lagrangian interpolation of the R^˜2(˜_i ) data.

Depending on the size of the spacecraft’s phase range noises ε_i^j, additional smoothing may be needed of the time series after M iterations.

O^˜M t_i ¹)+ ε¹ _i + ε² (69)

12 ^(˜i

In case additional smoothing is needed, a simple box car filter will suffice to attenuate the remaining ε_i^j noise. The derived time series in equation (69) is known as a function of the LGRS observation time t˜¹ _i . In the final step the t˜¹ _i time tags are changed to TDB time using the time correlation function described in section 4.4.6.

4.4 Measurements and Algorithm for determination of LGRS clock error model

This subsection describes the measurements and algorithm for the determination of the LGRS clock error model. The LGRS clock error models the observed LGRS clock drift after the clock drift due general relativity has been removed. As described in section 4.2 the LGRS clock error is modelled as a second order polynomial with TDB time as the independent variable. The second order polynomial will model a clock bias, clock drift due to a frequency offset and linear drift in the frequency. The measurements available to estimate the LGRS clock drift model parameters are shown in Table 3.

Table 3: Measurements available for LGRS clock error model determination