Technical note on the multi-GNSS, multi-frequency and near real-time ionospheric TEC monitoring system for South America

Development of regional services able to provide ionospheric total electron content (TEC) maps with a high spatial resolution, and in near real-time, are of high importance for applications and the research community. We provide here the methodologies, and a preliminary assessment, of such a system. The system relies on the public Global Navigational Satellite Systems (GNSS) infrastructure in South America, it incorporates data from multiple constellations (currently GPS, GLONASS, Galileo and BeiDou), it employs multiple frequencies, and it produces continental wide TEC maps with a latency of just few minutes. A year-round comparison of the produced maps with several products issued by the International GNSS Service (IGS) resulted in mean biases lower than 1 TEC units (TECU), whereas their evaluation against direct and independent GNSS-based slant TEC measurements resulted in accuracies of the same magnitude. 1 System Description


Data Preprocessing
Before any computation a per station preprocessing and data cleaning is performed. This includes the application of an appropriate time window, of an elevation cut off angle, of carrier phase wind-up corrections and the determination of phase-continuous intervals (i.e., with constant ambiguity). The data cleaning is non-parametric and consists in three steps. Firstly, pairs of satellite-receiver phase and code links, in two bands, are optimally selected according to the amount of available observations and the tracking modes, or channel attribute, employing the same default priorities as the ones defined by http://rtgpsout.unavco.org † Registration may be required in order to access same of the data streams.
Nischan [2016]. The corresponding, undifferenced, Melbourne-Wübbena (MW) linear combinations are screened for outliers and cycle slips ( Table 2). Any possible receiver's clock inconsistency between phase and code observations is handled as a cycle slip. Secondly, a per band, time-differenced, phase screening is performed, for unnoticed outliers and cycle slips. These screenings, performed both for-and backward in time, are repeated until no additional cycle slips, or outliers, are found. Furthermore, no attempt is made to correct any cycle slip. Finally, MW and (cuasi) ionosphere-free (IF) linear combinations, within each phase-continuous interval, are formed and modeled with low degree polynomials. Intervals resulting in residuals with root mean squared (RMS) greater than given thresholds are rejected. In total, less than 7 % of the original observations are generally left out, including those observations bellow the elevation cut off angle.
Thereafter, a clean and single set of undifferenced carrier phase φ ij,k (in cycles), code pseudorange C ij,k (in meters) and signal-to-noise ratio SN ij,k observables, between each pair of satellite i and receiver j, and for each tracked band k, is obtained. This data preprocessing is performed with the Fortran 2008 + OpenMP, in-house developed, software AGEO (library for Geodetic and Orbital Analysis or biblioteca de Análisis GEodésico y Orbital, in Spanish). In addition, it is externally parallelized, on a per station basis, by means of the GNU parallel software tool [Tange, 2011]. Table 2: Pairs of bands employed in the Melbourne-Wübbena (MW) data screening, with their corresponding wide-lane wavelengths λ WL (in meters). Those pairs of bands employed in the computation of geometry-free (GF) linear combinations and inter-frequency biases (IFBs) are also indicated. Observation codes according to RINEX version 3.03 [Gurtner and Estey, 2017]

Hardware Delays Calibration
Once per hour inter-frequency biases (IFBs) are estimated from carrier-to-code leveled geometry-free (GF) linear combinations [see, e.g., Spits, 2012]. These hardware delays are solved simultaneously with spherical harmonic (SH) coefficients of a single-layer VTEC representation [Schaer, 1999], assuming that all free electrons are constrained to an infinitesimally thin layer at a height of 450 km. Here only independent linear combinations are employed, between selected pairs of bands (Table 2). Hence, only independent IFBs are computed. Also, no closing restriction is imposed, resulting in satellite-receiver-, pair-of-bandsspecific IFBs. Moreover, different tracking modes or channel attributes, between each pair of satellitereceiver links, are taken into account. The estimation is made by means of a weighted least squares adjustment, performed also with the AGEO software, and executed in one single step, involving the most recent observations available, from all ground stations and all satellites, within the previews 24 hours (i.e., a 24 hours rolling-or moving-window).
In practice, after preprocessing the raw observations corresponding to each pair satellite i and receiver j, and for each pair of bands k and l, all possible carrier-to-code leveled GF linear combinations L GF,ij,kl (in meters) are computed by where L GF,ij,kl = λ k φ ij,k − λ l φ ij,l are the non-leveled GF linear combinations in phase (in meters) and C GF,ij,kl = C ij,l − C ij,k are the corresponding linear combinations in code pseudorange (in meters), being λ k and λ l the wavelengths of each band (in meters). Here the average is computed within each phasecontinuous interval, under the assumption of stable hardware delays. This commonly used methodology reduces the observations noise, from code to phase levels, and avoids the estimation of phase ambiguities, but it could also introduce some systematic errors [see, e.g., Ciraolo et al., 2007, Spits, 2012. In addition, each GF observation is weighted according to three factors: the instantaneous satellite elevation, the amount of observations employed during the carrier-to-code leveling (i.e., the length of each phasecontinuous interval) and the corresponding navigational system [GPS, GLONASS, Galileo or BeiDou, see Ren et al., 2016]. Then, the full set of GF observations is represented as where MF(z) is the Modified Single Layer Model (MSLM) mapping function [Schaer, 1999], being z the zenith distance of satellite i as seen by receiver j (in radians), IFB ij,kl are the corresponding specific IFBs (in meters), α kl is a proportionality constant (in meters per TECU, where 1 TECU is equivalent to 10 16 free e − per squared meter) being f k and f l the frequencies of each band (in Hertz), whereas the VTEC is expressed as a SH expansion in a sun-fixed frame VTEC(µ, t) = nmax n=0 n m=0 P nm (sin µ) a nm cos(m t) + b nm sin(m t) .
Here a nm and b nm are the coefficients of the SH expansion (in TECU), with maximum degree n max , whereas P nm are the corresponding Real Associated Legendre Functions [4π normalized, see for example Wieczorek and Meschede, 2018], t is the Local Time (LT, in radians) and µ is the modified dip latitude (also in radians). In this case the algorithms issued by the European GNSS (Galileo) Open Service [2016], together with the corresponding global grid, are employed for the computation of µ. As only regional observations are employed, the 24 hours time window helps into decoupling the hardware delays, from the ionospheric parameters, by always including in the adjustment observations spanning 24 hours of LT. Furthermore, a single set of constant coefficients a nm and b nm , loosely constrained to a zero ionosphere, is estimated for the entire time span, resulting in a mean (daily) VTEC representation. For the same reason a SH expansion of low degree is employed, in order to avoid ill conditioned normal equations [Haines, 1985], which in turn could produce mapping artifacts, particularly at the boundaries of the region. On the other hand, the IFB ij,kl are also parametrized as constants, and estimates with mean observational epoch at the middle of each moving-window are obtained. These hardware delays are also loosely constrained to their most recently estimated values. In fact, the main result of this hourly adjustment are precisely these decoupled IFB ij,kl estimates. However, they result systematically half a day old.

Near Real-Time TEC Mapping
Every 15 minutes, and also by means of a weighted least squares adjustment, both the IFBs and the SH coefficients for the regional VTEC representation are updated. In essence, the same software, methodology and parametrization described in the previews section are employed. However, in this case only the most recent observations available, within a one hour moving-window, are used. In addition, here four sets of pice-wise constant SH coefficients are estimated, each one valid for a quarter of an hour. Furthermore, the IFB ij,kl parameters are now actively constrained to their most recent, hourly, and decoupled estimates. In practice, this results in new hardware delays estimates that are simultaneously up-to-date (i.e., less than 30 minutes old) and decoupled from the coefficients of the SH expansion. Thereafter, tracks of instantaneous, VTEC estimates are obtained from the original GF observations by where ϕ and λ are the geographic latitude and longitude, respectively, of the ionospheric pierce points (IPPs), that is, the intersection point between the instantaneous satellite-receiver line-of-sight with the single layer of the model. Similarly, traces of slant TEC (STEC) estimates can be computed by At this point two representations of the current state of the regional ionospheric TEC are available. In one hand, an analytical representation, given by the coefficients of a low degree SH expansion in µ and t, with mean epoch at the middle of the latest 15 minutes of the observational window. On the other hand, a discreet and huge set of instantaneous VTEC ij,kl,ϕλ estimates, along the IPP tracks, during the same interval. In fact, the issued TEC product is obtained by mapping these tracks in the space domain. This postprocessing of the VTEC ij,kl,ϕλ estimates comprises three steps, all performed with the Generic Mapping Tools software package [GMT, Wessel et al., 2013]. Firstly, all available estimates are averaged within the cells of a uniform 0.5 × 0.5 degrees grid, previously discarding cells with very few observations, and effectively resulting in N space-and time-averaged VTEC p values, for p = 1, . . . , N . Secondly, this regular grid is approximated, using a generalized Green's function for continuous curvature spherical spline in tension [Wessel and Becker, 2008], by where ϕ and λ are arbitrary coordinates, M ≤ N is the number of employed coefficients, ϕ p and λ p are the coordinates of the corresponding cells, c 0 is the mean VTEC over all populated cells (in TECU), g is the generalized Green's function and c p are the spline coefficients (also in TECU), solved for by Singular-Value Decomposition (SVD) on the square linear system and retaining only those M eigenvalues whose ratios, to the largest, are greater than a given threshold. While we empirically determined optimal (fixed) values for both the tension and the threshold, searching over thousand of maps for minimization of the misfits, the number M of contributing eigenvalues is dynamically determined, every time, to accommodate the variance of the current data. That is, the more spatial variability in the regional VTEC the more eingenvalues are retained in the mapping procedure. Finally, the adjusted function is evaluated on a uniform 0.5×0.5 degrees grid and areas far away from IPP tracks are automatically masked out. The resulting grid constitutes the actual, near real-time, regional TEC map produced by the system, as no additional postprocessing (e.g., smoothing) is required nor performed.

Year-round Intercomparison with Global Products
To evaluate the quality of the produced maps, and particularly the possible presence of systematic biases, we compared them with several IGS final TEC products, provided in IONEX format, and computed  -Pajares et al., 1999. From UPC we employed both, their standard and their high rate products. We also included in the analysis TEC products from two additional IGS ACs: Natural Resources Canada [NRCAN, Canada; see Ghoddousi-Fard et al., 2011] and Wuhan University [WHU, China; see . These products are usually available with latencies of a few days or, at best, several hours. In addition, we also included in the intercomparison the (non-IGS) global and high resolution TEC products provided by the Massachusetts Institute of Technology [MIT Haystack Observatory, USA; see Rideout and Coster, 2006]. The comparison extends a full year, from June 1, 2017 to May 31, 2018, and it was performed on a map by map basis (i.e., epoch by epoch). In order to assess the expected differences we performed the same one-to-one comparison between pairs of global products. Although these maps have global coverage, the comparison was restricted to the area covered by our regional maps, that is, between 80 • S and 40 • N in latitude and 110 • W and 0 • E in longitude [see Figure 1 in Mendoza et al., 2019]. Also, no spatial or temporal interpolation was performed. Rather, only TEC samples at common epochs, and exactly the same reported locations, were differenced. For this reason, and before the comparisons, our high resolution maps were downsampled. Thus, the results were controlled by the standard 5 × 2.5 degrees spatial sampling (in longitude and latitude, respectively) of the IGS products or, alternatively, by the 1 × 1 degrees spatial sampling of the MIT products.
For this analysis, instead of the real-time data streams, we employed daily observational and navigational RINEX files available at the servers of the respective data providers. However, to reproduce exactly the results of the near real-time system, we only used data from those GNSS stations that are actually accessible in real-time, leaving all off-line stations out of the analysis. We also employed the very same broadcasted orbits and satellite clocks, and no other products. In addition, we followed exactly the same two-steps methodology previously described. That is, a first step resulted in IFBs estimates, from a 24 hours observational moving-window, while in a second and final step the TEC maps were produced, from a 15 minutes moving-window. To speed up this year-round analysis, and although some of the selected products are currently provided at a higher rate (e.g.,by CODE, NRCAN and particularly UPC and MIT), we computed maps with 2 hours of temporal sampling, following the classical IGS standard practice.
The year-round (and regional) comparison shows the existence of systematic differences between all the analyzed pairs of TEC products (Table 3). In average, our near real-time TEC maps show a very good agreement with the maps produced by ESA/ESOC, CODE, UPC (high rate) and especially NRCAN, resulting for all the cases in a mean bias lower or equivalent to 1 TECU, comparable with the differences found between pairs of global products. Nevertheless, while there seems to be no significantly biases between both, our products and the ones from NRCAN and between them and the ones from CODE, a small systematic bias do exists between the former and our products. The reason for this seemingly discordant results is simple: given their global coverage, the comparisons between pairs of IGS products span the entire area mapped [ Figure 1 in Mendoza et al., 2019], whereas those comparisons involving our regional product are mostly restricted to the land, leaving large portions of the oceans out of the analysis, and this is evident in the lower number of common TEC samples found (Table 3). This contributes also to the higher mean standard deviation encountered while comparing our maps with the other products, both in average and individually (Figure 1). Indeed, not only a smaller number of differences are averaged, also the smoothest areas of the IGS maps over the oceans, where no actual GNSS observations were available, are systematically left out of these comparisons. At the same time, all comparisons show, to a greater or lesser extent, smaller variance during the southern winter (i.e., June, July and August). This is probably due to the lower, regional, mean ionospheric TEC in that season.

Differential STEC Evaluation
In order to independently assess the accuracy of the produced TEC maps we applied a differential STEC (dSTEC) test developed by the IGS Ionosphere Working Group (IIWG) for the evaluation, and relative weighting, of their Global Ionosphere Maps (GIMs) [see, for example, , 2007, Roma-Dollase et al., 2018. In essence, the test is based on the ability to make highly accurate dSTEC measurements, on the order of 10 −2 TECU [Coster et al., 2013, Hernández-Pajares et al., 2017, and to compare them with synthetic (i.e., mapped) dSTEC values. In fact, we employed the very same Table 3: Year-round one-to-one comparison between selected (GNSS based) TEC products, from June 1, 2017 to May 31, 2018: codg (CODE), emrg (NRCAN), esag (ESA/ESOC), igsg (IGS, combination of codg and jplg), jplg (JPL/NASA), mapgps (MIT), upcg (UPC), uqrg (UPC, high rate), whug (WHU) and magn (our near real-time product). The mean difference x A−B and mean standard deviation σ A−B , over all compared maps, are expressed in TECU.
Products  Figure 1: Examples of mean differences and standard deviations, per map, resulting from the year-round one-to-one TEC products comparisons: between a final IGS and our near real-time product (codg and magn, respectively), between two final IGS products (codg and the high rate uqrg), between two final IGS products (codg and upcg, noting that no upcg IONEX files were available for October 8-13 and 21, 2017 and for January 28, 2018) and between a final and the combined IGS product (codg and igsg, respectively). The global K p index, provided by the GeoForschungsZentrum (GFZ, Germany), is also plotted.
implementation of the test as described in detail by Hernández-Pajares et al. [2017], the only difference being our extension of the test to the multi-frequency case. The analysis is performed on a per station basis, involving only GNSS stations that were not employed for the computation of the TEC maps being evaluated. Firstly, and after preprocessing the corresponding raw data (e.g., outliers rejection, phase wind-up correction, etc.), observed dSTEC o (in TECU) are obtained from (non-leveled) carrier phase GF linear combinations (in meters) by taking advantage of the total cancellation of the phase ambiguities within each phase-continuous interval.
Here t r (in hours) represents a reference epoch, when the satellite reaches its minimum zenith distance within each phase-continuous interval, whereas t s (in hours) are all other sample epochs within the same phase interval. Here we employed a sampling rate of 60 seconds and, following the convention stated by Hernández-Pajares et al. [2017], only GF observations no more than 900 seconds apart were differenced. This results in a maximum of 30 dSTEC o samples per phase-continuous interval, regardless of its total length. In addition, the corresponding zenith distances z r and z s (both in radians), with z r = z s , are stored for subsequent use. Secondly, and for each observed dSTEC o sample, synthetic dSTEC m values (in TECU) are computed by dSTEC m (t s ) = MF(z s ) VTEC(ϕ s , λ s , t s ) − MF(z r ) VTEC(ϕ r , λ r , t r ) where ϕ r , λ r and ϕ s , λ s (in degrees) are the coordinates of the corresponding IPPs. Here both VTEC(ϕ, λ, t) are obtained, following Schaer and Feltens [1998], by temporal interpolation between consecutive rotated TEC maps being T i and T i+1 the epochs of the corresponding maps (in hours), with T i < t < T i+1 , whereas the rotated longitudes λ i = λ + 15 (t − T i ) and λ i+1 = λ + 15 (t − T i+1 ) compensate the strong correlation between the ionospheric TEC and the (longitude of the) subsolar point. Within each map, the VTEC i and VTEC i+1 are spatially interpolated by a simple 4-point bilinear algorithm [see also Schaer and Feltens, 1998]. Finally, the observed minus computed ∆dSTEC (in TECU) are obtained In turn, the RMS of ∆dSTEC, per station, can be computed.  Table 4, and employed for the dSTEC evaluation of the near real-time TEC maps. For convenience the geomagnetic equator is also plotted.
In practice, we employed daily RINEX files from ten GNSS stations distributed over the study area ( Figure 2). In addition to files from the mentioned data providers we also employed observations from off-line GNSS stations supplied by the Centro Sismológico Nacional (CSN, Chile). The analysis was repeated in four independent days, during the years 2017 and 2018, near the ascending equinox, the descending equinox, the summer solstice and the winter solstice. Also, the TEC maps employed in this analysis are the very same produced for the year-round comparison with the IGS GIMs. However, for these four particular days, additional maps were produced in order to achieve the standard 15 minutes sampling rate of the monitoring system.
In summary, the observed dSTECs are fairly reproduced by the synthetic values, implying that the near real-time maps, in combination with the corresponding mapping function, are capable of representing the regional ionospheric TEC with an average accuracy better than 1 TECU (Table 4). Finally, the three stations leading to a total RMS > 0.7 TECU are located in areas where the IPPs coverage is, systematically, not optimal [especially near BOAV and PUMO, but to a lesser extent also near PISR, see Figure 1 in Mendoza et al., 2019], right). This suggests that the monitoring system could benefit from the use of additional data from GNSS stations located in these specific areas.

Conclusions and Outlook
A multi-GNSS, operational, high-rate and openly accessible ionospheric TEC monitoring system for South America has been successfully developed, tested and implemented [Mendoza et al., 2019]. Both the comparison of the produced maps against global products, including several final IGS GIMs, and also against independent and highly accurate dSTEC observations, resulted in no significant biases nor in unreasonable variances or RMS. On the oder hand, the monitoring system described here could certainly be improved in several ways. For example, the employment of observations at a higher rate could improve the data cleaning stage, resulting in less discarded observations. Data with a higher sampling rate could also open the possibility to reliably repair the cycle-slips, resulting in longer phase-continuous intervals, which in turn could reduce the systematic biases introduced during the carrier-to-code leveling process. Furthermore, the final spatial TEC mapping could be performed in a sun-fixed frame, or it could be parametrized in modified dip latitude instead of geographic latitude, or it could employ a different basis of analytical functions, just to mention some possible improvements. Moreover, while we are actually employing the classical double-frequency approach, to retrieve the ionospheric delay from the GNSS observations, a true triple-frequency methodology [see, e.g., Spits, 2012] could also be considered for those particular stations equipped with modern receivers. Finally, the (rather small) systematic biases found between our TEC maps and the ones produced by several IGS ACs should be addressed, looking at the differences in mathematical representation and parametrization and their impact on the resulting TEC and IFBs estimates. Comparisons with TEC measurements provided by other techniques or satellite missions could certainly help to identify the sources of these differences.
The near real-time TEC monitoring system described in this article is currently running on equipment acquired with funding by the Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, Argentina), trough the program Fondo para la Investigación Científica y Tecnológica (FONCyT), with grant No. PICT-2015PICT- -1776