Retrieving Precipitable Water Vapor Data Using GPS Zenith Delays and Global Reanalysis Data in China

Abstract

GPS has become a very effective tool to remotely sense precipitable water vapor (PWV) information, which is important for weather forecasting and nowcasting. The number of geodetic GNSS stations set up in China has substantially increased over the last few decades. However, GPS PWV derivation requires surface pressure to calculate the precise zenith hydrostatic delay and weighted mean temperature to map the zenith wet delay to precipitable water vapor. GPS stations without collocated meteorological sensors can retrieve water vapor using standard atmosphere parameters, which lead to a decrease in accuracy. In this paper, a method of interpolating NWP reanalysis data to site locations for generating corresponding meteorological elements is explored over China. The NCEP FNL dataset provided by the NCEP (National Centers for Environmental Prediction) and over 600 observed stations from different sources was selected to assess the quality of the results. A one-year experiment was performed in our study. The types of stations selected include meteorological sites, GPS stations, radio sounding stations, and a sun photometer station. Compared with real surface measurements, the accuracy of the interpolated surface pressure and air temperature both meet the requirements of GPS PWV derivation in most areas; however, the interpolated surface air temperature exhibits lower precision than the interpolated surface pressure. At more than 96% of selected stations, PWV differences caused by the differences between the interpolation results and real measurements were less than 1.0 mm. Our study also indicates that relief amplitude exerts great influence on the accuracy of the interpolation approach. Unsatisfactory interpolation results always occurred in areas of strong relief. GPS PWV data generated from interpolated meteorological parameters are consistent with other PWV products (radio soundings, the NWP reanalysis dataset, and sun photometer PWV data). The differences between them were approximately 1~3 mm at most at our selected stations, and GPS data processing is the main factor influencing the agreement of the GPS PWV results with those of other methods.

1. Introduction

Water vapor plays a key role in many weather and climate processes (e.g., see [1,2]). Many instruments can be employed to retrieve the atmospheric water vapor content, such as radiosondes, microwave radiometers, sun photometers, and ground-based GPS devices [3]. GPS meteorology, which was first proposed in 1992 by Bevis [4], has proven to be a very effective method to estimate precipitable water vapor (PWV) (e.g., see [5,6]). Information about atmospheric water vapor provided by ground-based GPS meteorology can be significant in many areas, such as numerical weather forecasting (NWP) (e.g., see [7,8]), global or regional atmospheric research (e.g., see [9,10]), and others (e.g., see [11]). In addition to numerous efforts to improve the accuracy of ground-based GPS meteorology (e.g., see [12,13]), the density of GPS sites has also increased in many countries, leading to a considerable increase in the PWV observation density (e.g., see [14,15]). These developments strongly promote the application of ground-based GPS PWV detection in the meteorological field.

A GPS signal is delayed when it propagates through the neutral atmosphere, and part of this delay is caused by atmospheric water vapor, which is mostly concentrated in the troposphere. To derive GPS PWV, tropospheric zenith total delay (ZTD) is estimated by first processing GPS data; then the zenith hydrostatic delay (ZHD) should be subtracted from ZTD to obtain the zenith wet delay (ZWD); finally, ZWD can be converted to PWV by a ratio value $\Pi$, which is related to the weighted mean temperature T_m. The Saastamoinen model [16] has been widely used to estimate ZHD because of its high accuracy of approximately 1–2 mm, and the surface barometric pressure P_s is one of the essential parameters in this model. T_m is often calculated from its linear relationship with the near-surface air temperature T_s. Therefore, to derive GPS PWV with high accuracy, site-specific surface pressure P_s and T_m are also essential, which are, respectively, used to calculate ZHD and $\Pi$ in combination with the GPS data. However, collecting such meteorological parameters requires collocated meteorological sensors, which are often not present at many geodetic GPS stations. To employ those sites without collocated meteorological sensors in GPS meteorology, different methods were studied in various areas. For example, Jade [17] firstly interpolated the 2.5° × 2.5° National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis data to obtain site-specific meteorological parameters for estimating GPS PWV. An interpolation of the NCEP pressure level data and NCEP surface data was explored. The experiments, conducted over a four-year period in India, indicated that the GPS PWV data derived using the interpolated surface pressure and weighted mean temperature from NCEP data were in strong agreement with those from real meteorological observations. However, bias between GPS PWV and horizontally interpolated NCEP PWV results apparently increased with highly undulated terrain. Vey [18] derived GPS PWV data that were obtained by using 141 global GPS stations over a 10-year period from 1994 to 2004, in which T_m data were obtained from the European Center for Medium-Range Weather Forecasts (ECMWF) and P_s data were interpolated from neighboring World Meteorological Organization (WMO) stations. In Karabatic’s study [19], P_s and T_s data of the nearest meteorological station were extrapolated to GPS stations in Austria, and they found that a PWV accuracy better than ±1 mm could be achieved when the distance between the two sites was within 20 km. In Poland, Bosy obtained specific GPS site P_s and T_s data by interpolating nearby meteorological data to the GPS site locations [20]. However, some systematic bias in the interpolation results appeared at some stations, which was attributed to unknown deficiencies in their interpolation procedure. The North American Regional Reanalysis (NARR) dataset for the year 2009 was used to estimate GPS PWV in California and Nevada by Means and Cayan [21]. The surface pressure and temperature were determined by the elevation of a specific station and the geopotential heights of the standard pressure levels of neighboring grids, and simple two-dimensional interpolation of the reanalysis gridpoint surface temperatures was performed. Means [22] also employed the same technique to obtain GPS PWV of over 500 sites from 2003 to 2009. Then, using these GPS PWV data, he studied the temporal and spatial extent of the North American monsoon in California and Nevada. Luo [12] presented a height-dependent linear ZHD correction model using freely accessible meteorological measurement data near GPS sites, and the mean ZHD bias was approximately 5 mm.

The number of GPS receivers deployed in China has greatly increased in recent years [23]. However, many stations are deployed for geodetic applications, leading to few collocated surface meteorological measurements, except at stations within the GPS meteorology network. In this paper, we employ a three-dimensional interpolation methodology and NCEP Final (NCEP FNL) Analysis data to generate surface pressure and temperature data for GPS stations over China. The goal of our study is to demonstrate that either the present available GPS networks or the historical data from previous GPS geodetic observations have great potential to reconstruct the precipitable water field with higher spatial and temporal resolution, which can be used to study the characteristics of water vapor (e.g., see [24]) or analyze the water vapor movement during severe weather (e.g., see [25]) in more detail.

In the following section, we describe the multi-source data and their processing strategy, including the three-dimensional interpolation algorithm applied in our research. Error analyses for the interpolated surface meteorological elements are presented in Section 3. Section 4 provides the results of the interpolation products, and the accuracy of GPS PWV from interpolation products is analyzed by comparison with the PWV results from other approaches. Conclusions are given in Section 5.

2. Methodology

2.1. GPS PWV Estimation

Basically, two solution approaches can be selected for GPS data processing: double-difference (DD) solution mode and precise point positioning (PPP) mode. For PPP, ZTDs are directly estimated together with other unknown parameters, such as the station coordinates or receiver clock-offsets. However, it requires a high-precision satellite clock and ephemeris product [26]. For the DD model, common errors between different receivers and satellites, e.g., the clock-offsets, can be eliminated or greatly reduced. However, to reduce the correlation of observations, it is necessary to introduce three or more distant stations into the network to form long baselines for estimating the absolute ZTD [27]. In fact, each ZTD can be regarded as an average of projection values in the zenith direction of the neutral atmospheric delays along the available slant signal paths from GPS satellites to the receiver during a short time period, which is dependent on the time interval for ZTD estimation. As previously mentioned, we applied the Saastamonien model to separate ZHD from ZTD. This model is expressed as a function of latitude $\theta$ in radians and height above the geoid H in km:

$$\text{ZHD}=[(0.0022768\pm 0.0000005){\text{m hPa}}^{-1}]\times \frac{P_0}{f(\theta ,H)}$$

(1)

where P₀ represents surface pressure in hPa and

$$f(\theta ,H)=1-0.00266\cdot \cos 2\theta -0.00028\cdot H$$

(2)

After the ZHD is calculated, ZWD can be generated by following equation:

$$\mathrm{Z}\mathrm{W}\mathrm{D}=\mathrm{Z}\mathrm{T}\mathrm{D}-\mathrm{Z}\mathrm{H}\mathrm{D}$$

(3)

Then, the ratio value Π can convert ZWD to PWV:

$$\mathrm{P}\mathrm{W}\mathrm{V}=\Pi \times \mathrm{Z}\mathrm{W}\mathrm{D}$$

(4)

The relational expression between Π and T_m (see [28]) is:

$$\Pi =\frac{{10}^6}{{\rho }_{\mathrm{w}}{R}_v(k_2^{\prime }+(k_3/T_m))}$$

(5)

where ${\rho }_w$ is the density of liquid water, R_v is the specific gas constant for water vapor, $k_2^{\prime }$ = (17 ± 10) K·mbar⁻¹, and $k_3$ = (3.776 ± 0.014) × 105 K²·mbar⁻¹. T_m is defined as:

$$T_m=\frac{{\displaystyle {\int }_{z_0}^{+\infty }(e/T)dz}}{{\displaystyle {\int }_{z_0}^{+\infty }(e/T^2)dz}}$$

(6)

where e and T represent vapor pressure and absolute temperature along the zenith path, respectively, and the integral intervals are from the earth surface to the atmospheric top. Simply, T_m can be estimated using the surface temperature measurement T_s:

$$T_m=a\times T_s+b$$

(7)

The frequently used values for coefficients a and b are, respectively, 0.72 and 70.2, proposed by Bevis [4]. However, their values can vary in different areas [29].

Numerous studies regarding ground-based GPS PWV estimations in different regions or under different weather conditions have indicated that the accuracy of the above derivation method reached 1~2 mm compared to other high-precision measuring techniques (e.g., see [30,31,32]). In this study, we used GPS data from the Crustal Movement Observation Network of China (CMONOC) and GAMIT/GLOBK software for GPS data processing. The CMONOC network was built to monitor crustal deformation, water vapor distribution, ionospheric variability, and other geophysical parameters. It consists of GNSS, VLBI, SLR, gravity, and level measurement sites. High spatial-temporal resolution observations are provided by 260 GNSS continuously operating stations, and the distribution of these reference stations is shown in Figure 1. Each station is equipped with a high-precision Trimble NetR8/R9 GNSS receiver, a meteorological sensor, and other auxiliary instruments. With this equipment, GPS data and surface meteorological parameters are continuously collected. GAMIT/GLOBK software was employed as a comprehensive suite of programs developed by MIT to process GPS data. The double-difference model is used in GAMIT software, so at least three long baselines should be included during the GPS data processing procedure. However, no external GPS data are required because of the large areas that CMONOC covers in our study.

2.2. Interpolation of NCEP FNL Global Analysis Data

Interpolating surface meteorological parameters P_s and T_s from the NWP data can be an effective complement for ground-based GPS meteorology as a solution for the absence of collocated meteorological sensors. In this study, we used the NCEP FNL Operational Global Analysis historical data to perform these interpolations. The NCEP FNL Operational Global Analysis data are produced from the same data assimilation and forecast system as the NCEP Global Forecast System (GFS) data. The difference between these two datasets is that approximately 10% more observations are assimilated into the initial condition for FNL than that for GFS. However, because it takes some time to wait for more observational data to be collected, the FNL analysis data are delayed by approximately 60–90 min compared with those from GFS analysis. These data can be downloaded freely from the website http://rda.ucar.edu/datasets/ds083.2/.

The NCEP FNL Analysis data are published every six hours and consist of 1° × 1° horizontal grids. The analyses provide various parameters on the surface, at 26 mandatory levels from 1000 millibars to 10 millibars, as well as at a few other levels. However, only several parameters are used in the interpolation. First, the four neighboring grid nodes are selected according to the latitude and longitude of specific stations. Assuming that the grid node number is i, we extracted pressure values and geopotential heights of the two mandatory levels nearest to the specific station’s height. Then, the scale height $H_p^i$ between these two levels is calculated as:

$$H_p^i=\frac{h_{upper}-h_{lower}}{\ln (p_{lower}/p_{upper})}$$

(8)

where h_lower and h_upper represent the upper and lower level geopotential height and p_lower and p_upper are the air pressure of the lower and upper level, respectively. Pressure $p_z^i$ at a GPS site’s height z is estimated by the following equation:

$$p_z^i=p^i\exp \left(-\frac{z-h^i}{H_p}\right)$$

(9)

where pⁱ is the pressure of height hⁱ and hⁱ is the height of one of the chosen mandatory levels or that of the surface, depending on which is the nearest level to the GPS station’s elevation.

After the vertical interpolations at the neighboring grid nodes, we horizontally interpolated the four pressure values to each GPS site’s position. The interpolation formula is expressed as:

$$p_s={\displaystyle \sum _{i=1}^n{w}^i{p}_z^i}$$

(10)

where p_s is the station’s surface pressure and wⁱ represents the interpolation coefficients of the vertical interpolated pressure $p_z^i$ at point i. wⁱ is determined by the following formula (see [17]):

$$w^i=\frac{{\left(R{\psi }^i\right)}^{-C}}{{\displaystyle \sum _{j=1}^n{\left(R{\psi }^j\right)}^{-C}}}$$

(11)

where R = 6378.17 km is the mean radius of the earth, C is the scale factor, which equals one in this study, and ${\psi }^i$ is the angular distance between the grid node i and the GPS receiver’s position (φ,θ):

$$\cos {\psi }_i=\sin {\phi }_i\cdot \sin \phi +\cos {\phi }_i\cdot \cos \phi \cdot \cos \left({\theta }_i-\theta \right)$$

(12)

In the above equation, φ and θ represent latitude and longitude, respectively. Because the reanalysis grids are definite, the interpolation coefficients in Equation (10) for the stationary GPS site are also fixed and can be stored as constants to avoid duplicate calculation.

Surface temperature generation also starts with vertical interpolation at neighboring points. The adiabatic lapse rate between two levels can be calculated simply from their temperatures and geopotential heights. We first estimated the mean adiabatic lapse rate between the 1000 hPa and 500 hPa pressure levels, and then, temperature $T_z^i$ at the height of the GPS station was given by:

$$T_z^i=T_s^i+\beta (h_s^i-z)$$

(13)

where $\beta$ is the computed adiabatic lapse rate, $T_s^i$ is the temperature at 2 m above the surface, and $h_s^i$ is the corresponding height. Finally, horizontal interpolation was employed to generate the temperature of the GPS station’s point. Equations (10)–(12) are used, while the pressure parameters in Equation (10) are replaced by temperature.

It is worth mentioning that the difference between the altitude height and geopotential height is neglected in our study (see [33]) because the top height we interpolated was below 30 km. However, only the geodetic height can be estimated from GPS data processing. Therefore, we used the EGM2008 model, which was developed by the U.S. National Geospatial-Intelligence Agency (NGA) to convert one GPS station’s geodetic height to the altitude height. Details of the EGM2008 model can be found at http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html. The accuracy of the EGM2008 model with 1° × 1° resolution can reach up to the sub-meter level over China [34]. The change of air pressure produced by several sub-meters of height difference can be ignored; we believe this sub-meter-level accuracy is acceptable for our purposes.

3. Interpolation Results

To assess the quality of interpolated surface pressure and temperature for GPS PWV derivation, we introduced two observation datasets for comparison: CMONOC meteorological measurements and the Integrated Surface Database (ISD), which are provided freely by the U.S. National Climatic Data Center (NCDC).

Specific altitude height is a key parameter during our interpolation procedure. From GAMIT/GLOBK processing, the geodetic height of each CMONOC GNSS station can be estimated. Then, conversion from the geodetic height to the altitude height has to be performed, and NCEP FNL data are interpolated to one site’s altitude height instead of at the geodetic height. As previously mentioned, geoid undulations of CMONOC GNSS sites estimated from the EGM2008 model were used to approximately adjust the geodetic heights to altitude heights. Their values are illustrated in Figure 2. It can be seen that they vary from approximately 20 m in eastern China to approximately −60 m in northwest China, which indicates that an apparent bias in the interpolation of air pressure would exist if no conversion between different height systems was performed.

ISD provides global hourly and synoptic observations of various parameters, such as temperature, dew point, sea level pressure, station pressure, and other observed elements (see [35]). Although ISD data originate from numerous sources, a single common ASCII format and common data model are still developed by ISD to facilitate the use of data. The input data sources compiled by ISD have to be previously processed through quality control. Currently, high spatial coverage is obtained by ISD stations in most parts of China, except some areas in western China. The number of active stations in China included in the ISD database is over 350, and their distribution in 2012 is shown in Figure 3. For some ISD sites, surface pressures have to be generated from sea level pressure.

3.1. Comparison between Interpolated and Observed P_s

Taking the first derivative of Equation (1) with respect to the surface pressure P₀, we can obtain the following formula:

$$\frac{\partial ZHD}{\partial P_0}=[(2.2779\pm \mathrm{0.0.0024}){\text{m hPa}}^{-1}]\times \frac{1}{f(\theta ,H)}$$

(14)

For a stable reference site, $\theta$ and H are constants, so the relationship between the errors of ZHD and P₀ is linear. From Equation (14), it can be inferred that the surface pressure error of 2.8 hPa produces an approximately 6.7 mm error in ZHD, and the equivalent error will be transferred to ZWD through the computation of Equation (3). This equates to a 1-mm error in PWV. The comparison between interpolated and measured surface pressure is based on the above analysis.

Considering that the CMONOC data and ISD data are derived from different sources and that most of the ISD observed data were likely previously assimilated into the NCEP FNL data, we compared the interpolated P_s with their corresponding measured values from two datasets. In our study, we make these comparisons using the full-year data of 2012 from both datasets. BIAS and RMSE are selected as the statistics, and their mathematical expressions are as follows:

$$\{\begin{array}{l}BIAS=\frac{{\displaystyle {\sum }_{i=1}^N\left(IN{T}_i-OB{S}_i\right)}}{N_{data}}\\ RMSE=\sqrt{\frac{{{\displaystyle {\sum }_{i=1}^N\left(IN{T}_i-OB{S}_i\right)}}^2}{N_{data}}}\end{array}$$

(15)

where INT is the interpolated value and OBS represents the observed value for the same element.

We calculated BIAS and RMSE between the interpolated and observed P_s at all of the selected CMONOC and ISD sites.

Table 1. Number of sites at which absolute BIAS or RMSE of interpolated surface pressure were within a given special value domain.

Value Domain	Absolute BIAS as Statistics		RMSE as Statistics
	ISD	CMONOC	ISD	CMONOC
≤1 hPa	356	200	318	158
1 ~ 2.8 hPa	17	36	55	76
2.8 ~ 5 hPa	2	3	2	5
≥5 hPa	1	1	1	1

gives the number of ISD and CMONOC sites at which absolute BIAS or RMSE were within the given value domain. This clearly indicates that both statistical values were below 1 hPa, that is, there was less than 0.5-mm error in the GPS PWV derivation, at most sites. Only at a few sites was the absolute BIAS or RMSE more than 2.8 hPa, which equates to more than a 1 mm GPS PWV error, as previously derived. Furthermore, it can be seen from Figure 4, which illustrates the BIAS and RMSE of all of the selected stations together with their locations on the map, that stations with a larger BIAS or RMSE are mainly distributed in southwest China or the northern part of Xinjiang Province. These sites had a common property in that they are all located in regions of higher relief. Taking the JIULONG station in Sichuan Province (ISD number: 564620) as an example, we found that it had a great difference in height with the four neighboring grids of the NCEP FNL model. The height of this station was 2994 m, whereas the lowest neighboring node’s height was 3703 m, so the altitude height differences were higher than 651 m, reflecting undulating terrain. Its interpolated and measured P_s time series is illustrated in Figure 5. There is an obvious system bias between the two series. The reason for such a system deviation may be that the four pressure values used for horizontal interpolation were all one-side vertical extrapolated instead of interpolated because the JIULONG station was lower than all of the neighboring nodes, and the accuracy would decrease substantially if the vertical extrapolation distances were too long. The other stations that exhibit large interpolation errors also faced the same problem, including the XJWQ and QHME stations in CMONOC. However, more than 97.5% of stations with GPS PWV derivation errors of interpolated P_s were below 1 mm. This demonstrates the effectiveness of our interpolation method in areas without significant relief.

Even under extreme weather conditions, the accuracy of the interpolated P_s did not decrease dramatically according to our experiment. Taking the Haikui typhoon, which landed in Ningbo City, Zhejiang Province on 8 August 2012, as an example, we compared the measured and interpolated P_s of the ZJZS station, which is the nearest CMONOC station to the typhoon landing position. As Figure 6 shows, the differences between the two P_s time series remain very small during the typhoon period from 5 to 11 August 2012. Their mean difference was −0.89 hPa, with a standard derivation of 0.37 hPa, and the largest difference was only −1.75 hPa. We also use the P_s series generated from the GPT model (see [36]) for comparison purposes; its difference from real observations could be as high as 20.51 hPa, which obviously cannot be accepted in GPS PWV retrieval.

3.2. Comparison between Interpolated and Observed T_s

The same statistics were employed to assess the quality of interpolated surface temperature T_s. We also take the first derivative of Equation (5) with respect to the weighted mean temperature T_m:

$$\frac{d\Pi }{d{T}_m}=\frac{{10}^6}{{\rho }_w{R}_v}\cdot \frac{k_3}{{(k_2^{\prime }+k_3/T_m)}^2{T}_m^2}$$

(16)

Figure 7 visualizes the relationship between Π and T_m. Π increases from 0.1363 to 0.1699 with the increase of T_m from 240 K to 300 K, while the first derivative of Π with respect to T_m only decreases from 5.6198 × 10⁻⁴ K⁻¹ to 5.5907 × 10⁻⁴ K⁻¹. The growth rate of Π is 20%, while the decline rate of Π’s derivative is merely 0.52%, which suggests that the derivative of Π with respect to T_m is not sensitive to changes in T_m and remains at a low value. When T_m equals 290 K and Π is 0.1643, even if T_m changed more than 5 K, the value of Π varied by only approximately 2.8 × 10⁻³, which represents a nearly 1.8% relative difference. The largest mean PWV value of approximately 50~60 mm occurs in Southern China in July, so such a relative difference is acceptable. This demonstrates that the impact of T_m error on GPS PWV derivation is substantially lower than that of surface pressure error because Π is less sensitive to T_m error.

To assess the quality of the interpolated T_s, we chose 5.0 K as the threshold value of absolute BIAS and RMSE between the interpolated and measured values based on the above analysis. During these comparisons, a few gross differences between the interpolated and observed temperatures were found. These errors were attributed to the obvious measurement errors. Using the CMONOC station QHBM as an example, abnormal measured temperature values of −167 °C or 145 °C occurred for unknown reasons. Therefore, it is necessary to perform a quality check before performing statistical work. We marked differences with absolute values larger than 20 °C as false and excluded them from the final statistical analysis; 0.003% of the data from ISD and 0.035% of the data from CMONOC were deleted. The corrected statistical results are given in

Table 2. Number of sites at which absolute BIAS or RMSE of interpolated surface air temperature were within a given special value domain.

Value Domain	Absolute BIAS as Statistics		RMSE as Statistics
	ISD	CMONOC	ISD	CMONOC
≤5.0 K	373	232	365	226
>5.0 K	3	8	11	14

and illustrated in Figure 8. Although the accuracy of the GPS-PWV derivation is not directly determined by T_s, as it is also affected by the precision of the T_m–T_s conversion formula, the quality of the interpolated T_s can also be a useful index to indicate the effectiveness of our interpolation scheme. The interpolation results for surface temperature at most stations showed the satisfactory accuracy of near-ground air temperatures for GPS-PWV derivation. However, sites with larger interpolation errors, which were mainly distributed in Southwestern and Northwestern China, also showed large surface temperature interpolation errors. The landforms of these regions included steep undulations, which likely led to poor outcomes. Overall, the accuracy of the interpolated T_s is slightly lower than that of the interpolated P_s, but their influences on the accuracy of the GPS-PWV derivation are of the same order because GPS-PWV is not as sensitive to T_s as to P_s.

4. Comparison of PWV Results

It is a difficult task to evaluate the precision of GPS PWV at every CMONOC station because there is a lack of alternative independent methods to measure PWV more precisely than the GPS measurements, which have a close temporal-spatial resolution. We attempted to use radiosonde data, which can be downloaded freely, to compare with the GPS PWV results. We selected 22 radiosonde stations from IGRA because there were CMONOC GNSS sites very close to them. The geodetic coordinates of the selected radiosonde stations and their nearest GNSS stations are given in

Table 3. Geodetic coordinates of RS sites and their nearby CMONOC GNSS sites.

GNSS Station	lat (°)	lon (°)	Height (m)	RS Station	lat (°)	lon (°)	Height (m)
HLAR	49.27	119.74	627.89	50527	49.22	119.75	611.00
HRBN	45.70	126.62	197.44	50953	45.68	126.62	143.00
JLYJ	42.87	129.50	283.79	54292	42.88	129.47	178.00
LNSY	41.83	123.58	69.02	54342	41.82	123.55	43.00
SDQD	36.08	120.30	12.78	54857	36.07	120.33	77.00
GDST	23.42	116.60	30.93	59316	23.35	116.68	3.00
HIHK	19.99	110.25	54.55	59758	20.03	110.35	24.00
YNMZ	23.36	103.40	1274.50	56985	23.38	103.38	1302.00
KMIN	25.03	102.80	1985.37	56778	25.02	102.68	1892.00
SCGZ	31.61	100.02	3352.28	56146	31.63	99.98	3394.00
GSPL	35.55	106.59	1408.69	53915	35.55	106.67	1348.00
AHAQ	30.62	116.99	57.75	58424	30.52	117.03	20.00
HBES	30.28	109.49	472.38	57447	30.27	109.48	458.00
NMEL	43.63	111.94	945.34	53068	43.65	112.00	966.00
NMEJ	41.96	101.06	888.42	52267	41.98	101.07	941.00
LHAS	29.66	91.10	3623.82	55591	29.70	91.13	3650.00
XZNQ	31.49	92.11	4572.55	55299	31.48	92.05	4508.00
XJYN	43.97	81.53	732.88	51431	43.95	81.33	664.00
XJKC	41.73	82.98	1028.08	51644	41.72	82.95	1100.00
XJAL	47.86	88.13	874.10	51076	47.73	88.08	737.00
GZGY	26.47	106.67	1093.13	57816	26.48	106.65	1222.00
SHAO	31.10	121.20	22.02	58362	31.40	121.47	4.00

. Their locations are shown in Figure 9. These stations are evenly distributed in China, and their altitudes ranged from several meters to higher than 4500 m. Therefore, they represent our derived GPS-PWV values under different conditions.

4.1. Results of Different T_m

We obtained the important parameter, the weighted mean temperature T_m, using three different approaches:

(1)

At all of the RS stations, we directly integrated the radiosonde data assuming that the balloon ascended along a vertical path. The following approximate formula is used:

$$T_m=\frac{{\displaystyle \sum \frac{(z_2-z_1)e}{T}}}{{\displaystyle \sum \frac{(z_2-z_1)e}{T^2}}}$$

(17)

where water vapor pressure e is calculated by e = e_s ×RH, RH represents the relative humidity, and saturation vapor pressure is generated using ITS-90 equations proposed by [37]. We denoted T_m by this method as T_{m_RS}.

(2)

At selected CMONOC stations in Table 3, T_m was estimated from T_s using the T_m–T_s linear equation of Equation (7). From Wang’s research [38], coefficients a and b are determined by the climatic region of every site. T_s was measured by the meteorological sensors collocated to the CMONOC GNSS sites. Then, T_{m_OW} could be obtained.

(3)

The same method as in Equation (2) was employed, except that T_s was obtained from the interpolation schemes described in Section 2. To differentiate their results, we denoted the results from this method as T_{m_IW}.

The former two methods are both based on real observations, whereas the NCEP FNL dataset was employed in the third scheme. The RMSE between the different T_m results are shown in

Table 4. RMSE between the T_m results from three different approaches.

RS Station Number	RMSE (K)
	Tm_OW–Tm_RS	Tm_IW–Tm_RS
50527	3.70	3.64
50953	4.74	4.61
54292	4.20	4.33
54342	3.63	3.69
54857	3.75	3.51
59316	2.25	2.10
59758	2.10	2.29
56985	2.74	3.07
56778	3.10	2.94
56146	2.58	2.28
53915	3.77	3.39
58424	3.33	3.24
57447	3.37	2.96
53068	4.86	4.27
52267	4.83	4.08
55591	3.41	2.58
55299	4.13	4.27
51431	4.01	3.25
51644	3.77	2.78
51076	3.89	3.68
57816	2.15	1.99
58362	3.48	3.66

. In the comparisons, radio sounding was considered to be the most precise measuring method, so the results of the other two approaches were both compared with T_{m_RS}. At each site, the RMSE of both T_{m_OW}s and T_{m_IW} to T_{m_RS} were always below 5 K, and their values were very close at the same site. This indicates that the accuracies of T_{m_OW} and T_{m_IW} are equally acceptable.

In areas without any radiosonde or surface meteorological stations, the accuracy of the NWP data decreased due to the lack of assimilation data. However, there are several CMONOC GNSS sites located in such regions that can provide high-precision surface meteorological observations. Because the accuracy of T_{m_OW} will not change under the above conditions, we can use T_{m_OW} to evaluate T_{m_IW}. We selected 10 CMONOC GNSS stations that met such conditions for further comparisons. The RMSE between the T_{m_IW} and T_{m_OW} of these stations are given in

Table 5. RMSE between T_{m_IW} and T_{m_OW} at CMONOC GNSS stations without nearby meteorological station.

Site Name	RMSE (K)
XZRT (Tibet)	3.29
XZGE (Tibet)	1.44
XZBG (Tibet)	1.76
XZGZ (Tibet)	1.98
WUSH (Xingjiang)	1.95
XJKE (Xingjiang)	1.39
XJDS (Xingjiang)	2.32
NMWL (Inner Mongolia)	2.63
QHWQ (Qinghai)	2.67
SCPZ (Sichuan)	5.69

. The values of the RMSE were still lower than 5 K at most of the selected sites, except SCPZ. This is attributed to, as is previously mentioned, the fact that the SCPZ site had a large height difference to its neighboring NWP grid nodes. Taking this factor out of consideration, we are confident that this dataset exhibits accuracy in terms of the T_{m_IW} that is acceptable, even if there was no precise meteorological data assimilated into the NCEP FNL dataset.

4.2. Results of Different PWV

At the stations given in Table 3 and Figure 9, different PWV results were generated in our study using different derivation methods:

(1)

At the 22 CMONOC GNSS stations, we first employed GAMIT software to process GPS data under ITRF2008. The cut-off elevation angle was 10°, and ZTD was estimated with the GMF mapping function at each station. Because of the long distances between our selected stations, absolute ZTD and atmospheric delay horizontal gradient values, at intervals of one hour and two hours respectively, could be estimated directly without introducing any other GPS sites into this network. Then, we calculated GPS PWV, as described in Section 2.1. Real surface pressure measurements of P_s and T_m were computed from real near-ground air temperature observations T_s. We denoted the PWV results as GPS_PWV_obs.

(2)

This scheme is similar to method (1), expect that P_s and T_s were generated from the interpolation of the NCEP FNL dataset, as proposed in Section 2.2. These GPS PWV results are referred to as GPS_PWV_NCEP.

(3)

PWV can be integrated from the vertical profile of several meteorological parameters using the following formula:

$$PWV={\displaystyle {\int }_0^{p\mathrm{s}}\frac{q}{{\rho }_{w}g}}dp,$$

(18)

where q is the specific humidity, ${\rho }_{\mathrm{w}}$ is the density of liquid water, g is gravitational acceleration, p is air pressure, and p_s represents surface pressure. Equation (18) can be approximated by:

$$PWV={{\displaystyle \sum _i\left(\frac{\overline{q}\Delta p}{{\rho }_{w}g}\right)}}_i.$$

(19)

The essential meteorological parameters can be obtained from the radio soundings. At the radiosonde stations listed in Table 3, the PWV integrated from the radiosonde data is denoted as RS_PWV.

(4)

PWV can also be integrated from NCEP FNL data using the same integral formula as Equation (19) and is referred to as NCEP_PWV.

First, we compared the GPS_PWV_obs and GPS_PWV_NCEP results. The only difference between them is the data source of P_s and T_s; therefore, such a comparison can directly indicate the feasibility of remedying the lack of P_s and T_s measurements in the NCEP FNL data interpolation. There were no meteorological data available at the XZNQ station, so it was excluded from the comparison. Figure 10 illustrates a scatter plot and several statistical and number density plots of differences between GPS_PWV_obs and GPS_PWV_NCEP. Simple linear regression shows that the two results are highly correlated, with a correlation coefficient of 0.9998 and a regression equation of y = 0.9983x − 0.1755, which is very close to y = x. All of the stations have BIAS within ±0.4 mm, and almost all of their absolute values are smaller than 0.2 mm; all of the RMSE values are smaller than 0.63 mm, with most below 0.4 mm. Deviations have a very low value in all of the PWV value fields. These comparisons clearly demonstrate that GPS_PWV_obs and GPS_PWV_NCEP are very similar. Therefore, the accuracy of deriving GPS PWV using P_s and T_s interpolated from NCEP reanalysis data is reliable compared with GPS PWV derived from real surface meteorological data.

We also compared GPS_PWV_NCEP with integrated PWV, including RS_PWV and NCEP_PWV. GPS_PWV_obs was not included in this comparison because its results are very close to those of GPS_PWV_NCEP. Figure 11 shows statistical and number density plots of differences between the three PWV result sets. High correlations exist between GPS_PWV_NCEP and integrated PWV, as demonstrated by the unary linear regressions, with correlation coefficients higher than 0.97. Overall, NCEP_PWV is more strongly correlated to GPS_PWV than RS_PWV. It is common that GPS_PWV_NCEP is more consistent with NCEP_PWV than with RS_PWV in our study. At most stations, the BIAS values between GPS_PWV and NCEP_PWV are within ±2 mm and the RMSEs are approximately 2~3 mm, whereas the BIAS values between GPS_PWV and RS_PWV are within ±3 mm and the RMSEs are approximately 2~4 mm. The number density plot of the differences between the three PWV results indicates that the PWVs were concentrated below 5 mm in our study, and the negative deviation between GPS_PWV and the other two PWVs was also mainly distributed within this PWV range. As the blue line indicates, GPS_PWV shows a dry bias compared with RS. However, this is probably because of the position offsets between the GPS and RS stations. The differences of their horizontal positions or height can also produce some PWV differences.

The higher RMSE at HIHK is attributed to the serious data loss and the very active ionosphere over this low latitude area. As Figure 12 shows, the ratio between L1 and L2 observations increased with increasing L2 data loss, which caused the difference between GPS_PWV_NCEP and NCEP_PWV to increase dramatically. Therefore, it is essential to perform quality control for GPS PWV retrieval in future work [39].

As previously mentioned, the accuracy of the NWP data will decline in regions without sufficient high-precision meteorological observations. Our NCWP_PWV is an indirect and integrated product of the NCEP FNL data; thus, its precision will also decrease. However, the derivation of GPS_PWV_NCEP only requires the interpolation of two levels of NCEP data, so changes in the reanalysis data exert less influence on it. For stations located in regions without intensive radio sounding observations, such as the LHAS and XZNQ stations, the differences between GPS_PWV_NCEP and RS_PWV are smaller than those between GPS_PWV_NCEP and NCEP_PWV.

PWV can also be retrieved by sun photometers [40], and in this paper, we denote this type of PWV product as SP_PWV. Fortunately, there was an active AERONET (Aerosol Robotic Network) site in 2012 located at Taihu, which is in close proximity to the SHAO GNSS station in Table 3. The distance between these stations was approximately 100 km. More details about AERONET and its products are presented on its website. To further assess our PWV products, we compared the GPS_PWV_NCEP of the SHAO site, RS_PWV of the 58362 radiosonde site and SP_PWV of the Taihu AERONET site. Figure 13 is a scatter plot of the three PWV datasets. Common observation times between GPS_PWV and SP_PWV occurred more frequently than between RS_PWV and SP_PWV. Statistics are given in

Table 6. Statistics of differences between SP_PWV, GPS_PWV_NCEP, and RS_PWV.

	BIAS(mm)	RMSE(mm)	Correlation
GPS_PWVNCEP vs. SP_PWV	0.3464	3.3439	0.9870
RS_PWV vs. SP_PWV	−1.7908	3.8969	0.9763

. Clearly, GPS_PWV_NCEP is highly consistent with SP_PWV, with a bias of only 0.3464 mm and RMSE of 3.3439 mm, and the correlation coefficient was higher than 0.98. Their agreement was even better than that of RS_PWV and SP_PWV.

The accuracy of PWV retrieved using our scheme under extreme weather conditions also required verification. Figure 14 illustrates the time series of RS_PWV and GPS_PWV_NCEP at the SHAO station from UTC time 00:00 28 July to 24:00 9 August in 2012. The time interval of GPS_PWV_NCEP was 6 h, whereas that of RS_PWV was 12 h. The SHAO site is located in Shanghai city. During this period, Shanghai suffered from two severe typhoons, Damrey and Haikui, in rapid succession. However, the two PWV time series still agreed very well, with a bias of −1.05 mm and RMSE of 2.83 mm. It can be seen that the PWV value fluctuated very dramatically in a short time period. Two rapid continuous PWV increase processes, that is, from 33 mm to 72 mm between 01.08.2012 06:00 and 02.08.2012 06:00 and from 47 mm to 73 mm between 06.08.2012 18:00 and 08.08.2012 12:00, were both reflected in the GPS_PWV_NCEP and RS_PWV time series.

5. Conclusions and Outlook

Employing as many GPS sites as possible to obtain GPS PWV can provide higher spatial and temporal resolution PWV data than other traditional methods, such as radio soundings or GPS precipitable water vapor products from limited GPS stations that equipped meteorological sensors. High spatiotemporal resolution water vapor information is beneficial for studying the details of vapor properties at a specific site location, or during a particular weather event. However, there are many geodetic GPS stations that do not have collocated meteorological sensors. Using the NWP dataset can be an alternative to obtaining the required meteorological parameters necessary for the traditional GPS PWV derivation. In this paper, we studied the scheme of interpolating NCEP FNL reanalysis data to generate surface pressure or temperature over China. Comparisons with real meteorological measurements show the accuracies of the interpolated surface pressure or surface temperature are acceptable for the GPS PWV derivation over most Chinese regions, except for some strong relief areas, even under particular synoptic cases, such as typhoons. Compared with radio sounding data, the weighted mean temperature can also be calculated using the surface temperature interpolated from the NCEP FNL dataset. Our study indicates that GPS PWV products that use interpolated surface meteorological elements were very close to those based on real meteorological observations, with differences within ±0.4 mm and a RMSE mainly below 0.6 mm. These GPS PWV results also strongly agreed with radio sounding observations, as demonstrated by their high correlations and low differences. The main factor influencing the accuracy of GPS PWV is GPS data processing. In our study, serious loss of GPS L2 data reduced the accuracy of the estimated GPS ZTD dramatically, and led to a further decrease of GPS PWV accuracy. However, the feasibility of applying such GPS PWV solutions to research long-term trends of water vapor over certain region requires further investigation. For example, experiments should be carried out over decades to confirm the long-term accuracy of our interpolation method. Another key issue is that numerous historical GPS data should be re-processed, which may involve many problems such as the difference between relative calibration and absolute calibration of antenna phase center variations [41].

Because interpolation of NWP data can be used to generate surface meteorological parameters, we demonstrate that it can serve as a very promising complement to the present ground-based GPS meteorological observation networks under appropriate conditions, such as in relatively flat regions. The time duration and spatial range of GPS PWV should be extended through future research endeavors. Details of some interesting climatic phenomena or extreme weather conditions can consequently be studied better. Another relevant research topic is to calculate essential meteorological elements for GPS PWV using forecast data instead of reanalysis data. This can help forecasters develop real-time, high-resolution, and high-precision precipitable water vapor maps. Such products are reconstructed much faster than data assimilation of the NWP system. If the accuracies are found to be suitable, these products would be helpful for rainstorm warnings and forecasting of flood hazards.