Degradation of Kinematic PPP of GNSS Stations in Central Europe Caused by Medium-Scale Traveling Ionospheric Disturbances During the St. Patrick’s Day 2015 Geomagnetic Storm

In solar cycle 24, the strongest geomagnetic storm took place on March 17, 2015, when the geomagnetic activity index was as high as -223 nT. To verify the impact that the storm had on the Global Navigation Satellite System (GNSS)’s positioning accuracy and precision, we used 30-s observations from 15 reference stations located in Central Europe. For each of them, we applied kinematic precise point positioning (PPP) using gLAB software for the day of the storm and, for comparison, for a selected quiet day (13 March 2015). Based on the conducted analyses, we found out that the position root mean square (RMS) values on the day of the geomagnetic storm were significantly high and amounted to several dozen centimeters. The average RMS for the altitude coordinates was 0.58 m between 12:00 and 24:00 (GPS time), and 0.37 and 0.26 m for directions North and East, respectively. The compromised accuracy level was caused by a sudden decrease in the number of satellites used for calculations. This was due to a high number of cycle slips (CSs) detected during this period. The occurrence of these effects was strictly correlated with the appearance of traveling ionospheric disturbances (TIDs). This was proven by analyzing changes in the total electron content (TEC) estimated for each station–satellite pair.


Introduction
The ionosphere is a layer of the Earth's atmosphere, wherein density distribution of free electrons is heterogeneous and is determined by the intensity of solar radiation, which stimulates ionization of atoms and gas particles. The density of free electrons and their variations significantly affects electromagnetic signals passing through the atmosphere. Along with the changes in the signal propagation velocity, irregularities of electron density may cause its diffraction and refraction. An assessment of the current state of the ionosphere can be made by analyzing the total electron content (TEC) parameter. The TEC, which represents the ionospheric refraction, depends on a given geographical location, the hour of a day and solar activity levels. Since the ionosphere is a dispersive region, refraction depends on the frequency of a signal. In case of the Global Navigation Satellite System (GNSS), 1 TECU (TEC units) is responsible for approximately a 0.16 m delay in L1 band. For single frequency receivers, ionospheric delay can be reduced only by modelling its value. Users can apply Klobuchar model coefficients [1] broadcasted with a GPS navigation message, which allows them to reduce the root mean square (RMS) ionospheric range error by about 50%. With the launch of the Galileo system, the NeQuick model was introduced [2], characterized by better quality; it can mitigate the ionospheric delay up to 68% [3]. Nonetheless, these values are insufficient for high precision applications. Thus, differential measurements can be used. With a vector length not exceeding 10 km, ionospheric refraction may be neglected [4]. However, long-range measurements are affected by a decrease in accuracy [5]. The solution is to use an ionospheric-free linear combination, which can eliminate 99% of first-order ionospheric delay but requires double frequency receivers. The second-and third-order delay can be efficiently modeled [6]. Ionospheric-free linear combination is applied particularly in the precise point positioning (PPP) approach [7,8], where only one receiver is used. However, this combination increases observation noise and makes resolving ambiguities difficult. It is worth noting that the ionosphere can significantly affect signal receiver acquisition. This may cause scintillation by rapid fluctuations in the signal amplitude and phase [9,10]. The accidental or rare phenomena, such as traveling ionospheric disturbances (TIDs) [11], may also be a problem. TIDs are sudden changes in plasma density in the ionosphere causing interference in a radio signal by changing the refraction index. They can cause sudden and notable changes in TEC, which can result in the hindered acquisition of GNSS signals and decrease positioning accuracy and precision [12]. TIDs can be generated not only by solar activity, but also by other sources, e.g., earthquakes [13] or intense phenomena in the troposphere [14]. Luo et al. (2018) investigated the accuracy of PPP positioning during a selected geomagnetic storm in solar cycle 24. Their findings clearly show that the RMS of a 3D position can reach over 1.5 m during superstorms [15]. This decrease of the accuracy was mainly caused by the occurrence of cycle-slip (CS) effects, which is a confirmation of previous works [16]. Moreno et al. (2010) found PPP large variations during the periods when the large rate of TEC (ROT) was observed at equatorial latitudes between April 18 and 22, 2004 [17]. Similar results were obtained by Jacobsen and Andalsvik who noticed the rapid increasing of vertical position errors together with increasing ROT index (ROTI) for both PPP and real-time kinematic (RTK) techniques [18]. Lu et al. (2020), based on data from Hong Kong, showed that ionospheric scintillation events vastly decreased the PPP accuracy due to the negative impact of cycle slips on convergence time [19]. Additionally, small-scale ionospheric irregularities of auroral origin can cause a PPP positioning error of about 0.5 m [20].
Because GNSS signals are sensitive to an ionospheric state, they can be successfully used for the purpose of investigating it. The research concerns global and local TEC mapping [21,22], investigating TIDs [23,24], modeling vertical profiles [25], estimating the height of ionospheric disturbances [26], etc.
In this paper, we present the negative effect of medium-scale TIDs that occurred during a severe geomagnetic storm on St. Patrick's Day in 2015 on kinematic PPP solutions. Our analyses were performed based on observations derived from a reference network in Poland. We focused on relative changes in TEC values and correlated them with the quality of the received signal. In Section 2, we briefly introduce the storm that took place on 17 March 2015. The data and the methodology are described in detail in Section 3. The results of the kinematic PPP positioning are show in Section 4. This is followed by a discussion in Section 5 where potential causes and effects of the obtained lower accuracy are analyzed. The article concludes in Section 6 with a summary of the obtained results.

The St. Patrick's Day 2015 Geomagnetic Storm
The source of the geomagnetic storm of 17-18 March 2015 (called the St. Patrick's Day event) was the ejection of coronal mass from the Sun on March 15. The initial velocity of the coronal masses' propagation in space was about 668 km/s, which made it possible for the magnetic cloud, formed during this phenomenon, to reach the area around the Earth on 17 March at 03:59 UT [27]. As a result of the interaction of echoes of coronal ejection with the magnetic field of the Earth, a geomagnetic storm was formed. The measurements conducted on that day showed a very low value of the disturbance storm time (Dst) parameter, which was -223 nT during the storm (Figure 1a). This allows it to be classified as a G4 (severe) geomagnetic storm [28]. This is also confirmed by a very strong total interplanetary magnetic field, which was over 35 nT about 13:30 UT on 15 March (Figure 1b). In Figure 1, we also present the planetary 3-h range Kp index to show the size of geomagnetic field Downloaded from mostwiedzy.pl variations caused by the system irregulars. All these parameters indicate that it was the most intense geomagnetic storm of the 24 th solar cycle [29].

Methodology and Data
In this study, we used 30-s GPS observations derived from 15 reference stations belonging to a national ASG-EUPOS network ( Figure 2a). Calculations were performed for two days, i.e., the one on which the geomagnetic storm disturbance occurred (17 March 2015) and a quiet day as a reference (March 13, 2015). The gLAB 5.4.4 software [30] was used for obtaining the kinematic PPP solution for each station with an interval equal to the data acquisition. The advantage of the gLAB software is epoch-by-epoch processing using the Kalman filter, which is very similar to the approach used in real-time processing. To ensure high accuracy and precision of the PPP solution, it is necessary that precise ephemeris and clocks of the satellites be used. In this study, products from the International GNSS Service were applied [31]. This allowed us to achieve a position with an accuracy of about 3 cm [32]. First-order ionospheric delays were reduced using dual-frequency observations, while tropospheric effects were eliminated by estimating correction to the a priori tropospheric delay. To reduce the multipath effect and remove signals passing low above the horizon, we adopted a 10° elevation mask. In addition, we considered satellites' and receivers' antenna phase center offsets and their variations, relativistic effects, as well as wind-up and solid-tide corrections. Cycle slips (CSs) were detected using two methods based on the Melbourne-Wübbena [33] and the geometry-free [34] linear combinations.
In order to determine the impact of TIDs on the position, we estimated vertical total electron content (VTEC) variations (ΔVTEC) based on observations from 650 reference stations (Figure 2b; coordinates presented in Table A1 in the Appendix). For each station-satellite pair, ΔVTEC time series were determined according to the methodology described by Nykiel et al. [24,26]. The use of a dense network of receivers allowed us to obtain maps for individual satellites. These maps are characterized by high 0.1° spatial resolution. The time resolution was consistent with the observation data (30 s), so it was possible to compare positioning results with ΔVTEC data explicitly. Examples of ΔVTEC maps showing the movement of TIDs determined from signals of GPS satellite number 24 are shown in Figure 2. As presented by Nykiel et al. [24], during the active phase of the storm (between 16:00 and 19:00 UTC) the RMS value of ΔVTEC maps amounted to more than 1 TECU, which is significantly higher than the average value estimated for quiet days (0.1 TECU). Moreover, at this time, the main direction of the motion of the inhomogeneities was from east to west with a velocity of about 1 km/s. Figure 3 shows ΔVTEC results of GPS/GLONASS satellites for selected stations analyzed in this study. The TIDs between 16 and 18 UT, causing VTEC disturbances of up to ±5 TECU, are clearly visible. They translate into signal delays of about 0.8 m on L1 frequency. Such sudden and considerable ΔVTEC jumps can affect the quality of satellite tracking and, consequently, positioning results.

Results
In this section, we present results of the coordinates' estimation conducted for the quiet (13 March) and the stormy (17 March) days. For the example of the CHNO (Figure 4, left) and KOSZ (Figure 4, right) stations, it can be observed how the position was changing depending on the ionospheric conditions. In both these cases, time series of topocentric coordinates clearly varied during the day with ionospheric disturbances (the red line) compared to the quiet day (the green line). The position variations can be observed in two main episodes. The first episode started around 16:00 and lasted until 19:00, while the second one started around 22:00 and lasted until the end of the day. The occurrence of position degradation is consistent with TIDs' appearance presented in Figure 3. In Figure 4, a slight delay between the disturbances in the KOSZ and CHNO stations is visible, which is due to the TIDs' movement.
The position degradation is evidenced in RMS values. During the quiet day, the RMS of horizontal coordinates was in the range of 0.01 and 0.02 m for the CHNO and KOSZ stations, respectively. At the same time, the RMS of the vertical component was only slightly higher and amounted to 0.03 and 0.04 m. Significantly higher values were obtained during the stormy day. For the CHNO station, they were 0.50, 0.32, and 0.63 m for North, East, and Up components, respectively, whereas for KOSZ they were 0.53, 0.53, and 0.99 m, respectively. The detailed information regarding the impact of ionospheric disturbances on the position quality can be found in Appendix A, which contains differences in North ( Figure A1), East ( Figure  A2), and Up ( Figure A3) components, between the stormy day and the reference day. To underline the impact of the geomagnetic storm itself, the time span of the analyzed data was reduced to the 12:00 -24:00 GPS period. For the North component, it was observed that ionospheric anomalies did not equally affect all stations. It is clearly visible that each station is represented by slightly different variations. Although, for most of them, the first episode of coordinate variations started around 16:00 and lasted three hours, for three stations (CCHN, ILAW, and WRKI) it started around 17:30 and lasted only half an hour. Similar discrepancies were found for the second episode. For most of the stations, it occurred between 22:00 and 24:00, while three stations (CCHN, GIZY, ILAW) were not clearly affected by the second ionospheric event. In the case of the East component, the distribution of the occurrence of both episodes was similar to the case of the North component differences. The exceptions here were the STRG and CHOJ stations, for which the variations obtained in the second episode were almost negligible. Variations in the vertical coordinates during both ionospheric episodes were similar to those found for horizontal components. The ionospheric storm had the least impact on the position quality at the ILAW and CCHN stations, while the WLAD station seems to have been affected the most.
The RMS positions obtained for the topocentric coordinates are summarized in Table 1. They were calculated for both the quiet day and the stormy day between 12:00 and 24:00 GPS time. While the quiet day was characterized by the average RMS equal to 0.02, 0.02, and 0.03 m, for the North, the East, and the Up components, respectively, the maximum RMS value did not exceed 0.05 m in any case. In contrast to this, much greater discrepancies between the coordinates took place in the occurrence of ionospheric disturbances. This can be observed, among others, in the average RMS that was equal to 0.37, 0.26, and 0.58 m for the North, the East, and the Up components, respectively.

Discussion
Above, we presented results for the kinematic PPP determination. Clearly, ionospheric disturbances had a considerable negative impact on solution accuracy. In this section, we discuss why this decreased accuracy is so significant. In Figure 5, the number of satellites used for position estimation is shown for each of the analyzed stations. It can be seen that, during the geomagnetic storm, the number of the used satellites falls significantly to 4, which directly translates to the value of the geometrical dilution of precision (GDOP) factor ( Figure 5, the magenta line). Its value is inversely proportional to the number of satellites. Here, we can observe high GDOP values, which at some moments were even greater than 14. This means that the estimated position was characterized by a low confidence level. High GDOP values are strictly correlated with compromised accuracy, which is clearly noticeable in the presented results ( Figure A1-3). Above, we show epochs when the positions were estimated. However, it is worth noting that, in some cases, the number of observed satellites was so low that the position equations could not be resolved. Such data breaks occurred most often for the WLAD (23 breaks) and KOSZ (17 breaks) stations. The average number of missing epochs was 10.
We determined the cause of these drops in the number of satellites. In Figure 6, we present the values of the estimated phase ambiguities. Clearly, no significant deterioration in ambiguity resolution was observed during calm ionospheric conditions (Figure 6, left). In contrast to this, the occurrence of ionospheric disturbances resulted in a clear degradation of ambiguity resolution at each of the analyzed stations ( Figure 6, right). In most cases, this took place during two episodes, namely, at 16:00 -19:00 and 22:00 -24:00 GPS time, exactly at the same time as in the case of the previously analyzed disorders. The effect characterized by a leap in the phase of the signal received from the satellite is called a cycle slip (CS). As a result of a CS, a given satellite cannot be used for solving the position. When this situation occurs for many satellites, determining the high accuracy position proves impossible. Moreover, here we can observe that ambiguities estimated during the main phase of the storm are very dispersed, which makes it difficult to apply CS repairing methods.  In Figure 7, the total number of CSs that occurred is shown for each of the analyzed stations. As can be seen, on the stormy day there were significantly more CSs than on the quiet day. The highest number of CSs occurred for the measurements at the WLAD station and amounted to 931 between 12:00 and 24:00 GPS time. This result is over twice as high as that observed for this station on the calm day. The second largest difference was recorded for the ILAW station, where the difference between the number of CSs for these two days was 416, meaning that the number of CSs observed on the stormy day was five times higher. Similar results were obtained for the CCHN station, where the number of the CSs observed during the stormy day was six times higher than that recorded on the calm day. The average number of CSs on the quiet day was 209, while on the stormy day it was estimated at 559, which is about 2.5 times more. To clearly show the impact of TIDs on discarded satellites from positioning, we estimated the average ΔVTEC for each station. Figure 8 presents the values of ΔVTEC for selected (included in positioning) and discarded satellites. In general, for discarded observations the ΔVTEC was about 0.95±0.49 TECU, which is twice as high as the value for the selected measurements (0.49±0.20 TECU). In most cases, the ΔVTEC for excluded observations is significantly higher. The highest value was observed at the CHOJ station and amounted to 2.29±1.48 TECU, which is over twice as high as the Downloaded from mostwiedzy.pl result for the selected observations (0.95±1.17 TECU). The three stations (CHNO, PPIL, and WRKI) are the exceptions. In these cases, the ΔVTEC of discarded observations is similar (WRKI and CHNO) or lower (PPIL) compared to the considered ones. This is due to the fact that the occurrence of CSs is determined not solely by the average value of ΔVTEC but also by the rapidity of its changes.

Conclusions
The analyzed geomagnetic storm was the result of coronal masses ejected from the Sun that reached the Earth's atmosphere on 17 March 2015. The ionospheric disturbances it caused clearly affected the quality of position determination based on the GNSS technique. Based on the measurements derived from the 15 reference stations of the ASG-EUPOS network, a position with an interval of 30 s was estimated using the kinematic PPP method.
In this paper, we showed the impact of the St. Patrick's Day storm on position. We observed a clear decrease in position accuracy in the main phase of the storm, during which high-intensity medium-scale TIDs occurred. Decreased accuracy was observed at all the analyzed stations and was caused by a fall in the number of satellites considered in positioning. We showed that the sudden drop in their number was caused by numerous occurrences of CSs, which made the satellites discarded in specific epochs. By analyzing the values of and changes in ΔVTEC over time for each of the receiver-satellite pairs, we proved that the occurrence of CSs was caused by TIDs. Moreover, in some cases, the decrease in the number of satellites was so severe that determining the position was impossible. On the day when the geomagnetic storm occurred, the RMS of the position was 0.37 ± 0.16, 0.26 ± 0.12, 0.58 ± 0.33 m on average for the North, East, Up topocentric components, respectively. For comparison, for the quiet day these values ranged from 0.02 to 0.03 m. We showed that, during the TIDs' transition, the CS occurrences were several times greater than during the quiet day. Furthermore, given the determined ambiguities, we claim that, during the storm, correct tracking of the phase of the signal proved troublesome for the receivers.

Conflicts of Interest:
The authors declare no conflict of interest. Downloaded from mostwiedzy.pl