Tropospheric refractivity and zenith path delays from least-squares collocation of meteorological and GNSS data

Abstract

Precise positioning requires an accurate a priori troposphere model to enhance the solution quality. Several empirical models are available, but they may not properly characterize the state of troposphere, especially in severe weather conditions. Another possible solution is to use regional troposphere models based on real-time or near-real time measurements. In this study, we present the total refractivity and zenith total delay (ZTD) models based on a numerical weather prediction (NWP) model, Global Navigation Satellite System (GNSS) data and ground-based meteorological observations. We reconstruct the total refractivity profiles over the western part of Switzerland and the total refractivity profiles as well as ZTDs over Poland using the least-squares collocation software COMEDIE (Collocation of Meteorological Data for Interpretation and Estimation of Tropospheric Pathdelays) developed at ETH Zürich. In these two case studies, profiles of the total refractivity and ZTDs are calculated from different data sets. For Switzerland, the data set with the best agreement with the reference radiosonde (RS) measurements is the combination of ground-based meteorological observations and GNSS ZTDs. Introducing the horizontal gradients does not improve the vertical interpolation, and results in slightly larger biases and standard deviations. For Poland, the data set based on meteorological parameters from the NWP Weather Research and Forecasting (WRF) model and from a combination of the NWP model and GNSS ZTDs shows the best agreement with the reference RS data. In terms of ZTD, the combined NWP-GNSS observations and GNSS-only data set exhibit the best accuracy with an average bias (from all stations) of 3.7 mm and average standard deviations of 17.0 mm w.r.t. the reference GNSS stations.

1 Introduction

The Global Navigation Satellite System (GNSS) signal propagating through the atmosphere is delayed due to the free electron content in the ionosphere and by the air density in the electrically neutral atmosphere. Both influences can be described by the refractive index n or the total refractivity $N_{\mathrm{tot}}$:

$$\begin{aligned} N_{\mathrm{tot}}=10^{6}(n-1) \text { }\left[ \frac{\text {mm}}{\text {km}}=\text {ppm}\right] . \end{aligned}$$

(1)

The refractivity of the troposphere is measured in GNSS meteorology by zenith total delay (ZTD) and troposphere gradients in north and east directions or as a function of meteorological parameters. Conversely, the zenith path delay can be calculated based on the refractivity values.

The tropospheric delay empirical models are usually functions of meteorological parameters (temperature, pressure and humidity). The application of standard atmosphere parameters or global models, such as the GPT (global pressure/temperature) model (Böhm et al. 2007) or the UNB3 (University of New Brunswick, version 3) model (Leandro et al. 2006), may not be sufficient, especially for positioning in non-standard weather conditions. In this study, we present a troposphere model that utilizes a collocation technique to reconstruct the troposphere conditions based on GNSS and meteorological observations. The goal is to obtain the model of troposphere parameters (i.e., total refractivity and ZTD), which can be used as an a priori model of troposphere or to constrain tropospheric estimates in positioning. This model can be used in various applications (not only in GNSS processing but also others, e.g., InSAR), but is mainly designed to provide troposphere estimates for Real-Time Kinematic Precise Point Positioning (RTK-PPP). For PPP, especially when processing in kinematic mode, the accuracy of estimated positions depends heavily on the applied a priori model (Jensen and Ovstedal 2008; Wielgosz et al. 2011). The main drawback of the PPP method is that a long interval of about 20-30 minutes is required for the solution convergence (Li et al. 2011; Dousa and Vaclavovic 2014). One of the reasons is the high correlation among the estimated parameters: troposphere delay, receiver clock offset and receiver height. A possible solution to efficiently de-correlate these parameters and shorten the convergence time is to introduce the external high-quality regional troposphere delay model to constrain the troposphere estimates (Hadaś 2015; Shi et al. 2014).

In previous investigations, several troposphere models have been incorporated into PPP software. Hadaś et al. (2013) have chosen two regional models: one based on GNSS data and one based on ground-based meteorological data to be applied into GNSS-WARP (Wroclaw Algorithms for Real-Time Positioning) software. The ZTD model based on the GNSS data exhibits −6.2 mm bias with a standard deviation of 8.8 mm w.r.t. the control solutions from the International GNSS Service (IGS) processing center—Military University of Technology Analysis Centre in Warsaw (MUT). The model from ground-based meteorological stations has various ZTD shifts for each station (from −100 to 10 mm). Application of the GNSS-based model improves the RMS by 33 % for all position components compared to positioning without any a priori model. The model based on GNSS data introduces the bias of 1 ± 7 cm RMS for the North component, 2 ± 10 cm for the East component and −5 ± 12 cm for the Up component w.r.t. the MUT final control solution. Jensen and Ovstedal (2008) have exploited a troposphere model based on NWP model High Resolution Limited Area Model (HIRLAM) along with standard atmosphere Saastamoinen and UNB3 models. The results for North and East components are similar for all models, but there are large differences in the biases of the Up component (the Saastamoinen model: −0.1 cm, the UNB3 model: −8.1 cm, the NWP model: −13.2 cm). The standard deviations for all approaches are similar and equal to about 13 cm. The results lead to conclusion that more work should be carried out with the NWP approach to improve its performance. The aim of this paper is to propose a more accurate model that can be an alternative for the aforementioned solutions based on meteorological data (both ground-based and NWP).

The researchers from the Geodesy and Geodynamics Lab at ETH Zürich have developed the software package COMEDIE (Collocation of Meteorological Data for Interpretation and Estimation of Tropospheric Pathdelays) to interpolate and extrapolate meteorological parameters from real measurements to arbitrary locations (Eckert et al. 1992a, b; Hurter 2014). The software allows collocation and interpolation of meteorological parameters (temperature, air pressure, humidity) as well as ZTD and tropospheric refractivity. The most recent studies on this topic from Hurter and Maier (2013) are related to the 4-D interpolation of wet refractivity, dew point temperature and relative humidity for a 3-year period (2009–2011) in Switzerland. Authors reconstruct the wet refractivity profiles from GNSS data, ground-based meteorological parameters and radio occultations at the location of the radiosonde (RS) station in Payerne and validate the results against RS observations. The best collocation solution (combined GNSS and ground-based meteorological data) results in 3-year root mean square (RMS) between 2 and 7 ppm (corresponding to 5–80 % relative wet refractivity difference) below the height of 2 km and 4 ppm (130 % relative difference) at the height of 4 km. The results are improved by adding radio occultations, but there are only 189 radio occultations within the study area during the whole 3-year period. Radio occultation data mostly improve the accuracy in the upper troposphere. Maximum median offsets have decreased from a 120 % relative error to 44 % at the height of 8 km. The tomographic approach has been applied for a 1-year period in Payerne by Perler et al. (2011). The validation with the RS results in standard deviations of ${\sim }10$ ppm at the ground, which decrease to ${\sim }5$ ppm at the height of 4.5 km. Hence, the wet refractivity profiles from COMEDIE are proven to have a comparable accuracy to the results from GNSS tomography and partly mitigate the problem that path delays from ground-based GNSS stations have very limited capability to recover vertical structures in the atmosphere above the station.

Fig. 1

Location of the stations used in the study: GNSS stations from the AGNES network (left) and ground-based meteorological stations from the SwissMetNet network (right). Location of the RS in Payerne is marked with a star. Original version: Hurter and Maier (2013)

The herein presented study is a continuation of the Swiss research. We combine the ZTD and total refractivity (instead of the wet refractivity) from ground-based meteorological stations with the horizontal gradients assuming that information on the azimuthal asymmetry contained in the gradients improves the interpolation of the total refractivity. Moreover, we implement the collocation model for Poland, wherein the number of meteorological stations is limited and the differences of station altitudes are smaller. Therefore, we take advantage of a different data source, mainly NWP Weather Research and Forecasting (WRF) model.

Refractivity studies have been presented in many papers, with the main emphasis on the wet part. Various approaches have been applied: GNSS tomography (Notarpietro et al. 2011; Troller et al. 2002), NWP tomography (Rohm and Bosy 2011), algebraic reconstruction techniques (ART) (Bender et al. 2011), radio occultations (Heise et al. 2006), ray-tracing (Gegout et al. 2011) and many others. The advantage of using the collocation technique is that it can provide a methodology to investigate individual instrument accuracies, because of a relatively easy implementation of additional data sources.

This introduction is followed by Sect. 2 which presents the data sets from both regions: ground-based meteorological and GNSS observations from Switzerland as well as NWP and GNSS data from Poland. Section 3 describes the collocation technique in more detail. Section 4 presents the collocation procedure results for both countries. Section 5 discusses the obtained outcomes between countries and Sect. 6 summarizes the study.

2 Data

The total refractivity profiles from COMEDIE are calculated for two regions: western part of Switzerland and Poland. For each case, different data sets are used, because the selected countries exhibit different terrain conditions. Switzerland is mountainous, so the ground-based meteorological stations are located at various altitudes, which allows the reconstruction of the total refractivity profiles from ground-based stations. On the other hand, Poland is located mostly on lowlands with mountains up to 2.5 km in altitude in the south of the country. The height distribution of ground-based meteorological stations is too flat to reconstruct the refractivity profiles with the collocation technique. Moreover, the horizontal resolution of the stations is very sparse (50–70 km) and the stations are highly inhomogeneous (Hadaś et al. 2013). Thus, the main data sources for Poland are the NWP model and the GNSS observations.

2.1 Swiss data

We process the observations acquired in the 3-year period (1.01.2009–31.12.2011) from two main data sources: meteorological ground-based observations and GNSS products. The meteorological observations are provided by permanent and automatic weather stations (AWS) from the part of SwissMetNet network (20 stations) of Swiss Federal Office of Meteorology and Climatology (MeteoSwiss¹). The stations measure air pressure, temperature and relative humidity with a 10-minute resolution. Values of the total refractivity from AWS and RS stations are calculated according to Eq. 2. The GNSS path delays and horizontal gradients are calculated (with 1-h resolution) for 18 stations of Automated GNSS Network for Switzerland (AGNES) deployed by the Swiss Federal Office of Topography (swisstopo²). The ZTDs and horizontal gradients are retrieved from the Bernese GNSS Software TROPO files (Dach et al. 2015). The processing carried out by swisstopo is based on the same procedure as described in Perler et al. (2011). The results of interpolation are compared with the reference RS station at the MeteoSwiss Regional Center of Payerne. The RS launches are performed twice a day at 0:00 and 12:00 UTC. Station locations of both networks are presented in Fig. 1. The stations are located at various altitudes: mostly from 300 m to ${\sim }2$ km above mean sea level (AMSL), with two stations at ${\sim }3$ km AMSL (one GNSS and one AWS) and two stations at nearly 4 km AMSL (one GNSS and one AWS). The height distribution of both AWS and GNSS stations is shown in Fig. 2.

Fig. 2

Height distribution of meteorological (green) and GNSS (red) stations shown in Fig. 1. Location of radiosonde in Payerne is marked in yellow

2.2 Polish data

The collocation procedure is performed for 5 months of data (11.04.2014–15.09.2014) from two sources: GNSS and NWP model WRF.³ The WRF model for Poland is calculated and provided by the Department of Climatology and Atmosphere Protection of the University of Wroclaw⁴ (Kryza et al. 2013). The WRF is a mesoscale numerical weather prediction system computed in Poland for two nested domains: 10 km $\times$ 10 km (for the whole country, 96 $\times$ 111 horizontal nodes) and 2 km $\times$ 2 km (for south-west Poland, 281 $\times$ 258 horizontal nodes). The 10 km $\times$ 10 km grid with the 34 vertical, unevenly spaced $\sigma$-type levels is chosen (up to ${\sim }31$ km AMSL). The initial and boundary conditions are taken from the Global Forecast System (GFS) $0.5^\circ \times 0.5^\circ$ model provided by National Centers for Environmental Information.⁵ Using the dynamical downscaling technique, the meteorological parameters are calculated on a denser WRF grid (Kryza et al. 2016). There are no data assimilated into the WRF model and, thus, the independence between WRF and reference data (RS and GNSS) can be assumed. For computational purposes, only WRF points within a horizontal distance of 20 km from each interpolation point are used (with the original spacing, thus, all vertical information for each point is taken into account). The 24-h forecasts (with 1-h resolution) of meteorological parameters are calculated once a day. Both the analysis at 0:00 UTC and the following 23-h predictions are used in the collocation procedure. The uncertainties in terms of mean absolute errors assigned to the NWP outputs are calculated from a reanalysis for 1981–2010 based on 3-h data for synoptic stations in Poland. For particular meteorological parameters (air pressure p, temperature T and relative humidity RH), the uncertainties are equal to: $\mathrm{d}p=1$ hPa, $\mathrm{d}T=1.66$ K, $d\mathrm{RH}=8.93~\%$ (Kryza et al. 2016).

The GNSS ZTD is the second data source included in the collocation procedure. The GNSS and Meteo working group researchers at the Institute of Geodesy and Geoinformatics⁶ (IGG, WUELS) are calculating the GNSS ZTD at about 120 GNSS stations of the European Position Determination System Active Geodetic Network (ASG-EUPOS⁷) in Poland and adjacent areas (Fig. 3). A model called IGGHZ-G (H means 1 h interval, Z is the abbreviation of Zenith, G stands for GNSS) is computed using the Bernese GNSS Software version 5.2 (Dach et al. 2015). A more comprehensive description of the data acquisition can be found in Bosy et al. (2012). About 15 stations used in the model are part of the EUREF Permanent Network (EPN) and provide also the ground-based meteorological parameters with 1-h resolution (Fig. 3, in red). These stations are used as a reference source for ZTD interpolation comparisons. The other reference source for validation of COMEDIE outputs are RS observations. Three RS stations: LEBA, LEGIONOWO and WROCLAW are located in Poland (Fig. 3). Values of meteorological parameters (air pressure, temperature and dew point temperature, which is converted to water vapor partial pressure), are given as a vertical profile with 30–70 different height levels. The number of levels varies on a daily basis with weather conditions and RS performance; usually, it is only 30–40 levels. The measurements are taken twice a day (at 0.00 and 12.00 UTC). The RS data are retrieved from the US NOAA (National Oceanic and Atmospheric Administration) Earth System Research Laboratory website.⁸

Fig. 3

Location of the Polish stations used in the study: GNSS stations of the ASG-EUPOS network (left) and RS stations (right). Original version: asgeupos.pl (left) and Wilgan et al. (2015) (right)

3 Collocation technique

The most convenient way to represent the GNSS environmental propagation effects is to introduce the atmosphere refractivity. In the neutral atmosphere, the total refractivity can be expressed as a function of meteorological parameters: air pressure p, temperature T and water vapor partial pressure e (Essen and Froome 1951):

$$\begin{aligned} N_{\mathrm{tot}}=N_{\mathrm{dry}}+N_{\mathrm{wet}}=k_1\cdot \frac{p-e}{T}+k_2\cdot \frac{e}{T}+k_3\cdot \frac{e}{T^2}, \end{aligned}$$

(2)

where $k_1=77.689 \text { K}\cdot \text {hPa}^{-1}$, $k_2=71.2952\text { K}\cdot \text {hPa}^{-1}$ and $k_3=375{,}463\text { K}^2\cdot \text {hPa}^{-1}$ are empirically determined coefficients. In this study, values given by Rüeger (2002) are adopted.

Table 1 Stochastic parameters applied in the least-squares collocation

The total propagation delay $\Delta ^{\text {PD}}$ can be expressed as an integral of the total refractivity $N_{\mathrm{tot}}$ along the propagation path s from receiver r to the satellite w:

$$\begin{aligned} \Delta ^{\text {PD}}=10^{-6}\int _r^w N_{\mathrm{tot}} \mathrm{d}s. \end{aligned}$$

(3)

For each pair of satellite-receiver the delay can be estimated individually, which is computationally very demanding. The usual approach is to calculate one delay in the zenith direction (ZTD) for each receiver and then to project this delay using mapping functions. The scope of this research covers only GNSS ZTDs, which can be determined as an integral of $N_{\mathrm{tot}}$ in the zenith direction:

$$\begin{aligned} \mathrm{ZTD}=10^{-6}\int _{\text {zenith direction}}N_{\mathrm{tot}} \text { ds}. \end{aligned}$$

(4)

It is not feasible to measure the refraction of the atmosphere at all points along the signal path; therefore, a method to infer these conditions to the arbitrary locations has to be developed. In this study, the total refractivity profiles and ZTDs are calculated by utilizing the least-squares collocation technique, which is based on the adjustments of the measurements to the deterministic part f(u, x, t) and to the stochastic parts s and $\epsilon$ (Troller 2004):

$$\begin{aligned} l=f(\mathbf {u},\mathbf {x},t)+s(\mathbf {C}_{\mathrm{ss}},\mathbf {x},t)+\epsilon , \end{aligned}$$

(5)

where l is the measurement, f($\mathbf {u},\mathbf {x},t$) is the function describing general field of measured values, $\mathbf {u}$ are the unknown parameters, $\mathbf {x},t$ are the coordinates in space and time, $s(\mathbf {C}_{\mathrm{ss}},\mathbf {x},t)$ is the stochastic parameter  $s\sim \mathcal {N}(0;\mathbf {C}_{\mathrm{ss}})$ (signal), $\epsilon$ is the stochastic parameter  $\epsilon \sim \mathcal {N}(0;\mathbf {C}_{\epsilon \epsilon })$ (noise).

The parameters of the deterministic and signal parts, estimated in the least-squares collocation procedure, allow the interpolation of both to the points where no measurements are available.

3.1 Collocation of ZTD

The least-squares collocation requires the deterministic model of the considered parameter to describe the general trends in the measurements. In this study, the following model of ZTD is utilized (modified after Hurter and Maier 2013):

$$\begin{aligned} \mathrm{ZTD}(x,y,z,t)= & {} [\mathrm{ZTD}_0+a_{\mathrm{ZTD}}(x-x_0)+b_{\mathrm{ZTD}}(y-y_0)\nonumber \\&+\,c_{\mathrm{ZTD}}(t-t_0)]\cdot e^{-\frac{z-z_0}{H_{\mathrm{ZTD}}}}, \end{aligned}$$

(6)

where $x_0, y_0, z_0=0, t_0$ are the coordinates of reference point and reference time, x, y, z, t are the coordinates and time of investigated point, $ZTD_0$ is the ZTD at a reference position and time, $H_{\mathrm{ZTD}}$ is the scale height, $a_{\mathrm{ZTD}}, b_{\mathrm{ZTD}}, c_{\mathrm{ZTD}}$ are the gradient parameters in  x, y and time, respectively.

The unknown parameters: $\mathrm{ZTD}_0$ at a reference position and time, scale height $H_{\mathrm{ZTD}}$, and gradients in x, y direction and time ($a_{\mathrm{ZTD}},b_{\mathrm{ZTD}},c_{\mathrm{ZTD}}$) are estimated for each time batch during the collocation procedure to allow the calculation of the deterministic part of the model during interpolation.

The signal part s of the collocation model is assumed to be of normal distribution with mean 0 and the covariance matrix $\mathbf {C_{\mathrm{ss}}}$. The matrix is described by a covariance function of the distances between the measurements, which shows the spatial and temporal dependencies:

$$\begin{aligned} C_{\mathrm{ss}}(\mathrm{ZTD}_k,\mathrm{ZTD}_l)=\frac{\sigma _{\mathrm{signal}}^2}{q}, \end{aligned}$$

(7)

where q is a scaling factor that increases the correlation length with height above the ground:

$$\begin{aligned} q= & {} 1+ \left[ \left( \frac{x_k-x_l}{\Delta x_0}\right) ^2+\left( \frac{y_k-y_l}{\Delta y_0}\right) ^2\right. \nonumber \\&\left. +\left( \frac{z_k-z_l}{\Delta z_0}\right) ^2+\left( \frac{t_k-t_l}{\Delta t_0}\right) ^2 \right] \cdot e^{-\frac{z_k+z_l}{2z_0}}, \end{aligned}$$

(8)

where $\sigma _{\mathrm{signal}}^2$ is the a priori covariance of signal, $x_k,y_k,z_k,t_k$ are the Cartesian coordinates and time of  observation  $k,~x_l,y_l,$ $z_l,t_l$ are the Cartesian coordinates and time of  observation  l, $z_0$ is the scale height modifying the correlation lengths, as a function of height, $\Delta x_0,\Delta y_0,\Delta z_0,\Delta t_0$ are the correlation lengths of space and time.

The stochastic parameters (correlation lengths) in Table 1 were empirically developed by Hirter (1996) and $\sigma _{\text {signal}}$ is an average formal uncertainty of ZTDs from GNSS processing of L1/L2 dual-frequency geodetic GNSS observations.

The stochastic parameter $\epsilon$ is described by the covariance matrix $\mathbf {C}_{\epsilon \epsilon }$ which consists of the noise of particular measurements on the diagonal. The off-diagonal elements are equal to zero. The uncertainties used to calculate the noise part are provided for each data source separately.

3.2 Collocation of ZTD and total refractivity

The total refractivity $N_{\mathrm{tot}}$ can be expressed as the derivative of ZTD in zenith direction. Thus, if the ZTD observations $l_{\mathrm{ZTD}}$ are described as:

$$\begin{aligned} l_{\mathrm{ZTD}}=f(\mathbf {u},\mathbf {x},t)+s(\mathbf {C}_{\mathrm{ss}},\mathbf {x},t)+\epsilon , \end{aligned}$$

(9)

then, the $N_{\mathrm{tot}}$ observations can be related to $l_{\mathrm{ZTD}}$ using the differential operator in zenith direction $D=-\frac{\partial }{\partial z}$:

$$\begin{aligned} l_{N_{\mathrm{tot}}}=D\left( f(\mathbf {u},\mathbf {x},t)+s(\mathbf {C}_{\mathrm{ss}},\mathbf {x},t)+\epsilon \right) . \end{aligned}$$

(10)

Applying the operator D to the deterministic part of the ZTD (Eq. 6) results in a deterministic model of $N_{\mathrm{tot}}$:

$$\begin{aligned} N_{\mathrm{tot}}(x,y,z,t)&=~D \mathrm{ZTD}(x,y,z,t)\nonumber \\&=\frac{1}{H_{\mathrm{ZTD}}}[\mathrm{ZTD}_0 +a_{\mathrm{ZTD}}(x-x_0)\nonumber \\&\quad +b_{\mathrm{ZTD}}(y-y_0)+c_{\mathrm{ZTD}}(t-t_0)]\cdot e^{-\frac{z}{H_{\mathrm{ZTD}}}}. \end{aligned}$$

(11)

To derive the covariance functions between $N_{\mathrm{tot}}$ and ZTD as well as between $N_{\mathrm{tot},k}$ and $N_{\mathrm{tot},l}$ the differential operator D is applied to the covariance function of ZTD (Eq. 7):

$$\begin{aligned} \begin{array}{ll} &{}\mathbf {C}_{\mathrm{ss}}\left( N_{\mathrm{tot}},\mathrm{ZTD}\right) = \mathbf {C}_{\mathrm{ss}}\left( \mathrm{ZTD},N_{\mathrm{tot}}\right) \\ &{}\quad = \frac{\sigma _{\mathrm{signal}}^2}{q^2}\left[ \frac{2\cdot (-z_{\mathrm{ZTD}}+z_{N_{\mathrm{tot}}})}{(\Delta z_0)^2}\cdot e^{-\frac{z_{N_{\mathrm{tot}}}+z_{\mathrm{ZTD}}}{2z_0}}+\frac{1-q}{2z_0}\right] , \end{array} \end{aligned}$$

(12)

$$\begin{aligned} \begin{array}{ll} &{}\mathbf {C}_{\mathrm{ss}}\left( N_{\mathrm{tot},k},N_{\mathrm{tot},l}\right) = \mathbf {C}_{\mathrm{ss}}({N_{\mathrm{tot},l}}, {N_{\mathrm{tot},k}})\\ &{}\quad = \frac{2\sigma _{\mathrm{signal}}^2}{q^2} \left[ \frac{e^{-\frac{z_k+z_l}{2z_0}}}{(\Delta z_0)^2}+\frac{(q-1)(q-2)}{8qz_0^2}- \frac{4(z_k-z_l)^2}{q(\Delta z_0)^4}\cdot e^{-\frac{z_k+z_l}{z_0}} \right] . \end{array} \end{aligned}$$

(13)

Under the influence of the differential operator D, the uncorrelated noise $\epsilon$ of ZTD becomes the uncorrelated noise of $N_{\mathrm{tot}}$.

3.3 Collocation of ZTD, total refractivity and horizontal gradients

The horizontal gradients of the ZTD distribution $\nabla \mathrm{ZTD} = \left( \frac{\partial \mathrm{ZTD}}{\partial x}, \frac{\partial \mathrm{ZTD}}{\partial y}\right)$ are introduced into the collocation procedure with the assumption that they can improve the interpolation as they contain information about the azimuthal asymmetry in the tropospheric delay. The $\nabla \mathrm{ZTD}$ can be related to the result of GNSS processing: atmospheric delay gradient $G=(G_E,G_N)$, with some assumptions about the vertical distribution of water vapor (Ruffini et al. 1999; Shoji 2013). For instance, an assumption of the exponential horizontal refractivity distribution with scale height $H_{\mathrm{ZTD}}$ implies:

$$\begin{aligned} G=H_{\mathrm{ZTD}}\cdot \nabla \mathrm{ZTD}. \end{aligned}$$

(14)

From this relation, the horizontal gradient models are derived as follows:

$$\begin{aligned} G_{E}=H_{\mathrm{ZTD}}\cdot \frac{\partial ZTD}{\partial x} = H_{\mathrm{ZTD}}\cdot a_{\mathrm{ZTD}} \cdot e^{-\frac{z}{H_{\mathrm{ZTD}}}}, \end{aligned}$$

(15)

$$\begin{aligned} G_{N}=H_{\mathrm{ZTD}}\cdot \frac{\partial \mathrm{ZTD}}{\partial y} = H_{\mathrm{ZTD}}\cdot b_{\mathrm{ZTD}} \cdot e^{-\frac{z}{H_{\mathrm{ZTD}}}}. \end{aligned}$$

(16)

In the data set with ZTDs, total refractivities and horizontal gradients, there are four different types of variables that are taken into consideration during the collocation procedure. Therefore, the covariance matrix consists of 4 by 4 segments, where each segment is described by a separate covariance function:

$$\begin{aligned} \mathbf {C}_{\mathrm{ss}}= \begin{bmatrix} C_{\mathrm{ss}}\left( \mathrm{ZTD},\mathrm{ZTD}\right)&C_{\mathrm{ss}}\left( \mathrm{ZTD}, N_{\mathrm{tot}}\right)&C_{\mathrm{ss}}\left( \mathrm{ZTD},G_E \right)&C_{\mathrm{ss}}\left( \mathrm{ZTD},G_N \right) \\ C_{\mathrm{ss}}\left( N_{\mathrm{tot}},\mathrm{ZTD}\right)&C_{\mathrm{ss}}\left( N_{\mathrm{tot}},N_{\mathrm{tot}}\right)&C_{\mathrm{ss}}\left( N_{\mathrm{tot}},G_E\right)&C_{\mathrm{ss}}\left( N_{\mathrm{tot}},G_N\right) \\ C_{\mathrm{ss}}\left( G_E,\mathrm{ZTD} \right)&C_{\mathrm{ss}}\left( G_E,N_{\mathrm{tot}}\right)&C_{\mathrm{ss}}\left( G_E,G_E \right)&C_{\mathrm{ss}}\left( G_E,G_N \right) \\ C_{\mathrm{ss}}\left( G_N,\mathrm{ZTD} \right)&C_{\mathrm{ss}}\left( G_N,N_{\mathrm{tot}}\right)&C_{\mathrm{ss}}\left( G_N,G_E \right)&C_{\mathrm{ss}}\left( G_N,G_N \right) \\ \end{bmatrix} \end{aligned}$$

(17)

The covariance functions for particular segments are derived from the covariance function of $\mathbf {C}_{\mathrm{ss}}\left( \mathrm{ZTD},\mathrm{ZTD}\right)$ (Eq. 7). Functions that include only ZTD and/or $N_{\mathrm{tot}}$ are described in Sect. 3.1 and 3.2. The remaining 7 functions are calculated analogically to Eqs. 12 or 13 and are described in Appendix A.

Fig. 4

Solution flowchart used in the framework of the study. AWS meteorological parameters are used for Switzerland while NWP parameters for Poland

3.4 Processing

The flowchart in Fig. 4 gives an overview of the work undertaken during this study. The parallelograms denote the data sets and the rectangles show the processing steps. Firstly, the combinations of data sets are established to be used in the collocation procedure. For Switzerland, three main data sets are considered:

• 1. ‘AWS only’ that includes the total refractivity calculated from ground-based AWS measurements using Eq. 2.

• 2. ‘AWS/GNSS’ that includes the total refractivity calculated from AWS measurements and ZTD from GNSS stations.

• 3. ‘AWS/GNSS/GRAD’ that includes the total refractivity calculated from AWS measurements, as well as ZTD and horizontal gradients of ZTD from GNSS stations.

For Poland, the data sets are different as the main data source is the NWP model WRF:

• 4. ‘WRF only’ that includes the total refractivity calculated from the WRF model using Eq. 2.

• 5. ‘WRF/GNSS’ that includes the total refractivity calculated from the WRF model and ZTD from GNSS stations.

• 6. ‘WRF/GNSS/AWS’ that includes the total refractivity calculated from the WRF model and from AWS that are a part of the EPN as well as ZTD from all GNSS stations.

• 7. ‘GNSS only’ that includes only ZTD from the GNSS network.

The collocation procedure is carried out for each combination of data sets 1–7. The deterministic models and covariance functions introduced in previous sections are used. When the input data are ZTDs, the deterministic model described by Eq. 6 is used. For the total refractivity, the model described by Eq. 11 is utilized and for horizontal gradients the model is represented by Eqs. 15 and 16. All models are employed with corresponding covariance matrices. During the procedure, the collocation parameters $\mathrm{ZTD}_0$, $H_{\mathrm{ZTD}}$, $a_{\mathrm{ZTD}}$, $b_{\mathrm{ZTD}}$ and $c_{\mathrm{ZTD}}$ are calculated for each time interval separately. To speed up the computation, each day is divided into 8-h batches with 1-h overlap to the previous and next batch for the Swiss data and in 12-h batches with 1-h overlap for the Polish data. Using the obtained collocation parameters, it is possible to interpolate the total refractivity and ZTD values to the locations of the RS and GNSS stations and validate the results from COMEDIE with the reference data. Interpolated total refractivity values are compared to RS, while the interpolated ZTDs are compared to both reference data sources.

Fig. 5

Biases and standard deviations of the total refractivity differences ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for two data sets: ‘AWS/GNSS’ (red) and ‘AWS/GNSS/GRAD’ (dashed blue). The data period is 1.01.2009–31.12.2011

Fig. 6

Biases and standard deviations of the total refractivity differences ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for two data sets: ‘AWS/GNSS’ (red) and ‘AWS only’ (green). The data period is 5.06.2011–31.12.2011

Fig. 7

Differences of the total refractivity ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for ‘WRF only’ data set. The data period is 11.04.2014–15.09.2014

4 Results

4.1 Switzerland

We reconstruct the total refractivity profiles at the RS station in Payerne according to the processing steps presented in Sect. 3.4. Firstly, we perform the collocation procedure on the data set ‘AWS/GNSS’. In the previous study from Hurter and Maier (2013), this data set exhibits the best accuracy w.r.t. the RS data in Payerne. Figure 5 presents biases and standard deviations from residuals $N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$ w.r.t. the heights of the RS. The comparisons are cut at 4 km AMSL, as the highest ground-based station is located at ${\sim }$3.6 km. Biases vary from −7 to 3 ppm with standard deviations from 3 to 9 ppm (averaged from the 3-year period). The differences in the residual distribution may be a result of the height resolution with very few stations above the height of 2 km.

In the next step, we include the horizontal gradients into the collocation procedure. The biases and standard deviations of residuals for the data set ‘AWS/GNSS/GRAD’ are also shown in Fig. 5. Unfortunately, including gradients into the collocation procedure does not improve the interpolation. The model with gradients exhibits larger biases of about 0.5 ppm than the model without gradients, with the largest difference of 2 ppm at the height of about 2.5 km. The deterioration of the total refractivity values when including the gradients raises the question of whether the GNSS ZTD information is necessary for the interpolation of the total refractivity field. Therefore, we investigate the total refractivity values from the data sets: ‘AWS/GNSS’ and ‘AWS only’. Figure 6 shows biases and standard deviations for residuals $N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$ for a shorter period 5.06.2011–31.12.2011 (due to the long computational time). For data set ‘AWS/GNSS’, the statistics for the 7-month period are consistent with the 3-year period. The results show that excluding the ZTD values from the model has a strong influence on the model accuracy. The absolute biases are larger for the data set ‘AWS only’ of more than 10 ppm above the height of 2 km, where there are only 2 ground-based meteorological stations. Therefore, including the GNSS ZTDs is necessary to obtain a good accuracy in the vertical direction of the model. We do not perform the collocation procedure in Switzerland for the data set ‘GNSS only’ because as shown in Hurter and Maier (2013) this data set has worse accuracy than ‘AWS/GNSS’ by about 5 ppm (for wet refractivity).

4.2 Poland

4.2.1 Interpolation of the total refractivity

We perform the least-squares collocation of the total refractivity based on different combinations of data sets and interpolate the obtained values to the locations of reference data as described in Sect. 3.4 for 5 months in 2014 (11.04.2014–15.09.2014). For Switzerland, the results from the 3-year period are consistent with the 7-month period; therefore, we decide to process a shorter period for Poland, where the collocation procedure is more computationally demanding. We use the summer months as the worst-case scenario, because the collocation results are likely to deviate stronger from the RS in summer than they do in winter. We employ four data sets for the Poland case study: ‘WRF only’, ‘WRF/GNSS’, ‘WRF/GNSS/AWS’ and ‘GNSS only’. We decide not to include horizontal gradients of ZTD into the collocation of the Polish data, because gradients in the Swiss data worsen the interpolation (Sect. 4.1).

In the first step, we calculate the total refractivity profiles from the data set ‘WRF only’. Figure 7 shows differences between profiles from RS station WROC and COMEDIE. For other RS stations, the results are similar and as such they are not shown. Near the Earth’s surface, the differences are larger, sometimes at the level of −10 ppm but they are decreasing with height. At middle levels (height of 5–15 km), the differences are much smaller, positive and higher in summer months (June – August) than in spring (April–May). For upper levels, the differences are again negative on the order of −2 ppm. It is worth noticing that the height of the model exceeds 30 km, because the NWP model reaches much higher altitudes than ground-based meteorological stations in Switzerland (where the comparisons are performed only up to 4 km).

Fig. 8

Differences of the total refractivity ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for ‘WRF/GNSS’ data set. The data period is 11.04.2014–15.09.2014

Fig. 9

Differences of the total refractivity ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for ‘GNSS only’. The data period is 11.04.2014–15.09.2014. Please note that the scale on the color-bar is different from the previous plots (Figs. 7, 8)

In the next step, we add the GNSS ZTDs to the collocation procedure. Figure 8 shows the differences between values from RS and calculated from the data set ‘WRF/GNSS’. There are only small differences between the data sets ‘WRF only’ and ‘WRF/GNSS’ at lower levels and also at the highest levels the data set ‘WRF/GNSS’ seems to perform slightly better. To test all the possible data sets, we also interpolate the total refractivity values from the ZTDs set ‘GNSS only’. As shown in Fig. 9, the collocation from ‘GNSS only’ gives much worse results than collocation from the sets that contain WRF data. At lower levels, the differences are equal to approximately 30 ppm. At middle levels, differences of 10 ppm are frequent. Only upper levels show slightly better agreement with a 5 ppm difference. The results are easy to predict, because we attempt to reconstruct a whole profile of refractivity from a single ZTD value. To achieve this goal, some more sophisticated techniques should be exploited (e.g., GNSS tomography). For the collocation technique, the vertical dense distribution of WRF data is necessary to achieve an accurate refractivity interpolation.

Fig. 10

Boxplots of the total refractivity differences $N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$ and relative differences $(N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}})/N_{\mathrm{COMEDIE}}*100~\%$ for chosen data sets: WRF only (top), WRF/GNSS (upper middle), WRF/GNSS/AWS (lower middle) and GNSS only (bottom). Boxes denote the 25th and 75th percentile. The median is marked inside the boxes. Lines show offsets from $q_{25~\%}-1.5(q_{75~\%}-q_{25~\%})$ to $q_{25~\%}+1.5(q_{75~\%}-q_{25~\%})$. The red dots denote outliers. The data period is 11.04.2014–15.09.2014. Please note that the scale for ‘GNSS only’ boxplots is different from the other plots

Fig. 11

Residuals of ZTDs ($\mathrm{ZTD}_{\mathrm{RS}}-\mathrm{ZTD}_{\mathrm{model}}$), based on one of the following models: ‘WRF only’ from COMEDIE (red), ‘WRF/GNSS’ from COMEDIE (green), ‘WRF/GNSS/AWS’ from COMEDIE (cyan) or direct GNSS observation for station WROC (black). Data period: 11.04.2014–15.09.2014

Table 2 Mean biases and standard deviations of residuals $\mathrm{ZTD}_{\mathrm{RS}}-\mathrm{ZTD}_{\mathrm{model}}$ (mm)

We compare all the data sets involved into the collocation procedure in Fig. 10. Boxplots show the total refractivity differences for sets: ‘WRF only’, ‘WRF/GNSS’, ‘WRF/GNSS/AWS’ and ‘GNSS only’. The advantage of using the data sets that include the WRF is evident. We intend to improve the refractivity interpolation even more by adding ground-based meteorological measurements from AWS, but data sets ‘WRF/GNSS’ and ‘WRF/GNSS/AWS’ show very similar accuracy w.r.t. reference RS observations. Thus, there is no clear advantage in including AWS information. Figures 7, 8 and 9 show very small differences in the upper levels, but it is important to note that the values of the total refractivity decrease exponentially with height and at the topmost level of the model (${\sim }30$ km) they are about 60 times smaller than at the bottommost level. We include the boxplots of fractional differences between RS and COMEDIE, which show that the relative differences are much larger at the topmost levels of the model, even though the absolute differences seem insignificant. To evaluate the accuracy of the model, we must consider residuals (absolute and relative) on every level.

4.2.2 Interpolation of ZTD

For the final assessment, we calculate the ZTD values as an integral from the total refractivities (Eq. 4) obtained using COMEDIE. Included data sets are ‘WRF only’, ‘WRF/GNSS’ and ‘WRF/GNSS/AWS’. Figure 11 shows differences between $\mathrm{ZTD}_{\mathrm{RS}}$ calculated from RS observations and $\mathrm{ZTD}_{\mathrm{model}}$, where ‘model’ is one of the COMEDIE data sets or direct GNSS observation at the station WROC (which at this point becomes another validation data source). Table 2 shows statistics for all aforementioned models.

The differences between the residuals from the ‘WRF only’ and ‘WRF/GNSS’ data sets indicate that GNSS information is important to improve the accuracy of the ZTD interpolation. The standard deviation from the ‘WRF only’ set is more than 15 mm larger than from all data sets containing GNSS results. The data sets ‘WRF/GNSS’ and ‘WRF/GNSS/AWS’ are again nearly identical, which additionally indicates that there is no need to include the ground-based meteorological stations. The standard deviation of the differences between the two reference data sources RS and GNSS is 12.8 mm. In our previous investigations (Wilgan et al. 2015), the standard deviation of residuals $\mathrm{ZTD}_{\mathrm{RS}}-\mathrm{ZTD}_{\mathrm{GNSS}}$ for station WROC was 9.9 mm, but the study was conducted during winter months, when the water vapor content is the smallest. The distance between RS WROCLAW and the closest GNSS station is ${\sim }$6 km; for RS LEGIONOWO, it is ${\sim }$9 km and for RS LEBA ${\sim }$40 km. The improvement after adding the GNSS observations into the collocation procedure is visible for all RS stations, but for LEBA the impact is much smaller.

It is important to acknowledge that the RS measurements are not error free. The $\mathrm{ZTD}_{\mathrm{RS}}$ error is on a similar level to the $\mathrm{ZTD}_{\mathrm{GNSS}}$ error. The RS needs time to ascent through the atmosphere (${\sim }$1 h to reach the highest altitude) and there might be some variations in the atmosphere during this time. Moreover, there are only 3 RS stations in Poland, too few to test the collocation procedure across the whole country. Thus, we interpolate the ZTD values to the locations of the GNSS stations and compare the outputs from COMEDIE with the reference GNSS data.

Fig. 12

Comparison of ZTD from GNSS station REDZ (black) with ZTDs from COMEDIE from 3 data sets: ‘WRF only’ (red), ‘WRF/GNSS’ (green) and ‘GNSS only’ (blue) (top) with corresponding residuals of $\mathrm{ZTD}_{\mathrm{GNSS}}-\mathrm{ZTD}_{\mathrm{COMEDIE}}$ (mm) for all data sets (bottom). Data period is 23–31.05.2014

Fig. 13

Comparison of ZTD from GNSS station KRAW (black) with ZTDs from COMEDIE from 3 data sets: ‘WRF only’ (red), ‘WRF/GNSS’ (green) and ‘GNSS only’ (blue) (top) with corresponding residuals of $\mathrm{ZTD}_{\mathrm{GNSS}}-\mathrm{ZTD}_{\mathrm{COMEDIE}}$ (mm) for all data sets (bottom). Data period is 23–31.05.2014

Table 3 Mean biases and standard deviations of residuals $\mathrm{ZTD}_{\mathrm{GNSS}}-\mathrm{ZTD}_{\mathrm{COMEDIE}}$ (mm)

Fig. 14

Boxplots of total refractivity differences for Poland (left) and Switzerland (right) up to 4 km AMSL height. Boxes denote the 25th and 75th percentile. The median is marked inside the boxes. Lines show offsets from $q_{25~\%}-1.5(q_{75~\%}-q_{25~\%})$ to $q_{25~\%}+1.5(q_{75~\%}-q_{25~\%})$. The red dots denote outliers. The data period is 11.04.2014–15.09.2014 for Poland and 1.01.2009–31.12.2011 for Switzerland

We choose 9 evenly distributed GNSS stations that are a part of the EPN. The station at which we interpolate the results is excluded from the collocation procedure. The chosen stations have the smallest biases and standard deviations w.r.t. the EPN final combined troposphere product presented in Bosy et al. (2012). We test three data sets to find the best solution for the interpolation of ZTDs: ‘WRF only’, ‘WRF/GNSS’ and ‘GNSS only’ for nine days in 2014 (23–31.05.2014), which contain a severe weather event (26–28.05.2014). The event is a torrential rain associated with strong movements of the ascending air within the large convection cells. In the days before the main event, the rainfalls were also recorded, but not as severe. Thus, we divide the solution into three 3-day periods: with the severe rainfall (‘heavy rainfall’), before the event (‘moderate rainfall’) and after the event (‘after rainfall’). Table 3 presents biases and standard deviations for all 9 stations and all 3 data sets with the division for 3 periods. Figures 12 and 13 show the comparisons of the interpolation results w.r.t. GNSS data for two sample stations: Redzikowo (REDZ) and Kraków (KRAW), respectively. The first station is located at the Baltic sea coast, with only a few GNSS stations nearby and the second one in the south of Poland (but not yet in the mountains) with many stations nearby, along with other EPN stations. In the case of KRAW and other stations surrounded by many GNSS stations, the best performing data sets are: ‘WRF/GNSS’ or ‘GNSS only’. On average, the ‘WRF/GNSS’ data set is better by less than a millimeter. Both sets that contain the GNSS data have smaller standard deviations than ‘WRF only’ set by about 14 mm on average. In case of REDZ, adding GNSS data is not improving the collocation results as strongly as for other stations. Moreover, the ‘GNSS only’ data set has worse accuracy than ‘WRF only’, especially during the rainfall. For all the stations, the ‘WRF only’ data set performs much better during the ‘heavy rainfall’ period compared to the ‘after rainfall’ period. After the rainfall, the $\mathrm{ZTD}_{\mathrm{WRF}}$ values still remain at a similar level as during the rainfall, but in reality the ZTD values drop significantly and this is visible in $\mathrm{ZTD}_{\mathrm{GNSS}}$. The reason for such behavior is that the humidity values provided by the WRF model after the rainfall are too high. We experienced similar problems in the past with another NWP model (COAMPS) producing too wet conditions (Wilgan et al. 2015). Moreover, for many stations, the WRF data are experiencing problems also during the ‘moderate rainfall’ period, but without a clear trend as in the ‘after rainfall’ period. Therefore, it is important to support the collocation procedure with the GNSS information during severe weather events.

5 Discussion

The terrain conditions of Poland and Switzerland are very different; therefore, various data sets are involved in the collocation procedure. In Switzerland, the vertical resolution of ground-based stations is diversified and the horizontal resolution is dense, so it is possible to reconstruct the state of the troposphere solely from the ground-based stations and GNSS data. In Poland, most of the ground-based meteorological stations are located in the lowlands with a horizontal resolution of 50–70 km. Thus, we employ the NWP data to reconstruct vertical profiles. Even though the terrain resolution and the data collection schemes are diverse, this paper shows that the same procedure can be applied for both countries to obtain ZTD and total refractivity models of similar accuracy. As an example, we compare the total refractivity results between the countries. Figure 14 shows the boxplots of residuals $N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$ for the best possible solutions: ‘AWS/GNSS’ for Switzerland and ‘WRF/GNSS’ for Poland w.r.t. heights of the RS. For Poland, the heights are cut at 4 km AMSL. Note that the sampling of the RS data is different for Poland and Switzerland, hence the different resolution of the plots. For the lowest levels, the differences for Switzerland are smaller than for Poland, but after the model reaches the level where there are not many ground-based stations (${\sim }$2 km), the results for Switzerland are worse than for Poland, where the median is close to 0. For both countries, most of the residual values are within ±5 ppm limits.

Utilizing the NWP model, we are not restricted by the number or location of the ground-based meteorological stations. In flat countries like Poland, an NWP model is a valuable data source containing the information on the vertical variability of the atmospheric parameters. Furthermore, we can reach much higher altitudes (here, ${\sim }$30 km AMSL) than in the case of ground-based stations in Switzerland, where we are constrained by the height of the highest station (${\sim }$4 km AMSL). Above that level, we need to extrapolate the results, which affect the accuracy significantly. Also, in countries where the AWS network is highly inhomogeneous and the horizontal resolution of stations is very sparse, the NWP model provides a good database of the atmospheric state. The disadvantage of using the NWP model is the time needed to compute the collocation, as there are much more data to process than in the case of ground-based meteorological stations. Therefore, our recommendation for providing a priori model for positioning (especially for the PPP strategy) is to use a model based on ‘AWS/GNSS’ for countries that can achieve the height variability solely from the ground-based meteorological stations and ‘NWP/GNSS’ for flat countries.

Fig. 15

Biases and standard deviations of the total refractivity differences ($N_{\mathrm{RS}}-N_{\mathrm{COMEDIE}}$) for data sets: ‘AWS/GNSS (ZTD)’ (red), ‘AWS/GNSS/GRAD (ZTD)’ (dashed blue) ‘AWS/GNSS (ZHD + ZWD)’ (green), and ‘AWS/GNSS/GRAD (ZHD + ZWD)’ (dashed cyan), where (ZHD + ZWD) means that the hydrostatic and wet parts of zenith delay were modeled separately. The data period is 1.01.2009–31.12.2011

The Switzerland case study of total refractivity has similar accuracy to the study on wet refractivity performed by Hurter and Maier (2013) presented in Sect. 1. For further improvements, we include the horizontal gradients, but unfortunately they worsen the collocation results. We attempt to improve the collocation with gradients. Firstly, we remove the gradients with relatively high errors: gradients that are smaller than their $3\sigma$ error from the Bernese log-files. From the original data set 28.45 % of gradients are removed, but there is almost no improvement after the reduction. In our next attempt, we remove gradients that do not follow some major patterns based on the hourly distribution of gradients. We remove gradients that do not satisfy $|\sphericalangle G - \text {median}|<\frac{\pi }{6}$. The elimination of some gradients improves the interpolation, but the data set with gradients is still worse than the one without gradients. A possible reason is that the gradient model may be wrong (the chosen model is very simple to hold the relation with ZTD and $N_{\mathrm{tot}}$ models) or the stochastic parameters shown in Table 1 are not adequate, since wet and dry gradients are likely to significantly differ in their correlation lengths. Therefore, we divide the ZTD model into the wet and hydrostatic parts ($\mathrm{ZTD} = \mathrm{ZHD} + \mathrm{ZWD}$), with the analogous model for ZWD and ZHD as the one for ZTD (Eq. 6). The divided model has 10 unknown collocation parameters (the wet and dry scale heights are considered separately). The models for gradients are also divided into the wet and dry parts. Figure 15 shows that the statistics for the divided model are only slightly better above 2.5 km than for the undivided model. Furthermore, we utilize only the horizontal gradients from GNSS processing, but there are many more techniques from which we can retrieve gradients like NWP, Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), Very Long Baseline Interferometry (VLBI) water vapor radiometer (WVR) and many others. The gradients retrieved for the same station by different techniques can differ significantly (e.g., Teke et al. (2013); Douša et al. (2016)); hence, the physical meaning of the gradients may not be understood well.

6 Summary

We investigated the total refractivity profiles calculated using the collocation software COMEDIE for two regions: western part of Switzerland and Poland. For Switzerland, three data sets were used: ‘AWS/GNSS’, ‘AWS/GNSS/GRAD’ and ‘AWS only’. The refractivity profiles were compared to the reference RS station in Payerne. The data set with the best performance of the total refractivity interpolation was ‘AWS/GNSS’ with the absolute biases from 0 to 7 ppm and standard deviations from 3 to 9 ppm. Excluding the ZTDs from the collocation results worsened the interpolation above a height of 2 km by more than 10 ppm. Introducing the horizontal gradients was presumed to improve the interpolation, but the data set with gradients had larger biases by about 0.5 ppm for most of the profile (up to 4 km) with the maximum offset of 2 ppm at 2.5 km height than the data set without gradients. For Poland, we considered four data sets: ‘WRF only’, ‘WRF/GNSS’, ‘WRF/GNSS/AWS’ and ‘GNSS only’. The data set with the best accuracy was the combined ‘WRF/GNSS’. The data set ‘WRF/GNSS/AWS’ exhibited very similar accuracy. Therefore, we see no benefit from including ground-based meteorological information from AWS into the collocation procedure. The data set ‘GNSS only’ showed much worse accuracy with the discrepancies at lower altitudes even at the level of −30 ppm. The data set ‘WRF only’ showed similar agreement with reference data as ‘WRF/GNSS’ in terms of the total refractivity, but for the interpolation of ZTDs from all sets, the standard deviations from residuals were almost two times larger for the ‘WRF only’ set than for all data sets containing GNSS results. We also performed the interpolation of ZTDs at the locations of GNSS stations for three data sets: ‘WRF only’, ‘WRF/GNSS’ and ‘GNSS only’ for nine days in May 2014 that contained a severe weather event. The best data sets with similar accuracy were: ‘WRF/GNSS’ and ‘GNSS only’ with average biases of 3.7 and 3.8 mm and average standard deviations of 16.7 and 17.2 mm, respectively.

We can conclude that the best troposphere models based on collocation can be obtained from the combination of meteorological (NWP or AWS) and GNSS data. Using NWP-only data biases, the troposphere delays in particular due to the overestimation of the humidity after rainfalls. Using GNSS-only data provides substantially larger differences of the total refractivity with respect to the RS measurements. Whereas, using a NWP/AWS-GNSS combination results in the smallest biases and the smallest residuals with respect to both, RS and GNSS data. We recommend to use these models as an a priori model of the troposphere for positioning.

Appendix A

In Sect. 3.3, we present the covariance matrix used in the collocation procedure for data set ‘AWS/GNSS/GRAD’ (Eq. 17). Three covariance functions are presented in Sect. 3.1 and 3.2. The remaining covariance functions for collocation with gradients are as follows:

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_E,\mathrm{ZTD} \right) = \mathbf {C}_{\mathrm{ss}}\left( \mathrm{ZTD},G_E\right) \nonumber \\&\quad = H_{\mathrm{ZTD}} \cdot \frac{\sigma _{\mathrm{signal}}^2}{q^2}\cdot \frac{2(x_{\mathrm{ZTD}}-x_{G_E})}{(\Delta x_0)^2}\cdot e^{-\frac{z_{G_E}+z_{\mathrm{ZTD}}}{2z_0}}, \end{aligned}$$

(18)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_N,\mathrm{ZTD} \right) = \mathbf {C}_{\mathrm{ss}}\left( \mathrm{ZTD},G_N\right) \nonumber \\&\quad = H_{\mathrm{ZTD}} \cdot \frac{\sigma _{\mathrm{signal}}^2}{q^2}\cdot \frac{2(y_{\mathrm{ZTD}}-y_{G_N})}{(\Delta y_0)^2}\cdot e^{-\frac{z_{G_N}+z_{\mathrm{ZTD}}}{2z_0}}, \end{aligned}$$

(19)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_E,G_N \right) = \mathbf {C}_{\mathrm{ss}}\left( G_N,G_E\right) = H_{\mathrm{ZTD}}^2 \nonumber \\&\quad \cdot \left( \frac{-8\sigma _{\mathrm{signal}}^2}{q^3} \right) \cdot \frac{x_{G_N}-x_{G_E}}{(\Delta x_0)^2} \cdot \frac{y_{G_N}-y_{G_E}}{(\Delta y_0)^2}\cdot e^{-\frac{z_{G_E}+z_{G_N}}{z_0}}, \end{aligned}$$

(20)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_{E,k},G_{E,l} \right) = \mathbf {C}_{\mathrm{ss}}\left( G_{E,l},G_{E,k} \right) \nonumber \\&\quad =H_{\mathrm{ZTD}}^2 \cdot \frac{2\sigma _{\mathrm{signal}}^2}{q^2} \cdot \frac{1}{(\Delta x_0)^2} \cdot e^{-\frac{z_{G_{E,k}}+z_{G_{E,l}}}{2z_0}} \nonumber \\&\qquad - \left( \frac{-8\sigma _{\mathrm{signal}}^2}{q^3} \right) \cdot \frac{(x_{G_{E,k}}-x_{G_{E,l}})^2}{(\Delta x_0)^4} \cdot e^{-\frac{z_{G_{E,k}}+z_{G_{E,l}}}{z_0}}, \end{aligned}$$

(21)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_{N,k},G_{N,l} \right) = \mathbf {C}_{\mathrm{ss}}\left( G_{N,l},G_{N,k} \right) \nonumber \\&\quad = H_{\mathrm{ZTD}}^2 \cdot \frac{2\sigma _{\mathrm{signal}}^2}{q^2} \cdot \frac{1}{(\Delta y_0)^2} \cdot e^{-\frac{z_{G_{N,k}}+z_{G_{N,l}}}{2z_0}} \nonumber \\&\qquad - \left( \frac{-8\sigma _{\mathrm{signal}}^2}{q^3} \right) \cdot \frac{(y_{G_{N,k}}-y_{G_{N,l}})^2}{(\Delta y_0)^4} \cdot e^{-\frac{z_{G_{N,k}}+z_{G_{N,l}}}{z_0}}, \end{aligned}$$

(22)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_E,N_{\mathrm{tot}} \right) = \mathbf {C}_{\mathrm{ss}}\left( N_{\mathrm{tot}},G_E \right) =\nonumber \\&\quad =H_{\mathrm{ZTD}} \cdot \frac{8\sigma _{\mathrm{signal}}^2}{q^3} \cdot \frac{(x_{N_{\mathrm{tot}}}-x_{G_E})(z_{N_{\mathrm{tot}}}-z_{G_E})}{(\Delta x_0)^2(\Delta z_0)^2} \nonumber \\&\qquad \cdot e^{-\frac{z_{G_E}+z_{N_{\mathrm{tot}}}}{z_0}} + \frac{\sigma _{\mathrm{signal}}^2(2-q)(x_{N_{\mathrm{tot}}}-x_{G_E})}{q^3z_0(\Delta x_0)^2} \nonumber \\&\qquad \cdot e^{-\frac{z_{G_E}+z_{N_{\mathrm{tot}}}}{2z_0}}, \end{aligned}$$

(23)

$$\begin{aligned}&\mathbf {C}_{\mathrm{ss}}\left( G_N,N_{\mathrm{tot}} \right) = \mathbf {C}_{\mathrm{ss}}\left( N_{\mathrm{tot}},G_N \right) \nonumber \\&\quad =H_{\mathrm{ZTD}} \cdot \frac{8\sigma _{\mathrm{signal}}^2}{q^3} \cdot \frac{(y_{N_{\mathrm{tot}}}-y_{G_N})(z_{N_{\mathrm{tot}}}-z_{G_N})}{(\Delta y_0)^2(\Delta z_0)^2} \nonumber \\&\qquad \cdot e^{-\frac{z_{G_N}+z_{N_{\mathrm{tot}}}}{z_0}} + \frac{\sigma _{\mathrm{signal}}^2(2-q)(y_{N_{\mathrm{tot}}}-y_{G_N})}{q^3z_0(\Delta y_0)^2} \nonumber \\&\qquad \cdot e^{-\frac{z_{G_N}+z_{N_{\mathrm{tot}}}}{2z_0}}. \end{aligned}$$

(24)