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

: 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 coefﬁcients 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 artiﬁcial 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 signiﬁcantly 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 modiﬁed wet MF, which is no longer based on the assumption of a spherically layered atmosphere. We show that with this modiﬁed 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 coefﬁcient of the developed closed-form expression. We provide some evidence that current NWMs are able to do so.


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

Tropospheric Delay
The tropospheric delay T is given through: 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: 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).

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: 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: Here the indices indicate the specific elevation and azimuth angle and m g denotes the gradient MF. m g (e) = 1 sin(e) tan(e) + C 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(e, a) ∼ T 0 (e) + m g (e)·[N cos(a) + E sin(a)] where T 0 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] [N, The entries of the design matrix Γ are given by the partial derivatives of the righthand side of Equation (7) with respect to the gradient coefficients and the weight matrix W reads as: W kl = sin(e k )sin(e l )δ 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.

Tropospheric Delays in the GNSS Analysis
In the GNSS analysis, the tropospheric delay is parameterized T(e, a) ∼ m h (e)·ZHD + m w (e)·ZWD + m g (e)·[N cos(a) + E sin(a)] 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: 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(e; a, b, c) = 1 + a 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.

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 orderionospheric 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 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 0 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 raytraced 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: [δr, ZWD 1 , . . . , ZWD 24 , N 1 , . . . , N 24 , E 1 , . . . , E 24 , δt 1 , . . . , 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: 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.

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

Modified Wet Mapping Function
The tropospheric delay is approximated as: 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

Modified Wet Mapping Function
The tropospheric delay is approximated as: 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 Remote Sens. 2021, 13, 3944 8 of 18 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 0 ; as Z 0 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 0 . Thus, we rewrite the tropospheric delay as: where the modified wet MF is defined as: and the fourth coefficient z w is given by: The coefficient Z 0 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 0 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 1 to Z 4 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 0 . 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.

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

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

Experiment 2
The left panel in Figure 3 shows the error in the ZTD and the error in the station upcomponent 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 upcomponent error. Again, the error in the ZTD reached 15 mm and the error in the upcomponent 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

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 stationspecific 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.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 18 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.

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.

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.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 18 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.

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.

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

Experiment 4
The left panel in Figure 5 shows the error in the ZTD and the error in the station upcomponent 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.

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

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

On the Relation between the Tropospheric Parameter Z0 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 zw is more important than the three wet MF coefficients aw, bw, and cw taken together. It is also important to note that the modified wet MF is reduced to the standard wet MF if the fourth coefficient zw, and hence the tropospheric parameter Z0, vanishes. Thus, it is worth studying the relation between the tropospheric parameter Z0 and the errors in the GNSS estimates in detail. The left panel in Figure 7 shows the tropospheric parameter Z0 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 Z0 are correlated. The error in the ZTD appears to be proportional to Z0. 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 Z0 for all stations and epochs analyzed. With the linear fit (indicated by the red line) we obtain the following relation:

On the Relation between the Tropospheric Parameter Z 0 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 0 , vanishes. Thus, it is worth studying the relation between the tropospheric parameter Z 0 and the errors in the GNSS estimates in detail. The left panel in Figure 7 shows the tropospheric parameter Z 0 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 0 are correlated. The error in the ZTD appears to be proportional to Z 0 . 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 0 for all stations and epochs analyzed. With the linear fit (indicated by the red line) we obtain the following relation: Another question is whether the tropospheric parameter Z 0 introduces a systematic or random error. We find that for the POTS station and the bulk of stations, the tropospheric parameter Z 0 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.
or random error. We find that for the POTS station and the bulk of stations, the tropospheric parameter Z0 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 firstorder 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 Z0. We have some evidence that current NWMs can predict the tropospheric parameter Z0. This is discussed in the next section.

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 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 0 . We have some evidence that current NWMs can predict the tropospheric parameter Z 0 . This is discussed in the next section.

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 0 . 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 0 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 0 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 0 . 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 0 . 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 Z0. 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 Z0 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 Z0 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 Z0. 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 Z0.  correspond to the mean and standard deviation between the GNSS and NWM ZTDs. The numbers with the yellow background correspond to the case when we corrected the GNSS ZTDs a posteriori. For details, refer to the text.

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: 1. 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 upcomponent 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. 2. 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,

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: 1.
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.

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