Keywords

1 Introduction

The connection between solar drivers, Earth’s magnetosphere, and Ionosphere-Thermosphere (IT) phenomena in the upper atmosphere is very complex and dependent on many processes, including energy-absorption, ionization, and dissociation of molecules due to variable X-ray and Extreme Ultra Violet (EUV) solar radiance (Calabia and Jin 2019). In addition, the variable solar wind plasma, combined with a favorable alignment of the Interplanetary Magnetic Field (IMF), can produce aurora particle precipitation at high latitudes, which results in chemical reactions and Joule heating through collisions between electrically-charged and neutral particles.

Consequences of upper-atmospheric conditions on human activity highlight the necessity to better understand and predict IT processes, potentially preventing detrimental effects to orbiting, aerial, and ground-based technologies. Charged particles (mostly free electrons in the ionosphere) are able to influence the propagation of electromagnetic radio waves in Global Navigation Satellite Systems (GNSS), satellite communication infrastructures, and Earth observation techniques. For instance, the ionospheric effects on GNSS result in a time-delay of the modulated components from which pseudo-range measurements are made, and a time-advance in carrier-phase measurements (Jin and Su 2020; Jin et al. 2019, 2020). These effects are directly related to the final accuracy and reliability of GNSS onboard satellites, aviation, and ground vehicles, providing unacceptable errors for single-frequency receivers (Hofmann-Wellenhof et al. 2008).

Existing ionospheric total electron content (TEC) models are designed to give a general description of the variations in the ionosphere. These include, for example, the Bent ionospheric model (Bent and Llewellyn 1973), the Parameterized Ionospheric Model (PIM) (Daniell et al. 1995), the NeQuick model (Radicella 2009), and the International Reference Ionosphere (IRI) model (Bilitza et al. 2011). However, the highly variable ionosphere and the complex processes involved in the IT coupling result in difficulties for accurate modeling, and the existing models only provide monthly averages of behavior for especially magnetically quiet conditions. For instance, IRI and NeQuick are only climatological models, and cannot accurately represent space weather events such as geomagnetic storms.

The variable aerodynamic-drag on LEO satellites is associated with the upper-atmospheric expansion/contraction in response to variable solar and geomagnetic activity. The drag makes tracking difficult, decelerates LEO orbits, reduces their altitude, and shortens the lifespan of satellite missions. In addition, the exponential increase of space debris (including the recent destructive events of Fengyun-1C, Iridium, and Mission Shakti) has recently highlighted the importance of orbital tracking and prediction of potential collisions with orbiting satellites. Unfortunately, the existing aerodynamic-drag models used in Precise Orbit Determination (POD) result in positioning errors, far away from the operational requirements (Anderson et al. 2009; Calabia et al. 2020), and several recent studies have exposed limitations to accurately predict the actual neutral density variability (e.g., Müller et al. 2009; Sutton et al. 2009; Doornbos et al. 2010; Emmert and Picone 2010; Liu et al. 2011; Lei et al. 2012; Chen et al. 2014; Cnossen and Förster 2016; Calabia and Jin 2016; Panzetta et al. 2018).

The first empirical models based on orbital decay were introduced by Harris and Priester (1962) and Jacchia (1964). Upper-atmospheric total mass density (TMD) models have since been improving with the use of new techniques, algorithms, and proxies, as for example, Bowman et al. (2008) with the JB2008, Bruinsma (2015) with the Drag Temperature Model (DTM), and Picone et al. (2002) with the US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter-radar Extended 2000 model (NRLMSISE-00).

The exact connection between the ionospheric TEC and thermospheric TMD models is still unclear, but physics-based models such as the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIEGCM) (Richmond et al. 1992), which generally use empirical parameterizations and boundary conditions to solve 3-dimensional fluid equations, can provide an approximated solution with various prognostic variables including, but not limited to, TEC, TMD, Joule heating, etc. Compared to empirical models, the physics-based models are far more complex, but can help to better understand the physical mechanisms responsible for the observed IT variability, and have a greater potential for predictions and projections of future occurrences.

Unfortunately, the existing upper-atmospheric models are incapable of predicting the IT variability as required, in spite of the efforts to model variations, anomalies, and climatology over the last half-century. This is largely due to the limited quality and quantity of observations used to better characterize the driver-response relationship of the IT variability, and the lack of comprehensive approaches for calibrating the models. In response to this situation, the international community sought to increase scientific research on upper-atmosphere modeling, develop safeguard schemes, and produce space weather products, centers, and services. As a result, on 22 April 2017, in Vienna, the Focus Area on Geodetic Space Weather Research (FA GSWR) was created within the structure of the Global Geodetic Observing System (GGOS) of the International Association of Geodesy (IAG). Since then, the GGOS FA GSWR has been defining its internal structure and main objectives with the purpose to collate existing initiatives and guide future activity in this area. Since many geodetic tasks depend on upper-atmospheric properties, the development of ionosphere and neutral density models as GGOS products for direct applications has become the main objective. The committee of the FA GSWR selected the electron and the neutral density as potential candidates for the list of essential geodetic variables (EGV).

In this paper, we aim at the first objective with an innovative study based on the Principal Component Analysis (PCA) of large data-sets of globally sparse observations. In this scheme, the main PCA modes derived from different physical variables will be investigated, modeled, and compared in a common space-time frame and the correlations, similarities, and differences will provide evidence regarding possible coupling and driving mechanisms. Previous studies using the PCA technique have mapped the data to an orthogonal basis composed of a number of spherical harmonics (SH) functions expressed in the Local Solar Time (LST) and magnetic dip latitude coordinate system (Matsuo and Forbes 2010; Wan et al. 2012). In our method, instead of employing a basis composed of a number of SH functions, we benefit from the full resolution provided by the initial variables and use of the geographical latitude and longitude coordinate system. This coordinate system will reveal possible geographical contributions, including the South Atlantic Anomaly (SAA), or the South magnetic dip, as well as simplify the modeling of the initial variables for practical applications. Other alternative techniques widely employed in spatiotemporal analyses include the use of neural networks (e.g., Gowtam et al. 2019) or wavelet analyses (e.g., Dabbakuti and Ratnam 2016).

Here we present the first results from comparing the main PCA components of TEC and TMD observations from a complete full solar cycle. Section 2 briefly introduces the data and methods employed in our research; PCA results are illustrated and compared in the third section; and finally, conclusions are given in the last section.

2 Data and Methods

2.1 Ionospheric Observations

The 16 year TEC time-series from IGS global ionosphere maps (GIMs) (2003–2018) have been downloaded from the Crustal Dynamics Data Information System (CDDIS) of National Aeronautics and Space Administration (NASA) (https://cddis.nasa.gov/index.html). The data is the result of an integrated set of GNSS dual-frequency tracking data from more than 350 permanent stations. The temporal resolution of each TEC GIM is provided at 2 h LST, and the spatial resolution is 2.5° by 5° in latitude by longitude, respectively.

2.2 Thermospheric Observations

We employ 13 years’ (2003–2015) of TMD time-series inferred from accelerometer and POD measurements made by the GRACE (Gravity Recovery and Climate Experiment) mission (Tapley et al. 2004). TMD estimates were obtained from the Information System and Data Center (ISDC) GeoForschungsZentrum (GFZ), computed by Calabia and Jin (2016, 2017), and provided at 3-min interval sampling in Calabia and Jin (2019). In order to deal with the different orbit locations and altitudes, the data have been normalized to 475 km altitude and interpolated between orbits. More details on the data processing can be found in Calabia and Jin (2016).

2.3 Space Weather and Geomagnetic Indices

Space weather and geomagnetic indices have been downloaded from the Low-Resolution OMNI (LRO) data set of NASA (http://omniweb.gsfc.nasa.gov/form/dx1.html) and from the International Service of Geomagnetic Indices (ISGI) website (http://isgi.unistra.fr/data_download.php). NASA/ESA (European Space Agency) Solar and Heliospheric Observatory (SOHO) satellite Solar Extreme-ultraviolet Monitor (SEM) measurements are provided by Bowman et al. (2008).

2.4 PCA Modeling

The aim of the PCA technique is to determine a new set of basis vectors that capture the most dominant structures in the data, based on eigenvalue decomposition of the covariance matrix. Detailed analyses and the selection of retained modes can be found in Preisendorfer (1988) and Wilks (1995), and a readily computable algorithm in Bjornsson and Venegas (1997). In our modeling scheme, given a time-series of estimates, the spatial patterns of its variability, temporal variation, and the measure of its importance are presented as a low-dimensional space spanned by a set of modes, which can be parameterized in terms of most representative proxies (e.g., solar flux index, annual cycles, etc.). The measure of importance of each component is provided by the analysis itself, and given in % to the total variance. This is calculated from the eigenvalues (contribution of each mode to the total variability). Our general method (Calabia and Jin 2016) is composed in three basic steps, which include (1) obtaining grids for given time moments; (2) arranging each grid in by columns and finding the eigenvectors (space-dependent components), eigenvalues (contribution of each mode to the total variability), and the projections of the initial matrix on each eigenvector (time-dependent components); and (3) analytical fitting of resulting modes in terms of most representative proxies, including a statistical assessment (i.e., standard deviation, Pearson correlation).

Fig. 1
figure 1

First spatial PCA component of (a) TEC and (b) TMD at 475 km altitude. All longitudes provide values at noon (12 h LST). Dimensionless quantities. Dip isoclinic lines are plotted in dash-dotted lines

Full size image

3 Results and Analysis

The results from both TEC and TMD analyses show that the dominant forcing is the solar flux cycle. Figures 1 and 2 show the space and time-dependent components of the first PCA mode, respectively. Secondary and above PCA modes are not included in this manuscript. Concerning the analysis of TMD, the first PCA mode accounts for 92% of the total variability, and its parameterization reaches 96% of the correlation. As for the TEC variability, the first leading mode accounts for 75% of the total variance at each LST case (12 PCA analyses, one each 2 h sampling), and its parameterization reaches 98% of correlation. For both cases, the high values of explained variance indicate marked patterns of variability. The high correlations in the fitting of the PCA modes indicate high accuracy to represent the actual TMD and TEC variability.

Each spatial PCA component is interpreted as a set of values in space that vary depending on a coefficient that varies on time (temporal PCA component). High PCA coefficient indicates high variability. Figure 1 shows the spatial components at noontime of the first PCA mode from both the 16 years’ time-series of TEC and the 13 years’ time-series of TMD (475 km altitude). The corresponding time-expansion coefficients (Fig. 2) are mainly related to the solar flux cycle (analysis presented in the next paragraph). In Fig. 1, we can observe that both values of TEC and TMD show a drop off along the magnetic equator, depicting the characteristic equatorial bulge with a two-crest shape, at about ±20° dip latitude for TEC (namely equatorial ionization anomaly or EIA), and about ±30° dip latitude for TMD (namely equatorial mass anomaly or EMA) (Liu et al. 2009). Minimum values along the magnetic equator are similarly located for both cases at 0°E and at about 90°E, and maximum values in the crests are located at the West, over the eastern Pacific Sea. However, distinct features show no clear spatial correlation between the TEC and TMD and suggest they originate from different physical mechanisms. For instance, a clear asymmetry is shown only in the global distribution of TMD, with higher values in the southern hemisphere than that shown in the northern hemisphere. The most characteristic discrepancy is the marked enhancement near the southern magnetic pole (dip = 90°). None of these structures are reflected by the ionospheric TEC distribution. In addition, we can observe wider and less marked equatorial crests of TMD when comparing to the TEC crests. The maximum peak of the northern TMD crest is displaced about 40° East from the maximum peak of the TEC crests. The northern TMD crest in the eastern longitudes is very weak when comparing to the southern crest. In general, the TMD crests are very asymmetric, while the TEC crests only show a minor asymmetry over the SAA. Although both TEC and TMD show a common trough along the magnetic equator, these distinct features show no clear spatial correlation.

Both TEC and TMD temporal PCA components show a clear common driver-response relationship to solar flux, plus a small annual modulation. The corresponding time-dependent components to Fig. 1 are presented in Fig. 2, in terms of the F10.7 solar radio flux and annual cycle. The annual variations shown in Fig. 2 have been normalized to common solar flux F10.7 = 110 sfu. Note as well that we only investigate the first PCA mode, and the North-South annual variation caused by the angle between the ecliptic and the equatorial plane (including also Earth-Sun distance) is usually reflected in the second PCA components (Calabia and Jin 2016, Fig. 6b). In both TMD and TEC cases, clear quadratic dependencies to the F10.7 solar radio flux are seen in Fig. 2a and c. Concerning the annual modulation of the first PCA component, both cases have shown to increase with solar activity (not shown in this manuscript), and during equinox seasons, with higher values in December than in June (Fig. 2b). However, no clear asymmetries are depicted between March and September for any of the both TEC or TMD analysis. The annual modulation in the first PCA components shows a larger contribution for TEC, by approximately one order of magnitude, while both have a very similar shape. The maximum peaks in equinox both show to be approximately 10 days delayed from March equinox (day-of-year ~90) and approximately 1 month from September equinox (day-of-year ~300). The minimum peaks in solstice both show to be delayed 20 days from June solstice (day-of-year ~191) and approximately 1 month from December solstice (day-of-year ~15).

4 Conclusions and Discussion

In this paper, we present our first results on the main PCA mode of TEC and TMD variability from 16 years of IGS TEC GIMs and from 13 years of TMD estimates at 475 km altitude inferred from GRACE accelerometers and POD. The sparse nature of the initial observations has been successfully synthesized by the PCA to derive the main spatiotemporal patterns that can be compared in a common frame. The dominant patterns of both TEC and TMD variations represented by the first PCA mode show a main dependence on the solar flux with a minor annual modulation. The effect of the varying Earth-Sun distance along the ecliptic plane for a seasonal asymmetry is clearly shown in Fig. 2b and d with higher values in December than in June. Although both TEC and TMD annual variations in the first PCA component show a very similar shape, the annual contribution in TEC is approximately one magnitude larger than that in the TMD variations. These results are in agreement with the common knowledge that neutral species of the upper atmosphere are heated and ionized by radiation from the Sun, resulting in the atmospheric expansion and the creation of electrically charged species, including free electrons and electrically charged atoms and molecules. However, as seen by the spatial analysis, although both EIA and EMA troughs show a clear alignment to the dip magnetic equator, there is no clear spatial correlation between the main spatial structures defined by TEC and TMD (Fig. 1). The EIA crests are clearly located at about ±20° dip latitude, but the EMA lacks such a clear two-crest structure. Moreover, a clear asymmetry is shown in the global distribution of TMD, with higher values in the southern hemisphere than that in the northern hemisphere. The most characteristic discrepancy is a clear TMD enhancement near the southern magnetic pole (dip = 90°). None of these other structures are reflected by the TEC distribution. These results might suggest the solar irradiance as the main driver of TEC variations, and that TMD variations are strongly driven by solar wind and magnetospheric forcing.

Fig. 2
figure 2

Solar flux (left panels) and annual (right panels) contributions to the first time-expansion PCA component of TEC (top panels) and TMD at 475 km altitude (bottom panels). Dimensionless quantities. Annual cycles normalized to F10.7 = 110 sfu. Fitting is shown in black solid lines

Full size image

The hemispheric asymmetry in TMD may be caused by a long-term forcing, related to processes under favorable solar wind input and particle precipitation over the southern magnetic dip. These processes would produce a higher thermospheric heating and subsequent TMD increase in the southern hemisphere (Huang et al. 2014). Previous results have shown that the southern hemisphere TMD is more susceptible to magnetospheric driving conditions (Ercha et al. 2012; Calabia and Jin 2019), and a larger conductance in the southern hemisphere also has been reported in many studies (e.g., Sheng et al. 2017; Wilder et al. 2013, Förster and Haaland 2015, Lu et al. 1994). Other studies on thermospheric temperature and winds (Maruyama et al. 2003; Lei et al. 2010; Drob et al. 2015) could suggest that the EMA trough along the magnetic equator is caused by temperature reduction, due to diverging meridional winds (adiabatic cooling due to upward vertical winds) associated with parallel ion drag, and that equatorward winds during elevated solar and magnetospheric activity could modulate ion drag with the consequent increase of temperature and TMD at the EMA. This could indicate that the long-term forcing related to magnetospheric activity and magnetic field configuration in the southern hemisphere could produce stronger disturbances to thermospheric TMD than that seen to ionospheric TEC.

By systematically improving the overall consistency between models and observations, future research might lead to new geophysical insights and better quantitative characterization of observations. Although the response of the IT system to space weather is still not well understood, the global nature of present satellite observations, and the ability of their advanced scientific instruments, provide the unprecedented opportunity to measure the upper-atmospheric conditions at a greater resolution and accuracy, enabling us to calibrate and constrain the existing models. Monitoring and predicting the Earth’s upper-atmospheric processes driven by solar activity is highly relevant to science, industry, and defense. These communities emphasize the need to increase efforts for a better understanding of the IT responses to highly variable solar conditions, as well as how detrimental space weather affects our life and society.