Effects of the 12 May 2021 Geomagnetic Storm on Georeferencing Precision

In this work, we present the positioning error analysis of the 12 May 2021 moderate geomagnetic storm. The storm happened during spring in the northern hemisphere (fall in the south). We selected 868 GNSS stations around the globe to study the ionospheric and the apparent position variations. We compared the day of the storm with the three previous days. The analysis shows the global impact of the storm. In the quiet days, 93% of the stations had 3D errors less than 10 cm, while during the storm, only 41% kept this level of accuracy. The higher impact was over the Up component. Although the stations have algorithms to correct ionospheric disturbances, the inaccuracies lasted for nine hours. The most severe effects on the positioning errors were noticed in the South American sector. More than 60% of the perturbed stations were located in this region. We also studied the effects produced by two other similar geomagnetic storms that occurred on 27 March 2017 and on 5 August 2019. The comparison of the storms shows that the effects on position inaccuracies are not directly deductible neither from the characteristics of geomagnetic storms nor from enhancement and/or variations of the ionospheric plasma.


Introduction
The tendency to incorporate autonomous or tele-operated systems in certain industries around the globe requires careful analysis over the performance of the technology that supports this activity. One of the key elements in those systems is the global navigation satellite system (GNSS), which allows the geolocation of a receiver almost anywhere in the world. However, GNSS as a radio-link-based system can be perturbed by variations in the ionosphere. In particular in Chile, mining, fishing, and agriculture activities are relevant economic sectors in Chile that are exploring the migration to a more teleoperated/autonomous operation. Due to this requirement, it is imperative to study the robustness of the GNSS systems along the country, in particular when geomagnetic storms are produced.
One of the most important drivers of ionospheric variations at a global scale is solar activity. During a geomagnetic storm, the ionosphere can be greatly perturbed through the Sun-solar wind-magnetosphere-ionosphere-atmosphere coupling process. Perturbations in the ionosphere, in particular electrons, can produce disturbances in space-based technology and their products such as communication disruptions and imprecision in the positioning estimation of global navigation satellite system (GNSS). The enhancement or oscillations in the ionospheric electron content influence the signal between a many underlying mechanisms needs a better understanding to obtain precise forecasting of its behavior based on geomagnetic storms parameters [20]. In addition, the effects of these physical processes on the ionosphere have also been reported to vary with solar activity, storm intensity, storm duration, season, location, local time, and altitude of the observing station, which increases the forecast uncertainties [6,[21][22][23][24][25][26][27]. For instance, Tsurutani et al. [28] and Mannucci et al. [29] have found that the response of the ionosphere to geomagnetic activity depends on the season of the year in which that portion of the ionosphere is located [30,31].
The ionospheric plasma density can be estimated by estimating the ionospheric electron density. Thus, it is the determining variable for investigation of the spatial and temporal variations in the ionosphere. The total electron content (TEC), which can be estimated from double-frequency GNSS receiver data, is used to study the ionospheric response during ionospheric storms [1,3,9,20,32,33]. The TEC is defined as the integral of the electron density from the ground height up to the ceiling height, i.e., the height of the transmitting satellite or infinity. Since the contribution to the TEC usually comes from low orbital altitudes (below 1000 km), then, vertical TEC (VTEC) is obtained from Slant TEC (STEC) at the ionospheric pierce point (IPP) at an altitude of 350 km. Another relevant ionospheric index used to study ionospheric variations is the rate of change of the TEC index (ROT). ROT is calculated in a temporal window. ROTI is defined as the standard deviation of ROT. It describes the small-scale variability of the line-of-sight electron content resulting from the ionosphere and plasmasphere. The total electron content (TEC) maps, together with other indexes derived from TEC, are used to estimate the locations and time where larger signal delays in GNSS receivers might produce higher positioning errors.
Although the position estimation of a GNSS receiver depends on TEC, for several reasons (location, correction algorithms, etc.), perturbations on TEC are not reflected in a simple manner over the position accuracy. Therefore, it is relevant to study the actual performance of the positioning estimation of GNSS receiver during geomagnetic storms. To evaluate the position error during geomagnetic storms, the precise point positioning (PPP) method is used. PPP is a method that performs efficient computation to determine high-quality coordinates. It uses a single receiver processing strategy for GNSS [50,51]. PPP does not require any additional data from a reference station and can provide a solution from a centimeter to decimeter level of positional accuracy both in static and kinematic modes [50,52]. For these reasons, PPP has become the predominant positioning technique [50]. Many free PPP online services are available, such as the Automatic Precise Positioning Service (APPS) of the Global Differential GPS (GDGPS) System (https://apps.gdgps.net/, accessed on 3 December 2021), the GNSS Analysis and Positioning Software (GAPS-PPP) (http://gaps.gge.unb.ca/, accessed on 3 December 2021), the magicGNSS solution (magicGNSS) (https://magicgnss.gmv.com/, accessed on 3 December 2021), and the Precise Point Positioning of Canadian Spatial Reference System (CSRS-PPP) (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php, accessed on 3 December 2021) [53][54][55]. The CSRS-PPP service is one of the most used PPP online services in the field, and we use this service in this work to obtain the PPP estimation.
More recent studies have presented TEC disturbances produced by geomagnetic storms including GNSS position errors (PPP). These latest studies of geomagnetic storms and errors in GNSS positioning have focused on storms of the Solar Cycle 24. A summary of these studies can be found in Table 1.  Table 1 shows that there are not many studies related to the impact over the GNSS positioning accuracy of moderate geomagnetic storms. In a recent study, it was shown that the storm of 5 August 2019 , which can be categorized as moderate (Dst peak = −53 nT, Kp = 5 + , AE∼1000 nT), had strong effects on TEC [61]. However, this work did not study the positioning accuracy of GNSS stations. The only study of the positioning error that considers a moderate geomagnetic storm was presented by Luo et al. [57] (see Table 1). Luo et al. analyzed three geomagnetic storms in Solar Cycle 24. For the analysis, they used 17 March 2017 as a reference day for all the analyzed storms, which is a quiet day (Kp = 1 + ). This work neither analyzes the conditions on the previous/posterior days for each storm nor removes the noise sources, such as other geophysical events (earthquakes) or interference from other radio emissions. In addition, in this work, the root means square (RMS) statistics per component (North, East, and Up) were obtained per latitude (high-, mid-, and low-). Thus, the reported 3D RMS was calculated with the components that combine different stations located at a similar latitude range.
In this work, we present the positioning error analysis of the 12 May 2021 geomagnetic storm (Dst peak = −61 nT, Kp = 7, AE∼1500 nT), which can be classified as moderate in terms of Dst but strong in terms of Kp. Although the study is at a global scale, unlike previous studies, we focus on the southern sector of America due to the large number of GNSS stations now available in this part of Latin America. We studied the spatial and temporal dependence of the higher errors estimated in the GNSS receivers for this storm. We also studied the time needed for some geodesic-quality GNSS stations to reduce the positioning error thanks to the algorithms that detect and correct the effects of TEC disturbances. By using three stations, one close to Madrid (Spain) and the other two in Chile (one close to Santiago and a second one 400 km south of it), we studied the effects on the positioning accuracy for two other moderate storms (27 March 2017 and 5 August 2019) to verify the methodology used in the 12 May 2021 study.

Materials and Methods
In this section, we describe the used data and processing procedure. The proposed procedure also includes analyzing the geophysical and geomagnetic conditions close to 12 May 2021.

Estimation of the Ionospheric Total Electron Content
The slant TEC (STEC) and vertical TEC (VTEC) data were obtained from GNSS measurements based on dual-frequency signals f 1 and f 2 . Then, STEC and VTEC were calculated using the program GPS-TEC from receiver independent exchange format (RINEX) files [62] (http://seemala.blogspot.com/, accessed on 3 December 2021). In this work, VTEC values corresponding to satellite cut-off elevation angle 30 o at 350 km altitude were selected to minimize possible errors. The temporal VTEC estimates were released every 30 s. The TEC values were corrected from the satellite and receiver bias using the data obtained from AIUB Data Center of Bern University in Switzerland (ftp://ftp.aiub.unibe.ch/CODE/, accessed on 3 December 2021).
The RINEX files were obtained from 868 GNSS stations (Figure 1), taked from: the International GNSS service (IGS) stations; the Chilean network of GNSS receivers operated by the National Seismological Center at University of Chile (CSN in Spanish); the Argentine Continuous Satellite Monitoring Network (RAMSAC in Spanish) [63]; the Brazilian Network for Continuous Monitoring of the Institute of Brazilian Geography and Statistics (IGBE in Portuguese); and the Crustal Dynamics Data Information System (CDDIS) of the National Aeronautics and Space Administration (NASA); University NAVSTAR Consortium (UNAVCO); the Geoscience Australia; and the African Geodetic Reference Frame (AFREF). In addition, the differential of VTEC (DVTE) in TECu and the percentage changes in VTEC (%DVTEC) were studied. These parameters are used in the analysis of ionospheric disturbances, defined as the relative variation of VTEC, epoch by epoch, with respect to the mean value (in time) of VTEC as shown in Equations (1) and (2) [3].
where t represents the epoch and VTEC is calculated by averaging the values of VTEC for the reference DoYs 129, 130, and 131.

Ionospheric TEC Maps
We used the ordinary Kriging interpolation technique [64] to map the estimation of the VTEC at each ionospheric pierce point (IPP). We selected this interpolation technique to fill in the data gap of the global ionosphere TEC map and the inhomogeneous sparsity of GNSS receivers [65][66][67]. We performed the ordinary Kriging interpolation in Python with the package Kriging. The documentation for the package can be found in https://github.com/ERSSLE/ordinary_kriging (accessed on 3 December 2021).

ROT and ROTI
Rate of change of the TEC Index (ROTI) is defined as the standard deviation of the rate of TEC (ROT). It is used in the detection and investigation of occurrences of ionospheric irregularities. ROTI was estimated by dual-frequency GNSS data with the time interval of 5 min by using Equation (3) [68][69][70][71][72]: where · represents the temporal average. The ROT and ROTI values are typically expressed in units of TECu/min. ROT is defined as the TEC variation rate of two successive epochs as stated in Equation (4) [ [68][69][70][71]: where k, t, and i represent the GPS time, the epoch, and the number of observation satellites, respectively. According to Liu et al. [73], the ROTI value can be divided into three groups: weak, if 0.25 ROTI < 0.5; moderate, if 0.5 ROTI < 1; and strong, if ROTI 1.

Apparent Position Variation Using Precise Point Positioning
The RINEX files of 868 GNSS stations were processed in the postprocessing kinematic precise point positioning with ambiguity resolution (PPP-AR) mode in the Precise Point Positioning of Canadian Spatial Reference System (CSRS-PPP) online service (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php, accessed on 3 December 2021) [53]. The CSRS-PPP provides centimeter-level estimations with converged float solutions [53][54][55]74]. The CSRS-PPP returns to the user a processing report via email.
Usually, the report provides a different reference value for different days. We process the data to have an equal reference for all the used data to facilitate the evaluation of the position-variation time series. At each of the 868 stations, we also apply the common noise filter to correct the time series of the North, East, and Up components (Equation (5)). Then, we estimate the 3D resultant.
where t is the epoch, CAPdoy t is the corrected apparent position, APdoy t is the apparent uncorrected position, and RP t is the reference position. We use AP from reference DoYs 129, 130, and 131 to calculate RP t . At each of the 868 stations, the maximum error was calculated for eight days and with time windows of 15 min. Subsequently, the percentage of stations in East, North, Up, and 3D position errors intervals was calculated.
On 12 May 2021 (DoY 132), a moderate geomagnetic storm (G3) took place with a Kp-index 4 from 6 UT. The maximum Kp-index of 7 was sustained between 12 and 15 UT. The Kp-index went below 3 after the 18 UT.
The geomagnetic storm started on the DoY 132. Its initial phase started at ∼6 UT and lasted ∼6 h. During this period, the Dst-index slowly decreased. The main phase of the storm started at ∼12.4 UT and lasted for ∼one hour. In the storm main phase, the Dst-peak of −61 nT was reached at 14 UT. At that time, the geomagnetic storm recovery phase began. The Dst-index went over −30 nT at ∼23 UT on the DoY 132 (see Figure 2b). We were unable to obtain the raw data of the AE-index, but it is possible to observe the behavior of this index in the graphical display of the data on the website of the WDC for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/ae_realtime/202105/index_20210512.html, accessed on 3 December 2021). In that image, it is possible to see an intensification of the AE-index with two peaks ∼1500 nT between 12 and 16 UT. The IMF Bz was below −10 nT between 12 and 15 UT on the day of the storm (DoY 132). The minimum Bz of −18.3 nT was reached at 13 UT. According to the Gonzalez and Tsurutani criteria [2], this event can be classified as a geomagnetic storm since it was caused by an interplanetary magnetic field Bz ≤ −10 nT that lasted more than 3 h (see Figure 2c). Tsurutani criteria [2], this event can be classified as a geomagnetic storm since it was caused by an interplanetary magnetic field Bz ≤ −10 nT that lasted more than 3 h (see Figure 2c).
The geomagnetic conditions were generally quiet between 9 and 11 of May 2021 (DoYs 129 and 131), where the Dst-index was predominantly positive, and the maximum Kp-index was 2 (see Figure 2a).

Possible Earthquakes Perturbations
We reviewed the occurrence of earthquakes (EQs) around the world with moment magnitudes over 5 Mw and depth less than 70 km between the 9 and 17 of May 2021 (DoYs 129 to 137). The data were obtained from https://earthquake.usgs.gov/earthquakes/ search/ (accessed on 3 December 2021). Five strong EQs (6-6.9 Mw) and twenty moderate EQs (5-5.9 Mw) occurred in this analyzed period. However, none of them produced noticeable effects on TEC during the geomagnetic storm period. The differential TEC was analyzed, taking care of the potential minor effects in the other days of the analyzed period.

Results
In this section, we present the main results obtained after applying the methodology described in the previous section to the 868 stations.
First, we isolated the period of time to analyze. We focused the study between 9 May and 17 May 2021 (DoY 129 to 137). Figure 2 presents an example of the signals analyzed for this work. We calculated VTEC, DVTEC, %DVTEC, and ROTI as described in the previous section. These variables can be seen in panels (g-j) of Figure 2, respectively, for a sample station, the Vegas de Itata station (known as vita station) located in Chile (36.42 • S, 72.86 • W). Then, by using the VTEC and the PPP-AR method (see Section 2.3), we were able to obtain the apparent position variation. In panels (d) to (f) in Figure 2, we present the root mean square (RMS) time series of the apparent position, after correcting using the common noise filter (see Equation (5)), per each component (East, North, and Up) also for the vita station. It is clear in the image that the period of the storm is the period with the larger uncertainties in the position estimation.
From each station data, we can estimate the VTEC over each station during the period of the geomagnetic storm (14.50-14.75 UT) on 12 May 2021 as shown in Figure 3a. By using the ordinary Kriging interpolation, as described in Section 2.1, it is possible to obtain filled-in VTEC maps for the geomagnetic storm day (Figure 3c) and for the averaged VTEC (VTEC). VTEC is calculated with the VTEC of the previous days to the geomagnetic storm, DoYs 129, 130, and 131 ( Figure 3b). Figure 3d,e show the differential VTEC (DVTEC) and perceptual DVTEC (%DVTEC) for the time of the geomagnetic storm, respectively, which are obtained by using Equations (1) and (2) with the data presented in (Figure 3b,c). From these figures, it is possible to notice that the South American sector is one of the most affected in terms of VTEC enhancement.
Since this study focuses on the positioning accuracy of the stations during a moderated storm, we estimate the PPP-AR as described in Section 2.3. Figure 4  After calculating the apparent position variation for each station in the time period of the geomagnetic storm (DoY 132, 14.50-14.75 UT), we quantified the number of stations around the world that had errors that fell in certain intervals (<5 cm, [5][6][7][8][9][10] cm, [10][11][12][13][14][15][16][17][18][19][20] cm,  cm,  cm,  cm and >100 cm). It is important to highlight that the PPP-AR procedure includes the use of a common noise filter (see Equation (5)). Without the filter, we had periods of time where certain GNSS stations consistently had very high errors even during quiet days. The results of the classification of the stations by their positioning errors is presented in Table 2. Each column on the table represents the percentage of stations with errors in a certain interval per each positioning component as well as the combination of these components in the 3D parameter. Each component uses two contiguous columns in Table 2, one with the percentages obtained for the average of the previous quiet days (e.g., North) and another with the percentage measured during the day and period of the geomagnetic storm (e.g., North gs). In this table, the percentage of the stations per interval for East, North, Up, and 3D position errors were obtained using the data from the total 868 GNSS stations around the globe we had available.
In Table 2, we present a snapshot of the errors focused on the geomagnetic storm time. Nevertheless, we also present a temporal evolution of the percentage of the stations with positioning errors over a certain threshold per component for the total of the available stations around the world in Figure 6. From this figure, it is possible to notice that the main errors are concentrated over the mid-day of the DoY 132, which is the geomagnetic storm period.
It is possible to notice from the TEC (e.g., Figure 3) and the derived PPP-AR (Figures 4 (left panel), and 5 (left panel)) data that there is a region where the errors are more severe during the geomagnetic storm (DoY 132, 14.50-14.75 UT). This region is South America. In South America are located 325 stations of the 868 total available GNSS stations (∼37%). In Table 3, we present similar results compared to in Table 2 but with the 325 stations of this region. However, there is no other localized area in which we can detect strong variations. The other perturbed stations are distributed around the world. For this reason, we concentrated this study in the South American sector.
The ROT index (ROTI) was calculated as described in Section 2.2 to study the relation of the variation of TEC in the positioning error for the 12 May 2021 geomagnetic storm. The images in the Figure 4 [57,61]. The August 2019 geomagnetic storm [61] was a moderate one that produced a extreme positive ionospheric storm with strong TEC perturbations over Europe. However, this work did not study positioning errors. Thus, we selected a station in Europe to be compared with the two Chilean stations, which are located in South America and in the area of interest for this work. In addition, the stations in Chile are both required since each of them was operative for different storms (see Table 4). Luo et al. [57] studied the position error produced by the moderated geomagnetic storm of March 2017 but in broad latitudinal ranges, calculating the maximum positioning errors per component for each of these ranges but not per station.  Table 4 summarized the positioning variations at each of the three stations produced by three different geomagnetic storms.

Discussion and Main Conclusions
In this section, we discuss the main findings regarding the geomagnetic storm of 12 May 2021.The main goal is to study the positioning errors of GNSS receivers caused by this storm, which can be classified as a moderate one. In addition, in order to verify our methodology and our findings we used previous equivalent geomagnetic storms as comparison. The comparison storms were the 27 March 2017 [61] and the 5 August 2019 [57].

12 May 2021 Geomagnetic Storm
The geomagnetic storm under study occurred in 12 May 2021. It can be classified as moderated in terms of Dst and strong in terms of Kp (Dst = −61 nT, Kp = 7, and AE = ∼1500 nT). This storm is the first strong storm, in terms of Kp, of the solar cycle 25. The main phase of this storm lasted one hour from 13 UT to 14 UT, with the Dst peak at 14 UT (see Figure 2). Since the storm occurred in May the southern and northern hemispheres were in fall and spring seasons, respectively. During the main phase of the storm the Sun moved from the (18.

The 27 March 2017 and 5 August 2019 Geomagnetic Storms
In Table 1 On the other hand, a recent study presented the case of a moderated geomagnetic storm occurring on 5 August 2019 (Dst peak = −53 nT, Kp = 5 + , AE ∼1000 nT) which produced a strong positive (decrease in electrons) ionospheric storm [61]. Although this geomagnetic storm was less intense than the May 2021 storm, its main phase was much longer (8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21). During this storm, the southern hemisphere was in winter, then, like in May 2021, the Sun was over the northern hemisphere (from (17.01 • N, 61.52 • E) to (16. 86 • N, 133.50 • W)). Thus, we used this storm as a comparison to the May 2021 storm. Unfortunately, the study of this storm did not analyze the positioning error caused by it. The study showed that the European sector was one of the most affected in terms of TEC variations. For this reason, the comparison was performed including a GNSS located in Europe (close to Madrid, Spain). The other two stations located in Chile were used, since that region was one of the most perturbed areas in terms of position during the May 2021 storm (see Figure 4g plot).

Ionospheric Effects
From Figure 3, we can see that the VTEC is perturbed around the globe, including the southern part of America during the main phase of the May 2021 geomagnetic storm. On the other hand, in Figure 4h, we can see an increment in ROTI, that starts at the polar regions, propagating later the increment toward the equator, agreeing with previous studies [77,78] (see Figure 5). In Figure 7, it is possible to see the VTEC and DVTEC perturbations of the GNSS stations over Madrid, Spain and Santiago, Chile (Vizcacha station) for the May 2021 storm. In comparison, we also present the VTEC and DVTEC perturbations on the same stations produced by the August 2019 storm. It is noticeable that TEC variations are similar between storms for the same station. In addition, the difference in the variations between the stations in southern and northern hemispheres is evident. The TEC variation seems stronger in the northern station. In Table 4, we can see that the VTEC is in percentage much higher over Madrid (Spain) than over Santiago (Chile) and Vega de Itata (Chile) for most of the storms, except for the March 2017 storm. We also calculated the ROT index for the May 2021 storm period. Five stations localized in high latitudes had values of maximum ROTI peak between 3.4 and 3.9 TECu/min, and average ROTI peak ∼1.1 TECu/min, and these values are similar to those presented by Kotulak et al. [77] for moderate geomagnetic storms.

Positioning Errors
The manner in which we quantified the positioning errors in this work was through the statistics of perturbed stations around the world and in particular in the South American sector. Our RMS position values for the quiet days are in accordance with the results presented by Katsigianni et al. [79]. They gathered that the performances of the kinematic postprocessed PPP-AR method are 0.8 and 2 cm for the horizontal components and the vertical component, respectively. Table 2 shows that for the quiet days, 71% of the stations had an error less than 5 cm in the 3D estimation. Table 2 also shows that for the May 2021 geomagnetic storm the main increment suffered by the Up component, passing from 27% of the total stations (868) with perturbations over 5 cm to 48%. Although the impact over the North and East components was less affected, they jumped from 1% to 13% and from 2% to 15%, respectively, for the stations with errors over 5 cm.  Figure 5g shows that the positioning errors were perceived mainly in stations in the South American sector. Table 3 shows that in the Up component, we went from 45% of the total station over South America (325) with an error over 5 cm to 87%. The similar increment can be perceived for the other components. In the North component, we passed from 2% to 28% of stations over an error of 5 cm. Similarly in the East component, the variation went from 2% to 30% for the same error range. Only  We compared the effects obtained over the positioning error due to the May 2021 storm, with the effects of the moderate geomagnetic storms that occurred on 27 March 2017 and 5 August 2019. Figure 7 compares the effects on the position in two different GNSS stations, one in Chile (South America sector) and one in Europe, where a previous storm produced strong ionospheric effects during the 5 August 2019 geomagnetic storm [16]. From the comparison, it is possible to notice that the perturbations in the position are much higher in the Vizcacha (Santiago, Chile) station and for the May 2021 storm. This is an intriguing result. Both storms, May 2021 and August 2019, are very similar, except for the duration of the main phases, but still, they produced very different effects on the position estimation (see Figure 7). Even more puzzling is the fact that the August 2019 storm has been reported to produce a strong ionospheric storm that affected the European sector [61]. For instance, the %DVTEC is almost double that in the Madrid station compared to the Santiago station (Vizcachas station), but we have found with our analysis on one station in the European sector that it did not impact the position estimation. Table 4 presents the positioning error and %DVTEC data for the three selected for comparison stations (one in Europe and two in South America) during three moderate geomagnetic storms. It is possible to notice the same pattern as that obtained between the May 2021 and August 2019 storms which was obtained between the March 2017 and the May 2021 storms. Thus, the effects over the position estimation were very severe during the May 2021 storm over the South America sector (Vega de Itata, vita-station), even though the %DVTEC was much higher for the March 2017 storm both for Spain (Madrid) and Chile (vita) stations.
The ROTI also has an intriguing behavior. We observe an increase in the number of stations with ROTI > 0.25 TECu/min ( Figure 5) especially in the poles and the South American region. In North America and northern Europe region, the number of stations with activity increased from 50% to 94%. In the Antarctic region, there was an increase from 12% to 62% of stations with activity. In South America, the percentage passed from 1.5% to 10.5% of the stations with activity. Nevertheless, we do not find a significant increase in the number of stations with higher positioning errors over the poles. Therefore, it is conclusive that the fast variations of TEC might be responsible for the variations over the South American sector.
Our results did not confirm that positioning errors increased rapidly with increasing ROTI (Figures 2d-f,j, 4g,h, 5 and 7; Table 4) as several studies suggest [56][57][58]60]. Therefore, our results suggested that position errors also occur, regardless of whether the ROTI has rapid variations, if it is ROTI < 0.25 TECu/min or no ROTI variations are appreciated.
In summary, the analyzed data show that the moderate 12 May 2021 geomagnetic storm strongly affected the GNSS precision for about an hour in several GNSS receivers, mainly on stations located in the South American sector. Comparison with previous moderate storms (27 March 2017 and 5 August 2019) showed that the effects in the position estimation are not directly deducible from the geomagnetic storms characteristics. The three analyzed storms were moderated, although the 12 May 2021 storm had a higher Kp, which tends to be a good indicator of midlatitude activity. By using three stations we compare the effects of the three storms, showing that the effects are stronger on the South American sector even though the ionospheric effects (DVTEC) are not severe.
The positioning error data show that the horizontal coordinates are more robust to TEC perturbations, although the error in the vertical component is still high. In the literature, the vertical coordinate tends to be neglected, but for current and future GNSS applications, it might be relevant. It could be important for autonomous aerial applications or for high-precision activities. The mining, agriculture, fishing, and disaster-control sectors in Chile are starting to adopt autonomous or tele-operated systems that might be sensitive to ionopsheric perturbations. The relevant height of the vehicles in these industries might impose a serious risk for infrastructure or people if the vertical error in the GNSS receiver skyrockets. For instance, for open-pit mines, a high vertical error may cause a failure in the estimation of the terrace in which a vehicle is with the consequent risk of falling. Potential risks could be reduced by stopping autonomous operations during these events. However, our results show that it is hard to predict when a storm will have serious effects over the position accuracy. False positives in the forecasting can be complex for these industries, since stopping operation even for a short time, such as an hour, could be prohibitively expensive. The new solar cycle is critical to unravel the relation between solar activity and positioning errors in the GNSS receivers thanks to the large number of GNSS stations available around the world, gathering relevant data. However, we need to process the data in a standardized manner that facilitates comparison between different cycles, paying attention to the storms that produce unexpected behaviors such as those in the 12 May 2021 geomagnetic storm.

Conflicts of Interest:
The authors declare no conflict of interest.