Comparison of accelerometer data calibration methods used in thermospheric neutral density estimation

. Ultra-sensitive space-borne accelerometers on-board of low-Earth orbit (LEO) satellites are used to measure non-gravitational forces acting on the surface of these satellites. These forces consist of the Earth radiation pressure, the solar radiation pressure and the atmospheric drag, where the ﬁrst two are caused by the radiation emitted from the Earth and the Sun, respectively, and the latter is related to the thermospheric density. On-board accelerometer measurements contain systematic errors, which need to be mitigated by applying a calibration before their use in gravity recovery or thermospheric neutral 5 density estimations. Therefore, we improve, apply, and compare three calibration procedures: (1) a multi-step numerical estimation approach, which is based on the numerical differentiation of the kinematic orbits of LEO satellites, (2) a calibration of accelerometer observations within the dynamic precise orbit determination procedure, and (3) a comparison of observed to modeled forces acting on the surface of LEO satellites. Here, accelerometer measurements obtained by the Gravity Recovery And Climate Experiment (GRACE) are used. Time series of bias and scale factor derived from the three calibration proce-10 dures are found to be different in time-scales of a few days to months. Results are more similar (statistically signiﬁcant) when considering longer time-scales, from which the results of approach (1) and (2) show better agreement to those of approach (3) during medium and high solar activity. Calibrated accelerometer observations are then applied to estimate thermospheric neutral densities. Differences between accelerometer-based density estimations and those from empirical neutral density models, e.g., NRLMSISE-00, are observed to be signiﬁcant during quiet periods, on average 22 % of the simulated densities (during 15 low solar activity), and up to 28 % during high solar activity. Therefore, daily corrections are estimated for neutral densities derived from NRLMSISE-00. Our results indicate that these corrections improve model-based density simulations in order to provide density estimates at locations outside the vicinity of the GRACE satellites, in particular during the period of high solar/magnetic activity, e.g., during St. Patrick’s Day storm on March 17 th , 2015.

1 Introduction Recent gravimetric satellites, for example, the satellite missions CHAllenging Minisatellite Payload (CHAMP) (Reigber et al., 2002) or Gravity Recovery And Climate Experiment (GRACE) (Tapley et al., 2004) are equipped with ultra-sensitive spaceborne accelerometers that allow for the measurement of non-gravitational forces acting on the surface of these satellites.These measurements reflect accelerations due to the atmospheric drag (the dominant component of the acceleration vector at the orbital altitude of these satellites), and thus enable studies of thermospheric neutral density and winds (e.g., Bruinsma et al., 2004;Liu et al., 2005;Sutton et al., 2007;Doornbos, 2012;Lei et al., 2012;Mehta et al., 2017).Non-gravitational forces also contain the effect of the Solar and Earth radiation.
Nevertheless, satellite accelerometer measurements need to be calibrated before being used in any applications such as solar terrestrial studies (e.g., Doornbos, 2012), gravity field recovery (e.g., Reigber et al., 2003) or precise orbit determination (e.g., Van Helleputte et al., 2009).This is mainly because of systematic errors that contaminate the sensor data.Calibrating the ultra-sensitive space-borne accelerometers before the satellite's launch is not possible, since gravity on the Earth's surface is too large and simulating the space environment is extremely difficult.Therefore, several studies have been developed during the last decade to ensure in-orbit calibration of accelerometer measurements.For example, Kim (2000); Tapley et al. (2004), and Klinger and Mayer-Gürr (2016) calibrate GRACE accelerometer observations within a gravity field recovery procedure.Bezděk (2010) and Calabia et al. (2015) apply numerical differentiation techniques, e.g., developed by Reubelt et al. (2003), to compute accelerations from precise kinematic orbits, which are then used to estimate calibration parameters and their uncertainties.Alternatively, calibration parameters can be estimated within the precise orbit determination procedure (Bettadpur, 2009;Van Helleputte et al., 2009;Visser and Van den IJssel, 2016).Each method mentioned above yields different calibration parameters and their uncertainty, as well as their influence on the final products such as thermospheric neutral density estimation has not yet been systematically investigated.
In order to better understand and reconcile differing results in the literature, in this study, three calibration procedures are applied to GRACE accelerometer measurements.The aim is to assess the impact of a specific calibration method on the estimation of global thermospheric neutral densities as will be discussed in what follows.
(1) The first approach is called here Multi-Step Numerical Estimation (MNE), which is based on the numerical differentiation of kinematic positions.The application of this method is similar to that of Bezděk (2010) with few differences concerning the orbit data and the stage in which calibration parameters are estimated.
(2) The second approach calibrates GRACE accelerometer measurements within the dynamic precise orbit determination procedure (Löcher, 2011).(3) Finally, calibration parameters are obtained by comparing the accelerometer measurements to modeled non-gravitational forces acting on the satellite.This procedure is commonly used to find initial calibration parameters in gravity recovery experiments (e.g., Kim, 2000;Van Helleputte et al., 2009;Klinger and Mayer-Gürr, 2016).
In recent decades, empirical and physical models of the atmosphere have gone through considerable development, while reflecting the range of density variability in response to solar and geomagnetic forcing.The Mass Spectrometer and Incoherent Scatter (MSIS) empirical models of the neutral atmosphere (Picone et al., 2002) have been developed since 1977.Other models such as Jacchia-Bowman (e.g., Bowman et al., 2008) have also been used in various satellite applications.The current version of the MSIS model, NRLMSISE-00, and Jacchia-Bowman 2008 are built from an extensive drag data set and they are parameterized in terms of solar and magnetic indices at daily and 3-hourly resolution, respectively.In this study, we show to what extent GRACE-derived calibrated accelerometer data affect the final estimation of atmospheric neutral density.Furthermore, as empirical thermospheric models fall short in simulating thermospheric neutral density (mostly at short-time scales and during high solar/magnetic activity, Bruinsma et al., 2004;Guo et al., 2007), daily empirical corrections are estimated, which can be applied to scale the outputs of these models, and therefore, improve their global performance particularly during high geomagnetic activity (see Sect. 4.3).
In the following, data sets and models are introduced in Sect.2, and the methodologies of calibration are discussed in Sect.
3. The results are presented in Sect. 4 and the study is concluded in Sect. 5.

GRACE data
We use GRACE Level-1B data (Case et al., 2010) 1 provided in the science reference frame (SRF) located at the centre of mass of each satellite.The axes of the SRF are parallel to the accelerometer frame (AF) and the satellite frame (SF, see Fig. 1).
Following this, the x-axis points to the phase centre of the K-Band instrument (along-track or anti-along-track depending on leading and trailing satellite), the z-axis is directed to the normal of the x-axis and the main equipment platform plane.The y-axis completes the right-handed triad.

Accelerometer
Space-borne capacitive accelerometers such as the SuperSTAR accelerometer on-board each GRACE twin satellites contain a proof mass, which is kept at the centre of mass of a satellite by compensating the non-gravitational forces with induced electrostatic forces.The measured accelerations of the proof mass of the SuperSTAR accelerometer, which are proportional to the voltage needed to generate the compensating electrostatic forces, are labeled as ACC1B within the GRACE Level-1B data.
This accelerometer has a resolution of 10 −10 ms −2 in the x-and z-directions, whereas the resolution of the y-axis is found to be one magnitude lower (Flury et al., 2008).The temporal resolution of the ACC1B data is 1 sec.Frommknecht (2007) showed that the noise level of the along-track and radial components of the accelerometer is 2-3 times higher than that specified in the handbook.However, the noise level is reported to be similar for both satellites.Possible reasons for systematic effects are axis non-orthogonality, displacement of the test mass with respect to the satellite's center of mass, and thermal effects (Frommknecht, 2007).

Star Camera
Two star cameras on-board each GRACE satellite provide its inertial orientation in terms of quaternions, labeled as SCA1B.
These data are used (in Sect.3.1) to transform measurements given in the SRF to the celestial reference frame (CRF).

Macro-model
In this study, a macro-model is required to model the non-gravitational accelerations acting on the surface of the satellite (Sect. 3.3).The geometry of the two identical GRACE satellites is represented in a macro-model (Bettadpur, 2012) including mass, surface area and material of each plane in terms of visible and infrared reflectivity coefficients for specular and diffuse reflection (see Tab. 1).The characteristics of the macro-model have been determined under laboratory conditions; their values may be different under space conditions.These values also change during the mission's life-time due to aging of surface coating under UV radiation (Vallado and Finkleman, 2014).Following this, we make use of the macro-model at hand, and an estimation of their uncertainty, as well as their impacts on the final thermospheric density estimations will be addressed in another study.

Precise orbits
Each GRACE satellite is equipped with three GPS receivers, whose data is used for precise orbit determination (POD) and which ensure the precise time tagging of other on-board measurements.Level-1B GRACE satellite orbits (GNV1B) are obtained from a dynamic POD procedure (Case et al., 2010).However, the choice of orbit data depends on their application.Since here we aim at calibrating accelerometer measurements for thermospheric density estimations, the kinematic precise orbits, which are free from accelerometer measurements are applied.For this, kinematic orbits are processed by Graz University of Technology (Zehentner and Mayer-Gürr, 2016) 2 .

Models
In this study, neutral thermospheric densities derived from two empirical density models JB2008 (Bowman et al., 2008) and NRLMSISE-00 (Picone et al., 2002) are compared to accelerometer derived densities from GRACE.The JB2008 model uses being a proxy for the solar electromagnetic radiation at a wavelength of 10.7 cm = 2.8 GHz.NRLMSISE-00 can be used for the entire space age (applied magnetic and solar indices start before 1957), and it predicts atmospheric temperature, density, and composition.JB2008 is constructed using a combination of solar-and geomagnetic proxies and indices that have been available since 1998, thus the model cannot be used before that year (Bruinsma et al., 2017).Both models account for spatial and temporal variations in the solar activity.Assessments of JB2008 and NRLMSISE-00 simulations indicate that JB2008 is closer to independent observations during average solar activity (Liu et al., 2017).The models show limited skill during the high solar/magnetic activity as model equations do not perfectly reflect timing of the heating transfer as demonstrated by Weimer et al. (2011).
As we will show in Sects.3.1 and 3.2, gravitational forces must be known while performing the two calibration procedures of the Multi-Step Numerical Estimation (Sect.3.1) and the Dynamic Estimation (Sect.3.2).Here, these forces are accounted for using background models as listed in Tab. 2.   (Petit and Luzum, 2010) ocean tides FES 2004 (Lyard et al., 2006) pole tides IERS Conventions 2010 (Petit and Luzum, 2010) pole ocean tides Desai 2004 (Petit and Luzum, 2010) 2009; Calabia et al., 2015).A parameterization using a fully populated scale factor matrix has been discussed by Klinger and Mayer-Gürr (2016), which is neglected here since its impact on final thermospheric density estimations is negligible.
In this study, the calibration equation is written as where a ng contains modeled non-gravitational accelerations, a raw the measured ones, and finally v represents errors.Since this parameterization yields highly (anti-)correlated calibration parameters, which cannot be physically interpreted, we follow an iterative estimation of the calibration parameters as recommended by Van Helleputte et al. (2009).During this iterative procedure, (1) daily calibration parameters are estimated following Eq.(1), then (2) the scales of step (1) are temporally averaged for each direction for the whole period of available data, and finally (3) daily biases are re-estimated with the constant scales computed in step (2).Then, the estimated scales are used to estimate daily biases for the three different approaches described in the following sections.

Multi-step numerical estimation (MNE)
This calibration method makes use of a two-fold numerical differentiation of kinematic orbit positions, a procedure that has been often applied and tested in gravity retrieval studies (e.g.Reubelt et al., 2003).The idea is that by applying a second numerical differentiation on kinematic orbit positions, one is able to obtain an estimate for the satellite's total acceleration a total .To obtain such derivatives, one option is to apply central differentiation operators, for example using seven-(or more) points for calculation (e.g., Reubelt et al., 2003).This implementation, however, introduces three main problems: (1) unwanted phase differences because of the averaging kernel in the differentiation operator and assuming linearity of the path, (2) ampli-fying the noise because of random errors in the original data points, and (3) amplifying temporally correlated noise because the differentiation operator convolves successive orbital positions within the filtering window.To mitigate these problems, the Savitzky-Golay filter (Savitzky and Golay, 1964) is recommended in Bezděk (2010), which combines smoothing and differentiation operators.Here, we perform a synthetic experiment, where different numerical derivatives are applied to orbit positions sampled from a true orbit defined through an analytical representation (see Appendix A, Eq. ( A1)).By definition, the main frequencies and amplitudes of the true orbit are known, and subsequently various numerical derivatives can be compared with the results of the analytical derivative.Our results (see Appendix A) confirm that the Savitzky-Golay filter is a suitable approach to estimate derivatives from kinematic orbits because the amplification of noise is limited and phase shifts are prevented.
The problem of designing the Savitzky-Golay filter is equivalent to finding coefficients c n,m , obtained from a least squares adjustment, which can be convolved with the time series of the kinematic positions to compute derivatives (Bezděk, 2010).
Thus, the derivative of the orbit r, shown by r, can be obtained as Here, r contains the total acceleration (a sum of both gravitation and non-gravitational forces).In Eq. ( 2), i is the epoch index to which the filter is applied; the width of the window is denoted by n and m expresses the order of the fitted polynomial.
From our numerical experiments (see Appendix A), combining a window length of n = 11 data points with a polynomial of degree m = 5 are found to be the best second derivative filter settings for the kinematic orbits used in this study, i.e. this filter adds less noise to the true derivative than others.The filtering requires equally spaced data, thus in this study the kinematic positions are interpolated to an interval of exactly 10 seconds using a cubic polynomial interpolation.We found that changing the interpolation technique to polynomial or a harmonic interpolation does not considerably change the final results.Data gaps are bridged by fitting a polynomial of degree 9 to one orbital revolution, where the chosen degree minimizes the difference between the true orbit and an interpolated orbit with simulated gaps.
The total acceleration r at a specific time t is related to the force f = ma acting on the satellite via the equation of motion The acting acceleration a depends on the time t, the satellite's position r, velocity ṙ, and force model parameters denoted by p consisting of a gravitational and a non-gravitational part.To remove gravitational accelerations a grav from Eq. ( 3), models of Tab. 2 are used and the desired non-gravitational accelerations acting on the spacecraft are estimated as In order to allow a physical interpretation of the non-gravitational accelerations, these are transformed from the celestial to the science reference frame by applying a rotation matrix that is derived from the entries of the star camera quaternions.
The resulting non-gravitational accelerations a ng can then be used to calibrate the on-board non-gravitational accelerometer measurements a raw .
Outliers in each direction of the non-gravitational accelerations a ng , which are caused by the constant time step interpolation of kinematic orbits, are detected and removed using the one-dimensional statistical test that considers the modified standard deviation σ a = ( i (a nga i − ãnga ) 2 )/(r − 1) for each axis (a = 1, 2, 3) with the number of data r based on the median ãnga instead of the mean ānga .Accelerometer observations a raw for the same epochs are discarded as well to keep the time domain of the data sets in agreement.Based on the data snooping procedure (Baarda , 1968), the outlier detection is iteratively repeated until the standard deviation is found to be smaller than the global standard deviation computed on the basis of days when no gap-filling interpolation was needed.The threshold in the along-track and cross-track directions is found to be 1 • 10 −5 ms −2 , whereas the threshold in the radial direction is equal to 2 • 10 −5 ms −2 .
Since errors of a ng are not Gaussian distributed (as a result of the numerical differentiation), an ordinary least squares estimation (OLS) does not provide the optimal solution.Therefore, we apply a generalized least squares (GLS) method (e.g., Rawlings et al., 2001), which allows us to jointly estimate the calibration parameters along with an iterative fit of an autoregressive process AR(q) to account for correlated errors of a ng .The optimal order of the AR process is determined from the Akaike Information Criterion (Akaike et al., 1988), which provides the relative quality of various AR models used in the GLS procedure.Our numerical experiments indicate that AR(q) with q = 20 removes the correlations sufficiently and avoids overparameterization.The AR coefficients are then used to decorrelate the input data of the OLS and the calibration parameters are estimated again.In an iterative procedure, an autoregressive process is fitted to residuals of the GLS again until convergence.
The MNE procedure applied here differs from the calibration procedure of Bezděk (2010) in the way it handles autocorrelation in a ng .Bezděk (2010) applies an inverse operator of the second derivative filter on a ng as well as on accelerometer measurements to recover the orbit and estimates the calibration parameters by comparing these orbits.Applying this procedure introduces new numerical errors while converting accelerations to orbital positions.Since the application of the GLS together with the AR process has already reduced the correlation errors, we directly find the calibration parameters from Eq. (1).

Dynamic estimation (DE)
Accelerometer calibration parameters can also be estimated within a dynamic precise orbit determination (POD) procedure.
In dynamic POD, orbits are estimated from observables, while accounting for the forces acting on the satellite including the non-gravitational forces from accelerometer measurements.In our implementation, kinematic orbits are the observables.
In the variational equation approach (e.g., Tapley, 1973), the dynamic POD is estimated by solving the equation of motion (see Eq. ( 3)) together with the associated variational equations.The force model parameters p in the equation of motion consist of gravitational and non-gravitational accelerations including accelerometer calibration parameters.The partial derivatives of the equation of motion with respect to the force model parameters p can be used to build the variational equations (Löcher, 2011) In Eqs. ( 5) and ( 6), a (t, r, ṙ, p) corresponds to the equation of motion introduced in Eq. (3).The system of equations is built using the partial derivatives of the equation of motion where r(t) is the given kinematic orbit and r(t) states the computed orbit using approximate values of p and x 0 (Löcher, 2011).According to Eq. ( 7), the parameters, which include the satellite's position r, velocity ṙ, and accelerometer calibration parameters b and S, are improved iteratively in a least-squares sense.Convergence is reached as soon as the starting position changes less than 1mm, since beyond this threshold the calibration parameters will not change significantly.

Empirical model approach (EMA)
Assuming that the non-gravitational accelerations are realistically modeled, e.g., using an empirical approach, calibration of the accelerometer measurements could be done by adjusting the on-board measurements to the modeled values (Sutton et al., 2007;Doornbos, 2012).For this, a ng in Eq. ( 1) is empirically modeled as In this formulation, we consider a ng to be consisted of atmospheric drag a drag as well as accelerations due to radiation pressure of the Sun a SRP and the Earth a ERP , which are presented in Sect.3.3.1,3.3.2,and 3.3.3,respectively.These non-gravitational forces also depend on the surface properties of the satellite, which are considered here using GRACE's satellite macro-model (see Sect. 2.1.3).

Atmospheric Drag
The atmospheric drag is caused by the interaction of the particles within the atmosphere with the surface of the satellite.The impact of drag on the satellite's surface is derived from as demonstrated by e.g., Bruinsma et al. (2004); Sutton et al. (2007); Doornbos (2012).Atmospheric drag depends mainly on the neutral density of the thermosphere ρ at the position of the satellite, as well as on the coefficient C D accounting for drag and lift coefficients.The drag coefficient is estimated following Doornbos (2012, Eqs. 3.55 -3.59), which is based on the original model by Sentman (1961) including modifications by Moe and Moe (2005).During the estimation of the drag coefficient, we choose an energy accommodation coefficient of 1 to account for diffuse reflection.In our EMA implementation, we use the empirical model NRLMSISE-00 (Picone et al., 2002) to determine the thermospheric neutral density ρ in Eq. ( 9).
Moreover, it is required to know the satellite's mass m, its areas A projected onto flight direction and the relative velocity v rel modeled by the satellite's initial velocity with respect to the velocity of the atmosphere, which is assumed to rotate with the Earth.Atmospheric winds are neglected as they have only a minor impact on the satellite's velocity and the estimated thermospheric densities (Sutton, 2008;Mehta et al., 2017).Unlike the calibration techniques of Sects.3.1 and 3.2, the modeled non-gravitational acceleration derived from the EMA depends on the model used to derive density at the position of the satellite, which mitigates the physical quality of this approach.This selection also has an influence of the thermospheric neutral density derived from GRACE accelerometer data, which will be discussed in Sect.4.2.

Solar radiation pressure (SRP)
Both visible and infrared radiation of the Sun interact with the surface of LEO satellites in terms of reflection or absorption, which accelerate the satellite due to the solar radiation pressure (SRP) (e.g.Sutton, 2008;Montenbruck and Gill, 2012).
Accounting for the interaction of photons with the satellite's surface requires information on the surface material.Moreover, the constellation of the Earth, the Sun and the satellite cause changes in the satellite's illumination.Ignoring the satellite's orientation, SRP is largest when the satellite is located in the direct sunlight, whereas it is lower in penumbra and it does not exist in umbra.Accelerations due to SRP acting on the GRACE satellite are modeled as in Doornbos (2012).During radiation pressure modeling in Sects.3.3.2and 3.3.3,the reflection model in Doornbos (2012, Eq. 3.48) is used together with the surface properties provided in Tab. 1 to account for diffuse and specular reflection, as well as for absorption.Testing other, potentially more realistic bidirectional reflectance distribution functions (e.g.Ashikhmin and Shirley, 2000) remains a subject of future research.

Earth radiation pressure (ERP)
The Earth emits thermal radiation and reflects a fraction of the incoming sunlight back into the space, where both radiations interact mainly with the nadir-pointing surface of the satellite in terms of reflection or absorption.This causes an acceleration due to the Earth radiation pressure (ERP), which decreases with an increasing distance of the satellite to the Earth, and cannot be neglected in force modeling for LEO satellites.
Accelerations due to ERP acting on the GRACE satellite are usually modeled following Knocke et al. (1988), where ERP is split into albedo (visible short wavelength radiation) and emission (infrared long wavelength radiation).Equations to estimate ERP acceleration on GRACE satellites can be found in Appendix B.
After several investigations, we conclude that the polynomial fit used in the Knocke et al. (1988) Research Center CERES ordering tool (http://ceres.larc.nasa.gov/).These data are used to estimate monthly global albedo and emission fields.Monthly albedo and emission fields are converted to equivalent spherical harmonic coefficients of low degree (<20) using a numerical integration as in Forootan et al. (2013).The temporal evolution of the albedo and emission coefficients are modeled by fitting cyclic functions of sine and cosine with annual and semi-annual frequencies (see also Rodríguez-Solano, 2009).

Thermospheric neutral density estimation with calibrated accelerometer data
Finally, the calibration parameters from the above methods are applied to the raw accelerometer measurements a raw to obtain calibrated accelerometer measurements according to Eq. ( 1) as Thermospheric neutral densities based on the calibrated accelerometer data a cal are estimated following (Sutton et al., 2005, Eq. ( 6)) where the equation of atmospheric drag a drag is solved for the density and the drag force is replaced by the relation of nongravitational forces acting on the satellite according to Eq. ( 8).In Eq. ( 11), r denotes the position of the satellite.

Calibration results
In the following, we compare the calibration parameters, represented in the SRF, obtained from the three calibration procedures during three individual months of the current (24 th ) solar cycle with varying solar activity.As already mentioned in Sect.3, the calibration parameters are estimated iteratively similar to Van Helleputte et al. (2009) to avoid highly (anti-)correlated biases and scales.To keep the calibration methods comparable, constant scales are assumed for all three approaches.We use the scales derived from the Dynamic Estimation (DE), since the (anti-)correlation dominates especially in the Multi- Step Numerical Estimation (MNE) and the Empirical Model Approach (EMA) causing unrealistic scales not close to 1, most  Bezděk (2010).In the along-track and cross-track directions, the scales of both satellites resulting from this study are about 0.03 (along-track) and 0.06 (cross-track) smaller than those from other methods.This variation is likely influenced by the period of data used to estimate the calibration parameters, which is at most seven years in other studies, whereas the scales here are derived from 14 years of data.When using data between August 2002 and March 2009, the along track scale of GRACE A obtained from DE is found to be 0.951, which is close to the scale of 0.960 published by Bettadpur (2009).We conclude that the scales of accelerometer measurements vary with the length of data used as well as with the method of estimation as already stated by Bettadpur (2009).
Daily biases for GRACE A during November 2008, February 2014 and March 2015 corresponding to low, medium and high solar activity obtained from the three calibration procedures using constant scales (Tab.3) are presented in Fig. 2. The magnitude of biases (Fig. 2) is confirmed to be different for the along-track, cross-track, and radial directions.Along-track biases vary between −1.4 • 10 −6 ms −2 and −1 • 10 −6 ms −2 , where biases from EMA differ from those of DE and MNE in particular during high and medium solar activity.In the cross-track direction, biases (2.7 • 10 −5 ms −2 to 3 • 10 −5 ms −2 ) are one order of magnitude larger than in the along-track direction.Cross-track biases of the three approaches are similar until the second decimal digit.Radial biases vary between −8•10 −7 ms −2 and −4•10 −7 ms −2 .Especially in the radial direction, biases show larger variations with time while comparing to the along-track and cross-track directions.Moreover, an offset is found between the calibration parameters obtained from MNE and other approaches, however, the reason for this difference remains unclear.
Based on the numerical results, one can see that calibration parameters obtained from MNE and DE are fairly similar in the along-track and cross-track directions.The offset between along-track calibration parameters obtained from EMA, compared to other methods, is caused by the dependency of its results on the density derived from empirical models (Doornbos, 2012).This empirical approach will likely become more realistic by introducing a horizontal wind model to the calculation of the satellite's relative velocities within drag estimations (Eq.( 9)) (Sutton, 2008).Moreover, Mehta et al. (2014) found that an improved understanding of the gas-surface interactions impacts the drag coefficient and hence yields more realistic neutral density for the GRACE satellites.The results obtained in this study confirm the order of magnitude of biases in all three directions, which are derived by Bezděk (2010, Fig. 14) for 1.5 years at the beginning of the mission.
The mean and standard deviation of the calibration parameters of March 2015 are provided in Tab. 4. In general, the standard deviation of biases in the along-track direction (5.95•10 −10 ms −2 to 1.15•10 −8 ms −2 ) are smaller for MNE and DE than those of the cross-track and radial directions.This is also observed for the scale factors with a standard deviation of 0.002 in the along-track direction (see Tab. 3).The relatively larger standard deviations of MNE biases may be more realistic compared to other techniques, since we estimate the impact of autocorrelation through the generalized least squares method.In the other techniques, autocorrelation is neglected, which results in lower error estimations.

Thermospheric neutral density estimation
By applying daily biases and the constant scale factors on raw accelerometer measurements, corresponding calibrated time series are computed.The calibrated accelerometer measurements obtained from the three applied calibration procedures (Sects. 3.1, 3.2, and 3.3) are used to compute thermospheric neutral density profiles in the along-track direction (Sect.3.4).In addition, the derived densities are compared to two different empirical models NRLMSISE-00 (Picone et al., 2002) and Jacchia-Bowman 2008 (Bowman et al., 2008).For evaluating the accelerometer-based densities resulting from this study, density sets by Sutton (2008), and more recently estimated densities from Mehta et al. (2017), both available between 2002 and 2010, are taken into account3 .Throughout this study, the densities are not normalized to a constant altitude, which is done for example in Lei et al. (2012), since such normalization needs an assumption about vertical changes of the thermospheric density.
Daily averages of along-track densities during three particular months with high, medium and low solar activity are presented in Fig. 3. Average densities during November 2008 vary about 2 • 10 −13 kg m −3 , which is more than one order of magnitude less than those during high (3 • 10 −12 kg m −3 ) or medium (2 • 10 −12 kg m −3 ) solar activity.Independent of the solar activity, empirical densities obtained from NRLMSISE-00 coincide best with the densities calculated from EMA.This is due to the fact that the NRLMSISE-00 model is already used for estimating the calibration parameters (see Sect.During this period, the DE approach is found to be the most suitable method to calibrate accelerometer measurements in order to derive thermospheric neutral densities.Since DE derived densities are closer to those of Sutton (2008) than to the recently computed densities by Mehta et al. (2017), we confirm that an improved drag estimation as performed by (Mehta et al., 2017) affects the density estimates.For example, using an effective energy accommodation coefficient is expected to increase the densities up to 20% (Mehta et al., 2013), and a better understanding of the gas-surface interaction will lead to further improvements (Mehta et al., 2017).
In addition to daily mean densities, the along-track densities on November 1 st , 2008, are presented in Fig. 4, whose results indicate that the general pattern is the same as the pattern of densities obtained in this study.Dominant peaks with a 1.5-hour period are found to be evident that is related to the orbital period of approximately 15 revolutions per day.DE derived densities are found to be closest to Sutton (2008)   densities.We observe that DE derived densities for this day are closer to those published by Mehta et al. (2017) than those in Sutton (2008).
On March 17 th and 18 th , 2015, thermospheric densities reach a maximum due to a strong solar event (St.Patrick's Day storm).Along-track densities during these days are presented in Fig. 5. Accelerometer based densities (MNE, DE, EMA) show stronger reactions to this storm of up to 1.2 • 10 −11 kg m −3 compared to empirically modeled densities of NRLMSISE-00 (0.6 • 10 −11 kg m −3 ) and JB2008 (0.9 • 10 −11 kg m −3 ) due to variations within the atmospheric composition affecting the atmospheric density and thus accelerometer measurements.However, JB2008 modeled densities seem to better represent the solar event than NRLMSISE-00 even though the maximum is delayed.The temporal resolution of input parameters such as F 10.7 index used in the empirical density models is daily, which leads to much weaker densities compared to accelerometer based densities during the St. Patrick's day storm.Our results confirm that the densities along the orbit resulting from calibrated accelerometer measurements are more reliable than empirically modeled densities, in particular during strong solar events.

Empirical density corrections for NRLMSISE-00
In the following, the accelerometer based densities ρ cal are used to estimate empirical corrections for NRLMSISE-00 densities and decrease afterwards due to the delayed and weakened maximum in the empirical thermospheric densities (see also Fig. 5).
The empirical model NRLMSISE-00 underestimates the thermospheric neutral density with values of up to 28% of simulated densities during high solar activity.It also overestimates neutral densities, which are found to be on average 22% of simulated densities during low solar activity.Only during medium solar activity, the empirical corrections are approximately close to 1. Overestimated model densities during low solar activity have also been reported in Doornbos (2012), and we confirm that empirical density models perform well only during moderate conditions.
Since the necessity of correcting empirical thermospheric neutral density models is evident, we provide global empirical corrections on a daily base, which can be used to scale model derived neutral density estimations for the altitude of 400km.
The corrections can likely be used for the whole altitude range of 300 -600km, since the altitude-dependent changes of neutral density is close to linear.Empirical corrections on March 2 nd , 2015, along the orbit of GRACE A are presented in  In order to derive global patterns of differences between GRACE densities and model output, as GRACE does not exactly repeat its daily tracks, daily density scales (s = ρ cal ρ −1 ) are first estimated by scaling the GRACE-derived densities along its orbit by the corresponding NRLMSISE-00 neutral density simulations.These daily fields are then converted to spherical harmonics coefficients, which allows spectral filtering to retain only low degree differences between GRACE and model estimations.To estimate global empirical correction maps, we consider a trade-off which resembles the problem of estimating time-variable gravity: Longer analysis intervals (currently one day) would allow to accumulate more data, with better spatial coverage, at the expense of temporal resolution.The spherical harmonic coefficients smooth out real variability at time scales below one day (see neighboring tracks in Fig. 8), but some signal may alias into the daily estimates.Shorter analysis intervals would enable better temporal resolution, but with less dense data coverage and thus a loss of spatial resolution.We find, in- deed, daily analysis an appropriate balance between temporal and spatial resolution.However, we have not carried out a formal optimization, since it was out of the scope of the present paper.
From the average cross-track spacing, we found that a spherical harmonic degree of 11 could be resolved, i.e. fixing the maximum degree to 10 is appropriate.The daily spherical harmonic coefficients are estimated using a least squares estimation.
Numerical problems are not expected, since the analysis of the normal equation matrix yields a condition number below 10 3 .
The scale correction fields are then synthesized on a 1 • × 1 • grid using the coefficients of up to degree and order 7, which are presented on the right column in Figure 8.The degree variance with a minimum of degree 7 suggests a synthesis of the same degree and an analysis of slightly higher degrees to avoid a loss of information.The pattern in the global corrections resulting from the three accelerometer calibration approaches is generally similar, however, the method of calibrating accelerometer data has a visible impact on the resulting densities.The global empirical corrections provide a measure of how much accelerometer data can improve empirically modeled densities obtained from NRLMSISE-00.According to the along track empirical corrections in the left column, the global expansions (right column) indicate that NRLMSISE-00 neutral density simulations need to be corrected by 22% on average during March 2 nd , 2015.

Conclusions
In this study, the measurements of the SuperSTAR accelerometers on-board the GRACE satellites are calibrated using three procedures.The Multi-Step Numerical Estimation approach is based on the numerical differentiation of kinematic orbits, where the main challenges are the noise amplification and the temporal correlation after applying a numerical differentiation operator.Here, similar to Bezděk (2010), the Savitzky-Golay filter is applied to mitigate the impact of noise and the temporal autocorrelation is reduced by fitting an autoregressive process within the generalized least squares estimation of the calibration parameters.In the Dynamic Estimation approach, the accelerometer measurements are calibrated within a dynamic POD based on the variational equation approach (Löcher, 2011).Finally, we apply an Empirical Model Approach, where the non- Furthermore, thermospheric neutral densities derived from calibrated accelerometer measurements along-track of GRACE are compared to densities obtained from the empirical models NRLMSISE-00 and Jacchia-Bowman 2008.The results suggest that accelerometer derived densities provide more reliable results, especially on short time scales and during strong solar events, for example during the St. Patrick's Day storm on March 17 th , 2015.Hence, accelerometer derived densities allow the improvement of empirical density models as already stated by Doornbos (2012), and ways to integrate these while retaining high temporal resolution should be found, e.g. by estimating 24-hour empirical correction fields at hourly or more frequent intervals.We conclude, from comparisons with densities from Sutton (2008) and recent results from Mehta et al. (2017), that densities estimated using the Dynamic Estimation fit better than those of Sutton (2008) but not as good as those of Mehta et al. (2017).
suggest that it is necessary to apply corrections to model densities depending on the solar activity.Due to overestimated empirical model densities during low solar activity, empirical corrections of 22% on average need to be applied on NRLMSISE-00 densities during quiet periods.In contrast, empirical corrections of up to 28% are required during high solar activity, since the model underestimates neutral densities.The spherical harmonic expansion of these corrections on a global grid provides a measure indicating to what extent GRACE-derived thermospheric density estimation can improve simulations of empirical density models on a daily base.These findings encourage the use of these factors to improve empirical density models.
Further efforts in satellite drag modeling will improve the Empirical Model Approach to calibrate accelerometer measurements, as well as the thermospheric neutral densities estimated from the three methods.Moreover, including a horizontal wind model in the Empirical Model Approach is expected to yield more realistic densities which might improve the consistency of the results.The Multi-Step Numerical Estimation method may be further developed through modeling the temporal correlations of accelerometer measurements.
In further studies, the empricial corrections derived from calibrated accelerometer measurements of GRACE A could be used to model densities in order to simulate non-gravitational accelerations acting on GRACE B, which contributes to filling data gaps during months where only one satellite provides accelerometer measurements.Other methods on transferring nongravitational accelerations of a satellite to a co-orbiting one are discussed in Kim and Tapley (2015).
The calibration procedures are applicable to other satellite missions carrying space-borne accelerometers as well.Combining the thermospheric neutral densities derived from different calibrated accelerometers allows further improvement of empirical density models.For example, the empirical density corrections at different altitudes can be used to obtain altitude profiles to correct empirical density models, which could then be used to derive accurate drag predictions for other satellites wich are not equipped with acceleromter, restricted to the period when the corrections are available.Besides, the assimilation of calibrated accelerometer measurements of various satellite missions into physical thermosphere/ionosphere models would likely enable an improved representation of physical processes in the atmosphere, e.g., following Matsuo et al. (2012) or Fedrizzi et al. (2012). of the analytical derivative.Random white noise of 2cm according to the standard deviation of the kinematic orbits has been added to the modeled orbit.The differences of analytical derivative to three different numerical derivatives are presented in The difference between the analytical derivative and the smoothing differentiation filter (black) shows that this filter introduces an unwanted phase shift.Therefore, the smoothing differentiation filter is not suitable.In comparison, the Savitzky-Golay filter (red and green) prevents phase shifts, and the application of different filter settings clarifies that difference between the analytical and the numerical derivative is minimized, i.e. the amplification of noise is limited, when using the settings n = 11 and m = 5 (red).The settings recommended in Bezděk (2010) (green) vary slightly due to the use of different kinematic orbits.
P u blis h e r s p a g e : h t t p:// dx.doi.o r g/ 1 0. 5 1 9 4/ a n g e o-3 6-7 6 1-2 0 1 8 Pl e a s e n o t e: C h a n g e s m a d e a s a r e s ul t of p u blis hi n g p r o c e s s e s s u c h a s c o py-e di ti n g, fo r m a t ti n g a n d p a g e n u m b e r s m a y n o t b e r efl e c t e d in t hi s v e r sio n.Fo r t h e d efi nitiv e v e r sio n of t hi s p u blic a tio n, pl e a s e r ef e r t o t h e p u blis h e d s o u r c e .You a r e a d vis e d t o c o n s ul t t h e p u blis h e r's v e r sio n if yo u wis h t o ci t e t hi s p a p er.This ve r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wi t h p u blis h e r p olici e s. S e e h t t p://o r c a .cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s.Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Figure 1 .
Figure 1.Satellite body-fixed reference frames (modified from Bettadpur, 2012).Here, SRF represents the science reference frame, AF and SF respectively stands for the accelerometer frame and the satellite frame, and SCF indicates the star camera frame.

3
Methods of calibrating accelerometer measurementsCalibration of the accelerometer measurements a raw requires an equation to link non-gravitational accelerations a ng with a set of calibration parameters.The parameterization is commonly formulated to estimate daily biases b = [b x , b y , b z ] T , and scale factors S = diag[s x , s y , s z ] for the x-, y-, and z-directions, respectively (e.g.Bettadpur, 2009;Van Helleputte et al.,

Figure 2 .
Figure 2. Comparison of daily biases of acceleration measurements of GRACE A during low (November 2008), high (February 2014), and medium (March 2015) solar activity.In these figures, MNE represents the Multi-Step Numerical Estimation (blue), DE indicates the Dynamic Estimation (red), and EMA stands for the Empirical Model Approach (green).Results are presented in the along-track (x), cross-track (y) and radial (z) directions of the SRF.Note that the scale of the y-axis differs due to varying biases along the three axes.

Figure 3 .
Figure 3. Daily average of along-track densities during November 2008 (low), February 2014 (medium) and March 2015 (high solar activity).Density obtained from calibrated accelerometer measurements: Multi-Step Numerical Estimation (MNE, blue), Dynamic Estimation (DE, red) and Empirical Model Approach (EMA, green).Density obtained from empirical density models are shown as NRLMSISE-00 (MSIS, black) and Jacchia-Bowman 2008 (JB, yellow).Density sets by Mehta et al. (2017) (light red) and Sutton (2008) (grey) are available in 2008.Note that the scale of the y-axis differs due to a strong variation of solar activity.
in terms of scale factors s = ρ cal ρ −1 similar toDoornbos et al. (2009).Mean empirical corrections during November 2008, February 2014 and March 2015 along the orbit of GRACE A are presented in Fig.6.Since densities derived from EMA coincide with NRLMSISE-00 densities (see Fig.3), empirical corrections are not estimated for EMA.In November 2008, the MNE based empirical corrections are less reliable due to less stable calibration parameters resulting from this approach during periods of a low signal-to-noise ratio.Due to the similarity of the calibration parameters derived from MNE and DE during high and medium solar activity, the density corrections obtained from both approaches are found to be similar with mean values of 1.11 (MNE) and 1.12 (DE) in March 2015, and in February 2014 the mean values are found to be 0.99 (MNE) and 0.97 (DE).During the St. Patrick's day storm (March 17 th and 18 th , 2015), the mean empirical corrections increase on March 17 th ,

Fig. 7 .
Fig. 7. Due to the similarity of the calibration parameters derived from MNE and DE, the density corrections obtained from both approaches are found to be similar on this day with a mean value of 1.21 (MNE) and 1.24 (DE).Additionally, a spatial representation of the corrections, which are estimated for the NRLMSISE-00 empirical model along the orbit of GRACE A on March 2 nd , 2015, are shown on the left column of Fig. 8.The empirical corrections derived from MNE and DE indicate that NRLMSISE-00 neutral densities need to be increased by 22% on average on this day.

Figure 6 .
Figure 6.NRLMSISE-00 daily mean empirical corrections of GRACE A during November 2008 (low), February 2014 (medium) and March 2015 (high solar activity).Densities obtained from calibrated accelerometer measurements of the Multi-Step Numerical Estimation (MNE, blue) and the Dynamic Estimation (DE, red).

Figure 7 .
Figure 7. NRLMSISE-00 empirical corrections along the orbit of GRACE A on March 2 nd , 2015.Densities obtained from calibrated accelerometer measurements of the Multi-Step Numerical Estimation (MNE, blue), the Dynamic Estimation (DE, red).

Figure 8 .
Figure 8. NRLMSISE-00 empirical corrections of GRACE A during March 2 nd , 2015.Left column: Corrections along the orbit.Right column: Corrections expanded on a global grid using a least squares adjustment to obtain spherical harmonic coefficients up to degree and order 10 (synthesized up to degree and order n = 7).Densities obtained from calibrated accelerometer measurements of the Multi-Step Numerical Estimation (MNE, top), the Dynamic Estimation (DE, bottom).

Table 1 .
(Bettadpur, 2012)s of the GRACE macro-model(Bettadpur, 2012)including coefficients for specular and diffuse reflection of visible and infrared radiation.the composition of the atmosphere provided by the Solar Maximum Mission.Both models are forced by a number of parameters, e.g. the geomagnetic planetary index K p accounting for variations of geomagnetic activity, and the F 10.7 index

Table 2 .
Gravitational force models and the partial derivatives of the equation of motion with respect to the initial state x 0 = [r 0 ṙ0 ] T as (Loeb et al., 2018)009)h is based on 48 monthly mean Earth radiation budget maps derived from satellite observations until 1979, does not sufficiently represent the spatial variability of albedo/emission (see alsoRodríguez-Solano, 2009).Hence a new ERP model is developed here, which considers spherical harmonics fit to remote sensing observations, including long wavelength and short wavelength flux at the top of atmosphere, as well as incoming solar radiation.As input data, we use the latest version of Cloud and the Earth's Radiant Energy System (CERES) data CERES EBAF_Ed4.0(Loebetal., 2018)obtained from the NASA Langley

Table 3 .
Scale factors S of accelerometer measurements of GRACE A and B obtained from different methods in the along-track (x), crosstrack (y) and radial (z) directions of the SRF.cross-track and radial directions (not shown).The mean and standard deviation of the scales of accelerometer measurements derived from DE using data between August 2002 and July 2016 are provided together with similar statistics from other studies for the GRACE satellites in Tab. 3. In general, the estimated scales for GRACE accelerometer are close to 1, which is in agreement with the expected behavior.Especially in the radial direction, the scale of 0.94 (GRACE A) and 0.93 Bettadpur (2009)d from DE is similar to the scales obtained during a POD byBettadpur (2009)and during the estimation similar to MNE by