Variations of TEC Over Iberian Peninsula in 2015 Due to Geomagnetic Storms and Solar Flares

The total electron content (TEC) over the Iberian Peninsula was studied using data from two locations obtained both by Global Navigation Satellite System receivers and an ionosonde. The principal component analysis applied to the TEC data allowed us to extract two main modes. Each mode is characterized by a daily TEC variation of a certain type (principal component [PC]) and its amplitude for each of the studied days (given by the empirical orthogonal functions [EOFs]). The variations of these modes as well as the original TEC data were studied in relation to four strongest geomagnetic storms of 2015 and three geomagnetic disturbances of lower amplitude observed during the same months. EOFs were found to correlate well with space weather parameters characterizing solar UV and XR fluxes, number of the solar flares, parameters of the solar wind, and geomagnetic indices. Multiple regression models were constructed to fit EOFs using combinations of the space weather parameters with a lag from 0 to 2 days. Combining the regression models of EOFs with the corresponding PCs, we reconstructed TEC variations as a function of space weather parameters observed in previous days. The possibility to use such reconstructions for the TEC forecasting was also studied.


Introduction
Modern society, industry, and science use Global Navigation Satellite System (GNSS)-based technologies more and more often: operation of unmanned air-born, floating and land vehicles, in-route and assisted landing procedures for commercial aviation, rescue operations, and so forth. The quality of the GNSS signal and therefore the reliability of the GNSS-based technological solutions depend on the conditions in the upper part of the Earth atmosphere-the partially ionized layer called ionosphere.
The main sources for the data on the ionospheric conditions are ionosondes installed at specific locations and networks of GNSS receivers. The ionosondes provide information about the altitude and peak electron density for different ionospheric layers (D, E, F1, and F2) allowing, consequently, to estimate a so-called ionospheric total electron content (iTEC) without plasmaspheric contribution. However, the number of ionosondes is limited (e.g., the Ebro Observatory, Spain, owns the only regularly working ionosonde on the Iberian Peninsula [IP]); they provide information for a certain area and with limited time resolution. As an alternative, the territories of heavily populated areas like the continental Europe are well covered by networks of GNSS receivers. These networks allow to monitor ionospheric conditions with high precision and to develop empirical models to predict ionospheric response to different external forcings, for example, solar flares and geomagnetic storms.
One of the widely used parameters characterizing ionospheric conditions that can be obtained from the GNSS data is the total electron content (TEC). TEC data are often used (Astafyeva et al., 2017;Goncharenko et al., 2013;Li et al., 2019) to analyze and model ionospheric disturbances caused by different events. In general, TEC response to space weather events (like flares or geomagnetic storms) consists of either changes of the amplitude and the shape of the regular daily TEC variation observed worldwide or a phenomenon called traveling ionospheric disturbances (TIDs) which are observed on smaller spatiotemporal scales. In particular, it was shown that at middle latitudes, ionospheric conditions (TEC, TIDs, and scintillation events) respond significantly to strong geomagnetic storms of recent years (Aa et al., 2018;Shim et al., 2018). These variations can be a reason for the GNSS signal degradation and a decrease of the precision of the positioning in the affected area.
The results of the data analysis can be used to develop regional empirical models allowing prediction of ionospheric response to different forcings, for example, geomagnetic storms. In particular, the method called principal component analysis (PCA) showed good potential for the analysis and modeling of ionospheric regular variations and disturbances (Chen et al., 2015;Li et al., 2019;Morozova et al., 2019). The advantage of the regional models (e.g., Hu & Zhang, 2018;Mukhtarov et al., 2018;Petry et al., 2014;Tebabal et al., 2019;Tsagouri et al., 2018) over the global ones is in their better reflection of specific local behavior of the ionosphere.
Large national and international research centers and companies can afford development and operation of complex physics-based models for the upper atmosphere and ionosphere. However, for smaller-scale private enterprises, such models are too expensive and demanding in human resources and computational costs. Thus, there is a demand from small enterprises providing services for local GNSS users for simplistic models forecasting ionospheric conditions for the following 1-2 days. Such models have to be cheap in maintenance, use reasonable amount of proprietary products, and provide reliable forecasts. One of the optimal candidates is the class of statistical models based on the regression of the ionospheric parameters on a set of parameters known to be forcings or precursors of ionospheric disturbances, such as the solar UV and XR irradiance, daily number of solar flares, parameters of the solar wind, and geomagnetic indices.
In this paper, we present results of the exploratory analysis of the TEC data obtained for a midlatitudinal region from both the ionosonde measurements and the data of GNSS receivers. The main goal of this study is to characterize the response of the ionosphere to such events like geomagnetic disturbances, solar flares, and variations of the solar UV and XR fluxes, test regression models based on different sets of regressors (solar [fluxes and flares], solar wind, and geomagnetic parameters) and with different time lags, and explore their forecasting abilities.
The paper is organized as follows. Section 1 gives an introduction. Sections 2 and 3 describe, respectively, the data in use and the methods. Section 4 gives the comparison of different TEC data series and their PCA reconstructions. Section 5 presents analysis of the TEC variations that were observed during March, June, October, and December of 2015 in relation to the geomagnetic disturbances, solar flares, and solar UV and XR flux variations occurring during these months. Section 6 presents the results of the multiple regression analysis of TEC variations on the solar and geomagnetic parameters. The final conclusions are given in section 7.

Data
In this work, we used the following data sets from different ground and space missions.

Total Electron Content
The TEC series are obtained for two locations on IP, on its west and east coasts, respectively, Lisbon (Portugal) and Ebro (Spain). Three sources for the vertical TEC data were used in this work.
The first data set is the GNSS TEC data from the Royal Observatory of Belgium (ROB) database available as vertical TEC maps in the IONEX format on a grid of 0.5°× 0.5°with 15-min time resolution (see also Bergeot et al., 2014). According to the information from the ROB website, the vertical TEC is estimated in near real time from the GPS data of the EUREF Permanent Network (EPN). The vertical TEC maps are produced from the slant TEC of each satellite-receiver pair as projections in the vertical TEC at the ionospheric piercing points (ionospheric shell at 450 km) and interpolated to a grid using splines. Here the TEC data for the grid points most close to Lisbon (39°N, 9°W) and Ebro (41°N, 0.5°E) were used, hereafter ROB-LIS and ROB-EBR, respectively.
Another series, hereafter SCI-LIS, is the GNSS TEC measured by the SCINDA receiver that has been active from November 2014 to July 2019 in the Lisbon Airport (38. 8°N, 9.1°W) in the frame of the ESA Small ARTES Apps project SWAIR (Space Weather and GNSS monitoring services for Air Navigation). The installed equipment is a NovAtel EURO4 receiver with a JAVAD Choke-Ring antenna. The data, originally of 1-min time resolution, were averaged to obtain the 1-h time resolution. The details of the SCINDA software functionalities can be found in Groves (2007), and Groves (2009a, 2009b). Since the data are provided in archived format (Barlyaeva et al., 2020a(Barlyaeva et al., , 2020b and contain some peculiarities (i.e., missed or erroneous blocks) to be taken into account and preprocessed, a software named "SCINDA-Iono" toolbox for MATLAB (Barlyaeva et al., 2020b) was developed by our group.
The calibration procedure was not performed during the installation of this receiver. To validate the noncalibrated SCI-LIS series, we present in section 4 a comparison of SCI-LIS with other TEC series.
The third data source is the iTEC (IONO-EBR) provided by the Ebro Observatory, Spain (40.8°N, 0.5°E, 50 m asl). The instrument currently installed at the Ebro Observatory is the DPS-4D ionospheric sounder, and the measured parameter is the critical frequencies of the ionospheric layers. The altitude profiles of electron density are calculated from the ionograms, and the integration of these electron profiles up to 1,000-km height gives the values of iTEC (Huang & Reinisch, 2001). The IONO-EBR data are of 1-h time resolution.
The 1-h series of TEC for the 4 months of 2015 (March, June, October, and December) are shown in the supporting information, Figure S1. The time is UT (for Lisbon UT ¼ LT, for Ebro UT ¼ LT-1); since the longitudinal distance between Lisbon and Ebro is~9.5°, the solar time difference is~40 min.

Space Weather Data
Five parameters were used to characterize geomagnetic field variations. The global Dst index have 1-h time resolution. The global Kp and ap indices and the local K index (K COI , calculated from the horizontal component of the geomagnetic field measured at the Coimbra Magnetic Observatory, COI, Coimbra, Portugal, 40.2°N, 8.4°W, 99 m asl) originally have 3-h time resolution. The fifth index is the auroral electrojet (AE) index characterizing the auroral activity in the polar regions. Weak and strong geomagnetic disturbances (with Dst ≤ −50 nT) were analyzed separately with a threshold between them set at −100 nT.
The data on the solar wind properties were obtained from the OMNI database. In this study, we used the following parameters with the original 1-h time resolution: interplanetary magnetic field (IMF) components as scalar B, Bx, GSM By and Bz (in nT), solar wind flow speed (v in km/s), proton density (n in n/cm 3 ), and flow pressure (p in nPa).
To parameterize the variations of the solar UV radiation, we used two proxies. The first one is the Mg II composite series (Snow et al., 2014), a proxy for the spectral solar irradiance variability in the spectral range from UV to EUV based on the measurements of the emission core of the Mg II doublet (280 nm). The second proxy is the F10.7 index from the OMNI data base. Both Mg II and F10.7 series have 1-day time resolution.
As a proxy for the variations of the solar XR flux, we used the data measured by the Solar EUV Experiment (SEE) for the NASA TIMED mission at the wavelength 0.5 nm with time resolution of 1 day (XR). Please note that for the plots, the UV and XR irradiance proxies were scaled individually for better visualization.
The information about the solar flares observed during the analyzed time interval was obtained through the NOAA National Geophysical Data Center (NGDC). The daily numbers of solar flares of classes C, M, and X separately as well as the total daily number of such flares (N) were calculated. Only flares that occurred during the local day time were considered.
The variations of the geomagnetic indices, the solar UV and XR fluxes, the number of solar flares, as well as parameters of the solar wind and the IMF for 2015 are presented in Figure S2 (daily means).

Methods
All series used in the presented work were averaged or interpolated to obtain series with 1-h and 1-day time resolution. A number of small gaps were linearly interpolated.
The analysis of the TEC series obtained from different sources was performed using the PCA. The input data set is used to construct a covariance matrix and calculate corresponding eigenvalues and eigenvectors. The eigenvectors are used to calculate principal components (PCs) and empirical orthogonal functions (EOFs). The combination of a PC and the corresponding EOF is called a "mode," and the eigenvalues allow us to estimate the explained variances of the extracted modes. PCs are orthogonal and conventionally nondimensional. The full descriptions of the method can be found in, for example, Bjornsson and Venegas (1997), Hannachi et al. (2007), and Shlens (2009).
Each of the TEC series obtained at different locations and by different instruments (ROB-LIS, ROB-EBR, SCI-LIS, and IONO-EBR) with 1-h time resolution was analyzed separately during each month-long time interval. The PCA input matrices were constructed in a way that each column contains 24 observations (every 1 h) for a specific day. Thus, PCA allows us to obtain daily variations of different types as PCs and the amplitudes of those daily variations for each day of a month as corresponding EOFs. Thus, we need no models for the daily variation of TEC or its seasonal changes as function of zenith angle or day of year (DOY): the amplitude and the phase (daily maximum local time) of the daily TEC variation are automatically extracted from the observational data. Further, the first and second PCA modes were used to reconstruct TEC variations.
Similarities between the variations of TEC and space weather parameters were analyzed using the correlation coefficients, r, that test linear relations between analyzed variables. The significance of the correlation coefficients was estimated using the Monte Carlo approach with artificial series constructed by the "phase randomization procedure" (Ebisuzaki, 1997). The obtained statistical significance (p value) considers the probability of a random series to have the same or higher absolute value of r as in the case of a tested pair of the original series.
The linear multiple regression models (MRMs) were constructed using TEC-related series (1-h and 1-day TEC, EOF1, and EOF2) as dependent variable and solar (i.e., UV and XR fluxes and the flare numbers), solar wind, and geomagnetic parameters as regressors. A "best subset" technique was chosen to ensure that only those regressors that are most influential for a particular TEC series were selected. The "best subset" was estimated using a so-called adjusted squared coefficient of multiple determination (R adj 2 ). The value R adj 2 ·100% shows the percent of the variations of the dependent parameter explained by the MRM in question.
To estimate the quality of the MRMs, we also used the following statistics: r, the correlation coefficient between the observed and modeled series; RMSE/σ, the root mean square error in units of the standard deviation (σ) (see Equation 1); and the χ 2 goodness-of-fit test (Equation 2).
where TEC is the observed and TEC MRM is the modeled TEC series and N is the length of the data series.

Comparison of TEC Series From Different Locations and Instruments and PCA Reconstructions
In this work, we used series of TEC obtained for two different locations on IP and by different instruments. Therefore, before further analysis, we want to discuss similarities and differences of these series during the analyzed time interval.
The correlation analysis of the TEC series from different locations shows high correlation: mean value for the four analyzed months r ¼ 0.93. There is also good agreement between the TEC values measured by the GNSS receivers and the ionosonde: mean r ¼ 0.95. The SCI-LIS series is well correlated both with other GNSS-based TEC series (mean r ¼ 0.95) and with iTEC (mean r ¼ 0.89). The fact that SCI-LIS is not calibrated does not affect the results of our analysis since we used only methods insensitive to scaling and shifting of the series. Please note also that for the plots, the SCI-LIS series was scaled to achieve better visualization (see Figures S1, S3, and S4). All correlation coefficients have p values ≤ 0.05. Individual correlation coefficients can be found in Table S1.
There is a systematic difference between the ROB and ionosonde series: the GNSS-based TEC values are higher than the iTEC series by~3-5 TECU, which, to our mind, reflect both the difference between the iTEC and TEC values (see discussion about plasmaspheric contribution in, e.g., Belehaki et al., 2003Belehaki et al., , 2004McKinnell et al., 2007) and the fact that the ROB data are interpolation of the actual observations to a regular grid.
Since the TEC series are highly correlated, we calculated the mean TEC series for each month using the normalized (varying between −1 and +1) individual TEC series. The mean TEC series are in arbitral units (a.u.).

Space Weather
The individual TEC series for each month were submitted to PCA as described in section 3. Here we will discuss only the first two PCA modes, Mode 1 and Mode 2. The mean PC1 and PC2 for both modes are shown in Figure 1 and PCs for individual series are presented in Figures S3 and S4. Corresponding mean EOF1 and EOF2 are shown in Figures 2 and 3, respectively, and the individual EOFs are shown in Figures S5 and S6.
Mode 1 explains 93-95%, 77-86%, 92-95%, and 87-94% of the TEC variations for March, June, October, and December, respectively. As one can see in Figure 1, PC of Mode 1 represents daily TEC variations caused by the regular changes of the insolation through a day. Daily minimum is observed during the local night and the daily maximum is in the afternoon hours. The differences between mean PCs for different months reflect seasonal changes of the ionospheric insolation and other climatologic conditions in the upper atmosphere.
Mode 2 of TEC series is defined as PC2/EOF2 for March and June for all TEC series and for October and December for the GNSS TEC series, whereas for the iTEC series, this mode rather consists of PC3/EOF3. Mode 2 represents daily TEC variations with a relatively shallow minimum around the noon and a maximum in the late afternoon (19-21 h). Mode 2 explains 2.4-2.9%, 6.1-8.4%, 1.5-3.0%, and 2.5-3.7% of the TEC variations for March, June, October, and December, respectively. The variations associated with Mode 2 are relatively well correlated between the locations and instruments (mean r ¼ 0.71, see also Table S1). Please note that for October and especially for December 2015, the iTEC Mode 2 shows lower correlation with other (GNSS) TEC series. The reason for this, probably, is that the variations associated with this mode for these months for the iTEC series are distributed between modes PC3/EOF3 (more) and PC2/EOF2 (less).
Overall, the first two PCA modes together explain 90-95% of the variations of the original series. Therefore, we used these modes to make the PCA reconstructions of TEC. The mean TEC Mode 1 and Mode 2 series can be found in Figures S7-S10.

TEC Variations in 2015
In   Also, on 20 March, a partial solar eclipse with the maximal obscuration of 60% at 09:00 UT (IP area) took place, resulting in the anomalously low amplitude of the TEC daily variations.

Geomagnetic Storms
The geomagnetic storm on 17-18 March 2015 (GS1, see Table S2) is characterized by Dst < −220 nT and Kp ¼ 8. This storm and its effect on the ionosphere is thoroughly analyzed by many authors (see, e.g., Balasis et al., 2018;Fagundes et al., 2016;Liu et al., 2015;Liu & Shen, 2017;Nava et al., 2016;Paul et al., 2018;Piersanti et al., 2017). The storm was triggered by two flares of C type on 14 and 15 March and two halo CMEs arrived at 1 AU on 17 March (Liu et al., 2015). The recovery phase of the storm started on 18 March and lasted for several days; Dst remained below −50 nT until the end of March.
The Earth's ionosphere responded to this storm by an increase of TEC during the main phase ("positive ionospheric storm," see Maruyama et al., 2009) on 17 March which was followed by a decrease of TEC ("negative ionospheric storm," Maruyama et al., 2009) during 18-19 March Fagundes et al., 2016;Liu & Shen, 2017). For Chinese middle latitudes, the decrease of TEC relative to the undisturbed level was 60-70% (Liu & Shen, 2017). In the European-African sector, the positive storm on 17 March was observed only in the Northern Hemisphere , and the daily TEC variations at 30-40°N were characterized by a second peak around 16-20 h LT Nava et al., 2016). The main mechanism behind the ionospheric perturbation during this time intervals is, most probably, the prompt penetration of electric field (PPEF) (Fagundes et al., 2016;Nava et al., 2016;Paul et al., 2018).
The TEC variations observed for the IP area (see Figures 2, 3, and S7) also show positive-negative ionospheric storm on 17-18 March. The change of the sign of the TEC response to the geomagnetic storm is most clearly seen in variations of EOF1 (Figures 2 and S5): the amplitude of EOF1 on 17 March reaches the highest value for this month, and the amplitude of EOF1 on 18 March is even smaller than the EOF1 monthly mean level. We also observe the second daily peak on 17 March, which is clearly seen in the variations of the EOF2 amplitude (Figures 3 and S7 [bottom]).
The geomagnetic storm on 22 June (GS2, see Table S2) is thoroughly described in, for example, Astafyeva et al. (2016Astafyeva et al. ( , 2017Astafyeva et al. ( , 2018 The ionospheric response to this storm also consisted of the positive (22 June) and negative (23-24 June) phases observed in all latitudinal sectors. In the middle latitudes of the European-African sector, TEC increased (relative to the undisturbed level) on 22 June during the evening, presunset hours, and even after the sunset (Astafyeva et al., 2017). On 23 June, TEC started to decrease on both dayside and nightside of the Earth (Astafyeva et al., 2017), and dual peaks in TEC diurnal variation were observed at the middle and low latitudes of the Northern Hemisphere (Paul et al., 2018). The main proposed mechanism for the ionospheric response to the geomagnetic storm was, again, PPEF (Astafyeva et al., 2016;Ngwira et al., 2018;Paul et al., 2018) and the disturbance dynamo electric field (DDEF) on 23 June (Paul et al., 2018).
Similar to the March event, the amplitude of the TEC daily cycle over IP significantly increases on 22 June and, again, on 25 June, during the second Dst decrease (see Figures 2, 3, and S8). On the following days (23 and 26 June), the amplitude of the daily variation decreased. These changes are reflected in the variations of the amplitude of Mode 1 (EOF1). Similar behavior but with smaller amplitude difference was also observed during the weak disturbance on 8-11 June. Mode 2 of the IP TEC variations shows more complex behavior. Similar to the March event, its amplitude is high and EOF2 is positive on the first day of the main storm (22 June) and negative on the second day of the main storm (23 June). On the contrary, for the second Dst decrease on 25 June, as well as for the weak disturbance on 8-11 June, the amplitude of Mode 2 is high but EOF2 is negative, changing to almost zero but positive on the following day.
On 7-8 October, a storm with Dst ¼ −124 nT and Kp ¼ 6 was registered (GS4 , Table S2). This storm was not associated with any CME and was probably caused by a coronal hole high-speed streams (Matsui et al., 2016;Pazos et al., 2019). No significant ionospheric response to this storm was previously reported.
The variations of TEC observed over IP (see Figures 2, 3, and S9) during this storm can be identified as a negative storm: the amplitude of the daily TEC cycle on the first day of the storm was lower than during previous and following days and they are most clearly seen in variations of Mode 1 (EOF1). Another feature of TEC variation during this storm is that Mode 2 increased on both days of the storm causing a second peak in TEC daily cycle on 7 October and an offset of the daily maximum on 8 October.
The last storm of 2015 started in the evening of 19 December and lasted until 22 December (GS3, see Table S2). Dst ¼ −155 nT with Kp ¼ 6 were observed on the night from 20 to 21 December. This was a two-step storm (Balasis et al., 2018) caused by CME that arrived on 19 December and related to a flare on 16 December (Balasis et al., 2018) resulting in an auroral activity all over the northern Europe (Cherniak & Zakharenkova, 2018). Detailed description of this storm can be found, for example, in During this geomagnetic storm, the positive ionospheric storm was observed in the middle latitudes of the Asian sector on 20-21 December (Paul et al., 2018) without the consequent negative phase. This storm was characterized by the increase of the amplitude of the TEC daily variation over IP (see Figures 2, 3, and S10), as well as Mode 1, on the first day (positive storm). On the second day, the amplitude of the daily variations was still high but lower than on the first day. Probably this ionospheric storm can be classified as a positive-negative storm as well, and the relatively high amplitude of the TEC variations on the second day is related to the two C-class flares that were observed on 21 December. Mode 2 variations were significantly amplified only on the first day of the storm. The small geomagnetic disturbance on 14-15 December (GD3, DOYs 348-349) seems to have no effect on the TEC variations.
Among the four strongest geomagnetic storms of 2015 (GS1 to GS4), only one (GS3) was characterized by just the negative ionospheric storm. Two geomagnetic storms (GS1 and GS2) caused ionospheric disturbances that can be classified as positive-negative storms. Lastly, the geomagnetic storm GS4 resulted in the ionospheric storm that is either positive or positive-negative. The types of the ionospheric storm defined by the original TEC values and by their Mode 1 are the same. Thus, the analysis of EOF1 allows easy classification of an ionospheric storm.
Besides the overall increase or decrease of the amplitude of the TEC variations described by Mode 1, there is another prominent feature of the TEC variations extracted from the data as Mode 2. Mode 2 is associated with appearance, as a rule, of a second daily peak or a sharp decrease in the TEC variations during the afternoon hours. Mode 2 was shown to intensify on the first day of the storm (causing second daily peak in TEC). This behavior is seen for all four storms of 2015. On the second day of the storms in March and June, Mode 2 was also high in amplitude but of the opposite sign, causing a sharp decrease of the TEC in the afternoon hours. During the storm in October, this mode was both significant in amplitude and of the same sign as for the previous days resulting, however, not in second peak but in the offset of the daily peak to the afternoon hours. During the December storm, the amplitude of Mode 2 on the second day was negligible.
Among the geomagnetic disturbances with lower amplitudes (GD1 to GD3), GD1 and GD2 caused positive-negative ionospheric disturbances. The amplitude of Mode 2 was high on the first day of the disturbances with positive EOF2 for GD2 and negative one for GD1. On the second day of both disturbances, EOF2s changed their signs to opposite. The GD3 in December showed no effect on the TEC variations.
The mean correlation coefficients between the EOF1s and EOF2s and the geomagnetic indices (Dst, Kp and K COI , ap, and AE) for each analyzed month are shown in Table 1. The correlation coefficients were calculated with lag ¼ 0-2 days (geomagnetic indices lead). EOF1s anticorrelate with variations of the K and AE indices with lag ¼ 1 day for all months except December (lag ¼ 0 day). EOF2s anticorrelate with variations of the K and AE indices for March and June (lag ¼ 1-2 days) and correlate for October and December (lag ¼ 0 day). For the Dst index, the signs of the correlation coefficients are opposite. The highest (in the absolute value) correlation coefficients were obtained, as a rule, for the Dst and AE indices. The change of the sign of the correlation coefficients between different events reflects, to our mind, differences in the geomagnetic conditions of the studied months, nonlinearity of the relation between the amplitude of Modes 1 and 2 and geomagnetic activity, and differences in the ionospheric response to different storms described above.
Previously, some conclusions were made about dependence of ionospheric response to geomagnetic storms on the season (Andonov et al., 2011;Kutiev et al., 2013 and references therein) or starting time (Borries et al., 2015;Yamamoto et al., 2000; see also Kutiev et al., 2013 and references therein). However, since only seven geomagnetic disturbances appeared during different seasons and started at different time were considered in this work, we can make no definite conclusion about dependence of TEC variations during geomagnetic storms on the season or starting time.
We also calculated correlation coefficients between TEC EOFs and the parameters of the solar wind (Table 1) with different lags (solar wind parameters lead). Both EOF1 and EOF2 show high and statistically significant correlations with solar wind parameters with highest (in the absolute value) correlation coefficients for B, Bx, Bz, v, and p.

Solar Flares and UV and XR Fluxes
The X-flare on 11 March had very weak effect on the TEC variations over IP: there is weak increase of TEC at 17 h (first hourly measurements after the flare). Neither EOF1 nor EOF2 show particular response to this flare (see Figures 2 and 3).
On the time scale of a month, there is dependence between the number of solar flares, solar UV and XR fluxes, and the amplitude of the daily TEC cycle and the amplitude of Mode 1 (or EOF1). The correlation coefficients (see Table 1 and its extended version Table S3) between the EOF1s and the solar UV and XR proxies are r ¼ 0.2-0.86 (depending on the month and the proxy). Generally, the UV flux proxies have higher correlation coefficients than the XR proxy: on average, r ¼ 0.57 for EOF1 versus UV and r ¼ 0.46 for EOF1 versus XR proxies. The lowest correlation is found for December. The correlation coefficients between the EOF2s and the solar UV and XR proxies are low and statistically insignificant. Thus, the overall increase Note. Only r ≥ 0.2 and p value ≤ 0.2 (in parentheses) are shown. The lag in days is shown in brackets (space weather parameters lead). Bold marks highest (in absolute values) correlation coefficients for each month for each group of parameters. The full version of Table 1 is Table S3. of the flare number and/or solar UV and XR fluxes results in the increase of the amplitude of the daily TEC variation and, consequently, EOF1. The dependence of EOFs on the number of flares is weaker and varies with month.

MRM of TEC Response to Space Weather Forcing
The results of the correlation analysis presented in section 5 show relatively high correlations between the TEC series and space weather parameters: the number of the solar flares, UV and XR fluxes, solar wind parameters, and geomagnetic indices. Consequently, linear multiple regression method can be used to reconstruct TEC variations using space weather parameters as regressors. Since the purpose of these MRMs is not to prove physical relations between the ionospheric and space weather parameters but to obtain best prediction quality possible, we can use as regressors parameters that could be related between each other (both through physical processes and statistically). Similar, but not exactly the same, approaches were taken in Chen et al. (2015) and Li et al. (2019) to model TEC variations over large areas in dependence on the solar (F10.7 index) and geomagnetic (Ap index) activity. Also, in the recently published paper by Razin and Voosoghi (2020), the first three PCA modes were used to train artificial neural networks to predict TEC variations.
The MRMs were constructed for EOF1 and EOF2 series as well as for the TEC daily means (1-day TEC series). Then, the MRM fits for the EOF1, EOF2, and 1-day TEC were used, together with previously obtained PC1 and PC2, to reconstruct TEC series (PCA-MRM reconstruction). The advantage of this method is that we do not need to parameterize daily TEC variations and their dependence on the season since these daily variations are already extracted from the data as corresponding PC1 and PC2. The MRMs were constructed using time lags from 0 (no lag) to 2 days (space weather parameters lead by 2 days). Using the lagged MRM, we test the possibility to use this method to forecast TEC variations.
We constructed MRMs both for the mean TEC series and for the series from the individual receiver (SCI-LIS). Since we make MRMs for the series with 1-day time resolution, we used daily mean series of space weather parameters too. Two types of MRM were constructed. The models of the first type use as regressors the solar (i.e., fluxes and flares) and geomagnetic parameters, and the second-type models use as regressors the solar and solar wind parameters. These MRMs were constructed using the "best subset" method to study which parameters are most/less frequently used for best fitting models (models with highest possible for the current set of regressors value of R adj 2 ). The R adj 2 parameters for MRMs with the solar wind parameters as regressors are shown in Table 2 (and for both sets of regression in Table S4) for the mean TEC series and in Table 3 for the SCI-LIS TEC series. The correlation coefficients between the original mean TEC series and the PCA-MRM reconstructions are also shown in Tables 2, 3, and S4. In most cases, the prediction quality of the MRMs using solar wind parameters as regressors is higher than the ones using geomagnetic indices. Solar UV flux proxy (Mg II) as well as the IMF components Bx and Bz are the parameters most often chosen (>43% of the MRMs use them) for "best set" of the regression parameters for the mean TEC series. For MRMs of the SCI-LIS TEC series, the most often selected regressors are solar UV flux proxy (Mg II) and the solar wind flow speed v The lag between the space weather and TEC series is between 0 and 2 days (space weather parameters lead).

Space Weather
(>75% of the MRMs use them), as well as the IMF components B, Bx, and By (58-67% of the MRMs use them). For most of the PCA-MRM reconstruction, the correlation coefficients between the original and reconstructed series are higher for the individual receiver than for the mean TEC series; however, in many cases, the MRM reconstructions of the daily mean TEC or EOFs for the mean TEC series are better than for the individual receiver (please compare Tables 2 and 3).
The PCA reconstructions using original EOFs for the SCI-LIS TEC series are shown in Figure 4, and the PCA-MRM reconstructions using regression on the solar wind parameters with lag ¼ 2 days are shown in Figure 5. The corresponding results for lag ¼ 0-1 days are in Figures S11 and S12. The two first PCA modes allow to reconstruct 90-95% of the original TEC series variability ( Figure 4). In case a higher reconstruction quality is needed, we suggest adding the third mode. The PCA-MRM reconstructions of the original TEC series ( Figure 5) have lower quality: 52-92% of the observed variations are reconstructed (Table 2), depending on the month and the lag.
To test forecasting abilities of the PCA-MRM reconstruction method, we used the data from the receiver SCI-LIS. First, we constructed MRMs for the daily mean TEC and EOFs using the data for a selected number of the days of a month (hereafter, "training days"). The minimal length of the "training days" is 16 for all months except December (due to statistical properties of the regressors' data sets, the minimal length of the "training days" for December is 20 days). We used these MRM reconstructions to forecast the TEC variations for the following day, and then, we gradually increased the length of the "training days." To make reconstruction technically easier, we did not apply the "best subset" selection; therefore, the constructed MRMs may not be the ones with the highest R adj 2 values. We made forecasts using the space weather parameters lagged by 1 and 2 days and also calculated their mean. Some examples of the forecasts based on different length of the "training days" for March 2015 are shown in Figure 6 and in Figures S13 to S15 for other months. To our mind, the use of the "best subset" and combination of different lags for different predictors (e.g., smaller lags for the solar flux proxies and larger for the solar wind parameters) could improve the forecast quality.
The quality of the PCA-MRM forecasting was estimated by calculation of the following statistics: correlation coefficients r, RMSE/σ, and χ 2 between the observed and forecasted series (see Figure 7 for March and

Space Weather
Figures S16, S17, and S18 for June, October, and December, respectively). In general, the quality of the TEC forecast grows with the number of the "training days." The forecasting quality (r, RMSE/σ, and χ 2 ) decreases for days of geomagnetic storms storm: for example, r (RMSE/σ) is significantly low ( The MRMs based on a set of days that already includes a storm have better forecasting abilities. Also, χ 2 gradually decreases to 1 when the number of the "training days" is increased. Averaging of the forecasts with different lags, in general, improves the forecast quality. We plan to test other regression methods, including neural networks based on the machine learning algorithms, especially for the storm days, to find an optimal set of the regressors and number of training days, to test different time lags and their combination, to study seasonal variations of the forecasting quality and make other optimizations. However, to our mind, the results of this exploratory work show good potential for the PCA-MRM forecast method. We also found that for a region as IP, the spatial TEC variations are not essential for forecasting daily TEC variations; therefore, both the TEC series from the individual receiver and mean TEC series can be used.

Conclusions
In this work, we analyzed variations of the TEC over the midlatitudinal area of the IP during seven geomagnetic disturbances of 2015 with Dst < −50 nT that occurred in March, June, October, and December. We also analyzed variations of TEC in relation to the daily total number of solar flares (C, M, and X classes) and the solar UV and XR fluxes.
The TEC data were obtained for two locations on the IP: Lisbon (Portugal) on the west coast and Ebro (Spain) on the north-east coast of the peninsula. Two types of TEC were used in the study: TEC calculated from GNSS receivers and the iTEC obtained from the ionosonde at the Ebro Observatory. All TEC series are in good agreement. The new data series, the noncalibrated TEC from the SCINDA GNSS receiver

Space Weather
installed between 2014 and 2019 in the Lisbon airport area, is found to be well correlated with other series which allowed us to use it for the analysis.
The TEC series were submitted to the PCA allowing to extract two main modes that together explain >90% of the TEC variability. Mode 1 represents daily TEC variation, and Mode 2 is associated with either second daily peak or a sharp decrease during the afternoon hours (19-21 h LT).
Similar changes in the amplitude of the TEC daily variation and EOF1were observed during all studied geomagnetic disturbances. Thus, Mode 1 (or corresponding EOF1) can be used for easy classification of an ionospheric storm.
Mode 2 of the TEC variations during most of the studied storms increased in amplitude with positive EOF2 values on the first day of the storm and negative (or almost zero in case of December storm) EOF2 values on the second day. Thus, the high positive Mode 2 (or high positive EOF2) is a good indicator of a second daily peak in the TEC daily variations.
The overall increase/decrease of the flare number as well as changes of the solar UV and XR fluxes during the analyzed months resulted in the increase/decrease of the amplitude of the TEC daily cycle (and therefore Mode 1). No relation between the amplitude of Mode 2 associated with the solar UV and XR fluxes was found.
The amplitudes of Mode 1 and Mode 2 (EOF1 and EOF2) were found to correlate with space weather parameters: solar UV and XR fluxes, daily numbers of the flares, parameters of the solar wind, and geomagnetic indices. Based on the correlation analysis, we constructed MRMs for the EOFs and the daily mean TEC values with space weather parameters as regressors. MRMs were made with time lags ¼ 0-2 days (space weather parameters lead). The MRM fits for EOFs and the daily mean TEC values were used to reconstruct TEC variations using corresponding PCs (PCA-MRM reconstructions).
We also explored a possibility to use the PCA-MRM reconstruction method to make forecasts of the daily TEC variations. In general, the predictive quality of the PCA-MRM forecasts was found good, but further studies are needed to improve their predictive quality.

Data Availability Statement
Data sets of ionospheric parameters are provided by SCINDA GNSS receiver from Lisbon airport area, Mendeley Data, v1 https://doi.org/10.17632/kkytn5d8yc.1 "SCINDA-Iono" toolbox for MATLAB by T. Barlyaeva is available online at https://www.mathworks.com/ matlabcentral/fileexchange/71784-scinda-iono_toolbox We acknowledge the mission scientists and principal investigators who provided the data used in this research.
The Mg II data are from Institute of Environmental Physics, University of Bremen available at http://www. iup.uni-bremen.de/gome/gomemgii.html (see also Snow et al., 2014 for more information).
The data on the variations of the solar XR flux are from the LASP Interactive Solar Irradiance Data Center (LISRD, http://lasp.colorado.edu/lisird/). LISIRD provides a uniform access interface to a comprehensive set of solar spectral irradiance (SSI) measurements and models from the soft X-ray (XUV) up to the near infrared (NIR), as well as total solar irradiance (TSI). The XR TIMED data are from the SEE measures the solar ultraviolet full-disk irradiance for the NASA TIMED mission. Level 3 data represent daily averages and are filtered to remove flares available at http://lasp.colorado.edu/lisird/data/timed_see_ssi_l3/