Estimating Satellite Orbital Drag During Historical Magnetic Superstorms

Understanding extreme space weather events is of paramount importance in efforts to protect technological systems in space and on the ground. Particularly in the thermosphere, the subsequent extreme magnetic storms can pose serious threats to low-Earth orbit (LEO) spacecraft by intensifying errors in orbit predictions. Extreme magnetic storms (minimum Dst $\leq$ --250 nT) are extremely rare: only 7 events occurred during the era of spacecraft with high-level accelerometers such as CHAMP (CHAllenge Mini-satellite Payload) and GRACE (Gravity Recovery And Climate experiment), and none with minimum Dst $\leq$ --500 nT, here termed magnetic superstorms. Therefore, current knowledge of thermospheric mass density response to superstorms is very limited. Thus, in order to advance this knowledge, four known magnetic superstorms in history, i.e., events occurring before CHAMP's and GRACE's commission times, with complete datasets, are used to empirically estimate density enhancements and subsequent orbital drag. The November 2003 magnetic storm (minimum Dst = --422 nT), the most extreme event observed by both satellites, is used as the benchmark event. Results show that, as expected, orbital degradation is more severe for the most intense storms. Additionally, results clearly point out that the time duration of the storm is strongly associated with storm-time orbital drag effects, being as important as or even more important than storm intensity itself. The most extreme storm-time decays during CHAMP/GRACE-like sample satellite orbits estimated for the March 1989 magnetic superstorm show that long-lasting superstorms can have highly detrimental consequences for the orbital dynamics of satellites in LEO.


Introduction
Magnetic storms are global phenomena that occur due to the interaction of solar perturbations with the Earth's magnetosphere (Gonzalez et al., 1994). The most intense and severe magnetic storms are commonly caused by coronal mass ejections (CMEs) (Balan et al., 2014;Daglis, Thorne, Baumjohann, & Orsini, 1999;Gonzalez et al., 1994). CMEs usually have a shock at their leading edge that is promptly followed by a sheath and a magnetic cloud (Balan et al., 2014;Gonzalez et al., 1994;Kilpua et al., 2019). Extreme magnetic storms are caused by the impact of extremely fast CMEs on the Earth's magnetosphere (Tsurutani & Lakhina, 2014), usually associated with highly depressed values of the southward component of the interplanetary magnetic field (Balan et al., 2014;Daglis et al., 1999;Gonzalez et al., 1994;Kilpua et al., 2019;Tsurutani & Lakhina, 2014).
Extreme space weather events like severe magnetic storms have been recognized by the U.S. Federal Government through the National Space Weather Strategy and Action Plan (National Science and Technology Council, 2015a, 2015b) as a natural hazard, and the need to establish benchmarks for extreme space weather events has also been recognized by the scientific community (e.g., Jonas, Murtagh, & Bonadonna, 2017;Lanzerotti, 2015;Riley et al., 2018). The intensity of magnetic storms is usually measured by depletions of the ground horizontal magnetic field component recorded by magnetometers located at mid-and lowlatitudes by means of the disturbance storm time (Dst) index (section 2.1). Extremely severe events, here termed magnetic superstorms, with minimum Dst ≤ -500 nT, are notably rare (Chapman, Horne, & Watkins, 2020;Cliver & Dietrich, 2013;Hayakawa, Ebihara, Willis, et al., 2019;Riley et al., 2018;Vennerstrøm et al., 2016). For instance, the March 1989 event, the only superstorm occurring during the space age (Meng, Tsurutani, & Mannucci, 2019), is well-known for the occurrence of low-latitude aurorae (Allen, Sauer, Frank, & Reiff, 1989;Pulkkinen, Bernabeu, Eichner, Beggan, & Thomson, 2012;Rich & Denig, 1992) and intense geomagnetically induced currents (GICs) which caused the blackout of the Hydro-Québec system in Canada for several hours, leading to serious economic losses (Bolduc, 2002;Kappenman, 2006;Pulkkinen et al., 2017). However, though arguably, the most extreme ground horizontal magnetic field perturbation (∼ -1600 nT) on record was recorded by the Colaba station during the Carrington event of September 1859 (Hayakawa, Ebihara, Willis, et al., 2019;Siscoe, Crooker, & Clauer, 2006;Tsurutani, Gonzalez, Lakhina, & Alex, 2003). Since that is the only known low-latitude data set available to date, a global analysis of that storm cannot be performed (Blake et al., 2020;Cliver & Dietrich, 2013;Hayakawa, Ebihara, Willis, et al., 2019;Siscoe et al., 2006). For this reason, the Carrington event is not addressed in this paper.
The most extreme magnetic storm experienced by CHAMP and GRACE took place in November 2003 with minimum Dst = -422 nT. Consequently, there are no assessments of satellite drag in LEO during magnetic superstorms inferred from high-accuracy accelerometer data. The orbital degradations of CHAMP and GRACE associated with the November 2003 event 60 hrs through stormy times were, respectively, ∼ -160 m and ∼ -71 m (Krauss, Temmer, Veronig, Baur, & Lammer, 2015;, much more severe than the natural drag caused by the quiet-time backgorund density estimated by Oliveira and Zesta (2019), namely -24.11 m and -6.86 m, respectively. Hence, these are the most extreme storm-time orbital decays measured with high-quality accelerometer data. In order to empirically estimate drag effects during magnetic superstorms, standard Dst data and ground magnetometer data of historical superstorms reconstructed from historical archives are used by a thermospheric empirical model (section 2.3) for density computations (section 2.4). These events occurred in March 1989 (Allen et al., 1989;Boteler, 2019), with the traditional Dst index available, September 1909Hayakawa, Ebihara, Cliver, et al. (2019Silverman (1995), May 1921(Hapgood, 2019Silverman & Cliver, 2001), and October/November 1903(Lockyer, 1903Ribeiro, Vaquero, Gallego, & Trigo, 2016), with an alternative version to the Dst index available. These four magnetic superstorms are here examined because they are the only events with known and complete magnetograms that satisfy the threshold Dst/Dstlike ≤ -500 nT. The main characteristics of these storms' effects will be presented in section 3.1. Effects of storm time duration associated with minimum values of Dst and Dst-like data will be estimated and compared. As a result, this effort will improve our understanding of severe satellite orbital drag effects in LEO caused by magnetic superstorms. 2 Data, model, and a framework for orbital drag estimations 2.1 Disturbance storm time indices In this study, magnetic activity is represented by the Dst index provided by the World Data Center for Geomagnetism, Kyoto, Nose, Iyemori, Sugiura, and Kamei (2015). This 1-hr-resolution index was defined in 1957, the International Geophysical Year (IGY), as described by Sugiura (1964). Specifically, Dst is computed by averaging latitudinally weighted horizontal magnetic field perturbations, with a background removal scheme, recorded by mid-and low-latitude stations with reasonably even longitudinal separation according to the expression where ∆H i is the horizontal magnetic perturbation of the i-th station, and Λ i is the contemporary magnetic latitude of the i-th station. The colored stars in Figure 1 show the stations, with their corresponding names, abbreviations, and geographic locations, used to compute standard Dst after the IGY.
where ∆H i is the horizontal magnetic perturbation of the i-th station, and Λ i is the contemporary magnetic latitude of the i-th station. The colored stars in Figure 1 show the stations, with their corresponding names, abbreviations, and geographic locations, used to compute standard Dst after the IGY. Additionally, recent efforts have been undertaken to provide alternative (but similar) versions to the standard Dst index for historical magnetic superstorms with archival material. The events took place in October/November 1903(Hayakawa, Ribeiro, et al., 2020), September 1909(Love et al., 2019b(Love et al., 2019a. This alternative index, also with resolution of 1 hr, was reconstructed with data obtained from four low/midlatitude stations, with the best possible longitudinal separation, and is represented here by Dst † (Dst "dagger"). The corresponding contemporary magnetic latitudes were computed by the authors. A background removal scheme similar to the one used to calculate Dst is used in the source papers as well. The stations used to compute Dst † used in this study are shown by the colored crosses in Figure 1. Therefore, the Dst † index is given by The Dst † data for the magnetic superstorms used here are available as supporting information provided by the respective references (Hayakawa, Ribeiro, et al., 2020

Neutral mass density data
CHAMP and GRACE neutral mass density (ρ) data obtained from their respective highaccuracy accelerometers are used in this work. CHAMP was launched in 2001 at the initial altitude 456 km and orbital inclination 87.25 • . It covered each 1 hr local time in 5.5 days with orbital period 90 min. The GRACE-A and -B spacecraft were launched in 2002 at the initial altitude 500 km and orbital inclination 89.5 • . The GRACE constellation covered each 1 hr local time in 6.7 days with orbital period 95 min. GRACE-A flew ∼220 km ahead of GRACE-B. As discussed in Oliveira and Zesta (2019), only GRACE-A data are used, henceforth GRACE data, because GRACE-A data show higher quality than GRACE-B data. CHAMP re-entered in 2010, while GRACE re-entered in 2018. Uncertainties and calibration techniques of both missions have been discussed by many papers (e.g., Bruinsma, Tamagnan, & Biancale, 2004;Doornbos & Klinkrad, 2006;Flury, Bettadpur, & Tapley, 2008).
The density data used in this study are normalized and intercalibrated as described in Oliveira, Zesta, Schuck, and Sutton (2017) and Zesta and Oliveira (2019). Basically, the Jacchia-Bowman 2008(hereafter JB2008, Bowman et al., 2008, see below) empirical model computes quiet-time densities (ρ 0 ) in order to obtain the background state for the quiet thermosphere. This approach ensures that the ratio and the difference between the stormtime and quiet-time densities are as close to one (ρ/ρ 0 ≈ 1) and zero (ρ -ρ 0 ≈ 0) as possible, respectively. As a result, storm-time density enhancements can be extracted more effectively ).

The Jacchia-Bowman 2008 (JB2008) empirical model
The first clear link between magnetic activity and satellite orbital drag effects was established by Jacchia (1959), who used Sputnik 1958δ1 data to discover that its altitude significantly decayed during an extreme magnetic storm. He correctly realized that this effect occurred due to augmented density levels at the satellite's altitude. Later on, this discovery led scientists to develop thermospheric empirical models such as the Jacchia 70 model (J70; Jacchia, 1970), the Mass Spectrometer Incoherent Scatter model series (MSIS; Hedin, 1987) which were the precursors to the Naval Research Laboratory Mass Spectrometer Incoherent Scatter Extended (NRLMSISE-00; Picone, Hedin, Drob, & Aikin, 2002), the Drag Temperature Model (DTM2013; Bruinsma, 2015), the High Accuracy Satellite Drag Model (HASDM, Storz, Bowman, Branson, Casalic, & Tobiska, 2005), and more importantly for this work, the JB2008 model. A description of the JB2008 model along with other popular thermospheric empirical models has recently been provided by He et al. (2018). The JB2008 empirical model computes thermospheric neutral mass density from a single parameter, the exospheric temperature (see equation 2 in . This temperature depends on several satellite parameters such as latitude, local time, and altitude. Additionally, JB2008 uses the solar radio flux at wavelength 10.7 cm, indicated by the F10.7 index, to account for thermospheric heating due to solar UV radiation (Bowman et al., 2008). Finally, a term that depends on Dst in the exospheric temperature represents the magnetic activity contribution, but JB2008 uses the 3-hr time resolution ap index for intervals when Dst > -75 nT. Dst and Dst † data of the historical magnetic superstorms recorded by the stations shown in Figure 1 will be used along with LEO satellite orbital data during the event of November 2003 to estimate the subsequent drag effects.

Orbital drag computations
Neutral mass densities are derived by high-accuracy accelerometers according to the drag equation (Prieto et al., 2014): where a d is the spacecraft acceleration caused by drag forces; ρ is the local thermospheric neutral mass density; C D is the drag coefficient; S/m is the area-to-mass ratio; and V is the relative velocity between the spacecraft velocity ( V s/c ) and the ambient neutral wind velocity ( V wind ). In this equation, all quantities are presumably known, and therefore it is solved for ρ in order to yield density. However, these parameters (particularly C D ) can introduce significant errors in density computations (Moe & Moe, 2005;Prieto et al., 2014;Zesta & Huang, 2016). In this study, drag coefficients computed with error mitigation methods by Sutton (2009) were used. Chen, Xu, Wang, Lei, and Burns (2012) provide the following expression for the computation of storm-time orbital decay rate: with a being the semi-major axis of the satellite orbit here replaced by the temporal Earth's radius plus satellite altitude , G = 6.67×10 −11 m 3 ·kg −1 ·s −2 the gravitational constant, M = 5.972×10 24 kg the Earth's mass, and ∆ρ the difference between the modeled storm-time and quiet-time densities. As outlined by Oliveira and Zesta (2019), the daily average of the semi-major axis a is represented by a . A comparison between the use of both a computation methods for a magnetically quiet day (not shown) reviewed a very minimal difference in da/dt. In addition, Krauss et al. (2015) and Oliveira and Zesta (2019) found the same results for the orbital decay of GRACE during the November 2003 storm.
Finally, the storm-time orbital decay (d(t)) is computed by the sum over all da/dt values along the satellite's path for any (t 1 , t 2 ) interval: where a (t) = da(t)/dt.

The selected magnetic superstorms
The benchmark event for the current study occurred in November 2003. That storm had minimum Dst = -422 nT, the most intense magnetic storm event with both CHAMP and GRACE neutral mass density data available. Ground magnetometer data and neutral mass density data for the GRACE satellite are shown in Figure 1 of Zesta and Oliveira (2019). The solar flux F10.7 index increased from 151 sfu (solar flux units) on 19 November to 175 sfu on 23 November. The Dst and F10.7 indices for that storm are shown in Figure 2. Figure 3 documents the orbits of CHAMP and GRACE in the time interval from 19 to 23 November 2003. The dial plots show orbits as a function of magnetic latitudes (MLATs) and magnetic local times (MLTs). The magnetic coordinate system used is the Altitude-Adjusted Corrected Geomagnetic Model (Laundal & Richmond, 2017;Shepherd, 2014, AACGM,). The left column shows altitudes for CHAMP, while the right column shows altitudes for GRACE. The top row indicates data for the northern hemisphere, while the bottom row indicates data for the southern hemisphere. The colorbars indicate altitudes for both satellites in the same periods.
CHAMP is in a near noon-midnight orbit. The orbit altitudes of CHAMP increased at high latitudes and at the magnetic poles of both hemispheres and decreased at midand low-latitudes. Similar behavior is shown by GRACE whose orbits were confined within the mid-noon/dusk and mid-midnight/dawn sectors. Therefore, both spacecraft provide reasonable coverage between the day and night sectors. The altitude variations shown in  CME leading edges are usually associated with the occurrence of positive jumps in the Dst index, while its sudden depression is associated with the arrival of CME magnetic material or sheaths (e.g., Gonzalez et al., 1994;Kilpua et al., 2019). The first perturbation, termed storm sudden commencement (SSC), is caused by the shock compression (e.g., Oliveira et al., 2018;Shi et al., 2019), while the second event, termed storm main phase, is associated with strong driving of the magnetosphere via magnetic reconnection (e.g., Daglis et al., 1999;Gonzalez et al., 1994;Kilpua et al., 2019). Examples of SSCs and storm main phases represented by the Dst and Dst † indices during magnetic superstorms caused by fast CMEs are illustrated in Figure 4.  (Hapgood, 2019). Given the similarities of UTs and GMTs, here they will be used interchangeably (Hapgood, 2019). The highlighted areas of each panel correspond to the time interval between SSC and minimum Dst/Dst † occurrences, which also mark the beginning of the storm recovery phase. This time interval will henceforth be referred to as the storm development duration time in this paper. Panels (a) and (b) show that the 1903 event is the weakest (minimum Dst † = -531 nT), whilst the 1921 event is the strongest (minimum Dst † = -907 nT) amongst all events. In contrast, the development duration times of both events are almost the same, ∼ 14 hr and ∼ 12 hr, respectively. Storm strengths can be estimated by computing how fast Dst (or Dst † ) is depressed during storm development.  peak at SSC compression by the development time. This provides a quantifiable measure of the impactfulness of the storm, meaning that storms with very low amplitude rates are commonly associated with high geomagnetic activity (e.g., Gonzalez et al., 1994). The estimated amplitude rates are -46.4 nT/hr and -80.0 nT/hr for the October/November 1903 and May 1921 events, respectively. These numbers explain why the effects of the 1921 event, such as equatorial extent of low-latitude aurorae (Chree, 1921;Silverman & Cliver, 2001), and GIC impacts on contemporary telegraph systems (Hapgood, 2019;Kappenman, 2006) were more severe than the effects of the 1903 event, mostly represented by mid-latitude aurorae (Hayakawa, Ribeiro, et al., 2020;Page, 1903), and local GIC impacts on contemporary telegraph systems in the United States and in the Iberian Peninsula (Hayakawa, Ribeiro, et al., 2020;Ribeiro et al., 2016 (Pulkkinen et al., 2012;Rich & Denig, 1992). Intense GICs occurred during both events, with several reports of geophysical disturbances on telegraph systems in 1909 (Hapgood, 2019;Hayakawa, Ebihara, Cliver, et al., 2019;Love et al., 2019b;Silverman, 1995), and on power transmission lines in 1989, particularly the power blackout in Québec, Canada (Allen et al., 1989;Boteler, 2019;Kappenman, 2006;Oliveira & Ngwira, 2017). During the 1989 event, the only event with satellite-based data amongst the four superstorms, the number of space objects "lost" in LEO increased dramatically around periods of maximum intensity due to errors introduced by storm heating effects into tracking systems (Allen et al., 1989;Burke, 2018;Joselyn, 1990). The left (not highlighted) part of Table 1 summarizes these storm properties. A comprehensive comparison of GIC effects caused by the superstorms on the contemporary ground infrastructure, i.e., telegraph systems and power grids, is a difficult task to be accomplished. However, the comparisons above show that the latitudinal extent of the auroral oval was more equatorward for the events with lower amplitude rates (May 1921 and September 1909 events). Next, the effects of these amplitude rates on storm-time orbital drag will be evaluated and compared for the four historical magnetic superstorms studied in this paper.

Storm-time orbital drag effects 3.2.1 The November 2003 extreme magnetic storm
Since the November 2003 magnetic storm is chosen in this work as the benchmark event, CHAMP and GRACE thermospheric neutral mass density response and the subsequent orbital drag effects for that storm are shown here, and an effort to compute errors associate with drag effects is performed. The orbital drag framework of Oliveira and Zesta (2019) summarized in section 2.4 is used for the drag computations. The Dst and F10.7 indices for the benchmark event are shown in Figure 2. Figure 5 documents density observed by CHAMP (upper left) and GRACE (upper right) along with JB2008 quiet-and storm-time density predictions for the benchmark storm. The dynamics of that storm orbital effects were discussed in detail by Oliveira and Zesta (2019) particularly for GRACE's case. In general, the predicted density dynamics follows CHAMP and GRACE observations quite well, but there are remarkable differences with respect to density values. Firstly, density for both satellites was highly underestimated during heating and cooling of the thermosphere, being more severe in CHAMP's case. Secondly, overestimations of JB2008 results for GRACE's orbit are higher than CHAMP's during thermospheric recovery times . This density dynamics is reflected on the observed/predicted density relative errors shown by the solid purple lines of Figure 5 in the lower left panel for CHAMP and in the lower right panel for GRACE. Figure 6 shows drag results for CHAMP's and GRACE's orbits, respectively. The grey highlighted areas in all panels indicate the storm development time (13 hrs) similarly to the ones shown in Figure 4. The odd rows of these figures show storm-time orbital decay rates (da/dt, equation 4), while the even rows show storm-time orbital decay (d, equation 5). The left column shows observation results, while the right column shows JB2008 results. In the even rows, the magenta lines indicate the "natural" orbital decay caused by the background density if there was no storm activity. The background density for storm-time drag computations was obtained by the method developed by .
As a result of the density dynamics shown in Figure 5, at t = 72 hrs, the storm-time orbital decays estimated for CHAMP and GRACE shown in Figure 6 are underestimated by 13.57% and overestimated by 16.32%, respectively. However, the uncertainties associated with the magnetic superstorms here investigated should differ from these uncertainties for two reasons: (i) the superstorms are more intense, and (ii) the superstorms had different development times and therefore different magnetic activity during different times. These uncertainties are obtained for the most extreme magnetic storm during both CHAMP and GRACE commission times, and therefore may represent an upper limit of JB2008 uncertainties for extreme magnetic storms with high-level thermosphere neutral mass density available.  Table  1). The same is shown in panels c1 (CHAMP) and c2 (GRACE) for the superstorms of March 1989 (red line) and September 1909 (blue line). In this case, the storms had very similar intensities, but different development durations (Table 1). The storm-time orbital degradation (equation 6), is shown for CHAMP (panels b1 and d1) and GRACE (panels b2 and d2). The same colors used to represent da/dt results in panels a1/c1 and a2/c2 above are used to represent d results in panels b1/d1 and b2/d2. Figure 7a1 shows that da/dt values during October/November 1903 for CHAMP were very close to zero before CME impact. On the other hand, da/dt values preceding the stormy period of May 1921 show some variations (meaning ∆ρ is not necessarily close to zero),  presumably linked to high magnetic activity shown by ground magnetometer data during the same pre-storm period (Hapgood, 2019;Love et al., 2019b). CHAMP da/dt values for the 1921 event decreased faster in comparison to minimum da/dt values for the 1903 event. Similar orbital drag dynamics is observed for GRACE (a2), but the absolute values of the drag response are smaller (Table 1) because GRACE operated at higher altitudes in comparison to CHAMP (Krauss, Temmer, & Vennerstrom, 2018;. The da/dt results for CHAMP and GRACE are summarized in Table 1.

The magnetic superstorms
For the same pair of storms, the storm-time orbital degradations of CHAMP (panel b1) at the end of 72 hr after CME impact were -91.23 m and -196.24 m for both events, respectively. The same estimated results for GRACE (b2) are -60.40 m (1903) and -142.09 m (1921). Comparatively, the percentual difference between drag effects during both superstorms for CHAMP (115.10%) are higher than the percentual difference of the superstorm intensites (70.81%) most likely because the magnetosphere was hit by another CME on 16 May 1921 (Love et al., 2019a, Figure 4;), leading to an additional magnetosphere energization during its recovery, which in turn impacted drag effects. Similarly, the orbital drag relative difference is higher in the case of GRACE (135.25%), when compared with the case of CHAMP. As suggested by   Figure 10), this is presumably due to the interplay between heating propagation from auroral-to-equatorial latitudes and (possibly) the direct uplift of neutrals at low and equatorial latitudes more evident at altitudes higher than 400 km Tsurutani et al. (2007).
In summary, the main features that arise from the comparison between these events are: (i) CHAMP and GRACE decayed faster during the most intense event (1921) due to its sharper negative excursion of the Dst † index and lower amplitude rate (Figure 4a and b; Table 1); and (ii) the relative differences between d for both events do not closely follow the relative differences between minimum Dst † values. This is likely the case because the magnetosphere was struck by another CME during its recovery, increasing the magnetospheric activity which in turn affected the subsequent orbital drag effects. Tables 1 and 2  May 1921 is 53.98% stronger a Percentual differences between more severe (March 1989) with respect to less severe (May 1921) drag effects summarize these results.
The comparisons between estimated drag effects for the March 1989 and September 1909 superstorms are remarkably different. These events had very similar strengths (similar minimum Dst and Dst † values), but their development times were quite distinct. Figure 7c1 shows that 1909 CHAMP da/dt values could have shown a very sharp negative excursion after CME impact, which follows very closely the same feature in the Dst † index ( Figure  4d). The minimum da/dt value (-285.14 m/day) for the September 1909 superstorm was reached shortly before minimum Dst † . On the other hand, the March 1989 drag effects are quite different, since da/dt decreased more slowly in comparison to the former case due to the differences in storm development amplitude rates. This is explained by the fact that the magnetosphere was most likely struck by multiple CMEs while the storm main phase was developing (Boteler, 2019;Fujii et al., 1992;Lakhina & Tsurutani, 2016). Similarly to the 1909 case, the minimum da/dt value (-621.29 m/day) occurred shortly before minimum Dst occurrence. The thermosphere recovery of the 1989 superstorm took longer than the thermosphere recovery of the 1909 superstorm, most likely because the magnetosphere was hit yet by more CMEs shortly after the beginning of the magnetosphere recovery ( Figure  4c). A similar behavior is shown by the GRACE results, panel c2, but with smaller absolute values due to higher GRACE altitudes. The relative differences between da/dt peak values of CHAMP and GRACE for both superstorms are 117.30% and 145.73%, even though both events had approximately the same minimum Dst and Dst † values and very different storm development durations and amplitude rates. Now the storm-time orbital degradations in both cases are evaluated. Figure 7d1 shows that CHAMP d decreased faster during the main phase of the 1909 event, reaching values near its minimum value around the beginning of storm recovery. This is a typical feature of drag effects triggered by a storm caused by an isolated CME (Krauss et al., 2018(Krauss et al., , 2015. Conversely, CHAMP's orbital degradation decreased more dramatically during the recovery of the 1989 superstorm. These drag effects correlate well with a very sharp negative excursion presented by the Dst index, which is also directly related with the occurrence of low-latitude aurorae and very intense GICs around the world (Allen et al., 1989;Hayakawa, Ebihara, Cliver, et al., 2019;Kappenman, 2006). This time also coincides with the loss of orbital control of several objects in LEO as shown by satellite-based data (Allen et al., 1989;Burke, 2018;Joselyn, 1990). The storm-time orbital decays for the 1909 and 1989 events are -96.61 m and -388.59 m for CHAMP and -62.14 m and -305.58 m for GRACE. Their relative differences are 302.22% and 391.76%, closely following the proportion of storm time developments in the case of CHAMP. Taking into consideration that both superstorms were almost equally intense, these results show that the storm time duration can play a major role in driving orbital drag effects. Note also that relative differences are higher in the case of GRACE, most likely explained by the reasons suggested by Oliveira and Zesta (2019) as mentioned before. Another striking difference concerning minimum Dst and Dst † values, storm development duration and subsequent amplitude rate impacts arises from the comparison between the May 1921 and March 1989 superstorms. The 1921 event was more than 50% stronger than the 1989 event, but active times during the latter lasted twice longer. The storm-time orbital decay for the March 1989 event was nearly twice more severe than the May 1921 event in both CHAMP's and GRACE's cases (Figure 7 and Tables 1 and 2). These results clearly reveal that a long-lasting magnetic superstorm can drive much more severe drag effects in comparison to a short-lasting, even stronger, superstorm. Tables 1 and 2 summarize the main results discussed in sections 3.1 and 3.2.
The results presented so far correspond to the storm-time orbital decay values estimated by JB2008. Furthermore, the uncertainties computed for CHAMP's and GRACE's orbital drag effects during November 2003 (section 3.2.1) can be used to obtain more realistic drag results. Results are shown in Table 3, where white cells show model results, whereas grey cells show corrected results. In these new computations, only assumptions on overall error levels (at t = 72 hrs) were used since realistic errors cannot be obtained for the different superstorms because there are neither CHAMP nor GRACE density data available during these superstorm times.
There are no solar wind nor interplanetary magnetic field data available for the magnetic superstorms discussed in this paper. Furthermore, it is important to emphasize that our statements concerning CME impacts are supported by our current knowledge of the underlying science: intense magnetic storms, particularly extreme events, are usually caused by CMEs (Balan et al., 2014;Daglis et al., 1999;Gonzalez et al., 1994;Kilpua et al., 2019;Lakhina & Tsurutani, 2016;Tsurutani & Lakhina, 2014).

Discussion and conclusion
Extreme magnetic storms (minimum Dst ≤ -250 nT) are very rare. Only 39 extreme events have taken place since the beginning of the space era (Meng et al., 2019), while only 7 extreme events were observed by CHAMP and GRACE . Additionally, only one magnetic superstorm (minimum Dst ≤ -500 nT) occurred since 1957, while none were ever observed by either CHAMP or GRACE. Therefore, current knowledge of thermospheric mass density response to magnetic supersotorms and the subsequent storm-time drag effects are very limited. Then, in order to estimate these effects, 4 historical magnetic superstorms with complete magnetograms were selected: one with standard Dst data (March 1989), and 3 with Dst † (Dst-like) data occurring on October/November 1903(Hayakawa, Ribeiro, et al., 2020), September 1909(Love et al., 2019b, and May 1921 (Love et al., 2019a). These Dst and Dst † data were used as input data for the JB2008 thermosphric empirical model for density computations. The extreme magnetic storm of November 2003 (minimum Dst = -422 nT), the most extreme event during CHAMP's and GRACE's commission times, at the altitudes ∼400 km and ∼490 km, respectively, was used as the benchmark event. The orbital drag framework provided by Oliveira and Zesta (2019) was used for drag estimations.
First, two events with different intensities but with approximately the same storm development times were compared (October/November 1903 andMay 1921). Although the 1921 superstorm was ∼ 70% stronger than the 1903 superstorm, the drag effects in the former were up to 135% more severe than the effects in the latter (GRACE's case). This is attributed to the likely impact of another CME during the recovery phase of the 1921 superstorm. Second, the other pair of superstorms, with very similar strengths, but with the September 1909 storm development being 3 times shorter than the March 1989 storm development, were compared. Results show that the relative difference of the storm-time orbital degradation for the 1989 event was about 400% higher than the 1909 event (GRACE's case). This is explained by the likely impacts of several CMEs on the magnetosphere during the main and recovery phases of the March 1989 superstorm (Boteler, 2019;Fujii et al., 1992;Lakhina & Tsurutani, 2016). Therefore, as opposed to latitudinal extent of aurorae, a superstorm with a smaller amplitude rate (absolute value) can cause more detrimental effects on orbital drag in comparison to an even stronger superstorm that develops faster (larger absolute value of amplitude rate). The CHAMP and GRACE storm-time orbital decays as predicted by JB2008 and corrected by errors obtained during the November 2003 event (Table 3) are much more severe than the orbital degradation due to the background densities during the benchmark storm shown in Figure 6 for CHAMP (-28.68 m) and GRACE (-9.59 m). For example, results for the March 1989 event show that the CHAMP storm-time orbital decay was estimated to be ∼-441.32 m: such value has never been measured by a LEO spacecraft with high-level accelerometers. Therefore, these results set a new basis for these effects. Despite the fact that these effects can have significant error levels particularly during the storm recovery phases due to the lack of nitric oxide cooling effects in the model (Bowman et al., 2008;Knipp et al., 2017;Mlynczak et al., 2003;, these results reveal the comparative roles of time durations and strengths of magnetic superstorms in controlling drag effects.
The results of this work clearly show that multiple CME impacts on the Earth's magnetosphere (as in the March 1989 superstorm), particularly occurring during active times, can largely enhance satellite orbital drag due to long and sustained storm times. These drag effects can be more severe when compared to drag effects during storms caused by a single CME leading to even more intense storms, but lasting shorter. Therefore, orbital drag forecasters should be aware of potential impacts of several CMEs on the terrestrial magnetosphere during ongoing magnetic storms (e.g., Zhao & Dryer, 2014, and many references therein). Additionally, different thermospheric empirical models should produce different results, with JB2008 outperforming NRLMSISE-00 and HASDM outperforming JB2008 Bowman et al. (2008), but with DTM2013 outperforming JB2008 (Bruinsma, 2015). In a future work, simulation results using different models of tens of historical severe and extreme magnetic storms, with minimum Dst ≤ -250 nT (Chapman et al., 2020;Meng et al., 2019;, will be statistically studied. 18