Operational Multi‐GNSS Global Ionosphere Maps at GFZ Derived From Uncombined Code and Phase Observations

Global navigation satellite system (GNSS) networks with multi‐frequency data can be used to monitor the activity of the Earth's ionosphere and to generate global maps of the vertical total electron content (VTEC). This article introduces and evaluates operational GFZ VTEC maps. The processing is based on a rigorous least squares approach using uncombined code and phase observations, and does not entail leveling techniques. A single‐layer model with a spherical harmonic VTEC representation is used. The solutions are generated in a daily post‐processing mode with GPS, GLONASS, and Galileo data, and are provided for the period since the beginning of 2000. A comparison of the GFZ VTEC maps with the final combined International GNSS Service (IGS) product shows a high consistency with the solutions of the IGS analysis centers. A validation with about four years of Jason‐3 altimetry‐derived VTEC data is provided, in which the GFZ solution has the smallest bias of 1.2 TEC units compared to the solutions of the IGS analysis centers, and with 3.0 TEC units one of the smallest standard deviations.


Introduction
Global navigation satellite systems (GNSS) are an effective tool to observe the Earth's ionosphere and measure its total electron content (TEC; e.g., Davies & Hartmann, 1997;Hernández-Pajares et al., 1999;Jakowski et al., 1999;Komjathy, 1997;Lanyi & Roth, 1988;Mannucci et al., 1998;Sardon et al., 1994;Schaer et al., 1996). GNSS-derived TEC data provide valuable information for atmospheric studies or can be used to assess and mitigate the impact of the ionosphere on GNSS-based services such as positioning or tropospheric sounding.
Due to the presence of the instrumental code biases and the phase ambiguities, the unbiased line-of-sight or slant TEC (STEC) of the receiver-satellite links cannot be directly observed with GNSS, but only through applying an external model describing the spatiotemporal TEC distribution. It is common practice to use the geometry-free or L4 combination of the GPS L1 and L2 signals, which represent the STEC biased by respectively code biases and ambiguities for code and phase observations (Schaer, 1999). The ambiguities are usually removed from the model through leveling the phase to the code observations (Ciraolo et al., 2007;Mannucci et al., 1998). The so obtained leveled phase observations are then used to subsequently estimate the model coefficients and a set of estimable receiver and satellite code biases. The same technique can also be applied to dual-frequency data of other systems. The disadvantages of this estimation strategy are that (a) it is limited to two frequencies, (b) it reduces the redundancy of the system, as one code and one phase observation are used in order to remove only one geometry parameter, and (c) it generally does not apply a weighting of the observations caused for instance by different elevation angles, see also the discussion in Khodabandeh and Teunissen (2016). Further, the temporal correlations introduced by the leveling operation seem to be generally ignored.
An alternative strategy for TEC estimation is to only use the geometry-free phase combinations and to estimate an ambiguity term for each receiver-satellite arc Hernández-Pajares et al., 1999;Krypiak-Gregorczyk & Wielgosz, 2018). The idea behind this strategy is to avoid leveling errors (Ciraolo et al., 2007) at the cost of additional ambiguity parameters. A least squares framework for multi-frequency TEC determination using uncombined code and phase observations of a single receiver or a local array of receivers avoiding all the above-mentioned shortcomings is presented in Teunissen (2016, 2017).
Abstract Global navigation satellite system (GNSS) networks with multi-frequency data can be used to monitor the activity of the Earth's ionosphere and to generate global maps of the vertical total electron content (VTEC). This article introduces and evaluates operational GFZ VTEC maps. The processing is based on a rigorous least squares approach using uncombined code and phase observations, and does not entail leveling techniques. A single-layer model with a spherical harmonic VTEC representation is used. The solutions are generated in a daily post-processing mode with GPS, GLONASS, and Galileo data, and are provided for the period since the beginning of 2000. A comparison of the GFZ VTEC maps with the final combined International GNSS Service (IGS) product shows a high consistency with the solutions of the IGS analysis centers. A validation with about four years of Jason-3 altimetry-derived VTEC data is provided, in which the GFZ solution has the smallest bias of 1.2 TEC units compared to the solutions of the IGS analysis centers, and with 3.0 TEC units one of the smallest standard deviations.
Global ionosphere maps (GIMs), containing a snapshot of the TEC in vertical direction (VTEC) on a spatial grid, are routinely generated by different institutes, for instance within the ionosphere working group of the International GNSS Service (IGS; see Hernández-Pajares et al., 2009or Roma-Dollase et al., 2018 for an overview).
The German Research Centre for Geosciences (GFZ) Potsdam offers a wide range of products for ground and space-based atmospheric sounding (Wickert et al., 2020). These efforts are now extended with operational ground-based GNSS GIMs, which are introduced and evaluated in this article. The main focus and novelty of the article is on the formulation of a rigorous least squares approach for the computation of multi-GNSS global VTEC maps based on the work in Teunissen (2016, 2017). The common single-layer ionospheric model with a spherical harmonic VTEC parameterization is employed. The GIMs are computed from GPS, GLONASS, and Galileo data for the period since 2000 and will be provided on an operational basis in the future. A high consistency with the GIMs of other analysis centers and the combined IGS solution is obtained. In a comparison with about four years of Jason-3 VTEC data, the GFZ solution has the smallest VTEC bias and one of the smallest standard deviations among the solutions by the IGS analysis centers.
In Section 2, the estimation framework and the ionospheric modeling approach are introduced. The processing settings and an analysis of the GFZ GIM product are presented in Section 3. A validation by means of a comparison with the IGS solution and Jason-3 VTEC data is given in Sections 4 and 5.

Methodology
In this section, the estimation framework for the GFZ GIMs is presented. We start with the GNSS observation equations and estimable parameters in Section 2.1 and the ionospheric modeling strategy in Section 2.2. The combination of these two parts yielding the unbiased TEC solutions is presented in Section 2.3. s r E d , which is only possible through an ionospheric model, such as discussed in Hernández-Pajares et al. (2011). We use the common single-layer model (Azpilicueta et al., 2006;, 2010Davies & Hartmann, 1997;Goss et al., 2019;Mannucci et al., 1998;Schaer, 1999), in which the TEC is assumed to be contained within an infinitesimal thin shell at a constant height, which we assume at 450km E above the mean Earth radius in agreement with the solutions provided within the IGS. The slant ionospheric delays ( )   This single-layer spherical harmonic model representation is a strong simplification of the physical reality and, therefore, a limiting factor for the quality of the TEC estimates. A higher resolution could for instance be obtained with regional models in areas with dense station networks. In order to reduce the error resulting from the constant thin layer height, the use of two and four layers has been proposed in Hernández-Pajares et al. (1997,2020).

Unbiased TEC Estimation
Estimating the unbiased STEC or, equivalently, the unbiased ionospheric slant delays ( ) be a vector containing the known values of the mapping function times the basis functions from Equations 4 and 5, and let ( ) c be a vector of the same size with the associated coefficients , ( ) w t exponential elevation dependent noise amplification factors from Euler and Goad (1991),  2 p E the zenith-referenced variance of the code observations, and E  a factor that accounts for the higher precision of the phase observations and is chosen as  0.0001 E  . Correlations between observations are assumed absent.
Through an invertible transformation the E T code observations , ( ) s r j E p t can equivalently be expressed by the time-averaged component: see Khodabandeh and Teunissen (2016). The same transformation is applied to the phase observations , the full-rank time-averaged and time-differenced observation equations according to Table 1 are given by: and with F E e an E F -vector of ones, The solution of the time-averaged ionospheric delay  ( ) s r E i t and its variance are given by: while the solution of the time-differenced ionospheric delays  i t r s ( ) 1 and their variances and covariances are given by: . While the precision of the time-averaged component   i t r s ( ) is driven by the code data, the time-differenced ionospheric delays   i t r s ( ) 1 are on the level of the phase data. Derivations are given in the Appendix and in Teunissen (2016, 2017).
The time-averaged and time-differenced estimable ionospheric delays are expressed and parameterized as: Since only the time-averaged ionospheric delay  ( ) s r E i t depends on the link-specific GF code bias ,GF is a free variate that does not contribute to the solution of the ionospheric model coefficients  E c and can simply be ignored (Teunissen, 2000), thereby also eliminating the bias parameter ,GF  Khodabandeh and Teunissen (2017).
The ionospheric observables   i t r s ( ) 1 of different arcs are uncorrelated, so that the normal equations from all arcs are added, where within each arc the inverse covariance matrix from Equations 16 and 17 is used for weighting.
We now consider the model with the receiver and satellite specific code biases, cf.
where the bias of the first receiver is chosen as reference to remove the inherent rank-deficiency. A different yet equivalent set of estimable biases is obtained with the commonly used zero-mean constraint on the GF satellite code biases of each constellation.
In the dual frequency case, the estimable parameters in Tables 1 and 2 are identical, so that ionospheric observables   i t r s ( ) can be obtained as described above. If more than two frequencies are used, the parame- and can no longer be ignored, except for the single-receiver case. We

Operational Global VTEC Maps
The above-described estimation strategy is implemented in GFZ's EPOS-P8 GNSS analysis software. Global VTEC maps with a temporal resolution of two hours, at even hours, and a spatial resolution of 2.5° in latitude and 5° in longitude are provided in the ionosphere map exchange format (IONEX, Schaer et al., 1998). They are generated for the period since the beginning of 2000, and will be provided on an operational basis in the future. Daily rapid solutions are combined with the two neighboring solutions on the normal equation level to a final solution, from which mainly the first and last maps that are contained twice benefit. The settings of the processing are in line with the GFZ activities for the third IGS reprocessing campaign (Männel et al., 2020), so that GLONASS is included since 2012 and Galileo since 2014. Dual-frequency GPS L1/ L2, GLONASS L1/L2, and Galileo E1/E5a data are used. The numbers of stations and satellites are shown in Figure 1. The decay of the number of stations from 2015 is caused by stations that become unavailable, while at the same time not as many new stations are added to the processing. This is a consequence of using an almost identical network definition as for the IGS reprocessing campaign, in which recently installed stations with short coordinate time series are given low priority. In 2021, around 250 globally distributed stations and 75 satellites are employed. The number of stations will be kept at this level in the future.
Daily receiver and satellite GF code biases/DCBs are estimated and provided for GPS to be compatible with the IGS solution. For GLONASS, the receiver-satellite decomposition of the DCBs is not valid due to the presence of inter-channel code biases caused by the FDMA scheme, so that the receiver DCBs are satellite dependent (Wang et al., 2016;X. Zhang et al., 2017). Furthermore, the impact of the receiver-satellite bias decomposition on the VTEC maps is negligible, as shown below. The biases of GLONASS and Galileo are therefore assumed arc-specific and DCBs are not estimated, cf. Section 2. This is conform with the definition of the future IONEX 1.1 format, in which DCBs are no longer supported, but should be provided in a separate bias file. To this end, DCBs between arbitrary signals of any system, in particular also for the ones that are not used to generate the ionospheric solution, can be estimated after correcting the ionospheric delays of  s r E i in Tables 1 and 2 Figure 3. The values strongly resemble the characteristics of the solar radio flux, as shown by the F10.7 index (Tapping, 2013). The daily equivalent planetary amplitude Ap as an index of the average level of geomagnetic activity   . With a perfect model, the estimated variance factor should be below one. In the time series in Figure 4, we see values of  5 40 E for the square root of the variance factor, meaning that the phase residuals are at the decimeter rather than the millimeter or centimeter level that would be possible from the observation precision. The close resemblance with the mean VTEC in Figure 3 shows that with increasing ionospheric activity also the modeling errors increase. The maximum corresponds to the Halloween solar storm at the end of October 2003 with one of the largest recorded solar flares, for which also a peak of the F10.7 and an Ap index of around 200 is observed, see Figure 3.
In order to assess the impact of assuming arc-specific code biases or the receiver-satellite bias decomposition, the year 2010 was processed with both models. The daily absolute biases between the global mean VTEC of both solutions are mostly below 4 10 TECU E and the RMS VTEC difference is around 0.006TECU E as shown in Figure 5. Both values are clearly below the resolution of 0.1TECU E of the IONEX format, so that the benefit of the bias decomposition on the VTEC estimation through increased redundancy is negligible. This shows that essentially only the time-differenced observation data is relevant for estimating the ionospheric model, cf. Section 2.3, when applying a rigorous least squares approach.
A description and comparison of the final GIMs of these seven IAACs can be found in Roma-Dollase et al. (2018).
The smoothed time series of the daily VTEC biases and standard deviations between the individual solutions and the combined solution IGSG are shown in Figure 6, starting from November 2002, which is when IGSG switched from 12 daily maps at odd hours to 13 daily maps at even hours. The biases and standard deviations are computed by means of integrating over the sphere as defined above. One has to be careful when interpreting these measures, as they merely indicate how close the solutions are to the IGS combination, but not how well they describe the reality. In addition, between end-2014 and mid-2019, the IGS combination is, with a few exceptions, exclusively based on CODG and JPLG, so that these two solutions are by design close to the combination.
The daily VTEC biases in Figure 6 show that all eight solutions generally agree very well with the IGS solution with absolute values of less than two TECU. The bias of each solution is largely consistent over time and does not depend on the ionospheric activity. The VTEC values of JPL are about two TECU larger compared to the other solutions, whereas the VTEC values of GFZ tend to be slightly smaller than the IGS combination. The daily VTEC standard deviations confirm an agreement of the individual solutions at the few TECU levels, with values of around two TECU or below during times of low ionospheric activity, cf. Figure 3, and a few TECU during high ionospheric activity. The GFZ solution is consistent with the solutions of the IAACs and the IGS combination.

Validation With Jason-3 Altimetry VTEC
The VTEC over the oceans can directly be observed from dual-frequency altimeters (5.3 and 13.6GHz E for the Jason satellites). Although the altimetry VTEC observations potentially suffer from instrumental biases , and do not capture the TEC of the upper ionosphere and plasmasphere above the orbital height ( 1.336km E for Jason-3), they are GNSS-independent and are well suited to validate the GNSS GIMs. Exemplary Jason-2 VTEC values for all satellite passes on the day of the St. Patrick's Day geomagnetic storm corresponding to Figure 2 are shown in Figure 7. We use the VTEC observations of the Jason-3 satellite that are provided by DGFI-TUM (Dettmering et al., 2011) from its launch in January 2016 until January 2020. Observations with active rain or ice flag are removed. The number of daily observations used to validate the GNSS GIMs is shown in Figure 8.  For every altimetry VTEC observation, the associated GNSS VTEC is computed through a temporal interpolation of the two closest VTEC maps, which are rotated to account for the correlation between the TEC and the Sun's position, and a spatial interpolation between the four closest grid points (Schaer et al., 1998).
Daily mean biases and standard deviations of the GNSS solutions compared to the Jason-3 VTEC are shown in Figure 9. All GNSS solutions show a positive VTEC bias, which is attributed to the situation mentioned above. The differences of the biases reflect the behavior of the VTEC biases in Figure 6. The GFZ solution has the smallest bias, whereas the bias of the JPL solution is the largest one. Due to the expected and observed bias of the Jason-3 VTEC, the more important measure to assess the quality of the GNSS GIMs is the standard deviation. The UPC solution shows the smallest values. The results of the GFZ, CAS, CODE, and IGS solutions are almost identical, and the remaining solutions have slightly larger standard deviations with the ESA solution having the largest values. Comparing Figure 9 to Figure 3 confirms that the modeling errors increase with increasing ionospheric activity. The agreement of the GNSS GIMs with the Jason-3 VTEC is at the few TECU levels.
The overall distribution of all VTEC differences between the GNSS GIMs and Jason-3 is presented in Figure 10 for the nine GNSS solutions. Each curve is normalized so as to represent a valid distribution function. The associated biases, standard deviations, and RMS differences are given in Table 3. These results confirm the findings of Figure 9. The GFZ solution has the smallest bias. The UPC solution is most peaked, but all solutions have a standard deviation of around three TECU. In terms of the RMS values, the GFZ solution is the best, caused by the smallest bias and one of the smallest standard deviations.

Conclusion
In this article, the operational GFZ GNSS GIMs were introduced and an initial assessment was provided.
A rigorous least squares approach was formulated to estimate the coefficients of the single-layer spherical harmonic ionospheric model from uncombined multi-frequency, multi-GNSS code and phase observations. The shortcomings of the commonly used geometry-free observations combined with a phase-to-code leveling procedure were avoided with this approach. In order to avoid a large-scale adjustment problem, a twostep procedure was formulated. Biased ionospheric delays and potentially code biases for the third frequency and beyond are determined on an arc-by-arc basis, for which closed form solutions exist, and serve as an input for estimating the ionospheric model. This solution is mathematically identical to a direct solution.
If the code biases are assumed arc-specific, only the time-differenced observation data are relevant, whereas with a receiver-satellite bias decomposition also the time-averaged observation data are used. It was demonstrated, however, that the impact of this bias decomposition is negligible for the estimation of the global ionospheric model, so that the time-averaged data only become relevant if one wants to provide a set of DCB parameters along with the ionospheric solution, and can be ignored otherwise.
The GFZ GIMs were demonstrated to be consistent with the solutions of the IGS IAACs and the combined IGS solution through a comparison of daily VTEC biases and standard deviations between the individual solutions and the combined solution for a time span of more than 18 yr.
An altimetry validation with around four years of Jason-3 VTEC data confirmed the quality of the GFZ GIMs, which had the smallest bias and one of the smallest standard deviations compared to the solutions of the IAACs.

Appendix A: Solution of the Time-Averaged System
In . These phase observations are therefore free variates that do not contribute to the solution of the remaining parameters (Teunissen, 2000) and are omitted from the model. The resulting system with E F unknowns and E F equations is solved via inversion as:  yields:

Data Availability Statement
The GFZ GIMs are available at ftp://isdcftp.gfz-potsdam.de/gnss/products/iono . The GNSS observation data are provided by the IGS (Johnston et al., 2017) at https://www.igs.org/data. The satellite orbits are taken from GFZ's contribution to the third IGS reprocessing campaign , and GFZ's sensor meta-information system (Bradke, 2020) is used for the station and satellite metadata. The solar flux F10.7 index is provided by the National Research Council of Canada at https://www.spaceweather.gc.ca, and the Ap index by GFZ Potsdam . The altimetry VTEC data are produced by DGFI-TUM (Dettmering et al., 2011) and distributed at https://openadb.dgfi.tum.de.