4.1 Static global-gridded T_s–T_m model

A linear formula $T_{\mathrm{m}}=a{T}_{\mathrm{s}}+b$ for the relation between T_m and T_s has been adopted in many studies. Based on the T_s and T_m products from the ERA-Interim data covering the years 2009 to 2012, we performed linear fittings of T_m versus T_s on each grid point. Then, the slope constant (a), the intercept constant (b) and the fitting root mean square error (RMSE) of each linear expression were calculated and contoured in Fig. 4. The a and b values are related to the latitude as well as the underlying surface (e.g., land, ocean). In the middle–high latitudes over the Northern Hemisphere, constant a value varies from 0.6 to 0.8, and constant b is approximately 100–50 over most of the continents. The constants in the Bevis equation are within these value ranges. Constant a is smaller (approximately 0.5–0.7) over land at the middle–high latitudes over the Southern Hemisphere. In particular, there are abrupt changes in the values of constants a and b from land to ocean at the middle–high latitudes due to the different feature variations of T_s and T_m (see Fig. 2). At low latitudes, the a value is smaller than over the other regions, because of the low variations of T_s and T_m. The fitting RMSEs are within 2–4 K over the middle–high latitude lands, and lower values are obtained over the oceans or at low latitudes. The reason for the low RMSE around the equator is the smaller fluctuation of T_m. Meanwhile, there is no RMSE larger than 4.5 K in the results of our model. As we did not perform any spatial or temporal smoothing of the data during the data processing, both the precision and resolution of our static model are better than other models (e.g., Lan et al., 2016).

Figure 3T_s–T_m scatter plots at two locations: (blue dots) 0.35^∘ N 180.00^∘ E and (red dots) 70.53^∘ N 180.00^∘ E. The magenta and green lines are their linear fitting curves. Temperature unit is Kelvin.

Download

4.2 Time-varying global-gridded T_s–T_m model

The time variation in the T_s–T_m relationship should also be considered in a precise T_s–T_m model. Therefore, a time-varying equation is applied for T_s–T_m regression at each grid node:

where “doy” represents the observed day of year and “h” is the observed hour in UTC time; (m₁, m₂), (n₁, n₂) and (p₁, p₂) are fitting coefficients. These equations can reflect the amplitudes of annual, semiannual and diurnal variations in our T_s–T_m models.

Figure 4Distributions of the (a) slope constant a, (b) intercept constant b and (c) RMSE of static linear T_s–T_m equations at ERA-Interim grid nodes. Temperature unit is Kelvin.

Download

Our new regression model found similar values for the coefficients a and b (of its static term) as for the static model in Sect. 4.1, except for some differences over the oceans. In Fig. 5, besides these constants a and b, we illustrate the amplitudes of annual, semiannual and diurnal terms. We can see that there are large annual variations (amplitude > 5 K) in the vast regions from Tibet to northern Africa, and in some places of the Siberia and Chile. Large diurnal variations (amplitude > 3 K) mainly occur over the midlatitude lands such as northeastern Asia or North America. Semiannual variations, however, are small in most areas except some high-latitude areas (amplitude > 3 K). All variations are smaller over the oceans due to the slower temperature changes over water than over land. The estimated T_m RMSE is also contoured in Fig. 5, and we can see that the RMSE dropped significantly in the regions with large annual or diurnal variations.

Figure 5(a) The slope constant a, (b) intercept constant b, amplitudes of T_m (c) annual, (d) semiannual and (e) diurnal terms in our time-varying global gridded T_s–T_m model, and (f) the model-estimated T_m RMSE distribution. Temperature unit is Kelvin.

Download

4.3 Assessments of T_s–T_m models

To further assess the precision of the T_s–T_m models using other independent data sources, we generated T_m and T_s from the radiosonde data at 723 radiosonde stations in the year 2016. These data are not assimilated into the 2009–2012 ERA-Interim data sets. As a result, we can regard them to be independent of our model. At each radiosonde site, different T_s–T_m models were employed to calculate T_m. In addition, we also estimated T_m using the 1^∘ × 1^∘ GPT2w model (Bohm et al., 2015), which is a global gridded T_m empirical model independent of the surface meteorological observation data. Then, these calculated T_m will be evaluated by being compared with the integrated T_m of radiosondes (denoted as T_{m_RS}) twice a day.

The model estimations of T_m are denoted by T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w} for the Bevis equation, the latitude-related model, our static global gridded model, time-varying global gridded model and the GPT2w model. When the global gridded models are employed, the radiosonde station may not be located at a grid node. Therefore, we interpolated the coefficients in the T_s–T_m equations from the neighboring grids to the radiosonde sites. The interpolation formula is expressed as follows (Jade and Vijayan, 2008):

C_site and $C_{\mathrm{site}}^i$ represent the coefficients in T_s–T_m equations at the radiosonde site location and its neighboring grids, respectively. wⁱ are the interpolation coefficients, which are determined using the equation

where R=6378.17 km is the mean radius of the earth, λ is the scale factor which equals one in our study, and ψⁱ is the angular distance between the ith grid node and the station's position. ψⁱ are computed using following formula (with latitude φ and longitude θ):

Considering the fact that the reanalysis grids are definite and every radiosonde site is in situ, we can compute the interpolation coefficients in Eq. (7) for all of the radiosonde stations. Then, these coefficients are stored as constants to avoid reduplicating the calculation.

Taking T_{m_RS} as the reference values, we calculated the biases and RMSEs of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w} at each radiosonde site. The results are illustrated in Fig. 6. Obviously, in many regions, the Bevis equation has bad precision with the absolute bias and RMSE both larger than 5 K. T_{m_LatR} can reduce the estimated biases in many areas, but the RMSEs remain large. Large biases still exist at quite a few radiosonde stations, e.g., in Africa or western Asia. T_{m_static} and T_{m_GPT2w} remove the large T_m biases at most of the radiosonde stations. T_{m_varying} performs significantly better over the world, especially in the Middle East, North America, Siberia, etc.

Detailed statistics of the distributions of the bias and RMSE using different models are shown in Fig. 7 and Table 2. At over 97.37 % of the radiosonde stations, the biases of T_{m_varying} are within −3–3 K. Large positive biases (>3 K) nearly disappear in T_{m_varying}. In contrast, there are significant large biases in T_{m_Bevis} and T_{m_LatR}. Improvements in RMSE are more evident. The RMSEs of T_{m_varying} are smaller than 4 K at over 91 % of the radiosonde sites, while few sites (<1 %) have RMSEs larger than 5 K. This is clearly better than the other models. In T_{m_Bevis} and T_{m_LatR}, there are more than 17 % of the radiosonde sites with RMSEs larger than 5 K. The overall performance of T_{m_GPT2w} is very close to T_{m_Bevis}, except that its absolute bias is smaller than the other T_s–T_m models.

Table 2Statistics of T_m estimates from different models. Reference data are the radiosonde T_m derivations.

Download Print Version | Download XLSX

Figure 6(a, c, e, g, i) The biases and (b, d, f, h, j) the RMSEs of the estimated T_m from (a, b) the Bevis equation, (c, d) the latitude-related model, (e, f) our static global gridded model, (g, h) our time-varying global gridded model and (i, j) the GPT2w model at each radiosonde station. Reference data are the radiosonde data of the year 2016. Temperature unit is Kelvin.

Figure 7The distributions of (1) the biases and (2) the RMSEs of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w} compared with the radiosonde data at 723 stations in the year 2016. Temperature unit is Kelvin.

Download

Table 3Number of radiosonde sites at which the five global applied T_m estimation models had superior performance.

Download Print Version | Download XLSX

Figure 8T_m series of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying}, T_{m_GPT2w} and T_{m_RS} at the IGRA station (a) no. 62378 and (b) no. 40841. Temperature unit is Kelvin.

Download

To identify the superior T_m estimation model at each radiosonde site, we employed the following statistical tests under the assumption of a normal distribution of the estimated T_m error:

1. First, Brown–Forsythe tests (Brown and Forsythe, 1974) of the equality of variances were carried out at each site for estimating the T_m errors from two different models, e.g., model A and B. The purpose of this step is to determine whether there are significant differences in the variances of the T_m results. If the test rejects the null hypothesis at a 5 % significance level and the errors of model A and B have the same variance, the model with the smaller sample variance is regarded as the better one. However, if the test does not reject the homogeneity of variances, analysis of variance (ANOVA) is performed in the next step.

2. ANOVA is a technique used to analyze the differences among group means (Hogg and Ledolter, 1987). It evaluates the null hypothesis that the samples all have the same mean against the alternative that the means are not the same. If the null hypothesis is rejected at a 5 % significance level, the T_m sample with smaller absolute mean value is believed to be better. Otherwise, we think that the two models perform almost as well at this radiosonde site.

3. After multiple tests and comparisons, the best model at each radiosonde station may be identified. However, at some sites no superior model can be confirmed. All the models are believed to have equivalent performance.

Finally, we counted the number of sites at which each T_m model performed the best. The results are given in Table 3. The time-varying global gridded model is superior to the others at 434 radiosonde stations (60.03 % of all sites), while the second-best estimation, T_{m_GPT2w}, is superior at only 12.86 % of the sites.

In Fig. 8 the T_m series at the IGRA station no. 62378 (29.86^∘ N 31.34^∘ E in Egypt) are given. We can see that large negative biases ($<-\mathrm{5}$ K) between T_{m_Bevis} (or T_{m_LatR}) and T_{m_RS} exist. T_{m_static} performs only slightly better from July to October. However, T_{m_varying} and T_{m_GPT2w} can eliminate most of the seasonal errors. Different properties of T_m series appear at another IGRA station no. 40841 (30.25^∘ N 56.97^∘ E in Iran). Some observation data are missing, but we can still see that there are large positive differences (>5 K) between T_{m_Bevis} (or T_{m_LatR}) and T_{m_RS} throughout the year. The biases of T_{m_static} are much smaller, but some large errors still appear in many months. The T_{m_varying}, however, performs as well as the T_m calculated from the radiosonde data, with small biases, and captures the variations well. The time series ofT_{m_GPT2w} are smoother and cannot capture the fluctuations of the T_m time series, causing an accuracy worse than T_{m_varying}.

On the other hand, even T_{m_varying} have large differences from T_{m_RS} at a few IGRA stations. This can be explained by the fact that our fitting analyses are based on the T_m values derived from ERA-Interim profiles. The quality of ERA-Interim data can be very poor in the regions with sparse observation data (Itterly et al., 2018).

Figure 9Variation in the uncertainty in Q with the value of T_m and the uncertainty in T_m.

Download

5.1 Theoretical analysis of the GPS–PWV uncertainty

Comprehensive research on the uncertainty in GPS–PWV has been carried out by Ning et al. (2016). The uncertainties in ZTD, ZHD and conversion factor Q have been studied in detail. The total uncertainty in GPS–PWV is as follows:

where σ_PWV, σ_ZTD, ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ and σ_Q are, respectively, the uncertainties of GPS–PWV, the ZTD estimation, the P_s observations and the conversion factor Q. σ_c=0.0015 denotes the uncertainty in constant C=2.2767 in Eq. (1), PWV is the value of GPS–PWV, and

where ${\mathit{\sigma }}_{k_{\mathrm{3}}}=\mathrm{0.012}\times {\mathrm{10}}^{\mathrm{5}}$ K² hPa⁻¹, ${\mathit{\sigma }}_{k_{\mathrm{2}}^{\prime }}=\mathrm{2.2}$ K hPa⁻¹ and ${\mathit{\sigma }}_{T_{\mathrm{m}}}$, respectively, denote the uncertainties of k₃, $k_{\mathrm{2}}^{\prime }$ and T_m in Eq. (4). The variation in σ_Q with the values of T_m and ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ is depicted in Fig. 9. Assuming that T_m is 280 K, we find that the σ_Q increases by over 60 % (from 0.069 to 0.112) as the ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ rises from 3.0 to 5.0 K. However, the σ_Q is less sensitive to the value of T_m. The σ_Q rises only by 17.96 % (about from 0.061 to 0.075) as the value of T_m drops from 300 to 270 K with ${\mathit{\sigma }}_{T_{\mathrm{m}}}=\mathrm{3.0}$ K.

Ning et al. (2016) assumed the T_m were obtained from NWP models so the uncertainty in T_m was set to be small (${\mathit{\sigma }}_{T_{\mathrm{m}}}=\mathrm{1.1}$ K). However, as shown in Sect. 4.3, the uncertainties of T_m from different T_m models are significantly larger at the radiosonde stations. For each radiosonde station, we calculated the mean value of T_m and assigned the ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ with the RMSEs of T_m given in Fig. 6. Then we obtained the σ_Q in Eq. (11). Our statistics indicate that the σ_Q using our varying T_s–T_m model decreases on average by 19.26 %, 17.77 %, 7.79 % and 18.67 % with respect to the σ_Q, respectively, using T_{m_Bevis}, T_{m_LatR}, T_{m_static} and T_{m_GPT2w}. For example, at the IGRA station no. 42724 (22.88^∘ N 91.25^∘ E in India), σ_Q drops by 53 % from 0.141 of the T_{m_Bevis} to 0.066 of the T_{m_varying}.

The uncertainty in Q will be propagated to the total uncertainty in GPS–PWV according to Eq. (10). We obtained the contributions of the different terms in Eq. (10) to the total GPS–PWV uncertainty. The contribution of one term is measured by the percentage it accounts for the total σ_PWV. The percentages are computed using the formulas

where p_ZTD, $p_{P_{\mathrm{s}}}$, p_C and p_Q, respectively, indicate the contributions of the uncertainties associated with ZTD, P_s, constant C and factor Q to the total σ_PWV. Following the summaries of Ning et al. (2016), we assumed that σ_ZTD=4 mm and σ_C=0.0015. T_m identically equals 280 K since the σ_Q is less sensitive to the value of T_m with respect to the ${\mathit{\sigma }}_{T_{\mathrm{m}}}$. Table 4 gives five sets of typical values which are assigned to the ${\mathit{\sigma }}_{P_{\mathrm{s}}}$, ${\mathit{\sigma }}_{T_{\mathrm{m}}}$, P_s and PWV in Eqs. (10)–(12).

Table 4Different typical values for ${\mathit{\sigma }}_{P_{\mathrm{s}}}$, ${\mathit{\sigma }}_{T_{\mathrm{m}}}$, P_s and PWV.

Download Print Version | Download XLSX

Figure 10Contributions of different terms to the total uncertainty in GPS–PWV with the typical values shown in Table 4.

Download

The ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ equals 0.2 hPa in Ning et al. (2016); however we enlarged its typical value to 0.5 hPa in consideration of the possible worse performance of the surface barometers. In Fig. 10, we illustrated the contributions of the terms in Eq. (12) based on the assumptions (a)–(e) in Table 4. Some feature variations of the contributions of different terms can be found from the comparisons between different subplots:

1. No significant difference exists between Fig. 10a and b. Because of the small value of σ_c in Eq. (10), the σ_PWV is not sensitive to the value of P_s. Meanwhile, the uncertainty associated with σ_c contributes less than 10 % of the σ_PWV.

2. With the typical values in Table 4 (values a and b), a reduction of ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ can reduce the p_Q significantly. For example, in Fig. 10a, the p_Q accounts for 69.54 % with ${\mathit{\sigma }}_{T_{\mathrm{m}}}=\mathrm{6}$ K, and it declines to 38.19 % with ${\mathit{\sigma }}_{T_{\mathrm{m}}}=\mathrm{3}$ K.

3. As Fig. 10c shows, the uncertainty associated with σ_ZTD accounts for the main part of σ_PWV when the values of PWV and ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ are not high. With the typical values in Table 4 (value c), the p_ZTD can be up to 74.21 % with ${\mathit{\sigma }}_{T_{\mathrm{m}}}=\mathrm{3}$ K. The p_Q, however, can drop from 26.76 % to 9.00 % as the ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ decreases from 6 to 3 K. Although the p_Q is not large under this situation, a smaller ${\mathit{\sigma }}_{T_{\mathrm{m}}}$can still reduce the contribution of σ_Q to the σ_PWV.

4. The uncertainty associated with ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ dominates the error budget of PWV when the ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ is large. In Fig. 10d and e, the $p_{P_{\mathrm{s}}}$ is over 80 % with ${\mathit{\sigma }}_{T_{\mathrm{m}}}<\mathrm{3}$ K and ${\mathit{\sigma }}_{P_{\mathrm{s}}}=\mathrm{5}$ hPa. In Fig. 10d, the p_Q increases from 7.55 % to 23.19 % as the ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ rises from 3 to 6 K. However, in Fig. 10e, the p_Q only grows from 1.29 % to 4.61 % with the same variation in ${\mathit{\sigma }}_{T_{\mathrm{m}}}$.

Theoretical analyses of σ_PWV were also carried out at two representative stations. At the IGRA station no. 42971 (20.25^∘ N 85.83^∘ E, in India), the mean value of PWV is 53.88 mm. The RMSEs of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w} are 4.30, 3.15, 2.41, 1.93 and 1.97 K. The ${\mathit{\sigma }}_{T_{\mathrm{m}}}$ in Eq. (11) was replaced by the calculated RMSEs, and the p_ZTD, $p_{P_{\mathrm{s}}}$, p_C and p_Q were generated with two typical values, 0.5 and 5 hPa, assigned to the ${\mathit{\sigma }}_{P_{\mathrm{s}}}$. With ${\mathit{\sigma }}_{P_{\mathrm{s}}}=\mathrm{0.5}$ hPa, the p_C accounts for around 7 %, while the $p_{P_{\mathrm{s}}}$ accounts for around 4 % of the total σ_PWV. By using different T_m estimations, the variations of p_C and $p_{P_{\mathrm{s}}}$ are both within 4 %. However, the p_Q varies more evidently. It accounts for averages of 55.69 %, 40.77 %, 30.70 %, 23.53 % and 24.11 % of the σ_PWV with the estimations of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w}, respectively. The p_ZTD rises with the reduction of p_Q, e.g., from 36.23 % of T_{m_Bevis} to 62.53 % of T_{m_varying}. On the other hand, with ${\mathit{\sigma }}_{P_{\mathrm{s}}}=\mathrm{5}$ hPa, the $p_{P_{\mathrm{s}}}$ accounts for more than 75 % of the σ_PWV, while the p_Q decreases from 14.21 % of T_{m_Bevis} to 3.9 % of T_{m_varying}.

At another representative station, the IGRA station no. 50557 (49.17^∘ N 125.22^∘ E, in northeastern China), the mean PWV is only 12.17 mm. The RMSEs of T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w} are 5.16, 3.94, 3.54, 2.99 and 5.10 K. We can see that the accuracy of T_m has been improved significantly. However, because of the low average value of PWV, the p_ZTD averagely contributes over 73.5 % of the σ_PWV, while the p_Q averagely contributes less than 10.5 % assuming ${\mathit{\sigma }}_{P_{\mathrm{s}}}=\mathrm{0.5}$ hPa and less than 1.5 % assuming ${\mathit{\sigma }}_{P_{\mathrm{s}}}=\mathrm{5}$ hPa. But such a discussion only concerns the average values. In fact, even at this station there are still some high values of PWV, for example at 12:00 UTC 22 July 2016, the PWV reached 48 mm. For the observations with high PWV, the improvement in the accuracy of T_m can still exert a significant positive impact on the reduction of p_Q.

Table 5Statistics on the relative errors of different PWV retrievals.

Download Print Version | Download XLSX

It is worth mentioning that the uncertainty in ZHD may be underestimated in some situations. There are two reasons for this. Firstly, the calculation of ZHD assumes that the water vapor does not contribute to the mass of the atmosphere. The ZHD error introduced by this assumption is often negligible. But in some very wet regions, the mass of water vapor could produce significant errors in the ZHD calculation. Secondly, and more importantly, the error of P_s in Eq. (1) can sometimes be very large. Small ${\mathit{\sigma }}_{P_{\mathrm{s}}}$ is reasonable when the surface barometer is calibrated routinely and equipped together with the GPS antenna. However, if there were significant height difference between the GPS antenna and the barometer, the error for ZHD would increase significantly. Snajdrova et al. (2006) found that 10 m of height difference approximately causes a difference of 3 mm in the ZHD. On the other hand, P_s can be generated from NWP data if there are no nearby barometers at GPS site. The error of P_s could be very large using this method (Means and Cayan, 2013; Jiang et al., 2016). In these cases, the GPS–PWV error reduction will be very limited due to the more precise T_m estimation.

5.2 Impact of real T_m estimation

To study the impact of T_m on the real GPS–PWV retrieval, we first downloaded GPS ZTD products (Byun and Bar-Sever, 2009) at 74 IGS sites in the year 2016 from the NASA Crustal Dynamics Data Information System (CDDIS) ftp address (ftp://cddis.gsfc.nasa.gov/pub/gps/products/troposphere/zpd, last access: 25 February 2019). These selected GPS sites were equipped with meteorological sensors so that the surface pressure and temperature measurements could also be obtained. ZHD was calculated using Eq. (1). It is subtracted from ZTD to obtain ZWD. Then, T_m was generated with six approaches: the first five T_m series were T_{m_Bevis}, T_{m_LatR}, T_{m_static}, T_{m_varying} and T_{m_GPT2w}. The sixth T_m was integrated from the ERA-Interim profiles and interpolated to each GPS site (Jiang et al., 2016; Wang et al., 2016). Finally, the GPS–PWV was generated from the ZWD and the six different T_m estimates leading to over 100 compared points for each GPS–PWV series. We denoted these GPS–PWV sets as PWV_BTm, PWV_LTm, PWV_STm, PWV_VTm, PWV_GTm and PWV_ETm. The only difference between these GPS and PWV estimations is the T_m estimation model; therefore, the impact of other errors is excluded.

Figure 11Relative RMSEs of PWV_BTm, PWV_STm, PWV_VTm and PWV_GTm compared with PWV_ETm at 74 IGS stations in the year 2016.

Download

Figure 12(a) PWV differences and (b) relative differences of PWV_BTm, PWV_LTm, PWV_STm, PWV_VTM and PWV_GTm compared with PWV_ETm at the ALIC station in the year 2016. PWV unit is millimeters.

Download

The T_m from ERA-Interim is believed to be the most accurate among our T_m estimates at the selected GPS sites. We therefore took the PWV_ETm as reference values to assess the other PWV. The relative RMSEs of PWV_BTm, PWV_LTm, PWV_STm, PWV_VTm and PWV_GTm at these selected stations were calculated and are illustrated in Fig. 11. Detailed statistics are given in Table 5. The mean relative error of all sites drops from 1.18 % of the PWV_BTm to 0.91 % of the PWV_VTm. PWV_VTm has the minimum mean relative errors at 51.35 % of the sites, while PWV_STm is superior at 27.03 % of the sites. PWV_STm and PWV_VTm obtain relative RMSEs smaller than 1.0 % at 55 sites, while only 28 sites of PWV_BTm, 31 sites of PWV_LTm and 22 sites of PWV_GTm perform similarly. For example, at the ALIC site (23.67^∘ S 133.89^∘ E, in Australia), with a mean PWV of approximately 23 mm, the relative RMSE dropped from 1.97 % of PWV_BTm to 1.10 % of PWV_VTm. The time series of the relative differences of PWV_BTm, PWV_LTm, PWV_STm, PWV_VTm and PWV_GTm are given in Fig. 12. We found that some relative RMSEs could be reduced by more than 2 % from PWV_BTm to PWV_VTm. Obviously, PWV_BTm and PWV_LTm have larger relative errors throughout the year, while the PWV differences are significantly larger only in the summer season (when the PWV values are highest). Apparently, the T_m variations in summer are not modeled well by either the Bevis model and the latitude-related model. PWV_STm eliminate those large differences but still retain some residual errors, which are removed by more than 0.5 mm in PWV_VTm. PWV_GTm has some large errors during the period from May to July. All of these results demonstrate that our time-varying model has a precision advantage.

Figure 13RMSEs of the PWV_BTm, PWV_STm, PWV_VTm, PWV_GTm and PWV_ETm compared with the PWV_RS at 11 IGS stations in 2016. PWV unit is millimeters.

Download

Figure 14PWV differences of the PWV_BTm, PWV_LTm, PWV_STm, PWV_VTm, PWV_GTm and PWV_ETm compared with the PWV_RS at the NIRL station in the year 2016. PWV unit is millimeters.

Download

5.3 Comparisons between GPS–PWV and radiosonde PWV

Among our selected 74 IGS sites, there are only 11 sites located within 5 km to a nearby IGRA radiosonde station. At these common stations, we generated PWV from the radiosonde data (PWV_RS) by adjusting the sounding profiles to the heights of IGS sites. It is worth noting that a geoid undulation correction should be carried out on each IGS site geoid height (Jiang et al., 2016). Then, we compared PWV_BTm, PWV_LTm, PWV_STm, PWV_VTm, PWV_GTm and PWV_ETm with PWV_RS. Figure 13 shows the statistics. The RMSEs of GPS–PWV are approximately 1–5 mm. Comparisons indicate that the RMSEs of different GPS–PWV retrievals are very close (differences < 0.2 mm) regardless of the applied T_m sources at most of the selected sites. This means that other errors (e.g., ZTD estimation errors or sounding sensors errors) instead of the T_m make up the bulk of the differences between the GPS–PWV and the radiosonde PWV. Actually, each sounding does not represent the vertical sounding centered at the radiosonde site because of the complex path of the balloon. And GPS–PWV represents the averaged value of the water vapor zenithal projection from all the slant signal paths during the observation period. Such differences can introduce significant uncertainty into our comparisons. However, we still found obvious gaps between PWV at the NRIL station (88.36^∘ N 69.36^∘ E, 4.1 km away from the IGRA station no. 23078 in Russia). The RMSE decreases from 2.29 mm of PWV_BTm to 1.84 mm of PWV_VTm and 1.42 mm of PWV_ETm. As shown in Fig. 14, the large PWV differences appear mainly from May to September. During those five months, the mean GPS–PWV difference to PWV_RS decreases by over 30 % from 2.52 mm of PWV_BTm to 1.67 mm of PWV_VTm, and the reductions of GPS–PWV error are mainly around 1–2 mm. This is attributed to the wetter atmosphere in these months. As indicated by the uncertainty analysis in Sect. 5.1, the improvement in the accuracy of T_m can be translated into more error reduction in the GPS–PWV retrieval with higher values of PWV.