Occurrence of GPS Loss of Lock Based on a Swarm Half-Solar Cycle Dataset and Its Relation to the Background Ionosphere

: This paper discusses the occurrence of Global Positioning System (GPS) loss of lock events obtained by considering total electron content (TEC) measurements carried out by the three satellites of the European Space Agency Swarm constellation from December 2013 to December 2020, which represents the longest dataset ever used to perform such an analysis. After describing the approach used to classify a GPS loss of lock, the corresponding occurrence is analyzed as a function of latitude, local time, season, and solar activity to identify well-deﬁned patterns. Moreover, the strict relation of the occurrence of the GPS loss of lock events with deﬁned values of both the rate of change of electron density index (RODI) and the rate of change of TEC index (ROTI) is highlighted. The scope of this study is, on one hand, to characterize the background conditions of the ionosphere for such events and, on the other hand, to pave the way for their possible future modeling. The results shown, especially the fact that GPS loss of lock events tend to happen for well-deﬁned values of both RODI and ROTI, are of utmost importance in the light of Space Weather effects mitigation.


Introduction
Electron density irregularities are at the base of sudden and rapid fluctuations characterizing both amplitude and phase of electromagnetic waves propagating through the ionospheric plasma, e.g., [1][2][3][4][5]. These fluctuations significantly affect the Global Navigation Satellite System (GNSS) and, in the worst case, can even cause a loss of lock (LoL) event, a condition for which a GNSS receiver can no longer track the signal sent by the satellite, with a consequent degradation of the positioning accuracy. Due to the relevance that these phenomena have in the framework of Space Weather, several works have recently focused on this matter. Buchert et al. [6] investigated the early LoL events affecting the European Space Agency (ESA) Swarm mission in January and February 2014 at low latitudes, and ascribed them to electron density gradients associated with equatorial plasma bubbles (EPBs) immersed in a very high electron density background. Using Global Positioning System (GPS) radio occultation observations of the Constellation Observing System for Meteorology, Ionosphere, and Climate, Yue et al. [7] analyzed the signal cycle slip occurrence from 2007 to 2011. They correlated it with the appearance of sporadic E layers, and more specifically with F-region irregularities and the equatorial ionospheric anomaly (EIA) at low and equatorial latitudes, and with the polar cap electron density gradients in the polar regions. Liu et al. [8] analyzed the LoL occurrence characterizing the GNSS data collected in Weipa (Australia; 12.45 • S, 130.95 • E) from 2011 to 2015, and highlighted how these events are closely related with the solar activity. Using GPS data recorded at Skibotn (Norway; 69.43 • N, 20.38 • E), Jin and Oksavik [9] found the occurrence of many LoL events connected with the St. Patrick storm that occurred in March 2015. On the contrary, Piersanti et al. [10], analyzing the geomagnetic storm occurred in August 2018, found no LoL events in the GPS data recorded by the ESA Swarm constellation, and attributed this to the fact that the storm was too weak to seriously degrade GPS signals. Damaceno et al. [11], using data from four GPS receivers located in a latitudinal belt ranging from about 25 • S to 2 • S in the eastern part of Brazil, performed a statistical analysis of LoL occurrence throughout the 24th solar cycle. Their results showed that LoL events occur mainly in the post-sunset hours, especially within the Southern crest of the EIA, i.e., where EPBs are most likely to occur. Zhang et al. [12] made a statistical study of the GPS LoL occurrence during 14 intense geomagnetic storms, defining at the same time two indices to identify, respectively, the number of LoL events affecting a single site and the global average rate of LoL. Xiong et al. [13,14] were the first that, using L1 and L2 GPS measurements acquired by the Swarm constellation from December 2013 to December 2016, showed a global view of the LoL occurrence.
This work is based on an extended dataset of LoL events recorded by the Swarm constellation from December 2013 to December 2020, the longest ever used, and discusses the corresponding occurrence as a function of latitude, local time, season, and solar activity. In addition, our analyses aim at finding a relation between the LoL occurrence and defined values of the following two ionospheric indices: the rate of change of electron density index (RODI) and the rate of change of total electron content (TEC) index (ROTI). This kind of investigation is done both to characterize the background conditions of the ionosphere for such events, and to start understanding whether these events can be somehow forecast, which would be very important for Space Weather effects mitigation purposes. Section 2 describes data and methods used for the analyses. Section 3 is devoted to the results and corresponding discussion, while the summary and conclusions are presented in Section 4.

Swarm Data
Swarm is an ESA constellation constituted by three low-Earth-orbit satellites launched at the end of 2013, and still operating, in a circular near-polar orbit [15]. Two of them (Swarm A and C) have the same orbit configuration: An inclination of 87.35 • , an initial altitude of about 460 km, and an initial east-west separation of about 1 • -1.5 • in longitude. The third (Swarm B) has a different orbital configuration: An inclination of 87.75 • and an initial altitude of about 520 km. For this study, we consider data recorded from all the three Swarm satellites for the seven-year period between 2 December 2013 and 31 December 2020.
We used Langmuir Probes (LPs) [16] and Precise Orbit Determination (POD) antennas [17] data. LPs provide measurements of in situ electron density (Ne) with a 2-Hz rate [18,19]. POD antennas are GPS receivers through which TEC values are estimated [20,21]. Level 2 (L2) TEC data contain time series of vertical and slant (both absolute and relative) TEC for each GPS satellite in view (at most eight due to the instrumentation design). TEC time series are at 0.1 Hz until 15 July 2014, and since then at 1 Hz. Swarm data are freely downloadable at ftp://swarm-diss.eo.esa.int (accessed on 31 May 2021). LP data are provided with different flags [18]. The ones we are interested in are two: Flags_LP related to the LP instruments; Flags_Ne related to electron density measurements. Only the most reliable LP data are considered for this study, namely, those with Flags_LP = 1 and Flags_Ne ≤ 29. Differently, no flag identifies TEC data, so all TEC data are considered.

Ionospheric Indices: RODI and ROTI
By means of Ne and TEC measured time series, it is possible to calculate indices, like RODI and ROTI, characterizing their fluctuations at different spatial and temporal scales. The Swarm constellation, due to its orbit configuration and to the high sampling rate of its measurements, is particularly suited for the calculation of these indices. Several studies, exploiting Swarm data, have used these indices to get very good results about the topside ionospheric dynamics characterizing the high latitudes [22,23], as well as the low latitudes [13,24,25].

RODI
RODI values [26] can be obtained according to the procedure recently proposed by Pignalberi [27] through his TITIPy (Topside Ionosphere Turbulence Indices with Python) tool.
Specifically, the rate of change of electron density (ROD) is the time derivative of the electron density calculated along the satellite's orbit using a 1st order finite difference scheme. Ne values measured by LPs onboard Swarm satellites in principle are continuous time series; however, missing measurements are possible. By indexing with k the contiguous measured values, the kth ROD value is defined as where Ne k is the electron density measured at time t k . ROD values are calculated only for time consecutive measurements, that is, when the condition (t k+1 − t k ) = δt = 0.5 s is met. Then, the kth RODI value is calculated as the standard deviation of ROD k in a window of N = 2j + 1 elements, i.e., where ROD k+i are the ROD values falling inside the window centered at the index k, whose width is N = (2j + 1), and ROD k is The choice of the width of the window used is crucial for the detection of fluctuations of different spatial and temporal scales. A good compromise is to select a 10-s window, as was done for instance by Jin et al. [22] and De Michelis et al. [23], which is precisely the width used in this study. This means that j = 10 and (at most) 21 values are used for RODI calculation.
RODI value is calculated only when at least half of the values are present in the window; otherwise, the corresponding value is considered not statistically significant and discarded. Then, in our case where ∆t = 10 s and δt = 0.5 s, at least 11 ROD values are needed to start calculating the index.

ROTI
ROTI value [28] is calculated as RODI but considering TEC values in place of Ne values. However, while for each instant of time a single value of Ne is available, several TEC values can be available, because Swarm satellites can track up to eight GPS satellites simultaneously.
Specifically, the rate of change of TEC (ROT) is the time derivative of TEC calculated along the path between a Swarm satellite and a specific GPS satellite. For the sake of simplicity, in the following, we simply refer to TEC values, having in mind that both slant TEC (sTEC) and vertical TEC (vTEC) values can be used in the ROT calculation. In this study, we use absolute sTEC values, then where s is the index running on the GPS satellites in view, identified by their Pseudo Random Noise (PRN) number. The kth ROTI value, for the GPS satellite s, is calculated as with The window used for ROTI calculation is as wide as that used for RODI calculation.

GPS Loss of Lock Identification Methodology
In our work, the GPS LoL is identified directly from the L2 sTEC time series calculated from ambiguity-corrected carrier phase observations, performed on L 1 and L 2 signals at frequencies that are respectively equal to 1575.42 MHz and 1227.60 MHz, contained in the receiver-independent exchange files generated by the POD antennas [20]. So, when L 1 or L 2 (or both) are unavailable, the corresponding sTEC value cannot be calculated. This can happen for two reasons: (1) when the GPS satellite is outside the POD antenna field of view; or (2) when a LoL is ongoing. In this work, LoL events are identified by looking for interruptions in the sTEC time series for a specific GPS satellite (identified by the corresponding PRN) following a 4-step procedure: sTEC time series of each specific GPS satellite (then, of each specific PRN) are extracted from the L2 TEC data files (which contain data of all the 32 GPS satellites); 2.
From the sTEC time series of each specific PRN, we identify the segments of orbit for which the GPS satellite is in the field of view; 3.
For each of these segments of orbit, interruptions in the sTEC time series ranging from 1 to 1200 s are identified. The choice of this range guarantees that the identified interruptions are actually LoL and not due to the fact that the satellite has gone out of the field of view. This means that if, for instance, a LoL begins just before the satellite is going outside the POD field of view, this event will not be considered in our dataset; 4.
Steps 1-3 loop for each of the 32 GPS satellites (then, for each PRN).
An example of a GPS LoL is given in Figure 1 for an event recorded by Swarm A on 17 March 2015, during the development of the well-known St. Patrick storm. Specifically, the event represented in Figure 1 relates to PRN = 9, that was in the Swarm A field of view between 13:46:35 and 14:06:50 universal time (UT) (ascending orbit), and is identified by the interruption in the sTEC time series. Measured sTEC values, and calculated RODI and ROTI values, are represented in geographic coordinates, and as plots latitude vs. value. RODI values are obtained through (1)-(3) from Ne measurements, and ROTI values are those calculated by applying (4)-(6) to PRN = 9 sTEC measurements. By looking at sTEC panels, it is clear how sTEC increases from South to North until the time series is interrupted (between dashed magenta curves) and then reappears after a definite time interval, whose duration in this case was of 191 s. The duration of the LoL event is calculated considering the difference between the timestamps corresponding to the magenta curves in Figure 1. The same behavior is also shown by ROTI values that are derived from sTEC. A quite similar behavior is displayed also by RODI values, with the difference that RODI can be calculated also during a LoL event, if Ne data are available and not flagged. The most important thing to notice is that both ROTI and RODI show an abrupt increase before the GPS signal is lost, and then they are both by far lower when the GPS signal reappears (the scale of both ROTI and RODI in Figure 1 is logarithmic).
Since RODI is an in situ index, being derived by in situ Ne values, while ROTI has no simple spatial characterization, being associated with the path linking Swarm and GPS satellites, the association of RODI and ROTI values to a specific LoL event is tricky. In fact, a LoL event is due to an ionospheric irregularity whose exact spatial location cannot be identified and depends on the geometric configuration of Swarm and GPS satellites. In Pignalberi [27], sTEC and corresponding ROTI values are associated with the Ionospheric Pierce Point (IPP), that is, the point where the link path between GPS and Swarm satellites pierces the spherical thin shell located 400 km above the Swarm satellite orbit. Depending on the elevation and azimuth angles between Swarm and GPS satellites, the IPP location can differ to some degrees in both latitude and longitude from the Swarm one, and this spatial difference results in a time difference between RODI and ROTI values. For these reasons, in this work, RODI and ROTI values associated with a specific LoL event will be the corresponding maximum values in a 10-s wide window centered at the LoL start timestamp (that in Figure 1 is the timestamp corresponding to the lower magenta curve).
The methodology just described has been applied to the dataset of POD sTEC data recorded by Swarm A, B, and C from 2 December 2013 to 31 December 2020. In this way, it has been possible to identify a number of LoL equal to 16,146 for Swarm A, 12,593 for Swarm B, and 15,992 for Swarm C. From here on, we talk always in terms of the cumulative dataset of LoL obtained from all the three satellites of the Swarm constellation. Figure 2 shows the distribution of the duration (in seconds), between 1 and 100 s, of the identified LoL events. The percentage of events whose duration is greater than 100 s is less than 5.5%. The largest fraction of LoL is made up of events whose duration is of 18 s (about 29% of the total) and 1 s (about 24% of the total). The other bins are unevenly populated. This distribution is consistent with that found by Xiong et al. [14], who analyzed data from December 2013 to November 2016. We investigated the reliability of very short LoL events. Specifically, we tried to verify whether LoL of 1 s are real or most likely due to hardware outages. To this aim, we plotted in geographic coordinates the LoL occurrences within each bin, for the following two cases: (1) Considering all the 44,731 LoL (top panel of Figure 3); (2) Discarding the 10,727 LoL whose duration is 1 s (bottom panel of Figure 3). This comparison shows that discarding LoL of 1 s allows to eliminate most of the events located at mid latitudes, which are somewhat unexpected, not affecting at all the patterns characterizing both low and high latitudes. In virtue of this result, the analyses that is shown in Section 3 is based on LoL events whose duration is greater than 1 s.

Loss of Lock Occurrence Distribution
This section discusses the distribution of the identified GPS LoL events as a function of different variables to highlight their spatial, diurnal, seasonal, and solar activity features. Figure 3 shows the geographical distribution of LoL occurrences in bins that are 2.5 • -wide in latitude and 5.0 • -wide in longitude. Most of the LoL events are not randomly distributed, follow a well-defined structure, and maximize along the crests of the EIA at low latitudes, and inside the auroral oval at high latitudes, in both hemispheres.
LoL events at low latitudes have to be ascribed mainly to the generation of EPBs [6], a phenomenon that, due to the peculiar joint geometric configuration of the geomagnetic field and the equatorial zonal electric field, is more frequent at equatorial and low latitudes in the dusk-midnight sector, e.g., [29][30][31][32]. EPBs are characterized by different shapes and sizes, their lifetime ranges from a few minutes to several hours, and can seriously compromise the quality of radio signals, e.g., [33,34]. Recently, De Michelis et al. [24], using data from the ESA Swarm mission, have highlighted the turbulent nature of EPBs, suggesting at the same time a link between the scaling features of electron density fluctuations characterizing EPBs and RODI.
Considering the strong link between the equatorial plasma density irregularities and LoL events, we can suggest that the low-latitude maximum occurrence of LoL between about 70 • W and 10 • E of longitude has to be ascribed to both the longitudinal variation of the vertical plasma drift, and the alignment between the sunset terminator and the geomagnetic field [35]. The latter holds, taking into account that, at low latitudes, most of the plasma irregularities form in the evening sector. There is in fact a general agreement about the higher occurrence of EPBs locations where the dusk terminator is parallel to the magnetic meridian [36], a condition met in this particular longitudinal sector. Tsunoda [37] was the first who suggested that this longitudinal variation could be explained in terms of the flux tube integrated Pedersen conductivity gradients near the dusk terminator. According to Tsunoda [37], flux tubes with declinations aligning with the solar terminator facilitate the onset of the Rayleigh-Taylor (RT) instability and the onset of turbulent irregularities.
Concerning the longitudinal variability of the vertical drift, it is worth highlighting that Abdu et al. [38] and Saito and Maruyama [39] showed that there are specific longitudinal sectors, like the Southeast Asian ones, where thermospheric meridional trans-equatorial winds decrease the effect of the pre-reversal enhancement (PRE), with a consequent suppression of EPBs.
There are also studies that, in order to justify the different low-latitude longitudinal occurrence of EPBs, suggest the importance of two additional phenomena: atmospheric gravity waves that would modulate the ionospheric instabilities trigger mechanisms, e.g., [40,41]; and the South Atlantic Anomaly that would enhance the vertical drift due to the low magnetic field [42].
With regard to the high latitudes, from a longitudinal point of view, the distribution of LoL is more uniform than that at low latitudes, and confined inside the auroral oval. The auroral oval is in fact associated with significant particle precipitations and field-aligned currents from the magnetosphere, both providing the necessary energy for the development of ionospheric irregularities [43][44][45][46].
The more uniform longitudinal high-latitude distribution of LoL, if compared to the low-latitude one, has to be ascribed to the fact that the types of irregularities characterizing the high latitudes are many (polar cap patches, auroral blobs, the equatorward wall of the ionospheric trough, polar holes, and auroral cavities), located everywhere (in the dayside cusp, in the polar cap, and in the nightside auroral oval), and due to different plasma processes (the Kelvin-Helmholtz instability, the gradient drift instability, the large-scale plasma convection, the high-speed plasma convection, and the ion frictional heating) [23,[47][48][49][50]. The maximum occurrence from 70 • E to 180 • E of geographic longitude characterizing the Southern hemisphere (SH) can be understood recalling that the geographical coordinates of the South magnetic pole are about 64.3 • S and 136.6 • E. So, those bins showing a high occurrence of LoL have to be thought around the South magnetic pole and inside the polar cap. Figure 4 shows the distribution of LoL occurrences as a function of the Quasi-Dipole (QD) latitude [51] and the day of the year, in bins 2.5 • -wide in QD latitude and 5-wide in days. Again, it is clear how LoL events are symmetric with respect to the magnetic equator. At low latitudes they cluster around equinoxes, while at high latitudes their occurrence is bounded between about September and May. For all latitudes, the minimum occurrence of LoL is in June, July, and August, independently of the season. The reason of the clustering around equinoxes at low latitudes has to be found again in the alignment between the solar terminator and the geomagnetic flux tubes [35,37], a condition that is satisfied especially during equinoxes. Moreover, equinoxes are characterized by a more ionized ionosphere, a phenomenon which is known as the semiannual anomaly, e.g., [52]. This causes an increase of the equatorial zonal electric field and, around sunset, a reinforcement of the PRE, with a consequent strengthening of the vertical plasma drift. This creates a favorable environment for the RT instability, which leads to the development of EPBs, e.g., [53]. Fejer et al. [54] have done one of the first studies showing the positive correlation between the strength of PRE and the EPBs occurrence. Scherliess and Fejer [55], on the other hand, showed that eastward electric fields are more intense near the west coast of South America during March and April, which supports also the geographical distribution displayed in Figure 3. Recent works by Cai et al. [56] and Karan et al. [57], based on Global-scale Observation of Limb and Disk mission data, have given interesting insights about the South American sector in terms of the EIA hemispheric asymmetry and EPB zonal drift velocities. Furthermore, additional works by Gentile et al. [58] and Lee et al. [59] support the clustering of LoL around equinoxes. The former, using Defense Meteorological Satellite Program observations, has shown that EPBs occur most frequently around equinoxes. The latter found that the most suitable conditions for the development of EPBs, due to the interlink between the PRE and the asymmetry of EIA, occur just in the equinoctial months.
The LoL distribution characterizing the high latitudes is consistent with the seasonal trend reflecting that of polar cap patches, both in the Northern hemisphere (NH), e.g., [21,60] and in the SH [61], and that of ground-based scintillations, e.g., [62,63]. In addition, Bjoland et al. [50], using data from the European Incoherent Scatter Svalbard Radar, found that electron density depletions seem to be clearer during equinox and winter, for moderate and high solar activity, than during the rest of the year. This seasonal trend can be explained in the NH in terms of the effect of solar EUV radiation that during summer ionizes the E region, which may significantly short-circuit the F-region irregularities. Instead, in the SH this explanation can no longer be considered valid to justify the clear damping of irregularities between June and August, and other mechanisms should be invoked. This non-seasonality characterizing the high-latitude LoL occurrence could be somehow linked to the ionospheric annual anomaly. The annual anomaly is related to the observation that, taking the world as a whole, the overall level of the maximum of electron density appears greater in December than in June, e.g., [64][65][66]. Although it has been known for more than 80 years [67], this electron density asymmetry is still unexplained. In principle, the simplest explanation of this anomaly would relate to the shape of the Earth's orbit and, in the first place, could be ascribed to the variation of the Sun-Earth distance. The electron density annual asymmetry between solstices is, however, of about 30% and greatly exceeds the ionization rate asymmetry due to the annual variation of Sun-Earth distance which is approximately equal to 7% [68]. This is a feature that needs further and more specific studies that are out of the scope of the paper. Figure 5 shows the distribution of LoL occurrences as a function of the QD latitude and Magnetic Local Time (MLT), in bins 2.5 • and 15 min wide, respectively, in QD latitude and MLT, both as global and polar projection. Specifically, the polar projections show the probability density, just to better highlight which are the regions where the LoL occurrence maximizes. As expected, at low latitudes the events cluster between 19 and 23 MLT, the time window for which EPBs are most likely to occur, e.g., [69,70].
Differently, at high latitudes the diurnal distribution of LoL is more homogeneous, with maxima located around the local noon and in the nightside sector, mainly between 18 and 00 MLT. This is rather expected because, contrary to low latitudes, at high latitudes the ionospheric irregularities do not show a preferential time window to occur. They take place indifferently in the sub-auroral regions, in the dayside cusp, in the polar cap, as well as in the nightside auroral oval and ionospheric trough, e.g., [43,50,71,72]. Maxima related to the local noon are linked to the irregularities characterizing the dayside cusp, e.g., [73,74], while those related to the nightside sector are related to the irregularities characterizing the nightside auroral oval dynamics, e.g., [74,75]. Figure 6 shows the LoL occurrences for each year separately, from 2014 to 2020, as a function of the QD latitude and the day of the year. The bottom panel of the figure also shows the daily values of solar index F 10.7 and the corresponding 81-days running mean. F 10.7 is the solar radio flux at 10.7 cm wavelength (2800 MHz) [76] and the corresponding data have been downloaded at the National Aeronautics and Space Administration OM-NIWeb Data Explorer website (https://omniweb.gsfc.nasa.gov/form/dx1.html, accessed on 31 May 2021). Figure 6 shows an evident dependence of the LoL occurrence on solar activity. In fact, LoL events maximize in 2014 and 2015, when solar activity is maximum, and they practically disappear afterwards. This confirms the findings of previous studies [35,61,77,78], namely, that ionospheric irregularities, independently of latitude, are directly proportional to solar activity and, thus, consequent phenomena such as LoL. On the other hand, for high solar activity, the physical processes at the base of the irregularities formation, like for instance the equatorial zonal electric field at low latitudes and particle precipitations at high latitudes, are much more intense, e.g., [53,79,80].

RODI and ROTI Values at GPS Loss of Lock
As previously shown in Figures 3-5, the LoL distribution presents a high-latitude asymmetry between North and South, with a higher occurrence characterizing the SH than the NH. The same asymmetry was observed analyzing RODI by Jin et al. [22] and De Michelis et al. [23], who found that this is more intense in the SH than in the NH. Jin et al. [22] also found that in the SH the distribution of high values of RODI extends over a larger area than in the NH, which reflects well the LoL distribution shown for instance in Figure 3. According to this, following the procedure described in Section 2.3, we show the binned median values of RODI and ROTI values calculated during LoL events.
The top panel of Figure 7 is equivalent to the bottom panel of Figure 3, but in this case, what is represented within each bin is no longer the LoL occurrence but the median of RODI values corresponding to all LoL events falling in that bin. The dataset, the binning procedure, as well as the graphical representation of Figure 7 are the same that have been applied to obtain Figure 3. The figure shows that low-latitude RODI values are higher than high-latitude ones. The values in a log 10 scale range between 5 and 6 (then, between 10 5 and 10 6 cm −3 /s) at low latitudes, and between 4.5 and 5.5 (then, between 10 4.5 and 10 5.5 cm −3 /s) at high latitudes. These values are rather consistent with those found by De Michelis et al. [24,81], who investigated the properties of electron density fluctuations in terms of the first-and second-order scaling exponents and of an intermittency coefficient. This suggests that at the base of LoL events most likely there are processes of turbulent nature. It is also important to point out that the RODI values found are by far higher than those, ranging in a log 10 scale between 2.9 and 3.8 (then, between 10 2.9 and 10 3.8 cm −3 /s), obtained considering all RODI values for the period December 2013-December 2020, and not only those corresponding to a LoL event, as highlighted by the bottom panel of Figure 7.   Figure 8, the median within each bin is calculated considering only ROTI values corresponding to LoL events and, differently from RODI, the obtained values are similar for low and high latitudes. Again, however, it is important to highlight how these values are by far higher than those displayed in the bottom panel of Figure 8, for which the median within each bin is calculated considering all ROTI values recorded from December 2013 to December 2020, and not only those corresponding to a LoL event. In a log 10 scale the values in the top panel of Figure 8 range between −1.1 and −0.3 (then, between 10 −1.1 and 10 −0.3 TECU/s), while those in the bottom panel of Figure 8 range between −2.1 and −1.9 (then, between 10 −2.1 and 10 −1.9 TECU/s). It is the intention of the authors to try to investigate also the properties of the TEC fluctuations, similarly to what was done for the electron density [24,77], to strengthen the assumption that at the base of a LoL event there are processes of turbulent nature.
Figures similar to Figure 8, but concerning seasonal and QD latitude-MLT distributions similarly to what has been shown in Figures 4 and 5 for LoL events, are visible in the Supplementary Materials as Figures S3 and S4. In the Supplementary Materials, the reader can also find corresponding RODI and ROTI versions of Figure 6 as Figures S5 and S6. Figures 7 and 8 highlight the fact that, in coincidence with the occurrence of a LoL event, there are high values of both RODI and ROTI, which are much higher than the background, as also highlighted by Figure 1. This means that a LoL event is most likely to occur inside regions characterized by large electron density gradients. From this perspective, with the aim of understanding if there is any value of RODI and ROTI at which it is possible to say that a LoL event has the greatest probability of occurrence, Figures 9-11 show the joint probability density distribution between the values of RODI and ROTI obtained in correspondence of LoL events. Figure 9 refers to all the LoL events recorded all over the globe, while Figures 10 and 11 consider only LoL events recorded at |QD latitude| > 45 • and |QD latitude| ≤ 45 • , respectively. Figure 9 shows that there is a well-defined family ranging between 4.5 and 5.5 values of RODI and between −0.5 and −1.0 values of ROTI (we recall that both scales are in log 10 ) for which most of the events tend to cluster. This family remains practically unchanged when considering LoL events occurred at |QD latitude| > 45 • , while it moves rightward when considering LoL events occurred at |QD latitude| ≤ 45 • . These results follow those recently found by De Michelis et al. [81], based only on RODI values, that have showed that high and mid-low latitudes are characterized by different families of RODI values, namely, by different kinds of physical processes responsible for the ionospheric irregularities formation. Further studies are needed, involving also additional data as the magnetic and electric field ones, to try to identify what are more specifically the processes behind the occurrence of a LoL event. This could pave the way for a future possible forecast of such events.

Summary and Conclusions
Considering data from the ESA Swarm constellation recorded from December 2013 to December 2020, we analyzed the occurrence of GPS LoL events and characterized them in terms of RODI and ROTI. This is the first time that such an analysis, based on such a large time window, is carried out globally, encompassing half a solar cycle, from the activity peak of 2013 to the minimum of 2020. The main outcomes of the study can be summarized as follows: 1.
LoL events are mainly located at low and high latitudes, for both hemispheres. Specifically, at low latitudes they maximize along the EIA crests between about 70 • W and 10 • E of longitude; 2.
The high-latitude LoL occurrence is higher in the Southern hemisphere than in the Northern hemisphere. Is this asymmetry real? This is a point that needs further analyses; 3.
At low latitudes, LoL events cluster around equinoxes; both high and low latitudes are characterized by a minimum of occurrence in June, July, and August, independently of the season. This non-seasonality of the LoL occurrence is intriguing and deserves additional analyses; 4.
At low latitudes, LoL events cluster between 19 and 23 MLT, while at high latitudes the diurnal distribution of LoL is more uniform, with maxima characterizing the local noon and the nighttime sector between 18 and 00 MLT; 5.
LoL events strongly depend on solar activity, maximizing for years of maximum solar activity and reducing significantly for years of minimum solar activity; 6.
LoL events are strictly connected with very high values of both RODI and ROTI, meaning that they are most likely to occur inside regions characterized by very large electron density gradients; 7.
Joint probability density distributions between RODI and ROTI values corresponding to LoL events showed that there is a well-defined family for which LoL events tend to cluster. Moreover, in terms of RODI, the values are somewhat consistent with those recently found by studies focused on the turbulent feature of ionospheric irregularities, which suggests that LoL events are likely caused by irregularities triggered by processes of turbulent nature. The position of this family changes moving from high to low latitudes. Does this movement imply that the physical processes responsible for the ionospheric irregularities formation at the base of LoL events are different? Moreover, is the second family visible for lower values of RODI (around 3) real? These questions deserve additional analyses.
Recalling that LoL events seriously degrade the accuracy of GNSS and GPS systems on which many of our infrastructures by now depend, the importance of these results can be appreciated. Trying to characterize these events in terms of ionospheric indices, as done in this work, could be valuable to lay the foundations for their possible future prediction, which would be of extreme importance in the framework of Space Weather effects mitigation. Data Availability Statement: Publicly available data were analyzed in this study. Elaborated data are available on request from the corresponding author.
Acknowledgments: Thanks to the European Space Agency for making Swarm data publicly available via ftp://swarm-diss.eo.esa.int (accessed on 31 May 2021), and for the considerable efforts made for the Langmuir probes and POD antennas data calibration and maintenance.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Acronyms
The