GRACE Follow‐On Accelerometer Data Recovery

In Gravity Recovery and Climate Experiment (GRACE) Follow‐on (GRACE‐FO) mission, similar to its predecessor GRACE, the twin satellites are equipped with three‐axis accelerometers, measuring the non‐gravitational forces. After 1 month in orbit, during the in‐orbit‐checkout phase, the noise on GRACE‐D accelerometer measurements elevated and resulted in systematical degradation of the data. For this reason, the GRACE‐D data need to be replaced by synthetic data, the so‐called transplant data, officially generated by the GRACE‐FO Science Data System (SDS). The SDS transplant data are derived from the GRACE‐C accelerometer measurements, by applying time and attitude corrections. Furthermore, model‐based residual accelerations due to thruster firings on GRACE‐D were added, proven to improve the data quality in gravity field recovery. However, preliminary studies of GRACE‐FO data during the single accelerometer months show that the low degree zonal harmonics, in particular C20 and C30, are sensitive to the current transplant approach. In this work, we present a novel approach to recover the GRACE‐D ACT1B data by incorporating non‐gravitational force models and analyze its impact on monthly gravity field solutions. The results show the improved ACT1B data not only contributed to a noise reduction but also improved the estimates of the C20 and C30 coefficients. The application of this new approach demonstrates that the offset between Satellite Laser Ranging (SLR) and GRACE‐FO derived C30 time series can be reduced by the use of the alternative accelerometer product.

sible because both satellites fly in the same orbit and have a time delay of 25-30 s. Therefore, the change in the non-gravitational accelerations during this time delay is very small, and both accelerometers would measure approximately the same signal. However, it is mandatory to apply time and attitude corrections to the available accelerometer measurements, to account for the variable separation between the two spacecrafts and the orientation differences relative to each other.
Over the final months of GRACE, after the GRACE-B accelerometer was turned off due to battery issues, the demand for synthetic ACC data raised once again. Bandikova et al. (2019) presented an improved ACC data transplant approach, which includes the modeling of residual linear accelerations due to thruster firings, in addition to the attitude and time correction. The modeling of thruster spikes has been studied before both on ACC1B (Meyer et al., 2011) and ACC Level-1A (ACC1A) (McCullough et al., 2015) using a statistical approach. However, Bandikova et al. (2019) show that by determining the dynamic system, which generates the impulse response of each thruster pair, one can improve the estimates of the spikes and reduce the noise in gravity field solutions. The improved ACC transplant contributes to the processing standards for final months of GRACE, and ACC transplant data were included in the official Level-1B products.
Based on the official Level-1B products, GRACE/GRACE-FO Level-2 RL06 products are provided by the SDS members. To provide long-term gravity field time series, the processing standards for GRACE-FO remained mostly consistent with GRACE. In addition to the official centers, various other research institutions are also preparing their processing setup for GRACE-FO gravity field solutions. Among these centers, the Institute of Geodesy at Graz University of Technology (TUG) routinely provides operational GRACE-FO products, consistent with ITSG-Grace2018 standards (Kvas et al., 2019), which are published by the International Center for Global Earth Models (ICGEM, Drewes et al., 2016).
Preliminary analysis of GRACE-FO data during the single accelerometer months proves that the derived low degree zonal harmonics, in particular C 20 and C 30 , are affected by the transplant approach, and, as recommended by SDS, are needed to be replaced by SLR-derived values (Loomis et al., 2020). Figure 1 shows the absolute difference between GRACE-FO-derived C 30 and SLR-derived time-series, provided in the Technical , and compares their variations with the β′ angle of the orbit. The β′ angle represents the angle between the orbital plane of the satellite and the Sun direction. It is visible that the peaks, indicating the largest deviations from the recommended values, are highly correlated with the periods with β′ angle close to zero. During these periods, the satellite is exposed to direct sunlight for approximately half of its revolution, and the rest of the time is passing through the Earth's shadow, causing large temperature differences. This orbital condition with 161-day cycle has an effect on the satellite's operation, specifically on the energy absorption by solar panels and on the thermal control due to the heating by the Sun. To further investigate the reason, the transplant approach can be evaluated in GRACE mid-timespan, when both accelerometers operated nominally and provided high-quality data. Meyer et al. (2011) compared the acceleration signal measured by each satellite in the same orbit position during a half-cycle of β′ angle between March and September 2007. They showed that a difference up to 3 nm/s 2 exists between the two measured signals. The most substantial difference happens for a short period of time when satellites are transiting through Earth's shadow, as the solar radiation pressure (SRP) changes for each satellite at the same position. Additionally, a long-term variation exists mainly during periods when β′ angle is close to zero. In other words, when satellites are under direct Sun exposure, due to high fluctuations in the atmosphere density, they experience different forces, which in this case is the drag force. Hence, the direct transplanting of ACC data can cause large errors under these conditions.
The main purpose of this paper is to present an improved method to recover the GRACE-D accelerometer data by incorporating non-gravitational force models and analyze its impact on the recovered gravity field solutions. Within the gravity field recovery, the transplanted GRACE-D accelerometer data show a significant impact on the low degree coefficients, in particular reduces the offset between SLR and GRACE-FO derived C 30 time series.
The paper is organized as follows: In Section 2, we provide a brief overview about the accelerometer data and focus on different ACC data products. In Section 3, we briefly describe the non-gravitational forces acting on the satellites. In Section 4, we present the alternative GRACE-D ACC recovery. And finally, in Section 5, we discuss the results by comparing our products with the official transplant data and their impact on the recovered gravity field solutions.

Accelerometer Data
The GRACE-FO accelerometer is a three-axis electrostatic accelerometer manufactured by the Office National d'Études et de Recherches Aérospatiales (ONERA) (Christophe et al., 2015). The accelerometer provides information about the linear and angular acceleration of the satellite. The sensor has two high-sensitive axes, the radial and along-track axes exhibiting a resolution better than 0.1 nm/( 2 s Hz ), and one less-sensitive axis, the cross-track axis with a precision of 1 nm/( 2 s Hz ).
The accelerometer is placed in the center of mass (CoM) of the satellite. The core of the sensor consists of a proof mass, surrounded by an electrode cage. The proof mass is suspended by electrostatic forces generated by the electrodes. The sensor measurement is determined from the usage of analog voltages, producing the electrostatic force. The electrostatic force is proportional to the sum of the non-gravitational forces acting on the satellite, including atmospheric drag, solar radiation pressure (SRP), earth radiation pressure (ERP), and other disturbances.
There are three types of science data products for the GRACE-FO accelerometer data, which in the following, will be briefly introduced.

ACC1A
The ACC Level-1A (ACC1A) data products include the 10-Hz linear acceleration measurements, given in accelerometer reference frame (AF). The origin of this frame is defined to be the center of mass of the proof mass and the frame axes are directed as shown by Figure 2. The time frame of the data is determined by the on-board computer (OBC).

ACT1A
To obtain better gravity field solutions, a series of calibration processes is applied to ACC1A, providing the calibrated accelerometer data (

ACT1B
During Level-1A to Level-1B data processing, the ACT1A linear accelerations are edited, time tagged, and low-pass filtered in order to reduce the high-frequent measurement noise (Wu et al., 2006). Figure 3 BEHZADPOUR ET AL.
10.1029/2020JB021297 4 of 17  describes the Level-1A to Level-1B processing for linear accelerations. The process is discussed in more details in Appendix A. The output of this process is 1-Hz ACT1B product, given in satellite reference frame (SRF) and GPS time.
As the unprocessed Level-1A data products of GRACE-FO are publicly available to all processing centers, an independent Level-1A to Level-1B data processing can be set up. At TUG, we have implemented standard algorithms for data screening, time synchronization and data rate reduction for the raw data, which produce in-house (i.e., independent of the officially released product) ACT1B data, serving as input for the gravity field recovery. To validate the implementation, we compare the TUG ACT1B from GRACE-C with the SDS product. Figure 4 shows time-series and power spectral densities of differences between the two products.

Table 1
Non-Gravitational Force Models due to Radiation Pressure Figure 5. Comparison between calibrated GRACE-C measurements (black) and simulated data on along-track (blue), cross-track (red), and radial (green) direction on December 05, 2019 (left) and September 12, 2019 (right).
The deviations are more than one order of magnitude smaller than the signal noise on each axis, which demonstrates the preciseness of our implementation.

Non-Gravitational Force Models
GRACE-FO is operating in a low-Earth orbit (LEO). The non-gravitational accelerations acting upon LEO satellites are mainly due to atmospheric drag, SRP, and ERP. The radiation pressure force models are summarized in Table 1. For the details on the formulae applied for the radiation pressure modeling, the reader is referred to the corresponding references. To model these effects, the satellites' geometry and surface properties are obtained from GRACE-FO macro model (Wen et al., 2019). As the drag force model plays an important role in recovering the GRACE-D missing measurements, the basic equations of respective accelerations are described in the following.

Atmospheric Drag
Aerodynamic force is the force acting on the satellite's surface caused by interchange of momentum with the atmosphere molecules. For LEO satellites, it is the dominant non-gravitational perturbation. The aerodynamic force is modeled as: depending on the aerodynamic coefficient C a , the atmospheric density ρ, the cross-sectional area A ref , the satellite mass m, and the true airspeed v TAS , that is, the velocity of the spacecraft relative to the atmosphere. Here, the cross-sectional area is obtained from the macro model, and the mass of the satellite is set to their launch mass values, as the MAS1B data product currently report the tank mass, and not the total mass. Therefore, for the purpose of simplification, this variation is being neglected.
For atmospheric density, there exist several models, such as the Jacchia-Bowman 2008 (JB2008; Bowman et al., 2008), and the Drag Temperature Model 2013 (DTM2013; Bruinsma, 2015), and the NRLMSISE-00 (Picone et al., 2002). Our initial analysis based on the DTM2013 and JB2008 models revealed that using either of these models leads to comparable results, as in our approach the density variations are corrected afterward based on GRACE-C real measurements. Nevertheless, due to the short delivery latency and the availability of real-time values, we obtain atmospheric density from the JB2008 model.
The v TAS is the relative velocity of the satellite with respect to the atmosphere, which is the sum of inertial velocity of the satellite  r, co-rotating atmosphere and atmospheric wind velocity v w :  where r is the satellite position and ω E is the angular velocity of the Earth sidereal rotation. The wind velocity is derived from Horizontal Wind Model 2014 (HWM14; Drob et al., 2015). The component of the aerodynamic force toward normal velocity direction ˆT AS v is referred to as drag and the component toward as lift, with n being the unit normal vector to the satellite plate. Hence, the aerodynamic coefficient can be expressed by: Therein, C D and C L are dimensionless drag and lift coefficients, respectively. As drag is the major component of the aerodynamic force acting on satellites, neglecting lift and referring to drag force instead of aerodynamic force are conventional. Consequently, the aerodynamic coefficient can be referred to as drag coefficient, which is set to a constant value of 2.4 for GRACE-FO. Figure 5 shows the summation of non-gravitational models for 2 days of mission at different β′ angles. This comparison reveals that at similar altitude, the maximum acceleration and deviation from model occur when satellite passes through the day-side under β′ ∼ 0° condition.
Due to uncertainties in the state and attitude of the satellite, interaction of the satellite's surface and atmosphere molecules affecting the drag coefficient, as well as uncertainties associated with atmospheric density models, it is not possible to model the drag force accurately. Therefore, the values used for density and drag coefficient play a significant role in drag model uncertainty (Moe & Moe, 2005;Prieto et al., 2014).
According to Equation 1, drag coefficient and density errors are not completely separable without a good knowledge of at least one of them. Hence, the model exhibits a multiplicative total error from these two sources. However, in case of co-orbiting missions, it is convenient to assume that both satellites experience the same environment, in which the drag model error for each satellite is approximately equivalent. Therefore, by having one actual accelerometer measurements, in this mission GRACE-C, one can estimate the model error, which then can be used to retrieve the missing measurements from the other spacecraft with a malfunctioning accelerometer, that is, GRACE-D.

GRACE-D ACC Recovery
In the following, the recovery processing for linear accelerations from GRACE-C ACT1A to GRACE-D ACT1B is described. Table 2 summarizes the required data products in this procedure. The input data of the process are the daily 10 Hz GRACE-C ACT1A. In the first step, we apply the standard Level-1A processing, as described in Appendix A, to convert the data to GPS time and SRF frame. Additionally, before applying the low-pass filter, we remove the thruster spikes using THR1B products with a margin of 1 s and fill the gaps with linear interpolation.

Calibration and Model Reduction
The next step is computing the simulated accelerations according to Section 3 for both satellites with 1-Hz sampling. Since orbit data and SCA1B are already created with the same sampling, no interpolation is needed. Using the simulated data, we calibrate the GRACE-C data, obtained from the first step, with daily constant bias on each axis. Lastly, we reduce the simulated data a model from the calibrated GRACE-C data a cal , creating the unmodeled acceleration signal Δa:

Time Correction
In this step, we transfer required GRACE-C Level-1B data, including orbit, SCA1B, and the unmodeled accelerations Δa to GRACE-D time frame by applying a transfer time correction. The input data for com-BEHZADPOUR ET AL.

10.1029/2020JB021297
8 of 17 puting this correction are linear interpolated GNV1B data of both satellites, which delivers the orbit positions and velocities in GPS time frame. By finding the distance between two satellites at each epoch and GRACE-D velocities, the transfer time correction can be estimated by the distance traveled divided by the velocity. Figure 6 shows the transfer time correction obtained for 01 July 2018.

Drag Model Correction
As mentioned in Section 3, the un-modeled acceleration signal contains an unknown multiplicative error S D = ρ ⋅ C D , that is, a scale factor, related to atmospheric density and aerodynamic coefficient. Within the recovery process, we estimate this scale factor in a least squares sense: daily interval (k = 1,440). Therefore, the total number of k + d = 1,443 parameters per interval are needed for the scale estimation. The outputs of this adjustment are: (a) the calculated scale, which will be applied to GRACE-D drag model, and (b) the residual signal, containing errors from other sources as well as the estimation process.

Attitude Correction
The residual signal needs to be transferred to GRACE-D frame by applying attitude correction, as the orientation of each satellite with respect to its velocity vector is different. The K-band ranging (KBR) measurement principle requires precise alignment of each satellite's KBR antenna toward each other, that is, in the direction of the so-called line of sight (LOS). Since the antenna is mounted at the front panel of each BEHZADPOUR ET AL.

10.1029/2020JB021297
10 of 17 satellite, the leading satellite is turned by 180° around its z-axis. Furthermore, both satellites fly with a pitch offset of approximately 1° with respect to the LOS.
According to Bandikova et al. (2019), for GRACE accelerometer transplant, it was not possible to correct the attitude with the actual attitude variations, due to high-frequency noise on the star camera measurement. Therefore, the attitude correction was approximated by a 180° yaw and an empirical (not physical) 3.2° pitch transform.
In GRACE-FO, the attitude sensors on-board each satellite consist of three star cameras and one angular rate sensing inertial measurement unit (IMU). These attitude data are combined by means of a Kalman filter to obtain an optimal attitude product (SCA1B) (Harvey & Sakumura, 2019). According to , this method reduces the noise level by approximately two orders of magnitude lower than GRACE. This allows us to directly use the GRACE-FO SCA1B to transfer residual accelerations from GRACE-C to GRACE-D frame.

Thruster Spikes
The final step in the recovery process is to add GRACE-D modeled thruster responses due to attitude thruster firings. This has been done for GRACE-D ACT1A based on regression of the available accelerometer data, as reported by the difference between two processed GRACE-D ACT1B time-series: one with modeled spikes and the other with removed spikes. For the latter, we remove the thruster events with a margin of 1 s and fill the gaps with linear interpolation during Level-1A processing. The resulting data only includes roll/pitch/yaw thruster firing spikes with 1 s sampling.

Results
In the following, we present the comparison of in-house and official ACT1B data sets in time and frequency domain as well as gravity field recovery based on these products. Figure 8 shows time-series and power spectral densities of differences of the official transplanted data, SDS ACT1B, and the TUG ACT1B. The differences between two time-series under indirect sunlight condition (β′ = 71.8°) are insignificant. Therefore, under this condition, our approach independently validates the official accelerometer products. A major difference is visible in radial component in frequency domain, in particular in the orbital frequency, which is approximately one cycle per 89 min or 1 cycle per revolution (cpr) frequency. As expected, higher frequencies are less affected by the recovery process, as we used the same thruster responses, which are dominant at frequencies over 3 mHz. Figure 9 reveals that the magni-BEHZADPOUR ET AL.

10.1029/2020JB021297
12 of 17 tude of the differences between the SDS and TUG data is dependent on orbit configuration with respect to Sun, that is, β′ angle. Here, for all components, major differences exist at 1-3 cpr frequencies. For further details, the differences are plotted as a function of argument of latitude and time in Figure 10. It is visible that the recovery process mainly affects radial and along-track components, when satellites are directly illuminated by the Sun (β′ ∼ 0°).

Drag Model Scale
As depicted in Section 4, we estimate drag model scale factors during the GRACE-D ACC recovery process. Figure 11 illustrates the time series of the estimated scale. The estimates show scattered variations comparable with β′ angle cycle. This behavior reflects the higher uncertainties in atmosphere density during periods with direct sunlight, that is, with β′ angle ranging from −20° to 30°. To some extent, the variations of the estimated values also depend on the temperature fluctuations in these periods, which affects the gas molecular behavior and therefore, the drag coefficient.
Additionally, the high-frequency perturbations of the atmospheric density during geomagnetic storms are also expected to be absorbed by the scale factors. However, our observations are mainly collected during solar minimum (2018-2020) and, therefore, sparse solar events are visible in the data. During this period, the severest geomagnetic storm happened on 26 August 2018, whose impact on GRACE-C accelerometer measurements has been studied by Krauss et al. (2020). Due to the gaps in GRACE-D data, this period is not included in the scale time series. Nevertheless, few minor geomagnetic storms have triggered atmospheric disturbances during 2019 and 2020 and subsequently affected the scale factors.
Among a variety of indices, which characterize the geomagnetic activities, it has been shown that the SYM-H index (Iyemori et al., 2010) has one of the highest correlation with the neutral atmospheric density, with nearly zero time delay at low earth orbiters (Krauss et al., 2015). SYM-H index describes the geomagnetic disturbances at mid-latitudes with a temporal resolution of 1 min. Figure 12 shows the time series of the estimated scale, compared to the variations of this index. For two selected events, a clear correlation exists between peaks in SYM-H and scale factor time series.
It is worth mentioning that in addition to the described effects, the estimates probably absorb other effects, which are not adequately modeled (e.g., variations in radiation pressure models).

Gravity Field
This section compares the monthly gravity field solutions based on the recovered GRACE-D accelerometer data and ITSG-Grace2018 scheme (Kvas et al., 2019). We recovered monthly solutions from July 2018 to August 2020 using (a) SDS ACT1B in preliminary results, denoted as ITSG-Grace2018p, and (b) TUG ACT1B in final GRACE-FO operational solutions, referred to as ITSG-Grace2018. The ITSG-Grace2018 solutions are publicly available and can be downloaded from ICGEM site (http://icgem.gfz-potsdam.de/or ifg.tugraz.at/ITSG-Grace2018). Figure 13 shows the degree variances for July 2020 (β′ = −0.4°) with respect to the static gravity field GO-CO06s (Kvas et al., 2021). Compared to the preliminary solution, ITSG-Grace2018 with in-house computed ACT1B performs better in terms of noise over all spherical harmonic degrees, with the largest differences in degrees 2 and 3. Figure 14 shows the power spectral densities of the post-fit range rate residuals from both solutions (ITSG-Grace2018p and ITSG-Grace2018). The post-fit residuals are the differences between adjusted and original range rate observations. The residuals indicate the total error sources, which are present within the gravity recovery process, as they are propagated errors of all involved instruments and modeling errors.
Here, strong 2 and 3 cpr signals, existing in the preliminary scenario, are mitigated using the GRACE-D TUG ACT1B. This clearly demonstrates the overall improvement of gravity field solutions and better estimation of low degree coefficients, particularly C 20 and C 30 .
For more details, we compare the monthly estimates of C 20 derived from GRACE-FO and SLR-derived values, provided in the TN-14, over the entire available period. Figure 15a shows that the bias in the estimates from ITSG-Grace2018 solutions has reduced with respect to the TN-14 values. However, similar to GRACE, the C 20 coefficients still exhibit a strong 161-day periodic signal and should be replaced by C 20 estimates from SLR data.
Additionally, Figure 15b demonstrates the C 30 time-series derived from GRACE-FO and SLR. The ITSG-Grace2018 estimates for the C 30 coefficients show a higher correlation with the SLR solution. This supports the hypothesis that by mitigating 2 and 3 cpr signals in transplanted data, resulting from alternative GRACE-D ACT recovery, the estimates of the C 30 coefficients are significantly improved. The improvement is mainly visible for the months with β′ ∼ 0°.

Conclusions
The results presented in this paper show the advantages of incorporating non-gravitational force models and applying drag model corrections within GRACE-D ACC recovery. The alternative ACT product contributes to an improved estimation of higher degrees of the recovered monthly gravity field solutions, as well as low zonal degrees 2 and 3. The estimates of the C 30 coefficients represents a significant improvement with respect to SLR-derived values, with which the GRACE-FO values are recommended to be substituted.
Additionally, we were able to prove that the non-gravitational forces, in particular the corrected drag model, recover part of the acceleration signal that is missing in the direct transplant approach in radial and alongtrack directions. The magnitude of this missing signal is β′ angle dependent and reaches its maximum during direct Sun exposure. Within the periods of passing the Earth's shadow, the modeled signal is approximately zero.
Furthermore, we show that the drag model scale correction mainly reflects the physical characterizations of changing thermospheric density. This includes long-term variations during periods with direct sunlight, as well as short-term fluctuations due to geomagnetic storms.
However, there are several issues, which have a potential for further improvements: in order to avoid the effects from outliers in GRACE-C measurements, the parameterization of the drag scale reflects smoothed variations longer than 1 min. This setup may result in losing high-frequency details in the satellites' changing environment. Moreover, thruster spikes, which also contributes to the high-frequency spectrum of the measurements, are currently modeled with simple impulse functions with constant magnitudes based on statistical estimation. As proved in GRACE mission , incorporating more realistic thruster responses in transplant data considerably reduces noise in degrees over 15, that is, the first orbital resonance. It is also worth mentioning that the alternative transplant data are based on the bias-corrected GRACE-C measurements; therefore, it still exhibits unknown instrument specific scales. To avoid a degra-BEHZADPOUR ET AL.

10.1029/2020JB021297
15 of 17 dation of the recovered monthly gravity field solutions, the accelerometer scale and bias need to be modeled and co-estimated during gravity field recovery.
In summary, this work contributes to an improvement of the GRACE-FO derived gravity field solutions by investigating the role of the missing GRACE-D accelerometer data. The efforts in constructing the constituents of the missing measurements not only increase the gravity field quality but also improve our knowledge about satellites' interaction with their environment.