Calibration of GRACE on-board accelerometers for thermosphere density derivation

ABSTRACT Low Earth Orbit satellite on-board accelerometers play an important role in improving our understanding of thermosphere density; however, the accelerometer-derived densities are subject to accelerometer calibration errors. In this study, two different dynamic calibration schemes, the accelerometer parameter-incorporated orbit fitting and precise orbit determination (POD), are investigated with the Gravity Recovery And Climate Experiment (GRACE) satellite accelerometers for thermosphere density derivation during years 2004–2007 (inclusive). We show that the GRACE accelerometer parametrization can be optimized by fixing scale coefficients and estimating biases every 60 min so that the orbit fitting and POD precision can be improved from 10 cm to 2 cm in the absence of empirical acceleration compensations and as a result the integrity of calibration parameters may be reserved. The orbit-fitting scheme demonstrates similar calibration precision with respect to POD. Their bias estimates in the along-track and cross-track components exhibit an offset within 0.1% and a standard deviation (STD) less than 0.3%. Correspondingly, a bias of 2.20% and a STD of 5.75% exists between their thermosphere density estimates. The orbit-fitting and POD-derived thermosphere densities are validated through the comparison against the results published by other institution. The comparison shows that either of them can achieve a precision level at 6%. To derive thermosphere density from the rapid-increasing amount of on-board accelerometer data sets, it is suggested to take full advantage of the orbit-fitting scheme due to its high efficiency as well as high precision.


Introduction
Accelerometers onboard the Low Earth Orbit (LEO) satellites, through providing high accuracy measurements of the nonconservative accelerations acting on the satellites, have been extensively investigated and substantially contributed to the derivation of the Earth gravity field (Reigber et al. 2002;Chen et al. 2005) and thermosphere total mass density (Bruinsma, Tamagnan, and Biancale 2004;Sutton, Forbes, and Nerem 2005;Bruinsma et al. 2006;Sutton, Nerem, and Forbes 2007;Doornbos et al. 2012;Jin, Calabia, and Yuan 2018). The precision of the accelerometers carried on the Gravity Recovery And Climate Experiment (GRACE) satellites can reach the level of 10 À 9 ms À 2 = ffi ffi ffi ffi ffi ffi Hz p for the less-sensitive axis and 10 À 10 ms À 2 = ffi ffi ffi ffi ffi ffi Hz p for the high-sensitive axes (Frommknecht 2007;Flury, Bettadpur, and Tapley 2008). However, due to variations of the proof mass voltage of the accelerometer as well as the variations of the on-board environments, the acceleration measurements usually require scale and bias calibration before application. For GRACE satellites, the scale errors are typically at 2-8% level while the bias errors can reach at 10 −5 m/s 2 magnitude and vary with time (Bettadpur 2003;Klinger and Mayer-Gürr 2016). These errors will contaminate the gravity and the thermosphere density signals represented by the accelerometer data. Thus, it is necessary to calibrate the onboard accelerometer measurements before applying them.
Many studies have been carried out on the onboard accelerometer calibration, among which the approaches incorporated accelerometer parameter into satellite dynamics are investigated (Bruinsma, Biancale, and Perosanz 2007;Van Helleputte and Visser 2008;Van Helleputte, Doornbos, and Visser 2009;Bezděk 2010;Visser and Vanden Ijssel 2016;Guo et al. 2018;Vielberg et al. 2018). In such calibration process, the nonconservative accelerations are represented by the accelerometer measurements and their calibration parameters. With the input of accelerometer measurements, the calibration parameters are estimated in the process of gravity field recovery or orbit determination (Reigber et al. 2002;. LEO-based accelerometer measurements have been widely investigated in gravity field recovery studies and require calibration as well. During such gravity estimation process, the LEO satellite accelerations are obtained from numerical differentiation of kinematic orbits that is fully free from any a priori force models. The accelerations are then divided into two parts, the gravitational part principally parameterized for gravity and the nongravitational part represented by the accelerometer measurements and their calibration parameters (Gerlach et al. 2003). However, gravity inversion is a very time-consuming task and brings about huge computation burden for those who only need calibrated accelerometer data for thermosphere studies.
Another typical implementation of the dynamic calibration method is to estimate calibration parameters within a reduced-dynamic precise orbit determination (POD) process (Van Helleputte and Visser 2008;Van Helleputte, Doornbos, and Visser 2009;Koch, Shabanloui, and Flury 2019). During POD, the nongravitational accelerations usually represented by the empirical models such as drag and solar radiation are replaced by the accelerometer measurements. The principle of this method is that only when the calibration parameters are accurately determined, the dynamic orbits calculated from integrations of the satellite accelerations are accurate. Due to the strong anticorrelations between scale and bias factors, Van Helleputte, Doornbos, and Visser (2009) applied very strong constraints on the bias parameters in radial and cross-track directions when calibrating the GRACE accelerometers. When accelerometer data were used and no empirical accelerations were estimated, the root mean square (RMS) value of the GRACE satellite orbit differences was just around 10 cm compared to the Jet Propulsion Laboratory (JPL) precision scientific orbit (PSO). To further improve the orbit precision, empirical acceleration models in along-track, crosstrack, and radial directions were introduced to absorb any un-/mis-modeled accelerations in their studies. By such dynamic modeling, the GRACE orbit precision could be improved to an RMS value of 3.5 cm. It is also found that the calibration factors in along-track component can be well determined due to its strong signal strength by atmospheric drag, but the radial component is barely meaningful. However, though the orbit precision reaches a higher level by estimating empirical accelerations, the estimated bias or scale parameters are no longer complete as the part of the calibration parameters can be absorbed in the empirical accelerations (personal communication, Bruinsma Sean). In addition, this method also requires the processing of raw GPS data sets, which increases complexity. Bezděk (2010) proposed another method for onboard accelerometer calibration that compares accelerometer data with the nongravitational accelerations obtained by subtracting the gravitational accelerations from total accelerations using background models. The satellite total accelerations were calculated from numerical differentiation of kinematic orbits. However, due to the errors in kinematic orbits as well as the error amplification and correlation through numerical differentiation, this method strongly depends on the following factors: (1) accurate calculation of satellite accelerations from kinematic orbits; (2) an autoregressive process to account for correlated errors; and (3) the accuracy of background gravitational models. Bezděk (2010) applied a second-derivative filter for orbit derivatives and the generalized least squares for estimation as it allows for joint estimation of calibration parameters and correlated errors, hence to produce realistic GRACE calibration parameters in agreement with Bettadpur (2003). Vielberg et al. (2018) reinvestigated this calibration method with the GRACE on-board data and refined it by using derived nongravitational accelerations as observations for direct estimation of calibration parameters as well as applying an autoregressive process to account for the correlation errors. In a similar approach, Wöske et al. (2019) calibrated the GRACE accelerometer straightforward by nongravitational force modeling using start-of-the-art models. In such procedure, nongravitational forces had to be handled very carefully, including the atmospheric drag forces and winds, as well as radiation forces due to solar radiation pressure, albedo, Earth infrared and satellite thermal radiation; and the high-precision satellitesurface model such as finite element model should be utilized for forces modeling. However, due to limited satellite drag model accuracy, this approach cannot produce precise along-track calibrations, e.g. for GRACE satellites the fitting residuals were around 10 nm/s 2 for total accelerations around ±70 nm/s 2 at low solar activity (Wöske et al. 2019); thus, it is not a suitable choice for thermosphere density retrieval that mainly depends on along-track accelerometer observations.
To date, a number of LEO missions are equipped with accelerometers on board, such as Gravity field and steady-state Ocean Circulation Explorer (GOCE), Swarm, GRACE follow-on, etc. The accuracy and efficiency of accelerometer calibration are both critical for accelerometer-based thermosphere density retrieval studies. An error of 10 nm/s 2 in the accelerometer measurements in the along-track component can lead to density retrieval errors at 2-3% level at the altitude of CHAMP (CHAllenging Minisatellite Payload) satellite (Doornbos 2012). Based on this, we explored another dynamic calibration scheme, which incorporates the accelerometer data into orbit fitting. In this method, the modeled gravitational accelerations and the accelerometer data are used in orbit integration to produce orbit that best fits the reference orbit, while the calibration parameters are estimated during the fitting process together with the satellite initial state vectors. In many orbit determination software such as PANDA (Liu and Ge 2003), orbit fitting is often used for preparing approximate initial state vectors and dynamic parameters. Thus, calibration through orbit fitting is simple to implement. In addition, the orbit-fitting method can achieve much higher efficiency than the methods proposed by Van Helleputte, Doornbos, and Visser (2009) since there is no need to process the on-board GPS data.
The GRACE L1B products from 2004 to 2007 (inclusive) are used for investigating the dynamic calibration methods. Mainly the POD and orbit-fitting calibration schemes are explored and evaluated. Data processing strategies are discussed and analyzed comprehensively, especially on the parameterization of calibration parameters. The two sets of calibration results are compared and their precision is assessed. The following sections are organized as below: the GRACE data and calibration methods are introduced in Section 2; the data processing strategy, especially the dynamic models, is discussed in Section 3; parameterization of the GRACE accelerometer measurements is analyzed in Section 4; in Section 5, the calibration results from POD and orbit fitting are presented and compared; in Section 6, the results of thermosphere density derivation using the two calibration schemes are presented and their precisions are validated; then the conclusions are given in Section 7.

Materials and methods
In this section, the data sources involved in the calibration procedure are introduced, and then the two dynamic calibration methods, orbit fitting and POD, are presented.

GRACE satellite data
The GRACE Level-1B (L1B) products from 2004 to 2007 are retrieved from the Information System and Data Center (ISDC) for this study (Case, Kruizinga, and Wu 2010). The L1B level products include measurements from different onboard sensors. In this study, the GPS data, accelerometer data, and star camera data are employed. The GRACE accelerometer data from L1B products is given in the satellite bodyfixed frame (SBF; its definition is depicted in Figure 1) (Case, Kruizinga, and Wu 2010;Dahle et al. 2013). The L1B star camera measurements, which describe the satellite attitude using quaternions, are used to convert coordinates between SBF and the celestial reference frame (CRF).
The GRACE GPS measurements are utilized during the POD calibration scheme. The GPS satellite orbit and clock products are also required in a GPS-based POD process. The 30 s GPS clock products that allow to reduce interpolation errors and the orbit products from JPL second IGS reprocessing campaign (repro2) are adopted (Desai et al. 2014). Since the JPL repro2 products are generated with the same configurations and models, they are more consistent than the legacy IGS products. However, it may increase orbit differences between the POD results and the GRACE PSO as they are generated using different GPS products.
The JPL GRACE PSO product (namely the GNV products) is obtained from typical reduced-dynamic POD procedure, which is represented by a set of background models and dynamic parameters. Fitting to the JPL PSO equals to fitting the accelerometer data backward to the dynamic models. In contrast, it is more appropriate to use kinematic orbits for orbit fitting as they are free from any a priori models. In our study, the GRACE kinematic orbits produced by Graz University of Technology are adopted for orbit fitting (Zehentner and Mayer-Gürr 2015). These orbits are calculated using raw undifferenced GPS observables with a kinematic orbit determination technique and are free from any dynamic models. They are considered to have an accuracy of 2-3 cm, which is sufficient for the calibration.
The orbit precision is regarded as a metric when evaluating the calibration performance. The GRACE orbit generated with calibrated acceleration observables is compared with the JPL PSO and validated with the satellite laser ranging (SLR) observations. The statistics of orbit differences and SLR residuals are then used for orbit precision assessment. The smaller the orbit discrepancy and SLR residuals are, the more accurate the calibration parameters are.

Accelerometer calibration by orbit fitting
The original accelerometer measurements refer to the accelerometer frame, whereas L1B accelerometer products are in SBF. For orbit calculations in the orbitfitting procedure, all satellite accelerations are represented in CRF. We convert the L1B accelerometer data from SBF to CRF using attitude transformation matrix calculated by the L1B attitude products where M is the attitude matrix. For calculating M with GRACE attitude data, please refer to Case, Kruizinga, and Wu (2010). b is the 3 × 1 bias vector consisting the accelerometer biases in SBF X, Y, Z components, while K is the scale factor matrix with a dimension of 3 × 3. a is the acceleration measurements in SBF, while e a CRF is the calibrated accelerations in CRF. Both a and e a CRF are three-dimensional vectors representing the acceleration components in X, Y, Z directions. To estimate the bias and scale factors, the partial derivatives of acceleration with respect to b and K need to be calculated when generating a priori dynamic orbit by numerical orbit integration, and they can be deduced from Equation (1) as below: In the orbit-fitting scheme, the positions from kinematic orbit r obs t at epoch t are treated as pseudoobservations for the satellite positions r t . Thus, the observation equation can be formed as In Equation (4), x 0 are the parameters needed to be estimated, including the initial satellite position and velocity at beginning epoch t 0 as well as the accelerometer calibration parameters b and K. r t j x 0 are the satellite positions taken from a priori dynamic orbit calculated with x 0 : v t is the observation error. Expanding r t in Equation (4) using the Taylor's theorem at t 0 , and then expressing them in a matrix form for subsequent adjustment, we have where Δx 0 is the correction vector for x 0 , A t is the design matrix, equaling to the product of the partial derivatives of satellite positions over satellite state at epoch t ( @r t @x t ) and the satellite state transition matrix @x t @x 0 , which can be calculated simultaneously during numerical orbit integration, and l t is the observed minus computed (OMC) residual vector. For each epoch, the observation equation can be established according to Equation (5). Normal equations can be formed by stacking all the observation equations epoch by epoch in a selected orbit arc. Thus, Δx 0 can be estimated using the least-square adjustment method by where A and L are the stacked design matrix and residual vectors, respectively. With the obtained corrections Δx 0 , the satellite initial state and accelerometer calibration parameters x 0 can then be updated. Since the observation equation is not a linear system, 3-4 times of iterations of above procedure is often required to achieve convergence. The final estimates x 0 represent the dynamic orbit that best fits the kinematic orbit.

Accelerometer calibration by POD
Accelerometer calibration in a GPS-based POD procedure differs with the orbit-fitting scheme in the respect of observation equation modeling. In the POD scheme, the GPS observations, including the pseudorange and carrier phase, are modeled with respect to the dynamic orbit. The estimated parameters include not only the initial satellite state and accelerometer calibration parameters, but also the receiver clock errors and ambiguities. In our study, the ionosphere-free pseudorange and carrier phase are utilized for observation equation modeling, which allows for elimination of the first-order ionosphere refraction delays. For a detailed description on GPSbased POD procedure, the readers are referred to Montenbruck and Gill (2000) and Jäggi et al. (2007).

Data processing strategy
The performances of orbit-fitting and POD calibration schemes both significantly depend on the satellite dynamic configurations, which means that the fitting or POD arc length and dynamic models can all affect the final calibration precision. In the GFZ GRACE Level 2 data processing standard documents, the detailed strategy for recovering Earth Gravity using GRACE data is described. In the early versions of these standards (0001/0002), the arc length was set 36 h while in more recent versions (0003/0004/0005) was 24 h (Dahle et al. 2013). In this study, the 24 h arc length is selected. The gravitational errors can be minimized by adopting background models as accurate as possible, such as utilizing the EIGEN-6C model for both static and time-varying gravity field (Foerste et al. 2011), FES2004 model for ocean tide (Lyard et al. 2006). For the GRACE satellites, the accelerometer scales and biases cannot be estimated simultaneously without applying constraints due to their strong correlation (Bruinsma, Biancale, and Perosanz 2007;Van Helleputte, Doornbos, and Visser 2009). To solve this problem, Van Helleputte and Visser (2008Visser ( , 2009 applied strong constraints on the a priori bias values in the cross-track and radial direction. Bruinsma, Biancale, and Perosanz (2007) found no significant variability of bias factors when strong constraints applied on scale factors, which supports the idea of fixing scales to constant values and estimating only biases. Such accelerometer parameterization is also consistent with the strategy proposed in Dahle et al. (2013).
The accelerometer biases are often estimated typically one set per single day (e.g. Van Helleputte and Visser 2008;Van Helleputte, Doornbos, and Visser 2009;Vielberg et al. 2018). But according to Van Helleputte andVisser (2008, 2009), the orbit precision estimated in such strategy can only reach at 10 cm level when no empirical accelerations are applied, which is significantly degraded compared to the GRACE POD results reported in other studies (Kang et al. 2006;Jäggi et al. 2007). Even though the precision is much improved after introducing the empirical accelerations, the bias estimations are subject to nontrivial deviation as the empirical accelerations and the biases are correlated. In the GFZ processing standard release 0005, it is suggested the scale factors are kept as constants while the biases are estimated every 60 min (Dahle et al. 2013). In addition, there could also be other couplings in the biases, e.g. a dependence on temperature. In our calibration strategy, the temperature-induced biases variations can be absorbed by the piecewise-constant biases and then represented in long-term trend. It should also be noted that the GRACE satellite thermal control during 2004-2007 was in a good condition and thus the thermal coupling effects were small. Thus, these technical couplings can be ignored at this point.
For the orbit-fitting scheme, the kinematic positions are treated as pseudo-observations. However, when performing POD calibration scheme, certain errors in GPS observations should be taken into consideration. Hence, the GPS observation model configurations can impact the POD calibration precision. To minimize the orbit discrepancies with respect to the JPL PSO, similar models are adopted to form the onboard GPS observation equation. Table 1 lists the detailed dynamic and observation models for both orbit-fitting and POD calibration schemes. The EIGEN-6C model is adopted to calculate the Earth's gravity, including both static and periodic terms. The solid Earth tide, pole tide are calculated using models recommended by IERS 2003 (McCarthy and Petit 2004). Ocean tide is computed using the FES 2004 model, while the atmosphere tide using GRACE AOD products with degree and order of 40. Particularly, the accelerometer data is first calibrated using the approximations provided in Bettadpur 2003, and then used in orbit integration. Different bias intervals for accelerometer are explored, and the results are given in Section 4. We use the ionospheric-free combination (PC: code combination, LC: phase combination) to eliminate the first-order ionosphere delay while the second and higher orders are neglected. The JPL repro2 precision GPS ephemeris and 30 s clock products are adopted for GPS satellite positions and clock errors. GPS satellite antenna phase center offset (PCO) and phase center variations (PCVs) are corrected using IGS antenna models (Schmid et al. 2005), while the receiver antenna PCO is corrected using the values provided by ISDC VOGPSN products (Case, Kruizinga, and Wu 2010).

Parameterization for accelerometer measurement
To evaluate the parameterization for accelerometer measurements, we processed the GRACE-A satellite data from day of year (DOY) 100 to 150 in 2006 (the time system is GPS Time) using the orbit-fitting calibration scheme with different estimation interval for bias factors. Dynamic models are adopted as listed in Table 1 with the estimation interval for the bias factor set as 45, 60 , 90, 180, and 360 min (respectively indicated as acc_45, acc_60, acc_90, acc_180, and acc_360 hereafter), in which 90 min is just about one orbit period of the GRACE satellites, 45, 180, and 360 min are the half, double, and quadruple of GRACE orbit period, respectively, while the 60 min interval is tested here as it is recommended in Dahle et al. 2013. After orbit fitting, the reduced-dynamic orbit generated from the fitted state vectors and calibration parameters is compared to the JPL PSO and validated by the SLR technique. The parameterization that gives highest precision statistics results will be considered as the best scheme against others. The 24 h estimation interval is not tried here as Van Helleputte and Visser (2008Visser ( , 2009) have already shown that the orbit precision can only reach 8.8 cm with respect to the JPL PSO. Figure 2 illustrates the RMS values of the orbit difference comparing to JPL PSO with different estimation intervals for each arc. The precision in the radial and along-track direction of all the five cases is very similar and their RMS series fluctuate steadily and smoothly except the results on DOY 112/2006 for cases acc_45 and acc_90. Generally, the radial and along-track precision is around 1.0 and 2.0 cm, respectively. However, the precision of cross-track component reaches around 1.5 cm for 45 min and 60 min interval, but degrades significantly to about 2.0-5.0 cm when the estimation interval is increased to 90,, 180, and 360 min. A very large orbit difference RMS for the arc in DOY 113 is noticed in all three components for case acc_45 and in cross-track component for case acc_90. This is possibly related to that there are no available kinematic orbit positions during 07:52:50-09:00:00 on 23 April 2006 (DOY 113/2006). In a kinematic POD process, orbit gaps occur quite frequently that is largely due to poor observation quality or data outage, and a variety of factors could induce such results such as strong ionosphere activity, signal loss, receiver reboot, etc. (Montenbruck et al. 2018). Table 2 shows the mean RMS value of the daily orbit differences for each component. Orbit precision with respect to JPL PSO is roughly all better than 5.0 cm (3D RMS) in all cases if the results on arc DOY 113 are excluded, while using the case of estimation interval 60 min achieves the best orbit precision against the others, which can reach 1.8, 1.4, and 1.0 cm in the along-track, cross-track, and radial directions, respectively. The precision in cross-track direction varies to somewhat larger extent among different cases compared to that in the along-track and radial directions, which might suggest it is more sensitive to the bias parameterization.

Orbit precision comparison for different bias intervals
The SLR fit is also used here for orbit precision validation. Table 2 shows the mean and STD of SLR residuals for the five cases. There is a bias about 0.2-0.5 cm in all the orbit solutions, and the STDs are all below 2.5 cm. The best precision (1.78 cm) is achieved by acc_60, which is consistent with the orbit comparison results.  The orbit precision in these five cases is very comparable to the result by Van Helleputte, Doornbos, and Visser (2009) that the empirical accelerations are estimated every 600 s in all three directions, but much better than that without estimating empirical accelerations. The 60 min interval gives best performance, which is exactly in agreement with the GFZ processing strategy. Figure 3 depicts the estimated bias series in X, Y, and Z directions (in SBF, denoted as Bx, By, and Bz hereafter) of the five cases. Since the actual attitudes of GRACE satellites are very close to their nominal ones, Bx, By, and Bz generally indicate biases in the along-track, cross-track and radial components, respectively. Generally, Bx, By, and Bz series are in the same magnitude among different cases. An outlier estimate in cross-track direction for cases 45 and 90 min is observed, which should be attributed to the missing positions from the kinematic ephemeris and is consistent with the orbit RMS statistics in Figure 1. During this period, the Bx, By, and Bz biases from different cases all fluctuate around −1.170 × 10 −6 , 2.886 × 10 −5 , −5.440 × 10 −7 m/s 2 , respectively, and show very consistent variation trends.

Comparison of estimated bias parameters
The mean and STD values of the estimated bias parameters are listed in Table 3. The differences of their bias mean values among different cases are quite small, mainly at the magnitude of 10 −9 m/s 2 . The Bx series generally exhibits the smallest amplitude  variation among the three components. In addition, its variations reduce significantly when the estimation interval is extended, which is revealed by the STD statistics as well. In the cases of acc_90, acc_180, and acc_360, the Bx estimates are almost constant during the testing period, and their STDs are 1.022 × 10 −9 , 6.510 × 10 −10 , and 4.910 × 10 −10 m/s 2 , respectively, while in cases acc_45 and acc_60 the Bx STDs are obviously larger, reaching 4.520 × 10 −9 and 2.078 × 10 −9 m/s 2 , respectively. On the contrary, the By and Bz series from different cases show obviously larger variations. Their STDs are generally at the level of (1-2) × 10 −8 m/s 2 and grow larger when estimation interval exceeds 60 min. This might be related to several reasons. On one hand, sensitivity of the accelerometer is weakest in the SBF-Y direction; on the other hand, dynamic modeling of GRACE-A satellite is incomplete in our study. The un-/mis-modeled gravitational forces including ignored high-order Earth gravity field (over order and degree of 130) and ocean tidal perturbations (over order and degree of 40) mainly act on the SBF-Z direction. This would result in more variable Bz estimates. According to Figure 3, in the cases of acc_45 and acc_60, the estimated Bx, By and Bz series mostly give signatures of high-frequency noise pattern. However, they tend to reveal periodical variations in very similar frequencies when estimating in larger intervals, as seen in the cases of acc_90, acc_180, and acc_360. The reason for these periodical signals is still under investigation.
When using accelerometer data for thermosphere derivation, the precision of calibration parameters in the along-track direction is much more important than that in other directions as the satellite drag acceleration mainly acts in this direction. The estimates for Bx from the above test cases all give very similar results except the case acc_45, which tends to overestimate the Bx parameters, likely due to over-parameterization.

Comparison of different dynamic calibration schemes
In this section, the GRACE data from 2004 and 2007 were processed and calibrated using the POD and orbit-fitting schemes, respectively. Their calibration results were then compared. According to the test results in the above section, the acc_60 case was selected as dynamic configuration for the GRACE accelerometer parameterization.

Orbit precision comparison
Similar to the above section, the precision of dynamic calibration schemes was firstly evaluated by the orbit precisions. For illustration, the GRACE-A satellite orbits from POD and orbit-fitting schemes were compared to the JPL PSO for each arc, and the RMS values of their orbit differences in along-track, cross-track, and radial components were calculated and are displayed in Figure 4. Generally, the orbit precisions of the two schemes can both reach cm level. Compared to JPL PSOs, the POD results give better agreement in the along-track component than the orbit-fitting method, while they agree in the cross-track direction. The radial component precision of orbit-fitting scheme demonstrates an obvious advantage over the POD scheme. The overall mean RMS values of orbit precisions with respect to JPL PSOs are listed in Table 4. The orbit precision from orbit fitting is 2.21, 1.55, and 1.27 cm in the along-track, cross-track, and radial component, respectively, while that from POD is 1.36, 1.49, and 1.65 cm in the three components. Their 3D position precisions are 2.98 and 2.61 cm, respectively. In general, the POD precision is slightly better than orbit fitting. This should be due to that the state transition constraints during a kinematic orbit determination process are much looser than those during a dynamic orbit determination since the satellite positions must be estimated for every epoch, which  can then result in worse orbit precision. As confirmation, we compared the Graz kinematic orbits and the JPL PSOs, which indicates that their differences are 2.31, 2.31, and 2.59 cm in the along-track, cross-track, radial components. This is even worse than the orbit precision of the orbit-fitting scheme. The radial precision from orbit-fitting scheme is better than POD, which might be attributed to that the kinematic positions involved in orbit fitting are vectors along the radial direction, thus making the constraints on this component stronger than others. In addition, the orbit differences show larger discrepancy during the first half year in 2004, which is in agreement with the results given in Zehentner and Mayer-Gürr (2015). The SLR observations are also taken in for orbit precision evaluation of the two schemes. The SLR OMC mean and STD of GRACE-A satellite for each arc are shown in Figure 5, while their mean values are listed in Table 4. For both schemes, the SLR OMC mean values fluctuate mostly in the ranges from −5 to 5 cm and their STDs are mainly below 5 cm. The SLR OMC mean and STD values of orbit-fitting scheme are −0.53 cm and 2.52 cm, which is larger than those of the POD scheme. This is consistent with the orbit comparison results.

Calibration parameters comparison
The bias parameters of GRACE-A satellite obtained using the orbit-fitting and POD methods are illustrated in Figure 6. Both schemes reveal noticeable gradual trending changes for the bias parameters especially in the SBF-X and SBF-Y directions, which fluctuate from −1.12 × 10 −6 m/s 2 to −1.20 × 10 −6 m/s 2 and from 2.8 × 10 −5 m/s 2 to 2.9 × 10 −5 m/s 2 , respectively, while the bias parameters in SBF-Z direction vary around −5 × 10 −7 nm/s 2 . The Bx and By parameter time series gives significant signatures of bias drifts, which is consistent with the previous studies (Bruinsma, Biancale, and Perosanz 2007;Van Helleputte, Doornbos, and Visser 2009;Klinger and Mayer-Gürr 2016). For both schemes, the bias drift in the SBF-X component is about 14 nm/s 2 per year while in the SBF-Y component 200 nm/s 2 per year. Such drifts could have resulted from mixed reasons. Parts of the drifts can be results of real satellite hardware changes, such as the spacecraft thermal variations in the long-term scale, which has been already reported (Doornbos 2012;Klinger and Mayer-Gürr 2016). Jumps of the bias parameters can be observed for both schemes in some arcs, which should be due to the effects of specific events, such as maneuvers, instrument reboots or switching from main to redundant board and vice versa. Comparing the calibration biases from the two schemes, excellent agreement can be observed in the SBF-X and SBF-Y components. However, there exists a systematic difference in the SBF-Z component between the two schemes. Since we applied the same dynamic models in the orbit-fitting and POD schemes, such systematic Bz differences might be due to the inconsistent handling of the GRACE-A satellite GPS data between our study and the Graz kinematic orbits.
The differences of the bias parameters between the orbit-fitting and POD schemes are illustrated in Figure 7. In comparison, the differences of the Bx parameters are the smallest, which mainly fluctuate within several nm/s 2 . However, the differences in the SBF-Y and SBF-Z components are at much larger magnitude about several tens of nm/s 2 while an additional offset can be found in the Bz differences. The mean and STD values of the differences are presented in Table 5. For Bx, the mean and STD values of the differences are only 0.52 ± 2.97 nm/s 2 , while they are −1.48 ± 31.82 nm/s 2 and 42.71 ± 22.63 nm/s 2 for By and Bz differences, respectively. The small differences  in Bx estimates indicate that both schemes can derive thermosphere densities at very similar precision. Treating the bias parameters from orbit-fitting scheme as references, the relative differences between the two schemes are also calculated. For Bx, the relative differences are very small, at −0.04% ± 0.26%, while the Bz parameter differences are as large as −7.25% ± 3.85%.
In many previous studies, the GRACE satellite accelerometer calibration parameters are usually estimated on a daily basis (e.g. Van Helleputte and Visser 2008;Van Helleputte, Doornbos, and Visser 2009;Vielberg et al. 2018). To facilitate comparison, we averaged the 60 min Bx, By, and Bz parameters from the orbit-fitting scheme into daily values. The results also demonstrate that for Bx the daily STD is under 4 nm/s 2 and is mainly less than 0.3% over its daily mean value; for By and Bz, their STDs both vary at a higher level and are typically under 20 nm/s 2 . These daily STD values can represent the bias estimation precision and are also consistent with the results in Vielberg et al. (2018).

Comparison of the calibration schemes on thermosphere density derivation
Thermosphere densities are important parameters in modeling the space environment. Deriving thermosphere densities through LEO-on-board accelerometers are attracting greatly increased attention in recent decades due to its good accuracy as well as huge data number. In this section, we investigated thermosphere density derivation using the GRACE accelerometer data by applying the obtained Bx parameters from the two calibration schemes described above. The derived two sets of densities were firstly compared to each other and then compared to the estimates by different approaches, including empirical models as well as densities calculated by other researchers.
To retrieve densities from accelerometer data, mainly the accelerometer measurements of the SBF-X components are explored. Detailed procedures on the density derivation can be found in Bruinsma, Tamagnan, and Biancale 2004;Sutton, Nerem, and Forbes 2007;Doornbos 2012;Vielberg et al. 2018. Firstly, the satellite aerodynamic accelerations are isolated by subtracting the solar radiation pressure (SRP) and Earth radiation pressure from the calibrated accelerometer data. Thus, we have a X;obs ¼ a X;cal À a srp À a alb À a IR (7) In Equation (7), a X;cal and a X;obs are the calibration accelerometer data in SBF-X and observed drag acceleration, respectively; the SRP accelerations a srp is calculated using the GRACE macro model (Bettadpur 2012); the Earth radiation accelerations including albedo a alb and infrared radiation a IR are calculated using grid model (Bhanderi and Bak 2005). In practice, a alb and a IR can be ignored as they are often 2-3 orders of magnitude lower than the aerodynamic accelerations and act mainly in the radial direction.
With the observed drag accelerations, atmosphere density can be calculated by where m is the satellite mass, v r is the relative velocity between satellite and atmosphere, and C A,X is the SBF-X component of satellite aerodynamic coefficient. As seen, the key factor in deriving the density from Equation (8) is the C A,X parameter. We selected the Sentman model for GRACE satellite aerodynamic coefficient calculation, which has been investigated in previous studies (Sutton, Nerem, and Forbes 2007;Doornbos 2012). In Sentman aerodynamic coefficient calculation, the atmosphere wind velocity is calculated using the HWM07 model (Drob et al. 2008), the atmosphere temperature is calculated using the NRLMSISE00 model (Picone et al. 2002), the satellite wall temperature is set as 273 K as recommended in Doornbos (2012), and the accommodation coefficient is set as 0.93 as recommended in Sutton, Nerem, and Forbes (2007).
To evaluate our density estimates, empirical models such as NRLMSISE00 and High Accuracy Satellite Drag Model (HASDM) (Storz et al. 2005) are used. The NRLMSISE00 model is developed based on the previous MSISE models and is widely adopted in atmosphere studies as well as LEO drag analysis. The typical errors of NRLMSISE00 are around 10-30% and  grow with increasing altitudes (Marcos, Bowman, and Sheehan 2006). At the altitude of GRACE satellite, its errors are mainly about 25%. The HASDM model is an empirical model based on the Jacchia-71 model. It employs drag data from 75-80 LEOs to estimate the inflection point temperature and the exosphere temperature every 3 h (Bowman and Storz 2003). The errors of HASDM are confirmed at 6-8% at the altitudes between 200 and 800 km (Storz et al. 2005). The GRACE accelerometer-based densities published by Doornbos (2012) are also taken for comparison purpose (referred to as Density_Direct hereafter).
In Figure 8, thermosphere densities from multisources are illustrated from 00:00 to 04:00 on DOY 002/2004, where specifically densities labeled as POD and orbit fitting are retrieved using the calibration parameters from the POD and orbit-fitting schemes in this study, respectively. The density variations from different sources coincide with each other in general, and their correlation coefficients between one and another all exceed 0.90. However, compared to the accelerometer-derived densities, the two empirical models both overestimate the densities to some extent. Specifically, the NRLMSISE00 model overestimates the densities during midnight while HASDM shows better agreement with POD, orbit fitting, and Density_Direct. The three accelerometer-based density data sets show excellent agreement with correlation coefficients above 0.99. The mean and STD values of the density ratio of NRLMSISE00, HASDM, and Density_Direct over POD and orbit fitting densities are listed in Table 6. Large bias of the NRLMSISE00 model is observed around 50% with respect to densities from both POD and orbit-fitting schemes and its STD is around 25%, which is in accord with its typical precision at this altitude. Compared to NRLMSISE00, the HASDM exhibits a much smaller bias at the order of 10% and a reduced STD value around 20%. Both the POD and orbit-fitting-derived densities demonstrate good agreement with the Density_Direct data, with bias below 5% and STD around 5%.
We compared all the POD and orbit-fitting-derived densities from 2004 to 2007 and obtained over 800 million density ratios. A log-normal distribution is applied to fit the density ratios to analyze the differences between these two density data sets. The frequency counts of the density ratios as well as the results of log-normal distribution fitting are illustrated in Figure 9. The two density data sets agree with each other with a very high correlation coefficient of 0.9976. A small bias of 2.20% is found between the orbitfitting and POD densities and their log-normal STD is 5.75%. These figures are obviously larger than the Bx differences given in Table 5, which should be attributed to the fact that the accelerometer-derived densities are not linearly related to the calibration parameters.
All the densities during 2004 and 2007 from POD and orbit-fitting schemes are compared to NRLMSISE00, HASDM, and Density_Direct for the precision evaluation. The log-normal distribution fit is applied to their density ratios; the fitting mean and   Figure 9. Frequency count of density ratios of orbit fitting over POD densities and its log-normal distribution fit. The bin interval for the frequency count is 0.01. Table 7. In general, the thermosphere densities generated by orbit fitting and POD are at very similar precision level, and they indicate very similar statistics when compared to the empirical models as well as the accelerometer-based densities. The NRLMSISE00 model indicates a bias of 33.7% and an STD of 0.262 with respect to the POD densities, as well as a bias of 29.9% and an STD of 27.2% with respect to the orbit-fitting densities. For HASDM, the bias is improved to 8.8% and 5.9% while the STD is reduced to 20.5% and 21.0% compared to the POD and orbit-fitting results, respectively. The density differences between POD, orbit fitting, and Density_Direct are the best among others. The bias and STD values between POD and Density_Direct are 1.6% and 5.7%, respectively, while those between orbit fitting and Density_Direct are −0.1% and 5.6%.

Conclusion
In this study, the GRACE on-board accelerometer calibration is investigated using the dynamic calibration method. Two different methods are implemented and compared. The first method is to estimate the calibration parameters during a POD procedure that combines the on-board GPS data, while the second during an orbit-fitting process that requires kinematic orbits. To confirm and validate the proposed calibration schemes, the GRACE data during 2004 and 2007 is analyzed.
The GRACE accelerometer calibration is investigated focusing on its parametrization. The scale coefficients of GRACE accelerometer are fixed to their a priori values while the bias parameters are estimated using piecewise-constant model. Different bias estimation intervals are tested in this study to optimize the orbit fit and orbit determination precision. The results indicate that estimating the bias parameters every 60 min can achieve orbit precision at 2 cm level, as confirmed by orbit comparison to the JPL PSOs as well as the SLR fitting residuals. With this processing strategy, the integrity of the accelerometer biases is preserved since no extra empirical dynamical compensations are estimated during the calibration. The estimated biases from the POD and orbit-fitting schemes show good agreement in the SBF-X and SBF-Y components with differences of −0.04% ± 0.26% and 0.00% ± 0.11%, respectively, while there is a systematic offset in the SBF-Z component at about −7.25% ± 3.85%.
We derived the densities from GRACE-A satellite accelerometer data using the Bx parameters from orbit-fitting and POD schemes, respectively. The density differences between these two schemes are very small and are only at 2.2%±5.75% level. The derived densities are also confirmed by comparison to empirical models such as NRLMSISE00 and HASDM as well as the accelerometer-based density products published by Doornbos (2012), finding that the differences are 1.6%±5.7% and −0.1%±5.6% for the POD and orbitfitting schemes, respectively. The density comparison of POD and orbit fitting over the NRLMSISE00 model reveals that the NRLMSISE00 exhibits a bias about 30% and a STD about 27% during 2004 and 2007 at the altitude of GRACE satellite, while that over the HASDM gives an improved bias at 7% level and a reduced STD at 21% level.
Considering the increasing LEO-based accelerometer measurements from the ongoing SWARM and GRACE follow-on missions, the efficiency in accelerometer calibration can be very critical in deriving the thermosphere density. The orbit-fitting scheme proposed in our study can avoid the complexity in processing the onboard GPS data and more importantly it can deliver calibration precision at similar level with respect to the POD scheme. Thus, it is very suitable in thermosphere density derivation from on-board accelerometers.