Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation

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


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][4][5]. Consequently, ZTDs are assimilated at 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 2 we show how ZTDs and horizontal delay gradients are derived from a NWM. In Section 3 we 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 4 we 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.

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 e and azimuth angle a, the tropospheric delay T reads as 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 m g (e) = 1 sin(e) tan(e) + K where the constant K equals 0.0032 [2]. The ZTD, denoted Z, is given by The horizontal delay gradient vector G reads as 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].

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 S is calculated 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 from Fermat's principle [14]. The index of refraction n is related to the refractivity Ψ through n = 10 −6 Ψ + 1 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) = m g (e) · [N · cos(a) + E · sin(a)] Given a set of tropospheric delays, denoted S and 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 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 Hence, there is no need to calculate tropospheric delays under the assumption of a spherically layered atmosphere.

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 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, t denotes 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 • √ 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 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 clock-and tropospheric parameter residuals epoch wise. The standard elevation angle dependent weighting is applied in the least-square fit. In short, let S denote ray-traced tropospheric delays and M denote the a priori tropospheric delays, then the solution vector V, which includes the coordinate residual (daily) and the clock-and tropospheric parameter residuals (epoch wise), is obtained through 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 W reads as W kl = sin(e k )sin(e l )δ 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.

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 east-and 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 east-and 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. The middle and lower panel shows the station specific root-mean-square deviation for the north-and 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).

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 and the wet delay Sw is defined as such that = + (15) 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 north-and 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).

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 and the wet delay S w is defined as Remote Sens. 2019, 11, 41 Therefore, the horizontal delay gradient vector can be written as where the hydrostatic horizontal delay gradient vector G h is given by and the wet horizontal delay gradient vector G w is given by 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 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. 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.
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. 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

Variational Data Assimilation
We consider the following inverse problem: given the observations y and the background refractivity field b determine 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 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 C must be determined iteratively since the forward operator H is 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 a can be written as Here H and H T 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 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.

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

Single Station
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 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.
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 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 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. However, there is no deviation close to Earth's surface. If ZTDs and horizontal delay gradients are assimilated (right panel) the refractivity profiles are affected accordingly. However, there is no deviation close to Earth's surface. If ZTDs and horizontal delay gradients are assimilated (right panel) the refractivity profiles are affected accordingly. Figure 8 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. Likewise, if we assimilate horizontal delay gradients (middle panel) we find a positive impact. The assimilation of horizontal delay gradients in addition to ZTDs (right panel) yields a small improvement in the refractivity field around 800 hPa. Again, this is in line with our expectation: there is no significant improvement when horizontal delay gradients are assimilated in addition to ZTDs because the horizontal delay gradients do not contain information that is not already contained in ZTDs from dense network of stations.    line with our expectation: there is no significant improvement when horizontal delay gradients are assimilated in addition to ZTDs because the horizontal delay gradients do not contain information that is not already contained in ZTDs from dense network of stations.

GNSS Analysis
The G-Nut/Tefnut software [25] was utilized for GNSS data processing in this study. The PPP method was used when supported with the precise orbits and clocks from the European Space Agency, ESA (http://navigation-office.esa.int/GNSS_based_products.html). The extended Kalman filter and backward smoothing methods [26] were applied for estimating station coordinates (in static mode), clock corrections, initial phase ambiguities and site-/epoch-wise tropospheric parameters. When using GPS and GLONASS observations together, inter-system bias was

GNSS Analysis
The G-Nut/Tefnut software [25] was utilized for GNSS data processing in this study. The PPP method was used when supported with the precise orbits and clocks from the European Space Agency, ESA (http://navigation-office.esa.int/GNSS_based_products.html). The extended Kalman filter and backward smoothing methods [26] were applied for estimating station coordinates (in static mode), clock corrections, initial phase ambiguities and site-/epoch-wise tropospheric parameters. When using GPS and GLONASS observations together, inter-system bias was

GNSS Analysis
The G-Nut/Tefnut software [25] was utilized for GNSS data processing in this study. The PPP method was used when supported with the precise orbits and clocks from the European Space Agency, ESA (http://navigation-office.esa.int/GNSS_based_products.html). The extended Kalman filter and backward smoothing methods [26] were applied for estimating station coordinates (in static mode), clock corrections, initial phase ambiguities and site-/epoch-wise tropospheric parameters. When using GPS and GLONASS observations together, inter-system bias was additionally estimated, while GLONASS inter-frequency biases were practically absorbed in initial phase ambiguities.
Both GPS and GLONASS pseudo range (carrier-phase) observations were utilized with a priori sigma of 2 m (0.02 m) and 3 m (0.03 m), respectively. The elevation angle cut-off for observations below 7 • was applied together with the standard elevation angle dependent data weighting. Applied models were compliant with the IERS 2010 conventions [27]. We also used the igs14_2013.ATX file for absolute antenna phase centre offsets and variations [28], that is consistent with the IGS2014 reference frame [29] used for precise satellite positions; both products provided by the International GNSS Service, IGS [30].
We estimated tropospheric parameters at a 5-min sampling rate as a stochastic process with a random walk of 6 mm/ √ h and 0.6 mm/ √ h for ZTDs and horizontal delay gradients, respectively. A priori values of ZTDs were derived from the GPT [20] and introduced into the adjustment model with a priori sigma of 10 cm. Null values were introduced as a priori values for horizontal delay gradients with an a priori sigma of 5 mm. The hydrostatic and wet MFs employed were those of the GMF [19] and the gradient MF employed was that of Chen and Herring [2].

Experiment Design
We consider only data from the station Potsdam. The simulated observations are replaced by real observations. The assimilation setup is the same as in the previous Section 4.1.1 with only one difference: we perform a bias correction for the ZTD. The bias in the ZTD is the mean difference between the GNSS and NWM ZTD for the considered time period. We suspect that this bias in the ZTD (5 mm) does not come from the GNSS but from the NWM [10]. No bias corrections are performed for the gradient coefficients as the mean differences between the GNSS and NWM gradient coefficients are negligible (< 0.1 mm). Figure 9 summarizes the statistics in observation space when real data from a single station is used. At first we note that the observed ZTDs are as accurate as the background ZTDs (the root mean square error is about 10 mm) and the observed horizontal delay gradients are as accurate as the background horizontal delay gradients (the root mean square error is about 0.5 mm). Therefore, it is not surprising that the results we obtain with real observations are comparable with the results we obtain for the simulated observations. In short, if we assimilate ZTDs only (left panel), the ZTDs but not the horizontal delay gradients are adjusted and if we assimilate horizontal delay gradients only (middle panel), the horizontal delay gradients but not the ZTDs are adjusted. If ZTDs and horizontal delay gradients are assimilated (right panel) both, ZTDs and horizontal delay gradients, are adjusted. For example, if we assimilate ZTDs the refractivity is improved by up to 6%. If we assimilate ZTDs and horizontal delay gradients the refractivity is further improved by up to 4%. We find that these values for the percentage improvement are completely consistent with the values for the percentage improvement we obtain when we use simulated observations    Figure 10 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. If we assimilate horizontal delay gradients (middle panel) the four refractivity profiles surrounding the station are equally affected. As expected, 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. Figure 11 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. The assimilation of ZTDs and horizontal delay gradients improves the GFS 24h forecast refractivity field (verified with its own analysis). The direct comparison of the left and right panels reveals that the assimilation of the horizontal delay gradients in addition to the ZTD has a positive impact.        In order to provide a value for the percentage improvement we average over the refractivity profiles surrounding the station. Figure 12 shows the root mean square error of the analysis refractivity as a function of the pressure when we assimilate ZTDs only (blue line) and when we assimilate both ZTDs and horizontal delay gradients (red line). For comparison the root mean square error of the background refractivity as a function of the pressure is shown as well (black line). For example, if we assimilate ZTDs the refractivity is improved by up to 6%. If we assimilate ZTDs and horizontal delay gradients the refractivity is further improved by up to 4%. We find that these values for the percentage improvement are completely consistent with the values for the percentage improvement we obtain when we use simulated observations. Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 16 for 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).

Conclusions
We developed the forward operator for horizontal delay gradients and its adjoint. The operator for the horizontal delay gradient is costly when compared to the operator for the ZTD. In fact, the calculation of a single horizontal delay gradient vector requires the calculation of 120 tropospheric delays. Reducing the computational cost must be subject to future work.
Armed with the developed operator we estimate the impact of horizontal delay gradients in variational data assimilation. 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. Our results are representative for meso-beta scale NWMs. The impact in a meso-gamma scale NWM remains to be studied.
We have shown that the assimilation of GNSS horizontal delay gradients in addition to the GNSS ZTDs improves the GFS 24h forecast refractivity field; for the considered station (Potsdam, Germany) and period (June and July, 2017) we find an improvement in the retrieved refractivity of up to 4%. We did not make any attempt to improve the GFS analysis refractivity field itself. Experiments with simulated observations followed by experiments with real observations should be carried out utilizing a sophisticated data assimilation system.

Conclusions
We developed the forward operator for horizontal delay gradients and its adjoint. The operator for the horizontal delay gradient is costly when compared to the operator for the ZTD. In fact, the calculation of a single horizontal delay gradient vector requires the calculation of 120 tropospheric delays. Reducing the computational cost must be subject to future work.
Armed with the developed operator we estimate the impact of horizontal delay gradients in variational data assimilation. 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. Our results are representative for meso-beta scale NWMs. The impact in a meso-gamma scale NWM remains to be studied.
We have shown that the assimilation of GNSS horizontal delay gradients in addition to the GNSS ZTDs improves the GFS 24h forecast refractivity field; for the considered station (Potsdam, Germany) and period (June and July, 2017) we find an improvement in the retrieved refractivity of up to 4%. We did not make any attempt to improve the GFS analysis refractivity field itself. Experiments with simulated observations followed by experiments with real observations should be carried out utilizing a sophisticated data assimilation system.

Appendix A
The horizontal delay gradient vector can be written as We now show that for the chosen observation geometry the second term in the equation above vanishes. To see this we note that Γ T S 0 = m g (e 1 ) cos(a 1 ) . . . m g (e n ) cos(a n ) m g (e 1 ) sin(a 1 ) . . . m g (e n ) sin(a n )