Impact of Tropospheric Mismodelling in GNSS Precise Point Positioning: A Simulation Study Utilizing Ray-Traced Tropospheric Delays from a High-Resolution NWM

Abstract

In GNSS analysis, the tropospheric delay is parameterized by applying mapping functions (MFs), zenith delays, and tropospheric gradients. Thereby, the wet and hydrostatic MF are derived under the assumption of a spherically layered atmosphere. The coefficients of the closed-form expression are computed utilizing a climatology or numerical weather model (NWM) data. In this study, we analyze the impact of tropospheric mismodelling on estimated parameters in precise point positioning (PPP). To do so, we mimic PPP in an artificial environment, i.e., we make use of a linearized observation equation, where the observed minus modelled term equals ray-traced tropospheric delays from a high-resolution NWM. The estimated parameters (station coordinates, clocks, zenith delays, and tropospheric gradients) are then compared with the known values. The simulation study utilized a cut-off elevation angle of 3° and the standard downweighting of low elevation angle observations. The results are representative of a station located in central Europe and the warm season. In essence, when climatology is utilized in GNSS analysis, the root mean square error (RMSE) of the estimated zenith delay and station up-component equal about 2.9 mm and 5.7 mm, respectively. The error of the GNSS estimates can be reduced significantly if the correct zenith hydrostatic delay and the correct hydrostatic MF are utilized in the GNSS analysis. In this case, the RMSE of the estimated zenith delay and station up-component is reduced to about 2.0 mm and 2.9 mm, respectively. The simulation study revealed that the choice of wet MF, when calculated under the assumption of a spherically layered troposphere, does not matter too much. In essence, when the ‘correct’ wet MF is utilized in the GNSS analysis, the RMSE of the estimated zenith delay and station up-component remain at about 1.8 mm and 2.4 mm, respectively. Finally, as a by-product of the simulation study, we developed a modified wet MF, which is no longer based on the assumption of a spherically layered atmosphere. We show that with this modified wet MF in the GNSS analysis, the RMSE of the estimated zenith delay and station up-component can be reduced to about 0.5 mm and 1.0 mm, respectively. In practice, its success depends on the ability of current (future) NWM to predict the fourth coefficient of the developed closed-form expression. We provide some evidence that current NWMs are able to do so.

1. Introduction

In GNSS analysis the signal travel time delay induced by the neutral atmosphere between the satellite and the station, known as tropospheric delay, is approximated by utilizing mapping functions (MFs), zenith delays, and tropospheric gradient components [1]. The so-called hydrostatic and wet MF (the ratio of slant and zenith delays) are derived under the assumption of a spherically layered atmosphere. In order to take into account the deviation from a spherically layered atmosphere, a dedicated gradient MF and gradient components are introduced into the tropospheric delay model. The parameters of the tropospheric delay model are derived from a climatology or numerical weather model (NWM). In one of the most demanding GNSS applications, known as precise point positioning (PPP) [2], the station coordinates are estimated with cm-level accuracy. To reach this level of accuracy in the positioning domain, the inaccuracy of the underlying climatology or NWM must be taken into account. To do so, some of the parameters in the tropospheric delay model, i.e., the zenith delay and gradient components, are estimated together with the station coordinates (the coefficients of the MFs are kept fixed). Thus, the output of PPP is two products: station coordinates and respective tropospheric parameters. This implies that PPP can be regarded as a tool for atmospheric remote sensing. The zenith delays can be related to the precipitable water vapor (PWV) above the considered station [3], and the gradient components can be roughly related to the (first-order) horizontal PWV gradient at the respective stations [4]. These relations explain the interest of meteorology and climate research in GNSS-based atmospheric data [5].

The GNSS is all-weather; i.e., the station coordinates and the respective tropospheric parameters are obtained under all weather conditions. This is a clear advantage over other atmospheric remote sensing instruments. However, the accuracy of GNSS estimates must be weather dependent to some extent; the reason being the inaccuracies in the tropospheric delay model. First, the functional form of the tropospheric delay model is inaccurate. This can be recognized when ray-traced tropospheric delays, especially those derived from a high-resolution NWM, are compared with parameterized tropospheric delays. Second, let us assume for a moment that the functional form of a tropospheric delays is correct; the parameters of the tropospheric delay model are derived from climatology or a NWM, and neither of them can be regarded as error free. As mentioned before, this error source is mitigated, since the zenith delays and gradient components are estimated together with the station coordinates in the GNSS analysis. Nevertheless, a priori zenith delays and the MFs are important error sources [6]. The most commonly used a priori zenith delay to date is probably the one from the global pressure and temperature model (GPT) [7], and the most commonly used MF to date is probably the global mapping function (GMF) [8]. The drawback of this MF is that it is based on climatology, and thus, limited in its ability to predict short-term variability of the atmosphere. In contrast, the NWMs capture the short-term variability of the atmosphere and, thus, are considered a valuable data source for MFs [9]. The most prominent example of an MF based on NWM data is the Vienna Mapping Functions 1 (VMF1) [10]. Recently, various other MFs based on NWM data have been successfully applied in practice [11,12,13,14]. This includes those MFs based on NWM data that are no longer based on the assumption of a spherically layered atmosphere, such as the Vienna Mapping Functions 3 (VMF3) [15] and variants thereof [16].

The larger the deviation of the climatology or the NWM from the true state of the atmosphere, the larger the errors in the GNSS estimates. This is obvious and suggested in many studies; however, it is not trivial to provide concrete numbers. For example, we may exchange one MF, say an MF based on a climatology, with another MF, say an MF based on an NWM, and study the impact of this exchange on the GNSS estimates. However, this exercise will not provide us with the absolute error of the tropospheric mismodelling of GNSS estimates. The purpose of our study is to estimate such numbers. To do so we mimic PPP in an artificial environment; i.e., we make use of a linearized observation equation, where the observed minus modelled term equals ray-traced tropospheric delays from a high-resolution NWM. The estimated parameters, station coordinates, zenith delays, and tropospheric gradients are then compared with the perfectly known values.

2. Materials and Methods

2.1. Tropospheric Delay

The tropospheric delay T is given through:

$$T={\displaystyle \int }n ds-g$$

(1)

where n denotes the index of refraction, s denotes the arclength of the ray-path, and g denotes the geometric distance. The ray-path follows Fermat’s principle. The index of refraction n is related to the refractivity Ψ through:

$$n={10}^{-6}\Psi +1$$

(2)

The refractivity field depends on the pressure, temperature, and humidity field [17]. The pressure, temperature, and humidity fields are unknown, and, hence, the true refractivity field is also unknown. In this study we assume that the output of a mesoscale NWM gives the true refractivity field. Specifically, we utilize the weather research and forecasting (WRF) model [18] to simulate the refractivity field of the atmosphere. We consider a limited area covering the central part of Europe (for details see below). The initial and boundary conditions for the limited area were the 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 models were applied: the Thompson scheme for the microphysics [19], the Kain–Fritsch scheme for the cumulus parameterization [20], the Yonsei University scheme for the planetary boundary layer [21], the RRTMG Short and Longwave scheme for the radiation [22], the Unified Noah Land Surface Model for the land surface [23], and the revised MM5 scheme for the surface layer [24]. The refractivity was calculated from the pressure, temperature, and humidity and 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 [25]. Then the algorithm by [26] allowed us to compute tropospheric delays with high speed and precision for any station–satellite link.

The artificial environment that we made use of in the following simulations has some limitations. We ran the WRF model with a limited horizontal resolution of 10 km. The simulation would be made more realistic by increasing the horizontal resolution. However, even with an horizontal resolution of say 2 km the mesoscale model would not resolve all tropospheric features. In particular, the mesoscale model does not explicitly resolve turbulence in the planetary boundary layer (the turbulence is parameterized).

2.2. Zenith Delays and Tropospheric Gradient

The tropospheric delay can be regarded as a function of the elevation angle e and the azimuth angle a. Thus, the tropospheric delay can be written as:

$$T=T\left(e,a\right)$$

(3)

For any station location, we consider k = 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°. With this set of tropospheric delays we are prepared to define the ZTD and the gradient components. Both the ZTD and the gradient components can be understood as specific linear combinations of tropospheric delays. Specifically, the ZTD is defined as:

$$ZTD=T\left({90}^{°},0^{°}\right)$$

(4)

The north gradient component N and the east gradient component E are defined as:

$$\begin{gathered} N=\frac{{\sum }_{i=1}^{\mathit{k}}{m}_g\left(e_i\right)sin{\left(e_i\right)}^{2}cos\left(a_i\right)T\left(e_i,a_i\right)}{{\sum }_{i=1}^{\mathit{k}}{m}_g{\left(e_i\right)}^{2}sin{\left(e_i\right)}^{2}cos{\left(a_i\right)}^2} \\ E=\frac{{\sum }_{i=1}^{\mathit{k}}{m}_g\left(e_i\right)sin{\left(e_i\right)}^{2}sin\left(a_i\right)T\left(e_i,a_i\right)}{{\sum }_{i=1}^{\mathit{k}}{m}_g{\left(e_i\right)}^{2}sin{\left(e_i\right)}^{2}sin{\left(a_i\right)}^2} \end{gathered}$$

(5)

Here the indices indicate the specific elevation and azimuth angle and m_g denotes the gradient MF.

$$m_g\left(e\right)=\frac{1}{sin\left(e\right)\tan \left(e\right)+C}$$

(6)

where C = 0.0031 [1].

The ZTD is defined as the tropospheric delay in the zenith direction and no further explanation is required. The definition for the north and east gradient components deserves some further explanation. To understand this definition let us assume that the tropospheric delay is approximated by:

$$T\left(e,a\right) ~ T_0\left(e\right)+m_g\left(e\right)·\left[N cos\left(a\right)+E sin\left(a\right)\right]$$

(7)

where T₀ denotes the tropospheric delay under the assumption of a spherically layered atmosphere. Given the set of tropospheric delays the gradient components are determined by a weighted least square fit [27]

$$\left[N,E\right]={\left({\mathbf{\Gamma }}^{\mathit{T}}\mathit{W}\mathbf{\Gamma }\right)}^{-1}{\mathbf{\Gamma }}^{\mathit{T}}\mathit{W} \left(T-T_0\right)$$

(8)

The entries of the design matrix Γ are given by the partial derivatives of the right- hand side of Equation (7) with respect to the gradient coefficients and the weight matrix W reads as:

$$W_{kl}=sin\left(e_k\right)sin\left(e_l\right){\delta }_{kl}$$

(9)

where the indices denote the specific elevation angles. After some algebra it is not difficult to verify that for the set of tropospheric delays Equation (8) reduces to Equation (5). For the details, the reader is referred to [28].

The definition of the ZTD, and in particular the definition of the tropospheric gradient components, appears arbitrary. However, as we will demonstrate in the following, these quantities can be estimated with a high accuracy using the GNSS.

2.3. Tropospheric Delays in the GNSS Analysis

In the GNSS analysis, the tropospheric delay is parameterized

$$T\left(e,a\right)~m_h\left(e\right)·ZHD+m_w\left(e\right)·ZWD+m_g\left(e\right)·\left[Ncos\left(a\right)+Esin\left(a\right)\right]$$

(10)

where ZHD denotes the Zenith Hydrostatic Delay, ZWD denotes the Zenith Wet Delay, m_h denotes the hydrostatic MF, and m_w denotes the wet MF. The ZTD is given by:

$$ZTD=ZHD+ZWD$$

(11)

Typically, the hydrostatic and wet MF are derived under the assumption of a spherically layered atmosphere. Under this assumption, it is recognized that the elevation angle dependency of the hydrostatic and wet MF is accurately described by the continued fraction form [29]

$$m\left(e;a,b,c\right)=\frac{1+\frac{a}{1+\frac{b}{1+c}}}{sin\left(e\right)+\frac{a}{sin\left(e\right)+\frac{b}{sin\left(e\right)+c}}}$$

(12)

where a, b, and c are called MF coefficients. The MF coefficients are computed utilizing a climatology or an NWM. Specifically, for some refractivity profiles (the assumption is that the refractivity is a function of the altitude only) the ratios of slant and zenith delays (mapping factors) are computed (the elevation angles are 3°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, and 90°) and the MF coefficients are determined by least-squares fitting. We can confirm that under the assumption of a spherically layered atmosphere the continued fraction form with three MF coefficients yields an exquisite level of precision. We call our realization of the MF, where all three MF coefficients are determined by least-squares fitting, the Potsdam mapping function (PMF). The most prominent example for a MF based on a climatology is the new mapping function (NMF) [29]. The most prominent example of a MF based on NWM data is the VMF1 [10]. The VMF1 is of particular interest because it is based on an efficient concept; the MF coefficients b and c come from a climatology, and the MF coefficient a is determined from a single mapping factor by inverting the continued fraction form. This is also why several other MFs based on the VMF1 concept exist, e.g., the UNB-VMF1 [30] and the GFZ-VMF1 [31]. The main difference with the original VMF1 is that they are based on different NWMs.

In all the above mentioned concepts, the assumption is that the atmosphere is spherically layered. Recently, attempts have been made to move away from this concept. For example, instead of utilizing a single refractivity profile above the station, we utilize the refractivity field above the station and one can compute 120 mapping factors for various elevation and azimuth angles (the elevation angles are 3°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, and 90° and the spacing in azimuth is 30°). Then, by averaging over the azimuth angles, 10 mapping factors are computed, and the three coefficients of the MF are determined by least-squares fitting. However, some caution is required here as the continued fraction form no longer gives a perfect fit with the mapping factors. The continued fraction form gives a perfect fit with the mapping factors provided that they are calculated under the assumption of a spherically layered atmosphere. This was also recognized in the recently developed VMF3 [15], and variants thereof [16]. In the present work, as a by-product of the simulation study, we will also develop a new MF concept. The details will be provided below.

2.4. Precise Point Positioning Simulation

We simulated PPP in an artificial environment, i.e., the refractivity field of the mesoscale NWM. The simulator that we implemented is very similar to the simulators utilized to study various other effects, such as multipath and geometry effects [32] or higher order-ionospheric effects [33]. The simulation was simplified by ignoring carrier-phase ambiguities in the observation equation. In essence, we utilized the following version of the linearized observation equation

$$T\left(e,a\right)-m_h\left(e\right)·ZHD=-\mathit{u}\left(e,a\right)·\delta \mathit{r}+c_0·\delta t+m_w\left(e\right)·ZWD+m_g\left(e\right)·\left[Ncos\left(a\right)+Esin\left(a\right)\right]$$

(13)

The left hand side of the equation represents the measured minus modelled term and the right hand side of the equation includes the parameters to be estimated. Here u denotes the tangent-unit vector of the station–satellite link, δr denotes the coordinate residual, δt denotes the clock residual, and c₀ denotes the vacuum speed of light. In the present study we consider 120 station satellite links for a single epoch, where azimuth angles are selected randomly and elevation angles are obtained through e = 90° − 87° √w, where w ∈ [0,1] is obtained from a random number generator. The set of station–satellite links mimics a realistic observation geometry with a cut-off elevation angle of 3°. The observation geometry is realistic in so far as the elevation angle dependency reflects the linearly increasing number of observations with decreasing elevation angles [25]. However, the simple azimuth angle dependency does not reflect the lack of observations to the North (South) for a station located in the Northern (Southern) Hemisphere. This might have an impact on the results, in particular the tropospheric gradients, and must be studied in future. The tropospheric delays that enter the left hand side of the equation are ray-traced tropospheric delays. For each day, we consider 24 epochs, i.e., every hour, we consider 120 station–satellite links. Then, we combine the 24 times 120 equations and obtain, using the least-square adjustment, the coordinate residual on a daily basis and the clock- and tropospheric parameter residuals epoch-wise. Standard elevation angle dependent weighting is applied in the least-square fit. In short, let T denote the ray-traced tropospheric delays and A denote the a priori tropospheric delays, i.e., the product of the hydrostatic MF and the ZHD, then the solution, which includes the coordinate residual (daily) and the clock- and tropospheric parameter residuals (epoch wise), is obtained through:

$$\left[\delta \mathit{r},ZW{D}_1,\dots ,ZW{D}_{24},N_1,\dots ,N_{24},E_1,\dots ,E_{24},\delta t_1,\dots ,\delta t_{24}\right]={\left({\mathit{Q}}^{\mathit{T}}\mathit{W}\mathit{Q}\right)}^{-1}{\mathit{Q}}^{\mathit{T}}\mathit{W} \left(T-A\right)$$

(14)

The design matrix Q is defined by the partial derivatives of the right-hand side of Equation (13), with respect to the coordinate residual, the clock residual, the ZWD, and the gradient components. The weight matrix equals the weight matrix defined in Equation (9) and, therefore, downweights the low elevation angle observations. One of the reasons to do so in practice is the known deficiencies of the tropospheric model at low elevation angles. The final ZTD is given by the sum of the a priori ZHD and the estimated ZWD.

Next, let ZTD_WRF denote the ZTD calculated from the NWM, following Equation (4), and let N_WRF and E_WRF denote the north and east gradient component calculated from the NWM, following Equation (5); then, the errors of the GNSS estimates are given by:

$$\begin{gathered} \Delta ZTD=\left(ZHD+ZWD\right)-ZT{D}_{WRF} \\ \Delta N=N-N_{WRF} \\ \Delta E=E-E_{WRF} \\ \Delta \mathit{r}=\delta \mathit{r} \\ \Delta t=\delta t \end{gathered}$$

(15)

and quantify the impact of the tropospheric mismodelling in the GNSS analysis. Note that ZTD_WRF, N_WRF, and E_WRF stand for the true ZTD, the true north, and the true east gradient component, since in the context of our simulation study the NWM represents the true state of the atmosphere.

2.5. Experiment Setup

The error of the GNSS estimates depends on the chosen ZHD, the chosen hydrostatic MF, and the chosen wet MF. In total, we ran five experiments, as summarized in

Table 1. Summary of the simulation experiments. The simulation experiments differed in the chosen ZHD, hydrostatic and wet MF. GMF stands for the global mapping function, GPT stands for the global pressure and temperature model, and PMF stands for the Potsdam mapping function, i.e., the MF calculated from the NWM. The subscripts h and w indicate hydrostatic and wet MF, respectively. The subscript z indicates modified wet MF. For details refer to the text.

Experiment	ZHD	mh	mw
#1	GPT	GMFh	GMFw
#2	NWM	GMFh	GMFw
#3	NWM	PMFh	GMFw
#4	NWM	PMFh	PMFw
#5	NWM	PMFh	PMFz

. In the first experiment, the ZHD came from the GPT [7] and the hydrostatic and the wet MF were taken from the GMF [8]. In essence all tropospheric parameters that enter the observation equation came from climatology. This approach is the most widely used in practice, due to its simplicity. There are no external data required in this approach. In the second experiment, the ZHD came from the NWM, but the hydrostatic and wet MF were based on the climatology. This approach is also considered practical, as it solely requires the pressure at the station, due to an existing relationship between the pressure and ZHD [34]. In the third experiment, the ZHD and the hydrostatic MF were derived from the NWM, whereas the wet MF was still based on the climatology. In the fourth experiment the ZHD, the hydrostatic and the wet MF were derived from the NWM. It is important to note that up to this point, both the hydrostatic and the wet MF were calculated under the assumption of a spherically layered atmosphere. Finally, in the fifth experiment, the ZHD, the hydrostatic MF and the wet MF were again derived from the NWM. The difference, however, was that the wet MF was no longer calculated on the assumption of a spherically layered atmosphere. The details on the derivation of this modified wet MF are provided in the next section.

In the simulation study we considered a limited area covering Germany, the Czech Republic, and parts of Poland and Austria, and a period of two months (May and June) in 2013. We utilized data from more than 400 stations. The locations of the stations are shown in Figure 1. This choice was motivated by an existing dense station network, which was utilized in the Benchmark campaign within the European COST Action ES1206 GNSS4SWEC (Advanced GNSS tropospheric products for monitoring severe weather and climate) [35]. The station locations and time period considered in the simulation study represent stations located in central Europe and the warm season. The station POTS (Potsdam, Germany), located in the center of the limited area model, was considered representative of the bulk of the stations. Thus, we will provide some more details for this station. For the rest of the stations we solely provide some statistics.

2.6. Modified Wet Mapping Function

The tropospheric delay is approximated as:

$$\begin{gathered} T\left(e,a\right)~m_h\left(e;a_h,b_h,c_h\right)·ZHD+m_w\left(e;a_w,b_w,c_w\right)·ZWD+\dots \\ {m}_g\left(e\right)·Z_0+\dots \\ {m}_g\left(e\right)·\left[Ncos\left(a\right)+Esin\left(a\right)\right]+\dots \\ {m}_g\left(e\right)·\left[Z_{1}cos\left(2·a\right)+Z_{2}sin\left(2·a\right)\right]+\dots \\ {m}_g\left(e\right)·\left[Z_{3}cos\left(3·a\right)+Z_{4}sin\left(3·a\right)\right]+\dots \end{gathered}$$

(16)

where the hydrostatic and wet MF are calculated as usual, under the assumption of a spherically layered atmosphere and the coefficients Z can be understood as higher-order tropospheric parameters. The proposed approximation for the tropospheric delay can be considered as follows: The difference between tropospheric delays and tropospheric delays calculated under the assumption of a spherically layered atmosphere can be regarded as a function of the elevation and azimuth angle. Given a set of tropospheric delays and a set of tropospheric delays calculated under the assumption of a spherically layered atmosphere (the elevation angles are 3°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, and 90° and the spacing in azimuth is 30°), then, for a specific elevation angle, the difference between the tropospheric delays can be expanded in a Fourier series. In general, the coefficients of the Fourier series expansion will be different for different elevation angles. Suppose we assume that the elevation angle dependency of these coefficients of the Fourier series expansion follows a simple law, namely that this elevation angle dependency can be represented by the gradient MF. In that case, we end up with the expression for the tropospheric delay provided above. In this interpretation, the north and east gradient components can be considered the second and third coefficients of the Fourier series expansion. The Z coefficients are discarded in the GNSS analysis. Obviously, the danger in the estimation of the ZTD, station up-component, and clock lies in the presence of Z₀; as Z₀ appears in a term that depends on the elevation angle only, and this term is present on the left-hand side of the observation equation, but there is no corresponding term on the right-hand side of the observation equation, and it will be absorbed by the estimated ZTD, station up-component, and clock (also see the linearized observation equation). One possible solution is to introduce a modified wet MF, which takes into account the coefficient Z₀. Thus, we rewrite the tropospheric delay as:

$$\begin{gathered} T\left(e,a\right)~mh\left(e;a_h,b_h,c_h\right)·ZHD+m_w^{*}\left(e;a_w,b_w,c_w,z_w\right)·ZWD+\dots \\ {m}_g\left(e\right)·\left[Ncos\left(a\right)+Esin\left(a\right)\right]+\dots \\ {m}_g\left(e\right)·\left[Z_{1}cos\left(2·a\right)+Z_{2}sin\left(2·a\right)\right]+\dots \\ {m}_g\left(e\right)·\left[Z_{3}cos\left(3·a\right)+Z_{4}sin\left(3·a\right)\right]+\dots \end{gathered}$$

(17)

where the modified wet MF is defined as:

$$m_w^{*}\left(e;a_w,b_w,c_w,z_w\right)=m_w\left(e;a_w,b_w,c_w\right)+m_g\left(e\right)·z_w$$

(18)

and the fourth coefficient z_w is given by:

$$z_w=\frac{Z_0}{ZWD}$$

(19)

The coefficient Z₀ is determined by a least-square fit. For the determination of the three coefficients a_w, b_w, and c_w, solely the refractivity profile above the station and 10 mapping factors are required, whereas for the determination of the fourth coefficient z_w the refractivity field above the station and 120 additional mapping factors are required. Clearly, this makes the determination of the modified wet MF more expensive. Adding the term containing Z₀ to the wet MF, and not to the hydrostatic MF, is motivated by the origin of this extra term, which is more likely to be the wet than the hydrostatic refractivity field. We make no attempt to estimate Z₁ to Z₄ in the GNSS analysis. We think that with the increasing number of parameters to be estimated besides the coordinates and clocks, the GNSS solution would probably get worse.

Equation (16) appears similar to the approximation proposed by [36], except for the extra term containing Z₀. This can be explained by the fact that in [36] the isotropic part of the tropospheric delay, i.e., the part of the tropospheric delay which depends solely on the elevation angle, already contains the average over various azimuth angles. The isotropic part in [36] is not calculated under the assumption of a spherically layered atmosphere. The hydrostatic and wet MF which enter our Equation (16) are calculated under the assumption of a spherically layered atmosphere.

3. Results

In this section we show the results of the simulation experiments listed in Table 1. The estimation of the station clock is part of the PPP simulation. However, we will only show and discuss the impact of the tropospheric mismodelling on the estimated station coordinates and tropospheric parameters, as they are the parameters of interest in precise positioning and atmospheric remote sensing.

3.1. Experiment 1

The error in the ZTD and the error in the station up-component for the station POTS as a function of time are shown in the left panel of Figure 2. The error in the ZTD reached 15 mm, and the error in the up-component reached 10 mm. The time series shows the known correlation between the ZTD and station up-component error. The step-like structure in the time series for the error of the station up-component is due to the fact that the coordinate residual was estimated on a daily basis. The noise-like structure in the time series for the error in the ZTD is due to the zenith delay residual being estimated epoch-wise. The error in the ZTD is the composite of two error sources; the climatology is not able to follow either the slow change in the hydrostatic portion of the refractivity field or the rapid change in the wet portion of the refractivity field. The right panel in Figure 2 shows the station-specific root mean square error (RMSE) for the ZTD, the gradient components, and the station up-component. With only a few exceptions, the RMSE is the same for the considered stations and, hence, it is meaningful to provide the average RMSE for the respective parameters. This is also why we considered the POTS station to represent the bulk of the stations. We obtained 2.9 mm for the ZTD, 0.11 mm for the north (east) gradient component, and 5.7 mm for the station up-component. The significant RMSE for the ZTD and the station up-component was owing to the hydrostatic MF; the wet MF and the ZHD are inaccurate because they are based on a climatology. The RMSE for the gradient components is not regarded as significant, as a deviation of 0.1 mm in the gradient component converts into a tropospheric delay deviation of about 10 mm at an elevation angle as low as 5°. The small deviations can be explained by the fact that the way the gradient components are defined and the way gradient components are estimated in the GNSS analysis are very similar. The remaining deviations were mainly due to the different geometries in their determination.

3.2. Experiment 2

The left panel in Figure 3 shows the error in the ZTD and the error in the station up-component for the station POTS as a function of time. It is not obvious that the introduction of the ZHD from the NWM yielded a significant reduction in the ZTD and station up-component error. Again, the error in the ZTD reached 15 mm and the error in the up-component reached 10 mm. However, a closer look reveals that the introduction of the ZHD from the NWM yielded a somewhat smaller deviation around zero. Due to the ZHD from the NWM, the slow change in the hydrostatic portion of the refractivity field was, to some extent, taken into account. The right panel in Figure 3 shows the station-specific RMSE for the ZTD, the gradient components, and the station up-component; we obtained 2.7 mm for the ZTD, 0.11 mm for the north (east) gradient component, and 5.2 mm for the station up-component. We can state that introducing the ZHD from the NWM yielded a small but consistent reduction of the error in the ZTD and station up-component. Regarding the gradient components, it appears that the chosen ZHD did not have an impact, since if we replace the ZHD from the climatology with the ZHD from the NWM the error for the gradient components remains unchanged.

3.3. Experiment 3

The error in the ZTD and the error in the station up-component for the POTS station as a function of time are shown in the left panel of Figure 4. The introduction of the ZHD and hydrostatic MF from the NWM yielded a significant reduction in the ZTD and station up-component error. The error in the station up-component remained below 6 mm and the error in the ZTD was typically below 10 mm. A large ZTD error, of about 15 mm around day 30, is still visible in the time series. Since the ZHD and hydrostatic MF come from the NWM, the slow change in the hydrostatic portion of the refractivity field was taken into account. The wet MF came from the climatology, and, hence, the rapid change in the wet portion of the refractivity field was not considered. The right panel in Figure 4 shows the station-specific RMSE for the ZTD, the gradient components, and the station up-component; we obtained 2.0 mm for the ZTD, 0.1 mm for the north (east) gradient component, and 2.9 mm for the station up-component. Regarding the gradient component, it appears that the chosen hydrostatic MF was not significant, as replacing the hydrostatic MF from the climatology with the hydrostatic MF from the NWM resulted in a difference in the error for the gradient components of only 0.01 mm.

3.4. Experiment 4

The left panel in Figure 5 shows the error in the ZTD and the error in the station up-component for the POTS station as a function of time. The introduction of the ZHD, hydrostatic, and wet MF from the NWM significantly reduced the ZTD and station up-component error. The error in the station up-component remained below 6 mm, and the error in the ZTD was typically below 10 mm. Again, a large ZTD error of about 15 mm around day 30 is visible in the time series. The ZTD and station up-component error reduction mainly came from the introduction of the NWM, ZHD, and NWM hydrostatic MF. The introduction of the NWM wet MF brought almost no reduction in the ZTD and station up-component error, and we can state that the errors of the GNSS estimates were insensitive to the chosen wet MF. However, we must be more specific, and better yet, we state that the errors of the GNSS estimates were insensitive to the chosen wet MF, as long as it was calculated under the assumption of a spherically layered atmosphere. The right panel in Figure 5 shows the station-specific RMSE for the ZTD, the gradient components, and the station up-component; we obtained 1.8 mm for the ZTD, 0.1 mm for the north (east) gradient component, and 2.4 mm for the station up-component. The statistics confirm that as long as the wet MF is calculated under the assumption of a spherically layered atmosphere, it does not have too much impact. In essence, if we replace the wet MF, which is based on the climatology, with the wet MF, which is based on the NWM, the RMSE is reduced by only 0.2 mm for the ZTD and 0.5 mm for the station up-component. Regarding the gradient components, it appears that the chosen wet MF did not have an impact. If we replace the wet MF from the climatology with the wet MF from the NWM, the error for the gradient components remains unchanged.

3.5. Experiment 5

The error in the ZTD and the error in the station up-component for the POTS station as a function of time are shown in the left panel of Figure 6. The introduction of the ZHD, hydrostatic, and modified wet MF from the NWM yielded a significant reduction in the ZTD and station up-component error. The error in the station up-component remained below 3 mm, and the error in the ZTD was typically below 2 mm. Since the ZHD and the hydrostatic MF came from the NWM, the slow change in the hydrostatic portion of the refractivity field was taken into account. In addition, the modified wet MF came from the NWM, so that the rapid change in the wet portion of the refractivity field was taken into account as well. In the previous section, we assumed that the errors of the GNSS estimates are not very sensitive to the chosen wet MF when the wet MF is calculated under the assumption of a spherically layered atmosphere. Here, we support this assertation; provided that the wet MF is no longer calculated under the assumption of a spherically layered atmosphere, the errors of the GNSS estimates can be reduced significantly. The right panel in Figure 6 shows the station-specific RMSE for the ZTD, the gradient components and the station up-component; we obtained 0.5 mm for the ZTD, 0.09 mm for the north (east) gradient component, and 1 mm for the station up-component. Regarding the gradient components, it appears that the chosen ZHD, hydrostatic, and modified wet MF from the NWM did not have a significant impact.

4. Discussion

4.1. On the Relation between the Tropospheric Parameter Z₀ and the Errors in the GNSS Estimates

The GNSS estimates are not very sensitive to the wet MF when the wet MF is calculated under the assumption of a spherically layered atmosphere. In other words, it does not matter whether the three coefficients of the wet MF are derived from the climatology or the NWM. The situation changes if the wet MF is not calculated under the assumption of a spherically layered atmosphere. The fourth wet MF coefficient z_w is more important than the three wet MF coefficients a_w, b_w, and c_w taken together. It is also important to note that the modified wet MF is reduced to the standard wet MF if the fourth coefficient z_w, and hence the tropospheric parameter Z₀, vanishes. Thus, it is worth studying the relation between the tropospheric parameter Z₀ and the errors in the GNSS estimates in detail. The left panel in Figure 7 shows the tropospheric parameter Z₀ for the Potsdam station as a function of time. A one-to-one comparison of the left panel in Figure 7 and the left panel in Figure 5 reveals that the error in the ZTD and the appearance of Z₀ are correlated. The error in the ZTD appears to be proportional to Z₀. This is confirmed when we take a look at all the stations and epochs analyzed. The right panel in Figure 7 shows the error in the ZTD as obtained from the fourth simulation experiment versus the tropospheric parameter Z₀ for all stations and epochs analyzed. With the linear fit (indicated by the red line) we obtain the following relation:

$$\Delta ZTD ~ 6.5·Z_0$$

(20)

Another question is whether the tropospheric parameter Z₀ introduces a systematic or random error. We find that for the POTS station and the bulk of stations, the tropospheric parameter Z₀ introduces a random error into the GNSS estimates. There are some stations where small systematic errors were introduced and those stations are located in complex terrain. Our interpretation is that in complex terrain, the topography disturbs the stratification of the atmosphere systematically, and this cannot be described with the first-order gradient in the PWV field alone, and, therefore, the second-order gradient in the PWV field is required. For example, the moist atmospheric boundary layer follows, to some extent, the topography, so that for a station located on a mountain top or in a valley the atmosphere (on average) will not appear spherically layered. From the station’s perspective, the surrounding PWV field looks either convex or concave (the second derivative in the PWV does not vanish). However, some caution is needed with this interpretation, as the limited horizontal resolution of 10 km that we utilized does not draw conclusive results for stations in complex terrain. Those stations should be analyzed case by case with high-resolution weather model data.

In short, the simulation study suggests that we have two options to reduce errors in the GNSS ZTD estimates: (1) replace the standard wet MF with the modified wet MF in the GNSS analysis, or (2) correct the estimated ZTD a posteriori by utilizing the correction provided in Equation (20). Clearly, in practice, success will depend on the ability of the current (future) NWM to predict the additional tropospheric parameter Z₀. We have some evidence that current NWMs can predict the tropospheric parameter Z₀. This is discussed in the next section.

4.2. Preliminary Results from GNSS Meteorology in Central Europe

In Europe the EUMETNET EIG GNSS water vapor programme (E-GVAP) is in charge of collecting operational GNSS tropospheric products for numerical weather prediction. GFZ is one of the operational E-GVAP Analysis Centers and processes more than 500 stations in near real-time (NRT). About 300 stations are located in central Europe (the area of interest). The Earth Parameter and Orbit determination System (EPOS) software, developed at the GFZ, is used to estimate tropospheric products from the GNSS carrier phase and code measurements in PPP mode [37]. The precise satellite orbits and clocks, as well as earth rotation parameters, are available from the International GNSS Service (IGS) analysis center at the GFZ. The station coordinates are estimated in a sliding window mode (24 h), station clock errors are estimated epoch by epoch, and ZTDs and tropospheric gradients are estimated every 15 min. The a prior ZHD comes from the GPT, and the hydrostatic and the wet MF is taken from the GMF. For details, the reader is referred to [38]. Notably, the ZTDs provided by GFZ are used by several European weather services, such as the Met Office, the United Kingdom’s national weather service; Météo-France, the French national meteorological service; and DWD, the German Weather service, for their day-by-day weather forecasts. Thus, the high quality of the tropospheric products must be ensured. As a part of the quality control at the GFZ, the tropospheric products are compared daily to the corresponding tropospheric parameters derived from the GFS of the NCEP. We utilized short range forecasts (the forecast length ranges from 6 to 11 h), available four times per day, to obtain refractivity fields valid for every hour. These refractivity fields are available with a horizontal resolution of 0.25°. Given these refractivity fields, we computed the station-specific ZHDs, ZWDs, the three hydrostatic (wet) MF coefficients, the two tropospheric gradient components, and, recently, also the tropospheric parameter Z₀. As an example, Figure 8 shows a one-to-one comparison of GNSS and NWM ZWDs for the area of interest on the 15 May 2021 at 11 UTC. An inspection by eye indicates that the GNSS and NWM ZWDs agree fairly well. We can also measure this agreement; e.g., the mean and standard deviation between the GNSS and NWM ZWDs. We find that these mean and standard deviations equal −5.0 mm and 6.2 mm, respectively. The deviations are a composite of the error of the GNSS ZWDs and the error of the NWM ZWDs. For example, we think that the origin of the negative bias is the NWM ZWDs [35]. The GNSS ZWDs include errors caused by the tropospheric mismodelling. For example, our simulation study suggested that the error in the GNSS ZWD caused by ignoring the tropospheric parameter Z₀ in the GNSS analysis is about 1 mm. Hence, we corrected the GNSS ZWDs a posteriori, utilizing the correction provided in Equation (20) and found that the mean and standard deviations equal −5.1 mm and 5.8 mm, respectively. The mean deviation is hardly affected. As mentioned before, the tropospheric parameter Z₀ introduces a random, and not systematic, error. The interesting point is the reduction in the standard deviation, which is small; however, for this particular epoch, it reached several percent. In order to show that the reduction in the standard deviation was not accidental, we show in the upper panel of Figure 9 the time series of the reduction in the standard deviation with a temporal resolution of 1 h for the month of May in 2021. The standard deviations are reduced by a small amount, but they are reduced systematically. For the considered period we found a reduction of the standard deviation for every hour. On average, the reduction in the standard deviation is about 5%. The lower panel of Figure 9 puts some weight on this statistic, as it shows the sample number that enters the statistic. There are some epochs with missing data, due to the NRT processing of the GNSS data, but in general the sample number per epoch is above 200. In conclusion, the small but systematic reduction in the standard deviation gives some evidence that current NWMs can predict the tropospheric parameter Z₀. This is not too surprising. A one-to-one comparison in Figure 8 shows that the NWM can predict the ZWD and, hence, the water vapor field fairly well; therefore, we have good reason to believe that the NWM can predict the tropospheric parameter Z₀.

5. Conclusions

We analyzed the impact of tropospheric mismodelling on the estimated parameters in PPP: station coordinates, clocks, zenith delays, and tropospheric gradients. The true state of the atmosphere is unknown; hence, we performed a simulation study. This was done by mimicking PPP in an artificial environment; i.e., we used a linearized observation equation, where the observed minus modelled term equals the tropospheric delays derived from a NWM. The model that we utilized had a horizontal resolution of 10 km and was thus on the edge between meso-beta and meso-gamma scale models. Thus, it can resolve some, but not all, small-scale atmospheric features. A more realistic description of the true state of the atmosphere is possible by increasing the horizontal resolution, and this is what we suggest for future studies. We utilized a limited area and time; i.e., our results are representative for central Europe and the warm season. Different setups, in particular in the tropics, are recommended for future simulation experiments. Stations in complex terrain should be carefully analyzed, with the respective high-resolution NWM data. Having these limitations in mind, we can draw the following conclusions:

• The quality of GNSS estimates is weather dependent. The reason being that the tropospheric delay model is inaccurate. This tropospheric delay model is based on a climatology or a NWM, and neither of these can be regarded as error free. The larger the deviation of the climatology or the NWM from the true state of the atmosphere, the larger the errors in the GNSS estimates. This is obvious and claimed in many studies; however, it is not trivial to provide some numbers. We provide such numbers. In essence, for the considered area and period, when the climatology was utilized in the tropospheric delay model, the error in the estimated ZTD and station up-component was about 2.9 mm and 5.7 mm, respectively. This error must be understood in a statistical sense; the individual errors can be roughly five times larger.
• The error in the GNSS estimates can be reduced significantly if the true ZHD and the true hydrostatic MF are utilized in the GNSS analysis; in this case the error in the estimated ZTD and station up-component will be reduced to about 2.0 mm and 2.9 mm respectively. The true ZHD and the true hydrostatic MF are unknown. However, there is good reason to believe that the ZHD and the hydrostatic MF from a NWM (analysis or short range forecast) are close to the true ZHD and the true hydrostatic MF. This also explains the success of the NWM based MFs that are currently recommended for analyzing space geodetic data [10]. An indication of the high quality of the NWMs in this respect is the agreement between them. For example, the two state of the art NWM based tropospheric delay models, VMF1 and UNB-VMF1, gave almost the same results in a geodetic analysis [14]. There were some differences for specific locations and times, but a clear advantage of one NWM over the other was not demonstrated in a statistical sense.
• The simulation study showed that the wet MF, when calculated under the assumption of a spherically layered troposphere, did not matter too much. In essence, when the true wet MF (calculated from the true refractivity profile above the station in question) is utilized in the GNSS analysis, the error in the estimated ZTD and station up-component remains about 1.8 mm and 2.4 mm, respectively. This is somewhat surprising; however, it is in line with the results from [9]. This study demonstrates that the improvement that may be realized for the wet MF is significantly less than for the hydrostatic MF. Reference [9] stated that: “… Contrary to the expectation that water vapor is the major source of error, the height error is dominated by the hydrostatic component, except possibly near the equator. Thus it is clear that an improvement in the hydrostatic component will have the larger impact on improving the measurement of height and atmosphere delay by VLBI and GPS”. However, the assumption is that the atmosphere is spherically layered. This brings us to the next item.
• We developed a modified wet MF, which is no longer based on the assumption of a spherically layered atmosphere. We have shown that with this modified wet MF the error in the estimated ZTD and station up-component can be reduced to about 0.5 mm and 1.0 mm, respectively. Its success, in practice, depends on the current (future) NWMs predicting the fourth coefficient of the developed closed form expression. We provided some evidence that current NWMs are able to do so. The determination of the modified wet MF is more expensive than the standard wet MF; the determination of the standard wet MF requires 10 tropospheric delays, whereas the determination of the modified wet MF requires an additional 120 tropospheric delays. The algorithm by [25], which calculates tropospheric delays with high speed and precision on an off-the-shelf PC also makes this feasible for real-time applications. We routinely calculate the additional 120 tropospheric delays, because they are required to determine the tropospheric gradients. The modified wet MF is a by-product of the simulation study (the primary purpose of our study was to quantify the tropospheric delay mismodelling). A more sophisticated MF concept could certainly be developed. In this respect, the newly developed VMF3 [15] sets a new standard, and we are going to analyze if we should follow this concept or refine our current concept in the future.