© L. Perrone & A. Mikhailov, Published by EDP Sciences 2025

1 Introduction

It is well known that the Earth’s ionosphere and neutral upper atmosphere (thermosphere) are heavily impacted by space weather events underlining the critical importance of accurate nowcasting and forecasting of the thermosphere’s state for the safety and operation of low-Earth orbiting (LEO) satellites. During geomagnetic storms, the increased neutral density can cause drag forces to rise by more than 100% in extreme cases. This creates strong orbit perturbations, leading to unpredictable changes in a satellite’s trajectory, which can degrade its orbit, cause it to re-enter the atmosphere prematurely, or in some cases, even lead to the complete loss of the satellite (e.g., Berger et al., 2023; Hapgood et al., 2022).

On February 3, 2022, at 18:13 UTC (12:45 LT), SpaceX launched 49 satellites into orbit at altitudes between 210 and 320 km from the Kennedy Space Centre in Florida (28.52°N, 80.67°W). This launch occurred during the recovery phase of a moderate geomagnetic storm, which nearly coincided with the peak in neutral density. Shortly after the launch, another moderate geomagnetic storm took place, resulting in the loss of 38 satellites due to increased atmospheric drag (Hapgood et al., 2022). This highlights the importance of monitoring the thermosphere. Several studies have analyzed these geomagnetic storms (Dang et al., 2022; Fang et al., 2022; Hapgood et al., 2022; Lin et al., 2022; Berger et al., 2023; He et al., 2023; Laskar et al., 2023) by using observational data, empirical models, and first-principles model simulations to study variations in neutral density during the event.

According to official SpaceX communications, “onboard GPS suggested that the storm’s speed and severity caused atmospheric drag to increase by up to 50% compared to previous launches” (He et al., 2023). Analysis of data from the Swarm-A satellite (440 km altitude) by He et al. (2023) and Lin et al. (2022) revealed a neutral density increase of approximately 90% at 11:00 UTC on February 3, and 110% at 21:00 UTC on February 4, compared to the reference baseline. Swarm-A recorded a gradual decrease in density, returning to pre-storm levels by February 7 during the recovery phase of the storm. Lin et al. (2022), using data from the Swarm-B satellite (515 km altitude), found a neutral density increase of 118% and 124% during the two storms, respectively. Dang et al. (2022), using accelerometer data from the Swarm-C satellite, observed density increases, with peaks at 12:00 UTC of 1.5 × 10⁻¹² kg/m³ on February 3 and 1.75 × 10⁻¹² kg/m³ on February 4.

He et al. (2023) and Berger et al. (2023) analyzed data from a pair of low-drag GRACE-FO satellites in polar orbits at altitudes of 482–535 km in the 06:00/18:00 LT sector. They reported neutral density increases of approximately 90–94% for the first storm and 120–163% for the second storm, respectively. The differences in results are attributed to the use of different reference levels in their analyses.

Empirical models such as NRLMSIS 2.0 (Laskar et al., 2023), HASDM (Berger et al., 2023), DTM2013 (He et al., 2023), and physics-based models like CTIPe (He et al., 2023), MAGE (Lin et al., 2022), WAM-IPE (Fang et al., 2022), GAIA (Kataoka et al., 2022), and TIGCM (Dang et al., 2022) generally underestimated the observed neutral density during this event.

The objectives of our analysis can be summarized as follows:

1. To describe the state of the ionosphere-thermosphere system in different longitudinal sectors at mid-latitudes during two geomagnetic storms in February 2022.

2. To retrieve neutral composition ([O], [N₂], [O₂]), exospheric temperature (Tex), and vertical plasma drift during the events using ground-based ionosonde data and Swarm-C satellite neutral density observations.

3. To compare the retrieved neutral density with modern empirical thermospheric models of neutral composition for the periods in question.

4. To quantify the role of individual thermospheric parameters in the increased neutral density observed during the February 2022 storms.

5. To explain the peculiarities of neutral composition and foF₂ variations in different longitudinal sectors during the storm periods.

2 Observations

Observations at nine ionosonde stations were used in our analysis. Six of them are located in the Northern and three in the Southern Hemispheres (Table 1). Observed bottom-side Ne(h) profiles along with Swarm-C neutral density observations allowed us to retrieve thermospheric parameters at these stations for our analysis.

Swarm is a constellation of three satellites (A, B and C) provide global measurements of magnetic field, electron and neutral density. Swarm-C has the same orbit as Swarm-A and currently it is the only Swarm satellite providing accelerometer-based density observations (Visser & van den Ijssel, 2016; van den Ijssel et al., 2020, Iorfida et al., 2023, Siemes et al., 2016).

3 Method to retrieve thermospheric parameters

The noontime mid-latitude ionospheric F-region manifests the state of the surrounding neutral atmosphere and the intensity of solar ionizing radiation. Therefore, it is possible to solve the inverse problem of aeronomy to retrieve basic aeronomic parameters responsible for the formation of the daytime mid-latitude F-region. Recently such a method THERION has been proposed by Perrone & Mikhailov (2018) to retrieve a self-consistent set of noon time neutral composition ([O], [N₂], [O₂] concentrations), temperature Tex along with vertical plasma drift W related to thermospheric winds, and total solar EUV ionizing flux.

The formation mechanism of the midlatitude daytime F-layer includes photoionization of neutral [O], [O₂], [N₂] species by solar EUV with λ < 105 nm, plasma transport by diffusion, and thermospheric winds and plasma recombination in the chain of ion-molecular reactions. The equations used in model calculations may be found in Mikhailov et al., (2012).

The input parameters include observed plasma frequency fp at 180 km for (10, 11, 12, 13, 14) LT, noontime foF₂ and satellite neutral gas density as fitted parameters. Daytime Swarm-C observations were used for our analysis and measured neutral densities were reduced to the ionosonde locations and 12 LT using the MSISE00 thermospheric model (Picone et al., 2002). The height of satellite orbit was kept unchanged during this reduction not to introduce an additional uncertainty related to the limitation of the MSISE00 model. The inclusion of neutral density into the retrieval process increases the reliability of the obtained results. In this case the inferred aeronomic parameters are found consistently with the observed neutral gas density. Standard indices of solar (F_10.7) and geomagnetic (Ap) activity are also used as the input parameters. Such indices are necessary to specify the starting state of the thermosphere and solar EUV in the fitting process. Neutral density from NRLMSIS 2.0 (Emmert et al., 2020) and JB2008 (Bowman et al., 2008) empirical climatological models of the thermosphere were used for a comparison.

Alongside accelerometer data, thermospheric density for all Swarm satellites has been derived using high-quality Swarm GPS data, employing a precise orbit determination approach (van den Ijssel et al., 2020). The resulting density data are included in the European Space Agency’s Level 2 Swarm products. A comparison of these two types of data for Rome station in the storm period manifests good coincidence (Table 2).

4 Results

The analyzed period was characterized by two moderate magnetic storms on February 03 with the minimum Dst = −66 nT and on February 04 with Dst = −62 nT and maximal 3-hour ap reaching 56 nT (Fig. 1). Solar activity was moderate during the analyzed period with F_10.7 ~ 127 sfu (1 sfu = 10⁻²² W m⁻² Hz⁻¹). After the first storm SpaceX launched 49 satellites at 18:13 UTC into the initial orbit around 210 km altitude. The 7th of February was the date when 38 failed Starlink satellites burned up on the re-entry.

Figure 1 Observed IMF-Bz, F10.7, Dst, and ap-3 hour variations for the February 01–09, 2022 period.

4.1 Electron density behaviour in the F₂ layer maximum during February geomagnetic storms

Stations were grouped in three longitudinal sectors – American (Millstone Hill, Eglin, Bahia Blanca), European (Juliusruh, Rome, Hermanus), and Japanese (Petropavlovsk, Wakkanai, Hobart) ones.

The terms “near-pole” and “far-from-pole” stations were introduced by Rishbeth & Müller-Wodarg (1999) to describe a station’s position relative to the magnetic pole. According to this classification, stations in the North American and Australian sectors are considered as “near-pole,” while those in the Eurasian sector are classified as “far-from-pole.” Rishbeth & Müller-Wodarg (1999) also noted that stations in these longitudinal sectors experience different levels of auroral heating due to their varying proximities to the magnetic pole (Fig. 2).

Figure 2 Magnetic inclination. Ionosonde stations used in the paper: near-pole stations (red asterisks), far-from-pole (blue asterisks).

The F₂-layer storm-time theory (Fuller-Rowell et al., 1994; Prölss, 1995; Forbes et al., 1996; Rishbeth, 1998; Field et al., 1998; Rishbeth & Müller-Wodarg, 1999; Rishbeth et al., 2000) suggests that the selected stations should exhibit different foF₂ storm-time variations based on factors such as season, storm phase, and their location within the “near-pole” or “far-from-pole” longitudinal sectors. Figures 3–5 display different foF₂ diurnal variations in accordance with this F₂-layer storm-time theory.

Figure 3 Diurnal variations of foF2 in the American longitudinal sector during the February 01–09, 2022 period. Solid pink curve – monthly median foF2. Vertical blue line – the moment of SpaceX launch. 3-hour ap indices (bottom panel) are given for the convenience of discussion.

A well-defined night-time negative foF₂ storm phase may be followed by positive daytime foF₂ deviations in the winter hemisphere at the “near-pole” stations. However, this daytime positive foF₂ effect is absent at the “far-from-pole” stations. Positive foF₂ deviations are not observed in the summer “near-pole” sector. Interestingly, Bahia Blanca, located in the summer hemisphere, exhibits very unusual foF₂ variations.

All of these foF₂ peculiarities require special analysis, which is provided in the Discussion section. Observed foF₂ variations should be analyzed in conjunction with the retrieved thermospheric parameters, given their close relationship.

4.2 Retrieved thermospheric parameters

The thermospheric parameters retrieved using our THERION method (Perrone & Mikhailov, 2018) are presented in Table 3. These results include exospheric temperature, concentrations of atomic and molecular oxygen, and molecular nitrogen at 250 km altitude, as well as vertical plasma drift related to thermospheric neutral winds, and neutral density at 250 km, denoted as ρ = m_o[O] + m_o2[O₂] + m_N2[N₂]. Normally, the contribution of molecular oxygen to total neutral density is small, but it may become more significant under storm conditions. Therefore, [O₂] is included in Tables 3 and 4.

Negative vertical plasma drift, W = V_nx sin I cos I, I is magnetic inclination, indicates that the effective meridional wind, V_nx, is poleward, while positive W corresponds to equatorward V_nx.

Two thermospheric empirical models were used for comparison with the retrieved ρ₂₅₀. NRLMSIS 2.0 (Emmert et al., 2020) is the latest version of the MSIS series, while JB2008 (Bowman et al., 2008) has shown the best results in comparison with CHAMP observations (Shim et al., 2012, Fig. 3).

5 Discussion

Our analysis of the Thermosphere-Ionosphere System’s response to the moderate geomagnetic storms in February 2022, which resulted in the loss of 38 SpaceX Starlink satellites, did not reveal any extraordinary features of this event. The observed storm-time variations in thermospheric and related ionospheric parameters are consistent with the current understanding of the F₂-layer storm mechanism (Fuller-Rowell et al., 1994; Prölss, 1995; Forbes et al., 1996; Rishbeth, 1998; Field et al., 1998; Rishbeth & Müller-Wodarg, 1999; Rishbeth et al., 2000).

However, our analysis also uncovered some interesting features in the storm-time variations of foF₂ that are not directly related to this mechanism. While these peculiarities in foF₂ and neutral composition variations have been discussed in previous literature (Fuller-Rowell et al., 1994; Prölss, 1995; Forbes et al., 1996; Rishbeth, 1998; Field et al., 1998; Rishbeth & Müller-Wodarg, 1999; Rishbeth et al., 2000) in the context of geomagnetic storms, our analysis shows that they also occurred during the analyzed period as a response of the thermosphere and F₂-layer to the geomagnetic storms.

In accordance with the storm mechanism, the mid-latitude daytime F₂-layer shows a dependence on the intensity of incident solar EUV radiation and the state of the surrounding thermosphere. During the period in question, variations in solar EUV were small (within 3%), according to the EUV model by Nusinov et al. (2021). This is also reflected in the variations of F_10.7 (Fig. 1); therefore, all observed effects should be attributed to geomagnetic activity. The latter is the most important factor, as neutral composition and temperature determine neutral density, which in turn affects satellite drag.

Storm-time thermospheric parameters depend on the intensity of the geomagnetic disturbance, the longitudinal sector (“near-pole” or “far-from-pole”), geomagnetic latitude, local time, and the prehistory of geomagnetic activity. The neutral density observed by satellites reflects a complex interplay between variations in neutral temperature which always increases during geomagnetic disturbances and neutral composition, which exhibits complex storm-time variations.

Neutral density is primarily composed of atomic oxygen and molecular nitrogen at the altitudes of SpaceX Starlink satellites. While molecular nitrogen, as an inertial species, generally follows variations in neutral temperature, atomic oxygen exhibits very complex spatial and temporal changes. With our method (Perrone & Mikhailov, 2018), we can track such storm-time daily variations of thermospheric parameters across different seasons and longitudinal sectors. The results are presented in Table 3.

Firstly, it should be noted that the period preceding the analyzed timeframe in the paper was not magnetically quiet, thus avoiding the recovery phase of the storm (Fig. 1) February 9 was the only day that could be used as a reference. We will discuss variations in thermospheric parameters at an altitude of 250 km. This altitude was selected for two reasons: on the one hand, it is close to the altitudes of SpaceX Starlink satellite orbits, and on the other hand, it is near the peak of the daytime F₂-layer, which is relevant for foF₂ variations in relation to the retrieved thermospheric parameters. Let us now consider these variations across different longitudinal sectors.

5.1. Northern Hemisphere-American longitudinal sector

Millstone Hill is a “near-pole” station located in the winter Northern Hemisphere, where the maximal neutral density (ρ) increase of 43.8% occurred on February 4 (Table 4). The previous day, February 3, was marked by a moderate 31% increase in ρ, in response to the first magnetic storm, which took place during the nighttime hours (Fig. 3).

The disturbed neutral composition, characterized by a decreased O/N₂ ratio, was transported from the auroral zone to the middle latitudes by both background and disturbed equatorward thermospheric wind circulation. This resulted in a negative NmF₂ storm phase that began after midnight and lasted throughout the entire day on February 3, lasting until sunrise on February 4 (Fig. 3). A 75.9% increase in [N₂]₂₅₀ (Table 4), with nearly unchanged [O]₂₅₀, indicates the O/N₂ decrease on February 3.

At Millstone Hill during daytime hours, we observe a strong poleward thermospheric wind (negative vertical plasma drift, Table 3) on all disturbed days. This indicates that the magnetic disturbances were not strong enough to override the solar-driven poleward thermospheric circulation during daytime hours. According to the previously mentioned F₂-layer storm mechanism, we would expect the interaction of storm-induced and background solar-driven circulations, leading to downwelling of neutral gas. This downwelling should result in an increase in atomic oxygen. In fact, we observe a 3.4–34% increase in [O]₂₅₀ from February 3 to 5 compared to the reference day (Table 4). Part of this increase in [O] can be attributed to the rise in exospheric temperature (Tex), and part to the downwelling process. However, we must remember that even this 3.4% increase in [O] occurred under elevated Tex = 1053 K, and we still observe a decreased O/N₂ ratio that was transported to middle latitudes during the nighttime hours on February 3, contributing to the negative NmF₂ storm phase on February 3 (Fig. 3).

The true downwelling effect took place later, on February 4–5, with [O]₂₅₀ increasing by 24–34% (Table 4), which led to daytime foF₂ values coinciding with the monthly median foF₂ (Fig. 3).

It is interesting to examine the contributions of [O] and [N₂] to the total neutral density at Millstone Hill throughout the storm. Table 4 shows that atomic oxygen makes the largest contribution to ρ variations, both during quiet and disturbed conditions. Only on February 3, after the first storm, do the contributions from [O] and [N₂] become comparable. This is also evident in the O/N₂ ratio on February 3, where the depletion of this ratio is due to the increase in [N₂].

On February 4 and 5, our results only partially align with TIMED/GUVI observations (Zhang et al., 2022), which found an O/N₂ depletion in the North American sector, associated with an increase in [N₂]. While the O/N₂ ratio did decrease on February 3–4, by February 5, due to downwelling, the ratio increased relative to the reference day (Table 3). It should be noted that TIMED/GUVI provides O/N₂ column density measurements, while we are analyzing neutral composition at a fixed altitude.

Finally, we can conclude that the 31% increase in ρ₂₅₀ on February 3 was due to equal contributions from [O] and [N₂], whereas the contribution from [O] was larger than that from [N₂] on February 4–5 (Table 4). Therefore, the role of [O] in the ρ₂₅₀ increase became more prominent as the storm progressed. Neutral density ρ₂₅₀ calculated with the empirical models JB2008 (Bowman et al., 2008) and NRLMSIS 2.0 (Emmert et al., 2020) showed smaller variations during the storm days compared to the neutral density retrieved by THERION, with results from JB2008 being closer to our findings.

Eglin, also located in the winter “near-pole” sector, largely mirrors the storm-time thermospheric and foF₂ variations seen at Millstone Hill, but with different magnitudes due to its lower magnetic latitude (Table 1). The maximal ρ₂₅₀ increase on February 4 was only 34.3%, and 24.1% on the previous day. Disturbed neutral composition with a low O/N₂ ratio could not reach the latitude of Eglin after the February 3 storm, and the daytime F₂-region showed a slight negative storm phase (Fig. 3). This only occurred on the night of February 3–4, when solar-driven (equatorward) thermospheric circulation, combined with the storm-induced circulation, carried the disturbed neutral composition to Eglin’s latitude, resulting in a well-pronounced negative storm phase (Fig. 3). Daytime northward circulation (Table 3) then shifted the disturbed neutral composition back to higher latitudes, and Eglin exhibited a positive foF₂ storm phase on February 4.

This process was described by Prölss (1995), who stated: “The idea that composition perturbations, once generated, ‘rotate’ with the Earth is only schematic. Actually, the disturbance bulge will be pushed around by winds and may move back and forth in latitude.” It should be noted that despite a strong downward plasma drift (Table 3), a positive daytime foF₂ storm phase occurred due to a large increase in [O] – 45.2% on February 4 and 29.4% on February 5 (Table 4). The increase in [O] was attributed to downwelling, until solar-driven (poleward) circulation interacted with storm-induced (equatorward) circulation.

After the storm ended, by February 6, both Millstone Hill and Eglin exhibited quiet-time foF₂ variations, with no positive daytime storm phase (Fig. 3).

Table 4 shows that atomic oxygen contributes 60–70% to ρ₂₅₀ during storm days. Contrary to TIMED/GUVI observations (Zhang et al., 2022), the O/N₂ ratio increases from 2.78 on February 3 to 4.0 on February 5 (Table 3). Finally, we can conclude that the ρ₂₅₀ increase from 24.1% on February 3 to 34.3% on February 4 was driven by a continuous increase in [O].

Neutral density ρ₂₅₀ calculated using the empirical models JB2008 (Bowman et al., 2008) and NRLMSIS 2.0 (Emmert et al., 2020) shows smaller variations during the storm days compared to the neutral density retrieved by THERION, except on February 5, when the neutral density calculated with NRLMSIS 2.0 is larger than the other values.

5.2 Southern Hemisphere-American longitudinal sector

Bahia Blanca, located at φ = 29.2°S, is situated on the polar slope of the equatorial anomaly (EA) crest. Formally, our method for retrieving aeronomic parameters may not be fully applicable here, as there is an additional plasma influx due to the “fountain effect” under the ExB drift. Observed foF₂ variations (Fig. 3) indicate this impact. According to Mikhailov et al. (1996), “…IMF-Bz turnings to a northward direction result in a decrease (or even reversal) in the normal Sq zonal component of the electric field (eastward during the daytime and westward at night). Southward IMF-Bz excursions enhance the normal Ey both during the day and night.” Figure 1 shows a northward IMF-Bz excursion on February 2 and a southward IMF-Bz excursion on February 3. This means a decrease in upward ExB drift on February 2 and an increase on February 3.

Correspondingly, the EA crest will move closer to the geomagnetic equator on February 2 and farther from it on February 3. In the first case, Bahia Blanca will be farther from the crest maximum, resulting in lower foF₂, while in the second case, it will be closer to the crest maximum, leading to larger foF₂ (Fig. 3). The effect of the increased post-sunset upward ExB drift, associated with the southward IMF-Bz excursion, is seen on February 4, when elevated foF₂ values are observed during the evening and night-time hours (Fig. 3).

All of this indicates that Bahia Blanca was strongly affected by equatorial ExB drifts, making our method unsuitable for this station. Unfortunately, no other mid-latitude stations with ionospheric observations are available in South America for the period in question

5.3 Northern Hemisphere-European longitudinal sector

The European longitudinal sector occupies an intermediate position between the “near-pole” North American and “far-from-pole” Japanese sectors. Therefore, by analogy with Millstone Hill and Eglin, the European stations, Juliusruh and Rome, show similar winter storm-time variations in thermospheric parameters, but with smaller magnitudes. A 35% storm-time ρ₂₅₀ increase at Juliusruh and approximately 22% ρ₂₅₀ increase at Rome were due to an increase in [O], which made the main contribution (54–69%) to neutral density (Table 4). The mechanism behind this [O] increase is the same: downwelling resulting from the interaction of solar-driven (poleward) and storm-induced (equatorward) thermospheric circulations.

In conclusion, we can say that the 16–35% storm-time increase in neutral density in the winter European sector was primarily due to the increase in atomic oxygen.

Neutral densityρ₂₅₀ calculated using the empirical model JB2008 shows values close to those from THERION in Rome, while in Juliusruh, there is an underestimation. In contrast, NRLMSIS 2.0 overestimates the neutral density at both stations.

However, the difference between longitudinal sectors in relation to the storm onset resulted in different foF₂ diurnal variations. Unlike the American sector, in Europe, the storm onset on February 3 occurred during daytime hours (Fig. 4). This period is referred to as the “forbidden time” for negative ionospheric storms (Prölss & von Zahn, 1978), when disturbed neutral composition with a low O/N₂ ratio is confined to high latitudes by the strong poleward solar-driven thermospheric circulation. Our calculations show this strong poleward wind (large downward W, Table 3). As a result, Juliusruh and Rome exhibited either positive foF₂ deviations during the daytime or no disturbances at night on February 3–4 (Fig. 4). The absence of disturbances at night is a clear indication that disturbed neutral composition with a low O/N₂ ratio was confined to higher latitudes. The second storm, on February 4, occurred during the night, and a negative foF₂ storm was observed at Juliusruh (a higher-latitude station) on February 5.

Figure 4 Same as Figure 3 but for the European longitudinal sector.

5.4 Southern Hemisphere-European longitudinal sector

Hermanus, located in the summer hemisphere, exhibits quite different storm-time variations of thermospheric parameters compared to stations in the winter hemisphere. The main difference lies in thermospheric circulation, which plays a significant role in neutral composition variations. As mentioned earlier, the first storm occurred during daytime hours on February 3 (Fig. 4). In the summer hemisphere, daytime solar-driven poleward circulation is weak (W > −10 m/s), and the magnetic disturbance was strong enough to reverse it to an equatorward direction during February 3–4 (Table 3). Disturbed neutral composition was moved from the auroral zone toward middle latitudes. The O/N₂ ratio decreased to 1.36 on February 3 and 1.50 on February 4, compared to the reference day, February 9. Table 4 shows a 13–24% decrease in atomic oxygen [O]₂₅₀, despite elevated neutral temperature, while [N₂], being an inertial species, increased by 12–52% following the rise in Tex. As a result, there was a modest storm-time increase of ≤15.3% in ρ₂₅₀ (Table 4), primarily due to the increase in [N₂] compared to the reference day. Therefore, Hermanus exhibits characteristics typical of a mid-latitude summer storm-time variations in thermospheric parameters, as expected for the intermediate longitudinal sector.

Neutral density ρ₂₅₀ calculated with the empirical model JB2008 shows an underestimation with respect to THERION values while NRLMSIS2.0 overestimates the neutral density at both stations.

5.5 Southern Hemisphere-Asia longitudinal sector

Hobart, located in the “near-pole” summer sector, should exhibit more pronounced storm-time effects compared to Hermanus. Let us consider the diurnal variations of foF₂ (Fig. 5). The first magnetic storm occurred during daytime hours on February 3, and despite the “forbidden time” for negative F₂-layer storms, thermospheric circulation was reversed (W = +11.1 m/s, Table 3). Disturbed neutral composition with a low O/N₂ ratio (0.59 compared to 1.34 on the reference day) was transferred from the auroral zone to middle latitudes, resulting in a negative F₂-layer storm phase. Equatorward circulation (W = +11.0 m/s) occurred on February 4, maintaining the O/N₂ ratio at a low level of 0.48. This resulted in a pronounced negative foF₂ storm. It should be noted that the decrease in the O/N₂ ratio was so significant that the equatorward wind could not compensate for in a comparison with the negative effect on foF₂.

Figure 5 Same as Figure 3 but for the Japanese longitudinal sector.

The second storm occurred during nighttime hours on February 4/5 (Fig. 5), and the disturbed neutral composition, characterized by an O/N₂ ratio of 0.55, was transported to the middle latitudes, producing a strong negative foF₂ disturbance during the daytime. This effect was enhanced by a negative W = −9.3 m/s (Table 3) as auroral activity decreased and normal solar-driven poleward circulation began to recover.

Figure 5 also shows sharp positive peaks in foF₂ on February 3–4, which are clearly related to the passage of Travelling Atmospheric Disturbances (TADs). We do not have AE-index variations for the Southern Hemisphere to verify these increases in auroral heating. Northern Hemisphere AE-indices may not be suitable for this test due to the North-South asymmetry of geomagnetic activity, which always exists, and the larger the level of magnetic activity, the larger the asymmetry (Makarov, 2016). Nevertheless, observed diurnal variations in hmF₂ at Hobart (not shown in the paper) display peaks in hmF₂ approximately one hour before the foF₂ peaks, clearly indicating the passage of TADs.

Table 4 shows that [N₂] (via Tex) provides the main contribution to neutral density. However, the large increases in [N₂] – approximately 43% on February 3 and 81% on February 4 – resulted in only modest increases in neutral density, with a 10.7% increase on February 3 and a 33.0% increase on February 4 in ρ₂₅₀, respectively. These small storm-time neutral density increases were puzzling at the beginning of the space era and even contradicted the predictions of the Jacchia-71 thermospheric model (Jacchia, 1971), where Δlog ρ is directly related to the Kp index. The explanation for this effect lies within the framework of the F₂-layer storm-time theory (Rishbeth et al., 2000, and references therein). In fact, the modest neutral density increase is the result of an interplay between storm-decreased atomic oxygen and storm-increased molecular nitrogen concentrations, which may partly cancel each other out, leading to relatively small changes in total neutral density at lower altitudes (around 250 km).

Thus, it can be concluded that the increase in ρ₂₅₀ from February 3 to 5 was entirely due to the increase in [N₂], which was associated with the strongly elevated Tex (Table 3). However, despite this Tex increase, the ρ₂₅₀ variations were very modest, with only an ~11% increase on February 3 and ~ 33% on February 4. These increases were primarily due to the anti-phase variations of [O] and [N₂], which largely compensated each other. It is important to note that both the NRLMSIS 2.0 and JB2008 models did not reproduce the observed ρ₂₅₀ variations (Table 3) and thus cannot be used for such an analysis. Additionally, neutral density ρ₂₅₀ calculated using the empirical models JB2008 and NRLMSIS 2.0 showed smaller variations during the storm days compared to the neutral density retrieved by THERION on February 4. However, the opposite trend was observed on February 3 and 5.

5.6 Northern Hemisphere-Asia longitudinal sector

Petropavlovsk and Wakkanai are located in the winter mid-latitude “far-from-pole” longitudinal sector. The first storm occurred during daytime hours at both stations (Fig. 5). According to the F₂-layer storm-time theory, in winter, the disturbed neutral composition is confined to higher latitudes by thermospheric circulation (Prölss, 1995), preventing its penetration to middle latitudes. The observed foF₂ variations indicate a slight increase above the monthly median level, which is likely due to downwelling, as discussed in relation to other Northern Hemisphere stations.

During night-time, when the thermospheric circulation becomes equatorward, disturbed neutral composition can be transported to middle latitudes, leading to negative foF₂ deviations during the daytime. This effect is more pronounced at the higher-latitude station, Petropavlovsk, which is closer to the auroral oval than Wakkanai. The second storm, on February 4, occurred during nighttime hours, and pronounced negative foF₂ deviations were observed at both stations.

On February 5, during the daytime, a strong (W ≈ −20 m/s) poleward wind shifted the disturbed neutral composition back to higher latitudes (as explained by Prölss, 1995), and foF₂ values returned to near their monthly medians. In the evening, the circulation reversed, and a strong negative foF₂ storm effect was again observed at Petropavlovsk, where observations were available (Fig. 5).

Table 4 shows a continuous increase in neutral density ρ₂₅₀ at Wakkanai, from 9.2% on February 3 to 42.2% on February 5. This increase was driven by the rise in both [O]₂₅₀ (from 11.7 to 15.8) × 10⁸ cm⁻³ and [N₂]₂₅₀ (from 4.4 to 5.4) × 10⁸ cm⁻³ (Table 3). Despite these increases, the contributions of [O] and [N₂]_to ρ₂₅₀ remained relatively stable during this period, with [O] contributing 60–65% and [N₂] contributing 34–39%. This effect may be linked to the increase in thermospheric temperature (Tex, Table 3) and the downwelling of atomic oxygen.

Petropavlovsk, located in the “far-from-pole” sector, exhibits relatively small day-to-day variations in [O] during daytime hours, while [N₂] shows a 24–52% increase on disturbed days (Tables 3 and 4), in response to variations in Tex. Since [O] contributes the main share (54–62%) of ρ₂₅₀, this results in a moderate 15–21% increase in neutral density (Table 4), which can also be attributed to the increase in Tex.

6 Summary results

In addition to Table 4, Figure 6 clearly summarizes the changes in [O]₂₅₀, [N₂]₂₅₀, and ρ₂₅₀ for February 2–5 with respect to the reference day, February 9. All stations show an increase in neutral density; however, this increase is less than 45%, which is much lower than the previously mentioned estimates of approximately 90–94% and 120–163% reported by He et al. (2023) and Berger et al. (2023), who observed values in the 06/18 LT sector at altitudes of 482–535 km. This difference can be partly explained by the variation in observation altitudes. At altitudes around 500 km, the contributions of molecular species such as N₂ and O₂ are absent, and [O] alone represents the neutral density. However, Figure 6 shows that the increase in [O] is less than 50% at all stations. Furthermore, Hobart and Hermanus show a decrease in [O] concentrations, which should result in negative deviations of neutral density at higher altitudes. Therefore, the discrepancy with the results of He et al. (2023) and Berger et al. (2023) can be attributed to differences in the baseline level used in their estimates.

Figure 6 Changes in % of [O]250, [N2]250, and ρ250 with respect to the reference February 09 day during the February 2022 storm period.

Molecular nitrogen, as an inertial species, shows an increase at all stations in response to the storm-time increase in Tex. The maximal [N₂] increase, up to 80%, occurs at the “near-pole” stations, Millstone Hill and Hobart, while the minimal increase, around 20%, is observed at the “far-from-pole” stations, Rome and Wakkanai (Fig. 6).

Retrieved neutral densities ρ₂₅₀ from Table 3 were used for a comparison with empirical NRLMSIS 2.0 and JB2008 models in three longitudinal sectors for quiet and disturbed days.

Both empirical models manifest larger MRD for quiet days compared to disturbed ones but the largest MRD manifests the new version of MSIS model. NRLMSIS 2.0 also exhibits the largest biases for quiet days. In principle, this seems strange. However, if we compare, for example, the neutral density observed at the 559 km satellite height at Millstone Hill with the MSISE neutral density, we find the same result: namely, a significant overestimation of the neutral density calculated by MSISE during quiet days compared to the neutral density observed by Swarm (Table 6).

This is due to a strong overestimation of [O] during quiet days. On average, the JB2008 model shows better agreement with the retrieved ρ₂₅₀_than NRLMSIS 2.0. This model demonstrated the best results when compared with CHAMP neutral density observations (Shim et al., 2012, Fig. 3).

7 Conclusions

Our recently developed method THERION (Perrone & Mikhailov, 2018) to extract neutral composition (O, O₂, N₂), temperature Tex, and vertical plasma drift related to thermospheric winds was applied to ground-based ionosonde and Swarm neutral density observations to analyze the Thermosphere-Ionosphere system during the February 2022 storm period when 38 of 49 SpaceX Starlink satellites were lost. The obtained results may be formulated as follows:

1. Day-time mid-latitude thermosphere has manifested very moderate <50% neutral density perturbations at 250 km height which is close to altitudes of SpaceX Starlink satellites. It should be noted that neutral temperature and density are maximal during analyzed day-time hours.

2. The largest perturbations of thermospheric and related ionospheric parameters took place (as this was expected) in North American and Australian “near-pole” longitudinal sectors. However, due to different seasons and different thermospheric circulation the ρ₂₅₀ increase was different in the two sectors. The formation mechanisms of these ρ₂₅₀ perturbations are also different in these sectors.

3. The largest atomic oxygen increase up to (30–50)% took place in the Northern winter Hemisphere where [O] provided the main contribution to the ρ₂₅₀ increase. On the contrary, [O] demonstrated a strong storm-time depression down to –(20–40)% in the summer Hemisphere (Hobart, Hermanus) and ρ₂₅₀ increase on February 03–05 was totally due to [N₂] increase related to elevated Tex.

4. The atomic oxygen increase due to downwelling and molecular nitrogen increase due to elevated Tex work in one direction in the Northern winter Hemisphere resulting in (35–45)% increase of ρ₂₅₀. On the contrary, [O] and [N₂] work in opposite directions in the Southern summer Hemisphere compensating to a great extent the contributions of each other to ρ₂₅₀. However depending on Tex the [N₂] increase may overpower the [O] decrease resulting in large ρ_250. This happened at Hobart on February 04 when due to large Tex = 1252 K the [N₂] contribution to ρ₂₅₀ was 75.5% while [O] contributed only 21%, the rest 3.5% of contribution was related to the [O₂] increase.

5. The European longitudinal sector manifested same features as in the American longitudinal sector but with less magnitude. A 35% storm-time ρ₂₅₀ increase at Juliusruh and ~22% ρ₂₅₀ increase at Rome were due to [O] increase which provided the main (54–69)% contribution to neutral density. The mechanism of this [O] increase is the same – downwelling resulted from the interaction of solar-driven (poleward) with storm-induced (equatorward) thermosheric circulations. Finally we may conclude that a (16–35)% storm-time increase of neutral density in the winter European sector was mainly related to the atomic oxygen increase.

6. Hermanus located in summer European longitudinal sector was subjected to the same [O]/[N₂] compensation mechanism. A 13% [O] decrease despite elevated neutral temperature and (40–52)% [N₂] increase resulted in a small (10–15)% ρ₂₅₀ increase on February 03–04. Therefore Hermanus exhibited a typical of mid-latitude summer storm-time variations of thermospheric parameters.

7. In general the “far-from-pole” winter longitudinal sector (Petropavlovsk, Wakkanai) manifested very moderate <21% neutral density perturbations mainly related to Tex increase. This was also seen in foF₂ diurnal variations. The main storm effects were seen during night-time hours while daytime foF₂ were close to monthly median values.

8. The obtained results have shown a moderate impact of two geomagnetic storms in February 2022 on the Thermosphere-Ionosphere system. However this impact was sufficient to result in loss of 38 satellites highlighting that the routine monitoring of the thermosphere, as the new Fast Track data from ESA and GOLD measurements, is crucial for understanding its dynamics and for operational needs.

Acknowledgments

The editor thanks Jaroslav Chum and an anonymous reviewer for their assistance in evaluating this paper.

Funding

This research was carried within the INGV Pianeta Dinamico Project (CUP D53J19000170001), Space weather effects on the South Atlantic anomaly Region (SESAR) – 2021, funded by MUR (law 145/2018).

Data availability statement

The European Space Agency provides Swarm data (https://earth.esa.int/web/guest/swarm/data-access, while Woods provides EUV observations (http://lasp.colorado.edu/lisird/). The Rome ionospheric data are provided by INGV (https://doi.org/10.13127/eswua/hf), and the Wakkanai ionospheric data are provided by NICT (https://wdc.nict.go.jp/Ionosphere/en/data). Ionosonde data are also available through the Lowell DIDBase via GIRO (http://giro.uml.edu/). Additional data sources include NOAA SWPC (https://www.swpc.noaa.gov/), GFZ Potsdam for the Kp index (https://www.gfz-potsdam.de/en/kp-index/), and the WDC for Geomagnetism, Kyoto, for the ap geomagnetic index (http://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html).

References

All Tables

All Figures

	Figure 1 Observed IMF-Bz, F10.7, Dst, and ap-3 hour variations for the February 01–09, 2022 period.
In the text

	Figure 2 Magnetic inclination. Ionosonde stations used in the paper: near-pole stations (red asterisks), far-from-pole (blue asterisks).
In the text

	Figure 3 Diurnal variations of foF2 in the American longitudinal sector during the February 01–09, 2022 period. Solid pink curve – monthly median foF2. Vertical blue line – the moment of SpaceX launch. 3-hour ap indices (bottom panel) are given for the convenience of discussion.
In the text

	Figure 4 Same as Figure 3 but for the European longitudinal sector.
In the text

	Figure 5 Same as Figure 3 but for the Japanese longitudinal sector.
In the text

	Figure 6 Changes in % of [O]250, [N2]250, and ρ250 with respect to the reference February 09 day during the February 2022 storm period.
In the text