Global‐Scale Ionospheric Tomography During the March 17, 2015 Geomagnetic Storm

The correct representation of global‐scale electron density is crucial for monitoring and exploring the space weather. This study investigates whether the ground‐based Global Navigation Satellite System (GNSS) tomography can be used to reflect the global spatial and temporal responses of the ionosphere under storm conditions. A global tomography of the ionosphere electron density is constructed based on data from over 2,700 GNSS stations. In comparison to previous techniques, advances are made in spatial and temporal resolution, and in the assessment of results. To demonstrate the capabilities of the approach, the developed method is applied to the March 17, 2015 geomagnetic storm. The tomographic reconstructions show good agreement with electron density observations from worldwide ionosondes, Millstone Hill incoherent scatter radar and in‐situ measurements from satellite missions. Also, the results show that the tomographic technique is capable of reproducing plasma variabilities during geomagnetically disturbed periods including features such as equatorial ionization anomaly enhancements and depletion. Validation results of this brief study period show that the accuracy of our tomography is better than the Neustrelitz Electron Density Model, which is the model used as background, and physics‐based thermosphere‐ionosphere‐electrodynamics general circulation model. The results show that our tomography approach allows us to specify the global electron density from ground to ∼900 km accurately. Given the demonstrated quality, this global electron density reconstruction has potential for improving applications such as assessment of the effects of the electron density on radio signals, GNSS positioning, computation of ray tracing for radio‐signal transmission, and space weather monitoring.

Besides tomography, data assimilation with numerical models can also be used to estimate the global electron density (e.g., Kodikara et al., 2021) but tomography is different from data assimilation techniques since tomography is not constrained by the internal physics of the model. Besides TOMION, She et al. (2017) developed a global tomography with a parameterized ionospheric electron density model based on eigenmode and leastsquare regression. This technique is often referred to as function-based tomography, where a set of spherical harmonic and empirical orthogonal functions are used to capture the ionospheric features. Another recent approach was developed by Cheng et al. (2021) using a voxel-based technique, in which the fundamental difference is that voxel-based technique does not depend on polynomial coefficients (Yao et al., 2013). Using ∼300 GNSS stations, Cheng et al. (2021) have reconstructed the global ionosphere with a resolution of 2.5°, 5°, 50 km, and 1 hr in latitude, longitude, altitude, and time, respectively. This non-constrained reconstruction technique was applied with a parallel computing method in order to characterize the ionospheric anomalies of a storm event that occurred on September 7-8, 2017. As mentioned by the authors, the use of global data can pose problems to invert the 3D electron density in the ionosphere, requiring the development of new techniques that are capable of reproducing the ionosphere with higher temporal and spatial resolutions.
Built on previous work, it is reasonable to say that four of the main challenges to improve global scale tomography are related to the (a) enhancements of the spatial and temporal resolution of the 3D grid, (b) incorporation of TEC data with better global coverage, (c) development of realistic constraints, and (d) advances on the system geometry, that is incorporation of TEC measurements with distinct elevation angles. Based on such motivations, we have developed a technique which makes advances in the first three points. In this study, the objective is to provide advances on the development of a voxel-based CIT technique for mapping the global ionosphere accurately. In comparison to previous global-scale ionospheric tomography methods, we use a denser network of GNSS stations (∼2700), higher spatial and temporal resolution (2°, 2°, 20 km, and 15 min in latitude, longitude, altitude, and time, respectively) and we incorporate optimized constraints to stratify the ionosphere in several layers. Additionally, the present tomography evaluation is considering a set of accurate in-situ electron density measurements in order to present a broader evaluation of the capabilities of the global-scale tomography for the digestion of the ionosphere under storm conditions. Despite the computational demands, we have incorporated constraints in the global tomography since it is an important step to stabilize the solution at locations with only a few GNSS measurements. The CIT technique is applied to the challenging case study of the geomagnetic storm on March 17, 2015, which is also known as the 2015 St. Patrick's Day storm. The evaluation of the reconstructed 3D electron density has been performed at distinct altitudes based on electron density measurements of ionosondes, incoherent scatter radars (ISRs), and in-situ measurements from the Defense Meteorological Satellite Program (DMSP), Gravity Recovery and Climate Experiment (GRACE), and Swarm missions. In addition to the model-data comparison, we highlight the capabilities of the CIT technique to represent ionosphere disturbances 10.1029/2021SW002889 3 of 21 during the geomagnetic storm. The occurrence of strong enhancements and depletion of ionization with respect to quiet times (positive and negative ionosphere storms) has been described in many studies (e.g., Borries et al., 2016;Nava et al., 2016;Nayak et al., 2016;. In our study, we pay special attention to the equatorial ionization anomaly (EIA) in order to understand whether the developed tomography can be used to reflect relevant spatial and temporal responses of the ionosphere during storm conditions.

Tomography Formulation
The basic quantity required to perform CIT is expressed by the integral of the electron density along radio frequency signal paths crossing the ionosphere. CIT is, therefore, mainly applied to solve an inverse problem using total electron content (TEC) measurements. A major source of TEC measurements nowadays is the growing number of ground-based GNSS receivers, which provide a substantially large spatial and temporal coverage of the ionosphere in comparison to several remote sensing instruments. The calibration of the slant TEC (STEC) estimated in this work is based on the formulation provided by . As mentioned by , the STEC calibration is implemented in two steps. The first step estimates the Differential Code Bias based on IGS Vertical TEC (VTEC) maps from the Center for Orbit Determination in Europe (CODE) solutions. The second step estimates ambiguities using the so-called phase leveling, where the ambiguity terms are determined by the mean difference of the phase and code measurements over one arc of data with no cycle slips. We have chosen this method since it has presented similar accuracy to ionosphere-free solutions when applied to GNSS precise point positioning.
In order to estimate a set of unknown electron density values based on the GNSS TEC measurements, the ionosphere is broken down into a grid of voxels with a specific spatial resolution and coverage. The equation to solve for the unknown electron density x in each voxel is then given by Bust and Mitchell (2008) as: where y is the vector of STEC measurements, E A is composed by the path length of each GNSS signal inside the voxels, and ε is the system and measurement noises.
The ground-based GNSS measurements only provide a limited number of viewing angles of the raypath. Due to the almost vertical GNSS geometry of the ground-based observations, the direct solution of Equation 1 is ill-conditioned to describe the ionosphere using several layers. Ionospheric imaging by CIT is therefore a challenge due to this incomplete geometrical coverage that renders an ill-posed system (e.g., Mitchell & Spencer, 2003;Pryse et al., 1998;Ssessanga et al., 2021;Wen et al., 2008;Yao et al., 2014;Zheng et al., 2018). Hernández-Pajares et al. (1999), for instance, use Equation 1 with only two layers in order to solve the ill-conditioned problem. The incorporation of radio-occultation measurements is a viable option to improve the GNSS geometry (Yin & Mitchell, 2005) and thus the number of layers. The number of low Earth orbit (LEO) satellites is still sparse for high vertical and temporal resolutions. Additionally, the length of the unknown vector E x in Equation 1 is directly proportional to the number of voxels. Considering the ionospheric dimensions, a high-resolution 3D grid estimation can, therefore, cost an intense computational effort, especially on a global scale. In this regard, approximations and constraints on Equation 1 are required to gain computational efficiency and solve the ill-conditioned system with high vertical resolutions.
First CIT applications to represent multiple layers of the ionosphere were based on 2D regional geometries using iterative algorithms, such as the Algebraic Reconstruction Technique (ART) (Austen et al., 1988) and the Multiplicative ART (MART) (Raymund et al., 1990). In principle, ART and MART provide fast updates of the ionosphere without applying constraints on the estimation, which poses a problem since 3D iterative tomographic algorithms are affected by the almost vertical geometry provided by the GNSS measurements. Some studies have shown that constrained techniques can avoid this problem and improve the overall specification of electron density (e.g., Debao et al., 2015;Liu et al., 2010;Wen et al., 2010;Yao et al., 2020). The constraint dependence can add to the computational burden due to the additional matrix multiplications. In addition to the matrix multiplications, some constrained techniques require a large data set with many intersections of the GNSS signals. Therefore, on a global scale, the large number of unknowns x together with several data gaps precludes the possibility of using the most sophisticated constrained techniques, such as the proposed by Seemala et al. (2014). In order to gain both computational efficiency and accurate specification of electron density, the CIT technique in where γ is the relaxing parameter, y i is the ith TEC observation, ij E A is the path length inside the voxel j for signal i, and N is the total number of voxels. The relaxing parameter usually varies from 0 < γ < 1 in order to guarantee the convergence of the estimation.
In this study, we use the relaxing parameter from Prol et al. (2021), which helps to stabilize the ill-conditioned solution with no relevant additional computational efforts, since it does not incorporate a set of new matrix multiplications. This relaxing parameter is reproduced here as follows: where 0 j E x is the electron density given by the background model and 0 max E x is the maximum electron density along the GNSS signal. The constant 0.2 controls the iteration convergence.
Equations 2 and 3 can be solely used to provide good specifications of the electron density ionosphere at regions with a dense GNSS network. However, a major issue on the global tomography relates to the sparse GNSS receiver locations, leading to several data gaps. In this study, similar to Prol et al. (2019) and further validated by Zheng et al. (2020), we use VTEC data together with STEC to increase the spatial coverage of data. To improve the solution over the oceans, we have incorporated VTEC observations obtained by the Global Ionospheric Maps (GIMs) from CODE, which provides a smooth interpolation of the data gaps. The VTEC values obtained from GIMs are first applied to Equations 2 and 3 in order to provide an initial correction on the background. Then, both equations are used again to update the 3D grid using the STEC observations of each local GNSS station. Afterwards, the regularization process of Heise et al. (2002) is used to calculate the electron density in the remaining empty voxels.
The background 3D electron density is specified by the Neustrelitz Electron Density Model (NEDM), which well describes climatological patterns of the ionosphere and plasmasphere (Jakowski & Hoque, 2019). The spatial distribution is defined with a 2° × 2° × 20 km resolution for latitude, longitude, and altitude, respectively, for both the background and tomography. The vertical distribution ranges from 50 to 20,000 km, covering from the bottomside ionosphere until the GNSS orbit height in the plasmasphere. The time resolution is 15 min, that is, the electron density in each grid cell is defined to remain constant over a time period of 15 min. In this regard, the selected spatial and temporal resolution can be considered high for global scale tomography since there is a great challenge on the usage of higher resolutions due to the number of unknowns in the imaging system. Although this resolution may not be enough to capture very fine regional structures of the ionosphere, such as plasma bubbles, the defined resolution is sufficient to capture quite a few signatures and disturbances of the analyzed geomagnetic storm.

Data and Case Study
In each reconstruction step, the input STEC data for the tomography scheme are obtained from eight GNSS networks with a total of ∼2700 GNSS stations. Figure 1 shows a map of these GNSS stations. We used data from the Global Positioning System (GPS) collected by the continuously operating GNSS networks of Africa Our method is validated using ionosonde measurements from the Global Ionosphere Radio Observatory (GIRO) network (Reinisch & Galkin, 2011) and reference electron density of ISR observations at Millstone Hill. To ensure a sufficient quality of the ionosonde measurements, we adopt the automatically scaled data from the GIRO network, ignoring values with a Confidence Score (C-Score) equal to zero or unknown precision. In addition, we remove all data from the ISR measurements with the standard deviation of the electron temperature larger than 500 K and only profiles at altitudes lower than 600 km (empirically observed cut off for outliers) are kept. For the topside ionosphere, we compare results with in-situ electron density measurements of DMSP, GRACE and Swarm missions. The DMSP satellite orbit is materialized at an altitude of ∼830 km above the Earth, GRACE is flying at an approximate altitude of 500 km and Swarm-A orbit altitude is around 462 km. As reference for the DMSP data, we are relying on the scintillation meter (SM) of the satellite number F15, which provides plasma density by a simple faraday cup (Hairston, 2018). We have manually checked the occasional SM glitches and removed them when necessary. GRACE instruments are equipped with the K-band ranging (KBR) system, which allows the derivation of high-precision electron density measurements (Xiong et al., 2015(Xiong et al., , 2021. In other words, no preprocessing has been performed. For the Swarm mission, data from the Langmuir probes of satellite A is used to measure the electron densities (Knudsen et al., 2017).   Zakharenkova et al. (2016) for this event and Lu et al. (2020) added modeling results to these TID studies. Astafyeva et al. (2015) and other publications reported hemispheric asymmetries in the storm evolution. Also, a superfountain effect has been observed for the St. Patrick's Day storm 2015 for example by Nayak et al. (2016) and Fagundes et al. (2016). This means that the two crests of the EIA, which are usually at about 15° to the North and South of the equator at quiet times, increase in electron density and shift in the direction towards mid-latitudes. This effect results from the magnetosphere driven penetration electric fields imposing on the normal background zonal electric fields (e.g., Kelley, 1989). The following results will mainly focus on the description of changes in the EIA and fountain effect. Figure 2 shows the VTEC maps from tomography from March 16 to 18, 2015 (DOY 75-77) as well as the VTEC by NEDM during DOY 75. The main phase of the storm occurred on DOY 76, and DOY 77 was part of the recovery phase. DOY 75 is the quiet pre-storm day. The x-axis is the geographical latitude, the right y-axis gives the longitude, and the left y-axis gives the universal time. Figure 2 shows these VTEC maps from 0 to 22 UT for every 2 hr during each day. A more detailed view of the electron density profiles obtained by tomography are shown in Figure 3. The slices in terms of altitude and geographic latitude refer to the longitude section 60° W. Figure 3 demonstrates the capability of the global tomography to represent vertical structures of the ionosphere. The tomography electron density profiles are shown in the top three panels at DOYs 75, 76, and 77. The bottom two panels show DOY 75 subtracted from DOYs 76 and 77. Each column from left to right shows the variation from 17 to 18.75 hr UT, which reflect the period of the instances of the highest variabilities in the VTEC maps. The variability in electron density on DOYs 75 and 77 is relatively low compared to DOY 76. At 17 UT on DOY 76, an enhancement in electron density in the equatorial region is visible. And then the enhancement intensifies on either side of the equator until 18 UT. It is apparent that the electron density in the Southern hemisphere enhances more than the Northern hemisphere after 18 UT on DOY 76. The level of enhancement in DOY 76 is around 9 × 10 11 el/m 3 , as shown in the fourth panel, which is more pronounced in the Southern hemisphere. An electron density decrease is then observed in the mid-latitudes of the Northern hemisphere around 3 × 10 11 el/m 3 . In the recovery phase of the storm, a decrease in electron density is observed in the Southern hemisphere, reaching a minimum decrease of 10 12 el/m 3 , as shown in the last panel. Therefore, the results show that electron density can vary significantly within a few hours during disturbed conditions and remain relatively steady under quiet or recovery conditions; however, the magnitude of the electron density decrease can be more intense in the recovery phase than the increase in the main phase.

Overview of the Main Storm Features Observed by Tomography
To explore the reconstructed EIA distributions in more detail, Figure 4 shows the electron density distribution over distinct altitudes at 18 UT. As it can be seen, the EIA crests are strongest and clearly distinguishable in the altitude range 290-300 km. With increasing altitude, the maximum electron density gets closer to the geomagnetic equator. While the differences between the electron density distributions of DOYs 75, 76 and 77 are significant in the altitude range 290-300 km, it is not so evident at 450 km. During DOY 76, the EIA crests have larger distance from the geomagnetic equator than during DOY 75 and 77. We can observe asymmetries between the Southern and Northern crests, presenting a higher level of the plasma density in the Northern hemisphere. During the recovery phase in DOY 77, the EIA crests are suppressed at all altitudes.
To illustrate a specific analysis on the EIA crests, Figure 5 shows the EIA development over time together with the SYM/H index. The top panel shows the electron density at the peak height (NmF2) in terms of latitude by universal time. We fixed the local time at 15 hr to represent the EIA dynamics close to the highest ionization locations. Our investigation also presents a segmentation of the Southern (blue line) and Northern (red line) EIA crests. The segmentation was conducted by computing the latitudinal location of the maximum NmF2 values for the South and North directions. The middle panel gives a closer look of the maximum NmF2 latitudes by the solid lines. The dashed line stands for the differences between DOY 76 and 75, as well as between DOY 77 and 75. Thus, the dashed lines start on zero. During DOY 75, the maximum NmF2 shows the well-known two crests, which have different intensities and the location varies between 0° and 25° in latitude during the quiet day. Since the crest is parallel to the geomagnetic equator, the latitude of the crest presented in geographic coordinates . This suggests the penetration of the eastward electric field in the daytime ionosphere, which increases the force driving the vertical drift of the fountain effect and changes the EIA configuration. Despite the higher displacement of the NmF2 peak at the South, the whole structure of the Northern crest is affected up to about 50° in latitude (top panel). The storm-induced disturbance on EIA continues severely affecting the Southern crest until ∼2 hr UT of DOY 77. However, after two inflection points of the SYM/H index (9 UT and 23 UT in DOY 76), the anomaly starts to be significantly suppressed with a sharp decrease in the displacement of the crests, reflecting a recovery phase. Since there are no evident EIA crests, the two crests of the NmF2 are merged into one. However, it is evident the overall capabilities of using such three-dimensional reconstructions for a global analysis of the storm induced effects in the ionosphere. A still missing point to verify is related to the accuracy of the electron density estimation with the developed technique, which will be discussed in the next section.

Electron Density Comparison to Observations and Models
The electron density assessment compares the tomography results with NEDM and TIE-GCM. We chose NEDM (briefly described in Jakowski & Hoque, 2019) as a representative empirical model. The TIE-GCM was selected for this comparison since it is a well-established physics-based model (Hsu et al., 2021;Kodikara et al., 2021;Qian et al., 2014). TIE-GCM uses a finite differencing technique to discretize the numerical solutions for the conservation of mass, energy, and momentum to model the coupled ionosphere-thermosphere (Qian et al., 2014;Richmond et al., 1992). In this study, we use the TIE-GCM  perturbations to the advective and diffusive transport. In order to compare with different data products presented in this study, we linearly interpolate the TIE-GCM estimates to the data location, first along the geographic latitude and longitude and then a logarithmic interpolation to the geometric height.
In order to evaluate the accuracy of the tomography electron density estimations in the peak height, Figure 6 shows the critical frequency (foF2) values retrieved by ionosondes, tomography, background NEDM, and TIE-GCM during the week of the St. Patrick's Day storm. Figure 7 shows a similar graph to Figure 6, but for the difference between the model results and the ionosonde data. Here 40 ionosondes are used to show the local time variations. Each latitude row in Figures 6 and 7 corresponds to a distinct ionosonde. As can be seen, the ionosonde data present typical distributions, such as higher ionization in low-latitudes at the daytime and lower foF2 values in mid-latitudes at the nighttime. In general, tomography reproduced the foF2 distributions well including storm signatures during DOY 76, which is mainly shown by electron density enhancements during the disturbed geomagnetic activity, and after 17 March, during the storm recovery phase. Tomography correctly reproduces the foF2 variations observed by the ionosondes with a small underestimation. The background NEDM remains relatively steady and hardly follows the observed foF2 variability. TIE-GCM is also capable of reproducing the foF2 enhancement and decrease in the main and recover phases, respectively. Indeed, TIE-GCM is driven with geophysical indices that represent the disturbed conditions. The bias in tomography range from −2 to 2 MHz, while bias in TIE-GCM (NEDM) is around 4 (3) MHz, mainly during the storm phase. A complete statistical view on the differences between the methods is presented in Table 1.
In order to assess the tomography capabilities in reproducing the electron density distribution at altitudes above the peak height, in-situ electron density measurements from GRACE, Swarm and DMSP satellites are used as reference. For each instrument, we discriminate the in-situ measurements in two geometries, referred to as rising and setting phases. The rising phase is referred to when the satellite is flying from the South to North direction, while the setting phase is related from the North to South direction. In the presented experiments, the rising phase of GRACE and Swarm occur during the daytime, around 18 LT in the equatorial region. The rising phase of DMSP occurs during the nighttime, around 2 LT in the equatorial region. The setting phases in the equatorial region occur around 06 LT for GRACE and SWARM satellites, and around 14 LT for DMSP. Figures 8-10 show the in-situ measurements of GRACE, Swarm and DMSP, respectively, for the week of interest as well as the corresponding tomography, background and TIE-GCM estimations. TIE-GCM is not compared with DMSP due to the model's limited vertical domain. The electron density distributions are estimated by the models along the trajectories of the satellites.
Assessing Figures 8 and 9, we can see two crests in the daytime equatorial region at the GRACE and Swarm altitudes, likely associated with the influence of the fountain effect. During the DOY 76, the main phase of the storm has contributed to an enhancement of electron density. Additionally, due to the disturbed geomagnetic fields, we can see a slight change in the EIA structure in DOY 76 in comparison to other DOYs. Although the background NEDM has not shown such variabilities in the daytime, the tomography and TIE-GCM show plasma enhancements around the equator but slightly less than as measured by GRACE and Swarm. During nighttime, low electron densities were observed for GRACE and Swarm. As shown in Figures 8 and 9, the tomography correctly maps the Ne enhancements in the equatorial region during the setting phase of both GRACE and Swarm on DOYs 76 through 78. A general overestimation from tomography is observed at nighttime. Overall, TIE-GCM corresponds well with GRACE and Swarm measurements; however, in Figure 8, TIE-GCM does not seem to capture the plasma enhancement in the low latitudes on DOYs 78 and 79 during nighttime. Swarm and TIE-GCM also show such enhancement during this nighttime period in Figure 9 but the TIE-GCM results are slightly different in intensity and spatial distribution. In case of the DMSP results, we have observed a better representation of the electron density by tomography in comparison to the background NEDM, being capable to show more correlated distributions as well as storm signatures in the electron density, such as the electron density  In order to show the tomography performance to represent electron density profiles, specific cases of the ISR instrument at Millstone Hill were used as reference. Figure 11 shows four representative examples of the electron density profiles seen by tomography, the background NEDM model and TIE-GCM. The first three columns show the electron density profiles at 19 UT (10.4 LT), when the electron density was considerably disturbed during DOYs 76, 77 and 78. Unfortunately, no data was recorded before DOY 76. Results at DOY 79/19 UT are equal to DOY 78 at 19 UT, so that the fourth column is provided as an example to show the profile estimation at a distinct hour of DOY 79. We chose an hour in the nighttime (21.6 LT) to show an example of a profile at DOY 79. In most cases, the tomography results show lowest deviation from the ISR measurements. In contrast, the deviations between the models (NEDM and TIE-GCM) are varying a lot. During the recovery days of 77, 78 and 79, tomography reproduced well the ISR electron density profiles. During DOY 76, the tomography profile is also closer to the ISR data than the models, but attention has to be paid to significant deviations from ISR measurements in altitudes below 350 km. During the storm day (DOY 76), the ISR data profile is uplifted in comparison to modeling profiles. Neither the background NEDM and tomography were capable of showing such vertical distribution. However, in any case, tomography has presented a better agreement in the electron density profile estimations. With regard to TIE-GCM, we observe a strong overestimation in the example of DOYs 76 and 79. In a few cases of the recovery phase (DOY 77), we observe that TIE-GCM to show best consistency with ISR data, overpassing the tomography performance. The background was in general quite far from the observed profile.
A broader comparison of the electron density profiles of the ISR instrument at Millstone Hill is shown in Figure 12. We have neglected the data of DOY 76 in this example since very few electron density measurements were collected in DOY 76, in such a way that we could not observe the daily evolution of electron density distributions. Overall, we can confirm the analysis made in the examples of Figure 11, in which tomography has clearly outperformed the background NEDM and TIE-GCM. Table 1 shows a statistical overview of the tomography performance in comparison to the background and TIE-GCM in terms of the average of the error (model -data), standard-deviation of the error, and percentage of the Pearson correlation coefficient between model and data. The average error was computed using the experimental data set of all measurements available in the period of interest from the ionosondes, GRACE, Swarm, DMSP, and ISR. A positive average error stands for an overestimation, while negative averages indicate underestimations. Overall, we can see a significant improvement of tomography in comparison to the background for all the instruments, in the average, standard deviation and correlation. The percentage of improvement considering the root mean square error (RMSE) was 45%, 31%, 30%, 7% and 55%, for ionosonde, GRACE, Swarm, DMSP, and ISR, respectively. There is no relevant improvement in the DMSP assessment because NEDM was developed using DMSP in-situ measurements to calibrate the topside model. The physics-based TIE-GCM shows the lowest error for GRACE and Swarm. However, the standard deviation for the same cases is larger than for the tomography results. Based on the analysis of Figures 8 and 9, we assume that the correct representation of the electron density magnitude for GRACE and Swarm, with no relevant under or overestimations, was the main reason TIE-GCM presents lower average errors in comparison to tomography. However, the TIE-GCM error dispersion and correlation is larger.

Discussion
In this work, we present a global 3D reconstruction of the ionosphere electron density capable of reproducing challenging conditions during ionosphere storms. One of the most recent advances in this field has been presented by Cheng et al. (2021). Although the objective of our work equals Cheng et al. (2021), the implementation and capabilities of the CIT studies differ. The spatial resolution of the tomography result presented here is 2° × 2° × 20 km for latitude, longitude, and altitude respectively, covering from 50 to 20,000 km in altitude with a temporal resolution of 15 min. According to the authors' knowledge, this work presents the highest resolution of global CIT to date. Our method also incorporates more GNSS receivers and utilizes optimized constraints on the reconstruction. While the number of used ionosphere measurements is significantly lower in Cheng et al. (2021), they have the advantage of data coverage also over the oceans, because they are not only using combined multi-GNSS (GPS, GLONASS, BDS, and Galileo) observations, but also radio occultation (RO) measurements from Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) satellites. Although more data and a higher spatial and temporal resolution is applied in the study presented here, the computation of the CIT still completes within 20 to 30 min. This computational time is close to Cheng et al. (2021) who achieve 10 min computation time. There is still room for improvements in our processing time since we have used a simple computer with no parallel computing. The processing machine uses a socket of six intel i7 cores from the ninth generation with 16 Gb of memory RAM.   The capability of the applied CIT to reproduce ionosphere storm conditions is assessed based on the example of the St. Patrick's Day storm 2015 and discussed here with respect to published ionosphere storm studies for this event. For example, Nava et al. (2016) describe that at the beginning of the storm, one or two periods of enhanced ionization with respect to quiet times (positive storm effects) were observed depending on the longitudinal zones. The authors showed that the American region exhibits the most remarkable increase in VTEC. Furthermore, a strong decrease in ionization (negative storm) is observed at all longitudes after the main phase of the storm, lasting several days (e.g., Fagundes et al., 2016;Nava et al., 2016). This decrease was reported to be the largest in the Asian sector. The remarkable ionization increase in the American region and a strong decrease in the Asian sector are also confirmed in the present work at Figure 2 at 18 hr UT in DOY 76 (America) and at 6 UT in DOY 77 (Asia).
The CIT results presented here indicate a hemispheric asymmetry of the ionosphere storm perturbations. In Figures 3 and 4, a larger plasma density enhancement has been shown in the Northern hemisphere compared to the Southern hemisphere, especially after 18 UT. This effect, with dependency on the longitudinal sector, has been discussed by Astafyeva et al. (2015). They explain this effect with composition changes and partly with the different non-dipolar portions of the geomagnetic field in the two hemispheres as well as with IMF by component variations. In addition, Liu et al. (2016) reported a hemispheric asymmetry also in the negative storm conditions close to the storm enhanced density on DOY 76, where the Northern hemisphere negative storm was stronger than in the Southern hemisphere. This effect is visible in Figure 3 (second panel row from the bottom).
In Section 4.1, special attention has been paid to the presentation of the EIA. The EIA is developed because the electric and magnetic fields are, in general, perpendicular to each other over the geomagnetic equator in the ionosphere. The resulting force in the daytime drives the electron density to higher altitudes with a drift velocity E × B, where E and B are the electric and geomagnetic field vectors (Anderson, 1981). As described in , the peak height and the whole profile distribution are proportionally uplifted depending on the intensity of the zonal electric fields at the equator. At higher altitudes, the uplifted plasma diffuses along with the geomagnetic field lines as a consequence of the pressure and gravity forces. The plasma diffusion along the geomagnetic field lines creates two crests at about 15° to the North and South. The two crests formed by the so-called fountain effect can be clearly observed in Figures 2 and 3. The EIA crests in DOY 75 represents a natural development during quiet times, showing daily variabilities due to the expected longitudinal asymmetries, which can be related to daily variations of the solar ionization, longitudinal differences of the zonal electric field and distinct ionization level due to the geomagnetic field configuration (Kil et al., 2011). In addition, we have observed that the EIA crests at DOY 76 after 14 UT have a higher ionization level and are more far apart than during the quiet and recovery days (Figures 3-5). These features also agree well with the observations of the superfountain (e.g., Fagundes et al., 2016;Nayak et al., 2016). Based on the electron density reconstructions, we could observe a displacement of the EIA crests induced by the storm by 9°-20° in latitude when compared to a quiet day. This effect relates to the penetration of an eastward electric field in the daytime ionosphere, as described by for example Hairston et al. (2016). This electric field increases the force driving the vertical drift of the fountain effect and changes the EIA configuration. Also, several other studies of the St. Patrick's Day storm 2015 (e.g., Astafyeva et al., 2015;Fagundes et al., 2016;Huba et al., 2017;Liu et al., 2016;Nayak et al., 2016;Yue et al., 2016) suggest the involvement of eastward penetration electric fields not only in the EIA enhancement but also more general in the electron density enhancements in the ionosphere during the main storm phase.
A significant decrease in the EIA intensity has been shown in Figure 5, starting at 2 UT on DOY 77 a few hours after the storm recovery started. A complete reversal in the longitudinal structure/pattern of EIA during the recovery phase of the storm as compared to the quiet day has been reported by Yadav et al. (2016). Also, in mid-latitudes, a strong negative storm is apparently visible in Figures 2 and 5. As shown by Yadav et al. (2016) and Borries et al. (2016), the suppression of EIA can be associated with the Disturbance Dynamo Electric Field (DDEF) and westward PPEF due to a northward turning point of the Bz component in the interplanetary magnetic field. A resulting counter equatorial electrojet (CEJ) imposes the westward electric field over the daytime ionosphere and generates a depression of TEC and EIA around 6-12 UT. In addition, a short partial recovery phase has been reported during 9-12 UT (e.g., Borries et al., 2016;Yadav et al., 2016), which is associated with reduced TEC enhancements with a certain delay to this period. A slight decrease in the intensity of the EIA on DOY 76 between 12 and 14 UT, visible in Figure 5, can be associated to this partial recovery phase. During the partial recovery phase a westward PPEF due to a northward turning point of the Bz component in the interplanetary magnetic field is suggested in Borries et al. (2016) and Yadav et al. (2016). A resulting CEJ imposes the westward electric field over the daytime ionosphere and generates a depression of TEC and EIA.
A significantly higher F2-layer peak height (hmF2) has been detected by the Millstone Hill ISR (Figure 11). Liu et al. (2016), who analyzed the Millstone Hill data in more detail, point out that while TEC was enhanced, a decrease in NmF2 and an increase in hmF2 took place. The plasma was distributed over a much larger range of altitude. An increase of the density of molecular neutrals, resulting in enhanced recombination, is responsible for the observed decrease in NmF2. These hmF2 changes during storm conditions can be caused by zonal electric fields and meridional thermosphere winds. Zhang, Erickson, et al. (2017) describe the contribution of the PPEF to the upward and poleward plasma transport for the Millstone Hill observatory in the association with the storm enhanced density during the St. Patrick's Day storm 2015. Also, Ratovsky et al. (2021) identify the zonal electric field as the main contributor to the hmF2 increase over the Irkutsk station during the St. Patrick's Days storm 2015. Figure 11 shows a limited capability of the tomography result to reproduce the hmF2 increase. This is mainly associated with the limitations of the background model and the complex electrodynamic processes generating this effect. Indeed, the developed tomography is highly dependent on the hmF2 values of the background. The background NEDM hmF2 was not capable of showing such vertical distribution, which is reflected in the tomography performance.
An extensive model-data comparison has been presented in Section 4.2, using ionosonde, incoherent scatter radar, Swarm, GRACE and DMSP electron density data. It has been shown that the tomography approach generates very good estimations of the 3D electron density and can provide improvements in comparison to NEDM and the numerical model TIE-GCM (lower RMSE and better correlation coefficients than the numerical model). The larger standard deviation of TIE-GCM is justified since an overall larger structure of the EIA was reproduced by TIE-GCM during the satellite's rising phase. After the storm onset, during DOY 76, higher electron densities are reproduced by TIE-GCM above ±40° in latitude, while the reference data and tomography are more concentrated to the low-latitudes. In terms of RMSE, tomography outperforms TIE-GCM by 57%, 30%, 11%, and 51% for ionosonde, GRACE, Swarm, and ISR data, respectively. However, there is still some potential for improvement, since some error estimates in Table 1 (bias with respect to GRACE and Swarm electron densities) are smaller for TIE-GCM than for the CIT result. Also, an overestimation of tomography electron densities above the peak height (measured with GRACE and Swarm, Figures 8 and 9) has been identified, which is expected to be a consequence of the even larger overestimation of the background NEDM. This shows that even after the tomographic update, there is still an important dependence on the background to regularize the solution. In greater altitudes, the error of the tomography results become rather small (Figures 10-12 and Table 1). Still, there is room for improvement at greater altitudes since the electron density magnitude during the DMSP setting phase is underestimated by the tomography. For instance, in Figure 10, tomography was not capable of reproducing the DMSP electron density storm enhancement over a wide range of latitudes in DOY 77. Analyzing the data set, it is plausible that this discrepancy may be due to the lack of GNSS data over the oceans. Indeed, over the oceans, the tomography performance mainly relies on the background and the CODE VTEC maps. A dedicated study may still be necessary to understand the performance of the tomography in these regions. Incorporation of radio-occultation-derived electron density in the tomography may help to improve the specification over larger areas. As discussed by Cheng et al. (2021), GNSS radio-occultation techniques can overcome the coverage limitation in ocean regions. Previous results, such as Prol et al. (2021), have presented a high correlation between DMSP electron densities and tomography using podTEC data, that is, topside TEC measurements by GNSS receivers onboard LEO satellites. Therefore, additional TEC data observed by the precise orbit determination could also improve the tomography estimations to fill the data gaps.
It remains to assess the capabilities of the developed tomography for distinct ionospheric responses to the storm event other than the ones discussed above. The large-scale effects like the ionosphere trough, BPD, large-scale TIDs and the TOI, reported for this event (c.f. Section 3), should be visible in the presented global tomography results. The amplitudes of the large-scale TIDs are rather small compared to other storm perturbations (e.g., trough and BPD) and unlikely to be present in our tomography. Still, Bolmgren et al. (2020) demonstrated that a 3D reconstruction of large-scale TID signatures is possible with tomography. A detailed discussion of all these large-scale storm effects is beyond the ambit of this study. Additional studies are necessary to assess the capabilities of the presented tomography with respect to these large-scale storm effects. There exists a high interest, to also study finer ionosphere perturbation structures with CIT reconstructions. For example, plasma bubbles are shown to be visible in regional high-resolution tomography results, as shown by  and medium-scale TIDs have been reconstructed with tomography by Chen et al. (2016). However, such effects cannot be discussed with the presented global tomography result, since the spatial resolution is not sufficient. In future work, it may be beneficial to investigate the performance and efficiency of the global tomography in comparison to regional reconstructions using dense GNSS networks, such as the regional methods applied over Japan (Saito et al., 2017). Data assimilation approaches applied on numerical models (e.g., Hsu et al., 2021;Kodikara et al., 2021) would also be an interesting option for comparison with the global-scale tomography results. Similarly, the high accuracy tomography outputs could be used to drive the numerical models in forecasting applications.

Conclusions
This study demonstrated improved capabilities in quality and performance of tomographic techniques to obtain reliable electron density profiles during the analyzed geomagnetic storm using large data sets of TEC estimates derived from ground-based GNSS receivers. Improved temporal and spatial resolution for global electron density reconstructions can be achieved in the reconstruction. To test the method, a global network of ionosondes, as well as ISR observations and in-situ measurements of DMSP, GRACE, and Swarm are used. The results indicate the inversion scheme as a capable method to improve the accuracy of the background specifications and common models. Additionally, many physical mechanisms of the observed ionospheric disturbances have been discussed, showing the capabilities of the CIT technique to reconstruct the EIA and electron density signatures during disturbed conditions. The electron density reconstructions show well the displacement of the EIA crests compared to a quiet day, which is associated with a PPEF effect during the storm event. We also observed a significant inhibition of the EIA formation in the recovery phase, associated with a DDEF in the daytime. Also positive and negative storm effects and hemispheric differences are captured with the CIT reconstruction. Changes in the F2-layer peak height are not yet captured sufficiently. Future research is still required to improve the temporal and spatial resolutions of the reconstructions to detect finer structures, such as equatorial plasma bubbles and traveling ionospheric disturbances. Further analysis should be conducted to evaluate the global tomography during longer periods in order to analyze the model accuracy in different seasons and solar activities. However, the overall capabilities of using such three-dimensional reconstructions for a global analysis of the storm induced effects in the ionosphere are evident.