Study on the Thermospheric Density Distribution Pattern during Geomagnetic Activity

Abstract

The atmospheric density of the thermosphere is a fundamental parameter for spacecraft launch and orbit control. Under magnetic storm conditions, the thermospheric atmospheric density experiences significant fluctuations, which have a negative impact on spacecraft control. Exploring thermospheric density during geomagnetic storms can help to mitigate the effects of such events. Research on the inversion of accelerometer measurements for different satellites and the variations of atmospheric density under extreme conditions is still in its infancy. In this paper, the distribution of atmospheric density during three geomagnetic storms is investigated from the inversion results of the Swarm-C accelerometer. Three major geomagnetic storms and their recovery phases are selected as case studies. The thermospheric density obtained by Swarm-C is separated into day and night regions. The empirical orthogonal function analysis method is used to study the spatiotemporal distribution of thermospheric density during geomagnetic storms. The results indicate that storms have a more significant impact on nighttime thermospheric density. The impact of magnetic storms on the temporal distribution of thermospheric density is considerable. The first-order empirical orthogonal function (EOF) time coefficient value on the day after the storm is the largest, reaching 2–3 times that before the magnetic storm. The impact of magnetic storms on atmospheric density is mainly reflected in the time distribution. The spatial distribution of atmospheric density is less affected by magnetic storms and is relatively stable in the short term. The impact of magnetic storms on the spatial distribution of nighttime thermospheric density is more significant than that of daytime regions, and the response of daytime regions to magnetic storms is slower.

1. Introduction

Existing atmospheric models are primarily classified into physical models and empirical models [1,2]. Both models have their own advantages and application areas. Physical models attempt to describe the fundamental physical characteristics through equations [3]. These models include thermospheric general circulation models (TGCM) and local high-resolution flow models [4], which simulate temperature, density, atmospheric composition, horizontal wind, and atmospheric responses to different energy inputs. Compared with the empirical model, the physical model requires far more calculations, so the operator must have a very professional knowledge of physics [5]. Empirical models consist of several relatively simple functions that generally exhibit excellent performance in representing long-term and large-scale atmospheric characteristic changes. However, these functions contain parameters that are fitted using existing observational data, and hence require significant observational data to improve model accuracy. The prediction accuracy of existing atmospheric models, whether physical or empirical, depends on the forecast accuracy of solar and geomagnetic activity indices.

Studying thermospheric density distribution under extreme atmospheric conditions (geomagnetic storms) can significantly prevent economic losses and disasters. Sur D. et al. researched the effects of solar activity on the thermospheric density over the South Atlantic Anomaly (SAA) [6]. The researchers found that, during periods of high solar activity, the density in the SAA increased, while during periods of low solar activity, the density decreased. Amaechi, P.O. et al. analyzed irregular changes in the Earth’s ionosphere caused by the 2013 St. Patrick’s Day geomagnetic storm [7]. They found that during storms, irregular changes in the upper atmosphere’s density can significantly impact the operation of satellites. This article utilizes the inversion results of the Swarm-C accelerometer [8,9,10] and employs the empirical orthogonal function analysis method to study the spatiotemporal distribution of atmospheric density during geomagnetic storms.

Scientific studies have used satellite spaceborne accelerometer data to detect atmospheric density [11,12,13,14,15]. The most popular satellites for studying thermospheric density are GRACE and CHAMP [16,17]. PM Mehta et al. [18] used the CHAMP and GRACE satellites’ accelerometer data inversion to obtain the latest atmospheric density dataset. They found that the difference between the new and existing datasets is 14–18% compared with CHAMP, and 10–24% compared with GRACE. Emmert et al. [19] calculated the measured results of the thermal layer’s atmospheric density and the empirical model’s average error, with an error range of approximately 15–20%. This error would increase during solar activity. Pieter Visser et al. [20] use the TDW (thermospheric density and wind) and GHOST (GPS high precision orbit determination software tools) software provided by NASA to invert the GRACE and CHAMP satellites’ thermal atmosphere and horizontal wind field.

Thermospheric density is mainly influenced by solar and geomagnetic activity cycles. Therefore, the thermospheric density is correlated with solar and geomagnetic activities. However, the change in geomagnetic activity is delayed by 5 min to 10 h compared to solar activity [21,22,23].

In relevant studies on space climate, the value of solar radiation index F10.7 helps evaluate solar activity. When there is no strong geomagnetic disturbance, the period of change of atmospheric density in the thermal layer is 27 days, the same as the sun’s rotation period [24]. A significant positive correlation exists between the daily change in atmospheric density and solar radiation index F10.7 [25]. However, once a geomagnetic storm occurs, the thermospheric atmosphere heats up rapidly. The density of the thermosphere near the equator increases by 50 percent or more in the hours after the storm begins [26].

In addition to the solar radiation index, the geomagnetic index is usually used to represent the intensity of geomagnetic activity [27,28,29]. Geomagnetic activity indices are used to measure the level of disturbance in the Earth’s magnetic field caused by solar activity [30,31,32,33]. There are several different types of geomagnetic activity indices, including:

K-index: This index measures the maximum deviation of the Earth’s magnetic field in a 3 h period. A-index: This index is based on the average of the K-indices measured at a single station over a 24 h period. Ap-index: This index is similar to the A-index, but it is based on a global average of K-indices measured at several stations around the world. Dst-index: This index measures the strength of the ring current in the Earth’s magnetosphere. AE-index: This index measures the level of geomagnetic activity caused by high-energy particles from the solar wind. PC-index: This index measures the level of geomagnetic activity caused by the interaction of the solar wind with the Earth’s magnetic field. In this study, the Dst index and Ap index were selected as the research objects. The choice of these two factors is based on their ability to provide complementary information on the state of the Earth’s magnetic field. The Dst index provides a measure of the strength of the ring current, while the Ap index provides a measure of the overall level of geomagnetic activity.

Geomagnetic storms can be divided into initial, main, and recovery phases. The initial phase is characterized by the Dst component decreasing by −20 to −50 nT in a few minutes. The initial phase with a rapid decrease in Dst is also referred to as a storm sudden commencement (SSC). However, not all storms have an initial phase, and not all SSCs are followed by a storm. The Dst index’s decrease to less than −50 characterizes geomagnetic storms’ main phase stage, scientists have determined after years of observation. The Dst index’s minimum value will be between −50 and −600 during the main phase of the geomagnetic storm, which usually lasts 2–8 h. This change in the Dst index characterizes the recovery phase of a geomagnetic storm. This change is from its minimum value to a quiet period of −20 to 20. The recovery phase may last more than 8 h, sometimes even as long as seven days [34,35,36].

Since the 20th century, scientists have recorded many geomagnetic storms that have caused substantial economic loss and infrastructure malfunction [37]. For example, on 29 October 2003, a geomagnetic storm led to the wide-area augmentation system (WAAS), operated by the Federal Aviation Administration (FAA), being unfunctional for about 30 h [38]. In addition, Japan’s ADEOS-2 satellite was seriously damaged, and the communication industry of many other satellites was interrupted by the geomagnetic storm [39].

This study selected two geomagnetic indices, the Dst index and the Ap index, to analyze the relationship between geomagnetic indices and thermospheric density in order to explore the relationship between geomagnetic storms and thermospheric atmospheric density. Finally, three recent typical geomagnetic storms were selected as research objects. Empirical orthogonal function analysis was used to perform a spatiotemporal decomposition of the latitude–time distribution of thermospheric density during the storms. The first-order time coefficients and EOF [40] (empirical orthogonal function) base functions obtained from the decomposition were used to study the temporal and spatial distribution of thermospheric density during the storms, and the spatiotemporal variation patterns of thermospheric density during geomagnetic storms were summarized and generalized. The study reveals a complex interplay between geomagnetic storms and the density of the thermospheric atmosphere. By using two different geomagnetic indices, this study can provide a more comprehensive understanding of the relationship between geomagnetic storms and thermospheric density. Additionally, by exploring the relationship between geomagnetic storms and thermosphere atmospheric density, this study can have practical applications in space weather forecasting and satellite operations.

2. Date and Methodology

2.1. Data Selection

This study focuses on the Swarm-C satellite accelerometers as the research object. It utilizes the inversion results to discuss the spatial and temporal distribution of thermospheric density during magnetic storms. The Swarm series of satellites were launched on 22 November 2013. Swarm C flies at an altitude of about 460 km (with an inclination of 87.4°, an east–west separation of 1°–1.5° in longitude). The GPS receiver provides location- and time-related information for geomagnetic field measurements. It can determine non-conservative forces such as thermospheric drag, solar radiation pressure, and the recovery of thermospheric density [41]. The accelerometers can directly measure these forces with much higher resolution than the GPS receiver. Additionally, by modeling calculations, thermospheric density can be extracted from the Swarm satellite accelerometers. The latitude range of thermospheric density measured by the satellite is from 87.4° S to 87.4° N. Since the Swarm-C satellite was fully launched in 2014, this study selected magnetic storm events occurring after 2015. All data, including Swarm accelerometer data and Dst and Ap indices, were provided by the European Space Agency at https://earth.esa.int/eogateway/missions/swarm (accessed on 23 March 2023).

Geomagnetic storms are divided into the initial, main, and recovery phases. The primary phase, also known as storm sudden commencements (SSC), lasts tens of minutes,. The main phase of a geomagnetic storm usually lasts 2–8 h. The recovery phase may last more than 8 h, sometimes as long as 7 days. Scientists have successively recorded many geomagnetic storms. Three geomagnetic storm events were selected from https://www.swpc.noaa.gov/products/real-time-solar-wind (accessed on 23 March 2023). The main phases of these events were on 23 June 2015, 7 October 2015, and 8 September 2017.

We obtained the solar wind parameters during the above three magnetic storms from https://omniweb.gsfc.nasa.gov/ow_min.html (accessed on 23 March 2023). Based on the moment when the SYM-H index became negative, we determined the start times of the storms according to the starting time of the magnetic obstacle, which were 23 June 2015 02:23, 6 October 2015 21:35, and 8 September 2017 00:19, respectively.

Figure 1 shows the Dst index and Ap index of the three geomagnetic storms on the day before, on, and after the occurrence of the storm’s main phase. Figure 1 shows that the Dst index’s minimum values for selected storms were less than −100.

They reached the level of intense geomagnetic storms (Dst index between −100 nT and −200 nT). Among them, the minimum Dst index of the storm on 23 June 2015 was less than −200, which reached a great geomagnetic storm level (Dst index between −200 nT and −300 nT). The overall trend of the Ap index change was similar to that of the Dst index. However, Dst and Ap characterize different current systems and different types of geomagnetic disturbances. The vertical dashed lines indicate SSC and the start times of the main phase (MP) and the recovery phase (RP).

Figure 1 records the variation trends of the Ap index and Dst index during the three geomagnetic storms. The process from the initial phase to the main phase and then to the recovery phase can be seen in the three pictures. This is a rapid process from the quiescent period to the initial phase. In this process, the Dst index of the three geomagnetic storms decreased sharply in a short time, and the average duration of the process was less than one hour. In contrast, the main phase of the geomagnetic storm lasted for a long time, and the Dst index did not always decrease during the main phase but fluctuated. On 8 September 2017, when the Dst index reached the minimum value during the main phase, it gradually increased. However, it declined from the 36 h until the 39 h, when there was another peak. Since the Swarm-C satellite was launched in 2014, this paper selects the geomagnetic storm events after 2015. The relevant data for the three geomagnetic storms will help study the relationship between geomagnetic storms and atmospheric density.

2.2. Swarm-C Accelerometer Data Calibration

In the Swarm-C satellite data obtained from the ESA, two types of problems need to be handled manually: the spike caused by the propulsion system in the accelerometer data and the abnormal data caused by error detection and correction (EDAC) events. During the experiment, we deleted the spike caused by the propulsion system in the accelerator data. Then, we filled the data gap by linear interpolation. For EDAC, we performed manual step correction on the step signal in the data.

Besides the manual correction, the calibration method used in this paper is based on the GRACE-A satellite accelerometer’s calibration method [42,43]. However, unlike the GRACE-A, the Swarm-C accelerometer corrects the data error caused by the accelerometer carrier’s temperature change. Therefore, a linear temperature correction will be introduced in this paper for this difference.

For the point of the satellite orbit, the calibration equation is defined as shown in Equation (1). In Equation, Q represents a linear temperature factor. $T(t+F)$ represents a temperature signal with a time shift of F. Parameter $\epsilon$ represents statistical noise and $a_{ACC}^{UNCAL}$ represents the uncalibrated acceleration data.

$$a_{NG}^{GPS}=B+S·a_{ACC}^{UNCAL}+Q·T\left(t+F\right)+\epsilon$$

(1)

This paper uses GPS-derived non-conservative force acceleration data as the calibration standard. It is projected in the direction of the calibrated accelerometer data component to obtain $a_{NG}^{GPS}$, the projected GPS-calibrated acceleration. The calibration parameters $\widehat{B}$, $\widehat{S}$, $\widehat{Q}$ are deduced by the least-squares method combined with Equation (1), where a variable cap with a variable represents the estimate of the linear calibration parameter. After obtaining the variable’s estimated value, the calibrated acceleration value $a_{ACC}^{CAL}$ can be calculated by Equation (2).

$$a_{ACC}^{CAL}=\widehat{B}+\widehat{S}·a_{ACC}^{UNCAL}+\widehat{Q}·T(t+F)$$

(2)

The non-conservative acceleration signal is the strongest in the satellite’s orbit direction for the Swarm-C satellite. The acceleration caused by the atmospheric resistance is dominant. Although the Swarm-C satellite accelerometer measures the acceleration of the three axes of the satellite reference system, this paper only calibrates and discusses the acceleration data along the satellite orbit direction. The acceleration data in other directions can be computed in the same way.

2.3. Decomposition of Atmospheric Density Based on EOF

This study refers to the Sentman model [44] for thermospheric density inversion based on accelerometer data. We also refer to the Swarm-C data provided by http://thermosphere.tudelft.nl/ (accessed on 23 March 2023).

The eigenvector analysis method is also known as the empirical orthogonal function analysis (EOF) method [45]. The empirical orthogonal function analysis method can decompose the meteorological variable field into two parts: function with time as the variable and function with the spatial position [46]. The variance of the first few components of the principal component decomposed by the empirical orthogonal function usually accounts for a large percentage of the total variance. Hence, the variation of the principal components can describe the variation of all components with time [47]. When the empirical orthogonal analysis method is applied to analyze the temporal and spatial distribution of thermospheric density in the thermosphere, the thermospheric density extracted by the Swarm-C satellite accelerometer should be smoothed to obtain the scatter points, whose abscissa is time and ordinate is latitude. Then, the scatter set X with the same time interval and latitude interval are selected. Each datum is represented by x(t, φ). According to the basic principle of empirical orthogonal decomposition, X can be decomposed into function V with time as a variable and function Z with space as a variable, i.e., X = VZ. The decomposition is shown in Equation (3).

$$x\left(t_i,{\phi }_i\right)={\sum }_{k=1}^P{v}^k\left({\phi }_i\right)z^k\left(t_j\right)$$

(3)

where $t_j(j=1,2,\cdots ,n)$ denotes the time corresponding to each lattice point, ${\phi }_i(i=1,2,\cdots ,m)$ is the latitude corresponding to each grid point. $v^k$ represents the k-th mode function, i.e., the variation of atmospheric density with latitude. $z^k$ is the k-th time coefficient, representing the change of atmospheric density with time. The results obtained by EOF analysis are evaluated in this paper using the storm data from 8 September 2017, the latest from the selected storms.

2.4. Correlation between Dst, Ap, and Atmospheric Density

The periodicity of solar activity determines the basic structure of the thermosphere and the regular long-term changes. In contrast, the influence of geomagnetic activity on thermospheric density is short-term and severe [48]. Therefore, the thermospheric density, Dst index, and Ap index obtained from the Swarm-C satellite accelerometer are compared to study the relationship between the intensity of geomagnetic activity and thermospheric density.

A point 72 h before the beginning of the geomagnetic storm is taken as the abscissa’s initial value. The density of the thermosphere and Dst is used as the ordinates to draw a biaxial map. The Ap index is recorded every 3 h. The atmospheric density of the thermosphere and Ap index are used as the ordinates to draw a biaxial map. The final results are shown in Figure 2.

Figure 2a,c,e show the relationship between thermospheric density and the Dst index. These three pictures show that the amplitude of thermospheric density and the magnitude of geomagnetic index Dst negatively correlate during geomagnetic storms. The thermospheric density value is in the trough when the Dst index peaks in the three geomagnetic storms. The peaks and valleys of the two curves are opposite. The correlation coefficient between the Dst index and thermospheric density value is −0.805 in the geomagnetic storm on 23 June 2015, −0.87 in the geomagnetic storm on 7 October 2015, and −0.795 in the geomagnetic storm on 8 September 2017. The absolute value of the correlation coefficient between the Dst index and thermospheric density is more than 0.8 in two storms, which has reached a negative significant correlation level.

Although the peak valley of thermospheric density is opposite to that of the Dst index, the maximum absolute value is not necessarily in the same position. For example, on 8 September 2017, peaks appeared at the 74 h and 90 h. However, the maximum absolute value of the Dst index appeared at the 74 h, and the maximum absolute value of thermospheric density appeared at the 90 h. Moreover, the peak value of the Dst index is not synchronous with the peak of thermospheric density; instead, there will be a specific time delay. Among the three geomagnetic storm events, the time delay between the Dst index peak and the peak of thermospheric density on 22 June 2015 is the largest.

Figure 2b,d,f show the relationship between thermospheric density and the Ap index. It can be seen from the three pictures that the Ap index shows a positive correlation with thermospheric density. The value of thermospheric density rises and falls with the rise and fall of the Ap index. The difference is that the maximum value of the Ap index appeared earlier than the maximum value of thermospheric density in the three geomagnetic storms. In the geomagnetic storm on 23 June 2015, the maximum value of the Ap index appeared about 10 h earlier than the maximum value of thermospheric density. In the geomagnetic storm on 7 October 2015, the maximum value of the Ap index appeared earlier than the maximum value of thermospheric density. In the geomagnetic storm on 8 September 2017, the maximum value of the Ap index appeared about 4 h earlier than the maximum value of thermospheric density. Moreover, unlike the thermospheric density curve showing a stable upward or downward trend, the curve fluctuation of the Ap index is more evident than that of the thermospheric density curve.

The Dst index lags behind thermospheric density, while the Ap index is ahead of thermospheric density. The different cadences may cause this difference. The change in thermospheric density and the Dst index always lags behind the Ap index for a while, so the Ap index is a more sensitive data choice when studying the change in thermospheric density. However, the disadvantage of the Ap index is that the correlation between the Ap index and thermospheric density is not as high as that between the Dst index and thermospheric density. Therefore, the correlation level between the Dst index and thermospheric density significantly negatively correlates most of the time. In contrast, the correlation between the Ap index and thermospheric density only reaches a highly positive correlation.

3. Results

3.1. Comparison of EOF Density Fields of Different Orders

The three-day thermospheric density data during storm events are selected to carry out empirical orthogonal decomposition to study the spatiotemporal decomposition of thermospheric density during a geomagnetic storm.

The specific steps for empirical orthogonal decomposition of thermospheric atmospheric density data are as follows:

Firstly, the measured thermospheric density field obtained by inverting the accelerometer data for the above three days is processed by subtracting the mean value, with a latitude interval of 5 degrees and a time interval of 1 h. This produces a matrix of anomalies in the measured thermospheric density field.

Next, using Equation (3), the matrix is multiplied by its transpose matrix to obtain a real symmetric matrix.

Finally, the eigenvalues and eigenvectors of the real symmetric matrix are calculated to obtain the characteristic vector fields of the accelerometer-measured thermospheric density variable field for each order.

Figure 3 shows the temporal and latitudinal distribution of thermospheric density in the thermosphere drawn by the characteristic vector fields of the thermospheric density of different orders obtained by empirical orthogonal decomposition.

As shown in Figure 3, on 8 September 2017, the atmospheric density increased due to the geomagnetic storm. As a result, the middle of the original measured atmospheric density field showed bright yellow. The distribution of the first-order atmospheric density eigenvector field obtained by EOF decomposition is similar to the original measured atmospheric density field distribution, but the details are different. The high-density area in the middle of the first-order EOF atmospheric density feature vector field is significantly larger than the original measured atmospheric density field. As a result, the color of the high-density area (orange–yellow) is brighter and more concentrated. Consequently, the color contrast with the surrounding area is more pronounced. Compared with the original measured atmospheric density field, the density range of the first-order EOF atmospheric density eigenvector field is smaller than that of the original measured atmospheric density field. After the second-order EOF, the atmospheric density eigenvector field is superimposed. As a result, the first two EOF atmospheric density eigenvector fields’ distribution is closer to the original measured atmospheric density field.

As a result, the continuous range of orange–yellow high-density areas in the figure is expanded. After superimposing the third-order EOF atmospheric density eigenvector field, the first three EOF atmospheric density eigenvector fields have almost identical characteristics to the original measured atmospheric density field.

The first-order thermospheric atmospheric density characteristic vector field represents significant large-scale changes in the measured atmospheric density field, with relatively less detailed information. As higher-order characteristic vector fields are gradually added, the short-term and local features of atmospheric density become apparent. The more orders of characteristic vector fields that are included, the more similar their distribution is to the original measured atmospheric density distribution, and the more comprehensive the description of the measured atmospheric density field becomes. The high-order characteristic vector fields often reveal features with small amplitudes and randomness, while the low-order vector fields typically reveal significant features with large amplitudes and long durations, whose physical meaning is clearer.

This paper explores statistics on the variance contribution rate of the first four order eigenvectors obtained after the EOF decomposition of atmospheric density on 8 September 2017.

Table 1. First-order EOF eigenvectors and their contribution statistics.

Order Number	The Eigenvalue (×10−20)	Variance Contribution Rate	Cumulative Variance Contribution Rate
First order	0.2408	0.9804	0.9804
Second order	0.0033	0.0135	0.9939
Third order	0.0010	0.0042	0.9980
Fourth order	0.0003	0.0014	0.9994

shows the eigenvectors obtained after the first four-order EOF decomposition and their statistical contribution results.

From the statistical data in Table 1, it can be seen that the variance contribution rate of the first-order EOF atmospheric density characteristic vector field reached 98.04%, while the variance contribution rate of the second-order EOF atmospheric density characteristic vector field dropped significantly to 1.35%, and from the second order onwards, the variance contribution rate of each order gradually decreases. In terms of cumulative variance contribution rate, the variance contribution rate of the first two EOF atmospheric density characteristic vector fields reached 99.39%, and the variance contribution rate of the first four orders reached 99.94%, which is already very similar to the original measured atmospheric density distribution field. Because the variance contribution rate of the first-order atmospheric density characteristic vector field has already exceeded 95%, it can be considered that the first-order atmospheric density characteristic vector field is sufficient to contain most of the information of the original measured atmospheric density field. Therefore, we assume the first-order EOF atmospheric density characteristic vector field is sufficient to highly fit and restore the measured atmospheric density distribution field. In the subsequent research, in order to filter out some incidental or even white noise features in the measured atmospheric density vector field and to avoid the contingency of research conclusions, the first-order EOF atmospheric density characteristic vector field is adopted as the research object in the data processing process to explore the long-term and large-scale changes of atmospheric density.

After the EOF decomposition is obtained, it is necessary to test the data obtained. Because the fitting degree of each order of the eigenvector field to the original measured density field has been analyzed in previous works, there is no composite analysis test here. Instead, the North test examines whether the vector fields of different orders are independent; that is to say, whether the vector fields of each order can be regarded as vector fields with independent characteristics. Sometimes the research data may be mixed with random or false data. After empirical orthogonal decomposition, these invalid and real data can be decomposed into several principal components and eigenvectors. Therefore, the North significance test is necessary.

In the North significance test, each order vector field’s eigenvalues need to be substituted into the function to obtain each eigenvalue’s lower and upper error limits. For example, suppose that the error of two adjacent eigenvalues intersects the lower and upper limits. In that case, the eigenvector field corresponding to the eigenvalue does not pass the significance test. Therefore, the vector field corresponding to the eigenvalue is invalid data and must be eliminated.

Table 2. Error upper and lower limits of the first five order eigenvector fields of EOF.

Order Number	The Lower Limit of Error	The Upper Limit of Error
1	1.84071 × 10−21	2.97603 × 10−21
2	2.52995 × 10−23	4.09037 × 10−23
3	7.81714 × 10−24	1.26386 × 10−23
4	2.57967 × 10−24	4.17076 × 10−24
5	6.23159 × 10−25	1.00751 × 10−24

shows the corresponding upper and lower error limits for the first five-order eigenvector fields obtained by EOF decomposition after the North test.

Table 2 shows that the error range of the characteristic root corresponding to the first-order atmospheric density eigenvector field is [1.84 × 10⁻²¹, 2.976 × 10⁻²¹]. The second-order’s error range of the characteristic root is [2.53 × 10⁻²³, 4.09 × 10⁻²³]. The error ranges of the third, fourth, and fifth orders are [1.264 × 10⁻²³, 7.817 × 10⁻²⁴], [4.170 × 10⁻²⁴, 2.580 × 10⁻²⁴], and [1.008 × 10⁻²⁴, 6.232 × 10⁻²⁵], respectively.

The error upper limit of the characteristic root of the second-order atmospheric density eigenvector field is less than that of the first-order atmospheric density eigenvector field. There is no intersection between the measurement field’s error ranges. Therefore, the first-order atmospheric density eigenvector field passes the modal significance test of EOF.

Similarly, according to the above steps, the EOF modal significance tests for the corresponding eigenvalues of the second, third, and fourth-order atmospheric density eigenvector fields show no intersection between the upper and lower limits of the error. Therefore, atmospheric density’s first, second, third, and fourth-order eigenvector fields are valuable signals and vector fields with independent characteristics. It should be noted that the eigenvector fields that do not pass the significance test are generally high-order eigenvector fields. Suppose a situation does not pass the significance test in the follow-up study. In that case, special instructions will be given when the elimination operation is needed.

3.2. Analysis of the Temporal Distribution of Atmospheric Density during Storms

We aim to analyze the temporal distribution of the thermospheric atmosphere density before and after a geomagnetic storm. Therefore, we obtained the measured atmospheric density data of the three geomagnetic storms the day before the storm, the day of the storm, and the day after.

The EOF daily spatiotemporal decomposition was performed for each data group separately, and only the first-order time coefficients after the decomposition were used. Figure 4 shows the first-order time coefficients of the day before, on, and after the three geomagnetic storms. The figures are furtherly separated into daytime and nighttime sub-figures for better presentation.

Data processing shows that the atmospheric density distribution in the daytime area differs from that in the nighttime because of radiation heating. Therefore, in the following research process, the atmospheric density in the daytime and the atmospheric density in the nighttime are processed separately. That is to say, the daily data of atmospheric density will be divided into two situations: daytime and nighttime.

The comparison between daytime and nighttime in Figure 4 shows the difference in the time coefficient. The EOF time coefficients of daytime and nighttime are very similar, but amplitude differences exist. In the daytime coefficients, the time coefficients of the day before the geomagnetic storm and the day after the geomagnetic storm are greater than those of the day before the geomagnetic storm, the day after the geomagnetic storm, and the day after the geomagnetic storm in the corresponding nighttime. In other words, when the thermosphere is in the daytime, its density is significantly higher than in the dark region. This phenomenon is true both before and after the geomagnetic storm.

The comparison in Figure 4 also shows the time coefficient’s difference before and after the geomagnetic storm. The time coefficient curve of the day before the occurrence of the geomagnetic storm (blue), the time coefficient curve of the day of the occurrence of the geomagnetic storm (red), and the time coefficient curve of the day after the occurrence of the geomagnetic storm (yellow) have obvious differences. The day’s curve is the most stable before the geomagnetic storm’s occurrence. The time coefficient of the day before the geomagnetic storm is similar to that during the occurrence. The difference is that the time coefficient value on the day of the occurrence peaks at the end. The peak value of the geomagnetic storm’s red curve on 23 June 2015 is the steepest and most apparent. The curve after the geomagnetic storm is at the top of the three geomagnetic storm events. That is to say, the value of the first-order time coefficient of the day after the geomagnetic storm is the largest. The first-order time coefficient’s peak value after the geomagnetic storm can reach 2–3 times that before the geomagnetic storm. All of them will fall back significantly. Therefore, the geomagnetic storm’s influence on the distribution of atmospheric density on the time scale is considerable. This effect will not be reflected immediately on the same day but the day’s coefficient after the geomagnetic storm.

3.3. Analysis of the Spatial Distribution of Atmospheric Density during Geomagnetic Storms

Similarly, using empirical orthogonal decomposition, Figure 5 shows the first-order EOF spatial distribution curves of the day before, on, and after the three geomagnetic storms.

From the longitudinal comparison in Figure 5, for the spatial distribution of thermospheric density in the daytime area, the spatial distribution characteristics of the thermospheric density in the daytime region of the three geomagnetic storms are similar. Their curves all present a bimodal structure. In other words, there are two atmospheric density peaks at 40° to 60° north–south latitude, and the minimum thermospheric density appears in two polar regions. The minimum values of the geomagnetic storms on 23 June 2015, and 7 October 2015, appeared near the south pole. The minimum values of the geomagnetic storm on 8 September 2017 appeared near the North Pole. Although there are peaks in both the northern and southern hemispheres, the two peaks are not symmetrical. The peak values of the southern hemisphere were higher during the geomagnetic storms of 23 June 2015 and 7 October 2015, and the peak values of the northern hemisphere were higher on 8 September 2017.

Another similarity is a trough near the equator for the spatial distribution of thermospheric density in the daytime area. Although the equatorial region’s atmospheric density is not the minimum value, it is slightly lower than nearby regions. Compared with the daytime region’s time coefficient curve, the atmospheric density’s spatial distribution curves before and after the geomagnetic storm are similar. There will not be a sudden change in one day, like the time coefficient, but instead there will be a gradual change.

There are some similarities in the spatial distribution of atmospheric density in the three geomagnetic storms’ thermospheres for the spatial distribution of atmospheric density in the nighttime. The curves show a bowl-like structure, and a depression appears in the middle of all the curves. This pattern means that the lowest value appears near the equator. Comparing the spatial distribution curve of atmospheric density on the day before the geomagnetic storm with the day after the storm, the curve before the storm is the gentlest, and the fluctuation of the atmosphere is the lowest. With the coming of the geomagnetic storm date, two asymmetric peaks of atmospheric density gradually appeared in the mid-latitude region of the northern and southern hemispheres.

From the horizontal comparison in Figure 5, the spatial distribution of the atmospheric density in the daytime and nighttime has a common point in the same geomagnetic storm. This is because a trough appears near the equator, a common feature of the three geomagnetic storms. The difference is that the spatial distribution curve of atmospheric density in the nighttime changes more obviously than in the daytime before and after the onset of the geomagnetic storm. That is to say, the influence of geomagnetic storms on the spatial distribution of atmospheric density in the night region is more significant than in the daytime region. As a result, daytime regions are slower to respond to geomagnetic storms.

4. Discussion

Based on the results obtained from the EOF decomposition experiment of atmospheric density during a geomagnetic storm, as well as the analysis of the spatial-temporal distribution of atmospheric density during a geomagnetic storm, we have reached the following conclusions:

(1) The first-order EOF characteristic vector field represents the most significant and large-scale changes in the original atmospheric density field while containing less detailed information. The addition of higher-order eigenvector fields gradually reveals the short-term and local characteristics of the thermosphere density and the more comprehensive description of the original atmospheric density distribution characteristics. The higher-order eigenvector field often reveals smaller amplitude characteristics, while the lower-order vector field reveals larger amplitude and longer duration features.

(2) The analysis of atmospheric density time distribution during a geomagnetic storm showed that geomagnetic storms significantly influence the temporal distribution of atmospheric density in the thermosphere. The time coefficient curve of the day after the geomagnetic storm differs from that of the other two days. In contrast, the influence of geomagnetic storms on the spatial distribution of atmospheric density in the thermosphere is smaller than that on the temporal distribution. The spatial distribution of atmospheric density in the thermosphere is generally stable during a certain period, especially without a geomagnetic storm. Even if a geomagnetic storm occurs, the disturbance to the spatial distribution of atmospheric density is relatively small and cannot fundamentally change its state. Although the nighttime response to a geomagnetic storm is obvious, the continuous rotation of the Earth results in alternating day and night, thus making the geomagnetic storm’s influence on atmospheric density’s spatial distribution relatively small. Consequently, the spatial distribution stability of the Earth’s thermosphere density is high, with slow short-term changes that accumulate over time. For the three geomagnetic storms, the spatial distribution curve of atmospheric density differs for each storm, resulting from long-term cumulative changes.

(3) In the analysis of the spatial distribution of atmospheric density during a geomagnetic storm, the spatial distribution of atmospheric density in the daytime and nighttime regions is asymmetric in the northern and southern hemispheres before and after the storm. This asymmetry is related to the intensity of the geomagnetic storm itself, the area of solar wind energy injection, and the penetration height of energy particles. When the geomagnetic storm begins, the geomagnetic index rises sharply, and the disturbance increases. The spatial distribution of EOF atmospheric density reflects the heat conduction process of the atmosphere.

5. Conclusions

The thermospheric density distribution during the geomagnetic storm was analyzed in this paper based on the Swarm-C accelerometer’s inversion results. First, three recent geomagnetic storms with more significant impacts were used as the research objects. Two geomagnetic indexes, Dst and Ap were selected to analyze the relationship between the geomagnetic index and the thermospheric density during the geomagnetic storm. After calculating the correlation between the geomagnetic index and the thermospheric density, it was found that the Ap index is advanced and can play a predictive role. Secondly, the measured atmospheric density data were separated according to the day and night regions. The empirical orthogonal analysis method decomposed the time distribution of the thermospheric density latitude and time during the three-field geomagnetic storm. The first-order time coefficient and the first-order EOF basis function were obtained to study the geomagnetic storm. Finally, we summarized temporal and spatial changes in atmospheric density during geomagnetic storms. We concluded that geomagnetic storms have a more significant impact on the temporal distribution of atmospheric density than on the spatial distribution. The spatial distribution of atmospheric density is relatively stable and changes slowly. Since this paper only takes three significant storms since 2015, s study using more storms would give a more robust result on the patterns of the thermosphere. There are also many factors related to the atmospheric density that this paper does not consider, such as thermospheric wind. For future studies, it is possible to consider more factors and use storms of different magnitudes to understand the thermosphere and geomagnetic activities better.