Unexpected space weather causing the reentry of 38 Starlink satellites in February 2022

– The accidental reentry of 38 Starlink satellites occurred in early February 2022, associated with the occurrence of moderate magnetic storms. A poorly understood structure of Coronal Mass Ejections (CMEs) caused the magnetic storms at unexpected timing. Therefore, a better understanding of minor CME structures is necessary for the modern space weather forecast. During this event, the “ up to 50% ” enhancement of air drag force was observed at ~200 km altitude, preventing the satellites ’ safety operations. Although the mass density enhancement predicted by the NRLMSIS2.0 empirical model is less than 25% under the present moderate magnetic storms, the real-time GAIA simulation showed a mass density enhancement of up to 50%. Further, the real-time GAIA simulation suggests that the actual thermospheric disturbances at 200 km altitude may occur with larger amplitude in a broader area than previously thought.


Introduction
Satellite drag is sensitive to space weather and space climate, and the social impact is high in modern society depending on thousands of satellites. There have been unfortunate accidents of satellite reentry due to the enhanced drag during large magnetic storms, including the old example of the Bastille event in July 2000 when the Japanese ASKA satellite orbiting at 440 km lost the attitude control and the solar power and finally reentered. The most recent example of Starlink satellites was somewhat surprising (Hapgood et al., 2022) because of the large number of lost satellites (38 out of 49) associated with moderate magnetic storms.
SpaceX launched the 49 Starlink satellites at 1813 UT on February 3. From the initial perigee of 210 km, the satellites were planned to raise the altitude to~340 km in the Low-Earth Orbit (LEO). The operation was a challenging attempt for SpaceX, targeting relatively low altitudes where the atmospheric drag is critical, partly to obtain future data for controlling space debris. In the press release on February 8 (https://www.spacex. com/updates/, February 8, 2022: GEOMAGNETIC STORM AND RECENTLY DEPLOYED STARLINK SATELLITES), SpaceX noted the space weather situation and operation as follows: ". . . onboard GPS suggests the escalation speed and severity of the storm caused atmospheric drag to increase up to 50% higher than during previous launches. The Starlink team commanded the satellites into a safe-mode where they would fly edge-on (like a sheet of paper) to minimize drag to effectively "take cover from the storm" and continued to work. . ." In general, Geomagnetic Disturbances (GMD) result in the Joule heating in the polar atmosphere, enhancing the mass density in the expanding thermosphere. The drag force of spacecraft is proportional to the mass density of the thermosphere. Therefore, unless we succeed in predicting the GMD, we cannot expect the satellite drag.
To investigate the thermospheric mass density variations, we consult the real-time results (Tao et al., 2020) of the Ground-to-topside model of Atmosphere and Ionosphere for Aeronomy (GAIA) (e.g., Jin et al., 2011). The GAIA can reproduce the thermospheric disturbance caused by atmospheric and magnetospheric phenomena (Miyoshi et al., , 2018Jin et al., 2012;Shinagawa et al., 2017). This paper examines our capability of the space weather forecast for this particular GMD event, mainly focusing on the puzzling parts for space weather forecasters and satellite operators. We hope this work contributes to accumulating our knowledge for future robust satellite operations.
2 Solar eruptions M1.1 class solar flare occurred at 2332 UT on January 29, 2022 (Fig. 1a), in the Active Region (AR) NOAA12936 located near the central meridian of the solar disc. Before the flare peak, a coronal dimming was observed by SDO/AIA (Pesnell et al., 2012, Lemen et al., 2012 at~2000 UT (Fig. 1b) near the northern envelope of the AR, together with several brightenings in the central part (Figs. 1c and 1d). Around 2200 UT a coronal loop in the central region expanded outward and eastward, and the M1.1-class Long Duration Event (LDE) flare occurred (Figs. 1e and 1f).
Associated with the LDE flare, a halo CME was observed at 2348 UT in SOHO/LASCO (Brueckner et al., 1995) images and 2353 UT in STEREO-A/SECCHI COR2 (Howard et al., 2008)  R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 41 images. The apparent propagation speeds of the CME in images of LASCO and COR2 were approximately 620 km/s and 440 km/s, respectively. Assuming constant speed propagation, the expected arrival time to 1 AU distance ranges from early February 2 to the end of February 2. The observed shock arrival time was 2137 UT on February 1, and most CME propagation models predicted the shock arrival time within the error range of several hours (https://kauai.ccmc.gsfc.nasa.gov/CMEscoreboard/). Note that STEREO-A was located at~35°longitude trailing behind the Earth. The CME edge is toward the northeast at LASCO, while it looks relatively isotropic at COR2, consistent with the observed early CME arrival at STEREO-A. Two CMEs appeared overlapped as marked by dashed circles in Figure 1h, one can be associated with the coronal dimming, and another can be associated with LDE flare. Although the "post-mortem" analysis can identify such a possibility, it was difficult for forecasters to utilize the information of possibly double CMEs for actual GMD forecast, as described in the following sections.

Geomagnetic disturbances
As shown in Figure 2, magnetic storms peaked at 1100 UT on February 3 (Dst = À77 nT) and 1100 UT on February 4 (Dst = À64 nT). The W-shape variation in the Dst index is characterized by the positive excursion in the middle of moderate storms at 0000 UT on February 4 (Dst peak = 4 nT). Note that the whole W-shape pattern took~2 days, essentially different from the typical two-step development of magnetic storms taking less than one day (Kamide et al., 1998).
Looking back at the Dst record for the last 20 years, we can identify similar W-shape Dst variations on July 9-10, 2005, and March 24-25, 2007. The July 2005 GMD event was caused by Figure 2. OMNI-2 hourly data of magnetic field strength (red curve is Bz in GSM coordinate system), solar wind speed V, dynamic pressure Pd, and Dst index (http://swdcwww.kugi.kyoto-u.ac.jp/dstae/index.html). The blue curve shows the result from the Burton model, where the input parameters are Bz, V, and Pd. The second panel from the bottom shows the ap index (https://www.gfz-potsdam.de/en/kp-index/), which is used as an input parameter for NRLMSIS2.0. Finally, the bottom panel shows the cross polar cap potential from the Weimer model, where the input parameters are the solar wind density, V, By, and Bz. R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 41 the arrival of CME just after a moderate CME storm, while March 2007 GMD event was caused by the arrival of the Corotating Interaction Region (CIR) after a CME storm. Therefore, the standard pictures of CME and CIR storms (Kataoka & Miyoshi, 2006) make these former examples understandable. The second Dst peak is expected for any forecasters in advance by the arrivals of shock and stream interface, respectively.
However, the February 2022 GMD event differs from the standard picture concerning the following two points. First, the leading edge of the first flux rope arrived~1 day later from the shock arrival at 2220 UT on February 1. Many experts expected that the oncoming flux rope would not likely hit the Earth anymore because the sheath duration is much longer than the average value of 10.7 h . After all, the first storm occurred with the very late arrival of the first flux rope on February 3 (yellow shaded area in Fig. 2). Another puzzling feature came next at the end of the first storm when the Earth exited the first flux rope. It seemed a reasonable (moderate and settling GMD) timing for the launch of Starlink satellites (1813 UT on February 3). Then, surprisingly, the second flux rope arrived (orange shaded area in Fig. 2), and the main phase of the second storm readily started on February 4, which must have confused the initial satellite operations. We address these two significant concerns documented above in more detail and better clarify by combining another in-situ solar wind observation data at STEREO-A at different longitude (~35°trailing behind in longitude from the DSCOVR), as shown in Section 4.
If we knew the arrival of two flux ropes in advance, the expectation and even rough prediction of these two moderate storms themselves would not have been likely a difficult task. For example, the empirical Burton model (Burton et al., 1975;O'Brien & McPherron, 2000) roughly works to predict the Dst variation (Fig. 2, blue curve). Further, suppose we knew the occurrence of these moderate storms. In that case, the thermospheric response (mass density enhancement at desired altitude) can also be roughly predictable by empirical models such as NRLMSIS2.0, as a function of the ap index. However, real-time GAIA simulation helps investigate detailed dynamic R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 41 response and somewhat significant mass density enhancement, as discussed in Section 5.

Interplanetary structure of coronal mass ejections
We can recognize the different looks of the overall shock-CME passage at different longitude by comparing the in-situ observation of the Interplanetary Magnetic Field (IMF), as shown in Figures 3 and 4. Interestingly, the STEREO-A data looks like a standard picture of shock-CME passage (Kataoka & Miyoshi, 2006), which was natural, especially at the headon location against the propagating CME. However, at the Earth position, as shown in Figure 4, the overall structure looks elongated and nearly doubled in time or space, which indicates the flank-side passage of CME (e.g., Marubashi & Lepping, 2007). Note also that the active region was located near the central meridian. Still, the CME was directed toward STEREO rather than the Earth and DSCOVR, as indicated in Figures 1g and 1h.
We separate the flux rope structure into two different flux ropes, as indicated by the yellow and orange shaded areas in Figure 3. The border of the first and second flux ropes is selected to be an apparent discontinuity in the magnetic field rotation.
Similarities also exist in the first flux rope. The West-North-East (and to South) rotation at STEREO-A and North-East-South (and to West) at DSCOVR are in the same helicity sense, just tilted by~90°. Such a relationship of~90°tilt of the flux rope is also consistent with the grim picture of head-on and flank-side passages of the same flux rope. For the second flux rope identified in DSCOVR data, we can identify a similar trailing part of the CME (orange shaded, small B field rotation merging to the Parker spiral) in STEREO-A data (Fig. 3). The orange shaded area in Figure 3 shows the magnetic field change from positive Bz to negative with negative By and positive Bx, similar to the orange shaded area in Figure 4. In Figure 4, at the back of the first flux rope, we observe a short duration kink of the magnetic field, which can also be related to the second flux rope at STEREO-A. This paper abandoned such a possible link simply because of the short duration at DSCOVR. Note, therefore, that this paper gives one of the possible scenarios.
The notable difference is the weak-B region (continued for more than half a day from 1200 UT on February 3) in DSCOVR data, separating the two flux ropes. Unfortunately, the appearance of the weak-B region is far from the standard picture, and it would be impossible for forecasters to expect the very late arrival of the second flux rope at the Earth.
Our interpretation is that there were originally two flux ropes, as indicated in Figure 1h and illustrated in Figure 5, which appeared to be stacked together at the STEREO-A position but separated at the Earth. The detailed modeling to examine several other possibilities is beyond the scope of this paper, although it would provide excellent, challenging material for future advanced modeling work. In future work, for example, we need additional sophisticated analysis and high-resolution MHD simulation to give a definitive 3D picture of the flux ropes, including a possibility of local kink appearance led by the interaction between CMEs and the ambient solar wind structure.

Thermospheric response to the geomagnetic disturbance
The GAIA covers the atmospheric regions from the ground surface to the upper thermosphere and the ionosphere without artificial boundaries. The GAIA consists of three parts: a whole atmosphere General Circulation Model (GCM) part, an iono-spheric part, and an electrodynamics part. Although Jin et al. (2011), Miyoshi et al. (2018, and their references explained the details of GAIA, we briefly describe GAIA as follows. The whole atmosphere GCM part consists of a global spectral model with a maximum horizontal wave number of 42, corresponding to a grid spacing of approximately 2.8°longitude by 2.8°latitude with 150 vertical layers with a vertical resolution of 0.2 scale height. The ionospheric part solves the dynamics and chemical processes of the major ion species such as O + , O þ 2 ; N þ 2 , and NO + (Shinagawa, 2009), including the coupling processes between the plasma and the neutrals. The horizontal resolution of the ionospheric part is 2.5°longitude by 1°latitude, and the vertical resolution is 10 km. The electrodynamics part solves the global distribution of the ionospheric currents and electric fields at each time step in a tilted dipole coordinate under the condition of equipotentiality along the geomagnetic field lines (Jin et al., 2008).
In this paper, we use the real-time GAIA simulation (Tao et al., 2020), coupled with the empirical high-latitude electric potential model (Weimer, 1995) and Kp-dependent auroral total power (Zhang & Paxton, 2008). The solar wind parameters (By, Bz, V, and density) obtained by DSCOVR (NOAA Space Weather Prediction Center, 2016) were used as inputs for the Weimer model and the Kp estimation referring to Luo et al. (2017). For the real-time GAIA simulation, at low altitude (<~40 km) in the GCM part, realistic meteorological forcing is also included by assimilating the meteorological reanalysis data set, the Japanese 55-year Reanalysis (JRA-55), provided by the Japan Meteorological Agency (Kobayashi et al., 2015. R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 41 Harada et al., 2016. As a proxy of the solar EUV input to the model, we used the latest daily flux value of F10.7, 128.6~129.6 (10 À22 W m À2 Hz À1 ) during the simulation period from February 1 to February 4, 2022 (https://www.spaceweather.gc. ca/forecast-prevision/solar-solaire/solarflux/sx-5-en.php). Figure 6 shows the overall distribution of thermospheric density variation plotted for a reference quiet day (February 1, 2022), before and during the second minor storm at 0000 UT and at 2100 UT on February 4, 2022. Here we evaluated the quiet background density of the thermosphere from the GAIA simulation results on February 1, 2022. Figures 6 and 7 show significant enhancements of the thermospheric mass density (up to 50% at both the 200 and 400 km altitudes) in the vast area during the second magnetic storm. The up to 50% enhancements of the thermospheric mass density indicate the enhancements of the drag force to the satellites passing through these altitude regions, which is consistent with the documentation of SpaceX. We summarized the time development of the simulation results in Supplemental Information Movies A1-A4.
The 50% drag increase was not surprisingly higher than expected. For example, at 200 km altitude, a~25% increase in enhanced geomagnetic activity can be expected from the empirical model, NRLMSIS2.0 (Emmert et al., 2020), as shown by Supplemental Information Movies B1-B4. Note that the NRLMSIS2.0 depends on the ap index (Fig. 2). Also, the actual density variation can largely deviate from the empirical model, such as the so-called "cellular structure" (Crowley et al., 1996), which can easily explain the additional 25% variation. Compared with those previous studies, the present real-time GAIA simulation shows the thermospheric mass density enhancements in the wider area. Note also that the density perturbations of up to 50% are observed at 400 km in the NRLMSIS-2.0 (Supplemental Information Movies B3 and B4). The values of (storm-quiet)/quiet in % are shown for 00:00 UT (before the second magnetic storm) and 21:00 UT (during the second magnetic storm) on February 4, 2022, in the upper and bottom panels, respectively. The selected quiet day is February 1, 2022.
From the comparison of Supplemental Information Movies A1-A4 and B1-B4, however, it is also clear that the empirical model of the thermosphere such as NRLMSIS-2.0 is insufficient to describe the dynamic response of the thermosphere during magnetic storms, highlighting the importance of physical models over empirical models for future satellites orbitography.
The high-to-middle latitude propagation of the thermospheric mass density enhancement can be explained by looking at the other basic parameters (Fig. 8) as follows: The successive magnetic storms caused the enhancements of the Joule heating in the polar regions (Fig. 8a). Next, the Joule heating caused the expanding atmosphere, i.e., mass density enhancement at higher altitudes (Fig. 8b). Finally, the mass density enhancement regions propagated toward middle latitudes by the equatorward flow (Fig. 8d). Although it is beyond the scope of this paper, a more detailed examination of this particular event, comparing different thermospheric models and observations, is ongoing and will be reported elsewhere.

Summary and conclusions
We showed that the occurrence of moderate storms at unexpected timings on February 3 and 4, 2022, possibly caused the change in the orbits of the Starlink satellites. The crucial lessons learned from this GMD event are: (1) In real-time, it was difficult to expect the arrivals of two separated flux ropes from the possible overlapped appearance of CMEs. (2) It was also difficult to accurately predict the solar wind profiles by the flank-side passage of CMEs. (3) The real-time simulation of the thermosphere was necessary to predict the 50% mass density enhancement. We are entering a new age to quantitatively address the hard-to-predict and poorly-understood minor CME structures by utilizing multi-spacecraft observations in the upstream solar wind, combined with realistic CME and solar wind simulations. For future satellite operation safety, we would also need a better understanding of the possible errors and limitations of cutting-edge simulations of the thermosphere. R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 41 Acknowledgements. Data use SDO/AIA is through the courtesy of NASA/SDO and the AIA science teams. SOHO is a project of international cooperation between ESA and NASA. DSCOVR Magnetometer 1 min averages are obtained from DSCOVR Space Weather Data Portal (NOAA Space Weather Prediction Center, 2016; https://www.ngdc.noaa.gov/dscovr/ portal). The STEREO-A IMPACT beacon data is provided from NASA STEREO Science Center (https://stereo-ssc.nascom.nasa.gov/data/beacon/ahead/impact/2022/02/). The Dst and ap indices are provided from Kyoto University, via NASA/OMNI website (https://omniweb.gsfc.nasa.gov/ow. html). RK thanks Okinawa Institute of Science and Technology for hosting his sabbatical visit. DS is supported by JSPS KAKENHI 21H04492. YM, HJ, CT, HS, HF are supported by JSPS KAKENHI 21H01150, JP20K04037. HJ and HS are supported by 20K04037. CT, YM, and HF are supported by JSPS KAKENHI 19K03942. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.
Movies A1-A4: Results from real-time GAIA simulation. The selected altitudes are 200 km and 400 km, and the selected days are February 3 and 4.
Movies B1-B4: Results from NRLMSIS2.0 model. The selected altitudes are 200 km and 400 km, and the selected days are February 3 and 4. Observed ap index was used every three hours.