Recovering Differential Forces From the GRACE-D Accelerometer

Abstract

The reduction of non-gravitational accelerations through an onboard accelerometer is a vital part of accurately estimating gravity fields for Gravity Recovery and Climate Experiment (GRACE) Follow-On (GRACE-FO). Current accelerometer quality, specifically a defect on the GRACE-D accelerometer, necessitates the use of extensive calibration techniques to optimally utilize the accelerometer measurements. This paper describes these techniques and correction algorithms as they are used for the official GRACE-FO accelerometer products, which rest on the information present in the partially functional GRACE-D accelerometer. First, transplanted data from GRACE-C is used to remove common effects between the spacecraft. Then, GRACE-D data is utilized to compute a correction to the transplanted data to account for solar radiation pressure (SRP) and albedo using a parity argument and to account for non-common drag and thruster leaks in the attitude control system. The SRP and albedo term predominantly improves C_3,0 while the drag and leaky thruster term improves the gravity field across the entire spectrum. Finally, some comments on alternative data usage, initial modifications to the algorithms to account for other spacecraft operational modes, and future outlook are provided.

Introduction

The Gravity Recovery and Climate Experiment (GRACE) Follow-On (GRACE-FO) mission (Kornfeld et al., 2019; Landerer et al., 2020), a successor to the GRACE mission (Tapley et al., 2004), recovers gravity by measuring inter-spacecraft acceleration between two spacecraft (GRACE-C and GRACE-D) formation flying in a common orbit, then subtracting off non-gravitational acceleration. A dual K-band receiver/transmitter (alternatively, a Laser-Ranging Interferometer (LRI)—Abich et al., 2019) measures the inter-spacecraft acceleration—in effect, echo-locating by measuring the round-trip travel time from GRACE-C to GRACE-D. Since non-gravitational acceleration can't currently be modeled at the necessary precision (except at times of low atmospheric drag), the non-gravitational forces are measured by an accelerometer on each spacecraft (Christophe et al., 2015).

As described in Harvey et al. (2022), noise levels on the GRACE-D accelerometer increased by orders of magnitude in June 2018, and have not recovered since. True signal is still recorded underneath the noise, leaving an accelerometer with clean readings along one linear axis and solid measurements of acceleration buried under jumpy, thermally correlated noise (thermals) along the other two axes. Note that throughout this paper, we use Radial to denote the body fixed Z-axis (rotation about this axis is yaw), Along to denote the body fixed X-axis (rotation about this axis is roll), and Cross to denote the body fixed Y-axis (rotation about this axis is pitch).

Previously, gravity field recoveries have ignored the GRACE-D accelerometer, and transplanted measurements from GRACE-C to GRACE-D by algorithms first used for the later part of the GRACE mission (Bandikova et al., 2019), sometimes augmented by non-gravitational force models (Behzadpour et al., 2021; Krauß et al., 2020).

However, even a single axis gives a great deal of information—and the other axes aren't a dead loss. Setting aside outgassing and thruster activity, the main sources of force on a LEO spacecraft are light—primarily solar radiation pressure (SRP), secondarily albedo (Earth radiation pressure) and thermal emission pressure—and matter—drag and lift. Both sources follow predictable physical laws. Drag hits along-track acceleration hardest, and impacts radial in a proportion predictable from spacecraft attitude relative to local atmospheric flow. SRP repeats at an orbital revolution, and when the attitude is well controlled most of the signal obeys symmetry laws around orbit angle zero, where the orbit angle is the angle in the orbit plane relative to the Sun. An orbit angle of zero has the Sun most directly overhead, noon for the spacecraft.

In this paper, we build a framework for modeling improved accelerometer data on GRACE-D. This framework is the foundation of the new ACH acclerometer data product released by the Level-1 processing performed at the Jet Propulsion Laboratory (JPL) (GRACE-FO, 2019b). This product is designed to be used as the nominal accelerometer product for GRACE-D; however, we have also included different components which make up the ACH product, as shown in Equation 1, in the public release (GRACE-FO, 2019b). These components are discussed in detail in this paper. [Image Omitted. See PDF]

Note that we do not discuss the AC0 component as this is the thruster model, which is unchanged from that given in Harvey et al. (2022).

As a preamble to the ACH product discussion, some important preliminary data analysis notes are provided in Section 2. The modeling framework, which can be divided into the different components shown in Equation 1, is then outlined as below, with relevant gravity field estimates provided to showcase their effect. Note that this modeling effort is specific to spacecraft operations in nominal Attitude and Orbit Control System pointing (or tight deadband mode).

1. RCYRPA. Jumpy, drifting, thermally correlated errors on GRACE-D, with two linear factors, one huge, one moderate—and relatively well behaved in late 2022/early 2023. These features are fully described by Harvey et al. (2022), but relevant features are recapped in Section 3.

2. SRP. SRP/albedo on GRACE-C and GRACE-D which basically obey parity rules. The theoretical motivation of the SRP component is discussed in Section 5 and the specifics of the algorithm for computing the ACS component of the ACH product are given in Section 6.

3. Drag. Drag on GRACE-C and GRACE-D which mostly cancels out in the differential drag combination, a transplant from GRACE-C to GRACE-D which comprises the ACU component of the ACH product, is discussed in Section 4. Additionally, within the drag component we account for scale factor differences between each accelerometer and smaller scale factor differences in each accelerometer axis, all within instrument tolerances.

4. Differential drag and leaky thrusters. Environmental drag differences between GRACE-C and GRACE-D and signal on both spacecraft due to low level thruster leaks, which are discussed in Section 7 and comprise the ACM component of the ACH product. Note that drag differences manifest mostly in along-track and the leaky thruster signal primarily occurs in cross-track, with a secondary effect in radial (driven by the thruster nozzle design). The leaky thruster signal manifests in the accelerometer data as a bias jump at each thruster firing, plus some drift between thruster firings (see Figure 1 for an example of bias jumps at thruster firings).

Finally, in Section 8 we briefly show that an alternative data observable (which has been used before in the context of GRACE-FO, centered on range acceleration measurements, as described in Allgeyer et al., 2022) which localizes the effect of errors, substantially increases science tolerance to accelerometer errors in recent higher drag months, and in Section 9 introduce recent work on modeling spacecraft acceleration differences when attitude control is loosened (wide deadband mode) in an effort to mitigate the effect of leaky thrusters.

Figure 1. Example of bias jumps at thruster firings due to the thruster leak on GF1.

Data Analysis and Frequency Ranges

Before beginning with the analysis of the accelerometer data, it is useful to define our data analysis techniques which were used to isolate particular frequency ranges for data analysis. The signal types, referenced throughout this paper are:

1. Low frequency signal (LFS)–Signal at frequencies equivalent to 1 cycle per orbital revolution (cpr) and lower.

2. High frequency repetitive signal (HFRS)–Signal at frequencies higher than 1 cpr that repeats from orbit to orbit.

3. High frequency non-repetitive signal (HFNRS)–Signal at frequencies higher than 1 cpr that does NOT repeat from orbit to orbit.

Note that before deriving these data types, there is a nominal starting point for the data stream, typically the approximately 10 Hz ACC1A data (GRACE-FO, 2019a). This starting point is achieved by applying the following steps:

1. Remove thruster firings and outliers as described by Harvey et al. (2022).

2. Compute 1 s averages. This is performed to decimate the data from 10 to 1 Hz and improve computational efficiency.

3. Compute a 101 s running median. The median filter is chosen to be resistant to outliers, which are frequent in the GRACE-D accelerometer data, and the 101 s filter width is empirically derived as a compromise between outlier removal and data attenuation.

With this as the common starting point, we can now more specifically define the signal types referenced throughout the paper.

LFS is derived by computing a running average equivalent to the orbital period. This procedure retains low frequency components (at 1 cpr and lower) while removing any higher frequency variations.

HFRS is derived to capture the high frequency signal (frequencies above 1 cpr) that repeats from orbit to orbit. More specifically, HFRS is derived as follows:

1. Subtract LFS from the time series.

2. Bin the data by orbit angle.

3. Compute the median per orbit angle.

4. Convert the orbit angle data into a time series. This time series will be identical for every orbital revolution.

HFNRS is now easily derived to capture the high frequency signal (frequencies above 1 cpr) that do NOT repeat from orbit to orbit. This is computed as follows:

1. Subtract LFS from the time series.

2. Subtract HFRS from the time series.

3. Fit a least squares spline with 60 s knots (note that 60 s was chosen to smooth out any high frequency outliers–an exhaustive search among different knot intervals was not performed).

Each of these time series highlights different aspects of the accelerometer signal, enabling detailed analysis and separation of various effects. These time series will be extensively used and referenced throughout the paper.

The GRACE-D Accelerometer

As on GRACE, the GRACE-FO accelerometers measure along-track, cross-track, and radial linear accelerations, as well as roll, pitch, and yaw angular accelerations. On GRACE-D, a noisy electrode plate pair, non-physical repetitive-per-orbital-revolution (also called per-rev) noise in roll, and a large cross-track sensitivity to thruster leaks make yaw, roll, and cross harder to use (Harvey et al., 2022), so in this note, we will focus on the more tractable radial, along, and pitch measurements. Gravity field estimation depends mostly on the accuracy of along-track and radial accelerations. In fact, the accelerometers were specifically designed to be more accurate in along-track and radial, less accurate in the less important cross-track direction (Christophe et al., 2015).

Both accelerometers can operate in different modes—usually, Normal Range Mode (NRM), or Large Range Mode (LRM). Note that the GRACE-D accelerometer was briefly operated in Bias Verification Mode (BVM), which is typically only used for ground testing, for some on-orbit diagnostics. For simplicity, we will discuss NRM operations, though the analysis also applies to LRM and BVM. For reference, the operational modes of the accelerometer are:

1. LRM1 = 28 May 2018. On initial activation, both accelerometers briefly in LRM.

2. NRM1 = 28 May 2018 to 20 June 2018. Both accelerometers in NRM.

3. LRM2 = 21 June 2018. Both accelerometers in LRM.

4. NRM2 = 22 June 2018 to 28 February 2019. From June 22 to February 7, both accelerometers in NRM. From February 7 to February 28, GRACE-C accelerometer in NRM, GRACE-D accelerometer inactive.

5. LRM3 = 28 February 2019 to 26 May 2020. GRACE-C accelerometer in NRM, GRACE-D accelerometer in LRM.

6. BVM = 26 May 2020 to 2 June 2020. GRACE-C accelerometer in NRM, GRACE-D accelerometer in BVM (a ground test mode that was used for diagnostic purposes).

7. NRM3 = 2 June 2020 to present. Both accelerometers in NRM.

GRACE-D accelerometer readings became much noisier after a mode transition a month into the mission with large thermal correlations and bias jumps. The noise (called RCYRPA = Roll/Cross/Yaw/Radial/Pitch/Along, since it hits all six acceleration measurements) has remained since, through a power cycle and multiple mode transitions. As shown by Harvey et al. (2022), the noise is largely common-mode in all directions, driven, in radial, along, and pitch, by two linear factors (see Equation 4). First-order corrections, which vary slightly from day-to-day, and noticeably over longer periods, but along the lines of: [Image Omitted. See PDF] [Image Omitted. See PDF]where $${\widehat{x}}_{\mathit{PITCH}}$$ is an estimate of pitch angular accelerations derived from a Kalman filtered/smoothed attitude based on star camera and IMU measurements (Goswami et al., 2021; Harvey et al., 2022; Harvey & Sakumura, 2019), reduce RCYRPA to 10⁻⁹ m/s/s scatter around an 10⁻⁸ m/s/s thermal, while a second-order combination with stable coefficients [Image Omitted. See PDF]cleans up what remains to within instrument requirements (Christophe et al., 2015). Note that more details on the first and second-order equations are provided in Harvey et al. (2022). Equation 4 provides a single clean accelerometer axis; however, since we require information from two axes (Radial, Along) at that level, the challenge is breaking the signal down to recover both.

Drag

GRACE-C and GRACE-D see a similar drag environment, separated by twenty-odd seconds and a few kilometers East-West. Limited evidence from minimizing a fit between GRACE-D data, before the accelerometer malfunctioned, or pitch-adjusted data (Equations 2 and 3) from after, with transplanted GRACE-C data suggest that a best-fit combination is given by [Image Omitted. See PDF] [Image Omitted. See PDF]where LagPeriod is the separation in seconds between C and D crossing the same line of latitude (in this context, a surrogate for attack angle or relative pitch between difference between the spacecraft) and the scale factors of 1.04 mostly represent an overall scale factor on the accelerometer, within instrument tolerances. Again, more recent data is currently leading us to adjust the scale factors in this decomposition, but for this analysis, this version suffices and is what is currently provided in the ACU product.

Note, that the fitted coefficient of 2.16 × 10⁻³ ⋅LagPeriod is comparable to the empirical 3.2° correction previously and independently derived for GRACE (Bandikova et al., 2019).

To get a better look at drag in isolation, we focus on the HFNRS signal in RAP (Equation 4), which ignores LFS (containing large thermals) and HFRS (which, under early mission conditions, includes a lot of SRP/albedo). A least squares fit minimizing the HFNRS differences between RAP (Equation 4) and a transplant of RAP from GRACE-C confirms the empirical fit. [Image Omitted. See PDF]

This combination (the difference between Equations 4 and 7) cancels out the bulk of drag from a common environment—although more recent data (March 2023) suggests that the fit to lower frequency signal (where higher atmospheric activity cause a noticeable increase in low frequency drag signal), a more accurate reflection of common environment response, is a bit larger (about half a percent). As a bonus, the same coefficients also significantly reduce SRP signal (and hence we see the repetitive signal—HFRS—from SRP mostly cancel out too).

Naturally, GRACE-C and GRACE-D do not see exactly the same atmospheric environment—they pass by at slightly different times, slightly different longitudes, and slightly different heights. In a low-drag environment, environmental differences are most visible at “bubbles” (so called by Gerard Kruizinga, JPL)—concentrations of denser atmosphere which usually show up near the poles. Bubbles were seen previously on the Challenging Mini-satellite Payload (Reigber, 2003) and GRACE.

Because drag primarily hits along-track, we might expect atmospheric fluctuations to affect along-track differences more than radial—the lift on the GFO satellites is non-negligible, but an order of magnitude lower than drag. Section 7 details how this differential drag is dealt with in the accelerometer data processing.

SRP and Albedo Parity

Acceleration contributions from light (SRP/albedo) are inherently easier to deal with than drag. SRP on GRACE-C and GRACE-D is repetitive and modelable, driven by a stable flux from the sun and spacecraft geometry. Albedo, with a dependence on Earth surface properties, might be tougher to model well, but is much smaller than SRP (by nearly an order of magnitude) to begin with, and cancels well between the spacecraft.

Our accelerometer treatment depends on a SRP parity decomposition. Here, with a simplified orbital model we show that SRP parity holds analytically.

For this SRP parity analysis, we assume an earth-centered, inertial cartesian coordinate system, where we assume both GRACE-FO satellites are orbiting in the x-y plane. We orient the x-y plane such that the sun-pointing direction vector $$\widehat{S}$$ is aligned with the x-axis, and is therefore expressed as the unit vector [Image Omitted. See PDF]where β is the angle between the orbital plane and the sun-pointing direction. We assume GRACE-C and GRACE-D are in identical circular orbits around the earth, where GRACE-D lags GRACE-C by a constant angle Δϕ so that the spacecraft pointing vectors from the origin to the spacecraft are given as a function of time t by [Image Omitted. See PDF] [Image Omitted. See PDF]where θ = θ(t) = ωt is the orbital angle of GRACE-C (ω is the orbital angular velocity). Note that θ = 0 corresponds to “orbital noon” and θ = π corresponds to “orbital midnight.” It is the parity of the differential SRP about θ = 0 that we are going to determine.

We assume that the orientation of the spacecraft is determined by a “nadir-pointing” condition (η = 0) or a “spacecraft-pointing” condition (η = 1), where “nadir-pointing” refers to the spacecraft radial axis pointing toward the center of the Earth and “spacecraft-pointing” refers to the spacecraft oriented to point at each other (introducing a difference between the spacecraft in angle of attack–or pitch angle–relative to the velocity direction). If we let $$\gamma =\eta \frac{{\Delta }\phi }{2}$$, then the spacecraft body fixed unit vectors, roughly in the radial and along-track directions are given as functions of time by [Image Omitted. See PDF] [Image Omitted. See PDF]and [Image Omitted. See PDF] [Image Omitted. See PDF]

We assume that the GRACE-C and GRACE-D have the same geometrical shape. Consider a flat surface element of each spacecraft U_i with unit normal $${\widehat{u}}_{0}=\widehat{u}(\theta =0,\eta =0)$$. Then the time evolution of that unit normal $$\widehat{u}(t)$$ on each spacecraft is given by [Image Omitted. See PDF] [Image Omitted. See PDF]where ξ(α) corresponds to rotation matrix about the z-axis by the angle α: [Image Omitted. See PDF]

We assume the radiation reflectivity and absorption coefficients for the surface element U_i are normalized such that the reflected SRP $${\vec{P}}_{R}$$ and absorption radiation pressure $${\vec{P}}_{A}$$ acting on the surface element are given by [Image Omitted. See PDF]and [Image Omitted. See PDF]

We define the differential SRP between the two spacecraft in the radial $$\widehat{r}$$ direction for element U_i to be [Image Omitted. See PDF]where the pressure and radial vector for GRACE-D is evaluated at advanced time t = t + Δt = t + Δϕ/ω so that the two spacecraft are at the same physical position when the pressure components are compared, since $${\widehat{T}}_{C}(t)={\widehat{T}}_{D}(t+{\Delta }t)$$. For the equation in the along-track, just replace $$\widehat{r}$$ with $$\widehat{l}$$.

As a simple specific example, we consider two elements with equal and opposite normals in the y-direction, U_y+ and U_y−: [Image Omitted. See PDF] [Image Omitted. See PDF]

The differential SRP, Equation 20, for the case of absorption pressure $${\vec{P}}_{A}$$ is [Image Omitted. See PDF]for the radial component, and [Image Omitted. See PDF]for the along-track component.

One can see from inspection the radial component of the differential SRP Equation 23 is antisymmetric in orbital angle θ, and that the along-track component of the differential SRP Equation 24 is symmetric in orbital angle θ.

It is algebraically laborious, though straightforward, to show that this radial/antisymmetric and along-track/symmetric result for both absorption SRP $${\vec{P}}_{A}$$ and reflective SRP $${\vec{P}}_{R}$$ holds in the general case, as long as each general surface element $${\widehat{u}}_{0}=\left[{u}_{x},{u}_{y},{u}_{z}\right]$$ corresponds to another u_y → −u_y surface element $${\widehat{u}}_{0}=\left[{u}_{x},-{u}_{y},{u}_{z}\right]$$.

The SRP signals themselves (without taking the differential) are demonstrably radial/symmetric and along-track/anti-symmetric under the simplifying assumptions applied here.

An Algorithm for SRP Differences

Concentrating on the period from December 2018 to September 2021—low outgassing, low atmosphere, negligible thruster leaks, we develop an algorithm to account for effects due to SRP and albedo.

Specifically, the process is:

1. Per day, break down the HFRS component of the RAP differential drag combination (the difference between Equations 4 and 7) into symmetric and anti-symmetric components. This can be done using a sinusoidal reconstruction, or Fourier transform.

2. Retain only components from twice-per-rev to twenty-times-per-rev. Twenty-times-per-rev is empirically chosen as the upper limit due to diminishing returns in the gravity field recovery.

3. Similarly, retain only the twice-per-rev to twenty-times-per-rev components of HFRS in pitch-adjusted-along and pitch-adjusted radial.

4. Over December 2018 to September 2021, least-squares fit scale factors on SRP in Along on GRACE-C and SRP in Radial on GRACE-C (model derived SRP data—for example, from a box-wing SRP model) to the difference between parity recoveries and the pitch-adjusted data. This step is computing a least squares scale factor of the differences between an SRP model and the difference between pitch-adjusted data and parity recoveries using RAP (this yields the portion of the SRP correction that does NOT directly obey the parity assumption due to suboptimal scaling when transplanting data from GRACE-C). We recover a factor of 1.1%.

5. Subtract 1.1% of Along and 1.1% of radial from differential drag RAP HFRS, and break down what remains into symmetric (Symm) and anti-symmetric (Anti) components, from once-per-rev to twenty-times-per-rev. Let SRPAlong = Symm, SRPRadial = Anti, plus, in each case, 1.1% of the box-wing model just mentioned.

This yields a correction to Equations 5 and 6, which now includes a correction for SRP/albedo, given by [Image Omitted. See PDF] [Image Omitted. See PDF]

Note that the terms SRPAlong and SRPRadial are the values provided in the ACS product. Additionally, note that we could pull low-degree per-revs directly from pitch-adjusted along and radial (smoothing out coefficients over several weeks)—excluding once-per-revs, for which thermals in the secondary RCYRPA noise dominate—but in the algorithm presented, we chose to use an SRP parity model, tuned to match pitch-adjusted along and radial. A typical recovery of the SRP difference using the algorithm above looks like Figure 2.

Figure 2. Estimated solar radiation pressure difference from transplant, on a low-drag day. Red line = along-track estimate, blue line = radial.

For comparison, using an SRP model (a simple box-wing model), instead of the GRACE-D data, with consistent transplant coefficients and scale assumptions appears in Figure 3.

Figure 3. Modeled solar radiation pressure difference from transplant. Red line = along-track estimate, blue line = radial.

The basic shapes are similar with abrupt changes at shadow entry/exit seen around 100 and 250°; the small night-time (orbit angle close to 180) slope present in the data but not the model might represent the most visible part of a thermally-driven once-per-rev accelerometer bias or an effect of the unmodeled thermal radiation pressure. As a point for future investigation, calibrating any thermal per-rev out might improve gravity estimation.

We estimate gravity fields with this correction (called ACH_SRP), and compared them to JPL's uncorrected, simple transplant fields (ACT), and fields generated with Graz's model-based TUG product (Behzadpour et al., 2021), which does well modeling SRP differences. Note that the ACH_SRP field is essentially a gravity field run with a product made up of ACU + ACS + AC0, in the nomenclature of Equation 1. We assessed recovery quality by comparing Delta Degree Variance (DDV) plots computed according to (Yuan, 2018) and recoveries of the C_3,0 harmonic term (as previously done in Behzadpour et al., 2021), for which there is an independent source of knowledge, based on Satellite-Laser-Ranging data described in, for instance, Lemoine et al. (2006), Loomis et al. (2020), Sun et al. (2023).

In this period, both Graz's model and our data driven recovery noticeably improve C_3,0 estimate consistency with Satellite Laser Ranging (SLR) (Figure 4 Table 1) as compared to a simple transplant, with a slight edge to the new algorithm. The effect on C_3,0, and other zonal coefficients, is substantial due to the repetitive nature of the correction with orbit angle.

Figure 4. C3,0 recoveries running from simple transplant (ACT), Graz accelerometer product (TUG), transplant with solar radiation pressure from parity (ACH_SRP). TN-14 is the recovery from Satellite Laser Ranging.

Table 1 RMS Difference From SLR-Based Estimate of C_3,0, by Average Monthly Beta Angle, in mm

	ACT	ACH_SRP	TUG
0° < |β| < 30°	0.74	0.14	0.16
30° < |β| < 60°	0.27	0.14	0.26
60° < |β| < 90°	0.19	0.21	0.22

Figure 5. Median Delta Degree Variance relative to GGM05C, December 2018 to September 2021, for high beta months (45° [less than] |β| [less than] 90°). Simple transplant (ACT), Graz accelerometer product (TUG), transplant with solar radiation pressure from parity (ACH_SRP).

Figure 6. Monthly gravity field Delta Degree Variance relative to GGM05C, December 2018, a high beta month (β = −46°). Simple transplant (ACT), transplant with solar radiation pressure from parity (ACH_SRP).

Figure 7. Monthly gravity field maps for December 2018, relative to GGM05C, a high beta month (β = −46°). Simple transplant (ACT), transplant with solar radiation pressure from parity (ACH_SRP).

Figure 8. Monthly gravity field Delta Degree Variance relative to GGM05C, April 2019, a low beta month (β = 8°). Simple transplant (ACT), transplant with solar radiation pressure from parity (ACH_SRP).

Figure 9. Monthly gravity field maps for April 2019, relative to GGM05C, a low beta month (β = 8°). Simple transplant (ACT), transplant with solar radiation pressure from parity (ACH_SRP).

Additionally, both SRP recovery strategies improve DDV, with the effect concentrated, as expected, at low beta angles (when the orbital plane is in line with the Earth-Sun line), where SRP variations are more significant (Figures 5–9).

For the future, a hybrid strategy, which tunes a model from empirical data, should work better than either. This will be particularly true as drag signal increases, since the algorithm presented assumes all of the HFRS is driven primarily by SRP/albedo; if drag begins to significantly contribute to HFRS, modifications will be required.

An Algorithm for Drag and Thruster Differences

In recent months—for this note, October 2021 to December 2022—drag has increased, as have thruster leaks. To simplify our analysis, we did not materially change our SRP analysis strategy when looking at this data (any recent parameter updates are not included), so solution quality is not optimized, but this analysis suffices to demonstrate that accelerometer data from GRACE-D can materially improve gravity field recovery in the presence of drag and leaky thruster differences. For differential drag and leaky thrusters:

1. Start with the HFNRS in RAP differential drag (difference between Equations 4 and 7). Call this time series δRAP.

2. Using the GRACE-C ACC1A data, compute the root mean square of jumps at all thruster firings. The jumps amongst all thruster firings on GRACE-C is used as a metric for how much variability is present in the thruster leak, called DMLTS (daily measured leaky thruster size). This metric assumes that the variability in the GRACE-C thruster leak is similar to that of GRACE-D.

3. Using the spline smoothed time series above, compute the running maximum of the time series over a 5 min window (300 s), $${\epsilon}\left({t}_{j}\right)=\max \left(\left\{\delta RAP\left({t}_{j}-{t}_{i}\right):{t}_{i}=-150,\text{\ldots },150\right\}\right)$$.

4. Compare $${\epsilon}\left({t}_{j}\right)$$ with DMLTS at each point in the running window to compute SOAD (size of atmospheric difference), $$\text{SOAD}\left({t}_{j}\right)=\sqrt{\max \left(\left\{{\epsilon}{\left({t}_{j}\right)}^{2}-DMLT{S}^{2},0\right\}\right)}$$.

5. Compute a weight based on SOAD and DMLTS, $$W\left({t}_{j}\right)=\frac{SOA{D}^{2}}{SOA{D}^{2}+DMLT{S}^{2}}$$.

6. Allocate δRAP to along and radial based on the weight as:

[Image Omitted. See PDF] [Image Omitted. See PDF]

Note that Equations 27 and 28 are the values provided in the ACM product.

Because we excluded the repetitive signal from processing here, adding this drag/thruster leak correction doesn't change C_3,0 much, in comparison to just using the SRP correction (Figure 10); we call this combined correction for SRP and drag ACH_ALL. Note that ACH_ALL is the full ACH product.

Figure 10. C3,0 recoveries running from simple transplant (ACT), Graz accelerometer product (TUG), transplant with solar radiation pressure (SRP) from parity (ACH_SRP), transplant with SRP and drag/thruster leak correction (ACH_ALL). TN-14 is the recovery from Satellite Laser Ranging.

While the C_3,0 results are noticeably closer to SLR than the simple transplant, and on par or slightly better than results processing with TUG data, they are not as good as during the earlier, low drag portion of the mission, due to limitations highlighted in the SRP algorithm description. In the future, applying an empirical SRP model, tuned to early-mission data, instead of recovering SRP locally in a high-drag environment, might help.

DDV plots, however, improve materially and systematically at all beta angles (Figures 11–13). The difference is clearly due to drag/leaky thruster correction from RAP; fields based on ACH_SRP or TUG don't improve nearly as much.

Figure 11. Median Delta Degree Variance, October 2021 to December 2022. Simple transplant (ACT), Graz accelerometer product (TUG), transplant with solar radiation pressure (SRP) from parity (ACH_SRP), transplant with SRP and drag/thruster leak correction (ACH_ALL).

Figure 12. Monthly gravity field Delta Degree Variance, May 2022. Simple transplant (ACT), transplant with solar radiation pressure (SRP) from parity (ACH_SRP), transplant with SRP and drag/thruster leak correction (ACH_ALL).

Figure 13. Monthly gravity field maps for May 2022. Simple transplant (ACT), transplant with solar radiation pressure (SRP) from parity (ACH_SRP), transplant with SRP and drag/thruster leak correction (ACH_ALL).

Localizing Effect of Error

While in principle there is no difference between range-rate and range-acceleration processing, in practice, the natural observable content in range-acceleration tend to localize acceleration error, by not integrating it through indefinitely. This can be advantageous when processing acceleration data contaminated by significant error sources.

In the GRACE-FO scenario, with prominent accelerometer errors, localizing the effect of these errors can provide a critical improvement, as described by Allgeyer et al. (2022). This kind of change has become increasingly important; we ran a year of data, in the recent high drag, high thruster leak period, and found a substantial and systematic improvement from range acceleration processing; for convenience, using LRI range rate and range acceleration data (Figure 14. Note that KBR could be used to provide identical results and the result is independent of KBR or LRI data).

Figure 14. Median Delta Degree Variance, LRI range acceleration (blue) versus LRI range rate (green) processing, late 2021 to late 2022.

First Steps in Loose Attitude Control Work

More recently (January/February 2023) the GRACE-FO mission has experimented with a looser attitude control (wide deadband mode), to cut down on thruster use and mitigate the leaking thrusters. Further operations in this mode are planned. While we do not yet have gravity field recovery results fully adapting our algorithms to this mode, preliminary results suggest that similar strategies, incorporating the GRACE-D data, can help there as well. Luckily, the pitch-adjusted Radial and Along data is particularly well behaved in this period—good enough that an accelerometer product based on this first-order corrected data suffices for decent science and is currently provided in the released ACH data product, only for months in wide deadband mode.

The pitch-adjusted Radial and Along data is also good enough for model tuning. Combined with high drag and large deviations in roll and yaw (and smaller, but not negligible, deviations in pitch), we can see that allowing drag coefficients to vary as a function of yaw-squared in Along, and a function of pitch and yaw-squared in radial, where pitch and yaw are defined with respect to the local flow of the atmosphere (including winds) relative to the spacecraft, accounts for most of the observed variation. [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where [Image Omitted. See PDF]and SRPModelAlongC is simply an SRP/albedo, which is removed from the raw GRACE-C accelerometer data. Note that an SRP/albedo model, instead of the algorithm described in Section 6, is used here because the atmospheric drag is large and the assumptions described in Section 5 are no longer valid in wide deadband mode; making adequate modeling of all forces other than drag (or tuning models using the GRACE-D accelerometer data) extremely invaluable. Finally, the HFNRS from differential RAP (difference between Equations 4 and 7) is utilized differently than presented in Section 7 by applying the entire time series to the along-track accelerations (denoted in Equation 29 by the δRAP term). The δRAP term is utilized in this manner based on the assumption that the wide deadband mode has mitigated the effect of the leaky thrusters and we now only need to account for drag effects.

We have begun estimating experimental gravity fields based on a pure empirical model for drag coefficients and local atmospheric flow recovery, derived from this data, with very encouraging results, shown in Figures 15–18.

Figure 15. Monthly gravity field Delta Degree Variance, January 2023. Nominal, publicly released ACH product (RL06.1), updated experimental ACH product for wide deadband mode (ACH_WDB).

Figure 16. Monthly gravity field maps for January 2023. Nominal, publicly released ACH product (RL06.1), updated experimental ACH product for wide deadband mode (ACH_WDB).

Figure 17. Monthly gravity field Delta Degree Variance, February 2023. Nominal, publicly released ACH product (RL06.1), updated experimental ACH product for wide deadband mode (ACH_WDB).

Figure 18. Monthly gravity field maps for February 2023. Nominal, publicly released ACH product (RL06.1), updated experimental ACH product for wide deadband mode (ACH_WDB).

The dramatic improvement at higher degrees, in the degree variance curves, and significantly reduced striping, in the global maps, are readily apparent.

We've also used the data from early 2023 to improve our parameter tuning in standard, tight attitude control—in the future, we will produce revised products that incorporate these improvements in combination with improved SRP/albedo modeling.

Conclusions and Future Work

Using data from the GRACE-D RAP axis substantially reduces error in GRACE-FO gravity field recoveries, as does a switch in observable to localize the effect of accelerometer error. It has been shown that the GRACE-D RAP data can be utilized to improve estimates of the C_3,0 coefficient and across the gravity field spectrum, indicating that there is significant signal still present underneath the noise in the GRACE-D accelerometer. Additionally, the GRACE-D accelerometer data has been shown to be a vital component in tuning empirical models. The results presented in this note for drag, SRP, and leaky thruster estimation can all be further extended. Better models or algorithms for true pitch recovery, thermals in RAP, phantom outlier removal, SRP, thruster leaks, drag coefficients, local wind recovery, or thermal emission pressure from the spacecraft would all improve solution quality further than is demonstrated here.

Acknowledgments

Research by Bertiger, Harvey, McCullough, Miller, and Yuan was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Research by Save was sponsored by JPL contract 1604489. On the GRACE-FO team, Vishala Arya, Srinivas Bettadpur, Eugene Fahnestock, Robert Gaston, Zhigui Kang, Felix Landerer, Peter Nagel, Meegyeong Paik, Chirag Patel, Athina Peidou, Nadege Pie, Steven Poole, Manoucher Shirbacheh, and Mark Tamisiea supported this investigation. Additionally, communication with Airbus Defense & Space, particularly conversations with Thomas Usbeck, has contributed significantly to our understanding of the thruster leak and thruster nozzle design. Finally, communication with ONERA, particularly Bruno Christophe, has been invaluable in our understanding of the accelerometer instrument.

Data Availability Statement

The Level-1 data used in this analysis is publicly available from PO. DAAC (https://podaac.jpl.nasa.gov/GRACE-FO?sections=data). Specifically, ACC1A data (raw acclerometer data) is available within the bundles at https://podaac.jpl.nasa.gov/dataset/GRACEFO_L1A_ASCII_GRAV_JPL_RL04 ACT1B data (calibrated accelerometer data for GRACE-C) is available within the noLRI bundles at https://podaac.jpl.nasa.gov/dataset/GRACEFO_L1B_ASCII_GRAV_JPL_RL04 and calibrated ACH1B data (calibrated accelerometer data for GRACE-D and its components) is available within the ACX bundles at https://podaac.jpl.nasa.gov/dataset/GRACEFO_L1B_ASCII_GRAV_JPL_RL04.