The Determination of Satellite Orbital Decay From POD Data During Geomagnetic Storms

In our previous study (Li et al., 2017, https://doi.org/10.1007/s11430-016-9052-1), we derived the satellite energy‐decay with a 20‐min resolution based on the Precise Orbit Determination (POD) data, using a proximate analytic approach to represent the time variation gravitational potential. In this follow‐up study, we improved the previous approach and calculated the POD‐based energy‐decays by using a numerical integration approach. Based on the precise energy‐decays, the orbital decays and decay rates with higher accuracy and resolution were further derived. The relative deviations of the orbital decays are generally less than 10% with respect to the accelerometer‐reference. The satellite orbital decays and decay rates derived from this approach were used to study the effects of geomagnetic activities and background density on the orbital changes. Our results show that, during the severe November 2003 storm, the storm‐induced orbital decay rate increased by a factor of 8 with respect to the quiet‐time reference. This POD‐based integration approach was also applied to study the orbital changes of multiple satellites at different altitudes during the September 2017 moderate storm. It is found that the storm‐induced orbital decay rates of Swarm‐B, Swarm‐A, and Gravity Recovery and Climate Experiment satellites increased by 100%–150% depending on their altitudes. Overall, the results suggest that our integration approach has better performance than the previous approach in deriving the orbital decay rate at solar minimum or at high altitude when the atmospheric density is relatively low.

. It was found that the satellite orbital decay during extreme Coronal Mass Ejections (CME)-driven storms changed dramatically (Knowles et al., 2001;Willis et al., 2005). Instead of CME-, the Corotating Interaction Regions (CIR)-driven storms could cause even larger orbital decay, since the CIR-driven storms last for a relatively longer period (Chen et al, 2012(Chen et al, , 2014Krauss et al., 2018). Based on the statistics study of the orbital decays of CHAMP and GRACE during hundreds of CMEs and CIRs, the orbital decay could be predicted in response to solar events (Krauss et al., 2020). Therefore, the orbital decays of more satellites from different altitudes during various solar conditions are in urgent need to increase the statistical samples and subsequently improve the accuracy of the orbit prediction.
The orbital decay of a satellite could be calculated approximately by integrating the drag force along the satellite trajectory (Anderson et al., 2009;Chen et al., 2012Chen et al., , 2014, which has been used to evaluate the satellite orbital decay during geomagnetic storm events (Krauss et al., 2020(Krauss et al., , 2018Oliveira et al., 2019). In this approach, the orbital density is the crucial parameter. The retrieval of the orbital density requires high precise onboard accelerometers as well as the exact satellite parameters such as shape, mass, and materials of satellite planes (March et al., 2019;Sutton et al., 2007). However, the orbital densities are not always available since only a few satellites are equipped with accelerometers. Therefore, the precise orbital decays used in the previous studies were mostly based on CHAMP and GRACE data. Since GPS receivers have become common space-borne instruments, satellite POD data are available for many satellites. The energy-decay of a satellite caused by non-conservative forces could be obtained from the POD data, which could be further used to derive the non-conservative acceleration and the orbital decay (Chen et al., 2012;.  derived the POD-based non-conservative acceleration of CHAMP with 20-min resolution by using an approximate analytic process to represent the time variation gravitational potential to further calculate the energy-decay and thermospheric density. Chen et al. (2012) used a similar analytic process to calculate satellite mean semi-major axis and further the orbital decay. However, the approximate analytic approach could induce large errors in the calculation of the energy-decay, especially when the satellite orbital density is relatively low, as to be discussed later.
In this study, we derived the POD-based satellite orbital decay and decay rate using a numerical integration approach in obtaining the time variation gravitational potential to calculate the energy-decay. The results were compared with those from the previous analytic method to evaluate the different accuracies of both approaches. Subsequently, the integration approach was utilized to study the storm effects on the variations of the satellite orbits. The orbital decays and decay rates of several satellites during two geomagnetic storm events with different magnitudes and under different geophysical conditions were analyzed.

Methods and Data
In this section, the database used in this study would be described first. Then we would introduce our optimized approach, which is also referred to as the POD-based integration approach. The derived satellite orbital decays and decay rates would be compared with the former analytic approach. The results from the accelerometer data would be regarded as the references.

Database
In this study, POD data from GRACE and Swarm constellation were used to analyze the orbital variations during geomagnetic storms. GRACE is a joint project of the Deutsches Zentrum für Luft-und Raumfahrt and the National Aeronautics and Space Administration, which was launched into a near-polar orbit with an initial altitude of about 500 km and operated from March 2002 to the end of 2017. GRACE consists of 2 satellites (A/B) and the twin satellites followed each other on the same orbital track with a distance of about 220 km. Therefore, only GRACE-A data was used in this study. Swarm constellation (A, B, and C) are also polar-orbiting satellites, which were launched on 22 November 2013 and operated by the European Space Agency (ESA). Swarm-A/C flew side-by-side at about 470 km and Swarm-B was located at a higher altitude of about 520 km. LI AND LEI 10.1029/2020SW002664 2 of 14

Description of the Approach
If the spherical Earth central gravity is the only force considered, the Earth-satellite system can be simplified as a two-body problem. In this case, satellite semi-major axis a is a constant and can be calculated from satellite mechanical energy (Montenbruck et al., 2000): where GM is the Earth gravitational constant. E represents satellite mechanical energy, which is the summation of the kinetic and the gravitational potential energy. v denotes the velocity of the satellite, and r is the distance between the satellite and the center of the Earth.
Considering the mechanical energy-decay induced by the non-conservative forces. It should be noted that the non-conservative forces consist of the atmospheric drag, the Solar Radiation Pressure (SRP), terrestrial infrared radiation, the equivalent force caused by the relativistic effect, and so on. Note that the atmospheric drag dominates at low altitudes. As the atmospheric drag decays exponentially with satellite altitude, the atmospheric drag could be on the same order of the SRP at high altitudes (Vielberg & Kusche, 2020). During a certain period, the mechanical energy of a satellite decreases under the effects of the non-conservative forces. According to Equation 1, the energy-decay ΔE between the initial time t 0 and the ending time t could be expressed as the change of the semi-major axis: Where 0 t a represents the semi-major axis of the satellite at the initial time. Due to the non-conservative forces, the semi-major axis changes to t a at the ending time. Consequently, the orbital decay of the satellite, which is defined as the difference between t a and 0 t a , is written as: Since the satellite orbital decay Δa could be calculated from the energy-decay ΔE, the mechanical energy-decay ΔE becomes the key parameter in the process of deriving satellite orbital decay.
In reality, the Earth-satellite system is not an ideal two-body system and the semi-major axis changes with positions. For instance, during the pre-storm time of 19 November 2003, according to Equation 1, GRACE instantaneous semi-major axis varied from 6843 to 6862 km, as shown in Figure 1a. This variation is mainly caused by the non-central gravitational potentials such as the non-spherical gravitational potential, the solar and lunar gravitational potentials. Therefore, we should consider all these conservative potentials and extract the exact energy-decay caused by the non-conservative forces to derive the orbital decay. The mechanical energy of a satellite is rewritten as follows (Jekeli et al., 1999): The terms on the right-hand side of Equation 5 represent the kinetic energy, the Earth gravitational potential V (including the central and the non-spherical gravitational potentials), the central gravitational potentials of the Sun V SUN and the Moon V MOON , respectively. The last term is the integration of time-variation gravitational potential, which is mainly caused by the Earth rotation and the time-variation Earth tides (solid and ocean tide).
LI AND LEI

10.1029/2020SW002664
3 of 14 Similarly to the two-body problem, the mechanical energy-decay ΔE should be calculated according to Equation 5 and subsequently to derive the orbital decay. Based on the precise orbit ephemeris, the kinetic energy of a satellite could be calculated directly, and the gravitational potentials as well as the solar and lunar potentials could be obtained using precise force models. However, the calculation of the time-variation gravitational potential (last term of Equation 5) is complicated, since it is difficult to be expressed explicitly.
Previous studies simplified this term by using a proximate analytic presentation (Chen et al., 2012;Jekeli et al., 1999;, while the simplification brings large relative error in estimating small-scale energy-decay such as the storm-induced energy-decay during solar minimum or satellites with high altitudes. To obtain the mechanical energy-decay, we proposed an optimized approach to calculate the time-variation gravitational potential using numerical integration, which is regarded as the integration method. Considering a short-arc of the satellite orbit, the change of the gravitational potential (including central, non-spherical, and time-variation gravitational potentials) between the initial and ending positions is obtained by calculating the work done by the Earth gravitational force along the satellite trajectory, using numerical integration approach: Where  ECI g is the Earth gravitational force in the Earth-Centered Inertial (ECI) coordinate. To guarantee the accuracy of the calculation, the satellite ephemeris should be interpolated first. In this study, an 8-point Lagrange-piecewise interpolation was utilized. The force models used to calculate satellite mechanical energy are as follow: 120-order EGM2008 gravity field model (Pavlis et al., 2008), the solar and the lunar positions are from DE405 ephemeris, the solid tide and the pole tide models are from IERS2010 (Petit et al., 2010), and the EOT ocean tide model (Rieser et al., 2012). Then the energy-decay caused by the non-conservative forces is expressed as: According to Equation 4, the satellite orbital decay could be derived from the energy-decay. The orbital decay of GRACE caused by the non-conservative forces on 19 November 2003 is shown in Figure 1b. We compared the orbital decay with the corresponding semi-major axis (Figure 1a) on the same day. It is found that the orbital decay of GRACE was about 20 m during one day, which was only one-thousandth as compared with the variation of the semi-major axis due to the conservative potentials (30 km). Overall, our approach is enabling to extract the tiny orbital variation of a satellite due to the non-conservative forces from the largescale variation of semi-major axis.
As mentioned above, previous studies used a simplified approach to calculate the time-variation gravitational potential, in which they assumed that the integration of the time-variation gravitational potential is equal to the earth rotation potential: where ω is the rotation speed of the Earth; x and y represent the location of the satellite in the ECI coordinate. With this approximation, the mechanical energy could be expressed in the analytic form. Therefore, the energy-decay ΔE is obtained by calculating the energy difference between the initial and ending position of the arc (ΔE = E t −E 0 ) without numerical integration (Chen et al., 2012;Jekeli et al., 1999;. This approach is referred to as the analytic method. This approach has the advantages of less computer time and convenient application. However, the approximation of the time-variation gravitational potential in Equation 9 would cause large uncertainty, especially when the satellite orbital density is relatively low, as shown later. In the next section, we would compare the two methods to illustrate their accuracies.

Method Verification
As mentioned above, the satellite energy-decays, associated with the non-conservative forces, are dominated by the atmospheric drag at low altitude and solar maximum. In this subsection, the energy-decays from both the analytic and integration methods are compared. Note that the energy-decay derived from the accelerometer data is used as the reference result, since GRACE accelerometer could measure the alongtrack acceleration with high accuracy, which could be applied to calculate very precise energy-decay. In this study, the calibrated accelerometer data was from Sutton et al. (2007) during 2002-2009and from Li and Lei (2021 during 2010-2016. The energy-decay from the accelerometer data is written as: Where  Acc represents satellite along-track acceleration from the accelerometer,  v is the satellite velocity, and dt is the time interval. Figure 2 shows a typical instance of the POD-based energy-decays from the analytic (gray line) and integration (blue line) methods on 19 November 2003. The energy-decay from the accelerometer data (red line) is shown as the reference. The time interval was 50 s and this 50-s time resolution is much higher than the 20-min resolution from the analytic method . It is seen that the energy-decay from the integration method is consistent with that from the accelerometer data. This demonstrates that the POD-based energy-decay from the integration method has the same accuracy with the accelerometer data. Note that the little difference between the two methods could be caused by the uncertainty of the POD data. However, LI AND LEI 10.1029/2020SW002664 5 of 14 large amplitude oscillations exist in the result from the previous analytic method, which is mainly caused by the simplification of the time-variation gravitational potential (Equation 9).
Next, we compared the POD-based orbital decays and decay rates from the analytic (gray line) and integration (blue line) methods. Since the accelerometer data is not available for Swarm, the results from the integration method are considered as the reference. During the pre-storm time of another case on 7 September 2017, GRACE flew at about 338 km with an orbital density of about 3 × 10 −12 kg/m 3 . Meanwhile, Swarm-B was at 515 km and the orbital density was less than 3 × 10 −13 kg/m 3 . Figure 3 shows the orbital decays and decay rates from those two methods at different orbital density levels. As illustrated in Figure 3a, the orbital decay rate of GRACE from the analytic method has a larger variation than the integration one. The difference of the decay rates between these two methods is about 100 m/d and the relative difference is about 50%. Although the difference of the decay rates between two methods is large, the orbital decays appear to be similar with about 10% relative difference (Figure 3b). This is due to the fact that the orbital decay is the integration of the decay rate and the large deviation of the decay rate may be partly canceled out. For Swarm-B at a higher altitude and with lower orbital density, the orbital decay rate is about 20 m/d, which is only one-tenth of that of GRACE. Therefore, although the absolute difference of the decay rates between the two methods is still about 100 m/d, the relative difference could be as large as more than 400% as shown in Figure 3c. Considering the orbital decay of Swarm-B, the orbital decay from the analytic method even turns to be positive as shown by Figure 3d, which is obviously incorrect. The comparisons of the orbital decays and the decay rates between the two methods suggest that the integration method has higher accuracy in calculating satellite orbital decay and the decay rate based on the POD data. However, as inferred by Figure 3b, the analytic method may also be used to estimate the orbital decay approximately when the orbital density is relatively large. LI AND LEI 10.1029/2020SW002664 6 of 14 The remaining question is when the analytic process could be used to estimate the orbital decay and how large is the deviation. We computed the daily orbital decays of GRACE from both analytic and integration methods from 2010 to 2016 (every 5 days) and compared with that from the accelerometer data. From 2010 to 2016, GRACE altitude decreased from about 480 to 360 km and the orbital densities varied from 1 × 10 −13 to 5 × 10 −12 kg/m 3 . Figure 4 shows histogram of the relative differences of the orbital decays from the two methods with respect to that from the accelerometer data under various orbital densities. The orbital densities used for statistical grouping are from Li and Lei (2021). It is seen that the relative deviations of the analytic method (red box) are more than 25% for the majority of the cases when the orbital densities are less than 5 × 10 −13 kg/m 3 . The deviations generally decrease to less than 25% when the orbital densities increase from 1 × 10 −12 to 5 × 10 −12 kg/m 3 . For the integration method (cyan box) the relative differences are mostly less than 10%. Overall, our integration method is more accurate than the previous analytic method. On the other hand, the analytic method can be used to estimate the orbital decay when the orbital density is larger than 10 −12 kg/m 3 . Therefore, the POD-process is an alternative way to derive satellite orbital decay and decay rate when the accelerometer data is not available. In the following section, two geomagnetic storms would be studied to evaluate the storm-effects on satellite orbits.

November 2003 Storm
The storm effect on satellite orbit was analyzed during the well-known November 2003 geomagnetic storm. The geomagnetic storm was one of the largest storm events since March 1989, and was triggered by the interaction of several magnetic clouds and terrestrial magnetosphere (Yushkov et al., 2005). Figure 5 shows the interplanetary and geophysical parameters during the storm. The November 2003 storm occurred during solar maximum with F10.7 at a high level of more than 140. During the main phase, the interplanetary magnetic field (IMF) B z component dropped to −47 nT and B y turned to be 49 nT. The corresponding AE and Kp indices increased to about 2000 nT and 9, respectively. The Dst minimum was −422 nT. Figure 6 shows the orbital decay and decay rate of GRACE during this event. Figure 6a   rates of GRACE from the analytic (gray line) and integration (blue line) methods based on POD data. The decay rate from the accelerometer data (red line) is shown for comparison. Again, the orbital decay rate from the integration method shows good agreement with that from the accelerometer data. The orbital oscillation is associated with the diurnal and latitudinal variations of the thermospheric density. However, large deviations are found in the result from the analytic method during the quiet time, as compared with the accelerometer data. During the main phase of the geomagnetic storm, when the thermospheric density was enhanced to more than 5 × 10 −12 kg/m 3 , the decay rate from the analytic method was generally consistent with that from the accelerometer data.
The orbital decay rate of GRACE during the quiet time was about 21 m/d as shown by the cyan line. During the main phase of the November 2003 storm, the orbital decay rate increased to more than 150 m/d, which was about 800% as compared with the quiet time value. This is expected since the thermospheric mass density increased by 300%-800% during the storm time due to the huge energy injection from the magnetosphere into the upper atmosphere (Bruinsma et al., 2006). The orbital change Δa due to the storm-induced density enhancement is given as (Chen et al., 2014): where  and  b represent the storm-time and background densities, respectively. a is the mean semi-major axis. C D , A, and m represent satellite drag coefficient, cross-section area, and mass respectively. Figure 6b shows the orbital decays from those three methods. To evaluate the extra orbital decay caused by the enhanced thermospheric density due to the geomagnetic storm, the background orbital decay has been removed, which is defined as the storm-induced orbital decay. The background orbital decay is regarded as the mean orbital decay during the first half day of 19 November when the geophysical condition was relatively quiet. As expected, the storm-induced orbital decays from the integration method (blue line) and the accelerometer data (red line) are almost the same. Since the orbital density was more than 10 −12 kg/m 3 , the analytic method also performs well and the relative deviation is about 10%, which is consistent with the conclusion in Section 2.3. Note that the storm-induced orbital decay was about 75 m, which contributed about 50% of the total decay during the selected 4 days (162 m). That is, although the geomagnetic storm lasted for about 1 day, the storm-induced orbital decay could be comparable with the 4-days orbital decay during the quiet time.

September 2017 Storm
The September 2017 geomagnetic storm was one of the largest events during the solar minimum of the 24th solar cycle. In this period, GRACE was at the end of its mission with an altitude of about 338 km. Swarm-A and Swarm-B flew at about 451 and 515 km, respectively. The simultaneous observations of the three satellites would reveal the various storm effects on satellites at different altitudes. Since the accelerometer data of Swarm-A/B were not available, the POD-based integration method was applied to investigate the storm effects on those satellites.
LI AND LEI 10.1029/2020SW002664 9 of 14  During this event, the orbital decays and decay rates of Swarm-A/B and GRACE are displayed in Figure 8. As shown in Figure 8a, the orbital decay rate of Swarm-B (∼515 km) increased from 7.7 m/d to more than 20 m/d from the quiet time to the storm time. Swarm-A flew at a lower altitude (∼451 km), and the decay rate increased from 20 to about 60 m/d (Figure 8b), which was greater than that of Swarm-B. With regard to GRACE, which flew at the lowest altitude (∼338 km), the decay rate changed from about 155 to 400 m/d (Figure 8c). Note that the positive values of Swarm-B decay rates were induced by the uncertainty of the POD data when the mass density was extremely low. The total orbital decays ( orbital decays were 7, 20, and 152 meters as shown by the solid lines. It is found that the storm induced a quarter of the total orbital decays, which appeared to be the same for those three satellites at different altitudes. It should be noted that the orbital decay of GRACE during the solar minimum storm (September 2017, Figure 8f) was much larger than that during the high solar activity storm (November 2003, Figure 6b), since the altitude of GRACE was about 150 km lower in September 2017 and consequently the orbital density was much higher. Figure 9 shows the relative changes of the orbital-averaged decay rates between the storm and the quiet time. The relative changes were about 100% for GRACE and about 150% for Swarm-A/B, since Swarm-A/B flew at higher altitudes than GRACE, which is associated with the expansion of the thermosphere and an increase of the scale height (Lei et al., 2010).
LI AND LEI 10.1029/2020SW002664 11 of 14 According to Equation 11, the relative change of the orbital decay rate represents the enhancement of the storm-time thermospheric density. It should be noted that the high-cadence thermospheric mass density could be further retrieved from the POD-based non-conservative force by precisely considering the solar radiation pressure, terrestrial infrared radiation, and the relativistic effect. The detailed process can be found in our recent work (Li & Lei, 2021).

Summary
In this study, it was found that the analytic presentation of the time-variation gravitational potential from Equation 9 could induce large uncertainties in the calculation of satellite orbital decay rate, and subsequently give unreliable satellite orbital decay when the orbital density is smaller (cf. less than 10 −12 kg/cm 3 ). Alternatively, the integration presentation of the time variation gravitational potential was used to calculate the orbital decay rate with a time resolution of less than 1 min, which is much higher than the previous analytic method. The accuracy of the corresponding orbital decay from the integration method was also high with relative deviations of less than 10% as compared with the accelerometer-reference. Subsequently, the November 2003 and September 2017 geomagnetic storms were considered to investigate the storm-induced drag effects on satellite orbits based on the integration method. In the severe November 2003 storm, the storm-induced orbital decay rate increased rapidly to ∼800% with respect to the quiet-time reference. The storm-induced orbital decay was comparable with the 4-days orbital decay during the quiet time. In the moderate September 2017 storm, the orbital variations of Swarm-A/B and GRACE satellites with different altitudes were analyzed. The orbital decay rates increased from 100% to 150% depending on satellite altitudes. The storm-induced orbital decays were also significant and contributed about one fourth of the total decays during the 3-days interval.