An empirical zenith wet delay correction model using piecewise height functions

Tropospheric delay is an important error source in space geodetic techniques. The temporal and spatial variations of the zenith wet delay (ZWD) are very large and thus limit the accuracy of tropospheric delay modelling. Thus, it is worthwhile undertaking research aimed at constructing a precise ZWD model. Based on the analysis of vertical variations of ZWD, we divided the troposphere into three height intervals (below 2 km, 2 to 5 km, and 5 to 10 km) and determined the fitting functions for the ZWD within these height intervals. The global empirical ZWD model HZWD, which considers the periodic variations of ZWD with a spatial resolution of $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, is established using the ECMWF ZWD profiles from 2001 to 2010. Validated by the ECMWF ZWD data in 2015, the precision of the ZWD estimation in the HZWD model over the three height intervals are improved by 1.4, 0.9, and 1.2 mm, respectively, compared to that of the currently best GPT2w model (23.8, 13.1, and 2.6 mm). The test results from ZWD data from 318 radiosonde stations show that the root mean square error (RMSE) in the HZWD model over the three height intervals was reduced by 2 % (0.6 mm), 5 % (0.9 mm), and 33 % (1.7 mm), respectively, compared to the GPT2w model (30.1, 15.8, and 3.5 mm) over the three height intervals. In addition, the spatial and temporal stabilities of the HZWD model are higher than those of GPT2w and UNB3m.

1 Introduction

The radio waves experience propagation delays when passing through the neutral atmosphere (primarily the troposphere), which are known as the tropospheric delays. The tropospheric delay is one of the main error sources in space geodetic techniques. In the processing of the space geodetic data, the tropospheric delay along the propagation path is generally expressed as the product of zenith tropospheric delay (ZTD) and mapping function (MF). The ZTD is divided into a zenith hydrostatic delay (ZHD) and a zenith wet delay (ZWD) (Davis et al., 1985), and the ZHD can be accurately determined using pressure observations. Unlike the ZHD, the ZWD is difficult to calculate accurately due to the high spatio-temporal variation in water vapour. Its spatial distribution is characterized with a near-zonal dependency, with values varying from about 2 cm at high latitudes to about 35 cm near the Equator (Fernandes et al., 2013). The temporal variation pattern of ZWD is mainly characterized by the seasonal variability, including annual and semi-annual components (Jin et al., 2007; Nilsson et al., 2008). The high variabilities in ZWD make it the main factor influencing tropospheric delay correction.

Various methods and models are developed to estimate the ZWD. Ray-tracing uses the observations from radiosonde profiles (Davis et al., 1985; Niell, 1996) or numerical weather models (Hobiger et al., 2008; Nafisi et al., 2012) to calculate the ZWD. It can provide the most accurate ZWD corrections. Models such as those developed by Bevis et al. (1992, 1994) make use of single-layer parameters from atmospheric models, such as total column water vapour (TCWV) and temperature. While Stum et al. (2011) proposed a model that only uses TCWV. These models provide similar results to the study of Davis et al. (1985) (that uses 3-D parameters) but only at the level of the model orography to which the meteorological parameters refer to. As this orography may depart significantly from the actual surface, and the vertical variation of the ZWD is not well known, at a different elevation they possess errors associated with the uncertainty in the modelling of the ZWD height variation (Fernandes et al., 2013, 2014; Vieira et al., 2018). The traditional Saastamoinen model (1972) and Hopfield model (1971) approximate the ZWD with surface observations as temperature and water vapour pressure observations. Without the information about the vertical distribution of water vapour, the stability and reliability of their ZWD estimates are poor. Moreover, both models are highly dependent on meteorological data. The aforementioned models have the limitations of application in wide area augmentation and real-time navigation and positioning. Therefore, the empirical climatological models were proposed as required practical conditions. The RTCA-MOPS (2016), designed by the US Wide Area Augmentation System (Collins et al., 1996), estimates ZWD by using the latitude band parameters table. The modified RTCA-MOPS model – called UNB3m (Leandro, 2006) – uses relative humidity as a parameter instead of the water vapour pressure to calculate the ZWD, effectively improving the precision of ZWD estimation to 5.5 cm (Möller et al., 2014), but the model deviation is increased when the height exceeds 2 km (Leandro, 2006). The TropGrid model (Krueger et al., 2004, 2005) provides the meteorological parameters needed to calculate tropospheric delay in the form of a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid. The improved TropGrid2 model (Schüler, 2014) enhances the efficiency of ZWD calculation by directly modelling ZWD with the exponential function. Based on the GPT2 model (Lagler et al., 2013), the GPT2w model (Böhm et al., 2015) adds weighted mean temperature and a vapour pressure decrease factor realised as a global grid to estimate ZWD by using the Askne and Nordius formula (Askne and Nordius, 1987). The GPT2w model has the best precision of ZWD estimation (3.6 cm) compared to other commonly used models (Möller et al., 2014).

The water vapour changes rapidly with respect to height, and the trends in water vapour at different heights vary, so the wet delay with direct relation to water vapour has complex spatio-temporal variations in the vertical direction. Kouba (2008) proposed an empirical exponential model to account for the height dependency of ZWD, but it will only be applicable within the height below 1000 m. The aforementioned empirical models are all based on a fixed height (average sea level or surface height) and use only a single decrease factor to describe the variation of water vapour or wet delay with respect to height, which makes it difficult to allow for the vertical distribution differences in water vapour (or wet delay) in the upper troposphere. In the course of aircraft dynamic navigation and positioning, the zenith delay error will result in twice as many errors in the station height estimate (Böhm and Schuh, 2013). Thus, it is necessary to correct the wet delay at different heights, which is clearly difficult for the aforementioned models. Based on the analysis of the characteristics of the ZWD profile, an empirical ZWD model, named HZWD, is established based on three functions applicable within corresponding height intervals, and the model precision is verified by European Centre for Medium-Range Weather Forecast (ECMWF) reanalysis data as well as radiosonde data.

Figure 1 Water vapour pressure profile (a) and ZWD profile (b) at a grid point (0^∘ N, 0^∘ E) at 12:00 UTC on 1 January 2010.

2 Vertical variations of ZWD

ZWD is defined as the integral of the wet refractivity along the vertical profile above the station.

$$\begin{array}{}\text{(1)}& \mathrm{ZWD}={\mathrm{10}}^{-\mathrm{6}}\underset{H}{\overset{\mathrm{\infty }}{\int }}{N}_{\mathrm{w}}\mathrm{d}h={\mathrm{10}}^{-\mathrm{6}}\underset{H}{\overset{\mathrm{\infty }}{\int }}(k_{\mathrm{2}}^{\prime }{\displaystyle \frac{e}{T}}+k_{\mathrm{3}}{\displaystyle \frac{e}{T^{\mathrm{2}}}})\mathrm{d}h\end{array}$$

In Eq. (1), N_w is the wet refractivity, e is the water vapour pressure in hectopascals (hPa), T is the temperature in Kelvin, $k_{\mathrm{2}}^{\prime }$ is 22.1 K hPa⁻¹, and k₃ is 373 900 K² hPa⁻¹ (Bevis et al., 1994). It can be seen from Eq. (1) that ZWD changes with height, vapour pressure, and temperature. The ZWD will decrease with increasing height due to the shorter integral length. With the profiles of water vapour pressure and temperature, one can obtain the accurate ZWD by the ray-tracing method. However, in practical applications (e.g. aircraft navigation and positioning, wide area augmentation), we usually use empirical models for ZWD corrections due to the unavailability of meteorological data profiles. Therefore, it is necessary to develop an empirical ZWD model with high precision. The temperature roughly decreases linearly with increasing height in the troposphere, while the change in water vapour is more variable, so the water vapour is the main determinant of vertical variation of ZWD. In the following content, we used the meteorological data profile of ERA-Interim pressure levels provided by ECMWF to analyse the vertical variation characteristics of ZWD and explore a suitable fitting function capable of describing the changes in ZWD with respect to height.

Figure 2 ZWD vertical gradients profile (a) and linear fit with height below 2 km (b) at a grid point (0^∘ N, 0^∘ E) at 12:00 UTC on 1 January 2010.

ERA-Interim can provide data at 00:00, 06:00, 12:00, and 18:00 UTC daily with a spatial resolution of not more than $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ and 37 pressure levels (Dee et al., 2011). The highest level data come from a height of approximately 50 km, covering almost the entire troposphere and stratosphere. We used the temperature, the geopotential height, and the specific humidity provided by the ERA-Interim pressure level data, and the discretised form of Eq. (1), to calculate the ZWD for each level height (Böhm and Schuh, 2013).

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}{e}_i=q_i\times P_i/(\mathrm{0.622}+\mathrm{0.378}\times q_i)\\ N_{{\mathrm{w}}_i}=k_{\mathrm{2}}^{\prime }{\displaystyle \frac{e_i}{T_i}}+k_{\mathrm{3}}{\displaystyle \frac{e_i}{T_i^{\mathrm{2}}}}\\ \mathrm{ZWD}={\mathrm{10}}^{-\mathrm{6}}\sum _i^{\mathrm{36}}{\displaystyle \frac{N_{{\mathrm{w}}_i}+N_{{\mathrm{w}}_{i+\mathrm{1}}}}{\mathrm{2}}}\cdot \left(h_{i+\mathrm{1}}-h_i\right)\end{array}\right.\end{array}$$

In Eq. (2), q is the specific humidity in grammes per gramme (g g⁻¹), P is the pressure in hPa, $k_{\mathrm{2}}^{\prime }$ and k₃ are empirical constants same as Eq. (1), and h is the geopotential height in metres. From Eq. (2), we can see that the ZWD at a specific level height is the sum of the ZWD portions in all layers above the specific level height. Figure 1 shows the water vapour pressure and ZWD profiles at a grid point (0^∘ N, 0^∘ E) at 12:00 UTC on 1 January 2010. From Fig. 1, it can be seen that the downward trend in the water vapour pressure varies significantly with height, and the decrease factor is different across different height intervals. The changes in ZWD with respect to height are similar to that of the water vapour pressure with respect to height: the decay is fastest up to a few kilometres height and slows down with increasing height; the ZWD values are close to zero after 10 km. Zhao et al. (2014) showed that about 50 % of the water vapour content is concentrated within 1.5 km of the surface and less than 10 % of the water vapour content remains above 5 km, leading to different ZWD decay rates within different height intervals. These results are basically consistent with our experiment results. Further, the derivative of the ZWD with respect to height (i.e. ZWD vertical gradient) is analysed to better understand the characteristic of the ZWD vertical distributions. Figure 2a shows the variation of ZWD vertical gradients with respect to height at the same grid point to Fig. 1. From Fig. 2a, it can be seen that the trends in ZWD vertical gradients at different height intervals are clearly different. Specifically, the linear fit of the ZWD gradients with height below 2 km shows a great agreement with an R² value of 0.99 (Fig. 2b). Thus we can come to a conclusion: ZWD gradients roughly change linearly below 2 km. From 2 to 5 km, and 5 to 10 km, the ZWD gradients vary non-linearly.

Figure 3 ZWD gradients profiles at grid points in different latitude bands (12:00 UTC, 1 January 2010).

Figure 4 ZWD gradients profiles at grid points in different latitude bands (12:00 UTC, 1 July 2010).

Figure 3 shows the ZWD vertical gradients with respect to height at grid points in different latitude bands. Figure 4 shows the similar ZWD vertical gradients as Fig. 3 but for a different season. The variations are similar to those in Fig. 2a, which show trend changes at about 2 and 5 km. It is worth noting that the ZWD gradients at low latitudes are much larger and water vapour is more variable than at high latitudes, resulting from the fact that the water vapour at low latitude is more variable. In addition, the ZWD gradient trends in the Southern Hemisphere are significant. In contrast, the ZWD gradients in the Northern Hemisphere are slightly complicated with respect to height: the reason for this may be that the Southern Hemisphere is mostly oceanic while the Northern Hemisphere has many sea coasts. The terrain complexity in the Northern Hemisphere contributes to the disturbances in the ZWD gradient in specific areas. According to the vertical variation characteristics of ZWD, we divided the troposphere into three height intervals (below 2 km, 2 to 5 km, and 5 to 10 km) and assumed 10 km to be the empirical tropopause beyond which the ZWD is assumed to be zero. For ZWD fitting with respect to height, TropGrid2 and GPT2w use exponential functions, while some scholars have also used a polynomial to describe the tropospheric delay with respect to height (Song et al., 2011). We used both polynomial and exponential functions to fit the variation trend of the ZWD with respect to height in the three selected intervals, respectively. The results showed that the quadratic polynomial used under 2 km and the exponential functions between 2 and 5 km and between 5 and 10 km gave the best fits. The combination of the quadratic polynomial and exponential functions for different height intervals are termed piecewise height functions. Table 1 summarises the global fitting statistics of different fit functions, demonstrating the superiority of piecewise height functions to the single polynomial function and single exponential function used for the whole troposphere.

Figure 5 Decadal time series and cycle fitting results of function coefficients z₁ (a), z₂ (b), and z₃ (c) at a grid point (0^∘ N, 0^∘ E) from 2001 to 2010.

Table 1 Fitting RMSE of piecewise height functions, single quadratic polynomial function, and single exponential function (unit: mm).

3 The HZWD model

From the above analysis of ZWD vertical variation and fitting, the piecewise height functions of the proposed HZWD model are:

$$\begin{array}{}\text{(3)}& \mathrm{ZWD}(B,L,H)=\left\{\begin{array}{l}{z}_{\mathrm{1}}+a_{\mathrm{1}}\cdot H+a_{\mathrm{2}}\cdot H^{\mathrm{2}}\\ H<\mathrm{2000}\mathrm{m}\\ z_{\mathrm{2}}\cdot \exp \mathit{\{}{\mathit{\beta }}_{\mathrm{2}}\cdot (H-\mathrm{2000}\left)\mathit{\right\}}\\ \mathrm{2000}\mathrm{m}\le H<\mathrm{5000}\mathrm{m}\\ z_{\mathrm{3}}\cdot \exp \mathit{\{}{\mathit{\beta }}_{\mathrm{3}}\cdot (H-\mathrm{5000}\left)\mathit{\right\}}\\ \mathrm{5000}\mathrm{m}\le H\le \mathrm{10}\mathrm{000}\mathrm{m}\\ \mathrm{0}\\ H>\mathrm{10}\mathrm{000}\mathrm{m}\end{array}\right.\end{array}$$

In Eq. (3), B is the latitude in degrees; L is the longitude in degrees; H is the height in metres; function coefficients z₁, z₂, and z₃ can be regarded as the ZWD at the height of 0, 2, and 5 km, respectively. We used the monthly mean profiles of ERA-Interim pressure levels from 2001 to 2010 with a horizontal resolution of $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ for ZWD modelling. The ZWD profiles calculated for each grid point are fitted by Eq. (3) to obtain the time series of the corresponding function coefficients: z₁, a₁, a₂, z₂, β₂, z₃, and β₃. It is worth noting that the ERA-Interim-derived ZWD data indicate that the averaged ZWD values at the three height intervals (i.e. below 2, 2 to 5 km, and 5 to 10 km) are 0.126, 0.0489, and 0.0111 m, respectively. Jin et al. (2007) found that the tropospheric delay has notable seasonal variations, mainly on annual and semi-annual cycles. Song et al. (2011) and Zhao et al. (2014) considered the temporal features of function coefficients in their troposphere models. We used the 10-year time series of those coefficients obtained to analyse their temporal variations. Figure 5 shows the time series and cycle fitting results of the function coefficients z₁, z₂, and z₃ at grid point (0^∘ N, 0^∘ E). Figure 5 shows that the time series of the function coefficients z₁, z₂, and z₃ have a significant characteristic annual cycle, and the semi-annual cycle is small but nevertheless evident.

Therefore, taking the annual and semi-annual cycles into consideration, we used Eq. (4) to fit the function coefficients derived from Eq. (3) to temporal parameters for each grid point (Böhm et al., 2015).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }r\left(t\right)=c_{\mathrm{0}}+c_{\mathrm{1}}\cos \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{2}\mathit{\pi }\right)+b_{\mathrm{1}}\sin \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{2}\mathit{\pi }\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }+c_{\mathrm{2}}\cos \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{4}\mathit{\pi }\right)+b_{\mathrm{2}}\sin \left({\displaystyle \frac{\mathrm{doy}}{\mathrm{365.25}}}\mathrm{4}\mathit{\pi }\right)\end{array}$$

In Eq. (4), c₀ is the annual mean, c₁ and b₁ are the annual cycle parameters, c₂ and b₂ are the semi-annual cycle parameters, and doy is the day of the year. The fittings show that the annual means and annual and semi-annual amplitudes of z₁, z₂, and z₃ are distinct. For instance, the cycle fitting results at a grid (0^∘ N, 0^∘ E) (Fig. 5) indicate that the temporal parameters (i.e. c₀, c₁, b₁, c₂, and b₂) of z₁ are 0.2911, 0.0237, 0.0312, −0.0006, and −0.0227 m, respectively; the temporal parameters of z₂ are 0.1215, 0.0118, 0.0203, 0.0004, and −0.0146 m, respectively; the temporal parameters of z₃ are 0.0255, 0.00031, 0.0070, −0.0019, and −0.0044 m, respectively. It should be noted that the fitting results of coefficients a₂, β₂, and β₃ show that all their annual means and annual and semi-annual amplitudes are small. However, below 2 km, the lack of cycle terms in a₂ would cause centimetre-level errors in the ZWD estimates, so these terms have been retained. For β₂ and β₃, ZWD itself is small at heights above 2 km, so the annual mean suffices for a desirable ZWD estimate. The experiment reveals that the loss of accuracy due to the lack of annual and semi-annual terms in β₂ and β₃ for the ZWD estimates is less than 0.1 mm. Therefore, only the annual means are retained for these two coefficients.

Figure 6 shows the global distributions of annual means of model coefficients z₁, z₂, and z₃. From Fig. 6 we can see that the extremum of ZWD annual means at 0 m height occur near the Equator and the maximum exceeds 0.36 m. The ZWD annual means decrease with increasing latitude. The distributions of ZWD annual means at 2 and 5 km heights are similar to that at 0 m, but the areas with the large values near the Equator decrease in extent and the ZWD distributions tend to be uniform, indicating that the water vapour content near the Equator is greater than that in other regions, and the ZWD value is also larger in low-altitude regions. As the height increases, the difference in water vapour content or ZWD between the Equator and other areas begins to decrease, but remains significant. Overall, there are some differences in the ZWD distribution at different heights, and it is necessary to model the spatio-temporal variations of ZWD at different heights.

Figure 6 Global distributions of annual means of HZWD model coefficients z₁ (a), z₂ (b), and z₃ (c).

After the fitting processes involving Eqs. (3) and (4), the global ZWD model HZWD, using piecewise height functions, is established. The spatial resolution of the HZWD model is $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$. Each grid point contains seven primary coefficients: z₁, a₁, a₂, z₂, β₂, z₃, and β₂. Among these coefficients, z₁, z₂, z₃, a₁, and a₂ are further expressed by Eq. (4) with five temporal parameters, respectively. Therefore, there are 27 parameters for each grid point and total 68 094 parameters for the HZWD model. As a comparison, the GPT2w model has a number of 77 760 parameters, which is 14 % more than that of the HZWD model. It is worth noting that the UNB3m model only has 50 parameters due to its coarse spatio-temporal resolution. When the HZWD model is applied, the four grid points surrounding the station are determined according to the horizontal position (latitude and longitude) of the station, and then the model coefficients of the corresponding height intervals at the four selected points are calculated according to Eq. (4). The ZWD of the four grid points are extrapolated to the station height by using Eq. (3), and finally the ZWD at the station location is obtained by using bilinear interpolation. The HZWD model only needs time, latitude, longitude, and height as input parameters. It can calculate ZWD without meteorological data and can provide wet delay correction products for navigation and positioning at different heights.

4 Validation and analysis of the HZWD model

To test the precision of the HZWD model and analyse the model correction performance compared to other troposphere models, we used the ERA-Interim pressure level data and radiosonde data from the year 2015 as external data sources and compared the results with the commonly used models UNB3m and GPT2w. The parameters used for the validation are bias and root mean square error (RMSE), expressed as follows.

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n({\mathrm{ZWD}}_i^M-{\mathrm{ZWD}}_i^{\mathrm{0}})\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n({\mathrm{ZWD}}_i^M-{\mathrm{ZWD}}_i^{\mathrm{0}}{)}^{\mathrm{2}}}\end{array}$$

In Eqs. (5) and (6), ZWD${}_i^M$ is the value estimated by the HZWD model developed in this study and ZWD${}_i^{\mathrm{0}}$ is the reference value.

For the UNB3m model, the ZWD at mean sea level (m.s.l.) is first calculated, then a vertical correction is applied to transform the ZWD to the target height. The formulae are as follows (Leandro et al., 2006):

$$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}{\mathrm{ZWD}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{6}}{\displaystyle \frac{(T_{\mathrm{m}}{k}_{\mathrm{2}}^{\prime }+k_{\mathrm{3}})R_{\mathrm{d}}}{g_{\mathrm{m}}(\mathit{\lambda }+\mathrm{1})-\mathit{\gamma }{R}_{\mathrm{d}}}}\cdot {\displaystyle \frac{e_{\mathrm{0}}}{T_{\mathrm{0}}}},\\ \mathrm{ZWD}={\mathrm{ZWD}}_{\mathrm{0}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\gamma }H}{T_{\mathrm{0}}}}\right)}^{\frac{(\mathit{\lambda }+\mathrm{1})g}{\mathit{\gamma }{R}_{\mathrm{d}}}-\mathrm{1}},\end{array}\right.\end{array}$$

where e₀, T₀, and ZWD₀ are the water vapour pressure, temperature, and ZWD at MSL, respectively; R_d is the specific gas constant for dry air (278.054 J kg⁻¹ K⁻¹); γ and λ are the temperature lapse rate and water vapour decrease factor, respectively; g_m is the gravity acceleration at the mass centre of the vertical column of the atmosphere and can be computed with geodetic latitude φ and height h by

$$\begin{array}{}\text{(8)}& g_{\mathrm{m}}=\mathrm{9.784}\left(\mathrm{1}-\mathrm{0.00266}\cos \mathrm{2}\mathit{\phi }-\mathrm{0.28}\cdot {\mathrm{10}}^{-\mathrm{6}}h\right);\end{array}$$

and T_m is the weighted mean temperature computed by

$$\begin{array}{}\text{(9)}& T_{\mathrm{m}}=T_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{\mathit{\gamma }{R}_{\mathrm{d}}}{g_{\mathrm{m}}(\mathit{\lambda }+\mathrm{1})}}\right).\end{array}$$

For the GPT2w model, the modelled meteorological parameters at the four grid points surrounding the target location are extrapolated vertically to the desired height, and then the Askne and Nordius formula (10) is used to calculate the wet delays at those base points: finally the wet delays are interpolated to the observation site in horizontal direction to get the target ZWD.

$$\begin{array}{}\text{(10)}& \mathrm{ZWD}={\mathrm{10}}^{-\mathrm{6}}\cdot (k_{\mathrm{2}}^{\prime }+{\displaystyle \frac{k_{\mathrm{3}}}{T_{\mathrm{m}}}})\cdot {\displaystyle \frac{R_{\mathrm{d}}e}{(\mathit{\lambda }+\mathrm{1})g_{\mathrm{m}}}}\end{array}$$

In GPT2w model, T_m is an empirical parameter modelled with seasonal components and g_m is simplified to a constant (9.80665 m s⁻²). It should be noted that the GPT2w model provides both $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ and $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ resolution versions. Since the horizontal resolution of the HZWD model is $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, we used the GPT2w model with the same resolution for validation.

4.1 Validation with ECMWF data

Modelling of the HZWD model is based on the monthly mean profiles of ERA-Interim pressure level data from 2001 to 2010, while we used the ERA-Interim pressure level data with the full time resolution of 6 h in 2015 for the model validation. This is to validate the model performance on the daily scale. Regarding the ZWD profiles calculated from these data as reference values, we calculated the global annual average bias and RMSE of the ZWD for three models (HZWD, GPT2w, and UNB3m) within the three height intervals: below 2 km, 2 to 5 km, and 5 to 10 km (Table 2).

Table 2 Error statistics for the three models compared to the 2015 ECMWF data (unit: mm).

From Table 2, it can be seen that the HZWD model is the most accurate model across all three intervals, followed by the GPT2w model, and the UNB3m model has the worst performance. The annual average biases of the HZWD model are lower than that of the GPT2w model and the UNB3m model, except below 2 km. Compared with the RMSEs in the GPT2w model, those of the HZWD model are decreased by 1.4, 0.9, and 1.2 mm within the three height intervals, corresponding to improvements of about 6 %, 6 %, and 32 %, respectively. The improvements of the HZWD model over GPT2w model will result in precision improvements of 2.8, 1.8, and 2.4 mm respectively in height estimates in real-time aircraft positioning. The correction performance improvement from 5 to 10 km height is particularly evident. Figure 7a shows the ECMWF ZWD profile and the ZWD profiles of the three models at 12:00 UTC on 1 January 2015 at a representative grid point (0^∘ N, 20^∘ E). More clearly, Fig. 7b shows the differences between the ZWD profiles of the three models and ECMWF ZWD profile at different heights. It can be seen that HZWD is the most stable model, showing the best agreement with the ECMWF ZWD data, which is superior to both the GPT2w and the UNB3m models.

Figure 7 The ZWD profiles (a) of ECMWF and the three models (HZWD, GPT2w, and UNB3m) and corresponding biases (b) at a grid point (0^∘ N, 20^∘ E) at 12:00 UTC on 1 January 2015.

The variation of the troposphere has a strong correlation with latitude. To analyse the correction performances of the three models in different regions around the world, we calculated the three models' errors in different latitude bands (10^∘ intervals). Figures 8 and 9 show the correction performances at different latitudes. It can be seen from Fig. 8 that the bias of the UNB3m model is basically positive in the three height intervals, indicating that its ZWD estimates are relatively large compared to the ECMWF data. Moreover, the bias in the Southern Hemisphere is significantly larger than that in the Northern Hemisphere, indicating systematic deviations in the Southern Hemisphere. Both the GPT2w model and the HZWD model have large biases in the low latitudes. The biases of the GPT2w model are positive from 2 to 5 km and 5 to 10 km height, indicating that the ZWD is overestimated by the GPT2w model with increasing height. For the HZWD model, the bias in each latitude band is relatively small with few exceptions, resulting in a global average bias close to zero (see Table 2).

Figure 8 Bias comparisons between the three models (HZWD, GPT2w, and UNB3m) in different latitude bands over the year 2015.

Figure 9 shows the RMSEs of the three models. It can be seen from Fig. 9 that the precision of the HZWD model is significantly better than that of the UNB3m model across the three height intervals and all latitude bands, which is better than the GPT2w model in general. The precision of the three models declines with decreasing latitude, because the active change in water vapour in these areas limits the precision of the model. Corresponding to Fig. 8, the errors in UNB3m are asymmetric: the main reason for this is that the meteorological parameters of UNB3m are interpolated from the coarse look-up table with a latitude interval of 15^∘, and UNB3m does not consider the longitudinal variations of any meteorological elements. It should be pointed out that the UNB3m model is based on the simple symmetric assumption of the Northern and Southern Hemispheres, and its modelling data source only comes from the atmospheric data collected over North America, which leads to poor precision in the Southern Hemisphere, especially in the high latitudes thereof.

Figure 9 RMSE comparisons between the three models (HZWD, GPT2w, and UNB3m) in different latitude bands over the year 2015.

Table 3 Error statistics for the three models validated by 2015 radiosonde data (unit: mm).

Summarising the distributions of bias and RMSE across different latitude bands, we can see that the HZWD model performs best with the ECMWF data as reference values. Compared with the models GPT2w and UNB3m, the HZWD model basically eliminates systematic error in the 5 to 10 km height interval and the correction performance is stable at all heights and regions. To investigate the model's performance over time, Fig. 10 shows the time series of RMSE for the three models at 6 h intervals throughout the year 2015 at grid point (0^∘ N, 20^∘ E). We can see that the HZWD model has the best overall performances within the three height intervals over the year 2015. We noticed the significantly large RMSE for all three models across all three height intervals around the doy 19 and doy 195 of 2015. This can be attributed to the sharp short-term ZWD variations in the Equator area. The short-term variations are hardly accounted for by all three models which only consider the seasonal variations of ZWD. Moreover, the GPT2w model has the worst performance from 5 to 10 km height, which is also identified by Figs. 8 and 9.

Figure 10 RMSEs in ZWD estimates of the three models (HZWD, GPT2w, and UNB3m) compared to the ECMWF data over the year 2015 at grid point (0^∘ N, 20^∘ E).

Figure 11 Uncertainty of ZWD with respect to height at station 01241 (63.70^∘ N, 9.60^∘ E at 10 m) at 11:00 UTC on 1 January 2015.

Figure 12 Global distributions of bias for the three models (HZWD, GPT2w, and UNB3m) compared to 2015 radiosonde data.

4.2 Validation with radiosonde data

A radiosonde is used in a sounding technique that regularly releases balloons to collect atmospheric meteorological data at different heights: it can obtain profiles of various meteorological data with high accuracy. At present, the Integrated Global Radiosonde Archive (IGRA) website (ftp://ftp.ncdc.noaa.gov/pub/data/igra/, last access: 1 November 2018) provides free downloads of global radiosonde data. We used radiosonde data from 318 stations collected in 2015 to test the HZWD model. After data pre-processing, the data with gross errors have been removed and a total of 163 671 radiosonde data epochs remained. With the provided profiles of geopotential height, temperature, and water vapour pressure, the data forms of the radiosonde data are very similar to the ECMWF pressure level data, and thus the radiosonde ZWDs can be calculated using the same method by Eq. (2). Before the validation, we conducted an assessment of the uncertainty of ZWD derived from radiosonde data. Rózsa (2014) showed that the uncertainty of ZWD is ±1.5 mm in the case of the Vaisala RS-92 radiosondes in central and eastern Europe. However, this uncertainty is only valid for the ZWD calculated from the height of the lowest layer and is limited to Europe. Using the same uncertainties of radiosonde meteorological data given by the technical specification of the radiosonde and the algorithm proposed by Rozsa (2014), we calculated the ZWD uncertainty for all heights in all radiosonde stations. Figure 11 shows the uncertainty of ZWD with respect to the height for radiosonde station 01241 located in Orland, Norway (63.70^∘ N, 9.60^∘ E at 10 m). We can see that the uncertainty of ZWD is less than ±1.5 mm near the height of 0 m and decrease quickly with increasing height. The global mean uncertainties of ZWD in all stations in the three height intervals are ±1.3, ±0.7, and ±0.2 mm, respectively, indicating the high accuracy of ZWD derived from radiosonde data.

Figure 13 Global distributions of RMSE for the three models (HZWD, GPT2w, and UNB3m) compared to 2015 radiosonde data.

Figure 14 RMSEs in ZWD estimates of the three models (HZWD, GPT2w, and UNB3m) for radiosonde station 01241 over the year 2015.

Taking the radiosonde ZWDs as reference ZWD values, we validated the ZWDs from the models HZWD, GPT2w, and UNB3m. Table 3 shows the statistical results of the three models. It can be seen from Table 3 that the HZWD model has the best overall stability of the average bias and RMSE, indicating the best precision, and the UNB3m model is the worst. Compared with the GPT2w model, the RMSEs in HZWD in the three height intervals are reduced by 0.6, 0.9, and 1.7 mm, which equates to precision improvements of 2 %, 5 %, and 33 %, respectively. Moreover, these improvements correspond to an error reduction of 1.2, 1.8, and 3.4 mm respectively in height estimates in geodetic techniques. Taking the uncertainty of radiosonde ZWD into account, the improvement of the HZWD model over GPT2w model below 2 km seems to be insignificant. Nevertheless, we can reasonably think that the ZWD predicted by HZWD is closer to true ZWD due to its smaller RMSE. It is worth noting that the bias and RMSE of the HZWD model and the GPT2w model are both larger than those of the results from ECMWF data in Table 2. The reason is that the HZWD model and the GPT2w model are based on ECMWF data, and thus the test results with radiosonde data are slightly worse than those using ECMWF data. On the contrary, the bias of the UNB3m model decreases, and the RMSE between 2 and 5 km and between 5 and 10 km are less than those in Table 2. It may be due to the fact that most of the radiosonde stations are in the Northern Hemisphere, accounting for more than 60 % (192∕318) of the total, which has a positive impact on the test results for the UNB3m model based on North American meteorological data.

Figure 12 shows the global distributions of bias for the three models within the three height intervals, and Fig. 13 shows the global distributions of RMSE for the three models. As can be seen from Fig. 12, the three models show a poorer performance in low-latitude areas than in mid- and high-latitude areas for all height intervals, similar to the results in Sect. 4.1. Within the 5 to 10 km interval, the bias of the GPT2w model is large and positive in the equatorial region, indicating that the ZWD of the GPT2w at this height is significantly overestimated, and the global bias of the UNB3m model in this height interval is positive, also indicating an overestimate of the ZWD in the UNB3m model. The bias of the HZWD model does not show obvious regional differences with respect to height, and the overall distribution of HZWD model bias has no tendency to either the positive or negative. Figure 13 further illustrates the precision of the HZWD model. The global RMSE distributions of the HZWD model are similar to that of GPT2w model below 2 km and between 2 and 5 km, but the precision of the HZWD model is slightly better. Combining this with the bias distribution of the GPT2w model in Fig. 12, the GPT2w model also has a large RMSE near the Equator in the 5 to 10 km interval, which shows that the GPT2w model is unstable at high height in low-latitude areas. The precision of the UNB3m model is poorer than that of both the HZWD and GPT2w models. Below 2 km, the UNB3m model reaches decimetre-level precision near the Equator, and even exceeds 12 cm in some areas: the distribution of north–south heterogeneity remains obvious.

These results validate the spatial stability of the precision of the HZWD model; furthermore the temporal stability of the model precision is verified next. Figure 14 shows the results of ZWD corrections of the three models for the radiosonde station 01241 for the whole of 2015. It can be seen from Fig. 14 that the HZWD model and the GPT2w model are relatively stable throughout the year, while the correction performance of the UNB3m model in 2015 is worse than those of the HZWD and GPT2w models. The probable reason for this is that the UNB3m model only takes into account the annual variations in the metrological elements with a fixed phase, resulting in precision instability throughout the year. The improved performance arising from use of the HZWD model, compared to that arising from use of the GPT2w model, is more apparent with increasing height: this shows that modelling ZWD piecewise with height can effectively approximate the real ZWD profile and improve the precision of ZWD estimation.

5 Conclusions

The complexity of spatio-temporal variations makes the modelling of tropospheric ZWD difficult. In this paper, the characteristics of vertical variation of wet delay are analysed. The troposphere is divided into three height intervals: below 2 km, 2 to 5 km, and 5 to 10 km according to different trends (10 km is assumed to represent the empirical tropopause). A quadratic polynomial and two exponential functions are used to describe the variation of wet delay within each of the three intervals. Based on the monthly mean data of ECMWF ZWD from 2001 to 2010, a global ZWD model with spatial resolution of $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ was established with height fitting followed by periodic fitting. Using the ECMWF ZWD data for 2015, the annual average RMSEs in the HZWD model are 23.8, 13.1, and 2.6 mm in the below 2 km, 2 to 5 km, and 5 to 10 km height intervals, respectively, which is far superior to the performance of the UNB3m model. Compared to the currently most accurate wet delay empirical model (the GPT2w model), the precisions within the three height intervals improved by 6 % (1.4 mm), 6 % (0.9 mm), and 32 % (1.2 mm), respectively. The improvements will result in precision improvements of 2.8, 1.8, and 2.4 mm respectively in height estimates in real-time aircraft positioning. The testing results of radiosonde data from 318 stations in 2015 show that the annual average RMSEs of the HZWD model are 30.1, 15.8, and 3.5 mm, which are 2 % (0.6 mm), 5 % (0.9 mm), and 33 % (1.7 mm) better than those of the GPT2w model, respectively, corresponding to height error reduction of 1.2, 1.8, and 3.4 mm in real-time aircraft positioning. Considering the ZWD fields (0.126, 0.0489, and 0.0111 m) in the three height intervals, the precision improvements at the top layer are especially significant, which accounts for about 15 % of the corresponding ZWD field. Moreover, compared with the GPT2w and UNB3m models, the HZWD model offers the highest spatio-temporal stability. With higher precision of ZWD estimates and fewer model parameters, the HZWD model is more efficient than the GPT2w model.

The HZWD model offers good precision stability in the vertical direction and can meet the requirements of ZWD correction at different heights within the troposphere; however, it can be seen that neither the HZWD nor the GPT2w models, i.e. those non-meteorological parameter-based models, performed well in the lowest region of the troposphere. In addition, compared with GPT2w, HZWD model is a closed model with a limitation to facilitate on-site meteorological observations. Further research is required to assess the variation and factors influencing the wet delay and to explore the possibility of incorporation of on-site meteorological data.