1 Introduction

The signal propagation delay due to atmospheric water vapour is an important error source in geodetic Very Long Baseline Interferometry (VLBI). This was realised already during the design phase of the Mark-III system (Shapiro 1976) aiming for centimetre accuracy for global baselines. Efforts were made in order to infer corrections for the wet delay from an independent stand-alone microwave radiometer, often referred to as a Water Vapour Radiometer (WVR), measuring the thermal emission from the sky at frequencies around the water vapour line at 22.2 GHz (Elgered et al. 1991; Kuehn et al. 1991; Emardson et al. 1999; Nilsson et al. 2017).

At that time it also became possible to estimate the atmospheric delays from the VLBI data themselves, as long as the observations were acquired at both high and low elevation angles (Herring et al. 1990; Davis et al. 1991). In spite of the fact that low elevation observations have a larger atmospheric delay a sequence of observations in different directions is essential to improve the geometry and decrease the correlation with the estimated vertical coordinate when solving for the atmospheric delay. There is an operational advantage of not having to keep a stand-alone WVR operating continuously and securing the data quality at each station in addition to the monitoring of the VLBI experiment. Furthermore, the WVR method has an insufficient accuracy during rain (Elgered et al. 1991).

Which of these two methods provides the lowest uncertainty of the geodetic results is an interplay between the accuracy of the WVR delays and the elevation cut-off angle of the VLBI observations (Herring 1986). Employing stand-alone WVRs instead of estimating the delays from the VLBI data, has resulted in small improvements in some cases but not always (Elgered et al. 1991; Kuehn et al. 1991; Emardson et al. 1999; Nilsson et al. 2017). Therefore, today the wet delays at all stations in an experiment are normally estimated from the VLBI data themselves, see e.g. Soja et al. (2015). Stand-alone WVR observations at specific sites are mainly used to verify these estimates, see e.g. Ning et al. (2012), Teke et al. (2013) and Klügel et al. (2019). The introduction of the VLBI Global Observing System (VGOS) (Niell et al. 2018) facilitates higher accuracy of estimated delays from the VLBI data because of an increased number of observations per time unit and thereby a better sampling of the atmosphere and its variability.

Validations of these delay estimates from a stand-alone WVR suffer from an important disadvantage, namely that the WVR and the VLBI telescope sample different atmospheric volumes. The atmospheric volumes observed by a telescope and a stand-alone WVR are illustrated in Fig. 1. A typical average profile of water vapour is exponential with a scale height of 2 km, meaning that approximately 90% of the water vapour exists below 5 km, which is in the near field of the VGOS telescope for elevation angles above around 10\(^\circ \). Line-of-sight radiometry directly with the telescope will sample the same air mass that also dictates the wet delay of the VLBI observations. The standard VGOS receiver has an upper frequency limit at 14 GHz. However, the use of higher frequencies has been proposed, where one motivation is to obtain information on the wet delay (Petrachenko et al. 2009). The idea of using an interferometer dish for sensing of the wet delay has been evaluated at the Very-Large Array (VLA) (Resch et al. 1984; Butler 2000).

We assess a method of using the wideband VLBI receiver, mounted in a particular VGOS telescope (at Onsala) to also function as a radiometer by sampling the sky brightness temperature at different frequencies in order to infer the wet delay. This method requires an extension of the frequency range beyond 14 GHz. It has previously been shown that frequency combinations within the interval 14–24 GHz can yield comparable results to that of a conventional 21/31 GHz channel WVR (Flygare et al. 2018). Recent developments in wideband feed horn and Low-Noise Amplifier (LNA) design has shown that 10:1 frequency ratios are possible (Flygare and Pantaleev 2020), with 6:1 feed having been successfully tested up to 50 GHz (Shi et al. 2017). Covering a 10:1 frequency ratio, a continuous band over 3–30 GHz is conceivable. Advantages of keeping the highest frequency as low as possible is the reduced level of tolerances needed in production, and the lower receiver noise expected. The suggested feed and LNA system used as an example in this study cover up to 24 GHz with successful manufacture and tests. The choice of this system is further discussed in Sect. 4.1.

Our assessment is carried out through simulations using data from the European Centre for Medium-Range Weather Forecasts (ECMWF). The simulations are described in Sect. 2. We address two different situations, for two different algorithm formulations: (1) when there is no liquid water in the atmosphere and (2) when there is cloud liquid water present. The implications of the expected performance for geodesy VLBI are discussed. In Sect. 4 we propose methods for the calibration of the sky brightness temperatures, inferred from the observations with the VGOS receiver. Finally, Sect. 5 contains our conclusions and suggestions for further work.

Fig. 1
figure 1

The concept of a stand-alone WVR and a VGOS telescope illustrating the different sampled air volumes. The WVR senses all the water vapour in its far field whereas the VGOS telescope has the entire troposphere in its near field when observing above a 10\(^\circ \) elevation angle. Note the different horizontal and vertical scales

Full size image
Fig. 2
figure 2

The distribution of observations with respect to the elevation angle for the VGOS schedule VO0076 run at Onsala on March 16, 2020 (left). In total there are 877 observations during this 24 h long experiment

Full size image

2 Simulated observations of sky brightness temperatures

The input data describing the atmospheric properties were taken from the ERA-Interim analysis of ECMWF (Dee et al. 2011) via its web interface (Berrisford et al. 2011). The quantities used are: pressure, temperature, humidity, and liquid water profiles. Data were downloaded for the years 2000–2003, for  00, 06, 12 and 18 UT, which gives 5840 atmospheric profiles, at the highest available height resolution. The position was selected to match the location of the Onsala Space Observatory meaning that the results are valid only for sites with similar weather conditions. It is known that the quality of the liquid water profiles in ERA-Interim is rather uncertain and difficult to model, see e.g. (Stengel et al. 2018). Furthermore, we have found that the Liquid Water Content (LWC) was significantly underestimated for the Onsala site by comparing these values to several years of independent WVR observations. It is of most importance that the true variability in the LWC is covered in the dataset used in the simulation. Therefore, the liquid water profiles in ERA-Interim were multiplied by a factor of 3, which is on the conservative side, not to underestimate the expected errors. Other VLBI sites have not been investigated in this study, but the method is applicable anywhere using ERA-Interim or any other global analysis.

The following assumptions were made about the calculated sky brightness temperatures: (i) a pencil beam (an infinitely narrow beam) calculation represents the antenna temperature; (ii) a monochromatic calculation represents the brightness temperature; (iii) the magnitude of the added uncorrelated noise is varied between 0.1 and 1.0 K (one standard deviation, SD); which is a reasonable interval taking thermal noise, gain variations, and ground noise pickup into account. (iv) instrument errors are assumed to be independent of the observed airmass.

VGOS observations cover a large range of elevation angles. Figure 2 depicts the distribution of observations at different elevations (airmasses). The so-called mapping functions, describing either the ratio between the wet delay and the zenith wet delay (ZWD) or the ratio between the sky brightness temperature and the zenith sky brightness temperature, are weather dependent and become less accurate for low elevation angles. Therefore, we limit our study to up to 6 airmasses, corresponding to an approximate elevation angle of 9.6\(^\circ \). Approximately 90% of all the observations in Fig. 2 are acquired at elevation angles above 10\(^\circ \). Given that radiometry for wet delay estimation is not accurate during rain and that the method is primarily proposed to validate the accuracy of the wet delays estimated from VLBI data, 90% of the observations is deemed a sufficiently large fraction. Furthermore, should the method, at some point, evolve into the primary source for wet delay estimates at sites with little rain observations at low elevation angles may not be used at all (Herring 1986).

Our simulated observations of the sky brightness temperature, using the ARTS software package v.2.3 (Buehler et al. 2018), covers the frequency range from 14 GHz to 24 GHz. The PWR98 model for water vapour attenuation was used (Rosenkranz 1998). See Moradi et al. (2020) for an assessment of water vapour attenuation models. The attenuation of liquid cloud water is taken from Ellison (2007).

We assume a horizontally layered atmosphere with no azimuthal dependence. The elevation angles of the simulated observations are equal to 90\(^\circ \), 30\(^\circ \), 19.5\(^\circ \), 14.5\(^\circ \), 11.5\(^\circ \), and 9.6\(^\circ \). Using the approximation of a flat Earth this corresponds to 1, 2, 3, 4, 5, and 6 airmasses, respectively.

The sky brightness temperature may be mapped from any elevation angle to an equivalent zenith temperature, similarly as is done for the wet delay, see e.g. Niell (1996) and Lagler et al. (2013) describe one of several mapping functions used in VLBI analyses. The difference between a flat earth mapping function and the mapping function for a spherical earth with an atmospheric thickness of 2 km is less than 1% for elevation angles above 14\(^\circ \) and increases to 2% for an elevation angle of 10\(^\circ \).

Fig. 3
figure 3

Simulated observations of the sky brightness temperature in the zenith direction (top) and in the direction of 6 airmasses (bottom). In the graphs to the left the liquid water content has been set to zero, and the graphs to the right include liquid water. The white lines depict the mean temperatures. We can see a saturation effect when airmass \(= 6\), and especially when liquid water is present

Full size image

Using the ARTS software we calculated the wet delay, the sky brightness temperature \(T_\mathrm{{b}}\), and the transmission (opacity, \(\tau \)) for the six different elevation angles. The brightness temperature is written

$$\begin{aligned} T_\mathrm{{b}} = T_\mathrm{{bg}}\, e^{-\tau _\infty } + \int \limits _0^\infty T(s) \, \alpha (s) \, e^{-\tau (s)} \hbox {d}s \end{aligned}$$
(1)

where \(T_\mathrm{{bg}}\) is the cosmic background temperature of 2.7 K, \(\tau _\infty \) is the total opacity through the atmosphere, \(\alpha \) is the frequency dependent total attenuation coefficient of the atmospheric constituents, and T is the physical temperature of the atmosphere. The opacity is defined as

$$\begin{aligned} \tau (s) = \int \limits _0^s \alpha (s') \hbox {d}s' \end{aligned}$$
(2)

Later we will prefer to use opacities instead of the observed brightness temperature in the algorithms for the wet delay (Eqs. 7, 8). By using the effective temperature of the atmosphere defined by

$$\begin{aligned} T_\mathrm{{eff}} = \frac{\int \nolimits _0^\infty T(s) \, \alpha (s) \, e^{-\tau (s)} \hbox {d}s}{\int \nolimits _0^\infty \alpha (s) \, e^{-\tau (s)} \hbox {d}s} \end{aligned}$$
(3)

we obtain the relation between the sky brightness temperature, and the total opacity in Eqs. (1) and (2):

$$\begin{aligned} T_\mathrm{{b}} = T_\mathrm{{bg}}\, e^{-\tau _\infty } + T_\mathrm{{eff}}\, \left( 1 - e^{-\tau _\infty }\right) \end{aligned}$$
(4)

Simulated values of \(T_\mathrm{{b}}\) and \(\tau _\infty \) together with Eq. (4) were used to calculate the effective temperatures at the frequencies and the elevation angles of interest.

Equation (4) and the definition \(\tau _\infty (\varepsilon ) = m(\varepsilon )\tau _z\), where \(\varepsilon \) is the elevation angle, m is the airmass factor, and \(\tau _z\) is the zenith opacity, gives:

$$\begin{aligned} \tau _z = -\frac{1}{m(\varepsilon )}\ln \left( \frac{T_\mathrm{{eff}} - T_\mathrm{{b}}(\varepsilon )}{T_\mathrm{{eff}} - T_\mathrm{{bg}}}\right) \end{aligned}$$
(5)

The dataset has 5840 profiles with a ZWD mean and a standard deviation of 94.3 mm and 45.3 mm, respectively. The mean and the standard deviation of the LWC is 0.16 mm and 0.31 mm, respectively. A second dataset was created with no liquid water in the atmosphere, simply by setting the liquid water profiles identical to zero. The simulated sky brightness temperatures for the two datasets are shown in Fig. 3. Comparing the left and right graphs we note that the contribution from liquid water can be dominant, but we are primarily interested in the contribution from water vapour. We also note that at 14 GHz the sensitivity in the brightness temperature, due to water vapour variability, is weak. This motivates our exclusion of lower frequencies in the following analyses.

3 Performance of wet delay algorithms

We have assessed whether to use either brightness temperatures or opacities in the algorithms for the wet delay. We note from Eq. (4) that when \(T_\mathrm{{b}} \ll T_\mathrm{{eff}}\) the opacity is approximately proportional to the sky brightness temperature, but for high opacities \(T_\mathrm{{b}}\) will saturate. However, when the water vapour in the atmosphere increase neither the wet delay nor the opacities will saturate. The use of opacities led to the best performance in terms of the smallest residuals so these will be used in the following. In a real observation only the observed brightness temperatures, at specific elevation angles, are available. Therefore, in order to use opacities we need a model for the effective temperature.

Table 1 Airmass factors at different elevation angles for the three frequency ranges
Full size table
Fig. 4
figure 4

The residuals for the model of the effective temperature using observations at all elevation angles. Left: without liquid water for frequencies 20–24 GHz. Right: with liquid water for frequencies 14–24 GHz

Full size image

3.1 Model for the effective temperature

We have tried several different formulations in order to estimate the effective temperature using the method of least-squares. The available input data consist of the air temperature at the ground, \(T_\mathrm{{g}}\), the relative humidity at the ground, \(r_\mathrm{{g}}\), the observed brightness temperature, \(T_\mathrm{{b}}\), and the airmass factor, m. We used a similar approach as Ingold et al. (1998) and the following formulation offered the best fit to the data:

$$\begin{aligned} T_\mathrm{{eff}} = a_0 + a_1 \cdot T_\mathrm{{g}} + a_2 \cdot r_\mathrm{{g}} + a_3 \cdot T_\mathrm{{b}} + \frac{a_4}{T_\mathrm{{b}}} + a_5 \cdot m \end{aligned}$$
(6)

where \(a_i\) are the estimated coefficients. The airmass, m, is defined as \(\tau _\infty (\varepsilon )/\tau _z\), where \(\varepsilon \) is the observed elevation and \(\tau _z\) the zenith opacity. For each atmospheric profile, ARTS calculates \(\tau _\infty (\varepsilon )\) for the six given \(\varepsilon \), including the zenith direction, \(\varepsilon \) = 90\(^\circ \). Using the ECMWF dataset for the Onsala site, mean values of m, for these elevations, can be estimated. Table 1 summarises the results for three different cases, one with no LWC and the frequency range 20–24 GHz and two with LWC and the frequency ranges 14–24 and 20–40 GHz. A very weak dependence on the frequency has been ignored. The airmass factors are quite close to the ones given by the flat Earth approximation, implying that \(m \approx 1/\sin (\varepsilon )\), since all elevation angles are above \(\approx \) 10\(^\circ \) and most of the water vapour and liquid water are found below a height of 5 km.

The coefficients in Eq. (6) are derived for the same three cases as for the airmass factors (see Table 1).

Examples of the fit to the model for the 20–24 GHz and 14–24 GHz frequency ranges are shown in Fig. 4 and the parameter values for all three frequency ranges are presented in Table 2. The rms error of the effective temperature estimations are below 3.1 K (or around 1%), compared to 3.5 K as reported by Ingold et al. (1998). Higher orders of Eq. (6) only gives a minor improvement of the fit. Note that the coefficients in Table 2 are tuned to the atmospheric conditions at the Onsala site. The coefficients for another location shall be obtained from ECMWF data for that site.

Table 2 Estimated coefficients for \(T_\mathrm{{eff}}\) in Eq. (6) for the three frequency ranges
Full size table

We will first assess an algorithm using observations in one frequency channel only, thereby optimised for the conditions with no liquid water in the atmosphere. Secondly we will do a similar assessment of a more general algorithm, based on observations at two frequency channels, which can provide reasonably accurate estimates of the wet delay also during cloudy conditions as long as the water drops are much smaller than the observed wavelength. This is the case also for a stand-alone WVR (Westwater and Guiraud 1980).

3.2 One-frequency wet delay algorithm

The one-frequency algorithm used for the ZWD, \(\ell _w\), is of the form

$$\begin{aligned} \ell _w= a_0 +a_1 \cdot p_\mathrm{{g}} +a_2 \cdot \tau _z +a_3 \cdot \tau _z^2 \end{aligned}$$
(7)

where \(p_\mathrm{{g}}\) is the surface pressure, \(\tau _z\) is the equivalent zenith opacity given by Eq. (5), and \(a_0\)–\(a_3\) are the coefficients determined by the method of least squares. An inclusion of the ground temperature was also tested but did not result in an improvement.

Fig. 5
figure 5

Simulations of the expected ZWD rms error (or SD, since no bias error exists) for the one frequency algorithm presented for the equivalent zenith direction for 1 airmass (left), 6 airmasses (middle), and 1–6 airmasses (right) and different receiver rms errors in the observed brightness temperatures. The lower plots zoom in on the frequency range giving the lowest SD. The circles mark the lowest SD at the optimal frequency

Full size image
Table 3 Expected error in the ZWD (SD) using the one-frequency algorithm at the optimal frequency for different airmasses and different errors in \(T_\mathrm{{b}}\)
Full size table

There is a need to have a sufficiently strong signal from the water vapour emission line compared to the receiver noise. Therefore, although a future receiver may cover a larger frequency band, we restrict our simulations to the range of 14–24 GHz.

The results for different frequencies and receiver noise are shown in Fig. 5 for observations acquired through 1, 6, and all (1–6) airmasses. We note that even for an rms error of 0.1 K and observations at 6 airmasses the water vapour signal becomes too weak below around 20 GHz. Please note the different y-scales in each graph. The optimal frequencies, and the resulting rms errors for these cases are presented in Table 3. The optimal frequency appears on one of the sides of the peak of the water vapour emission line. The reason is that most of the water vapour exists in the troposphere where the line is pressure broadened and contributions from the emission line to the sky brightness temperature shall match the contributions from the wet refractivity to the wet delay.

The increased scatter for a frequency towards the centre of the line is illustrated in Fig. 6. Note also that the optimal frequency appears closer to the line centre when the receiver noise increases, and for observations in the zenith direction only, since a stronger water vapour signal is then an advantage. The estimated, site dependent, coefficients for the case when observations are made at all airmasses 1–6 are presented in Table 4. We also note that the optimal frequency is not always on the same side of the emission line. The differences are however hardly significant, as seen in Fig. 5.

Fig. 6
figure 6

Illustrations of the fitted data to the one-frequency algorithm for three different frequencies, on and at each side of the emission line. The left graph refers to the blue line \((\mathrm{SD}\,\, T_b = 0.1 \,\mathrm{K})\) and the right graph the red line \((\mathrm{SD}\,\, T_b = 1.0 \,\mathrm{K})\) in the lower right graph in Fig. 5

Full size image
Table 4 Estimated coefficients for ZWD in Eq. (7) for 1–6 airmasses
Full size table

3.3 Two-frequency wet delay algorithm

The use of two frequencies is motivated by the need to separate the two unknown contributions to the opacity from water vapour and liquid water. Only the former has a large contribution to the wet delay. These shall therefore be chosen at significantly different frequencies to benefit from the different frequency dependence of water vapour and liquid water emissions (see e.g. Wu 1979).

The two-frequency algorithm for ZWD is formulated as

$$\begin{aligned} \ell _w= & {} a_0 +a_1 \cdot \tau _{z1} +a_2 \cdot \tau _{z1}^2 +a_3 \cdot \tau _{z2} \nonumber \\&+a_4 \cdot \tau _{z2}^2 +a_5 \cdot \tau _{z1} \cdot \tau _{z2} \end{aligned}$$
(8)

where again the coefficients \(a_{0-5}\) are determined by the method of least squares and \(\tau _{z1}\) and \(\tau _{z2}\) are the equivalent zenith opacities (Eq. 5) at the two frequencies.

Depending on the frequency band observed in future VGOS receivers we have derived two-frequency algorithms for two different cases: one where the upper frequency is set to 24 GHz, i.e. using K-band, and another where also the Ka-band (up to 40 GHz) is included.

3.3.1 Two-frequency wet delay algorithm K-band

The results for different frequency pairs and receiver noise are shown in Fig. 7, which is a 2-D version of Fig. 5, for observations acquired through 1, 6, and all (1–6) airmasses. The optimal frequency pairs, and the resulting rms errors for these cases are presented in Table 5. Compared to the one-frequency algorithm results, without any LWC, shown in Table 3, the rms errors for the two-frequency data, with LWC, are approximately twice as large. The estimated coefficients for the case when observations are made at all airmasses 1–6 are presented in Table 6. As in the one-frequency algorithm results, we note that the coefficients change significantly for the different levels of receiver noise because a noisier receiver implies that the higher of the two optimal frequencies is closer to the centre of the emission line but the differences in SD is almost insignificant. Again the higher of the two optimal frequencies is not always on the same side of the centre of the emission line.

The one-frequency algorithm improved when, apart from the opacity, also the surface pressure was included (see Eq. (7)). Such an improvement was not seen when testing the two-frequency algorithm. The reason is that the impact of the varying contribution from oxygen is suppressed by the linear combination of the opacities at the two frequencies. Therefore, the surface pressure is not included in Eq. (8).

A three-frequency algorithm was also tested. Adding a third channel gives a negligible improvement for the lowest receiver noise level (0.1 K) and no improvement was seen for the two noise levels of 0.5 K and 1.0 K.

Fig. 7
figure 7

Simulations of the expected ZWD rms error (SD) for the two-frequency algorithm for 1 airmass (left), 6 airmasses (middle) and 1–6 airmasses (right). The receiver noise is 0.1 K (top), 0.5 K (middle), and 1.0 K (bottom) for each row, respectively. The white areas correspond to rms errors larger than the upper limit of the scale and the circles mark the lowest rms error obtained for the optimal frequency pair

Full size image
Table 5 Expected error in the ZWD (SD) and the optimal frequency pair for the two-frequency algorithm using K-band, for different airmasses, and different errors in \(T_\mathrm{{b}}\)
Full size table
Table 6 Estimated coefficients for ZWD in Eq. (8) for 1–6 airmasses
Full size table
Table 7 Expected error in the ZWD (SD) and the optimal frequency pair for the two-frequency algorithm using K- and Ka-band, for different airmasses, and different \(T_\mathrm{{b}}\) errors
Full size table
Table 8 Estimated coefficients for ZWD in Eq. (8) for 1–6 airmasses
Full size table

3.3.2 Two-frequency wet delay algorithm K- and Ka-band

The results using the wider frequency band 20–40 GHz are presented in Tables 7 and 8, following the same structure as was used for the K-band algorithm. Many existing WVRs operate at about 21 and 31 GHz. As seen in column airmass 1–6 in Table 7, the optimal frequencies differ from these. Regarding the lower frequency channel, Fig. 5 shows that SD ZWD have local minima on both sides of the peak of the \(\hbox {H}_2\)O line. Depending on the receiver rms errors a frequency at about 23 GHz can be a slightly better choice than 21 GHz. Regarding the higher frequency channel, Table 7 gives somewhat higher frequencies than 31 GHz, but SD ZWD values are quite low in the whole 31–38 GHz frequency range. A motivation to observe around 31 GHz is that there is a protected frequency band for passive radiometry from 31.3 to 31.8 GHz.

The lower frequency channel (23–24 GHz) is most sensitive to water vapour and the higher frequency channel (31–37 GHz) is most sensitive to cloud water. When the receiver noise increases, the optimal frequencies are found in frequency ranges where the emission from the atmospheric water is stronger. This is clearly seen in the lower left sub-plot in Fig. 7 where one of the two frequencies is found at the peak of the water vapour emission line. Similarly, to have a better signal-to-noise ratio, the optimal frequency of the high frequency channel increases with increased receiver noise.

Finally, a comment on the empirical factor of 3 used to increase the LWC in the ERA-Interim dataset. The expected algorithm errors in Table 7 decrease by \(\approx \) 10% when the factor is not used. This is not a large difference given the other uncertainties. As mentioned in Sect. 2 the factor was introduced not to underestimate the algorithm errors.

3.4 Intermediate discussion

One obvious application of WVR data is to assess the accuracy of the estimated time series at the different VGOS sites. To put that into a context we need to discuss the ZWD uncertainty when it is estimated from the VLBI data. This is a complex issue where the result depends on the geometry of the observations and how the analysis is set up. Parameters to vary are for example the temporal resolution of the estimated time series of the ZWD and linear horizontal gradients (if they are estimated). In addition it will also depend on the constraints used for these parameters. As an example, in the CONT14 experiment estimates of the ZWD was made every 30 min and the gradients every 6 h, using the constraints for the variability of 15 mm \(\hbox {h}^{-1}\) for the ZWD rate segments and 2 mm \(\hbox {d}^{-1}\) for gradient rates. This resulted in a formal error of 1.7 mm for the ZWD (Elgered et al. 2019). The formal error is however in general always underestimating the real uncertainty. Heinkelmann et al. (2011) addressed this question by comparing the ZWD results from ten different analysis centres and concluded that overall there was an underestimation with approximately a factor of 2. Their results also show that the formal uncertainties of the ZWD from the different centres for the same experiment also vary with a factor of 2, depending on the assumptions made when setting up the analysis.

Because an estimated ZWD is strongly correlated with the estimated vertical coordinate of the VLBI station, an improvement in the ZWD will also result in a more accurate vertical coordinate. The correlation depends on the geometry of the VLBI network. An approximate relation assuming that the cutoff angle of observations is 10\(^\circ \) is that an error in the ZWD when estimated from the VLBI data is multiplied by a factor of 3 to give the associated error in the vertical coordinate. This factor is closer to 1 for the error of independent ZWD estimates, e.g. from a WVR, and no ZWD are estimated from the VLBI data, and observations below elevation angles of 10\(^\circ \) are ignored (Herring 1986). To summarize, we find that the expected ZWD uncertainties from the WVR and from VLBI data are comparable. The WVR data have the potential to improve the geodetic results during non-rainy conditions eliminating the need to solve for the ZWD and horizontal gradients in the VLBI data analysis. The first straight-forward application is, however, to use the WVR results to assess the quality of the atmospheric estimates from the VLBI data.

In order to estimate the turbulent atmosphere the temporal resolution of the estimated ZWD as well as the gradients should be as high as a few minutes. The WVR is very useful because it provides independent equivalent ZWD in the different directions without introducing any constraints of in terms of dependence on previous estimates. Such measurements can then be used to assess the corresponding estimates from the VLBI data.

This study used atmospheric data for the Onsala site. Algorithm errors depend only slightly on the weather conditions at the site, provided that a site specific algorithm is derived. The main difference of algorithms derived for dry or wet sites, e.g. in Alaska and in Florida, is the values of the optimal frequencies (Elgered 1993). In short, when there is a lot of water vapour in the atmosphere the emission is strong compared to the receiver noise. It is then an advantage to observe further above, or below, the peak of the line, where the emission from the water vapour profile correlates better with the profile of the wet refractivity.

Finally we note that the overall accuracy of WVR-determined ZWD is heavily depending on the uncertainty in the measured sky brightness temperatures. For the three simulated algorithms, including observations from 1 to 6 airmasses, the uncertainty is reduced by a factor \(\approx 3\) if the receiver rms error is decreased from 1.0 K to 0.1 K. The stability of the receiver is therefore of great interest and the next section discuss the state-of-the art of broadband receivers for radio astronomy telescopes.

4 Calibration of the sky brightness temperature

So far we have presented simulation results for the estimation of the wet delay based on one or two sky brightness temperatures observed with the receiver in a VGOS telescope. In this section we will discuss requirements on the receiver and the telescope in order to realise such an application.

A traditional stand-alone radiometer, collocated with a VLBI telescope, is typically equipped with a horn antenna, characterized by low sidelobes. Additionally, in order to infer the sky brightness temperature each channel has at least two reference loads with stable temperatures to determine the two unknowns, the gain and the system noise temperature. Furthermore, the absolute accuracy of the reference load that has a different temperature than ambient requires an external calibration because of temperature gradients in lossy waveguides which are difficult to model. This is typically accomplished through the use of elevation scans, also referred to as sky dips, or tip curves, see e.g. Elgered and Jarlemark (1998).

Fig. 8
figure 8

The wideband feed horn used to simulate ground-noise pickup when mounted in one of the VGOS antennas at Onsala

Full size image
Fig. 9
figure 9

The simulated antenna pattern at 21 GHz when the horn is mounted in one of the VGOS antennas at Onsala. A polar plot with an angle range of 0\(^\circ \)–105\(^\circ \) from the main lobe where the central part (\(\le 2\) \(^\circ \)) of the beam is not shown for clarity reasons (top) and a cut 270\(^\circ \)–90\(^\circ \) (left) and 180\(^\circ \)–0\(^\circ \) (right)

Full size image
Fig. 10
figure 10

The convolution between the VGOS antenna pattern and the sky brightness distribution. Left: the observed brightness temperatures from the VGOS antenna and a pencil beam. Right: the VGOS antenna ground pick-up and the difference between the VGOS antenna and a pencil beam

Full size image

4.1 Gain and system temperature of the receiver

The present standard VGOS receivers do not have two reference loads. Neither do they cover the frequency band around the water vapour emission line. To provide two reference temperatures for calibration, an internal noise injection and an external absorber load could be implemented. However, the difficulty to stabilise the temperature of the reference load, and to mechanically design and fit such a system on the telescope, makes the use of an external absorber less interesting. A more feasible solution would be to implement two separate internal references. Recent LNA designs have successfully included a noise-injection coupler directly within the wideband LNA (Pellegrini et al. 2021). With a conventional cryogenic system, a directional coupler with two injection points between the feed and the LNA could be used as reference loads. A stabilised noise diode would then be used as a noise source with a cold attenuator reference output. If a single injection point is used potentially a wideband ridge-waveguide switch system could be implemented.

We also note that with frequent changes of the observed extra-galactic sources tip-curve observations will follow automatically. A significant range of different elevation angles will be covered without having any impact on the VLBI observations. In fact, today’s scheduling procedures already require a sequence of different elevation angles in order to solve for the wet delay from the VLBI data. The observations of the reference loads is carried out during a few seconds and can be scheduled to occur during the time periods when the telescope is slewing between the sources.

A remaining question is how a VGOS telescope compares to the horn antennas of stand-alone WVRs. We have addressed this issue through simulations. We assume that the VGOS antenna is equipped with a quad-ridge flared horn (QRFH) feed identical to the one developed for the advanced instrumentation programme of the Square Kilometre Array (SKA) project (see Fig. 8) covering the frequency band 4.6–24.0 GHz. Details about the feed horn design and its measured performance is given by Dong et al. (2019). Simulations have shown that the horn is also a good candidate feed for the VGOS reflector geometry. Measured receiver noise of the feed with the LNA, cooled to 20 K, gives an estimated System Equivalent Flux Density (SEFD) comparable to that of current VGOS receiver systems (Flygare et al. 2018). Another reason why this frequency range can be a suitable candidate for future VGOS receivers is the allocation of 5G frequency bands for cellphone communication, potentially introducing more radio frequency interference (RFI). The frequency ranges of 3.4–3.8 GHz and 24.25–27.5 GHz are currently considered the main bands to be allocated in Europe, with similar bands considered in the rest of the world (International Telecommunication Union 2020). Together with already strong RFI sources below 3 GHz, the future potential of a 4.6–24 GHz receiver is interesting with the inclusion of line-of-sight radiometry.

4.2 Modelling input brightness temperature

For a VGOS telescope with a diameter of 13 m the half-power beamwidth at 21 GHz is \(\theta \sim 70\times \lambda /D=0.08^\circ \). The small angular spread for the telescope can be interpreted as the beam having a near-cylinder shape over the range, starting at the telescope and ending when the spread is equal to the antenna diameter. This means that we have a geometry as illustrated in Fig. 1. To accurately model the difference in observed brightness temperature between the assumed pencil-beam and the VGOS antenna pattern, the antenna temperature is calculated as the convolution of the antenna pattern with the surrounding brightness temperature distribution (Balanis 2005, p. 106). This distribution is strongly elevation dependent with varying contributions from cosmic sources, atmosphere and ground temperature (Cortés-Medellín 2007). For the frequency range of 14–24 GHz the cosmic sources are ignored in our calculations.

The antenna pattern for the VGOS antenna illuminated by the wideband QRFH, shown in Fig. 9 for 21 GHz, is normalised by integrating the power over the sphere. The polar plot shows the 0\(^\circ \)–105\(^\circ \) range from the main lobe. The central part of the beam is omitted to be able to resolve the far-out side-lobes in the 65\(^\circ \)–105\(^\circ \) range. The main and back lobes are seen in the two cuts below the polar plot covering the entire ±180\(^\circ \) range.

Figure 10 depicts the convolution between the VGOS antenna pattern and a sky brightness distribution at 21 GHz with an angular resolution comparable to the resolution of the antenna pattern. The result from a pencil beam is shown as comparison. In the zenith direction the VGOS antenna observes a larger antenna temperature than the pencil beam, which has no side lobes, since the VGOS antenna has contributions from lower elevation angles as well as the ground pick-up. The difference between the VGOS and the pencil beam decreases with decreasing elevation angle. At low elevation angles the VGOS antenna observes a lower brightness temperature than the pencil beam since the ground pick-up does not compensate for the contributions from the cold sky at high elevation angles. The curvatures at about 75\(^\circ \) elevation angle in the right plot is explained by the far-out side-lobes located 65\(^\circ \) to 105\(^\circ \) from the main lobe. As seen in Fig. 10, the simulated ground pick-up is quite stable in the 10\(^\circ \)–90\(^\circ \) elevation range, which is the interval for the simulated performance of the wet delay algorithms. It seems reasonable to assume that a correction to the pencil beam assumption (black curve in Fig. 10) can be modelled with an accuracy of a few tenths of a degree Kelvin.

5 Conclusions and suggestions for further work

We have described a technique providing independent ZWD using radiometry in a VGOS telescope. This could be considered if and when broadband receivers at higher frequencies are being implemented for geodetic VLBI. Three different algorithms have been studied: (i) a one-frequency algorithm based on radiometric observations up to 24 GHz to be used during cloud-free conditions, (ii) a two-frequency algorithm also using observations up to 24 GHz, and (iii) a two-frequency algorithm using observations up to 40 GHz. For the first two algorithms an existing wideband feed and LNA up to 24 GHz has been evaluated through simulation for the VGOS telescope. Results indicate that ground-noise pickup can be modeled with an accuracy of \(\approx 0.5~\hbox {K}\). The predicted formal error for these three cases, using optimal frequencies, and assuming receiver errors (SD) of 1 K averaged over observations spread over elevation angles between 10\(^\circ \) elevation and zenith, are 2.7 mm, 6.1 mm, and 3.8 mm, respectively. If the receiver errors are reduced to 0.1 K, these errors approximately decrease by a factor of 3, meaning that the receiver stability is likely to be the limiting factor.

The ZWD from the WVR cannot be the only method to handle the atmospheric calibration, unless the community accepts to ignore all observations acquired during rain. We think that the optimal use of the WVR method is to combine it with the present approach, i.e., the estimation of the ZWD and horizontal linear gradients from the VLBI data themselves.

A logic first application is to assess the accuracy of estimated ZWD and horizontal gradients. Compared to the stand-alone WVR, which has been used for this purpose at many sites and experiments (Teke et al. 2013; Nilsson et al. 2017), the main advantage of having the WVR in the telescope is that it will sample the same atmospheric volume that causes the delay of the VLBI signal. This is especially useful when assessing estimates of the horizontal gradients, as these quantities are often based on an overly simple model for a turbulent atmosphere, characterized by, e.g., convective processes. An outcome of such assessments is that one can determine optimal temporal resolutions and constraints for the estimated ZWD and horizontal gradients estimated from the VLBI data. Note that the outcome of the optimizations will be site dependent because of different atmospheric conditions.

Although the algorithms derived here for the wet delay and the effective temperature of the atmosphere are based on the weather conditions at the Onsala observatory, Sweden, the method can be used at any other site since the input data used have global coverage. The resulting ZWD uncertainties due to algorithm errors will also be similar provided that the optimal frequencies are determined for each site and its given receiver stability.