1 Determination of zenith hydrostatic delays and the development of new global long-term GNSS-derived precipitable water vapor

Surface pressure is a vital meteorological variable for the accurate determination of precipitable water vapor 10 (PWV) using Global Navigation Satellite Systems (GNSS). The lack of pressure observations is a big issue for the study of climate using historical GNSS observations, which is a relatively new area of GNSS applications in climatology. Hence the use of the surface pressure derived from either an empirical model (e.g. Global Pressure and Temperature 2 wet, GPT2w) or a global atmospheric reanalysis (e.g. ERA-Interim) 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 for the period 2000–2013. 15 Results show that a good accuracy is achieved from the GPT2w-derived pressure in the latitude band of –300~300 and the average value of Root-Mean-Square (RMS) errors across all the stations in this region is 2.4 mb. Correspondingly, an error of 5.6 mm and 1.0 mm in its resultant zenith hydrostatic delay (ZHD) and PWV is expected. In addition, GPT2w-derived pressure usually has a larger error in the cold season due to large diurnal ranges, which is not considered in the GPT2w model. The average value of the RMS errors of the ERA-Interim-derived pressure across all the 108 stations is 1.1 mb, 20 which will lead to an equivalent error of 2.5 mm and 0.4 mm in its resultant ZHD and PWV respectively. Our research also indicates that the ERA-Interim-derived pressure has the potential to be used as a useful meteorological data source to obtain high accuracy PWV on a global scale for climate studies and the GPT2w-derived pressure can be potentially used for climatology as well although it may be only suitable for the tropical regions.


Introduction
Water vapor (WV) as a principal atmospheric parameter is a central component in both earth's energy budget and water cycle.Accurate knowledge of WV is not only vital for weather forecasting but also an important independent data source for global climate studies.For the last decade, Global Navigation Satellite Systems (GNSS) have been used as an emerging and Atmos.Meas. Tech. Discuss., doi:10.5194/amt-2016-264, 2016 Manuscript under review for journal Atmos.Meas.Tech.Published: 16 August 2016 c Author(s) 2016.CC-BY 3.0 License.robust technique for remotely sensing precipitable WV (PWV) for the monitoring of the real-time PWV variations in the atmosphere (Schneider et al., 2010;Rohm et al., 2014;Zhang et al., 2015;Li et al., 2014;Li 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 and Elgered, 2012;Alshawaf et al., 2016) due to their 24-hour 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 PWV over a station can then be obtained by multiplying the ZWD with a conversion factor П which is a function of water-vaporweighted mean temperature  % over the station.The  % 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., 2016); 2) the relationship between surface temperature  & and the water-vapor-weighted mean temperature  % (Bevis et al., 1992;Bevis et al., 1994;Ross and Rosenfeld, 1997); and 3) an empirical model developed from atmospheric reanalysis products (Lagler et al., 2013;Böhm et al., 2014;Yao et al., 2012;Yao et al., 2015).
Motivated by our early research (Wang et al., 2016), 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 PWV 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.ERA-Interim) or an empirical model (e.g.Global Pressure and Temperature 2 wet, GPT2w).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 PWV determination has also been studied.The global atmospheric reanalysis datasets used are the ERA-Interim which is produced by European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011) and has been widely used for global climate studies.The empirical model selected is Global Pressure and Temperature 2 wet (GPT2w) (Böhm et al., 2014), which provides the mean value plus annual and semiannual amplitudes of the pressure.
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.Then 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 PWV is also assessed.In the determination of PWV, the ZTD reanalysis products provided by the Center for Orbit Determination in Europe (CODE), and the conversion factor П determined using the humidity and temperature profiles from the ERA-Interim are used.

Datasets
Recently, re-processing 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, the CODE has reprocessed GNSS observations at 371 global IGS (International GNSS services) 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 (SOPAC) (see: http://garner.ucsd.edu/pub/met/).In this study, the meteorological observations provided by the SOPAC and the ZTD results provided by the CODE are used to determine ZHD and PWV.
Since the time series of both the GNSS-derived ZTDs and the pressure measured contain gaps, the data missing 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, 06, 12, 18 UTC) of each day during the whole period of 2000-2013 is carried out.
The missing 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.Figure 2 indicates that the missing rate 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 PWV from historical GNSS observations.
As mentioned in the previous section, the ERA-Interim is also used to determine the conversion factor Π and compute surface pressure.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, 06, 12, 18 UTC of each day.Its spatial resolution is approximately 80 km in horizontal and at 60 vertical levels from the surface up to 0.1 hpa.

Methods for determining pressure from ERA-Interim
More details on the determination of the pressure from GPT2w can be found in (Böhm et al., 2014).The main focus of this section is the details of two computation methods used to obtain the pressure over GNSS stations from the ERA-Interim.

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  & for the station (Bosy et al., 2010;Karabatić et al., 2011;Wilgan et al., 2015) where  is the latitude of the GNSS station.

Four-point method
The four-point method (Schüler, 2001) uses the pressure values from two nearest pressure levels at four neighboring grid points to determine the pressure for the GNSS station ( & ,  & ,  & ).Two procedures need to be carried out: a vertical computation procedure and a horizontal computation procedure.The vertical computation uses the pressure values at these two pressure levels to determine the pressure at the height of the GNSS station at four nearby grid points.Let the horizontal coordinates of the four grids be  U ,  U ,  = 1, 2, 3, 4, the heights of the two nearest pressure levels at the  XY grid be  U Z and  U H , and the pressure values at these two levels be  Z and  H , respectively.The pressure value  U at the  XY grid and the height of the GNSS station (i.e. & ) can then be determined by, where  U Z and  U H are the two pressure values at the height of  & at the  XY grid calculated using  Z and  H , respectively, by, and  U Z and  U H in Eq. ( 3) are two weighting coefficients calculated by, After all the pressures  U ( = 1, 2, 3,4) in Eq. ( 3) are obtained, the next step is to horizontally interpolate the pressure value where  o is a weighting coefficient computed by, where  * U is a weighting coefficient for the  XY grid point calculated by, where  U is the spherical distance between the grid point and the GNSS station and  is the weighting power.As suggested by Schüler (2001),  can be set to a value between 1 and 1.5 for the pressure computation.In this study, comparisons among the pressures obtained from  = 1.0, 1.1 ,1.2, 1.3, 1.4, and 1.5 at all the 108 stations for the period 2000-2013 are conducted and the results indicate that the largest difference is between 1.0 and 1.5.The mean of the RMS values of the differences between these two sets of pressure values across all these 108 stations is about 0.05 millibar (mb), which will introduce a difference of only 0.12 mm in the resultant ZHD.In this study,  is set to 1.0 since the difference from different C values is very small.

Comparison and analysis
The pressure and ZHD at 108 stations for the period 2000-2013 determined from both GPT2w and the ERA-Interim using the aforementioned two computation methods are compared against the surface pressure measurements and their resultant To assess the impact of the errors in the three pressure sets, i.e.P_GPT2w, P_ERA1, and P_ERA4, on their resultant ZHD, the biases and RMS errors of the ZHDs with respect to ZHD values derived from the observed surface pressure are calculated.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 mm in its subsequent PWV determination, which can be neglected.
where  & is the total pressure (in mb) at the station's height (in kilometers) and accounts for the variation in the gravitational acceleration at the station with latitude  and height  above a reference ellipsoid.
For convenience, ZHDs resulting from P_GPT2w, P_ERA1 and P_ERA4 are named as ZHD_GPT2w, ZHD_ERA1, and ZHD_ERA4, respectively.As indicated in Fig. 4, the biases of ZHD_GPT2w (Fig. 4  P_GPT2w are latitude dependent, the RMS errors of its resultant ZHD i.e.ZHD_GPT2w also show a strong dependence on latitudes.For the stations located in the low-latitude belts, the RMS errors of ZHD_GPT2w are within 2.3-9.9 mm with a mean of 5.6 mm.For the stations located in the mid-latitude belts, the RMS errors of ZHD_GPT2w are in the range of 6.6-27.1 mm with a mean value of 16.8 mm.For the stations located in the high-latitude belts, the RMS errors of ZHD_GPT2w are in the range of 17.4-28.7 mm with a mean of 22.3 mm.As shown in Fig. 4 (b2) and Fig. 4 (b3), the RMS errors of ZHD_ERA1 and ZHD_ERA4 of the same station are similar and all in the range of 0.5-9.1 mm and the mean values of these two sets of RMS errors are about 2.5 mm.The percentages of the RMS errors of ZHD_ERA1 and ZHD_ERA4 at all the 108 stations less than 3 mm reach 78 % and 72 %, respectively.As suggested by previous studies (e.g.(Rohm and Bosy, 2007;Wilgan et al., 2015)), surface pressure at mountainous areas usually experiences rapid and frequent changes.The RMS errors of the ZHD_ERA1 and ZHD_ERA4 calculated at station WROC located in the Sudety Mountains (Poland) are 2.2 mm and 2.4 mm respectively.This indicates that the ERA-Interim can be used in the mountainous areas for the determination of ZHD and the two methods used in this study perform in a very similar manner.b2), 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 semi-annual amplitudes are much smaller than the annual amplitudes and thus are not discussed here.Apart from annual and semi-annual variations, obvious diurnal ranges, which are defined as the difference between the maximum and minimum values at the same day, are also found at most stations.As indicated by Fig. 6, the diurnal range of pressure and ZHD are related not only to the location of a station but also to the seasons.For example, the diurnal range in the northern polar region is obviously larger during the months from December to February than those from June to August.However, in the southern Polar Regions, the diurnal range is larger during June to August than that from December to February.It should be noted that, in the northern hemisphere the cold season (i.e.winter) is from December to February, whilst in the southern hemisphere the cold season is from June to August.Therefore, the diurnal range is usually larger in the cold season than that in the warm season.In addition, the differences between the diurnal ranges in cold and warm seasons are very small across the stations in the Tropics.

Comparisons of pressure and ZHD in temporal domain
Since the GPT2w model is built based on the monthly mean pressure level data of ERA-Interim and only takes into consideration the annual and semi-annual variations, it cannot reflect the diurnal variations of both pressure and ZHD.
Therefore, theoretically, the errors in the GPT2w-derived pressure and ZHD should be larger in the cold season than that in the warm season.This phenomenon was reaffirmed in our study.Figs.7 and 8 show the errors of the pressure and ZHD derived from the aforementioned three methods at the station SUTM (32.8⁰ S, 20.8⁰ E) located in the southern hemisphere and the station JPLM (34.2⁰N, 118.2⁰W) in the northern hemisphere, respectively.As shown in Fig. 6, the errors of the P_GPT2w and ZHD_GPT2w at SUTM in the southern hemisphere are obviously larger in June-August (winter in the southern hemisphere) than that in December-February.In contrast, for JPLM (Fig. 8) located in the northern hemisphere the errors of P_GPT2w and ZHD_GPT2w are noticeably larger in December-February (winter in the northern hemisphere) than that in June-August.
Fig. 9 shows the RMS errors of the GPT2w-derived ZHD at 81 GNSS stations with a time span of more than 3 years and the missing rate is smaller than 50 %.It should be noted that for long-term climate study, the length of the PWV data used should be as long as possible (e.g.>10 years).The results indicate that the magnitudes of these errors are usually larger in the cold season than that in the warm season.Similarly to Sect.4.2, the impact of the errors in the above three sets of pressure values on their resultant PWV is also analysed using the biases and RMS errors of the PWV.In the calculation of the PWV with Eq. ( 13), the ZTD provided by CODE and the conversion factor П computed from the temperature and humidity profiles from the ERA-Interim are adopted.

Comparisons of PWV
More details about the computation of П have been published in (Wang et al., 2016) and here is only a brief description of the calculation of П is given.As a function of  % , П can be calculated using, where  is the density of liquid water,  ‹ is the specific gas constant for water vapor, and m is the ratio of the molar masses of WV and dry air.The values of physical constants  Z = 70.60 ± 0.05   0Z ,  H = 70.4± 2.2   0Z , and  m = 3.739 ± 0.0012 10 l  H  0Z are from the widely used formula for atmospheric refractivity (Bevis et al., 1994) and the constant  H ˆ derived from Eq. ( 15) was set to 22.1±2.2  0Z as suggested by Bevis et al. (1994).
% in Eq. ( 14) is water-vapor-weighted mean temperature (Davis et al., 1985) and approximated as, where,  ‹ is the partial pressure (in hPa) of WV,  is the atmospheric temperature (in kelvin), and  is the i th pressure level.
‹ can be calculated using, where  & is the saturated vapor pressure, rh is the relative humidity  is the atmospheric temperature (in Celsius).

Comparison of PWV in spatial domain
The errors in the PWVs are investigated in spatial domain to study their spatial characteristics.PWVs obtained from surface pressure P_GPT2w, P_ERA1, and P_ERA4 are named as PWV_GPT2w, PWV_ERA1, and PWV_ERA4, respectively.As As suggested by The E-GVAP (http://egvap.dmi.dk/), the "breakthrough" accuracy of PWV for climate study should be 1.5 mm; the "goal" accuracy of PWV for climate study should be around 1 mm.As defined by the E-GVAP, the "goal" accuracy is an ideal requirement and no further improvements are needed if this requirement can be met.The "breakthrough" accuracy, if achieved, would result in a significant improvement for the climate study and may be considered as an optimum, from a cost-benefit point of view, when planning or designing observing systems.The statistic results of the RMS errors of PWV_GPT2w show that the errors are larger than 1.5 mm at about 91 % (73 of 80) stations located within the mid/highlatitude belts.Moreover, this is only the error in the PWV caused by the error in the pressure i.e. the errors in the GNSSestimated ZTD and the conversion factor П are not taken into account.Therefore, for the GPT2w, it cannot be used in the the mid to high latitude belts for computing the PWV with a "breakthrough" accuracy for the climate studies.On the contrary, the ERA-Interim can be used as a potential useful source of data for computing the PWV for the climate studies on a global scale since the percentages of the RMS errors of PWV_ERA1 and PWV_ERA4 at all the 108 stations that are less than 0.5 mm reach 78 % and 75 %, respectively.

Conclusions
The lack of meteorological data makes it very difficult to take the full advantage of historical GNSS data for climate studies.
This research is our continuous effort to extend our research presented in (Wang et al., 2016)  The results indicate that ERA-Interim has the potential to be used as a meteorological data source for obtaining high accuracy PWV results on a global scale for climate study and the GPT2w can be potentially used for a similar application but may be only suitable to be used in the tropical regions.
The biases of both ERA-Interim-derived and GPT2w-derived pressures are between -1 and 1 mb at most stations, whilst GPT2w has a far better performance in the tropical regions than both mid and high latitude regions.For those stations located in low latitude regions, the RMS errors of GPT2w-derived pressure are in the range of 1.0~4.3mb with the mean of 2.4 mb.However, for the stations located in the mid to high latitude belts, the mean value of the RMS errors is 7.4 mb, and for the stations located in the high-latitude belts, the mean of the RMS errors is 9.
ZHD.The comparison results are analysed in Sect.4.1 and 4.2.Since the ERA-Interim-derived pressure and ZHD are of very high accuracy on a global scale, the pressure and ZHD over 371 global stations at four epochs (00, 06, 12, 18 UTC) of each day during the whole period of 2000-2013 are derived from the ERA-interim and subsequently used to analyse their temporal variations.The analyses results are shown inSect.4.3.Section 4.4  shows the comparisons among the PWVs computed using the aforementioned four sets of pressure.in spatial domainTo investigate the performance of the GPT2w and ERA-Interim methods, the bias and RMS of the resultant pressures 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 the ERA-Interim with the four-point method are named as P_GPT2w, P_ERA1 and P_ERA4, respectively.As shown in Fig.3, the biases of P_GPT2w (Fig.3(a1)), P_ERA1 (Fig.3(a2)), P_ERA4(Fig. 3 (a3)) are between -1 and 1 mb at most stations.No obvious correlation between the magnitudes of the bias and the latitude is found.However, the RMS error of P_GPT2w (Fig.3 (b1)) is very latitude dependent and also larger than that of both P_ERA1 (Fig.3(b2)) and P_ERA4 (Fig.3(b3)) at the same stations.For the 28 stations located in the lowlatitude band of -30°~30°, the RMS errors of P_GPT2w are in the range of 1.0~4.3mb and the mean of these 28 RMS errors is 2.4 mb.For the 65 stations located in the mid-latitude bands of -30°~-60° and 30°~60°, the RMS errors of P_GPT2w are in the range of 2.8~12.0mb with the mean of 7.4 mb.For the stations located in the high-latitude belts of -60°~-90° and 60°~90°, the RMS errors are within 7.6-12.6mb and the mean value reaches 9.8 mb.Therefore, GPT2w has a far better performance in the tropical regions than the other regions in the determination of surface pressure.As shown inFig.3 (b2)      andFig.3 (b3), the results from the ERA-Interim using the two aforementioned computation methods have a similar accuracy.The RMS errors of both P_ERA1 and P_ERA4 are in the range of 0.2~4.0mb and the mean values of these two sets of RMS errors are both to be 1.1 mb.The percentages of those RMS errors of P_ERA1 and P_ERA4 across the 108 stations less than 1.5 mb reach 84 % and 80 %, respectively.4.2 Comparisons of ZHD in spatial domainAtmos.Meas.Tech.Discuss., doi:10.5194/amt-2016-264,2016 Manuscript under review for journal Atmos.Meas.Tech.Published: 16 August 2016 c Author(s) 2016.CC-BY 3.0 License.
indicated in Fig.10, the biases of PWV_GPT2w (Fig.10 (a1)), PWV_ERA1 (Fig.10 (a2)), and PWV_ERA4 (Fig.10 (a3)) are all in the range of -0.5~0.5 mm at most stations.The RMS errors of PWV_GPT2w (Fig.10 (b1)) are obviously smaller in low-latitude regions than that in high-latitude regions.For all the stations located in the low-latitude belts, the RMS errors of PWV_GPT2w are within 0.5~1.6 mm and the mean value of these RMS errors is 1.0 mm.For all the stations located in the mid-latitude belts, the RMS errors of PWV_GPT2w are in the range of 1.1~4.2mm with the mean of 2.6 mm.For all the stations located in the high-latitude belts, the RMS errors of PWV_GPT2w are in the range of 2.6-4.2mm with the mean of 3.4 mm.The RMS errors of PWV_ERA1 and PWV_ERA4 of the same stations are similar and all in the range of 0.1~1.5 mm with the mean of about 0.4 mm.Atmos.Meas.Tech.Discuss., doi:10.5194/amt-2016-264,2016 Manuscript under review for journal Atmos.Meas.Tech.Published: 16 August 2016 c Author(s) 2016.CC-BY 3.0 License.
through investigating the alternative methods for the determination of ZHD and developing of a new global long-term PWV time series, which is critical for an improved understanding of climate change using the state of the art GNSS technology.Our new PWV dataset across 371 stations for the period of 2000-2013 has been made available at https://doi.org/10.1594/PANGAEA.862525.In this study, the accuracy of the ERA-Interim-derived surface pressure (with two computation methods) and GPT2w-derived pressure for 108 GNSS stations during the period 2000-2013 has been investigated by comparing them against the observed pressure.
8 mb.The RMS errors of the two sets of ERA-Interim-derived pressure (with different computation methods) at all the 108 stations are both in the range of 0.2~4.0mb and the mean value is 1.1 mb.The biases and RMS errors of the ZHDs based on ERA-Interim-derived and GPT2w-derived pressure are calculated by using the ZHD values based on the observed pressure as a reference.The biases of the three sets of ZHD are all in the range Atmos.Meas.Tech.Discuss., doi:10.5194/amt-2016-264,2016 Manuscript under review for journal Atmos.Meas.Tech.Published: 16 August 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 1 .Figure 2 .
Figure 1.Missing rate of the surface pressure observations compared to the GNSS-derived ZTD at 371 IGS stations for the period of 2000-2013

Figure. 7 Figure. 8
Figure.7 Time series of errors in pressure (a) and ZHD (b) for station SUTM located in the southern hemisphere