Estimating the impact of Global Navigation Satellite System horizontal delay gradients in variational data assimilation

remote sensing Article Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation Florian Zus 1,*, Jan Douša 2, Michal Kaˇcmaˇrík3, Pavel Václavovic 2, Galina Dick 1 and Jens Wickert 1,4 1GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany; [email protected] (G.D.); [email protected] (J.W.) 2Geodetic Observatory Pecný, Research Institute of Geodesy, Topography and Cartography, Zdiby 250 66, Czech Republic; [email protected] (J.D.); [email protected] (P.V.) 3Department of Geoinformatics, VŠB – Technical University of Ostrava, Ostrava 708 33, Czech Republic; [email protected] 4Institute of Geodesy and Geoinformation Science, Technical University of Berlin, 10623 Berlin, Germany *Correspondence: [email protected]; Tel.: +49-331-288-1969 Received: 26 November 2018; Accepted: 20 December 2018; Published: 28 December 2018   Abstract: We developed operators to assimilate Global Navigation Satellite System (GNSS) Zenith Total Delays (ZTDs) and horizontal delay gradients into a numerical weather model. In this study we experiment with refractivity fields derived from the Global Forecast System (GFS) available with a horizontal resolution of 0.5 ◦ . We begin our investigations with simulated observations. In essence, we extract the tropospheric parameters from the GFS analysis, add noise to mimic observation errors and assimilate the simulated observations into the GFS 24h forecast valid at the same time. We consider three scenarios: (1) the assimilation of ZTDs (2) the assimilation of horizontal delay gradients and (3) the assimilation of both ZTDs and horizontal delay gradients. The impact is measured by utilizing the refractivity fields. We find that the assimilation of the horizontal delay gradients in addition to the ZTDs improves the refractivity field around 800 hPa. When we consider a single station there is a clear improvement when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients contain information that is not contained in the ZTDs. On the other hand, when we consider a dense station network there is not a significant improvement when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients do not contain information that is not already contained in the ZTDs. Finally, we replace simulated by real observations, that is, tropospheric parameters from a Precise Point Positioning solution provided with the G-Nut/Tefnut software, in order to show that the GFS 24h forecast is indeed improved when GNSS horizontal delay gradients are assimilated in addition to GNSS ZTDs; for the considered station (Potsdam, Germany) and period (June and July, 2017) we find an improvement in the retrieved refractivity of up to 4%. Keywords: GNSS; horizontal delay gradient; numerical weather prediction; data assimilation 1. Introduction Radio signals which are transmitted by the Global Navigation Satellite System (GNSS) constellation and received by a ground-based station allow the estimation of the Zenith Total Delay (ZTD) [ 1 ] and the horizontal delay gradient [ 2 ]. Such ZTDs and horizontal delay gradients are available from thousands of stations around the globe. Assimilation studies in meteorological models have shown that the usage of ZTDs improves the forecast skill [ 3 – 5 ]. Consequently, ZTDs are assimilated at Remote Sens. 2019,11, 41; doi:10.3390/rs11010041 www.mdpi.com/journal/remotesensing Remote Sens. 2019,11, 41 2 of 16 several weather agencies [ 6 ]. So far no assimilation studies including horizontal delay gradients have been published. To date the studies are concerned with the validation of horizontal delay gradients. In Reference [ 7 ] preliminary results were presented demonstrating that the horizontal delay gradients include real tropospheric features and are not artefacts caused by other errors such as multipath and orbital errors. However, only horizontal delay gradients from a single station were compared with horizontal delay gradients estimated from a collocated Water Vapor Radiometer. A number of studies followed. For example, Li et al. showed in Reference [ 8 ] that an improved observation geometry yields improved horizontal delay gradient estimates. Morel et al. realized in Reference [ 9 ] a comparison study on ZTDs and horizontal delay gradients from a dozen of stations located on an island in the Mediterranean sea (Corsica). Tropospheric parameters were estimated with two software packages (GAMIT and GIPSY). They found a good agreement between ZTDs (standard deviation of 6 mm) but the horizontal delay gradients were in less good agreement (standard deviation of 0.7 mm). Douša et al. presented in Reference [ 10 ] a comparison study on ZTDs and horizontal delay gradients for hundreds of stations located in central Europe. Again, tropospheric parameters were estimated with two software packages (EPOS and BERNESE). In addition, ZTDs and horizontal delay gradients were derived from two different Numerical Weather Models (NWMs); the re-analysis of the Integrated Forecast System (IFS) of the European Centre of Medium Range Weather Forecast (ECMWF) and the operational analysis of the Global Forecast System (GFS) of the National Centres for Environmental Prediction (NCEP). The comparison of GNSS and NWM horizontal delay gradients demonstrated a good agreement (standard deviation of 0.5 mm). The comparison of horizontal delay gradient maps (inspection by eye) provided a clear evidence that the GNSS horizontal delay gradients include real tropospheric features. In Reference [ 11 ] results for hundreds of stations distributed over Europe and a period of nearly two decades were presented. The comparison of GNSS and NWM ZTDs and horizontal delay gradients demonstrated a good agreement with standard deviations of 10 mm and 0.4 mm respectively. [ 12 ] provided more details for the same stations and time period. In fact, for several stations they find clear artificial signals in horizontal delay gradients which are caused by problems to receive low-elevation observations. They show that the standard deviations between GNSS and NWM ZTDs and horizontal delay gradients have a clear seasonal dependency, that is, the standard deviation is larger in summer than in winter. This can be explained by the fact that NWMs have difficulties to predict the high water vapour variability in summer. In addition, it appears that the standard deviation for the ZTD does not change over the years whereas the standard deviation for the horizontal delay gradient decreased over the years. This is an indication that the quality of GNSS horizontal delay gradients improved over the years. The GNSS ZTDs are of high quality since their introduction in meteorology in the early 90’s. In summary, recent studies begin to indicate that a single station does not only provide high-quality ZTD estimates but high-quality horizontal delay gradient estimates as well. The present study is concerned with one possible meteorological application: make use of horizontal delay gradients in numerical weather prediction. At first, we develop the operators that are required in a variational data assimilation context. Armed with these operators we estimate the impact of horizontal delay gradients in variational data assimilation. In this work we make use of NWM data from the GFS of the NCEP. The NWM data have a horizontal resolution of 0.5 ◦ on 31 pressure levels (the pressure levels range from 1000 hPa to 1 hPa). Thus, we study the impact in a meso-beta (>10 km horizontal resolution) scale NWM. The impact in a meso-gamma (<10 km horizontal resolution) scale NWM remains to be studied. The structure of this work is as follows. In Section 2we show how ZTDs and horizontal delay gradients are derived from a NWM. In Section 3we show that horizontal delay gradients can be related to ZTDs. This helps to better understand the impact when we consider observations from a sparse (dense) station network in variational data assimilation. In Section 4we show the principle behind variational data assimilation. We make use of simulated observations before we finally make use of real observations. The conclusion is given in Section 5. Remote Sens. 2019,11, 41 3 of 16 2. ZTDs and Horizontal Delay Gradients In the GNSS analysis the tropospheric delay, that is, the signal travel time delay induced by the neutral atmosphere, is parameterized utilizing Mapping Functions (MF). Specifically, for the elevation angle eand azimuth angle a, the tropospheric delay Treads as T(e,a)=mh(e)·Zh+mw(e)·Zw+mg(e)[cos(a)·N+sin(a)·E](1) where Z h denotes the Zenith Hydrostatic Delay (ZHD), Z w denotes the Zenith Wet Delay (ZWD), m h denotes the hydrostatic MF, m w denotes the wet MF, m g denotes the gradient MF, N denotes the north-gradient coefficient and E denotes the east-gradient coefficient in meter. The gradient MF reads as mg(e)=1 sin(e)tan(e)+K(2) where the constant Kequals 0.0032 [2]. The ZTD, denoted Z, is given by Z=Zh+ Zw(3) The horizontal delay gradient vector Greads as G = [E, N] (4) In the first sub-section we show how ZTDs and horizontal delay gradients are derived from a NWM. In order to make sure that what we model is what we measure, we show in the second sub-section how ZTDs and horizontal delay gradients are estimated with the Precise Point Positioning (PPP) method [13]. 2.1. ZTDs and Horizontal Delay Gradients Derived from the NWM ZTD and the horizontal delay gradients are derived from ray-traced tropospheric delays. Specifically, for a station satellite link, the tropospheric delay Sis calculated through S=Znds −g(5) where ndenotes the index of refraction, sdenotes the arclength of the ray-path and gdenotes the geometric distance. The ray-path follows from Fermat ' s principle [ 14 ]. The index of refraction nis related to the refractivity Ψthrough n=10−6Ψ+1 (6) Given the NWM pressure, temperature and humidity field the refractivity field can be computed [ 15 ]. Then the algorithm by [ 16 ] allows one to compute tropospheric delays with a high speed and precision for any station satellite link. This algorithm can be viewed as the forward operator for the tropospheric delay. The corresponding tangent-linear (adjoint) operator is constructed by application of the chain rule of differential calculus in forward (reverse) mode [17]. The ZTD is the tropospheric delay in the zenith direction. In essence, the ZTD is calculated by integrating the refractivity over height. The calculation of the horizontal delay gradients is based on the approximation for the azimuth dependency of the tropospheric delays (refer to Equation (1)) Γ(e, a) = mg(e) ·[N ·cos(a) + E ·sin(a)] (7) Remote Sens. 2019,11, 41 4 of 16 Given a set of tropospheric delays, denoted Sand a set of tropospheric delays computed under the assumption of a spherically layered atmosphere, denoted S 0 , the horizontal delay gradient vector is determined by G=(ΓTΓ)−1ΓT(S-S0) (8) where the design matrix Γ is defined by the partial derivatives of the right hand side of Equation (7) with respect to the gradient coefficients. We choose the set of tropospheric delays such that the computed horizontal delay gradient vector is independent of the actual observation geometry, that is, the set of tropospheric delays consists of tropospheric delays with elevation angles 3 ◦ , 5 ◦ , 7 ◦ , 10 ◦ , 15 ◦ , 20 ◦ , 30 ◦ , 50 ◦ , 70 ◦ , 90 ◦ and the spacing in azimuth is 30 ◦ . This set of tropospheric delays represents a sufficient discretization insofar as the refinement of azimuth (elevation) angle increments yields negligible differences in the estimated horizontal delay gradients. Another feature of this set of tropospheric delays is that the horizontal delay gradient vector actually simplifies to (refer to the Appendix A) G=(ΓTΓ)−1ΓTS (9) Hence, there is no need to calculate tropospheric delays under the assumption of a spherically layered atmosphere. 2.2. ZTDs and Horizontal Delay Gradients Estimated from the GNSS We mimic PPP. The simulation is simplified by ignoring initial carrier-phase ambiguities in the observation equation. In essence, we utilize the following simplified version of the linearized observation equation, see, for example Reference [18] S(e,a)−mh(e)·Zh=−u(e,a)·r+c·t+mw(e)·Zw+mg(e)[cos(a)·N+sin(a)·E](10) The left hand side of the equation represents the measured minus modelled term and the right hand side of the equation includes parameters to be estimated. Here u denotes the tangent-unit vector of the station-satellite link, r denotes the coordinate residual, tdenotes the clock residual and c denotes the vacuum speed of light. For a single epoch we consider 120 station satellite links where azimuth angles are selected randomly and elevation angles are obtained through e= 90 ◦ -87 ◦√ wwhere 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 tropospheric delays which enter the left hand side of the equation are ray-traced tropospheric delays. The hydrostatic and wet MF is taken from the Global Mapping Function (GMF) [ 19 ] and the a priori ZHD comes from the Global Pressure and Temperature model (GPT) [ 20 ]. For each day we consider four epochs (0, 6, 12, 18 UTC). Then we combine the four sets of equations and obtain the coordinate residuals by a least-square fit on a daily basis and the clockand tropospheric parameter residuals epoch wise. The standard elevation angle dependent weighting is applied in the least-square fit. In short, let Sdenote ray-traced tropospheric delays and Mdenote the a priori tropospheric delays, then the solution vector V, which includes the coordinate residual (daily) and the clockand tropospheric parameter residuals (epoch wise), is obtained through V=JTWJ−1JTW(S−M)(11) The design matrix J is defined by the partial derivatives of the right hand side of Equation (10) with respect to the coordinate residual, the clock residual, the ZWD and the gradient coefficients and the weight matrix Wreads as Wkl =sin(ek)sin(el)δkl (12) where the subscript indices denote the individual station satellite links. The ZTD is given by the sum of the a priori ZHD and the estimated ZWD. Remote Sens. 2019,11, 41 5 of 16 2.3. ZTD and Horizontal Delay Gradient Comparison We compare tropospheric parameters derived from the two methods described in Sections 2.1 and 2.2, in the following referred to as the NWM and the PPP method. The input data for both methods are ray-traced tropospheric delays. We select a few hundred stations, four epochs per day (0, 6, 12, 18 UTC) and a period of two months (June and July, 2017). The top panel in Figure 1 shows the station specific root-mean-square deviation between the ZTD derived from the NWM and the PPP method (the station locations are shown in the Figure 2). The deviations are about 3 mm. The deviations are due to the fact that the functional form of the tropospheric delay in PPP is inaccurate. In addition, the hydrostatic MF, the wet MF and the a priori ZHD are inaccurate, because they are based on a climatology. We also note that in PPP the tropospheric parameters are estimated together with the station coordinate and the station clock offsets. In particular, the ZTD is correlated with the station up-component and the station clock offset. The middle and lower panel in Figure 1 shows the station specific root-mean-square deviation for the eastand north-gradient coefficient, respectively. The deviations are below 0.1 mm. The small deviations can be explained by the fact that both methods determine the horizontal delay gradients in a very similar way (least-square fit). The remaining deviations are due to the imperfect observation geometry in PPP. With regards to the gradient coefficients the hydrostatic and wet MF does not matter. In essence, we replaced the hydrostatic and wet MF, which is taken from the GMF, by the hydrostatic and wet MF derived from the NWM and found negligible differences in the estimated gradient coefficients. However, we must agree on the gradient MF [ 21 ]. We used for both methods the same gradient MF [ 2 ]. In conclusion, the comparison shows that what we model is what we can measure. The forward operators for the ZTD and the horizontal delay gradients work as intended. Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 16 station specific root-mean-square deviation for the eastand north-gradient coefficient, respectively. The deviations are below 0.1 mm. The small deviations can be explained by the fact that both methods determine the horizontal delay gradients in a very similar way (least-square fit). The remaining deviations are due to the imperfect observation geometry in PPP. With regards to the gradient coefficients the hydrostatic and wet MF does not matter. In essence, we replaced the hydrostatic and wet MF, which is taken from the GMF, by the hydrostatic and wet MF derived from the NWM and found negligible differences in the estimated gradient coefficients. However, we must agree on the gradient MF [21]. We used for both methods the same gradient MF [2]. In conclusion, the comparison shows that what we model is what we can measure. The forward operators for the ZTD and the horizontal delay gradients work as intended. Figure 1. Tropospheric parameter differences between the NWM and PPP method. The top panel shows the station specific root-mean-square deviation for the ZTD. The middle and lower panel shows the station specific root-mean-square deviation for the northand east-gradient coefficients, respectively. We consider four epochs per day (0, 6, 12, 18 UTC) and a period of two months (June and July, 2017). 3. Relation between Horizontal Delay Gradients and ZTDs Let Ψh and Ψw denote the hydrostatic and wet refractivity, respectively. The hydrostatic delay Sh is defined as 𝑆=10𝛹𝑑𝑠+𝑠−𝑔 (13) and the wet delay Sw is defined as 𝑆=10𝛹𝑑𝑠 (14) such that 𝑆=𝑆+𝑆 (15) Therefore, the horizontal delay gradient vector can be written as 𝐺=𝐺+𝐺 (16) where the hydrostatic horizontal delay gradient vector Gh is given by 𝐺=(𝜞𝑻𝜞)𝜞𝑻𝑆 (17) and the wet horizontal delay gradient vector Gw is given by 𝐺=(𝜞𝑻𝜞)𝜞𝑻𝑆 (18) Figure 1. Tropospheric parameter differences between the NWM and PPP method. The top panel shows the station specific root-mean-square deviation for the ZTD. The middle and lower panel shows the station specific root-mean-square deviation for the northand east-gradient coefficients, respectively. We consider four epochs per day (0, 6, 12, 18 UTC) and a period of two months (June and July, 2017). 3. Relation between Horizontal Delay Gradients and ZTDs Let Ψh and Ψw denote the hydrostatic and wet refractivity, respectively. The hydrostatic delay S h is defined as Sh=10−6ZΨhds +s−g(13) and the wet delay Swis defined as Sw=10−6ZΨwds (14) Remote Sens. 2019,11, 41 6 of 16 such that S=Sh+Sw(15) Therefore, the horizontal delay gradient vector can be written as G=Gh+Gw(16) where the hydrostatic horizontal delay gradient vector Ghis given by Gh=ΓTΓ−1ΓTSh(17) and the wet horizontal delay gradient vector Gwis given by Gw=ΓTΓ−1ΓTSw(18) Let us assume that for some station we estimate the ZTD and the horizontal delay gradient with the GNSS and let us assume that the ZHD and the hydrostatic horizontal delay gradient are provided with high accuracy by a NWM. The following relation is exactly true Zw=Z−Zh(19) and the following relation is approximately true ∂Zw ∂x,∂Zw ∂y∼(G−Gh) F(20) where the factor Fcan be related to the scale height of the wet refractivity [ 22 ]. Roughly speaking, if we know the wet horizontal delay gradient at the station then we know the ZWD gradient at the station and vice versa. In order to demonstrate this, we chose some station and select stations in the vicinity of that station (stations within a radius of 1 ◦ ). Next, we collect the ZWDs from the stations, assume that the ZWD can be expanded in a Taylor-series and obtain by a least-square-fit the ZWD gradient at the chosen station. This ZWD gradient at the station is roughly proportional to the wet horizontal delay gradient at the station. We chose some epoch where large horizontal delay gradients are present (2 July 2017, 12 UTC) and chose the factor Fto be 3 km. The left panels of Figure 2show the wet east and north-gradient coefficient and the right panels of Figure 2show the approximation for the wet east and north-gradient coefficient. The gradient coefficients do not match perfectly but the validity of the approximation is obvious. The implication of this for variational data assimilation is the following. Let us consider a single station or a sparse station network. We expect to find a significant impact when the horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients contain information that is obviously not contained in the ZTDs. Next, let us consider a dense station network. In this case we do not expect to find a significant impact when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients do not contain information that is not already contained in the ZTDs. In essence, when we assimilate the ZTDs, we adjust the ZTDs, we adjust the ZTD gradients and hence we automatically adjust the horizontal delay gradients. However, this argument overlooks the scale height of the wet refractivity. In fact, if we know both the ZTDs and the horizontal delay gradients the scale height of the wet refractivity can be estimated [ 20 ]. This means when we assimilate both, ZTDs and the horizontal delay gradients, the scale height of the wet refractivity is adjusted as well. Remote Sens. 2019,11, 41 7 of 16 Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 16 Let us assume that for some station we estimate the ZTD and the horizontal delay gradient with the GNSS and let us assume that the ZHD and the hydrostatic horizontal delay gradient are provided with high accuracy by a NWM. The following relation is exactly true 𝑍=𝑍−𝑍 (19) and the following relation is approximately true 󰇣  , 󰇤∼()  (20) where the factor F can be related to the scale height of the wet refractivity [22]. Roughly speaking, if we know the wet horizontal delay gradient at the station then we know the ZWD gradient at the station and vice versa. In order to demonstrate this, we chose some station and select stations in the vicinity of that station (stations within a radius of 1°). Next, we collect the ZWDs from the stations, assume that the ZWD can be expanded in a Taylor-series and obtain by a least-square-fit the ZWD gradient at the chosen station. This ZWD gradient at the station is roughly proportional to the wet horizontal delay gradient at the station. We chose some epoch where large horizontal delay gradients are present (2 July 2017, 12 UTC) and chose the factor F to be 3 km. The left panels of Figure 2 show the wet east and north-gradient coefficient and the right panels of Figure 2 show the approximation for the wet east and north-gradient coefficient. The gradient coefficients do not match perfectly but the validity of the approximation is obvious. (a) (b) Figure 2. The left panels (a) show the wet east (E) and north-gradient (N) coefficient and the right panels (b) show the approximation for the wet east (E*) and north-gradient (N*) coefficient (2 July 2017, 12 UTC). The approximations for the wet gradient coefficients are obtained by utilizing ZWD gradients. For details refer to the text. The implication of this for variational data assimilation is the following. Let us consider a single station or a sparse station network. We expect to find a significant impact when the horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients contain information that is obviously not contained in the ZTDs. Next, let us consider a dense station network. In this case we do not expect to find a significant impact when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients do not contain information that is not already contained in the ZTDs. In essence, when we assimilate the ZTDs, we Figure 2. The left panels ( a ) show the wet east (E) and north-gradient (N) coefficient and the right panels ( b ) show the approximation for the wet east (E*) and north-gradient (N*) coefficient (2 July 2017, 12 UTC). The approximations for the wet gradient coefficients are obtained by utilizing ZWD gradients. For details refer to the text. 4. Variational Data Assimilation We consider the following inverse problem: given the observations yand the background refractivity field bdetermine the most probable refractivity field a. By taking into account Bayes’ theorem under the assumption that the background and observation errors have Gaussian statistics and are independent of each other, the most probable refractivity field is the one that minimizes the cost function [23] C[x]=1 2(x−b)TB−1(x−b)+1 2(y−H[x])TR−1(y−H[x]) (21) Here B and R denote the background and observation error covariance matrix, respectively. The diagonal elements of the matrix B store the squared refractivity error variances and the off-diagonal elements store the refractivity error covariances. Likewise, the diagonal elements of the matrix R store the squared observation error variances and the off-diagonal elements store the observation error covariances. The minimum of the cost function Cmust be determined iteratively since the forward operator His a nonlinear function of the refractivity field x. The first guess refractivity field is chosen to be the background refractivity field. In particular, if we perform a single iteration the most probable refractivity field acan be written as a=b+BHT[b]H[b]BHT[b]+R−1(y−H[b]) (22) Here H and HT denote the tangent-linear and adjoint operator, respectively. Hereinafter, we will call the most probable refractivity field simply the analysis refractivity field. Further, when we apply the ZTD forward operator to the background (analysis) refractivity field, we will call the resulting ZTD the background (analysis) ZTD. Likewise, when we apply the horizontal delay gradient forward Remote Sens. 2019,11, 41 8 of 16 operator to the background (analysis) refractivity field, we will call the resulting horizontal delay gradient the background (analysis) horizontal delay gradient. We investigate three scenarios: (1) the assimilation of ZTDs, (2) the assimilation of horizontal delay gradients and (3) the assimilation of both ZTDs and horizontal delay gradients. 4.1. Experiment with Simulated Observations In brief, the simulation works as follows: we chose some true refractivity field, extract error-free observations, add noise, assimilate the observations into some background refractivity field and finally check if the analysis refractivity field is closer to the true refractivity field than the background refractivity field. The details are as follows. The true refractivity field is derived from the GFS analysis. The background refractivity field is derived from the GFS 24h forecast valid at the same time. The background error covariance matrix is estimated by forming the difference between the GFS 48h and 24h forecast valid at the same time. This method is widely used in variational data assimilation [ 24 ] and also known as the National Meteorological Centre (NMC) method. We assume that refractivity at each grid point is correlated with the refractivity of its neighbouring grid points. Specifically, the vertical and horizontal error correlation length is chosen to be 500 m and 0.5 ◦ respectively. This defines the background error covariance matrix. The error-free observations are derived from the GFS analysis utilizing the developed forward operators. Gaussian errors with zero mean and a specific standard deviation are added to obtain observations. The standard deviation is chosen to be 10 mm for the ZTD and 0.5 mm for the gradient coefficients. This defines the observation error covariance matrix. The chosen standard deviations are rational insofar as they correspond to the typical standard deviations between GNSS and NWM tropospheric parameters [ 10 ]. We also note that with this choice the observed ZTD is approximately as accurate as the background ZTD and the observed horizontal delay gradient is approximately as accurate as the background horizontal delay gradient (see below), a prerequisite in variational data assimilation. Data are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). In the first sub-section we consider observations from a single station. The station is located in Potsdam (Germany). In the second sub-section we consider observations from a dense station network. The station network consists of 25 stations with a station separation of 0.5 ◦ where the centre station is the station Potsdam. The impact in observation space is measured by analysing the ZTDs and horizontal delay gradients at the station Potsdam. The impact in model space is measured by analysing the four refractivity profiles surrounding the station Potsdam. 4.1.1. Single Station Figure 3summarizes the statistics in observation space when a single station is used. If ZTDs are assimilated (left panel) the analysis ZTDs are more accurate than the background ZTDs. The analysis horizontal delay gradients are not more accurate than the background horizontal delay gradients. If horizontal delay gradients are assimilated (middle panel) the analysis horizontal delay gradients are more accurate than the background horizontal delay gradients. The analysis ZTDs are not more accurate than the background ZTDs. In short, if we assimilate ZTDs only, the ZTDs but not the horizontal delay gradients are adjusted and if we assimilate horizontal delay gradients only, the horizontal delay gradients but not the ZTDs are adjusted. Finally, if ZTDs and horizontal delay gradients are assimilated (right panel) both, ZTDs and horizontal delay gradients, are adjusted. Figure 4shows the mean and standard deviation between the analysis and background refractivity as a function of the pressure for the four refractivity profiles surrounding the station. If ZTDs are assimilated (left panel) the four refractivity profiles surrounding the station are not equally affected. The refractivity profile which is nearest to the station is most affected. If we assimilate horizontal delay gradients (middle panel) the four refractivity profiles surrounding the station are equally affected. The assimilation of the horizontal delay gradients yields no refractivity deviation at low Remote Sens. 2019,11, 41 9 of 16 and high altitudes. The largest refractivity deviation is around 800 hPa. If ZTDs and horizontal delay gradients are assimilated (right panel) the four refractivity profiles surrounding the station are affected accordingly. In essence, we find the combined impact of ZTDs and horizontal delay gradients. Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 16 observed horizontal delay gradient is approximately as accurate as the background horizontal delay gradient (see below), a prerequisite in variational data assimilation. Data are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). In the first sub-section we consider observations from a single station. The station is located in Potsdam (Germany). In the second sub-section we consider observations from a dense station network. The station network consists of 25 stations with a station separation of 0.5° where the centre station is the station Potsdam. The impact in observation space is measured by analysing the ZTDs and horizontal delay gradients at the station Potsdam. The impact in model space is measured by analysing the four refractivity profiles surrounding the station Potsdam. 4.1.1. Single Station Figure 3 summarizes the statistics in observation space when a single station is used. If ZTDs are assimilated (left panel) the analysis ZTDs are more accurate than the background ZTDs. The analysis horizontal delay gradients are not more accurate than the background horizontal delay gradients. If horizontal delay gradients are assimilated (middle panel) the analysis horizontal delay gradients are more accurate than the background horizontal delay gradients. The analysis ZTDs are not more accurate than the background ZTDs. In short, if we assimilate ZTDs only, the ZTDs but not the horizontal delay gradients are adjusted and if we assimilate horizontal delay gradients only, the horizontal delay gradients but not the ZTDs are adjusted. Finally, if ZTDs and horizontal delay gradients are assimilated (right panel) both, ZTDs and horizontal delay gradients, are adjusted. Figure 4 shows the mean and standard deviation between the analysis and background refractivity as a function of the pressure for the four refractivity profiles surrounding the station. If ZTDs are assimilated (left panel) the four refractivity profiles surrounding the station are not equally affected. The refractivity profile which is nearest to the station is most affected. If we assimilate horizontal delay gradients (middle panel) the four refractivity profiles surrounding the station are equally affected. The assimilation of the horizontal delay gradients yields no refractivity deviation at low and high altitudes. The largest refractivity deviation is around 800 hPa. If ZTDs and horizontal delay gradients are assimilated (right panel) the four refractivity profiles surrounding the station are affected accordingly. In essence, we find the combined impact of ZTDs and horizontal delay gradients. Figure 5 shows the root mean square error of the background and the analysis refractivity as a function of the pressure for the four refractivity profiles surrounding the station. If we assimilate ZTDs (left panel) we find a positive impact for the refractivity profile which is nearest to the station. If we assimilate horizontal delay gradients (middle panel) we find a positive impact for all four refractivity profiles surrounding the station. The assimilation of the horizontal delay gradients in addition to the ZTDs (right panel) improves the refractivity field around 800 hPa. This is in line with our expectation: there is an improvement when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients contain information that is not contained in the ZTDs. Figure 3. The root mean square errors of ZTDs and the horizontal delay gradients for the Background (B), the Observation (O) and the Analysis (A). Left panel: ZTDs are assimilated. Middle panel: Horizontal delay gradients are assimilated. Right panel: ZTDs and horizontal delay gradients Figure 3. The root mean square errors of ZTDs and the horizontal delay gradients for the Background (B), the Observation (O) and the Analysis (A). Left panel: ZTDs are assimilated. Middle panel: Horizontal delay gradients are assimilated. Right panel: ZTDs and horizontal delay gradients are assimilated. Data from the single station Potsdam are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 16 are assimilated. Data from the single station Potsdam are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). Figure 4. The mean (red) and one-sigma (black) refractivity deviation in percent between the background and the analysis as a function of the pressure. Left panel: ZTDs are assimilated. Middle panel: Horizontal delay gradients are assimilated. Right panel: ZTDs and horizontal delay gradients are assimilated. In each panel the four plots correspond to the grid points surrounding the station Potsdam. The labels (a), (b), (c) and (d) correspond to the grid points to the north-west, north-east, south-west and south-east respectively. Data from the single station Potsdam are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). Figure 5. The root mean square error of the background (analysis) refractivity in percent as a function of the pressure in black (red). Left panel: ZTDs are assimilated. Middle panel: horizontal delay gradients are assimilated. Right panel: ZTDs and Horizontal delay gradients are assimilated. In each panel the four plots correspond to the grid points surrounding the station Potsdam. The labels (a), (b), (c) and (d) correspond to the grid points to the north-west, north-east, south-west and south-east respectively. Data from the single station Potsdam are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). 4.1.2. Station Network Figure 6 summarizes the statistics in observation space when a network of stations is used. If we assimilate ZTDs (left panel) the analysis ZTDs are more accurate than the background ZTDs and the analysis horizontal delay gradients are more accurate than the background horizontal delay gradients as well. If we assimilate horizontal delay gradients (middle panel) the analysis horizontal delay gradients are more accurate than the background horizontal delay gradients and the analysis ZTDs are more accurate than the background ZTDs as well. In short, if we assimilate ZTDs only, both the ZTDs and the horizontal delay gradients are adjusted. Similarly, if we assimilate horizontal delay gradients only, both the horizontal delay gradients and the ZTDs are adjusted. Finally, if ZTDs and horizontal delay gradients are assimilated (right panel) both, ZTDs and horizontal delay gradients, are adjusted accordingly. In essence, when we assimilate the ZTDs, we do not only adjust the ZTDs but we automatically adjust the horizontal delay gradients. Figure 7 shows the mean and standard deviation between the analysis and background refractivity as a function of the pressure for the four refractivity profiles surrounding the station. If we assimilate ZTDs (left panel), the refractivity profiles are equally affected. Likewise, if we assimilate horizontal delay gradients (middle panel), the refractivity profiles are equally affected. Figure 4. The mean (red) and one-sigma (black) refractivity deviation in percent between the background and the analysis as a function of the pressure. Left panel: ZTDs are assimilated. Middle panel: Horizontal delay gradients are assimilated. Right panel: ZTDs and horizontal delay gradients are assimilated. In each panel the four plots correspond to the grid points surrounding the station Potsdam. The labels (a), (b), (c) and (d) correspond to the grid points to the north-west, north-east, south-west and south-east respectively. Data from the single station Potsdam are assimilated four times per day (0, 6, 12, 18 UTC) for a period of two months (June and July, 2017). Figure 5shows the root mean square error of the background and the analysis refractivity as a function of the pressure for the four refractivity profiles surrounding the station. If we assimilate ZTDs (left panel) we find a positive impact for the refractivity profile which is nearest to the station. If we assimilate horizontal delay gradients (middle panel) we find a positive impact for all four refractivity profiles surrounding the station. The assimilation of the horizontal delay gradients in addition to the ZTDs (right panel) improves the refractivity field around 800 hPa. This is in line with our expectation: there is an improvement when horizontal delay gradients are assimilated in addition to the ZTDs because the horizontal delay gradients contain information that is not contained in the ZTDs. Remote Sens. 2019,11, 41 16 of 16 8. Li, X.; Zus, F.; Lu, C.; Ning, T.; Dick, G.; Ge, M.; Wickert, J.; Schuh, H. Retrieving high resolution tropospheric gradients from multiconstellation GNSS observations. Geophys. Res. Lett. 2015,42, 4173–4181. [CrossRef] 9. Morel, L.; Pottiaux, E.; Durand, F.; Fund, F.; Boniface, K.; de Oliveira, P.S.; Van Baelen, J. Validity and behaviour of tropospheric gradients estimated by GPS in Corsica. Adv. Space Res. 2015 ,55, 135–149. [CrossRef] 10. Douša, J.; Dick, G.; Kaˇcmaˇrík, M.; Brožková, R.; Zus, F.; Brenot, H.; Stoycheva, A.; Möller, G.; Kaplon, J. Benchmark campaign and case study episode in central Europe for development and assessment of advanced GNSS tropospheric models and products. Atmos. Meas. Tech. 2016,9, 2989–3008. [CrossRef] 11. Pacione, R.; Araszkiewicz, A.; Brockmann, E.; Dousa, J. EPN-Repro2: A reference GNSS tropospheric data set over Europe. Atmos. Meas. Tech. 2017,10, 1689–1705. [CrossRef] 12. Douša, J.; Václavovic, P.; Eliaš, M. Tropospheric products of the second European GNSS reprocessing (1996–2014). Atmos. Meas. Tech. 2017,10, 3589–3607. [CrossRef] 13. Zumberge, J.F.; Heflin, M.B.; Jefferson, D.C.; Watkins, M.M.; Webb, F.H. Precise point positioning for the efficient and robust analysis of GPS data from large networks. J. Geophys. Res. 1997 ,102, 5005–5017. [CrossRef] 14. Zus, F.; Bender, M.; Deng, Z.; Dick, G.; Heise, S.; Shang-Guan, M.; Wickert, J. A methodology to compute GPS slant total delays in a numerical weather model. Radio Sci. 2012,47, RS2018. [CrossRef] 15. Thayer, G.D. An improved equation for the radio refractive index of air. Radio Sci. 1974 ,9, 803–807. [CrossRef] 16. Zus, F.; Dick, G.; Douša, J.; Heise, S.; Wickert, J. The rapid and precise computation of GPS slant total delays and mapping factors utilizing a numerical weather model. Radio Sci. 2014,49, 207–216. [CrossRef] 17. Zus, F.; Dick, G.; Heise, S.; Wickert, J. A forward operator and its adjoint for GPS slant total delays. Radio Sci. 2015,50, 393–405. [CrossRef] 18. Zus, F.; Deng, Z.; Wickert, J. The impact of higher-order ionospheric effects on estimated tropospheric parameters in Precise Point Positioning. Radio Sci. 2017,52, 963–971. [CrossRef] 19. Boehm, J.; Niell, A.; Tregoning, P.; Schuh, H. Global mapping function (GMF): A new empirical mapping function based on numerical weather model data. Geophys. Res. Lett. 2006,33, 943–951. [CrossRef] 20. Boehm, J.; Heinkelmann, R.; Schuh, H. Short Note: A global model of pressure and temperature for geodetic applications. J. Geod. 2007,81, 679–683. [CrossRef] 21. Kaˇcmaˇrík, M.; Douša, J.; Zus, F.; Václavovic, P.; Balidakis, K.; Dick, G.; Wickert, J. Sensitivity of GNSS tropospheric gradients to processing options. Ann. Geophys. Discuss. 2018. in review. [CrossRef] 22. Elosegui, P.; Davis, J.L.; Gradinarsky, L.P.; Elgered, G.; Johansson, J.M.; Tahmoush, D.A.; Rius, A. Sensing atmospheric structure using small-scale space geodetic networks. Geophys. Res. Lett. 1999 ,26, 2445–2448. [CrossRef] 23. Daley, R. Atmospheric Data Analysis; Cambridge Atmospheric and Space Science Series; Cambridge University Press: Cambridge, UK, 1991. 24. Parrish, D.; Derber, J. The National Meteorological Center’s spectral statistical analysis system. Mon. Weather Rev. 1992,120, 1747–1763. [CrossRef] 25. Douša, J.; Václavovic, P.; Zhao, L.; Kaˇcmaˇrík, M. New Adaptable All-in-One Strategy for Estimating Advanced Tropospheric Parameters and Using Real-Time Orbits and Clocks. Remote Sens. 2018,10, 232. [CrossRef] 26. Václavovic, P.; Douša, J. Backward smoothing for precise GNSS applications. Adv. Space Res. 2015 ,56, 1627–1634. [CrossRef] 27. Gérard, P.; Luzum, B. (Eds.) IERS Technical Note No. 36; Verlag des Bundesamts für Kartographie und Geodäsie: Frankfurt am Main, Germany, 2010; 179p, ISBN 3-89888-989-6. 28. Schmid, R.; Dach, R.; Collilieux, X.; Jaeggi, A.; Schmitz, M.; Dilssner, F. Absolute IGS antenna phase center model igs08.atx: Status and potential improvements. J. Geod. 2016,90, 343–364. [CrossRef] 29. Rebischung, P.; Altamimi, Z.; Ray, J.; Garayt, B. The IGS contribution to ITRF2014. J. Geod. 2016 ,90, 611–630. [CrossRef] 30. Dow, J.M.; Neilan, R.E.; Rizos, C. The International GNSS Service in a changing landscape of Global Navigation Satellite Systems. J. Geod. 2009,83, 191–198. [CrossRef] © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).