Gas-surface interactions modelling influence on satellite aerodynamics and thermosphere mass density

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 gas-surface interactions (GSI), which is part of the satellite aerodynamic modelling used in the acceleration to density data processing. Physics-based GSI models require in-situ atmospheric composition and temperature data that are not measured by any of the above-mentioned 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 self-consistency of the multi-mission thermosphere mass density data sets. We compare our results with those based on variable energy accommodation obtained by different studies and semi-empirical 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 semi-empirical models such as NRLMSISE-00 and DTM-2013 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 NRLMSISE-00 and DTM-2013, respectively. A comparison with the WACCM-X general circulation model is also performed. Similar to the other models, WACCM-X 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.


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 low-Earth orbit followed. The Gravity Recovery and Climate Experiment (GRACE), the Gravity field and steady-state Ocean Circulation Explorer (GOCE), Swarm, and the GRACE Follow-On (GRACE-FO) 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 before-mentioned 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 first-guess gas-surface 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 gas-dynamics modelling has been made in order to raise the level of accuracy of the density data sets (Pilinski et al., 2013a(Pilinski et al., , 2016Mehta 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 particle-surface collisions between atmospheric molecules and satellite surfaces is fundamental (March et al., 2019b).
The energy accommodation coefficient (a 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 (a E = 1 and a E = 0) describe particle-surface 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 in-situ 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, in-situ concentration and temperature observations are not available for any of the above-mentioned 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 physics-based 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. Mission-dependent biases, if not taken accurately into account, can cause problems in data assimilation, both in the empirical and physics-based 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 thermosphere-ionosphere system's energy budget (e.g., Wilson et al., 2006). Bias-free data will also be essential for investigating long-term trends over multiple solar cycles, induced by changes in solar activity and by cooling the upper-atmosphere due to increasing concentrations of greenhouse gases (e.g., Emmert et al., 2008;Qian et al., 2013). With the continuation of the Swarm and GRACE-FO 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 inter-mission biases are investigated and from which long-term trends and their uncertainties can be derived. Long-term trends strongly affect estimates of an in-orbit lifetime of future satellite missions, as well as the long-term 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 G. March et al.: J. Space Weather Space Clim. 2021, 11, 54 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 self-consistently 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 above-mentioned physics-based 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 "self-consistency", 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 self-consistency was achieved for a E = 0.85. Zonal winds within the analysed attitude manoeuvre had a lower variability under similar conditions for this new a 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/gocethermosphere-data?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 side-by-side 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 along-track 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 inter-compare 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 inter-compare 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 high-fidelity 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 Rarefied-gas Time-accurate 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.

Methodology
This study is based on the simulation output from the Stochastic PArallel Rarefied-gas Time-accurate Analyzer G. March et al.: J. Space Weather Space Clim. 2021, 11, 54 (SPARTA) software (Gallis et al., 2014) (https://sparta.github.io, August 2021) in combination with the new high-fidelity 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 particle-particle collisions as validation of such assumption. From technical drawings and pre-launch 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 NRLMSISE-00 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 gas-surface collisions are not assumed to be fully accommodated reflections (a 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 a 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: 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, 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 a 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 a E values (beyond the four simulated conditions) have been generated in the post-processing phase using a least-squares method. The processing was further validated with a few new SPARTA simulations generating results for a new set of a 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 a 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 in-situ 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 self-consistency 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 side-slip 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: along-track), horizontal sideways (Y: cross-track) and vertical downward (Z: anti-radial) 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 side-slip 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. 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 on-board accelerometers provided high-precision 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 GPS-derived accelerations was introduced as part of the official ESA documentation (van den IJssel et al., 2019). For CHAMP and GOCE the cross-track 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 platform-related issues) and especially for the small aerodynamic acceleration signal typical of high altitudes (around and over 500 km) (March et al., 2019b).

Gas-surface 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 side-slip 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 side-slip angles. The C x , C y , and C z components are available in the figure for the same range of a E values. During mission lifetime, the attack angle is mainly centered around the nominal flight configuration of 0°in attack angle, while the side-slip 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 a 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 a 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 a E range. When we inspect the differences between non-accommodated 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 cross-track 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 Z-component 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 anti-flight 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.
Appendix A shows the force coefficients for the CHAMP, GRACE, GOCE, and Swarm satellites for different attack and side-slip angles. For GOCE, large attitude manoeuvres were not performed. Therefore, the side-slip 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. G. March et al.: J. Space Weather Space Clim. 2021, 11, 54 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 side-slip 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 2002-10-27, 2005-05-15, and 2009-08-28, 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 a 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 NRLMSISE-00 within such activity levels. For GRACE is also illustrated the drag coefficient evolution for the energy accommodation coefficient of 0.85 in order to inter-compare 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.

Gas-surface interactions influence on neutral density
In this article, the optimal a 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 a 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 long-term averages. The differences with semi-empirical models are influenced by solar and geomagnetic activity because models perform differently for different geomagnetic and solar activity inputs (Emmert, 2015). The semi-empirical 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 physics-based 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 semi-empirical model. In Figure 7, the density ratio with respect to the NRLMSISE-00 model is shown for three different solar activity levels. The high activity is represented by 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 semi-empirical 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 NRLMSISE-00 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.
From this analysis, it is not straightforward to retrieve an optimal a 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 empirical-model (Pilinski et al., 2013b), which in this case adopts the aerodynamic information generated with the high-fidelity 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 HWM-07 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 low-frequency variations with respect to the new density. This is associated with the use of the same chemistry inputs from the NRLMSISE-00 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 a 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    March et al.: J. Space Weather Space Clim. 2021, 11, 54 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.
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 a 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 a 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 a E oscillates between the two fixed values of 0.85 and 0.93. In agreement with these results, interesting information about the a 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://iswat-cospar.org/) which involve part of the authors and will be further explained in a dedicated article.

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 Kalman-filter approach G. March et al.: J. Space Weather Space Clim. 2021, 11, 54 in the GPS-derived 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.
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 a 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 a 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 a 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 a 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 a E value of 0.93 (Fig. 13).
The results coincide with the values obtained from the accelerometer-derived wind analysis of the CHAMP and GOCE satellites in March et al. (2019b), which showed optimal a 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 a E value of 0.93 (available at https://swarm-diss.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 a E needs to be introduced in the calculations. In the next section, this will also be verified using longer time periods and inter-comparing the new densities with other sources and methods.

Comparisons with atmospheric and GSI models
In order to evaluate the new density data, further analysis is performed using the results from the NRLMSISE-00 and DTM-2013 semi-empirical 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, NRLMSISE-00 is one of the most widely used models, and DTM-2013 is one of the latest semiempirical models. A further comparison with the WACCM-X general circulation model (Liu et al., 2018) is also performed using a set of representative months. The WACCM-X model is a physics-based 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 NRLMSISE-00 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 WACCM-X upper pressure level altitude boundary extends to around 500 km altitude. The monthly time window was selected to limit the computational cost of the physics-based model simulations.
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 a E value of 0.85 was adopted for this comparison with WACCM-X 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 & Mayer-Gü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 NRLMSISE-00 and WACCM-X 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 WACCM-X than with NRLMSISE-00. For high activity, this is clearly visible for all the selected satellites. For GOCE, the density ratios show high agreement for both NRLMSISE-00 and the physics-based 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 WACCM-X 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 NRLMSISE-00 and DTM-2013 models have been computed (Fig. 15). The illustrated periods are between 2003-01-01 and 2005-01-01 for CHAMP and GRACE, 2011-01-01 and 2013-01-01 for GOCE, and 2014-07-19 and 2017 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.
The newly derived densities (labeled as "This work") are compared with the previous results from March et al. (2019a) for fully accommodated reflections (a 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 a E = 0.93 with new high-fidelity 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 NRLMSISE-00 model, and at the bottom the ones for DTM-2013. The figure shows the mean density ratio (l*) of the model and the observations as an index of their consistency. Under each group of bars, we provide the difference in percentage (D) between the maximum and minimum mean densities among the selected satellites. This is assumed to be a measure of self-consistency for multi-mission 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 semi-empirical models is not improving, the results are more consistent among the selected missions. Indeed, the scale with respect to the NRLMSISE-00 model is now more constant across the missions. Since our primary objective is a high self-consistency of the multi-mission density observations, we focus in the discussion of Figure 15 on the differences between the D values. The new a E bars from this study are within 4.13% of the variation among the satellites, whereas for the previous version of the data a E = 1.00 this value reached 8.26%. Introducing the a 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 NRLMSISE-00 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 a 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 self-consistency between the selected missions and both empirical models is a promising result and certainly a starting point for further improvements. Comparison of the density ratios between the newly derived density and the NRLMSISE-00 and WACCM-X models results. The densities are normalized at 400 km altitude except for GOCE, which is normalized at 250 km.

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 in-orbit lifetimes of space objects are longer, and this has adverse effects on the low orbit space debris population. Very recently, newly released or in-progress 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 long-term change effects.
The continued use of models based on older data in the acceleration-to-density 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 Fig. 15. Comparison between the mean density ratios (l*) between the NRLMSISE-00 model (top) and DTM-2013 (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).
G. March et al.: J. Space Weather Space Clim. 2021, 11, 54 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 a E = 0.93 and are presently available on the before-mentioned 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 high-fidelity 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 night-side, 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 in-situ 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.