Issue 
J. Space Weather Space Clim.
Volume 11, 2021



Article Number  54  
Number of page(s)  23  
DOI  https://doi.org/10.1051/swsc/2021035  
Published online  26 October 2021 
Research Article
Gassurface interactions modelling influence on satellite aerodynamics and thermosphere mass density
^{1}
Department of Astrodynamics and Space Missions, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, Delft, 2629 HS, The Netherlands
^{2}
Royal Netherlands Meteorological Institute, Utrechtseweg 297, De Bilt, 3731 GA, The Netherlands
^{3}
University of Colorado Boulder, 3775 Discovery Drive, Boulder, 80309 CO, USA
^{4}
RHEA for ESA – European Space Agency, Noordwijk, 2201 AZ, The Netherlands
^{*} Corresponding author: gunther.march@esa.int, g.march@hotmail.it
Received:
3
May
2021
Accepted:
7
September
2021
The satellite acceleration data from the CHAMP, GRACE, GOCE, and Swarm missions provide detailed information on the thermosphere density over the last two decades. Recent work on reducing errors in modelling the spacecraft geometry has greatly reduced scale differences between the thermosphere data sets from these missions. However, residual inconsistencies between the data sets and between data and models are still present. To a large extent, these differences originate in the modelling of the gassurface interactions (GSI), which is part of the satellite aerodynamic modelling used in the acceleration to density data processing. Physicsbased GSI models require insitu atmospheric composition and temperature data that are not measured by any of the abovementioned satellites and, as a consequence, rely on thermosphere models for these inputs. To reduce the dependence on existing thermosphere models, we choose a GSI model with a constant energy accommodation coefficient per mission, which we optimize exploiting particular attitude manoeuvres and wind analyses to increase the selfconsistency of the multimission thermosphere mass density data sets. We compare our results with those based on variable energy accommodation obtained by different studies and semiempirical models to show the principal differences. The presented comparisons provide novel opportunity to quantify the discrepancies between current GSI models. Among the presented data, density variations with variable accommodation are within ±10%, and peaks can reach up to 15% at the poles. The largest differences occur during low solar activity periods. In addition, we utilize a series of attitude manoeuvres performed in May 2014 by the Swarm A and C satellites, which are flying in close proximity, to evaluate the residual inconsistency of the density observations as a function of the energy accommodation coefficient. Our analysis demonstrates that an energy accommodation coefficient of 0.85 maximizes the consistency of the Swarm density observations during the attitude manoeuvres. Using such coefficient, for Swarm A and Swarm C, the new density would be lower in magnitude with a 4–5% difference. In recent studies, similar energy accommodation coefficients were retrieved for the CHAMP and GOCE missions by investigating thermospheric winds. These new values for the energy accommodation coefficient provide a higher consistency among different missions and models. A comparison of neutral densities between current thermosphere models and observations indicates that semiempirical models such as NRLMSISE00 and DTM2013 significantly overestimate the density, and that an overall higher consistency between the observations from the different missions can be achieved with the presented assumptions. The new densities from this work provide consistencies of 4.13% and 3.65% between the minimum and maximum mean ratios among the selected missions with NRLMSISE00 and DTM2013, respectively. A comparison with the WACCMX general circulation model is also performed. Similar to the other models, WACCMX seems to provide higher estimates of mass density especially under high and moderate solar activities. This work has the objective to guide density data users over the multiple data sets and highlight the remaining uncertainties associated with different GSI models.
Key words: thermosphere / mass density / satellite aerodynamics / gassurface interactions / atmospheric drag / accelerometer
© G. March et al., Published by EDP Sciences 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The launch of the Challenging Minisatellite Payload (CHAMP) satellite in 2000 marked a new era in which accelerometer measurements were used for producing highresolution and nearly continuous thermosphere density data sets. More satellite missions carrying precise accelerometers into a lowEarth orbit followed. The Gravity Recovery and Climate Experiment (GRACE), the Gravity field and steadystate Ocean Circulation Explorer (GOCE), Swarm, and the GRACE FollowOn (GRACEFO) missions were launched in 2002, 2009, 2013, and 2018, respectively. Though only the Swarm mission includes observing thermosphere density in its mission objectives, all of the beforementioned missions have provided valuable thermosphere density data sets. Among these missions, the initial altitude ranged between 270 km (for GOCE) and 530 km (for Swarm B). CHAMP and the other two Swarm satellites (Swarm A and Swarm C) were placed at 460 km and 470 km, respectively. The two GRACE satellites were initially around 505 km. Most of these missions contributed to creating density data sets which were initiated by Bowman et al. (2004) and Storz et al. (2005) and followed by Sutton (2008) (http://sisko.colorado.edu/sutton/data/ver2), Calabia & Jin (2016) (https://zenodo.org/record/4308315) and Doornbos (2012) and March et al. (2019a) (http://thermosphere.tudelft.nl).
Early generations of CHAMP and GRACE density data sets were based on simplified satellite geometry descriptions and idealised firstguess gassurface interaction (GSI) parameters (e.g., Bruinsma & Biancale, 2003; Sutton, 2008; Doornbos, 2011; Mehta et al., 2013). In particular, the energy accommodation coefficient model generated by Mehta et al. (2013) was used in Calabia & Jin (2016) to estimate densities. However, over the last years, a research effort towards improving geometry and rarefied gasdynamics modelling has been made in order to raise the level of accuracy of the density data sets (Pilinski et al., 2013a, 2016; Mehta et al., 2017; March et al., 2019a). New models of satellite surface geometries have been constructed to improve over the previously used models, which lacked the definition of geometry features such as baffles and antennas and physical features such as multiple reflections of atmospheric particles on the surfaces and shadowing of surfaces. These shortcomings had introduced large errors in the scaling of the individual density data sets (up to 32% in the case of Swarm), which led to readily apparent inconsistencies between data processed by different providers and between different missions (March et al., 2019a). However, the scale differences between data sets and atmospheric models are not exclusively due to geometry modelling. The characterization of particlesurface collisions between atmospheric molecules and satellite surfaces is fundamental (March et al., 2019b).
The energy accommodation coefficient (α_{E}) is an important parameter for GSI modelling. It describes the energy exchange between the atmospheric molecules and the satellite surfaces. Theoretically, the parameter can range between 1 and 0. The two extreme values (α_{E} = 1 and α_{E} = 0) describe particlesurface reflections with and without accommodation of the particle temperature to the temperature of the satellite surface, respectively. Depending on different level of adsorption of specific atmospheric constituents on satellite surfaces, GSI can drastically change (Pilinski et al., 2013a). The relative atmospheric concentrations of molecular nitrogen, atomic oxygen, and atomic helium play a relevant role in satellite aerodynamics (Mehta et al., 2019). The relative concentrations can vary significantly over the course of a solar cycle, and at shorter timescales, subject to solar and geomagnetic driving of the upper atmosphere.
Laboratory experiments under thermospheric conditions and insitu data for GSI in such regions are extremely sparse and often limited in practical usability due to the lack of auxiliary data and use of underlying assumptions. Currently, the most advanced approach to model GSI on satellites in a physically realistic way requires empirical relations between the atomic oxygen concentration and temperature of the gas and the energy accommodation coefficient (Pilinski et al., 2013a). However, insitu concentration and temperature observations are not available for any of the abovementioned satellites, so this approach in the data processing from acceleration to thermosphere densities relies completely on the use of existing thermosphere models for these input parameters. In addition, the parameters used in the aforementioned empirical relations are fitted to past satellite observations that were processed, making use of past thermosphere models as well. Therefore, accelerometer data processed with this physicsbased GSI approach depends on multiple previous thermosphere models and satellite data sets in complex ways, making it challenging to attribute and reduce the residual inconsistencies between the more modern data sets and models.
Reducing the current scale differences between the thermosphere data sets is crucial for many reasons. Missiondependent biases, if not taken accurately into account, can cause problems in data assimilation, both in the empirical and physicsbased models (e.g. Matsuo et al., 2012; Bruinsma, 2015; Mehta & Linares, 2018; Sutton, 2018). Even though in some cases bias differences are estimated and removed in the assimilation process, information on the absolute scale of the modeled thermosphere density cannot be recovered in this way, and inconsistencies with other thermosphere–ionosphere parameters could remain. Consequently, the biases also affect investigations of the thermosphereionosphere system’s energy budget (e.g., Wilson et al., 2006). Biasfree data will also be essential for investigating longterm trends over multiple solar cycles, induced by changes in solar activity and by cooling the upperatmosphere due to increasing concentrations of greenhouse gases (e.g., Emmert et al., 2008; Qian et al., 2013). With the continuation of the Swarm and GRACEFO missions, the data sets discussed in this paper will ultimately span more than two full solar cycles. This data is the beginning of large time series for which such intermission biases are investigated and from which longterm trends and their uncertainties can be derived. Longterm trends strongly affect estimates of an inorbit lifetime of future satellite missions, as well as the longterm evolution of the low Earth orbit space debris (Lewis et al., 2011). Moreover, orbit predictions of current and future missions will benefit from improved knowledge of satellite aerodynamics and upper atmospheric variability, which affects orbital lifetime and manoeuvre planning.
In this study, we investigate the influence of the GSI modelling on the consistency of thermosphere density data sets from the CHAMP, GRACE, GOCE, and Swarm missions. This research follows the previously published analyses on neutral wind products from CHAMP and GOCE that used a similar approach (March et al., 2019b). Similar to what was done for CHAMP and GOCE, our goal is to selfconsistently analyse and process these data sets and assess the results in combination with thermosphere models that were evaluated along the satellite trajectories. Due to limitations in the observation data and data/parameter relations for adsorbed constituents nearby satellite surfaces, the use of the abovementioned physicsbased variable accommodation GSI approach, based on oxygen concentration and temperature inputs, was not possible with our approach. For this reason, we chose a simpler GSI model with one constant energy accommodation coefficient per mission, which we treat as a free parameter to increase the consistency of the thermosphere density data sets. A comparison with the output of the variable accommodation coefficients models is performed to quantify the different results introduced by the two methods. When we use the term “selfconsistency”, we refer to studies of GSI effects on single mission products which turned out to provide a higher level of consistency for specific thermospheric products (e.g., density, wind) within short or long time windows (e.g., manoeuvres, seasons).
A similar approach to this study was selected in March et al. (2019b) for analysing the effect of GSI on thermosphere wind from the CHAMP and GOCE missions. Therein, the GSI model is based on diffuse reflections with incomplete accommodation (DRIA), which was adopted according to the current processing algorithms of the official GOCE and Swarm missions’ products. However, by assuming this constant incomplete accommodation of the reflected particles, we ignore possible variations of energy accommodation due to changes in atmospheric temperature and composition as a function of satellite altitude, solar and geomagnetic activity. This limitation, imposed by the currently available measurement data and inherent in the chosen approach, should be kept in mind by users of the data.
In particular, for CHAMP, an attitude manoeuvre performed in November 2002 provided detailed information on the energy accommodation coefficient. Studying the consistency of the thermosphere wind within similar orbits and magnetic field conditions, but using different satellite orientations, a higher selfconsistency was achieved for α_{E} = 0.85. Zonal winds within the analysed attitude manoeuvre had a lower variability under similar conditions for this new α_{E} value. Also, for GOCE, a study based on seasonal dependency resulted in a lower energy accommodation coefficient than the currently adopted value of 0.93 at TU Delft and for part of the official ESA mission products. Indeed, a greater consistency was achieved for a coefficient of 0.82. This value is currently adopted in the new release of GOCE+ data (https://earth.esa.int/eogateway/catalog/gocethermospheredata?text=goce+thermosphere+data+2.0; August, 2021).
Attitude manoeuvres are not common, and, in particular, the ones which guarantee a stable flight configuration without thruster activation, or with sufficient time windows between different attitude orientations are even more scarce. In order to provide reliable information, these manoeuvres need to be in periods of high to medium solar activity and at relatively low altitudes (i.e., below 450–500 km) to ensure a strong aerodynamic acceleration signal in relation to error sources such as radiation pressure and instrument bias mismodelling. Investigating these manoeuvres is particularly useful for thermosphere density and wind studies.
Neutral density is more sensitive to variations in the solar extreme ultraviolet (EUV) emissions. For this reason, studying densities to optimise the energy accommodation coefficient is more challenging and requires satellites with characteristics as similar as possible. The Swarm mission provides an excellent opportunity for such a comparison. Indeed, the Swarm A and C satellites are flying sidebyside at the same altitude (between 450 and 500 km) with up to 1.4° separation in longitude over the equator and 4–10 s (30–75 km) separation in alongtrack direction. This turns out to provide nearly identical density observations. If within a certain time window, these satellites perform attitude manoeuvres exposing a different side of their body to the atmospheric flow, it is possible to intercompare the two data sets and deduce information on GSI modelling parameters like the energy accommodation coefficient. This will be further explained in Section 3.3.
Beyond the introduction of the DRIA neutral densities, for the CHAMP and GRACE satellites, the results are compared with other GSI models based on variable accommodation coefficients. Among the external sources, the data from Mehta et al. (2017) and the SESAM empirical model Pilinski et al. (2013b) are used to highlight the differences between different methods. This additional study allows to enhance the understanding of the differences introduced by fixed accommodation, and offers the opportunity to intercompare different variable accommodation coefficient models. The two opportunities provide crucial input for GSI modelling studies, especially for improving thermosphere and space weather modelling capabilities.
Over the last years, numerous studies have been performed on satellite aerodynamics by Monte Carlo Test Particle (MCTP) techniques, and there is an increasing interest in processing satellite data with highfidelity geometries. The SPARCS software (Pilinski, 2011), based on the test particle technique, analyzes triangulated meshes to provide aerodynamic coefficients. In Mehta et al. (2017), a similar investigation was performed for the CHAMP and GRACE satellites assuming a variable energy accommodation coefficient. Throughout this work, drag coefficients are based on an extension of the Monte Carlo technique presented in March et al. (2019a) based on the Stochastic PArallel Rarefiedgas Timeaccurate Analyzer (SPARTA) simulator (Gallis et al., 2014). The density data resulting from this work have been made available to aid further development and validation of thermosphere models (e.g., NRLMSISE (Picone et al., 2001), DTM (Bruinsma, 2015)) and enhance space mission operations analysis and planning. The data sets are also made available at http://thermosphere.tudelft.nl (August 2021).
The adopted methodology is summarized in Section 2. The results about the GSI influence on satellite aerodynamics and the effects on mass density can be found in Section 3. Comparisons with external data sets and models are available in Section 4. Section 5 provides conclusions and an outlook for future work.
2 Methodology
This study is based on the simulation output from the Stochastic PArallel Rarefiedgas Timeaccurate Analyzer (SPARTA) software (Gallis et al., 2014) (https://sparta.github.io, August 2021) in combination with the new highfidelity satellite geometries from March et al. (2019a). At the analysed altitudes, the atmospheric flow is considered to be in free molecular regime due to the large mean free path. Therefore, collisions between particles are neglected because of the large distance between themselves (of the order of many meters). This allowed for a simplification resulting in faster convergence of the simulations, which were anyway tested with and without particleparticle collisions as validation of such assumption. From technical drawings and prelaunch pictures, triangulated geometries for all the selected satellites were created using CAD software. Afterwards, these new satellite geometries were given as input to SPARTA in order to calculate with a Monte Carlo approach the resulting forces under different conditions. In particular, the speed ratio was a fundamental parameter because it allowed to simulate different satellite velocities, but also different chemical compositions and atmospheric temperatures. From equation (1), it is possible to notice that this parameter is directly connected with the relative velocity between satellite and atmosphere (v_{inc}), local atmospheric temperature (T_{inc}), molecular mass (m) and the gas constant (R):
In the generated data sets, this ratio ranges between 0.5 and 14, and this interval fully describes all possible encountered mission scenarios. Furthermore, a wide range of attitude angles with respect to the incoming flow, satellite velocity, and chemical compositions was also simulated. Regarding the chemistry, most of the inputs were generated using the NRLMSISE00 model. An atmospheric temperature of 1000 K and a satellite surface temperature of 400 K were chosen among other initial inputs. The atmospheric composition was assumed to be 100% atomic oxygen. However, all the mentioned inputs were further modified according to the simulated speed ratios, which allowed to go beyond the mentioned inputs and span through different mass compositions and temperatures. The main simulation settings are fully described in March et al. (2019a). However, as a difference with respect to the previous work, the gassurface collisions are not assumed to be fully accommodated reflections (α_{E} = 1), which represent the case in which the temperature of the impinging particles on the satellite adjust to the surface temperature. Indeed, in this work, we allow an incomplete accommodation of the temperature, i.e., the energy accommodation coefficient α_{E} is allowed to differ from one. In this research, the energy accommodation coefficient is the key parameter to characterize the GSI effects. Indeed, this parameter provides tangible information on the energy exchange between atmospheric particles and satellite surfaces. Its value can be estimated from the incoming and reflected kinetic temperatures (T_{k,i} and T_{k,r}) and the satellite surface temperature (T_{wall}) according to:
(2)and may assume values ranging from 0 to 1. The incoming kinetic temperature is related to the molar mass (m) and the incoming velocity (v_{inc}) as,
(3)where k_{B} is Boltzmann’s constant. The kinetic temperature of the reflected particles can be determined from equation (2) by,
The lack of measurements or models for the wall temperature introduces uncertainty in calculations. However, since T_{wall} ≪ T_{k,i}, the sensitivity to this parameter is low and does not particularly affect aerodynamic computations. SPARTA did not have the capability to define the energy accommodation coefficient as defined in equation (2), therefore equation (4) was used to describe α_{E} through the modification of the wall temperature. Due to the expensive computational cost of the simulations, data sets were generated for four different energy accommodation coefficients (0.60, 0.80, 0.93, 1.00). Afterwards, the new data sets for the additional α_{E} values (beyond the four simulated conditions) have been generated in the postprocessing phase using a leastsquares method. The processing was further validated with a few new SPARTA simulations generating results for a new set of α_{E} values. The satellite surface is covered with a variety of materials, unfortunately, we could not associate different surface properties to different parts of the triangulated geometry due to software limitations. Therefore, GSI is assumed to be independent of the different surface materials.
From a few experimental observations (Hedin et al., 1973; Gregory & Peters, 1987; Moe et al., 1993), the α_{E} value is suggested to be closer to unity than zero. Unfortunately, the reliability of older observation methods is uncertain due to instrument limitations. The chemical constituents play a crucial role in the adsorption properties of the satellite surfaces. The amount of adsorbed chemical molecules on the satellite surface can be implemented in the Cercignani–Lampis–Lord (CLL) model (Cercignani & Lampis, 1971). This approach is widely used for GSI modelling (Walker et al., 2014; Mehta et al., 2017). The Response Surface Model (RSM) developed by Mehta et al. (2017) was designed to implement the CLL model applied to the CHAMP and GRACE satellites. In the RSM model, the energy accommodation coefficient is assumed to be variable and detailed geometry of the satellites was used to calculate the drag coefficients. However, the inputs for such models are difficult to estimate without insitu observations of the atmospheric composition. Therefore, in this article, we prefer to infer information on the energy accommodation coefficient by changing its value and analyzing the effect on the selfconsistency of the density data.
In the first part of the results (Sect. 3), the GSI modelling influence on the aerodynamic coefficients for CHAMP, GRACE, GOCE, and Swarm are introduced, and variations with attack and sideslip angles are examined. The starting point of this aerodynamic analysis is based on the vector expression for the aerodynamic acceleration, which is defined as follows:
According to this definition the drag coefficient is defined as the contribution along the relative velocity direction:
C_{x}, C_{y}, and C_{z} are components of the aerodynamic acceleration (Eq. (5)) computed along the axes of the satellite reference frame, which in nominal attitude corresponds to the longitudinal (X: alongtrack), horizontal sideways (Y: crosstrack) and vertical downward (Z: antiradial) directions, respectively. The selected aerodynamic coefficients are modeled in the satellite reference frame and have been normalized with a reference area set to 1 m^{2} for all missions. This fixed reference area does not depend on the attack and sideslip angles, and therefore variations of the true exposed satellite area do not need to be independently calculated, as they are included in the normalized force coefficients. A representation of the analysed angles is available in Figure 1.
Fig. 1 Representation of sideslip and attack angles for GRACE. 
The thermospheric density is retrieved using the previous aerodynamic coefficients with the algorithm discussed in Doornbos (2011) and March et al. (2019a). The methodology is based on processing the accelerations derived from satellite observations and using the previously mentioned drag equation (Eq. (5)). The CHAMP, GRACE, and GOCE onboard accelerometers provided highprecision information. Unfortunately, for the Swarm satellites, many anomalies occurred in the accelerometer instruments (Siemes et al., 2016), and for this reason, a density product based on GPSderived accelerations was introduced as part of the official ESA documentation (van den IJssel et al., 2019). For CHAMP and GOCE the crosstrack component has been used to derive wind products. However, for GRACE and Swarm, winds are difficult to retrieve due to lack of sufficiently precise accelerometer measurements (due to platformrelated issues) and especially for the small aerodynamic acceleration signal typical of high altitudes (around and over 500 km) (March et al., 2019b).
3 Results
3.1 Gassurface interactions influence on satellite aerodynamics
In this analysis, the GSI influence is investigated and illustrated for six different values of the energy accommodation coefficient (0.0, 0.2, 0.4, 0.6, 0.8, and 1.0). As previously mentioned, the value of 1 represents reflections with full accommodation to the spacecraft wall temperature, while the value 0 is for collisions without thermal energy exchange. The data sets are obtained as a function of attack and sideslip angles and a range of speed ratio values.
Figures 2–4 show the aerodynamic coefficients for the GRACE, CHAMP, and Swarm satellites for a wide range of sideslip angles. The C_{x}, C_{y}, and C_{z} components are available in the figure for the same range of α_{E} values. During mission lifetime, the attack angle is mainly centered around the nominal flight configuration of 0° in attack angle, while the sideslip angle is usually less stable and can vary over the full domain from approximately 0–180° during manoeuvres. The three plots of Figures 2–4 show that the coefficients are lower in magnitude when the collisions are closer to the fully accommodated mode. Moreover, looking at the constant step of 0.2 in α_{E}, it is clear that between 0.8 and 1.0 difference in the aerodynamic forces is much larger than between 0.0 and 0.2. A description of the maximum computed differences varying the attack and sideslip angles are available in Tables 1 and 2, respectively. The computed values are based on the full attitude ranges simulated and available on the aerodynamic tables. The full aerodynamic data sets are added as complementary material of this article. The higher sensitivity for coefficients nearby the fully diffusive mode was already observed in the zonal wind analysis by March et al. (2019b). Varying α_{E}, the shape of the aerodynamic force coefficient curves remain the same without relevant differences. However, the main change, as expected, is a bias between the different computed values within the selected α_{E} range. When we inspect the differences between nonaccommodated and fully diffusive modes, the percentage difference reaches up to 84.5% for CHAMP, 84.0% for GRACE, 52.0% for GOCE and 82.1% for Swarm (Tables 1 and 2). For the study of the attack angle variation, the illustrated coefficients are C_{x} and C_{z}. Indeed, the crosstrack component of the aerodynamic force (C_{y}) is negligible. The plots for the sideslip angle variation show the C_{x}, C_{y} coefficients. The C_{z} coefficients are mostly due to the inclined side panels of the satellites. The symmetrical shape of GOCE brings to a negligible contribution over the vertical direction. For this reason, only the most relevant aerodynamic contributions are shown in the enclosed plots. However, for the quantitative analysis of Table 2, the results for the Zcomponent are also provided. In Figure 2, the nearly symmetrical shape of the lobes for the C_{x}, C_{y}, and C_{z} coefficients can be observed for the GRACE satellite. Whereas, for CHAMP and Swarm, Figures 3 and 4 show the asymmetric shape of the C_{x} and C_{z} coefficient lobes, which is a consequence of the presence of booms (boom pointing into flight direction for CHAMP and antiflight direction for Swarm). A different sensitivity to the energy accommodation coefficient can also be observed depending on the nominal or backward orientation of the satellite. When CHAMP and Swarm have their boom pointing into flight direction, the drag coefficients are less dependent on the accommodation coefficient value. If the satellite exposes the large side to the incoming flow, the collisions play a crucial role in the force coefficients determination.
Fig. 2 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the GRACE satellite. 
Fig. 3 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the CHAMP satellite. 
Fig. 4 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the Swarm satellites. 
Maximum force coefficient differences for different α_{E} ranges along X and Z satellite reference frame axes for the simulated attack angles.
Maximum force coefficient differences for different α_{E} ranges along X, Y, and Z satellite reference frame axes for the simulated sideslip angles.
Appendix A shows the force coefficients for the CHAMP, GRACE, GOCE, and Swarm satellites for different attack and sideslip angles. For GOCE, large attitude manoeuvres were not performed. Therefore, the sideslip angle range is smaller (between −16° and 16°). All the figures from this section (including the previous polar plots) are obtained for a speed ratio of 7. A complete analysis showing the different coefficients for the analyzed range of speed ratios is provided in Appendix B. In this Appendix, the C_{x}, C_{y}, and C_{z} coefficients are shown in the 0.5–14 speed ratio range and for different sideslip angle configurations (i.e., 0°, 45°, and 90°).
As a further investigation, it is interesting to study how the latitude variations influence satellite aerodynamics. In Figure 5, this is shown for the CHAMP drag coefficient (C_{D}) for three different solar activity levels. The shaded areas represent the drag coefficient variability over the full day estimated over a resolution of 2° in the argument of latitude. With the argument of latitude, we refer to the angle along the orbit starting at the ascending node. In Figure 6, the same is produced for GRACE. The three days were already selected in March et al. (2019a). High, moderate and low activities correspond to 20021027, 20050515, and 20090828, respectively. For CHAMP, it is possible to notice that peaks in drag are reached at the equator for 0° and 180° in the argument of latitude. This is especially clear for high and low activity cases. Relevant differences in behaviour among different values of α_{E} are not present. Indeed, looking at the full day of observations, between different energy accommodation coefficients a mean bias can be identified as main effect on the satellite aerodynamics. For GRACE, the daily variations are limited within the selected days, this is most likely due to the simpler geometry and the less erratic reflections over satellite surfaces. However, in the proximity of the south pole, a larger variability can be noticed. This effect is especially visible for the low activity period, and it is probably connected to the more uncertain chemistry and wind inputs from NRLMSISE00 within such activity levels. For GRACE is also illustrated the drag coefficient evolution for the energy accommodation coefficient of 0.85 in order to intercompare similar conditions among the two missions. In both figures, the drag coefficient assumes an unitary reference area, therefore C_{D} already contains information about the satellite geometry.
Fig. 5 CHAMP aerodynamic drag coefficients as a function of the argument of latitude in three selected days describing high, moderate and low solar activity. 
Fig. 6 GRACE aerodynamic drag coefficients as a function of the argument of latitude in three selected days describing high, moderate and low solar activity. 
3.2 Gassurface interactions influence on neutral density
In this article, the optimal α_{E} values from March et al. (2019b) are implemented in the density processing to generate the newly derived data sets. When comparing the new and the previously adopted α_{E} value of 1.00 from March et al. (2019a), it is observed that the average difference between new and previous densities is around 6% for CHAMP and 11% for GOCE when looking at longterm averages. The differences with semiempirical models are influenced by solar and geomagnetic activity because models perform differently for different geomagnetic and solar activity inputs (Emmert, 2015). The semiempirical models are commonly closer to the observed density during periods of high solar activities and within this condition, the agreement with the new results is higher. However, when these models are introduced for comparisons, the results provide qualitative information that must be carefully interpreted. This is especially connected to the uncertainties introduced by the atmospheric models. The presented results aim to help reduce these uncertainties and improve current empirical or physicsbased models.
Analyzing long time periods, it is possible to investigate the sensitivity of the energy accommodation coefficient on the new density data with respect to the fixed output of a semiempirical model. In Figure 7, the density ratio with respect to the NRLMSISE00 model is shown for three different solar activity levels. The high activity is represented by F_{10.7} values between 130 and 300 solar flux units (sfu). The medium activity is in the range between 90 and 130 sfu, while the low solar activity is chosen for values of F_{10.7} under 90 sfu. The selected time periods are between January 2003 and January 2006 for CHAMP and GRACE, January 2010 and January 2013 for GOCE, and July 2014 and July 2017 for Swarm, respectively. Comparing the three subsets, it is possible to notice a lower agreement between different satellites for low values of F_{10.7}. Moreover, the average ratios move from a range of 0.73–0.99 to 0.61–0.91 when comparing the high and low F_{10.7} sets, respectively. From one side, this is due to the lower drag. Indeed, under such conditions, a larger influence of the error in solar radiation pressure modelling plays a crucial role in density estimation (van den IJssel et al., 2019). However, the lower performance of the semiempirical models during deep–low solar activity is also affecting the comparisons (Emmert, 2015). For the results for high activity, a clear optimal value of the accommodation coefficient which guarantees a higher consistency or lower variability among different missions cannot be identified. The large uncertainties introduced by NRLMSISE00 model may be the cause of such difficulties. However, out of this analysis, it is important to identify the evolution of the density ratios over different values of accommodation, quantifying the previously mentioned steeper variation nearby scenarios with accommodation coefficients close to one.
Fig. 7 Long period density ratios with respect to the NRLMSISE00 model versus the energy accommodation coefficient for different ranges of F_{10.7}. 
From this analysis, it is not straightforward to retrieve an optimal α_{E} coefficient. However, exploiting the particular manoeuvres of the Swarm constellation, as outlined in the next section, it is possible to further assess the neutral density and provide additional information on the energy accommodation coefficient determination.
Using a constant energy accommodation coefficient creates a constraint in specific features like differences in chemical composition during day–night transitions and along the orbit. To quantify the difference with other models like the CLL, a comparison with other densities is performed. For this study CHAMP was introduced because of its long mission lifetime with stable accelerometer performances and its representative altitudes among the selected satellites. In Figures 8–10, a comparison between different models is provided for high, moderate, and low solar activity, respectively. The new densities from this work are compared with the ones obtained by Mehta et al. (2017) and Doornbos (2011). Additionally, the density for CHAMP is also estimated using the SESAM empiricalmodel (Pilinski et al., 2013b), which in this case adopts the aerodynamic information generated with the highfidelity geometry from this work. Mehta et al. (2017) adopt an approach based on the CLL model, while Doornbos (2011) uses the DRIA method applied to a simplified macromodel. Only the first three hours for each selected day are shown to fully appreciate the variations within around two orbits. Comparing the new density with the ones derived by Mehta et al. (2017), larger variations can be noticed. The recorded discrepancies have a periodic behaviour associated with the latitude and orbit period. The greater peaks are over the poles and reach a maximum difference with respect to the new results of around 10–15%. Higher fluctuations in densities and their ratios are also localized nearby polar regions, this is mostly due to the complex dynamics in those areas. In particular, it is possible to find larger fluctuations in the horizontal winds, which are not included in the other processing schemes regarding the DRIA approaches. For completeness, in Figures 8–10, also the horizontal wind components estimated with the HWM07 model (along satellite’s longitudinal and cross directions) are introduced. Localized wiggles in the density ratios can be identified where these components are largely varying in magnitude. A more stable behaviour can be shown for the densities by Doornbos (2011), characterized by lowfrequency variations with respect to the new density. This is associated with the use of the same chemistry inputs from the NRLMSISE00 model and accelerometer calibration parameters from the same processing scheme. The additional output provided by the SESAM model covers the gap in knowledge about the differences between variable accommodation coefficients and the DRIA models. It is important to note that the SESAM density ratio highlights exclusively the difference in estimating α_{E} because the same accelerations and geometry are used to process the illustrated data. Analysing the variations between the two variable accommodation coefficients models, it is possible to identify a similar behaviour. However, a few differences can be found as well. Indeed, the SESAM model has an overall higher density ratio with respect to the RSM model by Mehta et al. (2017). The oscillations are similar, however, they do not seem to be solved by a pure bias, and therefore, there might be a chemical composition effect beyond the different geometry influence. This can be noticed by looking at the variations over the satellite argument of latitude. The largest differences appear during the low solar activity period (Fig. 10). Under this activity level, a good agreement among the variable accommodation models can be identified at the argument of latitudes nearby 180°. While the maximum discrepancy between the two CLL models can be found at 270° (South Pole). On the other side, at 270° it is possible to find the closest match with the DRIA densities. This needs a further investigation focused on the chemistry and solar activity inputs, chemical composition of nearby satellite surfaces, and the different collisions physics applied by the introduced models. A deeper investigation of the F_{10.7} input could enhance the quality of the data, however this remains out of the scope of this work.
Fig. 8 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under high solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 
Fig. 9 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under moderate solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 
Fig. 10 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under low solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 
The presented comparisons provide the novel opportunity to quantify the discrepancies between current CLL models, beyond the CLL and DRIA models. Among the presented data, density variations are within ±10%, however peaks can reach up to 15% over poles. The cyclic variations are expected with be associated to the different chemical constituents adsorbed over satellite surfaces along the orbit, which can be different also using similar CLL models with different inputs. Among the different features, it is interesting to note how the α_{E} varies with respect to the fixed value used in this work (0.85) and the one from Doornbos (2011) (0.93). Starting from the high activity (Fig. 8), it is possible to see that SESAM provides higher values for α_{E}. The minimum values over the orbit are close to the 0.93 value, and they seem to be localized over the poles. For moderate activity (Fig. 9), the accommodation coefficient used in Doornbos (2011) seems to agree with the average value from SESAM over the orbit. While for the low activity level (Fig. 10), SESAM’s α_{E} oscillates between the two fixed values of 0.85 and 0.93. In agreement with these results, interesting information about the α_{E} variability over the orbit can be retrieved and implemented in an updated version of the models. A possible way forward to process new density data could be based on combining the values retrieved with orbit manoeuvres with the periodical variation introduced in models such as SESAM. Possible studies over this opportunity are currently ongoing as part of the Committee on Space Research (COSPAR) International Space Weather Action Teams (ISWAT) activities (https://iswatcospar.org/) which involve part of the authors and will be further explained in a dedicated article.
3.3 Density consistency for Swarm A and C
Between the 13th and 14th of May 2014, the Swarm A and C satellites performed multiple attitude manoeuvres. The main objective of the four 90° yaw slew rotations was focused on investigating differences between the measurements of the vector and absolute magnetometers. However, it is possible to exploit these changes of orientation with respect to incoming atmospheric flow to retrieve information on the satellite aerodynamics.
Looking at a single satellite at a time, the thermospheric variability affects the density comparison in time. However, since Swarm A and C were flying next to each other, the simultaneous measurements of the two satellites can be compared. It is fair to assume that the extracted densities would be nearly identical. This assumption is especially justified during high solar activity, when the acceleration signal magnitude is high. Indeed, during high activity, the average ratio between the two densities is consistent to within 1%, and these small differences in density generally agree with the expected diurnal density gradient. For the current very low solar activity data, this density ratio can reach values up to 50%. If the experiment with the attitude manoeuvre would be repeated, it would therefore be necessary to wait for high solar activity conditions and/or significantly lower orbit altitudes.
In the two analyzed days in 2014, the average F_{10.7} value was around 165 sfu and represents a scenario of high solar activity. The absence of geomagnetic activity events provided a suitable output without undesired spikes. In Figure 11, the performed manoeuvre is represented by the absolute value of the yaw angles for both satellites. The coloured shaded areas highlight the orbits that are taken into account for the analysis. All orbits containing large slew rates have been discarded because of high thruster activity and uncertainty induced by the Kalmanfilter approach in the GPSderived accelerations. In the bottom part of Figure 11, the density ratio between Swarm A and C measurements is illustrated as a function of the energy accommodation coefficient. Each line corresponds to one of the selected orbits following the same color specified in the top part. Varying the accommodation coefficients, it is possible to notice a variation in the slope of the density ratios between Swarm A and C. With the green markers, the intersections between different periods are highlighted and appear to be concentrated around the value of 0.85.
Fig. 11 Yaw angle variations and density ratios between Swarm A and C measurements during the May 2014 manoeuvre. Selected orbit numbers are available in the top and bottom parts of the figure. 
Particular attention needs to be focused on the orbits, which, at the same time, are characterized by different satellite orientations with respect to the flow. Two of these orbits are numbered 11 and 17. For these two specific cases, the satellites are respectively in backward and sideways orientation with respect to the nominal flight. Having the capability to estimate the densities within the range of 0–1 for α_{E}, setting a different accommodation coefficient for each satellite, it is possible to create a map of density ratios within the two orbits. The two maps for the two selected orbits are available in Figure 12. In this representation of the α_{E} values for Swarm A and C, an area can be identified where the density ratio between the two satellites is close to one (in white). Highlighting the optimal ratios with dashed lines, the intersections for both orbits with the diagonal of the map coincide with an α_{E} value of 0.85. Similar conclusions can be drawn looking at the orbits with numbers 5 and 23. In this case, the two satellites are in the nominal and sideways configurations. For these two periods, an intersection is achieved for an α_{E} value of 0.81. Similar results can be found including all the remaining selected orbits, and all the retrieved values are nearby 0.85 and below the currently adopted α_{E} value of 0.93 (Fig. 13).
Fig. 12 Swarm A and C density ratios for periods 11 and 17 (on the left and right side, respectively). The optimal density ratios for which satellites are in better agreement are highlighted with dashed lines. 
Fig. 13 Swarm A and C optimal density ratios for May 2014 manoeuvres and selected periods. 
The results coincide with the values obtained from the accelerometerderived wind analysis of the CHAMP and GOCE satellites in March et al. (2019b), which showed optimal α_{E} of around 0.85 and 0.82 for the two missions respectively. Looking at the overall intersections available in Figure 11 and the other results from this section, a value of 0.85 is recommended for future Swarm density processing. This is especially recommended for the processing of Swarm A and C data. Comparing the current densities, processed with an α_{E} value of 0.93 (available at https://swarmdiss.eo.esa.int/, August 2021), the new density data would be slightly lower in magnitude with a 4–5% difference. Beyond the assumption to have a fixed energy accommodation coefficient, which likely does not reflect the collision physics in proximity of satellite surfaces, the results from the Swarm manoeuvre analysis further validate the outcome that, if the DRIA model is adopted, for the density (and wind) estimation, a lower value of α_{E} needs to be introduced in the calculations. In the next section, this will also be verified using longer time periods and intercomparing the new densities with other sources and methods.
4 Comparisons with atmospheric and GSI models
In order to evaluate the new density data, further analysis is performed using the results from the NRLMSISE00 and DTM2013 semiempirical models. The choice of these two semiempirical models was based on their wide use among scientific users and their applications to past and current TU Delft and ESA projects. In particular, NRLMSISE00 is one of the most widely used models, and DTM2013 is one of the latest semiempirical models. A further comparison with the WACCMX general circulation model (Liu et al., 2018) is also performed using a set of representative months. The WACCMX model is a physicsbased model developed at the High Altitude Observatory (HAO) of the National Center for Atmospheric Research (NCAR) (Liu et al., 2018). In order to have an overview of the newly generated neutral densities, the outputs are compared for specific conditions. The simulated scenarios represent three different periods with high, medium, and low solar activity. One month of data for each subset is investigated and shown in Figure 14. The months were selected according to acceleration data quality and availability and range of representative solar flux and geomagnetic inputs excluding months with large geomagnetic storms. The new results are normalized with the NRLMSISE00 model at the altitude of 400 km for the CHAMP and Swarm satellites. This procedure was already used with the DTM2000 model in Bruinsma et al. (2006). For GOCE, the normalization is performed at 250 km in order to provide a more representative altitude. This normalization was especially necessary for the Swarm satellites because the WACCMX upper pressure level altitude boundary extends to around 500 km altitude. The monthly time window was selected to limit the computational cost of the physicsbased model simulations.
Fig. 14 Comparison of the density ratios between the newly derived density and the NRLMSISE00 and WACCMX models results. The densities are normalized at 400 km altitude except for GOCE, which is normalized at 250 km. 
For the GRACE satellites, the relatively large density error due to the high altitude of this mission impedes the investigation of an optimal accommodation coefficient. An α_{E} value of 0.85 was adopted for this comparison with WACCMX and this needs to be taken into account when looking at the results. Additional studies would require a deeper investigation after reprocessing current GRACE data. This would be especially useful for the second phase of the mission, which was characterized by a dramatic degradation of accelerometer measurements (Klinger & MayerGürr, 2016).
As shown in Figure 14, which shows the Probability Density Function (PDF) of the density ratios between the new and the modelled densities, both models are currently providing higher estimates of the neutral density than the satellite data, especially under high and moderate solar activities. For low solar activity periods, the differences in the ratios between NRLMSISE00 and WACCMX compared with the new data are much larger. These comparisons are difficult to interpret because of the large errors, which allows a qualitative representation. Generally, the new data are in better agreement with WACCMX than with NRLMSISE00. For high activity, this is clearly visible for all the selected satellites. For GOCE, the density ratios show high agreement for both NRLMSISE00 and the physicsbased model, however for both, a bias of about 30–40% can be noticed. A similar behavior between the two introduced models and the new densities is verified for the moderate solar activity levels. This is valid for all the selected satellites. For Swarm, the newly derived density is in good agreement with the WACCMX model (similarly to CHAMP). However, if we look at the low activity periods, very different values can be detected for both models. This is especially the case for Swarm B, which is highly affected by the low signal magnitude in deep–low solar activity. Imperfections in solar radiation pressure modelling have a strong contribution. An improvement of this additional contribution would enhance the quality of derived neutral densities, which is currently under investigation.
In order to provide an overview over longer periods, the average density ratios with respect to the NRLMSISE00 and DTM2013 models have been computed (Fig. 15). The illustrated periods are between 20030101 and 20050101 for CHAMP and GRACE, 20110101 and 20130101 for GOCE, and 20140719 and 20170719 for the Swarm satellites. All these periods guarantee high solar activity conditions among the selected missions, which cover different periods. This enhances the accuracy of the analysis because of the higher aerodynamic signal level (e.g., lower impact of solar radiation pressure errors) and better atmospheric model performance during such high solar activity conditions.
Fig. 15 Comparison between the mean density ratios (μ*) between the NRLMSISE00 model (top) and DTM2013 (below) with different data sets for the CHAMP, GRACE, GOCE, and Swarm satellites. The percentage values indicate the maximum variability range among the satellites. Note that for GRACE and Swarm B the new density ratios are computed using a preliminary energy accommodation coefficient of 0.85 (which is not validated as for the other satellites). 
The newly derived densities (labeled as “This work”) are compared with the previous results from March et al. (2019a) for fully accommodated reflections (α_{E} = 1.00). The original densities computed with macromodel geometries and 0.93 as accommodation coefficient are available in the first column on the left side of the diagram in Figure 15. The densities estimated with α_{E} = 0.93 with new highfidelity geometries are included as well. In addition to the DRIA models, the estimated densities for CHAMP and GRACE are introduced from Mehta et al. (2017) and using the SESAM model. Currently, such information for the CLL models is only available for these two missions. Among the new densities, the bars for GRACE and Swarm B are shaded and with red dashed lines to highlight the preliminary assumption adopted in the illustration, which uses a 0.85 accommodation coefficient. Such value is not expected to reflect real conditions, however it is left in the figure to provide a comparison over similar GSI conditions. To identify a better value for the DRIA, it would be necessary to find other ways to mitigate the influence of low signal strength and quality for such satellites at high altitudes. This is out of the scope of the presented research, but it is recommended as future work to be done after an enhancement of accelerometer data quality (e.g., improvement of calibration, thruster activities filtering) and solar radiation pressure modelling.
Figure 15 illustrates the consistency of the density observation data sets with two representative thermosphere models. On the top, there are the results from the NRLMSISE00 model, and at the bottom the ones for DTM2013. The figure shows the mean density ratio (μ*) of the model and the observations as an index of their consistency. Under each group of bars, we provide the difference in percentage (Δ) between the maximum and minimum mean densities among the selected satellites. This is assumed to be a measure of selfconsistency for multimission observations. Figure 15 shows that the new densities are lower than those presented in March et al. (2019a). However, even if the agreement with the semiempirical models is not improving, the results are more consistent among the selected missions. Indeed, the scale with respect to the NRLMSISE00 model is now more constant across the missions. Since our primary objective is a high selfconsistency of the multimission density observations, we focus in the discussion of Figure 15 on the differences between the Δ values. The new α_{E} bars from this study are within 4.13% of the variation among the satellites, whereas for the previous version of the data α_{E} = 1.00 this value reached 8.26%. Introducing the α_{E} = 0.93 value, with the old panel method differences can reach 23.9%. For the same accommodation coefficient, the variability decreases to 5.67% when the new geometries and the SPARTA modelling are adopted. Looking at the differences between Swarm A and Swarm B, a lower accommodation coefficient for Swarm B, which flies at a higher altitude, might introduce a better match with the Swarm A density ratio.
Introducing a variable accommodation coefficient in the density processing, as in Mehta et al. (2017), gives differences around 11% for both density models. While, according to SESAM, these values are reduced to 2.48% and 3.34% for NRLMSISE00 and DTM, respectively. A clear difference between the two variable accommodation coefficient models can be found regarding the GRACE density. According to SESAM, it seems that the density ratio for GRACE is not dropping as for the data set from Mehta et al. (2017). From Figure 15, it seems that when the energy accommodation coefficient decreases in magnitude (for GRACE due to the higher altitude), SESAM provides a higher density with respect to the other CLL model. This could be already seen in Figure 10, where for low solar activity and in correspondence with the low peaks in α_{E} for SESAM, a similar behaviour could be captured for CHAMP. The different used geometries could not induce such a cyclic variation in the density ratio. For this reason, this feature could probably be related to different chemistry and adsorption assumptions introduced by the two models. Comparing the SESAM output with the new densities from this work, it is also intriguing to note that the two results are comparable. This is verified even though the GSI applied methods are very different between the two approaches. The low percentages achieved within this work need a careful interpretation due to the uncertainties associated with the atmospheric models. However, the high degree of selfconsistency between the selected missions and both empirical models is a promising result and certainly a starting point for further improvements.
5 Conclusions
The presented research investigates the influence of GSI on aerodynamic modelling and density data sets. New satellite aerodynamic models are made available to the scientific community as part of the Supplemental Material and on http://thermosphere.tudelft.nl (August 2021), including the used geometry models. The possibility to further validate this work with dedicated test campaigns is facilitated through the new aerodynamic coefficients described in Section 3.1. The figures and analyses of this section and Annexes Appendices A and B aim to provide insights for future aerodynamic studies and possible experimental campaigns investigating satellite aerodynamics. New neutral density data sets from the CHAMP, GRACE, GOCE, and Swarm satellites are obtained and investigated in view of the energy accommodation coefficient. This parameter plays a crucial role in GSI modelling. The previously adopted value of 0.93, which was used at TU Delft in the processing of the official ESA density data sets for GOCE and Swarm, is higher than the value of 0.85 that we found here by analyzing attitude manoeuvres of Swarm A and C. The value is also higher than the values of 0.82 and 0.85 for GOCE and CHAMP, respectively, which were derived from the study of the neutral winds presented in March et al. (2019b).
All of our comparisons between thermosphere models and observations show that the models have overestimated the mean density. Underestimating density means that inorbit lifetimes of space objects are longer, and this has adverse effects on the low orbit space debris population. Very recently, newly released or inprogress updates to thermosphere models (Emmert et al., 2020; Jackson et al., 2020) have started to reflect the lower thermospheric densities, although none of the models so far has incorporated a parametrization of longterm change effects.
The continued use of models based on older data in the accelerationtodensity data processing, knowing that they are biased, limits quantitative analyses of the energy accommodation coefficient. So far, such analyses are limited to studying attitude manoeuvres and exploiting synergies between satellites. This is demonstrated through the analysis of the Swarm A and C manoeuvres in May 2014. The two satellites provided the unique opportunity to compare simultaneous measurements at high solar activity with a large aerodynamic signal magnitude. Through investigating the density ratios, the optimal value of the energy accommodation coefficient for the Swarm satellites during the manoeuvre is found to be in the range of 0.80–0.90, indicating that the previously used value of 0.93 was too high. The newly found range is however in agreement with the previous analysis of thermosphere wind observed by the CHAMP and GOCE satellites. We find that applying an average value of 0.85 for the energy accommodation coefficient in future data processing will improve the Swarm density data set’s consistency, both between the three satellites, and with the other missions. Performing further attitude manoeuvres in future phases of the Swarm mission will enable the investigation of the dependence of the energy accommodation coefficient on, e.g., the altitude and solar activity level, which are both expected to affect the energy accommodation through the atomic oxygen concentration and temperature. The presented new Swarm densities are lower in magnitude than the ones obtained with α_{E} = 0.93 and are presently available on the beforementioned website as well as in ESA’s Swarm data archive. The difference in magnitude is expected to be around 5%. In the future, further exploitation of the Swarm A and C synergy is strongly encouraged. However, at the current (Mid 2021) low solar activity in combination with the high altitude, the aerodynamic signal magnitude is still too low to expect significant new insights.
The presented research would strongly benefit from an improved solar radiation pressure modelling, especially for the GRACE and Swarm satellites, which have spent a significant portion of their lifetime at relatively high altitude during solar minimum. Our efforts to augment the new highfidelity satellite geometries with surface properties for improving the solar radiation pressure models are currently ongoing.
We also see potential to improve the processing of GRACE accelerometer data, which would aid, among others GSI investigations as presented in this paper. In particular, the accelerometer data calibration for the last seven years of the mission has a great potential for improvement. Indeed, from April 2011 onward, the thermal control of the accelerometer was deactivated to save battery life. This resulted in significant perturbations related to the fluctuating instrument temperature.
The consequences of our assumption that the energy accommodation coefficient is constant for all solar and geomagnetic conditions need to be further investigated as well. This deserves a deeper investigation to characterize the impact of the chemical composition and temperature of the thermosphere, which significantly varies during the solar cycle and between the day and nightside, on GSI parameters and, thus, density observations. An improved level of accuracy of thermosphere data is expected when using GSI models based on variable accommodation coefficients with valid temperature and composition inputs. However, accurate new insitu measurements of the chemical composition are needed to reliably estimate the parameters of such models. In particular, additional data on the number density of nitrogen, helium, and atomic oxygen is crucial in this context. The need for new thermospheric temperature and composition data was highlighted by Emmert et al. (2020) as well. Dedicated satellite missions are strongly recommended to reduce current uncertainties. As shown in this research, GSI can presently only be investigated in an indirect and imperfect way, making use of scarce data from satellite constellations, manoeuvres and seasonal analysis. Nevertheless, some steps forward can be made with such analyses to improve the selfconsistency of the thermosphere data sets. More data of these types will certainly be helpful for further investigations. However, new experiments with more extensive instrumentation will be necessary to resolve open issues.
In the future, dedicated satellite thermosphere density and satellite drag experiments will need to measure not only the accelerations (using accelerometers, GNSS receivers, and star cameras as on the satellites analysed here) but also to independently and accurately measure the temperature, composition, wind and density on the same platform as input to the satellite aerodynamic model. Ideally, such experiments would eventually cover all possible temperature and composition environments, by spanning both solar minimum and maximum conditions, as well as a wide range of altitudes, including the nearly unexplored region below 200 km. The insights gained from such experiments could be retroactively applied to the data of the CHAMP, GRACE, GOCE, and Swarm missions.
Acknowledgments
This work was supported by the Netherlands Organisation for Scientific Research (NWO) project number ALWGO/1435. The authors would like to thank the European Space Agency for supporting the production of GOCE (Contract 18308/04/NL/MM) and Swarm (Contract 4000102140/10/NL/JA) thermospheric products. An additional acknowledgement goes to the High Altitude Observatory (HAO) from the National Center for Atmospheric Research (NCAR) and Dr. F. Gasperini for the possibility to work with the WACCMX model for part of this project. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.
Appendix A
New aerodynamics data sets
Fig. A.1 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the CHAMP satellite. 
Fig. A.2 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the CHAMP satellite. 
Fig. A.3 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the GRACE satellite. 
Fig. A.4 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the GRACE satellite. 
Fig. A.5 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the GOCE satellite. 
Fig. A.6 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the GOCE satellite. 
Fig. A.7 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the Swarm satellites. 
Fig. A.8 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the Swarm satellites. 
Appendix B
Speed ratio influence over satellite aerodynamics
The following figures show that an attitude variation affects the sensitivity of the aerodynamic forces to the speed ratio and energy accommodation coefficient. For low values of speed ratio (i.e., below 4) the force coefficients reach high magnitudes. Moreover, the force coefficients converge for very low speed ratios (s < 2), i.e. the drag coefficient is then less sensitive to the energy accommodation coefficient. The drag coefficient approaches a constant value towards higher speed ratios. However, due to the different collisions the convergence can vary depending on the satellite attitude. For GOCE, the Zcomponent is not included because of the symmetric shape of the satellite with respect to the X–Y plane.
Fig. B.1 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the CHAMP satellite for fixed sideslip angles. 
Fig. B.2 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the GRACE satellite for fixed sideslip angles. 
Fig. B.3 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the GOCE satellite for fixed sideslip angles. 
Fig. B.4 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the Swarm satellites for fixed sideslip angles. 
References
 Bowman BR, Kendra MJ, Marcos FA. 2004. A method for computing accurate daily atmospheric density values from satellite drag data. In: 14th AAS/AIAA Space Flight Mechanics Conference, Maui, Hawaii, February 8–12, 2004. [Google Scholar]
 Bruinsma S. 2015. The DTM2013 thermosphere model. J Space Weather Space Clim 5: A1. https://doi.org/10.1051/swsc/2015001. [CrossRef] [EDP Sciences] [Google Scholar]
 Bruinsma S, Biancale R. 2003. Total densities derived from accelerometer data. J Spacecr Rockets 40: 230–236. https://doi.org/10.2514/2.3937. [CrossRef] [Google Scholar]
 Bruinsma S, Forbes JM, Nerem RS, Zhang X. 2006. Thermosphere density response to the 20–21 November 2003 solar and geomagnetic storm from CHAMP and GRACE accelerometer data. J Geophys Res 111: A06303. https://doi.org/10.1029/2005JA011284. [Google Scholar]
 Calabia A, Jin S. 2016. New modes and mechanisms of thermospheric mass density variations from GRACE accelerometers. J Geophys Res: Space Phys 121(11): 11191–11212. https://doi.org/10.1002/2016JA022594. [Google Scholar]
 Cercignani C, Lampis M. 1971. Kinetic models for gassurface interactions. Transp Theory Stat Phys 1(2): 101–114. https://doi.org/10.1080/00411457108231440. [CrossRef] [Google Scholar]
 Doornbos E. 2011. Thermospheric density and wind determination from satellite dynamics. Ph.D. thesis, Delft University of Technology, Delft, The Netherlands. ISBN: 9789090260518. [Google Scholar]
 Doornbos E. 2012. Producing density and crosswind data from satellite dynamics observations. In: Thermospheric density and wind determination from satellite dynamics, Springer Berlin Heidelberg, Berlin, Heidelberg: 91–126. https://doi.org/10.1007/9783642251290. [CrossRef] [Google Scholar]
 Emmert J. 2015. Thermospheric mass density: A review. Adv Space Res 56: 773–824. https://doi.org/10.1016/j.asr.2015.05.038. [CrossRef] [Google Scholar]
 Emmert J, Picone JM, Meier RR. 2008. Thermospheric global average density trends, 1967–2007, derived from orbits of 5000 nearEarth objects. Geophys Res Lett 35(5): L05,101. https://doi.org/10.1029/2007GL032809. [CrossRef] [Google Scholar]
 Emmert JT, Drob DP, Picone JM, Siskind DE, Jones M, et al. 2020. NRLMSIS 2.0: A wholeatmosphere empirical model of temperature and neutral species densities. Earth Space Sci 8(3): e2020EA001321. https://doi.org/10.1029/2020EA001321. [Google Scholar]
 Gallis MA, Torczynski JR, Plimpton SJ, Rader DJ, Koehler T. 2014. Direct simulation Monte Carlo: The quest for speed. AIP Conf Proc 1628: 27–36. https://doi.org/10.1063/1.4902571. [CrossRef] [Google Scholar]
 Gregory JC, Peters PN. 1987. A measurement of the angular distribution of 5 eV atomic oxygen scattered off a solid surface in earth orbit. In: Proceedings of the 15th international symposium on rarefied gas dynamics, Vol. 1, Boffi V, Cercignani C (Eds.), B.G. Teubner, Stuttgart, pp. 644–656. [Google Scholar]
 Hedin AE, Hinton BB, Schmitt GA. 1973. Role of gassurface interactions in the reduction of OGO 6 neutral particle mass spectrometer data. J Geophys Res 78(22): 4651–4668. https://doi.org/10.1029/JA078i022p04651. [CrossRef] [Google Scholar]
 Jackson DR, Bruinsma S, Negrin S, Stolle C, Budd CJ, et al. 2020. The Space Weather Atmosphere Models and Indices (SWAMI) project: Overview and first results. J Space Weather Space Clim 10: 18. https://doi.org/10.1051/swsc/2020019. [CrossRef] [EDP Sciences] [Google Scholar]
 Klinger B, MayerGürr T. 2016. The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSGGrace2016. Adv Space Res 58(9): 1597–1609. https://doi.org/10.1016/j.asr.2016.08.007. [CrossRef] [Google Scholar]
 Lewis HG, Saunders A, Swinerd G, Newland RJ. 2011. Effect of thermospheric contraction on remediation of the nearEarth space debris environment. J Geophys Res: Space Phys 116(A2): n/a–n/a. https://doi.org/10.1029/2011JA016482. [CrossRef] [Google Scholar]
 Liu H, Bardeen CG, Foster BT, Lauritzen P, Liu J, et al. 2018. Development and validation of the Whole Atmosphere Community Climate Model With Thermosphere and Ionosphere Extension (WACCMX 2.0). J Adv Model Earth Syst 10(2): 381–402. https://doi.org/10.1002/2017MS001232. [CrossRef] [Google Scholar]
 March G, Doornbos E, Visser P. 2019a. Highfidelity geometry models for improving the consistency of CHAMP, GRACE, GOCE and Swarm thermospheric density data sets. Adv Space Res 63(1): 213–238. https://doi.org/10.1016/j.asr.2018.07.009. [CrossRef] [Google Scholar]
 March G, Visser T, Visser P, Doornbos E. 2019b. CHAMP and GOCE thermospheric wind characterization with improved gassurface interactions modelling. Adv Space Res 64(6): 1225–1242. https://doi.org/10.1016/j.asr.2019.06.023. [CrossRef] [Google Scholar]
 Matsuo T, Fedrizzi M, FullerRowell TJ, Codrescu MV. 2012. Data assimilation of thermospheric mass density. Space Weather 10(5). https://doi.org/10.1029/2012SW000773. [Google Scholar]
 Mehta P, Linares R, Sutton E. 2019. Datadriven inference of thermosphere composition during solar minimum conditions. Space Weather 17(9): 1364–1379. https://doi.org/10.1029/2019SW002264. [CrossRef] [Google Scholar]
 Mehta PM, Linares R. 2018. A new transformative framework for data assimilation and calibration of physical ionospherethermosphere models. Space Weather 34(27): 9. https://doi.org/10.1029/2018SW001875. [Google Scholar]
 Mehta PM, McLaughlin CA, Sutton EK. 2013. Drag coefficient modeling for grace using Direct Simulation Monte Carlo. Adv Space Res 52(12): 2035–2051. https://doi.org/10.1016/j.asr.2013.08.033. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Moe MM, Wallace SD, Moe K. 1993. Refinements in determining satellite drag coefficients: Method for resolving density discrepancies. J Guid Cont Dyn 16(3): 441–445. https://doi.org/10.2514/3.21029. [CrossRef] [Google Scholar]
 Picone JM, Hedin AE, Drob DP, Aikin AC. 2001. NRLMSISE–00 empirical model of the atmosphere: Statistical comparisons and scientific issues. J Geophys Res: Space Phys 107(A12): SIA 15–1–SIA 15–16. https://doi.org/10.1029/2002JA009430. [Google Scholar]
 Pilinski MD. 2011. Dynamic gassurface interaction modeling for satellite aerodynamic computations. Ph.D. thesis, University of Colorado at Boulder, Boulder, CO. https://scholar.colorado.edu/asen_gradetds/37/. [Google Scholar]
 Pilinski MD, Argrow BM, Palo SE. 2013a. Semiempirical model for satellite energyaccommodation coefficients. J Spacecr Rockets 47(6): 951–956. https://doi.org/10.2514/1.49330. [Google Scholar]
 Pilinski MD, Argrow BM, Palo SE, Bowman BR. 2013b. Semiempirical satellite accommodation model for spherical and randomly tumbling objects. J Spacecr Rockets 50(3): 556–571. https://doi.org/10.2514/1.A32348. [CrossRef] [Google Scholar]
 Pilinski MD, Bowman BA, Palo SE, Forbes JM, Davis BL, Moore RG, Koehler C, Sanders B. 2016. Comparative analysis of satellite aerodynamics and its application to spaceobject identification. Adv Space Res J 53(5): 1–11. https://doi.org/10.2514/1.A33482. [Google Scholar]
 Qian L, Marsh D, Merkel A, Solomon SC, Roble RG. 2013. Effect of trends of middle atmosphere gases on the mesosphere and thermosphere. J Geophys Res: Space Phys 118(6): 3846–3855. https://doi.org/10.1002/jgra.50354. [CrossRef] [Google Scholar]
 Siemes C, de Teixeira da Encarnação J, Doornbos E, Van den IJssel J, Kraus J, Pereštý R, Grunwaldt L, Apelbaum G, Flury J, Holmdahl Olsen PE. 2016. Swarm accelerometer data processing from raw accelerations to thermospheric neutral densities. Earth Planets Space 68(1): 92. https://doi.org/10.1186/s4062301604745. [CrossRef] [Google Scholar]
 Storz MF, Bowman BR, Branson MJI, Casali SJ, Tobiska WK. 2005. High accuracy satellite drag model (HASDM). Adv Space Res 36(12): 2497–2505. Space Weather, https://doi.org/10.1016/j.asr.2004.02.020. [Google Scholar]
 Sutton EK. 2008. Effects of solar disturbances on the thermosphere densities and winds from CHAMP and GRACE satellite accelerometer data. Ph.D. thesis, University of Colorado at Boulder, Boulder, CO, ISBN: 1243501634. [Google Scholar]
 Sutton EK. 2018. A new method of physicsbased data assimilation for the quiet and disturbed thermosphere. Space Weather 16(6): 736–753. https://doi.org/10.1002/2017SW001785. [CrossRef] [Google Scholar]
 van den IJssel J, Doornbos E, Iorfida E, March G, Montenbruck O. 2019. Thermosphere densities derived from Swarm GPS observations. Adv Space Res 65(7): 1758–1771. https://doi.org/10.1016/j.asr.2020.01.004. [CrossRef] [Google Scholar]
 Walker A, Mehta PM, Koller J. 2014. Drag coefficient model using the Cercignani–Lampis–Lord gassurface interaction model. J Spacecr Rockets 51(5): 1544–1563. https://doi.org/10.2514/1.A32677. [CrossRef] [Google Scholar]
 Wilson GR, Weimer DR, Wise JO, Marcos FA. 2006. Response of the thermosphere to Joule heating and particle precipitation. J Geophys Res 111: A10314. https://doi.org/10.1029/2005JA011274. [CrossRef] [Google Scholar]
Cite this article as: March G, van den IJssel J, Siemes C, Visser P.N.A.M, Doornbos E.N, et al. 2021. Gassurface interactions modelling influence on satellite aerodynamics and thermosphere mass density. J. Space Weather Space Clim. 11, 54. https://doi.org/10.1051/swsc/2021035.
Supplementary materials
Supplementary files provided by the authors. (Access here)
All Tables
Maximum force coefficient differences for different α_{E} ranges along X and Z satellite reference frame axes for the simulated attack angles.
Maximum force coefficient differences for different α_{E} ranges along X, Y, and Z satellite reference frame axes for the simulated sideslip angles.
All Figures
Fig. 1 Representation of sideslip and attack angles for GRACE. 

In the text 
Fig. 2 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the GRACE satellite. 

In the text 
Fig. 3 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the CHAMP satellite. 

In the text 
Fig. 4 Polar plots of the aerodynamic coefficients as a function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the Swarm satellites. 

In the text 
Fig. 5 CHAMP aerodynamic drag coefficients as a function of the argument of latitude in three selected days describing high, moderate and low solar activity. 

In the text 
Fig. 6 GRACE aerodynamic drag coefficients as a function of the argument of latitude in three selected days describing high, moderate and low solar activity. 

In the text 
Fig. 7 Long period density ratios with respect to the NRLMSISE00 model versus the energy accommodation coefficient for different ranges of F_{10.7}. 

In the text 
Fig. 8 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under high solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 

In the text 
Fig. 9 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under moderate solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 

In the text 
Fig. 10 Comparison between the new densities and the ones derived by Mehta et al. (2017), Doornbos (2011) and SESAM (Pilinski et al., 2013b) for the CHAMP satellite under low solar activity condition. Winds along longitudinal (XSC) and crosstrack (YSC) satellite directions are computed with HWM07 model. 

In the text 
Fig. 11 Yaw angle variations and density ratios between Swarm A and C measurements during the May 2014 manoeuvre. Selected orbit numbers are available in the top and bottom parts of the figure. 

In the text 
Fig. 12 Swarm A and C density ratios for periods 11 and 17 (on the left and right side, respectively). The optimal density ratios for which satellites are in better agreement are highlighted with dashed lines. 

In the text 
Fig. 13 Swarm A and C optimal density ratios for May 2014 manoeuvres and selected periods. 

In the text 
Fig. 14 Comparison of the density ratios between the newly derived density and the NRLMSISE00 and WACCMX models results. The densities are normalized at 400 km altitude except for GOCE, which is normalized at 250 km. 

In the text 
Fig. 15 Comparison between the mean density ratios (μ*) between the NRLMSISE00 model (top) and DTM2013 (below) with different data sets for the CHAMP, GRACE, GOCE, and Swarm satellites. The percentage values indicate the maximum variability range among the satellites. Note that for GRACE and Swarm B the new density ratios are computed using a preliminary energy accommodation coefficient of 0.85 (which is not validated as for the other satellites). 

In the text 
Fig. A.1 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the CHAMP satellite. 

In the text 
Fig. A.2 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the CHAMP satellite. 

In the text 
Fig. A.3 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the GRACE satellite. 

In the text 
Fig. A.4 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the GRACE satellite. 

In the text 
Fig. A.5 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the GOCE satellite. 

In the text 
Fig. A.6 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle and energy accommodation coefficient for the GOCE satellite. 

In the text 
Fig. A.7 Aerodynamic coefficients in the satellite reference frame as function of the attack angle and energy accommodation coefficient for the Swarm satellites. 

In the text 
Fig. A.8 Aerodynamic coefficients in the satellite reference frame as function of the sideslip angle (between 0° and 360°) and energy accommodation coefficient for the Swarm satellites. 

In the text 
Fig. B.1 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the CHAMP satellite for fixed sideslip angles. 

In the text 
Fig. B.2 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the GRACE satellite for fixed sideslip angles. 

In the text 
Fig. B.3 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the GOCE satellite for fixed sideslip angles. 

In the text 
Fig. B.4 Aerodynamic coefficients in the satellite reference frame as function of the speed ratio and energy accommodation coefficient for the Swarm satellites for fixed sideslip angles. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.