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

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.

1 Introduction

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 Data sets

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}{\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}{\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}{\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):

$$\begin{array}{}\text{(1)}& S(\mathit{ϵ},\mathit{\varphi })=m\left(\mathit{ϵ}\right)\left[Z+\text{cot}\left(\mathit{ϵ}\right)\left(G_{\mathrm{n}}\cos \left(\mathit{\varphi }\right)+G_{\mathrm{e}}\sin \left(\mathit{\varphi }\right)\right)\right],\end{array}$$

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.

Figure 1 Observations of the sky acquired from ONS1 for 0 to 24 h on 24 July 2019 for GPS (a), Glonass (b), and Galileo (c).

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.

Figure 2 Number of daily observations for each GNSS constellation acquired by ONS1 applying three different elevation cutoff angles (3, 10, and 15^∘).

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.

Figure 3 The water vapour radiometer (WVR) and GNSS stations (ONS1 and ONSA) at the Onsala Space Observatory.

Figure 4 One cycle of the WVR observations consists of 52 observations. The cycle is repeated every 5 min. Note that the zenith point is only measured once, the first time it passes over zenith.

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.

Figure 5 The time series of the ZWD estimated from the WVR data for each 5 min cycle. The radiometer was repaired from mid-August to early October.

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 Results

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

Table 1 The mean and standard deviations (SDs) of the east and the north gradients and the gradient amplitude, together with the mean and SD of their 1σ formal errors.

^a GPS + Galileo, ^b GPS + Glonass, ^c GPS + Glonass + Galileo.

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 The WRMS differences and correlations of the east and north gradients, obtained from different satellite constellations, relative to the WVR data.

^a GPS + Galileo, ^b GPS + Glonass, ^cGPS + Glonass + Galileo.

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.

Figure 6 Correlations between estimated gradients from the GNSS and the WVR data calculated for each month.

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.

Table 3 The WRMS differences and correlations of the east and the north gradients, obtained from different satellite constellations, relative to the WVR data for June and July 2019.

^a GPS + Galileo, ^b GPS + Glonass, ^c GPS + Glonass + Galileo.

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.

Figure 7 The mean code multipath RMS calculated from ONS1 (a) and ONSA (b).

Figure 8 Correlations and WRMS differences between estimated gradients from the GNSS and the WVR data using different effective temporal resolutions, Δt_eff.

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.

Figure 9 Correlations between estimated east gradients from the GNSS, given by the GRE solution with a 3^∘ elevation cutoff angle, and the WVR data. Four different values of Δt_eff are used: (a) 5 min, (b) 2 h, (c) 12 h, and (d) 24 h.

Figure 10 Correlations between estimated north gradients from the GNSS, given by the GRE solution with a 3^∘ elevation cutoff angle, and the WVR data. Four different values of Δt_eff are used: (a) 5 min, (b) 2 h, (c) 12 h, and (d) 24 h.

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.

Figure 11 ZWD (top) and gradient time series from the WVR and the GNSS solution using an elevation cutoff angle of 3^∘. The east and north gradients are shown in the second, third, and fourth row of graphs for Δt_eff of 5 min, 2 h, and 24 h, respectively. Note the different scales for the east and the north gradient graphs.

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.

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 $\mathrm{mm}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ relative to 0.3 $\mathrm{mm}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ (in Table 3).

^a GPS + Galileo, ^b GPS + Glonass, ^c GPS + Glonass + Galileo.

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.

Figure 12 Gradient time series from the WVR and from the GNSS (GRE) applying two different constraint values (0.3 and 1.0 $\mathrm{mm}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$) and using elevation cutoff angles of 3^∘ (a, b) and 10^∘ (c, d).

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}{\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}{\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}{\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}{\sqrt{\mathrm{h}}}^{-\mathrm{1}}$ are about twice as large compared to the ones obtained when using the constraint of 0.3 $\mathrm{mm}{\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.

4 Conclusions

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}{\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.