the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Tong Ning

### Gunnar Elgered

We have used 1 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 1 d. 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 of
detecting 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.

An accurate modelling of the atmospheric effects on GNSS observations is relevant for both geodetic and meteorological applications, in forecasting as well as in climate research. In geodetic applications the standard method is to estimate an equivalent zenith propagation delay together with a linear horizontal gradient. Early results showed an improved repeatability for the estimated coordinates when estimating gradients using GPS data (Bar-Sever et al., 1998; Meindl et al., 2004). A recent study (Zhou et al., 2017) found that estimating gradients with a temporal resolution of 1 h can achieve even better positioning performance than strategies where gradients are estimated with resolutions of many hours up to 1 d.

In meteorological applications the zenith total delay (ZTD) and horizontal gradients may be assimilated directly into the forecasting model; see e.g. (Zus et al., 2019). Inferred values of the zenith wet delay (ZWD) and the integrated water vapour (IWV) may be used to study long-term trends (Baldysz et al., 2018). Linear horizontal gradients estimated from the GNSS have been used to study specific meteorological conditions. For example, on the island of Corsica, where the physical meaning of gradients in coastal areas with a steep topography was studied (Morel et al., 2015), and in Texas, USA, where significant gradients during Hurricane Harvey were reported (Graffigna et al., 2019).

The quality of the estimated gradients has been assessed by comparisons to independent measurements, such as using a microwave radiometer (in the following this instrument is referred to as a water vapour radiometer, WVR), the space geodetic technique of very-long-baseline interferometry (VLBI), and numerical weather models.

Such an assessment was carried out by Elgered et al. (2019), where the GPS-derived gradients were compared with the ones obtained from WVR, VLBI,
and the European Centre for Medium-Range Weather Forecasts (ECMWF) analyses. The results show that the best agreement is obtained when an elevation
cutoff angle of 3^{∘} is applied in the GPS data processing in spite of the fact that the radiometer did not observe below 20^{∘}. They
also found that a homogeneous and frequent sampling of the sky is a critical parameter for gradient estimation. Using multi-GNSS observations instead
of GPS only, Li et al. (2015) found a significant increase in the correlation coefficient from below 0.5 to about 0.6 when compared to the
gradients computed from the ECMWF reanalysis product. The corresponding decrease in the root-mean-square (RMS) difference in the gradients was
25 %–35 % for multi-GNSS processing. The temporal resolutions of such comparisons are to our knowledge so far limited to 1 h for WVRs
(Lu et al., 2016), 2 h for VLBI (Steigenberger et al., 2007), and 1 h for numerical weather models (Kačmařík et al., 2019).

The aim of this study is to assess the quality of estimated gradients from multi-GNSS observations with temporal resolutions as high as 5 min using independent WVR data. Section 2 describes how the gradients are estimated from the GNSS and the WVR data. In Sect. 3.1 we present the gradients estimated for different GNSS constellations and different elevation cutoff angles. These are thereafter compared to the WVR data, first in Sect. 3.2 for the highest temporal resolution of 5 min and then in Sect. 3.3 over timescales up to 1 d, both for the entire data set and for a specific event of short-lived gradients associated with rapid changes in the ZWD. In Sect. 3.4 we also study the impact of using a weaker constraint for the random walk process of the GNSS gradient time series. Finally, Sect. 4 gives our conclusions.

## 2.1 GNSS

We have analysed 1 year (1 January–31 December 2019) of ground-based GNSS observations acquired from one station (ONS1) located at the Onsala Space Observatory, on the west coast of Sweden. For comparison purposes we also used GNSS data from June and July 2019 acquired at the collocated station ONSA. The data processing was carried out using GipsyX v.1.5 (https://gipsy-oasis.jpl.nasa.gov/gipsy/docs/releaseNotes-GipsyX-1.5.pdf, last access: 12 August 2021) with the precise point positioning (PPP) strategy (Zumberge et al., 1997). The input to the processing was ionosphere-free linear combinations formed by acquired GNSS phase-delay observations, while the output included station coordinates, clock biases, and tropospheric parameters. The final multi-GNSS orbit and clock products used were provided by the Center for Orbit Determination in Europe (CODE) (Prange et al., 2020). An ocean tide loading correction using the FES2004 model (Lyard et al., 2006) was applied, while no atmospheric pressure loading corrections were used. The absolute calibration of the phase centre variations (PCVs) for all antennas (from the file igs14_2101.atx) was implemented (Schmid et al., 2007). We used the Vienna Mapping Function 1 (VMF1) (Boehm et al., 2006) to map the zenith delay, and the gradient mapping function was the one presented by Bar-Sever et al. (1998).

The ZTD and linear horizontal delay gradients were estimated every 5 min using a random walk model with a standard deviation (SD) of 10 and 0.3 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$, respectively. No weighting was applied. The SD value used for the ZTD is given by Jarlemark et al. (1998), where they found a temporal variability in the wet delay, derived from 71 d of microwave radiometer measurements, varying in the interval 3–22 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$. Because our focus here is on high temporal resolution and because we expect wet gradients to sometimes be short-lived, we also use a weaker constraint (in Sect. 3.4) of 1.0 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ for the SD in the random walk model for the gradients.

The zenith hydrostatic delay (ZHD) was calculated using ground pressure measurements (Saastamoinen, 1973), and thereafter the ZWD was obtained by subtracting the ZHD from the ZTD. The gradients estimated from the GNSS data are total gradients, and they were interpreted as the sum of hydrostatic and wet components.

The model used for the gradient estimation is presented by Bar-Sever et al. (1998):

where *S*(*ϵ*,*ϕ*) is a slant delay for a certain elevation and azimuth angle; *Z* and *m*(*ϵ*) are the zenith delay and the elevation
mapping function; and *G*_{n} and *G*_{e} are the north and east horizontal gradient, respectively.

In order to compare to the wet component inferred by the WVR, we subtracted the hydrostatic component computed from the reanalysis product of the 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 timescales of minutes to hours (Elgered et al., 2019). The gradients were calculated from ERA5 by vertical integration of the horizontal refractivity gradients times the height. The profile closest to the site was used together with one profile to the east and one profile to the north to calculate the refractivity gradient profiles. The temporal resolution of the ERA5 gradients is 1 h, and we therefore interpolated them to a 5 min resolution in order to perform the subtraction for the GNSS data. For the same reason, we did not use wet gradients from ERA5 because we want to study the gradients with a temporal resolution down to 5 min.

The data processing was run for three different elevation cutoff angles (3, 10, and 15^{∘}). For each elevation cutoff angle 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
Fig. 1 for the ONS1 station. The three systems show a similar geometry for the observations. At this latitude 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.

Figure 2 depicts the number of daily observations given for each GNSS constellation and obtained for each elevation cutoff
angle. Some days have fewer 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 is 10 % and 22 %,
respectively, while for Glonass, the corresponding values are 8 % and 17 %, respectively.

## 2.2 Water vapour radiometer (WVR)

The WVR, shown in Fig. 3, is located close to the GNSS stations, 9 m from ONSA and 59 m from ONS1, 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).

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 Fig. 4. Data were acquired throughout 2019, except from mid-August to early
October because of a failure caused by a thunderstorm.

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. The jumps were later found to be caused by vibrations when the mechanical waveguide switch was activated at the beginning of each 5 min period. When such a jump was identified, one or several complete 5 min cycles were removed. The jumps were identified by viewing the ZWD during each day. The temporal resolution is then sufficient to identify a 5 min long group of data that is discontinuous with the adjacent 5 min periods. 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 synchronized with the estimation period. Smaller gain jumps may still have degraded the accuracy of the estimated gradients. Thereafter, when the model was applied, we required that at least 40 of the 52 observations were available in the 5 min cycle period.

There are 105 120 possible 5 min data points in 1 year. After removing data acquired during rain, data that were unstable because of the jumps, and all 5 min periods where there were less than 40 (of the 52 scheduled) observations (typically caused by large liquid water contents), we ended up with 56 612 data points. There were 14 236 periods of 5 min during the time when the WVR was in the lab, meaning that the gradients estimated from 62 % of the time when the WVR was operated were used to compare to the ones from the GNSS. 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 study does not need to assume that the WVR gradients are more accurate compared to the GNSS ones. The main advantage of the WVR gradients is that they are independent, and by comparing these to the gradients from different GNSS solutions we can assess the different GNSS processing methods. Furthermore, since we want to study the agreement with as high temporal resolution as possible, we do not apply constraints to the individual 5 min gradients in order to have them independent from adjacent estimates in terms of the atmospheric signals.

The estimated ZWD from the WVR data is shown in Fig. 5. The seasonal dependence is clearly visible as well as a large short-term variability. Because of the gain jumps this data set is not optimum for ZWD comparisons on an absolute scale, but our focus is on gradients.

## 3.1 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 summarizes the statistics of estimated gradients and their corresponding formal errors. A few features of the data are worth noting. The mean of and the variability (standard deviation) in the estimated gradient amplitude increase with increasing elevation cutoff angle together with its mean formal error. When the elevation cutoff angle increases, a smaller number of observations were included for the gradient estimation. For the GNSS, the geometry of the satellite constellation is also deteriorated for a larger elevation cutoff angle. As a result, the formal error in the estimated gradient increases as well as the variability. In addition, when using a lower elevation, the larger volume sensed by the GNSS introduces an averaging effect that reduces the mean amplitude of the estimated gradients (see Elgered et al., 2019).

As indicated by column 9 in Table 1, the gradient amplitudes estimated by the WVR (0.99 mm) are approximately twice as
large as the GNSS gradient amplitudes at 3^{∘} cutoff angle, i.e. 0.49 mm, for the GRE solution, but they decrease to around 50 % as
large for the cutoff angle of 15^{∘}, i.e. 0.69 mm, for the GRE solution. The GNSS gradient amplitudes are about twice as large as their
formal errors. The WVR gradient amplitudes are about 8 times larger than their mean formal errors, but the variability in the WVR formal errors is
significantly larger than those from the 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 the formal errors in the GNSS and the WVR gradients when calculating the weighted
root-mean-square (WRMS) differences and correlations. It should also be pointed out that the uncertainty (here represented by the formal error) in the
WVR gradients is scaled, meaning that if the true wet delays in the different directions have deviations from the linear gradient model the
uncertainties increase. Such deviations are common, e.g. during convection processes, and the assumption of linear changes in the wet refractivity in
a layered atmosphere will not be accurate. The formal error in the gradient given by GipsyX is, however, not scaled. Therefore, these uncertainties
are likely smaller than realistic values.

## 3.2 Comparison of gradients from the GNSS and WVR

We first carry out comparisons of the gradients estimated from the different GNSS constellations and using the three different 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.

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 in 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 Fig. 1). We also note that the GRE solution, in general, gives the best agreement with the WVR gradients.

For the GPS-only solution, the highest 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 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 is presented in Table 1, the estimated size of the gradients was in general increasing when
using a decreasing sky coverage of the observations.

Even though the WVR and the GNSS sample different parts of the sky, it is noted that the agreement becomes worse 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 Fig. 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. The GNSS gradients are
obtained using different constellations and a cutoff angle of 3^{∘} (see Fig. 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 Fig. 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 July only. Another reason for focusing on these 2 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 data set (see
Table 2). The two GNSS stations (ONS1 and ONSA) show very similar agreements with the WVR gradients. This is partly expected
since they are located close to each other, and therefore the gradients from the two stations 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
Fig. 3) do not have a significantly different impact on the estimated gradients. This can also be seen in Fig. 7, which shows the mean code multipath RMS values for both ONS1 and ONSA calculated using Anubis software (Václavovic et al., 2016) and an elevation cutoff
angle of 5^{∘}. Therefore, we continue to use only ONS1 GNSS data in the following.

## 3.3 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 timescales and temporal
resolutions. We therefore averaged the GNSS and the WVR gradients by applying a Gaussian window with different FWHM from ± 2.5 to
± 720 min for further comparisons. As described in Sect. 3.2, when an FWHM of ± 2.5 min is used, only the
WVR gradients were interpolated to the epochs in the GNSS time series. For all FWHMs larger than ± 2.5 min, both the GNSS and the WVR
gradients are interpolated to the epochs at 0, 5, 10, …, and 55 min after the hour. The 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 refer to and use the FWHM as an effective temporal resolution, Δ*t*_{eff}, from 5 to 1440 min (1 d), although the
time series will still have a value every 5 min.

The resulting WRMS differences and correlations are shown in Fig. 8. It is clear that the WRMS differences decrease when the
Δ*t*_{eff} increases, and more variations in the gradients are averaged out. Over all GNSS solutions, the highest correlation for 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 9 and 10 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 1 d, 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.

We also study a specific event of short-lived gradients, associated with rapid changes in the ZWD, starting from 00:00 UTC on 23 July (see Fig. 11). 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 the GNSS and WVR, but the amplitudes of the gradients from the GNSS are much smaller. For the north direction, there are also some large gradients, i.e. around 08:00 and 17:00 UTC on 23 July, which are detected by both the GNSS and WVR. Also in this case the amplitudes from the GNSS are smaller. In addition, the multi-GNSS solutions have a slightly higher possibility of capturing 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 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.

## 3.4 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 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ (see
Sect. 2.1). This value may be too small, especially when we have more observations from multi-GNSS constellations, to 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 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$. The changes in WRMS differences and correlations,
relative to the solution using the constraint value of 0.3 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{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 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ are about twice as large compared to the ones obtained when using the constraint of
0.3 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$.

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 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 Fig. 12, depicting the time series of the gradients for the 2.5 d, starting at
00:00 UTC on 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 2.5 d, i.e. east gradients at 21:00 UTC on 24 July and north gradients at 08:00 and 17:00 UTC on 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 volume of the atmosphere is more different compared to that
observed by the WVR.

We have estimated linear horizontal gradients using 1 year of data acquired from the GNSS station 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).

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
(Fig. 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 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 and 6 h, respectively
(Fig. 8). When these Δ*t*_{eff} values are applied, strong gradients of short duration detected by the WVR but not by
the GNSS are averaged out, and as a result the correlation increases. When estimating GNSS gradients the choice of Δ*t*_{eff} is a compromise
between getting a high correlation and losing 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 (Figs. 9,
10, and 11).

Furthermore, weakening the constraint used when estimating the GNSS gradients from 0.3 to 1.0 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ helps the GNSS data to track short-lived gradients, approaching a timescale of 5 min, however at the cost of increased formal errors (Table 4 and Fig. 12).

Possible improvements to study in similar future work would be to include BeiDou observations and use a WVR with better stability. It would also be of interest to carry out a similar study at low-latitude GNSS stations, where the sky coverage is better, and perhaps also the atmosphere is more variable. In addition, the role of the geometry of GNSS observations (see Fig. 1) can be further studied. For example, one can remove observations in a certain direction and investigate the change in the estimated gradients and their formal uncertainties for different observation geometries.

The input GNSS data, in RINEX format, are available from EUREF (https://igs.bkg.bund.de/, Federal Agency for Cartography and Geodesy, 2021). The ERA5 data are accessible from https://doi.org/10.5065/BH6N-5N20 (ECMWF, 2019). The estimated gradients from GNSS and WVR data have been registered and archived by the Swedish National Data Service (SND): https://doi.org/10.5878/fyt8-bs80 (Elgered and Ning, 2021).

The two authors (TN and GE) planned the work and the structure of the paper together. TN performed the GNSS data analyses and GE performed the WVR data analyses. Both contributed to the writing of the manuscript and approved it before the submission.

The authors declare that they have no conflict of interest.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank Tobias Nilsson for providing the hydrostatic gradients from the ERA5 data.

This paper was edited by Roeland Van Malderen and reviewed by two anonymous referees.

Baldysz, Z., Nykiel, G., Figurski, M., and Araszkiewicz, A.: Assessment of the impact of GNSS processing strategies on the long-term parameters of 20 years IWV time series, Remote Sens.-Basel, 10, 496, https://doi.org/10.3390/rs10040496, 2018. a

Bar-Sever, Y.-E., Kroger, P. M., and Börjesson, J. A.: Estimating horizontal gradients of tropospheric path delay with a single GPS receiver, J. Geophys. Res., 103, 5019–5035, https://doi.org/10.1029/97jb03534, 1998. a, b, c

Boehm, J., Werl, B. and Schuh, H.: Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data, J. Geophys. Res., 111, B02406, https://doi.org/10.1029/2005JB003629, 2006. a

Davis, J. L., Elgered, G., Niell, A. E., and Kuehn, C. E.: Ground-based measurement of gradients in the “wet” radio refractivity of air, Radio Sci., 28, 1003–1018, https://doi.org/10.1029/93RS01917, 1993. a

Elgered, G. and Ning, T.: High temporal resolution wet delay gradients estimated from multi-GNSS and microwave radiometer observations, Chalmers University of Technology, Space, Earth and Environmental Science, Swedish National Data Service [data set], Version 1, https://doi.org/10.5878/fyt8-bs80, 2021. a

Elgered, G., Ning, T., Forkman, P., and Haas, R.: On the information content in linear horizontal delay gradients estimated from space geodesy observations, Atmos. Meas. Tech., 12, 3805–3823, https://doi.org/10.5194/amt-12-3805-2019, 2019. a, b, c, d, e

European Centre for Medium-Range Weather Forecasts (ECMWF): ERA5 Reanalysis (0.25 Degree Latitude-Longitude Grid), updated monthly, Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory [data set], https://doi.org/10.5065/BH6N-5N20, 2019. a

Federal Agency for Cartography and Geodesy, Germany: https://igs.bkg.bund.de/, last access: 12 August 2021. a

Graffigna, V., Hernández-Pajares, M., Gende, M. A., Azpilicueta, F. J., and Antico, P. L.: Interpretation ofthe tropospheric gradients estimated with GPS during Hurricane Harvey, Earth Space Science, 6, 1348–1365, https://doi.org/10.1029/2018EA000527, 2019. a

Jarlemark, P. O. J., Emardson, T. R., and Johansson, J. M.: Wet delay variability calculated from radiometric measurements and its role in space geodetic parameter estimation, Radio Sci., 33, 719–730. https://doi.org/10.1029/98RS00551, 1998. a

Kačmařík, M., Douša, J., Zus, F., Václavovic, P., Balidakis, K., Dick, G., and Wickert, J.: Sensitivity of GNSS tropospheric gradients to processing options, Ann. Geophys., 37, 429–446, https://doi.org/10.5194/angeo-37-429-2019, 2019. a

Li, X., Zus, F., Lu, C., Ning, T., Dick, G., Ge, M., Wickert, J., and Schuh, H.: Retrieving high-resolution tropospheric gradients from multiconstellation GNSS observations, Geophys. Res. Lett., 42, 4173–4181, https://doi.org/10.1002/2015GL063856, 2015. a

Lu, C., Li, X., Li, Z., Heinkelmann, R., Nilsson, T., Dick, G., Ge, M., and Schuh, H.: GNSS tropospheric gradients with high temporal resolution and their effect on precise positioning, J. Geophys. Res.-Atmos., 121, 912–930, https://doi.org/10.1002/2015JD024255, 2016. a

Lyard, F., Lefevre, F., Letellier, T., and Francis, O.: Modelling the global ocean tides: Modern insights from FES2004, Ocean Dynam., 56, 394, https://doi.org/10.1007/s10236-006-0086-x, 2006. a

Meindl, M., Schaer, S., Hugentobler, U., and Beutler, G.: Tropospheric gradient estimation at CODE: results from global solutions, J. Meteorol. Soc. Jpn., 82, 331–338, https://doi.org/10.2151/jmsj.2004.331, 2004. a

Morel, L., Pottiaux, E., Durand, F., Fund, F., Boniface, K., de Oliveira, P. S., and Van Baelen, J.: Validity and behaviour of tropospheric gradients estimated by GPS in Corsica, Adv. Space Res., 55, 135–149, https://doi.org/10.1016/j.asr.2014.10.004, 2015. a

Ning, T., Haas, R., Elgered, G., and Willén, U.: Multi-technique comparisons of ten years of wet delay estimates on the west coast of Sweden, J. Geodesy, 86, 565–575, https://doi.org/10.1007/s00190-011-0527-2, 2012. a

Prange, L., Villiger, A., Sidorov, D., Schare, S., Beutler, G., Dach, R., and Jäggi, A.: Overview of CODE's MGEX solution with the focus on Galileo, Adv. Space Res., 66, 2786–2798, https://doi.org/10.1016/j.asr.2020.04.038, 2020. a

Saastamoinen, J.: Contributions to the theory of atmospheric refraction, B. Geod., 107, 13–34, https://doi.org/10.1007/BF02521844, 1973. a

Schmid, R., Steigenberger, P., Gendt, G., Ge, M., and Rothacher, M.: Generation of a consistent absolute phase center correction model for GPS receiver and satellite antennas, J. Geodesy, 81, 781–798, https://doi.org/10.1007/s00190-007-0148-y, 2007. a

Steigenberger, P., Tesmer, V., Krügel, M., Thaller, D., Schmid, R., Vey, S., and Rothacher, M.: Comparison of homogeneously reprocessed GPS and VLBI long time-series of tropospheric zenith delays and gradients, J. Geodesy, 81, 503–514, https://doi.org/10.1007/s00190-006-0124-y, 2007. a

Václavovic, P. and Douša, J.: G-Nut/Anubis: open-source tool for multi-GNSS data monitoring, in: IAG Symposia Series, Springer, 143, 775–782, https://doi.org/10.1007/1345_2015_97, 2016. a

Zhou F., Li, X., Li, W., Chen, W., Dong, D., Wickert, J., and Schuh, H.: The Impact of Estimating High-Resolution Tropospheric Gradients on Multi-GNSS Precise Positioning, Sensors, 17, 756, https://doi.org/10.3390/s17040756, 2017. a

Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M., and Webb, F. H.: Precise Point Positioning for the Efficient and Robust Analysis of GPS Data from Large Networks, J. Geophys. Res., 102, 5005–5017, https://doi.org/10.1029/96JB03860, 1997. a

Zus, F., Douša, J., Kačmařík, M., Václavovic, P., Dick, G. and Wickert, J.: Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation, Remote Sens.-Basel, 11, 41, https://doi.org/10.3390/rs11010041, 2019. a