Determination of zenith hydrostatic delay and its impact on GNSS-derived integrated water vapor

Surface pressure is a necessary meteorological variable for the accurate determination of integrated water vapor (IWV) using Global Navigation Satellite System (GNSS). The lack of pressure observations is a big issue for the conversion of historical GNSS observations, which is a relatively new area of GNSS applications in climatology. Hence the use of the surface pressure derived from either a blind model (e.g., Global Pressure and Temperature 2 wet, GPT2w) or a global atmospheric reanalysis (e.g., ERAInterim) becomes an important alternative solution. In this study, pressure derived from these two methods is compared against the pressure observed at 108 global GNSS stations at four epochs (00:00, 06:00, 12:00 and 18:00 UTC) each day for the period 2000–2013. Results show that a good accuracy is achieved from the GPT2w-derived pressure in the latitude band between −30 and 30 and the average value of 6 h root-mean-square errors (RMSEs) across all the stations in this region is 2.5 hPa. Correspondingly, an error of 5.8 mm and 0.9 kg m−2 in its resultant zenith hydrostatic delay (ZHD) and IWV is expected. However, for the stations located in the mid-latitude bands between −30 and −60 and between 30 and 60, the mean value of the RMSEs is 7.3 hPa, and for the stations located in the high-latitude bands from −60 to −90 and from 60 to 90, the mean value of the RMSEs is 9.9 hPa. The mean of the RMSEs of the ERAInterim-derived pressure across at the selected 100 stations is 0.9 hPa, which will lead to an equivalent error of 2.1 mm and 0.3 kg m−2 in the ZHD and IWV, respectively, determined from this ERA-Interim-derived pressure. Results also show that the monthly IWV determined using pressure from ERAInterim has a good accuracy − with a relative error of better than 3 % on a global scale; thus, the monthly IWV resulting from ERA-Interim-derived pressure has the potential to be used for climate studies, whilst the monthly IWV resulting from GPT2w-derived pressure has a relative error of 6.7 % in the mid-latitude regions and even reaches 20.8 % in the highlatitude regions. The comparison between GPT2w and seasonal models of pressure–ZHD derived from ERA-Interim and pressure observations indicates that GPT2w captures the seasonal variations in pressure–ZHD very well.


Introduction
Water vapor as a principal atmospheric parameter is a central component in both Earth's energy budget and water cycle.Accurate knowledge of water vapor is not only vital for weather forecasting but also an important independent data source for global climate studies.For the last decade, the Global Navigation Satellite System (GNSS) has been used as an emerging and robust tool for remotely sensing integrated water vapor (IWV) for the monitoring of the real-time IWV variations in the atmosphere (Schneider et al., 2010;Rohm et al., 2014;Zhang et al., 2015;Li et al., 2014Li et al., , 2015;;Guerova et al., 2016) or the studies of climate (Nilsson and Elgered, 2008;Jin and Luo, 2009;Vonder Haar et al., 2012;Ning Published by Copernicus Publications on behalf of the European Geosciences Union.and Elgered, 2012) due to its 24 h availability, high accuracy, global coverage, high resolution and low cost.The atmospheric parameter directly estimated from GNSS measurements is the GNSS signal's tropospheric zenith total delay (ZTD) which can be effectively divided into the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD).The ZHD can be accurately determined using the surface pressure observed by meteorological sensors.The GNSS-derived IWV over a station can then be obtained by multiplying the ZWD with a conversion factor which is a function of the water-vapor-weighted mean temperature T m over the station.The T m can be determined using one of the following three methods: (1) temperature and humidity profiles from either radiosonde observations or atmospheric reanalysis datasets (Wang et al., 2005(Wang et al., , 2016a)); (2) the relationship between surface temperature T s and the water-vapor-weighted mean temperature T m (Bevis et al., 1992(Bevis et al., , 1994;;Ross and Rosenfeld, 1997); and (3) a blind model developed from atmospheric reanalysis products (Yao et al., 2012(Yao et al., , 2015;;Böhm et al., 2015).
Motivated by our early research (Wang et al., 2016a), it is vital to assess the performance of different methods for determining the ZHD on a global scale, which is essential in the development of a reliable global long-term IWV time series for climate studies.Although the ZHD can be accurately obtained from surface pressure observations, few GNSS stations were installed with meteorological sensors back in the 1990s since these stations were established mainly for precise positioning and navigation applications.Therefore, the lack of meteorological data (i.e., pressure) at these stations is a serious issue for the use of these historical GNSS data for global climate studies.To address this issue, an alternative method is to use pressure derived from a global atmospheric reanalysis (e.g., European Centre for Medium-Range Weather Forecasts Reanalysis, ERA-Interim; Dee et al., 2011) or a blind model (e.g., Global Pressure and Temperature 2 wet, GPT2w; Böhm et al., 2015).In this study, the errors in the pressure derived from these two approaches are investigated and the impact of these errors on the subsequent ZHD and IWV determination has also been studied.
For the performance assessment of these two methods, the ERA-Interim-derived pressure (using two computation methods) and the GPT2w-derived pressure are compared against the surface pressure measured at the 108 global GNSS stations during the period 2000-2013.To determine the pressure at a specific station from the ERA-Interim, first it is necessary to vertically extrapolate the pressure at surrounding grids on the nearest pressure level to the height of the station.Then, the pressure at this station is either horizontally interpolated from the obtained four surrounding grids (fourpoint method) or set to the pressure at the nearest grid at the same height (one-point method, no horizontal interpolation).In this study, the impact of the error in these pressure values on the ZHD is evaluated by comparing the resultant ZHD against the ZHD derived from the surface pressure measured.Similarly, the impact of these errors on the IWV is also as-sessed.As stated in the product requirement document from the EIG EUMETNET GNSS Water Vapour Programme (Offiler, 2010), the accuracy of IWV for the Global Climate Observing System (GCOS) required is better than 3 kg m −2 .The "breakthrough" and "goal" accuracy are 1.5 and 1 kg m −2 , respectively.Since a 1 hPa error in pressure will lead to a 2.3 mm error in its resultant ZHD and about 0.35 kg m −2 in its resultant IWV, the errors of pressure and ZHD should be under 2.9 hPa and 6.6 mm, respectively, to ensure a better than 1 kg m −2 accuracy for their resultant IWV.Bock et al. (2013) also pointed out that for climate monitoring it is necessary to achieve a 3 % or better accuracy from GNSSderived IWV.Since for climate studies the accuracy of IWV (in terms of root-mean-square error, RMSE) can be obtained by averaging observations of a site over a long period (e.g., a month), the error in IWV is studied not only for the four epochs of each day but also for the monthly mean.
2 Datasets and methods

Datasets
In this study, surface pressure observations provided by the Scripps Orbit and Permanent Array Center (SOPAC), ERA-Interim from ECMWF and the GPT2w model are used to investigate the performance of GPT2w and ERA-Interim in the determination of pressure.To investigate the impact of the error in pressure on the resultant IWV, the IWVs at 108 stations are calculated and analyzed with the ZTD reanalysis products provided by the Center for Orbit Determination in Europe (CODE).

Surface pressure observations and GNSS-derived ZTD
Recently, reprocessing of historical GNSS data has emerged as a new initiative in the world GNSS community, driven by the high demand of various new scientific studies.For example, CODE has reprocessed GNSS observations at 371 global IGS (International GNSS Service) stations for the period of 1994-2013 and made the results public to the science community (see ftp://ftp.unibe.ch/aiub/REPRO_2013/).Some of these centers have also collected meteorological observations at these IGS stations, e.g., the Scripps Orbit and Permanent Array Center (see http://garner.ucsd.edu/pub/met/).In this study, the meteorological observations provided by the SOPAC and the ZTD results provided by CODE are used to determine ZHD and IWV.Surface pressure provided by SOPAC is used to validate the performance of ERA-Interim-and GPT2w-derived surface pressure.However, as suggested by previous studies (Wang et al., 2007;Heise et al., 2009;Bianchi et al., 2016), meteorological data provided by IGS need to be rigorously screened before they are used to calculate IWV.In this study, the pressure values from 131 stations are screened carefully based on the method suggested by Wang et al. (2007) to prevent those poor-quality surface pressure observations from being used.Firstly, the time series of pressure at all 131 stations are visually checked carefully to delete stations at which the pressure values have obvious large noises or offsets (which might be caused by the change of pressure sensors).This leads to 23 stations being deleted, and the remaining 108 stations (see the Supplement) are used in this study.Then, the pressure values at these 108 stations are further checked for detecting and excluding unrealistic values out of the range between 550 and 1100 hPa.The third and also the last step is identifying those pressure values that depart from the mean value of more than three standard deviations at each station, which leads to about 0.5 % of the data being detected and excluded from use.It also needs to be noted that the IGS surface pressure observations are not homogenized, and they are therefore not suitable for long-term climate studies.
Since the time series of both the GNSS-derived ZTDs and the pressure measured contain gaps, the missing-data rates of these time series need to be investigated.As a result, a comparison between the meteorological observations and ZTDs over 371 global stations at four epochs (00:00, 06:00, 12:00, 18:00 UTC) of each day during the whole period of 2000-2013 is carried out.The missing-data rates of pressure observations compared to the ZTDs are shown in Fig. 1, which indicates that at most stations only few meteorological observations were recorded.Statistical results indicate that the percentage of missing data smaller than 25 % happens only at 6 % of the stations and there were no pressure observations recorded at 67 % of the stations.Thus, a reliable alternative method for obtaining surface pressure is needed for the determination of IWV from historical GNSS observations.

ERA-Interim
The global atmospheric reanalysis datasets used are from ERA-Interim, which is produced by the European Centre for Medium-Range Weather Forecasts (ECMWF; Dee et al., 2011) and has been widely used for global climate studies.The ERA-Interim data cover the period from 1 January 1979 onwards, and near-real-time information is also generated.All these are available at 00:00, 06:00, 12:00 and 18:00 UTC of each day.The spatial resolution of the ERA-Interim is approximately 80 km (reduced Gaussian grid N128) in horizontal and at 37 vertical pressure levels.The ERA-Interim is used to determine the conversion factor and compute surface pressure in this study.

GPT2w
The blind model selected is GPT2w (Böhm et al., 2015), which provides the mean value plus annual and semiannual amplitudes of the pressure.Böhm et al. (2007) first developed a blind model, called GPT, based on 3-year monthlymean profiles of pressure and temperature from the ECMWF 40-year reanalysis data (ERA40) for providing pressure and temperature at any site in the vicinity of the Earth's surface.Lagler et al. (2013) then developed GPT2 by improving and combining GPT and Global Mapping Function using data of 10-year (2001-2010) 37-month-mean pressure levels from ERA-Interim.The latest GPT2w (Böhm et al., 2015) is an extension of GPT2 with an improved capability in the determination of zenith wet delays in the blind mode.The GPT2w model with the gridded input file is available at http:// ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/.To determine the pressure value at a GNSS station, the GPT2w model first determines the pressure values at four nearest grid points surrounding the station based on a 1 • × 1 • external grid file.Then, the pressure values at these four grids at the height of the station are calculated using an vertical exponential trend coefficient related to the inverse of the virtual temperature (Böhm and Schuh, 2013).Lastly, the pressure at the GNSS station is linearly interpolated from the above four pressure values X. Wang et al.: Determination of zenith hydrostatic delay 2.2 Methods for interpolating pressure data from ERA-Interim Since one nearest pressure level to a GNSS station is used to compute the pressure at this station, the vertical distance between pressure level and this GNSS station will probably affect the accuracy of the pressure derived from ERA-Interim due to its limited vertical resolution.We calculated the mean value of the height differences between the GNSS station and four surrounding grids at the nearest pressure level and studied the relationship between the errors in ERA-Interimderived pressure and this mean height difference.Results indicate that the difference between the GNSS station and its nearest pressure level is between 26 and 181 m at all 108 stations.Additionally, the error (bias and RMSE) in ERA-Interim-derived pressure does not show an obvious dependence on this height difference.For the height below 10 km, the height between two pressure levels in ERA-Interim is about 500 m, which means we can always find the nearest pressure level to a GNSS station with a height difference less than 250 m.Therefore, the vertical resolution of ERA-Interim seems sufficient for the pressure interpolation for the GNSS station.However, in mountainous areas where the pressure levels in ERA-Interim are below the actual terrain, the pressure values on these levels are extrapolated from the model and the results may not be in a good accordance with pressure observations.When ERA-Interim is used to obtain the pressure at a GNSS station, either the pressure from the nearest grid (one-point method) or four surrounding grids (four-point method) and also on the nearest pressure level is used.

One-point method
The one-point method uses meteorological data at an ERA-Interim grid that is closest to the GNSS station to compute the pressure P s for the station (ICAO, 1993) by where p c , T c and H c are the air pressure, temperature (in K) and height (in m) of the nearest ERA-Interim grid point (on the nearest pressure level), respectively; γ = 0.0065 K m −1 is the standard temperature lapse rate; H s is the station height; M = 0.0289644 kg mol −1 is the molar mass of dry air; R = 8.31432 (N m (mol K) −1 ) is the ideal gas constant and g is a gravitational parameter which can be determined by where ϕ is the latitude of the GNSS station.

Four-point method
The four-point method uses the pressure values from the nearest pressure level at four neighboring grid points to determine the pressure for the GNSS station (ϕ s , λ s , H s ).Two procedures need to be carried out: a vertical computation procedure and a horizontal computation procedure.The vertical computation uses pressure values at the nearest pressure level to determine four pressure values, p i (i = 1, 2, 3, 4), at four nearby grid points at the height of the GNSS station using Eqs.( 1) and (2).After p i (i = 1, 2, 3, 4) is obtained, the next step is to horizontally interpolate the pressure value at the GNSS station (ϕ s , λ s , H s ) by where w i is a weighting coefficient computed by where w * i is a weighting coefficient for the ith grid point calculated by where ψ i is the spherical distance between the grid point and the GNSS station.

Computation of ZHD and IWV
One of the most commonly used methods to obtain the ZHD is using the Saastamoinen formula, which is a function of surface pressure (Saastamoinen, 1972;Elgered et al., 1991;Niell et al., 2001).Davis et al. (1985) pointed out that the uncertainty in the ZHD obtained from the Saatamoinen formula was 0.5 mm, if uncertainties in the physical constants and the calculation of the mean value of gravity were taken into consideration.This magnitude of uncertainty will introduce an uncertainty of less than 0.1 kg m −2 in its subsequent IWV determination, which can be neglected.Equations ( 6) and ( 7), developed by Elgered et al. (1991) based on the Saastamoinen formula, are adopted in this study to obtain the ZHD (in mm).These formulae have been widely used in many studies (Bevis et al., 1992;Bock and Doerflinger, 2001;Bokoye et al., 2003;Kleijer, 2004;Musa et al., 2011;Kumar et al., 2013;Norazmi et al., 2015).
where P s is the total pressure (in hPa) at the station's height (in km) and accounts for the variation in the gravitational acceleration at the station with latitude ϕ and height H above a reference ellipsoid.
In the calculation of the IWV with Eq. ( 8), the ZTD provided by CODE and the conversion factor computed from the temperature and humidity profiles from the ERA-Interim are adopted.
More details about the computation of have been published in Wang et al. (2016a), and here only a brief description of the calculation of is given.As a function of T m , can be calculated using where ρ is the density of liquid water, R v is the specific gas constant for water vapor, and m is the ratio of the molar masses of water vapor and dry air.The values of physical constants k 1 = 77.60 ± 0.05 K mbar −1 , k 2 = 70.4± 2.2 K mbar −1 and k 3 = 3.739 ± 0.012 10 5 K 2 mbar −1 are from the widely used formula for atmospheric refractivity (Bevis et al., 1994), and the constant k 2 derived from Eq. ( 10) was set to 22.1±2.2K mbar −1 as suggested by Bevis et al. (1994).
T m in Eq. ( 9) is the water-vapor-weighted mean temperature (Davis et al., 1985) and is approximated as where P v is the partial pressure (in hPa) of WV, T is the atmospheric temperature (in K) and i is the ith pressure level.P v can be calculated using (Tetens, 1930) P si = 6.11 × 10 where P s is the saturated vapor pressure over water, rh is the relative humidity and T is the atmospheric temperature (in C).

Comparison and analysis
One possible error source for the determination of pressure from GPT2w and ERA-Interim is the representativeness error in this model due to the limited model resolution (Janjić and Cohn, 2006;Bock and Nuret, 2009).The representativeness error arises when the point observations can well represent small spatial scales but the model cannot, and this error may be extreme in complex mountainous terrain, where there is a mismatch between the model and actual terrain (Zhang et al., 2013).Although a previous study (Wilgan et al., 2015) indicated that the one-point method may have a better performance than the four-point method in mountainous areas, this phenomenon was not confirmed in our study.To investigate the performance of the aforementioned methods in the determination of pressure and its impact on the resultant ZHD and IWV, the pressure, ZHD and IWV at 108 stations for the period 2000-2013 resulting from both GPT2w and ERA-Interim using the aforementioned two methods are compared against surface pressure measurements and their resultant ZHD and IWV.

Comparisons of 6 h pressure data
To investigate the performance of the GPT2w and ERA-Interim methods, the biases and RMSEs of the resultant pressures, ZHD and IWV at 108 stations for the period of 2010-2013 (using observed pressure as a reference) are computed.
The pressure derived from the GPT2w model, the ERA-Interim with the one-point method and ERA-Interim with the four-point method are named as P_GPT2w, P_ERA1 and P_ERA4, respectively.Equation ( 14) shows the calculation of RMSE, which is a widely used measure of the differences between model-derived values and reference/true values such as the observed values.
where P _cal is the pressure derived from GPT2w or ERA-Interim, P _ref is the reference pressure from observations and n is the number of pressure observations.As indicated in Fig. 2, the RMSE of P_GPT2w is obviously latitude dependent and also obviously larger than that of both P_ERA1 and P_ERA4 at the same stations.For the 28 stations located in the low-latitude band between −30 and 30 • , the RMSEs of P_GPT2w are in the range of 1.4 to 4.3 hPa and the mean of these 28 RMSEs is 2.5 hPa.For the 59 stations located in the mid-latitude bands of −30 to −60 and 30 to 60 • , the RMSEs of P_GPT2w are in the range of 2.8 to 11.9 hPa with a mean of 7.3 hPa.For the 13 stations located in the high-latitude belts of −60 to −90 and 60 to 90 • , the RMSEs are in the range of 7.8 to 12.4 hPa and the mean value reaches 9.9 hPa.Therefore, GPT2w has a far better performance in the low-latitude regions than in the other regions in the determination of surface pressure, and the error of the GPT2w-derived 6 h pressure data may be dominated by the spatial and temporal representativeness error in GPT2w.
Although, the RMSEs of P_ERA1 and P_ERA4 across the 108 stations are mostly under 2 hPa, there are still 8 stations (JOZ2, WHIT, WROC, OHIG, GOPE, BOR1, SOFI and WUHN) that have an RMSE larger than 2 hPa, and 7 out of these 8 stations have a bias larger than 1 hPa or smaller than −1 hPa.The large difference between the observed surface pressure and the pressure derived from ERA-Interim is probably caused by poor-quality observations at these stations.Therefore, these eight stations are not used in the fol-  lowing sub-sections, i.e., only the remaining 100 stations (Table 1) are used to assess the performance of ERA-Interim and GPT2w in the determination of pressure.
It can be seen from Fig. 2a that there is no obvious correlation between the magnitudes of the biases (between −1 and 1 hPa at most stations) and the latitudes of stations for the pressure derived from GPT2w and ERA-Interim.Figure 2a also shows that the biases of the pressure obtained from all three methods are quite similar at most stations.As shown in Fig. 2b, the results from the ERA-Interim using the two aforementioned computation methods have a similar accuracy.The RMSEs of both P_ERA1 and P_ERA4 are in the range of 0.2 to 4.0 hPa, and the mean values of these two sets of RMSEs are both about 0.7 hPa in the low-latitude region and around 1.1 hPa in the other regions.Figure 3 indicates that the distribution of both bias and RMSEs in all three sets of pressure has no obvious correlation with stations' height.Therefore, the error in the ERA-Interim-derived pressure-  ZHD is not dominated by the spatial representativeness error in ERA-Interim.
As indicated by our results, except for the aforementioned eight stations that probably have poor-quality surface observations, the accuracies of the pressure derived from ERA-Interim at the remaining 100 stations are quite similar and they do not show any dependence on the latitude and height of stations.Therefore, the accuracy and spatial resolution of ERA-Interim seems sufficient for the determination of pressure on a global scale, and the model's representativeness error does not have an obvious impact on the resultant pressure.The similar accuracies of the P_ERA1 and P_ERA4 results indicate that if there are representativeness errors in ERA-Interim they tend to be similar in the one-point and four-point data.It is also likely that the differences between the ERA-Interim-derived pressure and observed pressure are dominated by the error in the observations rather than the error in the ERA-Interim model.However, the differences between the pressure derived from ERA-Interim and GPT2w are mainly caused by the differences of their spatial and temporal representativeness.

Comparison of 6 h and monthly pressure, ZHD and IWV
A 1 hPa error in determined pressure will lead to a 2.3 mm error in its resultant ZHD and about 0.35 kg m −2 in its re-sultant IWV.Therefore, the characteristics of the spatial distribution of the errors in ZHD and IWV are quite similar to that in pressure.Table 1 (more details can be found in Supplement 1) shows the statistical results of the bias and RMSE of the pressure, ZHD and IWV derived from the abovementioned three methods for three latitudinal bands.One can see that the accuracies of the ZHD and IWV derived from ERA-Interim with the two aforementioned methods are quite similar and noticeably better than the accuracy of that derived from GPT2w.As listed in Table 1, the biases of the ERA-Interim-derived ZHD are within −1 to 1 mm and the RMSE is about 2 mm at these 100 stations.In terms of IWV, the RMSE of the ERA-Interim-derived IWV is about 0.3 kg m −2 on a global scale.However, the RMSEs of the GPT2wderived ZHD and IWV reach 22.51 mm and 3.33 kg m −2 , respectively, in the high-latitude region, which are significantly larger than those in the low-latitude region.
Since for climate studies we are more interested in monthly mean data, the accuracy of monthly IWV is analyzed.The accuracy of monthly IWV calculated from GPT2w-and ERA-Interim-derived pressure is measured by the difference from the monthly IWV resulting from surface pressure observations (the latter is the reference).IWVs obtained from surface pressure P_GPT2w, P_ERA1 and P_ERA4 are named as IWV_GPT2w, IWV_ERA1 and IWV_ERA4, respectively.Table 2 (more details can be found in Supplement 2) shows the statistical results of the bias,  RMSE and relative error of the monthly IWV derived from the aforementioned three methods at 99 stations.The relative error is calculated using 100• RMS/IWV.It should be noted that the aforementioned eight stations with possible poor data quality and one station with a large amount of missing data, which therefore cannot be used to calculate the monthly mean result, are not used in Table 2.It can be seen that the error in both IWV_ERA1 and IWV_ERA4 is quite small, with an RMSE of about 0.23 kg m −2 on a global scale.The relative errors in both IWV_ERA1 and IWV_ERA4 in the low-latitude, mid-latitude and high-latitude regions are about 0.8, 1.9 and 2.5 %, respectively.The mean relative er-ror of both IWV_ERA1 and IWV_ERA4 across all 99 stations is about 1.6 %.However, for IWV_GPT2w, the relative error is as high as 6.7 % in the mid-latitude regions and can even be up to 20.8 % in the high-latitude regions.Therefore, the monthly IWV calculated from ERA-Interim-derived pressure has a good accuracy, especially in the low-and midlatitude regions; thus, it has the potential to be used for climate studies.Figure 4a1-a6 shows the scatter plots between the monthly IWV derived from surface pressure (x axis) and IWV_GPT2w, IWV_ERA1 and IWV_ERA4 (y axis) at six stations in three latitude belts.As shown in this figure, both monthly IWV_ERA1 and IWV_ERA4 do not contain obvious biases compared to the IWV derived from surface pressure observations.However, from the IWV_GPT2w result, obvious biases can be found at stations in mid-and highlatitude belts, especially at low IWV values.This is more obvious in the scatter plots for the 6 h IWV data as shown in Fig. 4b1-b6.The large error at low IWV values might be caused by the fact that GPT2w-derived pressure usually has a larger error in the cold season, when the IWV values are expected to be low.More discussion regarding this issue is presented in Sect.3.4.
As shown in Table 2, the accuracies of monthly P_ERA1 and P_ERA4 are not dependent on the horizontal location indicating that the monthly pressure data from ERA-Interim are not dominated by the spatial representativeness error in ERA-Interim.The comparison between the statistical results in Tables 1 and 2 also shows that the error in the GPT2wderived monthly pressure data is noticeably smaller than that in the 6 h pressure data derived from GPT2w.This implies that the error in the GPT2w-derived 6 h pressure is affected by the temporal representativeness error in GPT2w.Although the error in the ERA-Interim-derived monthly pressure data is also smaller than that in the 6 h pressure data, the improvements are not as significant as those of GPT2w.

Seasonal variation in pressure and ZHD
In this section, as shown in Eq. ( 15), seasonal models (annual and semiannual cycles) were fitted using the least squares estimation method from three sets of pressure (i.e., P_GPT2w, P_ERA1, P_ERA4 and surface pressure observations).
where P mean is the mean value (a 1 , b 1 ) and (a 2 , b 2 ) are the annual and semiannual amplitudes, respectively.Figure 5 gives an example of the comparison among the ZHD derived from GPT2w and seasonal models from ERA-Interim and pressure observations at site QAQ1.Then, the temporal characteristics of the errors in the pressure and ZHD time series reconstructed using Eq. ( 15) are studied.As expected, the mean values of pressure and ZHD at a GNSS station are closely related to the station height.However, the annual amplitudes of the pressure and ZHD do not show a very strong dependence on the station height but vary by locations.For example, the amplitudes seem larger over stations in eastern Asia than that in western Europe.The semiannual amplitudes are much smaller than the annual amplitudes and thus are not discussed here.Figure 6 shows the annual amplitudes of the ZHD calculated from GPT2w, ERA-Interim and pressure observations at 76 stations (15 stations in low-latitude, 49 in mid-latitude and 12 in high-latitude belts) with a time span of more than 3 years, and the percentage of missing data is less than 50 % in observed pressure.It shows that the annual amplitudes of these four sets of ZHDs show a very good agreement.
Next we compare the reconstructed seasonal models (including fitted annual and semiannual components) from models and observations.As shown in Fig. 7, the biases of pressure from the P_GPT2w and seasonal models of P_ERA1/P_ERA4 are in a range of −2 to 2 hPa, and the RMSE values are under 2 hPa at most stations.The biases of ZHD from the P_GPT2w and seasonal models of ZHD from ERA-Interim are in a range of −4.6 to 4.6 mm and the RMSE values are under 4.6 mm at most stations.Since the seasonal models estimated from P_GPT2w, P_ERA1 and P_ERA4 can be regarded as consistent in temporal scale, the differences among these models are mainly caused by the differences in their spatial resolutions.Since the accuracies of these three seasonal models are similar and agree well with the seasonal models estimated from the time series of pressure observations, the error in the seasonal models from GPT2w-derived pressure is not dominated by its spatial representativeness error.

Seasonal variations in the errors of the pressure and ZHD
To further study the characteristics of the error in pressure and ZHD derived from GPT2w and ERA-Interim, the amplitudes of the annual cycles in the errors of the P_GPT2w and P_ERA1-P_ERA4 time series were calculated based on the 6 h data at the aforementioned 76 stations, and the mean values of the annual amplitudes are around 0.5 hPa (P_GPT2w) and 0.2 hPa (P_ERA1-P_ERA4).Since the RMS errors shown in Table 1 are significantly larger than the amplitudes of the annual cycles in the errors, the variation of the errors in pressure-ZHD are not dominated by these seasonal cycles.
Figure 8 shows the time series of the errors in P_GPT2w and P_ERA1 and the seasonal cycles (Eq.15) estimated from the time series for station JPLM (34.2 • N, 118.2 • W; in the Northern Hemisphere) and station SUTM (32.8 • S, 20.8 • E; in the Southern Hemisphere).We can see that although the seasonal cycles in the time series of the errors in P_GPT2w are quite small, compared to the magnitude of the errors, the magnitudes of the errors at these two stations show a very obvious dependence on season.This shows that the errors in GPT2w-derived pressure-ZHD at JPLM are noticeably larger in December-February (i.e., winter in the Northern Hemisphere) than that in June-August.In contrast, the errors in the GPT2w-derived pressure-ZHD at SUTM are obviously larger in June-August (i.e., winter in the Southern Hemisphere) than that in December-February.This means that the RMSEs of the GPT2w-derived pressure-ZHD at these two stations are larger in the cold season than in the warm season.The results at all 76 stations show that the RM-SEs of P_GPT2w in the cold season are 1.7 times those in the warm season in the high-latitude belt, 1.8 times in midlatitude and 1.2 times in low-latitude.This also explains why

Conclusions
The lack of meteorological data makes it very difficult to take full advantage of historical GNSS data for climate studies.This research is part of our continuous effort to extend our research presented in Wang et al. (2016a) by investigating the alternative methods for the determination of ZHD and developing of a new global long-term IWV time series, which is critical for an improved understanding of climate change using the state of the art GNSS technology.In this study, the accuracies of the ERA-Interim-derived surface pressure (calculated from the nearest grid point and interpolated from four surrounding grid points) and GPT2wderived 6 h pressure have been investigated by comparing them against the observed pressure for 108 stations during the period 2000-2013.The biases of both ERA-Interimderived and GPT2w-derived pressures are between -1 and 1 hPa at most stations.For those stations located in lowlatitude regions, the RMSEs of GPT2w-derived 6 h pressure values are in the range of 1.0 to 4.3 hPa with a mean of 2.5 hPa.However, for the stations located in the mid-latitude belts, the mean value of the RMSEs is 7.3 hPa, and for the

Figure 1 .
Figure 1.Percentage of the missing surface pressure observations compared to the GNSS-derived ZTD at 371 IGS stations for the period of 2000-2013.Results indicate that no pressure observations are recorded at 67 % of stations, and the percentage of missing data is smaller than 25 % only at 6 % of stations.

Figure 2 .
Figure 2. Distribution of biases (a) and RMSEs (b) of 6 h P_GPT2w, P_ERA1 and P_ERA4 with latitudes across 108 GNSS stations for the period of 2000-2013.IWV is integrated water vapor.

Figure 3 .
Figure 3. Distribution of biases (a) and RMSEs (b) of 6 h P_GPT2w, P_ERA1 and P_ERA4 with heights across 108 GNSS stations for the period of 2000-2013.

Figure 4 .
Figure 4. Scatter plots of the monthly mean of IWV (a1-a6) and 6 h IWV (b1-b6) determined using surface pressure observations (x axis) and IWV determined using pressure from GPT2w and ERA-Interim (y axis) for the period of 2000-2013 (ZTD from CODE is used in the computation of IWV).

Figure 7 .
Figure 7. Bias of P_GPT2w and pressure derived from the seasonal models (annual and semiannual cycles) estimated from P_ERA1 and P_ERA4 time series (seasonal model derived from surface observations was used as a reference).

Figure 8 .
Figure 8.Time series of errors in pressure for stations JPLM (a 34.2 • N, 118.2 • W) and SUTM (b 32.8 • S, 20.8 • E) and the seasonal models estimated from the errors.

Table 2 .
Bias, RMSE and relative error of monthly IWV derived from GPT2w and ERA-Interim at 99 stations for the period 2000-2013.