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



Article Number  53  
Number of page(s)  12  
DOI  https://doi.org/10.1051/swsc/2021034  
Published online  19 October 2021 
Research Article
Uncertainty quantification of the DTM2020 thermosphere model
OMP/GETCNES, Space Geodesy Office, 14 Avenue E. Belin, 31401 Toulouse cedex 4, France
^{*} Corresponding author: claude.boniface@cnes.fr
Received:
11
June
2021
Accepted:
5
September
2021
Aims: The semiempirical Drag Temperature Models (DTM) calculate the Earth’s upper atmosphere’s temperature, density, and composition. They were applied mainly for spacecraft orbit computation. We developed an uncertainty tool that we implemented in the DTM2020 thermosphere model. The model is assessed and compared with the recently HASDM neutral density released publicly in 2020. Methods: The total neutral density dataset covers all highresolution CHAMP, GRACE, GOCE, and SWARM data spanning almost two solar cycles. We constructed the uncertainty model using statistical binning analysis and leastsquare fitting techniques, allowing the development of a global sigma error model to function the main variabilities driving the thermosphere state. The model is represented mathematically by a nonlinear manifold approximation in a 6D space parameter. Results: The results reveal that the altitude parameter presents the most notable error range during quiet and moderate magnetic activity (K_{p} ≤ 5). However, the most considerable uncertainty appears during severe or extreme geomagnetic activities. The comparison with density data provided by the SET HASDM database highlights some coherent features on the mechanisms occurring in the thermosphere. Moreover, it confirms the tool’s relevance to provide a qualitative database of neutral densities uncertainties values deduced from the DTM2020 model.
Key words: Thermosphere / DTM2020 / uncertainty / modeling / neutral density
© C. Boniface & S. Bruinsma, 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 thermosphere operational and research DTM2020 (Drag Temperature Model) models were developed in the framework of the Space Weather Atmosphere Models and Indices project (SWAMI, http://swamih2020.eu), which was a European Union Horizon 2020 project. Semiempirical thermosphere specification models are used to compute the atmospheric drag force in the orbit prediction of objects in Low Earth Orbit (LEO). Conjunction Assessment (CA) and subsequent satellite collision avoidance maneuver operational decisions presently are made without taking its utmost contributor to uncertainty, the thermosphere density estimation, properly into account in the orbit determination process (Bruinsma et al., 2021b; Hejduk & Snow, 2018). Until recently, propagation of the neutral density uncertainties to orbital state uncertainties has received little attention, but awareness for this theme is growing (e.g., Vallado & Finkleman, 2014; Emmert et al., 2017; Schiemenz et al., 2019; LopezJimenez et al., 2021). A key hurdle to including the neutral density uncertainty in a space object’s covariance is that thermosphere models do not deliver this information up to now. In this work, we developed a semianalytical model to characterize the uncertainty of atmospheric modeling for DTM2020 as a function of location (latitude, local time, altitude), season (day of the year), and two indices (F10.7 and K_{p}) describing the solar and geomagnetic activities. The model calculates the 1sigma uncertainty of Observed to Computed density ratio (“O/C”), which combines statistical data binning techniques and the leastsquares fitting procedure.
The drag force represents the largest source of uncertainty when extrapolating the orbit of spacecraft and space debris in LEO. Upper atmosphere density varies considerably in space and time due to highly variable solar and geomagnetic activity and seasonal effects (e.g., Guo et al., 2007; Lathuillère et al., 2008; Qian & Solomon, 2011; Emmert, 2015). Consequently, atmospheric drag on satellites varies significantly due to the proportional relation with thermospheric mass density. Other error sources in the drag computation, which are not addressed in this paper, are due to the unknown attitude of the object, approximate knowledge of its mass (Hejduk & Snow, 2018), and its aerodynamic drag coefficient (Mehta et al., 2014). The uncertainty in the thermospheric density modeling is generally 10–15% (Vallado & Finkleman, 2014) for a quiet thermosphere state. However, it can increase rapidly depending on the solar and geomagnetic activity, notably in the case of geomagnetic storms (e.g., Liu & Luhr, 2005; Forbes et al., 2005; Sutton et al., 2009; Emmert, 2015).
Provision of the model uncertainty was also one of the recommendations in COSPAR’s global roadmap 2015–2025 for Space Weather (Schrijver et al., 2015). The uncertainty model for DTM2020 is based on data in the altitude range from 200 to 500 km. Below and above that altitude interval, density data availability is too sparse due to too high and too low drag values, respectively. Neutral density decreases exponentially with height, and satellite drag is too large for satellites to maintain orbit without propulsion below approximately 350 km. In contrast, above approximately 500 km, satellite drag is usually too small to be accurately derived from orbital data. However, owing to the importance of reentry predictions on the one hand (Visser & Van den Ijssel, 2016; Hejduk & Snow, 2018) and the primary operating altitude of Earth Observation satellites at around 800 km, the uncertainties are provided for the altitude range 100–1000 km, by making presumptions described in the paper.
Section 2 introduces the datasets used in the model building, and Section 3 briefly describes the DTM2020 thermosphere model. The data analysis and modeling technique are presented in Sections 4 and 5, respectively, whereas the results are discussed in Section 6. A fully independent model validation using SET HASDM data (Tobiska et al., 2021) is specified in Section 7, and the paper is finally summarized in Section 8.
2 Neutral density data
2.1 Model development density datasets
Four density datasets from different satellite missions were used in the study to build the uncertainty model. These are only briefly described in this section.
CHAMP (Challenging Minisatellite Payload) was a small satellite mission dedicated to gravity and atmospheric research led by the GFZ institute (Reigber et al., 1996). The satellite was launched on 15 July 2000 in a nearcircular and nearly polar orbit with an inclination of 87.3° at an altitude of nearly 460 km. Densities inferred from accelerometer measurements, computed by TU Delft (Doornbos, 2011) using highfidelity satellite shape and aerodynamic coefficient models, are available for the entire 10year mission.
GRACE (Gravity Recovery and Climate Experiment) was a joint NASA and DLR mission to measure temporal variations in Earth’s gravity field (Tapley et al., 2004). The identical pair of spacecraft was launched on 17 March 2002 into a nearcircular and polar orbit. Densities were inferred for the entire mission from July 2002, at 500 km altitude, to December 2016, at 350 km due to natural decay.
GOCE (Gravity Field and SteadyState Ocean Circulation Explorer) was an ESA satellite mission to measure the Earth’s gravity field and model the geoid with high accuracy (Drinkwater et al., 2003). It was launched on 17 March 2009 in a 96.5° inclination, nearcircular orbit. Its science mission started in November 2009, at a mean altitude of 255 km. The orbit was lowered several times over the mission, and its lowest at 224 km was reached in May 2013. TU Delft derived densities under the ESA contract, using highfidelity satellite shape and aerodynamic coefficient models, from the thruster data from November 2009 to September 2013 (Bruinsma et al., 2014).
The three aforementioned satellites carried highresolution and highprecision accelerometers (Touboul et al., 1999), increasingly precise from CHAMP, GRACE to GOCE. The total densities are inferred from the 0.1 or 0.2 Hz measurements of nongravitational accelerations (even if indirectly in the case of GOCE). The satellite speed is approximately 8 km s^{−1}, which translates to an alongtrack resolution of approximately 80 or 40 km, respectively, between one and two orders of magnitude smaller than the thermosphere model resolution.
Finally, SWARM is an ESA Earth observation mission consisting of a constellation of three satellites launched in November 2013 (FriisChristensen et al., 2008). SWARM A and SWARM C are orbiting sidebyside at an altitude of 462 km at an inclination of 87.35°, whereas SWARM B is flying at a higher altitude of 511 km at the same inclination. The neutral density used in this study was deduced from GPS tracking data on SWARM A (Van den Ijssel et al., 2020) instead of the accelerometer due to severe problems with that instrument. At the altitude of about 460 km, neutral density can only be derived accurately for high and medium solar activity and with a spatial resolution of the order of thousands of kilometers with GPS data. In this study, SWARM A densities for F10.7 less than 100 sfu were rejected.
Figure 1 shows the time evolution of neutral density provided by the four satellites. The database 176 spans 2001 to 2019, i.e., covering almost two solar cycles (23 and 24).
Fig. 1 Neutral density derived from CHAMP, GRACE, GOCE, and SWARM satellites (2001–2019). Neutral density is in the log scale with the monthlyaveraged neutral density in black color. Solar Flux F10.7 is in red. 
2.2 Validation density dataset
Space Environment Technologies (SET) has recently released 20 years of neutral densities through the SET HASDM Database (Tobiska et al., 2021), which is derived from the US Space Force High Accuracy Satellite Drag Model (Storz et al., 2005). However, these densities are not observations but outputs of the HASDM model, which is a Jacchia model corrected in near realtime through data assimilation. The database provides neutral density every 3 h from 2000 to 2019 on 10° latitude × 1 h local time grids, and 25 km stepped altitude from 175 to 825 km. We validated the uncertainty model using the SET HASDM data on the altitudes 250, 325, 400, and 475 km for the entire 20year time span.
3 The DTM2020 thermosphere model
The DTM2020 model is briefly discussed in this section. A companion paper (Bruinsma & Boniface, 2021) details the improvements of this model and presents a complete assessment and comparison with the CIRA thermosphere models. The Drag Temperature Model (DTM) is a semiempirical model reproducing temperature, density, and composition of the thermosphere climatology. The total density data of the four satellites were assimilated in the DTM2020 thermosphere model. We developed a 1sigma uncertainty estimate of the neutral density for operational purposes. That is why the model is constructed based on the operational DTM2020, which uses the classical solar and geomagnetic drivers F10.7 and the planetary index K_{p}. Although the solar radio flux at 30 cm F30 (Dudok de Wit et al., 2014) and the highcadence Hpo geomagnetic indices (Matzka et al., 2019) represent variations in upper atmosphere heating accurately than F10.7 and K_{p}, they are not yet certified for operational use. We emphasize that the uncertainty model developed in this study only gives the uncertainty of the DTM2020 model density and not of atmospheric drag, for which more sophisticated tools like state vector covariance of space objects are required (Hilton et al., 2019).
4 Data processing and analysis method
Several works have studied uncertainty in the thermospheric density in orbit prediction, but their quantitative impact on the orbital dynamics has not been thoroughly investigated. Most of them concentrated on the errors in neutral density due to inaccurate space weather inputs and forecasts because these errors make a significant impact on orbit prediction (Anderson et al., 2009; Schiemenz et al., 2019). Therefore, density models must be adjusted, and the uncertainties in the model must be estimated. Several techniques exist for density model calibration (Doornbos et al., 2008; Chen et al., 2019). However, these techniques do not directly assess the uncertainty in the density.
Furthermore, longterm orbit prediction is computationally very timeconsuming (Petit & Lemaitre, 2016). Thus, the accurate prediction of neutral density and atmospheric drag for orbital uncertainty propagation generally considers simplified thermospheric mass density (Giza et al., 2009). Our goal is not to forecast the orbit of satellites but to provide an upstream estimate of neutral density uncertainties calculated with the DTM2020 thermosphere model for operational use. That is why a simplified statistical approach is preferred and described in the following section.
4.1 Data fusion
The neutral density datasets provided by the satellites GOCE, CHAMP, GRACE, and SWARM are combined. Data fusion by integrating the multiple data sources produces a complete neutral density database in the range of 200–500 km. Then, the new version of the semiempirical DTM model (see Sect. 3), using the classical K_{p} and F10.7 drivers, is used to compute the neutral densities to assess the uncertainty of the model’s prediction. The model and data are not compared by computing density residuals represented as absolute differences (observed minus computed). Indeed, an observedtocomputed density ratio (O/C) is preferred because it is a more appropriate quantity to express a model’s ability to reproduce the observations. Using density ratios comprehensively allows interpreting all features of the complex dynamic changes of neutral density as a function of geographical coordinates (latitude, altitude, local time) and solar and geomagnetic drivers (Bruinsma, 2015). Moreover, due to its skewed distribution, the mean density ratio is determined in log space. Equations (1) and (2) are used for the statistical analysis (Bruinsma et al., 2018):(1)(2)where 〈O/C〉 is the mean observedtocomputed density ratio and Std (O/C) is the standard deviation expressed in percent.
Several tests were performed to ensure the relevance of the data fusion process technique. Indeed, the neutral density provided by each satellite was analyzed separately. We observed similar 1sigma uncertainty trends. The histogram of the observed densities to the corresponding DTM2020 model densities (log scale) is illustrated in Figure 2. The neutral density distribution and a normal distribution N(μ = 0.02, σ = 0.2) present a relatively good agreement.
Fig. 2 Histogram observed densities to the corresponding DTM2020 model densities (log scale). Fusion data integrated GOCE, CHAMP, GRACE, and SWARM neutral densities datasets. The normal distribution N (μ = 0.02, σ = 0.2) is in red. 
4.2 Binning data
Statistical data binning is applied on the neutral density database to group numbers of more or less continuous values into a smaller number of “bins.” Then, the neutral density is averaged in each interval. This binning process on the density ratios is helpful to assess the quality of the model according to physical parameters. This technique allows the understanding of how the variability of physical parameters affects the density predictions and helps quantify the modeling errors (Bruinsma, 2015). Figure 3 illustrates the density ratios in 10° bins (left) and twohour local time bins (right) for DTM2020. The grey shaded bands represent the 1sigma error. The model has a higher error near the poles, where it can reach 20–25%, whereas it exhibits a semidiurnal residue (more than 20% at midnight).
Fig. 3 (Left) Density ratios binned in latitude of DTM2020, (Right) density ratios binned in local solar time of DTM2020. 
5 The 1sigma uncertainty semianalytical model
5.1 1sigma uncertainty polynomial functions
Density ratios were binned in several physical parameters of DTM2020 (latitude, altitude, day of the year, local time, solar flux, and geomagnetic index). The equations representing the best fit in the leastsquares sense were computed for each parameter. The leastsquare fitting model to data has polynomials function forms. Equation (3) summarizes the set of coupled equations:(3)where x_{alti}, x_{alat}, x_{hl}, x_{day}, x_{Fp}, x_{Kp} are altitude, latitude, local time, day of the year, solar flux, and geomagnetic index variables. σ_{x} is the 1sigma uncertainty for each parameter. Table 1 lists the coefficients {a_{i,x}} of each equation.
List of polynomial coefficient values for each parameter.
Figure 4 illustrates the 1sigma uncertainty (in percent) binned in latitude (top) and altitude (bottom) with the corresponding leastsquares fits to the data. The 1sigma error presents a quadratic evolution in latitude ranging from 22% in the tropics to 24–26% to the poles. Moreover, the uncertainty increases strongly linearly with the altitude (from 10% at 230 km to 30% at 500 km).
Fig. 4 (Top) 1σ uncertainty in percent binned in latitude and the model fitting to data in the leastsquares (LS) sense, (bottom) 1σ uncertainty in percent binned in altitude, and the model fitting to data in the leastsquares sense. Data is in blue, and the LS fitting in red color. 
We applied the exact process to the other parameters (see Fig. 5). The main trends are:
A quadratic evolution is noted in solar local time (from 20% at 15 h to 30% at midnight).
A constant error of 23% as a function of the day of the year.
The 1sigma error decreases linearly with the solar flux (from 25% at low solar activities to 20% at high solar activities).
The 1sigma uncertainty has a cubic dependency with the magnetic index (from 20% to 60% at high K_{p}Index). The uncertainty has the same order of values as the other parameters at low K_{p} (between 20 and 30% for low or moderate magnetic activities). A substantial increase in error appears for higher K_{p} (typically for severe geomagnetic storms).
Fig. 5 1σ uncertainty (in percent) leastsquare fitting: (topleft) 1σ uncertainty as a function of local time, (topright) 1σ uncertainty as a function of day of the year, (bottomleft) 1σ uncertainty as a function of solar flux, (bottomright) 1σ uncertainty as a function of K_{p}index. 
To summarize, the altitude parameter presents the largest uncertainty range during quiet and moderate geomagnetic activity (K_{p} ≤ 5). For higher K_{p}, the most significant uncertainty is mainly due to the geomagnetic activity (i.e., severe or extreme geomagnetic storms).
5.2 The leastsquare fit over manifolds
We performed a simplified procedure for fitting a set of leastsquares polynomial functions like those obtained in equation (3) to map a manifold structure (Sober et al., 2020; RubinDelanchy, 2021). The 1sigma uncertainty is deduced from the contribution of overall uncertainty described in equation (3). We assume that the 1sigma uncertainty is defined, in the space of N dimensions, by the following relations:(4)(5)(6)where {a_{i,k}} are the coordinateswise slopes, {x_{i}} the variables, and A is a constant. The coefficient k is related to the degree of the polynomial, and is the maximum degree of the polynomial for the variable x_{i}. Equation (6) is indeterminate, and the constant A must be approximated. For that, we postulate that A is inferred from equations (7) to (8):(7)(8)where is the mean value of each distribution of the variable except for x_{j}, chosen as a reference parameter. Then, the nonlinear manifoldlike equation (9) is solved stepbystep to estimate the 1sigma uncertainty:(9)
In this work, the 1sigma relation is simplified in 6D dimensions with the parameters latitude, altitude, local time, day of the year, solar and geomagnetic index.
The following equation approximates the 1sigma uncertainty by using the set of equations (3) and the process (4) → (9) for j = 0:(10)
The coefficients of the polynomial functions are related respectively to altitude, latitude, local time, day of the year, solar flux, and geomagnetic index. denotes the mean value of each physical parameter. The altitude is the reference variable in the equation. The averaged values are deduced from the neutral density distribution: = 415 km, = 0.15°, = 12 h, , .
The following remarks can be made here. First of all, a direct leastsquare approach on the overall neutral density dataset (i.e., without using binning process) may be possible. But due to the very large number of sparsely distributed data (not sparse and so much more for the HASDM dataset), the problem is too timeconsuming to solve. The main idea was to use the leastsquares mechanisms on a local rather than global level. This simplified approach can capture the local trends of the data to reconstruct better broader sets of functions like the one obtained in equation (3). Secondly, the first modeling performed to reduce the problem linearly was done using hyperplane (Reister, 1997) instead of manifold techniques. Indeed, density ratios were binned using piecewise linear functions deduced in the leastsquare sense. We observed that the two approaches lead to practically the same results (not presented in the paper). We adopted the manifold technique, allowing both a fitter model charting representation and a more efficient code. Finally, we highlight that the model calculates a mean 1sigma uncertainty of 22%, in ecellent agreement with the standard deviation deduced from the complete neutral density distribution (see Fig. 3).
5.3 1sigma error model
Equation (10) is well justified in the range 200–500 km. These presumptions were made to etend the altitude range to 100–1000 km due to the unavailability of data:
at a lower altitude (< 200 km), the worst case was chosen by fixing the evolution of 1sigma to their values at 200 km,
at a higher altitude (> 500 km), an optimistic case was assumed by keeping the evolution of 1sigma for an altitude of 500 km.
Therefore, the global model describing the 1uncertainty in the range [100 km–1000 km] is:(11)where is the function defined in equation (10).
Table 2 displays 1sigma error at different altitudes (from 100 to 1000 km) and for latitude = 0°, local time = 12 h, day of the year = 180, and two values of solar flux = 65 and 200 (low and high solar activity) and geomagnetic activity Kp = 4 and 8 (enhanced and high geomagnetic activity). The 1σ uncertainty is growing with altitude. The values are five times higher at high solar and low geomagnetic activities passing from 6.84% at 100 km to 30.60% at 1000 km. Likewise, the error is almost two times greater at low altitude (below 200 km) during low solar activity (6.84% at , against 11.37% at ). Furthermore, the 1σ uncertainty can have the same order of magnitude for two different altitudes due to the level of geomagnetic activity (34.19% at 100 km for against 35.13% at 1000 km for ).
Onesigma uncertainty as a function of altitude for different solar flux values and magnetic index K_{p} values.
6 Results and discussion
The model computes the 1sigma error in 6D parameter space, including geographical coordinates (latitude, altitude, solar local time), season (day of the year), solar and geomagnetic drivers (F10.7 proxy and Kp index). As discussed previously, altitude is one of the main drivers in error modeling during quiet and disturbed thermosphere states, and so this parameter is compared with the other physical quantities. The plots are limited to the range of 200–500 km to illustrate the results. Figure 6 shows 3D surface plotting of the 1sigma uncertainty as a function of latitude, local time, and altitude parameters. Curvilinear shapes are observed due to latitude (left side) and local time (right side) effects. The 1sigma uncertainty is minimum at the lowest altitude in the tropics (~5–10%), whereas a maximum of 30–35% is localized at the highest altitude (here 500 km) over the poles. It presents a predominantly diurnal variation with a peak uncertainty (~36%) at 500 km altitude at midnight. The error is minimum at the lowest altitude at 15 h solar local time (~6%). Moreover, the ratio between the uncertainty on the night side (22:00–05:00 LT) and dayside (10:00–17:00 LT) linearly decreases with altitude (1.2 at 500 km and 1.6 at 250 km). The latitude effects are also altitudedependant, with an uncertainty ratio pole/equator of 1.1 at 500 km against 1.4 at 250 km.
Fig. 6 3D plot of the 1σ uncertainty (in percent) as a function of space variabilities: (left) altitude and latitude, (right) altitude, and local time. 
The density error exhibits a solar cycle variation, driven by the intrinsic ∼11year variability of magnetic activity in the Sun (Qian & Solomon, 2011), as shown in Figure 7 (right frame). The maximum uncertainty is observed at a given altitude during low solar activity. For instance, it reached 22% in 2008 (minimum of solar cycle 23) while 16% error is deduced in 2002 at 400 km (high solar activity). A 5% uncertainty variation is observed between low and high solar activity in agreement with the range of solar flux variations estimated in Figure 9, left side. We remark that this modulation pattern matches the thermospheric variations due to the solar cycle (Solomon et al., 2011).
Fig. 7 3D plot of the 1σ uncertainty (in percent) as a function of altitude and time evolution: (left) day of the year (right) years. 
Surprisingly, the uncertainty does not exhibit seasonal variations, as illustrated in Figure 7 (left frame). However, as shown in Figure 8, different variations on the 1sigma error are perceived during low solar activity (minimum solar cycle in 2008) and high (2002).
Fig. 8 Density ratios 1σ uncertainty binned in the day of the year of DTM2020 with a 27day averaged in red color, for two periods (2001–2003 and 2007–2009, respectively maximum and minimum of solar cycle 23). 
The 1sigma uncertainty is minimum at the lowest altitude during high solar activities (~4% for at 200 km), while the maximum uncertainty of 30% is reached at the highest altitude during low solar activities (, as shown in Figure 9 (left frame). Whatever the altitude, we noted a 5% decrease in 1sigma from low to high solar activities. As a final point, the error is significantly enhanced by increasing geomagnetic activity to storm level, as illustrated in Figure 9 (right frame). Indeed, although the 1sigma uncertainty is minimum at the lowest altitude during quiet geomagnetic activities (~6% for K_{p} = 1 at 200 km), the maximum of 65% is obtained at 500 km altitude during extreme geomagnetic storms (K_{p} = 9 We observed a clear transition at K_{p} = 5 (low/medium geomagnetic activity K_{p} ≤ 5 and high geomagnetic activity K_{p} > 5). Indeed, independent of altitude, the error increases by 6% at low/medium geomagnetic activities, while we observe a 30% uncertainty increase from K_{p} = 5 to K_{p} = 9 (high/severe magnetic activities). It is wellknown that thermospheric density variations and uncertainties are enhanced in response to geomagnetic activity by the energy deposition from the solar wind and magnetosphere (Knipp et al., 2013; Bruinsma et al., 2021a).
Fig. 9 3D plot of the 1σ uncertainty (in percent) as a function of altitude and drivers : (left) solar flux, (right) magnetic index. 
We focus on the influence of geographical coordinates and solar/ geomagnetic activities in Figure 10. The 1sigma uncertainty is plotted as the latitude and local time function on the left side. A reverse saddlelike shape is observed. The maximum uncertainty occurs near the 24 LT in the pole region (~32%). In contrast, a minimum uncertainty of 20% is expected in the late afternoon near the equatorial zone. Therefore, the largest uncertainties are at high latitudes in the night/early morning (22:00–04:00 LT), whereas they are smallest in late morning and afternoon (10:00–18:00 LT) from the low to midlatitudes. Figure 10 (right frame) shows that the 1sigma uncertainty is minimum during high solar and low geomagnetic activity (~21% for F_{10.7} = 200 and K_{p} = 1). The maximum uncertainty of 60% is observed for F_{p} = 65 and K_{p} = 9. A clear transition always appears at K_{p} = 5. Moreover, independent of solar activity, an increase of 5% occurs at low geomagnetic activity K_{p} ≤ 5 whereas the uncertainty increases by a factor of 30% during high magnetic activities.
Fig. 10 3D surface plot of the 1σ uncertainty (in percent) (topleft) altitude and latitude, (left) geographical coordinates; (right) solar proxy and magnetic index. 
To conclude this section, we underline the chief impact of the altitude on the 1sigma uncertainty (ranging from 10 to 30%) compared to the other parameters excepting high geomagnetic activity, i.e., K_{p} > 5 (from 30 to 80%). Latitudinal and diurnal effects on the density uncertainties are observed. Finally, the variation in solar activity over a cycle impacts the neutral density uncertainty even if its amplitude is rather modest (variation in the uncertainty of 5%) and on a long time scale.
7 Validation using HASDM uncertainty model
HASDM densities at 250, 325, 400, and 475 km covering almost two solar cycles have been used to assess our uncertainty model developed with accelerometerinferred neutral density data (called UM in the following). For that, we used the same process to construct an uncertainty model but with HASDM (called HUM in the following) by binning the vast amount of data (more than 15 GB) to deduce a 6D polynomial function in the altitude range 200–500 km. The comparisons are shown in Figures 11–13 (left and right frames for the UM and HUM model, respectively). The two models present overall the same tendencies. For instance, they exhibit similar cubic shapes as a function of altitude and the local hour in Figure 11 (bottom), although a delay of two hours is observed on the maximum (15 LT for UM and 17 LT for HUM models). Besides, the longterm variation due to the solar cycle is well reproduced by the two models, as illustrated in Figure 12 (bottom frames). The 1sigma error, as a function of altitude and the magnetic index, presents the same shape with a transition at K_{p} = 5, as remarked in Figure 13 (bottom frames). One of the main features is an overall underestimation by nearly 5% of the HUM model. This difference increases significantly with the K_{p} index and reaches 20% at the highest altitude (Fig. 13 – bottom frames). We remind that the datasets have different sampling rates (and so also different spatial resolutions). The accelerometer measurements have a 10s temporal resolution, while HASDM’s is 3 h, which explains lower uncertainties due to the smoothing effect.
Fig. 11 Contour level plot of the 1σ uncertainty (in percent) as a function of geographical coordinates: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (Top) altitude and latitude, (bottom) altitude and local time. 
Fig. 12 Contour level plot of the 1σ uncertainty (in percent) as a function of temporal evolution: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (top) altitude and day of year, (bottom) altitude and years. 
Fig. 13 Contour level plot of the 1σ uncertainty (in percent) as a function of solar and magnetic drivers: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (Top) altitude and solar flux Fp, (bottom) altitude and magnetic index K_{p}. 
However, some differences need to be clarified. As shown in Figure 11 (top frames), latitude effects are not visible in the HUM model. This is due to the much lower spatial and temporal resolution of HASDM data, which does not capture the rapid and short spatialscale density variations, especially at auroral latitudes. Then, the HUM reveals seasonal effects not observed in the UM model (Fig. 10, top frames). Indeed, the UM model presents a constant 1sigma error at any given altitude, whereas the HUM uncertainty exhibits seasonal variation. The principle maximum is just after the June summer solstice (~24% at 450 km), and the secondary maximum is just after the winter solstice (~22% at 450 km). The absence of this variation in UM is most likely due to its much sparser database, but the amplitude of the variation around the mean uncertainty is only ±2%. Finally, a linear decrease as a function of solar flux is modeled with UM, whereas HUM is nearly constant on an altitude level (Fig. 13, top frames). We inferred that it is due to a DTM2020 biases effect during low solar activities (Bruinsma & Boniface, 2021) that overestimated the uncertainty value of UM model at low solar activities ().
To summarize, accelerometers deliver the most precise and highresolution measurements to infer neutral density. Comparison with HUM confirms that the UM model provides realistic estimates of the 1sigma uncertainty.
8 Conclusion
We developed a semianalytical model in 6D space parameters using the DTM2020 thermosphere model that characterizes the uncertainty of its neutral density estimates. The model is constructed with density data from 200 to 500 km altitude, but it was extended to cover 100–1000 km by making presumptions for altitudes without data coverage. The model is founded on a new combined statistical data binning and local least square fitting approach. Because the distribution of the highresolution density data is relatively sparse, the model was successfully validated in the 200–500 km altitude range using the recently released SET HASDM density data.
The principal results of the model are:
the main parameter for the uncertainty during low and moderate magnetic activity (K_{p} ≤ 5) is altitude;
the largest uncertainty is reached during severe or extreme geomagnetic storms;
the solar cycle causes a longterm modulation of uncertainty;
the model predicts latitudinal and diurnal effects on the density uncertainties.
Further study is needed to validate the uncertainty model in the extended range (100–1000 km), but that depends on future data availability.
Acknowledgments
The SWAMI project has been granted from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 776287. CNES APR equally supports SB and CB. The HASDM database is made available by SET (https://spacewx.com/hasdm/). GOCE and Swarm density data is made available by ESA (https://earth.esa.int/eogateway/). CB would like to acknowledge Julien LaurentVarin for their fruitful discussions on the mathematical method used in this work. The DTM models and the density uncertainty estimate model are available on Github via the SWAMI website (http://swamih2020.eu). The editor thanks David Vallado and an anonymous reviewer for their assistance in evaluating this paper.
References
 Anderson RL, Born GH, Forbes JM. 2009. Sensitivity of orbit predictions to density variability. J Spacecraft Rockets 46(6): 1214–1230. https://doi.org/10.2514/1.42138. [CrossRef] [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, Boniface C. 2021. The operational and research DTM2020 thermosphere models. J Space Weather and Space Clim 11: 47. https://doi.org/10.1051/swsc/2021032. [CrossRef] [EDP Sciences] [Google Scholar]
 Bruinsma SL, Doornbos E, Bowman BR. 2014. Validation of GOCE densities and thermosphere model evaluation. Adv Space Res 54: 576–585. https://doi.org/10.1016/j.asr.2014.04.008. [Google Scholar]
 Bruinsma S, Sutton E, Solomon SC, FullerRowell T, Fedrizzi M. 2018. Space Weather modeling capabilities assessment: Neutral density for orbit determination at low earth orbit. Space Weather 16(11): 1806–1816. https://doi.org/10.1029/2018SW002027. [CrossRef] [Google Scholar]
 Bruinsma S, Boniface C, Fedrizzi M, Sutton E. 2021a. Thermosphere modeling capabilities: geomagnetic storms. J Space Weather Space Clim 11: 12. https://doi.org/10.1051/swsc/2021002. [Google Scholar]
 Bruinsma S, Fedrizzi M, Lue J, Siemes C, Lemmens S. 2021b. Charting satellite courses in a crowded thermosphere. Eos 102: https://doi.org/10.1029/2021EO153475. [CrossRef] [Google Scholar]
 Chen J, Du J, Sang J. 2019. Improved orbit prediction of LEO objects with calibrated atmospheric mass density model. J Spatial Sci 64(1): 97–110. https://doi.org/10.1080/14498596.2017.1371089. [CrossRef] [Google Scholar]
 Doornbos E. 2011. Thermospheric density and wind determination from satellite dynamics. Ph.D. Dissertation. The University of Delft, 188 p. Available at http://repository.tudelft.nl/. [Google Scholar]
 Doornbos E, Klinkrad H, Visser P. 2008. Use of twoline element data for thermosphere neutral density model calibration. Adv Space Res 41(7): 1115–1122. https://doi.org/10.1016/j.asr.2006.12.025. [CrossRef] [Google Scholar]
 Drinkwater MR, Floberghagen R, Haagmans R, Muzi D, Popescu A. 2003. GOCE: ESA’s first Earth explorer core mission. Space Sci Rev 108: 419–432. https://doi.org/10.1007/9789401713337_36. [CrossRef] [Google Scholar]
 Dudok de Wit T, Bruinsma S, Shibasaki S. 2014. Synoptic radio observations as proxies for upper atmosphere modeling. J. Space Weather Space Clim 4: A06. https://doi.org/10.1051/swsc/2014003. [CrossRef] [EDP Sciences] [Google Scholar]
 Emmert J. 2015. Thermospheric mass density: A review. Adv Space Res. 56(5): 773–824. https://doi.org/10.1016/j.asr.2015.05.038. [CrossRef] [Google Scholar]
 Emmert JT, Warren HT, Segerman AM, Byers JM, Picone JM. 2017. Propagation of atmospheric density errors to satellite orbits. Adv Space Res 59(1): 147–165. https://doi.org/10.1016/j.asr.2016.07.036. [CrossRef] [Google Scholar]
 ESA. 1999. The four candidate Earth explorer core missions – Gravity field and steadystate ocean circulation explorer. European Space Agency Report SP1233(1). [Google Scholar]
 Forbes JM, Lu G, Bruinsma S, Nerem S, Zhang X. 2005. Thermospheric density variations due to the 15–24 April 2002 solar events from CHAMP/STAR accelerometer measurements. J Geophys Res 110: A12S27. https://doi.org/10.1029/2004JA010856. [Google Scholar]
 FriisChristensen E, Lühr H, Knudsen D, Haagmans R. 2008. Swarm – An earth observation mission investigating geospace. Adv Space Res 41(1): 210–216. https://doi.org/10.1016/j.asr.2006.10.008. [CrossRef] [Google Scholar]
 Giza D, Singla P, Jah M. 2009. An approach for nonlinear uncertainty propagation: Application to orbital mechanics. In: AIAA Guidance, Navigation, and Control Conference. AIAA 20096082. https://doi.org/10.2514/6.20096082. [Google Scholar]
 Guo J, Wan W, Forbes JM, Sutton E, Nerem RS, Woods TN, Bruinsma S, Liu L. 2007. Effects of solar variability on thermosphere density from CHAMP accelerometer data. J Geophys Res 112: A10308. https://doi.org/10.1029/2007JA012409. [CrossRef] [Google Scholar]
 Hejduk MD, Snow DE. 2018. The effect of neutral density estimation errors on satellite conjunction serious event rates. Space Weather 16: 849–869. https://doi.org/10.1029/2017SW00172. [CrossRef] [Google Scholar]
 Hilton S, Cairola F, Gardi A, Sabatini R, Pongsakornsathien N, Ezer N. 2019. Uncertainty quantification for Space Situational Awareness and Traffic Management. Sensors(Basel) 19(20): 4361. https://doi.org/10.3390/s19204361. [CrossRef] [Google Scholar]
 Knipp D, Kilcommons L, Hunt L, Mlynczak M, Pilipenko V, Bowman B, Deng Y, Drake K. 2013. Thermospheric damping response to sheathenhanced geospace storms. Geophys Res Lett 40: 1263–1267. https://doi.org/10.1002/grl.50197. [CrossRef] [Google Scholar]
 Lathuillère C, Menvielle M, Marchaudon A, Bruinsma S. 2008. A statistical study of the observed and modeled global thermosphere response to magnetic activity at mid and low latitudes. J Geophys Res 113: A07311. https://doi.org/10.1029/2007JA012991. [Google Scholar]
 Liu H, Luhr H. 2005. Strong disturbance of the upper thermospheric density due to magnetic storms: CHAMP observations. J Geophys Res 110: A09S29. https://doi.org/10.1029/2004JA010908. [Google Scholar]
 LopezJimenez S, Pastor A, Escobar D. 2021. Improving orbital uncertainty realism through covariance determination. Acta Astron 181: 679–693. https://doi.org/10.1016/j.actaastro.2020.09.026. [CrossRef] [Google Scholar]
 Matzka J, Stolle C, Kervalishvili G, Rauberg J, Yamazaki Y. 2019. The Hp geomagnetic index test dataset 2003, 2004, 2005, and 2017. GFZ Data Services. https://doi.org/10.5880/GFZ.2.3.2019.002. [Google Scholar]
 Mehta PM, Walker A, Lawrence E, Linares R, Higdon D, Koller J. 2014. Modeling satellite drag coefficient with response model. Adv Space Res 54(8): 1590–1607. https://doi.org/10.1016/j.asr.2014.06.033. [CrossRef] [Google Scholar]
 Petit A, Lemaitre A. 2016. The impact of the atmospheric model and of the space weather data on the dynamics of clouds of space debris. Adv Space Res. 57: 2245–2258. https://doi.org/10.1016/j.asr.2016.03.005. [CrossRef] [Google Scholar]
 Qian L, Solomon SC. 2011. Thermospheric Density: An overview of temporal and spatial variations. Space Sci Rev 168: 147–173. https://doi.org/10.1007/s112140119810z. [Google Scholar]
 Reigber C, Bock R, Förste C, Grunwaldt L, Jakowski N, Lühr H, Schwintzer P, Tilgner C. 1996. CHAMP Phase B: Executive Summary. Scientific Technical Report STR; 96/13. Potsdam: GFZ German Research Centre for Geosciences 24 p. https://doi.org/10.2312/gfz.b10396131. [Google Scholar]
 Reister DB. 1997. The leastsquares fit of a hyperplane to uncertain data. Robotica 15(4): 461–464. [CrossRef] [Google Scholar]
 RubinDelanchy P. 2021. Manifold structure in graph embeddings. arXiv:2006.05168v3. [Google Scholar]
 Schiemenz F, Utzmann J, Kayal H. 2019. Least squares orbit estimation including atmospheric density uncertainty consideration. Adv Space Res 63(12): 3916–3935. https://doi.org/10.1016/j.asr.2019.02.039. [CrossRef] [Google Scholar]
 Schrijver CJ, Kauristie K, Aylward AD, Denardini CM, Gibson SE, et al. 2015. Understanding space weather to shield society: A global road map for 2015–2025 commissioned by COSPAR and ILWS. Adv Space Res 55(12): 2745–2807. https://doi.org/10.1016/j.asr.2015.03.023. [Google Scholar]
 Sober B, Aizenbud Y, Levin D. 2020. Approximation of functions over Manifolds: a moving least LeastSquares Approach. arXiv:1711.00765v4. [Google Scholar]
 Solomon S, Quian L, Didkovsky L, Viereck R, Woods T. 2011. Causes of low thermospheric density during the 2007–2009 solar minimum. J Geophys Res 116: A00H07. https://doi.org/10.1029/2011JA016508. [Google Scholar]
 Storz MF, Bowman BR, Branson MJI, Casali SJ, Kent Tobiska W. 2005. High Accuracy satellite drag model (HASDM). Adv Space Res 36: 2497–2505. https://doi.org/10.1016/j.asr.2004.02.020. [CrossRef] [Google Scholar]
 Sutton EK, Forbes JM, Knipp DJ. 2009. Rapid response of the thermosphere to variations in Joule heating. J. Geophys Res 114: A04319. https://doi.org/10.1029/2008JA013667. [Google Scholar]
 SWAMI H2020 Program: Space Weather Models and Indexes. http://swamih2020.eu. [Google Scholar]
 Tapley BD, Bettadpur S, Watkins M, Reigber C. 2004. The gravity recovery and climate experiment: Mission overview and early results. Geophys Res Lett 31: L09607. https://doi.org/10.1029/2004GL019920. [Google Scholar]
 Tobiska W, Bowman B, Bouwer S, Cruz A, Wahl K, Pilinski M, Mehta P, Licata RJ. 2021. The SET HASDM density database. Space Weather 19: e2020SW002682. https://doi.org/10.1029/2020SW002682. [CrossRef] [Google Scholar]
 Touboul P, Willemenot E, Foulon B, Josselin V. 1999. Accelerometers for CHAMP, GRACE, and GOCE space missions: synergy and evolution. Boll Geof Teor Appl 40: 321–327. [Google Scholar]
 Vallado D, Finkleman D. 2014. A critical assessment of satellite drag and atmospheric density modeling. Acta Astron 95: 141–165. https://doi.org/10.1016/j.actaastro.2013.10.005. [CrossRef] [Google Scholar]
 Van den Ijssel J, Doornbos E, Iorfida E, March G, Siemes C, Montenbrück O. 2020. 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]
 Visser PNAM, Van den Ijssel JAA. 2016. Orbit determination and estimation of nongravitational accelerations for the GOCE reentry phase. Adv Space Res 58: 1840–1853. https://doi.org/10.1016/j.asr.2016.07.013. [CrossRef] [Google Scholar]
Cite this article as: Boniface C & Bruinsma S 2021. Uncertainty quantification of the DTM2020 thermosphere model. J. Space Weather Space Clim. 11, 53. https://doi.org/10.1051/swsc/2021034.
All Tables
Onesigma uncertainty as a function of altitude for different solar flux values and magnetic index K_{p} values.
All Figures
Fig. 1 Neutral density derived from CHAMP, GRACE, GOCE, and SWARM satellites (2001–2019). Neutral density is in the log scale with the monthlyaveraged neutral density in black color. Solar Flux F10.7 is in red. 

In the text 
Fig. 2 Histogram observed densities to the corresponding DTM2020 model densities (log scale). Fusion data integrated GOCE, CHAMP, GRACE, and SWARM neutral densities datasets. The normal distribution N (μ = 0.02, σ = 0.2) is in red. 

In the text 
Fig. 3 (Left) Density ratios binned in latitude of DTM2020, (Right) density ratios binned in local solar time of DTM2020. 

In the text 
Fig. 4 (Top) 1σ uncertainty in percent binned in latitude and the model fitting to data in the leastsquares (LS) sense, (bottom) 1σ uncertainty in percent binned in altitude, and the model fitting to data in the leastsquares sense. Data is in blue, and the LS fitting in red color. 

In the text 
Fig. 5 1σ uncertainty (in percent) leastsquare fitting: (topleft) 1σ uncertainty as a function of local time, (topright) 1σ uncertainty as a function of day of the year, (bottomleft) 1σ uncertainty as a function of solar flux, (bottomright) 1σ uncertainty as a function of K_{p}index. 

In the text 
Fig. 6 3D plot of the 1σ uncertainty (in percent) as a function of space variabilities: (left) altitude and latitude, (right) altitude, and local time. 

In the text 
Fig. 7 3D plot of the 1σ uncertainty (in percent) as a function of altitude and time evolution: (left) day of the year (right) years. 

In the text 
Fig. 8 Density ratios 1σ uncertainty binned in the day of the year of DTM2020 with a 27day averaged in red color, for two periods (2001–2003 and 2007–2009, respectively maximum and minimum of solar cycle 23). 

In the text 
Fig. 9 3D plot of the 1σ uncertainty (in percent) as a function of altitude and drivers : (left) solar flux, (right) magnetic index. 

In the text 
Fig. 10 3D surface plot of the 1σ uncertainty (in percent) (topleft) altitude and latitude, (left) geographical coordinates; (right) solar proxy and magnetic index. 

In the text 
Fig. 11 Contour level plot of the 1σ uncertainty (in percent) as a function of geographical coordinates: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (Top) altitude and latitude, (bottom) altitude and local time. 

In the text 
Fig. 12 Contour level plot of the 1σ uncertainty (in percent) as a function of temporal evolution: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (top) altitude and day of year, (bottom) altitude and years. 

In the text 
Fig. 13 Contour level plot of the 1σ uncertainty (in percent) as a function of solar and magnetic drivers: (left) Accelerometer uncertainty’s model, (right) HASDM uncertainty’s model. (Top) altitude and solar flux Fp, (bottom) altitude and magnetic index K_{p}. 

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.