Title: Assessing the performance of a Northeast Asia Japan-cantered 3-D ionosphere specification technique during the 2015 St. Patrick’s day solar storm

This paper demonstrates and assesses the capability of the advanced three-dimensional (3-D) ionosphere tomography technique, during severe conditions. The study area is northeast Asia and quasi-Japan-centred. Reconstructions are based on Total electron content data from a dense ground-based global navigation satellite system receiver network and parameters from operational ionosondes. We used observations from ionosondes, Swarm satellites and radio occultation (RO) to assess the 3-D picture. Specifically, we focus on St. Patrick’s day solar storm (17-19 March 2015), the most intense in solar cycle 24. During this event, the energy ingested into the ionosphere resulted in Dst and Kp and reaching values ~-223 nT and 8, respectively, and the region of interest, the East Asian sector, was characterized by a ~60% reduction in electron densities. Results show that the reconstructed densities follow the physical dynamics previously discussed in earlier publications about storm events. Moreover, even when ionosonde data were not available, the technique could still provide a consistent picture of the ionosphere vertical structure. Furthermore, analyses show that there is a profound agreement between the RO profiles/in-situ densities and the reconstructions. Therefore, the technique is a potential candidate for applications that are sensitive to ionospheric corrections.


Introduction
In trans-ionospheric radio communication, central to the precise performance is the ability to understand, characterize, and forecast the distribution of an irregular transiently changing dispersive terrestrial plasma (see for example Jakowski et al. 2011;Cannon et al. 2013;Kelly et al. 2014;Eastwood et al. 2017;Yasyukevich et al. 2018;Rovira-Garcia et al. 2020). Despite the complexity of this subject, traditionally in the past few years a simplistic horizontal distribution has been suffice, to a certain level of accuracy, to correct for the ionospheric effects on radio signals (e.g., Saito et al. 1998;Jakowski et al. 2012Jakowski et al. , 2011Ohashi et al. 2013). The horizontal distribution is mainly deducible from total electron content (TEC), a derivable quantity from abundantly available satellite-ground radio transmissions (Otsuka et al. 2002;Ma and Maruyama, 2003). However, as technological systems become more integrated, the demand for high precision on trans-ionospheric signals has increased. Thus, the need for a detailed three-dimensional (3-D) ionosphere picture. Computerized ionosphere tomography (CIT) is a way to make 3-D (or 2-D) distribution, i.e., an analysis technique to give height dimension from the organized multi-point measurements of TEC. Such advanced CIT is now called ionospheric imaging and is vital in radio communication since it provides information on the four-dimensional (time and space) evolution of the electron density structure. An excellent review and history on some of these techniques is found in Bust and Mitchell (2008).
The fidelity of the imaging technique depends on the accuracy of the measurements as well as the temporal and spatial distribution of the data used in the reconstruction.
Unfortunately, in most cases, measurements are imperfect and irregular in space and time. Besides, geometry constraints limit the observation information that can be included in the analysis. In return, imaging or tomography analysis (inverse problems) are generally under-determined and ill-posed (Yeh and Raymund, 1991;Raymund et al. 1994;Bust and Mitchell 2008). To have a consistent 3-D picture, a prior information from other sources such as models is needed. The accuracy of this information and the subject of how to best combine it with measurements has spawned an explosion of research and suggestions (e.g., Fremouw et al. 1992;Howe et al.1998;Rius et al. 1998;Bust et al. 2000;Tsai et al. 2002;Mitchell and Spencer, 2003;Ma et al. 2005;Schmidt et al. 2008;Okoh et al. 2010;Seemala et al. 2014;Ssessanga et al. 2015;Chen et al. 2016;Saito et al. 2017).
To maximize the advantages of imaging techniques, reconstructions are performed on a regional basis; wherein observations are abundant, and computation is tractable under high spatial grid resolution settings. Over the East-Asian sector, an imaging technique has been developed based on a regional network of dense GNSS (global navigation satellite system) receivers and ionosondes (Saito et al. 2017(Saito et al. , 2019Ssessanga et al. 2021). Figure 1 showcases the network and the horizontal resolution settings. The vertical dimension is not illustrated in the picture but extends to a GNSS altitude of 25,000 km (see Saito et al. 2017). The technique is a two-stage algorithm: First GNSS STEC (slant total electron content) are used to reconstruct the ionosphere electron density field based on constrained least-squares CIT. Second, to improve the exactness of the reconstructed picture, CIT densities are considered as background in an ionosonde data assimilation technique which assumes that ionosphere electron densities are better described by a log-normal distribution. Saito et al. (2017Saito et al. ( , 2019 and Ssessanga et al. 2021, have provided the mathematical derivation of this algorithm together with some preliminary results showing that the algorithm performs better than a stand-alone CIT technique. However, the authors did not specifically provide a performance analysis of the algorithm during geomagnetically disturbed conditions. Moreover, a wide spectrum of energy injected into the ionosphere during geomagnetic storms leads to a chaotic and non-linear ionosphere both in space and time. Although the main phase of a storm may last for a short period, the aftermath effects can last for many days before recovery (Prölss (1995) and Buonsanto (1999)).
Under such deviant conditions, most climate (such as international reference ionosphere(IRI)) or physics-based models used in forecasting and nowcasting ionospheric behaviour are rendered insufficient. Therefore, for applications that utilize trans-ionospheric signals, a need for a data-sensitive near-real-time observation driven algorithm is paramount. The goal of this paper is to show that the technique suggested by Ssessanga et al. 2021, can detect electron density dynamics even during severe storm conditions, and thus a potential candidate for near-real-time ionospheric corrections/applications. All hyperparameters are maintained the same as published in Saito et al. (2017Saito et al. ( , 2019 and Ssessanga et al. 2021. Section 2 gives a short mathematical description of the technique. Section 3 presents the reconstructionsmeasurement(ionosonde) comparisons illustrating the ionosphere electron density dynamics in connection to the ingested energy during the geomagnetic storm. It also analyses in detail the reconstructed densities in comparison to independent observations from Swarm satellites and radio occultation. Section 4 is the summary.

A brief description of the algorithm
Consider a required regional 3-D electron density field state ( n  ) and a column vector where A is an operator that maps the electron density field state into observation space, W is a weight matrix based on prior information and restrains the derived electron density from exceeding a certain value (Saito et al. (2017(Saito et al. ( , 2019). Elements superscripted with bc represent boundary conditions obtained from the NeQuick model (Nava et al. 2008) and where λ (> 0) is a regularization parameter. In the nutshell, the tomography analysis is based on the dense GPS-TEC data with model empirical data normalized to vertical profile peak densities and only used to constraint the solution (see Seemala et al. 2014). We developed a software system on a desktop PC, and realized the realtime tomography analysis from 1-second GEONET data from 200 stations over Japan. Performance of this tomography (Saito et al. (2017(Saito et al. ( , 2019) was generally good, but limited when the F region peak densities (NmF2) appeared below the ~270 km altitude range in October and November. Ssessanga et al. (2021) developed an advanced version of this tomography by incorporating two more features that provide a more realistic representation of the ionosphere. One is the inclusion of the ionosonde parameters h'F and foF2 from the nearby operational ionosondes. The other is the evaluation of data as logarithmic of the targeted densities; whereby all the elements in random vector n  are assumed to be drawn from log-normal distributions, and then apply Gaussian statics by taking the natural logarithm of each element ( ). That is A new cost function is then formulated as follows where R and B are data and background errors covariance matrices, and b X  is a vector of background log densities from the CIT solution.
The analysis log density a X  that minimizes (4), is derived as Where 1 − j A is a Jacobian evaluated at iteration j-1. At optimal, the real densities are computed as . 10 Colours green, blue and red represent the original tomography, improved/modified version, and ionosonde densities, respectively. In each subplot, the 270 km baseline is illustrated with a purple horizontal dashed line. Contrary to the original tomography, the modified version displays a distinct ability to vividly follow the measured data both vertically and horizontally. This is effective adaptiveness, and in the sequel, we further examine the penitential of the modified tomography under extremely chaotic ionosphere conditions. In fact, this work is part of a pre-analysis of the modified tomography version before advancing to a near-real-time on-line version (https://www.enri.go.jp/cnspub/tomo3/plotting.html).

Analysis of reconstructions during geomagnetic storm conditions
The geomagnetic solar storm analysed occurred on St. Patrick's day (day of year (DOY) 076-077) in 2015 and was among the most intense storms in the 24 th solar cycle (Astafyeva et al. 2015). We reconstruct the ionosphere vertical structure during this period (at a resolution of 15 minutes) and briefly analyse the variations in reference to the ionosphere disturbance indicator, Dst index. Specifically, the vertical structure is analysed at ionosonde locations for easy comparison with ground ionosonde observations. In-situ densities from the Swarm constellation and radio occultation density profiles are also utilized to assess the reconstructed F-region topside densities. Throughout the analysis, international reference ionosphere (IRI-2016, Bilitza et al. 2017) model densities are also presented as a form of reference.

Storm chronological highlights
A coronal mass ejection (CME) was observed erupting between ~ 00:30 UT -00:04 UT on DOY 074 and was predicted to encounter the Earth's magnetosphere on DOY 076. Figure 3 illustrates how the abundant energy injected into the magnetosphere-ionosphere system perturbed the terrestrial geomagnetic field leading to a Dst (Kp index) minimum (maximum) recording of approximately − 223 nT (8).
Kp as blue is scaled on the right and Dst as black is scaled on the left. Dashed vertical lines indicate the main events of the storm. Both indices can be obtained from the World Data Center for Geomagnetism, Kyoto http://wdc.kugi.kyoto-u.ac.jp or https://omniweb.gsfc.nasa.gov/form/dx1.html.
Before DOY 076, the geomagnetic indices Kp and Dst were relatively stable with magnitudes below 4 and 50 nT, respectively. On DOY 076, at 04:45 UT (sudden storm commencement, SSC), the CME hit the Earth, leading to an increase in Dst. On the same day, at ~07:30 UT the main phase of the storm commenced and the Dst decreased to a local minimum of ~ -80 nT. During this period the interplanetary magnetic field (IMF) Bz component was reported to have turned Southward (Cherniak et al. 2015). From about 9:30 UT to 12:20 UT there was a short-lived recovery in the Dst index to ~-50 nT; expected to be due to the IMF Bz component turning Northward (Astafyeva et al. 2015). From 12:20 UT the Dst continued with a gradually decreasing trend reaching a global minimum at ~22:00 UT on DOY 076.
The recovery phase continued through the following days with the Kp index ranging below 6.

Analysis of reconstructions at ionosonde locations
In Figure  In each image plate: for clarity, the altitude range is limited to the region of uttermost interest (F-region ~200 -600 km). Superimposed is the Dst index (ordinate axis on the right) in order to track the response of the ionosphere during the different storm dynamics. As indicated earlier, dashed vertical lines indicate the main events of the storm. The horizontal dashed line at 300 km is included as a baseline to vividly track the uplift of the plasma density. The time-altitude variation of the peak-density is represented by the purple solid-asterisk line. In the reconstructed, and ionosonde images, white spaces indicate periods when GPS-ground receiver ray links and ionosonde observations did not exist or intersect that particular location. At a particular time-stamp, if there exist any data gaps within the 200-500 km altitude range, the peak density is not determined.
Throughout the selected analysis period, at all stations, the IRI model only shows the general diurnal variation of the F region without detailed features. This is expected since IRI is a climatological model that represents monthly averages of terrestrial plasma densities/frequencies (Bilitza et al. 2017). Consequently, the IRI model performance is usually found lacking when ionosphere underlying driving mechanisms are deviant from the general trend.
Contrary to IRI, reconstructed results reflect changes in the plasma frequency following variations in Dst index as detailed: The onset of the storm occurred when the East Asian sector was on the day-side. Compared to reconstructed densities on the day before the storm (DOY 075); nearly just after the SSC, at OK426 station that is located geomagnetically in low latitudes, there is a slight noticeable instantaneous uplift (marked with a grey arrow) of the peak densities to altitudes above the 300 km baseline. The uplift reduces in amplitude towards the high latitudes (not noticeable at TO536) and is most probably due to geomagnetic storm induced eastward prompt penetration electric fields (PPEF), which end up modifying the equatorial and lowlatitude ionospheric phenomenology (Maruyama et al. 2004).
Approximately two hours after the SSC, there was a slight increase in F region plasma densities (~280 -340 km) particularly at stations located near the equatorial anomaly (OK426 and JJ433). In addition, there was a well defined wavy modulation of the peak density altitude at mid-latitude stations (JJ433, TO536 and IC437). The ~ 2-hour delay falls within the time frame the traveling atmospheric disturbance (TAD) created at high latitudes during geomagnetic storms, would take to propagate to the mid and low latitudes (for example see the works of Shiokawa et al. 2003 and references therein). Therefore, low latitude ionisation and wavy modulation of the mid-latitude F region could be related to an equator-wards propagating surge and plasma density enhancements in equatorial ionization anomaly (EIA) crests, following day-side PPEF.
By the time the East Asian sector enters the night side (~10:30 UT), there is a short-lived recovery in the plasma frequency between 10:30 -12:30 UT. Thereafter (~12:30 -24:00 UT), specifically at stations located in the mid-latitudes (JJ433, TO536 and IC437) and towards high latitudes (WK546), the plasma is uplifted to high altitudes (~500 km), gradually falls to ~320 km, rises again to about ~380 km and finally settles to an average altitude of ~300 km. Through this whole process, at all stations, the initial daytime plasma enhancements decrease gradually (to values nearly below 4 MHz), following the Dst index which hits a minimum at about 22:00 UT.
The decrease in F-region densities accompanied by descent in altitude of the peak densities could be attributed to an equator-ward wind surge and the night-side westward PPEF which has opposite effects (reduction in EIA ionisation) in reference to day-side PPEF.
On the day after the storm (DOY 077), at all stations, the plasma density (frequency) remained nearly ~60 % below the values observed on the day before the storm.
Interesting to note through this period, is the total or partial absence of data at all ionosonde stations. From the reconstructions, it is clear that these data gaps were due to a significant decrease in the F region plasma frequency (< 4 MHz) such that the ionosonde could not detect the reflected echoes. This is a well-known/typical negative ionospheric response following a major geomagnetic storm (e.g., Fuller-Rowell et al. 1994;Prölss, 1995;Tsagouri et al. 2000). A possible explanation for the intense plasma reduction is that the sub-storm activity on DOY 076 generated a composition bulge that entered the Asian sector, and persisted through the next day (DOY 077). In fact, Astafyeva et al. 2015 andNava et al. 2016, analysed the thermospheric column integrated O/N2 ratio changes measured by the GUVI (Global Ultraviolet Imager) instrument onboard the TIMED (Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics) satellite during the same geomagnetic storm (DOY 076-078, 2015), and found that the Asian sector had significant composition changes that could have led to the negative ionosphere response.
Latitudinal wise, stations towards the high latitudes (WK546) exhibit the highest reduction in plasma density. This enormous reduction in F plasma densities led to an increase in slab thickness, which is seen to increase polewards (extending >300 km in altitude), and most pronounced during the evening of DOY 076 and the early hours of DOY 077. Indeed, different studies have already shown that the slab thickness is dependent on diurnal, annual, latitudinal and storm-time variations (see Stankov and Warnant, 2009, and references therein). In the longitudinal perspective, stations TO536 and IC4337 that are almost on the same geomagnetic latitude (the difference is ~1°) show that the negative storm effect was greater on the Eastern side.
On the third day after the storm (DOY 078), except for OK426, all stations maintain a slow (< 6 MHz) gradual recovery towards the normal plasma frequency values observed during the quiet period. Nava et al. 2016 analysed the same storm and illustrated that the recovery process lasted more than 7 days. At OK426, reconstructions indicate that the low-latitudes had an undulating peak, accompanied by reinforced ionisation covering altitudes ~ 200-420 km. Surprisingly, these effects do not extend to mid-latitudes (see reconstructions at JJ433, TO536 and IC437). Nava et al. 2016 also generated regional TEC maps corresponding to periods, before, during and after the St. Patrick's day storm. Over the Asian sector, on the day after the storm, ionization was confined at low latitudes near the equator, and it almost disappeared at middle latitudes. This ambiguous or peculiar behaviour may be related to post-storm effects that might require further investigation (also see Kutiev et al. 2006).

GPS and ionosonde:
The significance of complementing GPS with ionosonde data is well illustrated at most stations; we can reconstruct the full extent of F region plasma structure beyond the bottomside limitations of the ionosondes. Moreover, as noted earlier, at points where the ionosondes failed to record any echoes due to the low plasma electron content, our technique was still able to provide reasonable plasma content estimations. This is a crucial result particularly for applications that utilize trans-ionospheric HF (high frequency) signals.

Comparison with Swarm densities
In-situ densities offer a good test of accuracy since they represent a specific point in space (ionosphere) at a particular time. The Swarm constellation consists of three identical polar-orbiting satellites that fly at two different altitudes ≤ 460 km (Alpha (A) and Charlie (C)) and ≤ 530 km (Bravo (B)). Each satellite has a Langmuir Probe that facilitates the measurement of in-situ densities. Part of the Swarm density results presented here are already presented in Ssessanga et al. 2021 (as preliminary results) to showcase the performance of the technique. We will add more analysis plots to cover the morning and evening sectors of the analysed period. The data used in the analysis are Level 1b electron density measurements at a 2 Hz rate, accessible at ftp://swarm-diss.eo.esa.int. For a further review on Swarm data see, for example, Olsen et al. (2013).
In Figure 5 In all sub-plots, the agreement between the reconstructed and in-situ densities is generally good. IRI on the other side exhibits a poor estimation of the densities. Also noticeable is that due to differences in orbital altitudes, there is a slight difference in the level of electron density distributions observed by satellites (A, C) and B. The difference between in-situ and reconstructed densities is on average less than 0.2 x 10 12 el/m 3 . This result offers confidence in the reconstructed topside F-region density structure, which is not attainable while using ground-based instruments such as ionosondes. A point of concern could arise from the poor density reconstructions at the awake of profiles in time range 9:00-11:00 UT (corresponding to satellites A and C on DOY 076 and 077). During this period Swarm A/C traces fall within or near the equatorial latitudes. Ssessanga et al. 2021, ascribed this to a poor grid specification over these latitudes. That is to say, the equatorial ionization anomaly region exhibits steep latitudinal densities gradients, yet, the currently utilized grid assumes a coarse horizontal spacing (5° × 5°) over this region (refer to Figure 1 and Saito et al. (2017; 2019)), hence reconstructions would not reflect the true ionospheric state and dynamics.
If we follow the traces of Swarm A and C, we observe that on DOY 076 when the storm commences (~10:00-12:00 UT) the topside densities range 0.5 ~ 2 × 10 12 el/ m 3 . On the day that follows, during the same period, when the negative storm effects are dominant, we observe that the densities remain nearly below 0.25 × 10 12 el/m 3 (~50% of the densities on DOY 076). The reduction in plasma densities is also observed at Swarm B altitudes and is maintained throughout the late hours of DOY 077. This result is consistent with the density reductions that were observed in the middle and bottom panels of Figure 5.
On DOY 078, when the F-region topside densities start returning to normal, the reconstructed densities relatively track well the in-situ measurements, but with a better performance at Swarm B altitudes. The discrepancy in performance may be related to the altitude level; Swarm (A and C) orbit at a lower altitude (≤ 460 km) than B, and at such altitudes the plasma densities are more likely to be influenced by different non-linear dynamical forces during the storm period, hence making the reconstruction acute. Certainly, the local time (~UT + 9) difference in Swarm B and A/C observations (different ionisation levels) might also influence the results.
Nevertheless, the reconstructed densities still outperform the IRI model estimations.

Comparison with radio occultation (RO)
In Figure 6, a set of RO electron density profiles (red) from COSMIC (Constellation Observing System for Meteorology, Climate, and Ionosphere; orbital altitude ~800 km) and GRACE (Gravity Recovery and Climate Experiment; orbital altitude ~490 km) constellations are plotted together with profiles from assimilated tomography (blue) and IRI (black) during the analysed period. RO data are accessible in level 2 format at https://cdaac-www.cosmic.ucar.edu/cdaac/. Comparison is limited to reconstructions away from the low latitudes and above the 200 km altitude mark, where RO density profiles are expected to have the best accuracy. That is to say, RO density profiles are an inversion of RO total electron content (ROTEC) using the so-called inverse Abel transform that assumes spherical symmetry in the ionosphere: For the most part, profiles from Abel transform have good accuracy, with the exception of the E-region (where rays have asymmetric contributions from the F region portions of the rays), and low latitudes (where large density gradients exist, Garcia-Fernandez et al. 2003;Wu et al. 2009;Yue et al. 2010). Rather than comparing the RO density profiles to a specific vertical profile within the grid, the comparison is performed at the location of tangent points (which contribute the most density along the RO ray path) shown on the maps in the top right corner of each subplot. Purple and green represent GRACE and COSMIC, respectively. Gaps in the profiles indicate instances when assimilated tomography reconstructions were not available. Fortunately, during the analysed period the GRACE constellation covered the ~ 200-400 km range, whereas COSMIC mostly covered the topside ~ 400-700 km.
This snip view gives us a chance to nearly observe and analyse the capabilities of the assimilated tomography in specifying the electron density field both in time and space (horizontally and vertically).
On the DOY 075 (top left corner subplot) before the storm, both assimilated tomography and IRI have a good estimation of the topside structure. However, on days that follow, assimilated tomography and RO profiles show a dramatic decrease in electron densities consequent to the geomagnetic storm. By contrast, the IRI model maintains a high-density output, with the largest deviation from the truth in the 200-400 km altitude range. This result is important because it gives a sample of what might be expected in applications that use models to correct for ionospheric effects during severe conditions. Despite the variability, reconstructed profiles on average adequately track the RO densities at all height ranges. Surprisingly, a combination of ionosonde densities and ground-based GNSS TEC is adequate to reconstruct a reliable topside structure (~400-700 km). Nonetheless, the noisy structure of some of the profiles motivates our suggestion in future to integrate RO data (both TEC and electron density profiles) in the analysis to further constrain the vertical structure.
Results presented in this study can be placed in context with previous performance analyses of the tomography technique during the early stages of development. In this sense, our work is a complement to studies by Seemala et al. 2014;Saito et al. 2017Saito et al. , 2019Ssessanga et al. 2021, who have already analysed the tomography reconstructions during the quiet period and found the technique to have good fidelity. Therefore, the technique seems fairly robust in handling different ionospheric conditions.

Summary
The material presented here has offered an insight into the ability of the tomography ionosonde data assimilated technique to characterize ionosphere plasma densities even under severe ionosphere conditions. The results clearly demonstrate that the proposed algorithm is a good candidate for a better specification of the regional 3-D ionosphere electron density field; reconstructed electron densities reflect changes following the energy sink into the ionosphere system during geomagnetic storms.
Moreover, the reconstructed densities are comparable to observations covering both bottom (ground-based ionosonde) and topside (space-borne satellites) ionosphere. Therefore, for a coherent 3-D ionosphere picture a joint use of both ground GNSS data and bottom side ground-based ionosondes should be emphasized. However, a more rigorous analysis should be performed to further restrain the variability in the reconstructed vertical structure. The introduction of RO data would significantly contribute to the structure stability. Indeed, then it would be of interest to extend our analysis to explore and quantify the impact of the plasmasphere on the stability/fidelity of the whole ionosphere 3-D structure.    The original tomography densities without the assimilation complexity are shown in green. Eminent is the robustness of the modified version (blue) to track the ionosonde measured (red) densities below the 270 km altitude range (purple line).    Comparison is performed along the trace of RO tangent points shown on maps in the upper right corner of each subplot (purple: GRACE, green: COSMIC). DOY 076 is not presented because the data points were few. Also, we limit the analysis to regions where RO density profiles have good accuracy (away from the low latitudes and altitudes above 200 km).