Development of time-varying global gridded Ts-Tm model for precise GPS-PWV retrieval

Water-vapor-weighted mean temperature, Tm, is the key variable for estimating the mapping factor between GPS zenith wet delay (ZWD) and precipitable water vapor (PWV). For the near real-time GPS-PWV retrieving, estimating Tm from 10 surface air temperature Ts is a widely used method because of its high temporal resolution and a fair degree of accuracy. Based on the estimations of Tm and Ts at each reanalysis grid node of the ERA-Interim data, we analyzed the relationship between Tm and Ts without data smoothing. The analyses demonstrate that the Ts–Tm relationship has significant spatial and temporal variations. Static and time-varying global gridded Ts–Tm models were established and evaluated by comparisons with the radiosonde data at 723 radiosonde stations in the Integrated Global Radiosonde Archive (IGRA). Results show that our global 15 gridded Ts–Tm equations have prominent advantages over the other globally applied models. At over 17% of the stations, errors larger than 5 K exist in the Bevis equation (Bevis et al., 1992) and in the latitude-related linear model (Yao et al., 2014b), while these large errors are removed in our time-varying Ts–Tm models. Multiple statistical tests at the 5 % significance level show that the time-varying global gridded model is superior to the other models at 60.03 % of the radiosonde sites. The second-best model is the 1o × 1o GPT2w model, which is superior at only 12.86 % of the sites. More accurate Tm can reduce the contribution 20 of the uncertainty associated with Tm to the total uncertainty of GPS-PWV, and the reduction augments with the growth of GPS-PWV. Our theoretical analyses with high PWV and small uncertainty in surface pressure indicate that the uncertainty associated with Tm can contribute more than 50 % of the total GPS-PWV uncertainty by using the Bevis equation, and it can decline to less than 25 % by using our time-varying Ts–Tm model. However, the uncertainty associated with surface pressure dominates the error budget of PWV (more than 75 %) when the surface pressure has error larger than 5 hPa. GPS-PWV 25 retrievals using different Tm estimates were compared at 74 International GNSS Service (IGS) stations. At 74.32% of the IGS sites, the relative differences of GPS-PWV are within 1 % by applying the static or the time-varying global gridded Ts–Tm equations, while the Bevis model, the latitude-related model and the GPT2w model perform the same at respectively 37.84 %, 41.89 % and 29.73 % of the sites. Compared with the radiosonde PWV, the error reduction in the GPS-PWV retrieval by using a more accurate Tm parameterization can be around 1~2 mm, which accounts for around 30 % of the total GPS-PWV error. 30

Water vapor is an important trace gas and one of the most variable components in the troposphere. The transport, concentration, and phase transition of water vapor are directly involved in the atmospheric radiation and hydrological cycle. It plays a key role in many climate changes and weather processes (Adler et al., 2016;Mahoney et al., 2016;Song et al., 2016).
However, water vapor has high spatial-temporal variability, and its content is often small within the atmosphere. It is a 35 challenge to measure water vapor content accurately and timely. For decades, several methods have been studied, such as radiosondes and water vapor radiometers, sun photometers, and GPS (Campmany et al., 2010;Ciesielski et al., 2010;Liu et al., 2013;Perez-Ramirez et al., 2014;Li et al., 2016). Compared with the traditional water vapor observations, ground-based GPS water vapor measurement has the advantages of high accuracy, high spatial-temporal resolution, all-weather availability, and low-cost (Haase et al., 2003;Pacione and Vespe, 2008;Lee et al., 2010;Means, 2013;Lu et al., 2015).  water vapor products, mainly including precipitable water vapor (PWV), are widely used in many fields such as real-time vapor monitoring, weather and climate research, and numerical weather prediction (NWP) (Van Baelen and Penide, 2009;Karabatic et al., 2011;Rohm et al., 2014;Adams et al., 2017).
GPS observations require some kinds of meteorological elements to estimate PWV. Zenith hydrostatic delay (ZHD) can be calculated using surface pressure P s by equation (Ning et al., 2013) Then, zenith wet delay (ZWD) is generated by subtracting ZHD from zenith total delay (ZTD). ZTD can be directly estimated from precise GPS data processing. Finally, a conversion factor Q, which is used to map ZWD onto PWV, is 50 determined by the water-vapor-weighted mean temperature T m over a GPS station. The mapping function from ZWD to PWV is expressed as (Bevis et al., 1992): where e and T respectively represent vapor pressure in hPa and temperature in Kelvin, i denotes the ith level and Δz i is the height difference of ith level . Vapor pressure e is calculated using equation e=e s RH; RH is the relative humidity, and the saturation vapor pressure e s can be estimated from the temperature observations using Goff-Gratch formula (Sheng et al., 2013).
There are three main approaches to estimate T m . They have respective advantages and disadvantages when they are applied for different purposes: 65 (1) The integration of vertical temperature and humidity profiles are believed to be the most accurate method. The profile data can be extracted from radio soundings or NWP datasets . However, some inconveniences have to be endured. It usually takes considerable amounts of time to acquire the NWP data, which is normally released with large volumes every 6 hours. This limits the use of NWP data in the near real-time GPS-PWV retrieving. Radiosonde data is another profile data source, but it has low spatial and temporal resolution. At most of the radiosonde sites, sounding balloons are daily cast at 70 00:00 UTC and 12:00 UTC. Furthermore, a large amount of GPS stations are not located close enough to the radio sounding sites. Therefore, such methods are appropriate for climate research or the study of long-term PWV trends, but do not meet the real-time requirements.
(2) Several global empirical models of T m are established based on the analyses of T m time series from NWP datasets or other sources (Yao et al., 2012;Chen et al., 2014;Bohm et al., 2015). T m at any time and any location can be estimated from 75 these models. They are often independent of the current meteorological observations which are required to be observed together with the GPS data. However, some important real variations, which may be dramatic during some extreme weather events, can be lost without the constraints of current real data (Jiang et al., 2016). Therefore, these modeled estimates are not accurate enough for high-precision meteorological applications, such as providing GPS-PWV estimates for weather prediction.
(3) Many studies indicated that the T m parameter has a relationship with some surface meteorological elements, such as 80 surface air temperature or surface air humidity (Bevis et al., 1992;Yao et al., 2014a). These surface meteorological parameters can be measured accurately and rapidly. T m is then estimated using these surface measurements. However, these studies also revealed that the relationships are often weak, except the T s -T m relationship. For example, Bevis et al. (1992) introduced the equation T m =0.72 T s +70.2 [K] after analyzing 8712 radiosonde profiles collected at 13 sites in the U.S. over two years. This equation has been widely used in many other studies. 85 According to Rohm et al. (2014), GPS-ZTD can be estimated very precisely by real-time GPS data processing. This means that T m is one of the key parameters in the near real-time GPS-PWV estimation. On the other hand, method (3) is the most suitable method for estimating T m in near real-time because of its balance between timeliness and accuracy. The T s -T m relationship has spatial-temporal variations. Several regional T s -T m equations were established using the profile data over corresponding fields (Wang et al., 2012). However, a T s -T m model without spatial variation is not good enough for a vast field, 90 e.g. the Indian region (Singh et al., 2014). Aside from this, some vast areas have no specific high-precision T s -T m model, for example over the oceans. In general, significant differences exist between oceanic and terrestrial atmospheric properties, especially near the surface layer and within the boundary layer. The change of T s from land to ocean may be very different from that of T m . Therefore it is necessary to model the T s -T m relationship over oceanic regions, since several ocean-based GPS meteorology experiments demonstrated the potential of such technique to retrieve PWV over the broad ocean (Rocken et al., 95 2005;Kealy et al., 2012). A global gridded T s -T m model has been established by Lan et al. (2016). In this model, the 2.0° × 2.5° T m data from "GGOS Atmosphere" and the 0.75° × 0.75° T s data from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data are both smoothed to the resolution of 4° × 5°. However, the T s -T m relationship is varying in time (Yao et al., 2014a), while the Lan et al. (2016) model is static.
The objective of this study is mainly to (1) develop global gridded T s -T m models without any smoothing of the data, then 100 assess their precision, and (2) study the performances of GPS-PWV retrievals using our T s -T m models. Table 1 lists the main differences between the T s -T m model developed in this study and the other global used T m models. In section 2, the data sources and determining methods of T m are introduced in detail. Then, in section 3 we analyze the T s -T m relationships and their variations on a global scale. Global-gridded T s -T m estimating models in different forms are established and evaluated in section 4. Section 5 assesses the accuracies of different PWV retrievals and section 6 presents conclusions based on our experiments. 105

Data Sources and Methodology
As the definition of T m in equation (5), e i parameter at the middle height of ith level is calculated by vertically exponential interpolation of the water vapor pressure of its two neighbor measurement points. The temperature is estimated by linear interpolation of the two neighbor temperatures. The integral intervals are from the earth surface to the top level of profile data. 110 The height of top level depends on the data sources we employed. The essential profile data, including the temperature, height and relative humidity values through the entire atmospheric column, can be obtained from the radiosondes or NWP datasets.
We employed radiosonde data from the Integrated Global Radiosonde Archive (IGRA, ftp://ftp.ncdc.noaa.gov/pub/data/igra) to calculate T m . Version 2.0 of the IGRA-derived sounding parameters provides pressure, geopotential height, temperature, saturation vapor pressure, and relative humidity observations at the observed levels. Bias 115 may be introduced if the integrals were terminated at lower levels (Wang et al., 2005), thus the integrations were performed up to the topmost valid radiosonde data. According to our quality control processes, some radiosonde profile data were rejected.
In each profile, the surface observations must be available and the top profile level should not be lower than 300 hPa standard level. Furthermore, the level number between the surface and the top level should be greater than 10 to avoid too sparse vertical profiles. At most of the radio sounding stations, sounding balloons are launched every 12 hours, and their ascending paths are 120 assumed to be vertical.
Profile data are usually provided by NWP products at certain vertical levels. The ERA-Interim product from ECMWF provides data on a regular 512 longitude by 256 latitude N128 Gaussian grid after the grid transforming performed by the NCAR Data Support Section (DSS). On each grid node of ERA-Interim, temperature, relative humidity and geopotential at 37 isobaric levels from 1000 hPa to 1 hPa can be obtained. Dividing the geopotential by constant gravitational acceleration value 125 (g≈9.80655 m/s 2 ), we can determine the geopotential heights of the surface and levels. Datasets are available at 00:00, 06:00, 12:00 and 18:00 UTC every day and have been covering a period from 1979.01 to present.
In theory, the computation of equation (5) should be integrated through the entire atmospheric column, and the geopotential height should be converted to the geometric height. However, water vapor is solely concentrated in the troposphere, and most of it is specifically located within the first 3 kilometers above sea-level. Moreover, in the two selected 130 datasets, the geopotential heights of top pressure levels are approximately 30~40 km. Geopotential height is very close to geometric height in such height ranges. According to our computation, the relative difference between them is only between 0.1 %~0.9 %. In fact, the height difference z  can be replaced by the geopotential height difference h  in equation (5), since the division can almost eliminate the difference between the two different height types. The T m value nearly has no change after such height replacement. For the convenience of calculation, we directly employed the geopotential height 135 variable. In this paper, we denoted the T m derived from ERA-Interim as T m_ERAI .
At each reanalysis grid node, the computation of equation (5) always starts from the surface height to the top pressure level. The pressure levels below surface height were rejected. T s is defined as the variable of "temperature at 2 meters above ground", and surface water vapor pressure can be derived from the "2 meter dewpoint temperature" variable in ERA-Interim.
These T s were also used in the regression analyses between T s and T m . 140

Correlation between Ts and Tm
Many studies have indicated the close relationship between T s and T m . However, T m is also found not being closely related to T s in some regions, e.g., in the Indian zone (Raju et al., 2007). Using the T m and T s generated from the global gridded reanalysis data, we are able to study the T s -T m relationship in detail.
We first carried out a linear regression analysis on four years of T s and T m data generated from the radiosonde data and 145 the global gridded ERA-Interim datasets, with data covering the period 2009\01 to 2012\12. The analysis results are shown in figure 1. Although the two datasets have different temporal resolutions (12 hours for the radiosonde data and 6 hours for the ERA-Interim data) and spatial resolutions, both analyses agree well with each other. This is expected because the radiosonde data have been assimilated into the ERA-Interim products. Our analyses also indicate that the T s -T m correlation coefficient is generally related to the latitude. The same conclusion has been drawn in other studies (Yao et al., 2014b). Significant positive 150 correlation coefficients can be found at mid-and high-latitudes and reach a maximum in the polar regions. The correlation coefficients drop dramatically at low latitudes. This is because T m is stable there, showing independency of the other parameters.
To study the variations of T s and T m , we illustrated the denary logarithm values of their standard deviations in figure 2. It is evident that T m varies to a lesser degree than T s at low latitudes. Aside from the latitude-related features, there are obvious differences of the T s -T m correlation coefficients between land and ocean. We even found that negative correlation coefficients 155 over certain oceans, e.g., low-latitude Western Pacific, Bay of Bengal or Arabian Sea (see figure 1). Unreliable regression analysis results may be derived when the T s and T m data both have small variations. In figure 3, scatter plots of T s and T m from ERA-Interim at two locations 0.35° N 180.00° E and 70.53° N 180.00° E are given. As the blue dots show, the T s -T m relationship is weak in the areas near the equator, because the entire variation ranges of T s and T m are both within 10 K. This results in a meaningless linear regression (see the magenta line). The T s -T m correlation coefficient is only -0.0893 there. Other 160 than the large spatial variations, studies have revealed that the T s -T m relationship also has temporal variations (Wang et al., 2005). Therefore, a good T s -T m model should take both the spatial and temporal variations into consideration, and this is the main aim in the following sections.

Development of global-gridded Ts-Tm models
Since the T s -T m relationship has large spatial variations, a global gridded T s -T m model is preferred for precise GPS-PWV 175 estimations. In this section, a static global gridded model and a time-varying global gridded model are established and assessed. Especially, there are abrupt changes in the values of constants a and b from land to ocean at the mid to high latitudes due to 185 the different variation features of T s and T m (see figure 2). At the 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 mid to high latitude lands, and lower values are obtained over the oceans or at the 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 is

Time-varying global-gridded Ts-Tm model
The time variation of 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 hr is the observed hour in UTC time; (m 1 , m 2 ), (n 1 , n 2 ) and (p 1 , p 2 ) are fitting coefficients. These equations can reflect the amplitudes of annual, semiannual and diurnal variations in our T s -T m models. Our new regression model found similar values for the coefficients a and b (of its static term) as for the static model in section 4.1, except for some differences over the oceans. In figure 5, besides these constants a and b, we also illustrate the amplitudes of annual, semiannual, and diurnal terms. We can see that there are large annual variations (amplitude > 5 K) in the 205 vast regions from Tibet to North Africa, and in some places of the Siberia and Chile. Large diurnal variations (amplitude > 3 K) mainly occur over the mid-latitude lands such as Northeast Asia or North America. Semiannual variations, however, are small in most of areas except some high-latitudes (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 figure 5, and we can see that the RMSE dropped significantly in the regions with large annual or diurnal variations.

Assessments of Ts-Tm 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 datasets. As a result, we can regard them as 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), 225 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 comparing them with the integrated T m of radiosondes (denoted as T m_RS ) twice a day.
The model estimations of T m are denoted as T m_Bevis , T m_LatR , T m_static , T m_varying , and T m_GPT2w from respectively 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 230 interpolated the coefficients in the T s -T m equations from the neighboring grids to the radiosonde sites. The interpolation formula is expressed as (Jade and Vijayan, 2008) : C site and C i site represent the coefficients in T s -T m equations at the radiosonde site location and its neighboring grids, respectively.
w i are the interpolation coefficients, which are determined using the equation: 235 where R=6378.17 km is the mean radius of the earth,  is the scale factor which equals one in our study, and i  is the 13 angular distance between the ith grid node and the station's position. i  are computed using following formula ( with latitude φ and longitude θ ):   cos sin sin cos cos cos Considering the fact that the reanalysis grids are definite, and every radiosonde site is in situ; we can compute the interpolation coefficients in equation (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 figure 6. Obviously, in many regions, the Bevis equation has a 245 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 West 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 region, etc.
Detailed statistics of the distributions of the bias and RMSE using different models are shown in figure 7 and table 2. At 250 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 have 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 255 than the other T s -T m models.   To identify the superior T m estimation model at each radiosonde site, we employed the following statistical tests under 270 the assumption of a normal distribution of the estimated T m error: (1) First, Brown-Forsythe tests (Brown and Forsythe, 1974) of 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 is significant variance difference between the T m results. If the test rejects the null hypothesis at a 5 % significance level that the errors of model A and B have the same variance, the model with the smaller sample variance is regarded as the better one. 275 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, 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 two models perform almost the same at this radiosonde site. 280 (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 performances. a) T m_Bevis 8% 26% 32% 25% 9% <-3 K -3~-1 K: Finally, we counted the number of sites at which each T m model respectively 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. 285 In figure 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 (< -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 290 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 capturing the variations well. The time series of T m_GPT2w are smoother and cannot capture the fluctuations of the T m time series, causing a worse accuracy 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 295 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).

GPS-PWV retrieving experiments
GPS-PWV has different error sources with different properties. It is complicated to evaluate the GPS-PWV uncertainty here due to the lack of collaborated additional independent techniques to monitor water vapor at the GPS site.

Theoretical analysis of the GPS-PWV uncertainty 305
A comprehensive research on the uncertainty of GPS-PWV has been carried out by Ning et al. (2016). The uncertainties of the ZTD, ZHD and conversion factor Q have been studied in detail. The total uncertainty of GPS-PWV is:  RMSEs of T m given in figure 6. Then we obtained the Q  in equation (11)  The uncertainty of Q will be propagated to the total uncertainty of GPS-PWV according to equation (10). We obtained the contributions of the different terms in equation (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 formulas: , , ,  Table   4 gives five sets of the typical values which are assigned respectively to the Ps  , Tm  , P s and PWV in equations (10)~(12).  (1) No significant difference exists between the figure 10(a) and 10(b). Because of the small value of  c in equation (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(a) and 4(b), a reduction of  Tm can reduce the Q p significantly. For example in 345 figure 10(a), the Q p accounts for 69.54 % with  Tm = 6 K, and it declines to 38.19 % with  Tm = 3 K.
(3) As figure 10(c) shows, the uncertainty associated with  ZTD accounts for the main part of  PWV when the values of PWV and  Ps are not high. With the typical values in table 4(c), the ZTD p can be up to 74.21 % with  Tm = 3 K. And the p Q , however, can drop from 26.76 % to 9.00 % as the  Tm decreases from 6 K to 3 K. Although the p Q is not large under this situation, a smaller  Tm can still reduce the contribution of  Q to the  PWV .

350
(4) The uncertainty associated with  Ps dominates the error budget of PWV when the  Ps is large. In figure 10(d~e), the Ps p is over 80 % with  Tm < 3 K and  Ps = 5 hPa. In figure 10(d), the Q p increases from 7.55 % to 23.19 % as the  Tm raises from 3 K to 6 K. However, in figure 10(e), the Q p only grows from 1.29 % to 4.61 % with the same variation of  Tm . Theoretical analyses on  PWV were also carried out at two representative stations. At the IGRA station No.42971 (20.25° 355 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 respectively 4.30 K, 3.15 K, 2.41 K, 1.93 K and 1.97 K. The  Tm in equation (11)  On the other hand, with  Ps = 5 hPa, the Ps p accounts for more than 75 % of the  PWV while the Q p 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 Northeast China), the mean PWV 365 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 respectively 5.16 K, 3.94 K, 3.54 K, 2.99 K 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 ZTD p averagely contributes over 73.5 % of the  PWV while the Q p averagely contributes less than 10.5 % assuming Ps  = 0.5 hPa and less than 1.5 % assuming Ps  = 5 hPa. But such discussion only concerns the average values. In It is worth mentioning that the uncertainty of ZHD may be underestimated in some situations. There are two reasons for this. Firstly, the calculation of ZHD assumes that the water vapor is not contributed 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 375 significant errors to the ZHD calculation. Secondly and more importantly, the error of P s in equation (1) can be very large sometime. Small Ps  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 of 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 were no nearby barometer to GPS site. The 380 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 due to the more precise T m estimation will be very limited.

Impact of real Tm 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 385 (ftp://cddis.gsfc.nasa.gov/pub/gps/products/troposphere/zpd). 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 equation (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 390 estimates leading to over one hundred 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-PWV estimations is the T m estimation model; therefore, the impact of other errors is excluded.
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 , 395 PWV VTm and PWV GTm at these selected stations were calculated and are illustrated in figure 11. The 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 RMSE 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 ALIC site (23.67° S 133.89° E, in Australia), with a mean PWV of 400 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 figure 12. We found that some relative RMSEs could reduce 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 both Bevis model and the latitude-related model. 405 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 timevarying model has precision advantages.

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 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 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 425 the water vapor zenithal projection from all the slant signal paths during the observation period. Such differences can introduce significant uncertainty to our comparisons. However, we still found obvious gaps between PWV at 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.84mm of PWV VTm and 1.42 mm of PWV ETm . As shown in figure 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 430 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 section 5.1, the improvement in the accuracy of T m can be translated in more error reduction in the GPS-PWV retrieval with higher value of PWV.

Summary and conclusion 440
We developed two global gridded T s -T m models which are respectively static and time-varying with a spatial resolution of 0.75° × 0.75°. The models are established by analyzing the ERA-Interim reanalysis datasets covering the year 2009~2012, which indicated the significant spatial-temporal variations in T s -T m relationship as well as the radiosondes covering the same period. The annual, semiannual, and diurnal variations in T s -T m relationship are considered in the time-varying model. The time-varying global gridded T s -T m model has a significant global precision advantage over the other global applied models, 445 including the Bevis equation, the latitude-related model and the GPT2w model. Average RMSE of T m reduces by approximately 1 K. At over 90 % of the radiosonde sites, our time-varying model has RMSE smaller than 4 K, while the RMSE larger than 5 K nearly disappear. On the other hand, in the Bevis model or in the latitude-related model, there are more than 17 % of the radiosonde sites having RMSE larger than 5 K. Multiple statistical tests at the 5 % significance level identified the significant superiority of our varying model at more than 60 % of the radiosonde sites. Analyses at the specific stations 450 demonstrate that the errors larger than 5 K in the estimated T m series can be eliminated by our varying T s -T m model.
More precise T m estimation can reduce around 20 % of the uncertainty in the conversion factor Q which maps GPS-ZWD to GPS-PWV, and the reduction can be even more than 50 % at some stations. The contribution of the uncertainty associated with Q to the total GPS-PWV uncertainty also declines by using a more precise T m model. The reduction is related to the value of PWV and the uncertainty of surface pressure. With GPS-PWV higher than 50 mm, the uncertainty associated with Q 455 contributes more than 55 % of the uncertainty of GPS-PWV by using the Bevis equation and less than 25 % by using our varying T s -T m model, assuming the ZTD and the surface pressure are measured accurately respectively with the uncertainties of 4 mm and 0.5 hPa. However, the uncertainty in ZTD or in surface pressure would dominate the error budget of GPS-PWV (> 70 %) if the value of GPS-PWV were small or the uncertainty of surface pressure were large. In these cases, the uncertainty associated with Q only contributes around 10 % of the GPS-PWV uncertainty or even smaller. Taking the GPS-PWV using 460 ERA-Interim T m estimates at 74 IGS sites as the references, we found that the GPS-PWV using our time-varying T s -T m model obtained the minimum mean relative error at 51.35 % of the sites, while the GPS-PWV using the static gridded T s -T m model is superior at only 27.03 % of the sites. The differences between GPS-PWV and radiosonde PWV are approximately 1~5 mm.
And our varying T s -T m model can reduce 30 % (around 1~2 mm) of the error in GPS-PWV retrieval with respect to the Bevis

equation. 465
According to our experiments, we are confident that the time-varying global gridded T s -T m models presented here will help us to retrieve GPS PWV more precisely and to study the precise PWV variations in high temporal resolution. Matlab array file consisting of the global gridded coefficients in our model, as well as codes to interpolate coefficients at any given location, are provided as the supplement of this study.