High temporal resolution wet delay gradients estimated from multi-GNSS and microwave radiometer observations

Abstract. We have used one year of multi-GNSS observations at the Onsala Space Observatory on the Swedish west coast to estimate the linear horizontal gradients in the wet propagation delay. The estimated gradients are compared to the corresponding ones from a microwave radiometer. We have investigated different temporal resolutions from 5 min to one day. Relative to the GPS-only solution and using an elevation cutoff angle of 10° and a temporal resolution of 5 min the improvement obtained for the solution using GPS, Glonass, and Galileo data is an increase in the correlation coefficient of 11 % for the east gradient and 20 % for the north gradient. Out of all the different GNSS solutions, the highest correlation is obtained for the east gradients and a resolution of 2 h, while the best agreement for the north gradients is obtained for 6 h. The choice of temporal resolution is a compromise between getting a high correlation and the possibility to detect rapid changes in the gradient. Due to the differences in geometry of the observations, gradients which happen suddenly, are either not captured at all or captured but with much less amplitude by the GNSS data. When a weak constraint is applied in the estimation of process, the GNSS data have an improved ability to track large gradients, however, at the cost of increased formal errors.


inferred by the WVR, we subtracted the hydrostatic component computed from the reanalysis product of the European Centre for Medium-Range Weather Forecasts (ECMWF), ERA5, from the total gradient to get the GNSS wet gradient. The hydrostatic gradients at the site are much less variable compared to the wet ones, and especially for time scales of minutes to hours (Elgered et al., 2019).
The data processing was run for three different elevation cutoff angles (3 • , 10 • and 15 • ). For each elevation cutoff angle 70 we used four different combinations of GNSS constellations in the processing: GPS only (G), GPS + Glonass (GR), GPS + Galileo (GE), and GPS + Glonass + Galileo (GRE). Due to limitations in the receiver capacity, not all BeiDou observations were recorded. Therefore we decided not to include BeiDou data in our GNSS data processing.
An example of the sky coverage of the observations for different GNSS constellations, applying an elevation cutoff angle of 3 • , is shown in Figure 1 for the ONS1 station. The three systems show a similar geometry for the observations. At this latitude 75 a significant part of the sky, just north of the zenith direction, is never sampled. The Glonass satellite orbits have a higher inclination angle implying that a smaller part of the sky is not sampled compared to GPS and Galileo. As a consequence there are more observations from Glonass to the north, especially below the elevation angle of 20 • . It is therefore interesting to study how this difference will affect the quality of the estimated gradients.  angle. Some days have less or no data due to receiver problems, especially from 13 to 15 July, where no Galileo observations were recorded. Averaged over the year, for the GPS system, the number of observations drops about 11 % and 23 % when the elevation cutoff angle changes from 3 • to 10 • and 15 • , respectively. For the Galileo system, the corresponding decrease in the number of observations are 10 % and 22 %, respectively, while for Glonass, the corresponding values are 8 % and 17 %, respectively.

Water Vapour Radiometer (WVR)
The WVR, shown in Figure 3, is located close to the GNSS sites, 9 m from ONSA and 59 m from ONS1, and with height differences of less than 2 m. The WVR was designed in order to provide independent estimates of the wet propagation delays for space geodetic applications. It measures the sky brightness temperature on and off the water vapour emission line at 22 GHz.
More detailed specifications are given by Elgered et al. (2019).

90
Starting in January 2019 the observations were scheduled in 5 min long cycles with the ambition to sample the whole atmosphere at elevation angles above 25 • , which is illustrated in Figure 4. Data were acquired throughout 2019, except from mid August to early October because of a failure caused by a thunder storm.
A four-parameter model was used to estimate the mean ZWD, a linear trend in the ZWD, and east and north linear horizontal gradients over 5 min (Davis et al., 1993). Before the model was applied all observations during rain and with a liquid water   content larger than 0.7 mm were removed. Thereafter, the quality of the data was assessed through manual editing. Gain jumps occurring sometimes at the beginning of a 5 min cycle were frequent during the whole year. When such a jump was identified one or several complete 5 min cycles were removed. Thereafter, when the model was applied, we required that at least 40 of the 52 observations were available in the 5 min cycle period. Thereby we eliminated the possibility that the large gain jumps could have an impact on the estimated ZWD trends and gradients for the 5 min cycle since it is not synchronised with the estimation 100 period. Smaller gain jumps may still have degraded the accuracy of the estimated gradients. An overall reason for applying this strict editing was that the primary goal was to have accurate gradients from the WVR rather than a statistical characterization of the specific atmospheric conditions at the site.
Finally it is noted that the WVR estimates are completely independent of the corresponding estimates from the GNSS data.
The estimated ZWD from the WVR data is shown in Figure 5. The seasonal dependence is clearly visible as well as a large 105 short term variability.

Estimated gradients and their formal errors
Before carrying out comparisons of the gradients estimated from different GNSS solutions with the WVR gradients we investigate the characteristics of the input data. Table 1   formal errors is significantly larger than those from GNSS. This is due to a varying uncertainty in the measured sky brightness temperatures. These variations are taken into account in the following comparisons by using of the formal errors of the GNSS and the WVR gradients when calculating the weighted root-mean-square (WRMS) differences and correlations.

Comparison of gradients from GNSS and WVR
We first carry out comparisons of the gradients estimated from the different GNSS constellations and using the three different 120 elevation cutoff angles, presented in Table 1, with the WVR gradients. Even though the gradients estimated from both the GNSS and the WVR data have a temporal resolution of 5 min, the estimates are not centered at exactly the same time epochs.
It is therefore necessary to synchronize the time series to compare the gradients. We first present results with the highest available temporal resolution of 5 min where the WVR gradients were interpolated to the epochs in the GNSS time series using the temporal Gaussian filter as described by Ning et al. (2012) with a full width at half maximum (FWHM) of ± 2.5 min.
125 Table 2 shows the WRMS differences and correlations of the east and the north gradients for the whole year of 2019.
As expected, using the data from multi-GNSS, we note a significant improvement (an increase in the correlation of up to 20 % and a maximum reduction of the WRMS difference of 11 %). Our interpretation is that the geometry of the observations is improved when observations from additional systems are added, especially in the south-north direction (see Figure 1). We also note that the GRE solution, in general, gives the best agreement with the WVR gradients.  For the GPS-only solution, the best correlation is obtained for the elevation cutoff angle of 3 • , especially for the north gradients. This is however not the case for the multi-GNSS solutions (GR, GE, and GRE) where the best correlation for the east gradient is obtained for the elevation cutoff angle of 10 • . For the GRE solutions, the correlations given by the 3 • and 10 • solutions are similar for the north gradients. These results indicate that the choice of elevation cutoff angle is a compromise between having a good geometry and avoiding elevation-angle-dependent systematic errors, e.g., multipath effects. Related to 135 this is the error introduced by the assumption that the turbulent atmosphere may be modelled with just a linear gradient when the elevation cutoff angle of 25 • has to be used for the WVR observations in order to avoid ground noise pickup. As was presented in Table 1 the estimated size of the gradients was in general increasing when using a decreasing sky coverage of the observations.
In spite of that the WVR and the GNSS sample different parts of the sky it is noted that the agreement becomes worse 140 for all GNSS solutions when the 15 • elevation cutoff angle is used. Our interpretation is that it is because many important observations are removed, especially in the north direction. We also note that even though there are slightly more observations from the Glonass system contributing to the south-north direction, especially below the elevation angle of 20 • (see Figure 1), the GR solution does not give a significantly better agreement with the WVR than the one for the GE solution.
In order to study any seasonal variability we compare the estimated gradients from the GNSS and the WVR for each month.

145
The GNSS gradients are obtained using different constellations and a cutoff angle of 3 • (see Figure 6). The change in the correlation is large from month to month and these changes seem to be related to the amplitude of the ZWD (see Figure 5). In general, a large ZWD variability results in a larger dynamic range for the gradients and consequently also a larger correlation. Figure 6 shows higher gradient correlations for June and July which is consistent with the results given by Elgered et al. (2019).
To investigate how well the GNSS data capture large gradients, we carried out the same comparison using data from June and 150 July only. Another reason for focusing on these two months is to include the ONSA GNSS station for validation purposes.
Multi-GNSS data from ONSA started to be acquired only in April 2019.  The results are summarized in Table 3. Because of the larger gradients the WRMS differences increase in all cases. However, at the same time we see higher correlations between GNSS and WVR gradients compared to the whole dataset (see Table 2).
The two GNSS stations (ONS1 and ONSA) show very similar agreements with the WVR gradients. This is partly expected 155 since they are located close to each other and therefore the gradients from the two sites are estimated based on the same observational directions and are affected by common error sources, such as orbit errors. However, it is of interest to note that the different antenna mountings (see Figure 3) do not have a significantly different impact on the estimated gradients.
Therefore, we continue to use only ONS1 GNSS data in the following.

Effective temporal resolutions from 5 min to 24 h
As mentioned in the introduction GNSS gradients have been compared to other independent estimates over different time scales and temporal resolutions. We therefore averaged the GNSS and the WVR gradients by applying a gaussian window with different FWHM, from ± 2.5 min to ± 720 min for further comparisons. As described in Section 3.2, when a FWHM of ± 2.5 min is used, only the WVR gradients were interpolated to the epochs in the GNSS time series. For all FWHM larger than ± 2.5 min both the GNSS and the WVR gradients are interpolated to the epochs at 0, 5, 10... 55 min after the hour. The 165 requirement to calculate a value at a specific epoch is that at least half of the original data points (with a 5 min resolution) exist within the FWHM. In the following we will refer to and use the FWHM as an effective temporal resolution, ∆t eff , from 5 min to 1440 min (one day), although the time series will still have a value every 5 min.
The resulting WRMS differences and correlations are shown in Figure 7. It is clear that the WRMS differences decrease when the ∆t eff increases and more variations of the gradients are averaged out. Over all GNSS solutions, the highest correlation for 170 the east gradients is obtained when a ∆t eff of 2 h is used, while for the north gradients, the best agreement is seen for a ∆t eff of 6 h. Figures 8 and 9 depict the gradients estimated from the GNSS data and the GRE solution for a 3 • elevation cutoff angle, against the gradients obtained from the WVR data for four different values for ∆t eff (5 min, 2 h, 12 h, and 24 h). It illustrates that even when the gradients are averaged over one day, i.e., applying ∆t eff of 24 h, there are substantial variations left which still give a clear correlation between the GNSS and the WVR data. 175 We also study a specific event of short lived gradients, associated with rapid changes in the ZWD, starting from 0 h, 23 July (see Figure 10). There is a passage of a warm front, indicated by a sudden increase in the ZWD, during the late evening of 24 July. As a result, we see a large gradient towards the west direction. They are detected by both GNSS and WVR but the amplitudes of the gradients from GNSS are much smaller. For the north direction, there are also some large gradients, i.e., at 8 and 17 h of 23 July, which are detected by both GNSS and WVR. Also in this case the amplitudes from GNSS are smaller.

180
In addition, the multi-GNSS solutions have a slightly higher possibility to capture sudden short lived gradients.
The results from the comparisons using different ∆t eff indicate that due to the poor geometry of the GNSS observations (especially in the south-north direction), if the gradient happens suddenly, the GNSS data do not capture the full picture of the gradients. When a small ∆t eff is used, all gradients are kept including the ones which are not correctly detected by the GNSS data. As a result, the correlation between the GNSS and the WVR data is deteriorated and this is the case when a ∆t eff of 5 min 185 is applied. However, when ∆t eff is too large (i.e., 24 h), the larger gradients which are captured correctly by both the GNSS and the WVR data will also be averaged out. The range of variations and the correlation decrease.   Table 4. The changes in WRMS differences and correlations of east and north gradients when a weaker constraint is applied in the GNSS data processing. The changes are the results from the constraint of 1.0 mm √ h −1 relative to the 0.3 mm √ h −1 (in Table 3).

Different constraints for gradient variability
All GNSS-derived gradients were so far estimated using a random walk model with a constraint value of 0.3 mm √ h −1 (see Section 2.1). This value may be too small, especially when we have more observations from multi-GNSS constellations, to 190 allow the GNSS data to detect sudden large gradients. In order to investigate this issue, we have processed the GNSS data from June and July again with the 5 min temporal resolution, applying a weaker constraint of 1.0 mm √ h −1 . The changes in WRMS differences and correlations, relative to the solution using the constraint value of 0.3 mm √ h −1 , are shown in Table 4 where the GNSS gradients are estimated applying the elevation cutoff angles of 3 • and 10 • . The formal errors obtained when using the weak constraint of 1.0 mm √ h −1 is about twice as large compared to the ones obtained when using the constraint of When a weak constraint is applied with an elevation cutoff angle of 3 • , the WRMS differences increase for both the east and the north gradients while the correlations are almost the same. We note that for the 10 • solution, the east gradients obtained from using the weak constraint values result in smaller WRMS differences compared to the WVR gradients. A slight improvement is also seen in the correlations for both the east and the north gradients. We interpret the result as the compromise between 200 capturing large gradients and including more noise. When a low elevation cutoff angle is used, the GNSS measurements will be more sensitive to the noise from the environment, i.e., multipath effects. In addition, the GNSS signals will have a lower signal-to-noise ratio due to the longer path through the atmosphere. When a weak constraint is applied, individual observations will have a larger influence on the estimated gradients.
More details are seen in Figure 11, depicting the time series of the gradients for the two and a half days, starting at 0 h, 205 23 July, from the WVR and the GNSS data based on the GRE solutions. There are several peaks, i.e., large gradients, shown for these two and a half days, i.e., east gradients at 21 h of 24 July and north gradients at 8 and 17 h of 23 July. When a weak constraint is applied, there is a clear improvement in tracking those larger gradients when the elevation cutoff angle of 10 • is used. This is not the case when the lower elevation cutoff angle of 3 • is used, possibly because the sampled atmosphere is more different compared to that observed by the WVR.  We have estimated linear horizontal gradients using one year of data acquired from the GNSS site ONS1 located on the Swedish west coast. The GNSS-derived gradients were compared to the ones obtained from a collocated WVR. Overall the multi-GNSS solutions, i.e., combinations of GPS, Glonass, and Galileo, show small but significant improvements with the WVR gradients compared to the GPS-only solution (Tables 2 and 3).

215
For the GPS-only solution, the best agreement, in terms of the correlation coefficient with the WVR gradients, is obtained when using an elevation cutoff angle of 3 • . For the multi-GNSS solution using all three constellations, the best agreement with the WVR data is obtained for the solution with an elevation cutoff angle of 10 • . The difference is largest for the east component which has the better sky coverage (Figure 1 and Tables 2 and 3). This indicates that if there are a sufficient number of observations, the low elevation observations are not that important. This is especially true when the comparison is made to a 220 WVR using observations evenly spread over the sky above an elevation angle of 25 • . It is also an indication that a linear model for horizontal variations in the wet refractivity does not describe the turbulent atmosphere well during all conditions.
We investigated different effective temporal resolutions, ∆t eff , of the compared time series. For all GNSS solutions, the highest correlations obtained for the east and the north gradients are for a ∆t eff of 2 h and 6 h, respectively (Figure 7). When these ∆t eff are applied, strong gradients of short duration detected by the WVR, but not by GNSS, are averaged out and as 225 a result the correlation increases. When estimating GNSS gradients the choice of ∆t eff is a compromise between getting a high correlation and loosing track of rapid gradient variations. However, when ∆t eff is even larger, e.g., 24 h, all gradients are further averaged and the dynamic range of gradient size and the correlation decreases (Figures 8, 9, and 10).
Furthermore, weakening the constraint used when estimating the GNSS gradients from 0.3 mm √ h −1 to 1.0 mm √ h −1 helps the GNSS data to track short-lived gradients, approaching a time scale of 5 min, however at the cost of increased formal errors 230 (Table 4 and Figure 11).
Possible improvements to study in similar future work would be to include BeiDou observations and use a WVR with a better stability. It would also be of interest to carry out a similar study at low latitude GNSS sites where the sky coverage is better and perhaps also the atmosphere is more variable.
Acknowledgement. We thank Tobias Nilsson for providing the ZHD gradients from the ERA5 data.
The ERA5 data are accessible from https://rda.ucar.edu/datasets/ds633.0/. The estimated gradients from the GNSS and the WVR data will be registered and archived by the Swedish National Data Service (SND) if and when the paper is accepted for publication. Before that these data are made available by the authors.