**Regular paper**
06 Nov 2018

**Regular paper** | 06 Nov 2018

# An empirical zenith wet delay correction model using piecewise height functions

YiBin Yao and YuFeng Hu

^{1,2,3}

^{1,4}

**YiBin Yao and YuFeng Hu**YiBin Yao and YuFeng Hu

^{1,2,3}

^{1,4}

^{1}School of Geodesy and Geomatics, Wuhan University, 129 Luoyu Road, Wuhan 430079, China^{2}Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, 129 Luoyu Road, Wuhan 430079, China^{3}Collaborative Innovation Center for Geospatial Technology, 129 Luoyu Road, Wuhan 430079, China^{4}College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, Shanxi, China

^{1}School of Geodesy and Geomatics, Wuhan University, 129 Luoyu Road, Wuhan 430079, China^{2}Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, 129 Luoyu Road, Wuhan 430079, China^{3}Collaborative Innovation Center for Geospatial Technology, 129 Luoyu Road, Wuhan 430079, China^{4}College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, Shanxi, China

**Correspondence**: YuFeng Hu (yfhu@whu.edu.cn)

**Correspondence**: YuFeng Hu (yfhu@whu.edu.cn)

Received: 22 May 2018 – Discussion started: 28 May 2018 – Revised: 04 Sep 2018 – Accepted: 18 Oct 2018 – Published: 06 Nov 2018

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.

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.

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

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^{−1}, and *k*_{3} is
373 900 K^{2} hPa^{−1} (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.

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

In Eq. (2), *q* is the specific humidity in grammes per gramme
(g g^{−1}), *P* is the pressure
in hPa, ${k}_{\mathrm{2}}^{\prime}$ and *k*_{3} 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*^{2} 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 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.

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

In Eq. (3), *B* is the latitude in degrees; *L* is the longitude in degrees;
*H* is the height in metres; function coefficients *z*_{1}, *z*_{2}, and
*z*_{3} 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*_{1}, *a*_{1}, *a*_{2}, *z*_{2}, *β*_{2},
*z*_{3}, and *β*_{3}. 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*_{1}, *z*_{2}, and *z*_{3} at grid point (0^{∘} N,
0^{∘} E). Figure 5 shows that the time series of the function
coefficients *z*_{1}, *z*_{2}, and *z*_{3} 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).

In Eq. (4), *c*_{0} is the annual mean, *c*_{1} and *b*_{1} are the annual
cycle parameters, *c*_{2} and *b*_{2} 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*_{1}, *z*_{2}, and *z*_{3} 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*_{0}, *c*_{1}, *b*_{1}, *c*_{2}, and *b*_{2}) of *z*_{1} are 0.2911,
0.0237, 0.0312, −0.0006, and −0.0227 m, respectively; the temporal
parameters of *z*_{2} are 0.1215, 0.0118, 0.0203, 0.0004, and −0.0146 m,
respectively; the temporal parameters of *z*_{3} 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*_{2}, *β*_{2}, and *β*_{3}
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*_{2} would cause
centimetre-level errors in the ZWD estimates, so these terms have been
retained. For *β*_{2} and *β*_{3}, 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 *β*_{2} and *β*_{3} 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*_{1}, *z*_{2}, and *z*_{3}. 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.

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*_{1}, *a*_{1}, *a*_{2}, *z*_{2}, *β*_{2}, *z*_{3}, and *β*_{2}. Among these coefficients, *z*_{1}, *z*_{2}, *z*_{3}, *a*_{1}, and *a*_{2} 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.

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.

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

where *e*_{0}, *T*_{0}, and ZWD_{0} 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^{−1} K^{−1}); *γ* 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

and *T*_{m} is the weighted mean temperature computed by

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.

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^{−2}). 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).

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.

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

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.

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

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.

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.

The ERA-Interim data can be downloaded from ECMWF (Dee et al., 2011). The radiosonde data can be downloaded from the IGRA (ftp://ftp.ncdc.noaa.gov/pub/data/igra/, last access: 1 November 2018) (NOAA, 2018).

YY and YH contributed to the conception of the study. YY contributed significantly to the data analysis and manuscript preparation. YH performed the model validation and wrote the manuscript.

The authors declare that they have no conflict of interest.

The authors would like to thank the ECMWF and IGRA
for providing relevant data. We thank the reviewers for their insightful
comments and constructive suggestions. This research was supported by the
National Key Research and Development Program of China (2016YFB0501803) and
the National Natural Science Foundation of China (41574028) and Key
Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan
University (16-02-03).

Edited by: Marc Salzmann

Reviewed by: Jakub Kalita and M. Joana Fernandes

Askne, J. and Nordius, H.: Estimation of tropospheric delay for microwaves from surface weather data, Radio Sci., 22, 379–386, 1987).

Bevis, M., Businger, S., Herring, T. A., Rocken, C., Anthes, R. A., and Ware, R. H.: GPS meteorology: Remote sensing of atmospheric water vapor using the Global Positioning System, J. Geophys. Res.-Atmos., 97, 15787–15801, 1992.

Bevis, M., Businger, S., Chiswell, S., Herring, T. A., Anthes, R. A., Rocken, C., and Ware, R. H.: GPS meteorology: Mapping zenith wet delays onto precipitable water, J. Appl. Meteorol., 33, 379–386, 1994.

Böhm, J. and Schuh, H. (Eds.): Atmospheric effects in space geodesy, Berlin, Springer, Vol. 5, 6–7, 2013.

Böhm, J., Möller, G., Schindelegger, M., Pain, G., and Weber, R.: Development of an improved empirical model for slant delays in the troposphere (GPT2w), GPS Solut., 19, 433–441, 2015.

Collins, P., Langley, R., and LaMance, J.: Limiting factors in tropospheric propagation delay error modelling for GPS airborne navigation, Proc. Inst. Navig. 52nd Ann. Meet, 3, http://gauss2.gge.unb.ca/papers.pdf/ion52am.collins.pdf, 1996.

Davis, J. L., Herring, T. A., Shapiro, I. I., Rogers, A. E. E., and Elgered, G.: Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length, Radio Sci., 20, 1593–1607, 1985.

Dee, D. P., Uppala, S. M., Simmons, A. J., et al.: The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, 2011.

Fernandes, M. J., Nunes, A. L., and Lázaro, C.: Analysis and inter-calibration of wet path delay datasets to compute the wet tropospheric correction for CryoSat-2 over ocean, Remote Sens., 5, 4977–5005, 2013.

Fernandes, M. J., Lázaro, C., Nunes, A. L., and Scharroo, R.: Atmospheric corrections for altimetry studies over inland water, Remote Sens., 6, 4952–4997, 2014.

Hobiger, T., Ichikawa, R., Koyama, Y., and Kondo, T.: Fast and accurate ray-tracing algorithms for real-time space geodetic applications using numerical weather models, J. Geophys. Res.-Atmos., 113, https://doi.org/10.1029/2008jd010503, 2008.

Hopfield, H. S.: Tropospheric effect on electromagnetically measured range: Prediction from surface weather data, Radio Sci., 6, 357–367, 1971.

Jin, S., Park, J. U., Cho, J. H., and Park, P. H.: Seasonal variability of GPS-derived zenith tropospheric delay (1994–2006) and climate implications, J. Geophys. Res.-Atmos., 112, https://doi.org/10.1029/2006jd007772, 2007.

Kouba, J.: Implementation and testing of the gridded Vienna Mapping Function 1 (VMF1), J. Geodesy, 82, 193–205, 2008.

Krueger, E., Schueler, T., Hein, G. W., Martellucci, A., and Blarzino, G.: Galileo tropospheric correction approaches developed within GSTB-V1, in: Proc. ENC-GNSS, Vol. 1619, 2004.

Krueger, E., Schueler, T., and Arbesser-Rastburg, B.:. The standard tropospheric correction model for the European satellite navigation system Galileo, Proc. General Assembly URSI, 2005.

Lagler, K., Schindelegger, M., Böhm, J., Krásná, H., and Nilsson, T.: GPT2: Empirical slant delay model for radio space geodetic techniques, Geophys. Res. Lett., 40, 1069–1073, 2013.

Leandro, R., Santos, M. C., and Langley, R. B.: UNB neutral atmosphere models: development and performance, Proc. ION NTM, 52, 564–573, 2006.

Möller, G., Weber, R., and Böhm, J.: Improved troposphere blind models based on numerical weather data, Navigation, 61, 203–211, 2014.

Nafisi, V., Madzak, M., Böhm, J., Ardalan, A. A., and Schuh, H.: Ray-traced tropospheric delays in VLBI analysis, Radio Sci., 47, 1–17, 2012.

Niell, A. E.: Global mapping functions for the atmosphere delay at radio wavelengths, J. Geophys. Res.-Sol. Ea., 101, 3227–3246, 1996.

Nilsson, T. and Elgered, G.: Long-term trends in the atmospheric water vapor content estimated from ground-based GPS data, J. Geophys. Res.-Atmos., 113, https://doi.org/10.1029/2008jd010110, 2008.

NOAA: Radiosonde data, available at: ftp://ftp.ncdc.noaa.gov/, last access: 1 November 2018.

Rózsa, S. Z.: Uncertainty considerations for the comparison of water vapour derived from radiosondes and GNSS[M]//Earth on the Edge: Science for a Sustainable Planet, Springer Berlin Heidelberg, 65–78, 2014.

RTCA: Minimum Operational Performance Standards For Global Positioning System/Satellite-Based Augmentation System Airborne Equipment, Radio Technical Commission for Aeronautics, SC-159, publication DO-229E, Washington, DC, 2016.

Saastamoinen, J.: Introduction to practical computation of astronomical refraction, Bull. Géodésique, 106, 383–397, 1972.

Schüler, T.: The TropGrid2 standard tropospheric correction model, GPS Solut., 18, 123–131, 2014.

Song, S., Zhu, W., Chen, Q., and Liou, Y.: Establishment of a new tropospheric delay correction model over China area, Sci. China Phys. Chem., 54, 2271–2283, 2011.

Stum, J., Sicard, P., Carrere, L., and Lambin, J.: Using objective analysis of scanning radiometer measurements to compute the water vapor path delay for altimetry, IEEE T. Geosci. Remote, 49, 3211–3224, 2011.

Vieira, T., Fernandes, M. J., and Lázaro, C.: Analysis and retrieval of tropospheric corrections for cryosat-2 over inland waters, Adv. Space Res., 62, 1479–1496, 2018.

Zhao, J. Y., Song, S. L., Chen, Q. M., Zhou, W. L., and Zhu, W. Y.: Establishment of a new global model for zenith tropospheric delay based on functional expression for its vertical profile, Chinese J. Geophys., 57, 3140–3153, 2014 (in Chinese).