Single-band VLBI Absolute Astrometry

The impact of the ionosphere dispersiveness on the fringe phase is reciprocal to the frequency in the first approximation. Therefore, the fringe phase in channel i in the presence of the ionosphere becomes

$$\begin{eqnarray}&&{\phi }_{i}=2\pi {\tau }_{p}\,{f}_{0}+{\tau }_{g}\,({f}_{i}-{f}_{0})-\displaystyle \frac{\alpha }{{f}_{i}},\end{eqnarray} \tag{ 1 }$$

where τ_p and τ_g are the phase and group delays, respectively, f_i is the frequency of the ith spectral channel, f₀ is the reference frequency, and

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{e}^{2}}{8\,{\pi }^{2}\,c\,{m}_{e}\,{\epsilon }_{o}}\left(\int {N}_{v}\,{{ds}}_{1}-\int {N}_{v}\,{{ds}}_{2}\right),\end{eqnarray} \tag{ 2 }$$

where N_v is the electron density, e is the charge of an electron, m_e is the mass of an electron, _o is the permittivity of free space, and c is the velocity of light in a vacuum. Integration is carried along the line of sight. Having substituted the values of the constants (Klobuchar 1996) and expressing the total electron content (TEC) along the line of sight ∫N_v ds in 1 · 10¹⁶ electrons/m² (so-called TEC units, or TECU), we arrive to α = 1.345 · 10⁹ s/TECU.

The phase and group delays are computed using the fringe phases ϕ_i with weights w_i using least squares. The result can be expressed analytically after some algebra

$$\begin{eqnarray}&&{\tau }_{\mathrm{gi}}={\tau }_{\mathrm{if}}+\displaystyle \frac{\alpha }{{f}_{e}^{2}}\,\mathrm{TEC},\end{eqnarray} \tag{ 3 }$$

where τ_if is the ionosphere-free group delay, TEC is ∫N_v ds expressed in TEC units, and f_e is the effective ionospheric frequency

$$\begin{eqnarray}&&{f}_{e}=\sqrt{\displaystyle \frac{{\sum }_{i}^{n}{w}_{i}\cdot {\sum }_{i}^{n}{w}_{i}{\left({f}_{i}-{f}_{0}\right)}^{2}-{\left({\sum }_{i}^{n}{w}_{i}({f}_{i}-{f}_{0})\right)}^{2}}{{\sum }_{i}^{n}{w}_{i}({f}_{i}-{f}_{0}){\sum }_{i}^{n}\tfrac{{w}_{i}}{{f}_{i}}-{\sum }_{i}^{n}{w}_{i}\cdot {\sum }_{i}^{n}{w}_{i}\tfrac{({f}_{i}-{f}_{0})}{{f}_{i}}}},\end{eqnarray} \tag{ 4 }$$

which depends on the weights of the spectral channels w_i . Typically, the effective ionospheric frequency is within several percents of the central frequency of the observing band.

The best way to mitigate the impact of the ionosphere on the group delay is to observe simultaneously in two or more widely separated frequency bands. Then the following linear combination of two group delays at the upper and lower bands, τ_u and τ_l , respectively, is ionosphere free

$$\begin{eqnarray}&&{\tau }_{{if}}=\displaystyle \frac{{f}_{u}^{2}}{{f}_{u}^{2}-{f}_{l}^{2}}\,{\tau }_{u}-\displaystyle \frac{{f}_{l}^{2}}{{f}_{u}^{2}-{f}_{l}^{2}}\,{\tau }_{l}.\end{eqnarray} \tag{ 5 }$$

Here f_u and f_l are the effective ionospheric frequencies at the upper and lower bands, respectively.

The residual contribution of the ionosphere in dual-band combinations is caused by systematic errors, namely (a) higher-order terms in the expansion of the dispersiveness on frequency (Hawarey et al. 2005), (b) the contribution of frequency-dependent source structure, and (c) the dispersiveness in the signal chain. These contributions affect the group delay at a level of several picoseconds and they are considered insignificant with respect to other systematic errors.

TEC maps (Hernández-Pajares et al. 2009) also known as global ionospheric models (GIMs) derived from the analysis of Global Navigation Satellite System (GNSS) data are used for the reduction of single-band observations. In particular, I used a CODE TEC time series (Schaer 1999) available at ftp.aiub.unibe.ch/CODE since 1995 January 01 with a spatial resolution of 5° × 25. The time resolution varied: it was 24 hr from 1995 January 01 to 1998 March 27; 2 hr from 1998 March 28 to 2014 October 18; and 1 hr after that date. The ionosphere is considered as a thin shell at a constant height H_i of 450 km above the mean Earth's radius. The ionospheric contribution is expressed via TEC as

$$\begin{eqnarray}&&{\tau }_{i}=\displaystyle \frac{\alpha }{{f}_{e}^{2}}\ M(e)\ {\rm{TEC}},\end{eqnarray} \tag{ 6 }$$

where M(e) is the so-called thin-shell ionospheric mapping function

$$\begin{eqnarray}&&M(e)=\displaystyle \frac{1}{\sqrt{1-{\left(\tfrac{{\bar{R}}_{\oplus }}{{R}_{\oplus }+{H}_{i}}\right)}^{2}\,{\cos }^{2}{e}_{\mathrm{gc}}}},\end{eqnarray} \tag{ 7 }$$

where R_⊕ is the mean Earth's radius and e_gc is the geocentric elevation angle with respect to the radius vector between the geocenter and the station. Computation of the TEC value at a given moment of time is reduced to computation of the position of the ionosphere piercing point at a given azimuth and elevation and interpolation of TEC at the piercing point at that latitude, longitude, and time.

TEC maps from GNSS provide a coarse model of the ionosphere. The errors of τ_i computed according to Equation (5) are greater than the residual ionosphere contribution of ionosphere-free linear combinations of group delays. Therefore, dual-band delay observables are preferred when they are available. However, there are two cases when they are not available: (a) dual-band observing sessions with some source detected only in one band and (b) single-band observing sessions. In these two cases we resort to computation of the ionospheric contribution to path delays using GNSS TEC maps and evaluation of the uncertainties of these contributions.

The simplest way to deal with a mixture of dual- and single-band data is to process experiments three times: (1) using dual-band data of those observations that detected a source in both bands, (2) using low-band data, and (3) using upper-band data and applying an ionospheric path delay computed from GNSS TEC maps. However, typically only a fraction, 2%–20%, of observations is detected in only one band; the rest of observations are detected in both bands. Therefore, we can use available dual-band observations at a given observing session to improve the TEC model.

I represent the ionospheric path delay at stations j, k as

$$\begin{eqnarray}\begin{array}{lcl}{\tau }_{i}(t) & = & {b}_{{\rm{j}}}({\rm{t}})-{b}_{k}(t)\\ & & +\,\displaystyle \frac{\alpha }{{f}_{e}^{2}}\left(\left({\mathrm{TEC}}_{{\rm{j}}}({\phi }_{{\rm{j}}},{\lambda }_{{\rm{j}},{\rm{t}}})+{{\rm{a}}}_{{\rm{j}}}(t)\right){\rm{M}}({{\rm{e}}}_{{\rm{j}}})\right.\\ & & -\,\left.\left({\mathrm{TEC}}_{{\rm{k}}}({\phi }_{{\rm{k}}},{\lambda }_{{\rm{k}}},{\rm{t}})+{{\rm{a}}}_{{\rm{k}}}({\rm{t}})\right){\rm{M}}({{\rm{e}}}_{{\rm{k}}}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where ${b}_{{\rm{j}}}({\rm{t}})={\sum }_{i}^{n}{B}_{i}^{0}(t){b}_{{ij}}$ is a delay bias expanded over the B-spline basis of the zeroth degree and ${a}_{{\rm{j}}}({\rm{t}})={\sum }_{i}^{n}{B}_{i}^{3}(t){a}_{{ij}}$ is the TEC bias expanded over the B-spline basis of the third degree. ϕ, λ are coordinates of the ionosphere piercing point, which depend on the positions of the observing stations as well as the azimuths and elevations of the observed sources.

The clock bias occurs due to path delay in VLBI hardware, which is different in different bands. This bias is constant for most of the experiments, however occasionally breaks may happen at some stations. The epochs of these breaks coincide with the epochs of breaks in the clock function. Expansion over the B-spline basis of the zeroth degree accounts for these breaks. (The B-spline of the zeroth degree is 1 within the range of knots [i, i+1] and 0 otherwise.)

I estimated the parameters a_j for all the stations and b_j for all the stations except the one taken as a reference using all the available dual-band observations of a given experiment using least squares with weights

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{1}{\sqrt{{y}^{2}+\tfrac{{f}_{u}^{4}{\sigma }^{2}({\tau }_{u})+{f}_{l}^{4}{\sigma }^{2}({\tau }_{l})}{{\left({f}_{u}^{2}-{f}_{l}\right)}^{2}}}},\end{eqnarray} \tag{ 9 }$$

where σ(τ_u ) and σ(τ_u ) are group delay uncertainties and y is the error floor, 12 ps, introduced to avoid observations with very high signal-to-noise ratios dominating the solution.

The time span of the B-spline knot sequence for the TEC bias in my solutions was 15 minutes. I applied constraints on the values of the B-spline coefficients and on the first and the second derivatives with the reciprocal weights 5 × 10⁻¹⁰ s, 4 × 10⁻¹⁴, and 2 × 10⁻¹⁸ s⁻¹, respectively. These constraints were introduced to ensure the continuity of biases and to prevent a singularity in rare cases when too few available observations at a given station could be used for bias estimation in a given spline segment. Figure 1 illustrates the estimates of the ionospheric bias.

Zoom In Zoom Out Reset image size

The resulting total electron content model TEC_j(t) + a_j(t) is more precise than the a priori TEC_j(t) model taken from GNSS maps because it uses additional information. Using estimates of the a_j and b_j spline coefficients, I compute τ_i (t) and its uncertainty according to the law of error propagation using the full variance-covariance matrix of the estimated spline coefficients. In order to evaluate the realism of these errors, I processed the trial data set and computed τ_i (t) using the estimates of the clock and TEC biases and compared them with the ionospheric contribution derived from dual-band observations. I removed clock biases from the VLBI dual-band ionospheric contributions τ_vi, formed the difference τ_i − τ_vi, and then divided them by σ(τ_i ).

The uncertainty σ(τ_i ) was derived from the variance-covariance matrix of a_j and b_j estimates. I generated a normalized histogram from a data set of 4,343,782 differences and computed the first two moments of the empirical distribution shown in Figure 2. The fitting parameters of the first and second moments of the distribution are 0.003 and 0.889, respectively. Two factors cause a deviation of the second moment from 1.0 in the opposite way: (a) TEC variations not accounted for by the parametric model and (b) statistical dependence of the estimates of a_j , b_j , and the VLBI path delay. After scaling the variance-covariance matrix by the square of 0.889, the distribution of the normalized residuals becomes close to Gaussian. The closeness of the empirical distribution to the normal distribution provides us with confidence that the extra noise introduced by the mismodeled ionospheric path delay after applying clock and TEC biases is properly considered.

Zoom In Zoom Out Reset image size

The closeness of the distribution of the normalized differences to a Gaussian with σ = 1 is encouraging, but it does not guarantee that the residual errors of the sum of the TEC from GNSS and the TEC bias adjustments cause no systematic errors in the estimates of source positions. To characterize the impact of the residual errors of the ionospheric contribution on source position, I ran solution XIA that had the following differences with respect to the reference solution: (1) it used X-band group delays; (2) data reduction for the ionosphere accounted for both a priori TEC from GNSS maps and ionosphere bias adjustment (expression in the parentheses in Equation (8)); and (3) the errors of the ionospheric biases $\sigma ={a}_{\mathrm{ij}}\,\alpha /{f}_{e}^{2}$ were added in quadrature to the reciprocal weights of the observables. It should be stressed that the parameterization, editing, and delay uncertainties were exactly the same as in the reference dual-band solution. In all trial solutions I only varied the input observables, changed the ionosphere-specific reduction model, and added noise to the delay uncertainties. Since I analyzed only the differences in the solutions, the contributions of other model deficiencies, such as the path delay in the neutral atmosphere, were canceled when I formed the differences.

I ran also solution SIA that differed from XIA by using S-band group delays and the S-band effective ionospheric frequency f_e . For control, I ran solutions XIN and SIN that used ionosphere-free combinations of the group delays and the same weights as used in the XIA and SIA solutions, respectively, and zero-mean Gaussian noise with $\sigma ={a}_{\mathrm{ij}}\,\alpha /{f}_{e}^{2}$ was added to each observable. The differences in the XIA and SIA solutions with respect to the reference solution provide us with a measure of the impact of the residual ionospheric errors on source position. The differences in the XIN and SIN solutions characterize the impact of the residual ionospheric errors if they were Gaussian and totally uncorrelated.

Figure 3 shows the distribution of the differences in the declinations of source position estimates from the XIA and XIN solutions. There is no noticeable deviation from the Gaussian shape. Table 1 lists first two moments of the distribution of the position differences in both right ascension and declination. The second moment from the SIN solution is close to 1.0, while the second moment from the XIN solution is 0.56. The errors of the ionospheric contribution in the S band dominate the error budget. These errors are 14 times less for the X band and are only a fraction of the overall group delay errors. The mean biases in right ascension and declination are negligible. The second moment of the position estimates from the XIA and SIA solutions is 15% higher than the moments from the positions from the XIN and SIN solutions. This increase occurs due to nonrandomness of the residual ionospheric errors and can be viewed as a measure of the unaccounted systematic errors in the source positions due to the ionosphere. Analysis of these trial solutions demonstrates that we are able to predict the impact of the ionospheric errors on source position with an accuracy of 15% when the TEC bias is estimated using dual-band observations.

Zoom In Zoom Out Reset image size

TEC biases cannot be computed when an entire session is observed only in one band. Therefore, we have to resort to deriving a regression model to provide estimates of these errors. In the past, Petrov et al. (2019) derived a regression against the so-called global TEC: the integral of the TEC over the entire Earth, following the ideas of Afraimovich et al. (2008). Later, Krásná & Petrov (2021) derived a regression against a function of the rms of the total ionospheric path delay from GNSS TEC. In this study I use the second approach with slight modifications. Following the general results of turbulence theory (see Tatarskii 1971), we can expect that fluctuations at scale x are related to fluctuations at scale y via a power law.

I processed the same data set of data from 263 24 hr VLBI experiments that were used in the previous section and computed the residual ionospheric path delay for each observation as

$$\begin{eqnarray}&&{\tau }_{r}=({\tau }_{\mathrm{gi}}-{\tau }_{\mathrm{vi}}-({b}_{j}-{b}_{k}))\,\tilde{M},\end{eqnarray} \tag{ 10 }$$

where τ_gi is the vertical ionospheric path delay from GNSS TEC maps, τ_vi is the vertical ionospheric path delay from VLBI, b_j − b_k are contributions of the clock bias, and $\tilde{M}$ is the averaged ionosphere mapping function between stations 1 and 2 of a given baseline: $\tilde{M}=(M({e}_{1})+M({e}_{2}))/2.0$. The clock biases are routinely adjusted during analysis of VLBI observations and therefore their contribution to VLBI results, such as source positions, is entirely eliminated. Subtracting them in Equation (10), I eliminate their impact on the statistics as well. I used only 24 hr VLBI experiments for deriving the statistics because the ionospheric path delay strongly depends on solar time, especially at low latitudes, and the statistics derived for shorter time intervals may not be representative.

Figure 4 shows the dependence of the rms of the residuals τ_r on the rms of the total ionospheric path delay from the GNSS TEC maps τ_gi. Each point on the plot corresponds to the rms for a given baseline and a given observing session. I confirm the early result of Krásná & Petrov (2021) but here I used a much larger data set. The result reported in Krásná & Petrov (2021) was slightly affected by an error in computing the ionospheric group delays for a case when some data are flagged because of radio interference. This error has been fixed and the affected experiments have been reprocessed from scratch. This dependence can be coarsely described as the square root of the rms of the total ionospheric path delay. For a better approximation I sought a regression in the form of an expansion over B-splines of the third degree. The spline coefficients computed using least squares are listed in Table 2.

Zoom In Zoom Out Reset image size

Using that regression, I developed the following algorithm for computing the errors of the ionospheric path delay from GNSS TEC maps. First, the coordinates of K points uniformly distributed over the sphere are computed using a random number generator. Then for each baseline and each time epoch, the azimuth and elevation angles of that point are computed at both stations of the baseline, and if the elevations above the horizon are greater than 5° at both stations, that point is selected for further computations. If not, the next point is drawn. Then the total ionospheric path delay τ_i (A₁, e₁, A₂, e₂) is computed using the GNSS TEC maps. It is worth mentioning here that unlike the troposphere path delay, τ_i (A, e) ≠ τ_i (A, π/2) × M(e), since the path delay depends on the positions of the ionosphere piercing points. It is not sufficient to compute the ionospheric path delay in the zenith direction and then map it via M(e): latitude and longitude of the piercing point can be as far as 1000 km from the station. Following this approach, we sample piercing points uniformly distributed within the mutual visibility zone. The process is repeated for 1440 time epochs that cover the time interval of the VLBI experiment under consideration with a step of one minute. Then for each baseline σ(τ_gt) is computed over a time series of 1440 τ_i values. Finally, the estimate of the rms of residual ionospheric path delay derived from the GNSS TEC maps is computed from this regression via the rms of the total ionospheric path delay as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{rr}}=\sum _{k}^{n}{B}_{k}^{3}(\sigma ({\tau }_{\mathrm{gt}}))\,\sqrt{{M}^{2}({e}_{1})+{M}^{2}({e}_{2})}.\end{eqnarray} \tag{ 11 }$$

Baseline-dependent data sets are considered independent for this computation: the mutual visibility at all the stations of the network at a given moment of time is not enforced. For several baselines longer than 96% of Earth's diameter, this algorithm has poor performance for selecting points above 5°. Therefore a minor modification is made for such an extreme case: the elevation angle is fixed to 5°, mutual visibility is not enforced, and azimuths are selected randomly within a range of [0, 2π] independently for both stations.

In order to evaluate the validity of this regression model of the residual ionospheric path delay errors, I computed τ_r from dual-band observations σ_gt following the algorithm described above for 263 24 hr experiments, and then computed histograms of the normalized residuals τ_r /σ_rr. One such histogram is presented in Figure 5. The first two moments of the distribution are −0.083 and 1.214, respectively. Since regression σ_rr was found using least squares, the number of observations with σ(τ_gr) less and greater than σ_rr for given σ(τ_gt) is approximately equal; the thick blue line cuts the cloud of green points in Figure 5 almost by half. However, the variance of the contribution of those points with σ(τ_r ) > σ_rr overweights the contribution of those points with σ(τ_r ) < σ_rr because variance quadratically depends on τ_r . This causes a positive bias. After multiplying σ_rr by 1.214, the distribution of normalized residuals becomes almost Gaussian.

Zoom In Zoom Out Reset image size

One can notice that $\sqrt{{M}^{2}({e}_{1})+{M}^{2}({e}_{2})}$ in Equation (11) is not the same as $\tilde{M}=(M({e}_{1})+M({e}_{2}))/2$ used for computation of the regression. I found that using $\tilde{M}$ instead of $\sqrt{{M}^{2}({e}_{1})+{M}^{2}({e}_{2})}$ decreases the second moment from to 1.214 to 1.196, which is negligible,

The distributions shown in Figures 2 and 5 are computed for the entire data set of 4.3 million path delays and they represent the general population over the interval of 23 yr. Statistics for an individual observing session may differ. In order to evaluate the scatter of the statistics, I computed a time series of second moments of the distribution of normalized residuals of ionospheric path delays and their uncertainties with and without TEC biases adjusted for each observing session separately. I divided the normalized residuals by scaling factors of 0.889 and 1.196, respectively. I computed the distribution of second-moment estimates and show it in Figure 6. The scatter of the second moments is small when TEC biases are adjusted. That means these statistics are robust. When TEC is not adjusted, the scatter is significantly larger, but even in that case 90% of the second-moment estimates deviate from 1.0 by no more than 30%. This provides us with a measure of the uncertainties in the computation of ionospheric path delay errors in individual single-band observing sessions.

Zoom In Zoom Out Reset image size

The error analysis presented in the previous section characterizes our ability to predict the first and second moments of the distribution, but it does not guarantee that the residual errors due to the ionosphere cause no systematic errors in the source positions. I ran trial solutions XIT and SIT that used X- and S-band group delays, respectively, applied ionospheric path delays derived from the GNSS TEC maps, and inflated the reciprocal weights of the observables by adding in quadrature the errors in the ionospheric path delays from the TEC maps.

The analysis of the differences in the source positions from the trial XIT and SIT solutions with respect to the reference dual-band solution revealed no peculiarities in right ascension, but revealed significant systematic errors in declination (see Figure 7). The pattern of systematic errors in declination from the S-band solution is similar but greater by a factor of ${f}_{x}^{2}/{f}_{s}^{2}\approx 14$. We cannot consider a solution with such errors as satisfactory. This was unexpected because the prior work of Sekido et al. (2003), Hobiger et al. (2006), Dettmering et al. (2011), and Motlaghzadeh et al. (2022) claimed a good agreement between the TEC derived from VLBI and GNSS. And indeed, the plot of the ionospheric contributions in the zenith direction from the VLBI after removal of the clock biases against the ionospheric contributions from GNSS (Figure 8) shows no peculiarities. The least-square approximation of this dependence to the straight line is τ_vi = −3.9 ps + 1.06 × τ_gi. Although the residuals of this dependence look random, they still cause systematic errors in source positions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

I made a number of trial solutions. The following leads turned out to be productive: (1) to modify the mapping function and (2) to scale the ionospheric path delays from the TEC. Schaer (1999) suggested the following modification of the ionospheric mapping function

$$\begin{eqnarray}&&M(e)=k\,\displaystyle \frac{1}{\sqrt{1-{\left(\tfrac{{\bar{R}}_{\oplus }}{{R}_{\oplus }+{H}_{i}+{\rm{\Delta }}H}\right)}^{2}\,{\cos }^{2}\alpha {e}_{\mathrm{gc}}}},\end{eqnarray} \tag{ 12 }$$

arguing that by varying parameters α and H_i one can account for a more realistic electron density distribution with height than a thin-shell model. Here ΔH is an increment in the ionosphere height and α is a fudge factor.

I ran 12 trial solutions with mapping functions with parameters (1) ΔH = 0.0 km, α = 1.0; (2) ΔH = 56.7 km, α = 0.9782; and (3) ΔH = 150.0 km, α = 0.9782. Variant 2 corresponds to the so-called JPL-modified single-layer ionospheric mapping function. It was used in a number of papers (Li et al. 2019; Xiang & Gao 2019; Wielgosz et al. 2021; Zhao et al. 2021), however no details of how these values of k, ΔH, and α were derived have been provided. I varied the mapping function scaling factor k, setting it to 0.7, 0.8, 0.9, and 1.0 and computed declination biases with respect to the reference solution. The results are presented in Figures 9–11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The declination bias is reduced when the parameter ionosphere height in the mapping function is increased. I interpret this result as a deficiency of the GNSS TEC maps and I associate the origin of this deficiency with an oversimplification of the mapping function that was used for derivation of the TEC maps from processing of GNSS observations.

Close analysis of these figures reveals that the declination bias has three components: (1) a constant; (2) a linear increase in the declination bias with a decrease of declination; and (3) the feature $\delta (x\,\sin \delta +y\,\cos \delta )$, with the minimum at approximately −12°, where δ is the declination. All three components depend on the parameters of the mapping function, ΔH and α, and on the scaling factor k. Selecting the optimal ΔH, α, k, we can reduce the declination biases. I selected ΔH = 56.7 km, α = 0.9782, and k = 0.85. This choice makes the weighted mean declination bias over all sources 0.013 mas. There is an element of subjectivity in this specific choice, since there exist another combinations of ΔH, α, k that provide a zero mean bias. Wielgosz et al. (2021) showed that the use of the modified mapping function results in a reduction of the mean rms of the difference between GNSS TEC maps and GNSS slant TEC in a range of 2 to 6%. That choice of ΔH and α led me to a selection of the specific scaling factor k.

Although scaling and modification of the ionospheric mapping function results in a substantial reduction of the declination bias, the remaining bias is still worrisome. To mitigate the bias even further, I introduce an ad hoc correction for the ionospheric bias in the data reduction model. I smoothed the biases D(δ) with a Gaussian filter with σ = 8°—these smoothed biases are shown in Figures 9–11—and expanded them over the basis of B-splines of the third degree with 13 equidistant knots in the range of [−45°, 90°]. The expansion coefficients are presented in the Appendix. I added the correction

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tau }{\partial \delta }\,D(\delta )/{f}_{e}^{2}\end{eqnarray} \tag{ 13 }$$

to the observables. Here f_e is the effective ionospheric path delay of a given observation. Figure 12 shows the effect of applying the debias correction. The bias has gone.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

I ran two solutions using X-band only data (XIB) and S-band only data (SIB), and applied both the a priori ionospheric path delays derived from GNSS TEC maps using the modified mapping function with parameters ΔH = 56.7 km, α = 0.9782, and k = 0.85 and the declination bias correction. The reciprocal weights of the observables were adjusted by adding in quadrature the residual ionospheric errors σ_rr modeled according to the regression in Equation (11). The first and the second moments of the normalized differences in sources positions are presented in Table 3. The second moments are less than 1.0 in the X-band solution. This indicates that the ionospheric contribution is not the dominant error source in these solutions.

Unfortunately, it does not look possible to fix the deficiency of the GNSS TEC maps without reprocessing the GNSS observations, which is well beyond the scope of this article, but nevertheless, it is still feasible to mitigate the impact of the imperfection of GNSS TEC maps on source position estimates.

A number of authors compared the ionospheric contribution from GNSS, VLBI, DORIS, radio occultation observations, and dual-band satellite altimeters (Sekido et al. 2003; Hobiger et al. 2006; Dettmering et al. 2011; Hernández-Pajares et al. 2017; Cokrlic et al. 2018; Li et al. 2018, 2019; Xiang & Gao 2019; Wielgosz et al. 2021; Zhao et al. 2021; Motlaghzadeh et al. 2022). A message these publications convey is there is a reasonable agreement between GNSS TEC maps and other techniques and there are no majors problems. Therefore, large systematic errors driven by mismodeling of the ionospheric contribution derived from GNSS maps came as a surprise. Do the results presented in this study contradict those of the prior publications?

Comparison in the past was often performed using very short continuous data sets of VLBI observations, from 5 to 15 days (Hobiger et al. 2005; Dettmering et al. 2011; Etemadfard et al. 2021; Motlaghzadeh et al. 2022). The level of agreement was characterized in terms of additive errors in the GNSS TEC model. Unfortunately, when data from a short time interval are analyzed, a distinction between additive and multiplicative errors becomes blurry. Characterizing the differences in terms of biases and rms of their scatter did not turn out productive in revealing systematic errors. Moreover, an unconscious bias of focusing a study on the assessment of an agreement rather than investigation of disagreements, which were noticed and just briefly reported, diverted attention from investigating the differences in depth.

However, reading between the lines of the published papers, we can find pieces of evidence supporting the findings in this work. The distribution of differences in the VLBI vertical TEC with respect to the TEC derived from GNSS TEC maps presented in Figure 6 in Hobiger et al. (2006) shows a very significant skew. The distribution has a negative mean and a much greater left tail than a right tail. The VLBI vertical TEC from that study appeared less than the TEC from the GNSS TEC maps, in agreement with what I have found. The authors did not investigate that stark deviation of the distribution of over one million differences from the Gaussian shape, only noting that the bias is less than 3 TECU. Considering the errors are normally distributed and additive, and considering that the derivation of the TEC from the VLBI does not introduce serious errors, one can expect to arrive to a Gaussian distribution of the residuals. For instance, the distribution of differences in Figure 5 indeed does not show any measurable deviation from a Gaussian shape. However, if the errors are multiplicative, the residuals will be non-Gaussian and their distribution will be skewed.

Satellite altimetry using Jason satellites (Vaze et al. 2010, and references therein) provides an independent way for assessing the level of disagreement between direct vertical TEC measurements and GNSS TEC maps. Li et al. (2018) showed that comparisons of the differences revealed significant systematic biases that depend on geomagnetic latitude. However, an attempt to characterize these additive biases in terms of the rms was not very productive. Liu et al. (2018) presented the spatial distribution of the differences (their Figures 6–7). That distribution strikingly reminds us of the average distribution of the TEC itself, suggesting the differences are multiplicative. The more recent paper of Dettmering & Schwatke (2022) revealed a highly significant scaling factor between GNSS TEC maps and four altimeter missions. The scaling factors estimated separately from the analysis of these four missions varied from 0.809 to 0.919, which are very close to what has been found in my analysis of astrometric observations of extragalactic radio sources.

Li et al. (2019) performed a comparison of TEC maps with Jason satellite altimetry and with ionospheric radio occultations (see the overview of this technique in Bonafoni et al. 2019). They characterized the differences as the superposition of an additive bias and a multiplicative scale factor. The scale factors defined as TEC${}_{\mathrm{Jason}/\mathrm{GIM}}$ and TEC_COSMIC/TEC${}_{\mathrm{GIM}}$ vary with time and latitude and stay in a range of 0.7–0.9 at night and 0.9–1.0 during the daytime. Since Jason orbits have an altitude of 1350 km and GNSS orbits have an altitude of about 20,000 km, Jason altimetry misses the contribution from the upper layers of the ionosphere and the plasmasphere. Yizengaw et al. (2008) analyzed the TEC between the Jason and GNSS satellites and found that the share of plasmasphere to the TEC is on average 15% and may reach 60% when the TEC is low. However, such a large share of the electron density at altitudes above 1350 km is inconsistent with the thin-shell model at altitudes of 450 or 506.7 km, which assumes no electron density above that height at all.

It is essential to note that the analysis of Jason and COSMIC data provides estimates of the vertical TEC, while the analysis of GNSS observations provides the slant TEC that is converted to the vertical TEC using a mapping function. The dependence of the declination biases on the mapping function found in this study suggests that a simple thin model may not be adequate. Schaer (1999) discussed the dependence of the effective ionosphere height on the solar zenith angle. Xiang & Gao (2019) studied it in more detail. The height of the peak electron density has annual and diurnal variations. The latter variations have an amplitude of about 100 km, being lower at daytime. They showed that the instantaneous mapping function that accounts for the height of the electron content maximum achieves a 8% reduction of mapping errors. It should be mentioned that if the effective height of the ionosphere is changed, the latitude and longitude of the ionosphere piercing point for an observation with a given elevation and azimuth will change as well, which will provide an additional change in the path delay.

These works strengthen argumentation in favor of that the thin-shell model used for the derivation of TEC maps since the mid-1990s is oversimplified. The effective height of the ionosphere changes with time and latitude and therefore, a realistic mapping function should also vary with time and latitude. Omission of this complexity results in a deficiency of GNSS TEC maps that manifests in multiplicative (scale) and additive (bias) errors that vary with time and latitude. The nonlinear dependence of the declination bias with declination in the form of $\delta (x\,\sin \delta +y\,\cos \delta )$ may be caused by the dependence of the GNSS TEC scaling factor on latitude, reported by Li et al. (2019). Radio waves from Southern Hemisphere sources observed at the array located in the Northern Hemisphere propagate through regions in the ionosphere with systematically lower latitudes than from Northern Hemisphere sources. Since a fixed mapping function is used for both computation of TEC maps and computation of the slant ionospheric path delay from these maps, no modification of that fixed mapping function for the data reduction of astronomical data is able to account for this kind of complexity, but as it was shown earlier, it is still possible to mitigate it. The remaining bias can be eliminated by applying the empirical debias correction.

Etemadfard et al. (2021) processed 60 days of VLBI data and estimated not only the TEC for each site using B-splines of the first degree, but also TEC derivatives over longitude and latitude that were considered constant over 24 hr periods. They claim that estimating TEC partial derivatives decreases the discrepancies of the TEC from VLBI to GNSS TEC maps by 36%. Inspired by this result, I introduced estimation of the latitude and longitude TEC gradients in the form of a B-spline in addition to TEC estimation using dual-band data. I considered the weighted rms (wrms) of the differences between the parameterized model of the TEC adjustment to the vertical TEC from dual-band data as a metric of improvement. I did not find a reduction of the rms greater than several percent and abandoned this approach. I should note this result does not disprove the findings of Etemadfard et al. (2021) since they used a different metric.

Expression for the ionospheric contribution (Equation (1)) is an approximation. Hawarey et al. (2005) considered the impact of a higher-order expansion on the group delay, namely proportional to f⁻³. They found that the maximum contribution of the third degree term varied from 3 to 9 ps at 8 GHz, depending on the baseline. This contribution is approximately one order of magnitude less than the contribution of the residual ionospheric errors after taking into account the GNSS TEC maps. It should be noted that the third degree term affects both single-band and ionosphere-free linear combinations of the group delays in two bands and therefore, cannot be retrieved in an analysis of the differences with respect to the dual-band reference solution.

As it was shown before, our ability to model the ionospheric contribution using single-band observations is limited. After applying the declination debias in the data reduction procedure, the declination bias is virtually eliminated. The residual ionospheric contribution causes additional random errors. In order to investigate their impact, I ran four trial solutions. I used ionosphere-free linear combinations of dual-band data and added to them the ionospheric contribution derived from VLBI data scaled to the specific frequency of the trial solution. I computed the contribution of the ionospheric path delay for modified observables using the GNSS TEC maps as if I processed single-band observations and then I estimated source positions. The differences in the source positions from these solutions with respect to the reference solutions are interpreted as the impact of the residual ionospheric errors at different frequencies. The declination dependence of the differences in right ascension and declination for four solutions at frequencies 2.3 GHz (S band), 4.3 GHz (C band), 8.4 GHz (X band), and 23.7 GHz (K band) is shown in Figure 13.

Zoom In Zoom Out Reset image size

These plots help to quantify the additional errors in source positions that would arise if the data set of data from 263 24 hr experiments used for this study would have been observed in one band only. We see that estimation of the TEC biases reduces the ionosphere-driven errors by a factor of 2–4. Still, even after estimation of the TEC biases, the declination dependence of the ionosphere-driven errors remains.

It is known that unaccounted source structure affects source position estimates. In general, the contribution of source structure at higher frequencies is less because jets become optically thin. Therefore, there was an expectation that observing at high frequencies, such as 22–24 GHz (K band), one may obtain more precise source positions. So far, observational evidence do not support that prediction (Lanyi et al. 2010; Charlot et al. 2020). A detailed comparison of the K-band absolute astrometry versus dual-band astrometry of Karbon & Nothnagel (2019) did not reveal an improvement. Therefore, it is instructive to see what is the contribution of the ionospheric errors in K band after applying data reduction from the GNSS TEC maps. Figure 14 shows the additional errors in right ascension and declination due to the residual contribution of the ionosphere. We see that the additional errors in right ascension are about 0.1 mas. Declination errors of Northern Hemisphere sources are at the same level. However, these errors grow with a decrease of declination for sources in the Southern Hemisphere approximately linearly and reach 0.3 mas at a decl. of −40°. Therefore, the unmodeled ionospheric contribution sets the error floor in position accuracy. This error floor is not accounted for in the ICRF3 catalog, and the unaccounted ionospheric contribution is 300% greater than the noise floor in K band adopted by Charlot et al. (2020) according to their Table 6. The potential of K-band astrometry cannot be utilized unless the accuracy of ionosphere modeling will be substantially improved.

Zoom In Zoom Out Reset image size

I investigated how the rms of the residual ionospheric contribution for 45 VLBA baselines from the data set of data from 263 24 hr observing sessions varied with time. I computed three statistics for each observing session: (1) τ_v − b_i , (2) τ_v − τ_g − b_i, and (3) τ_v − τ_g − a_i − b_i , and rescaled them to 8 GHz. Here τ_v and τ_g are the ionospheric path delays from VLBI and TEC maps respectively, a_i is the adjusted TEC bias, and b_i is the clock bias. The statistics are shown in Figure 15. In order to improve readability, the time series was smoothed using a Gaussian kernel with parameter σ = 1 yr. The first statistic characterizes the impact of the total ionosphere on the group delay. The second statistic characterizes the rms of the impact of the residual ionospheric errors on the group delay after applying the a priori GNSS TEC model. The final statistic characterizes the impact of the residual errors after applying the a priori GNSS TEC and adjusting the TEC biases using dual-band VLBI data.

Zoom In Zoom Out Reset image size

It would be instructive to explore to which extent the numerical values found in processing data from the VLBA network are representative to experiments conducted at other VLBI networks and to which extent they are specific to VLBA. I ran two solutions, XIS and XIR. Solution XIS used 63,984 X-band group delays from 36 VLBI experiments in 2016–2019 for absolute astrometry programs at core Southern Hemisphere stations hart15m, hartrao, hobart12, hobart26, kath12m, tidbin64, wark12m, and yarra12m with participation of other International VLBI Service for geodesy (IVS) stations. Solution XIR used 6,600,959 X-band group delays from 2154 geodetic VLBI experiments from IVS programs R1 and R4 (Lambert & Gontier 2006) for 2002–2022. Data from R1, R4, and Southern Hemisphere VLBI experiments are available at the NASA Crustal Dynamic Data Information System. ² In total, 36 globally distributed stations participated in experiments under R1 and R4 programs. These experiments were designed for Earth orientation parameter determination. I computed also reference solutions RIS and RIR that used ionosphere-free combinations of dual-band group delays. Solution parameterization and editing were identical for pairs of XIS/RIS and RIS/RIR solutions.

I compared the source positions from the XIS/RIS and XIR/RIR solutions. In both cases systematic declination errors are evident. I found that scaling that makes the mean declination bias of the XIS solution close to zero is k = 0.78, somewhat smaller than that from solution XIB based on data from the VLBA network. Moreover, estimates of the declinations from the XIS solution require a different debias correction than declinations from the XIB solution.

Figure 16 shows the differences in declination from position estimates of 638 sources with formal uncertainties <0.5 mas from the XIR solution with respect to the reference solution RIR based on data from the R1/R4 networks. Comparing solutions made with different values of the mapping function parameter k, I found that the mean declination bias vanishes when k = 0.75. The shape of the bias and its magnitude is noticeably different than the bias from declinations derived using the group delays from observations at the VLBA network (see Figure 7).

Zoom In Zoom Out Reset image size

The magnitude of the peak-to-peak declination bias from the solution that uses data from the R1/R4 network is 0.62 mas, while the magnitude of the declination bias from the solution that uses data from the VLBA network is 0.85 mas. The regional VLBA network stretches over 90° in longitude, but only within [18°, 48°] in latitude, while the global R1/R4 network stretches over all longitudes and over [−43°, +79°] in latitude. Observations of Southern Hemisphere sources with VLBA are possible in southern sectors only, and the further south a source, the lower the elevation that can be observed. At the same time, a given source can be observed at different azimuths and elevations in a global network.

Figure 17 shows the rms of the mean residual ionospheric contribution at R1/R4 baselines. Comparing it with a similar plot for the VLBA network (Figure 15), we see that the wrms of the ionospheric contribution observed in the VLBA network is a factor of 2–3 lower. This can be explained by both differences in the total ionospheric path delay sensed by the stations in the two networks, and by the amount of the ionospheric path delay that is absorbed in the estimates of the parameters, such as the source declination, causing biases. If a source is observed only at a given azimuth and elevation, most of the ionospheric path delay will be absorbed by the estimated parameters that reduce the residuals. The wider spread in azimuth and elevation a given source is observed, the smaller the share of the ionospheric path delay that will be absorbed by the estimated parameters and propagate to the residuals.

Zoom In Zoom Out Reset image size

We can estimate the scaling factor directly from single-band VLBI data using least squares. I ran twelve solutions using data from the VLBA, southern, and R1/R4 networks separately and using data from all three networks altogether and applying three different flavors of the ionospheric mapping function, ΔH = 0, α = 1.0; ΔH = 56.7, α = 0.9782; and ΔH = 150, α = 0.9782. Table 4 displays the estimates. The differences between the estimates of the scaling factors derived from the observations at the R1/R4 and Southern Hemisphere networks are statistically insignificant, while the scaling factor estimates from the VLBA network are systematically lower, and these differences with respect to the R1/R4 network are statistically significant. It is also remarkable that the scaling factors that minimize the postfit residuals differ from those that eliminate the mean declination bias. This implies that scaling does not eliminate systematic errors due to any deficiency of the GNSS TEC maps, but only alleviates their impact.

Negligible remaining biases and the realistic assessment of errors due to the residual ionospheric contribution allows us to use a mixture of dual-band and single-band group delays in a single least-square solution. Source positions derived from single-band data are in general less precise than those derived from dual-band data because an additional factor affects the uncertainties of group delays. These additional errors can be computed and accounted for in deriving the weights of observables. An increase in the uncertainties of group delay observables propagates to the uncertainties of source positions.

A combined analysis of a heterogeneous data set provides a significant advantage because all the data are used. We can fuse dual-band and single-band data in one data set and use it for estimation of source positions. A fused data set consists of observables of three types:

• 1.  
  Dual-band ionosphere-free linear combinations of group delays.
• 2.  
  Single-band delays of dual-band experiments detected in one band only. Data reduction accounts for the ionospheric path delay computed from the sum of the GNSS TEC and TEC bias adjustments. The uncertainty of the bias adjustment is added in quadrature to the group delay uncertainty for the reciprocal weights of such observables.
• 3.  
  Group delays from single-band experiments. Data reduction accounts for the ionospheric path delay computed from GNSS TEC and applies the debias correction. The uncertainty of the ionospheric model computed from the regression model is added in quadrature to the group delay uncertainty used for computation of the weights of such observables.

Solving for source coordinates using a fused data set provides estimates with minimum uncertainties and minimum correlations between right ascension and declination since it uses all available data provided the frequency-dependent position biases due to source structure are negligible. This condition may be violated for some strong sources with significant structure. An analysis of positional offsets between single-band and dual-band solutions will help to identify such sources (see, for example, Petrov et al. 2011), but these cases are encountered infrequently.

An improvement in modeling the ionospheric path delay when TEC biases are adjusted is quite impressive. However, this requires using additional information about the state of the ionosphere. Specifically, VLBI ionospheric path delays from dual-band data were used for computation of these biases. This information is not available for processing single-band observations. Ros et al. (2000) considered the use of dual-band GPS observations from collocated receivers for analysis of radio astronomy observations. They have shown that the GPS determination of the TEC from ground receivers alone without the use of GNSS TEC maps can be successfully applied to the astrometric analysis of VLBI observations.

There are ongoing efforts to install advanced GNSS receivers within a hundred meters of each VLBA antenna. The use of geometry-free pseudo-ranges at 1.2 and 1.5 GHz in a similar way as I used VLBI ionosphere-free group delays for adjustments of TEC biases promises a similar level of improvement.