A Simple Method for Correcting Empirical Model Densities During Geomagnetic Storms Using Satellite Orbit Data

Empirical models of the thermospheric density are routinely used to perform orbit maintenance, satellite collision avoidance, and estimate time and location of re‐entry for spacecraft. These models have characteristic errors in the thermospheric density below 10% during geomagnetic quiet time but are unable to reproduce the significant increase and subsequent recovery in the density observed during geomagnetic storms. Underestimation of the density during these conditions translates to errors in orbit propagation that reduce the accuracy of any resulting orbit predictions. These drawbacks risk the safety of astronauts and orbiting spacecraft and also limit understanding of the physics of thermospheric density enhancements. Numerous CubeSats with publicly available ephemeris in the form of two‐line element (TLEs) sets orbit in this region. We present the Multifaceted Optimization Algorithm (MOA), a method to estimate the thermospheric density by minimizing the error between a modeled trajectory and a set of TLEs. The algorithm first estimates a representative cross‐sectional area for several reference CubeSats during the quiet time 3 weeks prior to the storm, and then estimates modifications to the inputs of the NRLMSISE‐00 empirical density model in order to minimize the difference between the modeled and TLE‐provided semimajor axis of the CubeSats. For validation, the median value of the modifications across all CubeSats are applied along the Swarm spacecraft orbits. This results in orbit‐averaged empirical densities below 10% error in magnitude during a geomagnetic storm, compared to errors in excess of 25% for uncalibrated NLRMSISE‐00 when compared to Swarm GPS‐derived densities.


Background and Motivation
Earth's thermosphere is the region of the atmosphere between approximately 90 and 800 km, depending on solar conditions. Its middle and upper regions are the abode of numerous low-Earth orbiting (LEO) satellites that constitute billions of dollars in assets and are accompanied by some 20,000 currently known and trackable objects of space debris at least the size of a softball (Johnson, 1993). Understanding the behavior of the density of this region is vital to being able to accurately predict the orbits of these objects, as the amount of drag they experience is contingent on the magnitude of the local neutral density. Empirical models of the thermosphere like Jacchia-1972(Jacchia, 1979, DTM-2012, and NLRMSISE-00 (henceforth MSISE-00) (Picone et al., 2002) are used by the space situational awareness community to estimate the thermospheric density for orbit prediction. These models commonly exhibit errors in the density in excess of 20% during periods of high geomagnetic activity (Burke et al., 2007;. Density errors translate directly into orbit errors, jeopardizing the success of collision avoidance, re-entry prediction, and spacecraft maneuver planning (Bussy-Virat, Ridley, & Getchius, 2018; Doornbos & Klinkrad, 2006).
It is understood that the thermosphere's state is highly contingent on the energy input it receives in the form of solar extreme ultraviolet (EUV), Joule/frictional heating, and auroral precipitation, which all serve to control the thermospheric temperature. During geomagnetic quiet times, incoming solar EUV constitutes the largest part of the thermosphere's energy budget, but during a geomagnetic storm, up to two thirds of the energy budget can be comprised by Joule/frictional heating and auroral precipitation (Knipp et al., 2004). The direct dependence of the neutral density on the temperature results in the dynamics of the thermospheric density being influenced by diurnal tides, solar rotation, and the solar cycle (Forbes et al., 2012;Rhoden et al., 2000;Ruan et al., 2015;Vickers et al., 2014).
The second largest source of energy into the thermosphere is high-latitude Joule (or frictional) heating. This is due to collisions between ions and neutrals, since these populations have differing bulk velocities and temperatures. During quiet times, the Joule heating is relatively small compared to the solar EUV, but during geomagnetic storms, the heating at high latitudes can become the dominant source of energy into the thermosphere, causing up to an 800% increase in the density as up to 1,000 GW of energy is deposited during a storm Sutton et al., 2009;Vichare & Lakhina, 2005). Variations in the thermospheric density due to thermal expansion from EUV radiation and Joule heating affect orbiting satellites by changingthe degree of atmospheric drag they experience as energy from high latitudes is distributed globally via waves and pressure/temperature-gradient driven winds within hours (Burns et al., 1995;Mayr et al., 1984;Prolss, 1993;Zesta & Oliveria, 2019). There is thus a correlation between the rate of change in the semimajor axis of satellite orbits and geomagnetic activity, which can be visualized by observing an increase in spacecraft altitude decay during geomagnetic storms. This behavior has been investigated with the CHallenge Mini-satellite Payload (CHAMP) and the Gravity Recovery And Climate Experiment (GRACE) satellites, where it was found that storm time drag effects manifest first at high latitudes, that for extreme storms stronger orbital drag effects occur during the early main phase of the storm at low/equatorial latitudes, and that storm time orbital decay is directly proportional to storm intensity (Oliveira & Zesta, 2019). Drag modeling of both the CHAMP spacecraft and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite has also revealed that during geomagnetic storms, a 60% increase in orbit decay rate can occur due to the increased drag caused by thermospheric heating (Nwankwo et al., 2015). Figure 1 shows the relationship between geomagnetic activity and orbit decay with the strongly positive correlation between the rate of de-orbit of 20 identical 3U satellites launched by Planet Labs, Inc. and geomagnetic activity represented by the strength of the Earth's ring current (Dst). These Earth-observing 5-kg satellites are all sun-synchronous at altitudes in the vicinity of 450 km. The rate of change of altitude for these satellites was determined from two-line elements (TLEs) and averaged across all satellites in 6-hr windows. It was then compared to the 6-hr average of Kyoto Dst for the same time period and shifted back in time by 6 hr to account for the delay in behavior between Dst and orbital decay (Figure 1a). This yielded a peak correlation coefficient of 0.76 (Figure 1b).
For the majority of small satellites deployed in the thermosphere, records of their semimajor axes are available in the form of TLEs, a data product provided to the public by the the U.S. Air Force's 18th Space Control Squadron, located at the Combined Space Operations Center (CSpOC) at Vandenburg Air Force Base in California. They provide the mean orbital elements of a spacecraft at specific time and are typically Figure 1. The rate of change of the semimajor axis per year of 20 identical Flock 2 K satellites overplotted with the Dst during a geomagnetic storm (a). Both the semimajor axis per year and the Dst have been averaged in 6-hr time windows, and the value of the semimajor axis has been shifted forward in time by 6 hr to account for the characteristic time delay between initial storm onset and the resulting change in spacecraft altitude. Shown in panel (b) is the positive correlation between the rate of change of the semimajor axis and of Dst for the same time period.
reported once a day for LEO spacecraft. The orbit determination method of differential correction (DC) is used to generate TLEs and is essentially a multidimensional Newton-Raphson root solving method of y ¼ f ðxÞ with a least-squares statistical treatment of the known data ( y) provided by observations, either via GPS or visual/radar tracking (Vallado & Crawford, 2008). TLEs were designed to be used expressly for the purposes of orbit prediction with specific models, of which the most common is a set of Simplified General Perturbation (SGP) models referred to as SGP4 (Vallado et al., 2006). SGP4 is based off of theories of satellite motion described by Kozai (1959), Brouwer (1959, and Lyddane (1963), all of which neglected the effects of drag; SGP4 accounts for drag via power density functions (Hoots et al., 2004;Lane & Cranford, 1969;Lane & Hoots, 1979) that require a term that encapsulates the ballistic coefficient, B * , which can be used in the expression for the acceleration due to drag ; (1) where ρ is the local thermosphere density, ρ 0 is a reference air density given as 0.1570 kg/m 2 /R E , B is the ballistic coefficient in units of area per mass, v is the velocity of the object, C D is the drag coefficient, A is the cross-sectional area of the object as viewed from the ram direction, m is the mass of the object, and v m is the velocity of the medium through which the object is traveling. For LEO spacecraft, v m is representative only of the rotation speed of the atmosphere and usually neglects thermospheric winds. Unfortunately, the drag coefficient is often treated as constant for the given TLE, which is typically provided every 24 hr. This fails to capture how the changing composition of the thermosphere with altitude affects the gas-surface interactions on the spacecraft faces, leading to a variable drag coefficient. The Spacecraft Orbital Characterization Kit (SpOCK), an orbital propagator developed at the University of Michigan, does not rely on B * when using TLEs to perform orbit prediction . Instead, it allows the user to describe the spacecraft geometry with a CAD file or describe the orientation and area of all the spacecraft faces manually. It also restricts itself to using the mean orbital elements from the TLE, combined with an accommodation coefficient (α) for the spacecraft surfaces, permitting the calculation of a variable drag coefficient . It currently relies on the MSISE-00 empirical model for thermospheric density estimation.
Despite taking into account the spatial and temporal variations of the drag, modeled altitudes, specifically during geomagnetic storms, underpredict the semimajor axis decline (Figure 2). Immediately after the storm onset, the spacecraft's measured rate of decay increased, but SpOCK, relying on MSISE-00 for density estimation, failed to reproduce that increase. This shortcoming is partially attributable to the limitations of empirical models, which rely on parametric fits to a catalog of density measurements taken from a variety of sources, including sounding rockets and accelerometer data from spacecraft (Picone et al., 2002). During quiet times, these models typically exhibit density errors on the order of 10% (Bruinsma et al., 2014;Picone et al., 2002;Vallado & Finkleman, 2014). These uncertainties can increase to beyond 40% during geomagnetic storms for the reason that periods of intense geomagnetic activity are relatively infrequent and unpredictable, and their dynamical effects on the thermosphere are not understood well enough to provide for reliable density predictions over relatively short time scales (Eyiguler et al., 2019). These limitations are exacerbated by biases in TLEs themselves; as the orbital elements encoded in TLEs constitute mean Brouwer-Lyddane elements, the process of their calculation smooths out short-periodic affects in the elements that repeat on the order of a satellite's orbital period (Vallado et al., 2006). This study presents a new technique that addresses the density underestimation using TLEs from different satellites. The MSISE-00 density inaccuracy can be probed by using SpOCK to reproduce satellite trajectories during both periods of quiet and active geomagnetic activity. As SpOCK relies on both MSISE-00 and the TLEs for initialization and propagation, estimation of calibration factors to the geomagnetic inputs to MSISE-00 can yield a calibration method for thermospheric density models. A similar method has been shown to be successful by Doornbos et al. (2008), which involved the conversion of TLE data to drag data used in the daily adjustment of density model calibration parameters. This method found its inspiration in the Dynamics Calibration Atmosphere (DCA) used by the USAF's High-Altitude Satellite Drag Model (HASDM), which uses Space Surveillance Network observations of ∼75 orbiting spheres at various altitudes, performing a least squares DC across all satellites to solve for global density corrections to an empirical density model (Storz et al., 2005). Doornbos et al. (2008) relied on the techniques in Picone et al. (2005) to calculate the density directly from individual TLEs before performing a least-squares adjustment to minimize the difference between TLE-derived densities and those of empirical models. They calculated ballistic coefficients under the assumption that the long-term average of the ratio of observed over modeled density ratios should equal unity. The DCA in HASDM uses calibration satellites to solve for a global density correction field in terms of a set of coefficients to a local low-altitude and high-altitude temperature parameter. The DCA estimates a new set of zeroth degree coefficients for the low-altitude parameter every 3 hr, while the others are estimated every 12 hr (Storz et al., 2002). Gondelach and Linares (2020) demonstrated the power of a reduced-order model (ROM) that combines the predictive capabilities of physics-based models with the computational speed of empirical models in estimating the global density. Their ROM is a full data assimilation scheme that is capable of using a variety of data sets, including accelerometer-derived densities, nonlinear space weather model inputs, modified equinoctial elements that describe satellite orbits, and TLEs. A dynamic model is derived that retains the primary characteristics of the state space describing the thermospheric density but at lower dimensionality. The Dragster model developed and presented by Pilinski et al. (2016) operates by modifying the drag effects on spacecraft to minimize orbit error and uses three full-physics atmospheric models in tandem in an assimilative architecture. It also uses an ensemble Kalman filter to estimate solar and geomagnetic forcing parameters and density corrections (see Figure 3 in Pilinski et al. (2016) for more details). Dragster draws from a wide variety of high-latitude forcings, includinga p measurements, Defense Meteorological Satellite Program (DMSP) data, Super Dual Auroral Radar Network (SuperDARN) data, and Active Magnetosphere and Planetary Electrodynamics Response Experiment (AMPERE) data, as well as a wide variety of solar forcings that include F10.7, Solar Radiation Physical Model (SRPM) inputs, and Flare Irradiance Spectral Model (FISM) inputs. Dragster produces nowcast and 3-day forecasts of density, composition, and winds along satellite orbits in order to compute drag estimates.
While there are many different methods for computing the thermospheric density, we present a new one, the Multifaceted Optimization Algorithm (MOA). This method estimates the neutral density by minimizing the orbit error between several modeled spacecraft and a set of TLEs. MOA first estimates a representative cross-sectional area for each spacecraft during the geomagnetically quiet 3 weeks preceding a storm and then estimates modifications to F10.7 and a p inputs to MSISE-00 in order to minimize the RMS error between the modeled and TLE-provided semimajor axis of each spacecraft. The mean values of the modified F10.7 and a p values are then used to obtain corrected densities by using them to drive MSISE-00. Unlike Doornbos et al. (2008), MOA doesn't adjust temperature or scale factors to the density to generate new density estimations, is not validated by comparing its densities to TLE-derived densities, and it focuses on minimizing the orbit error between altitude changes from SpOCK and those from TLEs, by estimating modifications to geomagnetic indices that are inputs to MSISE-00. In contrast to Storz et al. (2002), and much more simply, MOA does not require as many calibration satellite to function and, instead of solving for a global density correction field in terms of coefficients related to temperatures, solves for a time series of median modifications to F10.7 and a p along the time period of interest. These modifications are then applied to their respective inputs to MSISE-00 to yield corrected densities. In contrast to Gondelach and Linares (2020), the method described here demonstrates a simple way of improving density modeling that relies on less information and less processing. This method requires only spacecraft TLEs from a small number of objects, an orbit propagator, and geomagnetic indices, and it demonstrates how improvements to storm-time density modeling can be achieved with limited information. In contrast to Pilinski et al. (2016), MOA does not use a Kalman filter, is not assimilative, only relies on orbit error minimization to estimate modifications to only two solar/geomagnetic model inputs to a single atmospheric model, and is primarily focused on improving storm-time modeling of neutral density. MOA is distinguished from Emmert et al. (2008) in that the technique presented uses a binary search algorithm to determine satellite cross-sectional areas before modifying geomagnetic indices, rather than following the procedures of Emmert et al. (2006) to calculate ballistic coefficients, and it does not determine corrected densities by first calculating density ratios and applying them to MSISE-00 densities. We compute modifications to MSISE-00 densities during the May 2017 geomagnetic storm across ten identical 3U CubeSats of Planet Labs, Inc. and apply these adjustments along the Swarm spacecrafts' orbits in order to directly compare the results with GPS-derived measurements of the thermospheric density.

Flock 3P and Swarm
The Flock 3P spacecraft represent the 10th "Flock" of Dove Earth Observation CubeSats launched by Planet Labs, Inc., an imaging company based in San Francisco, CA. Each CubeSat is approximately 5 kg in weight, conforms to the 3U form factor, and is equipped with a 9-mm diameter telescopic imager that can collect panchromatic, color, and near-infrared imagery of the Earth at a resolution of up to 3 m (see Foster et al., 2015& Safyan, 2015 for details). Each satellite has extendable solar panels that deploy from opposite sides, resulting in a maximum projected area of 1,950 cm 2 in one plane ( Figure 3). The satellites have a deployable cover over the telescopic imager and are nadir pointing, with the imager pointing directly downward to the Earth's surface. Flock 3P consisted of a total of 88 satellites, launched on 15 February 2017 into a sun-synchronous orbit at an altitude of ∼500 km by the Indian Space Research Organization's (ISRO) 39th flight of its Polar Satellite Launch Vehicle (PSLV) on the PSLV-C37 mission. The first 10 of these CubeSats were selected as calibration targets for use by MOA: the TLEs of the first 10 of these CubeSats were obtained from Spacetrack and used to generate adjustments to solar and geomagnetic model inputs.
Swarm is a European Space Agency (ESA) mission aimed at studying Earth's magnetic field. It consists of three satellites (A, B, and C) placed in two different orbits, two (Swarm-A and Swarm-C) flying alongside each other at 450 km in altitude at 87.4°inclination and the third (Swarm-B) at an altitude of 530 km at 88°inclination (Dunlop et al., 2015). The Swarm spacecraft, built largely by the former satellite manufacturer, Astrium, are all identical and carry onboard a vector field magnetometer, absolute scalar magnetometer, electric field instrument, accelerometer, and laser range reflector. The Swarm spacecraft were launched on 22 November 2013 from Plesetsk Cosmodrome in Russia on a Rockot vehicle and offer coverage through the auroral regions and across the high-latitude polar cap at a variety of local times similar to that of the CHAMP spacecraft (Reigber et al., 2002). GPS tracking of the Swarm spacecraft has allowed for the generation of a data set of GPS-derived thermospheric densities along the spacecraft orbits that were generated with the Precise Orbit Determination (POD) strategy, which converts range and phase information in the GPS measurements into nongravitational accelerations (van den Ijssel et al., 2016(van den Ijssel et al., , 2020. We compared MOA's densities to Swarm's GPS-derived densities in order to assess the quality of the former's performance.

SpOCK
SpOCK is an orbital propagator that simulates spacecraft's location given a series of inputs that may either be entirely user-supplied or provided by various scientific databases. SpOCK comprises a suite of C functions that require the user to supply a geometry file and a main input file. The geometry file describes each face of the spacecraft, including the unit vector describing the orientation, the surface area, the total surface area of any solar cells on that surface, drag coefficient or accommodation coefficient, and the solar radiation coefficient. SpOCK requires specification of the solar irradiance, as proxied by F10.7, and the planetary activity level, as specified by 3-hr a p . Both of these geomagnetic indices are available through either the NASA
Once SpOCK is commanded to be run with the appropriate initialization information, it obtains an estimation of the local spacecraft density using MSISE-00. After calculating an estimate of the local drag and other perturbing forces, such as higher order gravity terms, gravity due to the Sun and Moon, solar pressure, and albedo effects from sunlight reflecting off of the Earth, SpOCK then propagates the trajectory of the spacecraft for a single time step specified by the user in the main input file. SpOCK repeats this process until the given stopping point is reached. A large problem with techniques such as this is that the ballistic coefficient of the object is typically not known, for example. Objects such as collisional debris are highly uncertain in terms of their geometry, mass, drag coefficient, or orientation. Objects such as CubeSats, defunct satellites, or rocket stages may have information on the geometry and mass but may not have a good estimate on the drag coefficient or the orientation. Satellites such as CHAMP, GRACE, and GOCE, as well as the reference spheres, have a good estimate of all of the aforementioned characteristics, but this is quite rare. The ballistic coefficient can be derived from B * in the TLE, but this is specifically designed for SGP4, so it is ignored. Instead, MOA uses a series of TLEs during quiet geomagnetic conditions to estimate the surface area of the object, assuming a mass and accommodation coefficient, from which a drag coefficient is derived, following the work of Picone et al. (2005). This quiet period is set as the 3-week period before the time period being considered.

MOA
MOA first collects TLEs for a specific satellite for a user-specified interval of time. It uses the first TLE to initialize SpOCK. Subsequent TLEs are used to compare to the SpOCK semimajor axis predictions in an attempt to reproduce the orbital profile. MOA approximates the satellite geometry as a flat plate of an estimated cross-sectional area and the object's known mass. MOA uses this basic framework in three processes that use TLEs from selected satellites to obtain corrected model densities (Figure 4).
Orbit error minimization is performed using the root mean square error (RMSE). The RMSE (δ z ) between SpOCK altitude values (ξ i ) and the corresponding time steps closest to those of the TLE altitudes in the interval of the run (i ¼ 0:T) in question is calculated: ( 2) where N is the number of TLEs within the optimization interval. This method minimizes the RMSE between the SpOCK altitudes and TLE altitudes throughout the entire profile within the interval, using a binary search algorithm. This difference provides a more robust comparison as opposed to differencing the ending altitude.
The first process is the Area Optimization Algorithm (AROPT), which is a binary search algorithm that orbits a flat plate for a geomagnetic quiet period of 3 weeks prior to the user-defined time interval, varying the area of the flat plate in each iteration, searching for the orbit trajectory that best matches the behavior of the altitude specified in the series of TLEs. AROPT first computes the orbit error for the upper boundary (A U ), then the lower boundary (A L ), and finally the mean of both (A M ). These initial runs allow it to decide between which values the optimized area lies. The limits of the area search algorithm are twice the maximum possible projected area and half the minimum possible projected area for the given spacecraft. The projected area estimates are not obtained from TLE information but from the user-defined geometry file that specifies the orientation, material, and area of each face. The real unknown for the orbiting object is the ballistic coefficient, which is C D A m , but a mass and accommodation coefficient are approximated, implying that the real area  and derived area may be different, depending on how far the C D and mass estimates are off. The algorithm iteratively finds an area that allows the ballistic coefficient to minimize the error in the drag over the course of 2-3 days. This process assumes that the projected area is constant over the 2-3-day time period. This is most likely inaccurate, but it is permissible if the behavior of the object is repeating much faster than the minimization time period. For example, a tumbling object can be modeled with an average projected area, since it is most likely tumbling much more quickly than the optimization period. This assumption will fail if the object is systematically changing attitude for long periods of time. For example, the Cyclone Global Navigation Satellites (CYGNSS) satellites switch to a high drag mode for several days at a time to reduce the semimajor axis. This optimization scheme would work within the low-drag time and the high-drag time but would not come up with a proper area during the transition. The sun-pointed satellites used in this study have areas that change through an orbit, but that change is rapid compared to the minimization time-period (∼24 hr), allowing an average area to be deduced. The accommodation coefficient was assumed to be 0.9 for all spacecraft surfaces. This was used to calculate a variable drag coefficient using Equations 2 and 3 from Moe et al. (2004). The mass of each satellite was kept at 5 kg, which has been the case for every FLOCK spacecraft launched so far (eoPortal, 2017). MOA runs AROPT during the quiet three weeks prior to the storm, takes the 75th percentile of the areas it returns, and holds that area constant throughout the rest of the subprocesses. Areas returned by AROPT during this time period preserve the quiet-time bias of the empirical model SpOCK is using.
A single run of AROPT is shown in Figure 5 for one of the 3U "Flock 3P" satellites built by Planet Labs, Inc. The solid magenta line shows the altitude of the spacecraft as given by its TLEs, while the dashed lines show altitudes for the spacecraft simulated by SpOCK for different cross-sectional areas. The first three of those dashed lines indicate the lower, upper, and central boundaries used by AROPT, while the following six dashed lines are the results of AROPT running for different projected areas between the central and upper boundaries.
The same time period was run multiple times using a binary search algorithm until the RMS error was less than 1 m. This general search for the optimum area typically took around six iterations. The area optimization was conducted for subsequent intervals of time, creating a time series of optimized areas. After the completion of these runs, a histogram of the optimized areas was generated, along with a graph showing the behavior of the optimized areas over the entire period. A histogram of optimized areas for this spacecraftover a 3-week period is shown in Figure 6, where the dashed magenta line represents the 75th-percentile optimized area found by AROPT.
Using AROPT as the first step in MOA allows an inference of the orientation of the chosen spacecraft. For a spacecraft with known geometry, comparing the optimized areas returned by AROPT to the projected areas of the spacecraft from each of its sides allows for a coarse idea of its average orientation to be drawn. AROPT therefore represents a method of accounting for spacecraft variable geometry (deployable panels and antennae) and changes in attitude. AROPT selected a specified percentile from the distribution function of areas as the constant area to be used by the remaining subprocess(es) (e.g., in Figure 6, the 75th percentile is 817.83 cm 2 ). The rationale behind the selection of a specific quartile assumes the rate of de-orbit of the satellite in question will be overwhelmingly attributable to changes in the space environment captured by the behavior of geomagnetic indices. The optimized areas found by AROPT will be inextricably tied to the empirical model from which SpOCK receives an estimate of the density, as the finding of these areas essentially compensates for bias in the model. The optimized area obtained by AROPT will add bias in the downstream predictions, since it will be the static area for the F10.7 and 3-hr a p optimization. The goal of the area optimization is to allow a static area to be chosen, so that the thermospheric density can be altered away from the MSISE-00 predicted values. This assumes that, during the 3-week quiet interval selected for the area optimization, MSISE-00 predicted the correct mass density on average, but the values at time scales smaller than a couple of weeks may be incorrect. It is emphasized that variation in area ( Figure 6) could be due to (1) discrepancies between MSISE-00 and reality; (2) errors in TLEs; or (3) issues with F10.7 and a p describing the thermosphere during the optimization interval. This technique assumes that, on average, MSISE-00, the TLEs, and the models used to generate F10.7 and a p are unbiased such that the average of each of the those errors over a long time (over the ∼3-week long period being considered) is minimal. If this is not the case, then the bias extends into the storm period, so that the area bias compensates for the model bias.
It is important to emphasize that AROPT is not designed for the purpose of high-fidelity spacecraft geometry modeling. Studies that have focused on generating data sets of accelerometer-derived densities have had to contend with high-fidelity geometry modeling in order to minimize error in the generated densities. This is done to allow for data that has a temporal resolution of much less than an orbital period, even down to a few seconds. A major obstacle in achieving this resolution lies in uncertainty in the ballistic coefficient (Bhatnagar et al., 2005;Cook, 1965;McLaughlin et al., 2011), which can be quite accurately obtained from accelerometer data (Mehta et al., 2017) combined with models and algorithms to predict the drag coefficient using response surface modeling (Mehta et al., 2014), by even adjusting atomic oxygen accommodation (Bernstein et al., 2020), or derived directly from TLEs (Sang et al., 2013). MOA distinguishes itself from works like Mehta et al. (2017) in that it is not an algorithm that directly obtains densities using estimations of the ballistic coefficient. For each satellite, AROPT minimizes the error between the SpOCK-predicted semimajor axis and TLE-specified semimajor axis using an optimized cross-sectional area that changes every 2 days throughout the entire interval. The 75th percentile of all of the resulting optimized cross-sectional areas is then used by the resulting subprocesses. AROPT serves only as a means of coarsely removing contributions to changes in spacecraft altitude due to the drag before MOA estimates new densities by focusing entirely on contributions from geomagnetic indices. AROPT can be used for objects that either lack or don't have good CAD models or for objects that are rotating, since it does not attempt to have high temporal resolution. Instead, it focuses on longer time scales for individual objects. In short, it takes a group of objects that have some uncertainty and averages them together to obtain an acceptable result. It is a method that has significant error with only one object but reduces that error with many objects: 10, 20, or more.
While it is not the main purpose of MOA to obtain high-fidelity estimates of the satellites, the AROPT subprocess is capable of approximating the cross-sectional area of the satellites in a way that can account for variable orientation before estimating modifications to the geomagnetic indices, assuming that MSISE-00 is mostly unbiased during the time period of interest. Figure 7 demonstrates this for one of the CYGNSS spacecraft (Ruf et al., 2013). When in their nominal nadir pointing orientation, the projected area is 1,190.35 cm 2 , but occasionally, the spacecraft enter high-drag mode, when the projected area increases to approximately 7,784.58 cm 2 in order to more rapidly reduce their semimajor axis to reduce the relative velocity between the eight spacecraft. Figure 7 shows a transition from nominal nadir orientation on 27-28 July to a high-drag orientation for the rest of the time. AROPT was able to coarsely show the transition to high-drag by calculating a rise in the cross-sectional area above 6,000 cm 2 during the time of interest. While this was ultimately lower than the actual projected area, it shows that AROPT can retrieve approximate areas. In addition, this is assuming that MSISE-00 is completely unbiased during this time, and the difference in the areas could be due to such a bias. If the ultimate goal of MOA was to estimate the area perfectly, this would not be ideal, but because the end result is the density, any model bias in MSISE-00 is compensated for in the area optimization so that the resulting density would have significantly less bias.
Once a ballistic coefficient for an object is derived, the thermospheric density corrections can be derived. This is done across two time scales with two different indices: the F10.7 to represent ∼24-hr density adjustments (FOPT) and a p to represent 3-hr density adjustments (APOPT). By altering these indices, the global density can be altered, as opposed to the localized density. These alterations on the indices don't imply that the indices themselves are incorrect but are a way to alter the MSISE-00 densities without actually altering the MSISE-00 source code or using an arbitrary multiplicative factor. The FOPT and APOPT processes are binary search algorithms like AROPT, except they alter F10.7 and 3-hr a p , respectively. In order to set upper and lower boundaries, however, FOPT and APOPT set limitations to maximum and minimum values. For FOPT, the F10.7 maximum is set to 200 sfu, the standard reference value for F10.7 at solar maximum in literature (Husler et al., 2010;Tapping & DeTracey, 1990;Zhou et al., 2016) and the minimum set to 80% of the value of the smallest F10.7 value in the interval. The 80% value varied as the F10.7 for the days considered varied between 62 and 67 sfu. For APOPT, the 3-hr a p maximum is set to 400, the maximum possible, and the minimum to 0, the widest range of boundaries defined for this index (NCEI, 1999). MOA subsequently brackets between those limits to find the correction necessary to minimize orbit error. FOPT and APOPT do not apply adjustments to the solar and geomagnetic indices in a multiplicative manner but in an additive manner. The FOPT process is run after AROPT is complete.
APOPT only runs if MOA determines that a geomagnetic storm has occurred during the specified interval. It does this using data from the World Data Center for Geomagnetism and considers any geomagnetic disturbance in which Dst passes below −50 nT to be a storm (Syun-Ichi, 2018). In such a case, it will run APOPT for the 2 days following the date of initial storm onset, in order to account for the duration over which the impact of geomagnetic activity is felt globally. APOPT runs while holding the most recent correction to F10.7 obtained via FOPT constant throughout the duration of the storm. This is done for the reason that F10.7 not only varies on much longer time scales on the order of days (e.g., Wang et al., 2018) than 3-hr a p , which varies on the order of hours (e.g. Wrenn, 1987), but it also is only collected once a day. Therefore, any rapid changes in the density during stormtime will most strongly correlate with fluctuations in 3-hr a p and almost not at all with F10.7. Base values of F10.7 and a p data are taken either from NASA's Space Physics Data Facility OMNIWeb service or SWPC. This data source selection is held constant throughout all of the subprocesses.

Results
In order to demonstrate the MOA technique, the time period of 23 May 2017 to 6 June 2017 was explored and is presented here. During this time, 10 Flock 3P spacecraft launched by Planet Labs, Inc. were on orbit (Table 1), and a geomagnetic storm of moderate strength took place (Figure 8). The FLOCK 3P were in sun-synchronous orbits between 490 and 500 km in altitude. Multiple satellites were used to optimize the density corrections, which were then compared to GPS-derived densities from the Swarm spacecrafts. . TLE-derived rate of change of the semimajor axis for the Flock 3P satellites is in black (c). The rates of change for Flock 3P-1 (red) and Flock 3P-6 (blue) display unique behavior.
The same geometry file was used for each of these satellites: All of the CubeSats of the "Flock" series launched by Planet Labs, Inc. are identical and can be described as ordinary 3U CubeSats with deployed solar panels normal to the zenith direction possessing an area of approximately 1,950 cm 2 (Foster et al., 2015). The mass of each spacecraft was fixed at 5 kg, which has been the mass of each Flock spacecraft launched so far (eoPortal, 2017;Planet Labs, 2015). As they were Earth-imaging and sun-synchronous, their solar panels were always angled toward the Sun, and as a result the projected cross-sectional area was expected to vary around 1,000 cm 2 .
As described above, F10.7 adjustments are strongly influenced by the percentile of the optimized area distribution chosen. The adjustments to F10.7 in and of themselves are merely a means of correcting the density and are not the object of this study. We elected to compute adjustments using the 25th, 50th, and 75th percentiles of the optimized area distribution in order to determine its effect on the adjusted densities. For each of the 10 satellites, the following procedure was conducted: • The area distribution function was determined over the selected time period. • F10.7 adjustments were then calculated over that interval.
• 3-hr a p adjustments were finally calculated during the main phase of the storm.
Once this was completed, median daily F10.7 values were generated using results from all 10 satellites. These medians were calculated from constant F10.7 adjustments in each 24-hr interval and associated with noon of their respective days. Finally, median 3-hr a p adjustments were generated during the storm main phase. To validate the efficacy of the method, MSISE-00 was run at the locations of the three Swarm satellites using the modified median F10.7 and a p values. We focus on three major metrics to perform validation:  1. δ P : Percent difference between the peak orbit-averaged density between uncorrected MSISE-00 orbit-averaged densities, MOA orbit-averaged corrected densities, and orbit-averaged Swarm GPS-derived density data: where ρ N is either the MSISE-00 or MOA orbit-averaged density and ρ S is the Swarm GPS-derived orbit-averaged density. 2. η: Ratio of the peak orbit-averaged density magnitude to the 24-hr averaged orbit-averaged density prior to the peak density within the 24 hr immediately preceding the peak density. 3. ρ T : Total time-integrated density in kg · s km 3 during the main phase of the storm. In order to set the boundaries for calculating this integral, the following was done: (a) For each point, the arithmetic mean density and standard deviation of the density for the Swarm orbit-averaged density is calculated for the preceding 12 hr. (b) The lower bound of the integral is found when the density exceeds the mean + standard deviation at its associated time, and all of the density values for the next 12 hr satisfy that condition. (c) The upper bound of the integral is found using the same method but going backward from the density values at the end of the chosen time period. 4. t l : The time difference in hours between the peak in the MSISE-00 or MOA orbit-averaged densities and the peak in the Swarm GPS-derived orbit-averaged densities.
The two Swarm spacecraft chosen allow us to determine the utility of our adjustments to global geomagnetic indices along orbits at different altitudes: The Swarm-A spacecraft orbits at an altitude of ∼460 km and inclination of 87.4°, while Swarm-B orbits higher at ∼530 km in altitude and an inclination of 88°.

Optimized Areas and Adjustments to F10.7 and 3-hr a p
The geomagnetic indices during and surrounding the 27-28 May 2017 storm are shown in Figure 8. During this time, 3-hr a p surged sharply, peaking at ∼130 nT on 28 May, while Dst reached a minimum of −125 nT. This increased a p was associated with the density change that perturbed the orbits of the spacecraftand affected the baseline densities on which adjustments were performed. F10.7, in comparison, exhibited negligible variation.
The rate of change of the Flock 3P satellites' semimajor axis derived from the TLEs are shown in Figure 8c. A clear increase between 28 and 30 May corresponded to the main phase of the storm. The average maximum rate of change attained by the constellation during the storm main phase was over 20 km/year. This maximum occurred between 29 and 30 May, even though the storm's peak intensity, as indicated by 3-hr a p , occurred on the 28th. This delay in behavior is likely due both to the lower time resolution of TLEs and the fact that changes in the local density outside of high latitudes do not occur immediately in response to geomagnetic activity, but often hours later (Guo et al., 2010;Oliveira et al., 2017). There is a slight drop in the rate of change of the semimajor axis (dSMA) around 26 May just prior to the main phase of the storm. This drop is also observable in of all of the Flock 3P satellites, which may suggest a persistent structure in the thermospheric density response that has a consistent effect on all of the satellites. This indicates that the thermospheric density was most likely lower than expected on this day. It should be understood Figure 11. Linearly interpolated median F10.7 adjustments corresponding to the MOA runs where the 25th (black), 50th (blue), and 75th (red) percentiles were used. These results were from runs driven by OMNIWeb inputs. that the change in SMA is due to the integral of the density in the time prior to the measurement. These changes do not directly reflect changes in the instantaneous density.
The optimized area distributions found for each satellite clustered around 800 cm 2 (Figure 9). If 0°is considered to be parallel with the direction of travel, these results suggest that the Flock 3P satellites' largest face consisting of the solar panels was at an orbit-averaged angle of ∼20°. The high degree of overlap of the distributions suggests the orientations of the spacecraft were very similar during the time period. The notable exception was Flock 3P 1, which may have either have had its solar panels slightly closer to parallel or may have had panels that incompletely deployed, as suggested by its high count of values around 600 cm 2 .
Subsequently derived F10.7 adjustments corresponding to each percentile of the optimized areas all exhibit a drop around 26 May before they peaked during the main phase of the storm (Figures 10 and 11), similar to the behavior of the dSMA. They all consistently declined during the recovery phase. This prestorm drop may have been due to FOPT responding to the peak in F10.7 on the 27th, which occurred just over a day before the peak intensity of the storm on the 28th. This is clearest in Figure 12, where the peak negative F10.7 adjustments preceded the peak uncorrected F10.7 values by 1 day. There was a general trend of the F10.7 adjustments becoming less positive as a function of increasing percentile, while the overall behavior was preserved. The closeness of the lines corresponding to the 50th and 75th percentile suggests that this behavior tapered as the percentile increases. This is shown by the markedly smaller mean difference between the 75th percentile areas and 50th percentile areas (∼106.03 cm 2 ), compared with the 50th percentile areas and the 25th percentile areas (∼216.18 cm 2 ). As the cross-sectional area increased, the drag became an increasingly larger force, resulting in MOA lessening the contribution of increased F10.7.
Similar to the F10.7 adjustments, MOA's a p adjustments exhibit an increase beginning on or before 28 May, but this increase was much sharper, as the adjustments jump from −5 nT on the 28th to +80 nT on the 29th. This kept the 3-hr a p to MSISE-00 much higher after they reached their peak value of ∼140 nT on the 28th, before dropping sharply immediately after the start of the 29th. We will adjust this sharp drop in the future, allowing for a blending of the corrected 3-hr a p back down to the OMNIWeb values.
MOA's F10.7 and a p adjustments are used to drive the perturbations needed to get MSISE-00 to provide the best density predictions. The linearly interpolated median adjustments are applied to F10.7 inputs to MSISE-00 during the initial and recovery phases of the storm and to the a p inputs to MSISE-00 during the main phase of the storm, resulting in corrected densities along that satellite's orbit. The strength of this technique lies in that it derives adjustments to global drivers for MSISE-00, allowing a specification of the density at any other location in the thermosphere. This can be used not only to derive better densities for orbit propagation for other satellites but also to specify model biases for prediction.

Swarm Density Comparisons
In order to validate the technique, the MSISE-00 mass densities along the Swarm satellite orbit tracks were calculated using both the unperturbed drivers and the MOA-derived perturbed drivers. We first consider the effects of the selection of the quartile of the optimized area of the Flock CubeSats on the resulting orbit-averaged corrected densities along the trajectories of each Swarm spacecraft. This minimizes bias in the predicted drivers, since the absolute time-averaged projected area of each Flock satellite is not Figure 13. Orbit-averaged densities along-track Swarm-A (a) and Swarm-B (b) and for OMNIWeb inputs, with results shown corresponding to the percentiles of the optimized area. The gray area around the Swarm orbit-averaged densities denotes the ±1σ boundaries. publicly known. With enough data-model comparisons, an appropriate quartile can be selected and used for all future simulations. Figure 13 shows a comparison between Swarm-derived and MOA-derived densities with different quartiles of the area selected. The average altitudes of Swarm-A and Swarm-B during the time chosen were ∼452 and ∼515 km, respectively. From those, we derived a time resolution for the orbit-averaged densities of ∼93.6 and ∼94.9 min for each satellite, respectively.
During the main phase, the peak orbit-averaged densities returned by each MOA's percentiles were very close to the peak orbit-average densities for Swarm, with the 50th and 75th percentile cases being closest together. Along Swarm-A, percent differences between peak orbit-averaged MOA densities and those from Swarm were ∼4.5%, ∼5.8%, and ∼7.2% for the 25th, 50th, and 75th-percentile optimized areas, respectively ( Figure 13a). Note that increasing the percentile of the optimized area slightly reduced the accuracy in the peak orbit-averaged density along Swarm-A.
Along Swarm-B, this trend was reversed, with the 25th-percentile case yielding a peak orbit-averaged density error of ∼21.4% compared to ∼10.4% and ∼9.1% for the 50th-and 75th-percentile cases, respectively (Figure 13b). The rest of this study only considers the 75thpercentile case.
MOA alters the peak value a bit, but its main contribution is that it extends the period of high densities by a full day. While the peak is still too narrow, MOA does a much better job at capturing the second day response to the storm, due to the dramatic increase in a p that is added. Figure 14 overlays orbit-averaged densities from MSISE-00, MOA, and GPS-derived density data along Swarm-A and Swarm-B. It shows that the error in orbit-averaged density is reduced for MOA compared to uncorrected MSISE-00 along Swarm-B, but not Swarm-A, though the values of ρ T for MOA showed improvement over MSISE-00 along both orbits. Along Swarm-A, the integrated effect of the increased width of the storm can be quantified by the percent difference between MOA ρ T and Swarm ρ T , which was ∼4.0%, compared to a percent difference in integrated densities of ∼20.8% for MSISE-00. Along Swarm-B, the percent difference between MOA and Swarm was slightly lower, at only ∼3.7%, compared to ∼34.6% for MSISE-00. Visually, this translated to a widening of the peak in the orbit-averaged density. This behavior can be observed in Figure 14, where MOA orbit-averaged densities between 05-28 and 05-29 were noticeably higher than their uncorrected MSISE-00 counterparts. MOA matches Swarm quite well during the recovery of the storm, showing the same drop-off in density around the 29th. Along Swarm-A and Swarm-B, MOA attempted to recreate the second peak in the density occurring just before the 29th but was unable to reach the necessary amplitude to do so to the most accurately degree possible. This is especially obvious along Swarm-B. Finally, along both Swarm spacecraft, MOA placed the time of the peak orbit-averaged density as identical to that of MSISE-00. The times of the peak density, as shown in Tables 2 and  3, were an average of 4 hr prior to that of the Swarm data. This may be due to the fact that MSISE-00 and MOA are unable to account for the time delay between when geomagnetic indices peak and when local density at the spacecraft peaks. MSISE-00 and MOA apply Figure 14. Orbit-averaged densities along-track Swarm-A (a) and Swarm-B (b) for OMNIWeb inputs. MOA results corresponding to the 75th percentiles are in cyan, uncorrected MSISE-00 results are show in red, and Swarm data are shown in black. The gray vertical dashed lines represent the boundaries of the integral used to calculate ρ T across each data set, and the gray area around the Swarm orbit-averaged densities denotes the ±1σ boundaries.  (Bruinsma et al., 2006).
As shown in Table 2, usage of MOA's adjustments along Swarm-A resulted in a max ρ, δ P , η, and ρ T all closer to those of the Swarm data, compared to uncorrected MSISE-00. The closeness of ρ T in particular demonstrates that MOA greatly improved the width of the peak density during the main phase of the storm. Table 3 confirms these improvements along Swarm-B to a slightly greater degree, with the exception of max ρ, which MOA overestimated by 10.1%, resulting in a value of η different from that of Swarm by 11% compared to 4% for MSISE-00. The percent difference in ρ T along Swarm-A between MSISE-00 and Swarm was 4%, but this dropped slightly to 3.7% along Swarm-B. MOA's adjustments also predicted the time of the maximum orbit-averaged density more accurately along Swarm-A (early by 3.1 hr compared to Swarm), compared to Swarm-B (early by 4.8 hr compared to Swarm). MOA's adjustments thus performed better widening the peak density along Swarm-B but modeled the density time response better along Swarm-A.
Outside of the main phase, where only the adjustments to F10.7 were applied, MOA performed marginally better than MSISE-00 just before initial storm onset along Swarm-A, as well as just after the recovery phase along Swarm-B due to lowering the F10.7 just before the storm, and slightly raising it just after the storm. Both MSISE-00 and MOA overestimated the density during the initial and recovery phases along Swarm-A, though in the former case, MOA did so to a lesser degree, especially along Swarm-A. Along Swarm-B, the overestimation before the initial phase and after the recovery phase was much less noticeable, with MOA and MSISE-00 densities all residing within the 1σ boundaries of the Swarm densities. The adjustments to F10.7 were rather marginal compared to 3-hr a p and never exceeded ∼|17| sfu, compared to 3-hr a p adjustments, which grew to a maximum of +80 nT on the 29th.

Discussion and Conclusion
While widespread use of empirical density models like those of the MSISE-00 family are a testament to their efficacy and utility, this study highlights some of their limitations during high levels of geomagnetic activity. These limitations are in part largely owed to the models not being designed to fully capture the highly dynamic variations during these relatively rare events. Algorithms like HASDM (Storz et al., 2005) circumvent this problem by estimating corrections to the thermospheric temperature profile in order generate a corrected global density field. As biases are to be found in any model, the usage of correction factors to account for them is expected, and it is only a matter of how they are to be derived. Doornbos et al. (2008) used an algorithm in Picone et al. (2005) to derive thermospheric densities along satellite orbits using TLEs and then used a least-squares adjustment to estimate a set of calibration parameters to height-dependent scale factors of the densities and CIRA-72 temperatures in order to correct thermospheric densities. Pilinski et al. (2016) presented the Dragster Model, an assimilative algorithm that uses an ensemble of full-physics atmospheric models, a suite of high-latitude forcings and solar forcings, and an ensemble Kalman filter to specify real-time and 3-day forecasted densities, compositions, and winds along satellite orbits to compute drag estimates. Gondelach and Linares (2020) also presented a powerful method for both modeling storm time densities and performing density prediction by assimilating TLEs and historical empirical model density estimations into a ROM that allows for very accurate real-time density prediction. While MOA is similar to Doornbos et al. (2008) in that the TLE-derived SMA is the driving force behind the generation of the density adjustments, MOA differs in that changes in SMA are not used to derive density directly between subsequent TLEs but rather model input adjustments, which are in the form of median adjustments to F10.7 and a p , rather than height-dependent scale factors to the density and corrections to the thermospheric temperature profile. MOA is distinguished from Dragster in that it requires far fewer forcings to function, does not require an ensemble Kalman filter, does not rely on a super ensemble of physics-based atmospheric models, and does not predict thermospheric composition and winds. MOA differs from the ROM in that it is simpler to implement, not requiring the estimation of ballistic coefficients with the use of a Kalman filter, the modification of equinoctial elements, or the assimilation of accelerometer and GPS data from multiple sources. Additionally, MOA demonstrates the capacity to yield improved storm-time density predictions with the use of TLEs from only a few satellites, all which were in the same orbital plane. This differs also from the HASDM, which uses high temporal resolution drag information from the trajectories of ∼75 different calibration satellites at a variety of altitudes (Storz et al., 2005). From the close correspondence between orbit-averaged densities returned by MOA and GPS-derived data collected by the Swarm spacecraft, this preliminary study suggests that TLEs can be used effectively to correct empirical model densities through the process of orbit-error minimization in a relatively simple process.
The most prominent limitation of this technique is that its sole reliance on TLE-driven orbit propagation places a limit on the power of the obtained adjustments. Since TLEs are typically available once every day or two, the global adjustments to F10.7 and 3-hr a p obtained by MOA run the risk of "smoothing over" rapid changes the density. This is most clearly observed in two ways regarding our results: the width of the peak corrected orbit-averaged densities and the lack of distinct density features in the corrected orbit-averaged densities that are present in the orbit-averaged densities in Swarm's GPS-derived data, such as the second "peak" in the density during the main phase. While individual TLEs do not provide good temporal resolution, they are nevertheless available for a plethora of orbiting objects at a wide variety of altitudes. We aim to improve MOA by determining how storm-time density modeling can be improved with this simple method by using TLEs from many more objects. We also note that at present, MOA associates the corresponding adjustments to F10.7 and 3-hr a p found for each satellite with the either noon or the beginning of the day, respectively, and then estimates adjustments between those times with linear interpolation. We aim to see if improvements in density prediction may result from associating adjustments with the TLE epoch directly, taking the adjustments for each satellite, and then computing new adjustments by forming a univariate spline through the median adjustments found across all satellites combined. This may properly "fill in the gaps" as TLEs for each satellite are not reported at the same epoch each day, allowing for any resulting index adjustments to generate the distinct features in orbit-averaged density derived from GPS data.
Future work will involve a systematic study of MOA's capabilities across a series of cataloged geomagnetic storms that have occurred within the last decade. It will be necessary to determine the efficacy of this algorithm in handling anomalously large geomagnetic disturbances such as those characteristic of the St. Patrick's Day Storms of 2013 and 2015 and also minor disturbances such as those characteristic of 6 March 2016. It is additionally worth noting that MOA may perform differently when calculating adjustments to MSISE-00 or other models such as the those of the Jacchia family and the drag-temperature model. Achieving these milestones will elucidate the utility of using simple methods to address the common challenge of storm-time thermospheric density modeling and prediction.

Data Availability Statement
Data used to generate the figures in Sections 2 and 3, along with instructions on how to generate those figures are located in the University of Michigan's Deep Blue repository: https://deepblue.lib.umich.edu/data/ concern/data_sets/6q182k29p.