# New thermosphere neutral mass density and crosswind datasets from CHAMP, GRACE, and GRACE-FO

Christian Siemes, Claudia Borries, Sean Bruinsma, Isabel Fernandez-Gomez, Natalia Hladczuk, Josevan den Ijssel, Timothy Kodikara, Kristin Vielberg, Pieter Visser

► To cite this version:

Christian Siemes, Claudia Borries, Sean Bruinsma, Isabel Fernandez-Gomez, Natalia Hladczuk, et al. New thermosphere neutral mass density and crosswind datasets from CHAMP, GRACE, and GRACE-FO. Journal of Space Weather and Space Climate, 2023, 13, pp.16. 10.1051/swsc/2023014. hal-04831298

HAL Id: hal-04831298

https://hal.science/hal-04831298v1

Submitted on 12 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New thermosphere neutral mass density and crosswind datasets from CHAMP, GRACE, and GRACE-FO

Christian Siemes<fnref n="1" />*, Claudia Borries<fnref n="2" /> ID, Sean Bruinsma<fnref n="3" /> ID, Isabel Fernandez-Gomez<fnref n="2" /> ID,
Natalia Hładczuk<fnref n="1" />, Jose van den IJssel<fnref n="1" /> ID, Timothy Kodikara<fnref n="2" /> ID, Kristin Vielberg<fnref n="4" /> ID, and Pieter Visser<fnref n="1" /> ID

Received 3 March 2023 / Accepted 9 May 2023

**Abstract**—We present new neutral mass density and crosswind observations for the CHAMP, GRACE, and GRACE-FO missions, filling the last gaps in our database of accelerometer-derived thermosphere observations. For consistency, we processed the data over the entire lifetime of these missions, noting that the results for GRACE in 2011–2017 and GRACE-FO are entirely new. All accelerometer data are newly calibrated. We modeled the temperature-induced bias variations for the GRACE accelerometer data to counter the detrimental effects of the accelerometer thermal control deactivation in April 2011. Further, we developed a new radiation pressure model, which uses ray tracing to account for shadowing and multiple reflections and calculates the satellite’s thermal emissions based on the illumination history. The advances in calibration and radiation pressure modeling are essential when the radiation pressure acceleration is significant compared to the aerodynamic one above 450 km altitude during low solar activity, where the GRACE and GRACE-FO satellites spent a considerable fraction of their mission lifetime. The mean of the new density observations changes only marginally, but their standard deviation shows a substantial reduction compared to thermosphere models, up to 15% for GRACE in 2009. The mean and standard deviation of the new GRACE-FO density observations are in good agreement with the GRACE observations. The GRACE and CHAMP crosswind observations agree well with the physics-based TIE-GCM winds, particularly the polar wind patterns. The mean observed crosswind is a few tens of m·s⁻<fnref n="1" /> larger than the model one, which we attribute primarily to the crosswind errors being positive by the definition of the retrieval algorithm. The correlation between observed and model crosswind is about 60%, except for GRACE in 2004–2011 when the signal was too small to retrieve crosswinds reliably.

**Keywords:** Thermosphere / Neutral mass density observations / Neutral wind observations / Accelerometer data calibration / Radiation pressure modeling

## 1 Introduction

Accurate knowledge of the atmospheric density is critical for orbit predictions at altitudes below 600 km (Vallado & Finkleman, 2014) and for studying the coupling between the thermosphere and the ionosphere (Palmroth et al., 2021). Neutral mass density observations derived from precise accelerometer measurements on board satellites are the only current in-situ data source in the altitude range of 200–550 km with a high resolution along the orbit (Bruinsma et al., 2022). The accelerometer measurements allow deriving crosswind when the aerodynamic acceleration is large enough

(Sutton et al., 2007; Doornbos et al., 2010), which gives additional information on the dynamics of the thermosphere.

Since 2000, several satellites carrying precise accelerometers have been launched. The first was the CHallenging Mini satellite Payload (CHAMP) satellite, which was launched in 2000 and re-entered in 2010 (Reigber et al., 2002). The Gravity Recovery And Climate Experiment (GRACE) twin satellites followed in 2002 and provided data until decommissioning in October 2017, although they remained in orbit for another few months (Tapley et al., 2004). Then, the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite was launched in 2009 and re-entered in November 2013 (Floberghagen et al., 2011), the same month in which the three Swarm satellites were launched (Olsen et al., 2013). Finally, the twin

<footnote n="1">¹ Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629HS Delft, The Netherlands</footnote>

<footnote n="2">² Institute for Solar-Terrestrial Physics, German Aerospace Center, Kalkhorstweg 53, 17235 Neustrelitz, Germany</footnote>

<footnote n="3">³ Space Geodesy Office, CNES, 18 avenue Edouard Belin, 31401 Toulouse, France</footnote>

<footnote n="4">⁴ Rheinische Friedrich-Wilhelms-Universität Bonn, Nussallee 17, 53115 Bonn, Germany</footnote>

<footnote>*Corresponding author: c.siemes@tudelft.nl</footnote>
GRACE-Follow-On (GRACE-FO) satellites were launched in 2018 and provide measurements until today. All of these satellites were launched into circular, near-polar orbits with an initial altitude of roughly 500 km, except for the GOCE satellite, which was in a circular, sun-synchronous orbit at dawn/dusk local time at an altitude of 256 km for the largest fraction of its lifetime.

Early CHAMP and GRACE thermosphere density and crosswind datasets were produced by Bruinsma et al. (2004), Sutton (2008), and Doornbos (2011). They relied on panel models to derive the aerodynamic coefficient vector and evaluate the radiation pressure. More recently, Mehta et al. (2017) provided CHAMP and GRACE neutral mass density observations until 2010, and March et al. (2021) released new versions of the CHAMP and GRACE neutral mass density observations until 2009. Both authors used detailed geometry models of the satellites to simulate the aerodynamic coefficient vector by direct simulation Monte Carlo methods (Bird, 1994). The GOCE data was processed similarly by Doornbos et al. (2014) to obtain neutral mass density and crosswind observations. The Swarm satellites are a special case because the accelerometer measurements suffer from perturbations that require a substantial effort to correct (Siemes et al., 2016). Nevertheless, Swarm C neutral mass density observations have been produced with a nearly complete temporal coverage until today, and, recently, also Swarm A neutral mass density observations for 2014 have been released (Iorfida et al., 2022). Though some authors presented GRACE-FO neutral mass density observations, e.g., Krauss et al. (2020), we are unaware of systematically produced GRACE-FO neutral mass density datasets.

All these datasets have been produced with different methods for (a) the calibration of the accelerometer data, (b) the modeling of the radiation pressure to extract the aerodynamic acceleration, (c) accounting or not accounting for thermosphere wind in the neutral mass density retrieval, and (d) different assumptions on the gas-surface interaction, which affects the aerodynamic coefficient vector modeling and, hence, the scale of the neutral mass density observations. We aim to establish a database of thermosphere density and crosswind observations that is as complete, consistent, and accurate as possible. Such a database is the basis for thermosphere models like the Drag Temperature Model 2020 (Bruinsma & Boniface, 2021) and investigations of the gas-surface interaction (Bernstein & Pilinski, 2022).

This paper focuses on two critical aspects of neutral mass density and crosswind processing: the calibration of the accelerometer data and the modeling of radiation pressure. The neutral mass density and crosswind observations are derived from the along-track and cross-track accelerations, respectively (Doornbos, 2011). Therefore, the accuracy of the accelerometer data calibration and the radiation pressure modeling directly affects the accuracy of the neutral mass density and crosswind observations.

The data from electrostatic accelerometers are affected by biases and small instrument imperfections, which include a deviation of the scale factors from unity (Touboul et al., 2012, 2016). Thus, the accelerometer data need to be calibrated, and we refer to Kang et al. (2003) and Vielberg et al. (2018) for comparing the different calibration methods. The ones based on precise orbit determination (POD) give the best calibration results for the acceleration in the flight direction, although they perform less well in the other directions. van Helleputte et al.

(2009) reported an accuracy of the GRACE along-track and cross-track biases of 0.5 nm·s⁻² and 22 nm·s⁻², respectively. Further, the reported accuracy of along-track and cross-track scale factors is 1–3% and about 15%, respectively. The accuracy of the cross-track bias is of concern for the crosswind observations since Doornbos et al. (2010) demonstrated that a 10 nm·s⁻² cross-track acceleration error leads to a crosswind error of about 1 km·s⁻¹, which is larger than the expected average crosswind signal (Drob et al., 2015). A further point of attention is the thermally-induced bias variation of the GRACE accelerometer data, which is significant after the deactivation of the accelerometer thermal control in April 2011 (Klinger & Mayer-Gürr, 2016).

The size of the radiation pressure acceleration exceeds the size of the aerodynamic acceleration in the along-track direction when a satellite orbits at and above altitudes of 500 km during solar minimum, which was shown for the Swarm satellites by Van den IJssel et al. (2020). Since the GRACE and GRACE-FO satellites spent several years at such high altitudes during solar minimum, the accuracy of the radiation pressure modeling is essential for the neutral mass density calculation. Because the aerodynamic acceleration in the cross-track direction is much smaller than the one in the along-track direction, we expect that the accuracy of the radiation pressure modeling is relevant for the crosswind observations throughout the entire mission lifetime.

The key elements of the radiation pressure modeling are the outer satellite geometry, the satellite surface properties, the radiation flux, and the satellite thermal emissions. Solar radiation flux is near constant at a fixed distance to the Sun, with variations of only 0.1% (Dewitte & Clerbaux, 2017). Earth albedo and infrared radiation contribute less than 1 nm·s⁻² to the along-track and cross-track accelerations, which was shown by Vielberg & Kusche (2020) for the GRACE satellites. This is possibly a consequence of the symmetric satellite shape and Earth-pointing satellite attitude. Since the situation is very similar for the CHAMP and GRACE-FO satellites, we do not expect to gain much accuracy for the neutral mass density and crosswind observations by improving the radiation flux modeling.

The logical next step for advancing our processing is to augment the detailed satellite geometry models from March et al. (2019a) with thermo-optical surface properties to replace the presently used panel models. These models can then be used in ray-tracing simulations to determine the radiation pressure force coefficients. The key features of ray-tracing techniques for radiation pressure modeling have already been described in the 1990s (Klinkrad et al., 1991), which includes shadowing and multiple reflections.

We have two options for the implementation of the ray-tracing simulations. The first is to perform the simulations for each position along the orbit, where the calculations are accelerated by the graphics processing unit as suggested by Kenneally & Schaub (2020). Its advantage is the simplicity of calculating the amount of absorbed radiation per surface element, which we could use as an input to calculate the satellite's thermal emissions and account for moving parts such as steerable solar arrays. The disadvantage is the need to perform one simulation per incidence angle of incoming radiation, which quickly escalates when considering Earth's albedo and infrared radiation. More importantly, we intend to process the entire
CHAMP, GRACE, and GRACE-FO mission data, which renders this option infeasible due to the very long runtime of the calculations. Therefore, we perform the ray-tracing simulations in advance and store the results in a lookup table. Such an approach has been adopted by the University College London, which developed a tool for radiation pressure modeling based on ray tracing (Bhattarai et al., 2019). It has been extensively used for radiation pressure modeling within the context of precise orbit determination of GNSS satellites (Ziebart, 2004; Li et al., 2020; Bhattarai et al., 2022). Similarly, Darugna et al. (2018) used ray tracing to evaluate the radiation pressure for the first satellite of the Japanese Quasi-Zenith Satellite System, QZS-1. These methods have in common that they use detailed geometry models and generate grids of rays that illuminate the satellite, where one grid is generated for each incidence angle. Another example is provided by Wöske et al. (2019), who used ray tracing with a detailed, finite element model of the GRACE satellites to precalculate the radiation pressure for all incidence angles.

Wöske et al. (2019), as well as Vielberg & Kusche (2020), demonstrated for the GRACE satellites that the satellite thermal emissions are the second largest contributor to the radiation pressure acceleration, only surpassed by the solar radiation pressure. The thermal emissions are more substantial in the cross-track and radial directions due to the solar arrays on the sides, the satellites' top panels, and the radiator panels on the bottom. Thus, we need to model the thermal emissions accurately to obtain usable crosswind observations. Montenbruck et al. (2015) model thermal emissions as an instant reradiation of the absorbed radiation. However, the results of Wöske et al. (2019) show that thermal inertia is significant for the GRACE satellites. Therefore, we need to account for the history of absorbed radiation, implying that the assumption of an instant reradiation is inaccurate.

Though the methods to model radiation pressure have advanced significantly and detailed satellite geometry models were produced recently, our knowledge of the thermo-optical surface properties and thermal properties is not accurate for the satellites under consideration. For instance, the reflection and absorption coefficients for the GRACE and GRACE-FO satellites are specified identically (Bettadpur, 2012; Wen et al., 2019). In contrast, photos of the satellites clearly show that at least the solar arrays are different. Thus, empirically updating these properties based on the in-flight data was necessary.

The paper is structured as follows. Section 2 describes the methods we used to derive the neutral mass density and crosswind observations. We emphasize the accelerometer data calibration and the radiation pressure modeling, presented in Sections 2.1 and 2.4, respectively. The input data is specified in Section 3. We present and discuss the results of the accelerometer data calibration, radiation pressure modeling, and neutral mass density and crosswind retrieval in Section 4. We conclude with Section 5, which also provides an outlook for future improvements.

## 2 Methods

The method for deriving neutral mass density and crosswind observations is described in detail by Doornbos (2011), which we briefly summarize here. First, we preprocess the accelerometer data to remove the effects of thruster activations

and other artificial accelerations generated by the satellite, e.g., “twangs” in the case of the GRACE mission (Flury et al., 2008). The acceleration due to the thruster activations is stored since it is a perturbation when deriving density but a signal when estimating accelerometer data calibration parameters as described in Section 2.1.1.

Then, the measured acceleration vector $a_{meas}$ transforms into the calibrated acceleration vector $a_{cal}$ by

$$a_{cal} = S a_{meas} + b, \qquad (1)$$

where S is a diagonal scaling matrix and b is the bias vector. The calibrated acceleration vector is the sum of the aerodynamic acceleration vector $a_{aero}$ and the radiation pressure acceleration vector $a_{rp}$, i.e.

$$a_{cal} = a_{aero} + a_{rp}. \qquad (2)$$

It should be noted that this equation holds only when the accelerometer is placed into the satellite's center of mass so that accelerations due to rigid body dynamics and gravity gradients are zero by design. This is the case for CHAMP, GRACE, and GRACE-FO satellites.

The radiation pressure acceleration needs to be modeled to extract the aerodynamic acceleration. Therefore, we decompose the radiation pressure acceleration vector into several contributions:

$$a_{rp} = a_{srp} + a_{alb} + a_{ir} + a_{te}, \qquad (3)$$

where $a_{srp}$, $a_{alb}$, $a_{ir}$, and $a_{te}$ are the acceleration vectors due to solar radiation pressure, Earth's albedo, Earth's infrared radiation, and the satellite's thermal emission, respectively.

The aerodynamic acceleration vector is related to neutral mass density $\rho$ by

$$a_{aero} = \frac{v^2 \rho}{2m} C_{aero}, \qquad (4)$$

where v is the velocity of the atmosphere relative to the satellite, m is the satellite mass, and $C_{aero}$ is the aerodynamic coefficient vector, which is intrinsically multiplied by the cross-section area. The relative velocity v is the Euclidian norm of the relative velocity vector v, which is the sum of the satellite velocity vector $v_{sat}$, the velocity vector of the corotating atmosphere $v_{atm}$, and the wind vector $v_{wnd}$:

$$v = v_{sat} + v_{atm} + v_{wnd}. \qquad (5)$$

In a nominal situation, the satellite's x-axis is aligned with the flight direction within a few degrees so that the largest fraction of the drag signal is contained in the x-component of $a_{aero}$ (the orientation of the satellite reference frame is shown in Fig. 1). Further, the x-component is the most accurately calibrated component. Therefore, we derive the neutral mass density observations by projecting equation (4) onto the x-axis of the satellite reference frame and solving for $\rho$, which gives

$$\rho = \frac{2m}{v^2} \frac{a_{aero,x}}{C_{aero,x}}. \qquad (6)$$

The crosswind observations are also derived from equation (4), where the aerodynamic coefficient vector $C_{aero}$ is a function of the relative velocity v defined in equation (5). The algorithm of (Doornbos et al., 2010) retrieves crosswind observations by
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometry models of the CHAMP satellite (left), the GRACE satellites (center), and the GRACE-FO satellites (right). The colors correspond to the surface properties and the coordinate axes indicate the orientation of the satellite reference frames.</figcaption></figure>

<table><caption>**Table 1.** Gravitational force models used for accelerometer calibration.</caption><thead><tr><td>Force</td><td>Model</td></tr></thead><tbody><tr><td>Gravity field</td><td>GOCO05s to degree/order 200 including secular degree 2 terms</td></tr><tr><td>Solid Earth tides</td><td>IERS Conventions 2003 (McCarthy & Petit, 2004)</td></tr><tr><td>Third bodies</td><td>Sun, Moon, and planets (DE-405) (Standish, 1998)</td></tr><tr><td>Pole tide</td><td>IERS Conventions 2003 (McCarthy & Petit, 2004)</td></tr><tr><td>Ocean tides</td><td>GOT 4.7 20×20, update to (Ray, 1999)</td></tr><tr><td>Non-tidal atmosphere and ocean gravity field disturbances</td><td>AOD1B RL06 (Dobslaw et al., 2017)</td></tr></tbody></table>

iteratively adjusting the relative velocity vectors' directions such that the modeled aerodynamic acceleration in equation (4) aligns with the observed aerodynamic acceleration in equation (2). The observed wind vector $v_{obs}$ is the difference between the initial relative velocity vector and the adjusted one. We must interpret this difference as the (scalar) wind speed $v_{obs} = |v_{obs}|$ into the direction of the unit vector $u_{obs} = v_{obs}/v_{obs}$. Since the satellite velocity dominates the relative velocity vector, and the adjustments are small rotations about the satellite z-axis, which approximately points into the nadir direction, the observed wind direction $u_{obs}$ is roughly in the horizontal plane in the cross-track direction. Thence, we speak of crosswind observations.

We elaborate on the methods for the accelerometer data calibration and the radiation pressure modeling in Sections 2.1 and 2.4, respectively. We complement the information by describing the satellite geometry models, augmented by surface properties, in Section 2.2 and the simulation of the aerodynamic coefficient vector of the GRACE-FO satellites in Section 2.3.

## 2.1 Accelerometer data calibration

### 2.1.1 Estimation of calibration parameters

We used the method described in Visser & van den IJssel (2016) to calibrate the accelerometer data, which we employed earlier to calibrate the individual GOCE accelerometers. The method is based on dynamic orbit determination with the NASA Goddard Space Flight Center's GEODYN software (Pavlis et al., 2006), which uses an iterative Bayesian weighted least-squares method to solve the linearized observation equations. Table 1 gives an overview of the comprehensive gravitational force models. Since the accelerometer observations represent the total non-gravitational accelerations, we exclude non-gravitational models from the estimation. We processed the orbits in daily 30-h arcs, centered on a day with a 3-h overlap at the day boundaries, and based the satellite orientation on the star tracker observations. For efficient processing, we use precise orbits resulting from a GPS-based POD as input observations instead of the actual GPS observations.

To avoid calibration terms absorbing part of the aerodynamic signal, we estimate the limited set of calibration parameters in equation (1) in addition to the initial position and velocity, resulting in an estimation close to a purely dynamic orbit determination. Constant daily biases do not represent the sizable temperature-induced bias variations well. Therefore, we remove the thermally-induced bias variations using the accelerometer's temperature data as detailed in Section 2.1.2, after which a constant daily bias parameterization is suitable again. Since the z-component of the CHAMP accelerometer data is not usable (Perosanz et al., 2005), we estimate a set of constant and one-cycle-per-revolution empirical accelerations in this direction instead of using the actual measurements.

To reduce the high correlation between the scale factor and bias parameters, we use a multi-step approach for the calibration, based on a stacked normal matrices approach (Van Helleputte, 2011). The first step determines unconstrained daily calibration parameters for the entire mission. Although the estimated scale factors are very noisy, we identify periods in which they show stable behavior. For all missions, we expect hardly any variation of the scale factors over time. Therefore, we calculate the mean and standard deviation of the time series of daily scale factors, remove all scale factors from the time series that deviate by more than three times the standard deviation from the mean, and recalculate the mean scale factor in each direction. In the next step, we constrain the daily scale factors to these values and estimate unconstrained daily biases. In that estimation process, we store the part of the normal equations related to the scale factors. Then, we sum these parts of the normal equations to estimate new unconstrained scale factors (one set for the entire time series). This method has the advantage that days with a strong acceleration signal contribute more to the estimation. In the final step, we estimate constant daily biases using these newly estimated scale factor values. If necessary, we apply additional constraints to the bias estimation.

### 2.1.2 Modeling thermally-induced bias variations

The GRACE accelerometer biases are subject to a significant variation that correlates with the instrument's temperature.
McGirr et al. (2022) used the lowpass filtered cross-track acceleration as a proxy for the along-track bias variation. While that approach is effective in the gravity field model estimation context, it is not suited for our purposes since it implicitly removes low-frequency signals in the cross-track acceleration and, hence, the crosswind observations. Therefore, we developed a correction that accurately models the acceleration’s delayed response to changes in the measured temperature (Klinger & Mayer-Gürr, 2016). The reason for the delayed response is the thermistor’s location outside the sensor structure, whereas the bias is sensitive to the temperature inside the sensor structure. To model this delay, we assume a discretized equation for radiative heat transfer

$$U(t + \Delta t) = U(t) + \Delta t \cdot k_U (T^4(t) - U^4(t)), \quad (7)$$

where T is the absolute temperature measured by the thermistor on the sensor unit (variable TSU_Yp contained in the AHK1B files) and U is the modeled absolute temperature. 0 is the time step of the discretization, and $k_u$ is the model parameter that controls the delay in the response. To achieve a good fit, we assume a second heat transfer

$$V(t + \Delta t) = V(t) + \Delta t \cdot k_V (T^4(t) - V^4(t)), \quad (8)$$

where V is the second modeled absolute temperature. The thermally induced bias variations are

$$b_T(t) = S_U U(t) + S_V V(t), \quad (9)$$

where $S_U$ and $S_V$ are diagonal matrices with temperature sensitivity coefficients for the three accelerometer axes. Since we correct for the thermally induced bias variations prior to the accelerometer data calibration, equation (1) expands to

$$a_{cal} = S(a_{meas} + b_T) + b. \quad (10)$$

### 2.1.3 Cross-track acceleration data calibration

As discussed in the introduction, the accuracy of the calibration of the cross-track acceleration is too low to derive meaningful crosswind observations. Therefore, we transferred the high accuracy of the along-track to the cross-track acceleration as follows. From solving the x-component of equation (4) for neutral mass density and inserting the result into the y-component follows:

$$a_{aero,y} = \frac{C_{aero,y}}{C_{aero,x}} a_{aero,x}. \quad (11)$$

We used this equation to calculate the cross-track acceleration from the along-track acceleration and our knowledge of the aerodynamic coefficient vector. However, equation (11) requires knowledge of the crosswind because the aerodynamic coefficient vector is a function of the relative velocity vector, $C_{aero}(v)$. Changing the wind vector $v_{wd}$ changes the relative velocity vector $v$ and, hence, the aerodynamic coefficient vector. We used the Horizontal Wind Model 2007 (HWM-2007) Drob et al. (2008) to obtain the wind along the x- and z-axes of the satellite reference frame and set the wind along the y-axis to zero to avoid a biasing of the wind observations toward the model. With these assumptions about the wind, we first calculated the aerodynamic coefficient vector and then

$$\Delta a_{aero,y} = a_{aero,y} - \frac{C_{aero,y}}{C_{aero,x}} a_{aero,x}, \quad (12)$$

which is not zero, mainly because of the inaccurate y-axis bias in equation (10). Finally, we applied a lowpass filter with a cut-off frequency of 1 d⁻¹ to $\Delta a_{aero,y}$ to obtain a more accurate estimate of the bias, noting that the choice of the cut-off frequency is linked to the daily bias estimation.

To test the impact of setting the wind along the y-axis to zero, we calculated the aerodynamic coefficient vector using the HWM-2007 for the y-axis for GRACE A in 2014. Even though 2014 coincides with a solar maximum, it is a year with moderate geomagnetic activity, noting that the latter drives HWM-2007 while solar activity does not. Therefore, we selected it as a representative test case. After lowpass filtering, we found a root-mean-square (RMS) between the original and recalculated $\Delta a_{aero,y}$ of 0.15 nm·s⁻², which demonstrates that the assumption of zero wind along the y-axis has a negligible influence on the bias estimation. Before lowpass filtering, we found an RMS of 3.6 nm·s⁻², which means that $a_{aero,y}$ contains a significant signal due to crosswind at periods shorter than one day.

## 2.2 Satellite geometry models

The basis for the radiation pressure and aerodynamic simulations are the detailed geometry models for the CHAMP and GRACE satellites created by March et al. (2019a). The models for the CHAMP and GRACE satellites consist of 2248 and 5640 triangular facets, respectively. The model for the GRACE-FO satellites is new, consists of 23,746 facets, and was created specifically for this paper.

We augmented the models with thermo-optical surface properties for the radiation pressure modeling, as illustrated in Figure 1. In the augmentation process, we used the absorption, specular, and diffuse reflection coefficients for the visible and infrared wavelengths from the missions’ documentation (Bettadpur, 2012; Wen et al., 2019). We used the material properties provided by Fortescue et al. (2011) for the parts of the satellites for which the material properties are not specified in the documentation, e.g., the antennas.

## 2.3 Aerodynamic coefficient vector simulation

The simulation of the aerodynamic coefficient vectors using the detailed geometry model of the CHAMP satellite is described in detail by March et al. (2019b). Assuming diffuse reflections with incomplete energy accommodation as the model for gas-surface interaction (Sentman, 1961), the authors found an energy accommodation coefficient of 0.85 to be a suitable choice. This value was later adopted for the GRACE A satellite (March et al., 2021). For this paper, we used the same approach and value of the energy accommodation coefficient for the simulation of the aerodynamic coefficient vector of the GRACE-FO satellites.

## 2.4 Radiation pressure modeling

### 2.4.1 Radiation pressure coefficient vector simulation

We formulate the radiation pressure acceleration

$$a_{src}(\alpha, \beta) = \frac{P_{src}(\alpha, \beta)}{m} C_{w}(\alpha, \beta) \quad (13)$$
as the product of the radiation pressure coefficient vector $C_w$ and the radiation pressure $P_{src}$, divided by the satellite mass $m$. Subscript “src” denotes the external radiation sources in equation (3), i.e., solar radiation, Earth’s albedo, and Earth’s infrared radiation, and subscript “w” indicates the wavelength, which corresponds either to visible light or infrared radiation. The direction of incidence of the radiation is described by the angles $\alpha$ and $\beta$. They can be calculated from the unit vector

$$ \mathbf{u}_{\text{sat}} = \begin{bmatrix} u_x \\ u_y \\ u_z \end{bmatrix} = \begin{bmatrix} -\cos \alpha \cos \beta \\ \cos \alpha \sin \beta \\ \sin \alpha \end{bmatrix} \quad (14) $$

that is defined in the satellite reference frame and points from the radiation source to the satellite, i.e., in the direction of the ray propagation, by

$$ \alpha = \arcsin(u_z) \quad (15) $$

and

$$ \beta = \arctan_2(u_y, -u_x). \quad (16) $$

The definition of the angles follows from the setup of the ray tracing simulation illustrated in Figure 2, where the transformation matrix

$$ \begin{align} \mathbf{R}_{\text{sim-sat}} &= \begin{bmatrix} \cos \beta & \sin \beta & 0 \\ -\sin \beta & \cos \beta & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos \alpha & 0 & \sin \alpha \\ 0 & 1 & 0 \\ -\sin \alpha & 0 & \cos \alpha \end{bmatrix} \nonumber \\ &= \begin{bmatrix} \cos \alpha \cos \beta & \sin \beta & \sin \alpha \cos \beta \\ -\cos \alpha \sin \beta & \cos \beta & -\sin \alpha \sin \beta \\ -\sin \alpha & 0 & \cos \alpha \end{bmatrix} \tag{17} \end{align} $$

relates the reference frame of the simulation and the satellite. The rays are emitted in the direction

$$ \mathbf{u}_{\text{sim}} = \begin{bmatrix} -1 \\ 0 \\ 0 \end{bmatrix}, \quad (18) $$

such that $\mathbf{u}_{\text{sat}} = \mathbf{R}_{\text{sim-sat}}\mathbf{u}_{\text{sim}}$.
In the case of the solar radiation pressure, we have

$$ P_{\text{srp}} = \left( \frac{1\text{AU}}{r} \right)^2 P_{\text{1AU}} \quad (19) $$

where $r$ is the distance between the Sun and the satellite and $P_{1\text{AU}} = 4.56 \times 10^{-6}$ is the solar radiation pressure at a distance of 1 astronomical unit (AU), i.e., the mean distance between the Earth and the Sun (Montenbruck & Gill, 2012). We refer to Vielberg & Kusche (2020) for calculating the radiation pressure due to Earth’s albedo and infrared radiation. The advantage of the formulation in equation (13) is that we can use the same radiation pressure coefficient vectors for solar radiation pressure and Earth’s albedo.

We calculate the radiation pressure coefficient vector $C_w(\alpha, \beta)$ in a ray tracing simulation, where our method is very similar to the one described by Bhattarai et al. (2019). The simulation assumes that a radiation source emits parallel rays in the direction of the satellite and that the radiation source is larger than the satellite. With these assumptions, it is straightforward to simulate the radiation pressure that the source exerts on the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Grid configuration for the ray tracing simulation. The orientation of the grid’s axes is indicated by +X, +Y, and +Z along the edges of the box. The grid’s origin coincides with the satellite center of mass. The red line indicates a ray emitted from one of the grid cells on the +X side of the simulation box. The yellow patch indicates the cell area that is associated with the ray. The angles $\alpha$ and $\beta$ describe the satellite’s orientation to the radiation source.</figcaption></figure>

satellite. Conceptually, we place the satellite into a rectangular simulation box, which is just large enough to contain the satellite. Then, we subdivide one side of the box into a square grid, where the centers of the grid cells are the origins of the rays, which are emitted perpendicularly away from the grid into the simulation box. This configuration is illustrated in Figure 2.

We employ a simple geometric ray tracing algorithm to identify the nearest ray intersection with a satellite surface element. When the $i$th ray intersects the $j$th surface element, it adds an increment

$$ C_{w,i} = c_{w,a,j}C_{w,a,i} + c_{w,d,j}C_{w,d,i} + c_{w,s,j}C_{w,s,i} \quad (20) $$

to the radiation pressure coefficient vector $C_w$. The factors $c_{w,a,j}$, $c_{w,d,j}$, and $c_{w,s,j}$ are the coefficients of absorption, diffuse reflection, and specular reflection, respectively, of the $j$-th surface element. Subscript “w” denotes again the wavelength, i.e. visible light or infrared radiation. Note that these coefficients must sum to one. The coefficient vector increments on the right-hand side are defined as

$$ C_{w,a,i} = A_i r_i \quad (21) $$

for absorption,

$$ C_{w,d,i} = A_i \left( r_i - \frac{2}{3} n_j \right) \quad (22) $$

for diffuse reflection, and

$$ C_{w,s,i} = 2A_i(r_i \cdot n_j)n_j \quad (23) $$

for specular reflection, where $A_i$ is the size of a grid cell, $r_i$ is a unit vector pointing into the propagation direction of the $i$th ray, and $n_j$ is the outward normal of the $j$th surface element. Note that the distance between the radiation source and the satellite is irrelevant in the simulation. We loop over all grid nodes and sum all coefficient vector increments to obtain the radiation pressure coefficient vector:
$$C_w = \sum_i C_{w,i} \qquad (24)$$

The ray tracing simulation accounts not only for self-shadowing but also for multiple reflections. Each ray intersecting the surface spawns new rays if the surface element is reflective. The amount of radiation the initial ray carries is proportional to the cell size $A_i$. Unless the surface is a perfect mirror, a spawned ray with index $k$ will carry a smaller amount of radiation, which we consider by assigning a smaller area $A_k$. The direction of the spawned rays is determined from vectors $r_i$ and $n_j$ and the reflection mode, which can be diffuse and specular. In the case of a specular reflection, the new direction is

$$\mathbf{u}_k = \mathbf{r}_i + 2\mathbf{n}_j \qquad (25)$$

and its new area is

$$A_k = c_{w,s,j} A_i, \qquad (26)$$

where $c_{w,s,j}$ is the same specular reflection coefficient as in equation (20).

For diffuse reflections, we subdivide the half sphere centered around the surface normal into sections of size $\Delta\zeta \times \Delta\lambda$ as illustrated in Figure 3. The angles $\zeta$ and $\lambda$ are the local zenith angle and azimuth, respectively, where the surface normal points to the local zenith and zero azimuth are defined arbitrarily. We spawn a new ray in the center of each section, calculate the fraction of radiation that passes through the section according to Lambert's cosine law, and multiply the fraction by the amount of diffusely reflected radiation:

$$A_k = c_{w,d,j} \frac{\lambda_2 - \lambda_1}{2\pi} (\sin^2 \zeta_2 - \sin^2 \zeta_1) A_i \qquad (27)$$

where $\lambda_1, \lambda_2, \zeta_1$, and $\zeta_2$ define the corners of a section and $c_{w,d,j}$ is the same diffuse reflection coefficient as in equation (20).

### 2.4.2 Modeling the satellite's thermal emission

We need the temperature of the outer surfaces to calculate the satellite's thermal emission. Since no measurements of the surface temperatures are available at the time of writing, we rely on modeling to infer the surface temperatures. For that purpose, we use a thermal satellite model that consists of a set of panels and the inner body of the satellite, which is conceptually very similar to the model of Wöske et al. (2019). In our model, the panels gain heat by absorption of radiation, lose heat by diffusely emitting infrared radiation, and conductively exchange heat with the inner body. The latter will not only exchange heat with the panels but also generate heat $P_{gen}$ because a fraction of the electric power is converted to heat by the payload and other electric parts of the satellite. We select $P_{gen}$ such that the inner body's temperature is about 25 °C, noting that most instruments operate at room temperature.

The heat gain due to the absorption of radiation is

$$P_{\text{abs},j} = \Phi c_{a,j} A_j \cos \theta_j, \qquad (28)$$

where subscript $j$ indicates the panel, $\Phi$ is the radiation flux, $c_{a,j}$ is the coefficient of absorption of visible light or emissivity in the case of infrared radiation, $A_j$ is the area of the panel, and $\theta$ is the angle between the panel normal and the vector from the satellite to the radiation source. Note that we need to sum the heat gain due to the absorption of visible light and infrared radiation from all radiation sources and angles of incidence.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Simulating diffuse reflections by spawning new rays, marked by black arrows. Each ray represents the amount of radiation passing through a section of the half-sphere of size $\Delta\zeta \times \Delta\lambda$ by Lambert's cosine law. The red arrow indicates the surface normal.</figcaption></figure>

The heat loss due to the emission of infrared radiation follows the Stefan–Boltzmann law,

$$P_{\text{emit},j} = A_j \varepsilon_j \sigma T_j^4, \qquad (29)$$

where $P_{\text{emit},j}$ is the emitted power, $\varepsilon_j$ is the emissivity, $\sigma$ is the Stefan–Boltzmann constant, and $T_j$ is the absolute temperature of the panel.

The conductive heat exchange between the panels and the inner body is

$$P_{\text{cond},j} = k_j (T_j - T_{\text{body}}) \qquad (30)$$

where $k_j$ is the thermal conductivity and $T_{\text{body}}$ is the inner body's temperature.

Therefore, the total heat exchange for a panel is

$$P_j = P_{\text{abs},j} - P_{\text{emit},j} - P_{\text{cond},j}. \qquad (31)$$

For the inner body, it is

$$P_{\text{body}} = P_{\text{gen}} + \sum_j P_{\text{cond},j}. \qquad (32)$$

At each time step, we update the panels' and inner body's temperature based on their total heat exchange:

$$T_j(t + \Delta t) = T_j(t) + \frac{P_j}{C_j} \Delta t \qquad (33)$$

and

$$T_{\text{body}}(t + \Delta t) = T_{\text{body}}(t) + \frac{P_{\text{body}}}{C_{\text{body}}} \Delta t \qquad (34)$$

where $C_j$ and $C_{\text{body}}$ are the thermal capacity of the $j$th panel and the inner body, respectively. At the first epoch, we initialize the panels' temperatures with 273 and the body temperature with 298. The initial temperatures no longer influence the thermal model after about six hours. When calculating the temperatures, we use a time step of 1 for convenience because this is the sampling of the accelerometer data.

After calculating the temperatures, we insert them into equation (29) to obtain $P_{\text{emit},j}$ and calculate the radiation pressure due to the diffuse thermal emission by

$$\mathbf{a}_{\text{te}} = -\frac{2}{3} \sum_j \frac{P_{\text{emit},j}}{mc} \mathbf{n}_j \qquad (35)$$
<table><caption>**Table 2.** Characteristics of the orbits, acceleration data, and density and crosswind observations. LTAN is the acronym for the local time of the ascending node.</caption><thead><tr><th></th><th>CHAMP</th><th>GRACE A</th><th>GRACE B</th><th>GRACE C</th></tr></thead><tbody><tr><td>Initial orbit</td><td></td><td></td><td></td><td></td></tr><tr><td>Semi-major axis</td><td>6823 km</td><td>6869 km</td><td>6869 km</td><td>6870 km</td></tr><tr><td>Inclination</td><td>87.3°</td><td>89.0°</td><td>89.0°</td><td>88.9°</td></tr><tr><td>Eccentricity</td><td>0.0029</td><td>0.0021</td><td>0.0021</td><td>0.0011</td></tr><tr><td>LTAN progression</td><td>-5.44 min/d</td><td>-4.47 min/d</td><td>-4.47 min/d</td><td>-4.53 min/d</td></tr><tr><td>Initial altitude</td><td>461 km</td><td>506 km</td><td>506 km</td><td>507 km</td></tr><tr><td>Acceleration data</td><td></td><td></td><td></td><td></td></tr><tr><td>Instrument</td><td>STAR</td><td>SuperSTAR</td><td>SuperSTAR</td><td>SuperSTAR-FO</td></tr><tr><td>Nominal noise level</td><td>1 nm·s<sup>-1</sup></td><td>0.1 nm·s<sup>-1</sup></td><td>0.1 nm·s<sup>-1</sup></td><td>0.1 nm·s<sup>-1</sup></td></tr><tr><td>Data availability start</td><td>2000-07-29</td><td>2002-04-04</td><td>2002-04-04</td><td>2018-05-29</td></tr><tr><td>Data availability end</td><td>2010-09-04</td><td>2017-10-31</td><td>2017-05-22</td><td>-</td></tr><tr><td>Density and crosswind observations</td><td></td><td></td><td></td><td></td></tr><tr><td>Cadence</td><td>10 s</td><td>10 s</td><td>10 s</td><td>10 s</td></tr><tr><td>Resolution along orbit</td><td>152 km</td><td>422 km</td><td>422 km</td><td>152 km</td></tr></tbody></table>

where *c* is the speed of light. The thermal models for the CHAMP, GRACE, and GRACE-FO satellites are provided in the Appendix in Tables A1, A2, and A3, respectively.

# 3 Data

## 3.1 Input data for neutral mass density and crosswind data processing

The input to neutral mass density and crosswind data processing consists of acceleration, attitude, position, velocity, thruster, satellite mass, and accelerometer temperature data. We take these data from the Level 1 data archives of the CHAMP, GRACE, and GRACE-FO missions, all accessible from FTP ftp://isdcftp.gfz-potsdam.de. For the CHAMP, GRACE A, and GRACE C satellites, we used Level 1A data, allowing for cleaner removal of thrust events due to the higher temporal resolution compared to the Level 1B data. Only for the GRACE B satellite, we relied on Level 1B data because Level 1A data was unavailable to us.

The CHAMP acceleration data was measured by the Space Three-axis Accelerometer for Research (STAR) instrument developed by the Office National d'Études et de Recherches Aérospatiales (ONERA). The GRACE and GRACE-FO satellites were equipped with the more advanced SuperSTAR and SuperSTAR-FO accelerometers, respectively, described by Touboul et al. (2012) and Christophe et al. (2015). The nominal noise level of these accelerometers is reported in Table 2, noting that the value applies to the along-track acceleration measurement, whereas the cross-track acceleration measurement is about ten times worse by the design. In the case of the GRACE mission, perturbations stemming from the satellite platform limit the accuracy of the accelerometer data (Flury et al., 2008). The noise level of the acceleration data is about ten times higher than the nominal noise level of the accelerometers (Murböck et al., 2023), i.e., about 1 nm·s⁻², noting that this value depends on the satellite operations. As part of preprocessing the acceleration data, we applied a low-pass filter, whose cut-off frequency was carefully selected to suppress the perturbations. However, this low-pass filtering limits the resolution of the density and

crosswind observations. The resolution along the orbit is the orbital velocity divided by the cut-off frequency, multiplied by two following the Nyquist-Shannon sampling theorem. Table 2 reports relevant characteristics of the initial orbits, acceleration data, and density and crosswind observations.

The neutral mass density and crosswind data processing require modeling of the atmospheric composition and temperature, for which we use the NRLMSISE-00 model described in Section 3.3.1. Wind in the direction of the satellite's x-axis is accounted for by the HWM07 model (Drob et al., 2008).

## 3.2 Input data for accelerometer data calibration

The input to the calibration consists of the same satellite attitude and accelerometer data as mentioned above. However, thrust events are not removed from the accelerometer observations, as they contain an actual signal influencing the satellites' orbits. We use the GRACE and GRACE-FO missions' navigation orbits available in the GRACE Level 1 data as pseudo-observations (Case et al., 2010). For CHAMP, we use the precise orbit solutions determined at the Astronomical Institute of the University of Bern (AIUB) (Prange, 2010). This dataset is, unfortunately, not available for the early mission phase. For the period 2001–2002, therefore, the GeoForschungs Zentrum (GFZ) Rapid Science Orbits (RSO) are used as input (König et al., 2005). All input orbits are determined using a reduced-dynamic approach, which does not use any accelerometer data. We use the calibrated accelerations for thermospheric density retrieval, and therefore, ideally, the input orbits should not contain prior dynamical model information. However, when using only a minimal number of calibration parameters and with differences of only a few cm between both types of orbits, the impact of using reduced-dynamic instead of kinematic orbits on the calibration is small. To assess this, we also used kinematic orbits computed by the Graz University of Technology<fnref n="1" /> (Zehentner & Mayer-Gürr, 2016) as input observations for the CHAMP calibration for 2003–2008. The median RMS of the calibrated acceleration differences in flight direction is

<footnote n="1">¹ We received the orbit data on 9 June 2022 via personal communication with Barbara Suesser-Rechberger, Graz University of Technology.</footnote>
1.4 nm·s⁻², equivalent to about 0.5% of the aerodynamic acceleration and, therefore, negligible.

## 3.3 Thermosphere models

In-situ observations of neutral mass density above 200 km altitude are very sparse. While this emphasizes the importance of the observations presented in this paper, it also highlights the challenge of finding suitable data for validation. Therefore, we choose three thermosphere models to evaluate the new observations presented in this paper. The first is the Naval Research Laboratory Mass Spectrometer Incoherent Scatter radar Expanded model 2000 (NRLMSISE-00), a semi-empirical model with a wide user base. The second is the Drag Temperature Model 2020 (DTM-2020), one of the most recent semi-empirical models. The NRLMSISE-00 and DTM-2020 models output neutral mass density, but no wind. Finally, we use the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM), a physics-based model providing neutral mass density and wind. We provide some more information on these models is in the following sections.

### 3.3.1 NRLMSISE-00

The semi-empirical NRLMSISE-00 model is constructed using atmospheric composition from mass spectrometer data, thermosphere temperature derived from incoherent scatter radar data, and molecular oxygen (O₂) number density derived from solar ultraviolet occultation data amongst others (Picone et al., 2002). Further, it is based on neutral mass density datasets derived from accelerometer measurements and orbital decay data from the 1970s and 1980s. The model drivers are the daily F10.7 solar flux index and its 81-day centered mean F̄₁₀.₇ for solar activity, the 3-hourly ap index, and its daily average, the Ap index, for geomagnetic activity. There are two reasons to include the NRLMSISE-00 model in the evaluation. First, it is largely independent of the neutral mass density observations of the CHAMP, GRACE, and GRACE-FO satellites. The only minor dependency stems from using atmospheric composition and temperature as specified by the NRLMSISE-00 model in the calculation of the aerodynamic coefficient vector (Doornbos, 2011). Second, the model is widely used in science and operations.

### 3.3.2 DTM-2020

We use the operational version of the semi-empirical DTM-2020, driven by the F10.7 solar radio flux and geomagnetic Kp index (Bruinsma & Boniface, 2021). In sharp contrast to the NRLMSISE-00 model, it was constructed using neutral mass density observations derived from the accelerometer data of the GOCE, CHAMP, and GRACE satellites, GNSS tracking data of the Swarm A satellite, and laser ranging data on the Stella satellite. It is worthwhile to note that the neutral mass density observations presented in this paper were not used in the construction of the DTM-2020 model.

### 3.3.3 TIE-GCM

TIE-GCM is a three-dimensional, time-dependent, physics-based model of the Earth's upper atmosphere (Richmond et al.,

1992). The website, www.hao.ucar.edu/modeling/tgcm hosts the open-source TIE-GCM code. Kodikara (2019) and Qian et al. (2014) provide summaries of the recent developments in the model. This study uses TIE-GCM version 2.0 (released on 21 March 2016) with a horizontal resolution of 2.5° × 2.5° in geographic latitude and longitude, and a vertical resolution of 0.25 scale height. The solar irradiance input to the model is specified by the average of the daily solar flux F₁₀.₇, and its 81-day centered mean F̄₁₀.₇. The high latitude mean energy, energy flux, and electric potential are described by the default ion convection model and auroral particle precipitation scheme. The tidal forcing from the lower atmosphere is specified via perturbations from the Hagan et al. (2001) global scale wave model (GSWM). This study also uses the Qian et al. (2009) diurnal eddy diffusion coefficient to add perturbations to the advective and diffusive transport in the TIE-GCM.

For comparison with the satellites, the model neutral mass densities are interpolated in logarithmic space to the satellite location, first horizontally and then vertically to the satellite altitude. To minimize errors due to extrapolation, model estimates are not provided where the satellite altitude is 20 km above the model top. The zonal (east-west) and meridional (north-south) wind components of the model are only interpolated to the satellite altitude. In other words, horizontal winds are not extrapolated vertically upward beyond the model domain. To compare the model winds with the satellite crosswind speeds, the former is projected onto the observed crosswind directions.

# 4 Results and discussion

## 4.1 Accelerometer data calibration

To extract the thermally-induced bias variations, we calculate the residual acceleration

$$ \mathbf{a}_{\text{res}} = \mathbf{a}_{\text{meas}} - \mathbf{S}_{\text{aero,mod}}\mathbf{a}_{\text{aero,mod}} - \mathbf{a}_{\text{rp}} - \mathbf{b}_{\text{trend}} \quad (36) $$

by subtracting the scaled aerodynamic acceleration aₐero,mod based on the NRLMSISE-00 model, where Sₐero,mod is a diagonal matrix, the radiation pressure acceleration aᵣp from equation (3), and a linear trend bᵣtrend. Scaling of the modeled aerodynamic acceleration is necessary to compensate for biases in the thermosphere model, and removing the trend reduces the part of the bias unaffected by the temperature. The scaling matrix Sₐero,mod and the trend bᵣtrend are fitted to aᵣmeas in the intervals from 2007-01-01 00:00 UTC to 2007-01-14 21:20 UTC and 2007-01-24 03:33 UTC to 2007-02-01 00:00 UTC to fitting to the thermally-induced bias variations in-between.

Then, we fit the model for the thermally-induced bias variations in equation (9) to aᵣres to obtain the temperature sensitivity coefficients S_U and S_V and the temperature delay parameters k_U and k_V. The co-estimation of a linear trend ensures that the other parameters fit the temperature swing instead of the mean temperature. We illustrate the results in Figure 4, in which the left and center panels show the fit to aᵣres,x and aᵣres,y and the right panel illustrates the measured and modeled temperatures, T and U, respectively.

In the center panel, we see that the acceleration in the y-direction is most sensitive to temperature variations, showing a peak of −800 nm·s⁻² in response to a temperature swing of about 4.5 °C. Further, the measured temperature (red curve)
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Modeling GRACE A accelerometer bias variations for the along-track axis (left) and cross-track axis (center). Correction 1 uses the measured temperature, i.e. $b_T(t) = s_U T(t)$, and correction 2 the modeled temperature, i.e. $b_T(t) = s_U T(t)$. The RMS of the fit within the time window from 2007-01-17 00:00 UTC to 2007-01-21 00:00 UTC is reported in the brackets. The measured and modeled temperatures are shown in the right panel.</figcaption></figure>

<table><caption>**Table 3.** Model for thermally-induced bias variations for the GRACE accelerometers.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>GRACE A full mission</th><th>GRACE B until 2011-04-27</th><th>GRACE B 2011-04-27</th></tr></thead><tbody><tr><td>k<sub>U</sub></td><td>K<sup>-3</sup>·s<sup>-1</sup></td><td>3.1 × 10<sup>-13</sup></td><td>2.3 × 10<sup>-13</sup></td><td>2.3 × 10<sup>-13</sup></td></tr><tr><td>s<sub>U,x</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>-7.665</td><td>0.0</td><td>0.0</td></tr><tr><td>s<sub>U,y</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>177.7</td><td>69.3</td><td>69.3</td></tr><tr><td>s<sub>U,z</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>-0.2</td><td>-0.8</td><td>-0.8</td></tr><tr><td>k<sub>V</sub></td><td>K<sup>-3</sup>·s<sup>-1</sup></td><td>0.0</td><td>7.0 × 10<sup>-13</sup></td><td>7.0 × 10<sup>-13</sup></td></tr><tr><td>s<sub>V,x</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>0.0</td><td>-2.28</td><td>-3.8</td></tr><tr><td>s<sub>V,y</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>0.0</td><td>65.0</td><td>65.0</td></tr><tr><td>s<sub>V,z</sub></td><td>nm·s<sup>-2</sup>·K<sup>-1</sup></td><td>0.0</td><td>-0.3</td><td>-0.3</td></tr></tbody></table>

<table><caption>**Table 4.** Scale factors for the CHAMP, GRACE, and GRACE-FO accelerometer data.</caption><thead><tr><th>Scale factor</th><th>CHAMP full mission</th><th>GRACE A full mission</th><th>GRACE B until 2011-04-27</th><th>GRACE B after 2011-04-27</th><th>GRACE C full mission</th></tr></thead><tbody><tr><td>S<sub>x</sub></td><td>0.83 ± 2 × 10<sup>-5</sup></td><td>0.94 ± 2 × 10<sup>-5</sup></td><td>0.95 ± 2 × 10<sup>-5</sup></td><td>0.93 ± 3 × 10<sup>-5</sup></td><td>1.05 ± 2 × 10<sup>-4</sup></td></tr><tr><td>S<sub>y</sub></td><td>0.85 ± 4 × 10<sup>-4</sup></td><td>0.96 ± 4 × 10<sup>-4</sup></td><td>0.97 ± 6 × 10<sup>-4</sup></td><td>0.96 ± 6 × 10<sup>-4</sup></td><td>0.94 ± 9 × 10<sup>-4</sup></td></tr><tr><td>S<sub>y</sub></td><td>—</td><td>0.93 ± 1 × 10<sup>-4</sup></td><td>0.95 ± 5 × 10<sup>-4</sup></td><td>0.95 ± 9 × 10<sup>-4</sup></td><td>0.98 ± 7 × 10<sup>-4</sup></td></tr></tbody></table>

does not fit the residual acceleration (blue curve, which is behind the yellow curve). In contrast, the model temperature (yellow curve) fits almost perfectly, demonstrating the need to model the temperature delay as described in Section 2.1.2. The acceleration in the x-direction shows a response of 30 nm·s⁻² to the same temperature swing. Even though the acceleration in the x-direction is much less sensitive to temperature variations, we can still see that the model temperature fits better than the measured temperature. The acceleration in the z-direction (not shown) has an almost negligible temperature sensitivity of a few tenths of nm·s⁻².

We analyzed several similar temperature swings in the same way, and the results combined to obtain the values reported in Table 3. We modeled the thermally-induced bias variations with two temperatures for GRACE B to get a better fit. The need for two temperatures suggests that the heat is transferred through two paths. However, when summing the coefficients for GRACE B, e.g., $S_{Uy} + S_{Vy} = 69.3 \text{ nm} \cdot \text{s}^{-2} \cdot \text{K}^{-1} + 65.0 \text{ nm} \cdot \text{s}^{-2} \cdot \text{K}^{-1} = 134.3 \text{ nm} \cdot \text{s}^{-2} \cdot \text{K}^{-1}$, we obtain qualitatively similar values as for GRACE A, which is not surprising as the satellites and the instruments are identical.

Finally, we point out that the model for the thermally-induced bias variations acts as a lowpass filter on the temperature data. Even after deactivating the accelerometer's thermal control on the GRACE satellites in April 2011, the modeled bias variations at the orbital period have a tiny amplitude of only a few tens of pm·s⁻², which is negligible.

Table 4 summarizes the scale factors and their formal standard deviations determined with the stacked normal matrices approach for all satellites. Since the unconstrained estimation of daily scale factors for GRACE A has a stable behavior over time (not shown), we applied one set of constant scale factors for the entire lifetime. For GRACE B, we observed a small drop in the time series of daily scale factors after the accelerometer's thermal control deactivation (not shown). Therefore, we estimated two sets of constant scale factors for this satellite.

The scale factors' formal standard deviations are all smaller than $10^{-3}$, indicating a high estimation precision. Considering that the scale factors and the aerodynamic coefficient vectors in the x-direction both have a scaling effect on the density observations, we conclude that the scale factors' formal standard deviations are negligible compared to the uncertainty of at least a few percent in the aerodynamic coefficient vector
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** GRACE A accelerometer biases (top three panels) at different stages of the calibration. The bottom panel shows the RMS of the three-dimensional position differences resulting from the POD.</figcaption></figure>

(Mehta et al., 2022). Thence, we rounded the scale factors to two digits. For GRACE A, the standard deviations of the scale factors are $2 \times 10^{-5}$, $4 \times 10^{-4}$, and $1 \times 10^{-4}$ for $S_x$, $S_y$, and $S_z$, respectively, which indicates that especially in the x-direction, the scale factors are well-determined with the POD-based calibration. We obtained similar standard deviations for GRACE B and CHAMP. For GRACE C, the uncertainty in the scale factor calibration is substantially larger, with standard deviations of $2 \times 10^{-4}$, $9 \times 10^{-4}$, and $7 \times 10^{-4}$ for $S_x$, $S_y$, and $S_z$, respectively. The worse results are due to the shorter dataset and the small non-gravitational acceleration signal encountered at about 520 km altitude, the highest of all satellites considered in this paper (cf. Fig. 6, bottom panel). We expect to improve these results when more data at lower altitudes and higher solar activity becomes available.

After determining the scale factors, the next step in the accelerometer data calibration is the estimation of the biases. We illustrate the results of the bias estimation for GRACE A in Figure 5. The blue curve shows the unconstrained estimates of the daily biases without modeling the thermally-induced bias variations. As expected, the accelerometer biases have large offsets on the order of 1 for the x- and z-directions and 30 in the y-direction (Touboul et al., 2016). On top of the offsets, the biases show an exponential behavior until mid-2010, likely due to instrument aging. In mid-2010, we observed a step caused by a 5 °C drop in the sensor unit's temperature due to a change in the accelerometer thermal control. In April 2011, the accelerometer's thermal control was entirely deactivated, leading to temperature variations up to 10 °C, which has a significant impact on the bias estimates (cf. also Fig. 4).

The red curves show the unconstrained estimates of the daily biases when modeling the thermally-induced bias variations, noting that the red curves are covered by the yellow/purple curves when they are not visible. The large offsets to the unconstrained estimates without modeling the thermally-induced bias variations (blue curves) are due to using absolute temperatures in equation (9). The offsets are compensated by the daily biases and, thus, are of no concern. More importantly, when modeling the thermally-induced bias variations, the biases are much less affected by temperature variations, except for the wiggle in 2013–2014, where the GRACE A accelerometer temperature was unavailable. After experimenting, we found that using the GRACE B temperature was the best solution. The RMS in the bottom panel is improved, particularly from April 2011 onward. This reduction indicates that the calibration parameters are more accurate when modeling the thermally-induced bias variations (red dots are behind the yellow dots).

The yellow curve results from applying a constraint on the z-bias, suppressing artificial bias variations of about $200 \text{ nm} \cdot \text{s}^{-2}$ in that component. This suppression reduces the number of outliers in the x- and z-components and removes a slight artificial variation in the x-bias with an amplitude of $1 \text{ nm} \cdot \text{s}^{-2}$ as shown in the zoom box in the top panel.

The purple line is only relevant for the y-bias, which we adjusted as described in Section 2.1.3. The most significant adjustments have a size of $100 \text{ nm} \cdot \text{s}^{-2}$, which we highlight in the zoom box in the second panel from the top. We can also see that the purple line is much smoother than the yellow, demonstrating improved precision. After data gaps, the purple line shows a spike. The spikes are due to the imperfect
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Size of the radiation pressure acceleration relative to the aerodynamic acceleration, where the size is measured by the RMS of the acceleration within a sliding one-month window. The curves show the ratio of the RMS of the radiation pressure acceleration over that of the aerodynamic acceleration in the along-track (top) and cross-track direction (center). The bottom panel shows the altitude evolution of the satellites and the solar activity as indicated by the $F_{10.7}$ solar radio flux index.</figcaption></figure>

modeling of the increasing accelerometer sensor temperature after the instrument is switched on, noting that the y-axis is most sensitive to temperature (cf. Table 3). In that case, the heat source is the instrument's electronics, whereas the heat source is external to the instrument otherwise, which explains the poor modeling of that feature. We provide the crosswind observations with a flag that marks the affected periods as invalid.

We obtained very similar results for the GRACE B satellite. We also modeled the thermally-induced bias variations for the CHAMP satellite, though they had a much smaller impact than for the GRACE A satellite. In the case of the GRACE C satellite, the accelerometer's thermal control improved a lot, eliminating the need to model thermally-induced bias variations. Further, we did not apply z-bias constraints for the GRACE C satellite because they did not significantly affect the already smooth x-bias. The appendix provides figures of the estimated biases of the CHAMP, GRACE B, and GRACE C satellites.

## 4.2 Radiation pressure modeling

The radiation pressure acceleration is significant compared to the aerodynamic acceleration at altitudes above 450 km during low solar activity. Figure 6 compares the size of the radiation pressure acceleration to the size of the aerodynamic acceleration, where the RMS calculated within a one-month time window indicates the size. For the along-track direction, we see that the size of the radiation pressure acceleration is

about 100% of the aerodynamic acceleration in 2007–2009 for GRACE A and even 200% in 2018–2020 for GRACE C. Suppose that radiation pressure modeling has an uncertainty of 5%, which is a realistic assumption. Since the neutral mass density observations are approximately proportional to the along-track acceleration, the resulting uncertainty in the neutral mass density observations will be 5% for GRACE A and 10% for GRACE C. The situation is much more difficult for the cross-track accelerations, where the radiation pressure acceleration is 20 times larger than the aerodynamic acceleration for the GRACE A and GRACE C satellites. In this context, we point out that lift causes the most significant fraction of the cross-track aerodynamic acceleration, and only a small fraction is due to crosswind. Thus, accurately modeling the radiation pressure is crucial to obtain usable crosswind observations for the GRACE A and GRACE C satellites. The radiation pressure acceleration modeling is less critical for the CHAMP satellite because the aerodynamic acceleration is larger throughout the mission lifetime.

The CHAMP, GRACE, and GRACE-FO density and crosswind observations presented in this paper are available on the FTP ftp://thermosphere.tudelft.nl, designated as version 2 (V2). On the same FTP, we also provide the previous version presented by March et al. (2021), designated as version 1 (V1). We implemented two significant upgrades in V2 to improve the radiation pressure modeling compared to V1. First, we switched from commonly used panel models to augmented
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** GRACE A radiation pressure acceleration in the along-track (top) and cross-track (bottom) directions. Version 1 is the radiation pressure as applied to the previous neutral mass density dataset by March et al. (2021), i.e., based on a panel model and excluding the effect of the thermal emission. Version 2 is the radiation pressure as presented in this paper.</figcaption></figure>

<table><caption>**Table 5.** RMS of radiation pressure acceleration for GRACE A in 2008 in units of nm·s⁻².</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">Along-track</th><th colspan="3">Cross-track</th></tr><tr><th>Version 1</th><th>Version 2</th><th>Difference</th><th>Version 1</th><th>Version 2</th><th>Difference</th></tr></thead><tbody><tr><td><b>a</b><sub>rp</sub></td><td>15.5</td><td>13.9</td><td>1.9</td><td>23.1</td><td>27.3</td><td>4.4</td></tr><tr><td><b>a</b><sub>srp</sub></td><td>15.4</td><td>13.1</td><td>2.4</td><td>22.9</td><td>20.8</td><td>2.3</td></tr><tr><td><b>a</b><sub>alb</sub></td><td>0.2</td><td>0.2</td><td>0.0</td><td>0.2</td><td>0.2</td><td>0.0</td></tr><tr><td><b>a</b><sub>ir</sub></td><td>0.2</td><td>0.2</td><td>0.0</td><td>0.1</td><td>0.1</td><td>0.0</td></tr><tr><td><b>a</b><sub>te</sub></td><td>—</td><td>0.8</td><td>0.8</td><td>—</td><td>6.4</td><td>6.4</td></tr></tbody></table>

detailed geometry models described in Section 2.2 and fine-tuned the absorption and reflection coefficients of the solar array ($c_a = 0.65 \rightarrow 0.72$, $c_d = 0.30 \rightarrow 0.05$, $c_s = 0.05 \rightarrow 0.23$) and the nadir panel ($c_a = 0.12 \rightarrow 0.12$, $c_d = 0.20 \rightarrow 0.06$, $c_s = 0.68 \rightarrow 0.82$). We verified these changes by inspecting photos of the satellites, which confirmed that the surfaces were indeed more specular reflecting. Second, we model the radiation pressure due to the satellite thermal emission as described in Section 2.4.2, which is a new feature of the V2 dataset. We will elaborate on these changes and their effects in another paper and only present the total effect of the changes here.

Figure 7 shows the GRACE A radiation pressure acceleration in the along-track and cross-track directions as modeled in V1 and V2. The magnitude of the radiation pressure acceleration is about 30 nm·s⁻² in the along-track direction and 50 nm·s⁻² in the cross-track direction. The periodic behavior results from the beta angle progression, i.e., the angle between the orbital plane and the sun's direction, which is 360° in about 11 months. The differences between V1 and V2 amount to several nm·s⁻² in both directions, which is particularly significant for the cross-track direction. Table 5 reports the RMS of the radiation pressure acceleration for the contributions from solar radiation pressure, Earth's albedo, Earth's infrared radiation, and the satellite's thermal emission. The difference

between V1 and V2 radiation pressure acceleration $a_{rp}$ is 1.9 nm·s⁻² and 4.4 nm·s⁻² in the along-track and cross-track directions, respectively. These differences are 6% and 9% of the radiation pressure acceleration in along-track and cross-track directions, respectively. Given the relative size of the aerodynamic acceleration illustrated in Figure 6, we may expect differences of 6% in the neutral mass density observations and more than 200% in crosswind observations for GRACE during 2007–2010. Further, Table 5 demonstrates that the changes in solar radiation pressure and the newly implemented thermal emission cause the differences. This insight narrowed down the search space in fine-tuning the radiation pressure model.

### 4.3 Density observations

Neutral mass density observations $\rho_1$ are best compared to other observed or model densities $\rho_2$ in log space (Sutton, 2018). The mean density ratio

$$ \mu(\rho_1, \rho_2) = \exp \left( \frac{1}{N} \sum_{n=1}^{N} \ln \frac{\rho_{1,n}}{\rho_{2,n}} \right) \quad (37) $$

quantifies deviations in the scale of the density observations, which could be related to incorrect scale factors in the
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the old (V1) and new (V2) neutral mass density observations of GRACE A in terms of orbit mean µ(ρ<sub>V1</sub>, ρ<sub>V2</sub>) (top) and orbit standard deviation δσ(ρ<sub>1</sub>, ρ<sub>2</sub>) (middle). Peak values exceeding the limits of the vertical axis are as low as 0.4 and as large as 2.0 in the orbit mean and as large as 200% in the orbit standard deviation. The accelerometer temperature is illustrated for reference (bottom).</figcaption></figure>

accelerometer data calibration, mismodeling of the aerodynamic and radiation pressure coefficient vectors, or a biases in the density observations. The standard deviation of the density ratio

$$
\sigma(\rho_1, \rho_2) = \exp\left(\sqrt{\frac{1}{N-1} \sum_{n=1}^{N} \left(\ln\frac{\rho_{1,n}}{\rho_{2,n}} - \ln\mu(\rho_1, \rho_2)\right)^2}\right) \quad (38)
$$

measures if the variability in the densities is the same. It is
more convenient to express the deviation of σ(ρ₁, ρ₂) from
unity as a percentage

$$
\delta\sigma(\rho_1, \rho_2) = (\sigma(\rho_1, \rho_2) - 1) \times 100. \quad (39)
$$

Figure 8 compares the old (V1) and new (V2) neutral mass density observations of GRACE A, where we calculated the mean µ(ρ<sub>V1</sub>, ρ<sub>V2</sub>) and the standard deviation δσ(ρ<sub>1</sub>, ρ<sub>2</sub>) in a sliding window of the length of one orbital revolution. The mean is very close to 1 in the mission's first years but begins to deviate from 1 from 2006 onward, reaching values as low as 0.8 in 2009. The change in the radiation pressure modeling has a more considerable impact when the size of the radiation pressure acceleration approaches the size of the aerodynamic acceleration (cf. Fig. 6) and, thus, causes the deviation. The standard deviation also reflects this, which is about 10% in 2002–2003 and then gradually increases until 2009 when it reaches values up to 100%. Another cause of the large standard deviation in 2006–2010 is the improved accelerometer data preprocessing that reduces acceleration spikes due to attitude control thruster activations and other artifacts more effectively, which was not possible when producing the V1 dataset.

Further, several large spikes in the mean are associated with spikes in the temperature, which may occur, e.g., when the accelerometer is power-cycled. Figure 9 shows the V1 and V2 neutral mass density observations in April 2008, when one such temperature spike occurred (black line). The V1 density observations (blue line) show an artificial variation, even

dropping to negative values, which is not the case for the V2
density observations (red line). This comparison demonstrates
that modeling the temperature-induced bias variations (cf.
Sect. 2.1.2) improved the V2 dataset substantially.

We compare the neutral mass density observations to mod-
els in terms of the annual mean observation-model ratio in
Figure 10 and annual standard deviation in Figure 11. We will
not discuss the difference between the observations and the
models, for which we refer to Bruinsma et al. (2018). Instead,
we will focus on the difference between new (V1) and old
(V2) observations and how the previously unavailable GRACE
A and GRACE C observations (beyond 2011) compare to the
existing datasets (until 2010).

First, we notice that the mean of the new and old CHAMP
observations is practically identical. This agreement is not sur-
prising since the newly calibrated accelerometer data is almost
identical to the old version, the radiation pressure acceleration
is small compared to the aerodynamic acceleration (cf. Fig. 6)
and tends to average out to zero over long periods, and the aero-
dynamic model is the same for both versions. The standard
deviation of the new observations in 2005–2009 is about
0.5% smaller than that of the old version. For the
NRLMSISE-00 model, we notice that the standard deviation
of the new observations in 2001–2004 is higher than that of
the old version. However, this feature is not evident in the
comparison to the DTM-2020 and TIE-GCM models.

We find the best insights for the GRACE A satellite, for
which we have old and new observations for 2002–2009. We
see a decrease of the mean density ratio from 2006 to 2009
of about 10% for the DTM-2020 model and more than 20%
for the NRLMSISE-00 and TIE-GCM models (Fig. 10, middle
panel). We attribute this mainly to the models overestimating
density as the solar activity decreased towards the deep solar
minimum in 2008/2009. However, the decrease of the mean
density ratio is a few percent smaller for the new observations,
resulting in a more constant time series. This feature is most
evident for the DTM-2020 model, where the decrease of the
mean density ratio from 2005 to 2009 is about 6% smaller
<figure><img src="image_10.png" /><figcaption>Fig. 9. Old (V1) and new (V2) neutral mass density observations of GRACE A (top) and accelerometer temperature (bottom).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Observation-model mean ($\mu$) for the CHAMP (left), GRACE A (center), and GRACE C (right) satellites. GRACE B statistics are omitted because they are practically identical to GRACE A.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Observation-model standard deviation ($\delta\sigma$) for the CHAMP (left), GRACE A (center), and GRACE C (right) satellites. GRACE B statistics are omitted because they are practically identical to GRACE A.</figcaption></figure>

for the new density observations than the old version. We interpret this as a clear improvement in the radiation pressure modeling. From 2011 onward, the mean density ratio of the new observations and the DTM-2020 model varies between 1.0 and 1.2, which agrees well with the early years of the GRACE A satellite. The NRLMSISE-00 model shows a similar agreement, though the mean density ratio is lower, varying from 0.7 to 0.9 in 2011–2018 and from 0.8 to 0.85 in 2002–2005.

The standard deviation of the GRACE A observation-model density ratios (Fig. 11, middle panel) provides the strongest evidence of the improvements in the new observations. The standard deviation of the old observations increased from

2002 to 2009 by several tens of percent for all models, whereas the increase is 10–15% smaller for the new observations. In the case of the DTM-2020 model, the peak standard deviation in 2009 is 30%, which is only marginally larger than the standard deviation in 2004–2007, which is about 25%. Beyond 2010, the standard deviation for the DTM-2020 decreases to the lowest value of about 18%. Even though the standard deviation is 5% larger for the TIE-GCM model, the time series show the same variations. The standard deviation resembles the radiation pressure and aerodynamic acceleration ratio presented in Figure 6 (top panel, blue curve). This resemblance confirms that the lower standard deviation of the new observations is due to
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Crosswind speed observed by the GRACE A satellite (top), wind speed from the TIE-GCM model in the same direction (center), and the difference between the observed and model wind speed (bottom). The observed wind speed in the top panel is positive by definition. In contrast, negative values in the center panel indicate that the model wind flows in the opposite direction than the observed wind.</figcaption></figure>

the advanced radiation pressure modeling. However, reducing artifacts in the new preprocessing likely contributed to the improvement.

Since the GRACE C density observations are new, we compare them to the GRACE A observations focusing on the continuity from GRACE A to GRACE C, noting that the observation gap is only seven months. In the case of the DTM-2020 and the NRLMSISE-00 models, the mean density ratio connects well from GRACE A to GRACE C: The DTM-2020 model shows nearly identical values in 2017 and 2018, and the decrease in the NRLMSISE-00 mean density ratio from 2016 to 2017 continues to 2018. This agreement indicates that there are no significant offsets between GRACE A and GRACE C in the mean density ratio. We show the mean density ratio for the TIE-GCM model for completeness, noting that most GRACE C observations are above the upper model boundary and, therefore, are difficult to interpret.

We discussed above that the standard deviation is sensitive to the radiation pressure and aerodynamic acceleration ratio. There is a large discontinuity in the ratio of GRACE A and GRACE C because of the vast different altitudes near the re-entry of GRACE A and the launch of GRACE C

(cf. Fig. 6, bottom panel). Therefore, we compare the GRACE C standard deviation to that of GRACE A from 2008–2009 since the radiation pressure and aerodynamic acceleration ratios are closest then (cf. Fig. 6, top panel) due to comparable altitude and solar activity. The GRACE C standard deviation is approximately 25–31% in 2018–2021 (Fig. 11, right panel). These values are consistent with the GRACE A standard deviation for DTM-2020, which was 30% in 2008–2009. We observe a similar consistency for the NRLMSISE-00 and TIE-GCM models. In this context, we point out that the GRACE C accelerometer data contains substantially fewer artifacts and has a much lower measurement noise level than the GRACE A accelerometer data. Thus, we can only explain the consistency of the GRACE A and GRACE C standard deviation by the accuracy of the radiation pressure model. Therefore, we attribute the improvement observed for GRACE A mostly to the advanced radiation pressure model.

## 4.4 Crosswind observations

When comparing the crosswind observations with the TIE-GCM model, we project the model wind vector $v_{\text{mod}}$ onto
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Statistical comparison of the crosswind speed observed by the GRACE A satellite and wind speed from the TIE-GCM model in the same direction. The top, middle, and bottom panels show the mean, standard deviation, and Pearson's correlation coefficient calculated in a sliding one-day window.</figcaption></figure>

the direction of the observed crosswind $\mathbf{u}_{\text{obs}}$ (unit vector) to obtain the model wind

$$v_{\text{mod}} = \mathbf{v}_{\text{mod}} \cdot \mathbf{u}_{\text{obs}} \quad (40)$$

in the direction of the observed wind.

Figure 12 shows the GRACE A observed crosswind speed in July 2013 in the top panel, the TIE-GCM model wind in the direction of the observed wind in the center panel, and the difference between the observed and model wind speed in the bottom panel. The data is organized by time on the horizontal axis and argument of latitude on the vertical axis. The latter is the angle within the orbital plane, measured from the ascending equator crossing (Montenbruck & Gill, 2012). Each “column” in this plot represents one orbital revolution, starting from the ascending equator crossing at 0° argument of latitude, through the northernmost point in the orbit at 90° argument of latitude, the descending equator crossing at 180° argument of latitude, the southernmost point in the orbit at 270° argument of latitude, until the next ascending equator crossing at 360° = 0° argument of latitude.

We observe high wind speeds up to 800 m·s⁻¹ near the northern (90°) and southernmost (270°) points in orbit, which are associated with the auroral oval (Lühr et al., 2007). Outside of the polar regions, the observed wind speed ranges from 0 m·s⁻¹ to 200 m·s⁻¹. The almost horizontal stripe of fast wind speeds, starting at 180° argument of latitude on 1 July 2013 and lasting until 10 July, is an artifact of the eclipse entry, stemming inaccurate modeling of Earth's shadow. The pattern of the TIE-GCM model wind is very similar to the observed wind. In the polar regions, we see primarily positive model wind speeds, which means that the direction of the model and observed crosswind is the same, although the magnitude of the model wind is smaller than the observed one. Near the equator crossings (0° and 180°), the model wind is much smaller than the observed wind, which indicates that the inherent noise in the crosswind observations is about 200 m·s⁻¹. However, we will

demonstrate in the following that the noise in the crosswind observations is far from constant.

Figures 13 and 14 provide long-term comparisons of the GRACE A and CHAMP observed crosswinds to the TIE-GCM model wind. The top panels show the mean crosswind speed, the center panels show the standard deviation, and the bottom panels show Pearson's correlation coefficient, all calculated within a sliding window of one day. The GRACE A observed crosswind mean (Fig. 13, top panel, blue curve) shows a large variability during the mission lifetime. In 2002–2003 and 2012–2017, the observed crosswind mean agrees well with the model one, although it is slightly larger by a few tens of m·s⁻¹ as visible in the zoom on 2015–2016. This bias in the observed crosswind mean reflects the noise in the crosswind observations, which cannot assume negative values because the crosswind observations are positive by definition. The observed crosswind mean in 2004–2011 is unrealistically large compared to model one. The cause is remaining errors in the accelerometer data calibration, radiation pressure modeling errors, and cross-track acceleration measurement noise, which play a more considerable role when the radiation pressure to aerodynamic acceleration ratio is significant. Inspecting the ratio illustrated in Figure 6 (top panel, blue curve), we find that a ratio lower than 0.2 is necessary to derive meaningful GRACE A crosswind observations. Finally, the mean's magnitude follows a cyclic pattern of about 5.5 months, most apparent in 2012–2018. This pattern is caused by the beta angle progression of roughly 360° in 11 months, noting that the thermospheric wind is strongly organized by local time (Lühr et al., 2007).

The standard deviation of the GRACE A crosswind observations shows the same features as the mean, namely a good agreement with the TIE-GCM model in 2002–2003 and 2012–2018, unrealistically large values in 2004–2011, significant spikes after data gaps, and the cyclic pattern depending on the beta angle. The only difference to the mean is that the
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Statistical comparison of the crosswind speed observed by the CHAMP satellite and wind speed from the TIE-GCM model in the same direction. The top, middle, and bottom panels show the mean, standard deviation, and Pearson's correlation coefficient calculated in a sliding one-day window.</figcaption></figure>

standard deviations of the observations and the model have the same size.

The correlation between the GRACE A and TIE-GCM crosswinds varies between -0.2 and 0.8. When inspecting the zoom box, we notice that the correlation is low when the mean and the standard deviation are small. When that happens, it is impossible to reliably quantify the correlation, meaning we should not conclude those values. With that in mind, we state that the correlation between observed and model crosswinds is about 0.5 in 2002 and continuously decreases to 0.0 in 2007, indicating that observed crosswinds are essentially random. The decreasing aerodynamic acceleration signal ultimately caused the decrease in correlation. The zero correlation persisted until 2009 when the solar activity reached a minimum (cf. Fig. 6, bottom panel). In 2010–2011, the correlation increased to 0.3, and from 2012 onward to 0.6.

In summary, the GRACE A crosswind observations show a good agreement with the TIE-GCM model in 2002–2003 and 2012–2018. This agreement is a remarkable result considering the processing complexity required to mitigate the adverse effects of the deactivated accelerometer thermal control, noting that the cross-track acceleration measurements are most sensitive to temperature variations.

The CHAMP crosswind observations, compared to the TIE-GCM in Figure 14, behave more consistently than the GRACE A crosswind observations. For CHAMP, the mean, standard deviation, and correlation also show a cyclic behavior with a period of 4.5 months related to the beta angle progression of roughly 360° in 9 months. The mean of the CHAMP crosswind observations follows the same pattern as the TIE-GCM winds. However, the observed crosswind mean is about 100 m·s⁻¹, whereas the model crosswind is only about 50 m·s⁻¹ on average. The observed and model crosswinds' standard deviations are approximately 100 m·s⁻¹ and agree reasonably well. The correlation is 0.5–0.6 throughout the CHAMP mission lifetime. These consistently good results are a consequence

of the consistently large signal-to-noise ratio in the CHAMP aerodynamic accelerations.

## 5 Conclusions and outlook

We aim to provide the scientific community with the most complete, consistent, and accurate database of thermosphere density and crosswind observations derived from accelerometer measurements. The results presented here have filled the remaining gaps in our database, which is an important milestone in achieving this goal. In particular, we have produced the first GRACE-FO neutral mass density dataset extended the GRACE neutral mass density observations from 2009 until decommissioning in 2017 and generated a new GRACE crosswind dataset. The latter is published only for 2002–2003 and October 2011 to October 2017 due to the signal-to-noise ratio.

When producing the results of this study, we used a gas-surface interaction model with a constant value for the energy accommodation as described by March et al. (2021). This disputable choice affects the mean of the density observations as described by Mehta et al. (2022), with the largest effects at the highest altitude of 500 km during low solar activity. Incidentally, radiation pressure modeling is essential for deriving neutral mass density observations under the same conditions. Recent investigations by Bernstein et al. (2020) and Bernstein & Pilinski (2022) demonstrate that we need more accurate knowledge of the gas-surface interaction. The results of this study open up new possibilities for investigating gas-surface interaction modeling, which will benefit from the advances in radiation pressure modeling. After thoroughly investigating the model parameters and the resulting consistency between missions, we plan to employ a variable gas-surface interaction model. This upgrade will further improve the accuracy of the observations, especially those of the GRACE and GRACE-FO satellites, which have
been orbiting near 500 km altitude for a significant fraction of the mission lifetime.

The advances in radiation pressure modeling may also prove helpful for gravity field modeling. The accelerometer on the GRACE B satellite was switched off at the end of the mission because of battery issues, necessitating the use of the GRACE A accelerometer measurements for both satellites (Bandikova et al., 2019). Since the satellites flew in nearly identical orbits separated by approximately 220 km and in different attitudes, they experienced slightly different non-gravitational accelerations. Applying several corrections to the GRACE A accelerometer measurements to account for that difference was referred to as “transplanting” the GRACE A accelerometer measurements to the GRACE B satellite. In the case of the GRACE-FO mission, the GRACE D satellite suffered from an anomaly shortly after launch, after which the accelerometer performed much worse than expected (McCullough et al., 2019), requiring the same “transplantation” approach. Advanced radiation pressure modeling could increase the transplant’s accuracy. Similarly, it would be interesting to test whether the accelerometer calibration and the model for the thermally-induced bias variations can improve the retrieval of gravity field models in the case of the GRACE mission.

Although we have made considerable progress in radiation pressure modeling, there is still room for improvement. For example, we currently neglect the conversion of absorbed radiation into electricity by solar arrays. While implementing this process is straightforward, it requires non-trivial fine-tuning of the absorption and reflection coefficients, as well as the thermal capacitance and conductivity parameters, which we have optimized without considering this effect. This fine-tuning is underway and will be implemented in the next release of the neutral mass density and crosswind observations.

Finally, we plan to perform a comprehensive uncertainty quantification that takes into account all relevant sources of error and noise: accelerometer measurement noise, calibration errors in the accelerometer data, errors in the radiation pressure and aerodynamic models, including the gas-surface interaction model and thermosphere models providing atmospheric composition, temperature, and in-track wind (Bruinsma et al., 2022). The uncertainty quantification will facilitate the assimilation of neutral mass density observations into thermospheric models.

**Acknowledgements.** The Thermosphere Observations from Low Earth Orbiting Satellites (TOLEOS) project was funded by ESA via the Swarm Data, Innovation and Science Cluster (DISC), subcontract number SW-CO-DTU-GS-129. The GRACE Level 1A accelerometer data was kindly provided by the Center for Space Research, The University of Texas at Austin, US. The CHAMP Level 1 data was kindly provided by GFZ, Germany. The GRACE-FO geometry model was derived from a CAD model of the satellite, which was kindly provided by Airbus Defence and Space. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

## Data availability

This paper’s new density and crosswind observations are available in ASCII table format on the FTP server ftp://thermosphere.tudelft.nl, designated as version 2 (V2). On

the same server, we also publish the satellite geometry models, the aerodynamic coefficient vectors $C_{\text{aero}}$, and the radiation pressure coefficient vectors $C_w$. The new density and crosswind observations are also available in Common Data Format (CDF) on the ESA Swarm HTTP server http://swarm-diss.eo.esa.int and FTP server ftp://swarm-diss.eo.esa.int in the multi-mission directory.

## References

Bandikova T, McCullough C, Kruizinga GL, Save H, Christophe B. 2019. GRACE accelerometer data transplant. *Adv Space Res* **64(3)**: 623–644. https://doi.org/10.1016/j.asr.2019.05.021.

Bernstein V, Pilinski M. 2022. Drag coefficient constraints for space weather observations in the upper thermosphere. *Space Weather* **20**: e2021SW002977. https://doi.org/10.1029/2021SW002977.

Bernstein V, Pilinski M, Knipp D. 2020. Evidence for drag coefficient modeling errors near and above the oxygen-to-helium transition. *J Spacecr Rockets* **57(6)**: 1246–1263. https://doi.org/10.2514/1.A34740.

Bettadpur S. 2012. Gravity recovery and climate experiment: Product specification document. *Technical report GRACE 327-720*. Center for Space Research, The University of Texas at Austin.

Bhattarai S, Ziebart M, Allgeier S, Grey S, Springer T, Harrison D, Li Z. 2019. Demonstrating developments in high-fidelity analytical radiation force modelling methods for spacecraft with a new model for GPS IIR/IIR-M. *J Geod* **93**: 1515–1528. https://doi.org/10.1007/s00190-019-01265-7.

Bhattarai S, Ziebart M, Springer T, Gonzalez F, Tobias G. 2022. High-precision physics-based radiation force models for the Galileo spacecraft. *Adv Space Res* **69**: 4141–4154. https://doi.org/10.1016/j.asr.2022.04.003.

Bird GA. 1994. *Molecular gas dynamics and the direct simulation of gas flows*. Oxford Science Publication, New York.

Bruinsma S, Boniface C. 2021. The operational and research DTM-2020 thermosphere models. *J Space Weather Space Clim* **11**, 47. https://doi.org/10.1051/swsc/2021032.

Bruinsma S, Siemes C, Emmert JT, Mlynczak MG. 2022. Description and comparison of 21st century thermosphere data. *Adv Space Res*. https://doi.org/10.1016/j.asr.2022.09.038.

Bruinsma S, Sutton E, Solomon SC, Fuller-Rowell T, Fedrizzi M. 2018. Space weather modeling capabilities assessment: Neutral density for orbit determination at low Earth orbit. *Space Weather* **16**: 1806–1816. https://doi.org/10.1029/2018SW002027.

Bruinsma S, Tamagnan D, Biancale R. 2004. Atmospheric densities derived from CHAMP/STAR accelerometer observations. *Planet Space Sci* **52**: 297–312. https://doi.org/10.1016/j.pss.2003.11.004.

Case K, Kruizinga G, Wu S-C. 2010. GRACE level 1B data product user handbook. *Technical report JPL D-2207*. Jet Propulsion Laboratory.

Christophe B, Boulanger D, Foulon B, Huynh P-A, Lebat V, Liorzou F, Perrot E. 2015. A new generation of ultra-sensitive electrostatic accelerometers for GRACE Follow-on and towards the next generation gravity missions. *Acta Astronaut* **117**: 1–7. https://doi.org/10.1016/j.actaastro.2015.06.021.

Darugna F, Steigenberger P, Montenbruck O, Casotto S. 2018. Ray-tracing solar radiation pressure modeling for QZS-1. *Adv Space Res* **62**: 935–943. https://doi.org/10.1016/j.asr.2018.05.036.

Dewitte S, Clerbaux N. 2017. Measurement of the earth radiation budget at the top of the atmosphere – a review. *Remote Sens* **9**(11), 1143. https://doi.org/10.3390/rs9111143.
Dobslaw H, Bergmann-Wolf I, Dill R, Poropat L, Thomas M, Dahle C, Esselborn S, König R, Flechtner F. 2017. A new high-resolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06. *Geophys J Int* **211**: 263–269. https://doi.org/10.1093/gji/ggx302.

Doornbos E. 2011. Thermospheric density and wind determination from satellite dynamics. *Ph.D. Thesis*, Department of Astrodynamics and Satellite Missions, Delft University of Technology. Available at http://resolver.tudelft.nl/uuid:33002be1-1498-4bec-a440-4c90ec149aea.

Doornbos E, Bruinsma S, Fritsche B, Koppenwallner G, Visser P, Van den IJssel J, de Teixeira de Encarnação J. 2014. ESA contract 4000102847/NL/EL, GOCE+ Theme 3: Air density and wind retrieval using GOCE data – final report, Technical Report. TU Delft. Available at https://earth.esa.int/eogateway/documents/20142/1181177/GOCE-theme-3-final-report.pdf.

Doornbos E, van den IJssel J, Lühr H, Förster M, Koppenwallner G. 2010. Neutral density and crosswind determination from arbitrarily oriented multiaxis accelerometers on satellites. *J Spacecr Rockets* **47**(4): 580–589. https://doi.org/10.2514/1.48114.

Drob DP, Emmert JT, Crowley G, Picone JM, Shepherd GG, et al. 2008. An empirical model of the Earth's horizontal wind fields: HWM07. *J Geophys Res Space Phys* **113**, A12304. https://doi.org/10.1029/2008JA013668.

Drob DP, Emmert JT, Meriwether JW, Makela JJ, Doornbos E, et al. 2015. An update to the Horizontal Wind Model (HWM): The quiet time thermosphere. *Earth Space Sci* **2**: 301–319. https://doi.org/10.1002/2014EA000089.

Floberghagen R, Fehringer M, Lamarre D, Muzi D, Frommknecht B, Steiger C, Piñeiro J, da Costa A. 2011. Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission. *J Geod* **85**: 749–758. https://doi.org/10.1007/s00190-011-0498-3.

Flury J, Bettadpur S, Tapley BD. 2008. Precise accelerometry onboard the GRACE gravity field satellite mission. *Adv Space Res* **42**: 1414–1423. https://doi.org/10.1016/j.asr.2008.05.004.

Fortescue P, Swinerd G, Stark J. 2011. *Spacecraft systems engineering*, 4th edn. Wiley, Chichester. ISBN 978-0-470-75012-4.

Hagan ME, Roble RG, Hackney J. 2001. Migrating thermospheric tides. *J Geophys Res Space Phys* **106**(A7): 12739–12752. https://doi.org/10.1029/2000JA000344.

Iorfida E, Daras I, Haagmans R, Strømme A. 2022. Swarm A and C accelerometers: Data validation and scientific interpretation. *Earth Space Sci* **10**, e2022EA002458. https://doi.org/10.1029/2022EA002458.

Kang Z, Bettadpur S, Tapley B, Cheng M, Ries J. 2003. Determination of CHAMP accelerometer calibration parameters. In: *First CHAMP mission results for gravity, magnetic and atmospheric studies*, Reigber C, Lühr H, Schwintzer P, (Eds.), Springer, Berlin, Heidelberg. pp. 19–25. https://doi.org/10.1007/978-3-540-38366-6_3.

Kenneally PW, Schaub H. 2020. Fast spacecraft solar radiation pressure modeling by ray tracing on graphics processing unit. *Adv Space Res* **65**: 1951–1964. https://doi.org/10.1016/j.asr.2019.12.028.

Klinger B, Mayer-Gürr T. 2016. The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSG-Grace2016. *Adv Space Res* **458**: 1597–1609. https://doi.org/10.1016/j.asr.2016.08.007.

Klinkrad H, Koeck C, Renard P. 1991. Key features of a satellite skin force modelling technique by means of Monte-Carlo ray tracing. *Adv Space Res* **11**(6): 147–150. https://doi.org/10.1016/0273-1177(91)90244-E.

Kodikara T. 2019. Physical understanding and forecasting of the thermospheric structure and dynamics, *Ph.D. Thesis*, RMIT University, Melbourne, Australia. Available at https://researchrepository.rmit.edu.au/esploro/outputs/9921863942601341.

König R, Michalak G, Neumayer K, Schmidt R, Zhu S, Meixner H, Reigber C. 2005. Recent developments in CHAMP orbit determination at GFZ. In: *Earth observation with CHAMP results from three years in orbit*, Reigber C, Lühr H, Schwintzer P, Wickert J, (Eds.), Springer, Berlin, Heidelberg. pp. 65–70. https://doi.org/10.1007/3-540-26800-6_10.

Krauss S, Behzadpour S, Temmer M, Lhotka C. 2020. Exploring thermospheric variations triggered by severe geomagnetic storm on 26 August 2018 using GRACE follow-on data. *J Geophys Res Space Phys* **125**, e2019JA027731. https://doi.org/10.1029/2019JA027731.

Li Z, Ziebart M, Bhattarai S, Harrison D, Grey S. 2020. Fast solar radiation pressure modelling with ray tracing and multiple reflections. *Adv Space Res* **61**: 2352–2365. https://doi.org/10.1016/j.asr.2018.02.019.

Lühr H, Rentz S, Ritter P, Liu H, Häusler K. 2007. Average thermospheric wind patterns over the polar regions, as observed by CHAMP. *Ann Geophys* **25**: 1093–1101. https://doi.org/10.5194/angeo-25-1093-2007.

March G, Doornbos E, Visser P. 2019a. High-fidelity geometry models for improving the consistency of CHAMP, GRACE, GOCE and Swarm thermospheric density data sets. *Adv Space Res* **63**: 213–238. https://doi.org/10.1016/j.asr.2018.07.009.

March G, Van den IJssel J, Siemes C, Visser P, Doornbos E, Pilinski M. 2021. Gas-surface interactions modelling influence on satellite aerodynamics and thermosphere mass density. *J Space Weather Space Clim* **11**, 54. https://doi.org/10.1051/swsc/2021035.

March G, Visser T, Visser P, Doornbos E. 2019b. CHAMP and GOCE thermospheric wind characterization with improved gas-surface interactions modelling. *Adv Space Res* **64**: 1225–1242. https://doi.org/10.1016/j.asr.2019.06.023.

McCarthy D, Petit G. 2004. *IERS conventions (2003)*, IERS Technical Note 32. Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main. ISBN 3-89888-884-3.

McCullough CM, Harvey N, Save H, Bandikova T. 2019. Description of calibrated GRACE-FO accelerometer data products (ACT). *Technical report JPL D-103863*. NASA Jet Propulsion Laboratory, California Institute of Technology.

McGirr R, Tregoning P, Allgeyer S, McQueen H, Purcell A. 2022. Mitigation of thermal noise in GRACE accelerometer observations. *Adv Space Res* **69**: 386–401. https://doi.org/10.1016/j.asr.2021.10.055.

Mehta PM, Paul SN, Crisp NH, Sheridan PL, Siemes C, March G, Bruinsma S. 2022. Satellite drag coefficient modeling for thermosphere science and mission operations. *Adv Space Res*. https://doi.org/10.1016/j.asr.2022.05.064.

Mehta PM, Walker AC, Sutton EK, Godinez HC. 2017. New density estimates derived using accelerometers on board the CHAMP and GRACE satellites. *Space Weather* **15**(4): 558–576. https://doi.org/10.1002/2016SW001562.

Montenbruck O, Gill E. 2012. *Satellite orbits*. Springer. ISBN 978-3-540-67280-7. https://doi.org/10.1007/978-3-642-58351-3.

Montenbruck O, Steigenberger P, Hugentobler U. 2015. Enhanced solar radiation pressure modeling for Galileo satellites. *J Geod* **89**: 283–297. https://doi.org/10.1007/s00190-014-0774-0.

Murböck M, Abrykosov P, Dahle C, Hauk M, Pail R, Flechtner F. 2023. In-orbit performance of the GRACE accelerometers and microwave ranging instrument. *Remote Sens* **15**, 563. https://doi.org/10.3390/rs15030563.

Olsen N, Friis-Christensen E, Floberghagen R, Alken P, Beggan CD, et al. 2013. The Swarm Satellite Constellation Application and Research Facility (SCARF) and Swarm data products. *Earth Planet Space* **65**(11): 1189–1200. https://doi.org/10.5047/eps.2013.07.001.
Palmroth M, Grandin M, Sarris T, Doornbos E, Tourgaidis S, et al. 2021. Lower thermosphere – ionosphere (LTI) quantities: Current status of measuring techniques and models. *Ann Geophys* **39**: 189–237. https://doi.org/10.5194/angeo-39-189-2021.

Pavlis D, Poulouse S, McCarthy J. 2006. *GEODYN operations manual*. SGT Inc., Greenbelt.

Perosanz F, Biancale R, Lemoine J, Vales N, Loyer S, Bruinsma S. 2005. Evaluation of the CHAMP accelerometer on two years of mission. In: *Earth observation with CHAMP results from three years in orbit*, Reigber C, Lühr H, Schwintzer P, Wickert J, (Eds.), Springer, Berlin. pp. 77–82. https://doi.org/10.1007/3-540-26800-6_10.

Picone JM, Hedin AE, Drob DP, Aikin AC. 2002. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues. *J Geophys Res Space Phys* **107**(A12): 1468. https://doi.org/10.1029/2002JA009430.

Prange L. 2010. Global gravity field determination using the GPS measurements made onboard the low Earth orbiting satellite CHAMP. *Geodätisch-geophysikalische Arbeiten in der Schweiz*, **81**, Schweizerische Geodätische Kommission, Zürich, Switzerland. Available at http://www.sgc.ethz.ch/sgc-volumes/sgk-81.pdf.

Qian L, Burns AG, Emery BA, Foster B, Lu G, Maute A, Richmond AD, Roble RG, Solomon SC, Wang W. 2014. The NCAR TIE-GCM: A community model of the coupled thermosphere/ionosphere system. In: *Modeling the ionosphere thermosphere system*, Huba J, Schunk R, Khazanov G, (Eds.), John Wiley, Washington. pp. 73–83ISBN: 9781118704417. https://doi.org/10.1002/9781118704417.ch7.

Qian L, Solomon SC, Kane TJ. 2009. Seasonal variation of thermospheric density and composition. *J Geophys Res Space Phys* **114**(A1): A01312. https://doi.org/10.1029/2008JA013643.

Ray R. 1999. A global ocean tide model from TOPEX/POSEIDON altimetry: GOT99.2. *Technical Report NASA/TM1999-209478*. NASA Goddard Space Flight Center.

Reigber C, Lühr H, Schwintzer P. 2002. CHAMP mission status. *Adv Space Res* **30**(2): 129–134. https://doi.org/10.1016/S0273-1177(02)00276-4.

Richmond AD, Ridley EC, Roble RG. 1992. A thermosphere/ionosphere general circulation model with coupled electromagnics. *Geophys Res Lett* **19**(6): 601–604. https://doi.org/10.1029/92GL00401.

Sentman LH. 1961. Free molecule flow theory and its application to the determination of aerodynamic forces. *Technical Report*. Lockheed Missiles and Space Co Inc, Sunnyvale, CA.

Siemes C. 2020. Swarm satellite thermo-optical properties and external geometry. *Technical Report ESA-EOPG-MOM-MO-15*. European Space Agency. Available at https://earth.esa.int/eogate-way/documents/20142/37627/swarm-thermo-optical-properties-and-external-geometry.pdf.

Siemes C, de Teixeira J, da Encarnação E, Doornbos J, van den IJssel J, Kraus R, Pereštý L, Grunwaldt G, Apelbaum J Flury, Holmdahl Olsen PE. 2016. Swarm accelerometer data processing from raw accelerations to thermospheric neutral densities. *Earth Planet Space* **68**: 92. https://doi.org/10.1186/s40623-016-0474-5.

Standish E. 1998. *JPL planetary and lunar ephemerides, DE405/LE405*. JPL IOM 312.F-98-048. ftp://ssd.jpl.nasa.gov/pub/eph planets/ioms/de405.iom.pdf.

Sutton EK. 2008. Effects of solar disturbances on the thermosphere densities and winds from CHAMP and GRACE satellite accelerometer data. *Ph.D. Thesis*, Department of Aerospace Engineering Sciences, University of Colorado at Boulder. Available at https://www.proquest.com/docview/304639074.

Sutton EK. 2018. A new method of physics-based data assimilation for the quiet and disturbed thermosphere. *Space Weather* **16**: 736–753. https://doi.org/10.1002/2017SW001785.

Sutton EK, Nerem RS, Forbes JM. 2007. Density and winds in the thermosphere deduced from accelerometer data. *J Spacecr Rockets* **44**(6): 1210–1219. https://doi.org/10.2514/1.28641.

Tapley BD, Bettadpur S, Watkins M, Reigber C. 2004. The gravity recovery and climate experiment: Mission overview and early results. *Geophys Res Lett* **31**(9). L019920. https://doi.org/10.1029/2004GL019920.

Touboul P, Foulon B, Christophe B, Marque JP. 2012. CHAMP, GRACE, GOCE Instruments and Beyond. In: *Geodesy for Planet Earth*, vol 136, Kenyon S, Pacino M, Marti U, (Eds.), Springer. pp. 215–221. https://doi.org/10.1007/978-3-642-20338-1_26.

Touboul P, Métris G, Sélig H, Le Traon O, Bresson A, Zahzam N, Christophe B, Rodrigues M. 2016. Gravitation and geodesy with inertial sensors, from ground to space. *Aerospace Lab* **12**, AL12-11. https://doi.org/10.12762/2016.AL12-11.

Vallado DA, Finkleman D. 2014. A critical assessment of satellite drag and atmospheric density modeling. *Acta Astronaut* **95**: 141–165. https://doi.org/10.1016/j.actaastro.2013.10.005.

Van den IJssel J, Doornbos E, Iorfida E, March G, Siemes C, Montenbruck O. 2020. Thermosphere densities derived from Swarm GPS observations. *Adv Space Res* **65**: 1758–1771. https://doi.org/10.1016/j.asr.2020.01.004.

Van Helleputte T. 2011. The integration of spaceborne accelerometry in the precise orbit determination of low-flying satellites. *Ph.D. Thesis*, Department of Astrodynamics and Satellite Missions, Delft University of Technology. Available at http://resolver.tudelft.nl/uuid:c84217a2-967c-45c4-85cb-a1630590b926.

van Helleputte T, Doornbos E, Visser P. 2009. CHAMP and GRACE accelerometer calibration by GPS-based orbit determination. *Adv Space Res* **43**: 1890–1896. https://doi.org/10.1016/j.asr.2009.02.017.

Vielberg K, Forootan E, Lück C, Löcher A, Kusche J, Börger K. 2018. Comparison of accelerometer data calibration methods used in thermospheric neutral density estimation. *Ann Geophys* **36**(3): 761–779. https://doi.org/10.5194/angeo-36-761-2018.

Vielberg K, Kusche J. 2020. Extended forward and inverse modeling of radiation pressure accelerations for LEO satellites. *J Geod* **94**, 43. https://doi.org/10.1007/s00190-020-01368-6.

Visser P, van den IJssel J. 2016. Calibration and validation of individual GOCE accelerometers by precise orbit determination. *J Geod* **90**: 1–13. https://doi.org/10.1007/s00190-015-0850-0.

Wen HY, Kruizinga G, Paik M, Landerer F, Bertiger W, Sakumura C, Bandikova T, McCullough C. 2019. Gravity recovery and climate experiment follow-on (GRACE-FO) level-1 data product user handbook, Technical report JPL D-56935 (URS270772). NASA Jet Propulsion Laboratory, California Institute of Technology.

Wö ske F, Kato T, Rievers B, List M. 2019. GRACE accelerometer calibration by high precision non-gravitational force modeling. *Adv Space Res* **63**: 1318–1335. https://doi.org/10.1016/j.asr.2018.10.025.

Zehentner N, Mayer-Gürr T. 2016. Precise orbit determination based on raw GPS measurements. *J Geod* **90**: 275–286. https://doi.org/10.1007/s00190-015-0872-7.

Ziebart M. 2004. Generalized analytical solar radiation pressure modeling algorithm for spacecraft of complex shape. *J Spacecr Rockets* **41**(5): 840–848. https://doi.org/10.2514/1.13097.
# Appendix

## A.1 Satellite thermal models

<figure><img src="image_16.png" /><figcaption>Fig. A1. CHAMP accelerometer biases. We show only periods flagged as valid.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. A2. GRACE B accelerometer biases. We show only periods flagged as valid.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. A3.** GRACE C accelerometer biases. We show only periods flagged as valid.</figcaption></figure>

<table><caption>**Table A1.** CHAMP thermal radiation panel model. We used an internal heat generation $P_{gen} = 30$ W and a heat capacity $C_{body} = 100$ kJ·K⁻¹.</caption><thead><tr><th>Panel</th><th>A<sub>j</sub> (m<sup>2</sup>)</th><th>n<sub>x,j</sub> (-)</th><th>n<sub>y,j</sub> (-)</th><th>n<sub>z,j</sub> (-)</th><th>c<sub>vis,a</sub> (-)</th><th>c<sub>ir,a</sub> (-)</th><th>C<sub>j</sub> (J·K<sup>-1</sup>)</th><th>k<sub>j</sub> (W·K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Zenith</td><td>1.1500</td><td>0.0</td><td>0.0</td><td>-1.0</td><td>0.65</td><td>0.81</td><td>1500.0</td><td>0.01</td></tr><tr><td>Nadir</td><td>4.0000</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0.12</td><td>0.75</td><td>7000.0</td><td>0.5</td></tr><tr><td>Port (outer)</td><td>3.1000</td><td>0.0</td><td>-0.707107</td><td>-0.707107</td><td>0.65</td><td>0.81</td><td>5000.0</td><td>0.01</td></tr><tr><td>Starboard (outer)</td><td>3.1000</td><td>0.0</td><td>0.707107</td><td>-0.707107</td><td>0.65</td><td>0.81</td><td>5000.0</td><td>0.01</td></tr><tr><td>Front</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sandblasted Al</td><td>1.2199</td><td>0.3420</td><td>0.0</td><td>0.9397</td><td>0.40</td><td>0.00</td><td>2000.0</td><td>0.01</td></tr><tr><td>Sandblasted Al</td><td>0.0880</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.40</td><td>0.00</td><td>100.0</td><td>0.01</td></tr><tr><td>Gold foil</td><td>0.1314</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.00</td><td>200.0</td><td>0.01</td></tr><tr><td>Sandblasted Al</td><td>0.1500</td><td>0.7071</td><td>0.0</td><td>0.7071</td><td>0.34</td><td>0.00</td><td>200.0</td><td>0.01</td></tr><tr><td>Rear</td><td>0.8000</td><td>-0.9397</td><td>0.0</td><td>0.3420</td><td>0.34</td><td>0.62</td><td>1500.0</td><td>0.01</td></tr><tr><td>Boom</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Top (Kapton/SiO<sub>2</sub>)</td><td>0.5110</td><td>0.0</td><td>0.0</td><td>-1.0</td><td>0.34</td><td>0.62</td><td>800.0</td><td>0.0001</td></tr><tr><td>Top (Kapton/Al)</td><td>0.4190</td><td>0.0</td><td>0.0</td><td>-1.0</td><td>0.16</td><td>0.00</td><td>700.0</td><td>0.0001</td></tr><tr><td>Bottom</td><td>0.9300</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0.34</td><td>0.62</td><td>1400.0</td><td>0.0001</td></tr><tr><td>Right</td><td>0.9300</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>1400.0</td><td>0.0001</td></tr><tr><td>Left</td><td>0.9300</td><td>0.0</td><td>-1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>1400.0</td><td>0.0001</td></tr></tbody></table>
<table><caption>**Table A2.** GRACE thermal radiation panel model. We used an internal heat generation $P_{gen} = 70$ W and a heat capacity $C_{body} = 100$ kJ·K⁻¹.</caption><thead><tr><th>Panel</th><th>A<sub>j</sub> (m<sup>2</sup>)</th><th>n<sub>x,j</sub> (-)</th><th>n<sub>y,j</sub> (-)</th><th>n<sub>z,j</sub> (-)</th><th>c<sub>vis,a</sub> (-)</th><th>c<sub>ir,a</sub> (-)</th><th>C<sub>j</sub> (J·K<sup>-1</sup>)</th><th>k<sub>j</sub> (W·K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Zenith</td><td>2.1673620</td><td>0.0</td><td>0.0</td><td>-1.0</td><td>0.72</td><td>0.81</td><td>5000.0</td><td>0.1</td></tr><tr><td>Nadir</td><td>6.0711120</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0.12</td><td>0.75</td><td>10,000.0</td><td>0.5</td></tr><tr><td>Port (outer)</td><td>3.1554792</td><td>0.0</td><td>-0.766044</td><td>-0.642787</td><td>0.72</td><td>0.81</td><td>5000.0</td><td>0.1</td></tr><tr><td>Port (inner)</td><td>0.2282913</td><td>0.0</td><td>0.766044</td><td>0.642787</td><td>0.34</td><td>0.62</td><td>500.0</td><td>0.1</td></tr><tr><td>Starboard (outer)</td><td>3.1554792</td><td>0.0</td><td>0.766044</td><td>-0.642787</td><td>0.72</td><td>0.81</td><td>5000.0</td><td>0.1</td></tr><tr><td>Starboard (inner)</td><td>0.2282913</td><td>0.0</td><td>-0.766044</td><td>0.642787</td><td>0.34</td><td>0.62</td><td>500.0</td><td>0.1</td></tr><tr><td>Front (panel)</td><td>0.9551567</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>1000.0</td><td>0.1</td></tr><tr><td>Rear (panel)</td><td>0.9551567</td><td>-1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>1000.0</td><td>0.1</td></tr><tr><td>Front (antenna)</td><td>0.0461901</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Rear (antenna)</td><td>0.0461901</td><td>-1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Port (antenna)</td><td>0.0461901</td><td>0.0</td><td>-1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Starboard (antenna)</td><td>0.0461901</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr></tbody></table>

<table><caption>**Table A3.** GRACE-FO thermal radiation panel model. The absorption coefficients of the panels “Port (outer)” and “Starboard (outer)” are taken from the Swarm mission documentation (Siemes, 2020) because the solar arrays on these panels resemble those of the Swarm satellites rather than those of the GRACE satellites on photos. We used an internal heat generation $P_{gen} = 55$ W and a heat capacity $C_{body} = 100$ kJ·K⁻¹.</caption><thead><tr><th>Panel</th><th>A<sub>j</sub> (m<sup>2</sup>)</th><th>n<sub>x,j</sub> (-)</th><th>n<sub>y,j</sub> (-)</th><th>n<sub>z,j</sub> (-)</th><th>c<sub>vis,a</sub> (-)</th><th>c<sub>ir,a</sub> (-)</th><th>C<sub>j</sub> (J·K<sup>-1</sup>)</th><th>k<sub>j</sub> (W·K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Zenith</td><td>2.1673620</td><td>0.0</td><td>0.0</td><td>-1.0</td><td>0.65</td><td>0.81</td><td>12,000.0</td><td>0.1</td></tr><tr><td>Nadir</td><td>6.0711120</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0.12</td><td>0.75</td><td>10,000.0</td><td>0.5</td></tr><tr><td>Port (outer)</td><td>3.1554792</td><td>0.0</td><td>-0.766044</td><td>-0.642787</td><td>0.90</td><td>0.81</td><td>18,000.0</td><td>0.1</td></tr><tr><td>Port (inner)</td><td>0.2282913</td><td>0.0</td><td>0.766044</td><td>0.642787</td><td>0.34</td><td>0.62</td><td>800.0</td><td>0.1</td></tr><tr><td>Starboard (outer)</td><td>3.1554792</td><td>0.0</td><td>0.766044</td><td>-0.642787</td><td>0.90</td><td>0.81</td><td>18,000.0</td><td>0.1</td></tr><tr><td>Starboard (inner)</td><td>0.2282913</td><td>0.0</td><td>-0.766044</td><td>0.642787</td><td>0.34</td><td>0.62</td><td>800.0</td><td>0.1</td></tr><tr><td>Front (panel)</td><td>0.9551567</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>4000.0</td><td>0.1</td></tr><tr><td>Rear (panel)</td><td>0.9551567</td><td>-1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>4000.0</td><td>0.1</td></tr><tr><td>Front (antenna)</td><td>0.0461901</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Rear (antenna)</td><td>0.0461901</td><td>-1.0</td><td>0.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Port (antenna)</td><td>0.0461901</td><td>0.0</td><td>-1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr><tr><td>Starboard (antenna)</td><td>0.0461901</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.34</td><td>0.62</td><td>400.0</td><td>0.01</td></tr></tbody></table>

Cite this article as: Siemes C, Borries C, Bruinsma S, Fernandez-Gomez I, Hladczuk N, et al. 2023. New thermosphere neutral mass density and crosswind datasets from CHAMP, GRACE, and GRACE-FO. *J. Space Weather Space Clim.* **13**, 16. https://doi.org/10.1051/swsc/2023014.