Hourly PWV Dataset Derived from GNSS Observations in China

The rapid variation of atmospheric water vapor is important for a regional hydrologic cycle and climate change. However, it is rarely investigated in China, due to the lack of a precipitable water vapor (PWV) dataset with high temporal resolution. Therefore, this study focuses on the generation of an hourly PWV dataset using Global Navigation Satellite System (GNSS) observations derived from the Crustal Movement Observation Network of China. The zenith total delay parameters estimated by GAMIT/GLOBK software are used and validated with an average root mean square (RMS) error of 4–5 mm. The pressure (P) and temperature (T) parameters used to calculate the zenith hydrostatic delay (ZHD) and weighted average temperature of atmospheric water vapor (Tm) are derived from the fifth-generation reanalysis dataset of the European Centre for Medium-Range Weather Forecasting (ECMWF ERA5) products. The values of P and T at the GNSS stations are obtained by interpolation in the horizontal and vertical directions using empirical formulas. Tm is calculated at the GNSS stations using the improved global pressure and temperature 2 wet (IGPT2w) model in China with an RMS of 3.32 K. The interpolated P and T are validated by interpolating the grid-based ERA5 data into radiosonde stations. The average RMS and bias of P and T in China are 2.71/−1.11 hPa and 1.88/−0.51 K, respectively. Therefore, the error in PWV with a theoretical RMS of 1.85 mm over the period of 2011–2017 in China can be obtained. Finally, the hourly PWV dataset in China is generated and the practical accuracy of the generated PWV dataset is validated using the corresponding AERONET and radiosonde data at specific stations. Numerical results reveal that the average RMS values of the PWV dataset in the four geographical regions of China are less than 3 mm. A case analysis of the PWV diurnal variations as a response to the EI Niño event of 2015–2016 is performed. Results indicate the capability of the hourly PWV dataset of monitoring the rapid water vapor changes in China.


Introduction
Atmospheric water vapor is a key component of the troposphere, and it affects the global/regional water cycles and climate change [1]. Water vapor is also an important greenhouse gas, which is highly related to the heat feedback function. Therefore, knowledge of rapid water vapor variations is important for analyzing global/regional water vapor distributions [2]. However, the distribution of atmospheric water vapor with high temporal-spatial resolutions is still difficult to monitor, due to the shortcomings of traditional techniques for water vapor monitoring. Precipitable water vapor (PWV), which is used to reflect the content of atmospheric water vapor, is the value corresponding to the conversion of all atmospheric water vapor into liquid water in the air column of a unit area from the surface to the troposphere high. A number of techniques have been developed to monitor PWV, and they include radio sounding, water vapor radiometer (WVR), solar photometer, and remote

Interpolation Method
As mentioned previously, the grid-based P and T data derived from ERA5 must be interpolated into the RS stations prior to comparison; the interpolation can be completed on the basis of the following formula [26]: ERA-Interim and ERA5 are the fourth-and fifth-generation reanalysis datasets of the ECMWF that cover the reanalysis data from 1979 to two to three months before the present. Both datasets provide meteorological parameters, such as P, T, and PWV, with different temporal-spatial resolutions. The temporal resolution of the ERA-Interim dataset is 6 h, and that of the ERA5 dataset is 1 h. Therefore, the ERA5 dataset is expected to be widely used in the future and gradually replace the ERA-Interim dataset. Compared with the previous generations of data products, the ERA5 reanalysis dataset was generated by using the 4D-Var data assimilation in IFS Cycle 41R2 (CY41R2) of the ECMWF Integrated Prediction System. The parameters P and T are derived from the ERA5 dataset and interpolated into the GNSS station and height with a temporal resolution of 1 h (https://www.ecmwf.int/). Radiosonde data are also selected in the current work to validate the accuracy of P and T derived from ERA5 and PWV using GNSS and ERA5 data. The radiosonde data are derived from the Integrated Global Radiosonde Archive Version 2 (IGRA2) dataset. IGRA2 was generated in 2016 by the National Climate Data Center and contains more radiosonde stations than the previous dataset (IGRA1). This dataset includes P, T, relative humidity, and other parameters with temporal resolutions of twice or four times daily (ftp://ftp.ncdc.noaa.gov/pub/data/igra/). In this experiment, 87 radiosonde stations in China are selected (Figure 1), and 52 of them are collocated with the corresponding GNSS stations.
The PWV derived from AErosol RObotic NETwork (AERONET) with a temporal resolution of 1 h is also used to validate the generated PWV dataset at four stations during the period 2014-2017. AERONET, which is composed of a CE-318 solar photometer, mainly studies the optical and microphysical properties of global aerosols. The CE-318 solar photometer makes direct spectral measurements of solar radiation in channels between 340 and 1640 nm and in a water vapor band (936 nm). Multichannel measurements are used to retrieve aerosol optical thickness, and measurements at 936 nm are used to calculate PWV [21,22]. In the past 20 years, solar photometer-based PWV has been applied to many studies due to its high reliability and accuracy [23,24]. The three levels of PWV data provided by AERONET are raw data (Level 1.0), cloud-screened data (Level 1.5), and cloud-screened and quality-assured data (Level 2.0), respectively [25]. The corresponding PWV dataset can be obtained from the AERONET website (https://aeronet.gsfc.nasa.gov/). In this study, the PWV values derived from the Level 2.0 products are selected at the chosen four stations ( Figure 1). And we arranged the characteristics of the datasets selected for the whole experiment in Table 1. Here, China is divided into four regions according to the characteristics of geographical location, natural geography, and human geography.

Interpolation Method
As mentioned previously, the grid-based P and T data derived from ERA5 must be interpolated into the RS stations prior to comparison; the interpolation can be completed on the basis of the following formula [26]: where P 0 refers to the pressure at the grid point and h and h 0 are the heights at the station and grid point, respectively. dT and dh are the differences of the T and height between station and grid point, respectively. The height system should be unified due to the different height systems used for the ECMWF, RS, and GNSS data. The data provided by ECMWF and IGRA adopt geopotential height systems, whereas the GNSS data adopt the ellipsoidal height system. In this work, the normal height system is used. The specific steps to convert normal height from geopotential height to quasi-geoid height were presented in our previous study [27]. An empirical correction formula of PWV is also used herein to unify the PWV values of the collocated stations at different heights where PWV h1 and PWV h2 are the PWV values corresponding to the heights of h1 and h2, respectively.

Retrieval of PWV
In this work, ZTD is derived from CMONOC. ZHD at GNSS stations can be calculated using the interpolated P on the basis of Equation (3) [9]. Therefore, the ZWD used to calculate PWV can be obtained by extracting ZHD from ZTD. The ERA5-derived T at the GNSS stations are used to calculate the T m value on the basis of the IGPT2w model proposed by Huang et al. (2019a) [14]. PWV is then obtained with Equation (4): where ϕ and H are the latitude (rad) and orthometric height (km) of a GNSS station, respectively. ZHD is mainly affected by P, and 1 hPa of error in P leads to a 0.2 mm error in PWV [28]. PWV = 10 6 where K 2 , K 3 , and R V are constants with values of 16.48 K·hPa −1 , (3.776 ± 0.014) × 10 5 K 2 ·hPa −1 , and 461 J·kg −1 ·K −1 , respectively. ρ is the water vapor density [29].

Comparison between ERA5 and ERA-Interim Products
The ERA5-derived P and T products at grid points with temporal-spatial resolutions of 1 h and 0.5 • × 0.5 • , respectively, are first compared with those derived from ERA-Interim over China in the period of 2005-2017. China covers a vast territory with different characteristics in the east, west, north, and south regions. Hence, the country is divided into four geographical regions, namely, the northern, southern, northwest, and Qinghai Tibet regions, according to geographical location and the natural and human geographical characteristics of different regions ( Figure 1). Figure 2 presents the standard error (STD) and bias of P and T between ERA5 and ERA-Interim in China during the period 2005-2007. The accuracies of P and T derived from those two types of re-analysis data are similar, but the differences of P and T are relatively large in the northwest and Qinghai Tibet regions. This deviation may be related to the complex terrain and the lack of adequate ground-based meteorological observations. The statistical results in Table 2 reveal that the average STD and bias of P and T in China are 0.71/−0.36 hPa and 1.77/−0.33 K, respectively. Therefore, it can be concluded that the P and T data of ERA5 and ERA-Interim in most regions of China have good consistency. Because the time resolution of ERA5 data is higher than the ERA-interim data, the ERA5 data set is likely to replace ERA-interim data and widely used in the future.

Comparison of ERA5 Products with Radiosonde Data
The P and T derived from ERA5 are used to interpolate the corresponding values at the RS station and height during the period 2005-2017 in China. The interpolated P and T at the RS stations are validated with the RS-derived P and T, respectively. In the experiment, the differences of the ERA5-/RS-derived P and T larger than 3 are regarded as gross errors and thus removed. The accuracies of P and T derived from ERA5 can subsequently be calculated. Figure 3 presents the root mean square (RMS) and bias of P and T derived from ERA5 at 87 RS stations during the period 2005-2017 in China. Table 3 provides the specific statistical results of a number of RS stations and the average RMS and bias of the ERA5-derived P and T in the four regions of China. It can be found that the ERA5-derived P and T in northern and southern China are of good quality, whereas the accuracies of the values in the northwest and Qinghai Tibet regions are slightly poor. The statistical results reveal that the average RMS and bias of P and T in China are 2.71/−1.11 hPa and 1.88/−0.51 K, respectively. Such results are consistent with those in Section 3.1 and are thus enough to verify the good performance of the ERA5-derived P and T in China.

Comparison of ERA5 Products with Radiosonde Data
The P and T derived from ERA5 are used to interpolate the corresponding values at the RS station and height during the period 2005-2017 in China. The interpolated P and T at the RS stations are validated with the RS-derived P and T, respectively. In the experiment, the differences of the ERA5-/RS-derived P and T larger than 3σ are regarded as gross errors and thus removed. The accuracies of P and T derived from ERA5 can subsequently be calculated. Figure 3 presents the root mean square (RMS) and bias of P and T derived from ERA5 at 87 RS stations during the period 2005-2017 in China. Table 3 provides the specific statistical results of a number of RS stations and the average RMS and bias of the ERA5-derived P and T in the four regions of China. It can be found that the ERA5-derived P and T in northern and southern China are of good quality, whereas the accuracies of the values in the northwest and Qinghai Tibet regions are slightly poor. The statistical results reveal that the average RMS and bias of P and T in China are 2.71/−1.11 hPa and 1.88/−0.51 K, respectively. Such results are consistent with those in Section 3.1 and are thus enough to verify the good performance of the ERA5-derived P and T in China.

Analysis of GNSS-Derived ZTD from CMONOC
To validate the consistency of the GNSS-derived ZTD obtained on the basis of GAMIT/GLOBK (version 10.4) software with that processed by other softwares, the corresponding ZTD parameters calculated with the PPP technique using the positioning and navigation data analyst (PANDA) package are applied [30,31]. PANDA software has been independently developed by the research center of satellite navigation and positioning technology of Wuhan University since 2001, and has been adopted by many famous international research institutions. Its goal is to realize the postprocessing and real-time processing analysis of GNSS/SLR/DORIS/VLBI and other kinds of observation data [32].The double difference mode is used in the GAMMIT/GLOBK software while the non-difference mode is applied in the PANDA package. A total of 246 out of 260 GNSS stations are selected during the period 2011-2017 in China. The ZTD differences larger than 3 are removed as gross errors. Figure 4 presents the STD and bias distributions of the GNSS-derived ZTD differences using the GAMIT/GLOBK and PANDA software during the period 2011-2017. Table 4 provides the

Analysis of GNSS-Derived ZTD from CMONOC
To validate the consistency of the GNSS-derived ZTD obtained on the basis of GAMIT/GLOBK (version 10.4) software with that processed by other softwares, the corresponding ZTD parameters calculated with the PPP technique using the positioning and navigation data analyst (PANDA) package are applied [30,31]. PANDA software has been independently developed by the research center of satellite navigation and positioning technology of Wuhan University since 2001, and has been adopted by many famous international research institutions. Its goal is to realize the post-processing and real-time processing analysis of GNSS/SLR/DORIS/VLBI and other kinds of observation data [32].The double difference mode is used in the GAMMIT/GLOBK software while the non-difference mode is applied in the PANDA package. A total of 246 out of 260 GNSS stations are selected during the period 2011-2017 in China. The ZTD differences larger than 3σ are removed as gross errors. Figure 4 presents the STD and bias distributions of the GNSS-derived ZTD differences using the GAMIT/GLOBK and PANDA software during the period 2011-2017. Table 4 provides the statistical results of the STD and bias, as well as the number of RS stations, in four regions of China. As shown in Figure 4, the STD and  Table 4 indicates that the ZTD error is slightly large in southern China with an RMS of 5.5 mm; this region is mainly affected by the large atmospheric water vapor content. The RMS value is approximately 3.7 mm in northwest China, where the atmospheric water vapor is low. statistical results of the STD and bias, as well as the number of RS stations, in four regions of China. As shown in Figure 4, the STD and bias in mainland China are considerably small, and the statistical results reveal that their values are 4.6 and −0.4 mm, respectively. Table 4 indicates that the ZTD error is slightly large in southern China with an RMS of 5.5 mm; this region is mainly affected by the large atmospheric water vapor content. The RMS value is approximately 3.7 mm in northwest China, where the atmospheric water vapor is low.

Theoretical Error of PWV Calculated Using the Hourly P and T from ERA5
The theoretical error of PWV can be deduced using the error propagation law according to the accuracies of P, ZTD, and Tm. In our experiment, the accuracies of P, ZTD, and T are approximately 2.7 hPa, 5 mm, and 1.9 K, respectively. Therefore, the theoretical error of PWV in this work can be calculated on the basis of the following formula [33]: where

Theoretical Error of PWV Calculated Using the Hourly P and T from ERA5
The theoretical error of PWV can be deduced using the error propagation law according to the accuracies of P, ZTD, and T m . In our experiment, the accuracies of P, ZTD, and T are approximately 2.7 hPa, 5 mm, and 1.9 K, respectively. Therefore, the theoretical error of PWV in this work can be calculated on the basis of the following formula [33]: where σ v , σ ZTD , σ P 0 , and σ Q represent the errors in the ZTD, P, T, and conversion factor, respectively. V refers to PWV, and σ C = 0.0024, Q is a parameter for converting ZWD to PWV with an average value of 0.16. f (λ, H) accounts for the variation of gravitational acceleration at the ellipsoidal latitude ϕ and ellipsoidal height H of the GNSS station; it is generally considered to be equal to 1. In addition, σ Q is the effect of T m , which can be neglected for most cases. Therefore, the final theoretical error of PWV can be calculated using Equation (6) based on the obtained errors in P, T, ZTD and conversion factor. Table 5 shows the statistical mean values of the four regions. It can be found that the average theoretical error of PWV in the whole of China is about 1.85 mm. The maximum value occurs in the Qinghai-Tibet area, with the value being approximately 2 mm; whereas the minimum value is noted in northern China, with the value being about 1.3 mm.

T m Calculation Using Improved GPT2w (IGPT2w) Model
As described in Equation (2), T m is the only parameter affecting the accuracy of PWV when ZWD is known. Therefore, how to obtain a high precision T m using the interpolated T at a GNSS station is the key point in this section. The GPT2w model is one of the most widely used models for calculating T m , but it is prone to evident systematic bias [34]. The main reason is that the vertical correction of T m is not considered in the GPT2w model. To overcome this issue, Huang et al. (2019) proposed an IGPT2w model, which considers the height correction of T m ; such a model has greatly corrected the systematic error in the western region of China relative to the GPT2w model [18]. The corresponding vertical correction of T m is presented in Equation (7). Therefore, the IGPT2w model is used in this work to calculate the T m value of each GNSS station. Figure 5 presents the RMS distributions of the T m differences derived from the GPT2w and IGPT2w models at 87 radiosonde stations when the RS-derived T m is considered as the reference. The accuracy of the IGPT2w-derived T m in China, especially the southern region, shows an improvement. The statistical result reveals that the average improvement rate of T m using the IGPT2w model is approximately 10% in China. Figure 6 gives the comparison experiment of T m calculated by Bevis formula and IGPT2w model. It can be seen from the figure that the accuracy of T m derived from IGPT2w is higher than that from Bevis formula, so we chose the IGPT2w model to calculate T m in this paper.  (7) where T t m is the T m calculated by the GPT2w model at the target height while T r m is the T m calculated at the height of the grid point. δh t and δh r are the heights at the target and grid point in km, respectively. doy is day of year. α 0 represents the mean value of the decreasing rate of T m . (α 1 , α 2 ) and (α 3 , α 4 ) are the coefficients of the annual and semiannual periods, respectively. T t m is the T m calculated by the GPT2w model at the target height.

Comparison of Hourly PWV Dataset with AERONET Data
After the Tm is calculated, the hourly PWV dataset at each GNSS station can be obtained on the basis of the GNSS-derived ZTD and ERA5 data during the period 2011-2017. To validate the hourly PWV dataset, this study selects and compares four collocated stations between AERONET and GNSS. Figure 7 presents the time series of the hourly PWV derived from AERONET and GNSS at the four collocated stations. It can be observed from Figure 7 that the GNSS-derived PWV agrees well with that from AERONET at the four collocated stations. The calculated correlation coefficients between the GNSS-and AERONET-derived PWV are 0.97, 0.99, 0.99, and 0.93. A high correlation indicates that the generated hourly PWV dataset is of good quality. The statistical results reveal that the RMS values of the GNSS-derived PWV at the Beijing, Sacol, Xianghe, and Qoms_cas stations are 1.13, 2.04, 2.04, and 1.41 mm, respectively.

Comparison of Hourly PWV Dataset with AERONET Data
After the Tm is calculated, the hourly PWV dataset at each GNSS station can be obtained on the basis of the GNSS-derived ZTD and ERA5 data during the period 2011-2017. To validate the hourly PWV dataset, this study selects and compares four collocated stations between AERONET and GNSS. Figure 7 presents the time series of the hourly PWV derived from AERONET and GNSS at the four collocated stations. It can be observed from Figure 7 that the GNSS-derived PWV agrees well with that from AERONET at the four collocated stations. The calculated correlation coefficients between the GNSS-and AERONET-derived PWV are 0.97, 0.99, 0.99, and 0.93. A high correlation indicates that the generated hourly PWV dataset is of good quality. The statistical results reveal that the RMS values of the GNSS-derived PWV at the Beijing, Sacol, Xianghe, and Qoms_cas stations are 1.13, 2.04, 2.04, and 1.41 mm, respectively.

Comparison of Hourly PWV Dataset with AERONET Data
After the T m is calculated, the hourly PWV dataset at each GNSS station can be obtained on the basis of the GNSS-derived ZTD and ERA5 data during the period 2011-2017. To validate the hourly PWV dataset, this study selects and compares four collocated stations between AERONET and GNSS. Figure 7 presents the time series of the hourly PWV derived from AERONET and GNSS at the four collocated stations. It can be observed from Figure 7 that the GNSS-derived PWV agrees well with that from AERONET at the four collocated stations. The calculated correlation coefficients between the GNSS-and AERONET-derived PWV are 0.97, 0.99, 0.99, and 0.93. A high correlation indicates that the generated hourly PWV dataset is of good quality.

Comparison of Hourly PWV Dataset with RS Data
To further evaluate the accuracy of the hourly PWV dataset in China, this study selects the GNSS and radiosonde data from 52 out of the 87 collocated stations in China during the period 2011-2017. Herein, the horizontal distance less than 60 km in the latitudinal and longitudinal directions and the vertical difference between the GNSS and radiosonde stations less than 500 m are used as a principle to judge whether the radiosonde and GNSS stations are collocated or not. The empirical formula in Equation (3) is used to unify the PWV value from radiosonde height to GNSS height. Figure 8 presents the time series of the GNSS-and RS-derived PWV at the HLAR station (49.3° N, 119.7° E) during the period 2011-2017. A good consistency exists between the GNSS-and RS-derived PWV. Figure 9 presents the RMS and bias distributions of the PWV differences derived from GNSS and RS at 52 stations during the period 2011-2017. The statistical result reveals that the average RMS and bias are 2.25 and 1.57 mm, respectively. Table 6 provides the specific RMS and bias in the four regions of China. Excluding that of the Qinghai-Tibet region, the RMS values in the other three areas are less than 3 mm. Such a result obtained above indicates the potential application of the hourly PWV dataset generated in this work in weather nowcasting [35]. The average RMS is larger than the theoretical value, and this phenomenon is acceptable because of the uncertainties in radiosonde PWV data and the error introduced when using the empirical formula in Equation (3) [11].

Comparison of Hourly PWV Dataset with RS Data
To further evaluate the accuracy of the hourly PWV dataset in China, this study selects the GNSS and radiosonde data from 52 out of the 87 collocated stations in China during the period 2011-2017. Herein, the horizontal distance less than 60 km in the latitudinal and longitudinal directions and the vertical difference between the GNSS and radiosonde stations less than 500 m are used as a principle to judge whether the radiosonde and GNSS stations are collocated or not. The empirical formula in Equation (3) is used to unify the PWV value from radiosonde height to GNSS height. Figure 8 presents the time series of the GNSS-and RS-derived PWV at the HLAR station (49.3 • N, 119.7 • E) during the period 2011-2017. A good consistency exists between the GNSS-and RS-derived PWV. Figure 9 presents the RMS and bias distributions of the PWV differences derived from GNSS and RS at 52 stations during the period 2011-2017. The statistical result reveals that the average RMS and bias are 2.25 and 1.57 mm, respectively. Table 6 provides the specific RMS and bias in the four regions of China. Excluding that of the Qinghai-Tibet region, the RMS values in the other three areas are less than 3 mm. Such a result obtained above indicates the potential application of the hourly PWV dataset generated in this work in weather nowcasting [35]. The average RMS is larger than the theoretical value, and this phenomenon is acceptable because of the uncertainties in radiosonde PWV data and the error introduced when using the empirical formula in Equation (3) [11]. Sensors 2019, 19, x FOR PEER REVIEW 12 of 16

Comparison of PWV Image with ERA5
Apart from the comparison with the GNSS-derived PWV at collocated stations, the PWV image is also compared with that from the ERA5. A total of 246 GNSS stations are used to construct the Delaunay triangulation net in China, and the hourly PWVs of three points of a Delaunay triangle are used to interpolate the PWV value at the grid point with the spatial resolution of 0.5° × 0.5°. Figure  10 provides the mean PWV image derived from GNSS and ERA5 in different seasons in 2016 and their difference distribution. It can be observed that the PWV distribution derived from GNSS agrees well with that from ERA5 under different seasons, but the PWV difference is slightly larger in

Comparison of PWV Image with ERA5
Apart from the comparison with the GNSS-derived PWV at collocated stations, the PWV image is also compared with that from the ERA5. A total of 246 GNSS stations are used to construct the Delaunay triangulation net in China, and the hourly PWVs of three points of a Delaunay triangle are used to interpolate the PWV value at the grid point with the spatial resolution of 0.5° × 0.5°. Figure  10 provides the mean PWV image derived from GNSS and ERA5 in different seasons in 2016 and their difference distribution. It can be observed that the PWV distribution derived from GNSS agrees well with that from ERA5 under different seasons, but the PWV difference is slightly larger in

Comparison of PWV Image with ERA5
Apart from the comparison with the GNSS-derived PWV at collocated stations, the PWV image is also compared with that from the ERA5. A total of 246 GNSS stations are used to construct the Delaunay triangulation net in China, and the hourly PWVs of three points of a Delaunay triangle are used to interpolate the PWV value at the grid point with the spatial resolution of 0.5 • × 0.5 • . Figure 10 provides the mean PWV image derived from GNSS and ERA5 in different seasons in 2016 and their difference distribution. It can be observed that the PWV distribution derived from GNSS agrees well with that from ERA5 under different seasons, but the PWV difference is slightly larger in summer. The statistical results reveal that 91%, 81%, 90%, and 87% of the areas in China have PWV differences of less than 5 mm in four seasons, respectively.
summer. The statistical results reveal that 91%, 81%, 90%, and 87% of the areas in China have PWV differences of less than 5 mm in four seasons, respectively.

Analysis of Diurnal PWV Variations in China
Hourly PWV is generated in this work to calculate the diurnal variations of atmospheric water vapor, which is an important indicator for climate monitoring. Therefore, the diurnal anomalies of PWV in the four regions of China for four seasons are calculated. The hourly diurnal anomaly of PWV is calculated by removing the average PWV of the corresponding season. Spring, summer, autumn, and winter are from March to May, June to August, September to November, and December to February, respectively. Figure 11 presents the hourly diurnal anomalies of the PWV time series in the four regions of China under different seasons during the period 2015-2016. It can be found that the diurnal anomalies of the PWV in four seasons in 2016 were below and around 0 mm and increased to above 0 mm to a different extent gradually. In terms of season, this phenomenon is especially obvious in summer and autumn. Regionally speaking, the phenomenon is most obvious in the southern region. The main reason is the El Niño/Southern Oscillation (ENSO) event in 2016-2017. Heavy rainfall occurred in southern China, especially in summer and autumn. This condition led to a decrease in atmospheric water vapor content in the troposphere. Therefore, the hourly diurnal anomaly of PWV shows a decreasing trend with values below and around 0 mm. In 2017, the atmospheric water vapor returned to its normal average level, as reflected by the increasing diurnal anomaly of PWV.

Analysis of Diurnal PWV Variations in China
Hourly PWV is generated in this work to calculate the diurnal variations of atmospheric water vapor, which is an important indicator for climate monitoring. Therefore, the diurnal anomalies of PWV in the four regions of China for four seasons are calculated. The hourly diurnal anomaly of PWV is calculated by removing the average PWV of the corresponding season. Spring, summer, autumn, and winter are from March to May, June to August, September to November, and December to February, respectively. Figure 11 presents the hourly diurnal anomalies of the PWV time series in the four regions of China under different seasons during the period 2015-2016. It can be found that the diurnal anomalies of the PWV in four seasons in 2016 were below and around 0 mm and increased to above 0 mm to a different extent gradually. In terms of season, this phenomenon is especially obvious in summer and autumn. Regionally speaking, the phenomenon is most obvious in the southern region. The main reason is the El Niño/Southern Oscillation (ENSO) event in 2016-2017. Heavy rainfall occurred in southern China, especially in summer and autumn. This condition led to a decrease in atmospheric water vapor content in the troposphere. Therefore, the hourly diurnal anomaly of PWV shows a decreasing trend with values below and around 0 mm. In 2017, the atmospheric water vapor returned to its normal average level, as reflected by the increasing diurnal anomaly of PWV.

Conclusions
An hourly PWV dataset is generated at 246 GNSS stations during the period 2011-2017 in China using the GNSS observations from CMONOC. The meteorological parameters at the GNSS station and height are interpolated at the horizontal and vertical directions using the ERA5-derived hourly P and T data with temporal-spatial resolutions of 0.5° × 0.5° and 1 h, respectively. China is divided into four regions, and the quality of the ERA5-derived P and T is evaluated using the ERA-Interim and radiosonde data during the period 2006-2017. The comparison with 84 radiosonde stations reveals that the average RMS and bias of P/T in China are 2.71/−1.11 hPa and 1.88/−0.51 K, respectively. GNSS-derived ZTDs at 246 stations from CMONOC are analyzed on the basis of two processing modes, and the STD and bias of ZTD in China are about 4.7 and −0.47 mm, respectively. Therefore, the theoretical error of PWV caused by errors in ZTD, P, and T is approximately 1.9 mm according to the law of error propagation. A practical validation of the hourly PWV dataset at 264 GNSS stations is performed, and the results are compared with the AERONET-/RS-derived PWV data at collocated stations. The outcomes show that the average RMS and bias of the hourly PWV dataset are 2.25 and 1.57 mm, respectively. In addition, the PWV image is compared with that from ERA-Interim under different seasons in 2016. The comparison highlights the good performance of the generated PWV

Conclusions
An hourly PWV dataset is generated at 246 GNSS stations during the period 2011-2017 in China using the GNSS observations from CMONOC. The meteorological parameters at the GNSS station and height are interpolated at the horizontal and vertical directions using the ERA5-derived hourly P and T data with temporal-spatial resolutions of 0.5 • × 0.5 • and 1 h, respectively. China is divided into four regions, and the quality of the ERA5-derived P and T is evaluated using the ERA-Interim and radiosonde data during the period 2006-2017. The comparison with 84 radiosonde stations reveals that the average RMS and bias of P/T in China are 2.71/−1.11 hPa and 1.88/−0.51 K, respectively. GNSS-derived ZTDs at 246 stations from CMONOC are analyzed on the basis of two processing modes, and the STD and bias of ZTD in China are about 4.7 and −0.47 mm, respectively. Therefore, the theoretical error of PWV caused by errors in ZTD, P, and T is approximately 1.9 mm according to the law of error propagation. A practical validation of the hourly PWV dataset at 264 GNSS stations is performed, and the results are compared with the AERONET-/RS-derived PWV data at collocated stations. The outcomes show that the average RMS and bias of the hourly PWV dataset are 2.25 and 1.57 mm, respectively. In addition, the PWV image is compared with that from ERA-Interim under different seasons in 2016. The comparison highlights the good performance of the generated PWV dataset in the whole of China. The diurnal variation of the PWV time series is also analyzed under different seasons in 2016 and 2017, and the result indicates that the hourly PWV dataset is capable of reflecting the response of atmospheric water vapor to climate events, such as the ENSO event.
Author Contributions: Q.Z., P.Y., and W.Y. conceived the study and designed the review. Q.Z., and P.Y. conducted the literature review. Q.Z. wrote the first draft. Y.Y. performed the data collection. All authors have read and agreed to the published version of the manuscript.