An ionospheric index suitable for estimating the degree of ionospheric perturbations

This is anOp Abstract – Space weather can strongly affect trans-ionospheric radio signals depending on the used frequency. In order to assess the strength of a space weather event from its origin at the sun towards its impact on the ionosphere a number of physical quantities need to be derived from scientific measurements. These are for example the Wolf number sunspot index, the solar flux density F10.7, measurements of the interplanetary magnetic field, the proton density, the solar wind speed, the dynamical pressure, the geomagnetic indices Auroral Electrojet, Kp, Ap and Dst as well as the Total Electron Content (TEC), the Rate of TEC, the scintillation indices S4 and s(’) and the Along-Arc TEC Rate index index. All these quantities provide in combination with an additional classification an orientation in a physical complex environment. Hence, they are used for brief communication of a simplified but appropriate space situation awareness. However, space weather driven ionospheric phenomena can affect many customers in the communication and navigation domain, which are still served inadequately by the existing indices. We present a new robust index, that is able to properly characterize temporal and spatial ionospheric variations of small to medium scales. The proposed ionospheric disturbance index can overcome several drawbacks of other ionospheric measures and might be suitable as potential driver for an ionospheric space weather scale.


Introduction
Space weather scales as e.g.introduced by NOAA (National Oceanic and Atmospheric Administration, see http://www.swpc.noaa.gov/noaa-scales-explanation)describe the strength of selected observables of the space environment by numbers which are related to characteristic effects on people and technical systems at each level.Such a scale related to ionospheric perturbations has been requested by user groups in the domain of communication and navigation for a long time, but could not yet been realized satisfactorily due to the complexity of the ionospheric reaction with respect to space weather events.Ionospheric perturbations caused by space weather events are able to effectively disturb or even interrupt radio systems like the Global Navigation Satellite Systems (GNSS).Therefore a comprehensive ionospheric scale needs to cover the spatial and temporal ionospheric response with respect to such disturbances.The strong seasonal and regional dynamics due to solar inclination, solar cycle, day-night cycle, ionospheric anomalies like cusp, crest and trough in combination with the coupling to incoming space weather events makes the assessment of the ionospheric state in a certain region highly complicated.Furthermore, the huge diversity of different structures within the ionosphere like plasma bubbles, patches, gradients Pradipta & Doherty (2016), traveling ionospheric disturbances (TIDs) (Borries et al., 2009) influence high frequency and trans-ionospheric radio wave propagation.So the associated impact causes amplitude and phase scintillation of GNSS signals Basu & Basu (1981), Béniguel et al. (2009), Hlubek et al. (2014), Kriegel et al. (2017) and other ionospheric effects on space-based radar Xu et al. (2004).
A physical measure to drive an ionospheric scale must be based on a reasonable global measurement available in real time with a high temporal resolution and should be free from any model assumptions.Moreover, the driver must be userfriendly and should allow a significant assessment of the disturbance level faced at system and service level of the different main user groups in the area of satellite communication and navigation, positioning as well as high frequency propagation.
Following the general approach of the Disturbance Ionosphere Index by Jakowski et al. (2012) here, we present a specific version focusing on spatial gradients.We demonstrate that the Disturbance Ionosphere Index Spatial Gradient (DIXSG) is a comprehensive index for risk estimation and performance degradation for users in the navigation and communication domain and has the potential to be a valuable driver for an ionospheric warning scale.In the following we present the mathematical formulation and the application of the DIXSG to the intense geomagnetic/ionospheric storm event on March 17, 2015 ("St. Patrick's Day Storm").Furthermore, we compare the results of DIXSG with the ionospheric indices Rate of TEC (ROTI) and Along-Arc TEC Rate index (AATR).Finally we discuss how the DIXSG index can be of benefit for the definition of an ionospheric warning scale.

DIXSG method
The ionospheric disturbance index proposed in this study (DIXSG) is a modified and further developed version of the one presented in Jakowski et al. (2012) which itself is based on studies published in Jakowski et al. (2006).Especially in this version with its focus on spatial gradients an implicit elevation weighting, distance dependency, scalability of the sensitivity and the ability to simplify the index to a chose-able number of levels has been included whereas an explicit temporal term has been dropped.Like the former one it is based on a double difference method using cycle slip removed, dual frequency GNSS carrier phase measurements as input.The applied cycle slip methods detect outliers by monitoring the Total Electron Content (TEC) phase rate, the Melbourne-Wübbena combination, the signal to noise ratio and the residuals from fitting a third degree polynomial using 1 Hz GNSS data.After an offset is detected by at least one of the these methods the TEC phase measurements are properly readjusted.From the difference of two corresponding carrier phase measurements, e.g. on L1 and L2, from a given GNSS satellite to a ground receiver station, the non-calibrated slant TEC (STEC) can be derived, cf.e.g.Jakowski et al. (2012).To get the TEC-rate DSTEC the difference of two consecutive slant TEC measurement delayed by Dt (e.g. 30 s) has to be calculated.Since not only the time has changed during the two consecutive measurements by Dt but also the location, the corresponding distance Ds may be estimated.This can be done for example by determining the distance of the two ionospheric piercing points (IPP) on an assumed single layer ionosphere in the height of the center of gravity of a mean electron density profile, e.g.400 km.Now the absolute value of the quotient of DSTEC and (Dt • Ds) for a given GNSS satellite (k) to ground receiver (i respectively j) link can be computed as: This first term of the DIXSG has the advantage of being implicit elevation weighted by the Ds factor which is big at low elevation angles and increasingly smaller at higher elevations.
The final DIXSG follows as: where cROT (level) represents a selectable recognition level to adjust the sensitivity, d the distance of the corresponding IPP's of the links from satellite k to receiver i, respectively j and D the maximum allowed distance.Although the possible distances are bounded by the receiver network geometry, it is possible to choose a sub-set in order to filter the ionospheric disturbances by its scales.For the results presented d was limited to values between d min = 10 km and d max = 1000 km and consequently D = 1000 km.The limits d min and d max can freely, i.e.only bounded by the GNSS ground receiver network geometry, be chosen to extract gradients of a specific scale length or in order to optimize the efficiency for a given GNSS receiver network.
The location of a DIXSG value is chosen to be at the center point of the IPP pair, i.e., by half the distance at the big circle between both points.
If one now want to simplify the DIXSG to a number of distinct levels the corresponding number and magnitude of levels (cROT (level) (L) with L = 1, 2...n) are successively applied and each time simplified in the way that values higher than one are defined as one and zero else-wise.Finally the generated values are summed correspondingly.
Since the number and magnitude of the levels are freely selectable and only orientated on the user requirements, five level were chosen with the magnitudes of 50, 100, 150, 200 and 250 to give an example.The index is now calculated as follows: For the mapping of the DIXSG the maximum value inside a respective 1°Â 1°area has been chosen to represent the corresponding part of the ionosphere over the earth's surface at a given time.
The arithmetic mean of all of these mapped values may serve as a global DIXSG or DIXSGp: where N is the total number of all valid 1°Â 1°areas.
Considering equations ( 1)-( 4), this specific index approach differs from the general disturbance index approach presented in Jakowski et al. (2012) by focusing on spatial perturbations, by taking into account the covered distances Ds of individual piercing points, by including the distance of the matched IPPs and finally by introducing a concrete perturbation scale.The fixed scaling enables us to simultaneously compare different regions and also to combine regional indices to a global one.

Application of the DIXSG to ionospheric storms
In order to verify the above described technique, the DIXSG was applied to the GNSS data sets of 13 ionospheric storm events (three from solar cycle 23 using 1/30 Hz GNSS data and ten from solar cycle 24 using 1 Hz GNSS data).In all cases the storm induced ionospheric disturbances were successfully detected.Here we will show as an example the intense geomagnetic storm on March 17, 2015 (St.Patricks day storm) on solar cycle 24 with its strong ionospheric perturbations (Cherniak et al., 2015;Borries et al., 2016).
The storm, later categorized as "severe" on the NOAA geomagnetic storm scale, started after a coronal mass ejection hit the Earth's upper atmosphere at 04:45 UT on March 17.The Dst value reached a minimum of À223 nT at 23:00 UT on the same day, cf.  2, it is easy to detect an unusual expansion of the electron distribution especially to the South, which becomes even more evident when considering the difference between TEC and the median TEC taken over 27 days.To show in addition the temporal evolution of the ionospheric storm over Europe (i.e. at 0°l ongitude) in a compact way a time versus value plot over 5 days (March 16, 2015to March 20, 2015) was extracted from a time series of the corresponding "TEC minus TEC median" maps (cf.Fig. 3).It reveals a strong pattern correlated to the storm effect occurrence at the second half of March 17, 2015 and some weaker ones in the following days as a result of the recovery phase of the storm.
In a following step, ROTI values from IMPC (Ionospheric and Prediction Center, http://impc.dlr.de/) were analyzed with respect to the storm event.ROTI or the Rate of change of TEC Index is defined as the standard deviation of the TEC-rate ROT in units of TECU (1 TECU = 10 16 m À2 ) per minute: for a given GNSS satellite and ground receiver.
Confer Cherniak et al. (2015) and Jacobsen & Andalsvik (2016) for an extensive application of ROTI values gained from high latitude GNSS receiver stations on the ST.Patrick storm.Figure 4 shows three ROTI maps made from public available GNSS receiver stations data during the storm event.
Here the one hour maximum values are shown.Again the  corresponding values were converted to a time versus value plot taking this time mean values over Europe in an area of À40°to 50°longitude East during the period of the storm.A clear pattern appears during the main phase of the storm in latitudes between 55°and 75°North, cf. Figure 5. Comparing this with the time series of "TEC minus 27 day median TEC" Figure 3 one recognizes that the disturbance activity happens at the same time of the highest deviation from the "27 day median TEC" but at much higher latitudes.The difference of the actual TEC and the 27 day median TEC values reflect the deviation from the "normal" behavior of the ionosphere.For this reason indices which describe main features of ionospheric storms are valuable for GNSS users.
As next the AATR Sanz et al. ( 2014) is applied to the storm data.The AATR was developed to support SBAS algorithms design, risk analysis and performance qualification systematically and on a rational basis, a simplified criteria or parameter that allows characterizing ionospheric conditions under performance aspects.
The AATR index is defined by the computation of: where T is the observation time period, AATR T the along arc TEC rate, DSTEC the slant TEC rate between consecutive observations, Dt the time elapsed between consecutive observations (i.e. 30 s or 60 s according to Sanz et al., 2014) and M(e) the spherical thin shell mapping function defined as: with e representing the satellite elevation angle, R e the radius of the Earth and h I the approximate height of the center of gravity of the mean vertical profile of the electron density.Being, where N is the number of observations in 1 hour.Therefore, the AATR index is defined as the hourly RMS of the instantaneous values of AATR T computed from all measurements collected by a given receiver.The unit of the AATR is mm/s derived from TECU/s. Figure 6 shows the AATR values in mm/s of the GNSS receiver station in Hoefn/Iceland (hofn) (Lat.: þ64.2673°N,Lon.: À15.1979°E, Height: þ82.3 m) and the corresponding Dst values to allow an easy comparison.The AATR values reflect the variation of the geomagnetic activity very good.Taking all one hour AATR values from 25 selected public available stations in Europe it is possible to create a time versus value plot again, cf. Figure 7.The used stations are (from low to high latitudes) mas1 (Maspalomas/Gran Canaria/Spain), bshm (Haifa/Israel), rabt (Rabat/Morocco), nico (Nicosia/ Cyprus), pdel (Ponta Delgada/São Miguel Island-Acores/ Portugal), madr (Madrid/Spain), ebre (Roquetes/Spain), ajac (Ajaccio/Corsica/France), mars (Marseille/France), tlse (Toulouse/France), pado (Padova/Italy), bzrg (Bolzano/Italy), zim2 (Zimmerwald/Switzerland), brst (Brest/France), wtzr (Bad Kötzting/Germany), redu (Redu/Belgium), titz (Titz/ Germany), pots (Potsdam/Germany), warn (Rostock-Warnemünde/Germany), sass (Sassnitz/Germany), zwe2 (Zwenigorod/Russia), spt0 (Boras/Sweden), mar6 (Maartsbo/ Sweden), hofn (Hoefn/Iceland) and kiru (Kiruna/Sweden).
Since the AATR summarizes all measurements at an individual station location the resulting index is not Geolocated, i.e. provides only a rough picture of the ionospheric perturbation in the region surrounding the receiver site in a circle which depends on the elevation cut-off angle.Hence the corresponding map cannot provide a detailed image of the temporal and spatial storm features as confirmed in Figure 7.
Finally the derived five level DIXSG maps were analyzed.An example during the main phase of the storm is shown in Figure 8 globally and for Europe in Figure 9.The high level of the ionospheric disturbances in the Northern part of Europe is clearly reflected by the higher level of the index.This result could be achieved although the public available GNSS receiver station data have great lacks in the northern part of Europe.Like before the corresponding data from the time series of DIXSG maps are condensed in a time versus value plot showing the storm evolution within a five days period.Like in the case of ROTI the storm pattern becomes clearly visible but here even more pronounced in its characteristics, cf. Figure 10.
The large scale TEC increase during the positive storm phase (cf.Fig. 3), which correspond to the findings of Borries     caused by electron density patches or particle precipitation Borries et al. (2016).Thus, DIXSG is well suited to detect small and medium scale ionospheric irregularities.
Taking public available one second data (about 140 for the period of the discussed storm) for the calculation of a global DIXSG map in the described way (i.e.1°Â 1°area elements and etc.) it is possible to cover approximately 5% of the earth surface.From this data, a mean global DIXSG can be calculated like it is shown as red bars in Figure 1, while the blue line in the same plot indicates the Dst values to allow direct comparison.Although the Dst is a geomagnetic index and the DIXSG is purely an ionospheric index, they both show the same characteristic behavior with a high value of correlation.
To further test the ability of the DIXSG to correctly identify storm pattern in its temporal and spatial extent a number of storms were analyzed as mentioned before.

Conclusion
Since many years, the NOAAs Space Weather Scales serve as a guideline for risk assessment in case of space weather events.There exist three scales "Geomagnetic Storm", "Solar Radiation Storm" and "Radio Blackout", which inform users about strength, duration, impact on technical infrastructure, yearly probability of occurrence as well as on the underlying physical measure defining the scale.There have been small adaptations in wording and settings of the threshold in the past, but in general the overall system remained stable and clearly defined.However, an additional ionospheric scale has been requested by user groups in the domain of communication and navigation for a long time, but could not yet been realized sufficiently due to the complexity of the ionospheric reaction with respect to space weather events.The different scales and high variability of ionospheric disturbances require most probably different indices for a proper categorization of phenomena to the users operating a wide variety of  technological systems each with its own specific and requirements with respect to accuracy, precision, availability, integrity etc. and further will be met once the focus is moved from monitoring to prediction.
Nevertheless the authors believe that a TEC based index like the proposed ionospheric disturbance index can contribute to the definition of a global ionospheric disturbance scale.The introduced DIXSG has proven to be an easy to calculate, reliable and robust index with which it is possible to clearly identify ionospheric disturbances not only in time but also in terms of the affected area.It is scalable in its sensitivity to the strength of perturbations and it is possible to choose different level of recognition depending on the demands of a specific application.This makes it more suitable to a number of applications.In comparison with the AATR, which is available hourly at the location of the GNSS-receiver only, the DIXSG can cover large areas between the stations, depending on the cut-off elevation angle, with an update-rate of up to one minute.
The proposed DIXSG will be implemented as real time processor and the corresponding results will be offered to users for application-oriented utilization via IMPC.
Figure 1.Borries et al. (2016) identified four different phases of the storm with southward propagating ionospheric disturbances during the fourth phase, on the Northern hemisphere, originating at polar latitudes at the second half of March 17, 2015 in the European-African sector.To show the strength and expansion of the ionospheric disturbances due to the storm IGS (International GNSS Service) 2-hour TEC maps (Hernández-Pajares et al., 2009) were analyzed.In the IGS TEC map Figure

Fig. 1 .
Fig. 1.Global geomagnetic Dst value in the period of March 16 to 20, 2015 (blue line) and global (coverage approximately 5%) mean of the five level DIXSG values (red bars).The calculation is based on public available 1 Hz GNSS data.

Fig. 3 .
Fig. 3. Time series of 'TEC minus 27 day median TEC' at zero degree East longitude in TECU extracted from IGS 2 hour TEC maps taken from ftp://cdaweb.gsfc.nasa.gov/pub/data/gps/tec2hr_igs/.Shown is the period of March 16 to 20, 2015.Note the strong storm pattern at the second half of March 17, 2015.

Fig. 4 .
Fig. 4. Rate of change of TEC index (ROTI) in TECU per minute.Taken are the corresponding maximal values within an hour in an 2°Â 2°area element.The calculation is based on public available 1 Hz GNSS data.

Fig. 5 .
Fig. 5. Time series of ROTI in the period of March 16 to 20, 2015 over Europe, i.e. corresponding mean values were taken between longitude À40 to 50 degree East.The calculation is based on public available 1 Hz GNSS data.

Fig. 6 .
Fig. 6.Along Arc TEC Rate (AATR) at Hoefn/Iceland in mm/s (red bars) and the corresponding Dst values (blue line).For the calculation of the AATR a Dt of 30 s were chosen using 1 Hz GNSS data.

Fig. 8 .
Fig. 8. Global five level DIXSG at March 17, 2015 17:00 to 18:00 UT. are the corresponding maximal values within an hour in an 1°Â 1°a rea element.The calculation is based on public available 1 Hz GNSS data.

Fig. 7 .
Fig. 7. Time series of AATR values at 25 stations in the period of March 16 to 20, 2015 in Europe.The calculation is based on public available 1 Hz GNSS data.

Fig. 9 .
Fig. 9. Five level DIXSG over Europe.Taken are the corresponding maximal values within an hour in an 1°Â 1°area element.The calculation is based on public available 1 Hz GNSS data.
As a summary of their results the corresponding global DIXSG values of six of them are plotted exemplary in combination with the Dst values in Figures 1, 11-15.Figures 11 and 12 show two strong ionospheric storms following the "St.Patricks day storm" in the same year.Like in the former case the DIXSG changes reflect the Dst characteristics.The last three examples are very strong events from solar cycle 23, cf.Figures13-15.All cases again give evidence for the ability of the DIXSG to correctly trace the storm pattern, although only 30 s-data were used for this storm analysis while all previously discussed storm cases base on 1 s-GNSS-data.

Fig. 10 .
Fig. 10.Time series of the five level DIXSG (cf.Fig. 9) in the period of March 16 to 20, 2015 over Europe, i.e. corresponding maximal values were taken between longitude À40 to 50 degree East.The calculation is based on public available 1 Hz GNSS data.

Fig. 11 .
Fig. 11.Global geomagnetic Dst value in the period of June 22 to 27, 2015 (blue line) and global (coverage approximately 5%) mean of the five level DIXSG values (red bars).The calculation is based on public available 1 Hz GNSS data.

Fig. 13 .
Fig. 13.Global geomagnetic Dst value of the "Halloween Storm" in the period of Oct. 29 to Nov. 01, 2003 (blue line) and global mean of the five level DIXSG values (red bars).The calculation is based on public available 1/30 Hz GNSS data.

Fig. 12 .
Fig. 12. Global geomagnetic Dst value in the period of Dec. 19 to 24, 2015 (blue line) and global (coverage approximately 5%) mean of the five level DIXSG values (red bars).The calculation is based on public available 1 Hz GNSS data.

Fig. 14 .
Fig. 14.Global geomagnetic Dst value in the period of Nov. 19 to 22, 2003 (blue line) and global mean of the five level DIXSG values (red bars).The calculation is based on public available 1/30 Hz GNSS data.