The dense ground-based GPS provides a good tool to study water vapor distribution and multiscale variations, especially for linear trends on the interannual scale and short-term variations on the diurnal scale. It can also serve as an independent data source to evaluate performances of reanalyses. In this study, the 6-hourly precipitable water (PW) products at more than 260 GPS stations over China from 1999 to 2015 were analyzed and eight commonly used reanalyses, including 20CR version 2 (20CRv2), CFSR, ERA-Interim, JRA-25, JRA-55, MERRA, NCEP–NCAR, and NCEP–DOE AMIP-II, were evaluated. The climatological annual mean GPS PW distribution over China roughly shows a decreasing trend from southeast to northwest, with the largest annual and semiannual amplitudes in the lower reaches of the Yangtze River and mideastern China, respectively, and the smallest values in the Tibetan Plateau and southwestern China. All reanalyses (except for 20CRv2) can generally reproduce well the climatological annual mean PW (within 20%), annual amplitudes (within 20%), and semiannual amplitudes (within 20% except in the tropical monsoon region), but they all show wet biases in the Tibetan Plateau. Diurnal variation amplitudes reproduced by all reanalysis products are smaller than amplitudes estimated from GPS observations over China as a whole, and none of the reanalyses can capture the diurnal phases correctly. PW linear trends at most GPS stations in the recent 16 years are insignificant or with absolute values smaller than 0.10 mm yr−1. However, because of the assimilation of the unhomogenized radiosonde humidity data, most reanalyses show artificial decreasing PW trends (except in 20CRv2 and CFSR).
Water vapor, which is the most abundant greenhouse gas in the air, plays an important role in both the energy and hydrological cycle on Earth (Kiehl and Trenberth 1997). According to the Clausius–Clapeyron equation, a 1-K increase of temperature will induce around 7% increase of water vapor holding capacity of the air, and the increase of water vapor, as a greenhouse gas, will in turn enhance the warming. This strong positive feedback has notable influences on the climate system (Held and Soden 2000). Because of the high specific heat of water, fast variations of water vapor content associated with large changes of the latent heat always accompany the emergence, development, and vanishing of some extreme weather phenomena (e.g., Allan and Soden 2008). Therefore, monitoring and study of water vapor content in the air are indispensable to understand Earth’s climate and weather processes.
Traditionally, because they have the longest records and near-global coverage over the land, radiosondes are widely used to study the distribution, variations, and trends of water vapor content on both global and regional scales (e.g., Elliott and Gaffen 1991; Zhai and Eskridge 1996, 1997; Ross and Elliott 2001; Trenberth et al. 2005). However, the radiosonde measurements suffer from the well-known inhomogeneity issue, which is mainly caused by changes in the instrument type, the observation practice, and the analysis strategy (Garand et al. 1992; Zhai and Eskridge 1996). The inhomogeneity issue significantly limits applications of radiosonde records in the water vapor climatological study, especially in long-term trends study. In addition, radiosonde measurements are usually made twice daily, which makes it hard to capture diurnal and even subdiurnal variations of water vapor content.
Meteorological reanalysis products provide comprehensive global multidecadal records of historical atmosphere states using a single consistent numerical data assimilation scheme with various past observations. They have been used extensively in climate and weather research and applications. Although unified processing strategies and models are used in the reanalysis, biases or inhomogeneities existing in the assimilated data can still affect the reanalysis product. The reliability of reanalysis products is, therefore, a matter of great concern. Regarding the humidity fields, Trenberth et al. (2005) pointed out that both NCEP reanalyses and ERA-40 suffer from substantial problems over the ocean in terms of precipitable water (PW). Paltridge et al. (2009) and Dessler and Davis (2010) found that decadal trends of specific humidity in NCEP–NCAR reanalysis data are unrealistic, especially in the tropical middle and upper troposphere. At shorter time scales, verifications of NCEP reanalyses showed that global numerical weather prediction (NWP) models have difficulties in handling the diurnal cycle and moist processes on the synoptic scale (Dessler and Davis 2010).
The ground-based global positioning system (GPS) can utilize the delay of the radio signal transmitting from satellite to ground-based receiver caused by neutral atmosphere to infer the water vapor content in the atmosphere. Since 1990s, ground-based GPS has been proved to be a powerful tool in measuring water vapor content (Bevis et al. 1992) with advantages of low-cost, high accuracy, high temporal resolution, and all-weather operations. GPS data have not been assimilated into any of the reanalysis products, which provides an ideal way to study multiscale water vapor content variations on multiple time scales and to evaluate reanalyses as an independent data source.
China is located at the junction of the largest continental plate (the Eurasian plate) and the largest ocean (the Pacific Ocean). Distribution and variations of water vapor content over China are significantly affected by the land–sea discrepancy as well as terrain characteristics and climate types. As early as 1953, Xu (1958) used data from 33 radiosonde stations to calculate water vapor content in eastern China for January and July in 1957. Lu and Gao (1984) studied the climatological distribution of water vapor over China based on 1960–70 monthly mean radiosonde records. However, Zhai and Eskridge (1996) found that the radiosonde records used in Lu and Gao (1984) contain inhomogeneities caused by changes in sensor types and observation time. Later on, Zhai and Eskridge (1997) used relatively homogeneous radiosonde data from 1970 to 1990 to examine the water vapor climatology, trends, and variations over China. They concluded that the climatological distribution of water vapor over China is primarily dependent on latitude, topographical features, and monsoon conditions. Interannual variations and trends in water vapor content and surface temperature are closely correlated. Some of Chinese radiosonde stations have also been used in several previous studies of Northern Hemisphere water vapor trends (e.g., Ross and Elliott 2001; Durre et al. 2009; McCarthy et al. 2009). Recently, Zhao et al. (2012) analyzed the variations and trends of tropospheric humidity over China based on radiosonde dataset homogenized using the method proposed by Dai et al. (2011) and found that the PW increase between the surface and 300 hPa of about 2%–5% decade−1 over most of China in recent 40 years. Based on 28 GPS stations from 2004 to 2007, Jin et al. (2008) quantified the distribution and multiscale variations of PW over China, and the results show that the distribution and variations are mainly affected by latitude, topographical features, seasons, and monsoon conditions. A small number of GPS stations in China are also used in several published global water vapor studies (e.g., Wang and Zhang 2008; Vey et al. 2009).
Regarding the humidity fields in reanalyses over China, recently, Bao and Zhang (2013) evaluated NCEP–NCAR, CFSR, ERA-40, and ERA-Interim products using 11 radiosonde stations over the Tibetan Plateau and found that all the reanalysis products can produce temperature and horizontal winds well but relative humidity poorly. Lu et al. (2015) found that none of the commonly used reanalyses, including JRA-55, ERA-40, ERA-Interim, MERRA, NCEP-2, and ISCCP, are able to show the actual long-term trends and variations in PW over Tibetan Plateau. Based on the homogenized radiosonde dataset, evaluations show that few reanalyses can capture the observed long-term PW changes correctly, primarily resulting from the assimilation of inhomogeneous radiosonde humidity measurements (e.g., Zhao et al. 2015).
Because of the vital role of GPS in diverse research fields, GPS station networks around the world have been dramatically expanded in recent years, which provide invaluable dataset for water vapor study. The Crustal Movement Observation Network of China (CMONOC) includes a national and well-maintained GPS network that has been in operation since the late 1990s, but this dataset has not been explored extensively for water vapor study. In this paper, GPS PW products derived from CMONOC stations covering the period from 1999 to 2015 will be used to study the climatological distribution, multiscale variations, and linear trends of PW over China and to evaluate the commonly used global reanalysis products. The GPS and reanalysis products as well as the PW retrieval method used in this paper will be described in section 2. The PW variation extraction method and linear trend estimation method will be depicted in section 3. The distribution, multiscale variations, and linear trends of GPS PW over China from 1999 to 2015 will be examined and evaluations of reanalysis products will be carried out in section 4, followed by discussions and conclusions in section 5.
2. Data and methodology
a. GPS PW products
The 6-hourly GPS PW products were generated based on CMONOC GPS stations covering the period from 1 March 1999 to 30 April 2015. Details regarding the GPS data processing and PW retrievals can be found in Zhang et al. (2017). Here we will only give a brief description.
CMONOC project contains two phases. From 1999 to mid-2010, which is called phase I, there are about 28 permanent GPS stations over China (denoted by triangles in Fig. 1). Starting from the end of 2010 (phase II), the network was expanded and the number of GPS stations dramatically increased to about 265 by 2015 (Fig. 1). GPS data with a sampling interval of 30 s were processed by the precise point positioning (PPP) method (Zumberge et al. 1997) in a daily manner using the position and navigation data analyst (PANDA) software package (Shi et al. 2008) developed at Wuhan University.
The reprocessed satellite orbit and clock products provided by Massachusetts Institute of Technology (MIT), which are in the reference frame called IGS08 (Rebischung et al. 2012), were used to avoid inconsistences in satellite products regarding the reference frame. The absolute antenna phase center correction model (Schmid et al. 2007), the phase wind-up corrections (Wu et al. 1993), and the relativity corrections were applied in the GPS data analysis. Station coordinates were estimated as daily constants and receiver clock errors were solved every epoch regardless of temporal correlations. The zenith tropospheric delay (ZTD) can be divided into two components, namely the zenith hydrological delay (ZHD) and the nonhydrological delay, which is generally referred to as the zenith wet delay (ZWD). The empirical global pressure and temperature (GPT) model (Boehm et al. 2007) and relative humidity of 50% at mean sea level were used to estimate the a priori ZHD and a priori ZWD based on the Saastamoinen model (Saastamoinen 1972). Corrections to the a priori ZWD, which can also assimilate any mismodeling of the a priori ZHD, were then estimated as piecewise constant every 2 h with a power density of 20 mm h−1/2. The sum of the a priori ZHD, the a priori ZWD, and the estimated ZWD corrections was taken as the final ZTD. The global mapping function (GMF) (Boehm et al. 2006) was applied and tropospheric gradients in the north–south and east–west directions were estimated every 12 h to account for the azimuthally asymmetric of the atmosphere. Cutoff elevation angles of satellites were set to 7° and an elevation-dependent weighting strategy was applied to the measurements at low elevations (below 30°) (Gendt et al. 2003) to reduce the influences of multipath effects and the large uncertainties of mapping function at low elevations.
Air pressures at GPS stations were obtained from GPS collocated meteorological sensor measurements, nearby synoptic meteorological station records, or ERA-Interim products (Zhang et al. 2017). Precise values for ZHD were then calculated using the combined pressures and were subtracted from ZTD to derive ZWD. ZWD values were finally converted to PW based on water vapor weighted temperature Tm estimated from ERA-Interim. Error budget analysis showed that the uncertainty of the generated PW products is about 0.75 mm (Zhang et al. 2017).
b. Reanalysis products and PW retrievals
Eight global reanalysis products will be evaluated in this study, including the Twentieth Century Reanalysis, version 2c (20CRv2c), (Compo et al. 2011), Climate Forecast System Reanalysis (CFSR) (Saha et al. 2010), ERA-Interim (Dee et al. 2011), JRA-25 (Onogi et al. 2007), JRA-55 (Ebita et al. 2011), MERRA (Rienecker et al. 2011), NCEP–NCAR reanalysis (Kistler et al. 2001), and NCEP–DOE AMIP-II reanalysis (Kanamitsu et al. 2002). Main features of theses eight reanalysis products are summarized in Table 1.
The 6-hourly reanalysis grid products at pressure levels covering the same period as GPS data, namely from 1 March 1999 to 30 April 2015 (to 31 December 2014 for 20CRv2c and to 31 January 2014 for JRA-25), were used to derived PW at GPS stations. CFSR products were initially updated only to 31 December 2010, but then on 30 March 2011 NCEP upgraded their operational CFS to version 2 (CFSv2), which used the same model as CFSR to extend CFSR dataset continuously from 1 January 2011 to the present. Therefore, CFSR products used in this work are actually the combination of CFSR and CFSv2. Fields (pressure, geopotential, temperature, and specific humidity) at the four grid points nearest the GPS station are horizontally interpolated to the GPS station location at each pressure level. PW can be then calculated by integrating specific humidity from the GPS antenna level to the top level using
where ρw is the liquid water density; q denotes the specific humidity (kg kg−1), and P is the atmospheric pressure. The integral is carried out over height, taking into account of the difference between the pressure reference to mean sea level and the height reference to ellipsoidal height and considering that gravitational acceleration g varies as a function of height.
For GPS antenna levels above the lowermost reanalysis level, temperature and specific humidity at the antenna level are simply estimated by linear interpolation, while for GPS antenna levels below the lowermost reanalysis level the temperature at the antenna level is estimated by linear extrapolation, assuming a constant temperature lapse rate of −6.5 K km−1. The hydrostatic and ideal gas equations are used to adjust pressure to the GPS antenna height as described in Wang et al. (2007). Regarding the estimation of specific humidity for the antenna below the lowermost reanalysis level, the relative humidity at the lowermost reanalysis level is first calculated from the specific humidity, temperature, and pressure. Equivalent relative humidity between the GPS antenna and the lowermost reanalysis level is then assumed, and the specific humidity at the GPS antenna level is finally estimated on the basis of the relative humidity, temperature, and pressure. Detailed procedures for the specific humidity estimation are provided in the appendix.
3. Multiscale PW variation extraction method
a. PW fitting model
The power spectral density (PSD) of PW time series at each GPS station is estimated using the Welch method (Welch 1967) and is shown in Fig. 2. We can find that the most notable periodic signals are mainly annual and semiannual variation signals, which is consistent with conclusions in Jin et al. (2008) based on 28 stations over China with 3-yr PW time series. Accordingly, the following regression model is used to fit PW time series at each GPS station (Llamedo et al. 2017):
where y(t) denotes PW at time t (yr); A0 and A1 are the constant and trend coefficients; and A2, P1, A3, and P2 represent the annual amplitude, annual phase, semiannual amplitude, and semiannual phase, respectively. QBO1 and QBO2 are standardized time series representing the quasi-biennial oscillation (QBO) variations at 30- and 50-hPa height, respectively (http://www.cpc.ncep.noaa.gov/data/indices), and SF denotes the standardized time series representing the solar flux variations (http://www.esrl.noaa.gov/psd/data/correlation/solar.data). QBO1, QBO2, and SF are used to account for the influences of QBO and solar flux variations on the PW variations, and B1, B2, and B3 are the corresponding coefficients. The least squares method is used to estimate the coefficients (A0, A1, A2, P1, A3, P2, B1, B2, and B3) based on Eq. (2) at each GPS station. To avoid unrealistic estimates, stations with PW samples less than 600 will be excluded.
The geographic distribution of the fitting residual root-mean-square (RMS) is presented in Fig. 3, where we can find obvious correlations between RMS and geographic locations. Generally, RMS values are larger in eastern China than in western China, and the largest are found in the lower reaches of Yangtze River with typical values larger than 9 mm while the smallest are in the Tibetan Plateau (smaller than 3 mm). The average RMS at all stations over China is about 4.9 mm. Four stations [Fangshan (BJFS), Wuhan (WUHN), Urumqi (URUM), and Lhasa (LHAZ)] in different regions are selected and comparisons of raw PW time series and the modeled PW time series at these four stations are presented in Fig. 4. It is obvious that the residual RMS values are strongly correlated to PW values; that is, for stations BJFS and WUHN with larger PW, the RMS are larger than for stations URUM and LHAZ with relatively smaller PW, and the residuals usually maximize in summer when PW is also maximum during a year.
b. Diurnal variation extraction
In addition to the seasonal variations (annual and semiannual variations), the PW time series also exhibit a pronounced diurnal cycle although its amplitudes are generally one order of magnitude smaller than seasonal variations (Wang and Zhang 2009). Diurnal variations of water vapor are usually related to many weather processes, such as the diurnal variations of moist convection and precipitation, surface wind convergence, and surface evapotranspiration (Dai et al. 2002). Although they have considerable influences on the surface and atmospheric longwave radiation as well as the atmospheric absorption of solar radiation, the diurnal variations of water vapor are traditionally poorly studied due to the lack of high temporal resolution data, which can be improved by GPS observations (Wang and Zhang 2009).
In the derivations of PW products, pressure is obtained from GPS collocated sensors, nearby synoptic stations, or reanalysis products, with sampling intervals of 30 s, 3 h, and 6 h, respectively. The water vapor weighted temperatures are calculated from reanalysis with 6-h intervals. Although the ZTDs are with 2-h intervals, to avoid introducing additional diurnal errors into PW, no temporal interpolations were carried out for pressures and water vapor weighted temperatures. The 6-hourly PW products may not be sufficient enough to capture the accurate diurnal cycle amplitudes and phases, but they can still provide us with a general picture of PW diurnal variations over China. It is also fair to compare GPS PW diurnal variations with reanalysis products (generally with 6-h intervals for moisture fields) with the same temporal interval.
For the diurnal variation extraction, the daily mean PW is calculated and removed from the PW time series to get the PW diurnal anomalies from which the diurnal variation amplitudes and phases can then be analyzed.
c. Linear trend estimation method
For PW linear trend estimation, instead of directly using A1 in Eq. (2), we use Sen’s nonparametric method (Sen 1968), which is commonly adopted for trend estimations by the meteorological community. Sen’s method provides a more robust slope estimate than the least squares method since it is insensitive to outliers or extreme values (Fan and Yao 2003), particularly for relatively small datasets, and it can also reduce the impacts of autocorrelation in time series.
The monthly PW anomalies are first calculated as the deviations of the monthly PW from the monthly PW climatology (mean value of monthly PW over the whole period from 1999 to 2015). The 5-point-moving-average method is then applied to the monthly PW anomaly time series, and PW trends are then estimated as
4. PW distribution, variations, and trends
a. PW differences between GPS and reanalysis products
According to the climate features, the climate over China can be divided into five types (China Meteorological Administration 1979), namely the tropical monsoon (TM), subtropical monsoon (SM), midlatitude monsoon (MM), humid continental (HC), and plateau (PL) climate types. The CMONOC GPS stations located in regions with different climate types are indicated by different colors in Fig. 1. In this section, the reanalysis products are first evaluated by calculating PW differences between GPS and reanalysis products. The mean value and standard deviation (STD) of PW differences in different climate type regions covering the period from 1 March 1999 to 30 April 2015 (31 December 2014 for 20CRv2 and 31 January 2014 for JRA-25) are presented in Fig. 5. Hereinafter we will use the shortened acronyms 20CR, ERAI, NCEP1, and NCEP2 to denote 20CRv2, ERA-Interim, NCEP–NCAR, and NCEP–DOE AMIP-II, respectively. Among these eight reanalyses, we find that 20CR generally has the largest biases (always positive; i.e., wetter than GPS) and STD. This is related to the fact that humidity fields in 20CR over the land are not directly constrained by observations (Compo et al. 2011), which is different from the other seven reanalyses. The bias can reach about 4 mm in the MM region with STD of about 5 mm. On the other hand, ERAI agrees with GPS PW best with unobvious bias and STD of about 2.1 mm over China as a whole. The improvements of the newer generation of reanalysis (e.g., JRA-55) relative to the older one (e.g., JRA-25) provided by the same organization (i.e., JMA) are not obvious. Among five different climate type regions, all reanalyses show more abundant water vapor than GPS in the PL region, which is relevant to the well-known wet bias (significantly more precipitation than observations) existing in current global climate models (e.g., Gao et al. 2015). Based on nine GPS stations in southern Tibetan Plateau, Wang et al. (2017) also found similar evident PW wet biases (0.72–1.49 mm) in MERRA, ERAI, JRA-55, and NCEP final analysis. However, we also noticed that Liu et al. (2005) found an underestimation of PW in ECMWF analysis compared to three GPS stations in the Tibetan Plateau during the period from 2002 to 2003. This discrepancy in Liu et al. (2005) and our results may be due to the different GPS PW retrieval procedures, which needs further investigation. The STDs of absolute PW differences are roughly proportional to the mean PW (as will be shown later) in different regions for all reanalyses; that is, STDs are largest in the tropical monsoon region with the moistest air while smallest in the driest plateau region.
b. Climatological annual means
Average value of the modeled PW based on Eq. (2) from 1 January 1999 to 31 December 2015 is considered as the climatological annual mean PW at each GPS station. Figure 6 presents the geographic distribution of GPS annual mean PW as well as the differences of annual means between GPS and reanalyses. As can be seen from Fig. 6a, water vapor content in southeastern China is much more abundant than in western China, roughly showing a decreasing trend from southeast to northwest. The moistest place is found in southeastern China with mean PW ranging from 45 to 60 mm, while the driest region is located in the Tibetan Plateau with mean PW less than 15 mm due to the high elevation and the snow-covered surface, resulting in a colder atmosphere and thus less water. The average values of mean PW in different climate type regions are summarized in Table 2. We can find that the average climatological mean PW is more than 35 mm in the tropical monsoon region, which is more than 5 times the values in the plateau region. The average PW in the humid continental region is less than 10 mm, especially in northwestern China due to the natural barrier of the plateau and mountains to the south, which block the transport of moist air from the Bay of Bengal (Zhai and Eskridge 1997).
Regarding the relative differences of mean PW (relative to GPS) between GPS and reanalyses as shown in Figs. 6b–i and Table 2, all reanalyses agree best with GPS in the tropical monsoon or subtropical monsoon regions but worst in the Tibetan Plateau. Compared with other reanalyses, 20CR is distinctive in reproducing larger PW than GPS at most stations with relative differences of 20%–100%, except for stations in southeastern region (with unobvious biases). The other seven reanalyses can generally reproduce annual mean PW well, with relative differences from −20% to 20% at most stations (CFSR is slightly larger). All of ERAI, JRA-25, JRA-55, MERRA, NCEP1, and NCEP2 show smaller mean PW than GPS in the south of the middle and lower reaches of Yangtze River but larger mean PW in the region between Yangtze and Yellow Rivers. From Fig. 6 and Table 2, we can also easily find that all reanalyses overestimate mean PW of about 25%–60% in the Tibetan Plateau region. This finding is inconsistent with Zhao et al. (2015), who found underestimations of 40%–60% for mean PW relative bias (relative to radiosonde PW) over the western Tibetan Plateau in nine commonly used reanalysis products. The main reason is that Zhao et al. (2015) used radiosonde stations with sparse distribution over this area, and the PW values at grid points with horizontal resolution of 1° × 1° were interpolated from surrounding stations with lower elevations (thus with a deeper atmosphere and higher PW), resulting in artificially higher radiosonde PW over the Tibetan Plateau. This indicates that GPS data can be good complements to radiosonde observations in regions with sparse distributions (such as in the Tibetan Plateau) for PW climatological studies, especially as the number of GPS stations has dramatically increased in recent years.
c. Annual and semiannual variations
The geographic distribution of PW annual amplitudes [A2 in Eq. (2)] is presented in Fig. 7. Annual amplitudes in regions with monsoon climate types (>15 mm) are generally larger than regions with humidity continental (9.4 mm) and plateau climate types (7 mm). This is also evident in Fig. 8 where the seasonal variations of average climatological daily mean PW in five climate type regions are depicted. The most dramatic PW annual variations are found in the SM region as seen from Fig. 8 and Table 3, especially in the lower reaches of Yangtze River (with annual amplitudes larger than 25 mm) (as shown in Fig. 7). This is mainly due to the influences of monsoons. In the summer, the monsoon circulation (including the East Asian monsoon and the Indian monsoon) continuously brings a large amount of water vapor from the western Pacific Ocean, the South China Sea, and the Bay of Bengal (e.g., Ninomiya 1999; Ding and Sun 2001) to the monsoon region. Because of the encounter and stalemate of the cold and warm air masses, a quasi-stationary water vapor band is formed over the lower reaches of Yangtze River, which is well known as the mei-yu front (Lin 1979), resulting in continuously abundant water vapor content lasting for several weeks. Therefore, we can find the climatological daily mean PW in the SM region (near the lower reaches of Yangtze River) is almost equivalent to the values in the TM region (closer to the origins of the monsoon) in summer during days of year from about 170 to 230 as shown in Fig. 8. However, in winter the SM region is dominated by the winter monsoon with dry and cold air from the Siberia in the north (e.g., Huang et al. 2003), resulting in relatively low water vapor content (~10 mm). The smallest annual amplitudes are found in the PL region due to the dry and cold weather all the year around there as can be seen from the time series in Fig. 8.
All eight reanalyses can depict annual amplitudes fairly well with relative differences smaller than 20% as summarized in Table 3. 20CR (Fig. 7b) presents interesting spatial distribution characteristics where we can find negative differences in both the southeast and the northwest while positive differences between these two regions. This may indicate some spatial related systematic humidity errors existing in 20CR over China. CFSR agrees with GPS best overall, with average relative differences smaller than 10% in all subregions. Both ERAI and MERRA can reproduce PW annual amplitudes with small differences (from −5% to 5%) in the eastern half of China, but the differences can exceed 20% in the western part. For JRA-25 and JRA-55, annual amplitudes are larger than GPS at most stations and, surprisingly, the older reanalysis (JRA-25) has better performances than the newer generation of reanalysis (JRA-55) regarding the PW annual amplitudes. NCEP2 also reproduces larger amplitudes than GPS at most stations in China and almost in the same amplitudes with NCEP1 except in the Xinjiang Province located in northwestern China.
Compared to annual amplitudes, semiannual amplitudes are much smaller as shown in Fig. 9a and Table 4, with largest values of about 6–8 mm in mideastern China and smallest values less than 2 mm in southwestern China (the region with latitudes from 20° to 30°N and longitudes from 100° to 110°E). Peixoto et al. (1981) pointed out that the water vapor content in spring tends to resemble more closely to the low levels of winter, while the water vapor content in autumn is more closely approach the peaks of summer in the Northern Hemisphere due to the influences of the ocean’s thermal inertia. In the MM region where the largest semiannual amplitudes are found, for example, the peak PW value appears around day 200 as shown in Fig. 8. The day when PW reaches 20 and 30 mm from spring to summer is around days 150 and 170, respectively, namely about 50 and 30 days apart from day 200. However, from summer to autumn, the time when PW reaches 20 and 30 mm is about 66 and 38 days, respectively, apart from the day with the maximum PW value as shown in Fig. 8. The existence of the asymmetry in the annual cycle of PW is the main reason causing semiannual variations of PW. The distribution of semiannual amplitudes as presented in Fig. 9 also provides a good indicator to show influences of the ocean on moisture conditions in the atmosphere.
Relative differences of semiannual amplitudes of 20CR (Fig. 9b) show different distribution characteristics compared to the annual amplitudes (Fig. 7b), with systematic transitions from southwest to northeast (positive differences to negative and then back to positive). Overall, the performance of MERRA is best among these reanalyses regarding semiannual amplitudes. ERAI, JRA-25, NCEP1, and NCEP2 reproduce similar results as shown in Fig. 9 and the most notable differences with GPS are found in southwestern China (the region with smallest semiannual amplitudes) with negative relative differences of 30%–40%. Similar to the annual amplitudes, again, JRA-55 underestimates the semiannual amplitudes at most stations as shown in Fig. 9f.
d. Diurnal variations
The seasonal variations of PW diurnal anomalies over China and five subregions with different climate types are presented in Fig. 10. The diurnal amplitudes in the HC region (<0.5 mm) (with higher latitudes) are smaller than in other regions. Wang and Zhang (2009) also found that PW diurnal cycle amplitudes are generally smaller at higher latitudes based on International Global Navigation Satellite System (GNSS) Service (IGS) GPS stations. Compared to the HC region, the diurnal variations are evidently stronger in the plateau region (>1.0 mm in summer) as shown in Fig. 10e although these two regions have comparable mean PW values, which is consistent with conclusions in Jin et al. (2008). The diurnal amplitudes are generally largest in late spring and early summer in all regions. Regarding the diurnal phase, the time of PW minimum value varies with season as shown in Fig. 10. The minimum time is usually at 0800 Beijing time (BJT; =UTC + 8 h) in spring and summer (from March to August) but makes a transition to 0400 BJT in autumn (from September to November) except for the plateau region (near 1200 BJT), and then inverses and transforms to 1200–1400 BJT in winter (from December to February). On the other hand, the seasonal variations of PW peaks are not obvious (generally peaking at 2000 BJT) in all regions as shown in Fig. 10. Several previous studies, however–for example, by Dai et al. (2002) over North America, by Jin et al. (2008) over China, and by Wang and Zhang (2009) over the globe—found that PW generally peaks from noon to late evening, showing some seasonal variations. This discrepancy may be due to the relatively larger sampling intervals (6 h) of PW products we used in this work, compared to 30 min in Dai et al. (2002) and 2 h in Jin et al. (2008) and Wang and Zhang (2009).
For brevity, only the comparisons of PW diurnal variations over China as a whole between GPS and reanalysis products are depicted in Fig. 11. Results for individual subregion are provided in the supplemental material (Figs. S1–S5). Apparently, the diurnal variation amplitudes reproduced by all reanalysis products are smaller (less than 0.5 mm) than amplitudes estimated from GPS observations over China as a whole. For amplitudes in individual subregion, CFSR is the closest one to GPS in the TM and SM regions (Figs. S1 and S2), JRA-55 performs best in the MM (Fig. S3) region, and ERAI agrees well with GPS in the PL region (Fig. S5). Regarding the phases, the time of maximum PW for CFSR, JRA-25, and JRA-55 is almost invariant during one year (at 2000 BJT) as shown in Fig. 11, which is similar to GPS, whereas the peak PW time in all the other reanalyses changes with season. However, none of reanalyses can depict the seasonal variations of minimum PW time correctly. Among these five subregions, the largest phase differences are found in the HC region for all reanalyses (Fig. S4). In the TM region, 20CR is the only one showing diurnal variations with phase opposite to GPS in summer and autumn (Fig. S1), which is again partly due to the fact that no humidity observations were assimilated in 20CR.
e. Linear trends
PW linear trends at GPS stations are estimated using method described in section 3c based on PW monthly anomaly time series derived from GPS observations at stations operated since phase I (Fig. 12a) and eight reanalyses at all GPS stations from 1 March 1999 to 30 April 2015 (to 31 December 2014 for 20CR and to 31 January 2014 for JRA-25) (Figs. 12b–i). The absences of about 1.25 and 0.33 years of data for JRA-25 and 20CR are neglected in the trend comparisons. The correlation coefficient between reanalysis and GPS PW trends is also indicated in each panel. From Fig. 12a, we can find that PW linear trends at most GPS stations in recent 16 years are insignificant or with absolute values smaller than 0.10 mm yr−1. Only two stations show significant trends larger than 0.10 mm yr−1 in southeastern China. As for the reanalysis products, it is obvious that all reanalysis products, except 20CR and CFSR, reproduce significant decreasing trends at most stations, especially in the eastern half of China. This is mainly because of the assimilation of unhomogenized radiosonde humidity records in the reanalyses. The previous evaluations (Wang and Zhang 2008; Zhao et al. 2012; Zhang et al. 2017) show that the older radiosonde type GZZ2, which was widely used in China before the 2006, contains evident wet biases, and was then replaced gradually by the GTS1 model, with obvious dry biases, or by the latest generation of radiosondes, GTS1-1 and GTS1-2, with negligible biases. This will produce artificial decreasing PW trends for most radiosonde stations during the period from 1999 to 2015, and therefore result in negative correlation coefficients between GPS and reanalysis PW trends for ERAI, JRA-25, JRA-55, MERRA, NCEP1, and NCEP2. The spurious trends of PW in reanalysis products can project to trends of other related fields, such as precipitation (Gao et al. 2015). Among the eight reanalyses, 20CR performs best with correlation coefficient of about 0.07, which is partly due to the fact that 20CR assimilates only surface observations of synoptic pressure, monthly SST, and sea ice distribution (Compo et al. 2011), which avoids the inhomogeneity issue existing in radiosonde humidity data. 20CR presents significant decreasing trends between 0.10 and 0.15 mm yr−1 in southeastern China and northwest of Xinjiang Province while increasing trends at most stations in northeastern China. CFSR shows considerably different PW trend distributions compared to both GPS and the other reanalyses, with increasing trends at most stations, especially for stations located at 100°–110°E. This may be partly because that the CFSR is a fully coupled land–ocean–atmosphere reanalysis while all of the other reanalyses are forced with observed SSTs.
The actual climatic trends of PW are not simply linear and they can change dramatically over different periods. In this work, only the linear components of PW trends were estimated and compared, and readers should keep it in mind that estimations and comparisons in this section may not be valid for other different periods or different regions.
5. Discussion and conclusions
With the accumulation of GPS historical datasets and increasing number of operational GPS stations, the denser GPS station observations provide a valuable opportunity to study the water vapor climatology in regions with sparse radiosonde stations. In addition, GPS data can also serve as an independent data source to evaluate the performances of moisture fields in reanalysis products. Today reanalysis products are widely used in both climate and weather research. However, the reliabilities of different fields in reanalysis products need to be carefully evaluated before proper applications. This is the first attempt to evaluate, over China, the performances of eight commonly used reanalyses, 20CRv2c, CFSR, ERA-Interim, JRA-25, JRA-55, MERRA, NCEP–NCAR, and NCEP–DOE AMIP-II, in reproducing multiscale characteristics of PW based on a national GPS network over China covering more than 16 years.
PW differences between GPS and reanalyses show that 20CR generally has largest biases (wet biases) and STD among eight reanalyses that can reach 4 and 5 mm, respectively, in the midlatitude monsoon region. ERAI agrees with GPS PW best with negligible biases and STD of about 2.1 mm in China as a whole. However, all reanalyses reproduce more abundant water vapor than GPS in the Tibetan Plateau.
GPS climatological annual mean PW distribution roughly shows a decreasing trend from southeast to northwest, with more than 35 mm on average in the tropical monsoon region, which is more than 5 times the values in the plateau region. Relative differences of mean PW between reanalyses and GPS are small in the tropical monsoon region and subtropical monsoon region, but largest in the plateau region. All reanalyses except for 20CR can generally reproduce annual mean PW well, with relative differences smaller than 20% at most stations.
Annual amplitudes in regions with monsoon climate types (>15 mm) are generally larger than regions with humidity continental (9.4 mm) and plateau climate types (7 mm). On the other hand, the largest semiannual amplitudes are found in mideastern China of about 6–8 mm and smallest values (<2 mm) in southwestern China. All eight reanalyses can depict annual amplitudes fairly well with relative differences smaller than 20%. Overall, the performances of CFSR and MERRA are the best among these reanalyses regarding the annual and semiannual amplitudes.
The diurnal variation amplitudes in the humid continental (HC) region (<0.5 mm) are smaller than amplitudes in other regions based on GPS results. The PW generally peaks at 2000 BJT based on the 6-hourly PW products and the minimum time varies with season. Diurnal variation amplitudes reproduced by all reanalysis products are smaller (less than 0.5 mm) than amplitudes estimated from GPS observations over the whole of China and none of the reanalyses can capture the diurnal phases correctly.
PW linear trends at most GPS stations in recent 16 years are insignificant or with absolute values smaller than 0.10 mm yr−1. Only two stations show significant increasing trends larger than 0.10 mm yr−1 in southeastern China. However, because of the assimilation of unhomogenized radiosonde humidity data, most reanalyses show artificial decreasing PW trends except for 20CR and CFSR.
Results and conclusions in this work provide a more complete picture of the water vapor climatology over China in the last 16 years on a national scale as there were more GPS stations than radiosonde stations, especially in the Tibetan Plateau. It also gives us an idea about the reliabilities of water vapor derived from the widely used global reanalysis products over China on multiple time scales. The evaluation results of the current reanalysis products derived from this work can be potentially incorporated into the developments of next generation of reanalyses in the future.
The NCEP–NCAR, NCEP–DOE AMIP-II, 20CR, CFSR, JRA-55, and JRA-25 data were obtained from the CISL Research Data Archive (http://rda.ucar.edu/) at NCAR in Boulder, Colorado. The ERA-Interim data were obtained from the ECMWF Data Server (http://apps.ecmwf.int/datasets/). The MERRA data were obtained from the Goddard Earth Sciences Data and Information Services Center, Greenbelt, Maryland, from their website at http://disc.sci.gsfc.nasa.gov/mdisc. This work was supported by the National Key Research and Development Program of China (2016YFB0501800), the National Natural Science Foundation of China (41774036), and the Fundamental Research Funds for the Central Universities (2042017kf0179).
Specific Humidity Estimation at GPS Antenna Levels
For GPS antenna levels below the lowermost reanalysis level, the relative humidity (RH) at the lowermost reanalysis level is first estimated by
where es denotes the saturation water vapor pressure (hPa) that can be calculated by
where T is temperature (K). Also, e in Eq. (A1) denotes the water vapor pressure (hPa) that can be derived from specific humidity q (kg kg−1) as
where P denotes pressure (hPa).
After getting the temperature and pressure at the GPS antenna level following the method as described in section 2b, we can assume the equivalent relative humidity between the lowermost reanalysis level and the GPS antenna. The saturation water vapor pressure at the GPS antenna level is calculated on the basis of temperature through Eq. (A2), and the water vapor pressure can be estimated as
and the specific humidity at the GPS antenna level can be finally calculated by
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-17-0419.s1.