© F.S. Prol et al., Published by EDP Sciences 2021

1 Introduction

The Earth’s plasmasphere is an ionized region of the atmosphere sensitive to the geomagnetic activity. During quiet geomagnetic conditions, the plasmasphere is dominated by the dayside plasma that flows from the ionosphere upwards along the magnetic field lines (Belehaki et al., 2004). This dayside ionospheric filling in the presence of the Earth’s rotation produces a quite large torus of cold plasma, which becomes quite disturbed during active geomagnetic conditions. Numerous studies have analyzed the disturbed plasmasphere with particular interest in the internal dynamics (Singh et al., 2011), boundaries (Verbanac et al., 2015; Heilig & Lühr, 2018), erosion by enhanced magnetospheric convection (Carpenter, 1970; Zhelavskaya et al., 2017) and coupling processes with the magnetosphere (Borovsky & Valdivia, 2018) and ionosphere (Pierrard & Voiculescu, 2011). However, the electron density description during active geomagnetic conditions is still a relevant challenge for any imaging systems aiming describing the plasmasphere dynamics and structure.

Many imaging systems have been developed in the past few decades to describe the plasmasphere electron density (Gulyaeva et al., 2011; Katus et al., 2015). Most of the systems rely on whistler data (Carpenter & Anderson, 1992), climatological models (Gallagher et al., 2000; Webb & Essex, 2000) and in-situ observations (O’Brien & Moldwin, 2003; Huang et al., 2004; Liu et al., 2015; Zhelavskaya et al., 2017). More recently, the Global Navigation Satellite System (GNSS) has provided a unique opportunity to observe the electron density of the plasmasphere with global and continuous coverage. Considerable progress using GNSS observations has been achieved to represent the plasmasphere in Vertical TEC (VTEC) maps projected on a thin shell located at a specific altitude above the Earth’s surface (Chen et al., 2017). The developed VTEC calculations are useful tools for applications concerned with the horizontal morphology of the plasmasphere; however, since the GNSS signal travels at large altitude regions through the plasmasphere, tomographic algorithms allow the extraction of the vertical distribution, i.e., the exponential decay of the electron density with the increasing altitude.

Tomographic reconstruction techniques are widely applicable approaches to solve inverse problems and estimate the electron density based on TEC measurements. Despite the plasmasphere can contribute about 10–60% in the ground TEC measurements (up to GNSS height at 20,000 km), a few approaches were specifically developed to perform tomographic reconstruction in the region. Main difficulties are related to the low number of available satellites that can provide TEC observations, requiring severe underconstrained methods. Early studies using tomographic reconstructions and satellite-based TEC observations are given by Rius et al. (1997) and Hajj et al. (2000). Later, Heise et al. (2002) developed a technique which provides electron density estimations of 2D meridional cross-sections of each orbit plane of the CHAMP (Challenging Minisatellite Payload) satellite. Gerzen et al. (2015) evaluated the accuracy of these CHAMP reconstructions against IMAGE/RPI (Radio Plasmas Imager) data, showing improvements in comparison to the model used as the background. Spencer & Mitchell (2011) had also proposed a tomographic algorithm of the plasmasphere, in which they assumed the electron distribution as constant along the Earth’s magnetic field lines in the lower plasmasphere. Spencer & Mitchell (2011) used data of the Constellation Observing System for Meteorology, Ionosphere, and Climate/Formosa Satellite Mission 3 (COSMIC/FORMOSAT-3) and Pinto Jayawardena et al. (2016) have shown this method as an efficient tool to correctly represent some trends and patterns of plasmaspheric TEC distributions during quiet geomagnetic conditions. Another method was developed by Wu et al. (2015) based on a 3DVar approach, which has presented a general improvement in comparison to the global core plasma model (GCPM) when performing an internal quality check.

In this work, a new tomographic reconstruction method has been developed to describe the plasmasphere using METOP (METeorological OPerational) data. In comparison to previous methods, our developed method is the first based on a single METOP satellite and we are imposing a set of new constraints into the Simultaneous Iterative Reconstruction Technique (SIRT) (Pryse & Kersley, 1992). Previous methods were applied in quiet conditions; however, the main goal here is to present the potential of the proposed tomographic reconstruction for space weather monitoring during geomagnetic storms. In this direction, a total of four space weather events of the last solar cycle 24 were chosen. As input observations, we have used GNSS TEC data from METOP satellite. Since the satellite is flying at about 800 km, the topside TEC data already excludes a significant part of ionospheric contribution below the low Earth orbit (LEO). In-situ electron density observations from the Defense Meteorological Satellite Program (DMSP) were used as reference for validation purposes. Section 2 shows more details about the dataset. Section 3 presents the mathematical formulation of the developed method. Section 4 shows the experimental results and Section 5 presents the conclusions.

2 Dataset

Three primary datasets were used in this study: GNSS TEC data observed from METOP satellite, electron density measurements from DMSP satellite, and ground-based VTEC data from global ionospheric maps (GIMs). TEC measurements from the dual-frequency GNSS precise orbit determination (POD) receivers of METOP-A satellite were used as the main observations to perform the tomographic reconstruction. The METOP TEC data were processed by the University Corporation for Atmospheric Research (UCAR) and obtained through the COSMIC Data Analysis and Archive Center (CDAAC) via the portal http://cdaacwww.cosmic.ucar.edu/. A detailed review of the TEC estimation procedure as well as its quality evaluation is provided by Yue et al. (2011). All measurements with negative TEC values were neglected in the reconstructions. A cut-off elevation mask of 20° was defined. For the performance assessment of the reconstructed images, the DMSP electron density (Ne) data were used as reference values. The Ne measurements were processed by the Center for Space Sciences at the University of Texas (Dallas) and obtained through the Madrigal portal at http://cedar.openmadrigal.org/. A statistical review of the DMSP Ne quality is provided by Garner et al. (2010). GIMs were used to compute the VTEC contribution of the plasmasphere under geomagnetic storm events by comparing the ground TEC data with the onboard LEO data. The maps released by the Jet Propulsion Laboratory (JPL) were obtained at https://cdaweb.gsfc.nasa.gov/pub/data/gps/ with a spatial resolution of 5° longitude × 2.5° latitude and a temporal resolution of 2 h. JPL products provide smoother solutions in comparison to products with higher temporal resolution. The obtained tomographic reconstructions also provide smooth VTEC distributions due to the used regularization process, so that, the JPL distributions were selected in order to provide more uniform VTEC ratios.

Four weeks were chosen for probing the capabilities of the developed tomographic method during geomagnetic storm conditions. Figure 1 shows the Dst and Kp values during the four space weather events selected during the solar cycle 24. The two first events are related to the St. Patrick storms of 2013 and 2015 and the two last events are related to 2017 and 2018, where very low Dst values were observed during all events. Based on the Kp index, the storms can be classified as strong G3 (2013), severe G4 (2015), severe G4 (2017) and strong G3 (2018).

Fig. 1 Dst and Kp values observed in the selected storm events. Indices values were obtained from the OMNI data set provided by the Goddard Space Flight Center.

3 Tomographic reconstruction of the plasmasphere

The basic quantity required to perform plasmaspheric tomography based on GNSS data is the total electron content (TEC). Typically, TEC is defined as the integral of the electron density (Ne) along the path between the GNSS satellite and the receiving antenna, in a column with cross-sectional area of 1 m². In tomographic applications, however, the plasmasphere, or the ionosphere, is broken down into a 3D grid of numerous cells (voxels) and TEC is approximated by a finite sum. It can be written as:$$\mathrm{TE}{\mathrm{C}}_i=\sum _{j=1}^J\mathrm{}N{e}_j{d}_{\mathrm{ij}}$$(1)where Ne_j is the electron density in el/m³ for a cell j and d_ij is the path length in meters of each signal i inside the boundaries of the voxels.

The main goal of tomography is to estimate the electron density Ne based on TEC measurements. In this sense, a tomographic system is created as below:$$y=\mathrm{Ax}+ϵ$$(2)where y is a TEC vector of observations, x is the unknown vector of electron density, A is the design matrix composed of the path lengths inside the boundaries of the voxels, and ϵ consists of the measurement noises.

The solution of equation (2) involves many unknowns. In the case of high vertical resolutions, the problem is also ill-conditioned due to the limited number of GNSS viewing angles. To overcome this, a priori estimate is required as an initial solution known as a background. The background is usually obtained from a climatological model to fill in each voxel of the grid. Then, the final solution is obtained based on the distribution of the difference between the background TEC and the observed TEC into the resulting grid.

In this study, the background was obtained from a family of electron density models developed at the German Aerospace Center (DLR) in Neustrelitz. They are named Neustrelitz TEC Model (NTCM) (Jakowski et al., 2011), Neustrelitz Peak Density Model (NPDM) (Hoque & Jakowski, 2011), Neustrelitz Peak Height Model (NPHM) (Hoque & Jakowski, 2012) and Neustrelitz Plasmasphere Model (NPSM) (Jakowski & Hoque, 2018). The models are driven by the solar radio flux index F10.7 and a Chapman layer to describe the basic structure of the ionosphere and plasmasphere. The plasmasphere model (in this case NPSM) is superposed to the topside of the Chapman layer. Using such an approach we have introduced required variability in the topside scale height as a function of height, since it is an important point to be considered in Chapman layers, as demonstrated by (Olivares-Pulido et al., 2016). Despite a small number of coefficients and parameters, NTCM, NPDM, NPHM and NPSM provide a good estimate of the main features of the ionosphere and plasmasphere in terms of local time, geographic/geomagnetic location, solar irradiance and solar activity.

After defining the background, the tomographic reconstruction is performed using a 3D grid specifically built to fit into the orbital geometry of the METOP satellite. Figure 2, for instance, shows the coverage of the METOP-A satellite for day of year (DOY) 076 of 2013. Our investigations shows that when superimposing METOP data collected during the rising (or setting) orbit phase for an entire day, the measurements covers the entire globe with a very repetitive pattern in terms of local time (LT). This means that the local time of the observations repeats approximately at the same latitude for each revolution. In this regard, it is possible to superimpose the TEC data of the entire day and represent the system geometry with two 3D grids: one for satellite rising phase and another for setting phase.

Fig. 2 METOP-A orbit footprints for DOY 076 of 2013 in terms of latitude, longitude and local time. Left panel shows the satellite rising phase and the right panel shows the setting phase. The unit of the color bar is local time.

Both grids were defined with a temporal resolution of 1 day in order to cover the entire globe. The horizontal resolution of the 3D grids was defined as 2° × 15° in latitude and longitude, respectively, with a vertical resolution of 50 km, ranging from 800 to 20,000 km in order to cover the entire trajectory of the GNSS signals. The low longitude resolution of 15° was defined due to the longitudinal gaps between each METOP revolution, but the 15° resolution was enough to guarantee a reasonable level of intersections between signals from distinct orbital revolutions. A total of 880,425 cells are then created for each satellite revolution scenario.

Given the high amount of unknowns to cover the large region of the plasmasphere, a relatively light solution is required to solve the equation (2), such as the Algebraic Reconstruction Technique (ART) (Austen et al., 1988). ART is defined as a basis for the developed approach since it is very efficient for providing fast updates of the tomographic system. ART is based on the following iteration process:$$x_j^{k+1}=x_j^k+\gamma \frac{y_i-\sum _j^N\mathrm{}{A}_{\mathrm{ij}}{x}_j^k}{\sum _j^N\mathrm{}{A}_{\mathrm{ij}}^2}{A}_{\mathrm{ij}}$$(3)where γ is a relaxing parameter, y_i is the ith TEC observation, A_ij is the value of the corresponding i and j position in the A matrix, $x_j^0$ is the electron density given by the background, $x_j^k$ is the electron density at iteration k obtained after each update. Typically, ART produces minimum residuals after a few iterations, however the final solution is unlikelly to be stable due to the poor observation geometry. In order to stabilize the tomographic reconstruction, three main adaptations were used: (1) a simultaneous approach; (2) incorporation of VTEC values; and (3) a set of contraints.

ART updates the electron density voxels after each TEC observation is analyzed. This approach is affected by the way the observation vector is sorted. To overcome this problem, the Simultaneous Iterative Reconstruction Technique (SIRT) (Pryse & Kersley, 1992) is used in this work, where the final image is obtained by taking the average correction after the differences $y_i-{\sum }_j^N\mathrm{}{A}_{\mathrm{ij}}{x}_j^k$ of all i signals have been computed.

In addition to slant TEC, VTEC values are included as observations. Two values were therefore used in the iterative reconstruction: the slant TEC, related to the voxels of the real GNSS observation geometry, and the vertical TEC, related to the vertical column of voxels above and below the pierce points. As pointed out by Prol et al. (2019), the combination of VTEC and slant TEC is viable to increase the number of voxels with TEC measurements without significantly affecting the slant TEC estimations. The mapping function developed by Foelsche & Kirchengast (2002) is used in this study, with a single-shell height of 1300 km, since it has demostrated better results in the recent study by Zhong et al. (2016).

As reported in previous studies (Hobiger et al., 2008; Wen et al., 2010; Norberg et al., 2018), the use of constrained methods can efficiently stabilize the ill-conditioned solution in the presence of noises at altitudes of the ionosphere, from 50 to 2000 km. In addition, constrained methods can reduce excessive dependence on the background information, which helps to detect local and regional structures that are missing in the climatological models. The region of interest now, however, is much larger and composed by very few electron densities at high altitudes. Any small TEC error or mismodelling given by the background can effectively introduce artifacts in the reconstruction at high altitude regions. In order to stabilize the ill-conditioned geometry, the proposed algorithm employs a significant dependence on the background specifications in regions with very low electron density. This is performed using the following relaxing parameter:$${\gamma }_{\mathrm{ij}}=0.2(x_j^0/x_{\max }^0{)}^{w_{{\varphi }_m}}$$(4)where $x_j^0$ is the background electron density at cell j and $x_{\max }^0$ is the maximum electron density of the background in signal i. Differently from equation (3), γ varies with each signal i and voxel j and, now, the incorporation of the ratio ($x_j^0/x_{\max }^0$) keeps the profile shape of the final iteration very similar to the patterns provided by the background, mainly when the electron density $x_j^0$ is very small. Similar constant values of 0.2 were adopted by Pryse et al. (1998) and Prol et al. (2018) to guarantee a certain control in the iteration convergence. The weight $w_{{\varphi }_m}$ is a factor dependent on the geomagnetic latitude, and its use allows a reconstruction which relies more on the background in the high-latitudes. The weighting factor is defined as:$$w_{{\varphi }_m}=1+{\sin }^2\left({\varphi }_m\right)$$(5)where ϕ_m is the geomagnetic latitude.

It is worthy to mention that a similar regularization than proposed by Heise et al. (2002) is used for voxels without any ray passing through. Three gaussian functions are used to determine the amount of influence to the nearby voxels, where two equations are used to control the horizontal variability and the third one is used to keep the general characteristics of the vertical structure of the background. The Gaussian half widths are taken as σ_ϕ = 7° and σ_λ = 21° along the latitude and longitude direction, respectively.

In comparison to previous methods (Heise et al., 2002; Spencer & Mitchell, 2011; Wu et al., 2015), the proposed method has the following new features: (1) it uses observations only from single METOP satellites; (2) it performs tomographic reconstruction using two 3D grids distributed in local times for rising and setting phase of the satellite orbit; (3) it imposes a set of new constraints into SIRT; and (4) the number of signal cross-sections are increased since we include a whole day of data. Despite 3D grids are provided, our method was not developed to show the 3D distribution every hour. We can only show the global distributions of the fixed local times related to the rising and setting hours of the single METOP satellite.

4 Results and discussions

4.1 Tomographic responses to geomagnetic storms

Figure 3 demonstrates how well the developed tomography is adapting to the spatial distribution of the observations. Images of the top panels represent VTEC maps during the rising phase of the METOP-A satellite and images in the bottom panel represent corresponding VTEC maps during the setting phase. The rising phase occurred during the daytime around 9.5 h LT where the setting phase occurred around 21 LT, explaining the high level of ionization in the top panels. Despite the background is quite constant along the longitudinal VTEC distributions, there are specific variabilities that tomography has correctly updated, such as the larger TEC values in the west-sector (−180 to 0) and smaller TEC values in the east-sector (0–180) especially during the rising scenario. These results indicate that the developed algorithm is not so constrained to the background and capable to capture the general distributions of the observations, mainly at low- and mid-latitude regions.

Fig. 3 Comparison between VTEC maps obtained with the developed tomography, background, and used METOP observations (DOY 076/2015). Images of the top panel are related to the METOP orbit rising phase and in bottom panel are related to the METOP orbit setting phase. The unit of the color bar is TECU.

To investigate the capability of the developed method for capturing plasmaspheric and ionospheric responses to the geomagnetic storms, Figure 4 shows the tomographic results in terms of VTEC maps above LEO orbit height (up to GNSS satellite height) during selected space weather events. For comparisons, the VTEC images are given during a quiet day, the day when the main phase of the storm has perturbed the plasmasphere and the day after the main phase. The geometrical positions of the satellite are related to the rising phase of the MTEOP-A, where TEC observations were taken around 9 h LT in the low-latitude region of the plasmasphere. In general, the ionization is enhanced during the main phase of the storm and significantly reduced in the day during the recovery phase. This indicates a positive and negative TEC bias compared to quiet day values. It is worthy to mention that the background TEC maps during the storm phases are very similar to each other since the background is mainly driven by the F10.7 index and not geomagnetic or disturbance indices. However, the tomography has represented a subtle daily variation due to the geomagnetic events, meaning that the developed technique is capable of updating the background for displaying positive and negative TEC bias during the storm events. Additionally, we can observe longitudinal asymmetries in the reconstructions. They can be related to several reasons, such as: (1) the irregular geometrical distribution of TEC observations, (2) daily variabilities of the solar ionization, (3) longitudinal variations of the zonal electric field and (4) distinct ionization level due to the geomagnetic field configuration, i.e. magnetic declination, magnetic field strength, and displacement of the geographic and magnetic equators (Kil et al., 2011).

Fig. 4 Daily VTEC maps above LEO orbit height during analyzed storm days obtained with tomographic reconstruction based on METOP TEC measurements during satellite rising. The unit of the color bar is TECU and the given LT hour is computed in the equatorial region.

The positive and negative variabilities presented here corroborates to previous studies, when analyzing the same events of the ionosphere (Astafyeva et al., 2015; Fagundes et al., 2016; Yue et al., 2016; Krypiak-Gregorczyk, 2019). These studies, together with the work of Sharma et al. (2011), explain that the TEC enhancement in the ionosphere during storm involves the increases in the oxygen density, changes in the meridional winds, intensification of eastward electric fields, traveling ionospheric disturbances and penetration of electric fields. Negative biases are suggested to be caused by changes in neutral composition associated with the decrease in the O/N₂ density ratio. In the case of the plasmasphere, however, the ionospheric TEC variabilities are not a direct factor, but act as a precursor to the enhancement or decreasing of the dayside ionospheric filling related to the Earth’s rotation. Additionally, Lemaire et al. (2005) demonstrated the role of the Ring Current betatron mechanism for up-lifting the altitudes of the mirrors points and guiding centers of all ionospheric and magnetospheric particles when the equatorial magnetic field intensity decreases during the main phase of geomagnetic storms. Conversely, during storm recovery phases, the geomagnetic betatron mechanism moves the mirror points of all charged particles downwards back into the denser layers of the atmosphere. These upward and downward motions of ionospheric ions and electrons contribute respectively to the positive and the negative variation of the ionization density observed during geomagnetic storms in the topside ionosphere. During the successive main and recovery phases the kinetic energy of all charged particles is decreasing and increasing, as a result of the betatron mechanism which has been called the LBS mechanism, according to the initials of the authors who pointed it out. The LBS betatron mechanism adds on top of the various other geophysical mechanisms proposed previously for accounting of positive and negative ionospheric storm effects of the analyzed storms.

Figure 5 shows similar images as Figure 4, but the reconstructed maps are now showing the electron density at 800 km height. Figure 4 shows signatures of enhanced ionization at high-latitudes regions especially during the Saint Patrick storm days although corresponding VTEC maps in Figure 4 (main storm day at high-latitudes) do not show such enhancement. After the main storm phase, we can observe the negative bias related to the recovery phase. Similar analyses have been provided by numerous previous studies (Astafyeva et al., 2015; Fagundes et al., 2016; Yue et al., 2016; Krypiak-Gregorczyk, 2019). Our investigation, however, shows that the storm signatures at high-latitudes are more evident in the electron density level than those in VTEC level.

Fig. 5 Electron density distributions at height 800 km for analyzed storms obtained by tomographic reconstruction using METOP during satellite rising. The unit of the color bar is el/m3 and the given LT hour is computed in the equatorial region.

A clear indication of the positive and negative biases is provided by Figure 6. This figure shows the percentage of enhancement of the electron density in comparison to quiet days, defined here as DOYs 74/2013, 74/2015, 249/2017 and 236/2018. It is noticeable that the positive bias can lead to general electron density enhancements of 50–80% at all latitudes, depending on the solar activity, while there is a majority of negative bias in the recovery phase, reaching −40% during the daytime especially at mid- and high-latitudes.

Fig. 6 Percentage of enhancement of the electron density at 800 km during the days of the storm in comparison to quiet days (74/2013, 74/2015, 249/2017 and 236/2018). Units of the color bars are given in percentage.

4.2 Plasmaspheric contribution to ground-based TEC

In order to assess the percentage contribution of the plasmaspheric VTEC to the total VTEC observed from ground stations, daily mean VTEC values obtained from ground-based GIMs. For each storm event, a daily mean VTEC value of two different local times (9.3 and 21.0 LT) corresponding to the satellite raising and setting passages is computed. Figure 7 presents the evolution of the mean VTEC derived from the tomographic method during the period under investigation in absolute values (Fig. 7a) and its percentage variation relative to a quiet day (Fig. 7b). The global mean VTEC value was computed for each day considering latitudes between 40° S and 40° N. One can note that the model is capable capturing the daily variations of the quiet plasmasphere (first three days of each week) and also the VTEC depletion after the geomagnetic storms (last three days of each week).

Fig. 7 (a) Daily average of the VTEC derived from the tomography method at 9.3 and 21.0 LT. The average was computed based on data between 40° S and 40° N latitude. METOP data was missing to perform reconstructions at DOY 79/2013. (b) VTEC percentage relative to the first day of the time series, considered here as a quiet day (low Kp values). The dashed lines indicate the day of the storm onset. The VTEC unit is given in TECU.

In order to determine the daily mean percentage contribution of the plasmasphere, the global mean plasmaspheric VTEC from tomography is then divided by the global mean JPL GIM VTEC. It is important to highlight that since the GIMs and plasmaspheric model present different spatial and temporal resolution, a bilinear interpolation in time and space of the GIMs has been performed. Figure 8 shows the evolution of the mean percentage contribution of the plasmaspheric VTEC in comparison to the global VTEC values derived from the JPL GIMs. As expected, comparing plots at 9.3 and 21 LT we see that the plasmasphere contribution is higher during the nighttime hour in comparison to the day time hour. During daytime, the ionospheric electron content is rapidly enhanced due to the increasing of solar radiation, therefore reducing the percentage contribution of the plasmasphere. At night, the ionospheric electron content decreases considerably in one hand due to the absence of photo ionization on the other hand due to ongoing recombination (of free ions and electrons) and transport processes. Therefore, the percentage contribution of the plasmaspheric TEC to the GIM TEC is increased (Chen & Yao, 2015). In addition, one can note that the percentage contribution of the plasmaspheric TEC increases from 2013 to 2018. This is expected, since during the descending phase of the solar cycle the ionospheric TEC is reduced more than the plasmaspheric TEC.

Fig. 8 Percentage contribution of the plasmaspheric VTEC to the global VTEC derived from the JPL GIMs for 9.3 and 21.0 LT. The dashed lines indicate the day of the storm onset.

The presented analysis shows that the contribution of the plasmaspheric electron content to GPS TEC is in agreement with previous investigations suggesting that ground-based VTEC values may be impacted from about 10% up to 60% during the daytime and nighttime, respectively (Yizengaw et al., 2008; Chen & Yao, 2015). Figure 8 also shows that this contribution is considerably deviated during geomagnetic events. Indeed, one can note in all the investigated events that the mean percentage contribution of the plasmaspheric VTEC to the GIMs VTEC tends to reduce on the day following the main phase of the storm. The absolute values of the plasmaspheric VTEC (Fig. 7a) are likely reduced on the day following the storm onset due to the reduced plasma supply form the ionosphere (Cherniak et al., 2014). After the ionospheric recovery, the plasmasphere density starts to slowly rise and present an plasmaspheric contribution enhancement, such as shown in DOYs 78/2013, 253/2017 and 240/2018. These patterns are observed at 9.3 and 21.0 LT, except for the St. Patrick’s Day event of 2015. This is a more intense storm, with extreme lower Dst values during prolonged time. As result, we observe prolonged time with lower plasmasphere contribution in the daytime (9.3 LT), which is referred to a long period with reduced plasma supply form the ionosphere. During the nighttime of 2015, an enhancement of the percentage contribution on the day right after the storm main phase occurs due to a strong decrease on the global JPL GIM TEC mean, whereas the global plasmaspheric TEC mean reduces slightly.

4.3 Comparison with DMSP electron density data

In order to evaluate the retrieved tomographic images using DMSP data as reference, we present in Figure 9 six-hour averages of electron density values from DMSP, tomography and background. For better visualization, the units are given in log(Ne) and the color bar scales are not fixed. DMSP-18 was selected for years 2013 and 2015, and DMSP-17 was selected for years 2017 and 2018. Such satellite numbers were selected because their orbit is relatively coincident with the METOP orbits in the selected days. Indeed, the DMSP satellite geometry is very repetitive through the corresponding days. For the days used in the assessment, DMSP satellites are located near middle- and high-latitudes of the southern hemisphere at 0 h LT. The nighttime equator measurements are taken around 03 h LT. Several measurements are observed between 06 and 12 LT in the middle- and high-latitudes of the northern hemisphere. Measurements between 18 and 24 h LT are referred to the southern hemisphere. Due to the orbit geometry, just a few minutes are spent in the equatorial region, which creates near straight line of high ionization in the graphs.

Fig. 9 Electron density distributions of DMSP, tomography and background in terms of local time and days of year. White line indicates the beginning of the corresponding storm. DMSP-18 was selected for years 2013 and 2015 and DMSP-17 was selected for years 2017 and 2018. Color scale varies within each plot to adapt to a better visualization of the spatial distributions. Units of color bar are in log(Ne) with base 10.

The white line of Figure 9 indicates the beginning of each corresponding storm. We can see slight enhancements of the DMSP in-situ electron density measurements after the beginning of the four storms. Qualitatively, Figure 9 shows that tomography does an excellent job at replicating DMSP in-situ electron density and capturing enhanced structures after the storm trigger, especially considering that the background was the initial point. The background underestimates the electron density at middle- and high-latitudes and overestimates the nighttime equatorial region. It can be seen a great agreement of the tomography to represent almost all hours and not only during the storm events. The total RMSE considering all DOYs was 1.96 and 2.20 el/m³ for tomography and background, respectively. Therefore, the RMSE has revealed a general improvement of about 11% by the developed tomography. A major problem, however, was found to represent the equatorial region in the dusk (18 LT). As shown in Figure 3, the background VTEC maps are presenting higher values in comparison to the METOP TEC observations at low-latitudes. On the contrary, as shown in Figure 9, the background electron density is lower than the DMSP values for the same instances and locations. Since the electron density given by tomography is adapting to the METOP TEC observations, the tomographic results in the daytime equator are further underestimated when compared to DMSP. A plausible explanation is a misspecification of the background to describe the shape of the profiles of the topside ionosphere and plasmasphere. Indeed, the electron density at DMSP orbit high can be larger than the background while the METOP TEC presents lower values. The real electron density profile is likely uplifted in comparison to the background profile, despite the lower TEC from the real profile. At the dusk instances, the well-known pre-reverse enhancement of the ionosphere (Prol et al., 2019) intensifies the vertical drift of the equatorial ionosphere upward due to the resultant force of the electric field, increasing the peak height and displacing the plasma to higher altitudes. This strong upward vertical drift during the evening pre-reversal enhancement is difficult to be captured by any climatological model due to its relatively small spatial and temporal scales, which makes reasonable to obtain errors at such specific instants.

5 Conclusions

A tomographic technique to reconstruct the electron density based on TEC measurements of the METOP satellites is proposed. Although we have imposed strict constraints in the developed tomographic algorithm, the proposed method was capable to adapt well to the input TEC observations and describe the electron densities during geomagnetic storms. We have checked maximum electron density enhancements around 80% in comparison to quiet days and minimum negative bias around −40% in the recovery phase. In comparison to ground-based VTEC, the results suggested that the plasmaspheric VTEC contribution may vary during geomagnetic events, tending to reduce during the day following the Dst minima. A slow plasmaspheric contribution enhancement then occurs after the ionospheric recovery due to the increased plasma supply. In general, the tomographic method allowed verifying that the averaged plasmasphere VTEC contribution can vary from 10% to 60% in geomagnetic storm conditions. Additionally, the tomographic method was capable at replicating DMSP in-situ electron density and capturing enhanced structures after the storm trigger, revealing a general improvement of about 11% in comparison to background specifications. The proposed methodology, therefore, contributes to the few attempts to perform tomographic estimations in order to use METOP data and describe the electron density of such a large region as the plasmasphere during geomagnetic storms. In future, a key point would be to improve background representations, mainly regarding to the shape of the electron densities, since major problems were identified at the dusk period in the equatorial region, when typical vertical drifts of the ionosphere are effective to produce an intensified ionization at the altitudes of the DMSP orbit.

Acknowledgments

The authors would like to thank the University Corporation for Atmospheric Research, Boulder, CO, USA, and European Space Agency, for providing METOP TEC data (http://cdaac-www.cosmic.ucar.edu/). They would also like to thank the developers of the Defense Meteorological Satellite Program (DMSP), in which data was gathered from the Madrigal database (http://cedar.openmadrigal.org/). Ground-based VTEC data were provided by JPL and obtained at https://cdaweb.gsfc.nasa.gov/pub/data/gps/. This work has been funded by the Deutsche Forschungsgemeinschaft (DFG) under the Grant No. HO 6136/1-1 and by the “Helmholtz Pilot projects Information & Data Science II” (Grant support from the Initiative and Networking Fund of the Hermann von Helmholtz-Association Deutscher Forschungszentren e.V. (ZT-I-0022)) with the project named MAchine learning based Plasma density model (MAP). The editor thanks Joseph Lemaire and an anonymous reviewer for their assistance in evaluating this paper.

References

All Figures

	Fig. 1 Dst and Kp values observed in the selected storm events. Indices values were obtained from the OMNI data set provided by the Goddard Space Flight Center.
In the text

	Fig. 2 METOP-A orbit footprints for DOY 076 of 2013 in terms of latitude, longitude and local time. Left panel shows the satellite rising phase and the right panel shows the setting phase. The unit of the color bar is local time.
In the text

	Fig. 3 Comparison between VTEC maps obtained with the developed tomography, background, and used METOP observations (DOY 076/2015). Images of the top panel are related to the METOP orbit rising phase and in bottom panel are related to the METOP orbit setting phase. The unit of the color bar is TECU.
In the text

	Fig. 4 Daily VTEC maps above LEO orbit height during analyzed storm days obtained with tomographic reconstruction based on METOP TEC measurements during satellite rising. The unit of the color bar is TECU and the given LT hour is computed in the equatorial region.
In the text

	Fig. 5 Electron density distributions at height 800 km for analyzed storms obtained by tomographic reconstruction using METOP during satellite rising. The unit of the color bar is el/m3 and the given LT hour is computed in the equatorial region.
In the text

	Fig. 6 Percentage of enhancement of the electron density at 800 km during the days of the storm in comparison to quiet days (74/2013, 74/2015, 249/2017 and 236/2018). Units of the color bars are given in percentage.
In the text

Fig. 7 (a) Daily average of the VTEC derived from the tomography method at 9.3 and 21.0 LT. The average was computed based on data between 40° S and 40° N latitude. METOP data was missing to perform reconstructions at DOY 79/2013. (b) VTEC percentage relative to the first day of the time series, considered here as a quiet day (low Kp values). The dashed lines indicate the day of the storm onset. The VTEC unit is given in TECU.

In the text

	Fig. 8 Percentage contribution of the plasmaspheric VTEC to the global VTEC derived from the JPL GIMs for 9.3 and 21.0 LT. The dashed lines indicate the day of the storm onset.
In the text

Fig. 9 Electron density distributions of DMSP, tomography and background in terms of local time and days of year. White line indicates the beginning of the corresponding storm. DMSP-18 was selected for years 2013 and 2015 and DMSP-17 was selected for years 2017 and 2018. Color scale varies within each plot to adapt to a better visualization of the spatial distributions. Units of color bar are in log(Ne) with base 10.

In the text