Improving GNSS Zenith Wet Delay Interpolation by Utilizing Tropospheric Gradients: Experiments with a Dense Station Network in Central Europe in the Warm Season

Abstract

The Benchmark data set collected within the European COST Action ES1206 has aimed to support the development and validation of advanced Global Navigation Satellite System (GNSS) tropospheric products, in particular high-resolution zenith delays and tropospheric gradients. In this work we utilize this unique data set to show that the interpolation of GNSS Zenith Wet Delays (ZWDs) can be improved by utilizing tropospheric gradients. To do this we first prove the concept with simulated observations, that is, zenith delays and tropospheric gradients derived from a Numerical Weather Model. We show how tropospheric gradients can be converted to ZWD gradients. Then the ZWD gradients together with the ZWDs at selected reference stations are used in an inverse distance weighting interpolation scheme to estimate the ZWD at some target station. For a station configuration with an average station distance of 50 km in Germany and a period of two months (May and June 2013), we find an improvement of 20% in interpolated ZWDs when tropospheric gradients are taken into account. Next, we replace the simulated by real observations, that is, zenith delays and tropospheric gradients from a Precise Point Positioning (PPP) solution provided with the G-Nut/Tefnut analysis software. Here we find an improvement of 10% in interpolated ZWDs when tropospheric gradients are taken into account.

1. Introduction

Carrier phase (and code) observations from a ground-based Global Navigation Satellite System (GNSS) station allow the estimation of the Zenith Total Delay (ZTD) [1] and the tropospheric gradient [2]. From the ZTD, the Zenith Wet Delay (ZWD) can be retrieved. The tropospheric gradient at the station can be roughly related to the horizontal ZWD gradient at the station, that is, the partial derivative of the ZWD with respect to the longitude (latitude). In other words, a single station provides information on the ZWD and horizontal ZWD gradient. Thus it is natural to try to improve ZWD interpolation by utilizing tropospheric gradients. The ZWD can be converted to the Integrated Water Vapor (IWV), which is closely linked to precipitation. An improved ZWD map implies an improved IWV map, which is valuable in meteorological studies [3,4].

A prerequisite for an accurate interpolation is the quality of ZTDs and tropospheric gradients. To date most attention has been paid to ZTDs (see, e.g., [5] and references therein). Little attention has been paid to tropospheric gradients. Bar-Sever et al. presented in [2] preliminary results providing some evidence that the tropospheric gradients contain real tropospheric features. To what extent the tropospheric gradients contain useful tropospheric information was not exploited further by any means. A number of comparison studies followed. For example, Li et al. [6] showed that an improved observation geometry yields improved tropospheric gradient estimates. Morel et al. used in [7] different software packages and compared tropospheric gradients from a dozen of stations located on an island in the Mediterranean Sea (Corsica). Typically, such studies are restricted to a few stations. Dousa et al. presented in [8] a comparison for hundreds of stations including two GNSS estimation methods (double- and zero-differencing) and compared them with tropospheric gradients derived from Numerical Weather Models (NWMs). The comparison of tropospheric gradient maps, that is, the visual inspection, provided clear evidence that the GNSS tropospheric gradients include real tropospheric features. Although, the results were encouraging some points needed further investigation. For example, due to the low horizontal resolution of the underlying NWMs, the magnitude of NWM tropospheric gradients was underestimated with respect to the magnitude of the GNSS tropospheric gradients. In addition, the role of the GNSS data processing options and parameters remained unclear. For example, the impact of higher-order ionospheric corrections on estimated tropospheric gradients was analyzed in [9]. Clear artificial signals in tropospheric gradients were observed during the 2nd EUREF reprocessing [10] and identified as an absorption of asymmetric effect caused by problems related to low-elevation observations provided by progressively degrading instrumentation. Kacmarik et al. studied in [11] the sensitivity of tropospheric gradients to several processing options. Generally, all the solutions run in the post-processing mode provided a robust tropospheric gradient estimation with a clear relation to real weather conditions. The quality of tropospheric gradient estimates in real-time mode mainly depended on the actual quality of the real-time orbits and clocks. In addition, [11] showed that the higher horizontal resolution of the NWMs yields a better match between NWM and GNSS tropospheric gradient maps. As a matter of example, Figure 1 shows the GNSS and NWM tropospheric gradient map for an epoch where strong tropospheric gradients are present (10 May 2013, 6 UTC). Details on the GNSS and NWM tropospheric parameters are provided in Section 2 below. In order to show that tropospheric gradients are often accompanied by severe weather, we overlaid the tropospheric gradient maps with the radar image (instantaneous rain) provided by the Deutscher Wetter Dienst (DWD). We also added the station specific GNSS IWV values in order to show that the tropospheric gradients can be roughly related to horizontal IWV gradients; typically, the tropospheric gradients point from dry to moist areas. This was also recognized by [12] studying preliminary signs of deep convection. In fact, they used ZTDs and tropospheric gradients to derive so-called pseudo-ZTDs in order to improve IWV visualization. Similarly, Douša et al. [13] used the pseudo-ZTDs for improving a new two-stage tropospheric correction model. If precise tropospheric gradients are estimated along with ZTDs, the corresponding pseudo-ZTDs could be used to improve spatial interpolation.

The purpose of this study is to show that ZWD interpolation can be improved by utilizing tropospheric gradients. The algorithm developed in this study can be regarded an improved tropospheric correction model. Such a tropospheric correction model can be used for providing tropospheric augmentation corrections in positioning applications (see, e.g., [13] and references therein). The interpolated ZWD can be converted to IWV, which makes the tropospheric correction model valuable in meteorological applications as well.

The structure of this work is as follows. In Section 2 we introduce ZTDs and tropospheric gradients and describe the data sets used in this study. In Section 3 we show how tropospheric gradients can be converted to horizontal ZWD gradients and we describe the interpolation method. In Section 4 we present the results. At first, we make use of simulated observations before we make use of real observations. The purpose of utilizing simulated observations is to prove the concept in advance. In Section 5 we discuss the results. The conclusion is given in Section 6.

2. ZTDs and Tropospheric Gradients

In the GNSS analysis, the tropospheric delay is parameterized utilizing Mapping Functions (MF). Specifically, for the elevation angle e and azimuth angle a, the parameterized tropospheric delay T reads as

$$T\left(e,a\right)=m_h\left(e\right)\cdot ZHD+m_w\left(e\right)\cdot ZWD+m_g\left(e\right)\left[cos\left(a\right)\cdot N+sin\left(a\right)\cdot E\right]$$

(1)

where ZHD denotes the Zenith Hydrostatic Delay, ZWD denotes the Zenith Wet Delay, N denotes the north-gradient component, E denotes the east-gradient component, m_h denotes the hydrostatic MF, m_w denotes the wet MF, and m_g denotes the gradient MF [14]. The ZTD is given by

$$ZTD=ZHD+ZWD$$

(2)

Similarly to the ZTD, the tropospheric gradient can be written as the sum of a hydrostatic and wet contribution. In essence,

$$\begin{gathered} E=E_h+E_w \\ N=N_h+N_w \end{gathered}$$

(3)

where N_h (N_w) denotes the hydrostatic (wet) north-gradient component and E_h (E_w) denotes the hydrostatic (wet) east-gradient component. It is important to note that the GNSS analysis (see below) provides accurate ZTDs and (total) tropospheric gradients. The GNSS analysis does not provide ZWDs and wet gradients. In this study, the GNSS ZWDs and wet gradients are obtained from the GNSS ZTDs and tropospheric gradients by subtracting NWM ZHDs and hydrostatic gradients.

2.1. Benchmark Data Set

We base our study on the Benchmark data set which was collected within the European COST Action ES1206 GNSS4SWEC (Advanced GNSS tropospheric products for monitoring severe weather and climate). Among a broad set of various meteorological observations and products, it contains data from 430 GNSS reference stations in the area of central Europe covering Germany, the Czech Republic, and parts of Poland and Austria. The average station distance is about 50 km. The station locations are shown in Figure 1. All the data and some GNSS and NWM reference tropospheric products are available for the period of May and June 2013, when some significant rainfall episodes of various origin leading to extensive floods occurred in the selected area. For comprehensive and detailed information about the Benchmark campaign, the reader is referred to [8].

2.2. GNSS ZTDs and Tropospheric Gradients

We used the G-Nut/Tefnut software [15] for processing GNSS observations from GPS and GLONASS constellations. The software employed the Precise Point Positioning (PPP) method [16] and we utilized precise products for satellite orbits and clocks provided by the European Space Agency, ESA (http://navigation-office.esa.int/GNSS_based_products.html). We used pseudorange and carrier-phase observations with the 3° elevation angle cut-off and we applied an observation weighting based on the function 1/sin²(e). The estimated site-specific state vector comprised receiver coordinates, receiver clock corrections, ZTD, tropospheric gradient, float initial ambiguities of carrier-phase ionosphere-free linear combination, and Inter System Bias (ISB) for synchronizing different time frames of GPS and GLONASS. Extended Kalman filter supported by the backward smoother was applied for the state vector estimation [17]. We modelled GNSS observations according to the IERS 2010 conventions [18], and we used the file igs14_2013.ATX for deriving antenna phase center offsets and variations [19], which are consistent with the IGS2014 reference frame [20] used for precise satellite positions.

The station coordinates were estimated along with other parameters, considered constant over a day, that is, tightly constrained. The sigma of a prior coordinates were 100 m. The ambiguities were estimated as float values and thus not resolved to integer values. The station clock corrections were estimated without any dependency between consecutive epochs. We estimated the tropospheric parameters at a 5-min sampling rate as a random walk stochastic process with a noise of 6 mm/√h and 0.6 mm/√h for ZTD and tropospheric gradients, respectively. The ZHD, considered as the a priori value of ZTD, was calculated from the Global Pressure and Temperature model, GPT [21] and introduced into the GNSS adjustment model with an a priori sigma of 10 cm. Both hydrostatic and wet parts of ZTD were mapped from the zenith to particular satellite direction using hydrostatic and wet factors derived from the Global Mapping Function, GMF [22]. Tropospheric gradients were initialized by zero values, and Chen and Herring gradient MF [14] was employed for deriving slant delay contribution.

2.3. NWM ZTDs and Tropospheric Gradients

We utilized the Weather Research and Forecasting (WRF) model [23] to simulate the refractivity field of the atmosphere. The initial and boundary conditions for the limited area model were a Global Forecast System (GFS) analysis of the National Centers for Environmental Prediction (NCEP). The 24-h free forecasts start every day at 0 UTC. The following physics was applied: the Thompson scheme [24] for the microphysics, the Kain–Fritsch scheme [25] for the cumulus parameterization, the Yonsei University scheme [26] for the planetary boundary layer, the RRTMG Short and Longwave scheme [27] for the radiation, the Unified Noah Land Surface Model [28] for the land surface, and the revised MM5 scheme [29] for the surface layer. The refractivity, which was calculated from the pressure, temperature, and humidity [30], was available every hour with a horizontal resolution of 10 km on 50 vertical model levels (up to 50 hPa). The refractivity at an arbitrary point was obtained by interpolation [31].

For a given station location, the tropospheric gradients were derived from ray-traced tropospheric delays [32]. Firstly, a set of tropospheric delays was computed (it consists of 120 tropospheric delays where the elevation angles are 3°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, and 90° and the spacing in azimuth is 30°). Secondly, we computed azimuth independent tropospheric delays under the assumption of a spherically layered troposphere. Finally, the differences between the azimuth dependent tropospheric delays and the azimuth independent tropospheric delays were computed and then the gradient components were determined by a least-square fit [33]. The Chen and Herring gradient MF [14] was used. The hydrostatic delays were computed separately in order to calculate hydrostatic gradients. The wet gradients were obtained from the tropospheric gradients by subtracting the hydrostatic gradients. The ZTD was the tropospheric delay in the zenith direction. The ZHD was the hydrostatic delay in the zenith direction. The ZWD was obtained by subtracting the ZHD from the ZTD.

We estimated the accuracy of the NWM ZHD and hydrostatic gradient by comparing with the ZHD and hydrostatic gradient derived from the atmospheric reanalysis ERA5 [34]. The ERA5 is produced at the European Centre for Medium Range Weather Forecasts (ECMWF). The refractivity fields are available hourly with a horizontal resolution of about 31 km (T639 spectral triangular truncation). We considered the hourly data from the full benchmark period. On average (mean over all stations) the Root Mean Square Deviation (RMSD) was 3 mm for the ZHD and below 0.1 mm for both hydrostatic gradient components. For comparison, the RMSD was 12 mm for the ZWD and about 0.4 mm for both wet gradient components. In essence, the WRF model allows an accurate (short range) prediction of the ZHD and the hydrostatic gradient components. Hence, as far as one is concerned with the ZHD and the hydrostatic gradient components, the WRF model can be potentially used in real time applications as well.

3. Method

3.1. Convert Tropospheric Gradients to Horizontal ZWD Gradients

The following relations are approximately true, see, for example Reference [35]

$$\begin{gathered} E_w\sim C\frac{\partial ZWD}{\partial x}=\frac{C}{R}\frac{\partial ZWD}{\partial \lambda }\frac{1}{cos\varphi } \\ {N}_w\sim C\frac{\partial ZWD}{\partial y}=\frac{C}{R}\frac{\partial ZWD}{\partial \varphi } \end{gathered}$$

(4)

Here λ and ϕ denote the station longitude and latitude, R denotes Earth’s radius, and the factor C is related to the scale height of the wet refractivity gradient. Hence, if we know the wet gradient at the station then we know approximately the horizontal ZWD gradient at that station and vice versa [36]. In order to demonstrate this, we chose a station and selected neighboring stations within a radius of 1°. We collected the ZWDs from these stations, assumed that the ZWD can be expanded in a Taylor-series and obtained by a least square fit the horizontal ZWD gradient at the chosen station. This horizontal ZWD gradient at the station was roughly proportional to the wet gradient at the station. We repeated this procedure for each station in the benchmark data set. For some stations the complex terrain may play a role, however, so we left this investigation for future research. We chose an epoch where strong tropospheric gradients were present (31 May 2013, 18 UTC). We experimented with the factor C and found that an appropriate choice was 4 km. We utilized NWM data (similar results are obtained if we utilize GNSS data). The left panels of Figure 2 show the wet east and north-gradient component and the right panels of Figure 2 show the approximation for the wet east and north-gradient component. The gradient components do not match perfectly but the validity of the approximation is obvious. To some extent the differences between the left and right panel are due to the fact that C is chosen to be 4 km for any location. Nevertheless, it is clear that we have a tool with which we can at least approximately convert tropospheric gradients to horizontal ZWD gradients.

3.2. Interpolation Algorithm

The ZWD interpolation algorithm consists of two steps. At first, the ZWD and the wet gradient components at the reference stations are computed according to

$$\begin{gathered} ZWD=ZTD-ZHD \\ {E}_w=E-E_h \\ {N}_w=N-N_h \end{gathered}$$

(5)

Here the ZHD and the hydrostatic gradient components are derived from the NWM. Second, the ZWD at some target station (location), denoted $ZW{D}_t$, is obtained through an inverse distance weighting interpolation technique

$$ZW{D}_t=\sum w^i\cdot \left(ZW{D}^i+\frac{\partial ZW{D}^i}{\partial \lambda }\Delta {\lambda }^i+\frac{\partial ZW{D}^i}{\partial \varphi }\Delta {\varphi }^i\right)\cdot exp\left(-\frac{\Delta h^i}{H}\right)$$

(6)

where w denotes the weights of the inverse distance weighting scheme, ∆λ denotes the longitude difference, ∆ϕ denotes the latitude difference, and ∆h denotes the height difference between the target and the reference station respectively. The index i indicates the ith reference station. Following Equation (4), the partial derivative of the ZWD with respect to the longitude and latitude can be approximated utilizing the wet gradient components, which leads to

$$ZW{D}_t=\sum w^i\cdot \left(ZW{D}^i+\frac{E_w^i}{C}Rcos\left({\varphi }^i\right)\Delta {\lambda }^i+\frac{N_w^i}{C}R\Delta {\varphi }^i\right)\cdot exp\left(-\frac{\Delta h^i}{H}\right)$$

(7)

For any target station we consider reference stations within a radius of 1°. The factor C is chosen to be 4 km. The scale height of the ZWD, denoted H, is chosen to be 3 km. A more sophisticated vertical adjustment utilizing location and time dependent scale heights (see, e.g., [37]) remains to be implemented. We note that the interpolation algorithm above can be written as

$$ZW{D}_t=\sum w^i\cdot ZW{D}^i\cdot exp\left(-\frac{\Delta h^i}{H}\right)+\sum w^i\cdot \left(\frac{E_w^i}{C}Rcos\left({\varphi }^i\right)\Delta {\lambda }^i+\frac{N_w^i}{C}R\Delta {\varphi }^i\right)\cdot exp\left(-\frac{\Delta h^i}{H}\right)$$

(8)

where the first sum can be regarded as the standard interpolation algorithm and the second sum can be regarded as the correction to it.

Some end users are interested in the IWV and not the ZWD. The IWV at the target location, denoted $IW{V}_t$, is obtained through

$$IW{V}_t=\kappa \left(T_m\right)·ZW{D}_t$$

where the conversion factor κ depends on the weighted mean temperature Tm at the target location. The weighted mean temperature is related to the surface temperature [1]. We made use of the NWM to obtain the surface temperature at the target location.

We note that other more sophisticated interpolation schemes (for example kriging) are possible. However, we leave such investigations for future studies. The advantage of the chosen interpolation scheme is the ease of implementation.

4. Results

Each station in the benchmark data set was treated as a target station. The ZWD at the target station was regarded the true ZWD. Data from the surrounding stations, the reference stations, were used to estimate ZWD at the target station. This estimated ZWD at the target station was then compared with the true ZWD at the target station. It was thus straightforward to estimate the error of the interpolation strategy. We compared two interpolation strategies: (1) only ZWDs from the reference stations are used and (2) ZWDs and tropospheric gradients from the reference stations are used. Hereinafter the two methods are referred to as the first and second interpolation strategy respectively. In the first subsection we make use of simulated observations, that is, NWM ZTDs and tropospheric gradients. In the second subsection we replace the simulated observations by real observations, that is, GNSS ZTDs and tropospheric gradients.

4.1. Experiment with NWM Data

Differences between the two interpolation strategies can be expected for large tropospheric gradients. Therefore we first consider a single epoch where large tropospheric gradients are present. The scatter plot in Figure 3 shows the station specific ZWD differences between the interpolation strategies, that is, the correction to the standard interpolation algorithm, on 31 May 2013 at 18 UTC. The ZWD differences reach up to 25 mm. The ZWD differences show a specific pattern. A quick look at Figure 2 explains why ZWD differences are larger in some areas than in other areas. The tropospheric gradients show a specific pattern and the ZWD differences follow this pattern. Figure 3 does not provide information on how well the interpolation strategies perform. In Figure 4 the upper (lower) panel shows the ZWD error of the first (second) interpolation strategy. On average (mean over all stations) the Root Mean Square Error (RMSE) is reduced by more than 30% when tropospheric gradients are taken into account. Next, instead of a single epoch, we consider the hourly data from the full benchmark period of two months. In Figure 5, the upper (lower) panel shows the station specific ZWD RMSE of the first (second) interpolation strategy. On average the RMSE is reduced by about 20% when tropospheric gradients are taken into account.

4.2. Experiment with GNSS Data

Next we used GNSS ZTDs and tropospheric gradients. The scatter plot in Figure 6 shows the station specific ZWD differences between the interpolation strategies on 31 May 2013 at 18 UTC. Again, the ZWD differences reach up to 25 mm. The pattern in Figure 6 closely follows the pattern in Figure 3 due to the fact that the GNSS tropospheric gradients closely follow the NWM tropospheric gradients. In Figure 7 the upper (lower) panel shows the interpolation error of the first (second) strategy. On average (mean over all stations) the RMSE is reduced by about 20% when tropospheric gradients are taken into account. Figure 7 can be compared with Figure 4 provided above. The improvement when utilizing real observations is smaller compared to the improvement when utilizing simulated observations. Again, instead of a single epoch, we consider the hourly data from the benchmark period of two months. In Figure 8 the upper (lower) panel shows the station specific ZWD RMSE of the first (second) interpolation strategy. On average the RMSE is reduced by about 10% when tropospheric gradients are taken into account. Figure 8 can be compared to Figure 5 provided above. Clearly, the RMSE is reduced when tropospheric gradients are utilized, however, the improvement when utilizing real observations is smaller compared to the improvement when utilizing simulated observations.

The inclusion of the tropospheric gradients in the interpolation only makes a difference during severe weather. The tropospheric gradients are smaller during calm weather which is somewhat true for most of the days. As a result, the average impact is naturally small. It is thus recommended to take a closer look on specific epochs. As another case study we choose 10 May 2013, 6 UTC, the epoch used for illustration in the introduction (refer to Figure 1). However, instead of the interpolated ZWD we look at the interpolated IWV. The scatter plot in Figure 9 shows the station specific IWV differences between the interpolation strategies on 10 May 2013 at 6 UTC. The tropospheric gradients shown in Figure 1 explain why IWV differences are larger in some areas than in other areas. In south-west Germany, IWV differences reach up to 3 kg/m². In Figure 10 the upper (lower) panel shows the IWV error of the first (second) interpolation strategy. On average (mean over all stations) the Root Mean Square Error (RMSE) is reduced by about 10% when tropospheric gradients are taken into account.

5. Discussion

In the previous section we have shown that the ZWD interpolation can be improved by utilizing tropospheric gradients. The improvement when utilizing real observations is smaller compared to the improvement when utilizing simulated observations. The reason for that is twofold. First, the ZWD variability in reality is larger than the ZWD variability provided by a NWM with a horizontal resolution of 10 km. Thus, the interpolation errors in the real environment are larger than in an artificial environment. Second, the GNSS ZTDs and tropospheric gradients are not error free.

It is important to note that we have chosen final orbits and clocks in PPP. This means that we have shown that the ZWD interpolation can be improved by utilizing tropospheric gradients in post processing mode. For example, Shi et al. study in [38] the impact of Real Time (RT) satellite clock errors on the estimated ZTDs. In order to check if the ZWD interpolation can be improved by utilizing tropospheric gradients in RT mode we make use of RT orbits and clocks instead of final orbits and clocks. Specifically, we choose the IGS01 product from the IGS Real-Time Service (http://rts.igs.org). We also note that in RT mode we cannot apply backward smoothing. However, our experience is that the quality of RT orbits and clocks is more important than the backward smoothing for the tropospheric gradient estimation [11]. In essence, in RT mode we find that on average the RMSE is reduced by about 5% when tropospheric gradients are taken into account. We explain this reduction in the RMSE of (only) 5% instead of 10% by the fact that the GNSS ZTDs and tropospheric gradients are less accurate in RT than in post processing mode.

6. Conclusions

We have shown that the ZWD interpolation can be improved by utilizing tropospheric gradients. To do so we first showed how tropospheric gradients can be converted to horizontal ZWD gradients. We note that for this conversion ZHDs and hydrostatic gradients from an external source are required. We used NWM data for this purpose. Then the horizontal ZWD gradients together with the ZWDs at some reference stations were used to estimate the ZWDs at arbitrary locations. The interpolation scheme, which is essentially an inverse distance weighting scheme, remains to be refined in the future. In particular, an improved vertical adjustment utilizing location and time dependent scale heights should be considered.

At first, we tested the interpolation strategy with simulated observations, i.e., NWM zenith delays and tropospheric gradients. For a station configuration with an average station distance of about 50 km in Germany and a period of two months in the summer 2013, we found an improvement of 20% in the interpolated ZWD when tropospheric gradients were taken into account. Then, we replaced the simulated observations by real observations, i.e., zenith delays and tropospheric gradients from a PPP solution by the G-Nut/Tefnut software. We found an improvement of 10% in the interpolated ZWD when tropospheric gradients are taken into account. We found this improvement in post processing mode. In real time mode, we found an improvement of 5% in the interpolated ZWD when tropospheric gradients were taken into account. The fact that the ZWD interpolation is improved by utilizing tropospheric gradients can be regarded as another piece of evidence that the tropospheric gradients contain useful tropospheric information.

Though the results of this study are encouraging, it is important to keep in mind that they are based on a somewhat limited data set, that is, the Benchmark data set collected within the European COST action ES1206. Hence, the results are representative for a dense station network in central Europe in the warm season. Future studies should investigate other station networks in different climates.