Open AccessArticle

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

Florian Zus

^1,*,

Jan Douša

²,

Michal Kačmařík

³,

Pavel Václavovic

²,

Galina Dick

¹ and

Jens Wickert

^1,4

GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany

Geodetic Observatory Pecný, Research Institute of Geodesy, Topography and Cartography, Zdiby 250 66, Czech Republic

Department of Geoinformatics, VŠB – Technical University of Ostrava, Ostrava 708 33, Czech Republic

⁴

Institute of Geodesy and Geoinformation Science, Technical University of Berlin, 10623 Berlin, Germany

Author to whom correspondence should be addressed.

Remote Sens. 2019, 11(1), 41; https://doi.org/10.3390/rs11010041

Submission received: 26 November 2018 / Revised: 20 December 2018 / Accepted: 20 December 2018 / Published: 28 December 2018

(This article belongs to the Special Issue Selected Papers from IGL-1 2018 — First International Workshop on Innovating GNSS and LEO Occultations & Reflections for Weather, Climate and Space Weather)

Download

Browse Figures

Versions Notes

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

Graphical Abstract

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

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 e and azimuth angle a, the tropospheric delay T reads as

T (e, a) = m_{h} (e) \cdot Z_{h} + m_{w} (e) \cdot Z_{w} + m_{g} (e) [\cos (a) \cdot N + \sin (a) \cdot 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

m_{g} (e) = \frac{1}{\sin (e) \tan (e) + K}

(2)

where the constant K equals 0.0032 [2]. The ZTD, denoted Z, is given by

Z = Z_h + Z_w

(3)

The horizontal delay gradient vector G reads 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 S is calculated through

S = \int n d s - g

(5)

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

(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) = m_g (e) ⋅ [N ⋅ cos(a) + E ⋅ sin(a)]

(7)

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₀, the horizontal delay gradient vector is determined by

G = (Γ^TΓ)⁻¹ Γ^T (S-S₀)

(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)

G = (Γ^TΓ)⁻¹ Γ^T S

(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) - m_{h} (e) \cdot Z_{h} = - u (e, a) \cdot r + c \cdot t + m_{w} (e) \cdot Z_{w} + m_{g} (e) [c o s (a) \cdot N + s i n (a) \cdot 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, 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

V = {(J^{T} W J)}^{- 1} J^{T} W (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 W reads as

W_{k l} = s i n (e_{k}) s i n (e_{l}) δ_{k l}

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

2.3. ZTD and Horizontal Delay Gradient Comparison

We compare tropospheric parameters derived from the two methods described in Section 2.1 and Section 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.

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

S_{h} = 10^{- 6} \int Ψ_{h} d s + s - g

(13)

and the wet delay S_w is defined as

S_{w} = 10^{- 6} \int Ψ_{w} d s

(14)

such that

S = S_{h} + S_{w}

(15)

Therefore, the horizontal delay gradient vector can be written as

G = G_{h} + G_{w}

(16)

where the hydrostatic horizontal delay gradient vector G_h is given by

G_{h} = {(Γ^{T} Γ)}^{- 1} Γ^{T} S_{h}

(17)

and the wet horizontal delay gradient vector G_w is given by

G_{w} = {(Γ^{T} Γ)}^{- 1} Γ^{T} S_{w}

(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

Z_{w} = Z - Z_{h}

(19)

and the following relation is approximately true

[\frac{\partial Z_{w}}{\partial x}, \frac{\partial Z_{w}}{\partial y}] \sim \frac{(G - G_{h})}{F}

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

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.

4. 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 [x] = \frac{1}{2} {(x - b)}^{T} B^{- 1} (x - b) + \frac{1}{2} {(y - H [x])}^{T} R^{- 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 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

a = b + B H^{T} [b] {(H [b] B H^{T} [b] + R)}^{- 1} (y - H [b])

(22)

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.

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

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

4.2. Experiment with Real Observations

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

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

4.2.3. Results

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.

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.

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

Author Contributions

Conceptualization, F.Z.; Formal analysis, F.Z.; Investigation, F.Z. and M.K.; Resources, J.D. and P.V.; Software, F.Z., J.D. and P.V.; Visualization, M.K.; Writing—original draft, F.Z., J.D. and M.K.; Writing—review & editing, G.D. and J.W.

Funding

The mobility of M.K. at GFZ Potsdam was funded by the EU ESIF project No. CZ.02.2.69/0.0/0.0/16_027/0008463.

Acknowledgments

The Global Forecast System (GFS) data are provided by the National Centres for Environmental Prediction (NCEP). Precise GNSS orbits and clocks are provided by the European Space Agency (ESA). Absolute antenna phase centre offsets and variations are provided by the International GNSS Service (IGS). The four reviewers are gratefully acknowledged for their comments which helped to improve the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The horizontal delay gradient vector can be written as

G = {(Γ^{T} Γ)}^{- 1} Γ^{T} S - {(Γ^{T} Γ)}^{- 1} Γ^{T} S_{0}

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} = (\begin{matrix} m_{g} (e_{1}) \cos (a_{1}) & \dots & m_{g} (e_{n}) \cos (a_{n}) \\ m_{g} (e_{1}) \sin (a_{1}) & \dots & m_{g} (e_{n}) \sin (a_{n}) \end{matrix}) (\begin{matrix} S_{0} (e_{1}) \\ \dots \\ S_{0} (e_{n}) \end{matrix}) = (\begin{matrix} \sum_{j = 1}^{n} S_{0} (e_{j}) m_{g} (e_{j}) \cos (a_{j}) \\ \sum_{j = 1}^{n} S_{0} (e_{j}) m_{g} (e_{j}) \sin (a_{j}) \end{matrix})

where n denotes the number of tropospheric delays. For the considered 10 times 12 tropospheric delays (the elevation angles are 3°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, 90° and the spacing in azimuth is 30°) the two components of the vector can be written as

\sum_{j = 1}^{n} S_{0} (e_{j}) m_{g} (e_{j}) \cos (a_{j}) = \sum_{i = 1}^{10} \sum_{k = 1}^{12} S_{0} (e_{i}) m_{g} (e_{i}) \cos (a_{k}) = \sum_{i = 1}^{10} S_{0} (e_{i}) m_{g} (e_{i}) \sum_{k = 1}^{12} \cos (a_{k})

\sum_{j = 1}^{n} S_{0} (e_{j}) m_{g} (e_{j}) \sin (a_{j}) = \sum_{i = 1}^{10} \sum_{k = 1}^{12} S_{0} (e_{i}) m_{g} (e_{i}) \sin (a_{k}) = \sum_{i = 1}^{10} S_{0} (e_{i}) m_{g} (e_{i}) \sum_{k = 1}^{12} \sin (a_{k})

For the chosen observation geometry each tropospheric delay is accompanied by a tropospheric delay in opposite direction. Since cos(a) = −cos(a + 180°) and sin(a)= −sin(a + 180°) we have

\begin{matrix} \sum_{k = 1}^{12} \cos (a_{k}) = 0 \\ \sum_{k = 1}^{12} \sin (a_{k}) = 0 \end{matrix}

such that the two components of the vector are equal to zero.

References

Bevis, M.; Businger, S.; Herring, T.A.; Rocken, C.; Anthes, R.A.; Ware, R. GPS meteorology: Remote sensing of atmospheric water vapor using the Global Positioning System. J. Geophys. Res. 1992, 97, 15787–15801. [Google Scholar] [CrossRef]
Chen, G.; Herring, T.A. Effects of atmospheric azimuthal asymmetry on the analysis of space geodetic data. J. Geophys. Res. 1997, 102, 20489–20502. [Google Scholar] [CrossRef] [Green Version]
Poli, P.; Moll, P.; Rabier, F.; Desroziers, G.; Chapnik, B.; Berre, L.; Healy, S.B.; Andersson, E.; El Guelai, F.-Z. Forecast impact studies of zenith total delay data from European near real-time GPS stations in Meteo France 4DVAR. J. Geophys. Res. 2007, 112, D06114. [Google Scholar] [CrossRef]
Bennitt, G.V.; Jupp, A. Operational assimilation of GPS zenith total delay observations into the Met Office Numerical Weather Prediction models. Mon. Weather Rev. 2012, 140, 2706–2719. [Google Scholar] [CrossRef]
Mahfouf, J.-F.; Ahmed, F.; Moll, P.; Teferle, F.N. Assimilation of zenith total delays in the AROME France convective scale model: A recent assessment. Tellus A Dyn. Meteorol. Oceanogr. 2015, 67. [Google Scholar] [CrossRef]
Guerova, G.; Jones, J.; Douša, J.; Dick, G.; de Haan, S.; Pottiaux, E.; Bock, O.; Pacione, R.; Elgered, G.; Vedel, H.; et al. Review of the state of the art and future prospects of the ground-based GNSS meteorology in Europe. Atmos. Meas. Tech. 2016, 9, 5385–5406. [Google Scholar] [CrossRef] [Green Version]
Bar-Sever, Y.E.; Kroger, P.M.; Borjesson, J.A. Estimating horizontal gradients of tropospheric path delay with a single GPS receiver. J. Geophys. Res. 1998, 103, 5019–5035. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef]
Douša, J.; Dick, G.; Kačmaří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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
Thayer, G.D. An improved equation for the radio refractive index of air. Radio Sci. 1974, 9, 803–807. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
Boehm, J.; Heinkelmann, R.; Schuh, H. Short Note: A global model of pressure and temperature for geodetic applications. J. Geod. 2007, 81, 679–683. [Google Scholar] [CrossRef]
Kačmaří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. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
Daley, R. Atmospheric Data Analysis; Cambridge Atmospheric and Space Science Series; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
Parrish, D.; Derber, J. The National Meteorological Center’s spectral statistical analysis system. Mon. Weather Rev. 1992, 120, 1747–1763. [Google Scholar] [CrossRef]
Douša, J.; Václavovic, P.; Zhao, L.; Kačmaří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. [Google Scholar] [CrossRef]
Václavovic, P.; Douša, J. Backward smoothing for precise GNSS applications. Adv. Space Res. 2015, 56, 1627–1634. [Google Scholar] [CrossRef]
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. [Google Scholar]
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. [Google Scholar] [CrossRef]
Rebischung, P.; Altamimi, Z.; Ray, J.; Garayt, B. The IGS contribution to ITRF2014. J. Geod. 2016, 90, 611–630. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]

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

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.

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

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

Figure 6. Same as Figure 3 but data from the station network are assimilated.

Figure 7. Same as Figure 4 but data from the station network are assimilated.

Figure 8. Same as Figure 5 but data from the station network are assimilated.

Figure 9. 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. GNSS ZTDs and horizontal delay gradients from the G-Nut/Tefnut software 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).

Figure 10. 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. GNSS ZTDs and horizontal delay gradients from the G-Nut/Tefnut software 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).

Figure 11. 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. GNSS ZTDs and horizontal delay gradients from the G-Nut/Tefnut software 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).

Figure 12. The root mean square error of the refractivity as a function of the pressure for the background (black line), when we assimilate ZTDs only (blue line) and when we assimilate both ZTDs and horizontal delay gradients (red line). GNSS ZTDs and horizontal delay gradients from the G-Nut/Tefnut software 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). For details refer to the text.

© 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/).

Share and Cite

MDPI and ACS Style

Zus, F.; Douša, J.; Kačmařík, M.; Václavovic, P.; Dick, G.; Wickert, J. Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation. Remote Sens. 2019, 11, 41. https://doi.org/10.3390/rs11010041

AMA Style

Zus F, Douša J, Kačmařík M, Václavovic P, Dick G, Wickert J. Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation. Remote Sensing. 2019; 11(1):41. https://doi.org/10.3390/rs11010041

Chicago/Turabian Style

Zus, Florian, Jan Douša, Michal Kačmařík, Pavel Václavovic, Galina Dick, and Jens Wickert. 2019. "Estimating the Impact of Global Navigation Satellite System Horizontal Delay Gradients in Variational Data Assimilation" Remote Sensing 11, no. 1: 41. https://doi.org/10.3390/rs11010041

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu