A Globally Averaged Thermospheric Density Data Set Derived From Two‐Line Orbital Element Sets and Special Perturbations State Vectors

We describe a long‐term data set of global average thermospheric mass density derived from orbit data on ∼7,700 objects in low Earth orbit, via the effect of atmospheric drag. The data cover the years 1967–2019 and altitudes 250–575 km, and the temporal resolution is 3–4 days for most years. The data set is an extension and revision of a previous version. The most important change is the use of more precise orbit data: special perturbation state vectors are now used starting in 2001, instead of mean Keplerian orbital elements. The data are suitable for climatological studies of thermospheric variations and trends, and for space weather studies on time scales longer than 3 days.


Plain Language Summary
The region of Earth's atmosphere between about 90 and 500 km altitude, known as the thermosphere, exerts drag on the thousands of satellites and debris objects that orbit in this region. It is therefore important to understand and predict thermospheric variations in order to predict the future positions of these objects, for example, for collision avoidance and reentry prediction purposes. For the past five decades, the orbital trajectories of many of these objects have been routinely calculated from ground-based radar measurements. The orbits change in large part due to atmospheric drag, and this effect can be used to infer the density of the part of the atmosphere traversed by an object. This paper describes a long-term data set generated by applying the method to ∼7,700 objects to accurately estimate the global average atmospheric density as a function of time and altitude. The density data set covers the years 1967-2019 and altitudes 250-575 km, and it is suitable for climatological studies of thermospheric variations and trends, and for space weather studies on time scales longer than 3 days. The data set is a revision and extension of a previous version; the new version ingests more precise orbital trajectory data.
EMMERT ET AL.
Published 2021. This article is a U.S. Government work and is in the public domain in the USA. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. TLE format). For more information on this process and on orbit determination more generally, see Hejduk et al. (2013) and Vallado (2001), respectively.
Because the extrapolated SP ephemerides include the effects of forecasted thermospheric density, there is a strong potential for density derived from the resulting eGP elements to be contaminated by those forecasts. Therefore, we developed a new technique for deriving density from SP state vector data, which consist of the estimated position and velocity of an object at an epoch within the span of time ("fit span") covered by the observations used in the orbit determination. Such data have been operationally produced starting around 1999.
In this report, we describe a new 1967-2019 thermospheric density data set, based on TLE data through 2000 and SP state vectors starting in 2001. In the next section, we describe the new method for inferring density from SP state vectors. Section 3 describes the algorithm for producing the data set, which largely follows that of Emmert (2009, 2015a) (hereafter E09 andE15a, respectively). In Section 4, we present the new data set, compare with the current (E15a) data set, and give examples of scientific applications. Section 5 concludes the study.

Method for Deriving Density From SP State Vectors
SP state vectors are issued by the U.S. Space Force in the form of Vector Covariance Messages (VCMs) (National Research Council, 2012), which include the position and velocity (in both earth-centered inertial and earth-centered rotating frames) as well as the uncertainty covariance of the estimate and metadata such as the gravitational model used by the orbit propagator. From one state vector to the next, the change in translational kinetic energy per unit mass of a LEO object is the sum of the work done by gravity (mainly Earth's; here we also consider solar and lunar perturbations), atmospheric drag, and solar radiation pressure: Herein, all energy and work terms are expressed per unit mass. The work done by atmospheric drag between epochs t 1 and t 2 is (e.g., Picone et al., 2005): The work done by Earth's aspherical gravitational field is in general nonconservative because the field is rotating with respect to the Earth-centered inertial coordinate system. However, as derived in the supporting information (Text S1), it can be decomposed into a conservative potential term and a much smaller, path-dependent "rotational work" term that essentially consists of a torque applied by the rotating force field: The work done by solar and lunar gravitational perturbations is (Vallado, 2001): In our density estimates, we neglected the solar radiation pressure term. Below the 600 km maximum perigee height that we considered, we expect the net work done by solar radiation pressure to be small compared to the work done by drag, in part because drag always does negative work, whereas the positive and negative work done by radiation pressure tends to cancel out over the course of an orbit (cf., van den Ijssel et al., 2020, Figure 2 and p. 1766). Furthermore, as with the density data set derived by E09, we expect the large collection of objects used in our analysis to strongly reduce any individual-object errors that might be induced by the neglect of radiation pressure.
Inserting Equations 2 and 3 into Equation 1 and solving for the work done by drag, we obtain: Next, we define a ratio of the true weighted and integrated atmospheric density to a corresponding integrated model density: where we have introduced the superscripts T and M to denote true and modeled values, respectively. This ratio thus represents the ratio of the observed (as estimated from the state vectors) and expected (from an atmospheric density model) work by drag. Equation 6 is equivalent to Equation 22 of Picone et al. (2005), except that in the latter, the observed work by drag is computed from changes in the orbit-averaged elements contained in TLEs. Note that the orbit propagator used to produce LEO TLEs does not include azimuthal variations (i.e., sectoral and tesseral spherical harmonic components) in the geopotential field, nor does it include lunar or solar gravitational perturbations.
As a practical matter, the true ballistic coefficient B T is not known a priori for most objects. Therefore, it is convenient to define an atmospheric model-dependent ballistic coefficient B M as the value of B needed to reconcile the observed orbit changes with the atmospheric density model: The density ratio α can then be computed as once an estimate of the true ballistic coefficient is available (e.g., using the algorithm described in Section 3).
In our application of Equation 7, we calculate the kinetic energy change directly from the velocity vectors (in the Earth-centered inertial frame) of a pair of VCMs. The geopotential change ΔU is computed by evaluating the EGM-96 gravity model (Lemoine et al., 1998) at the position vectors (in the rotating frame) given in the two VCMs. EGM-96 is the gravity model used by the orbit propagator that produced the state vectors; in our calculation of EGM-96, we use the spherical harmonic truncation that is specified in the VCMs (typically degree and order 36 × 36 (National Research Council, 2012, p. 28). The calculation of the other four terms in Equation 7 (W rot , W S , W M , and the model density integral in the denominator) requires the evaluation of integrals along the trajectory of the object. For the trajectory, we use an ephemeris, at 1 min intervals, generated by the SGP4 propagator from the TLEs of a given object. An SGP4 ephemeris is computationally much faster and simpler than a more precise SP ephemeris derived from the VCMs, and we have found that the integral terms in Equation 7 are not sensitive to the errors in the SGP4 trajectories. This approach allows us to process a large number of objects while still exploiting the highly precise state vectors to calculate the changes in kinetic and potential energy, which are the dominant terms in the numerator of Equation 7 and thus the terms that require the most precision.
The method described above is a post-processing approach that applies conservation of energy to state vectors that were estimated via a separate orbit determination process. Similarly, the approach described by Picone et al. and used by E09 and E15a applies conservation of energy to existing orbital elements (which are essentially state vectors) contained in TLEs. In contrast, another approach is to tune drag accelerations as part of the orbit determination process; in this case, the orbit determination can be based on remote tracking observations (e.g., Storz et al., 2005;Tobiska et al., 2021), or on precise ephemerides derived from on-board GPS receivers (e.g., McLaughlin et al., 2011;van den IJssel, 2014). Note that the spatiotemporal resolution of both approaches is limited by the precision and cadence of the tracking observations or state vectors. In our case, the state vectors are derived using a typical fit span of 3 days, which limits the temporal resolution accordingly. Figure 1 illustrates the terms in Equation 5, for two example objects and for solar minimum and solar maximum years. The work terms were calculated using successive VCMs at least 3 days apart, and then normalized to the work done over 3 days (i.e., converted to average power with a time unit of 3 days). The kinetic and conservative potential (ΔU) terms are separately orders of magnitude larger than the other terms, so in Figure 1 they are combined (green traces). The expected drag term (left hand side of Equation 5; red traces) was computed using the NRLMSISE-00 density model (Picone et al., 2002), and the true ballistic coefficients of the two objects were estimated as described in Section 3. The amplitude of the rotational work term is ∼100 J/kg per 3 days, which is a non-negligible fraction of the work done by drag, especially at solar minimum and at the higher perigee altitude. The work done by solar and lunar gravity is typically less than 5 J/kg per 3 days, which even at solar minimum at the higher altitude is less than 5% of the work done by drag. The sum of all the terms on the right hand side of Equation 5 (black trace) closely follows the expected work done by drag (red trace), which gives confidence that we have accounted for all the major contributions to the kinetic energy change. The differences between the red and black traces can then be attributed, via Equation 6, to differences between the true and modeled atmospheric density. Figure S1 compares density ratios (α) computed from VCMs using Equations 7 and 8 with ratios computed from TLEs using the method of Picone et al. (2005), for four selected objects and two selected years. In cases where the drag is relatively strong (low perigee altitudes or solar maximum or large ballistic coefficient), the two methods agree well. However, the TLE results are noisier than the VCM results. In particular, at solar minimum for object 00063 (perigee ∼463 km), the VCM method produces density ratios that have a similar time dependence as two of the other objects (00060 and 00614), but the noise in the ratios produced by the TLE method is much larger than this signal. E09 found that such noise in the TLE-derived ratios could be greatly reduced by averaging over many objects. As demonstrated in the next section, combining VCM-derived density ratios from many objects produces a more precise estimate of the global average density ratio.

Data Processing Algorithm
Our production of a database of global average density largely followed the procedure of E09, except that beginning with the year 2001, we used VCMs and Equation 7 to calculate model-dependent ballistic coefficients (B M ). In this section, we outline the procedure and compare some interim processing quantities with the previous data set version (E15a).
Initial object selection was performed as described in E09 (Section 2.2), using TLEs for all years from 1967 to 2019. The purpose of this step is to remove objects that are maneuvering or whose orbital energy changes are dominated by sources other than atmospheric drag during a particular year. Although the eGP TLEs after about 2013 are not suitable for density calculations (see Section 1), we found that they are still reliable EMMERT ET AL.  and efficient for initial object screening. VCMs are not well suited for this purpose because they do not include the orbit-averaged mean motion that in TLEs effectively filter out aspherical gravitational effects. The reference set consists mainly of spherical or compact objects and includes the reference objects used in E09. This method requires a primary reference object whose B T value is known a priori. We used Starshine 1 (object 25769) for this purpose, with B T = 0.0105 m 2 /kg. This is the low end of the altitude-dependent range estimated by Pilinski et al. (2011), and is ∼8% larger than the value of 0.009744 m 2 /kg used by E09; the change is due in large part to the previous neglect of a launch vehicle adapter that made the satellite aspherical. As a result (via Equation 8) of the increased B T value, the new density values are smaller overall than those of E09 and E15a. We used the same set of reference objects used in E09, plus additional objects to cover recent years.
For the remaining objects, we estimated yearly B T values by applying the same method but pairing only with the reference objects, not with each other. Using the same screening procedure and criteria as E09, we removed objects with large temporal trends or interannual variability in their B T estimates. For each retained object, a single B T estimate was computed by averaging over the yearly values. Figure S2 shows yearly B T estimates as a function of time for a random selection of eight objects used in both the E15a data set and the new data set. In general, the old and new values agree closely and show the same interannual variability.
We then obtained a database of density ratios by applying Equation 8, retaining only those records for which the drag-weighted average altitude was between 150 and 600 km. Following Section 3.2 of E09, we computed an initial fit of the density ratios as a function of time and altitude, combining all remaining objects; we then applied a final quality control step by analyzing the residuals around the initial fit, removing outlying records and objects (using the criteria described in E09), and recomputing the fit.
The objects used in the final fit are tabulated in Data Set S1, along with their B T estimates, flags for the reference objects, and the number of observations used in each year. The total number of objects used in the final fit is 7,728, including 187 reference objects. Figure S3 shows the number of objects in the final fit as a function of year in five altitude bins (including all altitudes).
We computed estimated uncertainties of the time-altitude fit based on the object-to-object variance of the raw density ratios and the number of objects (i.e., the estimated error of the mean), following the procedure in Section 4 of E09. Figure S4 illustrates the density ratio variances (prior to removing outliers) of TLE-derived and VCM-derived density ratios, for selected years and altitude bins within the period of overlap (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) between the two methods. Especially under low drag conditions (higher altitudes and/or solar minimum), the object-to-object variance is smaller with the VCM-derived ratios, which therefore should provide a more precise estimate of atmospheric density than the TLE-derived ratios. There are numerous outliers in the VCM-derived ratios, but these were removed prior to computing the final time-altitude fit, following Section 3.2 of E09.
Finally, we applied the fitted density ratios and their uncertainties to gridded global averages of the reference model, NRLMSISE-00, following Section 5 of E09. The final density ratio and global average density products are described and discussed in the next section.

Results and Scientific Applications
Figure 2 compares final fitted density ratios (i.e., estimated true density divided by NRLMSISE-00 density) from E15a and the new VCM-derived ratios, for 2002 and 2008 (solar maximum and minimum, respectively) and at three selected altitudes. The new ratios are lower overall than the old values, as a result of the revised a priori B T value of the primary reference object, Starshine 1. The temporal dependence of the old and new values agrees very closely, except that the old values are noisier at the higher altitudes under solar minimum conditions, which is consistent with the object-to-object variance illustrated in Figure S4 and the typical density precision shown in Figure S5. Figure 3 shows annual average ratios of the new divided by old global average density values, at 250, 400, and 550 km. As expected, the new values are ∼8% lower overall due to the revised Starshine I B T value. However, the difference is larger (∼10%) at 250 km and smaller (∼4%-5%) at 550 km. This altitude-dependent sensitivity to the change in B T of the primary reference object possibly derives from the neglect of height-dependent changes in B T , which depends somewhat on atmospheric composition as well as temperature and orbital velocity (e.g., Doornbos, 2011;Pilinski et al., 2011) and is therefore expected to change with the decrease in perigee height over the orbital lifetime. The height dependence of the derived densities should therefore be interpreted with caution; however, we expect that the height dependence of variations in the density values are more reliable and more amenable to physical interpretation. Nonetheless, further development of the technique is warranted to estimate and mitigate possible long-term biases in the data set. For example, the calibration of the results could likely be improved by accounting for composition and temperature effects ( to apply to the non-reference objects, which are generally of unknown shape and orientation. Another possible improvement would be the application of theoretical B T values to multiple primary reference objects. At each altitude, Figure 3 indicates that the old and new results are temporally consistent to within 1%-2%. However, at 550 km, there is an upward shift of the new values (relative to the old values) of ∼3% from about 2004 to 2012, and there is a similar but smaller (∼1%) shift at 400 km. We believe that these shifts are due to the introduction of new objects, including new reference objects, for the years following 2007; the E15a data set used the same objects and B T estimates used in the 1967-2007 E09 data set. The anomalously low solar minimum period surrounding 2008 (e.g., Emmert et al., 2014;Solomon et al., 2011) may also have contributed. Whatever the cause, we expect that the new results are more accurate during this time period. Figure 4 summarizes the final global average log-density values as a function of time across all 53 years of the new data set. Daily, monthly running average, and yearly running average time series are shown for 250, 400, and 550 km altitudes. The data set now covers almost five full solar cycles. As with the E15a data set, the temporal resolution is 3-4 days.
The final fitted density ratios (relative to NRLMSISE-00) and their uncertainties are provided in Data Set S2 and in the repository listed in the acknowledgments, for ten altitudes from 250 to 575 km. The corresponding global average densities and uncertainties, computed as described in Section 5 of E09, are contained in Data Set S3. As illustrated in Figure S5, the precision of the new VCM-derived data (e.g., since 2001) is better than 1% above 400 km altitude and typically 1%-2% at lower altitudes. This precision is an improvement over the TLE-derived values, which in the 1990s is typically 1%-3% at all altitudes (some of the difference is due to the larger number of objects after 2000). We estimate that the overall accuracy of the data set remains 5%-10%, due to systematic uncertainty in the ballistic coefficient calibration.
Data Set S4 contains year-by-year plots of binned-average ratios compared to the final fitted ratios. Data Sets S5 and S6 are a similar set of plots showing the fitted ratios and their uncertainties (Data Set S5), and the log-density time series and their uncertainties (Data Set S6).
Orbit-derived mass density data are an important source of long-term information on the thermosphere and have been used in numerous scientific applications, including long-term trends (e.g., Emmert, 2015a), intra-annual variations (e.g., Weimer et al., 2018), global thermosphere response to stratospheric sudden warming events and the stratospheric quasi-biennial oscillation Liu, 2016;Yamazaki et al., 2015), tuning of empirical thermosphere models (e.g., Emmert et al., 2020;Jackson et al., 2020), the response of the thermosphere to solar extreme ultraviolet irradiance variations and coronal mass ejections (e.g., Dudok de Wit & Bruinsma 2017; Mendaza et al., 2017), and plasmasphere variations (e.g., Krall et al., 2016). Following up on Emmert et al., (2014), we are currently applying the data set to study the thermosphere's response to the current solar minimum, in comparison to previous solar minima, and we plan to update the long-term trend results of Emmert (2015a).

Conclusions
We derived a 53-year data set of global average, altitude-dependent thermospheric mass density from orbit data on ∼7,700 objects in low Earth orbit. The data set covers the years 1967-2019 and altitudes 250-575 km; the temporal resolution is typically 3-4 days (4-6 days during the first decade of the set). The data are available in the supporting information and at the repository listed in the acknowledgments.
The data set is a revision and extension of the one described in Emmert (2015a). The most important change to the methodology is that densities after the year 2000 are now derived from special perturbations state vectors (contained in vector covariance messages, or VCMs), instead of from Keplerian mean orbital elements (contained in two-line element sets, or TLEs, which are still used up to and including the year 2000). The state vectors are higher fidelity solutions to the orbit determination problem, and produce more precise estimates of thermospheric density. The new data set also incorporates a new set of reference objects, and it uses a higher a priori theoretical ballistic coefficient for the primary reference object (Starshine 1). As a result of the revised primary ballistic coefficient, the new densities are ∼8% lower overall than before. The difference between the new and old datasets depends on altitude, increasing from ∼4% at 550 km to ∼10% at 250 km.
The density data are suitable for climatological studies of global thermospheric variations and trends, and for space weather studies on time scales greater than 3-4 days. The height dependence of the absolute densities should be interpreted with caution, due to the assumption of constant ballistic coefficients, but we expect that the height dependence of density variations is reliable. Future improvements to the data set may include accounting for the effect of atmospheric composition and temperature on ballistic coefficients, as well as extending the data set further into the future.