Document [original]

1. 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 tropo-

spheric 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 re-

spectively 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 opera-

tion seem to be generally ignored.

An alternative strategy for TEC estimation is to only use the geometry-free phase combinations and to es-

timate an ambiguity term for each receiver-satellite arc (Brunini & Azpilicueta,2009; 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 mul-

ti-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 Khodabandeh and

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.

BRACK ET AL.

This is an open access article under

the terms of the Creative Commons

Attribution License, which permits use,

distribution and reproduction in any

medium, provided the original work is

properly cited.

Operational Multi-GNSS Global Ionosphere Maps at GFZ

Derived From Uncombined Code and Phase Observations

Andreas Brack1 , Benjamin Männel1 , Jens Wickert1,2 , and Harald Schuh1,2

1GFZ German Research Centre for Geosciences, Potsdam, Germany, 2Technische Universität Berlin, Institute of

Geodesy and Geoinformation Science, Berlin, Germany

Key Points:

• Operational GFZ global vertical total

electron content (VTEC) maps from

ground global navigation satellite

system (GNSS) data are introduced

• A rigorous least squares approach

with uncombined multi-GNSS,

multi-frequency code and phase data

is applied

• A validation by means of a

comparison with the solution of

the International GNSS Service

and Jason-3 altimetry VTEC data is

provided

Correspondence to:

A. Brack,

[email protected]

Citation:

Brack, A., Männel, B., Wickert, J.,

& Schuh, H. (2021). Operational

multi-GNSS global ionosphere maps

at GFZ derived from uncombined

code and phase observations. Radio

Science, 56, e2021RS007337. https://doi.

org/10.1029/2021RS007337

Received 16 JUL 2021

Accepted 2 SEP 2021

10.1029/2021RS007337

RESEARCH ARTICLE

1 of 14

Radio Science

BRACK ET AL.

10.1029/2021RS007337

2 of 14

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.,2009 or 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 opera-

tional 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 Khodabandeh and Teunissen(2016,2017). The com-

mon 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 pro-

cessing 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.

2. Methodology

In this section, the estimation framework for the GFZ GIMs is presented. We start with the GNSS obser-

vation equations and estimable parameters in Section2.1 and the ionospheric modeling strategy in Sec-

tion2.2. The combination of these two parts yielding the unbiased TEC solutions is presented in Section2.3.

2.1. GNSS Observation Equations and Estimable Parameters

The single-system GNSS code and phase observations

and



between receiver

and satellite

on fre-

quency

can be modeled as:

p t t it d d

rj j

, ,,

() () ()

() (







  















)) () ,





it a

(1)

where

refers to the time of reception. A geometry-free representation is used, in which no information on

the satellite-to-receiver geometry is employed. The use of a geometry-based or precise-point-positioning

model for ionospheric estimation is discussed in B. Zhang(2016). The link specific range parameters



contain all non-dispersive components such as the geometric ranges, the satellite and receiver clock offsets,

and the tropospheric line-of-sight delays, which cannot be estimated separately. The first order ionospheric

delays on the first frequency

are denoted by

and are multiplied with



1

. The receiver and sat-

ellite code biases are denoted as

,rj

and

, and the ambiguity parameter

 

 

, ,, ,

s ss

rj rj j j rj

Ea N

consists

of the receiver and satellite phase biases



,rj

and



and the integer ambiguity

, with



the wavelength

of frequency

. The integer ambiguities

are constant for each cycle-slip free arc, and the instrumental

code and phase biases are commonly assumed time-invariant (Hauschild,2017; Montenbruck & Haus-

child,2013). Therefore, the time dependency is omitted for

,rj

, and

. Phase wind-up, phase center

offset and variation, and other frequency specific effects that cannot be captured in



are corrected a priori.

We consider a network of receivers that observe the satellites on the frequencies

 1, ,Ej F

, with the

objective to retrieve information about the state of the ionosphere. The system of equations resulting from

Equation1 is rank deficient, implying that only certain linear combinations of the parameters can be unbi-

asedly estimated. In particular, the ionospheric slant delays

()

E it

cannot be separated from the code biases

,rj

and

and ambiguities

and can therefore not be determined in an absolute sense. We make use of

the invertible ionosphere-free (IF) and geometry-free (GF) decomposition of the code biases on the first two

frequencies (Khodabandeh & Teunissen,2015), shown for the receiver biases:





, ,IF ,GF ,2

rj r j r

dd d j

(2)

Radio Science

BRACK ET AL.

10.1029/2021RS007337

3 of 14

with

d dd

r rr

, ,,









21 12



(3)

The GF biases

,GFr

and

are usually referred to as differential code

biases (DCBs, e.g., Schaer,1999). The rank deficiencies can for instance

be removed by lumping the IF and GF code biases with



()

and

()

E it

respectively. The resulting estimable parameter combinations are given

in Tables1 and2. In Table1, the code biases

 

, ,,

rj rj j

Ed dd

are assumed

constant for each continuous satellite-receiver arc, so that receiver and satellite biases cannot be separated.

The estimable ionospheric delays

()

E it

are biased by the GF code biases. In Table2, the code biases are as-

sumed constant for the entire processing session, so that now receiver and satellite specific biases



,rj

and



are estimable, which increases the redundancy. The biases of the first receiver are chosen as reference

to remove the rank deficiencies. For now, this decomposition only affects the definition of the estimable

code biases



,rj

and



, but there are further implications when introducing an ionospheric model, so that

receiver and satellite GF code biases become estimable, see Section2.3.

Infinitely many other estimable combinations of the parameters can be formulated, leading to different

ionospheric parameters

()

E it

that can be estimated with different precision (Khodabandeh & Teunis-

sen,2016,2017). S-system theory tells us that there is a one-to-one correspondence between the parameters

of different estimable parameter sets, so that solutions from different choices can be transformed to one

another (Baarda,1973; Teunissen,1985) and can be used equivalently.

The extension to multi-GNSS is straightforward, since almost all parameters in Tables1 and2 are link spe-

cific. The only exception are the receiver and satellite code biases



,rj

and



in Table2, which are set up per

constellation, as the receiver biases are not identical across systems.

2.2. Global Ionospheric Modeling

In order to obtain estimates of the unbiased TEC, the ionospheric delays

()

E it

have to be separated from the

GF code biases

,GF

, which is only possible through an ionospheric model, such as discussed in Hernán-

dez-Pajares etal.(2011). We use the common single-layer model (Azpilicueta etal.,2006; Brunini & Az-

pilicueta,2009,2010; Davies & 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

450kmE

above the mean Earth radius in agreement with the solutions provided within the IGS.

The slant ionospheric delays

()

E it

are assumed to be related to the vertical ionospheric delays

()

E vt

in radial

direction via:





( ) ( ) ( ),

s ss

r rr

i t ME t v t

(4)

with the modified single-layer mapping-function



EM

(Schaer, 1999) depending on the elevation angle

()

E Et

. The global spatiotemporal distribution of the vertical ionospheric delays

()

E vt

is modeled with spher-

ical harmonic basis functions as (Schaer,1999):

Parameter Definition Condition

Ranges



()



,IF

()

E td

Ionospheric delays

()

E it

,GF

()

E it d

Code biases







, ,IF ,GF

ss s

rj r j r

Edd d

2Ej

Phase ambiguities

,





, ,IF ,GF

ss s

rj r j r

Ead d

Table 1

Estimable Parameter Combinations for Arc-Specific Code Biases

Parameter Definition Condition

Ranges



()





,IF IF

()

E td d

Ionospheric delays

()

E it



,GF GF

()

E it d d

Rec. code biases



,rj





1, 1,IF 1,GFrj r j r

Ed d d

1,Er



2Ej

Sat. code biases







, IF GF

ss s

E dd d



 

1, 1, IF 1,GFjj

E dd d

j > 2

Phase ambiguities

,



  

, ,IF IF ,GF GF

( )( )

ss s

rj r j r

Ea d d d d

Table 2

Estimable Parameter Combinations With a Receiver-Satellite Code Bias Decomposition;

   

() () ()

Radio Science

BRACK ET AL.

10.1029/2021RS007337

4 of 14



 



 



 





max

, () , ()

( ) sin ( ) cos ( ) sin ( ) ,

s sss

r nm r nm k t r nm k t r

vt P t c m t s m t

(5)

where

()

is the normalized associated Legendre function of degree

and order

, with a maximum de-

gree of



max 15En

. The latitude and longitude of the ionospheric pierce point (IPP) are denoted as



()

and



()

and are expressed in a solar-geomagnetic reference frame, in which the

-axis points toward the north

geomagnetic pole, the

-axis is in the plane containing the

-axis and the Sun, and the

-axis is positive

toward dusk and completes a right-handed system (Laundal & Richmond,2017). Freezing the position of

the Sun reduces the temporal variations of the ionospheric delays and facilitates the modeling. The model

coefficients are given by

, ()nm k t

and

, ()nm k t

, where

( ) {1, , }E kt K

represents the time interval that contains ,

that is, the observation span is partitioned into

intervals with constant coefficients.

The ionospheric delay

()

E vt

[m]E

is linked to the

VTEC ( )

[]

electrons/m2

as:

2

VTEC ( )

( ) 40.3 .

(6)

TEC values are usually given in TEC units (

1 1016 2

TECU electrons/m



). The equivalent expression holds

for

()

E it

and the STEC. When multiple systems are used to estimate the ionospheric model coefficients in

Equation5, it is important that all parameters

()

E vt

refer to the same reference frequency or are expressed

directly in TECU.

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 result-

ing from the constant thin layer height, the use of two and four layers has been proposed in Hernández-Pa-

jares etal.(1997,2020).

2.3. Unbiased TEC Estimation

Estimating the unbiased STEC or, equivalently, the unbiased ionospheric slant delays

()

E it

becomes possible

when including the ionospheric model from Section2.2 into the observation model from Section2.1. We

first focus on the model using arc-specific code biases

without the receiver-satellite decomposition.

The estimable biased ionospheric delay parameters

()

E it

for each arc and time instance in Table1 are re-

placed by the

common sets of ionospheric model coefficients



,nm

and



,nm

from Equation5 and the

GF code biases

,GF

for each arc. The so defined set of parameters can be estimated unbiasedly from the

network's observations using least squares.

Let

()

Etb

be a vector containing the known values of the mapping function times the basis functions from

Equations4 and5, and let

()kt

be a vector of the same size with the associated coefficients

, ()nm k t

and

, ()nm k t

so that

T

()

() ()

r r kt

E it tbc

. With the estimated coefficient vector

()

ˆkt

, the unbiased solution of the slant iono-

spheric delay

()

E it

is given by:

T

()

() () .

ˆˆ

r r kt

it tbc

(7)

By adapting the values of the basis functions through varying the latitude and longitude arguments, i.e., by

assuming artificial IPPs covering the entire sphere,

gridded global maps of the VTEC are generated from

the

sets of coefficients



Alternatively, we can first determine the estimates





()

of the biased ionospheric delays. These ionospheric

observables are then used in a second step to estimate the model coefficients



and the GF code biases

,GF

This solution is mathematically identical to the direct solution above, given that





()

are weighted with

their inverse covariance matrix. Since all parameters in Table1, including

s

, are link-specific, they can be

determined on an arc-by-arc basis, assuming no correlations between the observations of different links.

Considering daily GIMs that are computed from a network with hundreds of receivers, for which range pa-

rameters



()

have to be set up for each link and time instance and ambiguity parameters

,

for each link

and frequency, this is a suitable method to avoid the solution of a large-scale system of equations.

Radio Science

BRACK ET AL.

10.1029/2021RS007337

5 of 14

Let the continuous arc between receiver

and satellite

contain the

observation epochs

 1, ,Et T

. The

variance of the observations is assumed as:

 

 

 

 

D () () and D () () ,

ss s s

rj r p rj r p

p t wt t wt



(8)

with

()

E wt

exponential elevation dependent noise amplification factors from Euler and Goad(1991),



the

zenith-referenced variance of the code observations, and



a factor that accounts for the higher precision of

the phase observations and is chosen as

0.0001E



. Correlations between observations are assumed absent.

Through an invertible transformation the

code observations

,()

E pt

can equivalently be expressed by the

time-averaged component:





  1

( ) () ()

s ss s

rj r r rj

p t w wt p t

(9)

with





 

,1 1

1()

s Ts

r tr

E w wt

, and the

1ET

time-differenced components:



, ,,

(1 ) ( ) (1), 2, , ,

s ss

rj rj rj

p t pt p t T

(10)

see Khodabandeh and Teunissen(2016). The same transformation is applied to the phase observations



,()

, the range parameters



()

, and the estimable ionospheric delays

()

E it

With



 



,1 ,

() (), , ()

ss s

r r rF

E t pt p t

and









,1 ,

() (),…, ()

ss s

r r rF

Et t tφ

, the full-rank time-averaged and time-dif-

ferenced observation equations according to Table1 are given by:





  

    



  



    





  

    

    

  





()

E () ()

()

rF r

t it

p e μ Ed

eμ I

φa

(11)

and





  





  





  

  





(1 )

E (1 ) (1 ), 2, , ,

(1 )

t it t T

pe μ

eμ

(12)

with

-vector of ones,





 



1,,

 



,3 ,

[ ,, ]

ss s

r r rF

E ddd

, and

 



,1 ,

[ ,, ]

ss s

r r rF

E aaa

. Matrix

fol-

lows from removing the first two columns of the

-dimensional identity matrix

. The time-averaged and

time-differenced observations are uncorrelated and have mutually exclusive sets of parameters, so that both

systems can be solved separately.

The solution of the time-averaged ionospheric delay

()

E it

and its variance are given by:





it p t p t

() () ()







(13)

D



it w

()















(14)

while the solution of the time-differenced ionospheric delays



()

and their variances and covariances are

given by:





it Ftt t

() () () ( )[ ] ()

112 11 1 1

 





  











st()1





(15)

D



it Fw wt

() () () ()111

















(16)

D





i ti k F

w tk

( ), ( ) () ( ), ,11 11















(17)

with



 



22 2

(1 ) 4E







  1

, and





 





. While the precision of the time-aver-

aged component





()

is driven by the code data, the time-differenced ionospheric delays





()

are on the level

of the phase data. Derivations are given in the Appendix and in Khodabandeh and Teunissen(2016,2017).

The time-averaged and time-differenced estimable ionospheric delays are expressed and parameterized as:

Radio Science

BRACK ET AL.

10.1029/2021RS007337

6 of 14





 



,GF ( ) ,GF

( ) ( ) () ()

s s s ss s s

r r r r r r kt r

i t it d w wt t d

(18)

  

TT

( ) (1)

(1 ) (1 ) ( ) (1) .

ss s s

r r r kt r k

i t it tbc bc

(19)

Since only the time-averaged ionospheric delay

()

E it

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



and can

simply be ignored (Teunissen,2000), thereby also eliminating the bias parameter

,GF

. The model coeffi-

cients



are therefore estimated only from the time-differenced ionospheric observables





()

describing

the unbiased time-differenced ionospheric slant delays

(1 )

E it

, see also Khodabandeh and Teunissen(2017).

The ionospheric observables





()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. Table2. The ionospheric

model can again be directly estimated by replacing the biased ionospheric delay parameters

()

E it

by the

sets of model coefficients



and the GF code biases:



d ddr

d dd

, ,,

GF GF GF



 

(20)

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





()

can be obtained as described above. If more than two frequencies are used, the parame-

Figure 1. Numbers of stations (left) and satellites (right) used to generate the German Research Centre for

Geosciences (GFZ) global vertical total electron content maps.

2000 2005 2010 2015 2020

Time

# of satellites

GPS

GLONASS

Galileo

2000 2005 2010 2015 2020

Time

100

200

300

# of stations

Figure 2. Two exemplary global ionosphere maps showing the vertical total electron content in TECU for moderate ionospheric activity (left) and the St.

Patrick's Day geomagnetic storm (right), both at 8:00 UTC. The dashed lines mark the geomagnetic equator.

GFZ GIM on April 30, 2021

-180-90 090 180

Longitude [deg]

-45

Latitude [deg]

GFZ GIM on March 17, 2015

-180 -90090 180

Longitude [deg]

-45

Latitude [deg]

100

Radio Science

BRACK ET AL.

10.1029/2021RS007337

7 of 14

ters



,rj

and



are not link-specific, so that an arc-by-arc solution is no

longer possible. One can equivalently first determine the parameters of

Table1 and then jointly estimate the model coefficients



, the GF biases



,GFr

and



(Equation20), and the biases



,rj

and



(Table2) from the

resulting





()

and



,

Different from Equation18, the estimable time-averaged ionospheric de-

lay is now parameterized as:





 



1T

( ) ,GF GF

( ) () () ,

s ss s s

r r r r kt r

i t w wt t d d

(21)

and can no longer be ignored, except for the single-receiver case. We

therefore have two uncorrelated sets of observables, namely the time-av-

eraged block with





()

and possibly



s

, and the time-differenced block

with





()1

. The solution and variance of



s

and its covariance with





()

available through EquationA1. The normal equations of both blocks and

all arcs are added. It is shown later that the impact of the receiver-satellite

bias decomposition on the ionospheric model is negligible for GIMs. It

can, however, be used to obtain the GF biases or DCBs as a byproduct.

3. 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 lati-

tude 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 equa-

tion 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än-

nel 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 pro-

vided for GPS to be compatible with the IGS solution. For GLONASS, the

receiver-satellite decomposition of the DCBs is not valid due to the pres-

ence 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 de-

composition 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 cor-

recting the ionospheric delays of

s

in Tables1 and2 with the GIMs (Ke-

shin,2012; Montenbruck etal.,2014), and a complete set of DCBs can be

generated.

Figure 3. Time series of the German Research Centre for Geosciences

global navigation satellite systems-derived global mean vertical total

electron content values, of the F10.7 solar flux index in solar flux units

(SFU), and of the Ap geomagnetic index.

2000 2005 2010 2015 2020

Time

Global mean VTEC [TECU]

100

150

200

250

300

350

Solar flux [SFU] / Ap index

GNSS mean VTEC

Solar flux F10.7

Ap index

Figure 4. Time series of the square root of the estimated variance factor

of the global navigation satellite systems observables.

2000 2005 2010 201

020

Time

Sqrt. of variance factor

Radio Science

BRACK ET AL.

10.1029/2021RS007337

8 of 14

Two exemplary GIMs, one for April 30, 2021, a day with moderate ion-

ospheric activity, and the St. Patrick's Day geomagnetic storm on March

17, 2015, are shown in Figure2.

Global mean VTEC values are obtained by a weighted mean of all grid

values, where the weighting coefficient is proportional to the cosine of

the latitude to account for the surface area represented by each grid point.

Annual and seasonal variations of the GNSS derived daily global mean

VTEC, as well as the signature of the solar-cycle, are shown in 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

(Matzka, Stolle, etal.,2021) generally shows larger values for larger mean

VTEC values, and its maximum values coincide with peaks of the mean

VTEC and F10.7.

Due to modeling errors, the precision of the VTEC that is derived from the law of error propagation from the

estimated model coefficients



does generally not represent its accuracy. For the RMS maps that are provid-

ed along with the VTEC maps, a variance factor is estimated from the residuals relative to a unit standard

deviation of



1m

in zenith direction for the code observations and





1cm



for phase observations,

cf. Equation8. With a perfect model, the estimated variance factor should be below one. In the time series

in Figure4, we see values of

5 40E

for the square root of the variance factor, meaning that the phase re-

siduals 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 decompo-

sition, 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 TECUE

and the RMS VTEC difference is around

0.006TECUE

shown in Figure5. Both values are clearly below the resolution of

0.1TECUE

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.

Figure 5. Time series of daily absolute biases and RMS differences

between the German Research Centre for Geosciences global ionosphere

maps with daily and arc-specific code biases during the year 2010.

2010 2010.2 2010.42010.6 2010.

82011

Time

10-5

VTEC diff. [TECU]

RMS

Abs. bias

Figure 6. Smoothed time series of daily biases (left) and standard deviations (right) between the final global

ionosphere maps of the seven IAACs+GFZ and the final IGS solution.

2005 2010 2015 2020

Time

-5

VTEC bias [TECU]

VTEC bias wrt IGSG

2005 2010 2015 2020

Time

VTEC std [TECU]

VTEC std wrt IGSG

CAS

CODE

ESA

GFZ

JPL

NRCan

UPC

WHU

Radio Science

BRACK ET AL.

10.1029/2021RS007337

9 of 14

4. Comparison With the IGS GIMs

Operational GIMs are provided by several analysis centers based on dif-

ferent processing approaches. Within the IGS, the ionosphere working

group was established in 1998. Currently, seven IGS Ionosphere Associ-

ated Analysis Centers (IAACs) provide rapid and final GIMs, which are

combined to the rapid and final IGS solutions. We compare the final GFZ

solution and the final GIMs of the IAACs to the final IGS solution (IGSG,

Hernández-Pajares etal.,2009). The considered solutions are:

1. Chinese Academy of Sciences (CAS) with CASG (Li etal.,2015).

2. Center for Orbit Determination in Europe (CODE) with CODG

(Schaer etal.,1996).

3. European Space Agency/European Space Operations Center (ESA/

ESOC) with ESAG (Feltens,2007).

4. Jet Propulsion Laboratory (JPL) with JPLG (Mannucci etal.,1998).

5. Natural Resources Canada (NRCan) with EMRG (Ghoddousi-Fard,2014).

6. Universitat Politecnica de Catalunya (UPC-IonSAT) with UPCG (Hernández-Pajares etal.,1999).

7. Wuhan University (WHU) with WHUG (H. Zhang etal.,2013).

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 solu-

tions 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 de-

viations 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. Fig-

ure3, and a few TECU during high ionospheric activity. The GFZ solution is consistent with the solutions

of the IAACs and the IGS combination.

5. Validation With Jason-3 Altimetry VTEC

The VTEC over the oceans can directly be observed from dual-frequen-

cy altimeters (5.3 and

13.6GHzE

for the Jason satellites). Although the al-

timetry VTEC observations potentially suffer from instrumental biases

(Azpilicueta & Brunini,2009), and do not capture the TEC of the upper

ionosphere and plasmasphere above the orbital height (

1.336kmE

for Ja-

son-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 Fig-

ure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.

Figure 7. Jason-2 vertical total electron content in TECU during the day

of the St. Patrick's Day geomagnetic storm.

-180 -90090 180

Longitude [deg]

-45

Latitude [deg]

Jason-2 VTEC on March 17, 2015

100

Figure 8. Daily numbers of Jason-3 vertical total electron content

observations used for validation of the global navigation satellite systems

global ionosphere maps.

2017 20182019

2020

Time

# of Jason-3 VTEC obs.

Radio Science

BRACK ET AL.

10.1029/2021RS007337

10 of 14

For every altimetry VTEC observation, the associated GNSS VTEC is computed through a temporal inter-

polation 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 ob-

served 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 Fig-

ure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 con-

firm 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.

6. 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 coef-

ficients 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 two-

step procedure was formulated. Biased ionospheric delays and potential-

ly code biases for the third frequency and beyond are determined on an

Figure 9. Smoothed time series of daily biases (left) and standard deviations (right) between the final global navigation

satellite systems global ionosphere maps and the Jason-3 vertical total electron content data.

2017 2018 2019 2020

Time

VTEC bias [TECU]

VTEC bias wrt Jason-3

2017 20182019

2020

Time

2.5

3.5

4.5

VTEC std [TECU]

VTEC std wrt Jason-3

CAS

CODE

ESA

GFZ

JPL

NRCan

UPC

WHU

IGS

Figure 10. Distribution of the differences between the global navigation

satellite systems global ionosphere maps and Jason-3 vertical total electron

content data.

-5 05

VTEC GNSS - Jason-3 [TECU]

0.05

0.1

0.15

0.2

Distribution

CAS

CODE

ESA

GFZ

JPL

NRCan

UPC

WHU

IGS

Radio Science

BRACK ET AL.

10.1029/2021RS007337

11 of 14

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 demon-

strated, 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 Equation11, the

ambiguity parameters

s

only appear in the

time-averaged phase observations

()

Etφ

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

unknowns and

equations is solved via inversion as:















it t

()

() , , ()

e Ep























IIF

















st( ),

(A1)

with



 



IF 2 1

[ , ,0, ,0]

Eμ



 



[ 1,1,0, ,0]

, and

 

T TT

IF GF

()

EE E I eμ μμ

, cf.

Khodabandeh and Teunissen(2016). The variances of the estimates follow from applying the error propa-

gation law to Equation8.

Bias Std RMS

CAS 1.55 3.04 3.41

CODE 1.66 3.01 3.44

ESA 1.51 3.35 3.67

GFZ 1.22 3.02 3.25

JPL 3.83 3.10 4.92

NRCan 1.97 3.15 3.71

UPC 2.48 2.81 3.74

WHU 1.88 3.17 3.68

IGS 2.56 2.99 3.93

Note. The values of the GFZ solution are marked in bold.

Table 3

Biases, Standard Deviations, and RMS Differences Between Nine GNSS GIMs and the Jason-3 VTEC Observations; All

Values are in

[]TECU

Radio Science

BRACK ET AL.

10.1029/2021RS007337

12 of 14

Appendix B: Solution of the Time-Differenced System

With







(1 ) (1 ) , (1 )

s ss

r rr

E t tt

y pφ

, stacking the time-differenced observations for

 2, ,Et T

in Equation12

yields:

()





























 

















































()







































 





()













(B1)

The covariance matrix of

,td

follows from Equation8 as:





 





,td

r pr

Q DW D



(B2)

with the differencing operator





 



ED eI

and

-dimensional square matrix with the

values

()

E wt

on its diagonal. We further have:

























 

T ,1 T 1

,td 2

T ,1 T 1

,td 2

T ,1 T 1 2 2

,td 2

()

( ) ( ),

AQ A DW D

AQ B DW D

BQ B DW D



(B3)

where the last equality follows from





 

T 22

()EFμμ

The least squares solutions



(1 )

E it

in Equation15 are the rows of:





i BQ B BQ y

td T

td td





111

(B4)

where

 

  T ,1 1 T ,1

,td ,td

()

EB B AAQ A AQ B

follows from EquationB3 as





 



EB BA



. Using this expres-

sion and EquationB3, the inverse covariance matrix of



(1 )

E it

is derived as:

BQ B DW D

,()

()

(









22 2











 

))

(B5)

cf. Equations16 and17. Similar to EquationB3 we have:



 



 







 





T ,1 2 T 1 T T

,td

T ,1 2 T 1 T T

,td

(),

( ) ,.

r p r FF

r pr

AQ DWD e e

BQ DW D μ μ



(B6)

Using EquationsB5 andB6, the

( 1)Et

th row of EquationB4 reads:





it Ft

() () ()111

1111



















 















TT T

ep eeFr

































(),

(B7)

from which Equation15 is obtained by rearranging the terms.

Data Availability Statement

The GFZ GIMs are available at ftp://isdcftp.gfz-potsdam.de/gnss/products/iono (Brack etal.,2021). The

GNSS observation data are provided by the IGS (Johnston etal.,2017) at https://www.igs.org/data. The sat-

ellite orbits are taken from GFZ's contribution to the third IGS reprocessing campaign (Männel etal.,2021),

Radio Science

BRACK ET AL.

10.1029/2021RS007337

13 of 14

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.spaceweath-

er.gc.ca, and the Ap index by GFZ Potsdam (Matzka, Bronkalla, etal.,2021). The altimetry VTEC data are

produced by DGFI-TUM (Dettmering etal.,2011) and distributed at https://openadb.dgfi.tum.de.

References

Azpilicueta, F., & Brunini, C. (2009). Analysis of the bias between TOPEX and GPS vTEC determinations. Journal of Geodesy, 83(2),

121–127. https://doi.org/10.1007/s00190-008-0244-7

Azpilicueta, F., Brunini, C., & Radicella, S. M. (2006). Global ionospheric maps from GPS observations using modip latitude. Advances in

Space Research, 38(11), 2324–2331. https://doi.org/10.1016/j.asr.2005.07.069

Baarda, W. (1973). S-transformations and criterion matrices. Netherlands Geodetic Commission, Publications on Geodesy.

Brack, A., Männel, B., Bradke, M., Brandt, A., & Nischan, T. (2021). GFZ global ionosphere maps. GFZ Data Services. https://doi.

org/10.5880/GFZ.1.1.2021.006

Bradke, M. (2020). SEMISYS—Sensor Meta Information System. GFZ Data Services. https://doi.org/10.5880/GFZ.1.1.2020.005

Brunini, C., & Azpilicueta, F. (2010). GPS slant total electron content accuracy using the single-layer model under different geomagnetic

regions and ionospheric conditions. Journal of Geodesy, 84(5), 293–304. https://doi.org/10.1007/s00190-010-0367-5

Brunini, C., & Azpilicueta, F. J. (2009). Accuracy assessment of the GPS-based slant total electron content. Journal of Geodesy, 83(8),

773–785. https://doi.org/10.1007/s00190-008-0296-8

Ciraolo, L., Azpilicueta, F., Brunini, C., Meza, A., & Radicella, S. M. (2007). Calibration errors on experimental slant total electron content

(TEC) determined with GPS. Journal of Geodesy, 81(2), 111–120. https://doi.org/10.1007/s00190-006-0093-1

Davies, K., & Hartmann, G. K. (1997). Studying the ionosphere with the global positioning system. Radio Science, 32(4), 1695–1703. https://

doi.org/10.1029/97RS00451

Dettmering, D., Schmidt, M., Heinkelmann, R., & Seitz, M. (2011). Combination of different space-geodetic observations for regional

ionosphere modeling. Journal of Geodesy, 85(12), 989–998. https://doi.org/10.1007/s00190-010-0423-1

Euler, H. J., & Goad, C. C. (1991). On optimal filtering of GPS dual frequency observations without using orbit information. Bulletin Géodé-

sique, 65(2), 130–143. https://doi.org/10.1007/BF00806368

Feltens, J. (2007). Development of a new three-dimensional mathematical ionosphere model at European Space Agency/European Space

Operations Centre. Space Weather, 5(12). https://doi.org/10.1029/2006SW000294

Ghoddousi-Fard, R. (2014). GPS ionospheric mapping at Natural Resources Canada. In Proceeding of IGS Workshop 2014.

Goss, A., Schmidt, M., Erdogan, E., Görres, B., & Seitz, F. (2019). High-resolution vertical total electron content maps based on multi-scale

B-spline representations. Annales Geophysicae, 37(4), 699–717. https://doi.org/10.5194/angeo-37-699-2019

Hauschild, A. (2017). Basic observation equations. In P. J. Teunissen, & O. Montenbruck (Eds.), Springer handbook of global navigation

satellite systems (pp. 561–582). Springer. https://doi.org/10.1007/978-3-319-42928-1_19

Hernández-Pajares, M., Juan, J. M., & Sanz, J. (1997). Neural network modeling of the ionospheric electron content at global scale using

GPS data. Radio Science, 32(3), 1081–1089. https://doi.org/10.1029/97RS00431

Hernández-Pajares, M., Juan, J. M., & Sanz, J. (1999). New approaches in global ionospheric determination using ground GPS data. Journal

of Atmospheric and Solar-Terrestrial Physics, 61(16), 1237–1247. https://doi.org/10.1016/S1364-6826(99)00054-1

Hernández-Pajares, M., Juan, J. M., Sanz, J., Aragón-Àngel, À., García-Rigo, A., Salazar, D., & Escudero, M. (2011). The ionosphere:

Effects, GPS modeling and the benefits for space geodetic techniques. Journal of Geodesy, 85(12), 887–907. https://doi.org/10.1007/

s00190-011-0508-5

Hernández-Pajares, M., Juan, J. M., Sanz, J., Orus, R., Garcia-Rigo, A., Feltens, J., etal. (2009). The IGS VTEC maps: A reliable source of

ionospheric information since 1998. Journal of Geodesy, 83(3–4), 263–275. https://doi.org/10.1007/s00190-008-0266-1

Hernández-Pajares, M., Lyu, H., Garcia-Fernandez, M., & Orus-Perez, R. (2020). A new way of improving global ionospheric maps by ion-

ospheric tomography: Consistent combination of multi-GNSS and multi-space geodetic dual-frequency measurements gathered from

vessel-, LEO- and ground-based receivers. Journal of Geodesy, 94(8), 73. https://doi.org/10.1007/s00190-020-01397-1

Jakowski, N., Schlüter, S., & Sardon, E. (1999). Total electron content of the ionosphere during the geomagnetic storm on January 10, 1997.

Journal of Atmospheric and Solar-Terrestrial Physics, 61(3–4), 299–307. https://doi.org/10.1016/S1364-6826(98)00130-8

Johnston, G., Riddell, A., & Hausler, G. (2017). The international GNSS service. In P. J. Teunissen, & O. Montenbruck (Eds.), Springer

handbook of global navigation satellite systems (pp. 967–982). Springer. https://doi.org/10.1007/978-3-319-42928-1_33

Keshin, M. (2012). A new algorithm for single receiver DCB estimation using IGS TEC maps. GPS Solutions, 16(3), 283–292. https://doi.

org/10.1007/s10291-011-0230-z

Khodabandeh, A., & Teunissen, P. J. G. (2015). An analytical study of PPP-RTK corrections: Precision, correlation and user-impact. Journal

of Geodesy, 89(11), 1109–1132. https://doi.org/10.1007/s00190-015-0838-9

Khodabandeh, A., & Teunissen, P. J. G. (2016). Array-aided multifrequency GNSS ionospheric sensing: Estimability and precision analysis.

IEEE Transactions on Geoscience and Remote Sensing, 54(10), 5895–5913. https://doi.org/10.1109/TGRS.2016.2574809

Khodabandeh, A., & Teunissen, P. J. G. (2017). S-system theory applied to array-based GNSS ionospheric sensing. Studia Geophysica et

Geodaetica, 61(3), 429–452. https://doi.org/10.1007/s11200-016-1176-y

Komjathy, A. (1997). Global ionospheric total electron content mapping using the global positioning system (Dissertation). Department of

Geodesy and Geomatics Engineering Technical Report No. 188. University of New Brunswick.

Krypiak-Gregorczyk, A., & Wielgosz, P. (2018). Carrier phase bias estimation of geometry-free linear combination of GNSS signals for

ionospheric TEC modeling. GPS Solutions, 22(2), 45. https://doi.org/10.1007/s10291-018-0711-4

Lanyi, G. E., & Roth, T. (1988). A comparison of mapped and measured total ionospheric electron content using global positioning system

and beacon satellite observations. Radio Science, 23(4), 483–492. https://doi.org/10.1029/RS023i004p00483

Laundal, K. M., & Richmond, A. D. (2017). Magnetic coordinate systems. Space Science Reviews, 206, 27–59. https://doi.org/10.1007/

s11214-016-0275-y

Li, Z., Yuan, Y., Wang, N., Hernández-Pajares, M., & Huo, X. (2015). SHPTS: Towards a new method for generating precise global iono-

spheric TEC maps based on spherical harmonic and generalized trigonometric series functions. Journal of Geodesy, 89(4), 331–345.

https://doi.org/10.1007/s00190-014-0778-9

Appendix B: Solution of the Time-Differenced System

With





(1 ) (1 ) , (1 )

s ss

r rr

E t tty pφ

, stacking the time-differenced observations for

 2, ,Et T

in Equation12

yields:

()





























 

















































()







































 





()









(B1)

The covariance matrix of

,td

follows from Equation8 as:





  





,td ,

r pr

Q DW D I



(B2)

with the differencing operator





 



ED eI

and

-dimensional square matrix with the

values

()

E wt

on its diagonal. We further have:

























 

T ,1 T 1

,td 2

T ,1 T 1

,td 2

T ,1 T 1 2 2

,td 2

()

( ) ( ),

AQ A DW D

AQ B DW D

BQ B DW D



(B3)

where the last equality follows from





 

T 22

()EFμμ

The least squares solutions



(1 )

E it

in Equation15 are the rows of:





i BQ B BQ y

td T

td td





111

(B4)

where

 

  T ,1 1 T ,1

,td ,td

()

EB B AAQ A AQ B

follows from EquationB3 as





  

EB BA



. Using this expres-

sion and EquationB3, the inverse covariance matrix of



(1 )

E it

is derived as:

BQ B DW D

,()

()

(









22 2











 

))

(B5)

cf. Equations16 and17. Similar to EquationB3 we have:



 



  





 





T ,1 2 T 1 T T

,td

T ,1 2 T 1 T T

,td

(),

( ) ,.

r p r FF

r pr

AQ DWD e e

BQ DW D μ μ



(B6)

Using EquationsB5 andB6, the

( 1)Et

th row of EquationB4 reads:





it Ft

() () ()111

1111



















 







TT T

ep eeFr

































(),1

(B7)

from which Equation15 is obtained by rearranging the terms.

Data Availability Statement

The GFZ GIMs are available at ftp://isdcftp.gfz-potsdam.de/gnss/products/iono (Brack etal.,2021). The

GNSS observation data are provided by the IGS (Johnston etal.,2017) at https://www.igs.org/data. The sat-

ellite orbits are taken from GFZ's contribution to the third IGS reprocessing campaign (Männel etal.,2021),

Acknowledgments

Open access funding enabled and

organized by Projekt DEAL.

Männel, B., Brandt, A., Bradke, M., Sakic, P., Brack, A., & Nischan, T. (2020). Status of IGS reprocessing activities at GFZ. In IAG symposia.

https://doi.org/10.1007/1345_2020_98

Männel, B., Brandt, A., Bradke, M., Sakic, P., Brack, A., & Nischan, T. (2021). GFZ repro3 product series for the International GNSS Service

(IGS). GFZ Data Services. https://doi.org/10.5880/GFZ.1.1.2021.001

Mannucci, A. J., Wilson, B. D., Yuan, D. N., Ho, C. H., Lindqwister, U. J., & Runge, T. F. (1998). A global mapping technique for GPS-derived

ionospheric total electron content measurements. Radio Science, 33(3), 565–582. https://doi.org/10.1029/97RS02707

Matzka, J., Bronkalla, O., Tornow, K., Elger, K., & Stolle, C. (2021). Geomagnetic Kp index. GFZ Data Services. https://doi.org/10.5880/

Kp.0001

Matzka, J., Stolle, C., Yamazaki, Y., Bronkalla, O., & Morschhauser, A. (2021). The geomagnetic Kp index and derived indices of geomag-

netic activity. Space Weather, 19(5), e2020SW002641. https://doi.org/10.1029/2020SW002641

Montenbruck, O., & Hauschild, A. (2013). Code biases in multi-GNSS point positioning. Proceeding of ION-ITM 2013, 616–628.

Montenbruck, O., Hauschild, A., & Steigenberger, P. (2014). Differential code bias estimation using multi-GNSS observations and global

ionosphere maps. Navigation, 61(3), 191–201. https://doi.org/10.1002/navi.64

Roma-Dollase, D., Hernández-Pajares, M., Krankowski, A., Kotulak, K., Ghoddousi-Fard, R., Yuan, Y., etal. (2018). Consistency of sev-

en different GNSS global ionospheric mapping techniques during one solar cycle. Journal of Geodesy, 92(6), 691–706. https://doi.

org/10.1007/s00190-017-1088-9

Sardon, E., Rius, A., & Zarraoa, N. (1994). Estimation of the transmitter and receiver differential biases and the ionospheric total electron

content from global positioning system observations. Radio Science, 29(03), 577–586. https://doi.org/10.1029/94RS00449

Schaer, S. (1999). Mapping and predicting the Earth's ionosphere using the global positioning system (Dissertation). Astronomical Institute,

University of Bern.

Schaer, S., Beutler, G., Rothacher, M., & Springer, T. A. (1996). Daily global ionosphere maps based on GPS carrier phase data routinely

produced by the CODE analysis center. In Proceeding of IGS AC Workshop 1996.

Schaer, S., Gurtner, W. & Feltens, J. (1998). IONEX: The ionosphere map exchange format version 1. In Proceeding of IGS AC Workshop

1998.

Tapping, K. F. (2013). The 10.7 cm solar radio flux (F10.7). Space Weather, 11(7), 394–406. https://doi.org/10.1002/swe.20064

Teunissen, P. J. G. (1985). Generalized inverses, adjustment, the datum problem and S-transformations. In E. W. Grafarend, & F. Sansò

(Eds.), Optimization and design of geodetic networks (pp. 11–55). Springer. https://doi.org/10.1007/978-3-642-70659-2_3

Teunissen, P. J. G. (2000). Adjustment theory, an introduction. Delft Academic Press, Series on Mathematical Geodesy and Positioning.

Wang, N., Yuan, Y., Li, Z., Montenbruck, O., & Tan, B. (2016). Determination of differential code biases with multi-GNSS observations.

Journal of Geodesy, 90(3), 209–228. https://doi.org/10.1007/s00190-015-0867-4

Wickert, J., Dick, G., Schmidt, T., Asgarimehr, M., Antonoglou, N., Arras, C., etal. (2020). GNSS remote sensing at GFZ: Overview and recent

results. zfv—Zeitschrift für Geodäsie, Geoinformation und Landmanagement, 145(5), 266–278. https://doi.org/10.12902/zfv-0320-2020

Zhang, B. (2016). Three methods to retrieve slant total electron content measurements from ground-based GPS receivers and performance

assessment. Radio Science, 51(7), 972–988. https://doi.org/10.1002/2015RS005916

Zhang, H., Xu, P., Han, W., Ge, M., & Shi, C. (2013). Eliminating negative VTEC in global ionosphere maps using inequality-constrained

least squares. Advances in Space Research, 51(6), 988–1000. https://doi.org/10.1016/j.asr.2012.06.026

Zhang, X., Xie, W., Ren, X., Li, X., Zhang, K., & Jiang, W. (2017). Influence of the GLONASS inter-frequency bias on differential code bias

estimation and ionospheric modeling. GPS Solutions, 21(3), 1355–1367. https://doi.org/10.1007/s10291-017-0618-5

14 of 14

Radio Science 10.1029/2021RS007337

BRACK ET AL. 14 of 14