A New Method for Estimating Tropospheric Zenith Wet-Component Delay of GNSS Signals from Surface Meteorology Data

Abstract

A new concept is proposed for estimating the zenith wet delay (ZWD) and atmospheric weighted average temperature by inputting the temperature, total pressure, and specific humidity from surface weather data. In addition, a new ZWD integral method is described for highly accurate calculation of the ZWD from radiosonde observation. To evaluate the advantages of the new discrete integral formula, we utilized the 8-year radiosonde profiles of 85 stations in China from 2010 to 2017 to validate the accuracy of the radiosonde-derived ZWD. The results showed that the mean accuracy of the ZWD derived from radiosonde data was 4.28 mm. Next, the new ZWD model was assessed using two sets of reference values derived from radiosonde data and GNSS precise point positioning in China. The results confirm that the new development improved the accuracy of the estimation of the tropospheric wet delay from the surface meteorological data. The performance of this new model can be seen as an important step toward accurately correcting the tropospheric delay in Global Navigation Satellite System (GNSS) real-time navigation and positioning. It can also be used in GNSS meteorology for weather forecasting and climate research.

1. Introduction

The propagation delay induced in signals from the Global Navigation Satellite System (GNSS) through the neutral atmosphere is a major error source in GNSS real-time navigation and positioning [1]. Neutral-atmosphere propagation delay is usually called tropospheric delay and is considered to be composed of two components [2]. The hydrostatic term, which depends on the dry gases of the atmosphere, can be accurately estimated by empirical models [3]. The non-hydrostatic term (namely the wet term), caused by the water vapor pressure, is a relatively small component of the tropospheric delay, but is difficult to determine because of the challenge of accurately quantifying the distribution of water vapor pressure in the atmosphere [3,4].

To achieve a high-quality solution from GNSS observation in high-precision applications, the tropospheric delay must be compensated for appropriately. In the past few decades, many hydrostatic and non-hydrostatic delay models have been built to obtain a correction for the neutral-atmosphere propagation delay for geodetic space systems [5,6,7]. The hydrostatic component usually contributes more than 90% to the total tropospheric delay and varies smoothly both temporally and spatially as the dry gases are well mixed [5]. Many studies have confirmed that the zenith hydrostatic delay (ZHD) can be well estimated with an accuracy of a few millimeters by using surface meteorological data based on the Saastamoinen ZHD model [7,8,9,10].

In comparison with the hydrostatic delay, the non-hydrostatic term is highly variable in time and space, as the water vapor pressure is inhomogeneous. Hence, the delay at the zenith (zenith wet delay, ZWD) is very difficult to estimate accurately based on an empirical model. Therefore, since the 1970s, many ZWD models have been proposed to derive the ZWD from surface meteorology observations [3,11,12,13,14,15,16,17]. The performances of these empirical ZWD models have been comprehensively and systematically investigated in specific regions or globally using benchmark values derived from radiosonde data or products of the European Center for Medium-Range Weather Forecasts [7,18]. The results show that the performances of these models still have large inaccuracies (approximately 60 mm) in some areas, so it is necessary to establish a new model to compensate for the poor ZWD estimation accuracy [7].

This study proposes a new model to estimate the ZWD from meteorological observations. Assuming that the atmosphere satisfies the ideal-state equation and the specific humidity decreases with height at a much faster rate than does the total pressure, the wet-component refractivity can be integrated from the surface up to the height where the water vapor pressure is near to zero. With the surface temperature, surface specific humidity, and surface total pressure as the only inputs, the ZWD can be estimated. Then, using the radiosonde-derived ZWD and GNSS precise point positioning (PPP)-derived ZWD as references, the performance of the new ZWD model was evaluated in China.

In the following section, we introduce the method to derive the new ZWD model and describe a new integral model for estimating the ZWD from radiosonde observations. An introduction of the radiosonde measurement and derivation of the mixing ratio of the atmosphere are presented in the third section. The performance of the new ZWD model with respect to radiosonde data and in GNSS PPP is discussed in the fourth section. Conclusions and final remarks are given in the last section.

2. Methods

The ZWD is typically obtained by one of two methods: first, mathematical integration of radiosonde data, radio occultation data, numerical weather forecast products and other key data; or second, establishment of a high-precision ZWD model. In general, ZWD models can be assessed using benchmark values derived from the first method.

2.1. Empirical Zenith Wet Delay (ZWD) Models

Many empirical models are used to estimate the ZWD including those of Saastamoinen, Hopfield, Black, Baby, Goad–Goodman, UNB3–EGNOS, Berman, Callahan, Chao and Askne–Nordius [3,11,12,13,14,15,16,19,20,21,22]. These ZWD models have been assessed using benchmark values derived from radiosonde data, and the accuracy of their ZWD estimations has been confirmed to be within several centimeters [4,7]. However, the ZWD values derived from the Callahan and Askne–Nordius models are more accurate than those of the other models at low latitudes of China (especially at 15°N–25°N). If no surface meteorological observations are available in GNSS real-time PPP, the best choice to compute the ZWD in China is the Askne–Nordius model [7].

2.1.1. Callahan Model

This model, developed by Callahan (1973) [15], assumes that the water vapor pressure P_w can be expressed through an empirical exponential as a function of height:

$$P_w=P_w^s·exp\left(-0.248·h-0.048·h^2\right)$$

(1)

where $P_w^s$ is the surface water vapor pressure and h is the height, expressed in hPa and km, respectively. The complete derivation of the model is very complicated, but the simplified final form can be obtained by using a set of nominal values:

$$\mathrm{ZWD}=\frac{1035·P_w^s}{T_s^2}$$

(2)

From Equation (2), it can be seen that the ZWD can be obtained by inputting the surface temperature T_s and the water vapor pressure.

2.1.2. Askne and Nordius

Askne and Nordius (1987) [3] assumed that the change in water vapor pressure in the vertical direction could be approximately derived from the total pressure using an exponential decay model:

$$\left(\frac{P_w}{P_w^s}\right)={\left(\frac{P}{p_s}\right)}^{\lambda +1}$$

(3)

where P and $P_s$ are the total pressure and the total pressure at the surface, respectively, and λ is the exponential decay parameter, which is dependent on season and latitude. Therefore, the ZWD can be obtained from

$$\mathrm{ZWD}={10}^{-6}\left[\left(k_2-k_1\frac{R_d}{R_w}\right)+\frac{k_3}{T_m}\right]·\frac{R_d·P_w^s}{\left(\lambda +1\right)·g_s}$$

(4)

where $k_1$, $k_2$, and $k_3$ are refractivity constants; $R_d$ and $R_w$ are the mean specific gas constant for daytime air and water vapor, respectively; and $g_s$ is the acceleration due to gravity at the surface. The weighted mean temperature is calculated by

$$T_m=T_s·\left(1-\frac{\beta ·R_d}{g_s·\left(\lambda +1\right)}\right)$$

(5)

where $T_s$ is the surface absolute temperature, and β is the lapse rate of temperature with height.

2.2. New ZWD Model

The propagation of electromagnetic waves in the neutral atmosphere is mainly subject to the refractivity, which can be divided into two terms: the hydrostatic term and the wet term. The first term can be expressed as a product of the pressure, temperature, and empirical constants of the refractivity. The latter term can be expressed as a function of the water vapor pressure, temperature, and refractivity constants. A two-term expression for the wet-component refractivity N^W can be written as follows [9,22]:

$$N^{\mathrm{w}}=\left(k_2-k_1\frac{R_d}{R_w}\right)·\frac{P_w}{T}{+k}_3·\frac{P_w}{T^2}$$

(6)

where $k_1$, $k_2$, $k_3$, $R_d$, and $R_w$ are as defined in Equations (3) and (4), and T is the absolute temperature.

Then, the ZWD can be defined as an integral of N^W in a profile above the station [2,9]:

$${\mathrm{ZWD}=10}^{-6}\left[\left(k_2-k_1\frac{R_d}{R_w}\right)·{\int }_{h_s}^{\infty }\frac{P_w}{T}\mathrm{dh}+k_3·{\int }_{h_s}^{\infty }\frac{P_w}{T^2}\mathrm{dh}\right]$$

(7)

By introducing the mean temperature ($T_m$) of the water vapor

$$T_m=\frac{{\int }_{h_s}^{\infty }\frac{P_w}{T}\mathrm{dh}}{{\int }_{h_s}^{\infty }\frac{P_w}{T^2}\mathrm{dh}}$$

(8)

We may express the ZWD as

$$\mathrm{ZWD}={10}^{-6}·\left[\left(k_2-k_1\frac{R_d}{R_w}\right)+\frac{k_3}{T_m}\right]·{\int }_{h_s}^{\infty }\frac{P_w}{T}\mathrm{dh}$$

(9)

Under the condition of hydrostatic equilibrium, we can use the hydrostatic equation:

$$\frac{dP}{dh}=-{\rho }_m\left(h\right)·g\left(h\right)$$

(10)

where P is the total pressure; ${\rho }_m$ is the density of the air; g is the acceleration due to gravity; and h is the height of the observation station above the mean sea level.

The specific humidity q can be expressed as a function of the density of water vapor ${\rho }_v$ and the density of air ${\rho }_m$:

$$q=\frac{{\rho }_v}{{\rho }_m}$$

(11)

${\rho }_v$ from Equation (11) can usually be expressed as

$${\rho }_v=\frac{P_w}{R_w·T}$$

(12)

Under the assumption that the gravity is a constant, 𝑔_𝑠, combining Equation (10), (11), and (12), we get

$$\mathrm{ZWD}={10}^{-6}·R_w·\left[\left(k_2-k_1\frac{R_d}{R_w}\right)+\frac{k_3}{T_m}\right]{\int }_0^{P_s}\frac{q}{g_s}\mathrm{dP}$$

(13)

It can be seen from Equation (13) that the ZWD can be derived after finding the integral ${\int }_0^{P_s}\frac{q}{g_s}\mathrm{dP}$. The value of this integral can be determined as follows.

The specific humidity and pressure usually decrease with height, but the latter does much more rapidly. Empirically, this variation can be described as follows [3,23]:

$$\left(\frac{q}{q_s}\right)={\left(\frac{P}{p_s}\right)}^{\omega }$$

(14)

where $\omega$ is the mixing ratio of the atmosphere and $q_s$ and $P_s$ are the surface specific humidity and surface total pressure, respectively.

The relationship between q and P can be established based on Equation (14). Then, substituting q and P into Equation (13) and evaluating the integral, we obtain

$$\mathrm{ZWD}={10}^{-6}·R_w·\left[\left(k_2-k_1\frac{R_d}{R_w}\right)+\frac{k_3}{T_m}\right]·\frac{q_s·P_s}{\left(\omega +1\right)·g_s}$$

(15)

The hydrostatic equation may be expressed as

$$\frac{dP}{P}=-\frac{g_s}{T·R_d}\mathrm{dh}$$

(16)

From Equation (16), we can derive

$$\frac{P_w·\mathrm{dP}}{P}=-\frac{P_w·g_s}{T·R_d}\mathrm{dh}$$

(17)

By introducing the temperature into Equation (17), we get

$$\frac{P_w·\mathrm{dP}}{P·T}=-\frac{P_w·g_s}{T^2·R_d}\mathrm{dh}$$

(18)

Then, the integral of ${\int }_{h_s}^{\infty }\frac{P_w}{T^2}\mathrm{dh}$ can be expressed as

$${\int }_{h_s}^{\infty }\frac{P_w}{T^2}\mathrm{dh}=-\frac{R_d·P_s^{\frac{-\beta ·R_d}{g_s}}}{g_{s·}{T}_s}{\int }_{P_s}^0\frac{P_w}{P^{\left(1-\frac{\beta ·R_d}{g_s}\right)}}dP$$

(19)

Inserting Equation (18) into Equation (19) (for more specific steps, see Appendix A), we get

$${\int }_{h_s}^{\infty }\frac{P_w}{T^2}\mathrm{dh}=\frac{R_d·P_w^s}{g_{s·}{T}_s·\left(\left(\omega +1\right)+\frac{\beta ·R_d}{g_s}\right)}$$

(20)

Then, the T_m can finally be obtained (for more specific steps, see Appendix A):

$$T_m=\frac{R_w}{R_d}·\frac{q_s·P_s}{P_w^s}·T_s·\left(1+\frac{\beta ·R_d}{g_s·\left(\omega +1\right)}\right)$$

(21)

According to Equation (15) and Equation (21), the new ZWD model estimates the ZWD by inputting $q_s$, $P_s$, and $T_s$ after obtaining $\omega$. The parameter $\omega$ may be chosen to fit the site in question, and it can be derived based on the ZWD (see Section 3.3).

2.3. New ZWD Integral Model

The vertical decrease in temperature with height, assuming a constant β decay rate, can be expressed by

$$T=T_s+\beta \left(h-h_s\right)$$

(22)

Considering Equation (16) and integrating yields

$$\left(\frac{P}{P_s}\right)={\left(\frac{T}{T_s}\right)}^{\frac{-g_s}{\beta ·R_d}}$$

(23)

Inserting Equations (3), (22), and (23) into Equation (7) and integrating for the interval [h₁,h₂], the ZWD is expressed as follows (for more specific steps, see Appendix B):

$${\mathrm{ZWD}}_{h_1}^{h_2}={10}^{-6}·\left\{\frac{\left(k_2-k_1\frac{R_d}{R_w}\right)·P_w^1}{E^{1,2}·{\beta }^{1,2}}\left[{\left(\frac{T_2}{T_1}\right)}^{E^{1,2}}-1\right]+\frac{k_3·P_w^1}{\left(E^{1,2}-1\right)·T_1·{\beta }^{1,2}}\left[{\left(\frac{T_2}{T_1}\right)}^{E^{1,2}-1}-1\right]\right\}$$

(24)

where $E^{1,2}=\frac{-\left({\lambda }_{h_1}^{h_2}+1\right)·g_s}{{\beta }_{h_1}^{h_2}·R_d}$.

The ZWD values from the radiosonde measurements can be derived as

$$\begin{array}{c}\mathrm{ZWD}={\sum }_{i=1}^{n-1}\frac{{10}^{-6}·\left(k_2-k_1\frac{R_d}{R_w}\right)·P_w^i}{E^{i,i+1}·{\beta }^{i,i+1}}\left[{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}}-1\right]\\ +{\sum }_{i=1}^{n-1}\frac{{10}^{-6}·k_3·P_w^i}{\left(E^{i,i+1}-1\right)·T_i·{\beta }^{i,i+1}}\left[{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}-1}-1\right]\end{array}$$

(25)

where $E^{i,i+1}=\frac{-\left({\lambda }_{h_i}^{h_{i+1}}+1\right)·g_s}{{\beta }_{h_i}^{h_{i+1}}·R_d}$. The ZWD obtained using Equation (25) from the radiosonde data is usually used as a reference to evaluate the accuracy of the ZWD obtained by other means.

3. Deriving the Mixing Ratio of the Atmosphere

From Equation (15), it can been seen that the mixing ratio of the atmosphere $\omega$ is vital in obtaining the ZWD. This key variable in the new ZWD model needs to be determined in advance.

3.1. Data Collection

The radiosonde sounding technique is an important means for meteorological study, offering high accuracy, high vertical resolution, long-term stability, and all-weather capability for obtaining the water vapor density, pressure, and temperature in the troposphere and lower stratosphere. China has approximately 120 operable radiosonde stations. From 1950 to 2010, three generations of radiosonde sensors (type 49, type 59-7-1, and GTS1) were used at these stations [24,25]. The replacement of sensors introduced inhomogeneity into the historical records, leading to a systematic bias in the homogenized time series. To ensure consistency within the radiosonde measurements and ensure the quality of the ZWD derived from radiosonde, only radiosonde data after 2009 were used in our study. Radiosonde measurements are freely available online from the University of Wyoming (http://weather.uwyo.edu/upperair/sounding.html, last accessed: October 14, 2019), although only the measurements of 85 of the stations in China can be downloaded from the above URL (uniform resource locator). The distribution of the 85 stations across China is shown in Figure 1.

3.2. Uncertainty of the Radiosonde-Derived ZWD

To assess the uncertainty of the ZWD obtained from radiosonde measurements using the new integral model, we differentiated Equation (25):

$$d_{\mathrm{ZWD}}=A_1·d_{P_w}+A_2·d_T$$

(26)

where $d_{P_w}$ and $d_T$ are the derivatives of the water vapor pressure and temperature, respectively. $A_1$ and $A_2$ are given by

$$\begin{array}{c}{A}_1={10}^{-6}·\left(k_2{-k}_1\frac{R_d}{R_w}\right)·{\sum }_{i=1}^{n-1}\frac{1}{E^{i,i+1}·{\beta }^{i,i+1}}\left[{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}}-1\right]+{10}^{-6}·k_3·\\ {\sum }_{i=1}^{n-1}\frac{1}{\left(E^{i,i+1}-1\right)·T_i·{\beta }^{i,i+1}}\left[{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}-1}-1\right];\end{array}$$

$$\begin{array}{l}{A}_2={10}^{-6}·\left(k_2{-k}_1\frac{R_d}{R_w}\right)·{\displaystyle \sum _{i=1}^{n-1}}\frac{P_w^i}{{\beta }^{i,i+1}}·{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}-1}·\frac{T_i-T_{i+1}}{T_i^2}+{10}^{-6}·k_3\\ ·{\displaystyle \sum _{i=1}^{n-1}}\frac{P_w^i}{\left(E^{i,i+1}-1\right)·T_i^2·{\beta }^{i,i+1}}[\left(1-{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}-1}\right)+\left(E^{i,i+1}-1\right)\\ ·{\left(\frac{T_{i+1}}{T_i}\right)}^{E^{i,i+1}-2}·\frac{1}{T_i}·\left(T_i-T_{i+1}\right)]\end{array}$$

The water vapor pressure can be obtained from P, T, and the relative humidity RH:

$$\begin{array}{c}{P}_w=6.1121\left(1.0007+3.46\times {10}^{-6}·P\right)·RH·\hfill \\ \hfill e^{\left\{\frac{\left[18.729-\left(T-273.15\right)/227.3\right]\left(T-273.15\right)}{T-15.28}\right\}}\end{array}$$

(27)

which yields

$$d_{P_w}=B_1·d_P+B_2·d_T+B_3·d_{RH}$$

(28)

where d_P and d_RH are the derivatives of the total pressure and relative humidity; $B_1$, $B_2$, and $B_3$ are given by $B_1=21.147866\times {10}^{-6}·RH·e^{\theta }$, $B_2=6.1121\times RH·\left(1.0007+3.46\times {10}^{-6}·P\right)·e^{\theta }·\frac{18.729-2\times \left(T-273.15\right)/227.3-\theta }{T-15.28}$, $B_3=6.1121\times \left(1.0007+3.46\times {10}^{-6}·P\right)·e^{\theta }$, and $\theta =\frac{\left[18.729-\left(T-273.15\right)/227.3\right]·\left(T-273.15\right)}{T-15.28}$.

Assuming that the measurements of pressure, temperature, and relative humidity are uncorrelated, and considering Equation (28), we can rewrite Equation (26) as

$$d_{\mathrm{ZWD}}=D_1·d_P+D_2·d_T+D_3·d_{RH}$$

(29)

where $D_1$, $D_2$, and $D_3$ are given by $D_1=A_1·B_1$; $D_2=A_1·B_2+A_2$; and $D_3=A_1·B_3$.

Finally, we can obtain the uncertainty of the ZWD based on the error propagation law, which is denoted as ${\sigma }_{{\Delta }_{\mathrm{ZWD}}}$:

$${\sigma }_{{\Delta }_{\mathrm{ZWD}}}=\sqrt{D_1·{\sigma }_P^2+D_2·{\sigma }_T^2+D_3·{\sigma }_{RH}^2}$$

(30)

To calculate the empirical values of $D_1$, $D_2$, and $D_3$ for each station, a large volume of radiosonde observations was collected from the 85 available stations in China over an 8-year period from 2010 to 2017. Here, we only give their maximum, minimum, and mean values in

Table 1. Uncertainty values estimated from the radiosonde measurements.

Parameters	Max	Min	Mean
$D_1$	3.95 × 10−3	4.62 × 10−4	1.48 × 10−3
$D_2$	16.88	0.31	8.18
$D_3$	37.91	3.46	25.62
${\sigma }_{{\Delta }_{\mathrm{ZWD}}}$ (mm)	8.65	0.23	4.28

. The sensor type GTS1 met the accuracy requirements in the Commission for Instruments and Methods of Observations Guide. The specified accuracies of P, T, and $RH$ for the radiosonde sensors were ${\sigma }_P$ = 1 hPa, ${\sigma }_T$ = 0.5 K, and ${\sigma }_{RH}$ = 0.05, respectively [26]. By substituting these parameters into Equation (30), the ${\sigma }_{{\Delta }_{\mathrm{ZWD}}}$ values were calculated and are listed in Table 1. The statistical results showed that the uncertainty of the ZWD values derived from radiosonde using the new integral model was within 4.28 mm. Therefore, it is feasible to use these values as a reference to evaluate the new ZWD model.

3.3. The Determination of the Mixing Ratio of the Atmosphere (ω)

Equation (15) is difficult to solve directly using least squares. Instead, the best method to obtain the optimal $\omega$ is based on the following model:

$$\sqrt{\frac{{\sum }_{i=1}^{num}{\left({\mathrm{ZWD}}_i^R-{\mathrm{ZWD}}_i\right)}^2}{num}}=min$$

(31)

where ZWD^R is the ZWD obtained from radiosonde based on Equation (25), num is the number of radiosonde data points used to estimate the ZWD, and ZWD_i is the ZWD obtained from the new model. First, we give $\omega$ a search range in advance, and when a certain value for $\omega$ satisfies Equation (31), it is considered as the optimum mixing ratio of the atmosphere.

In our study, the ZWD^R values from 2010 to 2017 were obtained based on Equation (25), then inserted into Equation (31), where the search range of $\omega$ was set to [1,5] and the step size of $\omega$ was 0.01. When a certain $\omega$ value satisfies Equation (31), it is considered as the optimal $\omega$. The ZWD is especially sensitive to the variation in $\omega$, so the latter should be chosen to fit the location and season. The $\omega$ values obtained in different seasons (spring, summer, autumn, and winter) are shown in Figure 2. In addition, the distribution of $\omega$ across China was derived from these values based on the Gauss distance-weighting function [27,28], as shown in the four subgraphs in Figure 3.

Figure 2 and Figure 3 illustrate that the ω coefficients varied by more than ±2 in some areas, depending on the season. In addition, ω was larger in northwest than southeast China. The values of ω in winter were slightly higher than those in other seasons at most stations, especially in western China. Moreover, this parameter is especially sensitive to latitude. In addition, according to the seasonal variability of ω over China, the water vapor pressure decreased with height much more rapidly than the total pressure, and the difference was more obvious in winter. China covers a wide range of latitudes, from 15°N to 55°N. The statistical properties of ω in four latitudinal ranges, namely, 15–25°N, 25–35°N, 35–45°N, and 45–55°N, were studied separately.

Table 2. Summary of the empirical values of coefficient $\omega$ in China across different latitudinal ranges.

	Season				Annual Mean ω
	ω (Spring)	ω (Summer)	ω (Autumn)	ω (Winter)
15–25	3.12	2.57	2.51	2.80	2.76
25–35	2.90	2.64	2.59	3.01	2.80
35–45	2.92	2.89	2.84	3.29	3.01
45–55	2.91	3.11	3.12	2.88	3.02

lists the empirical coefficient ω in the four latitudinal ranges of China in different seasons, and the annual means. The annual mean of ω becomes continually larger as the latitude increases. The same trend of increasing ω with latitude is also followed in summer and autumn, and the values for those two seasons gradually become closer to the annual mean as the latitude increases. In spring, ω was highest in the latitudinal range of 15–25°N, and greater at 35–45°N than at either 25–35°N or 45–55°N.

4. Results, Validation, and Discussion

The radiosonde observations from all 85 of the available stations in China for the entire year of 2018 were utilized to evaluate the performance of the new ZWD model with statistical validity. In addition, we further evaluated the performance of the new model using the GNSS PPP technique. The PPP method allows for high-precision estimation of the zenith tropospheric delay (ZTD). By subtracting the zenith hydrostatic delay (ZHD), which can be accurately retrieved from models, the ZWD can be derived. The PPP is especially suitable for testing the performance of the new ZWD model as the ZWD value derived from each GNSS station is independent.

4.1. Assessment of ZWD Model with Radiosonde

The performance of the new ZWD model was evaluated by using as a reference the ZWD values retrieved using Equation (25) from the measurements collected at the 85 radiosonde stations during 2018. In addition, the accuracies of two empirical models (Askne–Nordius and Callahan) were assessed using benchmark values derived from radiosonde data based on Equation (25). In addition, the University of New Brunswick (UNB) series of models (UNB1–UNB4) has been extensively used in several regions around the world. The relatively recent version UNB3m has a higher overall accuracy than the UNB3 and UNB4 models [29,30,31]. Therefore, the performance of UNB3m was also investigated using radiosonde data and compared to that of the new model. Figure 4 shows the standard deviation (STD) and root mean square (RMS) errors of the differences between the ZWD derived from radiosonde measurements and those of both the new ZWD and empirical models.

In Figure 4, RAD means the Equation (25)-derived ZWD using radiosonde data, NZWD means the new ZWD model-derived ZWD, ANZWD means the Askne–Nordius model-derived ZWD, CAZWD means the Callahan model-derived ZWD, UNZWD means the UNB3m model-derived ZWD, and DZWD means the deviation of the models-derived ZWD from RAD. Figure 4 illustrates that the deviations between RAD and NZWD for the 85 stations were smaller than those between RAD and ANZWD and between RAD and CAZWD. The STD values were smaller in northwestern than southeastern China. In addition, the RMS of the deviations between RAD and NZWD was better than 3 cm at high latitudes in China, but poorer at low latitudes. This can be explained by the water vapor concentration being greater in the southeast than the northwest of China. Additionally, we investigated the performance of the ZWD models relative to the radiosonde values in the different latitudinal regions of China, as shown in Figure 5.

Figure 5 shows the RMS, STD, and mean ZWD values calculated from the radiosonde profiles and the three models for different latitudinal regions. It is obvious that the deviations of the ZWD models at high latitudes were less than at low latitudes. The comparatively poor performance at low latitudes can be explained by the fact that high-latitude regions contain less water vapor in the atmosphere. The RMS is less than 2.5 cm at latitudes above 35°. In addition, it can be seen that the best and worst-performing ZWD models were NZWD and CAZWD, respectively. China spans 15°N–55°N, and accurate ZWD models must be robust against considerable seasonal variations. Thus, a radiosonde-based evaluation of the ZWD models was performed for each individual season, as listed in

Table 3. Seasonal statistics of deviations between RAD and model-derived ZWD (Unit: cm).

Model		Spring	Summer	Autumn	Winter
RAD-NZWD	RMS	3.32	4.47	3.41	2.03
	STD	3.17	4.35	3.32	1.92
	Mean	−0.43	−0.94	−0.34	−0.07
RAD-ANZWD	RMS	3.83	5.21	4.32	2.56
	STD	3.67	5.05	3.92	2.51
	Mean	−0.68	−1.39	−0.42	−0.39
RAD-CAZWD	RMS	4.34	5.57	4.74	2.91
	STD	4.01	5.25	4.59	2.77
	Mean	−1.88	−2.02	−2.32	−0.93
RAD-UNZWD	RMS	4.23	5.34	4.45	2.67
	STD	3.88	5.18	4.02	2.44
	Mean	−1.76	−1.39	−2.11	−0.34

.

The statistical summary of the deviations between RAD and the models-derived ZWD in different seasons indicates that the mean RMS of the deviation between RAD and NZWD was better than 4 cm, and even better than 2.5 cm in winter. In addition, all of the ZWD models experienced the largest RMS and STD errors in summer. This can be explained by the fact that the water vapor concentration is usually highest in this season. In addition, the errors of the three models in spring were very similar to those in autumn. It is evident that the new ZWD model outperformed the other models. The overall statistical results showed that the ZWD accuracy obtained by the new model was 16.8% better than that obtained by the Askne–Nordius model and 22.3% better than that obtained by the Callahan model.

4.2. Validation of ZWD by GNSS PPP

In the previous paragraph (§ 4.1), the ZWD obtained from the radiosonde profiles was used as a reference to evaluate the performance of the new ZWD model. However, the accuracy of radiosonde products mainly depends on the quality of the radiosonde sensor. To explore the absolute accuracy of the new model, we further validated the ZWD models using the GNSS-derived ZWD as the benchmark value (as seen in Figure 6). GNSS observations from 28 stations in China were processed to estimate the ZTD with Bernese Software Version 5.2 at a data sampling rate of 30 s for the entire year of 2018 [32]. Then, the ZHD was obtained using the Saastamoinen ZHD model using observed surface pressure values, and the ZWD was computed by removing the ZHD from the ZTD.

In Figure 6, GZWD is the GNSS PPP-derived ZWD; the definitions of the other abbreviated variables (NZWD, ANZWD, CAZWD, UNZWD, and DZWD) were the same as before. Figure 6 gives the RMS of the deviations between GZWD and the model-derived ZWD. The RMS of the deviations between RAD and the models-derived ZWD decreased as the latitude increased. The RMS errors were smaller in northwestern than southeastern China. Moreover, the RMS of the differences between GZWD and NZWD was better than 4 cm at high latitudes and better than 5.5 cm at low latitudes. The errors of the three models at high-latitude areas were very close, but the new model was significantly more accurate than the other two traditional models at low latitudes. Additionally, we used the GNSS observations to separately investigate the performance of the ZWD models in each of our four latitudinal regions in China (see Figure 7).

Figure 7 gives the RMS, STD, and mean values calculated from GNSS PPP and the three models in the four latitudinal regions. It can be observed that the latitudinal trends of the deviations from the reference values were similar for all three ZWD models. In addition, the new model again achieved the best performance in all four regions. As each GNSS station also experiences considerable seasonal variations, the detailed statistics of the seasonal differences between the GZWD and the models-derived ZWD are listed in

Table 4. Seasonal statistics of deviations between GZWD and models-derived ZWD (Unit: cm).

Model		Spring	Summer	Autumn	Winter
GZWD-NZWD	RMS	3.37	4.32	3.55	1.95
	STD	3.29	4.27	3.46	1.85
	Mean	−0.52	0.79	0.29	0.04
GZWD-ANZWD	RMS	3.92	5.14	4.23	2.42
	STD	3.84	4.83	3.88	2.38
	Mean	−0.71	−1.73	−1.02	0.49
GZWD-CAZWD	RMS	4.39	5.27	4.75	2.35
	STD	4.21	5.16	4.56	2.28
	Mean	−0.98	−0.22	−1.32	−0.56
GZWD-UNZWD	RMS	4.38	5.08	4.39	2.63
	STD	4.08	4.65	3.81	2.45
	Mean	−1.22	−1.63	−0.97	−0.53

.

The statistical results showed that the RMS and STD of the three models in winter were significantly smaller than in the other seasons, whereas these error terms were largest in summer. This can be explained by the fact that summer and winter had the highest and lowest water vapor content, respectively. In addition, the RMS error of the new model was only 4.32 cm in summer, which was significantly better than the other two models. The overall statistical results showed that the ZWD accuracy obtained by the new model was 15.4% better than that obtained by the Askne–Nordius model and 21.5% better than that obtained by the Callahan model. Considering its performance, we concluded that the new model was the optimal ZWD model for China.

Finally, the four ZWD models were applied to correct the wet delays in the GNSS PPP data. The data for a total of 31 days in July 2018 were collected to conduct the PPP experiment based on the GAMP (GNSS Analysis software for multi-constellation and multi-frequency precise positioning) software in kinematic mode. Then, the accuracies of the kinematic solutions were evaluated using the precise coordinates provided by the International GNSS Service for four different stations (BJFS, TWTF, WUHN, and URUM). The statistical accuracies in the three coordinate directions are shown in Figure 8.

In Figure 8, the “New model” is the new ZWD model, the “AN model” is the Askne–Nordius model, the “CA model” is the Callahan model, and the “UN model” is the UNB3m model. The three sub-graphs in Figure 8 illustrate the RMS values for the east, north, and up components, respectively. With the use of the new ZWD model, the kinematic PPP accuracy was approximately 3–9 cm in the horizontal and 15–35 cm in the vertical direction for the four tested stations. The maximum deviation of the kinematic PPP results obtained by the four ZWD models was approximately 1.5 cm in the horizontal and 3.5 cm in the vertical directions. The new ZWD model outperformed the other three models.

5. Conclusions

In space geodesy, especially in precise GNSS real-time positioning as performed today, it is always important to minimize the impact of tropospheric delay to obtain the most scientifically accurate solutions. The tropospheric delay can be divided into two terms: the hydrostatic and wet terms. The former satisfies the hydrostatic equation and the zenith hydrostatic delay can be accurately corrected based on the Saastamoinen ZHD model. The latter mainly depends on the water vapor pressure, which is difficult to calibrate. In this study, we proposed a new ZWD model to improve the estimation of the ZWD from meteorological observations. In addition, a new ZWD integral model based on the radiosonde profile was suggested. The mean uncertainty of the ZWD obtained in this way from the radiosonde data was 4.28 mm when utilizing the 8-year radiosonde profiles of 85 stations in China from 2010 to 2017. Thus, the ZWD derived from radiosonde using the new ZWD integral model could be used to evaluate the performance of the new ZWD model.

The new model was finally assessed using radiosonde data and GNSS observations as references. The observations of 85 radiosonde stations and datasets of 28 GNSS stations in China were collected for the entire year of 2018. Using the ZWD from radiosonde products based on Equation (25) as references, the ZWD results obtained from the new method exhibited RMS values of 4.47 cm and 2.03 cm in summer and winter, respectively. In addition, using the GNSS PPP-derived ZWD as references, the RMS errors of the ZWD values from the new model were 4.32 cm and 1.95 cm in summer and winter, respectively. The statistical results confirm that this is the optimal model for estimating the ZWD in China. In addition, it also shows that the new model compares better to the radiosonde-derived ZWD than the Askne-Nordius, Callahan, and UNB3m models. The accuracy of the new model for the Earth as a whole has yet to be verified, which we will address in our future research.