# Comparison of Empirical and Theoretical Models of the Thermospheric Density Enhancement During the 3–4 February 2022 Geomagnetic Storm

Jianhui He, Elvira Astafyeva, Xinan Yue, Nicholas Pedatella, Dong Lin, Timothy Fuller-Rowell, Mariangel Fedrizzi, Mihail Codrescu, Eelco Doornbos, Christian Siemes, et al.

► To cite this version:

Jianhui He, Elvira Astafyeva, Xinan Yue, Nicholas Pedatella, Dong Lin, et al.. Comparison of Empirical and Theoretical Models of the Thermospheric Density Enhancement During the 3–4 February 2022 Geomagnetic Storm. Space Weather: The International Journal of Research and Applications, 2023, 21 (9), 10.1029/2023SW003521. [hal-04248456]

HAL Id: hal-04248456

https://hal.science/hal-04248456v1

Submitted on 18 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Space Weather®

## RESEARCH ARTICLE

10.1029/2023SW003521

**Key Points:**

* Observations by Swarm-A (~438 km) and Gravity Recovery and Climate Experiment Follow-On (~505 km) missions show 110% and 120% storm-time increase in the thermospheric density

* Both empirical and numerical models captured the observed storm-time density enhancements, however simulated results differ by up to 70%

* Simulations show that the storm-time density increase reached 35%–60% at the orbit altitude of the Starlink satellites, that is, 200–300 km

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

E. Astafyeva,
astafyeva@ipgp.fr

**Citation:**

He, J., Astafyeva, E., Yue, X., Pedatella, N. M., Lin, D., Fuller-Rowell, T. J., et al. (2023). Comparison of empirical and theoretical models of the thermospheric density enhancement during the 3-4 February 2022 geomagnetic storm. *Space Weather*, 21, e2023SW003521. https://doi.org/10.1029/2023SW003521

Received 2 APR 2023
Accepted 2 AUG 2023

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial-NoDerivs
License, which permits use and
distribution in any medium, provided the
original work is properly cited, the use is
non-commercial and no modifications or
adaptations are made.

# Comparison of Empirical and Theoretical Models of the Thermospheric Density Enhancement During the 3–4 February 2022 Geomagnetic Storm

<figure><img src="image_2.png" /></figure>

Jianhui He¹,²,³ Ⓓ, Elvira Astafyeva² Ⓓ, Xinan Yue¹,³ Ⓓ, Nicholas M. Pedatella⁴ Ⓓ, Dong Lin⁴ Ⓓ, Timothy J. Fuller-Rowell⁵,⁶ Ⓓ, Mariangel Fedrizzi⁵,⁶ Ⓓ, Mihail Codrescu⁵,⁶ Ⓓ, Eelco Doornbos⁷ Ⓓ, Christian Siemes⁸ Ⓓ, Sean Bruinsma⁹ Ⓓ, Frederic Pitout¹⁰ Ⓓ, and Adam Kubaryk⁵,⁶

¹Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China, ²Université Paris Cité, Institut de Physique du Globe de Paris (IPGP), CNRS UMR 7154, Paris, France, ³College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China, ⁴High Altitude Observatory, National Center for Atmospheric Research, Boulder, CO, USA, ⁵NOAA Space Weather Prediction Center, Boulder, CO, USA, ⁶Cooperative Institute for Research in Environmental Sciences, University of Colorado Boulder, Boulder, CO, USA, ⁷Koninklijk Nederlands Meteorologisch Instituut, De Bilt, The Netherlands, ⁸Delft University of Technology, Delft, The Netherlands, ⁹GET/CNES, Space Geodesy Office, Toulouse, France, ¹⁰Institut de Recherche en Astrophysique et Planétologie (IRAP), Toulouse, France

**Abstract** On 3 February 2022, at 18:13 UTC, SpaceX launched and a short time later deployed 49 Starlink satellites at an orbit altitude between 210 and 320 km. The satellites were meant to be further raised to 550 km. However, the deployment took place during the main phase of a moderate geomagnetic storm, and another moderate storm occurred on the next day. The resulting increase in atmospheric drag led to 38 out of the 49 satellites reentering the atmosphere in the following days. In this work, we use both observations and simulations to perform a detailed investigation of the thermospheric conditions during this storm. Observations at higher altitudes, by Swarm-A (~438 km, 09/21 Local Time [LT]) and the Gravity Recovery and Climate Experiment Follow-On (~505 km, 06/18 LT) missions show that during the main phase of the storms the neutral mass density increased by 110% and 120%, respectively. The storm-time enhancement extended to middle and low latitudes and was stronger in the northern hemisphere. To further investigate the thermospheric variations, we used six empirical and first-principle numerical models. We found the models captured the upper and lower thermosphere changes, however, their simulated density enhancements differ by up to 70%. Further, the models showed that at the low orbital altitudes of the Starlink satellites (i.e., 200–300 km) the global averaged storm-time density enhancement reached up to ~35%–60%. Although such storm effects are far from the largest, they seem to be responsible for the reentry of the 38 satellites.

**Plain Language Summary** On 3 February 2022, at 18:13 UTC, the 38 out of the 49 Starlink satellites reentered the atmosphere during a moderate geomagnetic storm that occurred on current and the next day. The storm-time enhanced thermospheric density is the cause of an enhanced atmospheric drag effect on the LEO satellites and space debris. In this work, the whole thermosphere neutral mass density conditions are investigated in detail by a combination of the observations and the model simulations. The results show that ~110%–120% storm-time increase in the upper thermospheric density (e.g., ~400–500 km), and the storm-time density increase reached ~35%–60% at the orbit altitude of the Starlink satellites, that is, 200–300 km. Although such storm effects are far from the largest, they seem to be responsible for the reentry of the 38 satellites.

## 1. Introduction

During geomagnetic storms, a large amount of energy in the form of Joule heating and particle precipitation is deposited in the high-latitude thermosphere, resulting in heating and expansion of the neutral gas (Mayr et al., 1978; Prölss, 1980). Storm-time high-latitude heating changes the pressure gradients in the *E* and *F* regions between the equatorial region and the poles, resulting in meridional winds perturbations that can quickly propagate to lower latitudes in a few hours. Then, the overall modification of the thermospheric composition and increase of the neutral mass density occur globally (Balan et al., 2011; Fuller-Rowell et al., 1996; Li et al., 2019; H. Liu & Lühr, 2005; Roble et al., 1982; Wang et al., 2021). The enhanced thermospheric density, in turn, is the
cause of an enhanced atmospheric drag effect on the low-Earth-orbit (LEO) satellites and space debris (Anderson et al., 2009; Chen et al., 2012; Wright, 2007). Even during geomagnetically quiet days, the long-term accumulated drag effect can also decrease the orbit altitude of LEO spacecraft due to the drag, which is proportional to the spacecraft area-over-mass ratio and the velocity of the spacecraft squared, as well as the neutral mass density (Emmert, 2015; Marcos, 2006; Sutton, 2018). The capability to accurately forecast variations of the thermospheric neutral density is crucial to the operations of spacecraft in LEO, notably for orbit maintenance and collision avoidance (Brandt et al., 2020; Bussy-Virat et al., 2018; Mehta et al., 2019; Sutton et al., 2021).

The storm-time thermosphere behavior and variations can be investigated by using the thermospheric neutral mass density data set from LEO satellites, for instance, the Challenging Mini-satellite Payload (CHAMP) (e.g., Reigber et al., 2002), Gravity Recovery and Climate Experiment (GRACE) and the GRACE follow-on mission (GRACE-FO) (Tapley et al., 2004) or the Swarm constellation mission (van den IJssel & Visser, 2007). The thermospheric mass density can be calculated from accelerometers onboard LEO satellites (Bruinsma et al., 2004; Doornbos, 2012; Villain, 1980), or from global navigation satellite system data used for precise orbit determination (POD) onboard LEO satellites (van den IJssel & Visser, 2007). The satellite-based thermospheric data helped to improve our understanding of the thermospheric state during storm times. Previous studies show that the relationship between the high-latitude input (e.g., Joule heating) and the thermosphere density enhancement driven by the storm is rather complex (Deng et al., 2011; Huang et al., 2012; Knipp et al., 2004, 2011; Lei et al., 2011; Lu et al., 2016). It has been well demonstrated that the neutral density increases at high latitudes shortly after the storm onset, generating large scale gravity waves, or traveling atmospheric disturbances (TAD), which take 4–6 hr to reach mid- and low latitudes (Bruinsma & Forbes, 2007; Oliveira et al., 2017; Zesta & Oliveira, 2019). Rapid variations of the heating cause oscillations of wind and temperature as well as large variations of the neutral composition and density. The observed storm-time variations of the thermospheric neutral mass density often show hemispheric/seasonal asymmetry, which is attributed to the prevailing summer-to-winter circulation (Astafyeva et al., 2017, 2020; Fuller-Rowell et al., 1996; H. Liu & Lühr, 2005; Sutton et al., 2005). Specifically, an easier equatorward propagation of the neutral density anomaly occurs in the summer hemisphere due to the superposition of the background and storm-driven circulations. However, the opposed circulation in the winter hemisphere hampers the propagation of the high-latitude perturbations toward low latitudes. In addition to the above-mentioned LEO satellites that measure density-related changes of satellite acceleration in situ, nowadays, the Global-scale Observations of the Limb and Disk (GOLD), the Ionospheric Connection Explorer (ICON) are also used to investigate the lower thermosphere responses during the storm times by remote sensing (Aa et al., 2022; Cai et al., 2021; Harding et al., 2022; Laskar et al., 2021b). However, because of the limitations of in situ LEO observations, notably due to the fixed orientation of the orbit plane with respect to the Sun for any given storm, it is not possible to infer the global storm-time response of the thermosphere solely based on observations.

Over the last decades, several general circulation models have been developed. General circulation models are based on the Navier-Stokes equation to self-consistently solve the density, velocity, and temperature of neutral components, and can reasonably provide time-dependent values for each parameter. These models include coupled ionosphere and thermosphere models or whole atmosphere models, such as the Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM, Richmond et al., 1992) and the Whole Atmosphere Community Climate Model with thermosphere and ionosphere eXtension (WACCM-X, H.-L. Liu et al., 2018), etc. In addition, several empirical thermosphere models have been developed based on a variety of data, including the Drag Temperature Model (DTM, Bruinsma, 2015), and the Mass Spectrometer and Incoherent Scatter model (MSIS, Emmert et al., 2020). Some of these empirical models are employed in (operational) orbit computation because of their easy implementation and use.

On 3 February 2022, at 18:13 UTC, a SpaceX Falcon 9 rocket launched 49 Starlink internet satellites from the Kennedy Space Center in Florida. Space-X initially deployed 49 Starlink satellites with elliptical orbits between 210 and 320 km altitude. The Local Time (LT) sampling from the ascending orbits at the equator were at ~12:45 LT. The satellites were meant to be further raised to 550 km. However, the launch and the satellite deployment took place during the main phase of a moderate geomagnetic storm that caused an increase in atmospheric drag. The official communication from SpaceX stated that “onboard GPS suggested the escalation speed and severity of the storm caused atmospheric drag to increase up to 50% higher than during previous launches” (https://www.spacex.com/updates/ accessed on 8 February 2022). SpaceX tried to put the satellites into a “safe mode,” turning them to fly edge-on to minimize drag. However, the second phase of the storm hit and further
increased the drag, causing 38 out of 49 satellites to re-enter in the following days. The remaining 11 satellites were successfully maneuvered into their operational orbits after these events.

Several research groups performed initial studies of the thermospheric conditions during this storm event (Dang et al., 2022; T.-W. Fang et al., 2022; Kataoka et al., 2022; Laskar et al., 2022; Lin et al., 2022; Zhang et al., 2022). Specifically, Laskar et al. (2022) used the empirical NRLMSIS2.0 and showed that the neutral mass density increased by only ~25% on 3–4 February 2022. Kataoka et al. (2022) used the real-time GAIA simulation and showed the mass density enhancement of up to 50% at 200 km. Dang et al. (2022) used TIE-GCM simulations and showed ~20% increase in atmospheric density at 210 km altitude on 3 February as compared to the quiet condition. These results suggest that the storm-time density increase might be smaller than the 50% stated by SpaceX. Previous results have demonstrated that the storm-time density enhancement may be not accurately estimated by a single model due to the discrepancies between different model simulations. In addition, so far, the detailed comparisons of different coupled (magnetosphere) ionosphere and thermosphere models and whole atmosphere models have not been conducted in the simulated upper and lower thermospheric changes.

In this study, we provide a detailed investigation of the upper and lower thermospheric changes due to the storm event. For this purpose, we analyze thermosphere mass density from the POD of the Swarm-A satellite (~438 km altitude) and the accelerometer of the Gravity Recovery and Climate Experiment Follow-On (GRACE-FO, ~505 km altitude). In addition, the observations are compared to densities from the following six empirical and physics-based models: DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model (CTIPe), Multiscale Atmosphere-Geospace Environment (MAGE), Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics (WAM-IPE), and SD-WACCMX. The main aims of the current work are: (a) based on both Swarm-A and GRACE-FO observations, qualitatively characterize the thermospheric disturbance structure and quantitatively estimate the upper thermosphere mass density changes due to this storm, (b) assess the capability of the empirical and physics-based models to simulate this storm's effect in both amplitude and pattern of the thermospheric disturbance, and (c) analyze storm-time density changes in the lower thermosphere (~200 and 300 km), where the Starlink satellites orbited after the launch.

## 2. Data and Models

### 2.1. Data

Our data set consists of GPS-derived thermosphere mass density from Swarm-A and accelerometer-derived mass density from GRACE-FO. During the storm, the Swarm-A orbital altitude varied between 426.8 and 456.3 km (the averaged altitude is 438 km) with 87.35° inclination. The LT sampling from the ascending and the descending orbits at the equator were at ~09:00 and 21:00 LT. The temporal resolution of the recovered aerodynamic accelerations is about 15–20 min in the GPS-derived accelerations, equivalent to 7,000–9,000 km along the orbit (van den IJssel et al., 2020).

The orbital altitude of the GRACE-FO varied between 483.8 and 532.7 km, and the average orbit height was about 505 km with an inclination of 89°. The LT samplings of the ascending and descending orbits were at ~06:00 and 18:00 LT. The time resolution is 10 s. The difference between the accelerometer data and the GPS-derived data are that the latter is a smoothed representation of the true accelerations, and therefore cannot be used to identify localized density variations and traveling atmosphere disturbances (Astafyeva et al., 2017; van den IJssel & Visser, 2007; van den IJssel et al., 2020). These two data sets can be downloaded from the websites https://swarm-diss.eo.esa.int and http://thermosphere.tudelft.nl.

We did not include data from Swarm-B because it flew in the same LT sector as Swarm-A and at roughly the same altitude as GRACE-FO. Swarm-C has the same orbit as Swarm-A, and it is currently the only Swarm satellite providing accelerometer-based densities. However, the Swarm accelerometers are suffering from many issues that affect the quality of their observations. Siemes et al. (2016) have shown that the Swarm accelerometers are perturbed by a number of anomalies, including spikes, sudden changes in the accelerometer bias, and large temperature-induced bias variations. This significantly affects their usefulness for density retrieval. In addition, accurate accelerometer data are not available in the directions perpendicular to the flight direction (van den IJssel et al., 2020).

In order to study the thermospheric response at lower altitudes, we used the GOLD level-2 disk temperatures (GOLD_Tndisk) data set. This remote sensing data set is available around the dayside over about ±70° latitude and 30°E–120°W longitude (~1/4 globe) between about 06:00 and 23:00 universal time (UT) each day. It can be
<table><caption>**Table 1**
The Summary of Models High Latitude Drivers and Resolutions</caption><thead><tr><td>Model</td><td>Setting up parameters (solar | geomagnetic)</td><td>Resolution (longitude, latitude, height, and time step)</td></tr></thead><tbody><tr><td>MSIS 2.0</td><td>F10.7 | ap (array of 7 values)</td><td>2.5°, 2.5°, 10 km, and 10 min</td></tr><tr><td>DTM2013</td><td>F30 | Kp</td><td>1.8°, 0.9°, 10 km, and 15 min</td></tr><tr><td>MAGE</td><td>Solomon and Qian bands | 1-min OMNI data</td><td>1.25°, 1.25°, 0.25 scale height, and 5 min</td></tr><tr><td>CTIPe</td><td>Solomon and Qian bands | Weimer and TIROS</td><td>18°, 2°, 1 scale height, and 1 min</td></tr><tr><td>WAM-IPE</td><td>F10.7 | Weimer, TIROS and OVATION</td><td>4°, 2°, 10 km, and 10 min</td></tr><tr><td>SD-WACCMX</td><td>Solomon and Qian bands | Heelis</td><td>2.5°, 1.9°, 0.25 scale height, and 1 hr</td></tr></tbody></table>

Note. Note that the resolution here means the used model output grid and rate in this study.

downloaded from the GOLD Science Data Center (http://gold.cs.ucf.edu). The GOLD_Tndisk is the integrated temperature accounting for weighting by the emissions at different heights, with the maximum weight at a height of approximately 160 km (Eastes et al., 2020). Therefore, it can be used to characterize the lower thermosphere variation (Laskar et al., 2022). To make the comparison between the observations and simulation results, the synthetic Tndisk from both the empirical and physics-based models are also calculated in the current work. In brief, the modeled Tndisk can be calculated as the integrated temperature accounting for the weighting by the emissions at different heights, and the maximum weight is at a height of approximately 160 km. The contribution function (weight) can be represented by the daytime N2 LBH band emissions, which can be found in more detail in Laskar et al. (2021a) and He et al. (2021).

## 2.2. Models

Two empirical and four physics-based models are used to evaluate the capability to reproduce the storm-time effects. The model resolutions and drivers used in this study are listed in Table 1. Note that all models were output on a regular grid in order to allow the same interpolation schemes and then to consistently make a comparison across the models.

1. The Naval Research Laboratory MSIS version 2.0 (MSIS2.0) is an empirical atmosphere model that extends from the ground to the exosphere. It describes the temperature, density, and composition of the Earth's thermosphere. The model inputs are the array of geomagnetic activity Ap index (i.e., daily and current time, 3, 6, and 9 hr before current time, and average values from 12 to 33 and 36 to 57 hr prior to current time), 81-day average and daily F10.7 index. MSIS2.0 is a reformulated upgrade of the previous version, which couples the neutral species densities to the entire vertical column using the effective mass of species in a geopotential height coordinate. The cooler stratospheric and mesospheric temperatures likely have an influence on the thermosphere neutral densities (Emmert et al., 2020). For model outputs, the latitude resolution is 2.5°. The longitude resolution is 2.5°. In the vertical direction, the step is 10 km. The output model parameters are calculated with a 10-min time step in this study.

2. The DTM 2013 (DTM2013) is a semi-empirical model reproducing the thermospheric temperature, density and composition in the altitude range from 120 to 1,500 km. DTM2013 assimilated high-resolution accelerometer-derived densities from CHAMP, GRACE, and low-altitude GOCE (Bruinsma, 2015). This model is run with the 30-cm solar radio flux (F30 index) and geomagnetic Kp index. In this study, the outputs of the model parameters have spatial resolutions of 0.9° in latitude, 1.8° in longitude, and 10 km in height, and temporal resolution of the output is 15 min.

3. The MAGE couples together the Grid Agnostic magnetohydrodynamic (MHD) with Extended Research Applications (GAMERA) model of the magnetosphere, the Rice Convection Model of the ring current, and the TIE-GCM. In MAGE, the high-latitude inputs to TIE-GCM are provided via a coupling layer (RE-developed Magnetosphere-Ionosphere Coupler/Solver). This provides a more dynamic temporal variation of high-latitude convection and auroral precipitation compared to using the empirical Heelis or Weimer convection models and Ray and Roble precipitation model (Lin et al., 2021; Pham et al., 2022). Solar inputs from F10.7 radio solar flux measurements are used as a proxy for XUV/EUV/FUV solar flux as described by Solomon and Qian (2005). In this study, the thermosphere part, that is, TIE-GCM, has a uniform resolution
of 1.25° in horizontal direction, and a vertical pressure grid resolution of 0.25 scale height. The temporal resolution of the output is 5 min.

4. The CTIPe consists of a global thermosphere, a high-latitude ionosphere, a mid- and low-latitude ionosphere/plasmasphere, and an electrodynamical calculation of the global dynamo electric field. The thermosphere code was originally developed by Fuller-Rowell (1981), which simulated the time-dependent global structure of the neutral thermosphere by numerically solving the non-linear primitive equations of momentum, energy, and continuity on a three-dimensional spherical polar grid rotating with the Earth. The model is driven by the high-latitude empirical Weimer model, and by the magnitude of the interplanetary magnetic field (IMF) in the y-z plane, together with the velocity and density of the solar wind. The auroral precipitation is specified by an empirical model based on TIROS-NOAA satellite observations. Solar heating, ionization and dissociation rates, and their variations with solar activity are specified by the Solomon and Qian (2005) solar EUV energy deposition scheme for upper atmospheric general circulation models. The resolutions of the model outputs are 2° in latitude, 18° in longitude. The temporal resolution of the output is 1 min.

Due to the limitation of the upper boundary of CTIPe being ~400 km for the time period of our study, that is, 15 levels in a vertical direction in a logarithm of pressure from a lower boundary of 1 Pa (~80 km altitude), we extrapolate the upper thermosphere mass density to make a comparison with the Swarm-A and GRACE-FO observations. The extrapolation is calculated using:

$$ \text{den}_i = \text{den}_0 * \exp\left(\frac{-(h_i - h_0) * (M_0 * g)}{R * T_0}\right) \quad (1) $$

where, denᵢ is the extrapolated density at the height of interest hᵢ, den₀, h₀, T₀, and M₀ are the density, height, temperature, and mean molecular mass at the upper boundary. g is the gravitational acceleration and R is the gas constant.

5. The coupled WAM-IPE includes two parts, which are focused on weather forecast timescales (Fuller-Rowell et al., 2008). The WAM is an extension of the Global Forecast System to model the upper atmospheric physics and chemistry. The IPE model (Maruyama et al., 2016) provides the plasma component of the atmosphere. The flux-tube solver is based on the Field Line Interhemispheric Plasma model and flux-tube structures are defined by the International Geomagnetic Reference Field. The thermospheric parameters calculated by WAM are fed into IPE for calculating the ionospheric responses. The space weather drivers of WAM-IPE include F10.7, Kp (estimated/forecast), 41-day backwards averaged F10.7 and 24 hr averaged Kp. This model uses the Weimer empirical ion convection model. The IMF and solar wind are used to derive the OVATION hemispheric power. This hemispheric power is then used to drive the TIROS model of the precipitation. The resolutions of the model outputs have 2° in latitude, 4° in longitude, 10 km in height. The temporal resolution of the output is 10 min.

6. The WACCM-X extends the lower atmospheric component of the Community Earth System Model into the thermosphere, with a model top boundary between 500 and 700 km (H.-L. Liu et al., 2018). Compared to the TIE-GCM, it has the advantage of self-consistently resolving lower atmospheric processes and therefore enables a more realistic simulation of upper atmospheric variability due to lower atmospheric forcing. In addition, the meteorology of the troposphere and stratosphere is constrained using the specified dynamics approach (Smith et al., 2017), namely, the “SD-WACCMX.” In this study, we use the empirical Heavis model to represent the magnetospheric energy input. The solar extreme-ultraviolet (EUV) variability and energy deposition scheme, including photoelectron effects, are parameterized as described in Solomon and Qian (2005). The model resolution is 2.5° in longitudinal direction, 1.9° in latitudinal direction and a vertical pressure grid resolution of 0.25 scale height. The temporal resolution of the output is 1 hr.

# 3. Results

## 3.1. Solar and Geophysical Conditions of the Geomagnetic Storm

Figure 1 shows the space environment conditions for 1–6 February 2022. Two stronger, but yet still moderate, geomagnetic disturbances occurred on 3 and 4 February under consideration (yellow shaded rectangles in Figure 1), for which we will refer to as the first and second storm. The solar wind speed increased at 00:00 UT on 3 February from 400 to 580 km/s. The IMF Bz turned southward at 01:00 UT on 3 February, and reached a minimum value of ~-18 nT at 09:00 UT. Correspondingly, the SYM-H index decreased, and a negative value of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Interplanetary magnetic field (IMF) and solar wind parameters and geomagnetic activity indices on 1–6 February 2022. From the top to the bottom are the solar wind velocity (km/s, 5 min) and proton number density (n/cm³, 5 min), IMF Bz components (nT, 5 min), geomagnetic symmetric (SYM) disturbances for H component (nT, 5 min) and Dst (nT, 1 hr) and Kp (3 hr).</figcaption></figure>

−79 nT was reached at 11:00 UT. We can also see that the minimum Dst is −66 nT, and the maximum Kp is 5 at this moment. On 4 February, the solar wind remained at around 500 km/s. The IMF Bz remained southward during the whole day, and minimum values of ~−10 nT can be seen at 09:00 and 17:00 UT. The SYM-H index reached a minimum of −61 nT at 10:00 UT and −70 nT around 21:00 UT. There exists a comparable Kp value on 3–4 February. Note that the SYM-H varied around ~−27 nT on 5–6 February, which means that there was a long-term recovery phase that lasted more than 2 days.

As quiet-time reference values, for both observations and modeling, we took the average of 1 and 2 February 2022. Both of these days are characterized by a minor geomagnetic disturbance (Figure 1), however, they did not have impact on the thermosphere behavior.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The comparison of the simulated thermosphere mass density for the DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics and SD-WACCMX with Swarm-A (average orbit height: ~438 km) and Gravity Recovery and Climate Experiment Follow-On (average orbit height: ~505 km) during the time period from 1 to 6 February.</figcaption></figure>

## 3.2. Thermosphere Mass Density Comparisons

Figure 2 shows variations of the thermosphere mass density from 1 to 6 February as observed by Swarm-A and GRACE-FO satellites. We also calculate neutral density from the model outputs along the satellite trajectories. Both satellites (blue lines) show two significant increases in the thermospheric mass density that occurred on 3–4 February 2022, which corresponds to the minimum of the SYM-H index. The maximum neutral mass density
during the second storm is higher than that during the first storm. For the Swarm-A, the density increased from ~0.8 × 10⁻¹² on 1 February to ~1.5 × 10⁻¹² kg/m³ around 11:00 UT during the first storm (i.e., 3 February) and to ~1.6 × 10⁻¹² kg/m³ around 21:00 UT during the second storm (i.e., 4 February). For the GRACE-FO, the density increased from ~0.3 × 10⁻¹² to ~0.5 × 10⁻¹² kg/m³ on 3 February and to ~0.6 × 10⁻¹² kg/m³ on 4 February. We also notice a small density increase on 2 February between 00:00 and 16:00 UT that corresponds to the very minor disturbance during this time. Similar to both observations, the two thermosphere density enhancements at around 11:00 on 3 February and around 21:00 UT on 4 February are captured by all six models (red lines in Figure 2). The enhancement on 2 February is also seen in the simulations. However, we also notice some discrepancies between the model results and observations. For instance, as compared to Swarm-A, DTM2013 overestimates the density during the recovery phase on 5–6 February, and during the minor geomagnetic disturbance on 2 February. While, MSIS2.0 underestimates the thermospheric density in response to the main phases on 3–4 February. MAGE overestimates by ~50% the value and by ~12 hr the duration of the thermospheric storm on 4 February. Finally, the largest discrepancy is observed for WAM-IPE, which overestimates the density for both quiet and disturbed conditions. Similar discrepancies between models and observations can also be confirmed by the observed-to-simulation density ratio in Figure S1 in Supporting Information S1.

However, it should be noted that the absolute values calculated from satellite observations can be affected by the methods of density estimation. The origin of systematic offsets between thermosphere neutral density data sets and empirical models stems from the selection of input data combined with uncertainties in satellite aerodynamic modeling. An additional factor, affecting comparisons of models derived from older data sets compared with newer data, stems from long-term change effects. This is discussed in more detail in Doornbos (2012). On the other hand, physics-based models contain parameters for sub-grid processes with an associated uncertainty, such as the eddy-diffusion coefficient, that are sometimes tuned toward specific data sets and/or empirical models. The discrepancies presented above are likely the result of the varying approaches used by the model developers.

To quantitatively describe the storm-time thermosphere mass density changes, we estimate the relative variations of the orbit-averaged density with respect to the quiet days, $\Delta\rho/\rho$, for both observations and models (Figure 3). Note that the orbit-averaged density is calculated over one orbital period for each time step. As a reference, the quiet and storm orbit-averaged densities for both observations and models are also provided in Figure S2 in Supporting Information S1. From Figure 3, for Swarm-A at an average altitude of ~438 km and 09/21 LT sector, the relative neutral mass density reached 90% at 11:00 UT on 3 February (the gray line in Figure 3), and ~110% at 21:00 UT on 4 February. GRACE-FO measurements at the average altitude of 505 km and in 06/18 LT sector showed ~90% and ~120% increases in neutral density for the first and second storms, respectively. For the simulation results (other colored lines in Figure 3) relative to Swarm-A on the first storm days (3 February), DTM2013 (the magenta line), CTIPe (the blue line), MAGE (the cyan line), and WAM-IPE (the red line) are more consistent with the observations. However, MSIS2.0 (the green line) and SD-WACCMX (the black line) underestimate the thermospheric density by ~50%. For the second disturbance on 4 February, CTIPe is much closer to the observations than other models. Meanwhile, we observe ~30% underestimation for DTM2013 and WAM-IPE, ~50%–70% underestimation for MSIS2.0 and SD-WACCMX models, and ~10% overestimation for MAGE. During the recovery phase on 5 February, the density values of DTM2013, MSIS2.0, WAM-IPE, and SD-WACCMX are close to the observations, while CTIPe and the MAGE show ~30% overestimation. On 6 February, the ~20% underestimation for both the WAM-IPE and the MAGE can be found as compared to the observations. The results are quite similar for the GRACE-FO observations (the bottom of Figure 3). Overall, we conclude that both empirical and theoretical models captured the enhanced storm-time density changes; however, different models exhibit positive and/or negative biases at different phases of this storm event.

To investigate the thermosphere mass density behavior due to the 3–4 February magnetic storm in more detail, we analyze the latitude-universal time (UT) distribution of the density enhancement (Figures 4–7). Note that we compare results measured and simulated at different local times, that is, ~09/21 LT for Swarm-A (Figures 4 and 5) and ~06/18 LT for GRACE-FO (Figures 6 and 7). An evident pattern is the hemispheric asymmetry, that is, a larger density enhancement existed in the northern hemisphere relative to the southern hemisphere. Specifically, during the main phase of the first storm on 3 February, the density from Swarm-A (the top panels of Figure 4) starts to increase at high latitudes of the northern hemisphere from 06:00 UT, and reaches a maximum of 150% at 11:00 UT. The enhancement propagates into the low latitudes around 13:00 UT, and lasts until ~21:00 UT. With the second storm development, the neutral density starts to increase at high and mid-latitudes in the northern hemisphere at 08:00 UT on 4 February, and at ~12:00 UT in the southern hemisphere. The maximum
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The relative changes comparison of the orbit averaged thermosphere mass density from simulations with the observed one from Swarm-A (average orbit height: ~438 km) (top panel) and Gravity Recovery and Climate Experiment Follow-On (average orbit height: ~505 km) (bottom panel) during the days from 3 to 6 February. The gray line indicates the observations. The magenta, green, blue, cyan, red, and black lines indicate simulation from DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics, and SD-WACCMX, respectively.</figcaption></figure>

enhancement of 150% is reached between 16:00 and 20:00 UT in both the northern and the southern hemispheres. From ~20:00 UT, the density enhancement drops to 30%–60% in both hemispheres, and it lasts during the next day at all latitudes. On 6 February, we observe 10%–30% enhancement at high latitudes, while at low and middle latitudes the neutral mass density returned to undisturbed levels. The simulation results in the 09:00 LT sector show quite similar patterns (Figure 4). The 3–4 February storm effects are also seen in the results of all models. However, MSIS2.0 and SD-WACCMX show an enhancement of the neutral mass density at high latitudes of the northern hemisphere, and only a weak disturbance propagates to the middle and low latitudes. Additionally, in MSIS2.0, the amplitude of the high-latitude neutral density disturbance is underestimated. DTM2013 reproduces the effects of both storms, but the amplitude of the thermospheric disturbance, especially at middle and low latitudes, is underestimated. In terms of the neutral density disturbance propagation from high latitudes to low latitudes, CTIPe, WAM-IPE, and MAGE show better agreement with observations at the altitude of Swarm-A and in the 09:00 LT sector. The recovery phase (5–6 February) reproduced by DTM2013 and WAM-IPE agree the most with the observations. CTIPe and MAGE show an obvious overestimation of the amplitude at all latitudes, especially on 5 February. In the 21:00 LT sector (Figure 5), the thermospheric storm effects are stronger and longer lasting than in the morning sector. They are underestimated by the empirical MSIS2.0 model, and by WAM-IPE and SD-WACCMX models. The agreement with observations is better for CTIPe and MAGE models; however, these models overestimate the amplitude of the thermosphere disturbance.

The GRACE-FO observations at 483.8 and 532.7 km and in 06/18 LT depict smaller-scale thermospheric disturbances structure due to the storm (Figures 6 and 7), because the resolution of the accelerometer-derived data set is higher than
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The latitude-universal time distribution of the relative changes comparison of the orbit thermosphere mass density from observed one around 09 local time provided by the Swarm-A with the simulated results from different models during the days of 3–6 February. From the top to bottom panels indicate Swarm-A, DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics, and SD-WACCMX, respectively.</figcaption></figure>

the GNSS-derived one. Similar to Swarm-A observations, we observe up to 150% density enhancement from the high to low latitudes during the storms on 3–4 February (the top panel of Figure 6). During both storms, the enhancements reach low latitudes due to large-scale TAD. The density enhancements are first seen at high latitudes in both the northern and southern hemisphere and then propagate to low latitude and equatorial regions after ~4 hr. DTM2013 predicts a smoothed density enhancement extending to the low latitudes from the high latitudes. The smooth enhancement in MSIS2.0 is mainly confined to the high latitudes, hardly extending to the low latitudes. The numerical models
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Figure 4, but for the 21 local time.</figcaption></figure>

can reproduce small-scale structures of the thermospheric response. CTIPe and MAGE agree best with observed densities. WAM-IPE underestimates the density enhancement in all latitudes on 3–4 February. SD-WACCMX shows a stronger enhancement at high latitudes of the northern hemisphere, and only a weak disturbance descends to the middle and low latitudes. A similar pattern can also be seen in the 18:00 LT sector (Figure 7). However, a smaller neutral density enhancement is seen around the low-latitude and equator regions on both 3 and 4 February. It suggests that thermospheric storm effects are weaker at 18:00 LT with respect to the 06:00 LT observations.

To study the response in the lower thermosphere, we analyze neutral temperatures (Tn) from GOLD disk measurements. Figure 8a presents the relative Tn changes (i.e., GOLD_Tndisk) at fixed 15:00 UT with respect to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The latitude-universal time distribution of the relative changes comparison of the orbit thermosphere mass density from the early morning (06 local time) observed one provided by the Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) with the simulated results from different models during the days of 3–6 February. From the top to bottom panels indicate GRACE-FO, DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics, and SD-WACCMX, respectively.</figcaption></figure>

the quiet day. Figure 8b shows relative changes in the synthetic temperature from all empirical and theoretical simulations during the days of 2–6 February. From the GOLD observations on 3 February, it is evident that a larger temperature enhancement occurred at higher latitudes in the northern hemisphere (~30%–40%), and there was no obvious enhancement at lower latitudes, which was around ~10%–15% in middle latitudes in the southern hemisphere in the American sector on 3 February. An overall larger temperature enhancement from the higher to lower latitudes for both the northern and southern hemispheres can be seen on 4 February. The temperature
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 6, but for the 18 local time.</figcaption></figure>

enhancement then decreased on 5–6 February, but it remained around 10% in the middle and low latitudes in the American sector. The observed GOLD_Tndisk results demonstrate that the neutral temperature changed in the lower altitudes due to this storm event, which is consistent with the observed mass density results. The simulations captured the disturbed lower thermosphere temperature during the storm period. However, the synthetic temperatures from the simulation result show clear deviations from the observations in both magnitudes and spatial patterns. The temperature enhancement of DTM2013 is underestimated (~5%) during the entire time period, and it only occurred at higher latitudes. A similar phenomenon can also be seen in MSIS2.0, that is, the underestimation (~5%) is confined in the higher latitudes. For CTIPe, it can generally represent the temperature changes in
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) The observed Global-scale Observations of the Limb and Disk neutral temperature enhancement at 15:00 universal time and (b) The synthetic neutral temperature enhancement from simulations during the days of 3–6 February. From the top to bottom panels of (b) indicate DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics, and SD-WACCMX, respectively.</figcaption></figure>
higher and middle latitudes around the American sector during the selected period. However, it slightly underestimates the enhancement relative to the observations. For example, in CTIPe, the enhancement in the northern high latitudes is around ~25%, and it is ~10% in the southern middle latitudes on 3 February. The underestimation of the temperature enhancement persists on 4–6 February, especially for the mid- and lower-latitudes in the American sector. For MAGE, it better represents the temperature changes during the selected period. However, an overall slight overestimation can be found during the storm's main phase (e.g., 3–4 February). For example, a larger temperature enhancement occurred in the northern higher latitudes on 4 February. WAM-IPE shows a realistic temperature enhancement pattern on 3 February, however, it underestimates the temperature changes at the selected times. SD-WACCMX reveals a smaller temperature enhancement value relative to the results observed by GOLD_Tndisk at all latitudes from 2 to 6 February.

We notice that before reentering the atmosphere, SpaceX initially launched and a short time later deployed 49 Starlink satellites at an orbit altitude between 210 and 320 km. The thermosphere mass density enhancement around 200 and 300 km should be investigated in detail. Figure 9 further gives the simulated zonal mean value of the thermosphere mass density enhancement at a height of 200 km during the days of 3–6 February. The reference day is again chosen as 1 and 2 February. With regard to the value of global mean density enhancement, DTM2013 predicts that the maximum density enhancement can be up to ~30%–35% on both 3 and 4 February and ~15% on 5–6 February. For MSIS2.0, it shows ~10%–15% enhancement during the two storm days, while no obvious enhancement is shown in the recovery phase. For CTIPe, enhancements of ~20% occurred on 3 February, ~30% on 4 February, and ~10%–20% on 5–6 February. For MAGE, it shows ~30% density enhancement occurred on 3 February in the middle and low latitudes, with enhancements of ~40% on 4 February, and ~10%–20% on 5–6 February. For WAM-IPE, it shows ~10%–15% density enhancement occurred on 3–4 February, and almost no obvious enhancement is shown on 5–6 February. For the SD-WACCMX, only ~10% density enhancement is seen on 4 February. It seems that there exists a reversed hemisphere asymmetry at 200 km compared to the higher altitudes (Figures 4–7), that is, the larger density enhancements are visible in the Southern Hemisphere at 200 km. The detailed reason for this will be given in discussion section.

The simulated zonal mean value of the thermosphere mass density enhancement at the altitude of 300 km is also shown in Figure 10. Similar to the results in the upper thermosphere (e.g., Figure 4), the majority of simulations show that two obvious density enhancements in the lower thermosphere occurred on 3–4 February, and a larger enhancement existed on 4 February. During the two storm days, the global averaged maximum density enhancement of DTM2013 can reach ~50% and ~45% on 3 and 4 February, respectively. They are ~20% and ~25% for MSIS2.0, ~41% and ~51% for CTIPe, ~57% and ~75% for MAGE, ~40% and ~40% for WAM-IPE, and ~15% and ~23% for SD-WACCMX.

Based on the thermosphere mass density comparison, during the two storms main phases, better density changes can be found in the theoretical models CTIPe and MAGE (seen in Figure 2). Meanwhile, according to the GOLD temperature comparison, better temperature changes can be found in the two theoretical models for the CTIPe and the MAGE, with a slight underestimation and overestimation for CTIPe and the MAGE, respectively (seen in Figure 8). Thus, we can estimate that the global averaged maximum density enhancement based on the average value of both the CTIPe and the MAGE models, which can be up to ~35% at 200 and ~60% at 300 km during the two storm days.

## 4. Discussion

The loss of 38 Starlink satellites in early February 2022 due to a geomagnetic storm is the perfect demonstration of the importance of reliable space weather forecasts. While space agencies and space weather prediction centers issued a warning of the possibility of a minor storm, its consequences on the orbits of the very low perigee Starlink satellites were apparently severely underestimated. The February 2022 storm was of a moderate intensity (the minimum SYM-H excursion of only ~79 nT, the minimum Dst of ~66 nT, and the maximum Kp of 5), and comprised two sub-events. Our results show that the thermospheric neutral mass density increased by 120%–150% at the altitudes of 400–500 km, and by 35%–60% at the orbital altitude of the Starlink satellites (210–320 km). Such effects are quite moderate as compared to previous observations for more intensive storms that can cause up to 900%–1100% increases in the thermospheric density (e.g., Astafyeva et al., 2017, 2020; Lei et al., 2011; Liu & Lühr, 2005) and, therefore, a significant increase in the atmospheric drag.

Both empirical and theoretical models captured the upper and lower thermospheric response to this storm; however, we found up to 70% discrepancy in the results of different models during different phases of this
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The simulated zonal mean value of the thermosphere mass density enhancement at a height of 200 km during the days of 03–06 February. From the top to bottom panels indicate DTM2013, MSIS2.0, Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, Whole Atmosphere Model-Ionosphere Plasmasphere Electrodynamics, and SD-WACCMX, respectively.</figcaption></figure>

storm. For the empirical models, both MSIS2.0 and DTM2013 use the same Bates temperature profile and a physical constraint of diffusive equilibrium for the vertical structure of the thermosphere (Emmert, 2015). The temperature changes are very similar, but the mass density enhancement variations show large differences in both upper and lower thermosphere. Since, in addition to the temperature, the density at a fixed height is dependent on the composition profile, the differences in the main thermosphere constituents and their response to geomagnetic
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Figure 9, but for 300 km.</figcaption></figure>

activity could be one possible explanation of the different response of the models. Another possible explanation is likely related to the assimilation of more recent total density data for the DTM2013 model (Bruinsma, 2015), especially the GOCE data at low altitudes (~224 km).

Generally, the Joule heating is the predominant driver of the increased temperatures in the high latitudes during the storm-time period. The time-dependent evolution of high-latitude Joule heating is adopted in the physics-based models, which is most calculated from the empirical Heelis and Weimer models (Heelis et al., 1982;
Weimer, 2005). However, for the empirical models, they are usually driven by the scalar geomagnetic indices (e.g., 3 hr Kp and Ap indices) to describe the varying heating distributions (Fedrizzi et al., 2012). Thus, both empirical models in this study poorly capture the inherent time dependence of the global atmospheric structure during geomagnetic activity. In addition, we note that the performance of the models, that is, the CTIPe, WAM-IPE and the MAGE, is better than the SD-WACCMX one, at least during the storm's main phase (seen in Figure 2). T.-W. Fang et al. (2022) indicated that WAM-IPE densities at perigee (200–250 km) agree well with the onboard calculated neutral density proxy provided by the Starlink satellites. According to T.-W. Fang et al. (2022), the neutral density proxy is obtained by dividing the drag acceleration by drag coefficient along the orbit for each satellite, which should be linearly proportional to the neutral density. Lin et al. (2022) compared the high latitude convection patterns between MAGE and standalone TIE-GCM during this storm, and showed that the TIE-GCM convection is generally weaker than DMSP measurements, while MAGE simulation results are more consistent with DMSP data. Besides, Zhu et al. (2022) also indicated the importance of accurate soft electron precipitation specifications in the whole auroral zone to improve the F-region neutral mass density estimations by the Global Ionosphere Thermosphere Model. The accurate specification of the energy input from the magnetosphere is crucial. Capturing the variabilities of electric field and electron precipitation through either coupled magnetosphere-ionosphere-thermosphere modeling or data-constrained modeling (e.g., AMIE) is therefore important to better simulate the thermosphere storm's response.

Figure 11 shows that the simulated integrated high latitude Joule heating rate poleward of 45° in both the northern hemisphere and southern hemisphere from three different physical models, that is, CTIPe, MAGE, and SD-WACCMX. Note that the WAM-IPE does not presently support output of its Joule heating calculations. For the northern hemisphere, CTIPe and MAGE show a larger Joule heating than other models during the whole period of this storm. Both values are close to each other on 3 February. However, a larger Joule heating occurred in the MAGE than that in the CTIPe around 06:00–14:00 UT on 4 February in both hemispheres. This is consistent with the bigger thermosphere mass density enhancement shown in the MAGE on 4 February. The SD-WACCMX shows a smallest Joule heating value, that is, the smallest high latitude energy input, and therefore shows a weaker thermosphere mass density enhancement during two main phases of this storm event. For the southern hemisphere, the same conclusion can also be conducted, while there exists a larger Joule heating than that in the northern hemisphere.

Besides the absolute value differences, we notice, for both observations and simulations, stronger density enhancements at 06/21 LT than other LT sectors. Such dependence can be explained by the LT asymmetry of the high-latitude input. The energy injection is the largest along the auroral oval, that is, it is shifted toward the night-side hemisphere. Therefore, the expansion of the upper atmosphere disturbance from the high latitude toward middle latitudes is essentially restricted to the midnight/early morning hours (Fuller-Rowell et al., 1996; Prölss, 1976, 1980). This can be confirmed by Figure 12 showing the simulated height integrated high latitude Joule heating rate in the northern hemisphere from CTIPe, MAGE, SD-WACCMX at 10:00 UT on 3 and 4 February. The selected UT corresponds to the minimum value of the IMF Bz (Figure 1). The simulated height integrated Joule heating changes during the whole days of 3 and 4 February are also given in Figures S4–S9 in Supporting Information S1. Even though the simulated absolute values differ, a larger Joule heating high-latitude energy input for all chosen models around the midnight/early morning hours.

Compared with other models, the SD-WACCMX can provide a more realistic representation of variability coming from the lower atmosphere. However, during the main phase of this storm, the SD-WACCMX has a worse performance than other physics-based models utilized in the study. The reason is that the density variation at LEO is mainly driven by geomagnetic and solar forcing. However, the values of the SD-WACCMX show much closer to that of both the Swarm and the GRACE-FO during the recovery phase (seen in Figure 3). This could be contributed to the storm-induced chemistry and dynamics changes of nitric oxide density, which will be investigated in a further study. This also suggests that a realistic lower atmosphere is also important for specifying thermosphere mass density even during the storm period. Even though this effect is not investigated in detail in this paper, according to Yue et al. (2022), at LEO altitudes, during solar minimum and geomagnetic quiet days, the impact from the lower atmosphere is much smaller compared to solar and geomagnetic effects, but is not negligible (5%–10% vs. 20%). Maute et al. (2022) showed that more realistic lower atmosphere variations improve especially the northern hemisphere neutral density by up to 15% compared to using a climatological representation of the lower atmosphere forcing. In terms of the thermosphere mass density enhancement during storm's main phase, the performance of SD-WACCMX is not as good as the other models most likely due to deficiencies in the
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The integrated high latitude Joule heating rate poleward of 45° in both the northern hemisphere (the top panel) and southern hemisphere (the bottom panel) simulated by the Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model (red), Multiscale Atmosphere-Geospace Environment (green), and SD-WACCMX (dark).</figcaption></figure>

Heelis empirical model. This again highlights the importance of accurate specification of the high latitude energy input for coupled ionosphere and thermosphere system. This also suggests that the Heelis model is less accurate to specify the high latitudes energy input at least during this moderate storm period.

The thermospheric behavior during the 3–4 February 2022 storm was hemispherically asymmetric around satel-
lite heights as seen in Figures 4–7, that is, the storm-time density enhancement in the northern hemisphere
(winter) is stronger than in the southern hemisphere (summer) in the upper thermosphere (i.e., >400 km). While
in the lower thermosphere (e.g., 200 km), as seen in Figure 9, the asymmetry is less pronounced or even reversed
on 4 February. Physically, neutral mass density changes at a fixed altitude are generally influenced by scale height
changes that are partly caused by storm-time neutral temperature enhancements due to the heating and partly by
changes in the mean molecular weight due to winds (Lei et al., 2010). To investigate the detailed reason, Figure 13
gives the mean molecular mass and temperature enhancement from MAGE during the storm events in the upper
(400 km) and lower thermosphere (200 km), respectively. A larger temperature enhancement still occurred in
the northern hemisphere at 200 and 400 km during this storm period. For the upper thermosphere, there exists a
comparable mean molecular mass enhancement in both the northern and southern hemispheres. Thus, there is a
stronger storm-time density enhancement at 400 km in the northern hemisphere (winter) compared to the south-
ern hemisphere (summer). However, for the lower thermosphere, a stronger mean molecular weight enhancement
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The height-integrated Joule heating in Northern Hemisphere simulated by the Coupled Thermosphere Ionosphere Plasmasphere Electrodynamics Model, Multiscale Atmosphere-Geospace Environment, and SD-WACCMX at 10:00 universal time on 03 February (the left panels) and on 04 February (the right panels). The patterns are plotted in latitude versus local time coordinates. The white dashed line represents the 45° latitude boundary in Northern Hemisphere.</figcaption></figure>

occurred in the northern hemisphere especially on 4 February. Those increases in mean molecular weight offset the effects of neutral temperature increases. Thus, we found a disappeared asymmetry or even a reversed asymmetry at 200 km. As a reference, Figure S3 in Supporting Information S1 also shows the background values of the thermosphere mass density, neutral temperature, and mean molecular mass for calculation of Figure 13.

The GOLD neutral temperature observations provide insight into the storm effects on the lower thermosphere. This type of observation shows that even during a moderate geomagnetic storm the lower thermospheric conditions can be perturbed significantly (e.g., 30%–40% seen in Figure 8). Previous studies have demonstrated that auroral heating effects on temperature, which include both particle precipitation and Joule heating, can penetrate down to the altitudes of ~105 km at high latitudes during storms (X. Fang et al., 2008; Xu et al., 2013). Joule heating values (Figure 11) suggest that a larger Joule heating can cause more obvious neutral temperature
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The temperature and mean molecular mass enhancement from Multiscale Atmosphere-Geospace Environment at 400 (the top panels) and 200 km (the bottom panels) during the time period from 3 to 6 February.</figcaption></figure>

enhancement in the lower thermosphere especially during the storm's main phase. Due to the exponential decrease of neutral density with altitude, the direct influence of Joule heating on temperature is largest in the *F* region. Deng et al. (2011) found that the energy deposited at high altitudes has a stronger influence on the neutral density variation at *F* region altitudes than the energy deposited at low altitudes. Their simulations show that heating below 150 km can cause a large buoyancy acceleration at 300 km. It mainly brings some temporal variation to the vertical velocity. However, the heating at high altitudes sets up a large vertical wind, which is the primary source for a large increase of neutral density at 300 km and higher altitudes.

**5. Conclusions**

In this work, we used observations and simulations of the thermosphere mass density to perform a detailed investigation of the thermospheric conditions during the moderate geomagnetic storm of 3–4 February 2022 that was responsible for the loss of 38 Starlink satellites.

Data from Swarm-A (~438 km) and GRACE-FO (~505 km) showed that the maximum thermosphere mass density enhancement due to this storm was ~110% for Swarm-A observations in the ~09/21 LT sector and ~120% for GRACE-FO observations in 06/18 LT. The simulations suggest that at altitudes of 200–300 km the global averaged density enhanced up to ~35%–60% during this storm event. Although this enhancement is modest when compared to other more intense storms, it was enough to deorbit 38 out of 49 Starlink satellites.

The observed storm-time enhancements extended to middle and low latitudes, and were stronger in the northern hemisphere. The enhancements show a LT difference, that is, a stronger one at 21:00 LT for the Swarm-A, and 06:00 LT for GRACE-FO than in the other LT sectors. Both empirical and theoretical models captured the upper and lower thermospheric response to this storm; however, we found up to 70% discrepancy in the results of different models during different phases of this storm.
This moderate geomagnetic storm also caused a larger (~30%–40%) neutral temperature enhancement in the lower thermosphere as seen in GOLD_Tndisk observations. The simulations show that the magnitude of the temperature enhancement is highly correlated with the Joule heating during the storm's main phase.

This study provides a comprehensive investigation of the whole thermosphere conditions during this moderate storm event. It shows the advantage of the multi-instrumental observations to simultaneously obtain insight into the atmosphere response to one particular storm event. Moreover, this work highlights the importance of both the empirical and theoretical models to better understand the complex storm processes in the global-scale geospace environment, and to satisfy a growing need in the space community to improve thermosphere density predictions. Furthermore, this study shows that even a moderate storm can significantly change the thermosphere conditions in upper and lower altitudes, and thereby cause a risk to the LEO satellites. Currently, with the greatly increasing human space activities, space weather events should be paid much attention to better alleviate the financial loss.

## Data Availability Statement

The Swarm-A and GRACE-FO thermosphere mass density data set can be downloaded from the ESA's Swarm website, https://swarm-diss.eo.esa.int. The GOLD level-2 disk temperatures can be downloaded from the GOLD Science Data Center, https://gold.cs.ucf.edu. The simulation results are available for free download from Open Science Framework (He, 2023).

## Acknowledgments

This work was supported by the Project of Stable Support for Youth Team in Basic Research Field, CAS (YSBR-018), the B-type Strategic Priority Program of the Chinese Academy of Sciences (GrantXDB41000000), the Chinese Meridian Project, the International Partnership Program Of Chinese Academy of Sciences (Grant 183311KYSB20200003), the National Natural Science Foundation of China (42104160), and the China Scholarship Council (CSC). This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #537 Data Assimilation in the Ionosphere and thermosphere. We thank the French National Space Agency (CNES), projects "Swarm" and "RealDetect," for the support. Christian Siemes acknowledges ESA for funding the generation of the Swarm data set through ESA Contract 4000102140/10/NL/JA and the GRACE-FO data set via the Swarm Data, Innovation and Science Cluster (DISC), Sub-contract SW-CO-DTU-GS-129. DL was supported by NASA DRIVE Science Center for Geospace Storms (CGS) under Grants 80NSSC22M0163, 80NSSC21K0008, 80NSSC21K1677, 80NSSC22K1635, and 80NSSC23K1055. This material is also based upon work supported by the National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the U.S. National Science Foundation under Cooperative Agreement No. 1852977.

## References

Aa, E., Zhang, S.-R., Wang, W., Erickson, P. J., Qian, L., Eastes, R., et al. (2022). Pronounced suppression and X-pattern merging of equatorial ionization anomalies after the 2022 Tonga volcano eruption. *Journal of Geophysical Research: Space Physics*, 127(6), e2022JA030527. https://doi.org/10.1029/2022JA030527

Anderson, R. L., Born, G. H., & Forbes, J. M. (2009). Sensitivity of orbit predictions to density variability. *Journal of Spacecraft and Rockets*, 46(6), 1214–1230. https://doi.org/10.2514/1.42138

Astafyeva, E., Bagiya, M. S., Förster, M., & Nishitani, N. (2020). Unprecedented hemispheric asymmetries during a surprise ionospheric storm: A game of drivers. *Journal of Geophysical Research: Space Physics*, 125(3), e2019JA027261. https://doi.org/10.1029/2019JA027261

Astafyeva, E., Zakharenkova, I., Huba, J. D., Doornbos, E., & van den IJssel, J. (2017). Global Ionospheric and thermospheric effects of the June 2015 geomagnetic disturbances: Multi-instrumental observations and modeling. *Journal of Geophysical Research: Space Physics*, 122(11), 11716–11742. https://doi.org/10.1002/2017JA024174

Balan, N., Yamamoto, M., Liu, J. Y., Liu, H., Lühr, H., & Lühr, H. (2011). New aspects of thermospheric and ionospheric storms revealed by CHAMP. *Journal of Geophysical Research*, 116(A7), A07305. https://doi.org/10.1029/2010JA016399

Brandt, D. A., Bussy-Virat, C. D., & Ridley, A. J. (2020). A simple method for correcting empirical model densities during geomagnetic storms using satellite orbit data. *Space Weather*, 18(12), e2020SW002565. https://doi.org/10.1029/2020SW002565

Bruinsma, S. (2015). The DTM-2013 thermosphere model. *Journal of Space Weather and Space Climate*, 5(27), A1. https://doi.org/10.1051/swsc/2015001

Bruinsma, S., Tamagnan, D., & Biancale, R. (2004). Atmospheric densities derived from CHAMP/STAR accelerometer observations. *Planetary and Space Science*, 52(4), 297–312. https://doi.org/10.1016/j.pss.2003.11.004

Bruinsma, S., & Forbes, J. M. (2007). Global observation of traveling atmospheric disturbances (TADs) in the thermosphere. *Geophysical Research Letters*, 34(14), L14103. https://doi.org/10.1029/2007GL030243

Bussy-Virat, C. D., Ridley, A. J., & Getchius, J. W. (2018). Effects of uncertainties in the atmospheric density on the probability of collision between space objects. *Space Weather*, 16(5), 519–537. https://doi.org/10.1029/2017SW001705

Cai, X., Burns, A. G., Wang, W., Qian, L., Solomon, S. C., Eastes, R. W., et al. (2021). Investigation of a neutral “tongue” observed by GOLD during the geomagnetic storm on May 11, 2019. *Journal of Geophysical Research: Space Physics*, 126(6), e2020JA028817. https://doi.org/10.1029/2020JA028817

Chen, G. M., Xu, J., Wang, W., Lei, J., & Burns, A. G. (2012). A comparison of the effects of CIR- and CME-induced geomagnetic activity on thermospheric densities and spacecraft orbits: Case studies. *Journal of Geophysical Research*, 117(A8), A08315. https://doi.org/10.1029/2012JA017782

Dang, T., Li, X., Luo, B., Li, R., Zhang, B., Pham, K., et al. (2022). Unveiling the space weather during the Starlink satellites destruction event on 4 February 2022. *Space Weather*, 20(8), e2022SW003152. https://doi.org/10.1029/2022SW003152

Deng, Y., Fuller-Rowell, T. J., Akmaev, R. A., & Ridley, A. J. (2011). Impact of the altitudinal Joule heating distribution on the thermosphere. *Journal of Geophysical Research*, 116(A5), A05313. https://doi.org/10.1029/2010JA016019

Doornbos, E. (2012). *Thermospheric density and wind determination from satellite dynamics*. Springer. https://doi.org/10.1007/978-3-642-25129-0

Eastes, R. W., McClintock, W. E., Burns, A. G., Anderson, D. N., Andersson, L., Aryal, S., et al. (2020). Initial observations by the GOLD mission. *Journal of Geophysical Research: Space Physics*, 125(7), e2020JA027823. https://doi.org/10.1029/2020ja027823

Emmert, J. T. (2015). Thermospheric mass density: A review. *Advances in Space Research*, 56(5), 773–824. https://doi.org/10.1016/j.asr.2015.05.038

Emmert, J. T., Drob, D. P., Picone, J. M., Siskind, D. E., Jones, M., Jr., Mlynczak, M. G., et al. (2020). NRLMSIS 2.0: A whole-atmosphere empirical model of temperature and neutral species densities. *Earth and Space Science*, 7(3), e2020EA001321. https://doi.org/10.1029/2020EA001321

Fang, T.-W., Kubaryk, A., Goldstein, D., Li, Z., Fuller-Rowell, T., Millward, G., et al. (2022). Space weather environment during the SpaceX Starlink satellite loss in February 2022. *Space Weather*, 20(11), e2022SW003193. https://doi.org/10.1029/2022SW003193

Fang, X., Randall, C. E., Lummerzheim, D., Solomon, S. C., Mills, M. J., Marsh, D. R., et al. (2008). Electron impact ionization: A new parameterization for 100 eV to 1 MeV electrons. *Journal of Geophysical Research*, 113(A9), A09311. https://doi.org/10.1029/2008JA013384
Fedrizzi, M., Fuller-Rowell, T. J., & Codrescu, M. V. (2012). Global Joule heating index derived from thermospheric density physics-based modeling and observations. *Space Weather*, 10(3), S03001. https://doi.org/10.1029/2011SW000724

Fuller-Rowell, T. J. (1981). A three-dimensional time-dependent global model of the thermosphere (PhD thesis). University College London.

Fuller-Rowell, T. J., Akmaev, R., Wu, F., Anghel, A., Maruyama, N., Anderson, D. N., et al. (2008). Impact of terrestrial weather on the upper atmosphere. *Geophysical Research Letters*, 35(9), L09808. https://doi.org/10.1029/2007GL032911

Fuller-Rowell, T. J., Codrescu, M. V., Rishbeth, H., Moffett, R. J., & Quegan, S. (1996). On the seasonal response of the thermosphere and ionosphere to geomagnetic storms. *Journal of Geophysical Research*, 101(A2), 2343–2353. https://doi.org/10.1029/95JA01614

Harding, B. J., Wu, Y.-J., Alken, P., Yamazaki, Y., Triplett, C. C., Immel, T. J., et al. (2022). Impacts of the January 2022 Tonga volcanic eruption on the ionospheric dynamo: ICON-MIGHTI and Swarm observations of extreme neutral winds and currents. *Geophysical Research Letters*, 49(9), e2022GL098577. https://doi.org/10.1029/2022GL098577

He, J. (2023). Comparison of empirical and theoretical models of the thermospheric density enhancement during the 3-4 February 2022 geomagnetic storm [Dataset]. OSF. https://doi.org/10.17605/OSF.I0/SDPC4

He, J., Yue, X., & Ren, Z. (2021). The impact of assimilating ionosphere and thermosphere observations on neutral temperature improvement: Observing system simulation experiments using EnKF. *Space Weather*, 19(10), e2021SW002844. https://doi.org/10.1029/2021SW002844

Heelis, R. A., Lowell, J. K., & Spiro, R. W. (1982). A model of the high-latitude ionospheric convection pattern. *Journal of Geophysical Research*, 87(A8), 6339–6345. https://doi.org/10.1029/JA087iA08p06339

Huang, Y., Richmond, A. D., Deng, Y., & Roble, R. (2012). Height distribution of Joule heating and its influence on the thermosphere. *Journal of Geophysical Research*, 117(A8), A08334. https://doi.org/10.1029/2012JA017885

Kataoka, R., Shiota, D., Fujiwara, H., Jin, H., Tao, C., Shinagawa, H., & Miyoshi, Y. (2022). Unexpected space weather causing the re-entry of 38 Starlink satellites in February 2022. *Journal of Space Weather and Space Climate*, 12, 41. https://doi.org/10.31223/X5GH0X

Knipp, D. J., Eriksson, S., Kilcommons, L., Crowley, G., Lei, J., Hairston, M., & Drake, K. (2011). Extreme Poynting flux in the dayside thermosphere: Examples and statistics. *Geophysical Research Letters*, 38(16), L16102. https://doi.org/10.1029/2011GL048302

Knipp, D. J., Tobiska, W. K., & Emery, B. A. (2004). Direct and indirect thermospheric heating sources for solar cycles 21–23. *Solar Physics*, 224(1–2), 495–505. https://doi.org/10.1007/s11207-005-6393-4

Laskar, F. I., Eastes, R. W., Codrescu, M. V., Evans, J. S., Burns, A. G. A. G., Wang, W., et al. (2021a). Response of GOLD retrieved thermospheric temperatures to geomagnetic activities of varying magnitudes. *Geophysical Research Letters*, 48(15), e2021GL093905. https://doi.org/10.1029/2021GL093905

Laskar, F. I., Pedatella, N. M., Codrescu, M. V., Eastes, R. W., Evans, J. S., Burns, A. G., & McClintock, W. (2021b). Impact of GOLD retrieved thermospheric temperatures on a whole atmosphere data assimilation model. *Journal of Geophysical Research: Space Physics*, 126(1), e2020JA028646. https://doi.org/10.1029/2020JA028646

Laskar, F. I., Sutton, E. K., Lin, D., Greer, K. R., Aryal, S., Cai, X., et al. (2022). Thermospheric temperature and density variability during 3 to 4 February 2022 minor geomagnetic storm: The SpaceX satellite loss event. *Space Weather*. https://doi.org/10.1002/essoar.10511973.1

Lei, J., Thayer, J. P., Burns, A. G., Lu, G., & Deng, Y. (2010). Wind and temperature effects on thermosphere mass density response to the November 2004 geomagnetic storm. *Journal of Geophysical Research*, 115(A5), A05303. https://doi.org/10.1029/2009JA014754

Lei, J., Thayer, J. P., Lu, G., Burns, A. G., Wang, W., Sutton, E. K., & Emery, B. A. (2011). Rapid recovery of thermosphere density during the October 2003 geomagnetic storms. *Journal of Geophysical Research*, 116(A3), A03306. https://doi.org/10.1029/2010JA016164

Li, J., Wang, W., Lu, J., Yue, J., Burns, A. G., Yuan, T., et al. (2019). A modeling study of the responses of mesosphere and lower thermosphere (MLT) winds to geomagnetic storms at middle latitudes. *Journal of Geophysical Research: Space Physics*, 124(5), 3666–3680. https://doi.org/10.1029/2019ja026533

Lin, D., Sorathia, K., Wang, W., Merkin, V., Bao, S., Pham, K., et al. (2021). The role of diffuse electron precipitation in the formation of subauroral polarization streams. *Journal of Geophysical Research: Space Physics*, 126(12), e2021JA029792. https://doi.org/10.1029/2021JA029792

Lin, D., Wang, W., Garcia-Sage, K., Yue, J., Merkin, V., McInnerney, J. M., et al. (2022). Thermospheric neutral density variation during the “SpaceX” storm: Implications from physics-based whole geospace modeling. *Space Weather*, 20(12), e2022SW003254. https://doi.org/10.1029/2022SW003254

Liu, H., & Lühr, H. (2005). Strong disturbance of the thermospheric density due to storms: CHAMP observations. *Journal of Geophysical Research*, 110, A09S29. https://doi.org/10.1029/2004JA01090

Liu, H.-L., Bardeen, C. G., Foster, B. T., Lauritzen, P., Liu, J., Lu, G., et al. (2018). Development and validation of the whole atmosphere community climate model with thermosphere and ionosphere extension (WACCM-X 2.0). *Journal of Advances in Modeling Earth Systems*, 10(2), 381–402. https://doi.org/10.1002/2017MS001232

Lu, G., Richmond, A. D., Lühr, H., & Paxton, L. (2016). High-latitude energy input and its impact on the thermosphere. *Journal of Geophysical Research: Space Physics*, 121(7), 7108–7124. https://doi.org/10.1002/2015JA022294

Marcos, F. A. (2006). New satellite drag modeling capabilities. In *44th AIAA aerospace sciences meeting and exhibit* (pp. 01731–3010).

Maruyama, N., Sun, Y.-Y., Richards, P., Middlecoff, J., Fang, T.-W., Fuller-Rowell, T. J., et al. (2016). A new source of the mid-latitude ionospheric peak density structure revealed by a new Ionosphere-Plasmasphere Model. *Geophysical Research Letters*, 43(6), 2429–2435. https://doi.org/10.1002/2015GL067312

Maute, A., Lu, G., Knipp, D. J., Anderson, B. J., & Vines, S. K. (2022). Importance of lower atmospheric forcing and magnetosphere-ionosphere coupling in simulating neutral density during the February 2016 geomagnetic storm. *Frontiers in Astronomy and Space Sciences*, 9, 932748. https://doi.org/10.3389/fspas.2022.932748

Mayr, H. G., Harris, I., & Spencer, N. W. (1978). Some properties of upper atmosphere dynamics. *Reviews of Geophysics*, 16(4), 539–565. https://doi.org/10.1029/RG016i004p00539

Mehta, P. M., Linares, R., & Sutton, E. K. (2019). Data-driven inference of thermosphere composition during solar minimum conditions. *Space Weather*, 17(9), 1364–1379. https://doi.org/10.1029/2019SW002264

Oliveira, D. M., Zesta, E., Schuck, P. W., & Sutton, E. K. (2017). Thermosphere global time response to geomagnetic storms caused by coronal mass ejections. *Journal of Geophysical Research: Space Physics*, 122(10), 10762–10782. https://doi.org/10.1002/2017JA024006

Pham, K. H., Zhang, B., Sorathia, K., Dang, T., Wang, W., Merkin, V., et al. (2022). Thermospheric density perturbations produced by traveling atmospheric disturbances during August 2005 storm. *Journal of Geophysical Research: Space Physics*, 127(2), e2021JA030071. https://doi.org/10.1029/2021JA030071

Prölss, G. W. (1976). On explaining the negative phase of ionospheric storms. *Planetary and Space Science*, 24(6), 607–609. https://doi.org/10.1016/0032-0633(76)90140-9

Prölss, G. W. (1980). Magnetic storm associated perturbation of the upper atmosphere: Recent results obtained by satellite-borne gas analyzers. *Reviews of Geophysics*, 18(1), 183–202. https://doi.org/10.1029/RG018i001p00183
Reigber, C., Lühr, H., & Schwintzer, P. (2002). CHAMP mission status. *Advances in Space Research*, 30(2), 129–134. https://doi.org/10.1016/S0273-1177(02)00276-4

Richmond, A. D., Ridley, E. C., & Roble, R. G. (1992). A thermosphere/ionosphere general circulation model with coupled electrodynamics. *Geophysical Research Letters*, 19(6), 601–604. https://doi.org/10.1029/92GL00401

Roble, R. G., Dickinson, R. E., & Ridley, E. C. (1982). Global circulation and temperature structures of thermosphere with high-latitude plasma convection. *Journal of Geophysical Research*, 87(A3), 1599–1614. https://doi.org/10.1029/ja087ia03p01599

Siemes, C., de Teixeira da Encarnação, J., Doornbos, E., van den IJssel, J., Kraus, J., & Pereštý, R. (2016). Swarm accelerometer data processing from raw accelerations to thermospheric neutral densities. *Earth Planets and Space*, 68(1), 1–16. https://doi.org/10.1186/s40623-016-0474-5

Smith, A. K., Pedatella, N. M., Marsh, D. R., & Matsuo, T. (2017). On the dynamical control of the mesosphere-lower thermosphere by the lower and middle atmosphere. *Journal of the Atmospheric Sciences*, 74(3), 933–947. https://doi.org/10.1175/jas-d-16-0226.1

Solomon, S. C., & Qian, L. (2005). Solar extreme-ultraviolet irradiance for general circulation models. *Journal of Geophysical Research*, 110(A10), A10306. https://doi.org/10.1029/2005JA011160

Sutton, E. K. (2018). A new method of physics-based data assimilation for the quiet and disturbed thermosphere. *Space Weather*, 16(6), 736–753. https://doi.org/10.1002/2017SW001785

Sutton, E. K., Forbes, J. M., & Nerem, R. S. (2005). Global thermospheric neutral density and wind response to the severe 2003 geomagnetic storms from CHAMP accelerometer data. *Journal of Geophysical Research*, 110(A9), A09S40. https://doi.org/10.1029/2004JA010985

Sutton, E. K., Thayer, J. P., Pilinski, M. D., Mutschler, S. M., Berger, T. E., Nguyen, V., & Masters, D. (2021). Toward accurate physics-based specifications of neutral density using GNSS-enabled small satellites. *Space Weather*, 19(6), e2021SW002736. https://doi.org/10.1029/2021SW002736

Tapley, B. D., Bettadpur, S., Watkins, M., & Reigber, C. (2004). The gravity recovery and climate experiment: Mission overview and early results. *Geophysical Research Letters*, 31(9), L09607. https://doi.org/10.1029/2004GL019920

van den IJssel, J., Doornbos, E., Iorfida, E., March, G., Siemes, C., & Montenbruck, O. (2020). Thermosphere densities derived from Swarm GPS observations. *Advances in Space Research*, 65(7), 1758–1771. https://doi.org/10.1016/j.asr.2020.01.004

van den IJssel, J., & Visser, P. (2007). Performance of GPS-based accelerometry: CHAMP and GRACE. *Advances in Space Research*, 39(10), 1597–1603. https://doi.org/10.1016/j.asr.2006.12.027

Villain, J. P. (1980). Traitement des données brutes de l'accéléromètre Cactus. Etude des perturbations de moyenne échelle de la densité thermosphérique. *Annales Geophysicae*, 36, 41–49.

Wang, W., Burns, A. G., & Liu, J. (2021). Upper thermospheric winds. In W. Wang, Y. Zhang, & L. J. Paxton (Eds.), *Upper atmosphere dynamics and energetics*. https://doi.org/10.1002/9781119815631.ch3

Weimer, D. R. (2005). Improved ionospheric electrodynamic models and application to calculating Joule heating rates. *Journal of Geophysical Research*, 110(A5), A05306. https://doi.org/10.1029/2004JA010884

Wright, D. (2007). Space debris. *Physics Today*, 60(10), 35–40. https://doi.org/10.1063/1.2800252

Xu, J. Y., Smith, A. K., Wang, W., Jiang, G., Yuan, W., Gao, H., et al. (2013). An observational and theoretical study of the longitudinal variation in neutral temperature induced by aurora heating in the lower thermosphere. *Journal of Geophysical Research: Space Physics*, 118(11), 7410–7425. https://doi.org/10.1002/2013JA019144

Yue, J., Yu, W., Pedatella, N. M., Bruinsma, S., Wang, N., & Liu, H. (2022). Contribution of the lower atmosphere to the day-to-day variation of thermospheric density. *Advances in Space Research*. https://doi.org/10.1016/j.asr.2022.06.011

Zesta, E., & Oliveira, D. M. (2019). Thermospheric heating and cooling times during geomagnetic storms, including extreme events. *Geophysical Research Letters*, 46(22), 12739–12746. https://doi.org/10.1029/2019GL085120

Zhang, Y., Paxton, L. J., Schaefer, R., & Swartz, W. H. (2022). Thermospheric conditions associated with the loss of 40 Starlink satellites. *Space Weather*, 20(10), e2022SW003168. https://doi.org/10.1029/2022SW003168

Zhu, Q., Deng, Y., Sheng, C., Anderson, P., & Bukowski, A. (2022). Impact of soft electron precipitation on the thermospheric neutral mass density during geomagnetic storms: GITM simulations. *Geophysical Research Letters*, 49(11), e2021GL097260. https://doi.org/10.1029/2021GL097260