Distinct thermospheric mass density variations following the September 2017 geomagnetic storm from GRACE and Swarm

Distinct thermospheric mass density variations following the September 2017 geomagnetic storm from GRACE and Swarm Liangliang Yuan 1, 2, Shuanggen Jin 1, Andres Calabia 3 1 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China 2 University of Chinese Academy of Sciences, Beijing 100049, China 3 Ann and J.H. Smead Department of Aerospace Engineering Sciences, University of Colorado, Boulder 80309, Colorado, USA Abstract: Geomorphologies and mechanisms of thermospheric mass density variations caused by geomagnetic storms are still a challenge due to limited observations and imprecise models. Recently, GRACE and SWARM satellites are unable to provide accurate and stable accelerometer measurements but new techniques to derive density estimates from Precise Orbit Determination (POD) are a potential tool to derive thermospheric density estimates. In this paper, we employ POD-based techniques to derive GRACE-A and SWARM-A mass density estimates and investigate the variability during the September 2017 geomagnetic storm. This unique opportunity can investigate the storm-time mass density variability in altitude because on September 2017 both missions have similar Local Solar Time (LST) location. We analyze the variability in latitude and altitude for this two-peak shaped storm, and our results reveal that the first peak affects upper altitudes more severely than the second peak, which enhances the lower thermosphere relatively. Furthermore, we investigate the coupling mechanism between thermosphere and ionosphere through the correlation between the mass density variability and the ionospheric [O/N2] ratio variability for the same storm. The Global Ultraviolet Imager (GUVI) [O/N2] ratio observations show for the first geomagnetic peak to be stronger in the southern hemisphere than in the northern hemisphere, and the contrary for the second geomagnetic peak. This feature shows a good consistence with our mass density variability with latitude. Keywords: Thermospheric density; POD; [O/N2]; Geomagnetic storm. 1. Introduction Main dynamic disturbances at Low Earth Orbit (LEO) are caused by changes of thermospheric mass density estates due to an atmospheric expansion caused by variations of solar activity. Therefore, the exact modeling of thermospheric mass density variations is a very important role in Precise Orbit Determination (POD) of debris, unmanned objects, and active satellites. In the thermosphere, the coupling between ion and neutral particles is mainly driven by the solar activity and Earth’s magnetic field, where the neutral density is very low and the dynamics are rather driven by the Extreme Ultra Violet (EUV) heating instead by the intermolecular interactions [Forbes et al., 1990]. In addition, rapid abrupt changes of thermospheric density and composition occur during geomagnetic storms [Forbes et al., 1996], and the involved geophysical processes are still not well understood nor modeled as required for POD. In the recent decades, accelerometer and POD have been the two major approaches to derive thermospheric total mass density estimates [Bruinsma et al., 2004; Mclaughlin et al., 2013;; Calabia and Jin, 2017], and widely employed to investigate thermospheric variations driven by Joule heating and air upwelling during geomagnetic storms [Forbes et al., 1996; Bruinsma et al., 2006; Calabia et al., 2016a]. Thermospheric mass density responses to geomagnetic storms are more apparent in the high-latitude regions, near the cusps, where the solar wind can access easily into the upper atmosphere through the Earth´s magnetic reconnection. Furthermore, in situ abrupt heating in the cusp can produce equator-ward gravity waves [Bruinsma et al., 2006]. During a magnetic storm, the energy in the thermosphere can be transferred from high to low latitudes through both gravity waves and meridional circulation [Richmond et al., 1979]. Moreover, disturbances in the high-latitude region propagate equator-ward through atmosphere circulation, reaching in few hours in the middle and lower latitudes Calabia et al. [2017] have recently provided the correlation index and time-delay for the most representative geomagnetic indices, as well as the time to reach the equatorial region. Studying the thermospheric density variability during geomagnetic storms is of high importance for Magnetosphere-Ionosphere-Thermosphere (MIT) coupling research and POD applications and for the improvement of current models used in POD. For instance, previous research (e.g., Calabia and Jin [2016b], Calabia and Jin [2017]) has proven that the latest empirical model NRLMSISE00 [Picone et al., 2002], standard used for POD, is currently unable to accurately reproduce and predict the mean value and LST amplitude of variability of the actual thermospheric density estates during geomagnetic storms. Burns et al. [2004] inferred that density variability is depends on the background quiet-time state and vary from season to season. Their results showed that the density variability in winter is higher than that in the summer under the same conditions. This feature is attributed to differences in solar heating. Moreover, based on Challenging Minisatellite Payload (CHAMP) accelerometer-based estimates, the results in Liu et al. [2005] showed small differences between the day-side and the night-side relative percentage changes, and that the response of thermospheric density during storm time is not similar from event to event. The westward equatorial flow following the diurnal variation of EUV heating is especially obvious in the night side, but it is difficult to investigate due to the large intervals between CHAMP recurrence sampling, and the convection driven by meridional and zonal winds is not sensitive to the current observations [Sutton et al., 2005]. Sutton et al. [2005] showed density enhancements of about 300–800% during the geomagnetic storm of November 2003.The time-delay with respect proxies showed shorter at high latitudes than at lower latitudes, being the later about 4-hour. Day side densities at 410 Km altitude in the high-latitude Southern hemisphere in Summer (i.e., -72° geomagnetic altitude) at noon LST are very responsive to increases in solar wind dynamic pressure, and even during periods when both Bz and By are near zero or positive [Bruinsma et al., 2006]. Finally, Jin et al. [2017] suggested that there is a strong relationship between thermospheric mass density variability and ionospheric [O]/[N2] variance, but more research studies are required to reveal the complex processes of coupling between ionosphere and neutrals in the upper atmosphere. This paper focuses on the detailed responses of thermospheric POD-based mass density during the geomagnetic storm of September 2017, and the relational variability with the ionospheric [O]/[N2] composition fluctuation. The following section describes the data and methods employed in this study. Section 3 analyzes in detail our analysis results for latitude and altitude variability dependences as well as the relational aspects with space weather indices and changes in ionospheric [O]/[N2] composition. 2. Data and Methods 2.1 Geomagnetic indices On September 8th, 2017, two peaks were detected in space weather indices with similar magnitude. The first peak reached the highest magnitude at 0 h UT with a Dst value of -140 nT, the latter peak reached a Dst value of -120 nT at 12 h UT. Figure 1 shows the time series of Space Weather measurements, including the interplanetary magnetic field (IMF) Bx, By, Bz, the Field Magnitude Average (FMA), the Plasma Speed (PS), and the Electric Field (EF), which have been downloaded from https://spdf.gsfc.nasa.gov/. The Dst index and the three-hour Ap index can be found at http://isgi.unistra.fr/data_download.php. The first and second Storm Sudden Commencement (SSC) [Araki, 1994] started at 21 h UT on 7 September 2017 and at 10 h UT on 8 September 2017, respectively. The minimum Bz value of IMF was -31 nT and -15 nT for each peak. Figure 1. Space weather indices during the 7-9 September 2017 geomagnetic storm events. SSCs occurred at 21 h UT UT on 7 September 2017 and at 10 h UT on 8 September 2017. 2.2 Thermospheric mass density retrieval In this work, POD-based thermospheric mass density estimates are derived from the Level 1B data of GRACE-A and SWARM-A missions, including accelerometer, navigation, thruster, and starcamera measurements, among others. GRACE's Level 1B data can be downloaded from http://isdcold.gfz-potsdam.de/ in binary big-endian format, and SWARM’s Level 1B data can be downloaded from ftp://swarm-diss.eo.esa.int in CDF (Common Data Format). More details and metadata can be found in the respective links. We employ the method developed in Calabia et al. [2015] and Calabia and Jin [2017] for accelerometer calibration and POD-based thermospheric mass density retrieval. Theoretical bases, specific details, and MATLAB® scripts can be found in Calabia [2017]. The method is based on the degravitation of numerically differentiated precise orbit velocities, where a first step employs the arcto-chord method in a piecewise interpolation scheme to minimize the arc-to-chord error committed in the numerical differentiation. In this scheme, POD-based total accelerations are derived from Precise Orbit velocities as follows: 𝑟󰇘𝑡0 = 𝑙𝑖𝑚 Δ 𝑡→0 𝑟𝑡2−2𝑟𝑡0+𝑟𝑡(−2) ( Δ 𝑡)2, where 𝑡(𝑖 + 1) − 𝑡(𝑖) = Δ 𝑡. (1) The following step is to de-gravitate or, in other words, remove from POD-based total accelerations the effects of the time-varying gravity. The time-varying gravity at the satellite location is accurately estimated through models, where the static average gravitational field is combined with the variable contributions, including the trends of the low-degree Stokes coefficients, solid tides, permanent tides, ocean tides, polar tides and others. The standard form to compute the variable Earth's gravitational potential is given by: 𝑉(𝑟, 𝜑, 𝜆) = 𝐺𝑀𝐸𝑎𝑟𝑡ℎ 𝑟∑(𝑎𝑒 𝑟) 𝑁 𝑛=0 𝑛∑[𝐶𝑛𝑚 𝑛 𝑚==0 𝑐𝑜𝑠( 𝑚𝜆) + 𝑆𝑛𝑚 𝑠𝑖𝑛( 𝑚𝜆)]𝑃 𝑛𝑚 𝑠𝑖𝑛( 𝜑) (2) The gravity model is calculated at the ITRS (International Terrestrial Reference System), but the non-gravitational accelerations are given at the SBS (Satellite Body System). Therefore, the timevariable gravity in the ITRS has to be transformed to the ICRF (International Celestial Reference System) through the star camera quaternion measurements. Moreover, the residual accelerations in the ICRF have to be transformed to the SBS (Satellite Body System) to remove the three-body perturbation and the relativistic effects. Non-gravitational accelerations are finally given in the SBS frame to be compared with accelerometer measurements if available. In order to derive the pure aerodynamic accelerations, direct Solar radiation pressure (direct and reflected) and Earth's albedo radiation need to be removed from the non-gravitational accelerations obtained in the previous step. Detailed theoretical basis and MATLAB algorithms for radiation pressure removal are provided in Calabia et al [2017]. The final sted is to obtain the air mass density along the orbital path through the aerodynamic drag equation 𝐹𝐷=1 2𝐶𝐴𝜌𝑣𝑟2 , (3) In this study, we set the drag coefficient C to 2.2. A is the cross-sectional area perpendicular to vr , ρ is the mass density, and vr is the relative velocity of the atmosphere with respect to the spacecraft, which includes the co-rotating atmosphere vC and horizontal winds vw. 𝑣𝑟= −𝑟󰇗𝑠𝑎𝑡 + 𝑣𝑐+ 𝑣𝑤 (4) 3. Results and Analyses 3.1 POD-based non-gravitational accelerations GRACE´s POD-based non-gravitational accelerations are first contrasted to accelerometer measurements. We employ the POD solution to calibrate accelerometer measurements. Figure 2 shows the comparison of our POD-based with the calibrated-accelerometer non-gravitational accelerations from GRACE-A on 8 September 2017. On the other side, in Figure 3 we compare SWARM´s POD-based non-gravitational accelerations computed through our numerical derivation with the ESA´s POD-based non-gravitational accelerations derived using a Kalman Filter on a least-squares POD scheme with the reduceddynamic solution as a reference. The low performance of SWARM-A accelerometers and the missing ESA´s POD-based non-gravitational accelerations 8 September 2017 forced us to compare our results for SWARM-C, and for a different date, on 15 August 2014. Figure 2. The comparison of our POD-based and calibrated-accelerometer non-gravitational accelerations from GRACE-A on September 8th , 2017. Figure 3. The comparison of POD-based non-gravitational accelerations from our results and from the official ESA products for SWARM-C on August 15th, 2014. 3.2 Thermospheric density variance The algorithm to derive mass density along the orbital path of GRACE and SWARM (Equation 3) has been used to estimate mass density estates from the POD-based non-gravitational accelerations during the September 8th , 2017 Geomagnetic Storm. We employ the surrounding quiet period to evaluate percentage changes during the storm, but accelerometer measurements of GRACE were not available during the adjacent periods. The accelerometer measurements of SWARM-C during this storm are also not provided by ESA. However, POD-based non-gravitational accelerations provide us the opportunity to derive and investigate mass density variations. Figure 4 pictures the mass density retrieved form POD, accelerometer, and from the NRLMSISE00 model. Clear differences are seen between the model and measurements, indicating the lack of detail from the model to reproduce mass density variations during the storm. Figure 4. The comparison of accelerometer-based thermospheric mass density from GRACE-A, our POD-based results, and the NRLMSISE00 estimates on September 8th, 2017. As it can be seen from the change in the IMF parameters in Figure 1, the peak of the two geomagnetic storms is located at 0 h UT and at 12 h UT on the 8th, respectively. SWARM and GRACE missions have a quasi-polar orbit and their latitude range from 80 ° S to 80 ° N, SWARM's Local Solar Time (LST) is located at about 10 h LST at the ascending orbit and 22 h LST at the descending orbit, and therefore divide our analysis into these two parts. GRACE’s LST at this period is located at about 21:30 h LST at the ascending orbit and about 9:30 h LST at the descending orbit, so this is an ideal configuration since SWARM has a similar LST location. Figures 5 and 7 pictures the mass densities along the orbital path of GRACE and SWARM, respectively. Figures 6 and 8 presents the same figures with normalized values at altitudes of 350 Km and 450 Km, respectively. The averaged atmospheric density at the day-side (10 h LST) shows higher values than the values at the night-side (22 h LST). It can be seen that two peaks are located at about 0 and 12 h UT during the geomagnetic storm, the average atmospheric density on the day side is significantly higher than that on the night side for both missions, and the short-term variations well correlated with the geomagnetic indices. Figure 5. POD-based thermospheric density variability from (a) ascending (21:30 h LST) and (b) descending orbits (9:30 h LST), from GRACE-A at orbital altitude during the geomagnetic storm of September 8th, 2017. Note time increases from right to left. Figure 6. Panels from Figure 5 normalized at 350 Km altitude. NRLMSISE00 is included in (b) and (d) for comparison. (a) and (b) are panels for day-time, and (c) and (d) panels for the night-time. Figure 7. POD-based thermospheric density variability from (a) ascending (10:00 h LST) and (b) descending orbits (22:00 h LST), from SWARM-A at orbital altitude during the geomagnetic storm of September 8th, 2017. Note time increases from right to left. Figure 8. Panels from Figure 7 normalized at 450 Km altitude. NRLMSISE00 is included in (b) and (d) for comparison. (a) and (b) are panels for day-time, and (c) and (d) panels for the night-time. 3.2 Density variability with latitude During the first phase of the geomagnetic storm, the GRACE-A results shown in Figure 6 reveal a slightly bigger response in the Northern hemisphere than in the Southern hemisphere. The global atmospheric density does not show an increase at this stage. On the other side, SWARM pictures a slightly bigger response in the Southern hemisphere. It is interesting to note that at Equinox periods density values are enhanced at the Southern latitude, suggesting a higher cusp precipitation probably related to the dipole-tilt angle variation [Calabia et al, 2017]. Moreover, from Figure 6 it seems that this effect is less pronounced at a lower altitude. NRLMSISE00 is unable to reproduce this variability in latitude. During the second geomagnetic peak, the global atmospheric density of SWARM rapidly increases. Figure 8 shows that the Southern hemisphere density values are bigger than those for the Northern hemisphere. This enhancement increases and rapidly covers more latitude range, especially in the daytime. On the other side for GRACE, Figure 6 shows equilibrium between both hemispheres.